macros/data/Fun4All_HF_SQM26.C

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