macros/leptoquarks_fun4all/G4Setup_sPHENIX.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_pstof = true,
0007         bool do_cemc = true,
0008         bool do_hcalin = true,
0009         bool do_magnet = true,
0010         bool do_hcalout = true,
0011         bool do_pipe = true,
0012         int n_TPC_layers = 40)
0013   {
0014
0015   // load detector/material macros and execute Init() function
0016
0017   if (do_pipe)
0018     {
0019       gROOT->LoadMacro("G4_Pipe.C");
0020       PipeInit();
0021     }
0022   if (do_svtx)
0023     {
0024       gROOT->LoadMacro("G4_Svtx_maps_ladders+intt_ladders+tpc_KalmanPatRec.C");
0025       SvtxInit(n_TPC_layers);
0026     }
0027
0028   if (do_pstof)
0029     {
0030       gROOT->LoadMacro("G4_PSTOF.C");
0031       PSTOFInit();
0032     }
0033
0034   if (do_cemc)
0035     {
0036       gROOT->LoadMacro("G4_CEmc_Spacal.C");
0037       CEmcInit(72); // make it 2*2*2*3*3 so we can try other combinations
0038     }
0039
0040   if (do_hcalin)
0041     {
0042       gROOT->LoadMacro("G4_HcalIn_ref.C");
0043       HCalInnerInit();
0044     }
0045
0046   if (do_magnet)
0047     {
0048       gROOT->LoadMacro("G4_Magnet.C");
0049       MagnetInit();
0050     }
0051   if (do_hcalout)
0052     {
0053       gROOT->LoadMacro("G4_HcalOut_ref.C");
0054       HCalOuterInit();
0055     }
0056
0057 }
0058
0059
0060 int G4Setup(const int absorberactive = 0,
0061         const string &field ="1.5",
0062         const EDecayType decayType = TPythia6Decayer::kAll,
0063         const bool do_svtx = true,
0064         const bool do_pstof = true,
0065         const bool do_cemc = true,
0066         const bool do_hcalin = true,
0067         const bool do_magnet = true,
0068         const bool do_hcalout = true,
0069         const bool do_pipe = true,
0070         const float magfield_rescale = 1.0) {
0071
0072   //---------------
0073   // Load libraries
0074   //---------------
0075
0076   gSystem->Load("libg4detectors.so");
0077   gSystem->Load("libg4testbench.so");
0078
0079   //---------------
0080   // Fun4All server
0081   //---------------
0082
0083   Fun4AllServer *se = Fun4AllServer::instance();
0084
0085   // read-in HepMC events to Geant4 if there is any
0086   HepMCNodeReader *hr = new HepMCNodeReader();
0087   se->registerSubsystem(hr);
0088
0089   PHG4Reco* g4Reco = new PHG4Reco();
0090   g4Reco->set_rapidity_coverage(1.1); // according to drawings
0091 // uncomment to set QGSP_BERT_HP physics list for productions
0092 // (default is QGSP_BERT for speed)
0093   //  g4Reco->SetPhysicsList("QGSP_BERT_HP");
0094   if (decayType != TPythia6Decayer::kAll) {
0095     g4Reco->set_force_decay(decayType);
0096   }
0097
0098   double fieldstrength;
0099   istringstream stringline(field);
0100   stringline >> fieldstrength;
0101   if (stringline.fail()) { // conversion to double fails -> we have a string
0102
0103     if (field.find("sPHENIX.root") != string::npos) {
0104       g4Reco->set_field_map(field, 1);
0105     } else {
0106       g4Reco->set_field_map(field, 2);
0107     }
0108   } else {
0109     g4Reco->set_field(fieldstrength); // use const soleniodal field
0110   }
0111   g4Reco->set_field_rescale(magfield_rescale);
0112
0113   double radius = 0.;
0114
0115   //----------------------------------------
0116   // PIPE
0117   if (do_pipe) radius = Pipe(g4Reco, radius, absorberactive);
0118
0119   //----------------------------------------
0120   // SVTX
0121   if (do_svtx) radius = Svtx(g4Reco, radius, absorberactive);
0122
0123   //----------------------------------------
0124   // PSTOF
0125
0126   if (do_pstof) radius = PSTOF(g4Reco, radius, absorberactive);
0127
0128   //----------------------------------------
0129   // CEMC
0130 //
0131   if (do_cemc) radius = CEmc(g4Reco, radius, 8, absorberactive);
0132 //  if (do_cemc) radius = CEmc_Vis(g4Reco, radius, 8, absorberactive);// for visualization substructure of SPACAL, slow to render
0133
0134   //----------------------------------------
0135   // HCALIN
0136
0137   if (do_hcalin) radius = HCalInner(g4Reco, radius, 4, absorberactive);
0138
0139   //----------------------------------------
0140   // MAGNET
0141
0142   if (do_magnet) radius = Magnet(g4Reco, radius, 0, absorberactive);
0143
0144   //----------------------------------------
0145   // HCALOUT
0146
0147   if (do_hcalout) radius = HCalOuter(g4Reco, radius, 4, absorberactive);
0148
0149   //----------------------------------------
0150   // BLACKHOLE
0151
0152   // swallow all particles coming out of the backend of sPHENIX
0153   PHG4CylinderSubsystem *blackhole = new PHG4CylinderSubsystem("BH", 1);
0154 blackhole->set_double_param("radius",radius + 10); // add 10 cm
0155
0156   blackhole->set_int_param("lengthviarapidity",0);
0157   blackhole->set_double_param("length",g4Reco->GetWorldSizeZ() - no_overlapp); // make it cover the world in length
0158   blackhole->BlackHole();
0159   blackhole->set_double_param("thickness",0.1); // it needs some thickness
0160   blackhole->SetActive(); // always see what leaks out
0161   blackhole->OverlapCheck(overlapcheck);
0162   g4Reco->registerSubsystem(blackhole);
0163
0164   //----------------------------------------
0165   // FORWARD BLACKHOLEs
0166   // +Z
0167   blackhole = new PHG4CylinderSubsystem("BH_FORWARD_PLUS", 1);
0168   blackhole->SuperDetector("BH_FORWARD_PLUS");
0169   blackhole->set_double_param("radius",0); // add 10 cm
0170   blackhole->set_int_param("lengthviarapidity",0);
0171   blackhole->set_double_param("length",0.1); // make it cover the world in length
0172   blackhole->set_double_param("place_z",g4Reco->GetWorldSizeZ()/2. - 0.1  - no_overlapp);
0173   blackhole->BlackHole();
0174   blackhole->set_double_param("thickness",radius - no_overlapp); // it needs some thickness
0175   blackhole->SetActive(); // always see what leaks out
0176   blackhole->OverlapCheck(overlapcheck);
0177   g4Reco->registerSubsystem(blackhole);
0178
0179   blackhole = new PHG4CylinderSubsystem("BH_FORWARD_NEG", 1);
0180   blackhole->SuperDetector("BH_FORWARD_NEG");
0181   blackhole->set_double_param("radius",0); // add 10 cm
0182   blackhole->set_int_param("lengthviarapidity",0);
0183   blackhole->set_double_param("length",0.1); // make it cover the world in length
0184   blackhole->set_double_param("place_z", - g4Reco->GetWorldSizeZ()/2. +0.1  + no_overlapp);
0185   blackhole->BlackHole();
0186   blackhole->set_double_param("thickness",radius - no_overlapp); // it needs some thickness
0187   blackhole->SetActive(); // always see what leaks out
0188   blackhole->OverlapCheck(overlapcheck);
0189   g4Reco->registerSubsystem(blackhole);
0190
0191   PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
0192   g4Reco->registerSubsystem(truth);
0193   se->registerSubsystem( g4Reco );
0194 }
0195
0196 void ShowerCompress(int verbosity = 0) {
0197
0198   gSystem->Load("libfun4all.so");
0199   gSystem->Load("libg4eval.so");
0200
0201   Fun4AllServer *se = Fun4AllServer::instance();
0202
0203   PHG4DstCompressReco* compress = new PHG4DstCompressReco("PHG4DstCompressReco");
0204   compress->AddHitContainer("G4HIT_PIPE");
0205   compress->AddHitContainer("G4HIT_SVTXSUPPORT");
0206   compress->AddHitContainer("G4HIT_CEMC_ELECTRONICS");
0207   compress->AddHitContainer("G4HIT_CEMC");
0208   compress->AddHitContainer("G4HIT_ABSORBER_CEMC");
0209   compress->AddHitContainer("G4HIT_CEMC_SPT");
0210   compress->AddHitContainer("G4HIT_ABSORBER_HCALIN");
0211   compress->AddHitContainer("G4HIT_HCALIN");
0212   compress->AddHitContainer("G4HIT_HCALIN_SPT");
0213   compress->AddHitContainer("G4HIT_MAGNET");
0214   compress->AddHitContainer("G4HIT_ABSORBER_HCALOUT");
0215   compress->AddHitContainer("G4HIT_HCALOUT");
0216   compress->AddHitContainer("G4HIT_BH_1");
0217   compress->AddHitContainer("G4HIT_BH_FORWARD_PLUS");
0218   compress->AddHitContainer("G4HIT_BH_FORWARD_NEG");
0219   compress->AddCellContainer("G4CELL_CEMC");
0220   compress->AddCellContainer("G4CELL_HCALIN");
0221   compress->AddCellContainer("G4CELL_HCALOUT");
0222   compress->AddTowerContainer("TOWER_SIM_CEMC");
0223   compress->AddTowerContainer("TOWER_RAW_CEMC");
0224   compress->AddTowerContainer("TOWER_CALIB_CEMC");
0225   compress->AddTowerContainer("TOWER_SIM_HCALIN");
0226   compress->AddTowerContainer("TOWER_RAW_HCALIN");
0227   compress->AddTowerContainer("TOWER_CALIB_HCALIN");
0228   compress->AddTowerContainer("TOWER_SIM_HCALOUT");
0229   compress->AddTowerContainer("TOWER_RAW_HCALOUT");
0230   compress->AddTowerContainer("TOWER_CALIB_HCALOUT");
0231   se->registerSubsystem(compress);
0232
0233   return;
0234 }
0235
0236 void DstCompress(Fun4AllDstOutputManager* out) {
0237   if (out) {
0238     out->StripNode("G4HIT_PIPE");
0239     out->StripNode("G4HIT_SVTXSUPPORT");
0240     out->StripNode("G4HIT_CEMC_ELECTRONICS");
0241     out->StripNode("G4HIT_CEMC");
0242     out->StripNode("G4HIT_ABSORBER_CEMC");
0243     out->StripNode("G4HIT_CEMC_SPT");
0244     out->StripNode("G4HIT_ABSORBER_HCALIN");
0245     out->StripNode("G4HIT_HCALIN");
0246     out->StripNode("G4HIT_HCALIN_SPT");
0247     out->StripNode("G4HIT_MAGNET");
0248     out->StripNode("G4HIT_ABSORBER_HCALOUT");
0249     out->StripNode("G4HIT_HCALOUT");
0250     out->StripNode("G4HIT_BH_1");
0251     out->StripNode("G4HIT_BH_FORWARD_PLUS");
0252     out->StripNode("G4HIT_BH_FORWARD_NEG");
0253     out->StripNode("G4CELL_CEMC");
0254     out->StripNode("G4CELL_HCALIN");
0255     out->StripNode("G4CELL_HCALOUT");
0256   }
0257 }