g4simulation/g4detectors/PHG4HcalSubsystem.cc

0001 #include "PHG4HcalSubsystem.h"
0002
0003 #include "PHG4HcalDetector.h"
0004 #include "PHG4HcalSteppingAction.h"
0005
0006 #include <g4main/PHG4HitContainer.h>
0007 #include <g4main/PHG4Utils.h>
0008
0009 #include <phool/PHCompositeNode.h>
0010 #include <phool/PHIODataNode.h>    // for PHIODataNode
0011 #include <phool/PHNode.h>          // for PHNode
0012 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0013 #include <phool/PHObject.h>        // for PHObject
0014 #include <phool/getClass.h>
0015
0016 #include <Geant4/G4String.hh>         // for G4String
0017 #include <Geant4/G4SystemOfUnits.hh>  // for cm
0018 #include <Geant4/G4Types.hh>          // for G4double
0019
0020 #include <cmath>     // for asin, cos, sin, sqrt, M_PI
0021 #include <cstdlib>   // for NULL, exit
0022 #include <iostream>  // for operator<<, basic_ostream
0023 #include <sstream>
0024
0025 class PHG4Detector;
0026 class PHG4SteppingAction;
0027
0028 //_______________________________________________________________________
0029 PHG4HcalSubsystem::PHG4HcalSubsystem(const std::string& na, const int lyr)
0030   : detector_(nullptr)
0031   , steppingAction_(nullptr)
0032   , radius(100)
0033   , length(100)
0034   , xpos(0)
0035   , ypos(0)
0036   , zpos(0)
0037   , lengthViaRapidityCoverage(true)
0038   , TrackerThickness(100)
0039   , material("G4_Fe")
0040   , _sciTilt(0)
0041   , _sciWidth(0.7)
0042   , _sciNum(256)
0043   , _sciPhi0(0)
0044   , active(0)
0045   , absorberactive(0)
0046   , layer(lyr)
0047   , detector_type(na)
0048   , superdetector("NONE")
0049   , light_scint_model_(true)
0050   , light_balance_(false)
0051   , light_balance_inner_radius_(0.0 * cm)
0052   , light_balance_inner_corr_(1.0)
0053   , light_balance_outer_radius_(10.0 * cm)
0054   , light_balance_outer_corr_(1.0)
0055 {
0056   // put the layer into the name so we get unique names
0057   // for multiple SVX layers
0058   // std::ostringstream nam;
0059   // nam << na << "_" << lyr;
0060   Name(na + "_" + std::to_string(lyr));
0061 }
0062
0063 //_______________________________________________________________________
0064 int PHG4HcalSubsystem::InitRun(PHCompositeNode* topNode)
0065 {
0066   // create hit list only for active layers
0067   PHNodeIterator iter(topNode);
0068   PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0069   // create detector
0070   detector_ = new PHG4HcalDetector(this, topNode, Name(), layer);
0071   detector_->SetRadius(radius);
0072   G4double detlength = length;
0073   if (lengthViaRapidityCoverage)
0074   {
0075     detlength = PHG4Utils::GetLengthForRapidityCoverage(radius + TrackerThickness) * 2;
0076   }
0077   detector_->SetLength(detlength);
0078   detector_->SetPosition(xpos, ypos, zpos);
0079   detector_->SetThickness(TrackerThickness);
0080   detector_->SetTilt(_sciTilt);
0081   detector_->SetScintWidth(_sciWidth);
0082   detector_->SetNumScint(_sciNum);
0083   detector_->SetScintPhi0(_sciPhi0);
0084   detector_->SetMaterial(material);
0085   detector_->SetActive(active);
0086   detector_->SetAbsorberActive(absorberactive);
0087   detector_->SuperDetector(superdetector);
0088   detector_->OverlapCheck(CheckOverlap());
0089   if (active)
0090   {
0091     std::string nodename;
0092     if (superdetector != "NONE")
0093     {
0094       nodename = "G4HIT_" + superdetector;
0095     }
0096     else
0097     {
0098       nodename = "G4HIT_" + detector_type + "_" + std::to_string(layer);
0099     }
0100     PHG4HitContainer* cylinder_hits = findNode::getClass<PHG4HitContainer>(topNode, nodename);
0101     if (!cylinder_hits)
0102     {
0103       dstNode->addNode(new PHIODataNode<PHObject>(cylinder_hits = new PHG4HitContainer(nodename), nodename, "PHObject"));
0104     }
0105     cylinder_hits->AddLayer(layer);
0106     if (absorberactive)
0107     {
0108       if (superdetector != "NONE")
0109       {
0110         nodename = "G4HIT_ABSORBER_" + superdetector;
0111       }
0112       else
0113       {
0114         nodename = "G4HIT_ABSORBER_" + detector_type + "_" + std::to_string(layer);
0115       }
0116       PHG4HitContainer* cylinder_hits_2 = findNode::getClass<PHG4HitContainer>(topNode, nodename);
0117       if (!cylinder_hits_2)
0118       {
0119         dstNode->addNode(new PHIODataNode<PHObject>(cylinder_hits_2 = new PHG4HitContainer(nodename), nodename, "PHObject"));
0120       }
0121       cylinder_hits_2->AddLayer(layer);
0122     }
0123     steppingAction_ = new PHG4HcalSteppingAction(detector_);
0124     steppingAction_->set_zmin(zpos - detlength / 2.);
0125     steppingAction_->set_zmax(zpos + detlength / 2.);
0126     if (light_balance_)
0127     {
0128       steppingAction_->SetLightCorrection(light_balance_inner_radius_,
0129                                           light_balance_inner_corr_,
0130                                           light_balance_outer_radius_,
0131                                           light_balance_outer_corr_);
0132       steppingAction_->SetLightScintModel(light_scint_model_);
0133     }
0134   }
0135
0136   return 0;
0137 }
0138
0139 //_______________________________________________________________________
0140 int PHG4HcalSubsystem::process_event(PHCompositeNode* topNode)
0141 {
0142   // pass top node to stepping action so that it gets
0143   // relevant nodes needed internally
0144   if (steppingAction_)
0145   {
0146     steppingAction_->SetInterfacePointers(topNode);
0147   }
0148   return 0;
0149 }
0150
0151 //_______________________________________________________________________
0152 PHG4Detector* PHG4HcalSubsystem::GetDetector() const
0153 {
0154   return detector_;
0155 }
0156
0157 //_______________________________________________________________________
0158 PHG4SteppingAction* PHG4HcalSubsystem::GetSteppingAction() const
0159 {
0160   return steppingAction_;
0161 }
0162
0163 void PHG4HcalSubsystem::Print(const std::string& what) const
0164 {
0165   detector_->Print(what);
0166   return;
0167 }
0168
0169 void PHG4HcalSubsystem::SetTiltViaNcross(const int ncross)
0170 {
0171   if (ncross == 0)
0172   {
0173     std::cout << Name() << " invalid crossing number " << ncross
0174               << " how do you think we can construct a meaningful detector with this number????" << std::endl;
0175     exit(1);
0176   }
0177   // The delta phi angle between 2 adjacent slats is just 360/nslats. If
0178   // there is just one crossing the tips of each slat can extend from
0179   // phi-delta-phi/2 to phi+delta-phi/2. G4 rotates around the center of a
0180   // slat so we are dealing in increments of just delta-phi/2. This
0181   // has to be multiplied by the number of crossings we want.
0182   //
0183   // To find this tilt angle we have to calculate a triangle
0184   // the long sides are the outer radius and the
0185   // inner radius + 1/2 tracker thickness
0186   // the angle between these sides is just (360/nslat)/2*(cross)
0187   // the we use a/sin(alpha) = b/sin(beta)
0188   // and c^2 = a^2+b^2 -2ab*cos(gamma)
0189   // the solution is not unique, we have to pick 180-beta
0190   // but the slat angle is 180-beta (it's on the other side of the triangle
0191   // just draw the damned thing if this is confusing)
0192   double sign = 1;
0193   if (ncross < 0)
0194   {
0195     sign = -1;
0196   }
0197   int ncr = std::abs(ncross);
0198   double thick = TrackerThickness;
0199   double c = radius + thick / 2.;
0200   double b = radius + thick;
0201
0202   double alpha = 0;
0203   alpha = ((360. / _sciNum * M_PI / 180.) / 2.) * ncr;
0204   double sinb = sin(alpha) * b / (sqrt(b * b + c * c - 2 * b * c * cos(alpha)));
0205   double beta = asin(sinb) * 180. / M_PI;  // this is already the slat angle
0206   _sciTilt = beta * sign;
0207   // print out triangle
0208   //   double tbeta = 180 - beta;
0209   //   double gamma = 180 - (alpha * 180. / M_PI) - tbeta;
0210   //   double a = b * sin(alpha) / sinb;
0211   //   std::cout << "triangle length: a: " << a
0212   //        << ", b: " << b << ", c: " << c << std::endl;
0213   //   std::cout << "triangle angle: alpha: " << alpha*180./M_PI
0214   //        << ", beta: " << tbeta
0215   //        << ", gamma: " << gamma
0216   //        << std::endl;
0217   std::cout << Name() << ": SetTiltViaNcross(" << ncross << ") setting slat angle to: " << _sciTilt << " degrees" << std::endl;
0218   return;
0219 }