calo/emc_pi0_tbt/Fun4All_pi0tbt_SIMPLE_EMBED.C

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