macros/TrackingProduction/Fun4All_PRDFReconstruction.C

0001 /*
0002  * This macro shows running the full event combining + tracking for
0003  * cosmics running.. There are some analysis modules run at the end
0004  * which package  hits, clusters, and clusters on tracks into trees
0005  * for analysis.
0006  */
0007 #include <GlobalVariables.C>
0008
0009 #include "../detectors/sPHENIX/G4Setup_sPHENIX.C"
0010
0011 #include <G4_ActsGeom.C>
0012 #include <G4_Magnet.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 <trackreco/PHActsSiliconSeeding.h>
0020 #include <trackreco/PHSiliconSeedMerger.h>
0021
0022 #include <mvtxrawhitqa/MvtxRawHitQA.h>
0023
0024 #include <inttrawhitqa/InttRawHitQA.h>
0025
0026 #include <tpcqa/TpcRawHitQA.h>
0027
0028 #include <trackingqa/InttClusterQA.h>
0029 #include <trackingqa/MicromegasClusterQA.h>
0030 #include <trackingqa/MvtxClusterQA.h>
0031 #include <trackingqa/SiliconSeedsQA.h>
0032 #include <trackingqa/TpcClusterQA.h>
0033 #include <trackingqa/TpcSeedsQA.h>
0034 #include <trackingqa/TpcSiliconQA.h>
0035 #include <trackingqa/TrackFittingQA.h>
0036
0037 #include <trackingdiagnostics/TrackResiduals.h>
0038 #include <trackingdiagnostics/TrkrNtuplizer.h>
0039
0040 #include <ffamodules/CDBInterface.h>
0041 #include <ffamodules/FlagHandler.h>
0042 #include <ffamodules/HeadReco.h>
0043 #include <ffamodules/SyncReco.h>
0044
0045 #include <fun4allraw/Fun4AllStreamingInputManager.h>
0046 #include <fun4allraw/InputManagerType.h>
0047 #include <fun4allraw/SingleGl1PoolInput.h>
0048 #include <fun4allraw/SingleInttPoolInput.h>
0049 #include <fun4allraw/SingleMicromegasPoolInput.h>
0050 #include <fun4allraw/SingleMvtxPoolInput.h>
0051 #include <fun4allraw/SingleTpcTimeFrameInput.h>
0052
0053 #include <fun4all/Fun4AllDstInputManager.h>
0054 #include <fun4all/Fun4AllDstOutputManager.h>
0055 #include <fun4all/Fun4AllInputManager.h>
0056 #include <fun4all/Fun4AllOutputManager.h>
0057 #include <fun4all/Fun4AllRunNodeInputManager.h>
0058 #include <fun4all/Fun4AllServer.h>
0059 #include <fun4all/Fun4AllUtils.h>
0060
0061 #include <phool/recoConsts.h>
0062
0063 R__LOAD_LIBRARY(libfun4all.so)
0064 R__LOAD_LIBRARY(libffamodules.so)
0065 R__LOAD_LIBRARY(libphool.so)
0066 R__LOAD_LIBRARY(libcdbobjects.so)
0067 R__LOAD_LIBRARY(libmvtx.so)
0068 R__LOAD_LIBRARY(libintt.so)
0069 R__LOAD_LIBRARY(libtpc.so)
0070 R__LOAD_LIBRARY(libmicromegas.so)
0071 R__LOAD_LIBRARY(libTrackingDiagnostics.so)
0072 R__LOAD_LIBRARY(libtrackingqa.so)
0073 R__LOAD_LIBRARY(libtrack_reco.so)
0074
0075 bool isGood(const std::string &infile);
0076 int getrunnumber(const std::string &listfile);
0077
0078 void Fun4All_PRDFReconstruction(
0079     const int nEvents = 0,
0080     const int /*runnumber_unused*/ = 53756,
0081     const std::string & /*dir*/ = "/.",
0082     const std::string &outfilename = "output_tracks",
0083     const std::string &input_gl1file = "gl1daq.list",
0084     const std::string &input_inttfile00 = "intt0.list",
0085     const std::string &input_inttfile01 = "intt1.list",
0086     const std::string &input_inttfile02 = "intt2.list",
0087     const std::string &input_inttfile03 = "intt3.list",
0088     const std::string &input_inttfile04 = "intt4.list",
0089     const std::string &input_inttfile05 = "intt5.list",
0090     const std::string &input_inttfile06 = "intt6.list",
0091     const std::string &input_inttfile07 = "intt7.list",
0092     const std::string &input_mvtxfile00 = "mvtx0.list",
0093     const std::string &input_mvtxfile01 = "mvtx1.list",
0094     const std::string &input_mvtxfile02 = "mvtx2.list",
0095     const std::string &input_mvtxfile03 = "mvtx3.list",
0096     const std::string &input_mvtxfile04 = "mvtx4.list",
0097     const std::string &input_mvtxfile05 = "mvtx5.list",
0098     const std::string &input_tpcfile00 = "tpc00_0.list",
0099     const std::string &input_tpcfile01 = "tpc01_0.list",
0100     const std::string &input_tpcfile02 = "tpc02_0.list",
0101     const std::string &input_tpcfile03 = "tpc03_0.list",
0102     const std::string &input_tpcfile04 = "tpc04_0.list",
0103     const std::string &input_tpcfile05 = "tpc05_0.list",
0104     const std::string &input_tpcfile06 = "tpc06_0.list",
0105     const std::string &input_tpcfile07 = "tpc07_0.list",
0106     const std::string &input_tpcfile08 = "tpc08_0.list",
0107     const std::string &input_tpcfile09 = "tpc09_0.list",
0108     const std::string &input_tpcfile10 = "tpc10_0.list",
0109     const std::string &input_tpcfile11 = "tpc11_0.list",
0110     const std::string &input_tpcfile12 = "tpc12_0.list",
0111     const std::string &input_tpcfile13 = "tpc13_0.list",
0112     const std::string &input_tpcfile14 = "tpc14_0.list",
0113     const std::string &input_tpcfile15 = "tpc15_0.list",
0114     const std::string &input_tpcfile16 = "tpc16_0.list",
0115     const std::string &input_tpcfile17 = "tpc17_0.list",
0116     const std::string &input_tpcfile18 = "tpc18_0.list",
0117     const std::string &input_tpcfile19 = "tpc19_0.list",
0118     const std::string &input_tpcfile20 = "tpc20_0.list",
0119     const std::string &input_tpcfile21 = "tpc21_0.list",
0120     const std::string &input_tpcfile22 = "tpc22_0.list",
0121     const std::string &input_tpcfile23 = "tpc23_0.list",
0122     const std::string &input_tpotfile = "tpot.list")
0123 {
0124   std::vector<std::string> gl1_infile;
0125   gl1_infile.push_back(input_gl1file);
0126
0127   // MVTX
0128   std::vector<std::string> mvtx_infile;
0129   mvtx_infile.push_back(input_mvtxfile00);
0130   mvtx_infile.push_back(input_mvtxfile01);
0131   mvtx_infile.push_back(input_mvtxfile02);
0132   mvtx_infile.push_back(input_mvtxfile03);
0133   mvtx_infile.push_back(input_mvtxfile04);
0134   mvtx_infile.push_back(input_mvtxfile05);
0135
0136   // INTT
0137   std::vector<std::string> intt_infile;
0138   intt_infile.push_back(input_inttfile00);
0139   intt_infile.push_back(input_inttfile01);
0140   intt_infile.push_back(input_inttfile02);
0141   intt_infile.push_back(input_inttfile03);
0142   intt_infile.push_back(input_inttfile04);
0143   intt_infile.push_back(input_inttfile05);
0144   intt_infile.push_back(input_inttfile06);
0145   intt_infile.push_back(input_inttfile07);
0146
0147   std::vector<std::string> tpc_infile;
0148   tpc_infile.push_back(input_tpcfile00);
0149   tpc_infile.push_back(input_tpcfile01);
0150   tpc_infile.push_back(input_tpcfile02);
0151   tpc_infile.push_back(input_tpcfile03);
0152   tpc_infile.push_back(input_tpcfile04);
0153   tpc_infile.push_back(input_tpcfile05);
0154   tpc_infile.push_back(input_tpcfile06);
0155   tpc_infile.push_back(input_tpcfile07);
0156   tpc_infile.push_back(input_tpcfile08);
0157   tpc_infile.push_back(input_tpcfile09);
0158   tpc_infile.push_back(input_tpcfile10);
0159   tpc_infile.push_back(input_tpcfile11);
0160   tpc_infile.push_back(input_tpcfile12);
0161   tpc_infile.push_back(input_tpcfile13);
0162   tpc_infile.push_back(input_tpcfile14);
0163   tpc_infile.push_back(input_tpcfile15);
0164   tpc_infile.push_back(input_tpcfile16);
0165   tpc_infile.push_back(input_tpcfile17);
0166   tpc_infile.push_back(input_tpcfile18);
0167   tpc_infile.push_back(input_tpcfile19);
0168   tpc_infile.push_back(input_tpcfile20);
0169   tpc_infile.push_back(input_tpcfile21);
0170   tpc_infile.push_back(input_tpcfile22);
0171   tpc_infile.push_back(input_tpcfile23);
0172
0173   // TPOT
0174   std::vector<std::string> tpot_infile;
0175   tpot_infile.push_back(input_tpotfile);
0176
0177   int runnumber = -99999;
0178   if (!gl1_infile.empty())
0179   {
0180     runnumber = getrunnumber(gl1_infile[0]);
0181   }
0182   else if (!mvtx_infile.empty())
0183   {
0184     runnumber = getrunnumber(mvtx_infile[0]);
0185   }
0186   else if (!intt_infile.empty())
0187   {
0188     runnumber = getrunnumber(intt_infile[0]);
0189   }
0190   else if (!tpc_infile.empty())
0191   {
0192     runnumber = getrunnumber(tpc_infile[0]);
0193   }
0194   else if (!tpot_infile.empty())
0195   {
0196     runnumber = getrunnumber(tpot_infile[0]);
0197   }
0198   if (runnumber == -99999)
0199   {
0200     std::cout << "could not extract run number from input files (all lists empty?)"
0201               << std::endl;
0202     gSystem->Exit(1);
0203   }
0204
0205   auto *se = Fun4AllServer::instance();
0206   se->Verbosity(1);
0207   auto *rc = recoConsts::instance();
0208   rc->set_IntFlag("RUNNUMBER", runnumber);
0209
0210   Enable::CDB = true;
0211   rc->set_StringFlag("CDB_GLOBALTAG", "newcdbtag");
0212
0213   rc->set_uint64Flag("TIMESTAMP", runnumber);
0214
0215   //! flags to set
0216   Enable::QA = true;
0217   TRACKING::tpc_zero_supp = true;
0218   TRACKING::pp_mode = true;
0219   G4TRACKING::convert_seeds_to_svtxtracks = false;
0220   G4TPC::REJECT_LASER_EVENTS = true;
0221
0222   Enable::MVTX_APPLYMISALIGNMENT = true;
0223   ACTSGEOM::mvtx_applymisalignment = Enable::MVTX_APPLYMISALIGNMENT;
0224   TpcReadoutInit(runnumber);
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   CDBInterface::instance()->Verbosity(1);
0232   // std::string geofile = CDBInterface::instance()->getUrl("Tracking_Geometry");
0233   // std::cout << "CDB tracking geometry file "<<geofile << std::endl;
0234   // Fun4AllRunNodeInputManager *ingeo = new Fun4AllRunNodeInputManager("GeoIn");
0235   // ingeo->AddFile(geofile);
0236   // se->registerInputManager(ingeo);
0237
0238   Enable::MVTX = true;
0239   Enable::INTT = true;
0240   Enable::TPC = true;
0241   Enable::MICROMEGAS = true;
0242
0243   G4Init();
0244   G4Setup();
0245
0246   // can use for zero field
0247   // double fieldstrength = 0.01;
0248   // G4MAGNET::magfield_tracking = "0.01";
0249   double fieldstrength = std::numeric_limits<double>::quiet_NaN();
0250   bool ConstField = isConstantField(G4MAGNET::magfield_tracking, fieldstrength);
0251
0252   if (ConstField && fieldstrength < 0.1)
0253   {
0254     G4MAGNET::magfield = "0.01";
0255     G4MAGNET::magfield_rescale = 1;
0256   }
0257
0258   TrackingInit();
0259
0260   int i = 0;
0261   int NumInputs = 0;
0262   Fun4AllStreamingInputManager *in = new Fun4AllStreamingInputManager("Comb");
0263
0264   for (const auto &iter : gl1_infile)
0265   {
0266     if (isGood(iter))
0267     {
0268       SingleGl1PoolInput *gl1_sngl = new SingleGl1PoolInput("GL1_" + std::to_string(i));
0269       //    gl1_sngl->Verbosity(3);
0270       gl1_sngl->AddListFile(iter);
0271       in->registerStreamingInput(gl1_sngl, InputManagerType::GL1);
0272       i++;
0273     }
0274   }
0275   NumInputs += i;
0276
0277   i = 0;
0278   for (const auto &iter : intt_infile)
0279   {
0280     if (isGood(iter))
0281     {
0282       SingleInttPoolInput *intt_sngl = new SingleInttPoolInput("INTT_" + std::to_string(i));
0283       intt_sngl->AddListFile(iter);
0284       in->registerStreamingInput(intt_sngl, InputManagerType::INTT);
0285       i++;
0286     }
0287   }
0288   NumInputs += i;
0289
0290   i = 0;
0291   for (const auto &iter : mvtx_infile)
0292   {
0293     if (isGood(iter))
0294     {
0295       SingleMvtxPoolInput *mvtx_sngl = new SingleMvtxPoolInput("MVTX_" + std::to_string(i));
0296       //    mvtx_sngl->Verbosity(3);
0297       mvtx_sngl->AddListFile(iter);
0298       in->registerStreamingInput(mvtx_sngl, InputManagerType::MVTX);
0299       i++;
0300     }
0301   }
0302   NumInputs += i;
0303
0304   i = 0;
0305
0306   for (const auto &iter : tpc_infile)
0307   {
0308     if (isGood(iter))
0309     {
0310       SingleTpcTimeFrameInput *tpc_sngl = new SingleTpcTimeFrameInput("TPC_" + std::to_string(i));
0311       //    tpc_sngl->Verbosity(2);
0312       //   tpc_sngl->DryRun();
0313       tpc_sngl->setHitContainerName("TPCRAWHIT");
0314       tpc_sngl->AddListFile(iter);
0315       in->registerStreamingInput(tpc_sngl, InputManagerType::TPC);
0316       i++;
0317     }
0318   }
0319   NumInputs += i;
0320   i = 0;
0321   for (const auto &iter : tpot_infile)
0322   {
0323     if (isGood(iter))
0324     {
0325       SingleMicromegasPoolInput *mm_sngl = new SingleMicromegasPoolInput("MICROMEGAS_" + std::to_string(i));
0326       //   sngl->Verbosity(3);
0327       mm_sngl->SetBcoRange(10);
0328       mm_sngl->SetNegativeBco(2);
0329       mm_sngl->SetBcoPoolSize(50);
0330       mm_sngl->AddListFile(iter);
0331       in->registerStreamingInput(mm_sngl, InputManagerType::MICROMEGAS);
0332       i++;
0333     }
0334   }
0335   NumInputs += i;
0336
0337   // if there is no input manager this macro will still run - so just quit here
0338   if (NumInputs == 0)
0339   {
0340     std::cout << "no file lists no input manager registered, quitting" << std::endl;
0341     gSystem->Exit(1);
0342   }
0343   se->registerInputManager(in);
0344
0345   SyncReco *sync = new SyncReco();
0346   se->registerSubsystem(sync);
0347
0348   HeadReco *head = new HeadReco();
0349   se->registerSubsystem(head);
0350
0351   FlagHandler *flag = new FlagHandler();
0352   se->registerSubsystem(flag);
0353
0354   Mvtx_HitUnpacking();
0355   Intt_HitUnpacking();
0356   Tpc_HitUnpacking();
0357   Micromegas_HitUnpacking();
0358
0359   Mvtx_Clustering();
0360
0361   Intt_Clustering();
0362
0363   Tpc_LaserEventIdentifying();
0364
0365   TPC_Clustering_run2pp();
0366   Micromegas_Clustering();
0367
0368   Reject_Laser_Events();
0369
0370   Tracking_Reco_TrackSeed_run2pp();
0371   Tracking_Reco_TrackMatching_run2pp();
0372
0373   Tracking_Reco_TrackFit_run2pp();
0374   // vertexing and propagation to vertex
0375   Tracking_Reco_Vertex_run2pp();
0376
0377   std::string residstring = outfilename + "_resid.root";
0378
0379   auto *resid = new TrackResiduals("TrackResiduals");
0380   resid->Verbosity(0);
0381   resid->outfileName(residstring);
0382   resid->alignment(false);
0383   // resid->hitTree();
0384   resid->convertSeeds(G4TRACKING::convert_seeds_to_svtxtracks);
0385
0386   if (ConstField && fieldstrength < 0.1)
0387   {
0388     resid->zeroField();
0389   }
0390   // se->registerSubsystem(resid);
0391
0392   // Fun4AllOutputManager *out = new Fun4AllDstOutputManager("out", "/sphenix/tg/tg01/hf/jdosbo/tracking_development/onlineoffline/hitsets.root");
0393   // se->registerOutputManager(out);
0394
0395   if (Enable::QA)
0396   {
0397     se->registerSubsystem(new MvtxRawHitQA);
0398     se->registerSubsystem(new InttRawHitQA);
0399     se->registerSubsystem(new TpcRawHitQA);
0400     se->registerSubsystem(new MvtxClusterQA);
0401     se->registerSubsystem(new InttClusterQA);
0402     se->registerSubsystem(new TpcClusterQA);
0403     se->registerSubsystem(new MicromegasClusterQA);
0404     se->registerSubsystem(new SiliconSeedsQA);
0405     se->registerSubsystem(new TpcSeedsQA);
0406     se->registerSubsystem(new TpcSiliconQA);
0407     se->registerSubsystem(new TrackFittingQA);
0408   }
0409
0410   se->run(nEvents);
0411   se->End();
0412   se->PrintTimer();
0413
0414   if (Enable::QA)
0415   {
0416     std::string qaOutputFileName = outfilename + std::to_string(runnumber) + "_qa.root";
0417     QAHistManagerDef::saveQARootFile(qaOutputFileName);
0418   }
0419   CDBInterface::instance()->Print();
0420   delete se;
0421   std::cout << "Finished" << std::endl;
0422   gSystem->Exit(0);
0423 }
0424
0425 bool isGood(const std::string &infile)
0426 {
0427   std::ifstream intest;
0428   intest.open(infile);
0429   bool goodfile = false;
0430   if (intest.is_open())
0431   {
0432     if (intest.peek() != std::ifstream::traits_type::eof())  // is it non zero?
0433     {
0434       goodfile = true;
0435     }
0436     intest.close();
0437   }
0438   return goodfile;
0439 }
0440
0441 int getrunnumber(const std::string &listfile)
0442 {
0443   if (!isGood(listfile))
0444   {
0445     std::cout << "listfile " << listfile << " is bad" << std::endl;
0446     gSystem->Exit(1);
0447   }
0448   std::ifstream ifs(listfile);
0449   std::string filepath;
0450   std::getline(ifs, filepath);
0451
0452   std::pair<int, int> runseg = Fun4AllUtils::GetRunSegment(filepath);
0453   int runnumber = runseg.first;
0454   //  int segment = abs(runseg.second);
0455   return runnumber;
0456 }