InttSeedTrackPerformance-Jingyu/macro/Fun4all_ProjectTruth.C

0001 //Truth Pro to Calo
0002 #include <G4_ActsGeom.C>
0003 #include <G4_Magnet.C>
0004 #include <GlobalVariables.C>
0005 #include <Trkr_Clustering.C>
0006 #include <Trkr_RecoInit.C>
0007 #include <Trkr_TpcReadoutInit.C>
0008 #include <Trkr_Reco.C>
0009 #include <G4_Global.C>
0010
0011 #include <caloreco/RawClusterBuilderTopo.h>
0012
0013 #include <ffamodules/CDBInterface.h>
0014 #include <fun4all/Fun4AllDstInputManager.h>
0015 #include <fun4all/Fun4AllDstOutputManager.h>
0016 #include <fun4all/Fun4AllInputManager.h>
0017 #include <fun4all/Fun4AllUtils.h>
0018 #include <fun4all/Fun4AllOutputManager.h>
0019 #include <fun4all/Fun4AllRunNodeInputManager.h>
0020 #include <fun4all/Fun4AllServer.h>
0021
0022 #include <phool/recoConsts.h>
0023
0024 #include <cdbobjects/CDBTTree.h>
0025
0026 #include <track_to_calo/TrackCaloMatch.h>
0027 #include <track_to_calo/TrackToCalo.h>
0028 #include <track_to_calo/CaloOnly.h>
0029 #include <track_to_calo/TrackOnly.h>
0030
0031 #include <trackreco/AzimuthalSeeder.h>
0032 #include <trackreco/PHActsTrackProjection.h>
0033 #include <trackingdiagnostics/TrackResiduals.h>
0034
0035 #include <trackbase_historic/SvtxTrack.h>
0036
0037 #include <stdio.h>
0038 #include <iostream>
0039 #include <filesystem>
0040
0041 #include "HFReco.C"
0042
0043 R__LOAD_LIBRARY(libfun4all.so)
0044 R__LOAD_LIBRARY(libffamodules.so)
0045 R__LOAD_LIBRARY(libphool.so)
0046 R__LOAD_LIBRARY(libcdbobjects.so)
0047 R__LOAD_LIBRARY(libmvtx.so)
0048 R__LOAD_LIBRARY(libintt.so)
0049 R__LOAD_LIBRARY(libtpc.so)
0050 R__LOAD_LIBRARY(libmicromegas.so)
0051 R__LOAD_LIBRARY(libTrackingDiagnostics.so)
0052 R__LOAD_LIBRARY(libtrack_reco.so)
0053 R__LOAD_LIBRARY(libtrack_to_calo.so)
0054 R__LOAD_LIBRARY(libcalo_reco.so)
0055 R__LOAD_LIBRARY(libkfparticle_sphenix.so)
0056
0057 using namespace std;
0058
0059 namespace fs = std::filesystem;
0060
0061 std::string GetFirstLine(std::string listname);
0062 bool is_directory_empty(const fs::path& dir_path);
0063
0064 void Fun4All_TrackAnalysis()
0065 {
0066   gSystem->Load("libg4dst.so");
0067
0068   int verbosity = 0;
0069
0070   //Create the server
0071   Fun4AllServer *se = Fun4AllServer::instance();
0072   se->Verbosity(1);
0073
0074   //Add all required input files
0075   for (unsigned int i = 0; i < myInputLists.size(); ++i)
0076   {
0077     Fun4AllInputManager *infile = new Fun4AllDstInputManager("DSTin_" + to_string(i));
0078     std::cout << "Including file " << myInputLists[i] << std::endl;
0079     //infile->AddFile(myInputLists[i]);
0080     infile->AddListFile(myInputLists[i]);
0081     se->registerInputManager(infile);
0082   }
0083
0084   Global_Reco();
0085
0086   auto projection = new PHActsTrackProjection("CaloProjection");
0087   float new_cemc_rad = 93.50;
0088   if (doEMcalRadiusCorr)
0089   {
0090     projection->setLayerRadius(SvtxTrack::CEMC, new_cemc_rad);
0091   }
0092   se->registerSubsystem(projection);
0093
0094   std::cout << "Building clusters" << std::endl;
0095   RawClusterBuilderTemplate *ClusterBuilder = new RawClusterBuilderTemplate("EmcRawClusterBuilderTemplate");
0096   ClusterBuilder->Detector("CEMC");
0097   ClusterBuilder->set_threshold_energy(0.030);  // for when using basic calibration
0098   std::string emc_prof = getenv("CALIBRATIONROOT");
0099   emc_prof += "/EmcProfile/CEMCprof_Thresh30MeV.root";
0100   ClusterBuilder->LoadProfile(emc_prof);
0101   ClusterBuilder->set_UseTowerInfo(1);  // to use towerinfo objects rather than old RawTower
0102   se->registerSubsystem(ClusterBuilder);
0103
0104 /*
0105   RawClusterBuilderTopo* ClusterBuilder1 = new RawClusterBuilderTopo("EMcalRawClusterBuilderTopo1");
0106   ClusterBuilder1->Verbosity(verbosity);
0107   ClusterBuilder1->set_nodename("TOPOCLUSTER_EMCAL");
0108   ClusterBuilder1->set_enable_HCal(false);
0109   ClusterBuilder1->set_enable_EMCal(true);
0110   //ClusterBuilder1->set_noise(0.0025, 0.006, 0.03);
0111   ClusterBuilder1->set_noise(0.01, 0.03, 0.03);
0112   ClusterBuilder1->set_significance(4.0, 2.0, 1.0);
0113   ClusterBuilder1->allow_corner_neighbor(true);
0114   ClusterBuilder1->set_do_split(true);
0115   ClusterBuilder1->set_minE_local_max(1.0, 2.0, 0.5);
0116   ClusterBuilder1->set_R_shower(0.025);
0117   se->registerSubsystem(ClusterBuilder1);
0118 */
0119
0120   //For particle flow studies
0121   RawClusterBuilderTopo* ClusterBuilder2 = new RawClusterBuilderTopo("EMcalRawClusterBuilderTopo2");
0122   ClusterBuilder2->Verbosity(verbosity);
0123   ClusterBuilder2->set_nodename("TOPOCLUSTER_HCAL");
0124   ClusterBuilder2->set_enable_HCal(true);
0125   ClusterBuilder2->set_enable_EMCal(false);
0126   //ClusterBuilder2->set_noise(0.0025, 0.006, 0.03);
0127   ClusterBuilder2->set_noise(0.01, 0.03, 0.03);
0128   ClusterBuilder2->set_significance(4.0, 2.0, 1.0);
0129   ClusterBuilder2->allow_corner_neighbor(true);
0130   ClusterBuilder2->set_do_split(true);
0131   ClusterBuilder2->set_minE_local_max(1.0, 2.0, 0.5);
0132   ClusterBuilder2->set_R_shower(0.025);
0133   se->registerSubsystem(ClusterBuilder2);
0134
0135   //CaloOnly *co = new CaloOnly("CaloOnly", outputRecoFile);
0136   //se->registerSubsystem(co);
0137
0138   //TrackOnly *to = new TrackOnly("TracksOnly", outputRecoFile);
0139   //se->registerSubsystem(to);
0140
0141   TrackCaloMatch *tcm = new TrackCaloMatch("Tracks_Calo_Match");
0142   tcm->SetMyTrackMapName("MySvtxTrackMap");
0143   tcm->writeEventDisplays(doEvtDisplay);
0144   tcm->setEventDisplayPath(outputJsonDir);
0145   tcm->setRunDate("2024-06-25");
0146   tcm->EMcalRadiusUser(doEMcalRadiusCorr);
0147   tcm->setEMcalRadius(new_cemc_rad);
0148   tcm->setdphicut(0.2);
0149   tcm->setdzcut(10000);
0150   tcm->setTrackPtLowCut(1.0);
0151   tcm->setEmcalELowCut(0.5);
0152   tcm->setnTpcClusters(30);
0153   tcm->setTrackQuality(100);
0154   tcm->setRawClusContEMName("CLUSTERINFO_CEMC");
0155   ttc->setDFNodeName("myFinder");
0156   //ttc->doTruthMatching();
0157   se->registerSubsystem(tcm);
0158
0159   // begin KFParticle
0160   PhotonConvKFPReco();
0161
0162   TrackToCalo *ttc = new TrackToCalo("Tracks_And_Calo", outputAnaFile);
0163   ttc->EMcalRadiusUser(doEMcalRadiusCorr);
0164   ttc->setEMcalRadius(new_cemc_rad);
0165   ttc->setKFPtrackMapName("PhotonConv_SvtxTrackMap");
0166   ttc->setKFPContName("PhotonConv_KFParticle_Container");
0167   //ttc->doTrkrCaloMatching();
0168   //ttc->anaTrkrInfo();
0169   //ttc->anaCaloInfo();
0170   ttc->setTrackPtLowCut(1.0);
0171   ttc->setEmcalELowCut(0.5);
0172   ttc->setnTpcClusters(30);
0173   ttc->setTrackQuality(100);
0174   ttc->doTrkrCaloMatching_KFP();
0175   ttc->setRawClusContEMName("CLUSTERINFO_CEMC");
0176   se->registerSubsystem(ttc);
0177
0178   if (Enable::DSTOUT)
0179   {
0180     Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", outputDstFile);
0181     //out->StripNode("RUN");
0182     //out->AddNode("Sync");
0183     out->SaveRunNode(0);
0184     se->registerOutputManager(out);
0185   }
0186
0187   se->run(nEvents);
0188   se->End();
0189   se->PrintTimer();
0190
0191   ifstream file_reco(outputRecoFile.c_str(), ios::binary | ios::ate);
0192   if (file_reco.good() && (file_reco.tellg() > 100))
0193   {
0194     string outputRecoDirMove = outDir + "/Reconstructed/" + to_string(runnumber) + "/";
0195     string makeDirectoryMove = "mkdir -p " + outputRecoDirMove;
0196     system(makeDirectoryMove.c_str());
0197     string moveOutput = "mv " + outputRecoFile + " " + outDir + "/Reconstructed/" + to_string(runnumber);
0198     std::cout << "moveOutput: " << moveOutput << std::endl;
0199     system(moveOutput.c_str());
0200   }
0201
0202   ifstream file_ana(outputAnaFile.c_str(), ios::binary | ios::ate);
0203   if (file_ana.good() && (file_ana.tellg() > 100))
0204   {
0205     string outputRecoDirMove = outDir + "/Reconstructed/" + to_string(runnumber) + "/";
0206     string makeDirectoryMove = "mkdir -p " + outputRecoDirMove;
0207     system(makeDirectoryMove.c_str());
0208     string moveOutput = "mv " + outputAnaFile + " " + outDir + "/Reconstructed/" + to_string(runnumber);
0209     std::cout << "moveOutput: " << moveOutput << std::endl;
0210     system(moveOutput.c_str());
0211   }
0212
0213   ifstream file_dst(outputDstFile.c_str(), ios::binary | ios::ate);
0214   if (file_dst.good() && (file_dst.tellg() > 100))
0215   {
0216     string outputRecoDirMove = outDir + "/Reconstructed/" + to_string(runnumber) + "/";
0217     string makeDirectoryMove = "mkdir -p " + outputRecoDirMove;
0218     system(makeDirectoryMove.c_str());
0219     string moveOutput = "mv " + outputDstFile + " " + outDir + "/Reconstructed/" + to_string(runnumber);
0220     std::cout << "moveOutput: " << moveOutput << std::endl;
0221     system(moveOutput.c_str());
0222   }
0223
0224   fs::path output_json_dir_path = outputJsonDir;
0225   if (!is_directory_empty(output_json_dir_path))
0226   {
0227     string outputJsonDirMove = outDir + "/Reconstructed/" + to_string(runnumber) + "/EvtDisplay/";
0228     string makeDirectoryMove = "mkdir -p " + outputJsonDirMove;
0229     system(makeDirectoryMove.c_str());
0230     string moveOutput = "mv " + outputJsonDir + "* " + outputJsonDirMove;
0231     std::cout << "moveOutput: " << moveOutput << std::endl;
0232     system(moveOutput.c_str());
0233   }
0234
0235   delete se;
0236   std::cout << "All done" << std::endl;
0237   gSystem->Exit(0);
0238
0239   return;
0240 }
0241
0242 std::string GetFirstLine(std::string listname)
0243 {
0244   std::ifstream file(listname);
0245
0246   std::string firstLine = "";
0247   if (file.is_open()) {
0248       if (std::getline(file, firstLine)) {
0249           std::cout << "First Line: " << firstLine << std::endl;
0250       } else {
0251           std::cerr << "Unable to read first line of file" << std::endl;
0252       }
0253       file.close();
0254   } else {
0255       std::cerr << "Unable to open file" << std::endl;
0256   }
0257   return firstLine;
0258 }
0259
0260 bool is_directory_empty(const fs::path& dir_path) {
0261     if (fs::exists(dir_path) && fs::is_directory(dir_path)) {
0262         return fs::directory_iterator(dir_path) == fs::directory_iterator();
0263     }
0264     return false;
0265 }