STrackMatchingComparator/macros/Fun4All_RunNewTruthMatcher.C

0001 #ifndef MACRO_FUN4ALLG4SPHENIX_C
0002 #define MACRO_FUN4ALLG4SPHENIX_C
0003
0004 #include <GlobalVariables.C>
0005
0006 #include <DisplayOn.C>
0007 #include <G4Setup_sPHENIX.C>
0008 #include <G4_Mbd.C>
0009 #include <G4_CaloTrigger.C>
0010 #include <G4_Centrality.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_KFParticle.C>
0017 #include <G4_ParticleFlow.C>
0018 #include <G4_Production.C>
0019 #include <G4_TopoClusterReco.C>
0020
0021 #include <Trkr_RecoInit.C>
0022 #include <Trkr_Clustering.C>
0023 #include <Trkr_LaserClustering.C>
0024 #include <Trkr_Reco.C>
0025 #include <Trkr_Eval.C>
0026 #include <Trkr_QA.C>
0027
0028 #include <Trkr_Diagnostics.C>
0029 #include <G4_User.C>
0030 #include <QA.C>
0031
0032 #include <ffamodules/FlagHandler.h>
0033 #include <ffamodules/HeadReco.h>
0034 #include <ffamodules/SyncReco.h>
0035 #include <ffamodules/CDBInterface.h>
0036
0037 #include <fun4all/Fun4AllDstOutputManager.h>
0038 #include <fun4all/Fun4AllOutputManager.h>
0039 #include <fun4all/Fun4AllServer.h>
0040
0041 #include <phool/PHRandomSeed.h>
0042 #include <phool/recoConsts.h>
0043
0044 R__LOAD_LIBRARY(libfun4all.so)
0045 R__LOAD_LIBRARY(libffamodules.so)
0046
0047 // For HepMC Hijing
0048 // try inputFile = /sphenix/sim/sim01/sphnxpro/sHijing_HepMC/sHijing_0-12fm.dat
0049
0050 int Fun4All_RunNewTruthMatcher(
0051     const int nEvents = 1,
0052     const int skip = 0,
0053     const string &outputFile = "G4sPHENIX.root",
0054     const string &inputFile = "https://www.phenix.bnl.gov/WWW/publish/phnxbld/sPHENIX/files/sPHENIX_G4Hits_sHijing_9-11fm_00000_00010.root",
0055     const string &embed_input_file = "https://www.phenix.bnl.gov/WWW/publish/phnxbld/sPHENIX/files/sPHENIX_G4Hits_sHijing_9-11fm_00000_00010.root",
0056     const string &outdir = ".")
0057 {
0058   Fun4AllServer *se = Fun4AllServer::instance();
0059   se->Verbosity(0);
0060
0061   //Opt to print all random seed used for debugging reproducibility. Comment out to reduce stdout prints.
0062   PHRandomSeed::Verbosity(1);
0063
0064   // just if we set some flags somewhere in this macro
0065   recoConsts *rc = recoConsts::instance();
0066   // By default every random number generator uses
0067   // PHRandomSeed() which reads /dev/urandom to get its seed
0068   // if the RANDOMSEED flag is set its value is taken as seed
0069   // You can either set this to a random value using PHRandomSeed()
0070   // which will make all seeds identical (not sure what the point of
0071   // this would be:
0072   //  rc->set_IntFlag("RANDOMSEED",PHRandomSeed());
0073   // or set it to a fixed value so you can debug your code
0074   //  rc->set_IntFlag("RANDOMSEED", 12345);
0075
0076
0077   //===============
0078   // Input options
0079   //===============
0080   // verbosity setting (applies to all input managers)
0081   Input::VERBOSITY = 0;
0082   // First enable the input generators
0083   // Either:
0084   // read previously generated g4-hits files, in this case it opens a DST and skips
0085   // the simulations step completely. The G4Setup macro is only loaded to get information
0086   // about the number of layers used for the cell reco code
0087   //  Input::READHITS = true;
0088   INPUTREADHITS::filename[0] = inputFile;
0089   // if you use a filelist
0090   // INPUTREADHITS::listfile[0] = inputFile;
0091   // Or:
0092   // Use particle generator
0093   // And
0094   // Further choose to embed newly simulated events to a previous simulation. Not compatible with `readhits = true`
0095   // In case embedding into a production output, please double check your G4Setup_sPHENIX.C and G4_*.C consistent with those in the production macro folder
0096   // E.g. /sphenix/sim//sim01/production/2016-07-21/single_particle/spacal2d/
0097   //  Input::EMBED = true;
0098   INPUTEMBED::filename[0] = embed_input_file;
0099   // if you use a filelist
0100   //INPUTEMBED::listfile[0] = embed_input_file;
0101
0102   Input::SIMPLE = true;
0103   // Input::SIMPLE_NUMBER = 2; // if you need 2 of them
0104   // Input::SIMPLE_VERBOSITY = 1;
0105
0106   // Enable this is emulating the nominal pp/pA/AA collision vertex distribution
0107   // Input::BEAM_CONFIGURATION = Input::AA_COLLISION; // Input::AA_COLLISION (default), Input::pA_COLLISION, Input::pp_COLLISION
0108
0109   //  Input::PYTHIA6 = true;
0110
0111   // Input::PYTHIA8 = true;
0112
0113   //  Input::GUN = true;
0114   //  Input::GUN_NUMBER = 3; // if you need 3 of them
0115   // Input::GUN_VERBOSITY = 1;
0116
0117   // Input::COSMIC = true;
0118
0119   //D0 generator
0120   //Input::DZERO = false;
0121   //Input::DZERO_VERBOSITY = 0;
0122   //Lambda_c generator //Not ready yet
0123   //Input::LAMBDAC = false;
0124   //Input::LAMBDAC_VERBOSITY = 0;
0125   // Upsilon generator
0126   //Input::UPSILON = true;
0127   //Input::UPSILON_NUMBER = 3; // if you need 3 of them
0128   //Input::UPSILON_VERBOSITY = 0;
0129
0130   //  Input::HEPMC = true;
0131   INPUTHEPMC::filename = inputFile;
0132   //-----------------
0133   // Hijing options (symmetrize hijing, add flow, add fermi motion)
0134   //-----------------
0135   //  INPUTHEPMC::HIJINGFLIP = true;
0136   //  INPUTHEPMC::FLOW = true;
0137   //  INPUTHEPMC::FLOW_VERBOSITY = 3;
0138   //  INPUTHEPMC::FERMIMOTION = true;
0139
0140
0141   // Event pile up simulation with collision rate in Hz MB collisions.
0142   //Input::PILEUPRATE = 50e3; // 50 kHz for AuAu
0143   //Input::PILEUPRATE = 3e6; // 3MHz for pp
0144
0145   // Enable this is emulating the nominal pp/pA/AA collision vertex distribution
0146   // for HepMC records (hijing, pythia8)
0147   //  Input::BEAM_CONFIGURATION = Input::AA_COLLISION; // for 2023 sims we want the AA geometry for no pileup sims
0148   //  Input::BEAM_CONFIGURATION = Input::pp_COLLISION; // for 2024 sims we want the pp geometry for no pileup sims
0149   //  Input::BEAM_CONFIGURATION = Input::pA_COLLISION; // for pAu sims we want the pA geometry for no pileup sims
0150
0151   //-----------------
0152   // Initialize the selected Input/Event generation
0153   //-----------------
0154   // This creates the input generator(s)
0155   InputInit();
0156
0157   //--------------
0158   // Set generator specific options
0159   //--------------
0160   // can only be set after InputInit() is called
0161
0162   // Simple Input generator:
0163   // if you run more than one of these Input::SIMPLE_NUMBER > 1
0164   // add the settings for other with [1], next with [2]...
0165
0166   // TEST [Derek, 09.05.2023]
0167   const UInt_t   nPerEvt(1);
0168   const TString  sParticle("pi-");
0169   const Double_t ptParMin(10.);
0170   const Double_t ptParMax(10.);
0171
0172   // TEST [Derek, 02.06.2024]
0173   const Bool_t   doGausVertex(false);
0174   const Double_t meanVtxXG(0.);
0175   const Double_t meanVtxYG(0.);
0176   const Double_t meanVtxZG(0.);
0177   const Double_t meanVtxXU(0.);
0178   const Double_t meanVtxYU(0.);
0179   const Double_t meanVtxZU(0.);
0180   const Double_t widthVtxXG(0.01);
0181   const Double_t widthVtxYG(0.01);
0182   const Double_t widthVtxZG(5.);
0183   const Double_t widthVtxXU(0.);
0184   const Double_t widthVtxYU(0.);
0185   const Double_t widthVtxZU(0.);
0186
0187   if (Input::SIMPLE)
0188   {
0189     INPUTGENERATOR::SimpleEventGenerator[0]->add_particles(sParticle.Data(), nPerEvt);
0190     if (Input::HEPMC || Input::EMBED)
0191     {
0192       INPUTGENERATOR::SimpleEventGenerator[0]->set_reuse_existing_vertex(true);
0193       INPUTGENERATOR::SimpleEventGenerator[0]->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
0194     }
0195     else
0196     {
0197       if (doGausVertex) {
0198         INPUTGENERATOR::SimpleEventGenerator[0] -> set_vertex_distribution_function(
0199           PHG4SimpleEventGenerator::Gaus,
0200           PHG4SimpleEventGenerator::Gaus,
0201           PHG4SimpleEventGenerator::Gaus
0202         );
0203         INPUTGENERATOR::SimpleEventGenerator[0] -> set_vertex_distribution_mean(meanVtxXG, meanVtxYG, meanVtxZG);
0204         INPUTGENERATOR::SimpleEventGenerator[0] -> set_vertex_distribution_width(widthVtxXG, widthVtxYG, widthVtxZG);
0205       } else {
0206         INPUTGENERATOR::SimpleEventGenerator[0] -> set_vertex_distribution_function(
0207           PHG4SimpleEventGenerator::Uniform,
0208           PHG4SimpleEventGenerator::Uniform,
0209           PHG4SimpleEventGenerator::Uniform
0210         );
0211         INPUTGENERATOR::SimpleEventGenerator[0] -> set_vertex_distribution_mean(meanVtxXU, meanVtxYU, meanVtxZU);
0212         INPUTGENERATOR::SimpleEventGenerator[0] -> set_vertex_distribution_width(widthVtxXU, widthVtxYU, widthVtxZU);
0213       }
0214     }
0215     INPUTGENERATOR::SimpleEventGenerator[0]->set_eta_range(-1, 1);
0216     INPUTGENERATOR::SimpleEventGenerator[0]->set_phi_range(-M_PI, M_PI);
0217     INPUTGENERATOR::SimpleEventGenerator[0]->set_pt_range(ptParMin, ptParMax);
0218   }
0219   // Upsilons
0220   // if you run more than one of these Input::UPSILON_NUMBER > 1
0221   // add the settings for other with [1], next with [2]...
0222   if (Input::UPSILON)
0223   {
0224     INPUTGENERATOR::VectorMesonGenerator[0]->add_decay_particles("e", 0);
0225     INPUTGENERATOR::VectorMesonGenerator[0]->set_rapidity_range(-1, 1);
0226     INPUTGENERATOR::VectorMesonGenerator[0]->set_pt_range(0., 10.);
0227     // Y species - select only one, last one wins
0228     INPUTGENERATOR::VectorMesonGenerator[0]->set_upsilon_1s();
0229     if (Input::HEPMC || Input::EMBED)
0230     {
0231       INPUTGENERATOR::VectorMesonGenerator[0]->set_reuse_existing_vertex(true);
0232       INPUTGENERATOR::VectorMesonGenerator[0]->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
0233     }
0234   }
0235   // particle gun
0236   // if you run more than one of these Input::GUN_NUMBER > 1
0237   // add the settings for other with [1], next with [2]...
0238   if (Input::GUN)
0239   {
0240     INPUTGENERATOR::Gun[0]->AddParticle("pi-", 0, 1, 0);
0241     INPUTGENERATOR::Gun[0]->set_vtx(0, 0, 0);
0242   }
0243
0244   // pythia6
0245   if (Input::PYTHIA6)
0246   {
0247     //! Nominal collision geometry is selected by Input::BEAM_CONFIGURATION
0248     Input::ApplysPHENIXBeamParameter(INPUTGENERATOR::Pythia6);
0249   }
0250   // pythia8
0251   if (Input::PYTHIA8)
0252   {
0253     //! Nominal collision geometry is selected by Input::BEAM_CONFIGURATION
0254     Input::ApplysPHENIXBeamParameter(INPUTGENERATOR::Pythia8);
0255   }
0256
0257   //--------------
0258   // Set Input Manager specific options
0259   //--------------
0260   // can only be set after InputInit() is called
0261
0262   if (Input::HEPMC)
0263   {
0264     //! Nominal collision geometry is selected by Input::BEAM_CONFIGURATION
0265     Input::ApplysPHENIXBeamParameter(INPUTMANAGER::HepMCInputManager);
0266
0267     // optional overriding beam parameters
0268     //INPUTMANAGER::HepMCInputManager->set_vertex_distribution_width(100e-4, 100e-4, 8, 0);  //optional collision smear in space, time
0269     //    INPUTMANAGER::HepMCInputManager->set_vertex_distribution_mean(0,0,0,0);//optional collision central position shift in space, time
0270     // //optional choice of vertex distribution function in space, time
0271     //INPUTMANAGER::HepMCInputManager->set_vertex_distribution_function(PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus);
0272     //! embedding ID for the event
0273     //! positive ID is the embedded event of interest, e.g. jetty event from pythia
0274     //! negative IDs are backgrounds, .e.g out of time pile up collisions
0275     //! Usually, ID = 0 means the primary Au+Au collision background
0276     //INPUTMANAGER::HepMCInputManager->set_embedding_id(Input::EmbedID);
0277     if (Input::PILEUPRATE > 0)
0278     {
0279       // Copy vertex settings from foreground hepmc input
0280       INPUTMANAGER::HepMCPileupInputManager->CopyHelperSettings(INPUTMANAGER::HepMCInputManager);
0281       // and then modify the ones you want to be different
0282       // INPUTMANAGER::HepMCPileupInputManager->set_vertex_distribution_width(100e-4,100e-4,8,0);
0283     }
0284   }
0285   if (Input::PILEUPRATE > 0)
0286   {
0287     //! Nominal collision geometry is selected by Input::BEAM_CONFIGURATION
0288     Input::ApplysPHENIXBeamParameter(INPUTMANAGER::HepMCPileupInputManager);
0289   }
0290   // register all input generators with Fun4All
0291   InputRegister();
0292
0293   if (! Input::READHITS)
0294   {
0295     rc->set_IntFlag("RUNNUMBER",1);
0296
0297     SyncReco *sync = new SyncReco();
0298     se->registerSubsystem(sync);
0299
0300     HeadReco *head = new HeadReco();
0301     se->registerSubsystem(head);
0302   }
0303 // Flag Handler is always needed to read flags from input (if used)
0304 // and update our rc flags with them. At the end it saves all flags
0305 // again on the DST in the Flags node under the RUN node
0306   FlagHandler *flag = new FlagHandler();
0307   se->registerSubsystem(flag);
0308
0309   // set up production relatedstuff
0310   //   Enable::PRODUCTION = true;
0311
0312   //======================
0313   // Write the DST
0314   //======================
0315
0316   //Enable::DSTOUT = true;
0317   Enable::DSTOUT_COMPRESS = false;
0318   DstOut::OutputDir = outdir;
0319   DstOut::OutputFile = outputFile;
0320
0321   //Option to convert DST to human command readable TTree for quick poke around the outputs
0322   //  Enable::DSTREADER = true;
0323
0324   // turn the display on (default off)
0325    //Enable::DISPLAY = true;
0326
0327   //======================
0328   // What to run
0329   //======================
0330
0331   // QA, main switch
0332   Enable::QA = false;
0333
0334   // Global options (enabled for all enables subsystems - if implemented)
0335   //  Enable::ABSORBER = true;
0336   //  Enable::OVERLAPCHECK = true;
0337   //  Enable::VERBOSITY = 1;
0338
0339   // Enable::MBD = true;
0340   // Enable::MBD_SUPPORT = true; // save hist in MBD/BBC support structure
0341   // Enable::MBDRECO = Enable::MBD && true;
0342   Enable::MBDFAKE = true;  // Smeared vtx and t0, use if you don't want real MBD/BBC in simulation
0343
0344   Enable::PIPE = true;
0345   Enable::PIPE_ABSORBER = true;
0346
0347   // central tracking
0348   Enable::MVTX = true;
0349   Enable::MVTX_CELL = Enable::MVTX && true;
0350   Enable::MVTX_CLUSTER = Enable::MVTX_CELL && true;
0351   Enable::MVTX_QA = Enable::MVTX_CLUSTER && Enable::QA && true;
0352
0353   Enable::INTT = true;
0354 //  Enable::INTT_ABSORBER = true; // enables layerwise support structure readout
0355 //  Enable::INTT_SUPPORT = true; // enable global support structure readout
0356   Enable::INTT_CELL = Enable::INTT && true;
0357   Enable::INTT_CLUSTER = Enable::INTT_CELL && true;
0358   Enable::INTT_QA = Enable::INTT_CLUSTER && Enable::QA && true;
0359
0360   Enable::TPC = true;
0361   Enable::TPC_ABSORBER = true;
0362   Enable::TPC_CELL = Enable::TPC && true;
0363   Enable::TPC_CLUSTER = Enable::TPC_CELL && true;
0364   Enable::TPC_QA = Enable::TPC_CLUSTER && Enable::QA && true;
0365
0366   Enable::MICROMEGAS = true;
0367   Enable::MICROMEGAS_CELL = Enable::MICROMEGAS && true;
0368   Enable::MICROMEGAS_CLUSTER = Enable::MICROMEGAS_CELL && true;
0369   Enable::MICROMEGAS_QA = Enable::MICROMEGAS_CLUSTER && Enable::QA && true;
0370
0371   Enable::TRACKING_TRACK = (Enable::MICROMEGAS_CLUSTER && Enable::TPC_CLUSTER && Enable::INTT_CLUSTER && Enable::MVTX_CLUSTER) && true;
0372   Enable::TRACKING_EVAL = Enable::TRACKING_TRACK && Enable::GLOBAL_RECO && true;
0373   Enable::TRACKING_QA = Enable::TRACKING_TRACK && Enable::QA && true;
0374
0375   // only do track matching if TRACKINGTRACK is also used
0376   Enable::TRACK_MATCHING = Enable::TRACKING_TRACK && true;
0377   Enable::TRACK_MATCHING_TREE = Enable::TRACK_MATCHING && true;
0378   Enable::TRACK_MATCHING_TREE_CLUSTERS = Enable::TRACK_MATCHING_TREE && true;
0379
0380   //Additional tracking tools
0381   //Enable::TRACKING_DIAGNOSTICS = Enable::TRACKING_TRACK && true;
0382   //G4TRACKING::filter_conversion_electrons = true;
0383
0384
0385   //  cemc electronics + thin layer of W-epoxy to get albedo from cemc
0386   //  into the tracking, cannot run together with CEMC
0387   //  Enable::CEMCALBEDO = true;
0388
0389   Enable::CEMC = true;
0390   Enable::CEMC_ABSORBER = true;
0391   Enable::CEMC_CELL = Enable::CEMC && true;
0392   Enable::CEMC_TOWER = Enable::CEMC_CELL && true;
0393   Enable::CEMC_CLUSTER = Enable::CEMC_TOWER && true;
0394   Enable::CEMC_EVAL = Enable::CEMC_G4Hit && Enable::CEMC_CLUSTER && false;
0395   Enable::CEMC_QA = Enable::CEMC_CLUSTER && Enable::QA && false;
0396
0397   Enable::HCALIN = true;
0398   Enable::HCALIN_ABSORBER = true;
0399   Enable::HCALIN_CELL = Enable::HCALIN && true;
0400   Enable::HCALIN_TOWER = Enable::HCALIN_CELL && true;
0401   Enable::HCALIN_CLUSTER = Enable::HCALIN_TOWER && true;
0402   Enable::HCALIN_EVAL = Enable::HCALIN_G4Hit && Enable::HCALIN_CLUSTER && false;
0403   Enable::HCALIN_QA = Enable::HCALIN_CLUSTER && Enable::QA && false;
0404
0405   Enable::MAGNET = true;
0406   Enable::MAGNET_ABSORBER = true;
0407
0408   Enable::HCALOUT = true;
0409   Enable::HCALOUT_ABSORBER = true;
0410   Enable::HCALOUT_CELL = Enable::HCALOUT && true;
0411   Enable::HCALOUT_TOWER = Enable::HCALOUT_CELL && true;
0412   Enable::HCALOUT_CLUSTER = Enable::HCALOUT_TOWER && true;
0413   Enable::HCALOUT_EVAL = Enable::HCALOUT_G4Hit && Enable::HCALOUT_CLUSTER && false;
0414   Enable::HCALOUT_QA = Enable::HCALOUT_CLUSTER && Enable::QA && false;
0415
0416   Enable::EPD = true;
0417   Enable::EPD_TILE = Enable::EPD && true;
0418
0419   Enable::BEAMLINE = true;
0420 //  Enable::BEAMLINE_ABSORBER = true;  // makes the beam line magnets sensitive volumes
0421 //  Enable::BEAMLINE_BLACKHOLE = true; // turns the beamline magnets into black holes
0422   Enable::ZDC = true;
0423 //  Enable::ZDC_ABSORBER = true;
0424 //  Enable::ZDC_SUPPORT = true;
0425   Enable::ZDC_TOWER = Enable::ZDC && true;
0426   Enable::ZDC_EVAL = Enable::ZDC_TOWER && true;
0427
0428   //! forward flux return plug door. Out of acceptance and off by default.
0429   //Enable::PLUGDOOR = true;
0430   Enable::PLUGDOOR_ABSORBER = true;
0431
0432   Enable::GLOBAL_RECO = (Enable::MBDFAKE || Enable::TRACKING_TRACK) && true;
0433   //Enable::GLOBAL_FASTSIM = true;
0434
0435   //Enable::KFPARTICLE = true;
0436   //Enable::KFPARTICLE_VERBOSITY = 1;
0437   //Enable::KFPARTICLE_TRUTH_MATCH = true;
0438   //Enable::KFPARTICLE_SAVE_NTUPLE = true;
0439
0440   Enable::CALOTRIGGER = Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER && false;
0441
0442   Enable::JETS = (Enable::GLOBAL_RECO || Enable::GLOBAL_FASTSIM) && false;
0443   Enable::JETS_EVAL = Enable::JETS && false;
0444   Enable::JETS_QA = Enable::JETS && Enable::QA && false;
0445
0446   // HI Jet Reco for p+Au / Au+Au collisions (default is false for
0447   // single particle / p+p-only simulations, or for p+Au / Au+Au
0448   // simulations which don't particularly care about jets)
0449   Enable::HIJETS = Enable::JETS && Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER && false;
0450
0451   // 3-D topoCluster reconstruction, potentially in all calorimeter layers
0452   Enable::TOPOCLUSTER = Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER && false;
0453   // particle flow jet reconstruction - needs topoClusters!
0454   Enable::PARTICLEFLOW = Enable::TOPOCLUSTER && false;
0455   // centrality reconstruction
0456   Enable::CENTRALITY = true;
0457
0458   // new settings using Enable namespace in GlobalVariables.C
0459   Enable::BLACKHOLE = true;
0460   //Enable::BLACKHOLE_SAVEHITS = false; // turn off saving of bh hits
0461   //Enable::BLACKHOLE_FORWARD_SAVEHITS = false; // disable forward/backward hits
0462   //BlackHoleGeometry::visible = true;
0463
0464   // run user provided code (from local G4_User.C)
0465   //Enable::USER = true;
0466
0467   //===============
0468   // conditions DB flags
0469   //===============
0470   Enable::CDB = true;
0471   // global tag
0472   rc->set_StringFlag("CDB_GLOBALTAG",CDB::global_tag);
0473   // 64 bit timestamp
0474   rc->set_uint64Flag("TIMESTAMP",CDB::timestamp);
0475   //---------------
0476   // World Settings
0477   //---------------
0478   //  G4WORLD::PhysicsList = "FTFP_BERT"; //FTFP_BERT_HP best for calo
0479   //  G4WORLD::WorldMaterial = "G4_AIR"; // set to G4_GALACTIC for material scans
0480
0481   //---------------
0482   // Magnet Settings
0483   //---------------
0484
0485   //  G4MAGNET::magfield =  string(getenv("CALIBRATIONROOT"))+ string("/Field/Map/sphenix3dbigmapxyz.root");  // default map from the calibration database
0486   //  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)
0487 //  G4MAGNET::magfield_rescale = 1.;  // make consistent with expected Babar field strength of 1.4T
0488
0489   //---------------
0490   // Pythia Decayer
0491   //---------------
0492   // list of decay types in
0493   // $OFFLINE_MAIN/include/g4decayer/EDecayType.hh
0494   // default is All:
0495   // G4P6DECAYER::decayType = EDecayType::kAll;
0496
0497   // Initialize the selected subsystems
0498   G4Init();
0499
0500   //---------------------
0501   // GEANT4 Detector description
0502   //---------------------
0503   if (!Input::READHITS)
0504   {
0505     G4Setup();
0506   }
0507
0508   //------------------
0509   // Detector Division
0510   //------------------
0511
0512   if ((Enable::MBD && Enable::MBDRECO) || Enable::MBDFAKE) Mbd_Reco();
0513
0514   if (Enable::MVTX_CELL) Mvtx_Cells();
0515   if (Enable::INTT_CELL) Intt_Cells();
0516   if (Enable::TPC_CELL) TPC_Cells();
0517   if (Enable::MICROMEGAS_CELL) Micromegas_Cells();
0518
0519   if (Enable::CEMC_CELL) CEMC_Cells();
0520
0521   if (Enable::HCALIN_CELL) HCALInner_Cells();
0522
0523   if (Enable::HCALOUT_CELL) HCALOuter_Cells();
0524
0525   //-----------------------------
0526   // CEMC towering and clustering
0527   //-----------------------------
0528
0529   if (Enable::CEMC_TOWER) CEMC_Towers();
0530   if (Enable::CEMC_CLUSTER) CEMC_Clusters();
0531
0532   //--------------
0533   // EPD tile reconstruction
0534   //--------------
0535
0536   if (Enable::EPD_TILE) EPD_Tiles();
0537
0538   //-----------------------------
0539   // HCAL towering and clustering
0540   //-----------------------------
0541
0542   if (Enable::HCALIN_TOWER) HCALInner_Towers();
0543   if (Enable::HCALIN_CLUSTER) HCALInner_Clusters();
0544
0545   if (Enable::HCALOUT_TOWER) HCALOuter_Towers();
0546   if (Enable::HCALOUT_CLUSTER) HCALOuter_Clusters();
0547
0548   // if enabled, do topoClustering early, upstream of any possible jet reconstruction
0549   if (Enable::TOPOCLUSTER) TopoClusterReco();
0550
0551   //--------------
0552   // SVTX tracking
0553   //--------------
0554   if(Enable::TRACKING_TRACK)
0555     {
0556       TrackingInit();
0557     }
0558   if (Enable::MVTX_CLUSTER) Mvtx_Clustering();
0559   if (Enable::INTT_CLUSTER) Intt_Clustering();
0560   if (Enable::TPC_CLUSTER)
0561     {
0562       if(G4TPC::ENABLE_DIRECT_LASER_HITS || G4TPC::ENABLE_CENTRAL_MEMBRANE_HITS)
0563     {
0564       TPC_LaserClustering();
0565     }
0566       else
0567     {
0568       TPC_Clustering();
0569     }
0570     }
0571   if (Enable::MICROMEGAS_CLUSTER) Micromegas_Clustering();
0572
0573   if (Enable::TRACKING_TRACK)
0574   {
0575     Tracking_Reco();
0576   }
0577
0578
0579
0580   if(Enable::TRACKING_DIAGNOSTICS)
0581     {
0582       const std::string kshortFile = "./kshort_" + outputFile;
0583       const std::string residualsFile = "./residuals_" + outputFile;
0584
0585       G4KshortReconstruction(kshortFile);
0586       seedResiduals(residualsFile);
0587     }
0588
0589
0590   //-----------------
0591   // Global Vertexing
0592   //-----------------
0593
0594   if (Enable::GLOBAL_RECO && Enable::GLOBAL_FASTSIM)
0595   {
0596     cout << "You can only enable Enable::GLOBAL_RECO or Enable::GLOBAL_FASTSIM, not both" << endl;
0597     gSystem->Exit(1);
0598   }
0599   if (Enable::GLOBAL_RECO)
0600   {
0601     Global_Reco();
0602   }
0603   else if (Enable::GLOBAL_FASTSIM)
0604   {
0605     Global_FastSim();
0606   }
0607
0608   //-----------------
0609   // Centrality Determination
0610   //-----------------
0611
0612   if (Enable::CENTRALITY)
0613   {
0614       Centrality();
0615   }
0616
0617   //-----------------
0618   // Calo Trigger Simulation
0619   //-----------------
0620
0621   if (Enable::CALOTRIGGER)
0622   {
0623     CaloTrigger_Sim();
0624   }
0625
0626   //---------
0627   // Jet reco
0628   //---------
0629
0630   if (Enable::JETS) Jet_Reco();
0631   if (Enable::HIJETS) HIJetReco();
0632
0633   if (Enable::PARTICLEFLOW) ParticleFlow();
0634
0635   //----------------------
0636   // Simulation evaluation
0637   //----------------------
0638   string outputroot = outputFile;
0639   string remove_this = ".root";
0640   size_t pos = outputroot.find(remove_this);
0641   if (pos != string::npos)
0642   {
0643     outputroot.erase(pos, remove_this.length());
0644   }
0645
0646   //if (Enable::TRACKING_EVAL) {
0647   if (true) {
0648
0649     build_truthreco_tables();
0650
0651     std::cout << "CHECK I am inside the svtx evaluator block!" << std::endl;
0652
0653     SvtxEvaluator* eval = new SvtxEvaluator(
0654       "SVTXEVALUATOR",
0655       outputroot + "_g4svtx_eval.root",
0656       "SvtxTrackMap",
0657       G4MVTX::n_maps_layer,
0658       G4INTT::n_intt_layer,
0659       G4TPC::n_gas_layer,
0660       G4MICROMEGAS::n_micromegas_layer
0661     );
0662     eval->do_cluster_eval(true);
0663     eval->do_g4hit_eval(false);
0664     eval->do_hit_eval(false);  // enable to see the hits that includes the chamber physics...
0665     eval->do_gpoint_eval(true);
0666     eval->do_vtx_eval_light(true);
0667     eval->do_eval_light(true);
0668     eval->do_track_eval(true);
0669     eval->do_gtrack_eval(true);
0670     eval->do_track_match(true);
0671     eval->set_use_initial_vertex(G4TRACKING::g4eval_use_initial_vertex);
0672
0673     bool embed_scan = true;
0674     if(TRACKING::pp_mode) {
0675       embed_scan = false;
0676     }
0677     eval->scan_for_embedded(embed_scan);   // take all tracks if false - take only embedded tracks if true
0678     eval->scan_for_primaries(embed_scan);  // defaults to only thrown particles for ntp_gtrack
0679     std::cout << "SvtxEvaluator: pp_mode set to " << TRACKING::pp_mode << " and scan_for_embedded set to " << embed_scan << std::endl;
0680     eval->Verbosity(Enable::TRACKING_VERBOSITY);
0681     se->registerSubsystem(eval);
0682   }  // end if (Enable::TRACKING_EVAL)
0683
0684   if (Enable::CEMC_EVAL) CEMC_Eval(outputroot + "_g4cemc_eval.root");
0685
0686   if (Enable::HCALIN_EVAL) HCALInner_Eval(outputroot + "_g4hcalin_eval.root");
0687
0688   if (Enable::HCALOUT_EVAL) HCALOuter_Eval(outputroot + "_g4hcalout_eval.root");
0689
0690   if (Enable::JETS_EVAL) Jet_Eval(outputroot + "_g4jet_eval.root");
0691
0692   if (Enable::DSTREADER) G4DSTreader(outputroot + "_DSTReader.root");
0693
0694
0695
0696   if (Enable::USER) UserAnalysisInit();
0697
0698   // Writes electrons from conversions to a new track map on the node tree
0699   // the ntuple file is for diagnostics, it is produced only if the flag is set in G4_Tracking.C
0700   if(G4TRACKING::filter_conversion_electrons) Filter_Conversion_Electrons(outputroot + "_secvert_ntuple.root");
0701
0702
0703   //======================
0704   // Run KFParticle on evt
0705   //======================
0706   if (Enable::KFPARTICLE && Input::UPSILON) KFParticle_Upsilon_Reco();
0707   if (Enable::KFPARTICLE && Input::DZERO) KFParticle_D0_Reco();
0708
0709   //----------------------
0710   // Standard QAs
0711   //----------------------
0712
0713   if (Enable::CEMC_QA) CEMC_QA();
0714   if (Enable::HCALIN_QA) HCALInner_QA();
0715   if (Enable::HCALOUT_QA) HCALOuter_QA();
0716
0717   if (Enable::JETS_QA) Jet_QA();
0718
0719   if (Enable::MVTX_QA) Mvtx_QA();
0720   if (Enable::INTT_QA) Intt_QA();
0721   if (Enable::TPC_QA) TPC_QA();
0722   if (Enable::MICROMEGAS_QA) Micromegas_QA();
0723   if (Enable::TRACKING_QA) Tracking_QA();
0724
0725   if (Enable::TRACKING_QA && Enable::CEMC_QA && Enable::HCALIN_QA && Enable::HCALOUT_QA) QA_G4CaloTracking();
0726
0727   //if (Enable::TRACK_MATCHING) {
0728   if (true) {
0729
0730     // initialize track matcher
0731     TrkrClusterIsMatcher* ismatcher = new TrkrClusterIsMatcher();
0732
0733     // These are the default values -- uncomment and change as desired
0734     //ismatcher->single_pixel_phi_MVTX = false ; // default to pitch*max(N_pixels_M,N_pixels_T)*tol_MVTX
0735     //ismatcher->single_pixel_phi_INTT = false ; // ... same as for MVTX
0736     //ismatcher->single_bin_phi_TPC    = true  ;  // default to pitch*tol_phi_TPC
0737     //ismatcher->single_pixel_z_MVTX = false ; // default to pitch*max(N_pixels_M,N_pixels_T)*tol_z_MVTX
0738     //ismatcher->single_pixel_z_INTT = false ; // ... same as for MVTX
0739     //ismatcher->single_bin_t_TPC    = true  ;  // default to pitch*tol_t_TPC
0740     //ismatcher-> tol_phi_MVTX = 0.5;
0741     //ismatcher-> tol_phi_INTT = 0.5;
0742     //ismatcher-> tol_phi_TPC  = 1.0;
0743     //ismatcher-> tol_z_MVTX  = 0.5;
0744     //ismatcher-> tol_t_TPC   = 1.0;
0745
0746     auto trackmatcher = new TruthRecoTrackMatching(ismatcher);
0747     trackmatcher->set_min_cl_match     (5);    // minimum number of matched clusters to make a matched track
0748     trackmatcher->set_min_cl_ratio     (0.1);  // at least 10% of truth clusters must be matched
0749     trackmatcher->set_cutoff_deta      (0.3);  // won't compare tracks with |Δeta|>0.3 away
0750     trackmatcher->set_cutoff_dphi      (0.3);  // won't compare tracks with |Δphi|>0.3 away
0751     trackmatcher->set_smallsearch_deta (0.05); // will first compare tracks within this |Δphi|
0752     trackmatcher->set_smallsearch_dphi (0.05); // will first compare tracks within this |Δeta|
0753     //trackmatcher->set_max_nreco_per_truth (4); // maximum reco tracks matched for any truth track
0754     //trackmatcher->set_max_ntruth_per_reco (4); // maximum truth tracks matched for any reco track
0755     trackmatcher->set_max_nreco_per_truth (1); // maximum reco tracks matched for any truth track // TEST [Derek, 10.13.2023]
0756     trackmatcher->set_max_ntruth_per_reco (1); // maximum truth tracks matched for any reco track // TEST [Derek, 10.13.2023]
0757
0758     // set verbosity
0759     int verbosity = std::max(Enable::VERBOSITY, Enable::TRACKING_VERBOSITY);
0760     trackmatcher->Verbosity(verbosity);
0761     se->registerSubsystem(trackmatcher);
0762
0763     if (Enable::TRACK_MATCHING_TREE) {
0764       auto treefiller = new FillClusMatchTree(ismatcher, outputroot + "_g4trackmatching.root");//, true, true, true, outputFile);
0765       treefiller->Verbosity(verbosity);
0766       if (Enable::TRACK_MATCHING_TREE_CLUSTERS) {
0767         treefiller->m_fill_clusters     = true;
0768         treefiller->m_fill_SvUnmatched  = true;
0769       } else {
0770         treefiller->m_fill_clusters     = false;
0771         treefiller->m_fill_SvUnmatched  = false;
0772       }
0773       treefiller->m_fill_clusverbose = false;
0774       se->registerSubsystem(treefiller);
0775     }
0776   }  // end if (ENABLE::TRACK_MATCHING)
0777
0778   //--------------
0779   // Set up Input Managers
0780   //--------------
0781
0782   InputManagers();
0783
0784   if (Enable::PRODUCTION)
0785   {
0786     Production_CreateOutputDir();
0787   }
0788
0789   if (Enable::DSTOUT)
0790   {
0791     string FullOutFile = DstOut::OutputDir + "/" + DstOut::OutputFile;
0792     Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", FullOutFile);
0793     if (Enable::DSTOUT_COMPRESS)
0794     {
0795       ShowerCompress();
0796       DstCompress(out);
0797     }
0798     se->registerOutputManager(out);
0799   }
0800   //-----------------
0801   // Event processing
0802   //-----------------
0803   if (Enable::DISPLAY)
0804   {
0805     DisplayOn();
0806
0807     gROOT->ProcessLine("Fun4AllServer *se = Fun4AllServer::instance();");
0808     gROOT->ProcessLine("PHG4Reco *g4 = (PHG4Reco *) se->getSubsysReco(\"PHG4RECO\");");
0809
0810     cout << "-------------------------------------------------" << endl;
0811     cout << "You are in event display mode. Run one event with" << endl;
0812     cout << "se->run(1)" << endl;
0813     cout << "Run Geant4 command with following examples" << endl;
0814     gROOT->ProcessLine("displaycmd()");
0815
0816     return 0;
0817   }
0818
0819   // if we use a negative number of events we go back to the command line here
0820   if (nEvents < 0)
0821   {
0822     return 0;
0823   }
0824   // if we run the particle generator and use 0 it'll run forever
0825   // for embedding it runs forever if the repeat flag is set
0826   if (nEvents == 0 && !Input::HEPMC && !Input::READHITS && INPUTEMBED::REPEAT)
0827   {
0828     cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
0829     cout << "it will run forever, so I just return without running anything" << endl;
0830     return 0;
0831   }
0832
0833   se->skip(skip);
0834   se->run(nEvents);
0835
0836   //-----
0837   // QA output
0838   //-----
0839
0840   if (Enable::QA) QA_Output(outputroot + "_qa.root");
0841
0842   //-----
0843   // Exit
0844   //-----
0845
0846 //  CDBInterface::instance()->Print(); // print used DB files
0847   se->End();
0848   std::cout << "All done" << std::endl;
0849   delete se;
0850   if (Enable::PRODUCTION)
0851   {
0852     Production_MoveOutput();
0853   }
0854
0855   gSystem->Exit(0);
0856   return 0;
0857 }
0858 #endif