tpc/TpcDVCalib/Fun4All_FieldOnAllTrackersCalos.C

0001 //Track + Calo matching for TPC Drift Velocity Calibration
0002 //authors: Xudong Yu <xyu3@bnl.gov>
0003
0004 #include <GlobalVariables.C>
0005
0006 #include <G4_ActsGeom.C>
0007 #include <G4_Magnet.C>
0008 #include <Trkr_Clustering.C>
0009 #include <Trkr_RecoInit.C>
0010 #include <Trkr_TpcReadoutInit.C>
0011 #include <Trkr_Reco.C>
0012 #include <G4_Global.C>
0013
0014 #include <caloreco/RawClusterBuilderTopo.h>
0015
0016 #include <tpcdvcalib/TrackToCalo.h>
0017
0018 #include <trackreco/AzimuthalSeeder.h>
0019 #include <trackreco/PHActsTrackProjection.h>
0020 #include <trackingdiagnostics/TrackResiduals.h>
0021
0022 #include <trackbase_historic/SvtxTrack.h>
0023
0024
0025 #include <ffamodules/CDBInterface.h>
0026
0027 #include <fun4all/Fun4AllDstInputManager.h>
0028 #include <fun4all/Fun4AllDstOutputManager.h>
0029 #include <fun4all/Fun4AllInputManager.h>
0030 #include <fun4all/Fun4AllUtils.h>
0031 #include <fun4all/Fun4AllOutputManager.h>
0032 #include <fun4all/Fun4AllRunNodeInputManager.h>
0033 #include <fun4all/Fun4AllServer.h>
0034
0035 #include <phool/recoConsts.h>
0036
0037 #include <iostream>
0038 #include <filesystem>
0039
0040 R__LOAD_LIBRARY(libfun4all.so)
0041 R__LOAD_LIBRARY(libffamodules.so)
0042 R__LOAD_LIBRARY(libmvtx.so)
0043 R__LOAD_LIBRARY(libintt.so)
0044 R__LOAD_LIBRARY(libtpc.so)
0045 R__LOAD_LIBRARY(libmicromegas.so)
0046 R__LOAD_LIBRARY(libTrackingDiagnostics.so)
0047 R__LOAD_LIBRARY(libtrack_reco.so)
0048 R__LOAD_LIBRARY(libTpcDVCalib.so)
0049 R__LOAD_LIBRARY(libcalo_reco.so)
0050
0051 std::string GetFirstLine(const std::string& listname);
0052 bool is_directory_empty(const std::filesystem::path& dir_path);
0053
0054 void Fun4All_FieldOnAllTrackersCalos(
0055     const int nEvents = 10,
0056     std::vector<std::string> myInputLists = {
0057         "run46730_0000_trkr.txt",
0058         "run46730_calo.list"},
0059     const std::string& outDir = "./",
0060     bool doTpcOnlyTracking = true,
0061     float initial_driftvelocity = 0.00710,
0062     bool doEMcalRadiusCorr = true,
0063     const bool convertSeeds = false)
0064 {
0065   int verbosity = 0;
0066
0067   gSystem->Load("libg4dst.so");
0068
0069   std::string firstfile = GetFirstLine(myInputLists[0]);
0070   if (*firstfile.c_str() == '\0') return;
0071   std::pair<int, int> runseg = Fun4AllUtils::GetRunSegment(firstfile);
0072   int runnumber = runseg.first;
0073   int segment = runseg.second;
0074
0075   auto *rc = recoConsts::instance();
0076   rc->set_IntFlag("RUNNUMBER", runnumber);
0077
0078   Enable::CDB = true;
0079   rc->set_StringFlag("CDB_GLOBALTAG", "ProdA_2024");
0080   rc->set_uint64Flag("TIMESTAMP", 6);
0081
0082   Enable::DSTOUT = false;
0083
0084   TpcReadoutInit( runnumber );
0085   std::cout<< " run: " << runnumber
0086        << " samples: " << TRACKING::reco_tpc_maxtime_sample
0087        << " pre: " << TRACKING::reco_tpc_time_presample
0088        << " vdrift: " << G4TPC::tpc_drift_velocity_reco
0089        << std::endl;
0090
0091   //Create the server
0092   Fun4AllServer *se = Fun4AllServer::instance();
0093   se->Verbosity(verbosity);
0094
0095   std::cout << "Including " << myInputLists.size() << " files." << std::endl;
0096
0097   //Add all required input files
0098   for (unsigned int i = 0; i < myInputLists.size(); ++i)
0099   {
0100     Fun4AllInputManager *infile = new Fun4AllDstInputManager("DSTin_" + std::to_string(i));
0101     std::cout << "Including file " << myInputLists[i] << std::endl;
0102     //infile->AddFile(myInputLists[i]);
0103     infile->AddListFile(myInputLists[i]);
0104     se->registerInputManager(infile);
0105   }
0106
0107   std::string outputAnaFileName = "TrackCalo_" + std::to_string(segment) + "_ana.root";
0108   std::string outputDstFileName = "TrackCalo_" + std::to_string(segment) + "_dst.root";
0109
0110   std::string outputRecoDir = outDir + "/inReconstruction/" + std::to_string(runnumber) + "/";
0111   std::string makeDirectory = "mkdir -p " + outputRecoDir;
0112   system(makeDirectory.c_str());
0113   std::string outputAnaFile = outputRecoDir + outputAnaFileName;
0114   std::string outputDstFile = outputRecoDir + outputDstFileName;
0115
0116   std::cout << "Reco ANA file: " << outputAnaFile << std::endl;
0117   std::cout << "Dst file: " << outputDstFile << std::endl;
0118
0119   std::string geofile = CDBInterface::instance()->getUrl("Tracking_Geometry");
0120
0121   Fun4AllRunNodeInputManager *ingeo = new Fun4AllRunNodeInputManager("GeoIn");
0122   ingeo->AddFile(geofile);
0123   se->registerInputManager(ingeo);
0124
0125   //G4TPC::tpc_drift_velocity_reco = (8.0 / 1000) * 107.0 / 105.0 * (1+ dz_separation_0 / 2. / 105.); //cm/ns
0126   G4TPC::tpc_drift_velocity_reco = initial_driftvelocity; //cm/ns
0127   G4TPC::ENABLE_MODULE_EDGE_CORRECTIONS = true;
0128   //to turn on the default static corrections, enable the two lines below
0129   //G4TPC::ENABLE_STATIC_CORRECTIONS = true;
0130   //G4TPC::DISTORTIONS_USE_PHI_AS_RADIANS = false;
0131   //G4MAGNET::magfield = "0.01";
0132   G4MAGNET::magfield_tracking = G4MAGNET::magfield;
0133   G4MAGNET::magfield_rescale = 1;
0134
0135   G4TRACKING::convert_seeds_to_svtxtracks = convertSeeds;
0136   std::cout << "Converting to seeds : " << G4TRACKING::convert_seeds_to_svtxtracks << std::endl;
0137
0138   TRACKING::pp_mode = false;
0139
0140   TrackingInit();
0141   if(doTpcOnlyTracking)
0142   {
0143     Tpc_HitUnpacking();
0144   }
0145   else
0146   {
0147     Mvtx_HitUnpacking();
0148     Intt_HitUnpacking();
0149     Tpc_HitUnpacking();
0150     Micromegas_HitUnpacking();
0151   }
0152
0153   Mvtx_Clustering();
0154   Intt_Clustering();
0155
0156   auto *tpcclusterizer = new TpcClusterizer;
0157   tpcclusterizer->Verbosity(verbosity);
0158   tpcclusterizer->set_do_hit_association(G4TPC::DO_HIT_ASSOCIATION);
0159   tpcclusterizer->set_rawdata_reco();
0160   tpcclusterizer->set_do_sequential(true);
0161   se->registerSubsystem(tpcclusterizer);
0162
0163   Micromegas_Clustering();
0164
0165   /*
0166    * Begin Track Seeding
0167    */
0168
0169   /*
0170    * Silicon Seeding
0171    */
0172   auto *silicon_Seeding = new PHActsSiliconSeeding;
0173   silicon_Seeding->Verbosity(verbosity);
0174   silicon_Seeding->searchInIntt();
0175   silicon_Seeding->setinttRPhiSearchWindow(0.4);
0176   silicon_Seeding->setinttZSearchWindow(1.6);
0177   silicon_Seeding->seedAnalysis(false);
0178   se->registerSubsystem(silicon_Seeding);
0179
0180   auto *merger = new PHSiliconSeedMerger;
0181   merger->Verbosity(verbosity);
0182   se->registerSubsystem(merger);
0183
0184   /*
0185    * Tpc Seeding
0186    */
0187   auto *seeder = new PHCASeeding("PHCASeeding");
0188   double fieldstrength = std::numeric_limits<double>::quiet_NaN();  // set by isConstantField if constant
0189   bool ConstField = isConstantField(G4MAGNET::magfield_tracking, fieldstrength);
0190   if (ConstField)
0191   {
0192     seeder->useConstBField(true);
0193     seeder->constBField(fieldstrength);
0194   }
0195   else
0196   {
0197     seeder->set_field_dir(-1 * G4MAGNET::magfield_rescale);
0198     seeder->useConstBField(false);
0199     seeder->magFieldFile(G4MAGNET::magfield_tracking);  // to get charge sign right
0200   }
0201   seeder->Verbosity(verbosity);
0202   seeder->SetLayerRange(7, 55);
0203   seeder->SetSearchWindow(2.,0.05); // z-width and phi-width, default in macro at 1.5 and 0.05
0204   seeder->SetClusAdd_delta_window(3.0,0.06); //  (0.5, 0.005) are default; sdzdr_cutoff, d2/dr2(phi)_cutoff
0205   //seeder->SetNClustersPerSeedRange(4,60); // default is 6, 6
0206   seeder->SetMinHitsPerCluster(0);
0207   seeder->SetMinClustersPerTrack(3);
0208   seeder->useFixedClusterError(true);
0209   seeder->set_pp_mode(TRACKING::pp_mode);
0210   se->registerSubsystem(seeder);
0211
0212   // expand stubs in the TPC using simple kalman filter
0213   auto *cprop = new PHSimpleKFProp("PHSimpleKFProp");
0214   cprop->set_field_dir(G4MAGNET::magfield_rescale);
0215   if (ConstField)
0216   {
0217     cprop->useConstBField(true);
0218     cprop->setConstBField(fieldstrength);
0219   }
0220   else
0221   {
0222     cprop->magFieldFile(G4MAGNET::magfield_tracking);
0223     cprop->set_field_dir(-1 * G4MAGNET::magfield_rescale);
0224   }
0225   cprop->useFixedClusterError(true);
0226   cprop->set_max_window(5.);
0227   cprop->Verbosity(verbosity);
0228   cprop->set_pp_mode(TRACKING::pp_mode);
0229   se->registerSubsystem(cprop);
0230
0231   /*
0232    * Track Matching between silicon and TPC
0233    */
0234   // The normal silicon association methods
0235   // Match the TPC track stubs from the CA seeder to silicon track stubs from PHSiliconTruthTrackSeeding
0236   auto *silicon_match = new PHSiliconTpcTrackMatching;
0237   silicon_match->Verbosity(verbosity);
0238   silicon_match->set_x_search_window(2.);
0239   silicon_match->set_y_search_window(2.);
0240   silicon_match->set_z_search_window(5.);
0241   silicon_match->set_phi_search_window(0.2);
0242   silicon_match->set_eta_search_window(0.1);
0243 //  silicon_match->set_use_old_matching(true); // deprecated
0244   silicon_match->set_pp_mode(true);
0245   se->registerSubsystem(silicon_match);
0246
0247   // Match TPC track stubs from CA seeder to clusters in the micromegas layers
0248   auto *mm_match = new PHMicromegasTpcTrackMatching;
0249   mm_match->Verbosity(verbosity);
0250   mm_match->set_pp_mode(TRACKING::pp_mode);
0251   mm_match->set_rphi_search_window_lyr1(3.);
0252   mm_match->set_rphi_search_window_lyr2(15.0);
0253   mm_match->set_z_search_window_lyr1(30.0);
0254   mm_match->set_z_search_window_lyr2(3.);
0255
0256   mm_match->set_min_tpc_layer(38);             // layer in TPC to start projection fit
0257   mm_match->set_test_windows_printout(false);  // used for tuning search windows only
0258   se->registerSubsystem(mm_match);
0259
0260   /*
0261    * End Track Seeding
0262    */
0263
0264   /*
0265    * Either converts seeds to tracks with a straight line/helix fit
0266    * or run the full Acts track kalman filter fit
0267    */
0268   if (G4TRACKING::convert_seeds_to_svtxtracks)
0269   {
0270     auto *converter = new TrackSeedTrackMapConverter;
0271     // Option to use TpcTrackSeedContainer or SvtxTrackSeeds
0272     // can be set to SiliconTrackSeedContainer for silicon-only track fit
0273     if (doTpcOnlyTracking)
0274     {
0275       converter->setTrackSeedName("TpcTrackSeedContainer");
0276     }
0277     else
0278     {
0279       converter->setTrackSeedName("SvtxTrackSeeds");
0280     }
0281     converter->setFieldMap(G4MAGNET::magfield_tracking);
0282     converter->Verbosity(verbosity);
0283     se->registerSubsystem(converter);
0284   }
0285   else
0286   {
0287     auto *deltazcorr = new PHTpcDeltaZCorrection;
0288     deltazcorr->Verbosity(verbosity);
0289     se->registerSubsystem(deltazcorr);
0290
0291     // perform final track fit with ACTS
0292     auto *actsFit = new PHActsTrkFitter;
0293     actsFit->Verbosity(verbosity);
0294     actsFit->commissioning(G4TRACKING::use_alignment);
0295     // in calibration mode, fit only Silicons and Micromegas hits
0296     if (!doTpcOnlyTracking)
0297     {
0298       actsFit->fitSiliconMMs(G4TRACKING::SC_CALIBMODE);
0299       actsFit->setUseMicromegas(G4TRACKING::SC_USE_MICROMEGAS);
0300     }
0301     actsFit->set_pp_mode(TRACKING::pp_mode);
0302     actsFit->set_use_clustermover(true);  // default is true for now
0303     actsFit->useActsEvaluator(false);
0304     actsFit->useOutlierFinder(false);
0305     actsFit->setFieldMap(G4MAGNET::magfield_tracking);
0306     se->registerSubsystem(actsFit);
0307   }
0308
0309   PHSimpleVertexFinder *finder = new PHSimpleVertexFinder;
0310   finder->Verbosity(verbosity);
0311   finder->setDcaCut(0.5);
0312   finder->setTrackPtCut(-99999.);
0313   finder->setBeamLineCut(1);
0314   finder->setTrackQualityCut(1000000000);
0315   if (!doTpcOnlyTracking)
0316   {
0317     finder->setRequireMVTX(true);
0318     finder->setNmvtxRequired(3);
0319   }
0320   else
0321   {
0322     finder->setRequireMVTX(false);
0323   }
0324   finder->setOutlierPairCut(0.1);
0325   se->registerSubsystem(finder);
0326
0327   Global_Reco();
0328
0329   auto *projection = new PHActsTrackProjection("CaloProjection");
0330   float new_cemc_rad = 100.70;//(1-(-0.077))*93.5 recommended cemc radius
0331   if (doEMcalRadiusCorr)
0332   {
0333     projection->setLayerRadius(SvtxTrack::CEMC, new_cemc_rad);
0334   }
0335   se->registerSubsystem(projection);
0336
0337   RawClusterBuilderTopo* ClusterBuilder1 = new RawClusterBuilderTopo("EMcalRawClusterBuilderTopo1");
0338   ClusterBuilder1->Verbosity(verbosity);
0339   ClusterBuilder1->set_nodename("TOPOCLUSTER_EMCAL");
0340   ClusterBuilder1->set_enable_HCal(false);
0341   ClusterBuilder1->set_enable_EMCal(true);
0342   //ClusterBuilder1->set_noise(0.0025, 0.006, 0.03);
0343   ClusterBuilder1->set_noise(0.01, 0.03, 0.03);
0344   ClusterBuilder1->set_significance(4.0, 2.0, 1.0);
0345   ClusterBuilder1->allow_corner_neighbor(true);
0346   ClusterBuilder1->set_do_split(true);
0347   ClusterBuilder1->set_minE_local_max(1.0, 2.0, 0.5);
0348   ClusterBuilder1->set_R_shower(0.025);
0349   se->registerSubsystem(ClusterBuilder1);
0350
0351   TrackToCalo *ttc = new TrackToCalo("Tracks_And_Calo", outputAnaFile);
0352   ttc->EMcalRadiusUser(true);
0353   ttc->setEMcalRadius(new_cemc_rad);
0354   se->registerSubsystem(ttc);
0355
0356   if (Enable::DSTOUT)
0357   {
0358     Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", outputDstFile);
0359     out->StripNode("TPC");
0360     out->StripNode("Sync");
0361     out->StripNode("MBD");
0362     out->StripNode("ZDC");
0363     out->StripNode("SEPD");
0364     out->StripNode("HCALIN");
0365     out->StripNode("HCALOUT");
0366     out->StripNode("alignmentTransformationContainer");
0367     out->StripNode("alignmentTransformationContainerTransient");
0368     out->StripNode("SiliconTrackSeedContainer");
0369     out->StripNode("RUN");
0370     out->SaveRunNode(0);
0371     se->registerOutputManager(out);
0372   }
0373
0374   se->run(nEvents);
0375   se->End();
0376   se->PrintTimer();
0377
0378   std::ifstream file_ana(outputAnaFile.c_str(), std::ios::binary | std::ios::ate);
0379   if (file_ana.good() && (file_ana.tellg() > 100))
0380   {
0381     std::string outputRecoDirMove = outDir + "/Reconstructed/" + std::to_string(runnumber) + "/";
0382     std::string makeDirectoryMove = "mkdir -p " + outputRecoDirMove;
0383     system(makeDirectoryMove.c_str());
0384     std::string moveOutput = "mv " + outputAnaFile + " " + outDir + "/Reconstructed/" + std::to_string(runnumber);
0385     std::cout << "moveOutput: " << moveOutput << std::endl;
0386     system(moveOutput.c_str());
0387   }
0388
0389   std::ifstream file_dst(outputDstFile.c_str(), std::ios::binary | std::ios::ate);
0390   if (file_dst.good() && (file_dst.tellg() > 100))
0391   {
0392     std::string outputRecoDirMove = outDir + "/Reconstructed/" + std::to_string(runnumber) + "/";
0393     std::string makeDirectoryMove = "mkdir -p " + outputRecoDirMove;
0394     system(makeDirectoryMove.c_str());
0395     std::string moveOutput = "mv " + outputDstFile + " " + outDir + "/Reconstructed/" + std::to_string(runnumber);
0396     std::cout << "moveOutput: " << moveOutput << std::endl;
0397     system(moveOutput.c_str());
0398   }
0399
0400   delete se;
0401   std::cout << "All done" << std::endl;
0402   gSystem->Exit(0);
0403
0404   return;
0405 }
0406
0407 std::string GetFirstLine(const std::string& listname)
0408 {
0409   std::ifstream file(listname);
0410
0411   std::string firstLine;
0412   if (file.is_open()) {
0413       if (std::getline(file, firstLine)) {
0414           std::cout << "First Line: " << firstLine << std::endl;
0415       } else {
0416           std::cerr << "Unable to read first line of file" << std::endl;
0417       }
0418       file.close();
0419   } else {
0420       std::cerr << "Unable to open file" << std::endl;
0421   }
0422   return firstLine;
0423 }
0424
0425 bool is_directory_empty(const std::filesystem::path& dir_path) {
0426     if (std::filesystem::exists(dir_path) && std::filesystem::is_directory(dir_path)) {
0427         return std::filesystem::directory_iterator(dir_path) == std::filesystem::directory_iterator();
0428     }
0429     return false;
0430 }