macros/gen/Fun4All_G4_Calo.C

0001 #ifndef MACRO_FUN4ALLG4CALO_C
0002 #define MACRO_FUN4ALLG4CALO_C
0003
0004 #include <GlobalVariables.C>
0005 #include <../../pi0Efficiency.h>
0006
0007 #include <DisplayOn.C>
0008 #include <G4Setup_sPHENIX.C>
0009 #include <G4_Bbc.C>
0010 #include <G4_CaloTrigger.C>
0011 #include <G4_DSTReader.C>
0012 #include <G4_Global.C>
0013 #include <G4_HIJetReco.C>
0014 #include <G4_Input.C>
0015 #include <G4_Jets.C>
0016 #include <G4_ParticleFlow.C>
0017 #include <G4_Production.C>
0018 #include <G4_TopoClusterReco.C>
0019 #include <G4_Tracking.C>
0020 #include <G4_User.C>
0021
0022 #include <ffamodules/FlagHandler.h>
0023 #include <ffamodules/XploadInterface.h>
0024
0025 #include <fun4all/Fun4AllDstOutputManager.h>
0026 #include <fun4all/Fun4AllOutputManager.h>
0027 #include <fun4all/Fun4AllServer.h>
0028
0029 #include <phool/PHRandomSeed.h>
0030 #include <phool/recoConsts.h>
0031
0032 R__LOAD_LIBRARY(libffamodules.so)
0033 R__LOAD_LIBRARY(libfun4all.so)
0034 R__LOAD_LIBRARY(libpi0Efficiency.so)
0035
0036 // For HepMC Hijing
0037 // try inputFile = /sphenix/sim/sim01/sphnxpro/sHijing_HepMC/sHijing_0-12fm.dat
0038
0039 int Fun4All_G4_Calo(
0040     const int nEvents = 1,
0041     const string &inputFile0 = "g4hits.list",
0042     //const string &inputFile1 = "DST_VERTEX_sHijing_0_12fm-0000000001-00000.root",
0043     const string &outputFile = "G4sPHENIX_calo.root",
0044     const string &outdir = ".")
0045 {
0046   Fun4AllServer *se = Fun4AllServer::instance();
0047   se->Verbosity(0);
0048
0049   //Opt to print all random seed used for debugging reproducibility. Comment out to reduce stdout prints.
0050   PHRandomSeed::Verbosity(0);
0051
0052   // just if we set some flags somewhere in this macro
0053   recoConsts *rc = recoConsts::instance();
0054   // By default every random number generator uses
0055   // PHRandomSeed() which reads /dev/urandom to get its seed
0056   // if the RANDOMSEED flag is set its value is taken as seed
0057   // You can either set this to a random value using PHRandomSeed()
0058   // which will make all seeds identical (not sure what the point of
0059   // this would be:
0060   //  rc->set_IntFlag("RANDOMSEED",PHRandomSeed());
0061   // or set it to a fixed value so you can debug your code
0062   //  rc->set_IntFlag("RANDOMSEED", 12345);
0063   Enable::XPLOAD = true;
0064   rc->set_StringFlag("XPLOAD_TAG",XPLOAD::tag);
0065   rc->set_StringFlag("XPLOAD_CONFIG",XPLOAD::config);
0066   rc->set_uint64Flag("TIMESTAMP",XPLOAD::timestamp);
0067
0068   //===============
0069   // Input options
0070   //===============
0071   // verbosity setting (applies to all input managers)
0072   Input::VERBOSITY = 1;
0073   // First enable the input generators
0074   // Either:
0075   // read previously generated g4-hits files, in this case it opens a DST and skips
0076   // the simulations step completely. The G4Setup macro is only loaded to get information
0077   // about the number of layers used for the cell reco code
0078   Input::READHITS = true;
0079   INPUTREADHITS::listfile[0] = inputFile0;
0080   //INPUTREADHITS::filename[0] = inputFile0;
0081   //INPUTREADHITS::filename[1] = inputFile1;
0082
0083   //-----------------
0084   // Initialize the selected Input/Event generation
0085   //-----------------
0086   // This creates the input generator(s)
0087   InputInit();
0088
0089   // register all input generators with Fun4All
0090   InputRegister();
0091
0092 // register the flag handling
0093   FlagHandler *flag = new FlagHandler();
0094   se->registerSubsystem(flag);
0095
0096   // set up production relatedstuff
0097    Enable::PRODUCTION = true;
0098
0099   //======================
0100   // Write the DST
0101   //======================
0102
0103   Enable::DSTOUT = true;
0104   Enable::DSTOUT_COMPRESS = false;
0105   DstOut::OutputDir = outdir;
0106   DstOut::OutputFile = outputFile;
0107
0108   //Option to convert DST to human command readable TTree for quick poke around the outputs
0109   //  Enable::DSTREADER = true;
0110
0111   // turn the display on (default off)
0112   Enable::DISPLAY = false;
0113
0114   //======================
0115   // What to run
0116   //======================
0117   // Global options (enabled for all enables subsystems - if implemented)
0118   //  Enable::ABSORBER = true;
0119   //  Enable::OVERLAPCHECK = true;
0120   //  Enable::VERBOSITY = 1;
0121
0122   Enable::CEMC = true;
0123   Enable::CEMC_CELL = Enable::CEMC && true;
0124   Enable::CEMC_TOWER = Enable::CEMC_CELL && true;
0125   Enable::CEMC_CLUSTER = Enable::CEMC_TOWER && true;
0126   Enable::CEMC_EVAL = Enable::CEMC_CLUSTER && true;
0127
0128   // Enable::HCALIN = true;
0129 //   Enable::HCALIN_CELL = Enable::HCALIN && true;
0130 //   Enable::HCALIN_TOWER = Enable::HCALIN_CELL && true;
0131 //   Enable::HCALIN_CLUSTER = Enable::HCALIN_TOWER && true;
0132 // //  Enable::HCALIN_EVAL = Enable::HCALIN_CLUSTER && true;
0133
0134 //   Enable::HCALOUT = true;
0135 //   Enable::HCALOUT_CELL = Enable::HCALOUT && true;
0136 //   Enable::HCALOUT_TOWER = Enable::HCALOUT_CELL && true;
0137 //   Enable::HCALOUT_CLUSTER = Enable::HCALOUT_TOWER && true;
0138 // //  Enable::HCALOUT_EVAL = Enable::HCALOUT_CLUSTER && true;
0139
0140
0141 //  Enable::EPD = false;
0142
0143   Enable::GLOBAL_RECO = true;
0144   //  Enable::GLOBAL_FASTSIM = true;
0145
0146 //  Enable::CALOTRIGGER = Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER && false;
0147
0148 //  Enable::JETS = true;
0149 //  Enable::JETS_EVAL = Enable::JETS && true;
0150
0151   // HI Jet Reco for p+Au / Au+Au collisions (default is false for
0152   // single particle / p+p-only simulations, or for p+Au / Au+Au
0153   // simulations which don't particularly care about jets)
0154 //  Enable::HIJETS = false && Enable::JETS && Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER;
0155
0156   // 3-D topoCluster reconstruction, potentially in all calorimeter layers
0157 //  Enable::TOPOCLUSTER = true && Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER;
0158   // particle flow jet reconstruction - needs topoClusters!
0159 //  Enable::PARTICLEFLOW = true && Enable::TOPOCLUSTER;
0160
0161   //---------------
0162   // Magnet Settings
0163   //---------------
0164
0165   //  const string magfield = "1.5"; // alternatively to specify a constant magnetic field, give a float number, which will be translated to solenoidal field in T, if string use as fieldmap name (including path)
0166   //  G4MAGNET::magfield = string(getenv("CALIBRATIONROOT")) + string("/Field/Map/sPHENIX.2d.root");  // default map from the calibration database
0167 //  G4MAGNET::magfield_rescale = 1;  // make consistent with expected Babar field strength of 1.4T
0168
0169   //---------------
0170   // Pythia Decayer
0171   //---------------
0172   // list of decay types in
0173   // $OFFLINE_MAIN/include/g4decayer/EDecayType.hh
0174   // default is All:
0175   // G4P6DECAYER::decayType = EDecayType::kAll;
0176
0177   // Initialize the selected subsystems
0178   G4Init();
0179
0180   //---------------------
0181   // GEANT4 Detector description
0182   //---------------------
0183   if (!Input::READHITS)
0184   {
0185     G4Setup();
0186   }
0187
0188   //------------------
0189   // Detector Division
0190   //------------------
0191
0192   if (Enable::BBC || Enable::BBCFAKE) Bbc_Reco();
0193
0194   if (Enable::MVTX_CELL) Mvtx_Cells();
0195   if (Enable::INTT_CELL) Intt_Cells();
0196   if (Enable::TPC_CELL) TPC_Cells();
0197   if (Enable::MICROMEGAS_CELL) Micromegas_Cells();
0198
0199   if (Enable::CEMC_CELL) CEMC_Cells();
0200
0201   if (Enable::HCALIN_CELL) HCALInner_Cells();
0202
0203   if (Enable::HCALOUT_CELL) HCALOuter_Cells();
0204
0205   //-----------------------------
0206   // CEMC towering and clustering
0207   //-----------------------------
0208
0209   if (Enable::CEMC_TOWER) CEMC_Towers();
0210   if (Enable::CEMC_CLUSTER) CEMC_Clusters();
0211
0212   //-----------------------------
0213   // HCAL towering and clustering
0214   //-----------------------------
0215
0216   if (Enable::HCALIN_TOWER) HCALInner_Towers();
0217   if (Enable::HCALIN_CLUSTER) HCALInner_Clusters();
0218
0219   if (Enable::HCALOUT_TOWER) HCALOuter_Towers();
0220   if (Enable::HCALOUT_CLUSTER) HCALOuter_Clusters();
0221
0222   // if enabled, do topoClustering early, upstream of any possible jet reconstruction
0223   if (Enable::TOPOCLUSTER) TopoClusterReco();
0224
0225   if (Enable::DSTOUT_COMPRESS) ShowerCompress();
0226
0227   //--------------
0228   // SVTX tracking
0229   //--------------
0230   if (Enable::MVTX_CLUSTER) Mvtx_Clustering();
0231   if (Enable::INTT_CLUSTER) Intt_Clustering();
0232   if (Enable::TPC_CLUSTER) TPC_Clustering();
0233   if (Enable::MICROMEGAS_CLUSTER) Micromegas_Clustering();
0234
0235   if (Enable::TRACKING_TRACK)
0236   {
0237     TrackingInit();
0238     Tracking_Reco();
0239   }
0240   //-----------------
0241   // Global Vertexing
0242   //-----------------
0243
0244   if (Enable::GLOBAL_RECO && Enable::GLOBAL_FASTSIM)
0245   {
0246     cout << "You can only enable Enable::GLOBAL_RECO or Enable::GLOBAL_FASTSIM, not both" << endl;
0247     gSystem->Exit(1);
0248   }
0249   if (Enable::GLOBAL_RECO)
0250   {
0251     Global_Reco();
0252   }
0253   else if (Enable::GLOBAL_FASTSIM)
0254   {
0255     Global_FastSim();
0256   }
0257
0258   //-----------------
0259   // Calo Trigger Simulation
0260   //-----------------
0261
0262   if (Enable::CALOTRIGGER)
0263   {
0264     CaloTrigger_Sim();
0265   }
0266
0267   //---------
0268   // Jet reco
0269   //---------
0270
0271   if (Enable::JETS) Jet_Reco();
0272   if (Enable::HIJETS) HIJetReco();
0273
0274   if (Enable::PARTICLEFLOW) ParticleFlow();
0275
0276   //----------------------
0277   // Simulation evaluation
0278   //----------------------
0279   string outputroot = outputFile;
0280   string remove_this = ".root";
0281   size_t pos = outputroot.find(remove_this);
0282   if (pos != string::npos)
0283   {
0284     outputroot.erase(pos, remove_this.length());
0285   }
0286
0287   //--------------
0288   // Set up Input Managers
0289   //--------------
0290
0291   InputManagers();
0292
0293   if (Enable::PRODUCTION)
0294   {
0295     Production_CreateOutputDir();
0296   }
0297
0298   if (Enable::DSTOUT)
0299   {
0300     string FullOutFile = DstOut::OutputFile;
0301     Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", FullOutFile);
0302     out->AddNode("Sync");
0303     out->AddNode("EventHeader");
0304     out->AddNode("TOWER_SIM_HCALIN");
0305     out->AddNode("TOWER_RAW_HCALIN");
0306     out->AddNode("TOWER_CALIB_HCALIN");
0307     out->AddNode("CLUSTER_HCALIN");
0308     out->AddNode("TOWER_SIM_HCALOUT");
0309     out->AddNode("TOWER_RAW_HCALOUT");
0310     out->AddNode("TOWER_CALIB_HCALOUT");
0311     out->AddNode("CLUSTER_HCALOUT");
0312     out->AddNode("TOWER_SIM_CEMC");
0313     out->AddNode("TOWER_RAW_CEMC");
0314     out->AddNode("TOWER_CALIB_CEMC");
0315     out->AddNode("CLUSTER_CEMC");
0316     out->AddNode("CLUSTER_POS_COR_CEMC");
0317 // leave the topo cluster here in case we run this during pass3
0318     out->AddNode("TOPOCLUSTER_ALLCALO");
0319     out->AddNode("TOPOCLUSTER_EMCAL");
0320     out->AddNode("TOPOCLUSTER_HCAL");
0321     out->AddNode("GlobalVertexMap");
0322     if (Enable::DSTOUT_COMPRESS) DstCompress(out);
0323     se->registerOutputManager(out);
0324   }
0325   //-----------------
0326   // Event processing
0327   //-----------------
0328   if (Enable::DISPLAY)
0329   {
0330     DisplayOn();
0331
0332     gROOT->ProcessLine("Fun4AllServer *se = Fun4AllServer::instance();");
0333     gROOT->ProcessLine("PHG4Reco *g4 = (PHG4Reco *) se->getSubsysReco(\"PHG4RECO\");");
0334
0335     cout << "-------------------------------------------------" << endl;
0336     cout << "You are in event display mode. Run one event with" << endl;
0337     cout << "se->run(1)" << endl;
0338     cout << "Run Geant4 command with following examples" << endl;
0339     gROOT->ProcessLine("displaycmd()");
0340
0341     return 0;
0342   }
0343
0344   // if we use a negative number of events we go back to the command line here
0345   if (nEvents < 0)
0346   {
0347     return 0;
0348   }
0349   // if we run the particle generator and use 0 it'll run forever
0350   if (nEvents == 0 && !Input::HEPMC && !Input::READHITS)
0351   {
0352     cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
0353     cout << "it will run forever, so I just return without running anything" << endl;
0354     return 0;
0355   }
0356   pi0Efficiency *eval = new pi0Efficiency("dummy", outputroot + "_pi0Efficiency.root");
0357   se->registerSubsystem(eval);
0358
0359   se->run(nEvents);
0360
0361   //-----
0362   // Exit
0363   //-----
0364
0365   //XploadInterface::instance()->Print(); // print used DB files
0366   se->End();
0367   //se->PrintTimer();
0368   std::cout << "All done" << std::endl;
0369   delete se;
0370   // if (Enable::PRODUCTION)
0371   // {
0372   //   Production_MoveOutput();
0373   // }
0374
0375   gSystem->Exit(0);
0376   return 0;
0377 }
0378 #endif