tutorials/sPHENIX_sims/Fun4All_G4_sPHENIX.C

0001
0002 //modified and simplified Fun4All macro (Simple inputgen + HCALOUT)
0003
0004 #ifndef MACRO_FUN4ALLG4SPHENIX_C
0005 #define MACRO_FUN4ALLG4SPHENIX_C
0006
0007 #include "Calib.h"
0008 #include <GlobalVariables.C>
0009
0010 #include <DisplayOn.C>
0011 #include <G4Setup_sPHENIX.C>
0012 #include <G4_Bbc.C>
0013 #include <G4_CaloTrigger.C>
0014 #include <G4_Centrality.C>
0015 #include <G4_DSTReader.C>
0016 #include <G4_Global.C>
0017 #include <G4_HIJetReco.C>
0018 #include <G4_Input.C>
0019 #include <G4_Jets.C>
0020 #include <G4_KFParticle.C>
0021 #include <G4_ParticleFlow.C>
0022 #include <G4_Production.C>
0023 #include <G4_TopoClusterReco.C>
0024 #include <G4_Tracking.C>
0025 #include <G4_User.C>
0026 #include <QA.C>
0027
0028 #include <ffamodules/FlagHandler.h>
0029 #include <ffamodules/HeadReco.h>
0030 #include <ffamodules/SyncReco.h>
0031 #include <ffamodules/XploadInterface.h>
0032
0033 #include <fun4all/Fun4AllDstOutputManager.h>
0034 #include <fun4all/Fun4AllOutputManager.h>
0035 #include <fun4all/Fun4AllServer.h>
0036
0037 #include <phool/PHRandomSeed.h>
0038 #include <phool/recoConsts.h>
0039
0040 R__LOAD_LIBRARY(libfun4all.so)
0041 R__LOAD_LIBRARY(libffamodules.so)
0042 R__LOAD_LIBRARY(libCalib.so.0.0.0)
0043
0044
0045 // For HepMC Hijing
0046 // try inputFile = /sphenix/sim/sim01/sphnxpro/sHijing_HepMC/sHijing_0-12fm.dat
0047
0048 int Fun4All_G4_sPHENIX(
0049     const int nEvents = 1,
0050     const string &inputFile = "https://www.phenix.bnl.gov/WWW/publish/phnxbld/sPHENIX/files/sPHENIX_G4Hits_sHijing_9-11fm_00000_00010.root",
0051     const string &outputFile = "G4sPHENIX.root",
0052     const string &embed_input_file = "https://www.phenix.bnl.gov/WWW/publish/phnxbld/sPHENIX/files/sPHENIX_G4Hits_sHijing_9-11fm_00000_00010.root",
0053     const int skip = 0,
0054     const string &outdir = ".")
0055 {
0056   Fun4AllServer *se = Fun4AllServer::instance();
0057   se->Verbosity(0);
0058
0059   //Opt to print all random seed used for debugging reproducibility. Comment out to reduce stdout prints.
0060   PHRandomSeed::Verbosity(1);
0061
0062   // just if we set some flags somewhere in this macro
0063   recoConsts *rc = recoConsts::instance();
0064
0065   //===============
0066   // Input options
0067   //===============
0068   // verbosity setting (applies to all input managers)
0069   Input::VERBOSITY = 0;
0070
0071   Input::SIMPLE = true;
0072
0073   //  Input::GUN = true;
0074
0075   //-----------------
0076   // Initialize the selected Input/Event generation
0077   //-----------------
0078   // This creates the input generator(s)
0079   InputInit();
0080
0081   //--------------
0082   // Set generator specific options
0083   //--------------
0084   // can only be set after InputInit() is called
0085
0086   // Simple Input generator:
0087   // if you run more than one of these Input::SIMPLE_NUMBER > 1
0088   // add the settings for other with [1], next with [2]...
0089
0090   if (Input::SIMPLE)
0091   {
0092     INPUTGENERATOR::SimpleEventGenerator[0]->add_particles("pi-", 1);
0093
0094     INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_function(PHG4SimpleEventGenerator::Uniform,
0095                                                                                 PHG4SimpleEventGenerator::Uniform,
0096                                                                                 PHG4SimpleEventGenerator::Uniform);                                      INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_mean(0., 0., 0.);
0097     INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_mean(140., 140., 10.);
0098     INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_width(0., 0., 0.);
0099     INPUTGENERATOR::SimpleEventGenerator[0]->set_eta_range(0, 0);
0100     INPUTGENERATOR::SimpleEventGenerator[0]->set_phi_range(0, 0);
0101     INPUTGENERATOR::SimpleEventGenerator[0]->set_pt_range(10., 10.);
0102   }
0103
0104   // particle gun
0105   // if you run more than one of these Input::GUN_NUMBER > 1
0106   // add the settings for other with [1], next with [2]...
0107   if (Input::GUN)
0108   {
0109     INPUTGENERATOR::Gun[0]->AddParticle("pi-", 0, 1, 0);
0110     INPUTGENERATOR::Gun[0]->set_vtx(0, 0, 0);
0111   }
0112
0113   //--------------
0114   // Set Input Manager specific options
0115   //--------------
0116   // can only be set after InputInit() is called
0117
0118  // register all input generators with Fun4All
0119   InputRegister();
0120
0121   if (! Input::READHITS)
0122   {
0123     rc->set_IntFlag("RUNNUMBER",1);
0124
0125     SyncReco *sync = new SyncReco();
0126     se->registerSubsystem(sync);
0127
0128     HeadReco *head = new HeadReco();
0129     se->registerSubsystem(head);
0130   }
0131
0132 // Flag Handler is always needed to read flags from input (if used)
0133 // and update our rc flags with them. At the end it saves all flags
0134 // again on the DST in the Flags node under the RUN node
0135   FlagHandler *flag = new FlagHandler();
0136   se->registerSubsystem(flag);
0137
0138   // set up production relatedstuff
0139   //   Enable::PRODUCTION = true;
0140
0141   //======================
0142   // Write the DST
0143   //======================
0144
0145   //Enable::DSTOUT = true;
0146   Enable::DSTOUT_COMPRESS = false;
0147   DstOut::OutputDir = outdir;
0148   DstOut::OutputFile = outputFile;
0149
0150   //Option to convert DST to human command readable TTree for quick poke around the outputs
0151   //  Enable::DSTREADER = true;
0152
0153   // turn the display on (default off)
0154    Enable::DISPLAY = false;
0155
0156   //======================
0157   // What to run
0158   //======================
0159
0160   // QA, main switch
0161   Enable::QA = false;
0162
0163   // Global options (enabled for all enables subsystems - if implemented)
0164   //  Enable::ABSORBER = true;
0165   //  Enable::OVERLAPCHECK = true;
0166   //  Enable::VERBOSITY = 1;
0167
0168   // Enable::BBC = true;
0169   // Enable::BBC_SUPPORT = true; // save hist in bbc support structure
0170   Enable::BBCFAKE = true;  // Smeared vtx and t0, use if you don't want real BBC in simulation
0171
0172   Enable::PIPE = true;
0173   Enable::PIPE_ABSORBER = true;
0174
0175   // central tracking
0176   Enable::MVTX = false;
0177   Enable::MVTX_CELL = Enable::MVTX && true;
0178   Enable::MVTX_CLUSTER = Enable::MVTX_CELL && true;
0179   Enable::MVTX_QA = Enable::MVTX_CLUSTER && Enable::QA && true;
0180
0181   Enable::INTT = false;
0182 //  Enable::INTT_ABSORBER = true; // enables layerwise support structure readout
0183 //  Enable::INTT_SUPPORT = true; // enable global support structure readout
0184   Enable::INTT_CELL = Enable::INTT && true;
0185   Enable::INTT_CLUSTER = Enable::INTT_CELL && true;
0186   Enable::INTT_QA = Enable::INTT_CLUSTER && Enable::QA && true;
0187
0188   Enable::TPC = false;
0189   Enable::TPC_ABSORBER = true;
0190   Enable::TPC_CELL = Enable::TPC && true;
0191   Enable::TPC_CLUSTER = Enable::TPC_CELL && true;
0192   Enable::TPC_QA = Enable::TPC_CLUSTER && Enable::QA && true;
0193
0194   Enable::MICROMEGAS = false;
0195   Enable::MICROMEGAS_CELL = Enable::MICROMEGAS && true;
0196   Enable::MICROMEGAS_CLUSTER = Enable::MICROMEGAS_CELL && true;
0197   Enable::MICROMEGAS_QA = Enable::MICROMEGAS_CLUSTER && Enable::QA && true;
0198
0199   Enable::TRACKING_TRACK = false;
0200   //Enable::TRACKING_EVAL = Enable::TRACKING_TRACK && true;
0201   Enable::TRACKING_QA = Enable::TRACKING_TRACK && Enable::QA && true;
0202
0203   //  cemc electronics + thin layer of W-epoxy to get albedo from cemc
0204   //  into the tracking, cannot run together with CEMC
0205   //  Enable::CEMCALBEDO = true;
0206
0207   Enable::CEMC = false;
0208   Enable::CEMC_ABSORBER = true;
0209   Enable::CEMC_CELL = Enable::CEMC && true;
0210   Enable::CEMC_TOWER = Enable::CEMC_CELL && true;
0211   Enable::CEMC_CLUSTER = Enable::CEMC_TOWER && true;
0212   //Enable::CEMC_EVAL = Enable::CEMC_CLUSTER && true;
0213   Enable::CEMC_QA = Enable::CEMC_CLUSTER && Enable::QA && true;
0214
0215   Enable::HCALIN = false;
0216   Enable::HCALIN_ABSORBER = true;
0217   Enable::HCALIN_CELL = Enable::HCALIN && true;
0218   Enable::HCALIN_TOWER = Enable::HCALIN_CELL && true;
0219   Enable::HCALIN_CLUSTER = Enable::HCALIN_TOWER && true;
0220   //Enable::HCALIN_EVAL = Enable::HCALIN_CLUSTER && true;
0221   Enable::HCALIN_QA = Enable::HCALIN_CLUSTER && Enable::QA && true;
0222
0223   Enable::MAGNET = false;
0224   Enable::MAGNET_ABSORBER = true;
0225
0226   Enable::HCALOUT = true;
0227   Enable::HCALOUT_ABSORBER = true;
0228   Enable::HCALOUT_CELL = Enable::HCALOUT && true;
0229   Enable::HCALOUT_TOWER = Enable::HCALOUT_CELL && true;
0230   Enable::HCALOUT_CLUSTER = Enable::HCALOUT_TOWER && true;
0231   //Enable::HCALOUT_EVAL = Enable::HCALOUT_CLUSTER && true;
0232   Enable::HCALOUT_QA = Enable::HCALOUT_CLUSTER && Enable::QA && true;
0233
0234   Enable::EPD = false;
0235
0236   Enable::BEAMLINE = true;
0237 //  Enable::BEAMLINE_ABSORBER = true;  // makes the beam line magnets sensitive volumes
0238 //  Enable::BEAMLINE_BLACKHOLE = true; // turns the beamline magnets into black holes
0239   Enable::ZDC = false;
0240 //  Enable::ZDC_ABSORBER = true;
0241 //  Enable::ZDC_SUPPORT = true;
0242   Enable::ZDC_TOWER = Enable::ZDC && true;
0243   //Enable::ZDC_EVAL = Enable::ZDC_TOWER && true;
0244
0245   //! forward flux return plug door. Out of acceptance and off by default.
0246   //Enable::PLUGDOOR = true;
0247   Enable::PLUGDOOR_ABSORBER = true;
0248
0249   //Enable::GLOBAL_RECO = (Enable::BBCFAKE || Enable::TRACKING_TRACK) && true;
0250   Enable::GLOBAL_FASTSIM = true;
0251
0252   //Enable::KFPARTICLE = true;
0253   //Enable::KFPARTICLE_VERBOSITY = 1;
0254   //Enable::KFPARTICLE_TRUTH_MATCH = true;
0255   //Enable::KFPARTICLE_SAVE_NTUPLE = true;
0256
0257   Enable::CALOTRIGGER = Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER && false;
0258
0259   //Enable::JETS = (Enable::GLOBAL_RECO || Enable::GLOBAL_FASTSIM) && true;
0260   //Enable::JETS_EVAL = Enable::JETS && true;
0261   Enable::JETS_QA = Enable::JETS && Enable::QA && true;
0262
0263   // HI Jet Reco for p+Au / Au+Au collisions (default is false for
0264   // single particle / p+p-only simulations, or for p+Au / Au+Au
0265   // simulations which don't particularly care about jets)
0266   //Enable::HIJETS = Enable::JETS && Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER && false;
0267
0268   // 3-D topoCluster reconstruction, potentially in all calorimeter layers
0269   Enable::TOPOCLUSTER = Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER && false;
0270   // particle flow jet reconstruction - needs topoClusters!
0271   Enable::PARTICLEFLOW = Enable::TOPOCLUSTER && true;
0272   // centrality reconstruction
0273   Enable::CENTRALITY = true;
0274
0275   // new settings using Enable namespace in GlobalVariables.C
0276   Enable::BLACKHOLE = true;
0277   //Enable::BLACKHOLE_SAVEHITS = false; // turn off saving of bh hits
0278   //Enable::BLACKHOLE_FORWARD_SAVEHITS = false; // disable forward/backward hits
0279   //BlackHoleGeometry::visible = true;
0280
0281   // run user provided code (from local G4_User.C)
0282   //Enable::USER = true;
0283
0284   //===============
0285   // conditions DB flags
0286   //===============
0287   //Enable::XPLOAD = true;
0288   // tag
0289   rc->set_StringFlag("XPLOAD_TAG",XPLOAD::tag);
0290   // database config
0291   rc->set_StringFlag("XPLOAD_CONFIG",XPLOAD::config);
0292   // 64 bit timestamp
0293   rc->set_uint64Flag("TIMESTAMP",XPLOAD::timestamp);
0294
0295   //---------------
0296   // World Settings
0297   //---------------
0298   //  G4WORLD::PhysicsList = "FTFP_BERT"; //FTFP_BERT_HP best for calo
0299   //  G4WORLD::WorldMaterial = "G4_AIR"; // set to G4_GALACTIC for material scans
0300
0301   //---------------
0302   // Magnet Settings
0303   //---------------
0304
0305   //  G4MAGNET::magfield =  string(getenv("CALIBRATIONROOT"))+ string("/Field/Map/sphenix3dbigmapxyz.root");  // default map from the calibration database
0306   //  G4MAGNET::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)
0307 //  G4MAGNET::magfield_rescale = 1.;  // make consistent with expected Babar field strength of 1.4T
0308
0309   //---------------
0310   // Pythia Decayer
0311   //---------------
0312   // list of decay types in
0313   // $OFFLINE_MAIN/include/g4decayer/EDecayType.hh
0314   // default is All:
0315   // G4P6DECAYER::decayType = EDecayType::kAll;
0316
0317   // Initialize the selected subsystems
0318   G4Init();
0319
0320   //---------------------
0321   // GEANT4 Detector description
0322   //---------------------
0323   if (!Input::READHITS)
0324   {
0325     G4Setup();
0326   }
0327
0328   //------------------
0329   // Detector Division
0330   //------------------
0331
0332   if (Enable::BBC || Enable::BBCFAKE) Bbc_Reco();
0333
0334   if (Enable::MVTX_CELL) Mvtx_Cells();
0335   if (Enable::INTT_CELL) Intt_Cells();
0336   if (Enable::TPC_CELL) TPC_Cells();
0337   if (Enable::MICROMEGAS_CELL) Micromegas_Cells();
0338
0339   if (Enable::CEMC_CELL) CEMC_Cells();
0340
0341   if (Enable::HCALIN_CELL) HCALInner_Cells();
0342
0343   if (Enable::HCALOUT_CELL) HCALOuter_Cells();
0344
0345   //-----------------------------
0346   // CEMC towering and clustering
0347   //-----------------------------
0348
0349   if (Enable::CEMC_TOWER) CEMC_Towers();
0350   if (Enable::CEMC_CLUSTER) CEMC_Clusters();
0351
0352   //-----------------------------
0353   // HCAL towering and clustering
0354   //-----------------------------
0355
0356   if (Enable::HCALIN_TOWER) HCALInner_Towers();
0357   if (Enable::HCALIN_CLUSTER) HCALInner_Clusters();
0358
0359   if (Enable::HCALOUT_TOWER) HCALOuter_Towers();
0360   if (Enable::HCALOUT_CLUSTER) HCALOuter_Clusters();
0361
0362   // if enabled, do topoClustering early, upstream of any possible jet reconstruction
0363   if (Enable::TOPOCLUSTER) TopoClusterReco();
0364
0365   //--------------
0366   // SVTX tracking
0367   //--------------
0368   if(Enable::TRACKING_TRACK)
0369     {
0370       TrackingInit();
0371     }
0372   if (Enable::MVTX_CLUSTER) Mvtx_Clustering();
0373   if (Enable::INTT_CLUSTER) Intt_Clustering();
0374   if (Enable::TPC_CLUSTER) TPC_Clustering();
0375   if (Enable::MICROMEGAS_CLUSTER) Micromegas_Clustering();
0376
0377   if (Enable::TRACKING_TRACK)
0378   {
0379     Tracking_Reco();
0380   }
0381   //-----------------
0382   // Global Vertexing
0383   //-----------------
0384
0385   if (Enable::GLOBAL_RECO && Enable::GLOBAL_FASTSIM)
0386   {
0387     cout << "You can only enable Enable::GLOBAL_RECO or Enable::GLOBAL_FASTSIM, not both" << endl;
0388     gSystem->Exit(1);
0389   }
0390   if (Enable::GLOBAL_RECO)
0391   {
0392     Global_Reco();
0393   }
0394   else if (Enable::GLOBAL_FASTSIM)
0395   {
0396     Global_FastSim();
0397   }
0398
0399   //-----------------
0400   // Centrality Determination
0401   //-----------------
0402
0403   if (Enable::CENTRALITY)
0404   {
0405       Centrality();
0406   }
0407
0408   //-----------------
0409   // Calo Trigger Simulation
0410   //-----------------
0411
0412   if (Enable::CALOTRIGGER)
0413   {
0414     CaloTrigger_Sim();
0415   }
0416
0417   //---------
0418   // Jet reco
0419   //---------
0420
0421   if (Enable::JETS) Jet_Reco();
0422   if (Enable::HIJETS) HIJetReco();
0423
0424   if (Enable::PARTICLEFLOW) ParticleFlow();
0425
0426   //----------------------
0427   // Simulation evaluation
0428   //----------------------
0429   string outputroot = outputFile;
0430   string remove_this = ".root";
0431   size_t pos = outputroot.find(remove_this);
0432   if (pos != string::npos)
0433   {
0434     outputroot.erase(pos, remove_this.length());
0435   }
0436
0437   if (Enable::TRACKING_EVAL) Tracking_Eval(outputroot + "_g4svtx_eval.root");
0438
0439   if (Enable::CEMC_EVAL) CEMC_Eval(outputroot + "_g4cemc_eval.root");
0440
0441   if (Enable::HCALIN_EVAL) HCALInner_Eval(outputroot + "_g4hcalin_eval.root");
0442
0443   if (Enable::HCALOUT_EVAL) HCALOuter_Eval(outputroot + "_g4hcalout_eval.root");
0444
0445   if (Enable::JETS_EVAL) Jet_Eval(outputroot + "_g4jet_eval.root");
0446
0447   if (Enable::DSTREADER) G4DSTreader(outputroot + "_DSTReader.root");
0448
0449   if (Enable::USER) UserAnalysisInit();
0450
0451
0452     gSystem->Load("libCalib.so.0.0.0");
0453     Calib *hcal_eval = new Calib();
0454     hcal_eval->set_filename("hcal_tower_energies.root");
0455     se->registerSubsystem(hcal_eval);
0456
0457
0458
0459
0460   //======================
0461   // Run KFParticle on evt
0462   //======================
0463   if (Enable::KFPARTICLE && Input::UPSILON) KFParticle_Upsilon_Reco();
0464   if (Enable::KFPARTICLE && Input::DZERO) KFParticle_D0_Reco();
0465
0466   //----------------------
0467   // Standard QAs
0468   //----------------------
0469
0470   if (Enable::CEMC_QA) CEMC_QA();
0471   if (Enable::HCALIN_QA) HCALInner_QA();
0472   if (Enable::HCALOUT_QA) HCALOuter_QA();
0473
0474   if (Enable::JETS_QA) Jet_QA();
0475
0476   if (Enable::MVTX_QA) Mvtx_QA();
0477   if (Enable::INTT_QA) Intt_QA();
0478   if (Enable::TPC_QA) TPC_QA();
0479   if (Enable::MICROMEGAS_QA) Micromegas_QA();
0480   if (Enable::TRACKING_QA) Tracking_QA();
0481
0482   if (Enable::TRACKING_QA && Enable::CEMC_QA && Enable::HCALIN_QA && Enable::HCALOUT_QA) QA_G4CaloTracking();
0483
0484   //--------------
0485   // Set up Input Managers
0486   //--------------
0487
0488   InputManagers();
0489
0490   if (Enable::PRODUCTION)
0491   {
0492     Production_CreateOutputDir();
0493   }
0494
0495   if (Enable::DSTOUT)
0496   {
0497     string FullOutFile = DstOut::OutputDir + "/" + DstOut::OutputFile;
0498     Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", FullOutFile);
0499     if (Enable::DSTOUT_COMPRESS)
0500     {
0501       ShowerCompress();
0502       DstCompress(out);
0503     }
0504     se->registerOutputManager(out);
0505   }
0506   //-----------------
0507   // Event processing
0508   //-----------------
0509   if (Enable::DISPLAY)
0510   {
0511     DisplayOn();
0512
0513     gROOT->ProcessLine("Fun4AllServer *se = Fun4AllServer::instance();");
0514     gROOT->ProcessLine("PHG4Reco *g4 = (PHG4Reco *) se->getSubsysReco(\"PHG4RECO\");");
0515
0516     cout << "-------------------------------------------------" << endl;
0517     cout << "You are in event display mode. Run one event with" << endl;
0518     cout << "se->run(1)" << endl;
0519     cout << "Run Geant4 command with following examples" << endl;
0520     gROOT->ProcessLine("displaycmd()");
0521
0522     return 0;
0523   }
0524
0525   // if we use a negative number of events we go back to the command line here
0526   if (nEvents < 0)
0527   {
0528     return 0;
0529   }
0530   // if we run the particle generator and use 0 it'll run forever
0531   // for embedding it runs forever if the repeat flag is set
0532   if (nEvents == 0 && !Input::HEPMC && !Input::READHITS && INPUTEMBED::REPEAT)
0533   {
0534     cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
0535     cout << "it will run forever, so I just return without running anything" << endl;
0536     return 0;
0537   }
0538
0539   se->skip(skip);
0540   se->run(nEvents);
0541
0542   //-----
0543   // QA output
0544   //-----
0545
0546   if (Enable::QA) QA_Output(outputroot + "_qa.root");
0547
0548   //-----
0549   // Exit
0550   //-----
0551
0552   se->End();
0553   std::cout << "All done" << std::endl;
0554   delete se;
0555   if (Enable::PRODUCTION)
0556   {
0557     Production_MoveOutput();
0558   }
0559
0560   gSystem->Exit(0);
0561   return 0;
0562 }
0563 #endif