tpc/TpcDVCalib/Fun4All_FieldOnAllTrackersCalos.C

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