g4simulation/g4tpc/PHG4TpcSubsystem.cc

0001 #include "PHG4TpcSubsystem.h"
0002 #include "PHG4TpcDetector.h"
0003 #include "PHG4TpcDisplayAction.h"
0004 #include "PHG4TpcSteppingAction.h"
0005
0006 #include <g4detectors/PHG4DetectorSubsystem.h>  // for PHG4DetectorSubsystem
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 #include <phool/RunnumberRange.h>
0021 #include <phool/recoConsts.h>
0022
0023 #include <iostream>  // for operator<<, basic_ost...
0024 #include <set>
0025
0026 class PHG4Detector;
0027
0028 //_______________________________________________________________________
0029 PHG4TpcSubsystem::PHG4TpcSubsystem(const std::string &name, const int lyr)
0030   : PHG4DetectorSubsystem(name, lyr)
0031 {
0032   InitializeParameters();
0033 }
0034
0035 //_______________________________________________________________________
0036 PHG4TpcSubsystem::~PHG4TpcSubsystem()
0037 {
0038   delete m_DisplayAction;
0039 }
0040
0041 //_______________________________________________________________________
0042 int PHG4TpcSubsystem::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 (detector adds its volumes to it)
0048   m_DisplayAction = new PHG4TpcDisplayAction(Name());
0049   // create detector
0050   m_Detector = new PHG4TpcDetector(this, topNode, GetParams(), Name());
0051   m_Detector->SuperDetector(SuperDetector());
0052   m_Detector->OverlapCheck(CheckOverlap());
0053   std::set<std::string> nodes;
0054   if (GetParams()->get_int_param("active"))
0055   {
0056     PHNodeIterator dstIter(dstNode);
0057     PHCompositeNode *DetNode = dstNode;
0058     if (SuperDetector() != "NONE" && !SuperDetector().empty())
0059     {
0060       PHNodeIterator iter_dst(dstNode);
0061       DetNode = dynamic_cast<PHCompositeNode *>(iter_dst.findFirst("PHCompositeNode", SuperDetector()));
0062       if (!DetNode)
0063       {
0064         DetNode = new PHCompositeNode(SuperDetector());
0065         dstNode->addNode(DetNode);
0066       }
0067     }
0068     std::string detector_suffix = SuperDetector();
0069     if (detector_suffix == "NONE" || detector_suffix.empty())
0070     {
0071       detector_suffix = Name();
0072     }
0073     m_HitNodeName = "G4HIT_" + detector_suffix;
0074     nodes.insert(m_HitNodeName);
0075     m_AbsorberNodeName = "G4HIT_ABSORBER_" + detector_suffix;
0076     if (GetParams()->get_int_param("absorberactive"))
0077     {
0078       nodes.insert(m_AbsorberNodeName);
0079     }
0080     for (const auto &nodename : nodes)
0081     {
0082       PHG4HitContainer *g4_hits = findNode::getClass<PHG4HitContainer>(topNode, nodename);
0083       if (!g4_hits)
0084       {
0085         g4_hits = new PHG4HitContainer(nodename);
0086         DetNode->addNode(new PHIODataNode<PHObject>(g4_hits, nodename, "PHObject"));
0087       }
0088     }
0089
0090     // create stepping action
0091     m_SteppingAction = new PHG4TpcSteppingAction(m_Detector, GetParams());
0092     m_SteppingAction->SetHitNodeName("G4HIT", m_HitNodeName);
0093     m_SteppingAction->SetHitNodeName("G4HIT_ABSORBER", m_AbsorberNodeName);
0094   }
0095   else
0096   {
0097     // if this is a black hole it does not have to be active
0098     if (GetParams()->get_int_param("blackhole"))
0099     {
0100       m_SteppingAction = new PHG4TpcSteppingAction(m_Detector, GetParams());
0101     }
0102   }
0103   return 0;
0104 }
0105
0106 //_______________________________________________________________________
0107 int PHG4TpcSubsystem::process_event(PHCompositeNode *topNode)
0108 {
0109   // pass top node to stepping action so that it gets
0110   // relevant nodes needed internally
0111   if (m_SteppingAction)
0112   {
0113     m_SteppingAction->SetInterfacePointers(topNode);
0114   }
0115   return 0;
0116 }
0117
0118 void PHG4TpcSubsystem::Print(const std::string &what) const
0119 {
0120   std::cout << Name() << " Parameters: " << std::endl;
0121   GetParams()->Print();
0122   if (m_Detector)
0123   {
0124     m_Detector->Print(what);
0125   }
0126   if (m_SteppingAction)
0127   {
0128     m_SteppingAction->Print(what);
0129   }
0130
0131   return;
0132 }
0133
0134 //_______________________________________________________________________
0135 PHG4Detector *PHG4TpcSubsystem::GetDetector() const
0136 {
0137   return m_Detector;
0138 }
0139
0140 void PHG4TpcSubsystem::SetDefaultParameters()
0141 {
0142   set_default_double_param("gas_inner_radius", 21.6);
0143   set_default_double_param("gas_outer_radius", 76.4);
0144   set_default_double_param("place_x", 0.);
0145   set_default_double_param("place_y", 0.);
0146   set_default_double_param("place_z", 0.);
0147   set_default_double_param("rot_x", 0.);
0148   set_default_double_param("rot_y", 0.);
0149   set_default_double_param("rot_z", 0.);
0150   set_default_double_param("tpc_length", 205.21);  // 2 * (maxdrift 102.325 + CM halfwidth 0.28) cm
0151
0152   set_default_double_param("steplimits", 1);  // 1cm by default
0153
0154   // material budget:
0155   // Cu (all layers): 0.5 oz cu per square foot, 1oz == 0.0347mm --> 0.5 oz ==  0.00347cm/2.
0156   // Kapton insulation 18 layers of * 5mil = 18*0.0127=0.2286
0157   // 250 um FR4 (Substrate for Cu layers)
0158   // HoneyComb (nomex) 1/2 inch=0.5*2.54 cm
0159   set_default_string_param("cage_layer_1_material", "G4_Cu");
0160   set_default_double_param("cage_layer_1_thickness", 0.00347 / 2.);
0161
0162   set_default_string_param("cage_layer_2_material", "FR4");
0163   set_default_double_param("cage_layer_2_thickness", 0.025);
0164
0165   set_default_string_param("cage_layer_3_material", "NOMEX");
0166   set_default_double_param("cage_layer_3_thickness", 0.5 * 2.54);
0167
0168   set_default_string_param("cage_layer_4_material", "G4_Cu");
0169   set_default_double_param("cage_layer_4_thickness", 0.00347 / 2.);
0170
0171   set_default_string_param("cage_layer_5_material", "FR4");
0172   set_default_double_param("cage_layer_5_thickness", 0.025);
0173
0174   set_default_string_param("cage_layer_6_material", "G4_KAPTON");
0175   set_default_double_param("cage_layer_6_thickness", 0.2286);
0176
0177   set_default_string_param("cage_layer_7_material", "G4_Cu");
0178   set_default_double_param("cage_layer_7_thickness", 0.00347 / 2.);
0179
0180   set_default_string_param("cage_layer_8_material", "G4_KAPTON");
0181   set_default_double_param("cage_layer_8_thickness", 0.05);  // 50 um
0182
0183   set_default_string_param("cage_layer_9_material", "G4_Cu");
0184   set_default_double_param("cage_layer_9_thickness", 0.00347 / 2.);
0185
0186   // Thomas K Hemmick <Thomas.Hemmick@stonybrook.edu>
0187   //  The total thickness along Zed would be 5.6 millimeters (+/- 2.8 mm around Zed=0).
0188   //  The outer surfaces would have 0.005 inches (125 um) FR4 coated with a negligible thickness of Al. (revised to Au as below)
0189   //  The interior would be some stiffener of either honeycomb or rohacell.  The range of radiation lengths for this material are:
0190   //  Large cell honeycomb:  1450 cm  (0.028 g/cm^3 density)
0191   //  rohacell:  760 cm (0.052 g/cm^3 density)
0192   //  Close cell honeycomb:  635 cm (0.064 g/cm^3 density)
0193   //  I think a calculation just for the rohacell would be more than sufficient.
0194   set_default_string_param("window_core_material", "ROHACELL_FOAM_51");
0195   set_default_double_param("window_thickness", 0.56);  // overall thickness
0196   // I just checked with PC manufacturers and we can get 8.9 micron thick copper in reasonably large sheets.
0197   //  At normal incidence, 8.9 microns is 0.06% of a radiation length.
0198   set_default_string_param("window_surface1_material", "G4_Cu");
0199   set_default_double_param("window_surface1_thickness", 8.9e-4);  // 8.9  um outter shell thickness be default
0200   // The FR4 should be either 5 or 10 mils thick.  10 mils is 254 microns and 5 mils is 0.127 microns.  I think either of these is mechanically fine...
0201   set_default_string_param("window_surface2_material", "FR4");
0202   set_default_double_param("window_surface2_thickness", 0.0127);  // 127  um 2nd shell thickness be default
0203
0204   // for geonode initialization
0205   set_default_double_param("drift_velocity_sim", 0.008);
0206
0207   set_default_int_param("ntpc_layers_inner", 16);
0208   set_default_int_param("ntpc_layers_mid", 16);
0209   set_default_int_param("ntpc_layers_outer", 16);
0210   set_default_int_param("tpc_minlayer_inner", 7);
0211
0212   set_default_double_param("tpc_minradius_inner", 31.105);  // 30.0);  // cm
0213   set_default_double_param("tpc_minradius_mid", 41.153);    // 40.0);
0214   set_default_double_param("tpc_minradius_outer", 58.367);  // 60.0);
0215
0216   set_default_double_param("tpc_maxradius_inner", 40.249);  // 40.0);  // cm
0217   set_default_double_param("tpc_maxradius_mid", 57.475);    // 60.0);
0218   set_default_double_param("tpc_maxradius_outer", 75.911);  // 77.0);  // from Tom
0219
0220   set_default_double_param("maxdriftlength", 102.325);       // cm
0221   set_default_double_param("CM_halfwidth", 0.28);       // cm
0222   recoConsts *rc = recoConsts::instance();
0223   int runnumber = rc->get_IntFlag("RUNNUMBER");
0224   if (runnumber < RunnumberRange::RUN2PP_FIRST)
0225   {
0226     // current default clock used in simulation. Need to decide how to handle long term
0227     // and ensure that the electron drift uses the same value
0228     set_default_double_param("tpc_adc_clock", 53.326184);  // ns, for 18.83 MHz clock
0229   }
0230   else if ( runnumber < RunnumberRange::RUN3_TPCFW_CLOCK_CHANGE)
0231   {
0232     set_default_double_param("tpc_adc_clock", 50.037280);  // ns, for 20 MHz clock
0233   }
0234   else
0235   {
0236     set_default_double_param("tpc_adc_clock", 56.881262);  // ns, for 17.5 MHz clock
0237   }
0238   if((runnumber <= RunnumberRange::RUN2AUAU_FIRST && runnumber >= RunnumberRange::RUN2PP_FIRST)
0239     || (runnumber > RunnumberRange::RUN3PP_FIRST))
0240   {
0241     // with drift length of 102cm and clock of 50 ns we get 542 time bins
0242     // to get to 1024 samples (time bins0) we therefore need 1024-542=482 additional time bins
0243     // and 482*50ns = 24100 ns. Add one some extra samples as they can always
0244     // be cut out later in the reconstruction
0245     set_default_double_param("extended_readout_time", 24800);  // ns, to account for 1024 samples
0246   }
0247   else
0248   {
0249     // with drift length of 102cm and clock of 56 ns we get 477 time bins
0250     // this already exceeds the 450 samples in the data
0251     set_default_double_param("extended_readout_time", 0);  // ns, to account for 450 samples
0252
0253   }
0254
0255   set_default_double_param("tpc_sector_phi_inner", 0.5024);  // 2 * M_PI / 12 );//sector size in phi for R1 sector
0256   set_default_double_param("tpc_sector_phi_mid", 0.5087);    // 2 * M_PI / 12 );//sector size in phi for R2 sector
0257   set_default_double_param("tpc_sector_phi_outer", 0.5097);  // 2 * M_PI / 12 );//sector size in phi for R3 sector
0258
0259   set_default_int_param("ntpc_phibins_inner", 1128);  // 94 * 12
0260   set_default_int_param("ntpc_phibins_mid", 1536);    // 128 * 12
0261   set_default_int_param("ntpc_phibins_outer", 2304);  // 192 * 12
0262
0263   set_default_double_param("TPC_gas_temperature", 15.0);  // in celcius
0264   set_default_double_param("TPC_gas_pressure", 1.0);      // in atmospheres
0265   set_default_double_param("Ne_frac", 0.00);
0266   set_default_double_param("Ar_frac", 0.75);
0267   set_default_double_param("CF4_frac", 0.20);
0268   set_default_double_param("N2_frac", 0.00);
0269   set_default_double_param("isobutane_frac", 0.05);
0270 }