g4simulation/g4detectors/PHG4CylinderSubsystem.cc

0001 #include "PHG4CylinderSubsystem.h"
0002 #include "PHG4CylinderDetector.h"
0003 #include "PHG4CylinderDisplayAction.h"
0004 #include "PHG4CylinderGeomContainer.h"
0005 #include "PHG4CylinderGeomv1.h"
0006 #include "PHG4CylinderSteppingAction.h"
0007
0008 #include <phparameter/PHParameters.h>
0009
0010 #include <g4main/PHG4DisplayAction.h>  // for PHG4DisplayAction
0011 #include <g4main/PHG4HitContainer.h>
0012 #include <g4main/PHG4SteppingAction.h>  // for PHG4SteppingAction
0013 #include <g4main/PHG4Utils.h>
0014
0015 #include <phool/PHCompositeNode.h>
0016 #include <phool/PHIODataNode.h>    // for PHIODataNode
0017 #include <phool/PHNode.h>          // for PHNode
0018 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0019 #include <phool/PHObject.h>        // for PHObject
0020 #include <phool/getClass.h>
0021 #include <phool/recoConsts.h>
0022
0023 #include <cmath>     // for NAN
0024 #include <iostream>  // for operator<<, basic_ostream, endl
0025 #include <sstream>
0026
0027 class PHG4CylinderGeom;
0028 class PHG4Detector;
0029
0030 //_______________________________________________________________________
0031 PHG4CylinderSubsystem::PHG4CylinderSubsystem(const std::string &na, const int lyr)
0032   : PHG4DetectorSubsystem(na, lyr)
0033 {
0034   m_ColorArray.fill(std::numeric_limits<double>::quiet_NaN());
0035   InitializeParameters();
0036 }
0037
0038 //_______________________________________________________________________
0039 PHG4CylinderSubsystem::~PHG4CylinderSubsystem()
0040 {
0041   delete m_DisplayAction;
0042 }
0043
0044 //_______________________________________________________________________
0045 int PHG4CylinderSubsystem::InitRunSubsystem(PHCompositeNode *topNode)
0046 {
0047   // create hit list only for active layers
0048   double detlength = GetParams()->get_double_param("length");
0049   if (!std::isfinite(detlength) && GetParams()->get_int_param("lengthviarapidity"))
0050   {
0051     GetParams()->set_double_param("length", PHG4Utils::GetLengthForRapidityCoverage(GetParams()->get_double_param("radius") + GetParams()->get_double_param("thickness")) * 2);
0052     detlength = GetParams()->get_double_param("length");
0053   }
0054   else
0055   {
0056     GetParams()->set_int_param("lengthviarapidity", 0);
0057   }
0058   // use world material if material was not set so far
0059   if (GetParams()->get_string_param("material") == "WorldMaterial")
0060   {
0061     recoConsts *rc = recoConsts::instance();
0062     GetParams()->set_string_param("material", rc->get_StringFlag("WorldMaterial"));
0063   }
0064   // create display settings before detector
0065   PHG4CylinderDisplayAction *disp_action = new PHG4CylinderDisplayAction(Name(), GetParams());
0066   if (std::isfinite(m_ColorArray[0]) &&
0067       std::isfinite(m_ColorArray[1]) &&
0068       std::isfinite(m_ColorArray[2]) &&
0069       std::isfinite(m_ColorArray[3]))
0070   {
0071     disp_action->SetColor(m_ColorArray[0], m_ColorArray[1], m_ColorArray[2], m_ColorArray[3]);
0072   }
0073   m_DisplayAction = disp_action;
0074
0075   // create detector
0076   m_Detector = new PHG4CylinderDetector(this, topNode, GetParams(), Name(), GetLayer());
0077   m_Detector->SuperDetector(SuperDetector());
0078   m_Detector->OverlapCheck(CheckOverlap());
0079   if (GetParams()->get_int_param("active"))
0080   {
0081     PHNodeIterator iter(topNode);
0082     PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0083     PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0084
0085     std::string nodename;
0086     std::string geonode;
0087     if (SuperDetector() != "NONE")
0088     {
0089       // create super detector subnodes
0090       PHNodeIterator iter_dst(dstNode);
0091       PHCompositeNode *superSubNode = dynamic_cast<PHCompositeNode *>(iter_dst.findFirst("PHCompositeNode", SuperDetector()));
0092       if (!superSubNode)
0093       {
0094         superSubNode = new PHCompositeNode(SuperDetector());
0095         dstNode->addNode(superSubNode);
0096       }
0097       dstNode = superSubNode;
0098       PHNodeIterator iter_run(runNode);
0099       superSubNode = dynamic_cast<PHCompositeNode *>(iter_run.findFirst("PHCompositeNode", SuperDetector()));
0100       if (!superSubNode)
0101       {
0102         superSubNode = new PHCompositeNode(SuperDetector());
0103         runNode->addNode(superSubNode);
0104       }
0105       runNode = superSubNode;
0106
0107       nodename = "G4HIT_" + SuperDetector();
0108       geonode = "CYLINDERGEOM_" + SuperDetector();
0109     }
0110
0111     else
0112     {
0113       nodename = "G4HIT_" + Name();
0114       geonode = "CYLINDERGEOM_" + Name();
0115     }
0116     PHG4HitContainer *cylinder_hits = findNode::getClass<PHG4HitContainer>(topNode, nodename);
0117     if (!cylinder_hits)
0118     {
0119       dstNode->addNode(new PHIODataNode<PHObject>(cylinder_hits = new PHG4HitContainer(nodename), nodename, "PHObject"));
0120     }
0121     cylinder_hits->AddLayer(GetLayer());
0122     PHG4CylinderGeomContainer *geo = findNode::getClass<PHG4CylinderGeomContainer>(topNode, geonode);
0123     if (!geo)
0124     {
0125       geo = new PHG4CylinderGeomContainer();
0126       PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(geo, geonode, "PHObject");
0127       runNode->addNode(newNode);
0128     }
0129     PHG4CylinderGeom *mygeom = new PHG4CylinderGeomv1(GetParams()->get_double_param("radius"), GetParams()->get_double_param("place_z") - detlength / 2., GetParams()->get_double_param("place_z") + detlength / 2., GetParams()->get_double_param("thickness"));
0130     geo->AddLayerGeom(GetLayer(), mygeom);
0131     auto *tmp = new PHG4CylinderSteppingAction(this, m_Detector, GetParams());
0132     tmp->HitNodeName(nodename);
0133     m_SteppingAction = tmp;
0134   }
0135   else if (GetParams()->get_int_param("blackhole"))
0136   {
0137     m_SteppingAction = new PHG4CylinderSteppingAction(this, m_Detector, GetParams());
0138   }
0139   if (m_SteppingAction)
0140   {
0141     (dynamic_cast<PHG4CylinderSteppingAction *>(m_SteppingAction))->SaveAllHits(m_SaveAllHitsFlag);
0142   }
0143   return 0;
0144 }
0145
0146 //_______________________________________________________________________
0147 int PHG4CylinderSubsystem::process_event(PHCompositeNode *topNode)
0148 {
0149   // pass top node to stepping action so that it gets
0150   // relevant nodes needed internally
0151   if (m_SteppingAction)
0152   {
0153     m_SteppingAction->SetInterfacePointers(topNode);
0154   }
0155   return 0;
0156 }
0157
0158 void PHG4CylinderSubsystem::SetDefaultParameters()
0159 {
0160   set_default_double_param("length", NAN);
0161   set_default_double_param("place_x", 0.);
0162   set_default_double_param("place_y", 0.);
0163   set_default_double_param("place_z", 0.);
0164   set_default_double_param("radius", NAN);
0165   set_default_double_param("steplimits", NAN);
0166   set_default_double_param("thickness", NAN);
0167   set_default_double_param("tmin", NAN);
0168   set_default_double_param("tmax", NAN);
0169   set_default_double_param("rot_x", 0.);
0170   set_default_double_param("rot_y", 0.);
0171   set_default_double_param("rot_z", 0.);
0172   set_default_double_param("start_phi_rad", 0.);
0173   set_default_double_param("delta_phi_rad", M_PI * 2);
0174   set_default_int_param("lengthviarapidity", 1);
0175   set_default_int_param("lightyield", 0);
0176   set_default_int_param("use_g4steps", 0);
0177
0178   // place holder, will be replaced by world material if not set by other means (macro)
0179   set_default_string_param("material", "WorldMaterial");
0180 }
0181
0182 PHG4Detector *
0183 PHG4CylinderSubsystem::GetDetector() const
0184 {
0185   return m_Detector;
0186 }
0187
0188 void PHG4CylinderSubsystem::Print(const std::string &what) const
0189 {
0190   std::cout << Name() << " Parameters: " << std::endl;
0191   if (!BeginRunExecuted())
0192   {
0193     std::cout << "Need to execute BeginRun() before parameter printout is meaningful" << std::endl;
0194     std::cout << "To do so either run one or more events or on the command line execute: " << std::endl;
0195     std::cout << "Fun4AllServer *se = Fun4AllServer::instance();" << std::endl;
0196     std::cout << "PHG4Reco *g4 = (PHG4Reco *) se->getSubsysReco(\"PHG4RECO\");" << std::endl;
0197     std::cout << "g4->InitRun(se->topNode());" << std::endl;
0198     std::cout << "PHG4CylinderSubsystem *cyl = (PHG4CylinderSubsystem *) g4->getSubsystem(\"" << Name() << "\");" << std::endl;
0199     std::cout << "cyl->Print()" << std::endl;
0200     return;
0201   }
0202   GetParams()->Print();
0203   if (m_SteppingAction)
0204   {
0205     m_SteppingAction->Print(what);
0206   }
0207   return;
0208 }