macros/common/G4_TrkrSimulation.C

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