TrackingProduction/run3auau/Fun4All_raw_hit_KFP.C

0001 /*
0002  * This macro shows a minimum working example of running the tracking
0003  * hit unpackers with some basic seeding algorithms to try to put together
0004  * tracks. There are some analysis modules run at the end which package
0005  * hits, clusters, and clusters on tracks into trees for analysis.
0006  */
0007
0008 // leave the GlobalVariables.C at the beginning, an empty line afterwards
0009 // protects its position against reshuffling by clang-format
0010 #include <GlobalVariables.C>
0011
0012 #include <G4_ActsGeom.C>
0013 #include <G4_Global.C>
0014 #include <G4_Magnet.C>
0015 #include <G4_Mbd.C>
0016 #include <QA.C>
0017 #include <Trkr_Clustering.C>
0018 #include <Trkr_LaserClustering.C>
0019 #include <Trkr_Reco.C>
0020 #include <Trkr_RecoInit.C>
0021 #include <Trkr_TpcReadoutInit.C>
0022
0023 #include <cdbobjects/CDBTTree.h>
0024
0025 #include <tpccalib/PHTpcResiduals.h>
0026
0027 #include <mvtxrawhitqa/MvtxRawHitQA.h>
0028
0029 #include <inttrawhitqa/InttRawHitQA.h>
0030
0031 #include <trackingqa/InttClusterQA.h>
0032 #include <trackingqa/MicromegasClusterQA.h>
0033 #include <trackingqa/MvtxClusterQA.h>
0034 #include <trackingqa/SiliconSeedsQA.h>
0035 #include <trackingqa/TpcClusterQA.h>
0036 #include <trackingqa/TpcSeedsQA.h>
0037 #include <trackingqa/TpcSiliconQA.h>
0038
0039 #include <tpcqa/TpcRawHitQA.h>
0040
0041 #include <trackingdiagnostics/KshortReconstruction.h>
0042 #include <trackingdiagnostics/TrackResiduals.h>
0043 #include <trackingdiagnostics/TrkrNtuplizer.h>
0044
0045 #include <kfparticle_sphenix/KFParticle_sPHENIX.h>
0046
0047 #include <ffamodules/CDBInterface.h>
0048 #include <ffamodules/FlagHandler.h>
0049
0050 #include <fun4all/Fun4AllDstInputManager.h>
0051 #include <fun4all/Fun4AllDstOutputManager.h>
0052 #include <fun4all/Fun4AllInputManager.h>
0053 #include <fun4all/Fun4AllOutputManager.h>
0054 #include <fun4all/Fun4AllRunNodeInputManager.h>
0055 #include <fun4all/Fun4AllServer.h>
0056 #include <fun4all/Fun4AllUtils.h>
0057
0058 #include <phool/recoConsts.h>
0059
0060 R__LOAD_LIBRARY(libkfparticle_sphenix.so)
0061 R__LOAD_LIBRARY(libfun4all.so)
0062 R__LOAD_LIBRARY(libffamodules.so)
0063 R__LOAD_LIBRARY(libphool.so)
0064 R__LOAD_LIBRARY(libcdbobjects.so)
0065 R__LOAD_LIBRARY(libmvtx.so)
0066 R__LOAD_LIBRARY(libintt.so)
0067 R__LOAD_LIBRARY(libtpc.so)
0068 R__LOAD_LIBRARY(libmicromegas.so)
0069 R__LOAD_LIBRARY(libTrackingDiagnostics.so)
0070 R__LOAD_LIBRARY(libtrackingqa.so)
0071 R__LOAD_LIBRARY(libtpcqa.so)
0072
0073 namespace HeavyFlavorReco
0074 {
0075   int VERBOSITY = 0;
0076
0077   std::string output_dir = "./";  // Top dir of where the output nTuples will be written
0078   std::string kfp_header = "outputKFParticle_";
0079   std::string processing_folder = "inReconstruction/";
0080   std::string trailer = ".root";
0081
0082   std::string pipi_decay_descriptor = "K_S0 -> pi^+ pi^-";  // See twiki on how to set this
0083   std::string pipi_reconstruction_name = "pipi_reco";       // Used for naming output folder, file and node
0084   std::string pipi_output_reco_file;
0085   std::string pipi_output_dir;
0086
0087   std::string ppi_decay_descriptor = "[Lambda0 -> proton^+ pi^-]cc";  // See twiki on how to set this
0088   std::string ppi_reconstruction_name = "ppi_reco";                   // Used for naming output folder, file and node
0089   std::string ppi_output_reco_file;
0090   std::string ppi_output_dir;
0091
0092   bool save_tracks_to_DST = false;
0093   bool dont_use_global_vertex = true;
0094   bool require_track_and_vertex_match = true;
0095   bool save_all_vtx_info = true;
0096   bool constrain_phi_mass = false;
0097   bool constrain_lambda_mass = false;
0098   bool constrain_D_mass = false;
0099   bool use_2D_matching = false;
0100   bool get_trigger_info = true;
0101   bool get_detector_info = true;
0102   bool get_dEdx_info = true;
0103   bool constrain_to_primary_vertex = true;
0104   bool use_pid = true;
0105   float pid_frac = 0.4;
0106 };  // namespace HeavyFlavorReco
0107
0108 //using namespace HeavyFlavorReco;
0109
0110 void create_hf_directories(const std::string& reconstruction_name, std::string &final_output_dir, std::string &output_reco_file)
0111 {
0112   std::string output_file_name = HeavyFlavorReco::kfp_header + reconstruction_name + HeavyFlavorReco::trailer;
0113   final_output_dir = HeavyFlavorReco::output_dir + reconstruction_name + "/";
0114   std::string output_reco_dir = final_output_dir + HeavyFlavorReco::processing_folder;
0115   output_reco_file = output_reco_dir + output_file_name;
0116
0117   std::string makeDirectory = "mkdir -p " + output_reco_dir;
0118   system(makeDirectory.c_str());
0119 }
0120
0121 void end_kfparticle(const std::string& full_file_name, const std::string& final_path)
0122 {
0123   std::ifstream file(full_file_name.c_str());
0124   if (file.good())
0125   {
0126     std::string moveOutput = "mv " + full_file_name + " " + final_path;
0127     system(moveOutput.c_str());
0128   }
0129 }
0130
0131 void Fun4All_raw_hit_KFP(
0132     const int nEvents = 10,
0133     const std::string& filelist = "filelist.list",
0134     const std::string& outfilename = "clusters_seeds",
0135     const bool convertSeeds = false,
0136     const int nSkip = 0,
0137     const bool doKFParticle = false)
0138 {
0139   auto *se = Fun4AllServer::instance();
0140   se->Verbosity(2);
0141   auto *rc = recoConsts::instance();
0142
0143   // input manager for QM production raw hit DST file
0144   std::ifstream ifs(filelist);
0145   std::string filepath;
0146
0147   int i = 0;
0148   int runnumber = std::numeric_limits<int>::quiet_NaN();
0149   int segment = std::numeric_limits<int>::quiet_NaN();
0150   bool process_endpoints = false;
0151
0152   while (std::getline(ifs, filepath))
0153   {
0154     std::cout << "Adding DST with filepath: " << filepath << std::endl;
0155     if (i == 0)
0156     {
0157       std::pair<int, int>
0158           runseg = Fun4AllUtils::GetRunSegment(filepath);
0159       runnumber = runseg.first;
0160       segment = runseg.second;
0161       rc->set_IntFlag("RUNNUMBER", runnumber);
0162       rc->set_uint64Flag("TIMESTAMP", runnumber);
0163     }
0164     if (filepath.find("ebdc") != std::string::npos)
0165     {
0166       if (filepath.find("_0_") != std::string::npos ||
0167           filepath.find("_1_") != std::string::npos)
0168       {
0169         process_endpoints = true;
0170       }
0171     }
0172     std::string inputname = "InputManager" + std::to_string(i);
0173     auto *hitsin = new Fun4AllDstInputManager(inputname);
0174     hitsin->fileopen(filepath);
0175     se->registerInputManager(hitsin);
0176     i++;
0177   }
0178
0179   rc->set_IntFlag("RUNNUMBER", runnumber);
0180   rc->set_IntFlag("RUNSEGMENT", segment);
0181
0182   Enable::QA = false;
0183   Enable::CDB = true;
0184   rc->set_StringFlag("CDB_GLOBALTAG", "newcdbtag");
0185   rc->set_uint64Flag("TIMESTAMP", runnumber);
0186
0187   std::stringstream nice_runnumber;
0188   nice_runnumber << std::setw(8) << std::setfill('0') << std::to_string(runnumber);
0189
0190   int rounded_up = 100 * (std::ceil((float) runnumber / 100));
0191   std::stringstream nice_rounded_up;
0192   nice_rounded_up << std::setw(8) << std::setfill('0') << std::to_string(rounded_up);
0193
0194   int rounded_down = 100 * (std::floor((float) runnumber / 100));
0195   std::stringstream nice_rounded_down;
0196   nice_rounded_down << std::setw(8) << std::setfill('0') << std::to_string(rounded_down);
0197
0198   std::stringstream nice_segment;
0199   nice_segment << std::setw(5) << std::setfill('0') << std::to_string(segment);
0200
0201   std::stringstream nice_skip;
0202   nice_skip << std::setw(5) << std::setfill('0') << std::to_string(nSkip);
0203
0204   HeavyFlavorReco::output_dir = "./";  // Top dir of where the output nTuples will be written
0205   HeavyFlavorReco::trailer = "_" + nice_runnumber.str() + "_" + nice_segment.str() + "_" + nice_skip.str() + ".root";
0206
0207   if (doKFParticle)
0208   {
0209     create_hf_directories(HeavyFlavorReco::pipi_reconstruction_name, HeavyFlavorReco::pipi_output_dir, HeavyFlavorReco::pipi_output_reco_file);
0210     create_hf_directories(HeavyFlavorReco::ppi_reconstruction_name, HeavyFlavorReco::ppi_output_dir, HeavyFlavorReco::ppi_output_reco_file);
0211   }
0212
0213   G4TRACKING::convert_seeds_to_svtxtracks = convertSeeds;
0214   std::cout << "Converting to seeds : " << G4TRACKING::convert_seeds_to_svtxtracks << std::endl;
0215
0216   std::cout << " run: " << runnumber
0217             << " samples: " << TRACKING::reco_tpc_maxtime_sample
0218             << " pre: " << TRACKING::reco_tpc_time_presample
0219             << " vdrift: " << G4TPC::tpc_drift_velocity_reco
0220             << std::endl;
0221
0222   TRACKING::pp_mode = false;
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
0230   TString outfile = outfilename + "_" + runnumber + "-" + segment + ".root";
0231   std::string theOutfile = outfile.Data();
0232
0233   FlagHandler *flag = new FlagHandler();
0234   se->registerSubsystem(flag);
0235
0236   std::string geofile = CDBInterface::instance()->getUrl("Tracking_Geometry");
0237
0238   Fun4AllRunNodeInputManager *ingeo = new Fun4AllRunNodeInputManager("GeoIn");
0239   ingeo->AddFile(geofile);
0240   se->registerInputManager(ingeo);
0241
0242   TpcReadoutInit(runnumber);
0243   G4TPC::REJECT_LASER_EVENTS = true;
0244   G4TPC::ENABLE_MODULE_EDGE_CORRECTIONS = true;
0245   // Flag for running the tpc hit unpacker with zero suppression on
0246   TRACKING::tpc_zero_supp = true;
0247
0248   // MVTX
0249   Enable::MVTX_APPLYMISALIGNMENT = true;
0250   ACTSGEOM::mvtx_applymisalignment = Enable::MVTX_APPLYMISALIGNMENT;
0251
0252   // to turn on the default static corrections, enable the two lines below
0253   G4TPC::ENABLE_STATIC_CORRECTIONS = true;
0254   G4TPC::USE_PHI_AS_RAD_STATIC_CORRECTIONS = false;
0255
0256   // to turn on the average corrections derived from simulation, enable the three lines below
0257   // note: these are designed to be used only if static corrections are also applied
0258   // G4TPC::ENABLE_AVERAGE_CORRECTIONS = true;
0259   // G4TPC::USE_PHI_AS_RAD_AVERAGE_CORRECTIONS = false;
0260   // G4TPC::average_correction_filename = std::string(getenv("CALIBRATIONROOT")) + "/distortion_maps/average_minus_static_distortion_inverted_10-new.root";
0261
0262   G4MAGNET::magfield_rescale = 1;
0263
0264   TrackingInit();
0265
0266   for (int felix = 0; felix < 6; felix++)
0267   {
0268     Mvtx_HitUnpacking(std::to_string(felix));
0269   }
0270   for (int server = 0; server < 8; server++)
0271   {
0272     Intt_HitUnpacking(std::to_string(server));
0273   }
0274   std::ostringstream ebdcname;
0275   for (int ebdc = 0; ebdc < 24; ebdc++)
0276   {
0277     if (!process_endpoints)
0278     {
0279       ebdcname.str("");
0280       if (ebdc < 10)
0281       {
0282         ebdcname << "0";
0283       }
0284       ebdcname << ebdc;
0285       Tpc_HitUnpacking(ebdcname.str());
0286     }
0287
0288     else if (process_endpoints)
0289     {
0290       for (int endpoint = 0; endpoint < 2; endpoint++)
0291       {
0292         ebdcname.str("");
0293         if (ebdc < 10)
0294         {
0295           ebdcname << "0";
0296         }
0297         ebdcname << ebdc << "_" << endpoint;
0298         Tpc_HitUnpacking(ebdcname.str());
0299       }
0300     }
0301   }
0302
0303   Micromegas_HitUnpacking();
0304
0305   Mvtx_Clustering();
0306
0307   Intt_Clustering();
0308
0309   Tpc_LaserEventIdentifying();
0310
0311   auto *tpcclusterizer = new TpcClusterizer;
0312   tpcclusterizer->Verbosity(0);
0313   tpcclusterizer->set_do_hit_association(G4TPC::DO_HIT_ASSOCIATION);
0314   tpcclusterizer->set_rawdata_reco();
0315   tpcclusterizer->set_reject_event(G4TPC::REJECT_LASER_EVENTS);
0316   se->registerSubsystem(tpcclusterizer);
0317
0318   Micromegas_Clustering();
0319
0320   Reject_Laser_Events();
0321   /*
0322    * Begin Track Seeding
0323    */
0324
0325   auto *silicon_Seeding = new PHActsSiliconSeeding;
0326   silicon_Seeding->Verbosity(0);
0327   silicon_Seeding->setStrobeRange(-5, 5);
0328   // these get us to about 83% INTT > 1
0329   silicon_Seeding->setinttRPhiSearchWindow(0.4);
0330   silicon_Seeding->setinttZSearchWindow(2.0);
0331   silicon_Seeding->seedAnalysis(false);
0332   se->registerSubsystem(silicon_Seeding);
0333
0334   auto *merger = new PHSiliconSeedMerger;
0335   merger->Verbosity(0);
0336   se->registerSubsystem(merger);
0337
0338   /*
0339    * Tpc Seeding
0340    */
0341   auto *seeder = new PHCASeeding("PHCASeeding");
0342   double fieldstrength = std::numeric_limits<double>::quiet_NaN();  // set by isConstantField if constant
0343   bool ConstField = isConstantField(G4MAGNET::magfield_tracking, fieldstrength);
0344   if (ConstField)
0345   {
0346     seeder->useConstBField(true);
0347     seeder->constBField(fieldstrength);
0348   }
0349   else
0350   {
0351     seeder->set_field_dir(-1 * G4MAGNET::magfield_rescale);
0352     seeder->useConstBField(false);
0353     seeder->magFieldFile(G4MAGNET::magfield_tracking);  // to get charge sign right
0354   }
0355   seeder->Verbosity(0);
0356   seeder->SetLayerRange(7, 55);
0357   seeder->SetSearchWindow(2., 0.05);           // z-width and phi-width, default in macro at 1.5 and 0.05
0358   seeder->SetClusAdd_delta_window(3.0, 0.06);  //  (0.5, 0.005) are default; sdzdr_cutoff, d2/dr2(phi)_cutoff
0359   // seeder->SetNClustersPerSeedRange(4,60); // default is 6, 6
0360   seeder->SetMinHitsPerCluster(0);
0361   seeder->SetMinClustersPerTrack(3);
0362   seeder->useFixedClusterError(true);
0363   seeder->set_pp_mode(true);
0364   se->registerSubsystem(seeder);
0365
0366   // expand stubs in the TPC using simple kalman filter
0367   auto *cprop = new PHSimpleKFProp("PHSimpleKFProp");
0368   cprop->set_field_dir(G4MAGNET::magfield_rescale);
0369   if (ConstField)
0370   {
0371     cprop->useConstBField(true);
0372     cprop->setConstBField(fieldstrength);
0373   }
0374   else
0375   {
0376     cprop->magFieldFile(G4MAGNET::magfield_tracking);
0377     cprop->set_field_dir(-1 * G4MAGNET::magfield_rescale);
0378   }
0379   cprop->useFixedClusterError(true);
0380   cprop->set_max_window(5.);
0381   cprop->Verbosity(0);
0382   cprop->set_pp_mode(true);
0383   se->registerSubsystem(cprop);
0384
0385   // Always apply preliminary distortion corrections to TPC clusters before silicon matching
0386   // and refit the trackseeds. Replace KFProp fits with the new fit parameters in the TPC seeds.
0387   auto *prelim_distcorr = new PrelimDistortionCorrection;
0388   prelim_distcorr->set_pp_mode(true);
0389   prelim_distcorr->Verbosity(0);
0390   se->registerSubsystem(prelim_distcorr);
0391
0392   /*
0393    * Track Matching between silicon and TPC
0394    */
0395   // The normal silicon association methods
0396   // Match the TPC track stubs from the CA seeder to silicon track stubs from PHSiliconTruthTrackSeeding
0397   auto *silicon_match = new PHSiliconTpcTrackMatching;
0398   silicon_match->Verbosity(0);
0399   silicon_match->set_pp_mode(TRACKING::pp_mode);
0400   if (G4TPC::ENABLE_AVERAGE_CORRECTIONS)
0401   {
0402     // for general tracking
0403     // Eta/Phi window is determined by 3 sigma window
0404     // X/Y/Z window is determined by 4 sigma window
0405     silicon_match->window_deta.set_posQoverpT_maxabs({-0.014, 0.0331, 0.48});
0406     silicon_match->window_deta.set_negQoverpT_maxabs({-0.006, 0.0235, 0.52});
0407     silicon_match->set_deltaeta_min(0.03);
0408     silicon_match->window_dphi.set_QoverpT_range({-0.15, 0, 0}, {0.15, 0, 0});
0409     silicon_match->window_dx.set_QoverpT_maxabs({3.0, 0, 0});
0410     silicon_match->window_dy.set_QoverpT_maxabs({3.0, 0, 0});
0411     silicon_match->window_dz.set_posQoverpT_maxabs({1.138, 0.3919, 0.84});
0412     silicon_match->window_dz.set_negQoverpT_maxabs({0.719, 0.6485, 0.65});
0413     silicon_match->set_crossing_deltaz_max(30);
0414     silicon_match->set_crossing_deltaz_min(2);
0415
0416     // for distortion correction using SI-TPOT fit and track pT>0.5
0417     if (G4TRACKING::SC_CALIBMODE)
0418     {
0419       silicon_match->window_deta.set_posQoverpT_maxabs({0.016, 0.0060, 1.13});
0420       silicon_match->window_deta.set_negQoverpT_maxabs({0.022, 0.0022, 1.44});
0421       silicon_match->set_deltaeta_min(0.03);
0422       silicon_match->window_dphi.set_QoverpT_range({-0.15, 0, 0}, {0.09, 0, 0});
0423       silicon_match->window_dx.set_QoverpT_maxabs({2.0, 0, 0});
0424       silicon_match->window_dy.set_QoverpT_maxabs({1.5, 0, 0});
0425       silicon_match->window_dz.set_posQoverpT_maxabs({1.213, 0.0211, 2.09});
0426       silicon_match->window_dz.set_negQoverpT_maxabs({1.307, 0.0001, 4.52});
0427       silicon_match->set_crossing_deltaz_min(1.2);
0428     }
0429   }
0430   se->registerSubsystem(silicon_match);
0431
0432   // Match TPC track stubs from CA seeder to clusters in the micromegas layers
0433   auto *mm_match = new PHMicromegasTpcTrackMatching;
0434   mm_match->Verbosity(0);
0435   mm_match->set_rphi_search_window_lyr1(3.);
0436   mm_match->set_rphi_search_window_lyr2(15.0);
0437   mm_match->set_z_search_window_lyr1(30.0);
0438   mm_match->set_z_search_window_lyr2(3.);
0439
0440   mm_match->set_min_tpc_layer(38);             // layer in TPC to start projection fit
0441   mm_match->set_test_windows_printout(false);  // used for tuning search windows only
0442   se->registerSubsystem(mm_match);
0443
0444   /*
0445    * End Track Seeding
0446    */
0447
0448   /*
0449    * Either converts seeds to tracks with a straight line/helix fit
0450    * or run the full Acts track kalman filter fit
0451    */
0452   if (G4TRACKING::convert_seeds_to_svtxtracks)
0453   {
0454     auto *converter = new TrackSeedTrackMapConverter;
0455     // Default set to full SvtxTrackSeeds. Can be set to
0456     // SiliconTrackSeedContainer or TpcTrackSeedContainer
0457     converter->setTrackSeedName("SvtxTrackSeedContainer");
0458     converter->setFieldMap(G4MAGNET::magfield_tracking);
0459     converter->Verbosity(0);
0460     se->registerSubsystem(converter);
0461   }
0462   else
0463   {
0464     auto *deltazcorr = new PHTpcDeltaZCorrection;
0465     deltazcorr->Verbosity(0);
0466     se->registerSubsystem(deltazcorr);
0467
0468     // perform final track fit with ACTS
0469     auto *actsFit = new PHActsTrkFitter;
0470     actsFit->Verbosity(0);
0471     actsFit->commissioning(G4TRACKING::use_alignment);
0472     // in calibration mode, fit only Silicons and Micromegas hits
0473     actsFit->fitSiliconMMs(G4TRACKING::SC_CALIBMODE);
0474     actsFit->setUseMicromegas(G4TRACKING::SC_USE_MICROMEGAS);
0475     actsFit->set_pp_mode(TRACKING::pp_mode);
0476     actsFit->set_use_clustermover(true);  // default is true for now
0477     actsFit->useActsEvaluator(false);
0478     actsFit->useOutlierFinder(false);
0479     actsFit->setFieldMap(G4MAGNET::magfield_tracking);
0480     se->registerSubsystem(actsFit);
0481
0482     auto *cleaner = new PHTrackCleaner();
0483     cleaner->Verbosity(0);
0484     cleaner->set_pp_mode(TRACKING::pp_mode);
0485     se->registerSubsystem(cleaner);
0486
0487     if (G4TRACKING::SC_CALIBMODE)
0488     {
0489       /*
0490        * in calibration mode, calculate residuals between TPC and fitted tracks,
0491        * store in dedicated structure for distortion correction
0492        */
0493       auto *residuals = new PHTpcResiduals;
0494       const TString tpc_residoutfile = theOutfile + "_PhTpcResiduals.root";
0495       residuals->setOutputfile(tpc_residoutfile.Data());
0496       residuals->setUseMicromegas(G4TRACKING::SC_USE_MICROMEGAS);
0497
0498       // matches Tony's analysis
0499       residuals->setMinPt(0.2);
0500
0501       // reconstructed distortion grid size (phi, r, z)
0502       residuals->setGridDimensions(36, 48, 80);
0503       se->registerSubsystem(residuals);
0504     }
0505   }
0506
0507   auto *finder = new PHSimpleVertexFinder;
0508   finder->Verbosity(0);
0509   finder->setDcaCut(0.5);
0510   finder->setTrackPtCut(0.3);
0511   finder->setBeamLineCut(1);
0512   finder->setTrackQualityCut(1000);
0513   finder->setNmvtxRequired(3);
0514   finder->setOutlierPairCut(0.1);
0515   se->registerSubsystem(finder);
0516
0517   if (!G4TRACKING::convert_seeds_to_svtxtracks)
0518   {
0519     // Propagate track positions to the vertex position
0520     auto *vtxProp = new PHActsVertexPropagator;
0521     vtxProp->Verbosity(0);
0522     vtxProp->fieldMap(G4MAGNET::magfield_tracking);
0523     se->registerSubsystem(vtxProp);
0524   }
0525
0526   TString residoutfile = theOutfile + "_resid.root";
0527   std::string residstring(residoutfile.Data());
0528
0529   auto *resid = new TrackResiduals("TrackResiduals");
0530   resid->outfileName(residstring);
0531   resid->alignment(false);
0532   resid->vertexTree();
0533   //   // adjust track map name
0534   //   if(G4TRACKING::SC_CALIBMODE && !G4TRACKING::convert_seeds_to_svtxtracks)
0535   //   {
0536   //     resid->trackmapName("SvtxSiliconMMTrackMap");
0537   //     if( G4TRACKING::SC_USE_MICROMEGAS )
0538   //     { resid->set_doMicromegasOnly(true); }
0539   //   }
0540
0541   //   resid->clusterTree();
0542   //   resid->hitTree();
0543   resid->convertSeeds(G4TRACKING::convert_seeds_to_svtxtracks);
0544
0545   //   resid->Verbosity(0);
0546   // se->registerSubsystem(resid);
0547
0548   // auto ntuplizer = new TrkrNtuplizer("TrkrNtuplizer");
0549   // se->registerSubsystem(ntuplizer);
0550
0551   /*
0552     // To write an output DST
0553     TString dstfile = theOutfile;
0554    std::string theDSTFile = dstfile.Data();
0555    Fun4AllOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", theDSTFile.c_str());
0556    out->AddNode("Sync");
0557    out->AddNode("EventHeader");
0558    out->AddNode("TRKR_CLUSTER");
0559    out->AddNode("SiliconTrackSeedContainer");
0560    out->AddNode("TpcTrackSeedContainer");
0561    out->AddNode("SvtxTrackSeedContainer");
0562    out->AddNode("SvtxTrackMap");
0563    out->AddNode("SvtxVertexMap");
0564    out->AddNode("MbdVertexMap");
0565    out->AddNode("GL1RAWHIT");
0566    se->registerOutputManager(out);
0567
0568   */
0569   if (Enable::QA)
0570   {
0571     se->registerSubsystem(new MvtxRawHitQA);
0572     se->registerSubsystem(new InttRawHitQA);
0573     se->registerSubsystem(new TpcRawHitQA);
0574     se->registerSubsystem(new MvtxClusterQA);
0575     se->registerSubsystem(new InttClusterQA);
0576     se->registerSubsystem(new TpcClusterQA);
0577     se->registerSubsystem(new MicromegasClusterQA);
0578
0579     auto *converter = new TrackSeedTrackMapConverter("SiliconSeedConverter");
0580     // Default set to full SvtxTrackSeeds. Can be set to
0581     // SiliconTrackSeedContainer or TpcTrackSeedContainer
0582     converter->setTrackSeedName("SiliconTrackSeedContainer");
0583     converter->setTrackMapName("SiliconSvtxTrackMap");
0584     converter->setFieldMap(G4MAGNET::magfield_tracking);
0585     converter->Verbosity(0);
0586     se->registerSubsystem(converter);
0587
0588     auto *finder_svx = new PHSimpleVertexFinder("SiliconVertexFinder");
0589     finder_svx->Verbosity(0);
0590     finder_svx->setDcaCut(0.1);
0591     finder_svx->setTrackPtCut(0.2);
0592     finder_svx->setBeamLineCut(1);
0593     finder_svx->setTrackQualityCut(500);
0594     finder_svx->setNmvtxRequired(3);
0595     finder_svx->setOutlierPairCut(0.1);
0596     finder_svx->setTrackMapName("SiliconSvtxTrackMap");
0597     finder_svx->setVertexMapName("SiliconSvtxVertexMap");
0598     se->registerSubsystem(finder_svx);
0599
0600     auto *siliconqa = new SiliconSeedsQA;
0601     siliconqa->setTrackMapName("SiliconSvtxTrackMap");
0602     siliconqa->setVertexMapName("SiliconSvtxVertexMap");
0603     se->registerSubsystem(siliconqa);
0604
0605     auto *convertertpc = new TrackSeedTrackMapConverter("TpcSeedConverter");
0606     // Default set to full SvtxTrackSeeds. Can be set to
0607     // SiliconTrackSeedContainer or TpcTrackSeedContainer
0608     convertertpc->setTrackSeedName("TpcTrackSeedContainer");
0609     convertertpc->setTrackMapName("TpcSvtxTrackMap");
0610     convertertpc->setFieldMap(G4MAGNET::magfield_tracking);
0611     convertertpc->Verbosity(0);
0612     se->registerSubsystem(convertertpc);
0613
0614     auto *findertpc = new PHSimpleVertexFinder("TpcSimpleVertexFinder");
0615     findertpc->Verbosity(0);
0616     findertpc->setDcaCut(0.5);
0617     findertpc->setTrackPtCut(0.2);
0618     findertpc->setBeamLineCut(1);
0619     findertpc->setTrackQualityCut(1000000000);
0620     // findertpc->setNmvtxRequired(3);
0621     findertpc->setRequireMVTX(false);
0622     findertpc->setOutlierPairCut(0.1);
0623     findertpc->setTrackMapName("TpcSvtxTrackMap");
0624     findertpc->setVertexMapName("TpcSvtxVertexMap");
0625     se->registerSubsystem(findertpc);
0626
0627     auto *tpcqa = new TpcSeedsQA;
0628     tpcqa->setTrackMapName("TpcSvtxTrackMap");
0629     tpcqa->setVertexMapName("TpcSvtxVertexMap");
0630     tpcqa->setSegment(rc->get_IntFlag("RUNSEGMENT"));
0631     se->registerSubsystem(tpcqa);
0632
0633     se->registerSubsystem(new TpcSiliconQA);
0634   }
0635
0636   if (doKFParticle)
0637   {
0638     // KFParticle dependancy
0639     Global_Reco();
0640
0641     // KFParticle setup
0642
0643     KFParticle_sPHENIX *kfparticle = new KFParticle_sPHENIX("pipi_reco");
0644     kfparticle->Verbosity(10);
0645     kfparticle->setDecayDescriptor("K_S0 -> pi^+ pi^-");
0646     // kfparticle->setDecayDescriptor("[K_S0 -> pi^+ pi^+]cc");
0647
0648     kfparticle->usePID(HeavyFlavorReco::use_pid);
0649     kfparticle->setPIDacceptFraction(HeavyFlavorReco::pid_frac);
0650     kfparticle->dontUseGlobalVertex(HeavyFlavorReco::dont_use_global_vertex);
0651     kfparticle->requireTrackVertexBunchCrossingMatch(HeavyFlavorReco::require_track_and_vertex_match);
0652     kfparticle->getAllPVInfo(HeavyFlavorReco::save_all_vtx_info);
0653     kfparticle->allowZeroMassTracks();
0654     kfparticle->use2Dmatching(HeavyFlavorReco::use_2D_matching);
0655     kfparticle->getTriggerInfo(HeavyFlavorReco::get_trigger_info);
0656     kfparticle->getDetectorInfo(HeavyFlavorReco::get_detector_info);
0657     kfparticle->saveDST(HeavyFlavorReco::save_tracks_to_DST);
0658     kfparticle->saveParticleContainer(false);
0659     kfparticle->magFieldFile("FIELDMAP_TRACKING");
0660
0661     // PV to SV cuts
0662     kfparticle->constrainToPrimaryVertex(HeavyFlavorReco::constrain_to_primary_vertex);
0663     kfparticle->setMotherIPchi2(100);
0664     kfparticle->setFlightDistancechi2(-1.);
0665     kfparticle->setMinDIRA(0.88);                    // was .95
0666     kfparticle->setDecayLengthRange(-0.1, FLT_MAX);  // was 0.1 min
0667
0668     kfparticle->setDecayLengthRange_XY(-10000, FLT_MAX);
0669     kfparticle->setDecayTimeRange_XY(-10000, FLT_MAX);
0670     kfparticle->setDecayTimeRange(-10000, FLT_MAX);
0671     kfparticle->setMinDecayTimeSignificance(-1e5);
0672     kfparticle->setMinDecayLengthSignificance(-1e5);
0673     kfparticle->setMinDecayLengthSignificance_XY(-1e5);
0674     kfparticle->setMaximumDaughterDCA_XY(100);
0675
0676     // Track parameters
0677     kfparticle->setMinimumTrackPT(0.0);
0678     kfparticle->setMinimumTrackIPchi2(-1.);
0679     kfparticle->setMinimumTrackIP(-1.);
0680     kfparticle->setMaximumTrackchi2nDOF(100.);
0681     kfparticle->setMinINTThits(0);
0682     kfparticle->setMinMVTXhits(0);
0683     kfparticle->setMinTPChits(20);
0684
0685     // Vertex parameters
0686     kfparticle->setMaximumVertexchi2nDOF(20);
0687     kfparticle->setMaximumDaughterDCA(0.5);  // 5 mm
0688
0689     // Parent parameters
0690     kfparticle->setMotherPT(0);
0691     kfparticle->setMinimumMass(0.40);  // Check mass ranges
0692     kfparticle->setMaximumMass(0.60);
0693     kfparticle->setMaximumMotherVertexVolume(0.1);
0694
0695     kfparticle->setOutputName(HeavyFlavorReco::pipi_output_reco_file);
0696
0697     se->registerSubsystem(kfparticle);
0698
0699     // Lambda reconstruction
0700     KFParticle_sPHENIX *kfparticleLambda = new KFParticle_sPHENIX("ppi_reco");
0701     kfparticleLambda->Verbosity(0);
0702     kfparticleLambda->setDecayDescriptor("[Lambda0 -> proton^+ pi^-]cc");
0703
0704     kfparticle->usePID(HeavyFlavorReco::use_pid);
0705     kfparticle->setPIDacceptFraction(HeavyFlavorReco::pid_frac);
0706     kfparticle->dontUseGlobalVertex(HeavyFlavorReco::dont_use_global_vertex);
0707     kfparticle->requireTrackVertexBunchCrossingMatch(HeavyFlavorReco::require_track_and_vertex_match);
0708     kfparticle->getAllPVInfo(HeavyFlavorReco::save_all_vtx_info);
0709     kfparticle->allowZeroMassTracks();
0710     kfparticle->use2Dmatching(HeavyFlavorReco::use_2D_matching);
0711     kfparticle->getTriggerInfo(HeavyFlavorReco::get_trigger_info);
0712     kfparticle->getDetectorInfo(HeavyFlavorReco::get_detector_info);
0713     kfparticle->saveDST(HeavyFlavorReco::save_tracks_to_DST);
0714     kfparticle->saveParticleContainer(false);
0715     kfparticle->magFieldFile("FIELDMAP_TRACKING");
0716
0717     // PV to SV cuts
0718     kfparticle->constrainToPrimaryVertex(HeavyFlavorReco::constrain_to_primary_vertex);
0719     kfparticle->setMotherIPchi2(100);
0720     kfparticle->setFlightDistancechi2(-1.);
0721     kfparticle->setMinDIRA(0.88);
0722     kfparticle->setDecayLengthRange(0.2, FLT_MAX);
0723
0724     // Track parameters
0725     kfparticle->setMinimumTrackPT(0.1);
0726     kfparticle->setMinimumTrackIPchi2(-1.);
0727     kfparticle->setMinimumTrackIP(-1.);
0728     kfparticle->setMaximumTrackchi2nDOF(100.);
0729     kfparticle->setMinTPChits(25);
0730
0731     // Vertex parameters
0732     kfparticle->setMaximumVertexchi2nDOF(20);
0733     kfparticle->setMaximumDaughterDCA(0.5);  // 5 mm
0734
0735     // Parent parameters
0736     kfparticle->setMotherPT(0);
0737     kfparticle->setMinimumMass(0.900);  // Check mass ranges
0738     kfparticle->setMaximumMass(1.300);
0739     kfparticle->setMaximumMotherVertexVolume(0.1);
0740
0741     kfparticle->setOutputName(HeavyFlavorReco::ppi_output_reco_file);
0742
0743     se->registerSubsystem(kfparticleLambda);
0744   }
0745
0746   se->skip(nSkip);
0747   se->run(nEvents);
0748   se->End();
0749   se->PrintTimer();
0750
0751   if (doKFParticle)
0752   {
0753     end_kfparticle(HeavyFlavorReco::pipi_output_reco_file, HeavyFlavorReco::pipi_output_dir);
0754     end_kfparticle(HeavyFlavorReco::ppi_output_reco_file, HeavyFlavorReco::ppi_output_dir);
0755   }
0756   if (Enable::QA)
0757   {
0758     std::string qaOutputFileName = theOutfile + "_qa.root";
0759     QAHistManagerDef::saveQARootFile(qaOutputFileName);
0760   }
0761   CDBInterface::instance()->Print();
0762   delete se;
0763   std::cout << "Finished" << std::endl;
0764   gSystem->Exit(0);
0765 }