SiCalo/gunmacro/Fun4All_singleParticle_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
0035 #include <ffamodules/FlagHandler.h>
0036 #include <ffamodules/HeadReco.h>
0037 #include <ffamodules/SyncReco.h>
0038 #include <ffamodules/CDBInterface.h>
0039 #include <fun4allutils/TimerStats.h>
0040 #include <fun4all/Fun4AllDstOutputManager.h>
0041 #include <fun4all/Fun4AllOutputManager.h>
0042 #include <fun4all/Fun4AllServer.h>
0043
0044 #include <phool/PHRandomSeed.h>
0045 #include <phool/recoConsts.h>
0046
0047 #include <simqa_modules/QAG4SimulationTracking.h>
0048 #include <qautils/QAHistManagerDef.h>
0049 #include <siliconseedsana/SiliconSeedsAna.h>
0050 #pragma GCC diagnostic push
0051
0052 #pragma GCC diagnostic ignored "-Wundefined-internal"
0053
0054
0055 #pragma GCC diagnostic pop
0056
0057 R__LOAD_LIBRARY(libfun4all.so)
0058 R__LOAD_LIBRARY(libffamodules.so)
0059 R__LOAD_LIBRARY(libsimqa_modules.so)
0060 R__LOAD_LIBRARY(libsiliconseedsana.so)
0061 R__LOAD_LIBRARY(libtrack_reco.so)
0062
0063 void ensure_dir(const std::string& path)
0064 {
0065   struct stat info;
0066
0067   if (stat(path.c_str(), &info) != 0)
0068   {
0069     std::cout << "Directory " << path << " does not exist. Creating..." << std::endl;
0070     if (mkdir(path.c_str(), 0777) != 0)
0071     {
0072       std::cerr << "Failed to create directory " << path << std::endl;
0073       exit(1);
0074     }
0075   }
0076   else if (!(info.st_mode & S_IFDIR))
0077   {
0078     std::cerr << "Path " << path << " exists but is not a directory!" << std::endl;
0079     exit(1);
0080   }
0081 };
0082
0083
0084 bool useTopologicalCluster = false;
0085 bool jpsiTodielectronOnly = true;
0086 int Fun4All_singleParticle_Silicon(std::string processID = "0")
0087 {
0088   const int nEvents = 10;
0089
0090   std::string particle_name = "J/psi";
0091 //  particle_name = "eta";
0092   std::string particle_name_tag = (particle_name == "J/psi") ? "jpsi" : particle_name;
0093   std::ostringstream pid;
0094   pid << std::setw(6) << std::setfill('0') << std::stoi(processID);
0095   std::string pid_str = pid.str();
0096
0097   char cwd[PATH_MAX];
0098   if (getcwd(cwd, sizeof(cwd)) == nullptr)
0099   {
0100     std::cerr << "Failed to get current working directory!" << std::endl;
0101     return 1;
0102   }
0103
0104 //  std::string outDir = "/sphenix/user/jaein213/tracking/SiliconSeeding/MC/macro/DST/" + particle_name_tag;
0105 //  std::string outDir2 = "/sphenix/user/jaein213/tracking/SiliconSeeding/MC/macro/ana/" + particle_name_tag;
0106   std::string baseDir(cwd);
0107   std::string outDir  = baseDir + "/DST_" + particle_name_tag;
0108   std::string outDir2 = baseDir + "/ana_" + particle_name_tag;
0109   std::string outDir3 = baseDir + "/jobtime_" + particle_name_tag;
0110   ensure_dir(outDir);
0111   ensure_dir(outDir2);
0112   ensure_dir(outDir3);
0113   ensure_dir(outDir + "/qa");
0114   Fun4AllServer *se = Fun4AllServer::instance();
0115   se->Verbosity(10);
0116
0117   // Opt to print all random seed used for debugging reproducibility. Comment out to reduce stdout prints.
0118   PHRandomSeed::Verbosity(1);
0119
0120   recoConsts *rc = recoConsts::instance();
0121   Input::VERBOSITY = 0;
0122   Input::SIMPLE = true;
0123   if(jpsiTodielectronOnly)
0124   {
0125     EVTGENDECAYER::DecayFile = "decayfile/JpsiDielectron.DEC";
0126   }
0127   InputInit();
0128   if (Input::SIMPLE)
0129   {
0130     INPUTGENERATOR::SimpleEventGenerator[0]->add_particles(particle_name, 1);
0131     INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_function(PHG4SimpleEventGenerator::Gaus,
0132         PHG4SimpleEventGenerator::Gaus,
0133         PHG4SimpleEventGenerator::Gaus);
0134     // INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_mean(0., 0., 0.);
0135     INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_width(0.01, 0.01, 5.);
0136     INPUTGENERATOR::SimpleEventGenerator[0]->set_vertex_distribution_width(0, 0, 0);
0137     INPUTGENERATOR::SimpleEventGenerator[0]->set_eta_range(-1.1, 1.1);
0138     INPUTGENERATOR::SimpleEventGenerator[0]->set_phi_range(-M_PI, M_PI);
0139     // INPUTGENERATOR::SimpleEventGenerator[0]->set_pt_range(0.1, 20.);
0140     INPUTGENERATOR::SimpleEventGenerator[0]->set_pt_range(0.3, 20.);
0141   }
0142
0143   InputRegister();
0144
0145   if (!Input::READHITS)
0146   {
0147     rc->set_IntFlag("RUNNUMBER", 1);
0148
0149     SyncReco *sync = new SyncReco();
0150     se->registerSubsystem(sync);
0151
0152     HeadReco *head = new HeadReco();
0153     se->registerSubsystem(head);
0154   }
0155   FlagHandler *flag = new FlagHandler();
0156   se->registerSubsystem(flag);
0157   TimerStats *ts = new TimerStats();
0158   std::string jobtimeFile = outDir3 + "/jobtime_" + processID + ".root";
0159   ts->OutFileName(jobtimeFile);
0160   se->registerSubsystem(ts);
0161   // Simulation setup
0162   Enable::MBDFAKE = true;
0163   Enable::PIPE = true;
0164   Enable::PIPE_ABSORBER = true;
0165   Enable::MVTX = true;
0166   Enable::INTT = true;
0167   Enable::TPC = true;
0168   Enable::MICROMEGAS = true;
0169
0170   Enable::MAGNET = true;
0171   Enable::MAGNET_ABSORBER = true;
0172
0173   Enable::PLUGDOOR_ABSORBER = true;
0174
0175   Enable::CEMC = true;
0176   Enable::CEMC_ABSORBER = true;
0177   Enable::CEMC_CELL = true;
0178   Enable::CEMC_TOWER = true;
0179   Enable::CEMC_CLUSTER = true;
0180   // Enable::CEMC_EVAL =  false;
0181   // Enable::CEMC_QA =  false;
0182
0183   Enable::HCALIN = true;
0184   Enable::HCALIN_ABSORBER = true;
0185   Enable::HCALIN_CELL =  true;
0186   Enable::HCALIN_TOWER =  true;
0187   // Enable::HCALIN_CLUSTER =  true;
0188   // Enable::HCALIN_EVAL =  true;
0189   // Enable::HCALIN_QA =  true;
0190
0191   Enable::HCALOUT = true;
0192   Enable::HCALOUT_ABSORBER = true;
0193   Enable::HCALOUT_CELL =  true;
0194   Enable::HCALOUT_TOWER =  true;
0195   // Enable::HCALOUT_CLUSTER = false;
0196   // Enable::HCALOUT_EVAL =  false;
0197   // Enable::HCALOUT_QA =  false;
0198
0199   // Tracking setup
0200
0201   Enable::CDB = true;
0202   rc->set_StringFlag("CDB_GLOBALTAG", CDB::global_tag);
0203   rc->set_uint64Flag("TIMESTAMP", CDB::timestamp);
0204
0205   MagnetInit();
0206
0207   G4Init();
0208   G4Setup();
0209
0210   Mbd_Reco();
0211   Mvtx_Cells();
0212   Intt_Cells();
0213   //TPC_Cells();
0214   //Micromegas_Cells();
0215
0216   // TrackingInit();
0217   ACTSGEOM::ActsGeomInit();
0218   Mvtx_Clustering();
0219   Intt_Clustering();
0220   // TPC_Clustering();
0221   // Micromegas_Clustering();
0222
0223   // // Calo setup
0224   CEMC_Cells();
0225   CEMC_Towers();
0226   CEMC_Clusters();
0227   TopoClusterReco();
0228   HCALInner_Cells();
0229   HCALOuter_Cells();
0230   HCALInner_Towers();
0231   // HCALInner_Clusters();
0232   HCALOuter_Towers();
0233   //  HCALOuter_Clusters();
0234
0235   auto silicon_Seeding = new PHActsSiliconSeeding;
0236   // silicon_Seeding->set_track_map_name("SvtxTrackSeedContainer");
0237   silicon_Seeding->Verbosity(0);
0238   silicon_Seeding->setinttRPhiSearchWindow(0.2);
0239   silicon_Seeding->setinttZSearchWindow(1.0);
0240   silicon_Seeding->seedAnalysis(false);
0241   se->registerSubsystem(silicon_Seeding);
0242
0243   auto merger = new PHSiliconSeedMerger;
0244   // merger->trackMapName("SvtxTrackSeedContainer");
0245   merger->Verbosity(0);
0246   se->registerSubsystem(merger);
0247
0248   auto converter = new TrackSeedTrackMapConverter;
0249   // SiliconTrackSeedContainer or TpcTrackSeedContainer
0250   converter->setTrackSeedName("SiliconTrackSeedContainer");
0251   converter->setFieldMap(G4MAGNET::magfield_tracking);
0252   converter->Verbosity(10);
0253   se->registerSubsystem(converter);
0254   PHSimpleVertexFinder *finder = new PHSimpleVertexFinder;
0255   finder->Verbosity(0);
0256   finder->setDcaCut(0.1);
0257   finder->setTrackPtCut(0.);
0258   finder->setBeamLineCut(1);
0259   finder->setTrackQualityCut(20);
0260   finder->setNmvtxRequired(3);
0261   finder->setOutlierPairCut(0.1);
0262   se->registerSubsystem(finder);
0263
0264   Global_Reco(); // Definition is below..
0265   /////////////////////
0266   // void Global_Reco()
0267   // {
0268   //   int verbosity = std::max(Enable::VERBOSITY, Enable::GLOBAL_VERBOSITY);
0269
0270   //   Fun4AllServer* se = Fun4AllServer::instance();
0271
0272   //   GlobalVertexReco* gblvertex = new GlobalVertexReco();
0273   //   gblvertex->Verbosity(verbosity);
0274   //   se->registerSubsystem(gblvertex);
0275
0276   //   return;
0277   // }
0278
0279   build_truthreco_tables(); // Definition is below..
0280   //////////////////////////////////
0281   // void build_truthreco_tables()
0282   // {
0283   //   int verbosity = std::max(Enable::VERBOSITY, Enable::TRACKING_VERBOSITY);
0284   //   Fun4AllServer* se = Fun4AllServer::instance();
0285
0286   //   // this module builds high level truth track association table.
0287   //   // If this module is used, this table should be called before any evaluator calls.
0288   //   // Removing this module, evaluation will still work but trace truth association through the layers of G4-hit-cluster
0289   //   SvtxTruthRecoTableEval *tables = new SvtxTruthRecoTableEval();
0290   //   tables->Verbosity(verbosity);
0291   //   se->registerSubsystem(tables);
0292
0293   //   return;
0294   // }
0295   /////////////////////////////////
0296   bool runTruth = false;
0297   auto ensure_dir = [](const std::string &path)
0298   {
0299     struct stat info;
0300     if (stat(path.c_str(), &info) != 0)
0301     {
0302       std::cout << "Directory " << path << " does not exist. Creating..." << std::endl;
0303       mkdir(path.c_str(), 0777);
0304     }
0305     else if (!(info.st_mode & S_IFDIR))
0306     {
0307       std::cerr << "Path " << path << " exists but is not a directory!" << std::endl;
0308       exit(1);
0309     }
0310   };
0311
0312
0313   std::string outputName = outDir + "/DST_SiliconOnly_single_" + particle_name_tag + "_vtxfixed_Caloevents_";
0314   if (runTruth)
0315     outputName += "truth";
0316   else
0317     outputName += "reconstructed";
0318   outputName += "Info_" + processID + ".root";
0319   std::string outputName2 = outDir2 + "/ana_" + processID + ".root";
0320
0321   Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", outputName);
0322   se->registerOutputManager(out);
0323
0324   auto projection = new PHActsTrackProjection("CaloProjection");
0325   se->registerSubsystem(projection);
0326   QAG4SimulationTracking *qa = new QAG4SimulationTracking();
0327   qa->Verbosity(5);
0328   se->registerSubsystem(qa);
0329   auto siana = new SiliconSeedsAna("SiliconSeedsAna");
0330   siana->setMC(true);
0331   siana->setVtxSkip(true);
0332   siana->setTopoCluster(useTopologicalCluster);
0333   siana->setOutputFileName(outputName2);
0334   int startnumber = nEvents * std::stoi(processID);
0335   siana->setStartEventNumber(startnumber);
0336   se->registerSubsystem(siana);
0337   se->run(nEvents);
0338
0339   std::string qaFile = outDir + "/qa/single_" + particle_name_tag + "_QA_VTXfixed" + processID + ".root";
0340   QAHistManagerDef::saveQARootFile(qaFile);
0341
0342   //-----
0343   // Exit
0344   //-----
0345
0346   //  CDBInterface::instance()->Print(); // print used DB files
0347   se->End();
0348   std::cout << "All done" << std::endl;
0349   delete se;
0350   if (Enable::PRODUCTION)
0351   {
0352     Production_MoveOutput();
0353   }
0354
0355   gSystem->Exit(0);
0356   return 0;
0357 }
0358 #endif