macros/TrackingProduction/Fun4All_TrackSeeding.C

0001 /*
0002  * This macro shows a working example of running TrackSeeding over the cluster DST
0003  * This has track residuals as default output but has KFParticle set up with a togglable flag
0004  * with the default set up for K Short reconstruction
0005  */
0006
0007 // GlobalVariables.C has to be first, an empty line afterwards
0008 // protects its position against reshuffling by clang-format
0009 #include <GlobalVariables.C>
0010
0011 #include <G4_ActsGeom.C>
0012 #include <G4_Global.C>
0013 #include <G4_Magnet.C>
0014 #include <QA.C>
0015 #include <Trkr_Clustering.C>
0016 #include <Trkr_Reco.C>
0017 #include <Trkr_RecoInit.C>
0018 #include <Trkr_TpcReadoutInit.C>
0019
0020
0021 #include <tpccalib/PHTpcResiduals.h>
0022
0023 #include <trackingqa/SiliconSeedsQA.h>
0024 #include <trackingqa/TpcSeedsQA.h>
0025 #include <trackingqa/TpcSiliconQA.h>
0026
0027 #include <trackingdiagnostics/TrackResiduals.h>
0028 #include <trackingdiagnostics/TrkrNtuplizer.h>
0029
0030 #include <kfparticle_sphenix/KFParticle_sPHENIX.h>
0031
0032 #include <ffamodules/CDBInterface.h>
0033
0034 #include <fun4all/Fun4AllDstInputManager.h>
0035 #include <fun4all/Fun4AllDstOutputManager.h>
0036 #include <fun4all/Fun4AllInputManager.h>
0037 #include <fun4all/Fun4AllOutputManager.h>
0038 #include <fun4all/Fun4AllRunNodeInputManager.h>
0039 #include <fun4all/Fun4AllServer.h>
0040 #include <fun4all/Fun4AllUtils.h>
0041
0042 #include <phool/recoConsts.h>
0043
0044
0045
0046 R__LOAD_LIBRARY(libfun4all.so)
0047 R__LOAD_LIBRARY(libffamodules.so)
0048 R__LOAD_LIBRARY(libphool.so)
0049 R__LOAD_LIBRARY(libcdbobjects.so)
0050 R__LOAD_LIBRARY(libTrackingDiagnostics.so)
0051 R__LOAD_LIBRARY(libtrackingqa.so)
0052 R__LOAD_LIBRARY(libkfparticle_sphenix.so)
0053
0054 void Fun4All_TrackSeeding(
0055     const int nEvents = 5,
0056     const std::string &clusterfilename = "DST_TRKR_CLUSTER_run2pp_ana517_2024p024_v001-00053877-02611.root",
0057     const std::string &dir = "/sphenix/lustre01/sphnxpro/production/run2pp/physics/ana517_2024p024_v001/DST_TRKR_CLUSTER/run_00053800_00053900/dst/",
0058     const std::string &outfilename = "clusters_seeds",
0059     const bool convertSeeds = false,
0060     const bool doKFParticle = false)
0061 {
0062   //std::string inputseedRawHitFile = dir + seedfilename;
0063   std::string inputclusterRawHitFile = dir + clusterfilename;
0064
0065   G4TRACKING::convert_seeds_to_svtxtracks = convertSeeds;
0066   std::cout << "Converting to seeds : " << G4TRACKING::convert_seeds_to_svtxtracks << std::endl;
0067   std::pair<int, int> runseg = Fun4AllUtils::GetRunSegment(clusterfilename);
0068   int runnumber = runseg.first;
0069   int segment = runseg.second;
0070
0071   auto *rc = recoConsts::instance();
0072   rc->set_IntFlag("RUNNUMBER", runnumber);
0073
0074   Enable::CDB = true;
0075   rc->set_StringFlag("CDB_GLOBALTAG", "newcdbtag");
0076   rc->set_uint64Flag("TIMESTAMP", runnumber);
0077   std::string geofile = CDBInterface::instance()->getUrl("Tracking_Geometry");
0078
0079   TpcReadoutInit(runnumber);
0080  // these lines show how to override the drift velocity and time offset values set in TpcReadoutInit
0081   // G4TPC::tpc_drift_velocity_reco = 0.0073844; // cm/ns
0082   // TpcClusterZCrossingCorrection::_vdrift = G4TPC::tpc_drift_velocity_reco;
0083   // G4TPC::tpc_tzero_reco = -5*50;  // ns
0084   std::cout << " run: " << runnumber
0085             << " samples: " << TRACKING::reco_tpc_maxtime_sample
0086             << " pre: " << TRACKING::reco_tpc_time_presample
0087             << " vdrift: " << G4TPC::tpc_drift_velocity_reco
0088             << std::endl;
0089
0090   std::string outDir = "myKShortReco/";
0091   std::string outputFileName = "outputKFParticle_" + std::to_string(runnumber) + "_" + std::to_string(segment) + ".root";
0092   std::string outputRecoDir = outDir + "inReconstruction/";
0093   std::string outputRecoFile = outputRecoDir + outputFileName;
0094
0095   if(doKFParticle){
0096     std::string makeMainDirectory = "mkdir -p " + outDir;
0097     system(makeMainDirectory.c_str());
0098     std::string makeDirectory = "mkdir -p " + outputRecoDir;
0099     system(makeDirectory.c_str());
0100   }
0101
0102   // distortion calibration mode
0103   /*
0104    * set to true to enable residuals in the TPC with
0105    * TPC clusters not participating to the ACTS track fit
0106    */
0107   G4TRACKING::SC_CALIBMODE = false;
0108   TRACKING::pp_mode = true;
0109
0110   Enable::MVTX_APPLYMISALIGNMENT = true;
0111   ACTSGEOM::mvtx_applymisalignment = Enable::MVTX_APPLYMISALIGNMENT;
0112
0113   std::string theOutfile = outfilename + "_" + std::to_string(runnumber) + "-" + std::to_string(segment) + ".root";
0114   auto *se = Fun4AllServer::instance();
0115   se->Verbosity(1);
0116
0117
0118   Fun4AllRunNodeInputManager *ingeo = new Fun4AllRunNodeInputManager("GeoIn");
0119   ingeo->AddFile(geofile);
0120   se->registerInputManager(ingeo);
0121
0122   G4TPC::REJECT_LASER_EVENTS=true;
0123   G4TPC::ENABLE_MODULE_EDGE_CORRECTIONS = true;
0124
0125   // to turn on the default static corrections, enable the two lines below
0126   G4TPC::ENABLE_STATIC_CORRECTIONS = true;
0127   G4TPC::USE_PHI_AS_RAD_STATIC_CORRECTIONS = false;
0128
0129   //to turn on the average corrections, enable the three lines below
0130   //note: these are designed to be used only if static corrections are also applied
0131   G4TPC::ENABLE_AVERAGE_CORRECTIONS = true;
0132   G4TPC::USE_PHI_AS_RAD_AVERAGE_CORRECTIONS = false;
0133    // to use a custom file instead of the database file:
0134   G4TPC::average_correction_filename = CDBInterface::instance()->getUrl("TPC_LAMINATION_FIT_CORRECTION");
0135
0136   G4MAGNET::magfield_rescale = 1;
0137   TrackingInit();
0138
0139   auto *hitsinclus = new Fun4AllDstInputManager("ClusterInputManager");
0140   hitsinclus->fileopen(inputclusterRawHitFile);
0141   se->registerInputManager(hitsinclus);
0142
0143   // reject laser events if G4TPC::REJECT_LASER_EVENTS is true
0144   Reject_Laser_Events();
0145
0146   /*
0147    * Begin Track Seeding
0148    */
0149
0150   Tracking_Reco_TrackSeed_run2pp();
0151   Tracking_Reco_TrackMatching_run2pp();
0152
0153   /*
0154    * Either converts seeds to tracks with a straight line/helix fit
0155    * or run the full Acts track kalman filter fit
0156    */
0157   if (G4TRACKING::convert_seeds_to_svtxtracks)
0158   {
0159     auto *converter = new TrackSeedTrackMapConverter;
0160     // Default set to full SvtxTrackSeeds. Can be set to
0161     // SiliconTrackSeedContainer or TpcTrackSeedContainer
0162     converter->setTrackSeedName("SvtxTrackSeedContainer");
0163     converter->setFieldMap(G4MAGNET::magfield_tracking);
0164     converter->Verbosity(0);
0165     se->registerSubsystem(converter);
0166   }
0167   else
0168   {
0169     Tracking_Reco_TrackFit_run2pp(theOutfile);
0170   }
0171
0172   //vertexing and propagation to vertex
0173   Tracking_Reco_Vertex_run2pp();
0174
0175   //run KFParticle
0176   if(doKFParticle){
0177      Global_Reco();
0178
0179   //KFParticle setup
0180
0181   KFParticle_sPHENIX *kfparticle = new KFParticle_sPHENIX("myKShortReco");
0182   kfparticle->Verbosity(1);
0183   kfparticle->setDecayDescriptor("K_S0 -> pi^+ pi^-");
0184
0185   //Basic node selection and configuration
0186   kfparticle->magFieldFile("FIELDMAP_TRACKING");
0187   kfparticle->getAllPVInfo(false);
0188   kfparticle->allowZeroMassTracks(true);
0189   kfparticle->useFakePrimaryVertex(false);
0190   kfparticle->getDetectorInfo(true);
0191
0192   kfparticle->constrainToPrimaryVertex(true);
0193   kfparticle->setMotherIPchi2(FLT_MAX);
0194   kfparticle->setFlightDistancechi2(-1.);
0195   kfparticle->setMinDIRA(-1.1);
0196   kfparticle->setDecayLengthRange(0., FLT_MAX);
0197   kfparticle->setDecayTimeRange(-1*FLT_MAX, FLT_MAX);
0198
0199   //Track parameters
0200   kfparticle->setMinMVTXhits(0);
0201   //kfparticle->setMinINTThits(0);
0202   kfparticle->setMinTPChits(20);
0203   kfparticle->setMinimumTrackPT(-1.);
0204   kfparticle->setMaximumTrackPTchi2(FLT_MAX);
0205   kfparticle->setMinimumTrackIPchi2(-1.);
0206   kfparticle->setMinimumTrackIP(-1.);
0207   //kfparticle->setMaximumTrackchi2nDOF(20.);
0208   kfparticle->setMaximumTrackchi2nDOF(300.);
0209
0210   //Vertex parameters
0211   kfparticle->setMaximumVertexchi2nDOF(50);
0212   kfparticle->setMaximumDaughterDCA(1.);
0213
0214   //Parent parameters
0215   kfparticle->setMotherPT(0);
0216   kfparticle->setMinimumMass(0.200);
0217   kfparticle->setMaximumMass(1.000);
0218   kfparticle->setMaximumMotherVertexVolume(0.1);
0219
0220   kfparticle->setOutputName(outputRecoFile);
0221
0222   se->registerSubsystem(kfparticle);
0223   }
0224
0225   std::string residstring = theOutfile + "_resid.root";
0226
0227   auto *resid = new TrackResiduals("TrackResiduals");
0228   resid->outfileName(residstring);
0229   resid->clusterTree();
0230   resid->alignment(false);
0231
0232   // adjust track map name
0233   if (G4TRACKING::SC_CALIBMODE && !G4TRACKING::convert_seeds_to_svtxtracks)
0234   {
0235     resid->trackmapName("SvtxSiliconMMTrackMap");
0236     if (G4TRACKING::SC_USE_MICROMEGAS)
0237     {
0238       resid->set_doMicromegasOnly(true);
0239     }
0240   }
0241
0242   resid->clusterTree();
0243   resid->convertSeeds(G4TRACKING::convert_seeds_to_svtxtracks);
0244   resid->Verbosity(0);
0245   se->registerSubsystem(resid);
0246
0247   if (Enable::QA)
0248   {
0249     se->registerSubsystem(new SiliconSeedsQA);
0250     se->registerSubsystem(new TpcSeedsQA);
0251     se->registerSubsystem(new TpcSiliconQA);
0252   }
0253   se->run(nEvents);
0254   se->End();
0255   se->PrintTimer();
0256   CDBInterface::instance()->Print();
0257   if (Enable::QA)
0258   {
0259     std::string qaOutputFileName = theOutfile + "_qa.root";
0260     QAHistManagerDef::saveQARootFile(qaOutputFileName);
0261   }
0262
0263   if(doKFParticle){
0264     std::ifstream file(outputRecoFile.c_str());
0265     if (file.good())
0266     {
0267       std::string moveOutput = "mv " + outputRecoFile + " " + outDir;
0268       system(moveOutput.c_str());
0269     }
0270   }
0271
0272   delete se;
0273   std::cout << "Finished" << std::endl;
0274   gSystem->Exit(0);
0275 }