SiCalo/PHYTIAmacro/Fun4All_PHYTIA_Silicon.C

0001 #ifndef MACRO_FUN4ALLG4SPHENIX_C
0002 #define MACRO_FUN4ALLG4SPHENIX_C
0003
0004 #include <fun4all/Fun4AllUtils.h>
0005 #include <G4_ActsGeom.C>
0006 #include <G4_Global.C>
0007 #include <G4_Magnet.C>
0008 #include <G4_Mbd.C>
0009 #include <GlobalVariables.C>
0010 #include <QA.C>
0011 #include <G4_TopoClusterReco.C>
0012 #include <Trkr_Clustering.C>
0013 #include <Trkr_Reco.C>
0014 #include <Trkr_RecoInit.C>
0015 #include <Trkr_TpcReadoutInit.C>
0016
0017 #include <ffamodules/FlagHandler.h>
0018 #include <ffamodules/HeadReco.h>
0019 #include <ffamodules/SyncReco.h>
0020 #include <ffamodules/CDBInterface.h>
0021
0022
0023 #include <fun4all/Fun4AllDstInputManager.h>
0024 #include <fun4all/Fun4AllDstOutputManager.h>
0025 #include <fun4all/Fun4AllInputManager.h>
0026 #include <fun4all/Fun4AllOutputManager.h>
0027 #include <fun4all/Fun4AllRunNodeInputManager.h>
0028 #include <fun4all/Fun4AllServer.h>
0029
0030 #include <phool/recoConsts.h>
0031
0032 #include <cdbobjects/CDBTTree.h>
0033
0034 #include <tpccalib/PHTpcResiduals.h>
0035
0036 #include <trackingqa/InttClusterQA.h>
0037 #include <trackingqa/MicromegasClusterQA.h>
0038 #include <trackingqa/MvtxClusterQA.h>
0039 #include <trackingqa/TpcClusterQA.h>
0040 #include <tpcqa/TpcRawHitQA.h>
0041 #include <trackingqa/TpcSeedsQA.h>
0042 #include <trackingqa/SiliconSeedsQA.h>
0043 #include <trackingqa/VertexQA.h>
0044 #include <trackingdiagnostics/TrackResiduals.h>
0045 #include <trackingdiagnostics/TrkrNtuplizer.h>
0046 //#include <trackingdiagnostics/KshortReconstruction.h>
0047 #include <siliconseedsana/SiliconSeedsAna.h>
0048
0049 #include <trackreco/PHActsTrackProjection.h>
0050
0051 // #include <caloreco/CaloGeomMapping.h>
0052 // #include <caloreco/CaloGeomMappingv2.h>
0053 // #include <caloreco/RawClusterBuilderTemplate.h>
0054 // #include <caloreco/RawClusterBuilderTopo.h>
0055
0056 #include <simqa_modules/QAG4SimulationTracking.h>
0057 #include <stdio.h>
0058 #include <sys/stat.h>   // stat, mkdir
0059 #include <limits.h>     // PATH_MAX
0060
0061 #pragma GCC diagnostic push
0062
0063 #pragma GCC diagnostic ignored "-Wundefined-internal"
0064
0065 #include <kfparticle_sphenix/KFParticle_sPHENIX.h>
0066 //#include <kfparticle_sphenix/KshortReconstruction_local.h>
0067
0068 #pragma GCC diagnostic pop
0069
0070 R__LOAD_LIBRARY(libkfparticle_sphenix.so)
0071 R__LOAD_LIBRARY(libfun4all.so)
0072 R__LOAD_LIBRARY(libffamodules.so)
0073 R__LOAD_LIBRARY(libphool.so)
0074 R__LOAD_LIBRARY(libcdbobjects.so)
0075 R__LOAD_LIBRARY(libmvtx.so)
0076 R__LOAD_LIBRARY(libintt.so)
0077 R__LOAD_LIBRARY(libtrack_reco.so)
0078   //R__LOAD_LIBRARY(libcalo_reco.so)
0079 R__LOAD_LIBRARY(libtpc.so)
0080 R__LOAD_LIBRARY(libmicromegas.so)
0081 R__LOAD_LIBRARY(libTrackingDiagnostics.so)
0082 R__LOAD_LIBRARY(libtrackingqa.so)
0083 R__LOAD_LIBRARY(libtpcqa.so)
0084 R__LOAD_LIBRARY(libsiliconseedsana.so)
0085 // For HepMC Hijing
0086 // try inputFile = /sphenix/sim/sim01/sphnxpro/sHijing_HepMC/sHijing_0-12fm.dat
0087 void ensure_dir(const std::string& path)
0088 {
0089   struct stat info;
0090
0091   if (stat(path.c_str(), &info) != 0)
0092   {
0093     std::cout << "Directory " << path << " does not exist. Creating..." << std::endl;
0094     if (mkdir(path.c_str(), 0777) != 0)
0095     {
0096       std::cerr << "Failed to create directory " << path << std::endl;
0097       exit(1);
0098     }
0099   }
0100   else if (!(info.st_mode & S_IFDIR))
0101   {
0102     std::cerr << "Path " << path << " exists but is not a directory!" << std::endl;
0103     exit(1);
0104   }
0105 };
0106
0107 const std::string trkrdir = "/sphenix/lustre01/sphnxpro/mdc2/js_pp200_signal/nopileup/trkrcluster/run0022/detroit/";
0108 const std::string calodir = "/sphenix/lustre01/sphnxpro/mdc2/js_pp200_signal/nopileup/calocluster/run0022/detroit/";
0109 const std::string truthdir = "/sphenix/lustre01/sphnxpro/mdc2/js_pp200_signal/nopileup/trkrhit/run0022/detroit/";
0110 const std::string trkrclusterfilename = "DST_TRKR_CLUSTER_pythia8_Detroit-0000000022-000000.root";
0111 const std::string calofilename = "DST_CALO_CLUSTER_pythia8_Detroit-0000000022-000000.root";
0112 const std::string truthfilename = "DST_TRUTH_pythia8_Detroit-0000000022-00000.root";
0113 int Fun4All_PHYTIA_Silicon(std::string processID = "0")
0114 {
0115   const int nEvents = 10;
0116
0117   std::ostringstream pid;
0118   pid << std::setw(6) << std::setfill('0') << std::stoi(processID);
0119   std::string pid_str = pid.str();
0120
0121   char cwd[PATH_MAX];
0122   if (getcwd(cwd, sizeof(cwd)) == nullptr)
0123   {
0124     std::cerr << "Failed to get current working directory!" << std::endl;
0125     return 1;
0126   }
0127
0128   std::string baseDir(cwd);
0129   std::string outDir  = baseDir + "/DST";
0130   std::string outDir2 = baseDir + "/ana";
0131 //  std::string outDir = "/sphenix/user/jaein213/tracking/SiliconSeeding/MC/PYTHIA/macro/DST";
0132 //  std::string outDir2 = "/sphenix/user/jaein213/tracking/SiliconSeeding/MC/PYTHIA/macro/ana";
0133   ensure_dir(outDir);
0134   ensure_dir(outDir + "/qa");
0135   ensure_dir(outDir2);
0136
0137   bool useTopologicalCluster = false;
0138
0139   std::string trkrclusterfilename = "DST_TRKR_CLUSTER_pythia8_Detroit-0000000022-" + pid_str + ".root";
0140   std::string calofilename = "DST_CALO_CLUSTER_pythia8_Detroit-0000000022-" + pid_str + ".root";
0141   std::string truthfilename = "DST_TRUTH_pythia8_Detroit-0000000022-" + pid_str + ".root";
0142   std::string inputclusterFile = trkrdir + trkrclusterfilename;
0143   std::string inputcaloclusterFile = calodir + calofilename;
0144   std::string inputtruthFile = truthdir + truthfilename;
0145
0146   // Enable::MVTX_APPLYMISALIGNMENT = true;
0147   // ACTSGEOM::mvtx_applymisalignment = Enable::MVTX_APPLYMISALIGNMENT;
0148   // ACTSGEOM::mvtxMisalignment = 100;
0149   // ACTSGEOM::inttMisalignment = 100.;
0150   // ACTSGEOM::tpotMisalignment = 100.;
0151   auto se = Fun4AllServer::instance();
0152   se->Verbosity(2);
0153   auto rc = recoConsts::instance();
0154
0155   Enable::CDB = true;
0156   rc->set_StringFlag("CDB_GLOBALTAG", "MDC2");
0157   rc->set_uint64Flag("TIMESTAMP", 6);
0158   rc->set_IntFlag("RUNNUMBER", 1);
0159
0160   SyncReco *sync = new SyncReco();
0161   se->registerSubsystem(sync);
0162
0163   HeadReco *head = new HeadReco();
0164   se->registerSubsystem(head);
0165   FlagHandler *flag = new FlagHandler();
0166   se->registerSubsystem(flag);
0167   // std::string geofile = CDBInterface::instance()->getUrl("Tracking_Geometry");
0168
0169   // Fun4AllRunNodeInputManager *ingeo = new Fun4AllRunNodeInputManager("GeoIn");
0170   // ingeo->AddFile(geofile);
0171   // se->registerInputManager(ingeo);
0172
0173   G4MAGNET::magfield_rescale = 1;
0174   TrackingInit();
0175
0176   auto hitsin = new Fun4AllDstInputManager("InputManager");
0177   hitsin->fileopen(inputclusterFile);
0178   se->registerInputManager(hitsin);
0179   auto hitsin_cluster = new Fun4AllDstInputManager("DSTin_cluster");
0180   hitsin_cluster->fileopen(inputcaloclusterFile);
0181   se->registerInputManager(hitsin_cluster);
0182   auto hitsin_truth = new Fun4AllDstInputManager("DSTin_truth");
0183   hitsin_truth->fileopen(inputtruthFile);
0184   se->registerInputManager(hitsin_truth);
0185
0186   // RawClusterBuilderTemplate *ClusterBuilder = new RawClusterBuilderTemplate("EmcRawClusterBuilderTemplate");
0187
0188   // ClusterBuilder->Detector("CEMC");
0189   // ClusterBuilder->set_threshold_energy(0.030); // This threshold should be the same as in CEMCprof_Thresh**.root file below
0190   // std::string emc_prof = getenv("CALIBRATIONROOT");
0191   // emc_prof += "/EmcProfile/CEMCprof_Thresh30MeV.root";
0192   // ClusterBuilder->LoadProfile(emc_prof);
0193   // ClusterBuilder->set_UseTowerInfo(1); // just use towerinfo
0194   // //    ClusterBuilder->set_UseTowerInfo(1); // to use towerinfo objects rather than old RawTower
0195   // se->registerSubsystem(ClusterBuilder);
0196   TopoClusterReco();
0197
0198   auto silicon_Seeding = new PHActsSiliconSeeding;
0199   // silicon_Seeding->set_track_map_name("SvtxTrackSeedContainer");
0200   silicon_Seeding->Verbosity(0);
0201   // silicon_Seeding->searchInIntt();
0202   // silicon_Seeding->setinttRPhiSearchWindow(0.4);
0203   // silicon_Seeding->setinttZSearchWindow(1.6);
0204   silicon_Seeding->setinttRPhiSearchWindow(0.2);
0205   silicon_Seeding->setinttZSearchWindow(1.0);
0206   silicon_Seeding->seedAnalysis(false);
0207   se->registerSubsystem(silicon_Seeding);
0208
0209   auto merger = new PHSiliconSeedMerger;
0210   // merger->trackMapName("SvtxTrackSeedContainer");
0211   merger->Verbosity(0);
0212   se->registerSubsystem(merger);
0213
0214   // Default set to full SvtxTrackSeeds. Can be set to
0215   // SiliconTrackSeedContainer or TpcTrackSeedContainer
0216   auto converter = new TrackSeedTrackMapConverter("SiliconSeedConverter");
0217   converter->setTrackSeedName("SiliconTrackSeedContainer");
0218   // converter->setTrackMapName("SiliconSvtxTrackMap");
0219
0220   converter->setFieldMap(G4MAGNET::magfield_tracking);
0221   converter->Verbosity(0);
0222   se->registerSubsystem(converter);
0223
0224   auto finder = new PHSimpleVertexFinder;
0225   finder->Verbosity(0);
0226   finder->setDcaCut(0.075);
0227   finder->setTrackPtCut(0.1);
0228   finder->setBeamLineCut(1);
0229   finder->setTrackQualityCut(300);
0230   finder->setNmvtxRequired(3);
0231   finder->setOutlierPairCut(0.075);
0232
0233   // finder->setTrackMapName("SiliconSvtxTrackMap");
0234   // finder->setVertexMapName("SiliconSvtxVertexMap");
0235   se->registerSubsystem(finder);
0236   // Tracking_Reco(); // Function used to setup the tracking reco. comment out here since we have setup already above
0237
0238   Global_Reco(); // Definition is below..
0239                  /////////////////////
0240                  // void Global_Reco()
0241                  // {
0242                  //   int verbosity = std::max(Enable::VERBOSITY, Enable::GLOBAL_VERBOSITY);
0243
0244   //   Fun4AllServer* se = Fun4AllServer::instance();
0245
0246   //   GlobalVertexReco* gblvertex = new GlobalVertexReco();
0247   //   gblvertex->Verbosity(verbosity);
0248   //   se->registerSubsystem(gblvertex);
0249
0250   //   return;
0251   // }
0252
0253   //  build_truthreco_tables(); // Definition is below..
0254   //////////////////////////////////
0255   // void build_truthreco_tables()
0256   // {
0257   //   int verbosity = std::max(Enable::VERBOSITY, Enable::TRACKING_VERBOSITY);
0258   //   Fun4AllServer* se = Fun4AllServer::instance();
0259
0260   //   // this module builds high level truth track association table.
0261   //   // If this module is used, this table should be called before any evaluator calls.
0262   //   // Removing this module, evaluation will still work but trace truth association through the layers of G4-hit-cluster
0263   //   SvtxTruthRecoTableEval *tables = new SvtxTruthRecoTableEval();
0264   //   tables->Verbosity(verbosity);
0265   //   se->registerSubsystem(tables);
0266
0267   //   return;
0268   // }
0269   /////////////////////////////////
0270   bool runTruth = false;
0271   // std::string outDir = "/sphenix/user/jaein213/tracking/SiliconSeeding/MC/macro/DST";
0272
0273   auto projection = new PHActsTrackProjection("CaloProjection");
0274   se->registerSubsystem(projection);
0275
0276   // Ensure output directories exist
0277
0278
0279
0280   std::string outputName = outDir + "/DST_SiliconOnly_PHYTIA_";
0281   outputName += "Info_" + processID + ".root";
0282   std::string outputName2 = outDir2 + "/ana_" + processID + ".root";
0283   // SecondaryVertexTagger* myTagger = new SecondaryVertexTagger();
0284   // myTagger->Verbosity(2);
0285   // if (runTruth)
0286   //{
0287   //   myTagger->truthMode();
0288   // }
0289   // else
0290   //{
0291   //   myTagger->truthRecoMatch();
0292   // }
0293   // myTagger->setOutputName(outputName.c_str());
0294   // myTagger->saveOutput();
0295   // se->registerSubsystem(myTagger);
0296   // outputName = processID + "test.root";
0297   Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", outputName);
0298   se->registerOutputManager(out);
0299   QAG4SimulationTracking *qa = new QAG4SimulationTracking();
0300   qa->Verbosity(0);
0301 //  se->registerSubsystem(qa);
0302   // outputName2 = processID + "ana.root";
0303   auto siana = new SiliconSeedsAna("SiliconSeedsAna");
0304   siana->setMC(true);
0305   siana->setTrackMapName("SvtxTrackMap");
0306   siana->setVertexMapName("SvtxVertexMap");
0307   siana->setTopoCluster(useTopologicalCluster);
0308   // siana->setVtxSkip(true);
0309   siana->setOutputFileName(outputName2);
0310   int startnumber = nEvents * std::stoi(processID);
0311   siana->setStartEventNumber(startnumber);
0312   siana->setEMCalClusterContainerName("CLUSTERINFO_CEMC");
0313   se->registerSubsystem(siana);
0314 //  se->skip(900);
0315   se->run(nEvents);
0316
0317   std::string qaFile = outDir + "/qa/singleElectron_QA_VTXfixed" + processID + ".root";
0318   QAHistManagerDef::saveQARootFile(qaFile);
0319
0320   //-----
0321   // Exit
0322   //-----
0323
0324   //  CDBInterface::instance()->Print(); // print used DB files
0325   se->End();
0326   std::cout << "All done" << std::endl;
0327   delete se;
0328
0329   gSystem->Exit(0);
0330   return 0;
0331 }
0332 #endif