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