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 #include <fun4all/Fun4AllUtils.h>
0008 #include <G4_ActsGeom.C>
0009 #include <G4_Global.C>
0010 #include <G4_Magnet.C>
0011 #include <GlobalVariables.C>
0012 #include <QA.C>
0013 #include <Trkr_Clustering.C>
0014 #include <Trkr_Reco.C>
0015 #include <Trkr_RecoInit.C>
0016 #include <Trkr_TpcReadoutInit.C>
0017
0018 #include <ffamodules/CDBInterface.h>
0019
0020 #include <fun4all/Fun4AllDstInputManager.h>
0021 #include <fun4all/Fun4AllDstOutputManager.h>
0022 #include <fun4all/Fun4AllInputManager.h>
0023 #include <fun4all/Fun4AllOutputManager.h>
0024 #include <fun4all/Fun4AllRunNodeInputManager.h>
0025 #include <fun4all/Fun4AllServer.h>
0026
0027 #include <phool/recoConsts.h>
0028
0029 #include <cdbobjects/CDBTTree.h>
0030
0031 #include <tpccalib/PHTpcResiduals.h>
0032
0033 #include <trackingqa/SiliconSeedsQA.h>
0034 #include <trackingqa/TpcSeedsQA.h>
0035 #include <trackingqa/TpcSiliconQA.h>
0036
0037 #include <trackingdiagnostics/TrackResiduals.h>
0038 #include <trackingdiagnostics/TrkrNtuplizer.h>
0039
0040 #include <kfparticle_sphenix/KFParticle_sPHENIX.h>
0041
0042 #include <stdio.h>
0043
0044 R__LOAD_LIBRARY(libkfparticle_sphenix.so)
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 void Fun4All_TrackSeeding(
0053     const int nEvents = 5,
0054     const std::string clusterfilename = "DST_TRKR_CLUSTER_run2pp_ana494_2024p021_v001-00053877-00000.root",
0055     const std::string dir = "/sphenix/lustre01/sphnxpro/production/run2pp/physics/ana494_2024p021_v001/DST_TRKR_CLUSTER/run_00053800_00053900/dst/",
0056     const std::string outfilename = "clusters_seeds",
0057     const bool convertSeeds = false,
0058     const bool doKFParticle = false)
0059 {
0060   //std::string inputseedRawHitFile = dir + seedfilename;
0061   std::string inputclusterRawHitFile = dir + clusterfilename;
0062
0063   G4TRACKING::convert_seeds_to_svtxtracks = convertSeeds;
0064   std::cout << "Converting to seeds : " << G4TRACKING::convert_seeds_to_svtxtracks << std::endl;
0065   std::pair<int, int>
0066       runseg = Fun4AllUtils::GetRunSegment(clusterfilename);
0067   int runnumber = runseg.first;
0068   int segment = runseg.second;
0069
0070   auto rc = recoConsts::instance();
0071   rc->set_IntFlag("RUNNUMBER", runnumber);
0072
0073   Enable::CDB = true;
0074   rc->set_StringFlag("CDB_GLOBALTAG", "newcdbtag");
0075   rc->set_uint64Flag("TIMESTAMP", runnumber);
0076   std::string geofile = CDBInterface::instance()->getUrl("Tracking_Geometry");
0077
0078   TpcReadoutInit(runnumber);
0079  // these lines show how to override the drift velocity and time offset values set in TpcReadoutInit
0080   // G4TPC::tpc_drift_velocity_reco = 0.0073844; // cm/ns
0081   // TpcClusterZCrossingCorrection::_vdrift = G4TPC::tpc_drift_velocity_reco;
0082   // G4TPC::tpc_tzero_reco = -5*50;  // ns
0083   std::cout << " run: " << runnumber
0084             << " samples: " << TRACKING::reco_tpc_maxtime_sample
0085             << " pre: " << TRACKING::reco_tpc_time_presample
0086             << " vdrift: " << G4TPC::tpc_drift_velocity_reco
0087             << std::endl;
0088
0089   string outDir = "myKShortReco/";
0090   string outputFileName = "outputKFParticle_" + to_string(runnumber) + "_" + to_string(segment) + ".root";
0091   string outputRecoDir = outDir + "inReconstruction/";
0092   string outputRecoFile = outputRecoDir + outputFileName;
0093
0094   if(doKFParticle){
0095     string makeMainDirectory = "mkdir -p " + outDir;
0096     system(makeMainDirectory.c_str());
0097     string makeDirectory = "mkdir -p " + outputRecoDir;
0098     system(makeDirectory.c_str());
0099   }
0100
0101   // distortion calibration mode
0102   /*
0103    * set to true to enable residuals in the TPC with
0104    * TPC clusters not participating to the ACTS track fit
0105    */
0106   G4TRACKING::SC_CALIBMODE = false;
0107   TRACKING::pp_mode = true;
0108
0109   Enable::MVTX_APPLYMISALIGNMENT = true;
0110   ACTSGEOM::mvtx_applymisalignment = Enable::MVTX_APPLYMISALIGNMENT;
0111
0112   TString outfile = outfilename + "_" + runnumber + "-" + segment + ".root";
0113   std::string theOutfile = outfile.Data();
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::ENABLE_MODULE_EDGE_CORRECTIONS = true;
0123
0124   // to turn on the default static corrections, enable the two lines below
0125   G4TPC::ENABLE_STATIC_CORRECTIONS = true;
0126   G4TPC::USE_PHI_AS_RAD_STATIC_CORRECTIONS = false;
0127
0128   //to turn on the average corrections, enable the three lines below
0129   //note: these are designed to be used only if static corrections are also applied
0130   G4TPC::ENABLE_AVERAGE_CORRECTIONS = true;
0131   G4TPC::USE_PHI_AS_RAD_AVERAGE_CORRECTIONS = false;
0132    // to use a custom file instead of the database file:
0133   G4TPC::average_correction_filename = CDBInterface::instance()->getUrl("TPC_LAMINATION_FIT_CORRECTION");
0134
0135   G4MAGNET::magfield_rescale = 1;
0136   TrackingInit();
0137
0138   auto hitsinclus = new Fun4AllDstInputManager("ClusterInputManager");
0139   hitsinclus->fileopen(inputclusterRawHitFile);
0140   se->registerInputManager(hitsinclus);
0141
0142   // reject laser events if G4TPC::REJECT_LASER_EVENTS is true
0143   Reject_Laser_Events();
0144
0145   /*
0146    * Begin Track Seeding
0147    */
0148
0149   /*
0150    * Silicon Seeding
0151    */
0152
0153   auto silicon_Seeding = new PHActsSiliconSeeding;
0154   silicon_Seeding->Verbosity(0);
0155   silicon_Seeding->setStrobeRange(-5,5);
0156   // these get us to about 83% INTT > 1
0157   silicon_Seeding->setinttRPhiSearchWindow(0.2);
0158   silicon_Seeding->setinttZSearchWindow(1.0);
0159   silicon_Seeding->seedAnalysis(false);
0160   se->registerSubsystem(silicon_Seeding);
0161
0162   auto merger = new PHSiliconSeedMerger;
0163   merger->Verbosity(0);
0164   se->registerSubsystem(merger);
0165
0166   /*
0167    * Tpc Seeding
0168    */
0169   auto seeder = new PHCASeeding("PHCASeeding");
0170   double fieldstrength = std::numeric_limits<double>::quiet_NaN();  // set by isConstantField if constant
0171   bool ConstField = isConstantField(G4MAGNET::magfield_tracking, fieldstrength);
0172   if (ConstField)
0173   {
0174     seeder->useConstBField(true);
0175     seeder->constBField(fieldstrength);
0176   }
0177   else
0178   {
0179     seeder->set_field_dir(-1 * G4MAGNET::magfield_rescale);
0180     seeder->useConstBField(false);
0181     seeder->magFieldFile(G4MAGNET::magfield_tracking);  // to get charge sign right
0182   }
0183   seeder->Verbosity(0);
0184   seeder->SetLayerRange(7, 55);
0185   seeder->SetSearchWindow(2.,0.05); // z-width and phi-width, default in macro at 1.5 and 0.05
0186   seeder->SetClusAdd_delta_window(3.0,0.06); //  (0.5, 0.005) are default; sdzdr_cutoff, d2/dr2(phi)_cutoff
0187   //seeder->SetNClustersPerSeedRange(4,60); // default is 6, 6
0188   seeder->SetMinHitsPerCluster(0);
0189   seeder->SetMinClustersPerTrack(3);
0190   seeder->useFixedClusterError(true);
0191   seeder->set_pp_mode(true);
0192   seeder->reject_zsize1_clusters(true);
0193   se->registerSubsystem(seeder);
0194
0195   // expand stubs in the TPC using simple kalman filter
0196   auto cprop = new PHSimpleKFProp("PHSimpleKFProp");
0197   cprop->set_field_dir(G4MAGNET::magfield_rescale);
0198   if (ConstField)
0199   {
0200     cprop->useConstBField(true);
0201     cprop->setConstBField(fieldstrength);
0202   }
0203   else
0204   {
0205     cprop->magFieldFile(G4MAGNET::magfield_tracking);
0206     cprop->set_field_dir(-1 * G4MAGNET::magfield_rescale);
0207   }
0208   cprop->useFixedClusterError(true);
0209   cprop->set_max_window(5.);
0210   cprop->Verbosity(0);
0211   cprop->set_pp_mode(true);
0212   cprop->set_max_seeds(5000);
0213   se->registerSubsystem(cprop);
0214
0215   // Always apply preliminary distortion corrections to TPC clusters before silicon matching
0216   // and refit the trackseeds. Replace KFProp fits with the new fit parameters in the TPC seeds.
0217   auto prelim_distcorr = new PrelimDistortionCorrection;
0218   prelim_distcorr->set_pp_mode(true);
0219   prelim_distcorr->Verbosity(0);
0220   se->registerSubsystem(prelim_distcorr);
0221
0222   /*
0223    * Track Matching between silicon and TPC
0224    */
0225   // The normal silicon association methods
0226   // Match the TPC track stubs from the CA seeder to silicon track stubs from PHSiliconTruthTrackSeeding
0227   auto silicon_match = new PHSiliconTpcTrackMatching;
0228   silicon_match->Verbosity(0);
0229   silicon_match->set_pp_mode(TRACKING::pp_mode);
0230   if(G4TPC::ENABLE_AVERAGE_CORRECTIONS)
0231   {
0232     // for general tracking
0233     // Eta/Phi window is determined by 3 sigma window
0234     // X/Y/Z window is determined by 4 sigma window
0235     silicon_match->window_deta.set_posQoverpT_maxabs({-0.014,0.0331,0.48});
0236     silicon_match->window_deta.set_negQoverpT_maxabs({-0.006,0.0235,0.52});
0237     silicon_match->set_deltaeta_min(0.03);
0238     silicon_match->window_dphi.set_QoverpT_range({-0.15,0,0}, {0.15,0,0});
0239     silicon_match->window_dx.set_QoverpT_maxabs({3.0,0,0});
0240     silicon_match->window_dy.set_QoverpT_maxabs({3.0,0,0});
0241     silicon_match->window_dz.set_posQoverpT_maxabs({1.138,0.3919,0.84});
0242     silicon_match->window_dz.set_negQoverpT_maxabs({0.719,0.6485,0.65});
0243     silicon_match->set_crossing_deltaz_max(30);
0244     silicon_match->set_crossing_deltaz_min(2);
0245
0246     // for distortion correction using SI-TPOT fit and track pT>0.5
0247     if (G4TRACKING::SC_CALIBMODE)
0248     {
0249       silicon_match->window_deta.set_posQoverpT_maxabs({0.016,0.0060,1.13});
0250       silicon_match->window_deta.set_negQoverpT_maxabs({0.022,0.0022,1.44});
0251       silicon_match->set_deltaeta_min(0.03);
0252       silicon_match->window_dphi.set_QoverpT_range({-0.15,0,0}, {0.09,0,0});
0253       silicon_match->window_dx.set_QoverpT_maxabs({2.0,0,0});
0254       silicon_match->window_dy.set_QoverpT_maxabs({1.5,0,0});
0255       silicon_match->window_dz.set_posQoverpT_maxabs({1.213,0.0211,2.09});
0256       silicon_match->window_dz.set_negQoverpT_maxabs({1.307,0.0001,4.52});
0257       silicon_match->set_crossing_deltaz_min(1.2);
0258     }
0259   }
0260   se->registerSubsystem(silicon_match);
0261
0262   // Match TPC track stubs from CA seeder to clusters in the micromegas layers
0263   auto mm_match = new PHMicromegasTpcTrackMatching;
0264   mm_match->Verbosity(0);
0265   mm_match->set_pp_mode(TRACKING::pp_mode);
0266   mm_match->set_rphi_search_window_lyr1(3.);
0267   mm_match->set_rphi_search_window_lyr2(15.0);
0268   mm_match->set_z_search_window_lyr1(30.0);
0269   mm_match->set_z_search_window_lyr2(3.);
0270
0271   mm_match->set_min_tpc_layer(38);             // layer in TPC to start projection fit
0272   mm_match->set_test_windows_printout(false);  // used for tuning search windows only
0273   se->registerSubsystem(mm_match);
0274
0275   /*
0276    * End Track Seeding
0277    */
0278
0279   /*
0280    * Either converts seeds to tracks with a straight line/helix fit
0281    * or run the full Acts track kalman filter fit
0282    */
0283   if (G4TRACKING::convert_seeds_to_svtxtracks)
0284   {
0285     auto converter = new TrackSeedTrackMapConverter;
0286     // Default set to full SvtxTrackSeeds. Can be set to
0287     // SiliconTrackSeedContainer or TpcTrackSeedContainer
0288     converter->setTrackSeedName("SvtxTrackSeedContainer");
0289     converter->setFieldMap(G4MAGNET::magfield_tracking);
0290     converter->Verbosity(0);
0291     se->registerSubsystem(converter);
0292   }
0293   else
0294   {
0295     auto deltazcorr = new PHTpcDeltaZCorrection;
0296     deltazcorr->Verbosity(0);
0297     se->registerSubsystem(deltazcorr);
0298
0299     // perform final track fit with ACTS
0300     auto actsFit = new PHActsTrkFitter;
0301     actsFit->Verbosity(0);
0302     actsFit->commissioning(G4TRACKING::use_alignment);
0303     // in calibration mode, fit only Silicons and Micromegas hits
0304     actsFit->fitSiliconMMs(G4TRACKING::SC_CALIBMODE);
0305     actsFit->setUseMicromegas(G4TRACKING::SC_USE_MICROMEGAS);
0306     actsFit->set_pp_mode(TRACKING::pp_mode);
0307     actsFit->set_use_clustermover(true);  // default is true for now
0308     actsFit->useActsEvaluator(false);
0309     actsFit->useOutlierFinder(false);
0310     actsFit->setFieldMap(G4MAGNET::magfield_tracking);
0311     se->registerSubsystem(actsFit);
0312
0313     auto cleaner = new PHTrackCleaner();
0314     cleaner->Verbosity(0);
0315     cleaner->set_pp_mode(TRACKING::pp_mode);
0316     se->registerSubsystem(cleaner);
0317
0318     if (G4TRACKING::SC_CALIBMODE)
0319     {
0320       /*
0321        * in calibration mode, calculate residuals between TPC and fitted tracks,
0322        * store in dedicated structure for distortion correction
0323        */
0324       auto residuals = new PHTpcResiduals;
0325       const TString tpc_residoutfile = theOutfile + "_PhTpcResiduals.root";
0326       residuals->setOutputfile(tpc_residoutfile.Data());
0327       residuals->setUseMicromegas(G4TRACKING::SC_USE_MICROMEGAS);
0328
0329       // matches Tony's analysis
0330       residuals->setMinPt(0.2);
0331
0332       // reconstructed distortion grid size (phi, r, z)
0333       residuals->setGridDimensions(36, 48, 80);
0334       se->registerSubsystem(residuals);
0335     }
0336   }
0337
0338   auto finder = new PHSimpleVertexFinder;
0339   finder->Verbosity(0);
0340
0341   //new cuts
0342   finder->setDcaCut(0.05);
0343   finder->setTrackPtCut(0.1);
0344   finder->setBeamLineCut(1);
0345   finder->setTrackQualityCut(300);
0346   finder->setNmvtxRequired(3);
0347   finder->setOutlierPairCut(0.10);
0348
0349   se->registerSubsystem(finder);
0350
0351   // Propagate track positions to the vertex position
0352   auto vtxProp = new PHActsVertexPropagator;
0353   vtxProp->Verbosity(0);
0354   vtxProp->fieldMap(G4MAGNET::magfield_tracking);
0355   se->registerSubsystem(vtxProp);
0356
0357   //run KFParticle
0358   if(doKFParticle){
0359      Global_Reco();
0360
0361   //KFParticle setup
0362
0363   KFParticle_sPHENIX *kfparticle = new KFParticle_sPHENIX("myKShortReco");
0364   kfparticle->Verbosity(1);
0365   kfparticle->setDecayDescriptor("K_S0 -> pi^+ pi^-");
0366
0367   //Basic node selection and configuration
0368   kfparticle->magFieldFile("FIELDMAP_TRACKING");
0369   kfparticle->getAllPVInfo(false);
0370   kfparticle->allowZeroMassTracks(true);
0371   kfparticle->useFakePrimaryVertex(false);
0372   kfparticle->getDetectorInfo(true);
0373
0374   kfparticle->constrainToPrimaryVertex(true);
0375   kfparticle->setMotherIPchi2(FLT_MAX);
0376   kfparticle->setFlightDistancechi2(-1.);
0377   kfparticle->setMinDIRA(-1.1);
0378   kfparticle->setDecayLengthRange(0., FLT_MAX);
0379   kfparticle->setDecayTimeRange(-1*FLT_MAX, FLT_MAX);
0380
0381   //Track parameters
0382   kfparticle->setMinMVTXhits(0);
0383   //kfparticle->setMinINTThits(0);
0384   kfparticle->setMinTPChits(20);
0385   kfparticle->setMinimumTrackPT(-1.);
0386   kfparticle->setMaximumTrackPTchi2(FLT_MAX);
0387   kfparticle->setMinimumTrackIPchi2(-1.);
0388   kfparticle->setMinimumTrackIP(-1.);
0389   //kfparticle->setMaximumTrackchi2nDOF(20.);
0390   kfparticle->setMaximumTrackchi2nDOF(300.);
0391
0392   //Vertex parameters
0393   kfparticle->setMaximumVertexchi2nDOF(50);
0394   kfparticle->setMaximumDaughterDCA(1.);
0395
0396   //Parent parameters
0397   kfparticle->setMotherPT(0);
0398   kfparticle->setMinimumMass(0.200);
0399   kfparticle->setMaximumMass(1.000);
0400   kfparticle->setMaximumMotherVertexVolume(0.1);
0401
0402   kfparticle->setOutputName(outputRecoFile);
0403
0404   se->registerSubsystem(kfparticle);
0405   }
0406
0407   TString residoutfile = theOutfile + "_resid.root";
0408   std::string residstring(residoutfile.Data());
0409
0410   auto resid = new TrackResiduals("TrackResiduals");
0411   resid->outfileName(residstring);
0412   resid->alignment(false);
0413
0414   // adjust track map name
0415   if (G4TRACKING::SC_CALIBMODE && !G4TRACKING::convert_seeds_to_svtxtracks)
0416   {
0417     resid->trackmapName("SvtxSiliconMMTrackMap");
0418     if (G4TRACKING::SC_USE_MICROMEGAS)
0419     {
0420       resid->set_doMicromegasOnly(true);
0421     }
0422   }
0423
0424   resid->clusterTree();
0425   resid->convertSeeds(G4TRACKING::convert_seeds_to_svtxtracks);
0426   resid->Verbosity(0);
0427   se->registerSubsystem(resid);
0428
0429   if (Enable::QA)
0430   {
0431     se->registerSubsystem(new SiliconSeedsQA);
0432     se->registerSubsystem(new TpcSeedsQA);
0433     se->registerSubsystem(new TpcSiliconQA);
0434   }
0435   se->run(nEvents);
0436   se->End();
0437   se->PrintTimer();
0438   CDBInterface::instance()->Print();
0439   if (Enable::QA)
0440   {
0441     TString qaname = theOutfile + "_qa.root";
0442     std::string qaOutputFileName(qaname.Data());
0443     QAHistManagerDef::saveQARootFile(qaOutputFileName);
0444   }
0445
0446   if(doKFParticle){
0447     ifstream file(outputRecoFile.c_str());
0448     if (file.good())
0449     {
0450       string moveOutput = "mv " + outputRecoFile + " " + outDir;
0451       system(moveOutput.c_str());
0452     }
0453   }
0454
0455   delete se;
0456   std::cout << "Finished" << std::endl;
0457   gSystem->Exit(0);
0458 }