g4simulation/g4detectors/PHG4OuterHcalSubsystem.cc

0001 #include "PHG4OuterHcalSubsystem.h"
0002
0003 #include "PHG4HcalDefs.h"
0004 #include "PHG4OuterHcalDetector.h"
0005 #include "PHG4OuterHcalDisplayAction.h"
0006 #include "PHG4OuterHcalSteppingAction.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
0014 #include <phool/PHCompositeNode.h>
0015 #include <phool/PHIODataNode.h>    // for PHIODataNode
0016 #include <phool/PHNode.h>          // for PHNode
0017 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0018 #include <phool/PHObject.h>        // for PHObject
0019 #include <phool/getClass.h>
0020
0021 #include <iostream>  // for operator<<, basic_ostream
0022 #include <limits>
0023 #include <set>       // for set
0024 #include <sstream>
0025
0026 class PHG4Detector;
0027
0028 //_______________________________________________________________________
0029 PHG4OuterHcalSubsystem::PHG4OuterHcalSubsystem(const std::string &name, const int lyr)
0030   : PHG4DetectorSubsystem(name, lyr)
0031 {
0032   InitializeParameters();
0033 }
0034
0035 //_______________________________________________________________________
0036 PHG4OuterHcalSubsystem::~PHG4OuterHcalSubsystem()
0037 {
0038   delete m_DisplayAction;
0039 }
0040
0041 //_______________________________________________________________________
0042 int PHG4OuterHcalSubsystem::InitRunSubsystem(PHCompositeNode *topNode)
0043 {
0044   PHNodeIterator iter(topNode);
0045   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0046
0047   // create display settings before detector
0048   m_DisplayAction = new PHG4OuterHcalDisplayAction(Name());
0049
0050   // create detector
0051   m_Detector = new PHG4OuterHcalDetector(this, topNode, GetParams(), Name());
0052   m_Detector->SuperDetector(SuperDetector());
0053   m_Detector->OverlapCheck(CheckOverlap());
0054
0055   if (GetParams()->get_int_param("active"))
0056   {
0057     std::set<std::string> nodes;
0058     PHNodeIterator dstIter(dstNode);
0059     PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstIter.findFirst("PHCompositeNode", SuperDetector()));
0060     if (!DetNode)
0061     {
0062       DetNode = new PHCompositeNode(SuperDetector());
0063       dstNode->addNode(DetNode);
0064     }
0065     std::string nodename;
0066     if (SuperDetector() != "NONE")
0067     {
0068       nodename = "G4HIT_" + SuperDetector();
0069     }
0070     else
0071     {
0072       nodename = "G4HIT_" + Name();
0073     }
0074     nodes.insert(nodename);
0075     if (GetParams()->get_int_param("absorberactive"))
0076     {
0077       if (SuperDetector() != "NONE")
0078       {
0079         nodename = "G4HIT_ABSORBER_" + SuperDetector();
0080       }
0081       else
0082       {
0083         nodename = "G4HIT_ABSORBER_" + Name();
0084       }
0085       nodes.insert(nodename);
0086     }
0087     for (const auto &node : nodes)
0088     {
0089       PHG4HitContainer *g4_hits = findNode::getClass<PHG4HitContainer>(topNode, node);
0090       if (!g4_hits)
0091       {
0092         g4_hits = new PHG4HitContainer(node);
0093         DetNode->addNode(new PHIODataNode<PHObject>(g4_hits, node, "PHObject"));
0094       }
0095     }
0096     // create stepping action
0097     m_SteppingAction = new PHG4OuterHcalSteppingAction(m_Detector, GetParams());
0098     m_SteppingAction->InitWithNode(topNode);
0099   }
0100   else
0101   {
0102     if (GetParams()->get_int_param("blackhole"))
0103     {
0104       m_SteppingAction = new PHG4OuterHcalSteppingAction(m_Detector, GetParams());
0105       m_SteppingAction->InitWithNode(topNode);
0106     }
0107   }
0108
0109   return 0;
0110 }
0111
0112 //_______________________________________________________________________
0113 int PHG4OuterHcalSubsystem::process_event(PHCompositeNode *topNode)
0114 {
0115   // pass top node to stepping action so that it gets
0116   // relevant nodes needed internally
0117   if (m_SteppingAction)
0118   {
0119     m_SteppingAction->SetInterfacePointers(topNode);
0120   }
0121   return 0;
0122 }
0123
0124 void PHG4OuterHcalSubsystem::Print(const std::string &what) const
0125 {
0126   std::cout << "Outer Hcal Parameters: " << std::endl;
0127   GetParams()->Print();
0128   if (m_Detector)
0129   {
0130     m_Detector->Print(what);
0131   }
0132   return;
0133 }
0134
0135 //_______________________________________________________________________
0136 PHG4Detector *PHG4OuterHcalSubsystem::GetDetector() const
0137 {
0138   return m_Detector;
0139 }
0140
0141 void PHG4OuterHcalSubsystem::SetLightCorrection(const double inner_radius, const double inner_corr, const double outer_radius, const double outer_corr)
0142 {
0143   set_double_param("light_balance_inner_corr", inner_corr);
0144   set_double_param("light_balance_inner_radius", inner_radius);
0145   set_double_param("light_balance_outer_corr", outer_corr);
0146   set_double_param("light_balance_outer_radius", outer_radius);
0147   return;
0148 }
0149
0150 void PHG4OuterHcalSubsystem::SetDefaultParameters()
0151 {
0152   set_default_double_param("inner_radius", 183.3);
0153   set_default_double_param("light_balance_inner_corr", std::numeric_limits<double>::quiet_NaN());
0154   set_default_double_param("light_balance_inner_radius", std::numeric_limits<double>::quiet_NaN());
0155   set_default_double_param("light_balance_outer_corr", std::numeric_limits<double>::quiet_NaN());
0156   set_default_double_param("light_balance_outer_radius", std::numeric_limits<double>::quiet_NaN());
0157   set_default_double_param("phistart", std::numeric_limits<double>::quiet_NaN());
0158   set_default_double_param("scinti_eta_coverage_neg", 1.1);
0159   set_default_double_param("scinti_eta_coverage_pos", 1.1);
0160   // some math issue in the code does not subtract the magnet cutout correctly
0161   // (maybe some factor of 2 in a G4 volume creation)
0162   // The engineering drawing values are:
0163   //  set_default_double_param("magnet_cutout_radius", 195.31);
0164   //  set_default_double_param("magnet_cutout_scinti_radius", 195.96);
0165   // seting this to these values results in the correct edges
0166   // (verified by looking at the G4 hit coordinates of the inner edges)
0167   set_default_double_param("magnet_cutout_radius", 195.72);
0168   set_default_double_param("magnet_cutout_scinti_radius", 197.04);
0169   set_default_double_param("outer_radius", 264.71);
0170   set_default_double_param("place_x", 0.);
0171   set_default_double_param("place_y", 0.);
0172   set_default_double_param("place_z", 0.);
0173   set_default_double_param("rot_x", 0.);
0174   set_default_double_param("rot_y", 0.);
0175   set_default_double_param("rot_z", 0.);
0176   set_default_double_param("tmin", -20.);
0177   set_default_double_param("tmax", 60.);
0178   set_default_double_param("dt", 100.);
0179   set_default_double_param("scinti_eta_coverage", 1.1);
0180   set_default_double_param("scinti_gap", 0.85);
0181   set_default_double_param("scinti_gap_neighbor", 0.1);
0182   set_default_double_param("scinti_inner_radius", 183.89);
0183   // some math issue in the code subtracts 0.1mm+ so the scintillator
0184   // does not end at 263.27 as per drawing but at 263.26
0185   // adding 0.125mm compensates for this (so 263.2825 gives the desired 263.27
0186   set_default_double_param("scinti_outer_radius", 263.2825);
0187   set_default_double_param("scinti_tile_thickness", 0.7);
0188   set_default_double_param("size_z", 304.91 * 2);
0189   set_default_double_param("steplimits", std::numeric_limits<double>::quiet_NaN());
0190   set_default_double_param("tilt_angle", -11.23);  // engineering drawing
0191                                                    // corresponds very closely to 4 crossinge (-11.7826 deg)
0192
0193   set_default_int_param("field_check", 0);
0194   set_default_int_param("light_scint_model", 1);
0195   set_default_int_param("magnet_cutout_first_scinti", 8);  // tile start at 0, drawing tile starts at 1
0196   set_default_int_param("etabins", 24);
0197   set_default_int_param("saveg4hit", 1);
0198
0199   // if ncross is set (and tilt_angle is NAN) tilt_angle is calculated
0200   // from number of crossings
0201   set_default_int_param("ncross", 0);
0202   set_default_int_param("n_towers", 64);
0203   set_default_int_param(PHG4HcalDefs::scipertwr, 5);
0204   set_default_int_param("n_scinti_tiles", 12);
0205
0206   set_default_string_param("material", "Steel_1006");
0207   std::string defaultmapfilename;
0208   const char *Calibroot = getenv("CALIBRATIONROOT");
0209   if (Calibroot)
0210   {
0211     defaultmapfilename = Calibroot;
0212     defaultmapfilename += "/HCALOUT/tilemap/oHCALMaps092021.root";
0213   }
0214   set_default_string_param("MapFileName", defaultmapfilename);
0215   set_default_string_param("MapHistoName", "hCombinedMap");
0216 }