macros/TrackingProduction/Fun4All_raw_hit_TPC_Matched_reco.C

0001 /*
0002  * This macro shows a minimum working example of running the
0003  * TPC track reconstruction starting with hit unpackers
0004  */
0005
0006 // leave the GlobalVariables.C at the beginning, an empty line afterwards
0007 // protects its position against reshuffling by clang-format
0008 #include <GlobalVariables.C>
0009
0010 #include <G4_ActsGeom.C>
0011 #include <G4_Global.C>
0012 #include <G4_Magnet.C>
0013 #include <G4_Mbd.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 <cdbobjects/CDBTTree.h>
0022
0023 #include <ffamodules/CDBInterface.h>
0024 #include <ffamodules/FlagHandler.h>
0025
0026 #include <fun4all/Fun4AllDstInputManager.h>
0027 #include <fun4all/Fun4AllDstOutputManager.h>
0028 #include <fun4all/Fun4AllInputManager.h>
0029 #include <fun4all/Fun4AllOutputManager.h>
0030 #include <fun4all/Fun4AllRunNodeInputManager.h>
0031 #include <fun4all/Fun4AllServer.h>
0032 #include <fun4all/Fun4AllUtils.h>
0033
0034 #include <fun4all/Fun4AllReturnCodes.h>
0035 #include <fun4all/SubsysReco.h>
0036
0037 #include <phool/recoConsts.h>
0038
0039 #include <tpctrackreco/TpcCrossingFinder.h>
0040 #include <tpctrackreco/TpcPolyClusterTrkrClusterConverter.h>
0041 #include <tpctrackreco/TpcPolyTrackSeedConverter.h>
0042 #include <tpctrackreco/Tpc_AssembledTrackReco.h>
0043 #include <tpctrackreco/Tpc_ModuleTrackReco.h>
0044 #include <tpctrackreco/Tpc_PolyClusterizer.h>
0045 #include <tpctrackreco/Tpc_PolyTrackReco.h>
0046 #include <tpctrackreco/Tpc_PolyTrackVertexer.h>
0047
0048 #include <trackingdiagnostics/Tpc_AssembledTrackDisplay.h>
0049 #include <trackingdiagnostics/Tpc_ModuleTrackDisplay.h>
0050 #include <trackingdiagnostics/Tpc_PolyClusterDisplay.h>
0051 #include <trackingdiagnostics/Tpc_PolyClusterResiduals.h>
0052 #include <trackingdiagnostics/TrackResiduals.h>
0053
0054 #include <format>
0055
0056 R__LOAD_LIBRARY(libfun4all.so)
0057 R__LOAD_LIBRARY(libffamodules.so)
0058 R__LOAD_LIBRARY(libphool.so)
0059 R__LOAD_LIBRARY(libcdbobjects.so)
0060 R__LOAD_LIBRARY(libmvtx.so)
0061 R__LOAD_LIBRARY(libintt.so)
0062 R__LOAD_LIBRARY(libtpc.so)
0063 R__LOAD_LIBRARY(libmicromegas.so)
0064 R__LOAD_LIBRARY(libPHGarfield.so)
0065 R__LOAD_LIBRARY(libtpctrackreco.so)
0066 R__LOAD_LIBRARY(libTrackingDiagnostics.so)
0067
0068 class SkipFirstN : public SubsysReco
0069 {
0070  public:
0071   explicit SkipFirstN(int n)
0072     : SubsysReco("SkipFirstN")
0073     , target_(n)
0074   {
0075   }
0076   int process_event(PHCompositeNode * /*unused*/) override
0077   {
0078     if (count_ < target_)
0079     {
0080       ++count_;
0081       return Fun4AllReturnCodes::ABORTEVENT;
0082     }
0083     return Fun4AllReturnCodes::EVENT_OK;
0084   }
0085
0086  private:
0087   int target_ = 0;
0088   int count_ = 0;
0089 };
0090 // 'Golden' store of p+p
0091 //.x Fun4All_raw_hit_TPC_reco.C(2, 79513, 0, ".", 0, "run3pp", "ana532_nocdbtag_v001","HITS_ppFieldOn")
0092 //.x Fun4All_raw_hit_TPC_reco.C(2, 79516, 0, ".", 0, "run3pp", "ana532_nocdbtag_v001","HITS_ppFieldOn")
0093
0094 // 6x6       75570-75573 - 1mrad; 75574 - 0mrad
0095 //.x Fun4All_raw_hit_TPC_reco.C(10, 75570, 0, ".", 0, "run3auau", "ana514_nocdbtag_v001","HITS_AuAu_ZeroField_1mrad")
0096 //.x Fun4All_raw_hit_TPC_reco.C(10, 75574, 0, ".", 0, "run3auau", "ana514_nocdbtag_v001","HITS_AuAu_ZeroField_0mrad")
0097
0098 // 111x111 Au+Au Field On; 75391 - close to end of store, 75405 - beginning of the next store
0099 //.x Fun4All_raw_hit_TPC_reco.C(2, 75391, 0, ".", 0, "run3auau", "ana514_nocdbtag_v001","HITS_AuAuFieldOn")
0100 //.x Fun4All_raw_hit_TPC_reco.C(2, 75405, 0, ".", 0, "run3auau", "ana514_nocdbtag_v001","HITS_AuAuFieldOn")
0101
0102 // 111x111 O+O
0103 //.x Fun4All_raw_hit_TPC_reco.C(2, 82626, 0, ".", 0, "run3oo", "ana537_nocdbtag_v001","HITS_OO")
0104
0105 void Fun4All_raw_hit_TPC_Matched_reco(
0106     const int nEvents = 10,
0107     const int runnumber = 79513,
0108     const int segment = 0,
0109     const std::string &outdir = ".",
0110     const int nSkip = 0,
0111     const std::string &collision = "run3pp",
0112     const std::string &production = "ana532_nocdbtag_v001",
0113     const std::string &outfilename = "ppFieldOn",
0114     const std::string &datatype = "physics")
0115 {
0116   const bool convertSeeds = false;
0117   auto *se = Fun4AllServer::instance();
0118   se->Verbosity(1);
0119   auto *rc = recoConsts::instance();
0120   se->registerSubsystem(new SkipFirstN(nSkip));
0121
0122   const std::string dsttype = "STREAMING_EVENT";
0123   const std::string dsttype_to_save = "TPC";
0124
0125   G4TPC::sampa_tzero_bias = 0;
0126
0127   // First order corrections will be applied from PHGarfield
0128   G4TPC::ENABLE_MODULE_EDGE_CORRECTIONS = false;
0129   G4TPC::ENABLE_STATIC_CORRECTIONS = false;
0130   G4TPC::ENABLE_AVERAGE_CORRECTIONS = false;
0131
0132   std::ostringstream runstr;
0133   runstr << std::setw(8) << std::setfill('0') << runnumber;
0134
0135   std::ostringstream segstr;
0136   segstr << std::setw(5) << std::setfill('0') << segment;
0137
0138   const int runbase = (runnumber / 100) * 100;
0139   const int runnext = runbase + 100;
0140
0141   std::ostringstream rundir;
0142   rundir << "run_" << std::setw(8) << std::setfill('0') << runbase << "_" << std::setw(8) << std::setfill('0') << runnext;
0143
0144   std::vector<std::string> streams;
0145
0146   // TPC EBDC streams
0147   for (int ebdc = 0; ebdc < 24; ++ebdc)
0148   {
0149     for (int endpoint = 0; endpoint < 2; ++endpoint)
0150     {
0151       std::ostringstream s;
0152       s << "ebdc" << std::setw(2) << std::setfill('0') << ebdc << "_" << endpoint;
0153       streams.push_back(s.str());
0154     }
0155   }
0156
0157   // TPOT
0158   streams.emplace_back("ebdc39");
0159
0160   // INTT streams
0161   for (int server = 0; server < 8; ++server)
0162   {
0163     std::ostringstream s;
0164     s << "intt" << server;
0165     streams.push_back(s.str());
0166   }
0167
0168   // MVTX streams
0169   for (int felix = 0; felix < 6; ++felix)
0170   {
0171     std::ostringstream s;
0172     s << "mvtx" << felix;
0173     streams.push_back(s.str());
0174   }
0175
0176   int i = 0;
0177   std::stringstream nice_runnumber;
0178   nice_runnumber << std::setw(8) << std::setfill('0') << std::to_string(runnumber);
0179
0180   int rounded_up = 100 * (std::ceil((float) runnumber / 100));
0181   std::stringstream nice_rounded_up;
0182   nice_rounded_up << std::setw(8) << std::setfill('0') << std::to_string(rounded_up);
0183
0184   int rounded_down = 100 * (std::floor((float) runnumber / 100));
0185   std::stringstream nice_rounded_down;
0186   nice_rounded_down << std::setw(8) << std::setfill('0') << std::to_string(rounded_down);
0187
0188   for (const auto& stream : streams)
0189   {
0190     std::string filename = std::format(
0191     "DST_{}_{}_{}_{}-{}-{}.root",
0192     dsttype, stream, collision, production, runstr.str(), segstr.str());
0193
0194     std::string filepath = std::format(
0195     "/sphenix/lustre01/sphnxpro/production/{}/{}/{}/DST_{}_{}/run_{}_{}/{}",
0196     collision, datatype, production, dsttype, stream,
0197     nice_rounded_down.str(), nice_rounded_up.str(), filename);
0198     std::cout << "Adding DST: " << filepath << std::endl;
0199     if (i == 0)
0200     {
0201       rc->set_IntFlag("RUNNUMBER", runnumber);
0202       rc->set_uint64Flag("TIMESTAMP", runnumber);
0203     }
0204
0205     std::string inputname = "InputManager" + std::to_string(i);
0206     auto *hitsin = new Fun4AllDstInputManager(inputname);
0207     hitsin->fileopen(filename);
0208     se->registerInputManager(hitsin);
0209
0210     ++i;
0211   }
0212
0213   rc->set_IntFlag("RUNNUMBER", runnumber);
0214   rc->set_IntFlag("RUNSEGMENT", segment);
0215
0216   Enable::QA = false;
0217   Enable::CDB = true;
0218   rc->set_StringFlag("CDB_GLOBALTAG", "newcdbtag");
0219   rc->set_uint64Flag("TIMESTAMP", runnumber);
0220
0221   G4TRACKING::convert_seeds_to_svtxtracks = convertSeeds;
0222
0223   G4MAGNET::magfield_rescale = 1;
0224
0225   std::cout << " run: " << runnumber
0226             << " samples: " << TRACKING::reco_tpc_maxtime_sample
0227             << " pre: " << TRACKING::reco_tpc_time_presample
0228             << " vdrift: " << G4TPC::tpc_drift_velocity_reco
0229             << std::endl;
0230
0231   TRACKING::streaming_mode = true;
0232
0233   FlagHandler *flag = new FlagHandler();
0234   se->registerSubsystem(flag);
0235
0236   std::string geofile = CDBInterface::instance()->getUrl("Tracking_Geometry");
0237
0238   Fun4AllRunNodeInputManager *ingeo = new Fun4AllRunNodeInputManager("GeoIn");
0239   ingeo->AddFile(geofile);
0240   se->registerInputManager(ingeo);
0241
0242   TpcReadoutInit(runnumber);
0243   G4TPC::REJECT_LASER_EVENTS = true;
0244   // Flag for running the tpc hit unpacker with zero suppression on
0245   TRACKING::tpc_zero_supp = true;
0246
0247   TrackingInit();
0248
0249   for (int felix = 0; felix < 6; felix++)
0250   {
0251     Mvtx_HitUnpacking(std::to_string(felix));
0252   }
0253   for (int server = 0; server < 8; server++)
0254   {
0255     Intt_HitUnpacking(std::to_string(server));
0256   }
0257   std::ostringstream ebdcname;
0258   for (int ebdc = 0; ebdc < 24; ebdc++)
0259   {
0260     for (int endpoint = 0; endpoint < 2; endpoint++)
0261     {
0262       ebdcname.str("");
0263       if (ebdc < 10)
0264       {
0265         ebdcname << "0";
0266       }
0267       ebdcname << ebdc << "_" << endpoint;
0268       Tpc_HitUnpacking(ebdcname.str());
0269     }
0270   }
0271   Micromegas_HitUnpacking();
0272
0273   //==============================================================
0274
0275   Mvtx_Clustering();
0276   Intt_Clustering();
0277   Micromegas_Clustering();
0278
0279   //==============================================================
0280
0281   Tpc_LaserEventIdentifying();
0282   Reject_Laser_Events();
0283
0284   //==============================================================
0285
0286   Tracking_Reco_SiliconSeed_run2pp();
0287   auto *converter = new TrackSeedTrackMapConverter("SiliconSeedToSvtxTrackMap");
0288   converter->setTrackSeedName("SiliconTrackSeedContainer");
0289   converter->setTrackMapName("SiliconSvtxTrackMap");
0290   converter->setClusterMapName("TRKR_CLUSTER");
0291   se->registerSubsystem(converter);
0292
0293   auto *finder_svx = new PHSimpleVertexFinder("SiliconVertexFinder");
0294   finder_svx->Verbosity(0);
0295   finder_svx->setDcaCut(0.1);
0296   finder_svx->setTrackPtCut(0.2);
0297   finder_svx->setBeamLineCut(1);
0298   finder_svx->setTrackQualityCut(500);
0299   finder_svx->setNmvtxRequired(3);
0300   finder_svx->setOutlierPairCut(0.1);
0301   finder_svx->setTrackMapName("SiliconSvtxTrackMap");
0302   finder_svx->setVertexMapName("SiliconSvtxVertexMap");
0303   se->registerSubsystem(finder_svx);
0304
0305   //==============================================================
0306
0307   se->registerSubsystem(new Tpc_ModuleTrackReco());     // makes TPC_MODULETRACKS
0308   se->registerSubsystem(new Tpc_AssembledTrackReco());  // makes TPC_ASSEMBLEDTRACKS
0309
0310   auto *crossingFinder = new TpcCrossingFinder();
0311   crossingFinder->Verbosity(0);
0312   crossingFinder->setInputNodeName("TPC_ASSEMBLEDTRACKS");
0313   crossingFinder->setOutputNodeName("TPC_CROSSING_DECISIONS");
0314   crossingFinder->setVertexMapNodeName("SiliconSvtxVertexMap");  // optional, configurable
0315   se->registerSubsystem(crossingFinder);
0316
0317   auto *cluster = new Tpc_PolyClusterizer();  // makes TPC_POLYCLUSTERS
0318   cluster->setUseSurveyGeometry(false);
0319   cluster->setKEffSide0(1.00);  // OO 82626 - 4.5, AuAu 6x6 76905 -0, pp 79513 - 1.0, 75391 5.8 75405 4.8
0320   cluster->setKEffSide1(1.60);  // OO 82626 - 5.0, AuAu 6x6 76905 -0, pp 79513 - 1.6, 75391 5.6 75408 4.8
0321   cluster->setMaxAcceptedTier(2);
0322   se->registerSubsystem(cluster);
0323
0324   se->registerSubsystem(new Tpc_PolyTrackReco());      // makes TPC_POLYTRACKS
0325   se->registerSubsystem(new Tpc_PolyTrackVertexer());  // makes TPC_POLYTRACKVERTICES
0326
0327   se->registerSubsystem(new TpcPolyTrackSeedConverter());           // converts TPC_POLYTRACKS to TpcTrackSeed
0328   se->registerSubsystem(new TpcPolyClusterTrkrClusterConverter());  // converts TPC_POLYCLUSTERS to TRKR_CLUSTER
0329
0330   auto *silicon_match = new PHSiliconTpcTrackMatching;
0331   silicon_match->Verbosity(0);
0332   silicon_match->set_pp_mode(TRACKING::streaming_mode);
0333   silicon_match->set_cluster_map_name("TRKR_CLUSTER");
0334   silicon_match->window_deta.set_posQoverpT_maxabs({-0.014, 0.0331, 0.48});
0335   silicon_match->window_deta.set_negQoverpT_maxabs({-0.014, 0.0331, 0.48});
0336   silicon_match->set_deltaeta_min(0.03);
0337   silicon_match->window_dphi.set_QoverpT_range({-0.15, 0, 0}, {0.15, 0, 0});
0338   silicon_match->window_dx.set_QoverpT_maxabs({3.0, 0, 0});
0339   silicon_match->window_dy.set_QoverpT_maxabs({3.0, 0, 0});
0340   silicon_match->window_dz.set_posQoverpT_maxabs({1.138, 0.3919, 0.84});
0341   silicon_match->window_dz.set_negQoverpT_maxabs({1.138, 0.3919, 0.84});
0342   silicon_match->set_crossing_deltaz_max(10);
0343   silicon_match->set_crossing_deltaz_min(0);
0344   silicon_match->set_test_windows_printout(false);
0345   silicon_match->set_max_crossing_diff(10); // good for poly seeding case
0346   // these are for testing. and default to false. The seed matcher will choose the crossing that works best
0347   //  silicon_match->set_use_tpc_crossing_only(false);  // use crossing information from TPC SA seed
0348   //  silicon_match->set_use_silicon_crossing_only(false);  // use crossing information from silicon seed
0349   se->registerSubsystem(silicon_match);
0350
0351   auto *deltazcorr = new PHTpcDeltaZCorrection;
0352   deltazcorr->Verbosity(0);
0353   se->registerSubsystem(deltazcorr);
0354
0355   auto *actsFit = new PHActsTrkFitter;
0356   actsFit->Verbosity(0);
0357   actsFit->commissioning(G4TRACKING::use_alignment);
0358   actsFit->setTrkrClusterContainerName("TRKR_CLUSTER");
0359   // in calibration mode, fit only Silicons and Micromegas hits
0360   actsFit->fitSiliconMMs(G4TRACKING::SC_CALIBMODE);
0361   actsFit->set_pp_mode(TRACKING::streaming_mode);
0362   actsFit->setUseMicromegas(false);
0363   actsFit->set_use_clustermover(false);  // default is true for now
0364   actsFit->useActsEvaluator(false);
0365   actsFit->useOutlierFinder(false);
0366   actsFit->setFieldMap(G4MAGNET::magfield_tracking);
0367   se->registerSubsystem(actsFit);
0368
0369   auto *cleaner = new PHTrackCleaner();
0370   cleaner->Verbosity(0);
0371   cleaner->set_quality_cut(1000);
0372   cleaner->set_pp_mode(TRACKING::streaming_mode);
0373   se->registerSubsystem(cleaner);
0374
0375   auto *finder = new PHSimpleVertexFinder("SvtxVertexFinder");
0376   finder->Verbosity(0);
0377   finder->setDcaCut(0.1);
0378   finder->setTrackPtCut(0.2);
0379   finder->setBeamLineCut(1);
0380   finder->setTrackQualityCut(500);
0381   finder->setNmvtxRequired(3);
0382   finder->setOutlierPairCut(0.1);
0383   finder->setTrackMapName("SvtxTrackMap");
0384   finder->setVertexMapName("SvtxVertexMap");
0385   se->registerSubsystem(finder);
0386
0387   auto *resid = new TrackResiduals("TrackResiduals");
0388   resid->outfileName(outdir + "/TrackResiduals_" + outfilename + "_" + std::to_string(runnumber) + "_" + std::to_string(segment) + ".root");
0389   resid->alignment(false);
0390   resid->clusterTree();
0391   resid->vertexTree();
0392   resid->hitTree();
0393   resid->convertSeeds(G4TRACKING::convert_seeds_to_svtxtracks);
0394   resid->Verbosity(0);
0395   se->registerSubsystem(resid);
0396
0397   Fun4AllOutputManager *out = new Fun4AllDstOutputManager("out", std::format("{}/output_DST/DST_{}_{}_{}-{}-{}.root",
0398             outdir,
0399             dsttype_to_save,
0400             collision,
0401             production,
0402             runnumber,
0403             segment));
0404
0405   out->AddNode("Sync");
0406   out->AddNode("EventHeader");
0407   out->AddRunNode("TPCGEOMCONTAINER");
0408   // out->AddNode("TRKR_HITSET");
0409   // out->AddNode("TPC_MODULETRACKS");
0410   // out->AddNode("TPC_ASSEMBLEDTRACKS");
0411   out->AddNode("TPC_CROSSING_DECISIONS");
0412   out->AddNode("TPC_POLYCLUSTERS");
0413   out->AddNode("TPC_POLYTRACKS");
0414   out->AddNode("TPC_POLYTRACKVERTICES");
0415   out->AddNode("TRKR_CLUSTER");
0416
0417   // se->registerOutputManager(out);
0418
0419   se->run(nEvents + nSkip);
0420   se->Print("NODETREE");
0421   se->End();
0422   se->PrintTimer();
0423
0424   CDBInterface::instance()->Print();
0425   delete se;
0426   std::cout << "Finished" << std::endl;
0427   gSystem->Exit(0);
0428 }