macros/TrackingProduction/Fun4All_TrackFitting.C

0001 /*
0002  * This macro shows a minimum working example of running track fitting over
0003  * the production cluster and track seed DSTs.. There are some analysis
0004  * modules run at the end which package clusters, and clusters on tracks
0005  * into trees for analysis.
0006  */
0007
0008 #include <fun4all/Fun4AllUtils.h>
0009 #include <G4_ActsGeom.C>
0010 #include <G4_Global.C>
0011 #include <G4_Magnet.C>
0012 #include <GlobalVariables.C>
0013 #include <QA.C>
0014 #include <Trkr_Clustering.C>
0015 #include <Trkr_Reco.C>
0016 #include <Trkr_RecoInit.C>
0017 #include <Trkr_TpcReadoutInit.C>
0018
0019 #include <ffamodules/CDBInterface.h>
0020
0021 #include <fun4all/Fun4AllDstInputManager.h>
0022 #include <fun4all/Fun4AllDstOutputManager.h>
0023 #include <fun4all/Fun4AllInputManager.h>
0024 #include <fun4all/Fun4AllOutputManager.h>
0025 #include <fun4all/Fun4AllRunNodeInputManager.h>
0026 #include <fun4all/Fun4AllServer.h>
0027
0028 #include <phool/recoConsts.h>
0029
0030 #include <cdbobjects/CDBTTree.h>
0031
0032 #include <tpccalib/PHTpcResiduals.h>
0033
0034 #include <trackingqa/SiliconSeedsQA.h>
0035 #include <trackingqa/TpcSeedsQA.h>
0036 #include <trackingqa/TpcSiliconQA.h>
0037
0038 #include <trackingdiagnostics/TrackResiduals.h>
0039 #include <trackingdiagnostics/TrkrNtuplizer.h>
0040
0041 #include <stdio.h>
0042
0043 R__LOAD_LIBRARY(libfun4all.so)
0044 R__LOAD_LIBRARY(libffamodules.so)
0045 R__LOAD_LIBRARY(libphool.so)
0046 R__LOAD_LIBRARY(libcdbobjects.so)
0047 R__LOAD_LIBRARY(libTrackingDiagnostics.so)
0048 R__LOAD_LIBRARY(libtrackingqa.so)
0049 void Fun4All_TrackFitting(
0050     const int nEvents = 10,
0051     const std::string seedfilename = "/sphenix/lustre01/sphnxpro/production/run2pp/physics/ana494_2024p021_v001/DST_TRKR_SEED/run_00053800_00053900/dst/DST_TRKR_SEED_run2pp_ana494_2024p021_v001-00053877-00000.root",
0052     const std::string clusterfilename = "/sphenix/lustre01/sphnxpro/production/run2pp/physics/ana494_2024p021_v001/DST_TRKR_CLUSTER/run_00053800_00053900/dst/DST_TRKR_CLUSTER_run2pp_ana494_2024p021_v001-00053877-00000.root",
0053     const std::string outfilename = "clusters_seeds",
0054     const bool convertSeeds = false)
0055 {
0056
0057   G4TRACKING::convert_seeds_to_svtxtracks = convertSeeds;
0058   std::cout << "Converting to seeds : " << G4TRACKING::convert_seeds_to_svtxtracks << std::endl;
0059   std::pair<int, int>
0060       runseg = Fun4AllUtils::GetRunSegment(seedfilename);
0061   int runnumber = runseg.first;
0062   int segment = runseg.second;
0063
0064   std::cout << " run: " << runnumber
0065             << " samples: " << TRACKING::reco_tpc_maxtime_sample
0066             << " pre: " << TRACKING::reco_tpc_time_presample
0067             << " vdrift: " << G4TPC::tpc_drift_velocity_reco
0068             << std::endl;
0069
0070   // distortion calibration mode
0071   /*
0072    * set to true to enable residuals in the TPC with
0073    * TPC clusters not participating to the ACTS track fit
0074    */
0075   G4TRACKING::SC_CALIBMODE = false;
0076   Enable::MVTX_APPLYMISALIGNMENT = true;
0077   ACTSGEOM::mvtx_applymisalignment = Enable::MVTX_APPLYMISALIGNMENT;
0078   TRACKING::pp_mode = true;
0079
0080   TString outfile = outfilename + "_" + runnumber + "-" + segment + ".root";
0081
0082   std::string theOutfile = outfile.Data();
0083
0084   auto se = Fun4AllServer::instance();
0085   se->Verbosity(1);
0086
0087   auto rc = recoConsts::instance();
0088   rc->set_IntFlag("RUNNUMBER", runnumber);
0089
0090   Enable::CDB = true;
0091   rc->set_StringFlag("CDB_GLOBALTAG", "newcdbtag");
0092   rc->set_uint64Flag("TIMESTAMP", runnumber);
0093   std::string geofile = CDBInterface::instance()->getUrl("Tracking_Geometry");
0094
0095   Fun4AllRunNodeInputManager *ingeo = new Fun4AllRunNodeInputManager("GeoIn");
0096   ingeo->AddFile(geofile);
0097   se->registerInputManager(ingeo);
0098
0099   TpcReadoutInit( runnumber );
0100   // these lines show how to override the drift velocity and time offset values set in TpcReadoutInit
0101   // G4TPC::tpc_drift_velocity_reco = 0.0073844; // cm/ns
0102   // TpcClusterZCrossingCorrection::_vdrift = G4TPC::tpc_drift_velocity_reco;
0103   // G4TPC::tpc_tzero_reco = -5*50;  // ns
0104
0105   G4TPC::ENABLE_MODULE_EDGE_CORRECTIONS = true;
0106
0107   // to turn on the default static corrections, enable the two lines below
0108   G4TPC::ENABLE_STATIC_CORRECTIONS = true;
0109   G4TPC::USE_PHI_AS_RAD_STATIC_CORRECTIONS = false;
0110
0111   //to turn on the average corrections, enable the three lines below
0112   //note: these are designed to be used only if static corrections are also applied
0113   G4TPC::ENABLE_AVERAGE_CORRECTIONS = true;
0114   G4TPC::USE_PHI_AS_RAD_AVERAGE_CORRECTIONS = false;
0115    // to use a custom file instead of the database file:
0116   G4TPC::average_correction_filename = CDBInterface::instance()->getUrl("TPC_LAMINATION_FIT_CORRECTION");
0117   G4MAGNET::magfield_rescale = 1;
0118   TrackingInit();
0119
0120   auto hitsinseed = new Fun4AllDstInputManager("SeedInputManager");
0121   hitsinseed->fileopen(seedfilename);
0122   se->registerInputManager(hitsinseed);
0123
0124   auto hitsinclus = new Fun4AllDstInputManager("ClusterInputManager");
0125   hitsinclus->fileopen(clusterfilename);
0126   se->registerInputManager(hitsinclus);
0127
0128
0129
0130   /*
0131    * Track Matching between silicon and TPC
0132    */
0133   // The normal silicon association methods
0134   // Match the TPC track stubs from the CA seeder to silicon track stubs from PHSiliconTruthTrackSeeding
0135   auto silicon_match = new PHSiliconTpcTrackMatching;
0136   silicon_match->Verbosity(0);
0137   silicon_match->set_pp_mode(TRACKING::pp_mode);
0138   if(G4TPC::ENABLE_AVERAGE_CORRECTIONS)
0139   {
0140     // for general tracking
0141     // Eta/Phi window is determined by 3 sigma window
0142     // X/Y/Z window is determined by 4 sigma window
0143     silicon_match->window_deta.set_posQoverpT_maxabs({-0.014,0.0331,0.48});
0144     silicon_match->window_deta.set_negQoverpT_maxabs({-0.006,0.0235,0.52});
0145     silicon_match->set_deltaeta_min(0.03);
0146     silicon_match->window_dphi.set_QoverpT_range({-0.15,0,0}, {0.15,0,0});
0147     silicon_match->window_dx.set_QoverpT_maxabs({3.0,0,0});
0148     silicon_match->window_dy.set_QoverpT_maxabs({3.0,0,0});
0149     silicon_match->window_dz.set_posQoverpT_maxabs({1.138,0.3919,0.84});
0150     silicon_match->window_dz.set_negQoverpT_maxabs({0.719,0.6485,0.65});
0151     silicon_match->set_crossing_deltaz_max(30);
0152     silicon_match->set_crossing_deltaz_min(2);
0153
0154     // for distortion correction using SI-TPOT fit and track pT>0.5
0155     if (G4TRACKING::SC_CALIBMODE)
0156     {
0157       silicon_match->window_deta.set_posQoverpT_maxabs({0.016,0.0060,1.13});
0158       silicon_match->window_deta.set_negQoverpT_maxabs({0.022,0.0022,1.44});
0159       silicon_match->set_deltaeta_min(0.03);
0160       silicon_match->window_dphi.set_QoverpT_range({-0.15,0,0}, {0.09,0,0});
0161       silicon_match->window_dx.set_QoverpT_maxabs({2.0,0,0});
0162       silicon_match->window_dy.set_QoverpT_maxabs({1.5,0,0});
0163       silicon_match->window_dz.set_posQoverpT_maxabs({1.213,0.0211,2.09});
0164       silicon_match->window_dz.set_negQoverpT_maxabs({1.307,0.0001,4.52});
0165       silicon_match->set_crossing_deltaz_min(1.2);
0166     }
0167   }
0168   se->registerSubsystem(silicon_match);
0169
0170   // Match TPC track stubs from CA seeder to clusters in the micromegas layers
0171   auto mm_match = new PHMicromegasTpcTrackMatching;
0172   mm_match->Verbosity(0);
0173   mm_match->set_pp_mode(TRACKING::pp_mode);
0174   mm_match->set_rphi_search_window_lyr1(3.);
0175   mm_match->set_rphi_search_window_lyr2(15.0);
0176   mm_match->set_z_search_window_lyr1(30.0);
0177   mm_match->set_z_search_window_lyr2(3.);
0178
0179   mm_match->set_min_tpc_layer(38);             // layer in TPC to start projection fit
0180   mm_match->set_test_windows_printout(false);  // used for tuning search windows only
0181   se->registerSubsystem(mm_match);
0182
0183
0184   /*
0185    * Either converts seeds to tracks with a straight line/helix fit
0186    * or run the full Acts track kalman filter fit
0187    */
0188   if (G4TRACKING::convert_seeds_to_svtxtracks)
0189   {
0190     auto converter = new TrackSeedTrackMapConverter;
0191     // Default set to full SvtxTrackSeeds. Can be set to
0192     // SiliconTrackSeedContainer or TpcTrackSeedContainer
0193     converter->setTrackSeedName("SvtxTrackSeedContainer");
0194     converter->setFieldMap(G4MAGNET::magfield_tracking);
0195     converter->Verbosity(0);
0196     se->registerSubsystem(converter);
0197   }
0198   else
0199   {
0200     auto deltazcorr = new PHTpcDeltaZCorrection;
0201     deltazcorr->Verbosity(0);
0202     se->registerSubsystem(deltazcorr);
0203
0204     // perform final track fit with ACTS
0205     auto actsFit = new PHActsTrkFitter;
0206     actsFit->Verbosity(0);
0207     actsFit->commissioning(G4TRACKING::use_alignment);
0208     // in calibration mode, fit only Silicons and Micromegas hits
0209     actsFit->fitSiliconMMs(G4TRACKING::SC_CALIBMODE);
0210     actsFit->setUseMicromegas(G4TRACKING::SC_USE_MICROMEGAS);
0211     actsFit->set_pp_mode(TRACKING::pp_mode);
0212     actsFit->set_use_clustermover(true);  // default is true for now
0213     actsFit->useActsEvaluator(false);
0214     actsFit->useOutlierFinder(false);
0215     actsFit->setFieldMap(G4MAGNET::magfield_tracking);
0216     se->registerSubsystem(actsFit);
0217
0218     auto cleaner = new PHTrackCleaner();
0219     cleaner->Verbosity(0);
0220     cleaner->set_pp_mode(TRACKING::pp_mode);
0221     se->registerSubsystem(cleaner);
0222
0223     if (G4TRACKING::SC_CALIBMODE)
0224     {
0225       /*
0226        * in calibration mode, calculate residuals between TPC and fitted tracks,
0227        * store in dedicated structure for distortion correction
0228        */
0229       auto residuals = new PHTpcResiduals;
0230       const TString tpc_residoutfile = theOutfile + "_PhTpcResiduals.root";
0231       residuals->setOutputfile(tpc_residoutfile.Data());
0232       residuals->setUseMicromegas(G4TRACKING::SC_USE_MICROMEGAS);
0233
0234       // matches Tony's analysis
0235       residuals->setMinPt(0.2);
0236
0237       // reconstructed distortion grid size (phi, r, z)
0238       residuals->setGridDimensions(36, 48, 80);
0239       se->registerSubsystem(residuals);
0240     }
0241   }
0242
0243
0244   auto finder = new PHSimpleVertexFinder;
0245   finder->Verbosity(0);
0246
0247   //new cuts
0248   finder->setDcaCut(0.05);
0249   finder->setTrackPtCut(0.1);
0250   finder->setBeamLineCut(1);
0251   finder->setTrackQualityCut(300);
0252   finder->setNmvtxRequired(3);
0253   finder->setOutlierPairCut(0.10);
0254
0255   se->registerSubsystem(finder);
0256
0257   // Propagate track positions to the vertex position
0258   auto vtxProp = new PHActsVertexPropagator;
0259   vtxProp->Verbosity(0);
0260   vtxProp->fieldMap(G4MAGNET::magfield_tracking);
0261   se->registerSubsystem(vtxProp);
0262
0263
0264
0265   TString residoutfile = theOutfile + "_resid.root";
0266   std::string residstring(residoutfile.Data());
0267
0268   auto resid = new TrackResiduals("TrackResiduals");
0269   resid->outfileName(residstring);
0270   resid->alignment(false);
0271
0272   // adjust track map name
0273   if (G4TRACKING::SC_CALIBMODE && !G4TRACKING::convert_seeds_to_svtxtracks)
0274   {
0275     resid->trackmapName("SvtxSiliconMMTrackMap");
0276     if (G4TRACKING::SC_USE_MICROMEGAS)
0277     {
0278       resid->set_doMicromegasOnly(true);
0279     }
0280   }
0281
0282   resid->clusterTree();
0283   resid->convertSeeds(G4TRACKING::convert_seeds_to_svtxtracks);
0284   resid->Verbosity(0);
0285   se->registerSubsystem(resid);
0286
0287   if (Enable::QA)
0288   {
0289     se->registerSubsystem(new SiliconSeedsQA);
0290     se->registerSubsystem(new TpcSeedsQA);
0291     se->registerSubsystem(new TpcSiliconQA);
0292   }
0293   se->run(nEvents);
0294   se->End();
0295   se->PrintTimer();
0296
0297   std::cout << "CDB Files used:" << std::endl;
0298   CDBInterface::instance()->Print();
0299
0300   if (Enable::QA)
0301   {
0302     TString qaname = theOutfile + "_qa.root";
0303     std::string qaOutputFileName(qaname.Data());
0304     QAHistManagerDef::saveQARootFile(qaOutputFileName);
0305   }
0306
0307   delete se;
0308   std::cout << "Finished" << std::endl;
0309   gSystem->Exit(0);
0310 }