macros/TrackingProduction/Fun4All_FullReconstruction.C

0001 /*
0002  * This macro shows a minimum working example of running the tracking
0003  * hit unpackers with some basic seeding algorithms to try to put together
0004  * tracks. There are some analysis modules run at the end which package
0005  * hits, clusters, and clusters on tracks into trees for analysis.
0006  */
0007
0008 // leave the GlobalVariables.C at the beginning, an empty line afterwards
0009 // protects its position against reshuffling by clang-format
0010 #include <GlobalVariables.C>
0011
0012 #include <G4_ActsGeom.C>
0013 #include <G4_Global.C>
0014 #include <G4_Magnet.C>
0015 #include <G4_Mbd.C>
0016 #include <QA.C>
0017 #include <Trkr_Clustering.C>
0018 #include <Trkr_LaserClustering.C>
0019 #include <Trkr_Reco.C>
0020 #include <Trkr_RecoInit.C>
0021 #include <Trkr_TpcReadoutInit.C>
0022
0023 #include <tpccalib/PHTpcResiduals.h>
0024
0025 #include <tpcqa/TpcRawHitQA.h>
0026
0027 #include <trackingqa/InttClusterQA.h>
0028 #include <trackingqa/MicromegasClusterQA.h>
0029 #include <trackingqa/MvtxClusterQA.h>
0030 #include <trackingqa/TpcClusterQA.h>
0031 #include <trackingqa/TpcSeedsQA.h>
0032
0033 #include <trackingdiagnostics/TrackResiduals.h>
0034 #include <trackingdiagnostics/TrkrNtuplizer.h>
0035
0036 #include <ffamodules/CDBInterface.h>
0037
0038 #include <fun4all/Fun4AllDstInputManager.h>
0039 #include <fun4all/Fun4AllDstOutputManager.h>
0040 #include <fun4all/Fun4AllInputManager.h>
0041 #include <fun4all/Fun4AllOutputManager.h>
0042 #include <fun4all/Fun4AllRunNodeInputManager.h>
0043 #include <fun4all/Fun4AllServer.h>
0044 #include <fun4all/Fun4AllUtils.h>
0045
0046 #include <phool/recoConsts.h>
0047
0048 R__LOAD_LIBRARY(libfun4all.so)
0049 R__LOAD_LIBRARY(libffamodules.so)
0050 R__LOAD_LIBRARY(libphool.so)
0051 R__LOAD_LIBRARY(libcdbobjects.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(libtrackingqa.so)
0058 R__LOAD_LIBRARY(libtpcqa.so)
0059
0060 void Fun4All_FullReconstruction(
0061     const int nEvents = 10,
0062     const std::string &filelist = "filelist.list",
0063     const std::string &outfilename = "clusters_seeds",
0064     const bool convertSeeds = false)
0065 {
0066   G4TRACKING::convert_seeds_to_svtxtracks = convertSeeds;
0067   std::cout << "Converting to seeds : " << G4TRACKING::convert_seeds_to_svtxtracks << std::endl;
0068
0069   auto *se = Fun4AllServer::instance();
0070   se->Verbosity(2);
0071   auto *rc = recoConsts::instance();
0072
0073   std::ifstream ifs(filelist);
0074   std::string filepath;
0075   int runnumber = std::numeric_limits<int>::quiet_NaN();
0076   int segment = std::numeric_limits<int>::quiet_NaN();
0077   int i = 0;
0078   while (std::getline(ifs, filepath))
0079   {
0080     std::cout << "Adding DST with filepath: " << filepath << std::endl;
0081     if (i == 0)
0082     {
0083       std::pair<int, int> runseg = Fun4AllUtils::GetRunSegment(filepath);
0084       runnumber = runseg.first;
0085       segment = runseg.second;
0086       rc->set_IntFlag("RUNNUMBER", runnumber);
0087       rc->set_uint64Flag("TIMESTAMP", runnumber);
0088     }
0089     std::string inputname = "InputManager" + std::to_string(i);
0090     auto *hitsin = new Fun4AllDstInputManager(inputname);
0091     hitsin->fileopen(filepath);
0092     se->registerInputManager(hitsin);
0093     i++;
0094   }
0095
0096   std::cout << " run: " << runnumber
0097             << " samples: " << TRACKING::reco_tpc_maxtime_sample
0098             << " pre: " << TRACKING::reco_tpc_time_presample
0099             << " vdrift: " << G4TPC::tpc_drift_velocity_reco
0100             << std::endl;
0101
0102   TRACKING::pp_mode = true;
0103
0104   Enable::MVTX_APPLYMISALIGNMENT = true;
0105   ACTSGEOM::mvtx_applymisalignment = Enable::MVTX_APPLYMISALIGNMENT;
0106
0107   // distortion calibration mode
0108   /*
0109    * set to true to enable residuals in the TPC with
0110    * TPC clusters not participating to the ACTS track fit
0111    */
0112   G4TRACKING::SC_CALIBMODE = false;
0113
0114   std::string theOutfile = outfilename + "_" + std::to_string(runnumber) + "-" + std::to_string(segment) + ".root";
0115
0116   Enable::CDB = true;
0117   rc->set_StringFlag("CDB_GLOBALTAG", "newcdbtag");
0118   rc->set_uint64Flag("TIMESTAMP", runnumber);
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   TpcReadoutInit(runnumber);
0126
0127   G4TPC::REJECT_LASER_EVENTS = true;
0128   G4TPC::ENABLE_MODULE_EDGE_CORRECTIONS = true;
0129   // Flag for running the tpc hit unpacker with zero suppression on
0130   TRACKING::tpc_zero_supp = true;
0131
0132   // to turn on the default static corrections, enable the two lines below
0133   G4TPC::ENABLE_STATIC_CORRECTIONS = true;
0134   G4TPC::USE_PHI_AS_RAD_STATIC_CORRECTIONS = false;
0135
0136   // to turn on the average corrections, enable the three lines below
0137   // note: these are designed to be used only if static corrections are also applied
0138   G4TPC::ENABLE_AVERAGE_CORRECTIONS = true;
0139   G4TPC::USE_PHI_AS_RAD_AVERAGE_CORRECTIONS = false;
0140   // to use a custom file instead of the database file:
0141   G4TPC::average_correction_filename = CDBInterface::instance()->getUrl("TPC_LAMINATION_FIT_CORRECTION");
0142
0143   G4MAGNET::magfield_rescale = 1;
0144   TrackingInit();
0145
0146   for (int felix = 0; felix < 6; felix++)
0147   {
0148     Mvtx_HitUnpacking(std::to_string(felix));
0149   }
0150   for (int server = 0; server < 8; server++)
0151   {
0152     Intt_HitUnpacking(std::to_string(server));
0153   }
0154   std::ostringstream ebdcname;
0155   for (int ebdc = 0; ebdc < 24; ebdc++)
0156   {
0157     ebdcname.str("");
0158     if (ebdc < 10)
0159     {
0160       ebdcname << "0";
0161     }
0162     ebdcname << ebdc;
0163     Tpc_HitUnpacking(ebdcname.str());
0164   }
0165
0166   Micromegas_HitUnpacking();
0167
0168   Mvtx_Clustering();
0169
0170   Intt_Clustering();
0171
0172   Tpc_LaserEventIdentifying();
0173
0174   TPC_Clustering_run2pp();
0175
0176   Micromegas_Clustering();
0177
0178   Reject_Laser_Events();
0179   /*
0180    * Begin Track Seeding
0181    */
0182
0183   /*
0184    * Silicon Seeding
0185    */
0186
0187   /*
0188   auto silicon_Seeding = new PHActsSiliconSeeding;
0189   silicon_Seeding->Verbosity(0);
0190   silicon_Seeding->searchInIntt();
0191   silicon_Seeding->setinttRPhiSearchWindow(0.4);
0192   silicon_Seeding->setinttZSearchWindow(1.6);
0193   silicon_Seeding->seedAnalysis(false);
0194   se->registerSubsystem(silicon_Seeding);
0195   */
0196
0197   auto *silicon_Seeding = new PHActsSiliconSeeding;
0198   silicon_Seeding->Verbosity(0);
0199   silicon_Seeding->setStrobeRange(-5, 5);
0200   // these get us to about 83% INTT > 1
0201   silicon_Seeding->setinttRPhiSearchWindow(0.4);
0202   silicon_Seeding->setinttZSearchWindow(2.0);
0203   silicon_Seeding->seedAnalysis(false);
0204   se->registerSubsystem(silicon_Seeding);
0205
0206   auto *merger = new PHSiliconSeedMerger;
0207   merger->Verbosity(0);
0208   se->registerSubsystem(merger);
0209
0210   /*
0211    * Tpc Seeding
0212    */
0213   auto *seeder = new PHCASeeding("PHCASeeding");
0214   double fieldstrength = std::numeric_limits<double>::quiet_NaN();  // set by isConstantField if constant
0215   bool ConstField = isConstantField(G4MAGNET::magfield_tracking, fieldstrength);
0216   if (ConstField)
0217   {
0218     seeder->useConstBField(true);
0219     seeder->constBField(fieldstrength);
0220   }
0221   else
0222   {
0223     seeder->set_field_dir(-1 * G4MAGNET::magfield_rescale);
0224     seeder->useConstBField(false);
0225     seeder->magFieldFile(G4MAGNET::magfield_tracking);  // to get charge sign right
0226   }
0227   seeder->Verbosity(0);
0228   seeder->SetLayerRange(7, 55);
0229   seeder->SetSearchWindow(2., 0.05);           // z-width and phi-width, default in macro at 1.5 and 0.05
0230   seeder->SetClusAdd_delta_window(3.0, 0.06);  //  (0.5, 0.005) are default; sdzdr_cutoff, d2/dr2(phi)_cutoff
0231   // seeder->SetNClustersPerSeedRange(4,60); // default is 6, 6
0232   seeder->SetMinHitsPerCluster(0);
0233   seeder->SetMinClustersPerTrack(3);
0234   seeder->useFixedClusterError(true);
0235   seeder->set_pp_mode(true);
0236   se->registerSubsystem(seeder);
0237
0238   // expand stubs in the TPC using simple kalman filter
0239   auto *cprop = new PHSimpleKFProp("PHSimpleKFProp");
0240   cprop->set_field_dir(G4MAGNET::magfield_rescale);
0241   if (ConstField)
0242   {
0243     cprop->useConstBField(true);
0244     cprop->setConstBField(fieldstrength);
0245   }
0246   else
0247   {
0248     cprop->magFieldFile(G4MAGNET::magfield_tracking);
0249     cprop->set_field_dir(-1 * G4MAGNET::magfield_rescale);
0250   }
0251   cprop->useFixedClusterError(true);
0252   cprop->set_max_window(5.);
0253   cprop->Verbosity(0);
0254   cprop->set_pp_mode(true);
0255   se->registerSubsystem(cprop);
0256
0257   // Always apply preliminary distortion corrections to TPC clusters before silicon matching
0258   // and refit the trackseeds. Replace KFProp fits with the new fit parameters in the TPC seeds.
0259   auto *prelim_distcorr = new PrelimDistortionCorrection;
0260   prelim_distcorr->set_pp_mode(true);
0261   prelim_distcorr->Verbosity(0);
0262   se->registerSubsystem(prelim_distcorr);
0263
0264   /*
0265    * Track Matching between silicon and TPC
0266    */
0267   // The normal silicon association methods
0268   // Match the TPC track stubs from the CA seeder to silicon track stubs from PHSiliconTruthTrackSeeding
0269   auto *silicon_match = new PHSiliconTpcTrackMatching;
0270   silicon_match->Verbosity(0);
0271   silicon_match->set_pp_mode(TRACKING::pp_mode);
0272   if (G4TPC::ENABLE_AVERAGE_CORRECTIONS)
0273   {
0274     // for general tracking
0275     // Eta/Phi window is determined by 3 sigma window
0276     // X/Y/Z window is determined by 4 sigma window
0277     silicon_match->window_deta.set_posQoverpT_maxabs({-0.014, 0.0331, 0.48});
0278     silicon_match->window_deta.set_negQoverpT_maxabs({-0.006, 0.0235, 0.52});
0279     silicon_match->set_deltaeta_min(0.03);
0280     silicon_match->window_dphi.set_QoverpT_range({-0.15, 0, 0}, {0.15, 0, 0});
0281     silicon_match->window_dx.set_QoverpT_maxabs({3.0, 0, 0});
0282     silicon_match->window_dy.set_QoverpT_maxabs({3.0, 0, 0});
0283     silicon_match->window_dz.set_posQoverpT_maxabs({1.138, 0.3919, 0.84});
0284     silicon_match->window_dz.set_negQoverpT_maxabs({0.719, 0.6485, 0.65});
0285     silicon_match->set_crossing_deltaz_max(30);
0286     silicon_match->set_crossing_deltaz_min(2);
0287
0288     // for distortion correction using SI-TPOT fit and track pT>0.5
0289     if (G4TRACKING::SC_CALIBMODE)
0290     {
0291       silicon_match->window_deta.set_posQoverpT_maxabs({0.016, 0.0060, 1.13});
0292       silicon_match->window_deta.set_negQoverpT_maxabs({0.022, 0.0022, 1.44});
0293       silicon_match->set_deltaeta_min(0.03);
0294       silicon_match->window_dphi.set_QoverpT_range({-0.15, 0, 0}, {0.09, 0, 0});
0295       silicon_match->window_dx.set_QoverpT_maxabs({2.0, 0, 0});
0296       silicon_match->window_dy.set_QoverpT_maxabs({1.5, 0, 0});
0297       silicon_match->window_dz.set_posQoverpT_maxabs({1.213, 0.0211, 2.09});
0298       silicon_match->window_dz.set_negQoverpT_maxabs({1.307, 0.0001, 4.52});
0299       silicon_match->set_crossing_deltaz_min(1.2);
0300     }
0301   }
0302   se->registerSubsystem(silicon_match);
0303
0304   // Match TPC track stubs from CA seeder to clusters in the micromegas layers
0305   auto *mm_match = new PHMicromegasTpcTrackMatching;
0306   mm_match->Verbosity(0);
0307   mm_match->set_pp_mode(TRACKING::pp_mode);
0308   mm_match->set_rphi_search_window_lyr1(3.);
0309   mm_match->set_rphi_search_window_lyr2(15.0);
0310   mm_match->set_z_search_window_lyr1(30.0);
0311   mm_match->set_z_search_window_lyr2(3.);
0312
0313   mm_match->set_min_tpc_layer(38);             // layer in TPC to start projection fit
0314   mm_match->set_test_windows_printout(false);  // used for tuning search windows only
0315   se->registerSubsystem(mm_match);
0316
0317   /*
0318    * End Track Seeding
0319    */
0320
0321   /*
0322    * Either converts seeds to tracks with a straight line/helix fit
0323    * or run the full Acts track kalman filter fit
0324    */
0325   if (G4TRACKING::convert_seeds_to_svtxtracks)
0326   {
0327     auto *converter = new TrackSeedTrackMapConverter;
0328     // Default set to full SvtxTrackSeeds. Can be set to
0329     // SiliconTrackSeedContainer or TpcTrackSeedContainer
0330     converter->setTrackSeedName("SvtxTrackSeedContainer");
0331     converter->setFieldMap(G4MAGNET::magfield_tracking);
0332     converter->Verbosity(0);
0333     se->registerSubsystem(converter);
0334   }
0335   else
0336   {
0337     auto *deltazcorr = new PHTpcDeltaZCorrection;
0338     deltazcorr->Verbosity(0);
0339     se->registerSubsystem(deltazcorr);
0340
0341     // perform final track fit with ACTS
0342     auto *actsFit = new PHActsTrkFitter;
0343     actsFit->Verbosity(0);
0344     actsFit->commissioning(G4TRACKING::use_alignment);
0345     // in calibration mode, fit only Silicons and Micromegas hits
0346     actsFit->fitSiliconMMs(G4TRACKING::SC_CALIBMODE);
0347     actsFit->setUseMicromegas(G4TRACKING::SC_USE_MICROMEGAS);
0348     actsFit->set_pp_mode(TRACKING::pp_mode);
0349     actsFit->set_use_clustermover(true);  // default is true for now
0350     actsFit->useActsEvaluator(false);
0351     actsFit->useOutlierFinder(false);
0352     actsFit->setFieldMap(G4MAGNET::magfield_tracking);
0353     se->registerSubsystem(actsFit);
0354
0355     auto *cleaner = new PHTrackCleaner();
0356     cleaner->Verbosity(0);
0357     cleaner->set_pp_mode(TRACKING::pp_mode);
0358     se->registerSubsystem(cleaner);
0359
0360     if (G4TRACKING::SC_CALIBMODE)
0361     {
0362       /*
0363        * in calibration mode, calculate residuals between TPC and fitted tracks,
0364        * store in dedicated structure for distortion correction
0365        */
0366       auto *residuals = new PHTpcResiduals;
0367       std::string tpc_residoutfile = theOutfile + "_PhTpcResiduals.root";
0368       residuals->setOutputfile(tpc_residoutfile);
0369       residuals->setUseMicromegas(G4TRACKING::SC_USE_MICROMEGAS);
0370
0371       // matches Tony's analysis
0372       residuals->setMinPt(0.2);
0373
0374       // reconstructed distortion grid size (phi, r, z)
0375       residuals->setGridDimensions(36, 48, 80);
0376       se->registerSubsystem(residuals);
0377     }
0378   }
0379
0380   auto *finder = new PHSimpleVertexFinder;
0381   finder->Verbosity(0);
0382   finder->setDcaCut(0.5);
0383   finder->setTrackPtCut(-99999.);
0384   finder->setBeamLineCut(1);
0385   finder->setTrackQualityCut(1000000000);
0386   finder->setNmvtxRequired(3);
0387   finder->setOutlierPairCut(0.1);
0388   se->registerSubsystem(finder);
0389
0390   // Propagate track positions to the vertex position
0391   auto *vtxProp = new PHActsVertexPropagator;
0392   vtxProp->Verbosity(0);
0393   vtxProp->fieldMap(G4MAGNET::magfield_tracking);
0394   se->registerSubsystem(vtxProp);
0395
0396   std::string residstring = theOutfile + "_resid.root";
0397
0398   auto *resid = new TrackResiduals("TrackResiduals");
0399   resid->outfileName(residstring);
0400   resid->alignment(false);
0401
0402   // adjust track map name
0403   if (G4TRACKING::SC_CALIBMODE && !G4TRACKING::convert_seeds_to_svtxtracks)
0404   {
0405     resid->trackmapName("SvtxSiliconMMTrackMap");
0406     if (G4TRACKING::SC_USE_MICROMEGAS)
0407     {
0408       resid->set_doMicromegasOnly(true);
0409     }
0410   }
0411
0412   resid->clusterTree();
0413   resid->hitTree();
0414   resid->convertSeeds(G4TRACKING::convert_seeds_to_svtxtracks);
0415
0416   resid->Verbosity(0);
0417   se->registerSubsystem(resid);
0418
0419   // auto ntuplizer = new TrkrNtuplizer("TrkrNtuplizer");
0420   // se->registerSubsystem(ntuplizer);
0421
0422   // Fun4AllOutputManager *out = new Fun4AllDstOutputManager("out", "/sphenix/tg/tg01/hf/jdosbo/tracking_development/Run24/Beam/41626/hitsets.root");
0423   // se->registerOutputManager(out);
0424   if (Enable::QA)
0425   {
0426     se->registerSubsystem(new TpcRawHitQA);
0427     se->registerSubsystem(new MvtxClusterQA);
0428     se->registerSubsystem(new InttClusterQA);
0429     se->registerSubsystem(new TpcClusterQA);
0430     se->registerSubsystem(new MicromegasClusterQA);
0431     auto *tpcqa = new TpcSeedsQA;
0432     tpcqa->setTrackMapName("TpcSvtxTrackMap");
0433     tpcqa->setVertexMapName("TpcSvtxVertexMap");
0434     tpcqa->setSegment(rc->get_IntFlag("RUNSEGMENT"));
0435     se->registerSubsystem(tpcqa);
0436   }
0437   se->run(nEvents);
0438   se->End();
0439   se->PrintTimer();
0440   CDBInterface::instance()->Print();
0441   if (Enable::QA)
0442   {
0443     std::string qaOutputFileName = theOutfile + "_qa.root";
0444     QAHistManagerDef::saveQARootFile(qaOutputFileName);
0445   }
0446
0447   delete se;
0448   std::cout << "Finished" << std::endl;
0449   gSystem->Exit(0);
0450 }