tpc/TpcDVCalib/Fun4All_TrackAnalysis.C

0001 /*
0002  * Track + Calo matching for TPC Drift Velocity Calibration
0003  * using the production track cluster and track seed DSTs
0004  */
0005
0006 #include <GlobalVariables.C>
0007
0008 #include <G4_ActsGeom.C>
0009 #include <G4_Magnet.C>
0010 #include <Trkr_Clustering.C>
0011 #include <Trkr_RecoInit.C>
0012 #include <Trkr_TpcReadoutInit.C>
0013 #include <Trkr_Reco.C>
0014 #include <G4_Global.C>
0015
0016 #include <caloreco/RawClusterBuilderTopo.h>
0017
0018
0019 #include <tpcdvcalib/TrackToCalo.h>
0020
0021 #include <cdbobjects/CDBTTree.h>
0022
0023 #include <tpccalib/PHTpcResiduals.h>
0024
0025 #include <trackingdiagnostics/TrackResiduals.h>
0026 #include <trackingdiagnostics/TrkrNtuplizer.h>
0027
0028 #include <trackreco/AzimuthalSeeder.h>
0029 #include <trackreco/PHActsTrackProjection.h>
0030 #include <trackbase_historic/SvtxTrack.h>
0031
0032 #include <ffamodules/CDBInterface.h>
0033
0034 #include <fun4all/Fun4AllDstInputManager.h>
0035 #include <fun4all/Fun4AllDstOutputManager.h>
0036 #include <fun4all/Fun4AllInputManager.h>
0037 #include <fun4all/Fun4AllUtils.h>
0038 #include <fun4all/Fun4AllOutputManager.h>
0039 #include <fun4all/Fun4AllRunNodeInputManager.h>
0040 #include <fun4all/Fun4AllServer.h>
0041
0042 #include <phool/recoConsts.h>
0043
0044 #include <iostream>
0045 #include <filesystem>
0046
0047 R__LOAD_LIBRARY(libfun4all.so)
0048 R__LOAD_LIBRARY(libffamodules.so)
0049 R__LOAD_LIBRARY(libphool.so)
0050 R__LOAD_LIBRARY(libcdbobjects.so)
0051 R__LOAD_LIBRARY(libg4dst.so)
0052 R__LOAD_LIBRARY(libmvtx.so)
0053 R__LOAD_LIBRARY(libintt.so)
0054 R__LOAD_LIBRARY(libtpc.so)
0055 R__LOAD_LIBRARY(libmicromegas.so)
0056 R__LOAD_LIBRARY(libTrackingDiagnostics.so)
0057 R__LOAD_LIBRARY(libtrack_reco.so)
0058 R__LOAD_LIBRARY(libTpcDVCalib.so)
0059 R__LOAD_LIBRARY(libcalo_reco.so)
0060
0061 std::string GetFirstLine(const std::string& listname);
0062 bool is_directory_empty(const std::filesystem::path& dir_path);
0063
0064 void Fun4All_TrackAnalysis(
0065     const int nEvents = 10,
0066     std::vector<std::string> myInputLists = {
0067         "run46730_0000_trkr_seed.txt",
0068         "run46730_0000_trkr_cluster.txt",
0069         "run46730_calo.list"},
0070     const std::string& outDir = "./",
0071     bool doTpcOnlyTracking = true,
0072     float initial_driftvelocity = 0.00710,
0073     bool doEMcalRadiusCorr = true,
0074     const bool convertSeeds = false)
0075 {
0076   int verbosity = 0;
0077
0078   G4TRACKING::convert_seeds_to_svtxtracks = convertSeeds;
0079   std::cout << "Converting to seeds : " << G4TRACKING::convert_seeds_to_svtxtracks << std::endl;
0080   std::string firstfile = GetFirstLine(myInputLists[0]);
0081   if (*firstfile.c_str() == '\0') return;
0082   std::pair<int, int> runseg = Fun4AllUtils::GetRunSegment(firstfile);
0083   int runnumber = runseg.first;
0084   int segment = runseg.second;
0085
0086   TpcReadoutInit(runnumber);
0087   std::cout << " run: " << runnumber
0088             << " samples: " << TRACKING::reco_tpc_maxtime_sample
0089             << " pre: " << TRACKING::reco_tpc_time_presample
0090             << " vdrift: " << G4TPC::tpc_drift_velocity_reco
0091             << std::endl;
0092
0093   // distortion calibration mode
0094   /*
0095    * set to true to enable residuals in the TPC with
0096    * TPC clusters not participating to the ACTS track fit
0097    */
0098   G4TRACKING::SC_CALIBMODE = false;
0099
0100   ACTSGEOM::mvtxMisalignment = 100;
0101   ACTSGEOM::inttMisalignment = 100.;
0102   ACTSGEOM::tpotMisalignment = 100.;
0103
0104   auto *se = Fun4AllServer::instance();
0105   se->Verbosity(1);
0106   auto *rc = recoConsts::instance();
0107   rc->set_IntFlag("RUNNUMBER", runnumber);
0108
0109   Enable::CDB = true;
0110   rc->set_StringFlag("CDB_GLOBALTAG", "ProdA_2024");
0111   rc->set_uint64Flag("TIMESTAMP", runnumber);
0112   std::string geofile = CDBInterface::instance()->getUrl("Tracking_Geometry");
0113
0114   Fun4AllRunNodeInputManager *ingeo = new Fun4AllRunNodeInputManager("GeoIn");
0115   ingeo->AddFile(geofile);
0116   se->registerInputManager(ingeo);
0117
0118   G4TPC::tpc_drift_velocity_reco = initial_driftvelocity; //cm/ns
0119   G4TPC::ENABLE_MODULE_EDGE_CORRECTIONS = true;
0120   // to turn on the default static corrections, enable the two lines below
0121   G4TPC::ENABLE_STATIC_CORRECTIONS = true;
0122   // G4TPC::DISTORTIONS_USE_PHI_AS_RADIANS = false;
0123
0124   G4MAGNET::magfield_rescale = 1;
0125   TrackingInit();
0126
0127   //Add all required input files
0128   for (unsigned int i = 0; i < myInputLists.size(); ++i)
0129   {
0130     Fun4AllInputManager *infile = new Fun4AllDstInputManager("DSTin_" + std::to_string(i));
0131     std::cout << "Including file " << myInputLists[i] << std::endl;
0132     infile->AddListFile(myInputLists[i]);
0133     se->registerInputManager(infile);
0134   }
0135
0136   /*
0137    * Either converts seeds to tracks with a straight line/helix fit
0138    * or run the full Acts track kalman filter fit
0139    */
0140   if (G4TRACKING::convert_seeds_to_svtxtracks)
0141   {
0142     auto *converter = new TrackSeedTrackMapConverter;
0143     // Option to use TpcTrackSeedContainer or SvtxTrackSeeds
0144     // can be set to SiliconTrackSeedContainer for silicon-only track fit
0145     if (doTpcOnlyTracking)
0146     {
0147       converter->setTrackSeedName("TpcTrackSeedContainer");
0148     }
0149     else
0150     {
0151       converter->setTrackSeedName("SvtxTrackSeeds");
0152     }
0153     converter->setFieldMap(G4MAGNET::magfield_tracking);
0154     converter->Verbosity(0);
0155     se->registerSubsystem(converter);
0156   }
0157   else
0158   {
0159     auto *deltazcorr = new PHTpcDeltaZCorrection;
0160     deltazcorr->Verbosity(0);
0161     se->registerSubsystem(deltazcorr);
0162
0163     // perform final track fit with ACTS
0164     auto *actsFit = new PHActsTrkFitter;
0165     actsFit->Verbosity(0);
0166     actsFit->commissioning(G4TRACKING::use_alignment);
0167     // in calibration mode, fit only Silicons and Micromegas hits
0168     if (!doTpcOnlyTracking)
0169     {
0170       actsFit->fitSiliconMMs(G4TRACKING::SC_CALIBMODE);
0171       actsFit->setUseMicromegas(G4TRACKING::SC_USE_MICROMEGAS);
0172     }
0173     actsFit->set_pp_mode(TRACKING::pp_mode);
0174     actsFit->set_use_clustermover(true);  // default is true for now
0175     actsFit->useActsEvaluator(false);
0176     actsFit->useOutlierFinder(false);
0177     actsFit->setFieldMap(G4MAGNET::magfield_tracking);
0178     se->registerSubsystem(actsFit);
0179   }
0180
0181   PHSimpleVertexFinder *finder = new PHSimpleVertexFinder;
0182   finder->Verbosity(verbosity);
0183   finder->setDcaCut(0.5);
0184   finder->setTrackPtCut(-99999.);
0185   finder->setBeamLineCut(1);
0186   finder->setTrackQualityCut(1000000000);
0187   if (!doTpcOnlyTracking)
0188   {
0189     finder->setRequireMVTX(true);
0190     finder->setNmvtxRequired(3);
0191   }
0192   else
0193   {
0194     finder->setRequireMVTX(false);
0195   }
0196   finder->setOutlierPairCut(0.1);
0197   se->registerSubsystem(finder);
0198
0199   Global_Reco();
0200
0201   auto *projection = new PHActsTrackProjection("CaloProjection");
0202   float new_cemc_rad = 100.70;//(1-(-0.077))*93.5 recommended cemc radius
0203   if (doEMcalRadiusCorr)
0204   {
0205     projection->setLayerRadius(SvtxTrack::CEMC, new_cemc_rad);
0206   }
0207   se->registerSubsystem(projection);
0208
0209   RawClusterBuilderTopo* ClusterBuilder1 = new RawClusterBuilderTopo("EMcalRawClusterBuilderTopo1");
0210   ClusterBuilder1->Verbosity(verbosity);
0211   ClusterBuilder1->set_nodename("TOPOCLUSTER_EMCAL");
0212   ClusterBuilder1->set_enable_HCal(false);
0213   ClusterBuilder1->set_enable_EMCal(true);
0214   //ClusterBuilder1->set_noise(0.0025, 0.006, 0.03);
0215   ClusterBuilder1->set_noise(0.01, 0.03, 0.03);
0216   ClusterBuilder1->set_significance(4.0, 2.0, 1.0);
0217   ClusterBuilder1->allow_corner_neighbor(true);
0218   ClusterBuilder1->set_do_split(true);
0219   ClusterBuilder1->set_minE_local_max(1.0, 2.0, 0.5);
0220   ClusterBuilder1->set_R_shower(0.025);
0221   se->registerSubsystem(ClusterBuilder1);
0222
0223   std::string outputAnaFileName = "TrackCalo_" + std::to_string(segment) + "_ana.root";
0224   std::string outputRecoDir = outDir + "inReconstruction/" + std::to_string(runnumber) + "/";
0225   std::string makeDirectory = "mkdir -p " + outputRecoDir;
0226   system(makeDirectory.c_str());
0227   std::string outputAnaFile = outputRecoDir + outputAnaFileName;
0228   std::cout << "Reco ANA file: " << outputAnaFile << std::endl;
0229
0230   TrackToCalo *ttc = new TrackToCalo("Tracks_And_Calo", outputAnaFile);
0231   ttc->EMcalRadiusUser(true);
0232   ttc->setEMcalRadius(new_cemc_rad);
0233   se->registerSubsystem(ttc);
0234
0235   std::ifstream file_ana(outputAnaFile.c_str(), std::ios::binary | std::ios::ate);
0236   if (file_ana.good() && (file_ana.tellg() > 100))
0237   {
0238     std::string outputRecoDirMove = outDir + "Reconstructed/" + std::to_string(runnumber) + "/";
0239     std::string makeDirectoryMove = "mkdir -p " + outputRecoDirMove;
0240     system(makeDirectoryMove.c_str());
0241     std::string moveOutput = "mv " + outputAnaFile + " " + outDir + "Reconstructed/" + std::to_string(runnumber);
0242     std::cout << "moveOutput: " << moveOutput << std::endl;
0243     system(moveOutput.c_str());
0244   }
0245
0246   se->run(nEvents);
0247   se->End();
0248   se->PrintTimer();
0249
0250   delete se;
0251   std::cout << "Finished" << std::endl;
0252   gSystem->Exit(0);
0253 }
0254
0255 std::string GetFirstLine(const std::string& listname)
0256 {
0257   std::ifstream file(listname);
0258
0259   std::string firstLine;
0260   if (file.is_open()) {
0261       if (std::getline(file, firstLine)) {
0262           std::cout << "First Line: " << firstLine << std::endl;
0263       } else {
0264           std::cerr << "Unable to read first line of file" << std::endl;
0265       }
0266       file.close();
0267   } else {
0268       std::cerr << "Unable to open file" << std::endl;
0269   }
0270   return firstLine;
0271 }
0272
0273 bool is_directory_empty(const std::filesystem::path& dir_path) {
0274     if (std::filesystem::exists(dir_path) && std::filesystem::is_directory(dir_path)) {
0275         return std::filesystem::directory_iterator(dir_path) == std::filesystem::directory_iterator();
0276     }
0277     return false;
0278 }