macros/leptoquarks_fun4all/G4Setup_EICDetector.C

0001
0002 double no_overlapp = 0.0001; // added to radii to avoid overlapping volumes
0003 bool overlapcheck = false; // set to true if you want to check for overlaps
0004
0005 void G4Init(bool do_svtx = true,
0006             bool do_cemc = true,
0007             bool do_hcalin = true,
0008             bool do_magnet = true,
0009             bool do_hcalout = true,
0010             bool do_pipe = true,
0011             bool do_FGEM = true,
0012             bool do_EGEM = true,
0013             bool do_FEMC = true,
0014             bool do_FHCAL = true,
0015             bool do_EEMC = true,
0016             bool do_DIRC = true,
0017             bool do_RICH = true,
0018             bool do_Aerogel = true,
0019             int n_TPC_layers = 40)
0020 {
0021
0022   // load detector/material macros and execute Init() function
0023
0024   if (do_pipe)
0025     {
0026       gROOT->LoadMacro("G4_Pipe.C");
0027       PipeInit();
0028     }
0029   if (do_svtx)
0030     {
0031       gROOT->LoadMacro("G4_Svtx_maps_ladders+intt_ladders+tpc_KalmanPatRec.C");
0032       SvtxInit(n_TPC_layers);
0033     }
0034
0035   if (do_cemc)
0036     {
0037       gROOT->LoadMacro("G4_CEmc_Spacal.C");
0038       CEmcInit(72); // make it 2*2*2*3*3 so we can try other combinations
0039     }
0040
0041   if (do_hcalin)
0042     {
0043       gROOT->LoadMacro("G4_HcalIn_ref.C");
0044       HCalInnerInit();
0045     }
0046
0047   if (do_magnet)
0048     {
0049       gROOT->LoadMacro("G4_Magnet.C");
0050       MagnetInit();
0051     }
0052   if (do_hcalout)
0053     {
0054       gROOT->LoadMacro("G4_HcalOut_ref.C");
0055       HCalOuterInit();
0056     }
0057
0058   if (do_FGEM)
0059     {
0060       gROOT->LoadMacro("G4_FGEM_EIC.C");
0061       FGEM_Init();
0062     }
0063
0064   if (do_EGEM)
0065     {
0066       gROOT->LoadMacro("G4_EGEM_EIC.C");
0067       EGEM_Init();
0068     }
0069
0070   if (do_FEMC)
0071     {
0072       gROOT->LoadMacro("G4_FEMC.C");
0073       FEMCInit();
0074     }
0075
0076   if (do_FHCAL)
0077     {
0078       gROOT->LoadMacro("G4_FHCAL.C");
0079       FHCALInit();
0080     }
0081
0082   if (do_EEMC)
0083     {
0084       gROOT->LoadMacro("G4_EEMC.C");
0085       EEMCInit();
0086     }
0087
0088   if (do_DIRC)
0089     {
0090       gROOT->LoadMacro("G4_DIRC.C");
0091       DIRCInit();
0092     }
0093
0094   if (do_RICH)
0095     {
0096       gROOT->LoadMacro("G4_RICH.C");
0097       RICHInit();
0098     }
0099
0100   if (do_Aerogel)
0101     {
0102       gROOT->LoadMacro("G4_Aerogel.C");
0103       AerogelInit();
0104     }
0105
0106
0107 }
0108
0109
0110 int G4Setup(const int absorberactive = 0,
0111             const string &field ="1.5",
0112             const EDecayType decayType = TPythia6Decayer::kAll,
0113             const bool do_svtx = true,
0114             const bool do_cemc = true,
0115             const bool do_hcalin = true,
0116             const bool do_magnet = true,
0117             const bool do_hcalout = true,
0118             const bool do_pipe = true,
0119             const bool do_FGEM = true,
0120             const bool do_EGEM = true,
0121             const bool do_FEMC = false,
0122             const bool do_FHCAL = false,
0123             const bool do_EEMC = true,
0124             const bool do_DIRC = true,
0125             const bool do_RICH = true,
0126             const bool do_Aerogel = true,
0127             const float magfield_rescale = 1.0) {
0128
0129   //---------------
0130   // Load libraries
0131   //---------------
0132
0133   gSystem->Load("libg4detectors.so");
0134   gSystem->Load("libg4testbench.so");
0135
0136   //---------------
0137   // Fun4All server
0138   //---------------
0139
0140   Fun4AllServer *se = Fun4AllServer::instance();
0141
0142   // read-in HepMC events to Geant4 if there is any
0143   HepMCNodeReader *hr = new HepMCNodeReader();
0144   se->registerSubsystem(hr);
0145
0146   PHG4Reco* g4Reco = new PHG4Reco();
0147   g4Reco->set_rapidity_coverage(1.1); // according to drawings
0148 // uncomment to set QGSP_BERT_HP physics list for productions
0149 // (default is QGSP_BERT for speed)
0150   //  g4Reco->SetPhysicsList("QGSP_BERT_HP");
0151
0152   if (decayType != TPythia6Decayer::kAll) {
0153     g4Reco->set_force_decay(decayType);
0154   }
0155
0156   double fieldstrength;
0157   istringstream stringline(field);
0158   stringline >> fieldstrength;
0159   if (stringline.fail()) { // conversion to double fails -> we have a string
0160
0161     if (field.find("sPHENIX.root") != string::npos) {
0162       g4Reco->set_field_map(field, 1);
0163     } else {
0164       g4Reco->set_field_map(field, 2);
0165     }
0166   } else {
0167     g4Reco->set_field(fieldstrength); // use const soleniodal field
0168   }
0169   g4Reco->set_field_rescale(magfield_rescale);
0170
0171   double radius = 0.;
0172
0173   //----------------------------------------
0174   // PIPE
0175   if (do_pipe) radius = Pipe(g4Reco, radius, absorberactive);
0176
0177   //----------------------------------------
0178   // SVTX
0179   if (do_svtx) radius = Svtx(g4Reco, radius, absorberactive);
0180
0181   //----------------------------------------
0182   // CEMC
0183   //
0184   if (do_cemc) radius = CEmc(g4Reco, radius, 8, absorberactive);
0185   //  if (do_cemc) radius = CEmc_Vis(g4Reco, radius, 8, absorberactive);// for visualization substructure of SPACAL, slow to render
0186
0187   //----------------------------------------
0188   // HCALIN
0189
0190   if (do_hcalin) radius = HCalInner(g4Reco, radius, 4, absorberactive);
0191
0192   //----------------------------------------
0193   // MAGNET
0194
0195   if (do_magnet) radius = Magnet(g4Reco, radius, 0, absorberactive);
0196
0197   //----------------------------------------
0198   // HCALOUT
0199
0200   if (do_hcalout) radius = HCalOuter(g4Reco, radius, 4, absorberactive);
0201
0202   //----------------------------------------
0203   // Forward tracking
0204
0205   if ( do_FGEM )
0206     FGEMSetup(g4Reco);
0207
0208   if ( do_EGEM )
0209     EGEMSetup(g4Reco);
0210
0211   //----------------------------------------
0212   // FEMC
0213
0214   if ( do_FEMC )
0215     FEMCSetup(g4Reco, absorberactive);
0216
0217   //----------------------------------------
0218   // FHCAL
0219
0220   if ( do_FHCAL )
0221     FHCALSetup(g4Reco, absorberactive);
0222
0223   //----------------------------------------
0224   // EEMC
0225
0226   if ( do_EEMC )
0227     EEMCSetup(g4Reco, absorberactive);
0228
0229   //----------------------------------------
0230   // PID
0231
0232   if ( do_DIRC )
0233     DIRCSetup(g4Reco);
0234
0235   if ( do_RICH )
0236     RICHSetup(g4Reco);
0237
0238   if ( do_Aerogel )
0239     AerogelSetup(g4Reco);
0240
0241   // sPHENIX forward flux return(s)
0242   PHG4CylinderSubsystem *flux_return_plus = new PHG4CylinderSubsystem("FWDFLUXRET", 0);
0243   flux_return_plus->set_int_param("lengthviarapidity",0);
0244   flux_return_plus->set_double_param("length",10.2);
0245   flux_return_plus->set_double_param("radius",2.1);
0246   flux_return_plus->set_double_param("thickness",263.5-5.0);
0247   flux_return_plus->set_double_param("place_z",335.9);
0248   flux_return_plus->set_string_param("material","G4_Fe");
0249   flux_return_plus->SetActive(false);
0250   flux_return_plus->SuperDetector("FLUXRET_ETA_PLUS");
0251   flux_return_plus->OverlapCheck(overlapcheck);
0252   g4Reco->registerSubsystem(flux_return_plus);
0253
0254   PHG4CylinderSubsystem *flux_return_minus = new PHG4CylinderSubsystem("FWDFLUXRET", 0);
0255   flux_return_minus->set_int_param("lengthviarapidity",0);
0256   flux_return_minus->set_double_param("length",10.2);
0257   flux_return_minus->set_double_param("radius",90.0);
0258   flux_return_minus->set_double_param("place_z",-335.9);
0259   flux_return_minus->set_double_param("thickness",263.5-5.0 - (90-2.1));
0260   flux_return_minus->set_string_param("material","G4_Fe");
0261   flux_return_minus->SetActive(false);
0262   flux_return_minus->SuperDetector("FLUXRET_ETA_MINUS");
0263   flux_return_minus->OverlapCheck(overlapcheck);
0264   g4Reco->registerSubsystem(flux_return_minus);
0265
0266   //----------------------------------------
0267   // BLACKHOLE
0268
0269   // swallow all particles coming out of the backend of sPHENIX
0270   PHG4CylinderSubsystem *blackhole = new PHG4CylinderSubsystem("BH", 1);
0271   blackhole->set_double_param("radius",radius + 100); // add 100 cm
0272
0273   blackhole->set_int_param("lengthviarapidity",0);
0274   blackhole->set_double_param("length",g4Reco->GetWorldSizeZ() - no_overlapp); // make it cover the world in length
0275   blackhole->BlackHole();
0276   blackhole->set_double_param("thickness",0.1); // it needs some thickness
0277   blackhole->SetActive(); // always see what leaks out
0278   blackhole->OverlapCheck(overlapcheck);
0279   g4Reco->registerSubsystem(blackhole);
0280
0281   //----------------------------------------
0282   // FORWARD BLACKHOLEs
0283   // +Z
0284   blackhole = new PHG4CylinderSubsystem("BH_FORWARD_PLUS", 1);
0285   blackhole->SuperDetector("BH_FORWARD_PLUS");
0286   blackhole->set_double_param("radius",0); // add 10 cm
0287   blackhole->set_int_param("lengthviarapidity",0);
0288   blackhole->set_double_param("length",0.1); // make it cover the world in length
0289   blackhole->set_double_param("place_z",g4Reco->GetWorldSizeZ()/2. - 0.1  - no_overlapp);
0290   blackhole->BlackHole();
0291   blackhole->set_double_param("thickness",radius - no_overlapp); // it needs some thickness
0292   blackhole->SetActive(); // always see what leaks out
0293   blackhole->OverlapCheck(overlapcheck);
0294   g4Reco->registerSubsystem(blackhole);
0295
0296   blackhole = new PHG4CylinderSubsystem("BH_FORWARD_NEG", 1);
0297   blackhole->SuperDetector("BH_FORWARD_NEG");
0298   blackhole->set_double_param("radius",0); // add 10 cm
0299   blackhole->set_int_param("lengthviarapidity",0);
0300   blackhole->set_double_param("length",0.1); // make it cover the world in length
0301   blackhole->set_double_param("place_z", - g4Reco->GetWorldSizeZ()/2. +0.1  + no_overlapp);
0302   blackhole->BlackHole();
0303   blackhole->set_double_param("thickness",radius - no_overlapp); // it needs some thickness
0304   blackhole->SetActive(); // always see what leaks out
0305   blackhole->OverlapCheck(overlapcheck);
0306   g4Reco->registerSubsystem(blackhole);
0307
0308   PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
0309   g4Reco->registerSubsystem(truth);
0310   se->registerSubsystem( g4Reco );
0311 }
0312
0313
0314 void ShowerCompress(int verbosity = 0) {
0315
0316   gSystem->Load("libfun4all.so");
0317   gSystem->Load("libg4eval.so");
0318
0319   Fun4AllServer *se = Fun4AllServer::instance();
0320
0321   PHG4DstCompressReco* compress = new PHG4DstCompressReco("PHG4DstCompressReco");
0322   compress->AddHitContainer("G4HIT_PIPE");
0323   compress->AddHitContainer("G4HIT_SVTXSUPPORT");
0324   compress->AddHitContainer("G4HIT_CEMC_ELECTRONICS");
0325   compress->AddHitContainer("G4HIT_CEMC");
0326   compress->AddHitContainer("G4HIT_ABSORBER_CEMC");
0327   compress->AddHitContainer("G4HIT_CEMC_SPT");
0328   compress->AddHitContainer("G4HIT_ABSORBER_HCALIN");
0329   compress->AddHitContainer("G4HIT_HCALIN");
0330   compress->AddHitContainer("G4HIT_HCALIN_SPT");
0331   compress->AddHitContainer("G4HIT_MAGNET");
0332   compress->AddHitContainer("G4HIT_ABSORBER_HCALOUT");
0333   compress->AddHitContainer("G4HIT_HCALOUT");
0334   compress->AddHitContainer("G4HIT_BH_1");
0335   compress->AddHitContainer("G4HIT_BH_FORWARD_PLUS");
0336   compress->AddHitContainer("G4HIT_BH_FORWARD_NEG");
0337   compress->AddCellContainer("G4CELL_CEMC");
0338   compress->AddCellContainer("G4CELL_HCALIN");
0339   compress->AddCellContainer("G4CELL_HCALOUT");
0340   compress->AddTowerContainer("TOWER_SIM_CEMC");
0341   compress->AddTowerContainer("TOWER_RAW_CEMC");
0342   compress->AddTowerContainer("TOWER_CALIB_CEMC");
0343   compress->AddTowerContainer("TOWER_SIM_HCALIN");
0344   compress->AddTowerContainer("TOWER_RAW_HCALIN");
0345   compress->AddTowerContainer("TOWER_CALIB_HCALIN");
0346   compress->AddTowerContainer("TOWER_SIM_HCALOUT");
0347   compress->AddTowerContainer("TOWER_RAW_HCALOUT");
0348   compress->AddTowerContainer("TOWER_CALIB_HCALOUT");
0349
0350   compress->AddHitContainer("G4HIT_FEMC");
0351   compress->AddHitContainer("G4HIT_ABSORBER_FEMC");
0352   compress->AddHitContainer("G4HIT_FHCAL");
0353   compress->AddHitContainer("G4HIT_ABSORBER_FHCAL");
0354   compress->AddCellContainer("G4CELL_FEMC");
0355   compress->AddCellContainer("G4CELL_FHCAL");
0356   compress->AddTowerContainer("TOWER_SIM_FEMC");
0357   compress->AddTowerContainer("TOWER_RAW_FEMC");
0358   compress->AddTowerContainer("TOWER_CALIB_FEMC");
0359   compress->AddTowerContainer("TOWER_SIM_FHCAL");
0360   compress->AddTowerContainer("TOWER_RAW_FHCAL");
0361   compress->AddTowerContainer("TOWER_CALIB_FHCAL");
0362
0363   compress->AddHitContainer("G4HIT_EEMC");
0364   compress->AddHitContainer("G4HIT_ABSORBER_EEMC");
0365   compress->AddCellContainer("G4CELL_EEMC");
0366   compress->AddTowerContainer("TOWER_SIM_EEMC");
0367   compress->AddTowerContainer("TOWER_RAW_EEMC");
0368   compress->AddTowerContainer("TOWER_CALIB_EEMC");
0369
0370   se->registerSubsystem(compress);
0371
0372   return;
0373 }
0374
0375 void DstCompress(Fun4AllDstOutputManager* out) {
0376   if (out) {
0377     out->StripNode("G4HIT_PIPE");
0378     out->StripNode("G4HIT_SVTXSUPPORT");
0379     out->StripNode("G4HIT_CEMC_ELECTRONICS");
0380     out->StripNode("G4HIT_CEMC");
0381     out->StripNode("G4HIT_ABSORBER_CEMC");
0382     out->StripNode("G4HIT_CEMC_SPT");
0383     out->StripNode("G4HIT_ABSORBER_HCALIN");
0384     out->StripNode("G4HIT_HCALIN");
0385     out->StripNode("G4HIT_HCALIN_SPT");
0386     out->StripNode("G4HIT_MAGNET");
0387     out->StripNode("G4HIT_ABSORBER_HCALOUT");
0388     out->StripNode("G4HIT_HCALOUT");
0389     out->StripNode("G4HIT_BH_1");
0390     out->StripNode("G4HIT_BH_FORWARD_PLUS");
0391     out->StripNode("G4HIT_BH_FORWARD_NEG");
0392     out->StripNode("G4CELL_CEMC");
0393     out->StripNode("G4CELL_HCALIN");
0394     out->StripNode("G4CELL_HCALOUT");
0395
0396     out->StripNode("G4HIT_FEMC");
0397     out->StripNode("G4HIT_ABSORBER_FEMC");
0398     out->StripNode("G4HIT_FHCAL");
0399     out->StripNode("G4HIT_ABSORBER_FHCAL");
0400     out->StripNode("G4CELL_FEMC");
0401     out->StripNode("G4CELL_FHCAL");
0402
0403     out->StripNode("G4HIT_EEMC");
0404     out->StripNode("G4HIT_ABSORBER_EEMC");
0405     out->StripNode("G4CELL_EEMC");
0406   }
0407 }