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