macros/TrackingProduction/Fun4All_Cosmics.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_Magnet.C>
0014 #include <QA.C>
0015 #include <Trkr_Clustering.C>
0016 #include <Trkr_RecoInit.C>
0017 #include <Trkr_Reco_Cosmics.C>
0018 #include <Trkr_TpcReadoutInit.C>
0019
0020 #include <trackingqa/InttClusterQA.h>
0021 #include <trackingqa/MicromegasClusterQA.h>
0022 #include <trackingqa/MvtxClusterQA.h>
0023 #include <trackingqa/TpcClusterQA.h>
0024 #include <trackingqa/TpcSeedsQA.h>
0025
0026 #include <trackingdiagnostics/TrackResiduals.h>
0027 #include <trackingdiagnostics/TrkrNtuplizer.h>
0028
0029 #include <cdbobjects/CDBTTree.h>
0030
0031 #include <ffamodules/CDBInterface.h>
0032
0033 #include <fun4all/Fun4AllDstInputManager.h>
0034 #include <fun4all/Fun4AllDstOutputManager.h>
0035 #include <fun4all/Fun4AllInputManager.h>
0036 #include <fun4all/Fun4AllOutputManager.h>
0037 #include <fun4all/Fun4AllRunNodeInputManager.h>
0038 #include <fun4all/Fun4AllServer.h>
0039 #include <fun4all/Fun4AllUtils.h>
0040
0041 #include <phool/recoConsts.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(libmvtx.so)
0048 R__LOAD_LIBRARY(libintt.so)
0049 R__LOAD_LIBRARY(libtpc.so)
0050 R__LOAD_LIBRARY(libmicromegas.so)
0051 R__LOAD_LIBRARY(libTrackingDiagnostics.so)
0052 R__LOAD_LIBRARY(libtrackingqa.so)
0053
0054 void Fun4All_Cosmics(
0055     const int nEvents = 0,
0056     const std::string &filelist = "filelist.list",
0057     const std::string &dir = "./",
0058     const std::string &outfilename = "cosmics")
0059 {
0060   TRACKING::tpc_zero_supp = true;
0061   TRACKING::tpc_baseline_corr = true;
0062   Enable::MVTX_APPLYMISALIGNMENT = true;
0063   ACTSGEOM::mvtx_applymisalignment = Enable::MVTX_APPLYMISALIGNMENT;
0064   Enable::QA = false;
0065
0066   auto *se = Fun4AllServer::instance();
0067   se->Verbosity(1);
0068   se->VerbosityDownscale(1000);
0069   auto *rc = recoConsts::instance();
0070
0071   rc->set_StringFlag("CDB_GLOBALTAG", "ProdA_2024");
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   int nTpcFiles = 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
0091     if (filepath.find("ebdc") != std::string::npos)
0092     {
0093       if (filepath.find("ebdc39") == std::string::npos)
0094       {
0095         nTpcFiles++;
0096       }
0097     }
0098     auto *hitsin = new Fun4AllDstInputManager(inputname);
0099     hitsin->fileopen(filepath);
0100     se->registerInputManager(hitsin);
0101     i++;
0102   }
0103
0104   TpcReadoutInit(runnumber);
0105   std::cout << " run: " << runnumber
0106             << " samples: " << TRACKING::reco_tpc_maxtime_sample
0107             << " pre: " << TRACKING::reco_tpc_time_presample
0108             << " vdrift: " << G4TPC::tpc_drift_velocity_reco
0109             << std::endl;
0110
0111   Enable::CDB = true;
0112   rc->set_uint64Flag("TIMESTAMP", runnumber);
0113   std::string geofile = CDBInterface::instance()->getUrl("Tracking_Geometry");
0114
0115   Fun4AllRunNodeInputManager *ingeo = new Fun4AllRunNodeInputManager("GeoIn");
0116   ingeo->AddFile(geofile);
0117   se->registerInputManager(ingeo);
0118
0119   // can use for zero field
0120   // double fieldstrength = 0.01;
0121   // G4MAGNET::magfield_tracking = "0.01";
0122   double fieldstrength = std::numeric_limits<double>::quiet_NaN();
0123   bool ConstField = isConstantField(G4MAGNET::magfield_tracking, fieldstrength);
0124   std::cout << "const field " << ConstField << std::endl;
0125   if (ConstField && fieldstrength < 0.1)
0126   {
0127     G4MAGNET::magfield = "0.01";
0128     G4MAGNET::magfield_rescale = 1;
0129   }
0130   TrackingInit();
0131
0132   for (int felix = 0; felix < 6; felix++)
0133   {
0134     Mvtx_HitUnpacking(std::to_string(felix));
0135   }
0136   for (int server = 0; server < 8; server++)
0137   {
0138     Intt_HitUnpacking(std::to_string(server));
0139   }
0140   std::ostringstream ebdcname;
0141   for (int ebdc = 0; ebdc < 24; ebdc++)
0142   {
0143     if (nTpcFiles == 24)
0144     {
0145       ebdcname.str("");
0146       if (ebdc < 10)
0147       {
0148         ebdcname << "0";
0149       }
0150       ebdcname << ebdc;
0151       Tpc_HitUnpacking(ebdcname.str());
0152     }
0153
0154     else if (nTpcFiles == 48)
0155     {
0156       for (int endpoint = 0; endpoint < 2; endpoint++)
0157       {
0158         ebdcname.str("");
0159         if (ebdc < 10)
0160         {
0161           ebdcname << "0";
0162         }
0163         ebdcname << ebdc << "_" << endpoint;
0164         Tpc_HitUnpacking(ebdcname.str());
0165       }
0166     }
0167     else
0168     {
0169       std::cout << "Wrong number of tpc files input " << nTpcFiles << "! Exiting now." << std::endl;
0170       gSystem->Exit(1);
0171     }
0172   }
0173
0174   Micromegas_HitUnpacking();
0175
0176   Mvtx_Clustering();
0177   Intt_Clustering();
0178
0179   Tpc_LaserEventIdentifying();
0180
0181   auto *tpcclusterizer = new TpcClusterizer;
0182   tpcclusterizer->Verbosity(0);
0183   tpcclusterizer->set_reject_event(true);
0184   tpcclusterizer->set_do_hit_association(G4TPC::DO_HIT_ASSOCIATION);
0185   tpcclusterizer->set_rawdata_reco();
0186   se->registerSubsystem(tpcclusterizer);
0187
0188   Micromegas_Clustering();
0189
0190   Tracking_Reco_TrackSeed();
0191
0192   TrackSeedTrackMapConverter *converter = new TrackSeedTrackMapConverter();
0193   // Default set to full SvtxTrackSeeds. Can be set to
0194   // SiliconTrackSeedContainer or TpcTrackSeedContainer
0195   converter->setTrackSeedName("SvtxTrackSeedContainer");
0196   converter->Verbosity(0);
0197   converter->cosmics();
0198   converter->setFieldMap(G4MAGNET::magfield_tracking);
0199   se->registerSubsystem(converter);
0200
0201   std::string residstring = dir + outfilename;
0202
0203   auto *resid = new TrackResiduals("TrackResiduals");
0204   resid->Verbosity(0);
0205   resid->outfileName(residstring);
0206   resid->alignment(false);
0207   resid->clusterTree();
0208   // resid->hitTree();
0209   resid->convertSeeds(true);
0210
0211   if (ConstField && fieldstrength < 0.1)
0212   {
0213     std::cout << "zero field" << std::endl;
0214     resid->zeroField();
0215   }
0216   resid->setSegment(segment);
0217   se->registerSubsystem(resid);
0218
0219   // Fun4AllOutputManager *out = new Fun4AllDstOutputManager("out", "/sphenix/tg/tg01/hf/jdosbo/tracking_development/onlineoffline/hitsets.root");
0220   // se->registerOutputManager(out);
0221
0222   if (Enable::QA)
0223   {
0224     se->registerSubsystem(new MvtxClusterQA);
0225     se->registerSubsystem(new InttClusterQA);
0226     se->registerSubsystem(new TpcClusterQA);
0227     se->registerSubsystem(new MicromegasClusterQA);
0228     se->registerSubsystem(new TpcSeedsQA);
0229   }
0230
0231   se->run(nEvents);
0232   se->End();
0233   se->PrintTimer();
0234
0235   if (Enable::QA)
0236   {
0237     TString qaname = dir + outfilename + "_qa.root";
0238     std::string qaOutputFileName(qaname.Data());
0239     QAHistManagerDef::saveQARootFile(qaOutputFileName);
0240   }
0241   delete se;
0242   std::cout << "Finished" << std::endl;
0243   gSystem->Exit(0);
0244 }