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