macros/TrackingProduction/Fun4All_PRDFReconstruction.C

0001 /*
0002  * This macro shows running the full event combining + tracking for
0003  * cosmics running.. There are some analysis modules run at the end
0004  * which package  hits, clusters, and clusters on tracks into trees
0005  * for analysis.
0006  */
0007
0008 #include <G4_ActsGeom.C>
0009 #include <G4_Magnet.C>
0010 #include <GlobalVariables.C>
0011 #include <QA.C>
0012 #include <Trkr_Clustering.C>
0013 #include <Trkr_Reco.C>
0014 #include <Trkr_RecoInit.C>
0015 #include <Trkr_TpcReadoutInit.C>
0016
0017 #include <intt/InttOdbcQuery.h>
0018 #include <ffamodules/CDBInterface.h>
0019 #include <fun4allraw/Fun4AllStreamingInputManager.h>
0020 #include <fun4allraw/InputManagerType.h>
0021 #include <fun4allraw/SingleGl1PoolInput.h>
0022 #include <fun4allraw/SingleInttPoolInput.h>
0023 #include <fun4allraw/SingleMicromegasPoolInput.h>
0024 #include <fun4allraw/SingleMvtxPoolInput.h>
0025 #include <fun4all/Fun4AllDstInputManager.h>
0026 #include <fun4all/Fun4AllDstOutputManager.h>
0027 #include <fun4all/Fun4AllInputManager.h>
0028 #include <fun4all/Fun4AllOutputManager.h>
0029 #include <fun4all/Fun4AllRunNodeInputManager.h>
0030 #include <fun4all/Fun4AllServer.h>
0031
0032 #include <phool/recoConsts.h>
0033 #include <ffamodules/HeadReco.h>
0034 #include <ffamodules/FlagHandler.h>
0035 #include <ffamodules/SyncReco.h>
0036
0037 #include <fun4allraw/SingleTpcTimeFrameInput.h>
0038
0039 #include <cdbobjects/CDBTTree.h>
0040
0041 #include <mvtxrawhitqa/MvtxRawHitQA.h>
0042 #include <inttrawhitqa/InttRawHitQA.h>
0043 #include <tpcqa/TpcRawHitQA.h>
0044 #include <trackingqa/MvtxClusterQA.h>
0045 #include <trackingqa/InttClusterQA.h>
0046 #include <trackingqa/TpcClusterQA.h>
0047 #include <trackingqa/TpcSeedsQA.h>
0048 #include <trackingqa/SiliconSeedsQA.h>
0049 #include <trackingqa/TpcSiliconQA.h>
0050
0051 #include <trackingqa/MicromegasClusterQA.h>
0052
0053 #include <trackingqa/MvtxClusterQA.h>
0054 #include <trackreco/PHActsSiliconSeeding.h>
0055 #include <trackreco/PHSiliconSeedMerger.h>
0056
0057 #include <trackingqa/TpcClusterQA.h>
0058 #include <trackingqa/TpcSeedsQA.h>
0059
0060 #include <trackingdiagnostics/TrackResiduals.h>
0061 #include <trackingdiagnostics/TrkrNtuplizer.h>
0062
0063 #include <fun4all/Fun4AllUtils.h>
0064
0065 #include <stdio.h>
0066
0067 R__LOAD_LIBRARY(libfun4all.so)
0068 R__LOAD_LIBRARY(libffamodules.so)
0069 R__LOAD_LIBRARY(libphool.so)
0070 R__LOAD_LIBRARY(libcdbobjects.so)
0071 R__LOAD_LIBRARY(libmvtx.so)
0072 R__LOAD_LIBRARY(libintt.so)
0073 R__LOAD_LIBRARY(libtpc.so)
0074 R__LOAD_LIBRARY(libmicromegas.so)
0075 R__LOAD_LIBRARY(libTrackingDiagnostics.so)
0076 R__LOAD_LIBRARY(libtrackingqa.so)
0077 R__LOAD_LIBRARY(libtrack_reco.so)
0078
0079 bool isGood(const string &infile);
0080
0081 void Fun4All_PRDFReconstruction(
0082     const int nEvents = 0,
0083     const int runnumber = 53756,
0084     const std::string &dir = "/.",
0085     const std::string &outfilename = "output_tracks",
0086     const std::string &input_gl1file = "gl1daq.list",
0087     const string &input_inttfile00 = "intt0.list",
0088     const string &input_inttfile01 = "intt1.list",
0089     const string &input_inttfile02 = "intt2.list",
0090     const string &input_inttfile03 = "intt3.list",
0091     const string &input_inttfile04 = "intt4.list",
0092     const string &input_inttfile05 = "intt5.list",
0093     const string &input_inttfile06 = "intt6.list",
0094     const string &input_inttfile07 = "intt7.list",
0095     const string &input_mvtxfile00 = "mvtx0.list",
0096     const string &input_mvtxfile01 = "mvtx1.list",
0097     const string &input_mvtxfile02 = "mvtx2.list",
0098     const string &input_mvtxfile03 = "mvtx3.list",
0099     const string &input_mvtxfile04 = "mvtx4.list",
0100     const string &input_mvtxfile05 = "mvtx5.list",
0101     const string &input_tpcfile00 = "ebdc00.list",
0102     const string &input_tpcfile01 = "ebdc01.list",
0103     const string &input_tpcfile02 = "ebdc02.list",
0104     const string &input_tpcfile03 = "ebdc03.list",
0105     const string &input_tpcfile04 = "ebdc04.list",
0106     const string &input_tpcfile05 = "ebdc05.list",
0107     const string &input_tpcfile06 = "ebdc06.list",
0108     const string &input_tpcfile07 = "ebdc07.list",
0109     const string &input_tpcfile08 = "ebdc08.list",
0110     const string &input_tpcfile09 = "ebdc09.list",
0111     const string &input_tpcfile10 = "ebdc10.list",
0112     const string &input_tpcfile11 = "ebdc11.list",
0113     const string &input_tpcfile12 = "ebdc12.list",
0114     const string &input_tpcfile13 = "ebdc13.list",
0115     const string &input_tpcfile14 = "ebdc14.list",
0116     const string &input_tpcfile15 = "ebdc15.list",
0117     const string &input_tpcfile16 = "ebdc16.list",
0118     const string &input_tpcfile17 = "ebdc17.list",
0119     const string &input_tpcfile18 = "ebdc18.list",
0120     const string &input_tpcfile19 = "ebdc19.list",
0121     const string &input_tpcfile20 = "ebdc20.list",
0122     const string &input_tpcfile21 = "ebdc21.list",
0123     const string &input_tpcfile22 = "ebdc22.list",
0124     const string &input_tpcfile23 = "ebdc23.list",
0125     const string &input_tpotfile = "tpot.list"
0126 )
0127 {
0128  vector<string> gl1_infile;
0129   gl1_infile.push_back(input_gl1file);
0130
0131 // MVTX
0132   vector<string> mvtx_infile;
0133   mvtx_infile.push_back(input_mvtxfile00);
0134   mvtx_infile.push_back(input_mvtxfile01);
0135   mvtx_infile.push_back(input_mvtxfile02);
0136   mvtx_infile.push_back(input_mvtxfile03);
0137   mvtx_infile.push_back(input_mvtxfile04);
0138   mvtx_infile.push_back(input_mvtxfile05);
0139
0140 // INTT
0141   vector<string> intt_infile;
0142   intt_infile.push_back(input_inttfile00);
0143   intt_infile.push_back(input_inttfile01);
0144   intt_infile.push_back(input_inttfile02);
0145   intt_infile.push_back(input_inttfile03);
0146   intt_infile.push_back(input_inttfile04);
0147   intt_infile.push_back(input_inttfile05);
0148   intt_infile.push_back(input_inttfile06);
0149   intt_infile.push_back(input_inttfile07);
0150
0151   vector<string> tpc_infile;
0152   tpc_infile.push_back(input_tpcfile00);
0153   tpc_infile.push_back(input_tpcfile01);
0154   tpc_infile.push_back(input_tpcfile02);
0155   tpc_infile.push_back(input_tpcfile03);
0156   tpc_infile.push_back(input_tpcfile04);
0157   tpc_infile.push_back(input_tpcfile05);
0158   tpc_infile.push_back(input_tpcfile06);
0159   tpc_infile.push_back(input_tpcfile07);
0160   tpc_infile.push_back(input_tpcfile08);
0161   tpc_infile.push_back(input_tpcfile09);
0162   tpc_infile.push_back(input_tpcfile10);
0163   tpc_infile.push_back(input_tpcfile11);
0164   tpc_infile.push_back(input_tpcfile12);
0165   tpc_infile.push_back(input_tpcfile13);
0166   tpc_infile.push_back(input_tpcfile14);
0167   tpc_infile.push_back(input_tpcfile15);
0168   tpc_infile.push_back(input_tpcfile16);
0169   tpc_infile.push_back(input_tpcfile17);
0170   tpc_infile.push_back(input_tpcfile18);
0171   tpc_infile.push_back(input_tpcfile19);
0172   tpc_infile.push_back(input_tpcfile20);
0173   tpc_infile.push_back(input_tpcfile21);
0174   tpc_infile.push_back(input_tpcfile22);
0175   tpc_infile.push_back(input_tpcfile23);
0176
0177 // TPOT
0178   vector<string> tpot_infile;
0179   tpot_infile.push_back(input_tpotfile);
0180
0181  auto se = Fun4AllServer::instance();
0182   se->Verbosity(2);
0183   auto rc = recoConsts::instance();
0184   rc->set_IntFlag("RUNNUMBER", runnumber);
0185
0186   Enable::CDB = true;
0187   rc->set_StringFlag("CDB_GLOBALTAG", "ProdA_2024");
0188
0189   rc->set_uint64Flag("TIMESTAMP", runnumber);
0190
0191   //! flags to set
0192   Enable::QA = true;
0193   TRACKING::tpc_zero_supp = true;
0194   TRACKING::pp_mode = true;
0195   G4TRACKING::convert_seeds_to_svtxtracks = false;
0196
0197   Enable::MVTX_APPLYMISALIGNMENT = true;
0198   ACTSGEOM::mvtx_applymisalignment = Enable::MVTX_APPLYMISALIGNMENT;
0199   TpcReadoutInit( runnumber );
0200   std::cout<< " run: " << runnumber
0201        << " samples: " << TRACKING::reco_tpc_maxtime_sample
0202        << " pre: " << TRACKING::reco_tpc_time_presample
0203        << " vdrift: " << G4TPC::tpc_drift_velocity_reco
0204        << std::endl;
0205
0206   CDBInterface::instance()->Verbosity(1);
0207   std::string geofile = CDBInterface::instance()->getUrl("Tracking_Geometry");
0208   std::cout << "CDB tracking geometry file "<<geofile << std::endl;
0209   Fun4AllRunNodeInputManager *ingeo = new Fun4AllRunNodeInputManager("GeoIn");
0210   ingeo->AddFile(geofile);
0211   se->registerInputManager(ingeo);
0212
0213   CDBInterface *cdb = CDBInterface::instance();
0214
0215   // can use for zero field
0216   //double fieldstrength = 0.01;
0217   //G4MAGNET::magfield_tracking = "0.01";
0218   double fieldstrength = std::numeric_limits<double>::quiet_NaN();
0219   bool ConstField = isConstantField(G4MAGNET::magfield_tracking,fieldstrength);
0220
0221   if(ConstField && fieldstrength < 0.1)
0222   {
0223     G4MAGNET::magfield = "0.01";
0224     G4MAGNET::magfield_rescale = 1;
0225   }
0226
0227   TrackingInit();
0228
0229   int i= 0;
0230   int NumInputs = 0;
0231   Fun4AllStreamingInputManager *in = new Fun4AllStreamingInputManager("Comb");
0232
0233   for (auto iter : gl1_infile)
0234   {
0235     if (isGood(iter))
0236     {
0237       SingleGl1PoolInput *gl1_sngl = new SingleGl1PoolInput("GL1_" + to_string(i));
0238       //    gl1_sngl->Verbosity(3);
0239       gl1_sngl->AddListFile(iter);
0240       in->registerStreamingInput(gl1_sngl, InputManagerType::GL1);
0241       i++;
0242     }
0243   }
0244   NumInputs += i;
0245
0246   i = 0;
0247   for (auto iter : intt_infile)
0248   {
0249     if (isGood(iter))
0250     {
0251     SingleInttPoolInput *intt_sngl = new SingleInttPoolInput("INTT_" + to_string(i));
0252     //    intt_sngl->Verbosity(3);
0253     InttOdbcQuery query;
0254     bool isStreaming = false;
0255     if(runnumber != 0)
0256       {
0257     query.Query(runnumber);
0258     isStreaming = query.IsStreaming();
0259       }
0260     intt_sngl->streamingMode(isStreaming);
0261     intt_sngl->AddListFile(iter);
0262     in->registerStreamingInput(intt_sngl, InputManagerType::INTT);
0263     i++;
0264     }
0265   }
0266   NumInputs += i;
0267
0268   i = 0;
0269   for (auto iter : mvtx_infile)
0270   {
0271     if (isGood(iter))
0272     {
0273     SingleMvtxPoolInput *mvtx_sngl = new SingleMvtxPoolInput("MVTX_" + to_string(i));
0274     //    mvtx_sngl->Verbosity(3);
0275     mvtx_sngl->AddListFile(iter);
0276     in->registerStreamingInput(mvtx_sngl, InputManagerType::MVTX);
0277     i++;
0278     }
0279   }
0280   NumInputs += i;
0281
0282   i = 0;
0283
0284
0285
0286   for (auto iter : tpc_infile)
0287   {
0288     if (isGood(iter))
0289     {
0290
0291       SingleTpcTimeFrameInput *tpc_sngl = new SingleTpcTimeFrameInput("TPC_" + to_string(i));
0292 //    tpc_sngl->Verbosity(2);
0293       //   tpc_sngl->DryRun();
0294       tpc_sngl->setHitContainerName("TPCRAWHIT");
0295       tpc_sngl->AddListFile(iter);
0296       in->registerStreamingInput(tpc_sngl, InputManagerType::TPC);
0297       i++;
0298     }
0299   }
0300   NumInputs += i;
0301     i = 0;
0302   for (auto iter : tpot_infile)
0303   {
0304     if (isGood(iter))
0305     {
0306     SingleMicromegasPoolInput *mm_sngl = new SingleMicromegasPoolInput("MICROMEGAS_" + to_string(i));
0307     //   sngl->Verbosity(3);
0308     mm_sngl->SetBcoRange(10);
0309     mm_sngl->SetNegativeBco(2);
0310     mm_sngl->SetBcoPoolSize(50);
0311     mm_sngl->AddListFile(iter);
0312     in->registerStreamingInput(mm_sngl, InputManagerType::MICROMEGAS);
0313     i++;
0314     }
0315   }
0316   NumInputs += i;
0317
0318   // if there is no input manager this macro will still run - so just quit here
0319   if (NumInputs == 0)
0320   {
0321     std::cout << "no file lists no input manager registered, quitting" << std::endl;
0322     gSystem->Exit(1);
0323   }
0324   se->registerInputManager(in);
0325
0326   SyncReco *sync = new SyncReco();
0327   se->registerSubsystem(sync);
0328
0329   HeadReco *head = new HeadReco();
0330   se->registerSubsystem(head);
0331
0332   FlagHandler *flag = new FlagHandler();
0333   se->registerSubsystem(flag);
0334
0335
0336   Mvtx_HitUnpacking();
0337   Intt_HitUnpacking();
0338   Tpc_HitUnpacking();
0339   Micromegas_HitUnpacking();
0340
0341   MvtxClusterizer* mvtxclusterizer = new MvtxClusterizer("MvtxClusterizer");
0342   mvtxclusterizer->Verbosity(0);
0343   se->registerSubsystem(mvtxclusterizer);
0344
0345   Intt_Clustering();
0346
0347   Tpc_LaserEventIdentifying();
0348
0349   auto tpcclusterizer = new TpcClusterizer;
0350   tpcclusterizer->Verbosity(0);
0351   tpcclusterizer->set_do_hit_association(G4TPC::DO_HIT_ASSOCIATION);
0352   tpcclusterizer->set_rawdata_reco();
0353   tpcclusterizer->set_reject_event(G4TPC::REJECT_LASER_EVENTS);
0354   se->registerSubsystem(tpcclusterizer);
0355
0356   Micromegas_Clustering();
0357
0358
0359   auto silicon_Seeding = new PHActsSiliconSeeding;
0360   silicon_Seeding->Verbosity(0);
0361   // these get us to about 83% INTT > 1
0362   silicon_Seeding->setinttRPhiSearchWindow(0.4);
0363   silicon_Seeding->setinttZSearchWindow(2.0);
0364   silicon_Seeding->setStrobeRange(-5,5);
0365   silicon_Seeding->seedAnalysis(false);
0366   se->registerSubsystem(silicon_Seeding);
0367
0368   auto merger = new PHSiliconSeedMerger;
0369   merger->Verbosity(0);
0370   se->registerSubsystem(merger);
0371
0372   /*
0373    * Tpc Seeding
0374    */
0375   auto seeder = new PHCASeeding("PHCASeeding");
0376  if (ConstField)
0377   {
0378     seeder->useConstBField(true);
0379     seeder->constBField(fieldstrength);
0380   }
0381   else
0382   {
0383     seeder->set_field_dir(-1 * G4MAGNET::magfield_rescale);
0384     seeder->useConstBField(false);
0385     seeder->magFieldFile(G4MAGNET::magfield_tracking);  // to get charge sign right
0386   }
0387   seeder->Verbosity(0);
0388   seeder->SetLayerRange(7, 55);
0389   seeder->SetSearchWindow(2.,0.05); // z-width and phi-width, default in macro at 1.5 and 0.05
0390   seeder->SetClusAdd_delta_window(3.0,0.06); //  (0.5, 0.005) are default; sdzdr_cutoff, d2/dr2(phi)_cutoff
0391   //seeder->SetNClustersPerSeedRange(4,60); // default is 6, 6
0392   seeder->SetMinHitsPerCluster(0);
0393   seeder->SetMinClustersPerTrack(3);
0394   seeder->useFixedClusterError(true);
0395   seeder->set_pp_mode(true);
0396   se->registerSubsystem(seeder);
0397
0398   // expand stubs in the TPC using simple kalman filter
0399   auto cprop = new PHSimpleKFProp("PHSimpleKFProp");
0400   cprop->set_field_dir(G4MAGNET::magfield_rescale);
0401   if (ConstField)
0402   {
0403     cprop->useConstBField(true);
0404     cprop->setConstBField(fieldstrength);
0405   }
0406   else
0407   {
0408     cprop->magFieldFile(G4MAGNET::magfield_tracking);
0409     cprop->set_field_dir(-1 * G4MAGNET::magfield_rescale);
0410   }
0411   cprop->useFixedClusterError(true);
0412   cprop->set_max_window(5.);
0413   cprop->Verbosity(0);
0414   cprop->set_pp_mode(true);
0415   se->registerSubsystem(cprop);
0416
0417   // Always apply preliminary distortion corrections to TPC clusters before silicon matching
0418   // and refit the trackseeds. Replace KFProp fits with the new fit parameters in the TPC seeds.
0419   auto prelim_distcorr = new PrelimDistortionCorrection;
0420   prelim_distcorr->set_pp_mode(true);
0421   prelim_distcorr->Verbosity(0);
0422   se->registerSubsystem(prelim_distcorr);
0423
0424   /*
0425    * Track Matching between silicon and TPC
0426    */
0427   // The normal silicon association methods
0428   // Match the TPC track stubs from the CA seeder to silicon track stubs from PHSiliconTruthTrackSeeding
0429   auto silicon_match = new PHSiliconTpcTrackMatching;
0430   silicon_match->Verbosity(0);
0431   silicon_match->set_pp_mode(TRACKING::pp_mode);
0432   if(G4TPC::ENABLE_AVERAGE_CORRECTIONS)
0433   {
0434     // for general tracking
0435     // Eta/Phi window is determined by 3 sigma window
0436     // X/Y/Z window is determined by 4 sigma window
0437     silicon_match->window_deta.set_posQoverpT_maxabs({-0.014,0.0331,0.48});
0438     silicon_match->window_deta.set_negQoverpT_maxabs({-0.006,0.0235,0.52});
0439     silicon_match->set_deltaeta_min(0.03);
0440     silicon_match->window_dphi.set_QoverpT_range({-0.15,0,0}, {0.15,0,0});
0441     silicon_match->window_dx.set_QoverpT_maxabs({3.0,0,0});
0442     silicon_match->window_dy.set_QoverpT_maxabs({3.0,0,0});
0443     silicon_match->window_dz.set_posQoverpT_maxabs({1.138,0.3919,0.84});
0444     silicon_match->window_dz.set_negQoverpT_maxabs({0.719,0.6485,0.65});
0445     silicon_match->set_crossing_deltaz_max(30);
0446     silicon_match->set_crossing_deltaz_min(2);
0447
0448     // for distortion correction using SI-TPOT fit and track pT>0.5
0449     if (G4TRACKING::SC_CALIBMODE)
0450     {
0451       silicon_match->window_deta.set_posQoverpT_maxabs({0.016,0.0060,1.13});
0452       silicon_match->window_deta.set_negQoverpT_maxabs({0.022,0.0022,1.44});
0453       silicon_match->set_deltaeta_min(0.03);
0454       silicon_match->window_dphi.set_QoverpT_range({-0.15,0,0}, {0.09,0,0});
0455       silicon_match->window_dx.set_QoverpT_maxabs({2.0,0,0});
0456       silicon_match->window_dy.set_QoverpT_maxabs({1.5,0,0});
0457       silicon_match->window_dz.set_posQoverpT_maxabs({1.213,0.0211,2.09});
0458       silicon_match->window_dz.set_negQoverpT_maxabs({1.307,0.0001,4.52});
0459       silicon_match->set_crossing_deltaz_min(1.2);
0460     }
0461   }
0462   se->registerSubsystem(silicon_match);
0463
0464   // Match TPC track stubs from CA seeder to clusters in the micromegas layers
0465   auto mm_match = new PHMicromegasTpcTrackMatching;
0466   mm_match->Verbosity(0);
0467   mm_match->set_rphi_search_window_lyr1(3.);
0468   mm_match->set_rphi_search_window_lyr2(15.0);
0469   mm_match->set_z_search_window_lyr1(30.0);
0470   mm_match->set_z_search_window_lyr2(3.);
0471
0472   mm_match->set_min_tpc_layer(38);             // layer in TPC to start projection fit
0473   mm_match->set_test_windows_printout(false);  // used for tuning search windows only
0474   se->registerSubsystem(mm_match);
0475
0476   /*
0477    * End Track Seeding
0478    */
0479
0480
0481   /*
0482    * Either converts seeds to tracks with a straight line/helix fit
0483    * or run the full Acts track kalman filter fit
0484    */
0485   if (G4TRACKING::convert_seeds_to_svtxtracks)
0486   {
0487     auto converter = new TrackSeedTrackMapConverter;
0488     // Default set to full SvtxTrackSeeds. Can be set to
0489     // SiliconTrackSeedContainer or TpcTrackSeedContainer
0490     converter->setTrackSeedName("SvtxTrackSeedContainer");
0491     converter->setFieldMap(G4MAGNET::magfield_tracking);
0492     converter->Verbosity(0);
0493     se->registerSubsystem(converter);
0494   }
0495   else
0496   {
0497     auto deltazcorr = new PHTpcDeltaZCorrection;
0498     deltazcorr->Verbosity(0);
0499     se->registerSubsystem(deltazcorr);
0500
0501     // perform final track fit with ACTS
0502     auto actsFit = new PHActsTrkFitter;
0503     actsFit->Verbosity(0);
0504     actsFit->commissioning(G4TRACKING::use_alignment);
0505     // in calibration mode, fit only Silicons and Micromegas hits
0506     actsFit->fitSiliconMMs(G4TRACKING::SC_CALIBMODE);
0507     actsFit->setUseMicromegas(G4TRACKING::SC_USE_MICROMEGAS);
0508     actsFit->set_pp_mode(TRACKING::pp_mode);
0509     actsFit->set_use_clustermover(true);  // default is true for now
0510     actsFit->useActsEvaluator(false);
0511     actsFit->useOutlierFinder(false);
0512     actsFit->setFieldMap(G4MAGNET::magfield_tracking);
0513     se->registerSubsystem(actsFit);
0514
0515     auto cleaner = new PHTrackCleaner();
0516     cleaner->Verbosity(0);
0517     cleaner->set_pp_mode(TRACKING::pp_mode);
0518     //se->registerSubsystem(cleaner);
0519   }
0520
0521
0522
0523   auto finder = new PHSimpleVertexFinder;
0524   finder->Verbosity(0);
0525   finder->setDcaCut(0.5);
0526   finder->setTrackPtCut(-99999.);
0527   finder->setBeamLineCut(1);
0528   finder->setTrackQualityCut(1000000000);
0529   finder->setNmvtxRequired(3);
0530   finder->setOutlierPairCut(0.1);
0531   se->registerSubsystem(finder);
0532
0533   // Propagate track positions to the vertex position
0534   auto vtxProp = new PHActsVertexPropagator;
0535   vtxProp->Verbosity(0);
0536   vtxProp->fieldMap(G4MAGNET::magfield_tracking);
0537   se->registerSubsystem(vtxProp);
0538
0539   TString residoutfile = outfilename + "_resid.root";
0540   std::string residstring(residoutfile.Data());
0541
0542   auto resid = new TrackResiduals("TrackResiduals");
0543   resid->Verbosity(0);
0544   resid->outfileName(residstring);
0545   resid->alignment(false);
0546   //resid->hitTree();
0547   resid->convertSeeds(G4TRACKING::convert_seeds_to_svtxtracks);
0548
0549
0550   if(ConstField && fieldstrength < 0.1)
0551   {
0552     resid->zeroField();
0553   }
0554   //se->registerSubsystem(resid);
0555
0556   // Fun4AllOutputManager *out = new Fun4AllDstOutputManager("out", "/sphenix/tg/tg01/hf/jdosbo/tracking_development/onlineoffline/hitsets.root");
0557   // se->registerOutputManager(out);
0558
0559   if (Enable::QA)
0560   {
0561     se->registerSubsystem(new MvtxRawHitQA);
0562     se->registerSubsystem(new InttRawHitQA);
0563     se->registerSubsystem(new TpcRawHitQA);
0564     se->registerSubsystem(new MvtxClusterQA);
0565     se->registerSubsystem(new InttClusterQA);
0566     se->registerSubsystem(new TpcClusterQA);
0567     se->registerSubsystem(new MicromegasClusterQA);
0568     se->registerSubsystem(new SiliconSeedsQA);
0569     se->registerSubsystem(new TpcSeedsQA);
0570     se->registerSubsystem(new TpcSiliconQA);
0571
0572   }
0573
0574   se->run(nEvents);
0575   se->End();
0576   se->PrintTimer();
0577
0578   if(Enable::QA)
0579     {
0580       TString qaname = outfilename + runnumber +"_qa.root";
0581       std::string qaOutputFileName(qaname.Data());
0582       QAHistManagerDef::saveQARootFile(qaOutputFileName);
0583     }
0584   CDBInterface::instance()->Print();
0585   delete se;
0586   std::cout << "Finished" << std::endl;
0587   gSystem->Exit(0);
0588 }
0589
0590
0591 bool isGood(const string &infile)
0592 {
0593   ifstream intest;
0594   intest.open(infile);
0595   bool goodfile = false;
0596   if (intest.is_open())
0597   {
0598     if (intest.peek() != std::ifstream::traits_type::eof()) // is it non zero?
0599     {
0600       goodfile = true;
0601     }
0602     intest.close();
0603   }
0604   return goodfile;
0605 }