rawYieldReco/macro/Fun4All_HF.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.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"};
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(const int nEvents = 500,                                                                                                                                                                                            //
0119                 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", //
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     // se->VerbosityDownscale(500);
0127     auto rc = recoConsts::instance();
0128
0129     if (inputDST.find(".list") != std::string::npos)
0130     {
0131         IsListFile = true;
0132         std::cout << "Input is list file" << std::endl;
0133     }
0134     if (inputDST.find(".root") != std::string::npos)
0135     {
0136         IsDstFile = true;
0137         std::cout << "Input is dst file" << std::endl;
0138     }
0139     if (!IsListFile && !IsDstFile)
0140     {
0141         std::cout << "Check your input! Exit" << std::endl;
0142         gSystem->Exit(0);
0143     }
0144
0145     int runnumber = std::numeric_limits<int>::quiet_NaN();
0146     int segment = std::numeric_limits<int>::quiet_NaN();
0147     std::string runspecies = "";
0148
0149     if (IsListFile)
0150     {
0151         auto *hitsin = new Fun4AllDstInputManager("InputManager");
0152         hitsin->AddListFile(inputDST);
0153         se->registerInputManager(hitsin);
0154
0155         std::ifstream ifs(inputDST);
0156         std::string filepath;
0157         int i = 0;
0158         while (std::getline(ifs, filepath))
0159         {
0160             std::cout << "Adding DST with filepath: " << filepath << std::endl;
0161             if (i == 0)
0162             {
0163                 std::pair<int, int> runseg = Fun4AllUtils::GetRunSegment(filepath);
0164                 runnumber = runseg.first;
0165                 segment = runseg.second;
0166                 runspecies = GetRunInfo(filepath);
0167                 CheckDstType(filepath);
0168                 break;
0169             }
0170             /*
0171             if (!DoSeeding && (get_dEdx_info || get_detector_info))
0172             {
0173               std::pair<int, int> runseg = Fun4AllUtils::GetRunSegment(filepath);
0174               int segmentDST = runseg.second;
0175               std::stringstream nice_segment;
0176               nice_segment << std::setw(5) << std::setfill('0') << to_string(segmentDST);
0177
0178               std::stringstream nice_runnumber;
0179               nice_runnumber << std::setw(8) << std::setfill('0') << to_string(runnumber);
0180
0181               int rounded_up = 100 * (std::ceil((float)(runnumber+0.1) / 100));
0182               std::stringstream nice_rounded_up;
0183               nice_rounded_up << std::setw(8) << std::setfill('0') << to_string(rounded_up);
0184
0185               int rounded_down = 100 * (std::floor((float)(runnumber+0.1) / 100));
0186               std::stringstream nice_rounded_down;
0187               nice_rounded_down << std::setw(8) << std::setfill('0') << to_string(rounded_down);
0188               //unsigned begin = inputDST.find("ana");
0189               //unsigned end = inputDST.find("-");
0190               std::string clus_anacdbver = "ana532_2025p009_v001"; //inputDST.substr(begin, end - begin);
0191               std::string clus_file = "";
0192               if (runspecies == "run3pp" || runspecies == "run3auau")
0193               {
0194                 if (runnumber == 79516)
0195                 {
0196                   clus_file = "/sphenix/lustre01/sphnxpro/production/" + runspecies + "/physics/ana538_2025p011_v001/DST_TRKR_CLUSTER/run_" + nice_rounded_down.str() + "_" + nice_rounded_up.str() + "/DST_TRKR_CLUSTER_" + runspecies + "_" + clus_anacdbver + "-" + nice_runnumber.str() + "-" + nice_segment.str() + ".root";
0197                 }
0198                 else
0199                 {
0200                   clus_file = "/sphenix/lustre01/sphnxpro/production/" + runspecies + "/physics/" + clus_anacdbver + "/DST_TRKR_CLUSTER/run_" + nice_rounded_down.str() + "_" + nice_rounded_up.str() + "/DST_TRKR_CLUSTER_" + runspecies + "_" + clus_anacdbver + "-" + nice_runnumber.str() + "-" + nice_segment.str() + ".root";
0201                 }
0202               }
0203               if (runspecies == "run2pp")
0204               {
0205                 clus_file = "/sphenix/lustre01/sphnxpro/production/" + runspecies + "/physics/" + clus_anacdbver + "/DST_TRKR_CLUSTER/run_" + nice_rounded_down.str() + "_" + nice_rounded_up.str() + "/dst/DST_TRKR_CLUSTER_" + runspecies + "_" + clus_anacdbver + "-" + nice_runnumber.str() + "-" + nice_segment.str() + ".root";
0206               }
0207               std::cout << "Input cluster DST: " << clus_file << std::endl;
0208
0209               std::string inputnameCluster = "InputManagerCluster" + std::to_string(i);
0210               auto hitsinclus = new Fun4AllDstInputManager(inputnameCluster);
0211               hitsinclus->fileopen(clus_file);
0212               se->registerInputManager(hitsinclus);
0213             }
0214             */
0215             i++;
0216         }
0217     }
0218
0219     if (IsDstFile)
0220     {
0221         std::pair<int, int> runseg = Fun4AllUtils::GetRunSegment(inputDST);
0222         runnumber = runseg.first;
0223         segment = runseg.second;
0224         runspecies = GetRunInfo(inputDST);
0225
0226         CheckDstType(inputDST);
0227
0228         auto *hitsin = new Fun4AllDstInputManager("InputManager");
0229         hitsin->fileopen(inputDST);
0230         se->registerInputManager(hitsin);
0231     }
0232
0233     rc->set_IntFlag("RUNNUMBER", runnumber);
0234     rc->set_IntFlag("RUNSEGMENT", segment);
0235     rc->set_uint64Flag("TIMESTAMP", runnumber);
0236
0237     Enable::CDB = true;
0238     rc->set_StringFlag("CDB_GLOBALTAG", "newcdbtag");
0239     rc->set_uint64Flag("TIMESTAMP", runnumber);
0240
0241     std::stringstream nice_runnumber;
0242     nice_runnumber << std::setw(8) << std::setfill('0') << to_string(runnumber);
0243
0244     int rounded_up = 100 * (std::ceil((float)(runnumber+0.1) / 100));
0245     std::stringstream nice_rounded_up;
0246     nice_rounded_up << std::setw(8) << std::setfill('0') << to_string(rounded_up);
0247
0248     int rounded_down = 100 * (std::floor((float)(runnumber+0.1) / 100));
0249     std::stringstream nice_rounded_down;
0250     nice_rounded_down << std::setw(8) << std::setfill('0') << to_string(rounded_down);
0251
0252     std::stringstream nice_segment;
0253     nice_segment << std::setw(5) << std::setfill('0') << to_string(segment);
0254
0255     std::stringstream nice_skip;
0256     nice_skip << std::setw(5) << std::setfill('0') << to_string(nSkip);
0257
0258     if (!DoSeeding && (get_dEdx_info || get_detector_info) && IsDstFile)
0259     {
0260         //unsigned begin = inputDST.find("ana");
0261         //unsigned end = inputDST.find("-");
0262         std::string clus_anacdbver = "ana532_2025p009_v001"; //inputDST.substr(begin, end - begin);
0263         std::string clus_file = "";
0264         if (runspecies == "run3pp" || runspecies == "run3auau")
0265         {
0266             clus_file = "/sphenix/lustre01/sphnxpro/production/" + runspecies + "/physics/" + clus_anacdbver + "/DST_TRKR_CLUSTER/run_" + nice_rounded_down.str() + "_" + nice_rounded_up.str() + "/DST_TRKR_CLUSTER_" + runspecies + "_" + clus_anacdbver + "-" + nice_runnumber.str() + "-" + nice_segment.str() + ".root";
0267         }
0268         if (runspecies == "run2pp")
0269         {
0270             clus_file = "/sphenix/lustre01/sphnxpro/production/" + runspecies + "/physics/" + clus_anacdbver + "/DST_TRKR_CLUSTER/run_" + nice_rounded_down.str() + "_" + nice_rounded_up.str() + "/dst/DST_TRKR_CLUSTER_" + runspecies + "_" + clus_anacdbver + "-" + nice_runnumber.str() + "-" + nice_segment.str() + ".root";
0271         }
0272         std::cout << "Input cluster DST: " << clus_file << std::endl;
0273
0274         auto hitsinclus = new Fun4AllDstInputManager("ClusterInputManager");
0275         hitsinclus->fileopen(clus_file);
0276         se->registerInputManager(hitsinclus);
0277     }
0278
0279     if (get_calo_info)
0280     {
0281         int ratio_nCalo_over_nTrkr = 10; // number of events per calo DST divided by number of events per tracking DST
0282         std::stringstream nice_segment_calo;
0283         nice_segment_calo << std::setw(5) << std::setfill('0') << to_string(segment / ratio_nCalo_over_nTrkr);
0284
0285         std::string calo_anacdbver = "ana509_2024p022_v001";
0286         std::string calo_file = "";
0287         calo_file = "/sphenix/lustre01/sphnxpro/production2/" + runspecies + "/physics/calofitting/" + calo_anacdbver + "/run_" + nice_rounded_down.str() + "_" + nice_rounded_up.str() + "/DST_CALOFITTING_" + runspecies + "_" + calo_anacdbver + "-" + nice_runnumber.str() + "-" + nice_segment_calo.str() + ".root";
0288         std::cout << "Input calo DST: " << calo_file << std::endl;
0289
0290         auto hitsincalo = new Fun4AllDstInputManager("CaloInputManager");
0291         hitsincalo->fileopen(calo_file);
0292         se->registerInputManager(hitsincalo);
0293     }
0294
0295     std::string geofile = CDBInterface::instance()->getUrl("Tracking_Geometry");
0296
0297     Fun4AllRunNodeInputManager *ingeo = new Fun4AllRunNodeInputManager("GeoIn");
0298     ingeo->AddFile(geofile);
0299     se->registerInputManager(ingeo);
0300
0301     // set flags
0302     TRACKING::pp_mode = true;
0303
0304     Enable::MVTX_APPLYMISALIGNMENT = true;
0305     ACTSGEOM::mvtx_applymisalignment = Enable::MVTX_APPLYMISALIGNMENT;
0306
0307     // distortion calibration mode
0308     /*
0309      * set to true to enable residuals in the TPC with
0310      * TPC clusters not participating to the ACTS track fit
0311      */
0312     G4TRACKING::SC_CALIBMODE = false;
0313
0314     TpcReadoutInit(runnumber);
0315     // these lines show how to override the drift velocity and time offset values set in TpcReadoutInit
0316     // G4TPC::tpc_drift_velocity_reco = 0.0073844; // cm/ns
0317     // TpcClusterZCrossingCorrection::_vdrift = G4TPC::tpc_drift_velocity_reco;
0318     // G4TPC::tpc_tzero_reco = -5*50;  // ns
0319     std::cout << " run: " << runnumber << " samples: " << TRACKING::reco_tpc_maxtime_sample << " pre: " << TRACKING::reco_tpc_time_presample << " vdrift: " << G4TPC::tpc_drift_velocity_reco << std::endl;
0320
0321     G4TPC::REJECT_LASER_EVENTS = true;
0322     G4TPC::ENABLE_MODULE_EDGE_CORRECTIONS = true;
0323     // Flag for running the tpc hit unpacker with zero suppression on
0324     TRACKING::tpc_zero_supp = true;
0325
0326     // to turn on the default static corrections, enable the two lines below
0327     G4TPC::ENABLE_STATIC_CORRECTIONS = true;
0328     G4TPC::USE_PHI_AS_RAD_STATIC_CORRECTIONS = false;
0329
0330     // to turn on the average corrections, enable the three lines below
0331     // note: these are designed to be used only if static corrections are also applied
0332     G4TPC::ENABLE_AVERAGE_CORRECTIONS = true;
0333     G4TPC::USE_PHI_AS_RAD_AVERAGE_CORRECTIONS = false;
0334     // to use a custom file instead of the database file:
0335     G4TPC::average_correction_filename = CDBInterface::instance()->getUrl("TPC_LAMINATION_FIT_CORRECTION");
0336
0337     G4MAGNET::magfield_rescale = 1;
0338     TrackingInit();
0339
0340     output_dir = "./"; // Top dir of where the output nTuples will be written
0341     trailer = "_" + nice_runnumber.str() + "_" + nice_segment.str() + "_" + nice_skip.str() + ".root";
0342
0343     if (DoUnpacking)
0344     {
0345         for (int felix = 0; felix < 6; felix++)
0346         {
0347             Mvtx_HitUnpacking(std::to_string(felix));
0348         }
0349         for (int server = 0; server < 8; server++)
0350         {
0351             Intt_HitUnpacking(std::to_string(server));
0352         }
0353         std::ostringstream ebdcname;
0354         if (runspecies == "run3pp" || runspecies == "run3auau")
0355         {
0356             for (int ebdc = 0; ebdc < 24; ebdc++)
0357             {
0358                 for (int endpoint = 0; endpoint < 2; endpoint++)
0359                 {
0360                     ebdcname.str("");
0361                     if (ebdc < 10)
0362                     {
0363                         ebdcname << "0";
0364                     }
0365                     ebdcname << ebdc << "_" << endpoint;
0366                     Tpc_HitUnpacking(ebdcname.str());
0367                 }
0368             }
0369         }
0370         if (runspecies == "run2pp")
0371         {
0372             for (int ebdc = 0; ebdc < 24; ebdc++)
0373             {
0374                 ebdcname.str("");
0375                 if (ebdc < 10)
0376                 {
0377                     ebdcname << "0";
0378                 }
0379                 ebdcname << ebdc;
0380                 Tpc_HitUnpacking(ebdcname.str());
0381             }
0382         }
0383
0384         Micromegas_HitUnpacking();
0385
0386         //se->registerSubsystem(new MvtxRawHitQA);
0387         //se->registerSubsystem(new InttRawHitQA);
0388         //se->registerSubsystem(new TpcRawHitQA);
0389     }
0390
0391     if (DoClustering)
0392     {
0393         Mvtx_Clustering();
0394
0395         Intt_Clustering();
0396
0397         Tpc_LaserEventIdentifying();
0398
0399         TPC_Clustering_run2pp();
0400
0401         Micromegas_Clustering();
0402
0403         Reject_Laser_Events();
0404
0405         //se->registerSubsystem(new MvtxClusterQA);
0406         //se->registerSubsystem(new InttClusterQA);
0407         //se->registerSubsystem(new TpcClusterQA);
0408         //se->registerSubsystem(new MicromegasClusterQA);
0409     }
0410
0411     if (DoSeeding)
0412     {
0413         Tracking_Reco_TrackSeed_run2pp();
0414
0415         auto converter = new TrackSeedTrackMapConverter("SiliconSeedConverter");
0416         // Default set to full SvtxTrackSeeds. Can be set to
0417         // SiliconTrackSeedContainer or TpcTrackSeedContainer
0418         converter->setTrackSeedName("SiliconTrackSeedContainer");
0419         converter->setTrackMapName("SiliconSvtxTrackMap");
0420         converter->setFieldMap(G4MAGNET::magfield_tracking);
0421         converter->Verbosity(0);
0422         se->registerSubsystem(converter);
0423
0424         auto finder = new PHSimpleVertexFinder("SiliconVertexFinder");
0425         finder->Verbosity(0);
0426         finder->setDcaCut(0.1);
0427         finder->setTrackPtCut(0.1);
0428         finder->setBeamLineCut(1);
0429         finder->setTrackQualityCut(1000000000);
0430         finder->setNmvtxRequired(3);
0431         finder->setOutlierPairCut(0.1);
0432         finder->setTrackMapName("SiliconSvtxTrackMap");
0433         finder->setVertexMapName("SiliconSvtxVertexMap");
0434         se->registerSubsystem(finder);
0435
0436         auto siliconqa = new SiliconSeedsQA;
0437         siliconqa->setTrackMapName("SiliconSvtxTrackMap");
0438         siliconqa->setVertexMapName("SiliconSvtxVertexMap");
0439         se->registerSubsystem(siliconqa);
0440
0441         auto convertertpc = new TrackSeedTrackMapConverter("TpcSeedConverter");
0442         // Default set to full SvtxTrackSeeds. Can be set to
0443         // SiliconTrackSeedContainer or TpcTrackSeedContainer
0444         convertertpc->setTrackSeedName("TpcTrackSeedContainer");
0445         convertertpc->setTrackMapName("TpcSvtxTrackMap");
0446         convertertpc->setFieldMap(G4MAGNET::magfield_tracking);
0447         convertertpc->Verbosity(0);
0448         se->registerSubsystem(convertertpc);
0449
0450         auto findertpc = new PHSimpleVertexFinder("TpcSimpleVertexFinder");
0451         findertpc->Verbosity(0);
0452         findertpc->setDcaCut(0.5);
0453         findertpc->setTrackPtCut(0.2);
0454         findertpc->setBeamLineCut(1);
0455         findertpc->setTrackQualityCut(1000000000);
0456         // findertpc->setNmvtxRequired(3);
0457         findertpc->setRequireMVTX(false);
0458         findertpc->setOutlierPairCut(0.1);
0459         findertpc->setTrackMapName("TpcSvtxTrackMap");
0460         findertpc->setVertexMapName("TpcSvtxVertexMap");
0461         se->registerSubsystem(findertpc);
0462
0463         //auto tpcqa = new TpcSeedsQA;
0464         //tpcqa->setTrackMapName("TpcSvtxTrackMap");
0465         //tpcqa->setVertexMapName("TpcSvtxVertexMap");
0466         //tpcqa->setSegment(rc->get_IntFlag("RUNSEGMENT"));
0467         //se->registerSubsystem(tpcqa);
0468     }
0469
0470     if (DoFitting)
0471     {
0472         Tracking_Reco_TrackMatching_run2pp("TRKR_CLUSTER_SEED");
0473
0474         G4TRACKING::convert_seeds_to_svtxtracks = convertSeeds;
0475         std::cout << "Converting to seeds : " << G4TRACKING::convert_seeds_to_svtxtracks << std::endl;
0476         /*
0477          * Either converts seeds to tracks with a straight line/helix fit
0478          * or run the full Acts track kalman filter fit
0479          */
0480         if (G4TRACKING::convert_seeds_to_svtxtracks)
0481         {
0482             auto *converter = new TrackSeedTrackMapConverter;
0483             // Default set to full SvtxTrackSeeds. Can be set to
0484             // SiliconTrackSeedContainer or TpcTrackSeedContainer
0485             converter->setTrackSeedName("SvtxTrackSeedContainer");
0486             converter->setFieldMap(G4MAGNET::magfield_tracking);
0487             converter->Verbosity(0);
0488             se->registerSubsystem(converter);
0489         }
0490         else
0491         {
0492             Tracking_Reco_TrackFit_run2pp("outfile.root", "TRKR_CLUSTER_SEED");
0493             //Tracking_Reco_TrackFit_run2pp();
0494         }
0495
0496         // vertexing and propagation to vertex
0497         //Tracking_Reco_Vertex_run2pp();
0498         Tracking_Reco_Vertex_run2pp("TRKR_CLUSTER_SEED");
0499
0500         //se->registerSubsystem(new TpcSiliconQA);
0501         //se->registerSubsystem(new TrackFittingQA);
0502         //se->registerSubsystem(new VertexQA);
0503     }
0504
0505     if (get_calo_info)
0506     {
0507         std::cout << "Begin my calo reco" << std::endl;
0508
0509         Process_Calo_Calib();
0510
0511         Global_Reco();
0512
0513         bool doEMcalRadiusCorr = true;
0514         auto projection = new PHActsTrackProjection("CaloProjection");
0515         if (doEMcalRadiusCorr)
0516         {
0517             projection->setLayerRadius(SvtxTrack::CEMC, cemc_proj_radius);
0518         }
0519         se->registerSubsystem(projection);
0520
0521         CaloGeomMapping *cgm = new CaloGeomMapping();
0522         cgm->set_detector_name("CEMC");
0523         cgm->set_UseDetailedGeometry(true);
0524         se->registerSubsystem(cgm);
0525
0526         //////////////////
0527         // Clusters
0528         std::cout << "Building clusters" << std::endl;
0529         RawClusterBuilderTemplate *ClusterBuilder = new RawClusterBuilderTemplate("EmcRawClusterBuilderTemplate");
0530         ClusterBuilder->Detector("CEMC");
0531         ClusterBuilder->set_UseDetailedGeometry(true);
0532         ClusterBuilder->set_threshold_energy(0.070); // for when using basic calibration
0533         std::string emc_prof = getenv("CALIBRATIONROOT");
0534         emc_prof += "/EmcProfile/CEMCprof_Thresh30MeV.root";
0535         // ClusterBuilder->set_UseAltZVertex(3); //0: GlobalVertexMap, 1: MbdVertexMap, 2: Nothing, 3: G4TruthInfo
0536         ClusterBuilder->LoadProfile(emc_prof);
0537         ClusterBuilder->set_UseTowerInfo(1); // to use towerinfo objects rather than old RawTower
0538         se->registerSubsystem(ClusterBuilder);
0539     }
0540
0541     output_dir = outDir;
0542
0543     if (run_pipi_reco)
0544         create_hf_directories(pipi_reconstruction_name, pipi_output_dir, pipi_output_reco_file);
0545     if (run_Kpi_reco)
0546         create_hf_directories(Kpi_reconstruction_name, Kpi_output_dir, Kpi_output_reco_file);
0547     if (run_KK_reco)
0548         create_hf_directories(KK_reconstruction_name, KK_output_dir, KK_output_reco_file);
0549     if (run_ppi_reco)
0550         create_hf_directories(ppi_reconstruction_name, ppi_output_dir, ppi_output_reco_file);
0551     if (run_ee_reco)
0552         create_hf_directories(ee_reconstruction_name, ee_output_dir, ee_output_reco_file);
0553     if (run_Lambdapi_reco)
0554         create_hf_directories(Lambdapi_reconstruction_name, Lambdapi_output_dir, Lambdapi_output_reco_file);
0555     if (run_LambdaK_reco)
0556         create_hf_directories(LambdaK_reconstruction_name, LambdaK_output_dir, LambdaK_output_reco_file);
0557
0558     if (run_pipi_reco || run_Kpi_reco || run_KK_reco || run_ppi_reco || run_ee_reco || run_Lambdapi_reco || run_LambdaK_reco)
0559         init_kfp_dependencies();
0560
0561     if (run_pipi_reco)
0562         reconstruct_pipi_mass();
0563     if (run_Kpi_reco)
0564         reconstruct_Kpi_mass();
0565     if (run_KK_reco)
0566         reconstruct_KK_mass();
0567     if (run_ppi_reco)
0568         reconstruct_ppi_mass();
0569     if (run_ee_reco)
0570         reconstruct_ee_mass();
0571     if (run_Lambdapi_reco)
0572         reconstruct_Lambdapi_mass();
0573     if (run_LambdaK_reco)
0574         reconstruct_LambdaK_mass();
0575
0576     se->skip(nSkip);
0577     se->run(nEvents);
0578     se->End();
0579     se->PrintTimer();
0580
0581     /*
0582     TString qaname = output_dir + "qaOut/HIST_" + nice_runnumber.str() + "_" + nice_segment.str() + "_" + nice_skip.str() + "_QA.root";
0583     std::string makeDirectory = "mkdir -p " + output_dir + "qaOut";
0584     if (run_pipi_reco)
0585     {
0586         qaname = pipi_output_dir + "qaOut/HIST_" + nice_runnumber.str() + "_" + nice_segment.str() + "_" + nice_skip.str() + "_QA.root";
0587         makeDirectory = "mkdir -p " + pipi_output_dir + "qaOut";
0588     }
0589     else if (run_Kpi_reco)
0590     {
0591         qaname = Kpi_output_dir + "qaOut/HIST_" + nice_runnumber.str() + "_" + nice_segment.str() + "_" + nice_skip.str() + "_QA.root";
0592         makeDirectory = "mkdir -p " + Kpi_output_dir + "qaOut";
0593     }
0594     else if (run_KK_reco)
0595     {
0596         qaname = KK_output_dir + "qaOut/HIST_" + nice_runnumber.str() + "_" + nice_segment.str() + "_" + nice_skip.str() + "_QA.root";
0597         makeDirectory = "mkdir -p " + KK_output_dir + "qaOut";
0598     }
0599     else if (run_ppi_reco)
0600     {
0601         qaname = ppi_output_dir + "qaOut/HIST_" + nice_runnumber.str() + "_" + nice_segment.str() + "_" + nice_skip.str() + "_QA.root";
0602         makeDirectory = "mkdir -p " + ppi_output_dir + "qaOut";
0603     }
0604     else if (run_ee_reco)
0605     {
0606         qaname = ee_output_dir + "qaOut/HIST_" + nice_runnumber.str() + "_" + nice_segment.str() + "_" + nice_skip.str() + "_QA.root";
0607         makeDirectory = "mkdir -p " + ee_output_dir + "qaOut";
0608     }
0609     else if (run_Lambdapi_reco)
0610     {
0611         qaname = Lambdapi_output_dir + "qaOut/HIST_" + nice_runnumber.str() + "_" + nice_segment.str() + "_" + nice_skip.str() + "_QA.root";
0612         makeDirectory = "mkdir -p " + Lambdapi_output_dir + "qaOut";
0613     }
0614     else if (run_LambdaK_reco)
0615     {
0616         qaname = ee_output_dir + "qaOut/HIST_" + nice_runnumber.str() + "_" + nice_segment.str() + "_" + nice_skip.str() + "_QA.root";
0617         makeDirectory = "mkdir -p " + LambdaK_output_dir + "qaOut";
0618     }
0619     std::cout << "Output QA file: " << qaname << std::endl;
0620     system(makeDirectory.c_str());
0621     std::string qaOutputFileName(qaname.Data());
0622     //QA_Output(qaOutputFileName);
0623     */
0624
0625     if (run_pipi_reco)
0626         end_kfparticle(pipi_output_reco_file, pipi_output_dir);
0627     if (run_Kpi_reco)
0628         end_kfparticle(Kpi_output_reco_file, Kpi_output_dir);
0629     if (run_KK_reco)
0630         end_kfparticle(KK_output_reco_file, KK_output_dir);
0631     if (run_ppi_reco)
0632         end_kfparticle(ppi_output_reco_file, ppi_output_dir);
0633     if (run_ee_reco)
0634         end_kfparticle(ee_output_reco_file, ee_output_dir);
0635     if (run_Lambdapi_reco)
0636         end_kfparticle(Lambdapi_output_reco_file, Lambdapi_output_dir);
0637     if (run_LambdaK_reco)
0638         end_kfparticle(LambdaK_output_reco_file, LambdaK_output_dir);
0639
0640     delete se;
0641
0642     std::cout << "Finished" << std::endl;
0643     gSystem->Exit(0);
0644 }