SiCalo/PYHTIAGen/Fun4All_PYHTIAGen_Silicon.C

0001 #ifndef MACRO_FUN4ALLG4SPHENIX_C
0002 #define MACRO_FUN4ALLG4SPHENIX_C
0003
0004 #include <sys/stat.h>
0005 #include <limits.h>     // PATH_MAX
0006 #include <unistd.h>
0007 #include <GlobalVariables.C>
0008
0009 #include <DisplayOn.C>
0010 #include <G4Setup_sPHENIX.C>
0011 #include <G4_Mbd.C>
0012 #include <G4_CaloTrigger.C>
0013 #include <G4_Centrality.C>
0014 #include <G4_DSTReader.C>
0015 #include <G4_Global.C>
0016 #include <G4_HIJetReco.C>
0017 #include <G4_Input.C>
0018 #include <G4_Jets.C>
0019 #include <G4_KFParticle.C>
0020 #include <G4_ParticleFlow.C>
0021 #include <G4_Production.C>
0022 #include <G4_TopoClusterReco.C>
0023
0024 #include <Trkr_RecoInit.C>
0025 #include <Trkr_Clustering.C>
0026 #include <Trkr_Reco.C>
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 #include <trackreco/PHActsTrackProjection.h>
0034 #include <phpythia8/PHPy8ParticleTrigger.h>
0035 #include <phpythia8/PHPy8JetTrigger.h>
0036
0037 #include <ffamodules/FlagHandler.h>
0038 #include <ffamodules/HeadReco.h>
0039 #include <ffamodules/SyncReco.h>
0040 #include <ffamodules/CDBInterface.h>
0041 #include <fun4allutils/TimerStats.h>
0042 #include <fun4all/Fun4AllDstOutputManager.h>
0043 #include <fun4all/Fun4AllOutputManager.h>
0044 #include <fun4all/Fun4AllServer.h>
0045
0046 #include <phool/PHRandomSeed.h>
0047 #include <phool/recoConsts.h>
0048
0049 #include <simqa_modules/QAG4SimulationTracking.h>
0050 #include <qautils/QAHistManagerDef.h>
0051 #include <siliconseedsana/SiliconSeedsAna.h>
0052 #pragma GCC diagnostic push
0053
0054 #pragma GCC diagnostic ignored "-Wundefined-internal"
0055
0056
0057 #pragma GCC diagnostic pop
0058
0059 R__LOAD_LIBRARY(libfun4all.so)
0060 R__LOAD_LIBRARY(libffamodules.so)
0061 R__LOAD_LIBRARY(libsimqa_modules.so)
0062 R__LOAD_LIBRARY(libsiliconseedsana.so)
0063 R__LOAD_LIBRARY(libtrack_reco.so)
0064
0065 void ensure_dir(const std::string& path)
0066 {
0067   struct stat info;
0068
0069   if (stat(path.c_str(), &info) != 0)
0070   {
0071     std::cout << "Directory " << path << " does not exist. Creating..." << std::endl;
0072     if (mkdir(path.c_str(), 0777) != 0)
0073     {
0074       std::cerr << "Failed to create directory " << path << std::endl;
0075       exit(1);
0076     }
0077   }
0078   else if (!(info.st_mode & S_IFDIR))
0079   {
0080     std::cerr << "Path " << path << " exists but is not a directory!" << std::endl;
0081     exit(1);
0082   }
0083 };
0084
0085
0086 bool useTopologicalCluster = false;
0087 bool jpsiTodielectronOnly = true;
0088 int Fun4All_PYHTIAGen_Silicon(std::string processID = "0")
0089 {
0090   const int nEvents = 20;
0091
0092   std::ostringstream pid;
0093   pid << std::setw(6) << std::setfill('0') << std::stoi(processID);
0094   std::string pid_str = pid.str();
0095
0096   char cwd[PATH_MAX];
0097   if (getcwd(cwd, sizeof(cwd)) == nullptr)
0098   {
0099     std::cerr << "Failed to get current working directory!" << std::endl;
0100     return 1;
0101   }
0102
0103   //  std::string outDir = "/sphenix/user/jaein213/tracking/SiliconSeeding/MC/macro/DST/" + particle_name_tag;
0104   //  std::string outDir2 = "/sphenix/user/jaein213/tracking/SiliconSeeding/MC/macro/ana/" + particle_name_tag;
0105   std::string baseDir(cwd);
0106   std::string outDir = baseDir + "/DST";
0107   std::string outDir2 = baseDir + "/ana";
0108   std::string outDir3 = baseDir + "/jobtime";
0109   ensure_dir(outDir);
0110   ensure_dir(outDir2);
0111   ensure_dir(outDir3);
0112   ensure_dir(outDir + "/qa");
0113   Fun4AllServer *se = Fun4AllServer::instance();
0114   se->Verbosity(10);
0115
0116   // Opt to print all random seed used for debugging reproducibility. Comment out to reduce stdout prints.
0117   PHRandomSeed::Verbosity(1);
0118
0119   recoConsts *rc = recoConsts::instance();
0120   Input::VERBOSITY = 0;
0121   //===============
0122   // Input options
0123   //===============
0124   Input::PYTHIA8 = true;
0125   PYTHIA8::config_file = "PYHTIA8_JPSI_DielectronOnly.cfg";
0126
0127   InputInit();
0128
0129   // can only be set after InputInit() is called
0130   // pythia6
0131   if (Input::PYTHIA6)
0132   {
0133     //! Nominal collision geometry is selected by Input::BEAM_CONFIGURATION
0134     Input::ApplysPHENIXBeamParameter(INPUTGENERATOR::Pythia6);
0135   }
0136   // pythia8
0137   if (Input::PYTHIA8)
0138   {
0139     //! Nominal collision geometry is selected by Input::BEAM_CONFIGURATION
0140
0141     PHPy8ParticleTrigger *p8_trigger = new PHPy8ParticleTrigger();
0142     p8_trigger->AddParticles(443); // J/psi
0143     p8_trigger->SetStableParticleOnly(false);
0144
0145     INPUTGENERATOR::Pythia8->register_trigger(p8_trigger);
0146     INPUTGENERATOR::Pythia8->set_trigger_AND();
0147
0148     /////////////////////// Part to Tune the Phytia beam parameters started ////////////////////
0149     // Option 1) You can use some default configuration provided by sPHENIX
0150     Input::ApplysPHENIXBeamParameter(INPUTGENERATOR::Pythia8, Input::pp_COLLISION);
0151     //    Input::ApplysPHENIXBeamParameter(INPUTGENERATOR::Pythia8,Input::pp_ZEROANGLE);
0152     //    Input::ApplysPHENIXBeamParameter(INPUTGENERATOR::Pythia8,Input::AA_COLLISION);
0153
0154     // Option 2) You can set the beam crossing or set the vertex_distribution_width at the INPUT generator level
0155     if (false)
0156     {
0157       Input::beam_crossing = 1.;
0158       double localbcross = Input::beam_crossing / 2. * 1e-3;
0159       localbcross = Input::beam_crossing / 2. * 1e-3;
0160       //  Xing angle is split among both beams, means set to 0.5 mRad
0161       INPUTGENERATOR::Pythia8->set_beam_direction_theta_phi(localbcross, 0, M_PI - localbcross, 0); // 1.5mrad x-ing of sPHENIX
0162       INPUTGENERATOR::Pythia8->set_vertex_distribution_width(
0163         120e-4,         // approximation from past PHENIX data
0164         120e-4,         // approximation from past PHENIX data
0165         16,             // measured in 2024 for 1.5mrad Xing angle
0166         20 / 29.9792);  // 20cm collision length / speed of light in cm/ns
0167     }
0168     /////////////////////// Part to Tune the Phytia beam parameters done ////////////////////
0169   }
0170
0171   if (Input::PILEUPRATE > 0)
0172   {
0173     //! Nominal collision geometry is selected by Input::BEAM_CONFIGURATION
0174     Input::ApplysPHENIXBeamParameter(INPUTMANAGER::HepMCPileupInputManager);
0175   }
0176   // register all input generators with Fun4All
0177
0178   InputRegister();
0179
0180   if (!Input::READHITS)
0181   {
0182     rc->set_IntFlag("RUNNUMBER", 1);
0183
0184     SyncReco *sync = new SyncReco();
0185     se->registerSubsystem(sync);
0186
0187     HeadReco *head = new HeadReco();
0188     se->registerSubsystem(head);
0189   }
0190   FlagHandler *flag = new FlagHandler();
0191   se->registerSubsystem(flag);
0192   TimerStats *ts = new TimerStats();
0193   std::string jobtimeFile = outDir3 + "/jobtime_" + processID + ".root";
0194   ts->OutFileName(jobtimeFile);
0195   se->registerSubsystem(ts);
0196   // Simulation setup
0197   Enable::MBDFAKE = true;
0198   Enable::PIPE = true;
0199   Enable::PIPE_ABSORBER = true;
0200   Enable::MVTX = true;
0201   Enable::INTT = true;
0202   Enable::TPC = true;
0203   Enable::MICROMEGAS = true;
0204
0205   Enable::MAGNET = true;
0206   Enable::MAGNET_ABSORBER = true;
0207
0208   Enable::PLUGDOOR_ABSORBER = true;
0209
0210   Enable::CEMC = true;
0211   Enable::CEMC_ABSORBER = true;
0212   Enable::CEMC_CELL = true;
0213   Enable::CEMC_TOWER = true;
0214   Enable::CEMC_CLUSTER = true;
0215   // Enable::CEMC_EVAL =  false;
0216   // Enable::CEMC_QA =  false;
0217
0218   Enable::HCALIN = true;
0219   Enable::HCALIN_ABSORBER = true;
0220   Enable::HCALIN_CELL = true;
0221   Enable::HCALIN_TOWER = true;
0222   // Enable::HCALIN_CLUSTER =  true;
0223   // Enable::HCALIN_EVAL =  true;
0224   // Enable::HCALIN_QA =  true;
0225
0226   Enable::HCALOUT = true;
0227   Enable::HCALOUT_ABSORBER = true;
0228   Enable::HCALOUT_CELL = true;
0229   Enable::HCALOUT_TOWER = true;
0230   // Enable::HCALOUT_CLUSTER = false;
0231   // Enable::HCALOUT_EVAL =  false;
0232   // Enable::HCALOUT_QA =  false;
0233
0234   // Tracking setup
0235
0236   Enable::CDB = true;
0237   rc->set_StringFlag("CDB_GLOBALTAG", CDB::global_tag);
0238   rc->set_uint64Flag("TIMESTAMP", CDB::timestamp);
0239
0240   MagnetInit();
0241
0242   G4Init();
0243   G4Setup();
0244
0245   Mbd_Reco();
0246   Mvtx_Cells();
0247   Intt_Cells();
0248   // TPC_Cells();
0249   // Micromegas_Cells();
0250
0251   // TrackingInit();
0252   ACTSGEOM::ActsGeomInit();
0253   Mvtx_Clustering();
0254   Intt_Clustering();
0255   // TPC_Clustering();
0256   // Micromegas_Clustering();
0257
0258   // // Calo setup
0259   CEMC_Cells();
0260   CEMC_Towers();
0261   CEMC_Clusters();
0262   TopoClusterReco();
0263   HCALInner_Cells();
0264   HCALOuter_Cells();
0265   HCALInner_Towers();
0266   // HCALInner_Clusters();
0267   HCALOuter_Towers();
0268   //  HCALOuter_Clusters();
0269
0270   auto silicon_Seeding = new PHActsSiliconSeeding;
0271   // silicon_Seeding->set_track_map_name("SvtxTrackSeedContainer");
0272   silicon_Seeding->Verbosity(0);
0273   silicon_Seeding->setinttRPhiSearchWindow(0.2);
0274   silicon_Seeding->setinttZSearchWindow(1.0);
0275   silicon_Seeding->seedAnalysis(false);
0276   se->registerSubsystem(silicon_Seeding);
0277
0278   auto merger = new PHSiliconSeedMerger;
0279   // merger->trackMapName("SvtxTrackSeedContainer");
0280   merger->Verbosity(0);
0281   se->registerSubsystem(merger);
0282
0283   auto converter = new TrackSeedTrackMapConverter;
0284   // SiliconTrackSeedContainer or TpcTrackSeedContainer
0285   converter->setTrackSeedName("SiliconTrackSeedContainer");
0286   converter->setFieldMap(G4MAGNET::magfield_tracking);
0287   converter->Verbosity(10);
0288   se->registerSubsystem(converter);
0289   PHSimpleVertexFinder *finder = new PHSimpleVertexFinder;
0290   finder->Verbosity(0);
0291   finder->setDcaCut(0.1);
0292   finder->setTrackPtCut(0.);
0293   finder->setBeamLineCut(1);
0294   finder->setTrackQualityCut(20);
0295   finder->setNmvtxRequired(3);
0296   finder->setOutlierPairCut(0.1);
0297   se->registerSubsystem(finder);
0298
0299   Global_Reco(); // Definition is below..
0300   /////////////////////
0301   // void Global_Reco()
0302   // {
0303   //   int verbosity = std::max(Enable::VERBOSITY, Enable::GLOBAL_VERBOSITY);
0304
0305   //   Fun4AllServer* se = Fun4AllServer::instance();
0306
0307   //   GlobalVertexReco* gblvertex = new GlobalVertexReco();
0308   //   gblvertex->Verbosity(verbosity);
0309   //   se->registerSubsystem(gblvertex);
0310
0311   //   return;
0312   // }
0313
0314   build_truthreco_tables(); // Definition is below..
0315   //////////////////////////////////
0316   // void build_truthreco_tables()
0317   // {
0318   //   int verbosity = std::max(Enable::VERBOSITY, Enable::TRACKING_VERBOSITY);
0319   //   Fun4AllServer* se = Fun4AllServer::instance();
0320
0321   //   // this module builds high level truth track association table.
0322   //   // If this module is used, this table should be called before any evaluator calls.
0323   //   // Removing this module, evaluation will still work but trace truth association through the layers of G4-hit-cluster
0324   //   SvtxTruthRecoTableEval *tables = new SvtxTruthRecoTableEval();
0325   //   tables->Verbosity(verbosity);
0326   //   se->registerSubsystem(tables);
0327
0328   //   return;
0329   // }
0330   /////////////////////////////////
0331   bool runTruth = false;
0332   auto ensure_dir = [](const std::string &path)
0333   {
0334     struct stat info;
0335     if (stat(path.c_str(), &info) != 0)
0336     {
0337       std::cout << "Directory " << path << " does not exist. Creating..." << std::endl;
0338       mkdir(path.c_str(), 0777);
0339     }
0340     else if (!(info.st_mode & S_IFDIR))
0341     {
0342       std::cerr << "Path " << path << " exists but is not a directory!" << std::endl;
0343       exit(1);
0344     }
0345   };
0346
0347   std::string outputName = outDir + "/DST_SiliconOnly_PHYTIAGen_";
0348   if (runTruth)
0349     outputName += "truth";
0350   else
0351     outputName += "reconstructed";
0352   outputName += "Info_" + processID + ".root";
0353   std::string outputName2 = outDir2 + "/ana_" + processID + ".root";
0354
0355   Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", outputName);
0356   se->registerOutputManager(out);
0357
0358   auto projection = new PHActsTrackProjection("CaloProjection");
0359   se->registerSubsystem(projection);
0360   QAG4SimulationTracking *qa = new QAG4SimulationTracking();
0361   qa->Verbosity(5);
0362   se->registerSubsystem(qa);
0363   auto siana = new SiliconSeedsAna("SiliconSeedsAna");
0364   siana->setMC(true);
0365   siana->setVtxSkip(true);
0366   siana->setTopoCluster(useTopologicalCluster);
0367   siana->setOutputFileName(outputName2);
0368   int startnumber = nEvents * std::stoi(processID);
0369   siana->setStartEventNumber(startnumber);
0370   se->registerSubsystem(siana);
0371   se->run(nEvents);
0372
0373   std::string qaFile = outDir + "/qa/PHYTIAGen_JPSI" + processID + ".root";
0374   QAHistManagerDef::saveQARootFile(qaFile);
0375
0376   //-----
0377   // Exit
0378   //-----
0379
0380   //  CDBInterface::instance()->Print(); // print used DB files
0381   se->End();
0382   std::cout << "All done" << std::endl;
0383   delete se;
0384   if (Enable::PRODUCTION)
0385   {
0386     Production_MoveOutput();
0387   }
0388
0389   gSystem->Exit(0);
0390   return 0;
0391 }
0392 #endif