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