Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2026-04-07 08:10:32 

0001 #include <fun4all/Fun4AllUtils.h>
0002 #include <G4_ActsGeom.C>
0003 #include <G4_Magnet.C>
0004 #include <GlobalVariables.C>
0005 #include <Trkr_Clustering.C>
0006 #include <Trkr_LaserClustering.C>
0007 #include <Trkr_Reco.C>
0008 #include <Trkr_RecoInit.C>
0009 #include <Trkr_TpcReadoutInit.C>
0010 #include <QA.C>
0011 
0012 #include <ffamodules/CDBInterface.h>
0013 
0014 #include <fun4all/Fun4AllDstInputManager.h>
0015 #include <fun4all/Fun4AllDstOutputManager.h>
0016 #include <fun4all/Fun4AllInputManager.h>
0017 #include <fun4all/Fun4AllOutputManager.h>
0018 #include <fun4all/Fun4AllRunNodeInputManager.h>
0019 #include <fun4all/Fun4AllServer.h>
0020 
0021 #include <phool/recoConsts.h>
0022 
0023 #include <mvtxrawhitqa/MvtxRawHitQA.h>
0024 #include <inttrawhitqa/InttRawHitQA.h>
0025 #include <tpcqa/TpcRawHitQA.h>
0026 #include <trackingqa/InttClusterQA.h>
0027 #include <trackingqa/MicromegasClusterQA.h>
0028 #include <trackingqa/MvtxClusterQA.h>
0029 #include <trackingqa/TpcClusterQA.h>
0030 #include <trackingqa/SiliconSeedsQA.h>
0031 #include <trackingqa/TpcSeedsQA.h>
0032 #include <trackingqa/TrackFittingQA.h>
0033 #include <trackingqa/TpcSiliconQA.h>
0034 #include <trackingqa/VertexQA.h>
0035 
0036 #include <kfparticle_sphenix/KFParticle_sPHENIX.h>
0037 
0038 #include <trackingdiagnostics/TrkrNtuplizer.h>
0039 #include <dedxfitter/dEdxFitter.h>
0040 
0041 #include <stdio.h>
0042 
0043 #include "HF_selections_OO.C"
0044 
0045 using namespace HeavyFlavorReco;
0046 
0047 R__LOAD_LIBRARY(libfun4all.so)
0048 R__LOAD_LIBRARY(libffamodules.so)
0049 R__LOAD_LIBRARY(libphool.so)
0050 R__LOAD_LIBRARY(libcdbobjects.so)
0051 R__LOAD_LIBRARY(libdedxfitter.so)
0052 
0053 bool IsListFile = false;
0054 bool IsDstFile = false;
0055 bool DoUnpacking = false;
0056 bool DoClustering = false;
0057 bool DoSeeding = false;
0058 bool DoFitting = false;
0059 
0060 void CheckDstType(const std::string inputDST)
0061 {
0062   if (inputDST.find("DST_STREAMING_EVENT") != std::string::npos) 
0063   {
0064     DoUnpacking = true;
0065     DoClustering = true;
0066     DoSeeding = true;
0067     DoFitting = true;
0068   }
0069   else if (inputDST.find("DST_TRKR_CLUSTER") != std::string::npos) 
0070   {
0071     DoUnpacking = false;
0072     DoClustering = false;
0073     DoSeeding = true;
0074     DoFitting = true;
0075   }
0076   else if (inputDST.find("DST_TRKR_SEED") != std::string::npos) 
0077   {
0078     DoUnpacking = false;
0079     DoClustering = false;
0080     DoSeeding = false;
0081     DoFitting = true;
0082   }
0083   else if (inputDST.find("DST_TRKR_TRACKS") != std::string::npos)
0084   {
0085     DoUnpacking = false;
0086     DoClustering = false;
0087     DoSeeding = false;
0088     DoFitting = false;
0089   }
0090   return;
0091 }
0092 
0093 void Fun4All_HF_OO(
0094     const int nEvents = 500,//
0095     const std::string inputDST = "/sphenix/lustre01/sphnxpro/production/run2pp/physics/ana506_2024p023_v001/DST_TRKR_TRACKS/run_00053800_00053900/dst/DST_TRKR_TRACKS_run2pp_ana506_2024p023_v001-00053877-00000.root",
0096     const std::string inputDir = "",
0097     const int nSkip = 0,//
0098     const bool convertSeeds = false)//
0099 {
0100   bool runTrkrNtuplizer = true;
0101 
0102   output_dir = "/sphenix/tg/tg01/hf/mjpeters/LightFlavorProduction/data/OO/";//outDir;
0103 
0104   auto se = Fun4AllServer::instance();
0105   se->Verbosity(1);
0106 
0107   if (inputDST.find(".list") != std::string::npos) 
0108   {
0109     IsListFile = true;
0110     std::cout << "Input is list file" << std::endl;
0111   }
0112   if (inputDST.find(".root") != std::string::npos)
0113   {
0114     IsDstFile = true;
0115     std::cout << "Input is dst file" << std::endl;
0116   }
0117   if (!IsListFile && !IsDstFile)
0118   {
0119     std::cout << "Check your input! Exit" << std::endl;
0120     gSystem->Exit(0);
0121   }
0122 
0123   int runnumber = std::numeric_limits<int>::quiet_NaN();
0124   int segment = std::numeric_limits<int>::quiet_NaN();
0125 
0126   if (IsListFile)
0127   {
0128     std::ifstream ifs(inputDST);
0129     std::string filepath;
0130     int i = 0;
0131     while (std::getline(ifs, filepath))
0132     {
0133       std::cout << "Adding DST with filepath: " << filepath << std::endl;
0134       if (i == 0)
0135       {
0136         std::pair<int, int> runseg = Fun4AllUtils::GetRunSegment(filepath);
0137         runnumber = runseg.first;
0138         segment = runseg.second;
0139         CheckDstType(filepath);
0140       }
0141       std::string inputname = "InputManager" + std::to_string(i);
0142       auto *hitsin = new Fun4AllDstInputManager(inputname);
0143       hitsin->fileopen(filepath);
0144       se->registerInputManager(hitsin);
0145       i++;
0146     }
0147   }
0148 
0149   if (IsDstFile)
0150   {
0151     std::pair<int, int> runseg = Fun4AllUtils::GetRunSegment(inputDST);
0152     runnumber = runseg.first;
0153     segment = runseg.second;
0154 
0155     CheckDstType(inputDST);
0156 
0157     std::string inputTrackFile = inputDir + "/" + inputDST;
0158 
0159     auto *hitsin = new Fun4AllDstInputManager("InputManager");
0160     hitsin->fileopen(inputTrackFile);
0161     se->registerInputManager(hitsin);
0162   }
0163 
0164   std::stringstream nice_runnumber;
0165   nice_runnumber << std::setw(8) << std::setfill('0') << to_string(runnumber);
0166 
0167   int rounded_up = 100*(std::ceil((float) runnumber/100));
0168   std::stringstream nice_rounded_up;
0169   nice_rounded_up << std::setw(8) << std::setfill('0') << to_string(rounded_up);
0170 
0171   int rounded_down = 100*(std::floor((float) runnumber/100));
0172   std::stringstream nice_rounded_down;
0173   nice_rounded_down << std::setw(8) << std::setfill('0') << to_string(rounded_down);
0174 
0175   std::stringstream nice_segment;
0176   nice_segment << std::setw(5) << std::setfill('0') << to_string(segment);
0177 
0178   std::stringstream nice_skip;
0179   nice_skip << std::setw(5) << std::setfill('0') << to_string(nSkip);
0180 /*
0181   if (get_trigger_info)
0182   {
0183     std::string gl1_file = "/sphenix/lustre01/sphnxpro/production/run2pp/physics/delme_ana479_nocdbtag_v001/DST_STREAMING_EVENT_INTT0/run_"
0184                          + nice_rounded_down.str() + "_" + nice_rounded_up.str()
0185                          + "/dst/DST_STREAMING_EVENT_INTT0_run2pp_ana479_nocdbtag_v001-" + nice_runnumber.str() + "-" + nice_segment.str() + ".root";
0186 
0187     auto hitsintrig = new Fun4AllDstInputManager("TriggerInputManager");
0188     hitsintrig->fileopen(gl1_file);
0189     se->registerInputManager(hitsintrig);
0190   }
0191 */
0192   if (get_dEdx_info || get_detector_info)
0193   {
0194     std::string clus_file = "/sphenix/lustre01/sphnxpro/production/run3oo/physics/ana534_2025p009_v001/DST_TRKR_CLUSTER/run_"
0195                           + nice_rounded_down.str() + "_" + nice_rounded_up.str()
0196                           + "/DST_TRKR_CLUSTER_run3oo_ana534_2025p009_v001-" + nice_runnumber.str() + "-" + nice_segment.str() + ".root";
0197 
0198     auto hitsinclus = new Fun4AllDstInputManager("ClusterInputManager");
0199     hitsinclus->fileopen(clus_file);
0200     se->registerInputManager(hitsinclus);
0201   }
0202 
0203   auto rc = recoConsts::instance();
0204   rc->set_IntFlag("RUNNUMBER", runnumber);
0205   rc->set_uint64Flag("TIMESTAMP", runnumber);
0206 
0207   Enable::CDB = true;
0208   rc->set_StringFlag("CDB_GLOBALTAG", "newcdbtag");
0209   rc->set_uint64Flag("TIMESTAMP", runnumber);
0210   std::string geofile = CDBInterface::instance()->getUrl("Tracking_Geometry");
0211 
0212   Fun4AllRunNodeInputManager *ingeo = new Fun4AllRunNodeInputManager("GeoIn");
0213   ingeo->AddFile(geofile);
0214   se->registerInputManager(ingeo);
0215 
0216   // set flags
0217   TRACKING::pp_mode = true;
0218 
0219   Enable::MVTX_APPLYMISALIGNMENT = true;
0220   ACTSGEOM::mvtx_applymisalignment = Enable::MVTX_APPLYMISALIGNMENT;
0221 
0222   // distortion calibration mode
0223   /*
0224    * set to true to enable residuals in the TPC with
0225    * TPC clusters not participating to the ACTS track fit
0226    */
0227   G4TRACKING::SC_CALIBMODE = false;
0228 
0229   TpcReadoutInit(runnumber);
0230   // these lines show how to override the drift velocity and time offset values set in TpcReadoutInit
0231   // G4TPC::tpc_drift_velocity_reco = 0.0073844; // cm/ns
0232   // TpcClusterZCrossingCorrection::_vdrift = G4TPC::tpc_drift_velocity_reco;
0233   // G4TPC::tpc_tzero_reco = -5*50;  // ns
0234   std::cout << " run: " << runnumber
0235             << " samples: " << TRACKING::reco_tpc_maxtime_sample
0236             << " pre: " << TRACKING::reco_tpc_time_presample
0237             << " vdrift: " << G4TPC::tpc_drift_velocity_reco
0238             << std::endl;
0239 
0240   G4TPC::REJECT_LASER_EVENTS = true;
0241   G4TPC::ENABLE_MODULE_EDGE_CORRECTIONS = true;
0242   // Flag for running the tpc hit unpacker with zero suppression on
0243   TRACKING::tpc_zero_supp = true;
0244 
0245   // to turn on the default static corrections, enable the two lines below
0246   G4TPC::ENABLE_STATIC_CORRECTIONS = true;
0247   G4TPC::USE_PHI_AS_RAD_STATIC_CORRECTIONS = false;
0248 
0249   // to turn on the average corrections, enable the three lines below
0250   // note: these are designed to be used only if static corrections are also applied
0251   G4TPC::ENABLE_AVERAGE_CORRECTIONS = true;
0252   G4TPC::USE_PHI_AS_RAD_AVERAGE_CORRECTIONS = false;
0253   // to use a custom file instead of the database file:
0254   G4TPC::average_correction_filename = CDBInterface::instance()->getUrl("TPC_LAMINATION_FIT_CORRECTION");
0255 
0256   G4MAGNET::magfield_rescale = 1;
0257   TrackingInit();
0258 
0259   //output_dir = "./"; //Top dir of where the output nTuples will be written
0260   trailer = "_" + nice_runnumber.str() + "_" + nice_segment.str() + "_" + nice_skip.str() + ".root";
0261 
0262   if (DoUnpacking)
0263   {
0264     for (int felix = 0; felix < 6; felix++)
0265     {
0266       Mvtx_HitUnpacking(std::to_string(felix));
0267     }
0268     for (int server = 0; server < 8; server++)
0269     {
0270       Intt_HitUnpacking(std::to_string(server));
0271     }
0272     std::ostringstream ebdcname;
0273     for (int ebdc = 0; ebdc < 24; ebdc++)
0274     {
0275       ebdcname.str("");
0276       if (ebdc < 10)
0277       {
0278         ebdcname << "0";
0279       }
0280       ebdcname << ebdc;
0281       Tpc_HitUnpacking(ebdcname.str());
0282     }
0283 
0284     Micromegas_HitUnpacking();
0285 
0286     se->registerSubsystem(new MvtxRawHitQA);
0287     se->registerSubsystem(new InttRawHitQA);
0288     se->registerSubsystem(new TpcRawHitQA);
0289   }
0290 
0291   if (DoClustering)
0292   {
0293     Mvtx_Clustering();
0294 
0295     Intt_Clustering();
0296 
0297     Tpc_LaserEventIdentifying();
0298 
0299     TPC_Clustering_run2pp();
0300 
0301     Micromegas_Clustering();
0302 
0303     Reject_Laser_Events();
0304 /*
0305     se->registerSubsystem(new MvtxClusterQA);
0306     se->registerSubsystem(new InttClusterQA);
0307     se->registerSubsystem(new TpcClusterQA);
0308     se->registerSubsystem(new MicromegasClusterQA);
0309 */
0310   }
0311 
0312   if (DoSeeding)
0313   {
0314     Tracking_Reco_TrackSeed_run2pp();
0315 
0316     auto converter = new TrackSeedTrackMapConverter("SiliconSeedConverter");
0317     // Default set to full SvtxTrackSeeds. Can be set to
0318     // SiliconTrackSeedContainer or TpcTrackSeedContainer
0319     converter->setTrackSeedName("SiliconTrackSeedContainer");
0320     converter->setTrackMapName("SiliconSvtxTrackMap");
0321     converter->setFieldMap(G4MAGNET::magfield_tracking);
0322     converter->Verbosity(0);
0323     se->registerSubsystem(converter);
0324 
0325     auto finder = new PHSimpleVertexFinder("SiliconVertexFinder");
0326     finder->Verbosity(0);
0327     finder->setDcaCut(0.1);
0328     finder->setTrackPtCut(0.1);
0329     finder->setBeamLineCut(1);
0330     finder->setTrackQualityCut(1000000000);
0331     finder->setNmvtxRequired(3);
0332     finder->setOutlierPairCut(0.1);
0333     finder->setTrackMapName("SiliconSvtxTrackMap");
0334     finder->setVertexMapName("SiliconSvtxVertexMap");
0335     se->registerSubsystem(finder);
0336 
0337     auto siliconqa = new SiliconSeedsQA;
0338     siliconqa->setTrackMapName("SiliconSvtxTrackMap");
0339     siliconqa->setVertexMapName("SiliconSvtxVertexMap");
0340     se->registerSubsystem(siliconqa);
0341 
0342     auto convertertpc = new TrackSeedTrackMapConverter("TpcSeedConverter");
0343     // Default set to full SvtxTrackSeeds. Can be set to
0344     // SiliconTrackSeedContainer or TpcTrackSeedContainer
0345     convertertpc->setTrackSeedName("TpcTrackSeedContainer");
0346     convertertpc->setTrackMapName("TpcSvtxTrackMap");
0347     convertertpc->setFieldMap(G4MAGNET::magfield_tracking);
0348     convertertpc->Verbosity(0);
0349     se->registerSubsystem(convertertpc);
0350 
0351     auto findertpc = new PHSimpleVertexFinder("TpcSimpleVertexFinder");
0352     findertpc->Verbosity(0);
0353     findertpc->setDcaCut(0.5);
0354     findertpc->setTrackPtCut(0.2);
0355     findertpc->setBeamLineCut(1);
0356     findertpc->setTrackQualityCut(1000000000);
0357     //findertpc->setNmvtxRequired(3);
0358     findertpc->setRequireMVTX(false);
0359     findertpc->setOutlierPairCut(0.1);
0360     findertpc->setTrackMapName("TpcSvtxTrackMap");
0361     findertpc->setVertexMapName("TpcSvtxVertexMap");
0362     se->registerSubsystem(findertpc);
0363 
0364     auto tpcqa = new TpcSeedsQA;
0365     tpcqa->setTrackMapName("TpcSvtxTrackMap");
0366     tpcqa->setVertexMapName("TpcSvtxVertexMap");
0367     tpcqa->setSegment(rc->get_IntFlag("RUNSEGMENT"));
0368     se->registerSubsystem(tpcqa);
0369   }
0370 
0371   if (DoFitting)
0372   {
0373     Tracking_Reco_TrackMatching_run2pp();
0374 
0375     G4TRACKING::convert_seeds_to_svtxtracks = convertSeeds;
0376     std::cout << "Converting to seeds : " << G4TRACKING::convert_seeds_to_svtxtracks << std::endl;
0377     /*
0378      * Either converts seeds to tracks with a straight line/helix fit
0379      * or run the full Acts track kalman filter fit
0380      */
0381     if (G4TRACKING::convert_seeds_to_svtxtracks)
0382     {
0383       auto *converter = new TrackSeedTrackMapConverter;
0384       // Default set to full SvtxTrackSeeds. Can be set to
0385       // SiliconTrackSeedContainer or TpcTrackSeedContainer
0386       converter->setTrackSeedName("SvtxTrackSeedContainer");
0387       converter->setFieldMap(G4MAGNET::magfield_tracking);
0388       converter->Verbosity(0);
0389       se->registerSubsystem(converter);
0390     }
0391     else
0392     {
0393       Tracking_Reco_TrackFit_run2pp();
0394     }
0395 
0396     //vertexing and propagation to vertex
0397     Tracking_Reco_Vertex_run2pp();
0398 /*
0399     se->registerSubsystem(new TpcSiliconQA);
0400     se->registerSubsystem(new TrackFittingQA);
0401     se->registerSubsystem(new VertexQA);
0402 */
0403   }
0404 
0405   string nTuplizer_dir = output_dir + "/nTP/";
0406   string NTP_output_reco_file = nTuplizer_dir + processing_folder + "nTuplizer_tracks" + trailer;
0407 
0408   if (runTrkrNtuplizer)
0409   {
0410     std::string makeDirectory = "mkdir -p " + nTuplizer_dir + processing_folder;
0411     system(makeDirectory.c_str());
0412 
0413     TrkrNtuplizer* eval = new TrkrNtuplizer("TrkrNtuplizer", NTP_output_reco_file, "SvtxTrackMap", 3, 4, 48);
0414     int do_default = 0;
0415     eval->Verbosity(0);
0416     eval->runnumber(runnumber);
0417     eval->segment(segment);
0418     eval->do_hit_eval(false);
0419     eval->do_cluster_eval(false);
0420     eval->do_clus_trk_eval(false);
0421     eval->do_vertex_eval(true);
0422     eval->do_track_eval(true);
0423     eval->do_tpcseed_eval(false);
0424     eval->do_siseed_eval(false);
0425     se->registerSubsystem(eval);
0426   }
0427 
0428   dEdxFitter* fitter = new dEdxFitter("dEdxFitter");
0429   fitter->set_filename((output_dir + "/dedx/dedx_" + nice_runnumber.str() + "_" + nice_segment.str() + trailer).c_str());
0430   fitter->Verbosity(1);
0431   se->registerSubsystem(fitter);
0432 
0433   if (run_pipi_reco) create_hf_directories(pipi_reconstruction_name, pipi_output_dir, pipi_output_reco_file);
0434   if (run_Kpi_reco) create_hf_directories(Kpi_reconstruction_name, Kpi_output_dir, Kpi_output_reco_file);
0435   if (run_KK_reco) create_hf_directories(KK_reconstruction_name, KK_output_dir, KK_output_reco_file);
0436   if (run_ppi_reco) create_hf_directories(ppi_reconstruction_name, ppi_output_dir, ppi_output_reco_file);
0437 
0438   if (run_pipi_reco || run_Kpi_reco || run_KK_reco || run_ppi_reco) init_kfp_dependencies();
0439 
0440   if (run_pipi_reco) reconstruct_pipi_mass();
0441   if (run_Kpi_reco) reconstruct_Kpi_mass();
0442   if (run_KK_reco) reconstruct_KK_mass();
0443   if (run_ppi_reco) reconstruct_ppi_mass();
0444 
0445   
0446 
0447 /*
0448   string decayTag = "200 GeV O+O";
0449 
0450   event_display_maker *myEvtDisplays = new event_display_maker();
0451   myEvtDisplays->setDecayTag(decayTag.c_str());
0452   myEvtDisplays->setKFParticleContainerName("pipi_reco");
0453   myEvtDisplays->setEventDisplayPath("/sphenix/user/cdean/software/TrackingAnalysis/Physics_Val_TF/data/macro/displays/");
0454   myEvtDisplays->plotsApproved(false);
0455   myEvtDisplays->addCaloInfo();
0456   myEvtDisplays->setMaxEvtDisplays(5);
0457   se->registerSubsystem(myEvtDisplays);
0458 */
0459   se->skip(nSkip);
0460   se->run(nEvents);
0461   se->End();
0462   se->PrintTimer();
0463 
0464   TString qaname = output_dir + "qaOut/HIST_" + nice_runnumber.str() + "_" + nice_segment.str() + "_" + nice_skip.str() + "_QA.root";
0465   std::string makeDirectory = "mkdir -p " + output_dir + "qaOut";
0466   if (run_pipi_reco) 
0467   {
0468     qaname = pipi_output_dir + "qaOut/HIST_" + nice_runnumber.str() + "_" + nice_segment.str() + "_" + nice_skip.str() + "_QA.root";
0469     makeDirectory = "mkdir -p " + pipi_output_dir + "qaOut";
0470   }
0471   else if (run_Kpi_reco)
0472   {
0473     qaname = Kpi_output_dir + "qaOut/HIST_" + nice_runnumber.str() + "_" + nice_segment.str() + "_" + nice_skip.str() + "_QA.root";
0474     makeDirectory = "mkdir -p " + Kpi_output_dir + "qaOut";
0475   }
0476   else if (run_KK_reco)
0477   {
0478     qaname = KK_output_dir + "qaOut/HIST_" + nice_runnumber.str() + "_" + nice_segment.str() + "_" + nice_skip.str() + "_QA.root";
0479     makeDirectory = "mkdir -p " + KK_output_dir + "qaOut";
0480   }
0481   else if (run_ppi_reco)
0482   {
0483     qaname = ppi_output_dir + "qaOut/HIST_" + nice_runnumber.str() + "_" + nice_segment.str() + "_" + nice_skip.str() + "_QA.root";
0484     makeDirectory = "mkdir -p " + ppi_output_dir + "qaOut";
0485   }
0486   std::cout << "Output QA file: " << qaname << std::endl;
0487   system(makeDirectory.c_str());
0488   std::string qaOutputFileName(qaname.Data());
0489   //QA_Output(qaOutputFileName);
0490 
0491   if (run_pipi_reco) end_kfparticle(pipi_output_reco_file, pipi_output_dir);
0492   if (run_Kpi_reco) end_kfparticle(Kpi_output_reco_file, Kpi_output_dir);
0493   if (run_KK_reco) end_kfparticle(KK_output_reco_file, KK_output_dir);
0494   if (run_ppi_reco) end_kfparticle(ppi_output_reco_file, ppi_output_dir);
0495   if (runTrkrNtuplizer) end_kfparticle(NTP_output_reco_file, nTuplizer_dir);
0496 
0497   delete se;
0498 
0499   std::cout << "Finished" << std::endl;
0500   gSystem->Exit(0);
0501 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!