macros/common/G4_TrkrSimulation.C

0001 #ifndef MACRO_G4TRKRSIM_C
0002 #define MACRO_G4TRKRSIM_C
0003
0004 #include <GlobalVariables.C>
0005
0006 #include <G4_ActsGeom.C>
0007 #include <G4_TrkrVariables.C>
0008
0009 #include <g4detectors/PHG4CylinderSubsystem.h>
0010
0011 #include <g4mvtx/PHG4MvtxDefs.h>
0012 #include <g4mvtx/PHG4MvtxDigitizer.h>
0013 #include <g4mvtx/PHG4MvtxHitReco.h>
0014 #include <g4mvtx/PHG4MvtxSubsystem.h>
0015
0016 #include <g4intt/PHG4InttDefs.h>
0017 #include <g4intt/PHG4InttDigitizer.h>
0018 #include <g4intt/PHG4InttHitReco.h>
0019 #include <g4intt/PHG4InttSubsystem.h>
0020
0021 #include <g4tpc/PHG4TpcCentralMembrane.h>
0022 #include <g4tpc/PHG4TpcDigitizer.h>
0023 #include <g4tpc/PHG4TpcDirectLaser.h>
0024 #include <g4tpc/PHG4TpcDistortion.h>
0025 #include <g4tpc/PHG4TpcElectronDrift.h>
0026 #include <g4tpc/PHG4TpcEndCapSubsystem.h>
0027 #include <g4tpc/PHG4TpcPadPlane.h>
0028 #include <g4tpc/PHG4TpcPadPlaneReadout.h>
0029 #include <g4tpc/PHG4TpcSubsystem.h>
0030
0031 #pragma GCC diagnostic push
0032 #pragma GCC diagnostic ignored "-Wundefined-internal"
0033 #include <tpc/TpcClusterZCrossingCorrection.h>
0034 #pragma GCC diagnostic pop
0035
0036 #include <g4micromegas/PHG4MicromegasDigitizer.h>
0037 #include <g4micromegas/PHG4MicromegasHitReco.h>
0038 #include <g4micromegas/PHG4MicromegasSubsystem.h>
0039
0040 #include <g4main/PHG4Reco.h>
0041
0042 #include <ffamodules/CDBInterface.h>
0043
0044 #include <cdbobjects/CDBTTree.h>
0045
0046 #include <TFile.h>
0047
0048 #include <cmath>
0049 #include <vector>
0050
0051 R__LOAD_LIBRARY(libg4tpc.so)
0052 R__LOAD_LIBRARY(libg4mvtx.so)
0053 R__LOAD_LIBRARY(libg4intt.so)
0054 R__LOAD_LIBRARY(libg4micromegas.so)
0055
0056 void MvtxInit()
0057 {
0058   BlackHoleGeometry::max_radius = std::max(BlackHoleGeometry::max_radius, 12.);
0059   BlackHoleGeometry::min_z = std::min(BlackHoleGeometry::min_z, -165.);
0060   BlackHoleGeometry::max_z = std::max(BlackHoleGeometry::max_z, 24.);
0061 }
0062
0063 double Mvtx(PHG4Reco* g4Reco, double radius,
0064             const int /*supportactive*/ = 0)
0065 {
0066   bool maps_overlapcheck = Enable::OVERLAPCHECK || Enable::MVTX_OVERLAPCHECK;
0067   int verbosity = std::max(Enable::VERBOSITY, Enable::MVTX_VERBOSITY);
0068   bool SupportActive = Enable::SUPPORT || Enable::MVTX_SUPPORT;
0069
0070   PHG4MvtxSubsystem* mvtx = new PHG4MvtxSubsystem("MVTX");
0071   mvtx->Verbosity(verbosity);
0072
0073   for (int ilayer = 0; ilayer < G4MVTX::n_maps_layer; ilayer++)
0074   {
0075     double radius_lyr = PHG4MvtxDefs::mvtxdat[ilayer][PHG4MvtxDefs::kRmd];
0076     //    mvtx->set_double_param(ilayer, "layer_z_offset", G4MVTXAlignment::z_offset[ilayer]);
0077     if (verbosity)
0078     {
0079       std::cout << "Create Maps layer " << ilayer << " with radius " << radius_lyr << " mm." << std::endl;
0080     }
0081     radius = radius_lyr / 10.;
0082   }
0083   //  mvtx->set_string_param(PHG4MvtxDefs::GLOBAL, "alignment_path",  G4MVTXAlignment::alignment_path);
0084   mvtx->set_string_param(PHG4MvtxDefs::GLOBAL, "stave_geometry_file", std::string(getenv("CALIBRATIONROOT")) + std::string("/Tracking/geometry/mvtx_stave.gdml"));
0085
0086   mvtx->SetActive();
0087   if (SupportActive)
0088   {
0089     mvtx->SetSupportActive();
0090   }
0091   std::cout << "PHG4MvtxSubsystem: Apply misalignment? Enable::MVTX_APPLYMISALIGNMENT=" << Enable::MVTX_APPLYMISALIGNMENT << std::endl;
0092   mvtx->Apply_Misalignment(Enable::MVTX_APPLYMISALIGNMENT);
0093   if (Enable::MVTX_APPLYMISALIGNMENT)
0094   {
0095     std::string file = CDBInterface::instance()->getUrl("MVTX_ALIGNMENT");
0096     std::cout << "applying with file " << file << std::endl;
0097     mvtx->MisalignmentFile(file);
0098   }
0099   mvtx->OverlapCheck(maps_overlapcheck);
0100   g4Reco->registerSubsystem(mvtx);
0101   radius += G4MVTX::radius_offset;
0102   return radius;
0103 }
0104
0105 void Mvtx_Cells()
0106 {
0107   int verbosity = std::max(Enable::VERBOSITY, Enable::MVTX_VERBOSITY);
0108   Fun4AllServer* se = Fun4AllServer::instance();
0109   // new storage containers
0110   PHG4MvtxHitReco* maps_hits = new PHG4MvtxHitReco("MVTX");
0111   maps_hits->Verbosity(verbosity);
0112
0113   double maps_readout_window = 9900.0;  // ns
0114   double extended_readout_time = 0.0;
0115   if (TRACKING::pp_mode) extended_readout_time = TRACKING::pp_extended_readout_time;
0116   // override the default timing window - default is +/- 5000 ns
0117   maps_hits->set_double_param("mvtx_tmin", -maps_readout_window);
0118   maps_hits->set_double_param("mvtx_tmax", maps_readout_window + extended_readout_time);
0119   std::cout << "PHG4MvtxHitReco: readout window is from " << -maps_readout_window << " to " << maps_readout_window + extended_readout_time << std::endl;
0120   std::cout << "PHG4MvtxHitReco: Mvtx set to streaming in simulation? (Enable::MVTX_STREAMING) " << ((Enable::MVTX_STREAMING) ? "True" : "False") << std::endl;
0121   maps_hits->set_int_param("mvtx_in_sphenix_srdo", Enable::MVTX_STREAMING);
0122   se->registerSubsystem(maps_hits);
0123
0124   PHG4MvtxDigitizer* digimvtx = new PHG4MvtxDigitizer();
0125   digimvtx->Verbosity(verbosity);
0126   // energy deposit in 25 microns = 9.6 KeV = 1000 electrons collected after recombination
0127   // digimvtx->set_adc_scale(0.95e-6);  // default set in code is 0.95e-06, which is 99 electrons
0128   se->registerSubsystem(digimvtx);
0129
0130   return;
0131 }
0132
0133 void InttInit()
0134 {
0135   BlackHoleGeometry::max_radius = std::max(BlackHoleGeometry::max_radius, 20.);  // estimated from display, can be made smaller but good enough
0136   BlackHoleGeometry::max_z = std::max(BlackHoleGeometry::max_z, 410. / 2.);
0137   BlackHoleGeometry::min_z = std::min(BlackHoleGeometry::min_z, -410. / 2.);
0138   // the mvtx is not called if disabled but the default number of layers is set to 3, so we need to set it
0139   // to zero
0140   if (!Enable::MVTX)
0141   {
0142     G4MVTX::n_maps_layer = 0;
0143   }
0144   if (!Enable::INTT_USEG4SURVEYGEOM)
0145   {
0146     G4INTT::sensor_radius[0] = 7.188 - 36e-4;
0147     G4INTT::sensor_radius[1] = 7.732 - 36e-4;
0148     G4INTT::sensor_radius[2] = 9.680 - 36e-4;
0149     G4INTT::sensor_radius[3] = 10.262 - 36e-4;
0150   }
0151 }
0152
0153 double Intt(PHG4Reco* g4Reco, double radius,
0154             const int /*absorberactive*/ = 0)
0155 {
0156   std::cout << "G4_TrkrSimulation::Intt" << std::endl;
0157   int verbosity = std::max(Enable::VERBOSITY, Enable::INTT_VERBOSITY);
0158   bool intt_overlapcheck = Enable::OVERLAPCHECK || Enable::INTT_OVERLAPCHECK;
0159
0160   // instantiate the INTT subsystem and register it
0161   // We make one instance of PHG4INTTSubsystem for all four layers of tracker
0162   // dimensions are in mm, angles are in radians
0163
0164   // PHG4InttSubsystem creates the detetor layer using PHG4InttDetector
0165   // and instantiates the appropriate PHG4SteppingAction
0166
0167   // The length of vpair is used to determine the number of layers
0168   std::vector<std::pair<int, int>> vpair;  // (sphxlayer, inttlayer)
0169   for (int i = 0; i < G4INTT::n_intt_layer; i++)
0170   {
0171     // We want the sPHENIX layer numbers for the Intt to be from n_maps_layer to n_maps_layer+n_intt_layer - 1
0172     vpair.emplace_back(G4MVTX::n_maps_layer + i, i);  // sphxlayer=n_maps_layer+i corresponding to inttlayer=i
0173     if (verbosity) std::cout << "Create strip tracker layer " << vpair[i].second << " as  sphenix layer  " << vpair[i].first << std::endl;
0174   }
0175
0176   PHG4InttSubsystem* sitrack = new PHG4InttSubsystem("INTT", vpair);
0177   sitrack->Verbosity(verbosity);
0178   sitrack->SetActive(1);
0179   sitrack->OverlapCheck(intt_overlapcheck);
0180   std::cout << "PHG4InttSubsystem: Use survey geometry? Enable::INTT_USEG4SURVEYGEOM=" << Enable::INTT_USEG4SURVEYGEOM << std::endl;
0181   sitrack->SetSurveyGeometry(Enable::INTT_USEG4SURVEYGEOM);
0182   if (Enable::INTT_ABSORBER || Enable::ABSORBER)
0183   {
0184     sitrack->SetAbsorberActive();
0185   }
0186   if (Enable::INTT_SUPPORT || Enable::SUPPORT)
0187   {
0188     sitrack->set_int_param(PHG4InttDefs::SUPPORTPARAMS, "supportactive", 1);
0189   }
0190   g4Reco->registerSubsystem(sitrack);
0191
0192   // Set the laddertype and ladder spacing configuration
0193
0194   std::cout << "Intt has " << G4INTT::n_intt_layer << " layers with layer setup:" << std::endl;
0195   for (int i = 0; i < G4INTT::n_intt_layer; i++)
0196   {
0197     std::cout << " Intt layer " << i << " laddertype " << G4INTT::laddertype[i] << " nladders " << G4INTT::nladder[i]
0198               << " sensor radius " << G4INTT::sensor_radius[i] << std::endl;
0199     sitrack->set_int_param(i, "laddertype", G4INTT::laddertype[i]);
0200     sitrack->set_int_param(i, "nladder", G4INTT::nladder[i]);
0201     sitrack->set_double_param(i, "sensor_radius", G4INTT::sensor_radius[i]);  // expecting cm
0202   }
0203
0204   // outer radius marker (translation back to cm)
0205   radius = G4INTT::intt_radius_max * 0.1;
0206   return radius;
0207 }
0208
0209 // Central detector cell reco is disabled as EIC setup use the fast tracking sim for now
0210 void Intt_Cells()
0211 {
0212   int verbosity = std::max(Enable::VERBOSITY, Enable::INTT_VERBOSITY);
0213   Fun4AllServer* se = Fun4AllServer::instance();
0214
0215   // new storage containers
0216   PHG4InttHitReco* reco = new PHG4InttHitReco();
0217   reco->setHotStripMaskFile("INTT_HotMap");
0218
0219   // The timing window defaults are set in the INTT ladder model, they can be overridden here
0220   double extended_readout_time = 0.0;
0221   if (TRACKING::pp_mode) extended_readout_time = TRACKING::pp_extended_readout_time;
0222   reco->set_double_param("tmax", 80.0 + extended_readout_time);
0223   reco->set_double_param("tmin", -20.0);
0224   std::cout << "INTT readout window is set to -20 to " << 80.0 + extended_readout_time << std::endl;
0225   reco->Verbosity(verbosity);
0226   se->registerSubsystem(reco);
0227
0228   // Intt digitization
0229   //===========
0230   // these should be used for the Intt
0231   /*
0232      How threshold are calculated based on default FPHX settings
0233      Four part information goes to the threshold calculation:
0234      1. In 320 um thick silicon, the MIP e-h pair for a nominally indenting tracking is 3.87 MeV/cm * 320 um / 3.62 eV/e-h = 3.4e4 e-h pairs
0235      2. From DOI: 10.1016/j.nima.2014.04.017, FPHX integrator amplifier gain is 100mV / fC. That translate MIP voltage to 550 mV.
0236      3. From [FPHX Final Design Document](https://www.phenix.bnl.gov/WWW/fvtx/DetectorHardware/FPHX/FPHX2_June2009Revision.doc), the DAC0-7 setting for 8-ADC thresholds above the V_ref, as in Table 2 - Register Addresses and Defaults
0237      4, From [FPHX Final Design Document](https://www.phenix.bnl.gov/WWW/fvtx/DetectorHardware/FPHX/FPHX2_June2009Revision.doc) section Front-end Program Bits, the formula to translate DAC setting to comparitor voltages.
0238      The result threshold table based on FPHX default value is as following
0239      | FPHX Register Address | Name            | Default value | Voltage - Vref (mV) | To electrons based on calibration | Electrons | Fraction to MIP |
0240      |-----------------------|-----------------|---------------|---------------------|-----------------------------------|-----------|-----------------|
0241      | 4                     | Threshold DAC 0 | 8             | 32                  | 2500                              | 2000      | 5.85E-02        |
0242      | 5                     | Threshold DAC 1 | 16            | 64                  | 5000                              | 4000      | 1.17E-01        |
0243      | 6                     | Threshold DAC 2 | 32            | 128                 | 10000                             | 8000      | 2.34E-01        |
0244      | 7                     | Threshold DAC 3 | 48            | 192                 | 15000                             | 12000     | 3.51E-01        |
0245      | 8                     | Threshold DAC 4 | 80            | 320                 | 25000                             | 20000     | 5.85E-01        |
0246      | 9                     | Threshold DAC 5 | 112           | 448                 | 35000                             | 28000     | 8.18E-01        |
0247      | 10                    | Threshold DAC 6 | 144           | 576                 | 45000                             | 36000     | 1.05E+00        |
0248      | 11                    | Threshold DAC 7 | 176           | 704                 | 55000                             | 44000     | 1.29E+00        |
0249      DAC0-7 threshold as fraction to MIP voltage are set to PHG4InttDigitizer::set_adc_scale as 3-bit ADC threshold as fractions to MIP energy deposition.
0250      */
0251   std::vector<double> userrange;  // 3-bit ADC threshold relative to the mip_e at each layer.
0252   userrange.push_back(0.0584625322997416);
0253   userrange.push_back(0.116925064599483);
0254   userrange.push_back(0.233850129198966);
0255   userrange.push_back(0.35077519379845);
0256   userrange.push_back(0.584625322997416);
0257   userrange.push_back(0.818475452196383);
0258   userrange.push_back(1.05232558139535);
0259   userrange.push_back(1.28617571059432);
0260
0261   // new containers
0262   PHG4InttDigitizer* digiintt = new PHG4InttDigitizer();
0263   digiintt->Verbosity(verbosity);
0264   digiintt->LoadBadChannelMap("INTT_Hotmap");
0265   for (int i = 0; i < G4INTT::n_intt_layer; i++)
0266   {
0267     digiintt->set_adc_scale(G4MVTX::n_maps_layer + i, userrange);
0268   }
0269   se->registerSubsystem(digiintt);
0270
0271   return;
0272 }
0273
0274 void TPCInit()
0275 {
0276   std::cout << "G4_TrkrSimulation::TpcInit" << std::endl;
0277   BlackHoleGeometry::max_radius = std::max(BlackHoleGeometry::max_radius, 82.);  // radius of the tpc hanger + .3 cm
0278
0279   if (Enable::TPC_ENDCAP)
0280   {
0281     BlackHoleGeometry::max_z = std::max(BlackHoleGeometry::max_z, 130.);
0282     BlackHoleGeometry::min_z = std::min(BlackHoleGeometry::min_z, -130.);
0283   }
0284   else
0285   {
0286     BlackHoleGeometry::max_z = std::max(BlackHoleGeometry::max_z, 211. / 2.);
0287     BlackHoleGeometry::min_z = std::min(BlackHoleGeometry::min_z, -211. / 2.);
0288   }
0289   // the mvtx is not called if disabled but the default number of layers is set to 3,
0290   // so we need to set it to zero
0291   if (!Enable::MVTX)
0292   {
0293     G4MVTX::n_maps_layer = 0;
0294   }
0295   // same for the INTT
0296   if (!Enable::INTT)
0297   {
0298     G4INTT::n_intt_layer = 0;
0299   }
0300
0301   // Set the (static) drift velocity in the cluster Z crossing correction module
0302   TpcClusterZCrossingCorrection::_vdrift = G4TPC::tpc_drift_velocity_reco;
0303 }
0304
0305 //! TPC end cap, wagon wheel, electronics
0306 void TPC_Endcaps(PHG4Reco* g4Reco)
0307 {
0308   bool OverlapCheck = Enable::OVERLAPCHECK || Enable::TPC_OVERLAPCHECK;
0309   bool AbsorberActive = Enable::ABSORBER || Enable::TPC_ABSORBER;
0310
0311   PHG4TpcEndCapSubsystem* tpc_endcap = new PHG4TpcEndCapSubsystem("TPC_ENDCAP");
0312   tpc_endcap->SuperDetector("TPC_ENDCAP");
0313
0314   if (AbsorberActive) tpc_endcap->SetActive();
0315   tpc_endcap->OverlapCheck(OverlapCheck);
0316
0317   //  tpc_endcap->set_int_param("construction_verbosity", 2);
0318
0319   g4Reco->registerSubsystem(tpc_endcap);
0320
0321   return;
0322 }
0323
0324 double TPC(PHG4Reco* g4Reco,
0325            double radius)
0326 {
0327   std::cout << "G4_TrkrSimulation::TPC" << std::endl;
0328   bool OverlapCheck = Enable::OVERLAPCHECK || Enable::TPC_OVERLAPCHECK;
0329   bool AbsorberActive = Enable::ABSORBER || Enable::TPC_ABSORBER;
0330   double drift_vel = G4TPC::tpc_drift_velocity_sim;
0331   if (G4TPC::TPC_GAS_MIXTURE == "NeCF4")
0332   {
0333     drift_vel = G4TPC::NeCF4_drift_velocity;
0334   }
0335   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4")
0336   {
0337     drift_vel = G4TPC::ArCF4_drift_velocity;
0338   }
0339   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4N2")
0340   {
0341     drift_vel = G4TPC::ArCF4N2_drift_velocity;
0342   }
0343   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4Isobutane")
0344   {
0345     drift_vel = G4TPC::ArCF4Isobutane_drift_velocity;
0346   }
0347
0348   PHG4TpcSubsystem* tpc = new PHG4TpcSubsystem("TPC");
0349   tpc->SetActive();
0350   tpc->SuperDetector("TPC");
0351   tpc->set_double_param("steplimits", 1);                  // 1cm steps
0352   tpc->set_double_param("drift_velocity_sim", drift_vel);  // this is the only place that drift_velocity_sim is set
0353   tpc->set_double_param("tpc_length", G4TPC::maxDriftLength * 2 + G4TPC::CM_halfwidth * 2);
0354   tpc->set_double_param("maxdriftlength", G4TPC::maxDriftLength);
0355   tpc->set_double_param("CM_halfwidth", G4TPC::CM_halfwidth);
0356   tpc->set_int_param("tpc_minlayer_inner", G4MVTX::n_maps_layer + G4INTT::n_intt_layer);
0357   tpc->set_int_param("ntpc_layers_inner", G4TPC::n_tpc_layer_inner);
0358   tpc->set_int_param("ntpc_phibins_inner", G4TPC::tpc_layer_rphi_count_inner);
0359
0360   if (AbsorberActive)
0361   {
0362     tpc->SetAbsorberActive();
0363   }
0364
0365   double extended_readout_time = 0.0;
0366   if (TRACKING::pp_mode)
0367   {
0368     extended_readout_time = TRACKING::pp_extended_readout_time;
0369   }
0370
0371   tpc->set_double_param("extended_readout_time", extended_readout_time);
0372
0373   // Note that we default to 75:20:05 Ar:CF4:i-C4H10
0374   if (G4TPC::TPC_GAS_MIXTURE == "NeCF4")
0375   {
0376     tpc->set_double_param("Ne_frac", G4TPC::NeCF4_Ne_frac);
0377     tpc->set_double_param("Ar_frac", G4TPC::NeCF4_Ar_frac);
0378     tpc->set_double_param("CF4_frac", G4TPC::NeCF4_CF4_frac);
0379     tpc->set_double_param("N2_frac", G4TPC::NeCF4_N2_frac);
0380     tpc->set_double_param("isobutane_frac", G4TPC::NeCF4_isobutane_frac);
0381   }
0382   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4")
0383   {
0384     tpc->set_double_param("Ne_frac", G4TPC::ArCF4_Ne_frac);
0385     tpc->set_double_param("Ar_frac", G4TPC::ArCF4_Ar_frac);
0386     tpc->set_double_param("CF4_frac", G4TPC::ArCF4_CF4_frac);
0387     tpc->set_double_param("N2_frac", G4TPC::ArCF4_N2_frac);
0388     tpc->set_double_param("isobutane_frac", G4TPC::ArCF4_isobutane_frac);
0389   }
0390   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4N2")
0391   {
0392     tpc->set_double_param("Ne_frac", G4TPC::ArCF4N2_Ne_frac);
0393     tpc->set_double_param("Ar_frac", G4TPC::ArCF4N2_Ar_frac);
0394     tpc->set_double_param("CF4_frac", G4TPC::ArCF4N2_CF4_frac);
0395     tpc->set_double_param("N2_frac", G4TPC::ArCF4N2_N2_frac);
0396     tpc->set_double_param("isobutane_frac", G4TPC::ArCF4N2_isobutane_frac);
0397   }
0398   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4Isobutane")
0399   {
0400     tpc->set_double_param("Ne_frac", G4TPC::ArCF4Isobutane_Ne_frac);
0401     tpc->set_double_param("Ar_frac", G4TPC::ArCF4Isobutane_Ar_frac);
0402     tpc->set_double_param("CF4_frac", G4TPC::ArCF4Isobutane_CF4_frac);
0403     tpc->set_double_param("N2_frac", G4TPC::ArCF4Isobutane_N2_frac);
0404     tpc->set_double_param("isobutane_frac", G4TPC::ArCF4Isobutane_isobutane_frac);
0405   }
0406   else
0407   {
0408     std::cout << "Your gas mixture, " << G4TPC::TPC_GAS_MIXTURE << ", is unknown. Defaulting to ArCF4Isobutane" << std::endl;
0409   }
0410
0411   tpc->OverlapCheck(OverlapCheck);
0412
0413   std::cout << "PHG4TpcSubsystem z geometry parameters: " << std::endl;
0414   std::cout << "    drift_velocity_sim " << drift_vel << std::endl;
0415   std::cout << "    max_driftlength " << G4TPC::maxDriftLength << std::endl;
0416   std::cout << "    CM_halfwidth " << G4TPC::CM_halfwidth << std::endl;
0417   std::cout << "    pp_extended_readout_time " << TRACKING::pp_extended_readout_time << std::endl;
0418
0419   g4Reco->registerSubsystem(tpc);
0420
0421   if (Enable::TPC_ENDCAP)
0422   {
0423     TPC_Endcaps(g4Reco);
0424   }
0425
0426   radius = G4TPC::tpc_outer_radius;
0427
0428   radius += no_overlapp;
0429
0430   return radius;
0431 }
0432
0433 void TPC_Cells()
0434 {
0435   int verbosity = std::max(Enable::VERBOSITY, Enable::TPC_VERBOSITY);
0436   auto* se = Fun4AllServer::instance();
0437   /* // commented out later on
0438     double drift_vel = G4TPC::tpc_drift_velocity_sim;
0439     if (G4TPC::TPC_GAS_MIXTURE == "NeCF4")
0440     {
0441      drift_vel = G4TPC::NeCF4_drift_velocity;
0442     }
0443     else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4")
0444     {
0445      drift_vel = G4TPC::ArCF4_drift_velocity;
0446     }
0447     else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4N2")
0448     {
0449      drift_vel = G4TPC::ArCF4N2_drift_velocity;
0450     }
0451     else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4Isobutane")
0452     {
0453      drift_vel = G4TPC::ArCF4Isobutane_drift_velocity;
0454     }
0455     else
0456     {
0457     }
0458   */
0459   // central membrane G4Hit generation
0460   if (G4TPC::ENABLE_CENTRAL_MEMBRANE_HITS)
0461   {
0462     auto* centralMembrane = new PHG4TpcCentralMembrane;
0463     centralMembrane->setCentralMembraneDelay(0);
0464     centralMembrane->setCentralMembraneEventModulo(5);
0465     se->registerSubsystem(centralMembrane);
0466   }
0467
0468   // direct laser G4Hit generation
0469   if (G4TPC::ENABLE_DIRECT_LASER_HITS)
0470   {
0471     auto* directLaser = new PHG4TpcDirectLaser;
0472
0473     // setup phi and theta steps
0474     /* use 5deg steps */
0475     static constexpr double deg_to_rad = M_PI / 180.;
0476     directLaser->SetPhiStepping(144, 0 * deg_to_rad, 360 * deg_to_rad);
0477     directLaser->SetThetaStepping(36, 0 * deg_to_rad, 90 * deg_to_rad);
0478     // directLaser->SetArbitraryThetaPhi(50*deg_to_rad, 145*deg_to_rad);
0479     directLaser->SetDirectLaserAuto(true);
0480     //__Variable stepping: hitting all of the central membrane____________
0481     // directLaser->SetDirectLaserPatternfromFile( true );
0482     // directLaser->SetFileStepping(13802);
0483     //___________________________________________________________________
0484
0485     //    directLaser->set_double_param("drift_velocity", drift_vel);
0486     se->registerSubsystem(directLaser);
0487   }
0488
0489   //=========================
0490   // setup Tpc readout for filling cells
0491   // g4tpc/PHG4TpcElectronDrift uses
0492   // g4tpc/PHG4TpcPadPlaneReadout
0493   //=========================
0494
0495   auto* padplane = new PHG4TpcPadPlaneReadout;
0496   padplane->Verbosity(verbosity);
0497
0498   padplane->set_int_param("ntpc_phibins_inner", G4TPC::tpc_layer_rphi_count_inner);
0499   //  padplane->SetDriftVelocity(drift_vel);
0500
0501   auto* edrift = new PHG4TpcElectronDrift;
0502   edrift->Detector("TPC");
0503   edrift->Verbosity(verbosity);
0504   if (G4TPC::ENABLE_STATIC_DISTORTIONS || G4TPC::ENABLE_TIME_ORDERED_DISTORTIONS)
0505   {
0506     auto* distortionMap = new PHG4TpcDistortion;
0507
0508     distortionMap->set_read_phi_as_radians(G4TPC::DISTORTIONS_USE_PHI_AS_RADIANS);
0509
0510     distortionMap->set_do_static_distortions(G4TPC::ENABLE_STATIC_DISTORTIONS);
0511     if (isRootFile(G4TPC::static_distortion_filename))
0512     {
0513       distortionMap->set_static_distortion_filename(G4TPC::static_distortion_filename);
0514     }
0515     else
0516     {
0517       distortionMap->set_static_distortion_filename(CDBInterface::instance()->getUrl(G4TPC::static_distortion_filename));
0518     }
0519     if (isRootFile(G4TPC::time_ordered_distortion_filename))
0520     {
0521       distortionMap->set_time_ordered_distortion_filename(G4TPC::time_ordered_distortion_filename);
0522     }
0523     else
0524     {
0525       distortionMap->set_time_ordered_distortion_filename(CDBInterface::instance()->getUrl(G4TPC::time_ordered_distortion_filename));
0526     }
0527     distortionMap->set_do_ReachesReadout(G4TPC::ENABLE_REACHES_READOUT);
0528     distortionMap->Init();
0529     edrift->setTpcDistortion(distortionMap);
0530   }
0531
0532   edrift->set_double_param("added_smear_trans", 0.085);
0533   edrift->set_double_param("added_smear_long", 0.105);
0534
0535   // Note that we default to 75:20:05 Ar:CF4:i-C4H10
0536   if (G4TPC::TPC_GAS_MIXTURE == "NeCF4")
0537   {
0538     edrift->set_double_param("added_smear_trans", 0.085);
0539     edrift->set_double_param("added_smear_long", 0.105);
0540     edrift->set_double_param("diffusion_long", G4TPC::NeCF4_diffusion_long);
0541     edrift->set_double_param("diffusion_trans", G4TPC::NeCF4_diffusion_trans);
0542     edrift->set_double_param("Ne_frac", G4TPC::NeCF4_Ne_frac);
0543     edrift->set_double_param("Ar_frac", G4TPC::NeCF4_Ar_frac);
0544     edrift->set_double_param("CF4_frac", G4TPC::NeCF4_CF4_frac);
0545     edrift->set_double_param("N2_frac", G4TPC::NeCF4_N2_frac);
0546     edrift->set_double_param("isobutane_frac", G4TPC::NeCF4_isobutane_frac);
0547   }
0548   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4")
0549   {
0550     edrift->set_double_param("added_smear_trans", 0.085);
0551     edrift->set_double_param("added_smear_long", 0.105);
0552     edrift->set_double_param("diffusion_long", G4TPC::ArCF4_diffusion_long);
0553     edrift->set_double_param("diffusion_trans", G4TPC::ArCF4_diffusion_trans);
0554     edrift->set_double_param("Ne_frac", G4TPC::ArCF4_Ne_frac);
0555     edrift->set_double_param("Ar_frac", G4TPC::ArCF4_Ar_frac);
0556     edrift->set_double_param("CF4_frac", G4TPC::ArCF4_CF4_frac);
0557     edrift->set_double_param("N2_frac", G4TPC::ArCF4_N2_frac);
0558     edrift->set_double_param("isobutane_frac", G4TPC::ArCF4_isobutane_frac);
0559   }
0560   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4N2")
0561   {
0562     edrift->set_double_param("added_smear_trans", G4TPC::tpc_added_smear_trans);
0563     edrift->set_double_param("added_smear_long", G4TPC::tpc_added_smear_long);
0564     edrift->set_double_param("diffusion_long", G4TPC::ArCF4N2_diffusion_long);
0565     edrift->set_double_param("diffusion_trans", G4TPC::ArCF4N2_diffusion_trans);
0566     edrift->set_double_param("Ne_frac", G4TPC::ArCF4N2_Ne_frac);
0567     edrift->set_double_param("Ar_frac", G4TPC::ArCF4N2_Ar_frac);
0568     edrift->set_double_param("CF4_frac", G4TPC::ArCF4N2_CF4_frac);
0569     edrift->set_double_param("N2_frac", G4TPC::ArCF4N2_N2_frac);
0570     edrift->set_double_param("isobutane_frac", G4TPC::ArCF4N2_isobutane_frac);
0571   }
0572   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4Isobutane")
0573   {
0574     edrift->set_double_param("added_smear_trans", G4TPC::tpc_added_smear_trans);
0575     edrift->set_double_param("added_smear_long", G4TPC::tpc_added_smear_long);
0576     edrift->set_double_param("diffusion_long", G4TPC::ArCF4Isobutane_diffusion_long);
0577     edrift->set_double_param("diffusion_trans", G4TPC::ArCF4Isobutane_diffusion_trans);
0578     edrift->set_double_param("Ne_frac", G4TPC::ArCF4Isobutane_Ne_frac);
0579     edrift->set_double_param("Ar_frac", G4TPC::ArCF4Isobutane_Ar_frac);
0580     edrift->set_double_param("CF4_frac", G4TPC::ArCF4Isobutane_CF4_frac);
0581     edrift->set_double_param("N2_frac", G4TPC::ArCF4Isobutane_N2_frac);
0582     edrift->set_double_param("isobutane_frac", G4TPC::ArCF4Isobutane_isobutane_frac);
0583   }
0584   else
0585   {
0586   }
0587
0588   // fudge factors to get drphi 150 microns (in mid and outer Tpc) and dz 500 microns cluster resolution
0589   // They represent effects not due to ideal gas properties and ideal readout plane behavior
0590   // defaults are 0.085 and 0.105, they can be changed here to get a different resolution
0591
0592   edrift->registerPadPlane(padplane);
0593   se->registerSubsystem(edrift);
0594
0595   // Tpc digitizer
0596   //=========
0597   PHG4TpcDigitizer* digitpc = new PHG4TpcDigitizer();
0598   digitpc->SetTpcMinLayer(G4MVTX::n_maps_layer + G4INTT::n_intt_layer);
0599   double ENC = 670.0;  // standard
0600   digitpc->SetENC(ENC);
0601   double ADC_threshold = 4.0 * ENC;
0602   digitpc->SetADCThreshold(ADC_threshold);  // 4 * ENC seems OK
0603   digitpc->Verbosity(verbosity);
0604   std::cout << " Tpc digitizer: Setting ENC to " << ENC << " ADC threshold to " << ADC_threshold
0605             << " maps+Intt layers set to " << G4MVTX::n_maps_layer + G4INTT::n_intt_layer << std::endl;
0606   digitpc->set_skip_noise_flag(false);
0607   se->registerSubsystem(digitpc);
0608 }
0609
0610 void MicromegasInit()
0611 {
0612   if (!Enable::MVTX)
0613   {
0614     G4MVTX::n_maps_layer = 0;
0615   }
0616   // same for the INTT
0617   if (!Enable::INTT)
0618   {
0619     G4INTT::n_intt_layer = 0;
0620   }
0621   if (!Enable::TPC)
0622   {
0623     G4TPC::n_gas_layer = 0;
0624   }
0625   BlackHoleGeometry::max_radius = std::max(BlackHoleGeometry::max_radius, 88.);
0626   BlackHoleGeometry::max_z = std::max(BlackHoleGeometry::max_z, 220. / 2.);
0627   BlackHoleGeometry::min_z = std::min(BlackHoleGeometry::min_z, -220. / 2.);
0628 }
0629
0630 void Micromegas(PHG4Reco* g4Reco)
0631 {
0632   bool overlapcheck = Enable::OVERLAPCHECK || Enable::MICROMEGAS_OVERLAPCHECK;
0633   int verbosity = std::max(Enable::VERBOSITY, Enable::MICROMEGAS_VERBOSITY);
0634   bool SupportActive = Enable::SUPPORT || Enable::MICROMEGAS_SUPPORT;
0635   const int mm_layer = G4MVTX::n_maps_layer + G4INTT::n_intt_layer + G4TPC::n_gas_layer;
0636   auto* mm = new PHG4MicromegasSubsystem("MICROMEGAS", mm_layer);
0637   mm->Verbosity(verbosity);
0638   if (SupportActive)
0639   {
0640     mm->SetSupportActive();
0641   }
0642   mm->OverlapCheck(overlapcheck);
0643   mm->SetActive();
0644   g4Reco->registerSubsystem(mm);
0645 }
0646
0647 void Micromegas_Cells()
0648 {
0649   // the acts geometry needs to go here since it will be used by the PHG4MicromegasHitReco
0650   ACTSGEOM::ActsGeomInit();
0651   auto* se = Fun4AllServer::instance();
0652   int verbosity = std::max(Enable::VERBOSITY, Enable::MICROMEGAS_VERBOSITY);
0653   // micromegas
0654   auto* reco = new PHG4MicromegasHitReco;
0655   reco->Verbosity(verbosity);
0656   double extended_readout_time = 0.0;
0657   if (TRACKING::pp_mode) extended_readout_time = TRACKING::pp_extended_readout_time;
0658
0659   reco->set_double_param("micromegas_tmax", 800.0 + extended_readout_time);
0660   se->registerSubsystem(reco);
0661
0662   se->registerSubsystem(new PHG4MicromegasDigitizer);
0663 }
0664
0665 #endif