TrackingProduction/run3auau/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_0 = "tpc00_0.list",
0102     const string &input_tpcfile00_1 = "tpc00_1.list",
0103     const string &input_tpcfile01_0 = "tpc01_0.list",
0104     const string &input_tpcfile01_1 = "tpc01_1.list",
0105     const string &input_tpcfile02_0 = "tpc02_0.list",
0106     const string &input_tpcfile02_1 = "tpc02_1.list",
0107     const string &input_tpcfile03_0 = "tpc03_0.list",
0108     const string &input_tpcfile03_1 = "tpc03_1.list",
0109     const string &input_tpcfile04_0 = "tpc04_0.list",
0110     const string &input_tpcfile04_1 = "tpc04_1.list",
0111     const string &input_tpcfile05_0 = "tpc05_0.list",
0112     const string &input_tpcfile05_1 = "tpc05_1.list",
0113     const string &input_tpcfile06_0 = "tpc06_0.list",
0114     const string &input_tpcfile06_1 = "tpc06_1.list",
0115     const string &input_tpcfile07_0 = "tpc07_0.list",
0116     const string &input_tpcfile07_1 = "tpc07_1.list",
0117     const string &input_tpcfile08_0 = "tpc08_0.list",
0118     const string &input_tpcfile08_1 = "tpc08_1.list",
0119     const string &input_tpcfile09_0 = "tpc09_0.list",
0120     const string &input_tpcfile09_1 = "tpc09_1.list",
0121     const string &input_tpcfile10_0 = "tpc10_0.list",
0122     const string &input_tpcfile10_1 = "tpc10_1.list",
0123     const string &input_tpcfile11_0 = "tpc11_0.list",
0124     const string &input_tpcfile11_1 = "tpc11_1.list",
0125     const string &input_tpcfile12_0 = "tpc12_0.list",
0126     const string &input_tpcfile12_1 = "tpc12_1.list",
0127     const string &input_tpcfile13_0 = "tpc13_0.list",
0128     const string &input_tpcfile13_1 = "tpc13_1.list",
0129     const string &input_tpcfile14_0 = "tpc14_0.list",
0130     const string &input_tpcfile14_1 = "tpc14_1.list",
0131     const string &input_tpcfile15_0 = "tpc15_0.list",
0132     const string &input_tpcfile15_1 = "tpc15_1.list",
0133     const string &input_tpcfile16_0 = "tpc16_0.list",
0134     const string &input_tpcfile16_1 = "tpc16_1.list",
0135     const string &input_tpcfile17_0 = "tpc17_0.list",
0136     const string &input_tpcfile17_1 = "tpc17_1.list",
0137     const string &input_tpcfile18_0 = "tpc18_0.list",
0138     const string &input_tpcfile18_1 = "tpc18_1.list",
0139     const string &input_tpcfile19_0 = "tpc19_0.list",
0140     const string &input_tpcfile19_1 = "tpc19_1.list",
0141     const string &input_tpcfile20_0 = "tpc20_0.list",
0142     const string &input_tpcfile20_1 = "tpc20_1.list",
0143     const string &input_tpcfile21_0 = "tpc21_0.list",
0144     const string &input_tpcfile21_1 = "tpc21_1.list",
0145     const string &input_tpcfile22_0 = "tpc22_0.list",
0146     const string &input_tpcfile22_1 = "tpc22_1.list",
0147     const string &input_tpcfile23_0 = "tpc23_0.list",
0148     const string &input_tpcfile23_1 = "tpc23_1.list",
0149     const string &input_tpotfile = "tpot.list"
0150 )
0151 {
0152  vector<string> gl1_infile;
0153   gl1_infile.push_back(input_gl1file);
0154
0155 // MVTX
0156   vector<string> mvtx_infile;
0157   mvtx_infile.push_back(input_mvtxfile00);
0158   mvtx_infile.push_back(input_mvtxfile01);
0159   mvtx_infile.push_back(input_mvtxfile02);
0160   mvtx_infile.push_back(input_mvtxfile03);
0161   mvtx_infile.push_back(input_mvtxfile04);
0162   mvtx_infile.push_back(input_mvtxfile05);
0163
0164 // INTT
0165   vector<string> intt_infile;
0166   intt_infile.push_back(input_inttfile00);
0167   intt_infile.push_back(input_inttfile01);
0168   intt_infile.push_back(input_inttfile02);
0169   intt_infile.push_back(input_inttfile03);
0170   intt_infile.push_back(input_inttfile04);
0171   intt_infile.push_back(input_inttfile05);
0172   intt_infile.push_back(input_inttfile06);
0173   intt_infile.push_back(input_inttfile07);
0174
0175   vector<string> tpc_infile;
0176   tpc_infile.push_back(input_tpcfile00_0);
0177   tpc_infile.push_back(input_tpcfile01_0);
0178   tpc_infile.push_back(input_tpcfile02_0);
0179   tpc_infile.push_back(input_tpcfile03_0);
0180   tpc_infile.push_back(input_tpcfile04_0);
0181   tpc_infile.push_back(input_tpcfile05_0);
0182   tpc_infile.push_back(input_tpcfile06_0);
0183   tpc_infile.push_back(input_tpcfile07_0);
0184   tpc_infile.push_back(input_tpcfile08_0);
0185   tpc_infile.push_back(input_tpcfile09_0);
0186   tpc_infile.push_back(input_tpcfile10_0);
0187   tpc_infile.push_back(input_tpcfile11_0);
0188   tpc_infile.push_back(input_tpcfile12_0);
0189   tpc_infile.push_back(input_tpcfile13_0);
0190   tpc_infile.push_back(input_tpcfile14_0);
0191   tpc_infile.push_back(input_tpcfile15_0);
0192   tpc_infile.push_back(input_tpcfile16_0);
0193   tpc_infile.push_back(input_tpcfile17_0);
0194   tpc_infile.push_back(input_tpcfile18_0);
0195   tpc_infile.push_back(input_tpcfile19_0);
0196   tpc_infile.push_back(input_tpcfile20_0);
0197   tpc_infile.push_back(input_tpcfile21_0);
0198   tpc_infile.push_back(input_tpcfile22_0);
0199   tpc_infile.push_back(input_tpcfile23_0);
0200   tpc_infile.push_back(input_tpcfile00_1);
0201   tpc_infile.push_back(input_tpcfile01_1);
0202   tpc_infile.push_back(input_tpcfile02_1);
0203   tpc_infile.push_back(input_tpcfile03_1);
0204   tpc_infile.push_back(input_tpcfile04_1);
0205   tpc_infile.push_back(input_tpcfile05_1);
0206   tpc_infile.push_back(input_tpcfile06_1);
0207   tpc_infile.push_back(input_tpcfile07_1);
0208   tpc_infile.push_back(input_tpcfile08_1);
0209   tpc_infile.push_back(input_tpcfile09_1);
0210   tpc_infile.push_back(input_tpcfile10_1);
0211   tpc_infile.push_back(input_tpcfile11_1);
0212   tpc_infile.push_back(input_tpcfile12_1);
0213   tpc_infile.push_back(input_tpcfile13_1);
0214   tpc_infile.push_back(input_tpcfile14_1);
0215   tpc_infile.push_back(input_tpcfile15_1);
0216   tpc_infile.push_back(input_tpcfile16_1);
0217   tpc_infile.push_back(input_tpcfile17_1);
0218   tpc_infile.push_back(input_tpcfile18_1);
0219   tpc_infile.push_back(input_tpcfile19_1);
0220   tpc_infile.push_back(input_tpcfile20_1);
0221   tpc_infile.push_back(input_tpcfile21_1);
0222   tpc_infile.push_back(input_tpcfile22_1);
0223   tpc_infile.push_back(input_tpcfile23_1);
0224
0225 // TPOT
0226   vector<string> tpot_infile;
0227   tpot_infile.push_back(input_tpotfile);
0228
0229  auto se = Fun4AllServer::instance();
0230   se->Verbosity(2);
0231   auto rc = recoConsts::instance();
0232   rc->set_IntFlag("RUNNUMBER", runnumber);
0233
0234   Enable::CDB = true;
0235   rc->set_StringFlag("CDB_GLOBALTAG", "ProdA_2024");
0236
0237   rc->set_uint64Flag("TIMESTAMP", runnumber);
0238
0239   //! flags to set
0240   Enable::QA = true;
0241   TRACKING::tpc_zero_supp = true;
0242   G4TRACKING::convert_seeds_to_svtxtracks = false;
0243
0244   Enable::MVTX_APPLYMISALIGNMENT = true;
0245   ACTSGEOM::mvtx_applymisalignment = Enable::MVTX_APPLYMISALIGNMENT;
0246   TpcReadoutInit( runnumber );
0247   std::cout<< " run: " << runnumber
0248        << " samples: " << TRACKING::reco_tpc_maxtime_sample
0249        << " pre: " << TRACKING::reco_tpc_time_presample
0250        << " vdrift: " << G4TPC::tpc_drift_velocity_reco
0251        << std::endl;
0252
0253   CDBInterface::instance()->Verbosity(1);
0254   std::string geofile = CDBInterface::instance()->getUrl("Tracking_Geometry");
0255   std::cout << "CDB tracking geometry file "<<geofile << std::endl;
0256   Fun4AllRunNodeInputManager *ingeo = new Fun4AllRunNodeInputManager("GeoIn");
0257   ingeo->AddFile(geofile);
0258   se->registerInputManager(ingeo);
0259
0260   CDBInterface *cdb = CDBInterface::instance();
0261
0262   // can use for zero field
0263   //double fieldstrength = 0.01;
0264   //G4MAGNET::magfield_tracking = "0.01";
0265   double fieldstrength = std::numeric_limits<double>::quiet_NaN();
0266   bool ConstField = isConstantField(G4MAGNET::magfield_tracking,fieldstrength);
0267
0268   if(ConstField && fieldstrength < 0.1)
0269   {
0270     G4MAGNET::magfield = "0.01";
0271     G4MAGNET::magfield_rescale = 1;
0272   }
0273
0274   TrackingInit();
0275
0276   int i= 0;
0277   int NumInputs = 0;
0278   Fun4AllStreamingInputManager *in = new Fun4AllStreamingInputManager("Comb");
0279
0280   for (auto iter : gl1_infile)
0281   {
0282     if (isGood(iter))
0283     {
0284       SingleGl1PoolInput *gl1_sngl = new SingleGl1PoolInput("GL1_" + to_string(i));
0285       //    gl1_sngl->Verbosity(3);
0286       gl1_sngl->AddListFile(iter);
0287       in->registerStreamingInput(gl1_sngl, InputManagerType::GL1);
0288       i++;
0289     }
0290   }
0291   NumInputs += i;
0292
0293   i = 0;
0294   for (auto iter : intt_infile)
0295   {
0296     if (isGood(iter))
0297     {
0298     SingleInttPoolInput *intt_sngl = new SingleInttPoolInput("INTT_" + to_string(i));
0299     //    intt_sngl->Verbosity(3);
0300     InttOdbcQuery query;
0301     bool isStreaming = false;
0302     if(runnumber != 0)
0303       {
0304     query.Query(runnumber);
0305     isStreaming = query.IsStreaming();
0306       }
0307     intt_sngl->streamingMode(isStreaming);
0308     intt_sngl->AddListFile(iter);
0309     in->registerStreamingInput(intt_sngl, InputManagerType::INTT);
0310     i++;
0311     }
0312   }
0313   NumInputs += i;
0314
0315   i = 0;
0316   for (auto iter : mvtx_infile)
0317   {
0318     if (isGood(iter))
0319     {
0320     SingleMvtxPoolInput *mvtx_sngl = new SingleMvtxPoolInput("MVTX_" + to_string(i));
0321     //    mvtx_sngl->Verbosity(3);
0322     mvtx_sngl->AddListFile(iter);
0323     in->registerStreamingInput(mvtx_sngl, InputManagerType::MVTX);
0324     i++;
0325     }
0326   }
0327   NumInputs += i;
0328
0329   i = 0;
0330
0331
0332
0333   for (auto iter : tpc_infile)
0334   {
0335     if (isGood(iter))
0336     {
0337
0338       SingleTpcTimeFrameInput *tpc_sngl = new SingleTpcTimeFrameInput("TPC_" + to_string(i));
0339 //    tpc_sngl->Verbosity(2);
0340       //   tpc_sngl->DryRun();
0341       tpc_sngl->setHitContainerName("TPCRAWHIT");
0342       tpc_sngl->AddListFile(iter);
0343       in->registerStreamingInput(tpc_sngl, InputManagerType::TPC);
0344       i++;
0345     }
0346   }
0347   NumInputs += i;
0348     i = 0;
0349   for (auto iter : tpot_infile)
0350   {
0351     if (isGood(iter))
0352     {
0353     SingleMicromegasPoolInput *mm_sngl = new SingleMicromegasPoolInput("MICROMEGAS_" + to_string(i));
0354     //   sngl->Verbosity(3);
0355     mm_sngl->SetBcoRange(10);
0356     mm_sngl->SetNegativeBco(2);
0357     mm_sngl->SetBcoPoolSize(50);
0358     mm_sngl->AddListFile(iter);
0359     in->registerStreamingInput(mm_sngl, InputManagerType::MICROMEGAS);
0360     i++;
0361     }
0362   }
0363   NumInputs += i;
0364
0365   // if there is no input manager this macro will still run - so just quit here
0366   if (NumInputs == 0)
0367   {
0368     std::cout << "no file lists no input manager registered, quitting" << std::endl;
0369     gSystem->Exit(1);
0370   }
0371   se->registerInputManager(in);
0372
0373   SyncReco *sync = new SyncReco();
0374   se->registerSubsystem(sync);
0375
0376   HeadReco *head = new HeadReco();
0377   se->registerSubsystem(head);
0378
0379   FlagHandler *flag = new FlagHandler();
0380   se->registerSubsystem(flag);
0381
0382
0383   Mvtx_HitUnpacking();
0384   Intt_HitUnpacking();
0385   Tpc_HitUnpacking();
0386   Micromegas_HitUnpacking();
0387
0388   MvtxClusterizer* mvtxclusterizer = new MvtxClusterizer("MvtxClusterizer");
0389
0390   se->registerSubsystem(mvtxclusterizer);
0391
0392   Intt_Clustering();
0393
0394   Tpc_LaserEventIdentifying();
0395
0396   auto tpcclusterizer = new TpcClusterizer;
0397   tpcclusterizer->Verbosity(0);
0398   tpcclusterizer->set_do_hit_association(G4TPC::DO_HIT_ASSOCIATION);
0399   tpcclusterizer->set_rawdata_reco();
0400   tpcclusterizer->set_reject_event(G4TPC::REJECT_LASER_EVENTS);
0401   se->registerSubsystem(tpcclusterizer);
0402
0403   Micromegas_Clustering();
0404
0405
0406   auto silicon_Seeding = new PHActsSiliconSeeding;
0407   silicon_Seeding->Verbosity(0);
0408   // these get us to about 83% INTT > 1
0409   silicon_Seeding->setinttRPhiSearchWindow(0.4);
0410   silicon_Seeding->setinttZSearchWindow(2.0);
0411   silicon_Seeding->setStrobeRange(-5,5);
0412   silicon_Seeding->seedAnalysis(false);
0413   se->registerSubsystem(silicon_Seeding);
0414
0415   auto merger = new PHSiliconSeedMerger;
0416   merger->Verbosity(0);
0417   se->registerSubsystem(merger);
0418
0419   /*
0420    * Tpc Seeding
0421    */
0422   auto seeder = new PHCASeeding("PHCASeeding");
0423  if (ConstField)
0424   {
0425     seeder->useConstBField(true);
0426     seeder->constBField(fieldstrength);
0427   }
0428   else
0429   {
0430     seeder->set_field_dir(-1 * G4MAGNET::magfield_rescale);
0431     seeder->useConstBField(false);
0432     seeder->magFieldFile(G4MAGNET::magfield_tracking);  // to get charge sign right
0433   }
0434   seeder->Verbosity(0);
0435   seeder->SetLayerRange(7, 55);
0436   seeder->SetSearchWindow(2.,0.05); // z-width and phi-width, default in macro at 1.5 and 0.05
0437   seeder->SetClusAdd_delta_window(3.0,0.06); //  (0.5, 0.005) are default; sdzdr_cutoff, d2/dr2(phi)_cutoff
0438   //seeder->SetNClustersPerSeedRange(4,60); // default is 6, 6
0439   seeder->SetMinHitsPerCluster(0);
0440   seeder->SetMinClustersPerTrack(3);
0441   seeder->useFixedClusterError(true);
0442   seeder->set_pp_mode(true);
0443   se->registerSubsystem(seeder);
0444
0445   // expand stubs in the TPC using simple kalman filter
0446   auto cprop = new PHSimpleKFProp("PHSimpleKFProp");
0447   cprop->set_field_dir(G4MAGNET::magfield_rescale);
0448   if (ConstField)
0449   {
0450     cprop->useConstBField(true);
0451     cprop->setConstBField(fieldstrength);
0452   }
0453   else
0454   {
0455     cprop->magFieldFile(G4MAGNET::magfield_tracking);
0456     cprop->set_field_dir(-1 * G4MAGNET::magfield_rescale);
0457   }
0458   cprop->useFixedClusterError(true);
0459   cprop->set_max_window(5.);
0460   cprop->Verbosity(0);
0461   cprop->set_pp_mode(true);
0462   se->registerSubsystem(cprop);
0463
0464   // Always apply preliminary distortion corrections to TPC clusters before silicon matching
0465   // and refit the trackseeds. Replace KFProp fits with the new fit parameters in the TPC seeds.
0466   auto prelim_distcorr = new PrelimDistortionCorrection;
0467   prelim_distcorr->set_pp_mode(true);
0468   prelim_distcorr->Verbosity(0);
0469   se->registerSubsystem(prelim_distcorr);
0470
0471   /*
0472    * Track Matching between silicon and TPC
0473    */
0474   // The normal silicon association methods
0475   // Match the TPC track stubs from the CA seeder to silicon track stubs from PHSiliconTruthTrackSeeding
0476   auto silicon_match = new PHSiliconTpcTrackMatching;
0477   silicon_match->Verbosity(0);
0478   silicon_match->set_pp_mode(TRACKING::pp_mode);
0479   if(G4TPC::ENABLE_AVERAGE_CORRECTIONS)
0480   {
0481     // for general tracking
0482     // Eta/Phi window is determined by 3 sigma window
0483     // X/Y/Z window is determined by 4 sigma window
0484     silicon_match->window_deta.set_posQoverpT_maxabs({-0.014,0.0331,0.48});
0485     silicon_match->window_deta.set_negQoverpT_maxabs({-0.006,0.0235,0.52});
0486     silicon_match->set_deltaeta_min(0.03);
0487     silicon_match->window_dphi.set_QoverpT_range({-0.15,0,0}, {0.15,0,0});
0488     silicon_match->window_dx.set_QoverpT_maxabs({3.0,0,0});
0489     silicon_match->window_dy.set_QoverpT_maxabs({3.0,0,0});
0490     silicon_match->window_dz.set_posQoverpT_maxabs({1.138,0.3919,0.84});
0491     silicon_match->window_dz.set_negQoverpT_maxabs({0.719,0.6485,0.65});
0492     silicon_match->set_crossing_deltaz_max(30);
0493     silicon_match->set_crossing_deltaz_min(2);
0494
0495     // for distortion correction using SI-TPOT fit and track pT>0.5
0496     if (G4TRACKING::SC_CALIBMODE)
0497     {
0498       silicon_match->window_deta.set_posQoverpT_maxabs({0.016,0.0060,1.13});
0499       silicon_match->window_deta.set_negQoverpT_maxabs({0.022,0.0022,1.44});
0500       silicon_match->set_deltaeta_min(0.03);
0501       silicon_match->window_dphi.set_QoverpT_range({-0.15,0,0}, {0.09,0,0});
0502       silicon_match->window_dx.set_QoverpT_maxabs({2.0,0,0});
0503       silicon_match->window_dy.set_QoverpT_maxabs({1.5,0,0});
0504       silicon_match->window_dz.set_posQoverpT_maxabs({1.213,0.0211,2.09});
0505       silicon_match->window_dz.set_negQoverpT_maxabs({1.307,0.0001,4.52});
0506       silicon_match->set_crossing_deltaz_min(1.2);
0507     }
0508   }
0509   se->registerSubsystem(silicon_match);
0510
0511   // Match TPC track stubs from CA seeder to clusters in the micromegas layers
0512   auto mm_match = new PHMicromegasTpcTrackMatching;
0513   mm_match->Verbosity(0);
0514   mm_match->set_rphi_search_window_lyr1(3.);
0515   mm_match->set_rphi_search_window_lyr2(15.0);
0516   mm_match->set_z_search_window_lyr1(30.0);
0517   mm_match->set_z_search_window_lyr2(3.);
0518
0519   mm_match->set_min_tpc_layer(38);             // layer in TPC to start projection fit
0520   mm_match->set_test_windows_printout(false);  // used for tuning search windows only
0521   se->registerSubsystem(mm_match);
0522
0523   /*
0524    * End Track Seeding
0525    */
0526
0527
0528   /*
0529    * Either converts seeds to tracks with a straight line/helix fit
0530    * or run the full Acts track kalman filter fit
0531    */
0532   if (G4TRACKING::convert_seeds_to_svtxtracks)
0533   {
0534     auto converter = new TrackSeedTrackMapConverter;
0535     // Default set to full SvtxTrackSeeds. Can be set to
0536     // SiliconTrackSeedContainer or TpcTrackSeedContainer
0537     converter->setTrackSeedName("SvtxTrackSeedContainer");
0538     converter->setFieldMap(G4MAGNET::magfield_tracking);
0539     converter->Verbosity(0);
0540     se->registerSubsystem(converter);
0541   }
0542   else
0543   {
0544     auto deltazcorr = new PHTpcDeltaZCorrection;
0545     deltazcorr->Verbosity(0);
0546     se->registerSubsystem(deltazcorr);
0547
0548     // perform final track fit with ACTS
0549     auto actsFit = new PHActsTrkFitter;
0550     actsFit->Verbosity(0);
0551     actsFit->commissioning(G4TRACKING::use_alignment);
0552     // in calibration mode, fit only Silicons and Micromegas hits
0553     actsFit->fitSiliconMMs(G4TRACKING::SC_CALIBMODE);
0554     actsFit->setUseMicromegas(G4TRACKING::SC_USE_MICROMEGAS);
0555     actsFit->set_pp_mode(TRACKING::pp_mode);
0556     actsFit->set_use_clustermover(true);  // default is true for now
0557     actsFit->useActsEvaluator(false);
0558     actsFit->useOutlierFinder(false);
0559     actsFit->setFieldMap(G4MAGNET::magfield_tracking);
0560     se->registerSubsystem(actsFit);
0561
0562     auto cleaner = new PHTrackCleaner();
0563     cleaner->Verbosity(0);
0564     cleaner->set_pp_mode(TRACKING::pp_mode);
0565     //se->registerSubsystem(cleaner);
0566   }
0567
0568
0569
0570   auto finder = new PHSimpleVertexFinder;
0571   finder->Verbosity(0);
0572   finder->setDcaCut(0.5);
0573   finder->setTrackPtCut(-99999.);
0574   finder->setBeamLineCut(1);
0575   finder->setTrackQualityCut(1000000000);
0576   finder->setNmvtxRequired(3);
0577   finder->setOutlierPairCut(0.1);
0578   se->registerSubsystem(finder);
0579
0580   TString residoutfile = outfilename + "_resid.root";
0581   std::string residstring(residoutfile.Data());
0582
0583   auto resid = new TrackResiduals("TrackResiduals");
0584   resid->Verbosity(0);
0585   resid->outfileName(residstring);
0586   resid->alignment(false);
0587   resid->clusterTree();
0588   //resid->failedTree();
0589   //resid->hitTree();
0590   resid->noEventTree();
0591   resid->convertSeeds(G4TRACKING::convert_seeds_to_svtxtracks);
0592
0593
0594   if(ConstField && fieldstrength < 0.1)
0595   {
0596     resid->zeroField();
0597   }
0598   se->registerSubsystem(resid);
0599
0600   // Fun4AllOutputManager *out = new Fun4AllDstOutputManager("out", "/sphenix/tg/tg01/hf/jdosbo/tracking_development/onlineoffline/hitsets.root");
0601   // se->registerOutputManager(out);
0602
0603   if (Enable::QA)
0604   {
0605     se->registerSubsystem(new MvtxRawHitQA);
0606     se->registerSubsystem(new InttRawHitQA);
0607     se->registerSubsystem(new TpcRawHitQA);
0608     se->registerSubsystem(new MvtxClusterQA);
0609     se->registerSubsystem(new InttClusterQA);
0610     se->registerSubsystem(new TpcClusterQA);
0611     se->registerSubsystem(new MicromegasClusterQA);
0612     se->registerSubsystem(new SiliconSeedsQA);
0613     se->registerSubsystem(new TpcSeedsQA);
0614     se->registerSubsystem(new TpcSiliconQA);
0615
0616   }
0617
0618   se->run(nEvents);
0619   se->End();
0620   se->PrintTimer();
0621
0622   if(Enable::QA)
0623     {
0624       TString qaname = outfilename + runnumber +"_qa.root";
0625       std::string qaOutputFileName(qaname.Data());
0626       QAHistManagerDef::saveQARootFile(qaOutputFileName);
0627     }
0628   CDBInterface::instance()->Print();
0629   delete se;
0630   std::cout << "Finished" << std::endl;
0631   gSystem->Exit(0);
0632 }
0633
0634
0635 bool isGood(const string &infile)
0636 {
0637   ifstream intest;
0638   intest.open(infile);
0639   bool goodfile = false;
0640   if (intest.is_open())
0641   {
0642     if (intest.peek() != std::ifstream::traits_type::eof()) // is it non zero?
0643     {
0644       goodfile = true;
0645     }
0646     intest.close();
0647   }
0648   return goodfile;
0649 }