ChargedHadronMultiplicity-SmallSystem/simulation/Fun4All_G4_sPHENIX.C

0001 #ifndef MACRO_FUN4ALLG4SPHENIX_C
0002 #define MACRO_FUN4ALLG4SPHENIX_C
0003
0004 #include <GlobalVariables.C>
0005
0006 #include "G4Setup_sPHENIX.C"
0007 #include "RecoScanConfig.h"
0008
0009 #include <DisplayOn.C>
0010 #include <G4_Mbd.C>
0011 #include <G4_CaloTrigger.C>
0012 #include <G4_Centrality.C>
0013 #include <G4_DSTReader.C>
0014 #include <G4_Global.C>
0015 #include <G4_HIJetReco.C>
0016 #include <G4_Input.C>
0017 #include <G4_Jets.C>
0018 #include <G4_KFParticle.C>
0019 #include <G4_ParticleFlow.C>
0020 #include <G4_Production.C>
0021 #include <G4_TopoClusterReco.C>
0022
0023 #include <Trkr_RecoInit.C>
0024 #include <Trkr_Clustering.C>
0025 #include <Trkr_LaserClustering.C>
0026 #include "./Trkr_Reco.C" // for testing purpose
0027 #include <Trkr_Eval.C>
0028 #include <Trkr_QA.C>
0029
0030 #include <Trkr_Diagnostics.C>
0031 #include <G4_User.C>
0032 #include <QA.C>
0033
0034 #include <ffamodules/FlagHandler.h>
0035 #include <ffamodules/HeadReco.h>
0036 #include <ffamodules/SyncReco.h>
0037 #include <ffamodules/CDBInterface.h>
0038
0039 #include <fun4all/Fun4AllDstOutputManager.h>
0040 #include <fun4all/Fun4AllOutputManager.h>
0041 #include <fun4all/Fun4AllServer.h>
0042
0043 #include <calotrigger/MinimumBiasClassifier.h>
0044 #include <centrality/CentralityReco.h>
0045 #include <g4centrality/PHG4CentralityReco.h>
0046 #include <trackreco/PHTruthVertexing.h>
0047
0048 #include <caloreco/CaloTowerBuilder.h>
0049 #include <caloreco/CaloTowerCalib.h>
0050 #include <globalvertex/GlobalVertexReco.h>
0051 #include <mbd/MbdReco.h>
0052 #include <zdcinfo/ZdcReco.h>
0053
0054 #include <phool/PHRandomSeed.h>
0055 #include <phool/recoConsts.h>
0056
0057 #include <Rtypes.h> // resolves R__LOAD_LIBRARY for clang-tidy
0058 #include <TROOT.h>
0059
0060 #include <trackreco/PHCASiliconSeeding.h> // CA silicon seeding
0061
0062 #include <g4eval/SvtxTruthRecoTableEval.h>
0063
0064 #include <vertexcompare/VertexCompare.h>
0065 #include <limits>
0066 #include <sstream>
0067
0068 R__LOAD_LIBRARY(libfun4all.so)
0069 R__LOAD_LIBRARY(libffamodules.so)
0070
0071 R__LOAD_LIBRARY(libmbd.so)
0072 R__LOAD_LIBRARY(libglobalvertex.so)
0073 R__LOAD_LIBRARY(libVertexCompare.so)
0074 // R__LOAD_LIBRARY(libfun4allraw.so)
0075 // R__LOAD_LIBRARY(libffarawmodules.so)
0076 R__LOAD_LIBRARY(libcentrality_io.so)
0077 R__LOAD_LIBRARY(libcentrality.so)
0078 R__LOAD_LIBRARY(libg4centrality.so)
0079 R__LOAD_LIBRARY(libcalotrigger.so)
0080 // R__LOAD_LIBRARY(libcentralityvalid.so)
0081 R__LOAD_LIBRARY(libzdcinfo.so)
0082
0083 R__LOAD_LIBRARY(libg4eval.so)
0084
0085 // For HepMC Hijing
0086 // try inputFile = /sphenix/sim/sim01/sphnxpro/sHijing_HepMC/sHijing_0-12fm.dat
0087
0088 int Fun4All_G4_sPHENIX(const int nEvents = 1,                                             //
0089                        const string generator = "HIJING",                                 //
0090                        const std::string &seedingalgo = "ACTS",                           //
0091                        const std::string &outputFile = "test.root",                       //
0092                        const int process = 0,                                             //
0093                        const std::string &scanName = "",                                  //
0094                        const double scanValue = std::numeric_limits<double>::quiet_NaN(), //
0095                        const std::string &scanTag = ""                                    //
0096 )
0097 {
0098     bool UseActsSiliconSeeding = true;
0099     if (seedingalgo == "CA")
0100     {
0101         UseActsSiliconSeeding = false;
0102     }
0103     else if (seedingalgo == "ACTS")
0104     {
0105         UseActsSiliconSeeding = true;
0106     }
0107     else
0108     {
0109         std::cout << "Seeding algorithm " << seedingalgo << " not currently supported, exiting now and please check input." << std::endl;
0110         return -1;
0111     }
0112     const std::string productionTag = "OO_0_15fm";
0113     const std::string species = "OO";
0114     const int runNumber = 1;
0115     std::string compareOutputFile = outputFile;
0116
0117     gRecoScanConfig = RecoScanConfig{};
0118     if (!scanName.empty() && !applyScanOverride(gRecoScanConfig, scanName, scanValue))
0119     {
0120         std::cout << "Scan variable " << scanName << " not currently supported or has an invalid value, exiting now and please check input." << std::endl;
0121         return -1;
0122     }
0123
0124     std::string effectiveTag = scanTag;
0125     if (effectiveTag.empty() && !scanName.empty())
0126     {
0127         effectiveTag = RecoScan::buildTagToken(scanName, scanValue);
0128     }
0129
0130     std::string defaultSuffix = buildRecoScanDefaultSuffix(scanName, gRecoScanConfig);
0131
0132     if (!defaultSuffix.empty())
0133     {
0134         if (!effectiveTag.empty())
0135         {
0136             effectiveTag += "_" + defaultSuffix;
0137         }
0138         else
0139         {
0140             effectiveTag = defaultSuffix;
0141         }
0142     }
0143
0144     if (!effectiveTag.empty())
0145     {
0146         const auto slashPos = compareOutputFile.find_last_of("/");
0147         if (slashPos == std::string::npos)
0148         {
0149             compareOutputFile = effectiveTag + "/" + compareOutputFile;
0150         }
0151         else
0152         {
0153             compareOutputFile = compareOutputFile.substr(0, slashPos + 1) + effectiveTag + compareOutputFile.substr(slashPos);
0154         }
0155     }
0156
0157     Fun4AllServer *se = Fun4AllServer::instance();
0158     se->Verbosity(0);
0159
0160     //===============
0161     // Input options
0162     //===============
0163     // verbosity setting (applies to all input managers)
0164     Input::VERBOSITY = 0;
0165     // First enable the input generators
0166     // Either:
0167     // read previously generated g4-hits files, in this case it opens a DST and skips
0168     // the simulations step completely. The G4Setup macro is only loaded to get information
0169     // about the number of layers used for the cell reco code
0170     Input::READHITS = true;
0171     bool fromG4Hits = false;
0172     if (generator == "HIJING")
0173     {
0174         // G4 hits files (Pass1)
0175         if (fromG4Hits)
0176         {
0177             // Run35, OO 0-15fm, no pileup
0178             INPUTREADHITS::filename[0] = "/sphenix/lustre01/sphnxpro/mdc2/shijing_hepmc/OO_0_15fm/g4hits/run0035/G4Hits_sHijing_OO_0_15fm-0000000035-" + std::string(Form("%06d", process)) + ".root";
0179             // INPUTREADHITS::filename[0] = "/sphenix/user/adeebsaed/macros/detectors/sPHENIX/outputfiles_fun4all/G4sPHENIX-0000000001-" + std::string(Form("%06d", process)) + ".root";
0180         }
0181         else
0182         {
0183             // TRKR G4HIT (Pass2)
0184             // Run34, OO 0-15fm, pileup, 220kHz
0185             INPUTREADHITS::filename[0] = "/sphenix/lustre01/sphnxpro/mdc2/shijing_hepmc/OO_0_15fm/pileup/run0034/220kHz/DST_TRKR_G4HIT_sHijing_OO_0_15fm_220kHz_bkg_0_15fm-0000000034-" + std::string(Form("%06d", process)) + ".root";
0186
0187             // TRUTH_G4HIT (Pass2)
0188             // Run34, OO 0-15fm, pileup, 220kHz
0189             INPUTREADHITS::filename[1] = "/sphenix/lustre01/sphnxpro/mdc2/shijing_hepmc/OO_0_15fm/pileup/run0034/220kHz/DST_TRUTH_G4HIT_sHijing_OO_0_15fm_220kHz_bkg_0_15fm-0000000034-" + std::string(Form("%06d", process)) + ".root";
0190
0191             // MBD/sEPD (Pass3, MBD/EPD)
0192             // Run34, OO 0-15fm, pileup, 220kHz
0193             INPUTREADHITS::filename[2] = "/sphenix/lustre01/sphnxpro/mdc2/shijing_hepmc/OO_0_15fm/mbdepd/run0034/220khz/DST_MBD_EPD_sHijing_OO_0_15fm_220kHz_bkg_0_15fm-0000000034-" + std::string(Form("%06d", process)) + ".root";
0194
0195             // TRUTH (Pass3 Track)
0196             // Run34, OO 0-15fm, pileup, 220kHz
0197             INPUTREADHITS::filename[3] = "/sphenix/lustre01/sphnxpro/mdc2/shijing_hepmc/OO_0_15fm/trkrhit/run0034/220khz/DST_TRUTH_sHijing_OO_0_15fm_220kHz_bkg_0_15fm-0000000034-" + std::string(Form("%06d", process)) + ".root";
0198
0199             // TrkrHits (Pass3 Track)
0200             // Run34, OO 0-15fm, pileup, 220kHz
0201             INPUTREADHITS::filename[4] = "/sphenix/lustre01/sphnxpro/mdc2/shijing_hepmc/OO_0_15fm/trkrhit/run0034/220khz/DST_TRKR_HIT_sHijing_OO_0_15fm_220kHz_bkg_0_15fm-0000000034-" + std::string(Form("%06d", process)) + ".root";
0202
0203             // TrkrClusters (Pass4 Job0)
0204             // Run34, OO 0-15fm, pileup, 220kHz
0205             INPUTREADHITS::filename[5] = "/sphenix/lustre01/sphnxpro/mdc2/shijing_hepmc/OO_0_15fm/trkrcluster/run0034/220khz/DST_TRKR_CLUSTER_sHijing_OO_0_15fm_220kHz_bkg_0_15fm-0000000034-" + std::string(Form("%06d", process)) + ".root";
0206
0207             // TRACKSEEDS (Pass4 JobA)
0208             // Run34, OO 0-15fm, pileup, 220kHz
0209             // INPUTREADHITS::filename[6] = "/sphenix/lustre01/sphnxpro/mdc2/shijing_hepmc/OO_0_15fm/nopileup/trackseeds/run0035/DST_TRACKSEEDS_sHijing_OO_0_15fm-0000000035-" + std::string(Form("%06d", process)) + ".root";
0210
0211             // GLOBAL (Pass5 Global) --> disable, because we want to build vertices from silicon instead of using the ones in DSTs which are built from the full tracks
0212             // Run34, OO 0-15fm, pileup, 220kHz
0213             // INPUTREADHITS::filename[7] = "/sphenix/lustre01/sphnxpro/mdc2/shijing_hepmc/OO_0_15fm/global/run0034/220khz/DST_GLOBAL_sHijing_OO_0_15fm_220kHz_bkg_0_15fm-0000000034-" + std::string(Form("%06d", process)) + ".root";
0214
0215             // TRUTH RECO (Pass5 Truth Track) --> disable, build the reco-truth matching map/container from SvtxTrackEval
0216             // Run34, OO 0-15fm, pileup, 220kHz
0217             // INPUTREADHITS::filename[8] = "/sphenix/lustre01/sphnxpro/mdc2/shijing_hepmc/OO_0_15fm/truthreco/run0034/220khz/DST_TRUTH_RECO_sHijing_OO_0_15fm_220kHz_bkg_0_15fm-0000000034-" + std::string(Form("%06d", process)) + ".root";
0218         }
0219     }
0220     else
0221     {
0222         std::cout << "Generator " << generator << " not currently supported, exiting now and please check input." << std::endl;
0223         return -1;
0224     }
0225
0226     //======================
0227     // What to run
0228     //======================
0229
0230     // QA, main switch
0231     // Enable::QA = true;
0232     // Global options (enabled for all enables subsystems - if implemented)
0233     //  Enable::ABSORBER = true;
0234     //  Enable::OVERLAPCHECK = true;
0235     //  Enable::VERBOSITY = 1;
0236
0237     Enable::MBD = true;
0238     // Enable::MBD_SUPPORT = true; // save hist in MBD/BBC support structure
0239     Enable::MBDRECO = Enable::MBD && true;
0240     // Enable::MBDFAKE = true; // Smeared vtx and t0, use if you don't want real MBD/BBC in simulation
0241
0242     Enable::PIPE = true;
0243     Enable::PIPE_ABSORBER = true;
0244
0245     // central tracking
0246     Enable::MVTX = true;
0247     Enable::INTT = true;
0248     Enable::TPC = true;
0249     Enable::MICROMEGAS = true;
0250     Enable::BEAMLINE = true;
0251
0252     //===============
0253     // conditions DB flags
0254     //===============
0255     recoConsts *rc = recoConsts::instance();
0256     Enable::CDB = true;
0257     // global tag
0258     rc->set_StringFlag("CDB_GLOBALTAG", CDB::global_tag);
0259     rc->set_uint64Flag("TIMESTAMP", CDB::timestamp);
0260     rc->set_IntFlag("RUNNUMBER", runNumber);
0261
0262     // Initialize the selected subsystems
0263     G4Init();
0264     if (!fromG4Hits) G4Setup(); // disable this if run from g4hits
0265
0266     Mbd_Reco();
0267     if (fromG4Hits)
0268     {
0269         Mvtx_Cells();
0270         Intt_Cells();
0271     }
0272
0273     TrackingInit();
0274
0275     if (fromG4Hits)
0276     {
0277         Mvtx_Clustering();
0278         Intt_Clustering();
0279     }
0280
0281     // Silicon Seeding
0282     bool doSiliconSeeding = true;
0283     Enable::TRACKING_VERBOSITY = 10;
0284     if (doSiliconSeeding)
0285     {
0286         if (UseActsSiliconSeeding)
0287         {
0288             Tracking_Reco_SiliconSeed_run2pp(); // default Acts silicon seeding, iterative
0289         }
0290         else // CA silicon seeding, one pass
0291         {
0292             auto silicon_Seeding = new PHCASiliconSeeding();
0293             silicon_Seeding->Verbosity(0);
0294             // silicon_Seeding->SetTrackMapName("SiliconTrackSeedContainer_CA");
0295             silicon_Seeding->SetMVTXStrobeIDRange(-1, 1);
0296             silicon_Seeding->SetLayerRange(0, 6);         // MVTX + INTT
0297             silicon_Seeding->SetSearchWindow(1.5, 0.2);   // SetSearchWindow(float eta_allowance, float phi_width) // default: float neighbor_phi_width = .02; float eta_allowance = 1.1
0298             silicon_Seeding->SetAlgoUseBestTriplet(true); // default: true
0299             // silicon_Seeding->SetRequireINTTConsistency(bool req); // default: true
0300             silicon_Seeding->SetMinMVTXClusters(2);
0301             silicon_Seeding->SetMinINTTClusters(1);
0302             silicon_Seeding->SetMinClustersPerTrack(3); //
0303             se->registerSubsystem(silicon_Seeding);
0304
0305             auto merger = new PHSiliconSeedMerger;
0306             merger->Verbosity(0);
0307             // merger->trackMapName("SiliconTrackSeedContainer_CA");
0308             se->registerSubsystem(merger);
0309         }
0310
0311         auto converter = new TrackSeedTrackMapConverter;
0312         // Default set to full SvtxTrackSeeds. Can be set to SiliconTrackSeedContainer or TpcTrackSeedContainer
0313         converter->setTrackSeedName("SiliconTrackSeedContainer");
0314         converter->setFieldMap(G4MAGNET::magfield_tracking);
0315         converter->Verbosity(0);
0316         se->registerSubsystem(converter);
0317
0318         // PHSimpleVertexFinder to find primary vertex using silicon-seeded tracks
0319         auto finder = new PHSimpleVertexFinder;
0320         finder->Verbosity(0);
0321         // finder->setTrackMapName("SvtxTrackMap_Acts");
0322         applyRecoScanConfig(finder, gRecoScanConfig);
0323         se->registerSubsystem(finder);
0324     }
0325
0326     // if (Enable::GLOBAL_RECO) // This does GlobalVertexReco (G4_Global.C)
0327     bool doGlobalReco = true;
0328     if (doGlobalReco) // run global reco regardless of the flag setting since we want to compare the vertexing performance with and without silicon seeding, and global reco is needed for that comparison
0329     {
0330         Global_Reco();
0331
0332         /*
0333         auto get_species_enum = [](const string &species_str) -> MinimumBiasInfo::SPECIES
0334         {
0335             if (species_str == "pp")
0336                 return MinimumBiasInfo::SPECIES::PP;
0337             else if (species_str == "OO")
0338                 return MinimumBiasInfo::SPECIES::OO;
0339             else if (species_str == "AuAu")
0340                 return MinimumBiasInfo::SPECIES::AUAU;
0341             else
0342             {
0343                 std::cerr << "Unknown species string: " << species_str << ", defaulting to PP" << std::endl;
0344                 return MinimumBiasInfo::SPECIES::PP;
0345             }
0346         };
0347
0348         MinimumBiasClassifier *mb = new MinimumBiasClassifier();
0349         mb->Verbosity(0);
0350         mb->setSpecies(get_species_enum(species));
0351         mb->setIsSim(true);
0352         mb->setOverwriteScale("/sphenix/user/dlis/Projects/centrality/cdb/calibrations/scales/cdb_centrality_scale_1.root");
0353         mb->setOverwriteVtx("/sphenix/user/dlis/Projects/centrality/cdb/calibrations/vertexscales/cdb_centrality_vertex_scale_1.root");
0354         se->registerSubsystem(mb);
0355         */
0356     }
0357
0358     // truth vertexing
0359     {
0360         auto vtxing = new PHTruthVertexing;
0361         vtxing->associate_tracks(false);
0362         se->registerSubsystem(vtxing);
0363     }
0364
0365     // Truth evaluator for reco-truth matching
0366     {
0367         auto truthRecoEval = new SvtxTruthRecoTableEval;
0368         truthRecoEval->Verbosity(0);
0369         se->registerSubsystem(truthRecoEval);
0370     }
0371
0372     // ntuplizer
0373     VertexCompareVerbosity::fillSilconSeed = 0;
0374     VertexCompareVerbosity::fillCluster = 0;
0375     VertexCompareVerbosity::fillTruthParticle = 0;
0376     bool doTruthMatching = true;
0377     auto compare = new VertexCompare();
0378     compare->IsSimulation();
0379     compare->setOutputName(compareOutputFile.c_str());
0380     if (doTruthMatching) compare->doTruthMatching();
0381     se->registerSubsystem(compare);
0382
0383     //--------------
0384     // Set up Input Managers
0385     //--------------
0386
0387     InputManagers();
0388
0389     // if we use a negative number of events we go back to the command line here
0390     if (nEvents < 0)
0391     {
0392         return 0;
0393     }
0394
0395     // se->skip(skip);
0396     se->run(nEvents);
0397     //  se->PrintTimer();
0398
0399     //-----
0400     // Exit
0401     //-----
0402
0403     CDBInterface::instance()->Print(); // print used DB files
0404     se->End();
0405     std::cout << "All done" << std::endl;
0406     delete se;
0407
0408     gSystem->Exit(0);
0409     return 0;
0410 }
0411 #endif