RICHAnalysis_macros/macros_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   PHG4Reco* g4Reco = new PHG4Reco();
0143   g4Reco->set_rapidity_coverage(1.1); // according to drawings
0144 // uncomment to set QGSP_BERT_HP physics list for productions
0145 // (default is QGSP_BERT for speed)
0146   //  g4Reco->SetPhysicsList("QGSP_BERT_HP");
0147
0148   if (decayType != TPythia6Decayer::kAll) {
0149     g4Reco->set_force_decay(decayType);
0150   }
0151
0152   double fieldstrength;
0153   istringstream stringline(field);
0154   stringline >> fieldstrength;
0155   if (stringline.fail()) { // conversion to double fails -> we have a string
0156
0157     if (field.find("sPHENIX.root") != string::npos) {
0158       g4Reco->set_field_map(field, 1);
0159     } else {
0160       g4Reco->set_field_map(field, 2);
0161     }
0162   } else {
0163     g4Reco->set_field(fieldstrength); // use const soleniodal field
0164   }
0165   g4Reco->set_field_rescale(magfield_rescale);
0166
0167   double radius = 0.;
0168
0169   //----------------------------------------
0170   // PIPE
0171   if (do_pipe) radius = Pipe(g4Reco, radius, absorberactive);
0172
0173   //----------------------------------------
0174   // SVTX
0175   if (do_svtx) radius = Svtx(g4Reco, radius, absorberactive);
0176
0177   //----------------------------------------
0178   // CEMC
0179   //
0180   if (do_cemc) radius = CEmc(g4Reco, radius, 8, absorberactive);
0181   //  if (do_cemc) radius = CEmc_Vis(g4Reco, radius, 8, absorberactive);// for visualization substructure of SPACAL, slow to render
0182
0183   //----------------------------------------
0184   // HCALIN
0185
0186   if (do_hcalin) radius = HCalInner(g4Reco, radius, 4, absorberactive);
0187
0188   //----------------------------------------
0189   // MAGNET
0190
0191   if (do_magnet) radius = Magnet(g4Reco, radius, 0, absorberactive);
0192
0193   //----------------------------------------
0194   // HCALOUT
0195
0196   if (do_hcalout) radius = HCalOuter(g4Reco, radius, 4, absorberactive);
0197
0198   //----------------------------------------
0199   // Forward tracking
0200
0201   if ( do_FGEM )
0202     FGEMSetup(g4Reco);
0203
0204   if ( do_EGEM )
0205     EGEMSetup(g4Reco);
0206
0207   //----------------------------------------
0208   // FEMC
0209
0210   if ( do_FEMC )
0211     FEMCSetup(g4Reco, absorberactive);
0212
0213   //----------------------------------------
0214   // FHCAL
0215
0216   if ( do_FHCAL )
0217     FHCALSetup(g4Reco, absorberactive);
0218
0219   //----------------------------------------
0220   // EEMC
0221
0222   if ( do_EEMC )
0223     EEMCSetup(g4Reco, absorberactive);
0224
0225   //----------------------------------------
0226   // PID
0227
0228   if ( do_DIRC )
0229     DIRCSetup(g4Reco);
0230
0231   if ( do_RICH )
0232     RICHSetup(g4Reco);
0233
0234   if ( do_Aerogel )
0235     AerogelSetup(g4Reco);
0236
0237   // sPHENIX forward flux return(s)
0238   PHG4CylinderSubsystem *flux_return_plus = new PHG4CylinderSubsystem("FWDFLUXRET", 0);
0239   flux_return_plus->set_int_param("lengthviarapidity",0);
0240   flux_return_plus->set_double_param("length",10.2);
0241   flux_return_plus->set_double_param("radius",2.1);
0242   flux_return_plus->set_double_param("thickness",263.5-5.0);
0243   flux_return_plus->set_double_param("place_z",335.9);
0244   flux_return_plus->set_string_param("material","G4_Fe");
0245   flux_return_plus->SetActive(false);
0246   flux_return_plus->SuperDetector("FLUXRET_ETA_PLUS");
0247   flux_return_plus->OverlapCheck(overlapcheck);
0248   g4Reco->registerSubsystem(flux_return_plus);
0249
0250   PHG4CylinderSubsystem *flux_return_minus = new PHG4CylinderSubsystem("FWDFLUXRET", 0);
0251   flux_return_minus->set_int_param("lengthviarapidity",0);
0252   flux_return_minus->set_double_param("length",10.2);
0253   flux_return_minus->set_double_param("radius",90.0);
0254   flux_return_minus->set_double_param("place_z",-335.9);
0255   flux_return_minus->set_double_param("thickness",263.5-5.0 - (90-2.1));
0256   flux_return_minus->set_string_param("material","G4_Fe");
0257   flux_return_minus->SetActive(false);
0258   flux_return_minus->SuperDetector("FLUXRET_ETA_MINUS");
0259   flux_return_minus->OverlapCheck(overlapcheck);
0260   g4Reco->registerSubsystem(flux_return_minus);
0261
0262   //----------------------------------------
0263   // BLACKHOLE
0264
0265   // swallow all particles coming out of the backend of sPHENIX
0266   PHG4CylinderSubsystem *blackhole = new PHG4CylinderSubsystem("BH", 1);
0267   blackhole->set_double_param("radius",radius + 100); // add 100 cm
0268
0269   blackhole->set_int_param("lengthviarapidity",0);
0270   blackhole->set_double_param("length",g4Reco->GetWorldSizeZ() - no_overlapp); // make it cover the world in length
0271   blackhole->BlackHole();
0272   blackhole->set_double_param("thickness",0.1); // it needs some thickness
0273   blackhole->SetActive(); // always see what leaks out
0274   blackhole->OverlapCheck(overlapcheck);
0275   g4Reco->registerSubsystem(blackhole);
0276
0277   //----------------------------------------
0278   // FORWARD BLACKHOLEs
0279   // +Z
0280   blackhole = new PHG4CylinderSubsystem("BH_FORWARD_PLUS", 1);
0281   blackhole->SuperDetector("BH_FORWARD_PLUS");
0282   blackhole->set_double_param("radius",0); // add 10 cm
0283   blackhole->set_int_param("lengthviarapidity",0);
0284   blackhole->set_double_param("length",0.1); // make it cover the world in length
0285   blackhole->set_double_param("place_z",g4Reco->GetWorldSizeZ()/2. - 0.1  - no_overlapp);
0286   blackhole->BlackHole();
0287   blackhole->set_double_param("thickness",radius - no_overlapp); // it needs some thickness
0288   blackhole->SetActive(); // always see what leaks out
0289   blackhole->OverlapCheck(overlapcheck);
0290   g4Reco->registerSubsystem(blackhole);
0291
0292   blackhole = new PHG4CylinderSubsystem("BH_FORWARD_NEG", 1);
0293   blackhole->SuperDetector("BH_FORWARD_NEG");
0294   blackhole->set_double_param("radius",0); // add 10 cm
0295   blackhole->set_int_param("lengthviarapidity",0);
0296   blackhole->set_double_param("length",0.1); // make it cover the world in length
0297   blackhole->set_double_param("place_z", - g4Reco->GetWorldSizeZ()/2. +0.1  + no_overlapp);
0298   blackhole->BlackHole();
0299   blackhole->set_double_param("thickness",radius - no_overlapp); // it needs some thickness
0300   blackhole->SetActive(); // always see what leaks out
0301   blackhole->OverlapCheck(overlapcheck);
0302   g4Reco->registerSubsystem(blackhole);
0303
0304   PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
0305   g4Reco->registerSubsystem(truth);
0306   se->registerSubsystem( g4Reco );
0307 }
0308
0309
0310 void ShowerCompress(int verbosity = 0) {
0311
0312   gSystem->Load("libfun4all.so");
0313   gSystem->Load("libg4eval.so");
0314
0315   Fun4AllServer *se = Fun4AllServer::instance();
0316
0317   PHG4DstCompressReco* compress = new PHG4DstCompressReco("PHG4DstCompressReco");
0318   compress->AddHitContainer("G4HIT_PIPE");
0319   compress->AddHitContainer("G4HIT_SVTXSUPPORT");
0320   compress->AddHitContainer("G4HIT_CEMC_ELECTRONICS");
0321   compress->AddHitContainer("G4HIT_CEMC");
0322   compress->AddHitContainer("G4HIT_ABSORBER_CEMC");
0323   compress->AddHitContainer("G4HIT_CEMC_SPT");
0324   compress->AddHitContainer("G4HIT_ABSORBER_HCALIN");
0325   compress->AddHitContainer("G4HIT_HCALIN");
0326   compress->AddHitContainer("G4HIT_HCALIN_SPT");
0327   compress->AddHitContainer("G4HIT_MAGNET");
0328   compress->AddHitContainer("G4HIT_ABSORBER_HCALOUT");
0329   compress->AddHitContainer("G4HIT_HCALOUT");
0330   compress->AddHitContainer("G4HIT_BH_1");
0331   compress->AddHitContainer("G4HIT_BH_FORWARD_PLUS");
0332   compress->AddHitContainer("G4HIT_BH_FORWARD_NEG");
0333   compress->AddCellContainer("G4CELL_CEMC");
0334   compress->AddCellContainer("G4CELL_HCALIN");
0335   compress->AddCellContainer("G4CELL_HCALOUT");
0336   compress->AddTowerContainer("TOWER_SIM_CEMC");
0337   compress->AddTowerContainer("TOWER_RAW_CEMC");
0338   compress->AddTowerContainer("TOWER_CALIB_CEMC");
0339   compress->AddTowerContainer("TOWER_SIM_HCALIN");
0340   compress->AddTowerContainer("TOWER_RAW_HCALIN");
0341   compress->AddTowerContainer("TOWER_CALIB_HCALIN");
0342   compress->AddTowerContainer("TOWER_SIM_HCALOUT");
0343   compress->AddTowerContainer("TOWER_RAW_HCALOUT");
0344   compress->AddTowerContainer("TOWER_CALIB_HCALOUT");
0345
0346   compress->AddHitContainer("G4HIT_FEMC");
0347   compress->AddHitContainer("G4HIT_ABSORBER_FEMC");
0348   compress->AddHitContainer("G4HIT_FHCAL");
0349   compress->AddHitContainer("G4HIT_ABSORBER_FHCAL");
0350   compress->AddCellContainer("G4CELL_FEMC");
0351   compress->AddCellContainer("G4CELL_FHCAL");
0352   compress->AddTowerContainer("TOWER_SIM_FEMC");
0353   compress->AddTowerContainer("TOWER_RAW_FEMC");
0354   compress->AddTowerContainer("TOWER_CALIB_FEMC");
0355   compress->AddTowerContainer("TOWER_SIM_FHCAL");
0356   compress->AddTowerContainer("TOWER_RAW_FHCAL");
0357   compress->AddTowerContainer("TOWER_CALIB_FHCAL");
0358
0359   compress->AddHitContainer("G4HIT_EEMC");
0360   compress->AddHitContainer("G4HIT_ABSORBER_EEMC");
0361   compress->AddCellContainer("G4CELL_EEMC");
0362   compress->AddTowerContainer("TOWER_SIM_EEMC");
0363   compress->AddTowerContainer("TOWER_RAW_EEMC");
0364   compress->AddTowerContainer("TOWER_CALIB_EEMC");
0365
0366   se->registerSubsystem(compress);
0367
0368   return;
0369 }
0370
0371 void DstCompress(Fun4AllDstOutputManager* out) {
0372   if (out) {
0373     out->StripNode("G4HIT_PIPE");
0374     out->StripNode("G4HIT_SVTXSUPPORT");
0375     out->StripNode("G4HIT_CEMC_ELECTRONICS");
0376     out->StripNode("G4HIT_CEMC");
0377     out->StripNode("G4HIT_ABSORBER_CEMC");
0378     out->StripNode("G4HIT_CEMC_SPT");
0379     out->StripNode("G4HIT_ABSORBER_HCALIN");
0380     out->StripNode("G4HIT_HCALIN");
0381     out->StripNode("G4HIT_HCALIN_SPT");
0382     out->StripNode("G4HIT_MAGNET");
0383     out->StripNode("G4HIT_ABSORBER_HCALOUT");
0384     out->StripNode("G4HIT_HCALOUT");
0385     out->StripNode("G4HIT_BH_1");
0386     out->StripNode("G4HIT_BH_FORWARD_PLUS");
0387     out->StripNode("G4HIT_BH_FORWARD_NEG");
0388     out->StripNode("G4CELL_CEMC");
0389     out->StripNode("G4CELL_HCALIN");
0390     out->StripNode("G4CELL_HCALOUT");
0391
0392     out->StripNode("G4HIT_FEMC");
0393     out->StripNode("G4HIT_ABSORBER_FEMC");
0394     out->StripNode("G4HIT_FHCAL");
0395     out->StripNode("G4HIT_ABSORBER_FHCAL");
0396     out->StripNode("G4CELL_FEMC");
0397     out->StripNode("G4CELL_FHCAL");
0398
0399     out->StripNode("G4HIT_EEMC");
0400     out->StripNode("G4HIT_ABSORBER_EEMC");
0401     out->StripNode("G4CELL_EEMC");
0402   }
0403 }