Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-03 08:20:27 

0001 #ifndef MACRO_TRKRRECO_C
0002 #define MACRO_TRKRRECO_C
0003 
0004 #include <G4_TrkrVariables.C>
0005 
0006 #include <trackingdiagnostics/TrackSeedTrackMapConverter.h>
0007 
0008 #include <trackreco/MakeActsGeometry.h>
0009 #include <trackreco/PHActsSiliconSeeding.h>
0010 #include <trackreco/PHActsTrackProjection.h>
0011 #include <trackreco/PHActsTrkFitter.h>
0012 #include <trackreco/PHActsVertexPropagator.h>
0013 #include <trackreco/PHCASeeding.h>
0014 #include <trackreco/PHGenFitTrkFitter.h>
0015 #include <trackreco/PHMicromegasTpcTrackMatching.h>
0016 #include <trackreco/PHSiliconHelicalPropagator.h>
0017 #include <trackreco/PHSiliconSeedMerger.h>
0018 #include <trackreco/PHSiliconTpcTrackMatching.h>
0019 #include <trackreco/PHSimpleKFProp.h>
0020 #include <trackreco/PHSimpleVertexFinder.h>
0021 #include <trackreco/PHTpcDeltaZCorrection.h>
0022 #include <trackreco/PHTrackCleaner.h>
0023 #include <trackreco/PHTrackSeeding.h>
0024 #include <trackreco/PrelimDistortionCorrection.h>
0025 #include <trackreco/SecondaryVertexFinder.h>
0026 #include <trackreco/TrackingIterationCounter.h>
0027 
0028 #include <tpc/TpcLoadDistortionCorrection.h>
0029 #include <tpc/LaserEventRejecter.h>
0030 
0031 #include <tpccalib/PHTpcResiduals.h>
0032 #include <tpccalib/TpcSpaceChargeReconstruction.h>
0033 
0034 #include <trackermillepedealignment/HelicalFitter.h>
0035 #include <trackermillepedealignment/MakeMilleFiles.h>
0036 
0037 #include <fun4all/Fun4AllServer.h>
0038 
0039 #include <trackingdiagnostics/TrackContainerCombiner.h>
0040 
0041 #include <string>
0042 
0043 R__LOAD_LIBRARY(libTrackingDiagnostics.so)
0044 R__LOAD_LIBRARY(libtrack_reco.so)
0045 R__LOAD_LIBRARY(libtpccalib.so)
0046 R__LOAD_LIBRARY(libtpc.so)
0047 R__LOAD_LIBRARY(libtrackeralign.so)
0048 
0049 void convert_seeds()
0050 {
0051   auto se = Fun4AllServer::instance();
0052   int verbosity = std::max(Enable::VERBOSITY, Enable::TRACKING_VERBOSITY);
0053 
0054   auto converter = new TrackSeedTrackMapConverter;
0055   // Default set to full SvtxTrackSeeds. Can be set to
0056   // SiliconTrackSeedContainer or TpcTrackSeedContainer
0057   converter->setTrackSeedName("SvtxTrackSeedContainer");
0058   converter->setFieldMap(G4MAGNET::magfield_tracking);
0059   converter->Verbosity(verbosity);
0060   se->registerSubsystem(converter);
0061 }
0062 void Tracking_Reco_TrackSeed_ZeroField()
0063 {
0064   // set up verbosity
0065   int verbosity = std::max(Enable::VERBOSITY, Enable::TRACKING_VERBOSITY);
0066 
0067   // get fun4all server instance
0068   auto se = Fun4AllServer::instance();
0069 
0070  auto silicon_Seeding = new PHActsSiliconSeeding;
0071   silicon_Seeding->Verbosity(verbosity);
0072   se->registerSubsystem(silicon_Seeding);
0073 
0074   auto merger = new PHSiliconSeedMerger;
0075   merger->Verbosity(verbosity);
0076   se->registerSubsystem(merger);
0077 
0078   auto seeder = new PHCASeeding("PHCASeeding");
0079   double fieldstrength = std::numeric_limits<double>::quiet_NaN();  // set by isConstantField if constant
0080   bool ConstField = isConstantField(G4MAGNET::magfield_tracking, fieldstrength);
0081   if (ConstField)
0082   {
0083     seeder->useConstBField(true);
0084     seeder->constBField(fieldstrength);
0085   }
0086   else
0087   {
0088     seeder->set_field_dir(-1 * G4MAGNET::magfield_rescale);
0089     seeder->useConstBField(false);
0090     seeder->magFieldFile(G4MAGNET::magfield_tracking);  // to get charge sign right
0091   }
0092   seeder->Verbosity(verbosity);
0093   seeder->SetLayerRange(7, 55);
0094   seeder->SetSearchWindow(1.5, 0.05);  // (z width, phi width)
0095   seeder->SetMinHitsPerCluster(0);
0096   seeder->SetMinClustersPerTrack(3);
0097   seeder->useFixedClusterError(true);
0098 
0099   if (G4TPC::TPC_GAS_MIXTURE == "NeCF4")
0100   {
0101     seeder->setNeonFraction(G4TPC::NeCF4_Ne_frac);
0102     seeder->setArgonFraction(G4TPC::NeCF4_Ar_frac);
0103     seeder->setCF4Fraction(G4TPC::NeCF4_CF4_frac);
0104     seeder->setNitrogenFraction(G4TPC::NeCF4_N2_frac);
0105     seeder->setIsobutaneFraction(G4TPC::NeCF4_isobutane_frac);
0106   }
0107   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4")
0108   {
0109     seeder->setNeonFraction(G4TPC::ArCF4_Ne_frac);
0110     seeder->setArgonFraction(G4TPC::ArCF4_Ar_frac);
0111     seeder->setCF4Fraction(G4TPC::ArCF4_CF4_frac);
0112     seeder->setNitrogenFraction(G4TPC::ArCF4_N2_frac);
0113     seeder->setIsobutaneFraction(G4TPC::ArCF4_isobutane_frac);
0114   }
0115   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4N2")
0116   {
0117     seeder->setNeonFraction(G4TPC::ArCF4N2_Ne_frac);
0118     seeder->setArgonFraction(G4TPC::ArCF4N2_Ar_frac);
0119     seeder->setCF4Fraction(G4TPC::ArCF4N2_CF4_frac);
0120     seeder->setNitrogenFraction(G4TPC::ArCF4N2_N2_frac);
0121     seeder->setIsobutaneFraction(G4TPC::ArCF4N2_isobutane_frac);
0122   }
0123   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4Isobutane")
0124   {
0125     seeder->setNeonFraction(G4TPC::ArCF4Isobutane_Ne_frac);
0126     seeder->setArgonFraction(G4TPC::ArCF4Isobutane_Ar_frac);
0127     seeder->setCF4Fraction(G4TPC::ArCF4Isobutane_CF4_frac);
0128     seeder->setNitrogenFraction(G4TPC::ArCF4Isobutane_N2_frac);
0129     seeder->setIsobutaneFraction(G4TPC::ArCF4Isobutane_isobutane_frac);
0130   }
0131   else
0132   {
0133   }
0134 
0135   seeder->set_pp_mode(true);
0136   se->registerSubsystem(seeder);
0137 
0138   // expand stubs in the TPC using simple kalman filter
0139   auto cprop = new PHSimpleKFProp("PHSimpleKFProp");
0140   cprop->set_field_dir(G4MAGNET::magfield_rescale);
0141   if (ConstField)
0142   {
0143     cprop->useConstBField(true);
0144     cprop->setConstBField(fieldstrength);
0145   }
0146   else
0147   {
0148     cprop->magFieldFile(G4MAGNET::magfield_tracking);
0149     cprop->set_field_dir(-1 * G4MAGNET::magfield_rescale);
0150   }
0151   cprop->useFixedClusterError(true);
0152   cprop->set_max_window(5.);
0153   cprop->Verbosity(verbosity);
0154 
0155   if (G4TPC::TPC_GAS_MIXTURE == "NeCF4")
0156   {
0157     cprop->setNeonFraction(G4TPC::NeCF4_Ne_frac);
0158     cprop->setArgonFraction(G4TPC::NeCF4_Ar_frac);
0159     cprop->setCF4Fraction(G4TPC::NeCF4_CF4_frac);
0160     cprop->setNitrogenFraction(G4TPC::NeCF4_N2_frac);
0161     cprop->setIsobutaneFraction(G4TPC::NeCF4_isobutane_frac);
0162   }
0163   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4")
0164   {
0165     cprop->setNeonFraction(G4TPC::ArCF4_Ne_frac);
0166     cprop->setArgonFraction(G4TPC::ArCF4_Ar_frac);
0167     cprop->setCF4Fraction(G4TPC::ArCF4_CF4_frac);
0168     cprop->setNitrogenFraction(G4TPC::ArCF4_N2_frac);
0169     cprop->setIsobutaneFraction(G4TPC::ArCF4_isobutane_frac);
0170   }
0171   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4N2")
0172   {
0173     cprop->setNeonFraction(G4TPC::ArCF4N2_Ne_frac);
0174     cprop->setArgonFraction(G4TPC::ArCF4N2_Ar_frac);
0175     cprop->setCF4Fraction(G4TPC::ArCF4N2_CF4_frac);
0176     cprop->setNitrogenFraction(G4TPC::ArCF4N2_N2_frac);
0177     cprop->setIsobutaneFraction(G4TPC::ArCF4N2_isobutane_frac);
0178   }
0179   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4Isobutane")
0180   {
0181     cprop->setNeonFraction(G4TPC::ArCF4Isobutane_Ne_frac);
0182     cprop->setArgonFraction(G4TPC::ArCF4Isobutane_Ar_frac);
0183     cprop->setCF4Fraction(G4TPC::ArCF4Isobutane_CF4_frac);
0184     cprop->setNitrogenFraction(G4TPC::ArCF4Isobutane_N2_frac);
0185     cprop->setIsobutaneFraction(G4TPC::ArCF4Isobutane_isobutane_frac);
0186   }
0187   else
0188   {
0189   }
0190 
0191   cprop->set_pp_mode(true);
0192   se->registerSubsystem(cprop);
0193 
0194 }
0195 void Tracking_Reco_TrackSeed()
0196 {
0197   // set up verbosity
0198   int verbosity = std::max(Enable::VERBOSITY, Enable::TRACKING_VERBOSITY);
0199 
0200   // get fun4all server instance
0201   auto se = Fun4AllServer::instance();
0202 
0203   // Assemble silicon clusters into track stubs
0204 
0205   auto silicon_Seeding = new PHActsSiliconSeeding;
0206   silicon_Seeding->Verbosity(verbosity);
0207 
0208   // modify strobe range
0209   silicon_Seeding->setStrobeRange(0,2);
0210 
0211   se->registerSubsystem(silicon_Seeding);
0212 
0213   auto merger = new PHSiliconSeedMerger;
0214   merger->Verbosity(verbosity);
0215   se->registerSubsystem(merger);
0216 
0217   auto seeder = new PHCASeeding("PHCASeeding");
0218   double fieldstrength = std::numeric_limits<double>::quiet_NaN();  // set by isConstantField if constant
0219   bool ConstField = isConstantField(G4MAGNET::magfield_tracking, fieldstrength);
0220   if (ConstField)
0221   {
0222     seeder->useConstBField(true);
0223     seeder->constBField(fieldstrength);
0224   }
0225   else
0226   {
0227     seeder->set_field_dir(-1 * G4MAGNET::magfield_rescale);
0228     seeder->useConstBField(false);
0229     seeder->magFieldFile(G4MAGNET::magfield_tracking);  // to get charge sign right
0230   }
0231   seeder->Verbosity(verbosity);
0232   seeder->SetLayerRange(7, 55);
0233   seeder->SetSearchWindow(1.5, 0.05);  // (z width, phi width)
0234   seeder->SetMinHitsPerCluster(0);
0235   seeder->SetMinClustersPerTrack(3);
0236   seeder->useFixedClusterError(true);
0237 
0238   if (G4TPC::TPC_GAS_MIXTURE == "NeCF4")
0239   {
0240     seeder->setNeonFraction(G4TPC::NeCF4_Ne_frac);
0241     seeder->setArgonFraction(G4TPC::NeCF4_Ar_frac);
0242     seeder->setCF4Fraction(G4TPC::NeCF4_CF4_frac);
0243     seeder->setNitrogenFraction(G4TPC::NeCF4_N2_frac);
0244     seeder->setIsobutaneFraction(G4TPC::NeCF4_isobutane_frac);
0245   }
0246   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4")
0247   {
0248     seeder->setNeonFraction(G4TPC::ArCF4_Ne_frac);
0249     seeder->setArgonFraction(G4TPC::ArCF4_Ar_frac);
0250     seeder->setCF4Fraction(G4TPC::ArCF4_CF4_frac);
0251     seeder->setNitrogenFraction(G4TPC::ArCF4_N2_frac);
0252     seeder->setIsobutaneFraction(G4TPC::ArCF4_isobutane_frac);
0253   }
0254   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4N2")
0255   {
0256     seeder->setNeonFraction(G4TPC::ArCF4N2_Ne_frac);
0257     seeder->setArgonFraction(G4TPC::ArCF4N2_Ar_frac);
0258     seeder->setCF4Fraction(G4TPC::ArCF4N2_CF4_frac);
0259     seeder->setNitrogenFraction(G4TPC::ArCF4N2_N2_frac);
0260     seeder->setIsobutaneFraction(G4TPC::ArCF4N2_isobutane_frac);
0261   }
0262   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4Isobutane")
0263   {
0264     seeder->setNeonFraction(G4TPC::ArCF4Isobutane_Ne_frac);
0265     seeder->setArgonFraction(G4TPC::ArCF4Isobutane_Ar_frac);
0266     seeder->setCF4Fraction(G4TPC::ArCF4Isobutane_CF4_frac);
0267     seeder->setNitrogenFraction(G4TPC::ArCF4Isobutane_N2_frac);
0268     seeder->setIsobutaneFraction(G4TPC::ArCF4Isobutane_isobutane_frac);
0269   }
0270   else
0271   {
0272   }
0273 
0274   seeder->set_pp_mode(TRACKING::pp_mode);
0275   se->registerSubsystem(seeder);
0276 
0277   // expand stubs in the TPC using simple kalman filter
0278   auto cprop = new PHSimpleKFProp("PHSimpleKFProp");
0279   cprop->set_field_dir(G4MAGNET::magfield_rescale);
0280   if (ConstField)
0281   {
0282     cprop->useConstBField(true);
0283     cprop->setConstBField(fieldstrength);
0284   }
0285   else
0286   {
0287     cprop->magFieldFile(G4MAGNET::magfield_tracking);
0288     cprop->set_field_dir(-1 * G4MAGNET::magfield_rescale);
0289   }
0290   cprop->useFixedClusterError(true);
0291   cprop->set_max_window(5.);
0292   cprop->Verbosity(verbosity);
0293 
0294   if (G4TPC::TPC_GAS_MIXTURE == "NeCF4")
0295   {
0296     cprop->setNeonFraction(G4TPC::NeCF4_Ne_frac);
0297     cprop->setArgonFraction(G4TPC::NeCF4_Ar_frac);
0298     cprop->setCF4Fraction(G4TPC::NeCF4_CF4_frac);
0299     cprop->setNitrogenFraction(G4TPC::NeCF4_N2_frac);
0300     cprop->setIsobutaneFraction(G4TPC::NeCF4_isobutane_frac);
0301   }
0302   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4")
0303   {
0304     cprop->setNeonFraction(G4TPC::ArCF4_Ne_frac);
0305     cprop->setArgonFraction(G4TPC::ArCF4_Ar_frac);
0306     cprop->setCF4Fraction(G4TPC::ArCF4_CF4_frac);
0307     cprop->setNitrogenFraction(G4TPC::ArCF4_N2_frac);
0308     cprop->setIsobutaneFraction(G4TPC::ArCF4_isobutane_frac);
0309   }
0310   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4N2")
0311   {
0312     cprop->setNeonFraction(G4TPC::ArCF4N2_Ne_frac);
0313     cprop->setArgonFraction(G4TPC::ArCF4N2_Ar_frac);
0314     cprop->setCF4Fraction(G4TPC::ArCF4N2_CF4_frac);
0315     cprop->setNitrogenFraction(G4TPC::ArCF4N2_N2_frac);
0316     cprop->setIsobutaneFraction(G4TPC::ArCF4N2_isobutane_frac);
0317   }
0318   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4Isobutane")
0319   {
0320     cprop->setNeonFraction(G4TPC::ArCF4Isobutane_Ne_frac);
0321     cprop->setArgonFraction(G4TPC::ArCF4Isobutane_Ar_frac);
0322     cprop->setCF4Fraction(G4TPC::ArCF4Isobutane_CF4_frac);
0323     cprop->setNitrogenFraction(G4TPC::ArCF4Isobutane_N2_frac);
0324     cprop->setIsobutaneFraction(G4TPC::ArCF4Isobutane_isobutane_frac);
0325   }
0326   else
0327   {
0328   }
0329 
0330   cprop->set_pp_mode(TRACKING::pp_mode);
0331   se->registerSubsystem(cprop);
0332 
0333   if (TRACKING::pp_mode)
0334   {
0335     // for pp mode, apply preliminary distortion corrections to TPC clusters before crossing is known
0336     // and refit the trackseeds. Replace KFProp fits with the new fit parameters in the TPC seeds.
0337     auto prelim_distcorr = new PrelimDistortionCorrection;
0338     prelim_distcorr->set_pp_mode(TRACKING::pp_mode);
0339     prelim_distcorr->Verbosity(verbosity);
0340 
0341     if (G4TPC::TPC_GAS_MIXTURE == "NeCF4")
0342     {
0343       prelim_distcorr->setNeonFraction(G4TPC::NeCF4_Ne_frac);
0344       prelim_distcorr->setArgonFraction(G4TPC::NeCF4_Ar_frac);
0345       prelim_distcorr->setCF4Fraction(G4TPC::NeCF4_CF4_frac);
0346       prelim_distcorr->setNitrogenFraction(G4TPC::NeCF4_N2_frac);
0347       prelim_distcorr->setIsobutaneFraction(G4TPC::NeCF4_isobutane_frac);
0348     }
0349     else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4")
0350     {
0351       prelim_distcorr->setNeonFraction(G4TPC::ArCF4_Ne_frac);
0352       prelim_distcorr->setArgonFraction(G4TPC::ArCF4_Ar_frac);
0353       prelim_distcorr->setCF4Fraction(G4TPC::ArCF4_CF4_frac);
0354       prelim_distcorr->setNitrogenFraction(G4TPC::ArCF4_N2_frac);
0355       prelim_distcorr->setIsobutaneFraction(G4TPC::ArCF4_isobutane_frac);
0356     }
0357     else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4N2")
0358     {
0359       prelim_distcorr->setNeonFraction(G4TPC::ArCF4N2_Ne_frac);
0360       prelim_distcorr->setArgonFraction(G4TPC::ArCF4N2_Ar_frac);
0361       prelim_distcorr->setCF4Fraction(G4TPC::ArCF4N2_CF4_frac);
0362       prelim_distcorr->setNitrogenFraction(G4TPC::ArCF4N2_N2_frac);
0363       prelim_distcorr->setIsobutaneFraction(G4TPC::ArCF4N2_isobutane_frac);
0364     }
0365     else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4Isobutane")
0366     {
0367       prelim_distcorr->setNeonFraction(G4TPC::ArCF4Isobutane_Ne_frac);
0368       prelim_distcorr->setArgonFraction(G4TPC::ArCF4Isobutane_Ar_frac);
0369       prelim_distcorr->setCF4Fraction(G4TPC::ArCF4Isobutane_CF4_frac);
0370       prelim_distcorr->setNitrogenFraction(G4TPC::ArCF4Isobutane_N2_frac);
0371       prelim_distcorr->setIsobutaneFraction(G4TPC::ArCF4Isobutane_isobutane_frac);
0372     }
0373     else
0374     {
0375     }
0376 
0377     se->registerSubsystem(prelim_distcorr);
0378   }
0379 
0380   std::cout << "Tracking_Reco_TrackSeed - Using stub matching for Si matching " << std::endl;
0381   // The normal silicon association methods
0382   // Match the TPC track stubs from the CA seeder to silicon track stubs from PHSiliconTruthTrackSeeding
0383   auto silicon_match = new PHSiliconTpcTrackMatching;
0384   silicon_match->Verbosity(verbosity);
0385   silicon_match->set_pp_mode(TRACKING::pp_mode);
0386   std::cout << "PHSiliconTpcTrackMatching pp_mode set to " << TRACKING::pp_mode << std::endl;
0387   if (G4TRACKING::SC_CALIBMODE)
0388   {
0389     // search windows for initial matching with distortions
0390     // tuned values are 0.04 and 0.008 in distorted events
0391     silicon_match->set_phi_search_window(0.04);
0392     silicon_match->set_eta_search_window(0.008);
0393   }
0394   else
0395   {
0396     // after distortion corrections and rerunning clustering, default tuned values are 0.02 and 0.004 in low occupancy events
0397     silicon_match->set_phi_search_window(0.03);
0398     silicon_match->set_eta_search_window(0.005);
0399   }
0400   silicon_match->set_test_windows_printout(false);  // used for tuning search windows
0401   se->registerSubsystem(silicon_match);
0402 
0403   // Associate Micromegas clusters with the tracks
0404   if (Enable::MICROMEGAS)
0405   {
0406     std::cout << "Tracking_Reco_TrackSeed - Using Micromegas matching " << std::endl;
0407 
0408     // Match TPC track stubs from CA seeder to clusters in the micromegas layers
0409     auto mm_match = new PHMicromegasTpcTrackMatching;
0410     mm_match->Verbosity(verbosity);
0411     mm_match->set_pp_mode(TRACKING::pp_mode);
0412     mm_match->set_rphi_search_window_lyr1(0.2);
0413     mm_match->set_rphi_search_window_lyr2(13.0);
0414     mm_match->set_z_search_window_lyr1(26.0);
0415     mm_match->set_z_search_window_lyr2(0.2);
0416 
0417     mm_match->set_min_tpc_layer(38);             // layer in TPC to start projection fit
0418     mm_match->set_test_windows_printout(false);  // used for tuning search windows only
0419     se->registerSubsystem(mm_match);
0420   }
0421 }
0422 
0423 void Tracking_Reco_TrackSeed_pass1()
0424 {
0425   Fun4AllServer* se = Fun4AllServer::instance();
0426   int verbosity = std::max(Enable::VERBOSITY, Enable::TRACKING_VERBOSITY);
0427 
0428   TrackingIterationCounter* counter = new TrackingIterationCounter("TrkrIter1");
0429   /// Clusters already used are in the 0th iteration
0430   counter->iteration(0);
0431   se->registerSubsystem(counter);
0432 
0433   PHActsSiliconSeeding* silseed = new PHActsSiliconSeeding("PHActsSiliconSeedingIt1");
0434   silseed->Verbosity(verbosity);
0435   silseed->searchInIntt();
0436   silseed->iteration(1);
0437   silseed->set_track_map_name("SiliconTrackSeedContainerIt1");
0438   se->registerSubsystem(silseed);
0439 
0440   PHSiliconSeedMerger* merger = new PHSiliconSeedMerger("SiliconSeedMargerIt1");
0441   merger->Verbosity(verbosity);
0442   merger->clusterOverlap(2);
0443   merger->searchIntt();
0444   merger->trackMapName("SiliconTrackSeedContainerIt1");
0445   se->registerSubsystem(merger);
0446 
0447   TrackContainerCombiner* combiner = new TrackContainerCombiner;
0448   combiner->newContainerName("SiliconTrackSeedContainer");
0449   combiner->oldContainerName("SiliconTrackSeedContainerIt1");
0450   combiner->Verbosity(verbosity);
0451   se->registerSubsystem(combiner);
0452 }
0453 
0454 void vertexing()
0455 {
0456   Fun4AllServer* se = Fun4AllServer::instance();
0457   int verbosity = std::max(Enable::VERBOSITY, Enable::TRACKING_VERBOSITY);
0458 
0459   auto vtxfinder = new PHSimpleVertexFinder;
0460   vtxfinder->Verbosity(verbosity);
0461   se->registerSubsystem(vtxfinder);
0462 }
0463 
0464 void Tracking_Reco_TrackFit()
0465 {
0466   int verbosity = std::max(Enable::VERBOSITY, Enable::TRACKING_VERBOSITY);
0467   auto se = Fun4AllServer::instance();
0468 
0469   // correct clusters for particle propagation in TPC
0470   auto deltazcorr = new PHTpcDeltaZCorrection;
0471   deltazcorr->Verbosity(verbosity);
0472   se->registerSubsystem(deltazcorr);
0473 
0474   if (G4TRACKING::use_genfit_track_fitter)
0475   {
0476     // perform final track fit with GENFIT
0477     auto genfitFit = new PHGenFitTrkFitter;
0478     genfitFit->Verbosity(verbosity);
0479     genfitFit->set_fit_silicon_mms(G4TRACKING::SC_CALIBMODE);
0480     se->registerSubsystem(genfitFit);
0481 
0482     if (G4TRACKING::SC_CALIBMODE)
0483     {
0484       // Genfit based Tpc space charge Reconstruction
0485       auto tpcSpaceChargeReconstruction = new TpcSpaceChargeReconstruction;
0486       tpcSpaceChargeReconstruction->set_use_micromegas(G4TRACKING::SC_USE_MICROMEGAS);
0487       tpcSpaceChargeReconstruction->set_outputfile(G4TRACKING::SC_ROOTOUTPUT_FILENAME);
0488       // reconstructed distortion grid size (phi, r, z)
0489       tpcSpaceChargeReconstruction->set_grid_dimensions(36, 48, 80);
0490       se->registerSubsystem(tpcSpaceChargeReconstruction);
0491     }
0492   }
0493   else
0494   {
0495     // perform final track fit with ACTS
0496     auto actsFit = new PHActsTrkFitter;
0497     actsFit->Verbosity(verbosity);
0498     actsFit->commissioning(G4TRACKING::use_alignment);
0499     // in calibration mode, fit only Silicons and Micromegas hits
0500     actsFit->fitSiliconMMs(G4TRACKING::SC_CALIBMODE);
0501     actsFit->setUseMicromegas(G4TRACKING::SC_USE_MICROMEGAS);
0502     actsFit->set_pp_mode(TRACKING::pp_mode);
0503     actsFit->set_use_clustermover(true);  // default is true for now
0504     actsFit->useActsEvaluator(false);
0505     actsFit->useOutlierFinder(false);
0506     actsFit->setFieldMap(G4MAGNET::magfield_tracking);
0507     se->registerSubsystem(actsFit);
0508 
0509     if (G4TRACKING::SC_CALIBMODE)
0510     {
0511       /*
0512        * in calibration mode, calculate residuals between TPC and fitted tracks,
0513        * store in dedicated structure for distortion correction
0514        */
0515       auto residuals = new PHTpcResiduals();
0516       residuals->setOutputfile(G4TRACKING::SC_ROOTOUTPUT_FILENAME);
0517       residuals->setUseMicromegas(G4TRACKING::SC_USE_MICROMEGAS);
0518       // reconstructed distortion grid size (phi, r, z)
0519       residuals->setGridDimensions(36, 48, 80);
0520       residuals->Verbosity(verbosity);
0521       se->registerSubsystem(residuals);
0522     }
0523   }
0524 
0525   if (!G4TRACKING::SC_CALIBMODE)
0526   {
0527     /*
0528      * in full tracking mode, run track cleaner, vertex finder,
0529      * propagete tracks to vertex
0530      * propagate tracks to EMCAL
0531      */
0532 
0533     if (!G4TRACKING::use_full_truth_track_seeding)
0534     {
0535       // Choose the best silicon matched track for each TPC track seed
0536       /* this breaks in truth_track seeding mode because there is no TpcSeed */
0537       auto cleaner = new PHTrackCleaner;
0538       cleaner->Verbosity(verbosity);
0539       // cleaner->set_quality_cut(30.0);
0540       se->registerSubsystem(cleaner);
0541     }
0542 
0543     vertexing();
0544 
0545     if (!G4TRACKING::use_genfit_track_fitter)
0546     {
0547       // Propagate track positions to the vertex position
0548       auto vtxProp = new PHActsVertexPropagator;
0549       vtxProp->Verbosity(verbosity);
0550       vtxProp->fieldMap(G4MAGNET::magfield_tracking);
0551       se->registerSubsystem(vtxProp);
0552 
0553       // project tracks to EMCAL
0554       auto projection = new PHActsTrackProjection;
0555       projection->Verbosity(verbosity);
0556       double fieldstrength = std::numeric_limits<double>::quiet_NaN();
0557       if (isConstantField(G4MAGNET::magfield_tracking, fieldstrength))
0558       {
0559         projection->setConstFieldVal(fieldstrength);
0560       }
0561       se->registerSubsystem(projection);
0562     }
0563   }
0564 }
0565 
0566 void Tracking_Reco_CommissioningTrackSeed()
0567 {
0568   // set up verbosity
0569   int verbosity = std::max(Enable::VERBOSITY, Enable::TRACKING_VERBOSITY);
0570 
0571   // get fun4all server instance
0572   auto se = Fun4AllServer::instance();
0573 
0574   auto silicon_Seeding = new PHActsSiliconSeeding;
0575   silicon_Seeding->Verbosity(verbosity);
0576   silicon_Seeding->sigmaScattering(50.);
0577   silicon_Seeding->setinttRPhiSearchWindow(2.);
0578   silicon_Seeding->helixcut(0.01);
0579   se->registerSubsystem(silicon_Seeding);
0580 
0581   auto merger = new PHSiliconSeedMerger;
0582   merger->Verbosity(verbosity);
0583   se->registerSubsystem(merger);
0584 
0585   // Assemble TPC clusters into track stubs
0586   auto seeder = new PHCASeeding("PHCASeeding");
0587   seeder->set_field_dir(G4MAGNET::magfield_rescale);                // to get charge sign right
0588   double fieldstrength = std::numeric_limits<double>::quiet_NaN();  // set by isConstantField if constant
0589   bool ConstField = isConstantField(G4MAGNET::magfield_tracking, fieldstrength);
0590   if (!ConstField)
0591   {
0592     seeder->magFieldFile(G4MAGNET::magfield_tracking);
0593     seeder->set_field_dir(-1 * G4MAGNET::magfield_rescale);
0594   }
0595   seeder->Verbosity(verbosity);
0596   seeder->SetLayerRange(7, 55);
0597   seeder->SetSearchWindow(1.5, 0.05);  // (z width, phi width)
0598   seeder->SetMinHitsPerCluster(0);
0599   seeder->SetMinClustersPerTrack(3);
0600   seeder->useConstBField(false);
0601   seeder->useFixedClusterError(true);
0602 
0603   if (G4TPC::TPC_GAS_MIXTURE == "NeCF4")
0604   {
0605     seeder->setNeonFraction(G4TPC::NeCF4_Ne_frac);
0606     seeder->setArgonFraction(G4TPC::NeCF4_Ar_frac);
0607     seeder->setCF4Fraction(G4TPC::NeCF4_CF4_frac);
0608     seeder->setNitrogenFraction(G4TPC::NeCF4_N2_frac);
0609     seeder->setIsobutaneFraction(G4TPC::NeCF4_isobutane_frac);
0610   }
0611   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4")
0612   {
0613     seeder->setNeonFraction(G4TPC::ArCF4_Ne_frac);
0614     seeder->setArgonFraction(G4TPC::ArCF4_Ar_frac);
0615     seeder->setCF4Fraction(G4TPC::ArCF4_CF4_frac);
0616     seeder->setNitrogenFraction(G4TPC::ArCF4_N2_frac);
0617     seeder->setIsobutaneFraction(G4TPC::ArCF4_isobutane_frac);
0618   }
0619   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4N2")
0620   {
0621     seeder->setNeonFraction(G4TPC::ArCF4N2_Ne_frac);
0622     seeder->setArgonFraction(G4TPC::ArCF4N2_Ar_frac);
0623     seeder->setCF4Fraction(G4TPC::ArCF4N2_CF4_frac);
0624     seeder->setNitrogenFraction(G4TPC::ArCF4N2_N2_frac);
0625     seeder->setIsobutaneFraction(G4TPC::ArCF4N2_isobutane_frac);
0626   }
0627   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4Isobutane")
0628   {
0629     seeder->setNeonFraction(G4TPC::ArCF4Isobutane_Ne_frac);
0630     seeder->setArgonFraction(G4TPC::ArCF4Isobutane_Ar_frac);
0631     seeder->setCF4Fraction(G4TPC::ArCF4Isobutane_CF4_frac);
0632     seeder->setNitrogenFraction(G4TPC::ArCF4Isobutane_N2_frac);
0633     seeder->setIsobutaneFraction(G4TPC::ArCF4Isobutane_isobutane_frac);
0634   }
0635   else
0636   {
0637   }
0638 
0639   seeder->set_pp_mode(TRACKING::pp_mode);
0640   se->registerSubsystem(seeder);
0641 
0642   // expand stubs in the TPC using simple kalman filter
0643   auto cprop = new PHSimpleKFProp("PHSimpleKFProp");
0644   cprop->set_field_dir(G4MAGNET::magfield_rescale);
0645   if (!ConstField)
0646   {
0647     cprop->set_field_dir(-1 * G4MAGNET::magfield_rescale);
0648   }
0649   cprop->useConstBField(false);
0650   cprop->useFixedClusterError(true);
0651   cprop->set_max_window(5.);
0652   cprop->Verbosity(verbosity);
0653 
0654   if (G4TPC::TPC_GAS_MIXTURE == "NeCF4")
0655   {
0656     cprop->setNeonFraction(G4TPC::NeCF4_Ne_frac);
0657     cprop->setArgonFraction(G4TPC::NeCF4_Ar_frac);
0658     cprop->setCF4Fraction(G4TPC::NeCF4_CF4_frac);
0659     cprop->setNitrogenFraction(G4TPC::NeCF4_N2_frac);
0660     cprop->setIsobutaneFraction(G4TPC::NeCF4_isobutane_frac);
0661   }
0662   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4")
0663   {
0664     cprop->setNeonFraction(G4TPC::ArCF4_Ne_frac);
0665     cprop->setArgonFraction(G4TPC::ArCF4_Ar_frac);
0666     cprop->setCF4Fraction(G4TPC::ArCF4_CF4_frac);
0667     cprop->setNitrogenFraction(G4TPC::ArCF4_N2_frac);
0668     cprop->setIsobutaneFraction(G4TPC::ArCF4_isobutane_frac);
0669   }
0670   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4N2")
0671   {
0672     cprop->setNeonFraction(G4TPC::ArCF4N2_Ne_frac);
0673     cprop->setArgonFraction(G4TPC::ArCF4N2_Ar_frac);
0674     cprop->setCF4Fraction(G4TPC::ArCF4N2_CF4_frac);
0675     cprop->setNitrogenFraction(G4TPC::ArCF4N2_N2_frac);
0676     cprop->setIsobutaneFraction(G4TPC::ArCF4N2_isobutane_frac);
0677   }
0678   else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4Isobutane")
0679   {
0680     cprop->setNeonFraction(G4TPC::ArCF4Isobutane_Ne_frac);
0681     cprop->setArgonFraction(G4TPC::ArCF4Isobutane_Ar_frac);
0682     cprop->setCF4Fraction(G4TPC::ArCF4Isobutane_CF4_frac);
0683     cprop->setNitrogenFraction(G4TPC::ArCF4Isobutane_N2_frac);
0684     cprop->setIsobutaneFraction(G4TPC::ArCF4Isobutane_isobutane_frac);
0685   }
0686   else
0687   {
0688   }
0689 
0690   cprop->set_pp_mode(TRACKING::pp_mode);
0691   se->registerSubsystem(cprop);
0692 
0693   if (TRACKING::pp_mode)
0694   {
0695     // for pp mode, apply preliminary distortion corrections to TPC clusters before crossing is known
0696     // and refit the trackseeds. Replace KFProp fits with the new fit parameters in the TPC seeds.
0697     auto prelim_distcorr = new PrelimDistortionCorrection;
0698     prelim_distcorr->set_pp_mode(TRACKING::pp_mode);
0699     prelim_distcorr->Verbosity(verbosity);
0700 
0701     if (G4TPC::TPC_GAS_MIXTURE == "NeCF4")
0702     {
0703       prelim_distcorr->setNeonFraction(G4TPC::NeCF4_Ne_frac);
0704       prelim_distcorr->setArgonFraction(G4TPC::NeCF4_Ar_frac);
0705       prelim_distcorr->setCF4Fraction(G4TPC::NeCF4_CF4_frac);
0706       prelim_distcorr->setNitrogenFraction(G4TPC::NeCF4_N2_frac);
0707       prelim_distcorr->setIsobutaneFraction(G4TPC::NeCF4_isobutane_frac);
0708     }
0709     else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4")
0710     {
0711       prelim_distcorr->setNeonFraction(G4TPC::ArCF4_Ne_frac);
0712       prelim_distcorr->setArgonFraction(G4TPC::ArCF4_Ar_frac);
0713       prelim_distcorr->setCF4Fraction(G4TPC::ArCF4_CF4_frac);
0714       prelim_distcorr->setNitrogenFraction(G4TPC::ArCF4_N2_frac);
0715       prelim_distcorr->setIsobutaneFraction(G4TPC::ArCF4_isobutane_frac);
0716     }
0717     else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4N2")
0718     {
0719       prelim_distcorr->setNeonFraction(G4TPC::ArCF4N2_Ne_frac);
0720       prelim_distcorr->setArgonFraction(G4TPC::ArCF4N2_Ar_frac);
0721       prelim_distcorr->setCF4Fraction(G4TPC::ArCF4N2_CF4_frac);
0722       prelim_distcorr->setNitrogenFraction(G4TPC::ArCF4N2_N2_frac);
0723       prelim_distcorr->setIsobutaneFraction(G4TPC::ArCF4N2_isobutane_frac);
0724     }
0725     else if (G4TPC::TPC_GAS_MIXTURE == "ArCF4Isobutane")
0726     {
0727       prelim_distcorr->setNeonFraction(G4TPC::ArCF4Isobutane_Ne_frac);
0728       prelim_distcorr->setArgonFraction(G4TPC::ArCF4Isobutane_Ar_frac);
0729       prelim_distcorr->setCF4Fraction(G4TPC::ArCF4Isobutane_CF4_frac);
0730       prelim_distcorr->setNitrogenFraction(G4TPC::ArCF4Isobutane_N2_frac);
0731       prelim_distcorr->setIsobutaneFraction(G4TPC::ArCF4Isobutane_isobutane_frac);
0732     }
0733     else
0734     {
0735     }
0736 
0737     se->registerSubsystem(prelim_distcorr);
0738   }
0739 
0740   // match silicon track seeds to TPC track seeds
0741 
0742   // The normal silicon association methods
0743   // Match the TPC track stubs from the CA seeder to silicon track stubs from PHSiliconTruthTrackSeeding
0744   auto silicon_match = new PHSiliconTpcTrackMatching;
0745   silicon_match->Verbosity(verbosity);
0746   silicon_match->set_pp_mode(TRACKING::pp_mode);
0747 
0748   silicon_match->set_phi_search_window(0.2);
0749   silicon_match->set_eta_search_window(0.015);
0750   silicon_match->set_x_search_window(std::numeric_limits<double>::max());
0751   silicon_match->set_y_search_window(std::numeric_limits<double>::max());
0752   silicon_match->set_z_search_window(std::numeric_limits<double>::max());
0753 
0754   silicon_match->set_test_windows_printout(false);  // used for tuning search windows
0755   se->registerSubsystem(silicon_match);
0756 
0757   // Associate Micromegas clusters with the tracks
0758   if (Enable::MICROMEGAS)
0759   {
0760     // Match TPC track stubs from CA seeder to clusters in the micromegas layers
0761     auto mm_match = new PHMicromegasTpcTrackMatching;
0762     mm_match->Verbosity(verbosity);
0763 
0764     mm_match->set_rphi_search_window_lyr1(0.4);
0765     mm_match->set_rphi_search_window_lyr2(13.0);
0766     mm_match->set_z_search_window_lyr1(26.0);
0767     mm_match->set_z_search_window_lyr2(0.2);
0768 
0769     mm_match->set_min_tpc_layer(38);             // layer in TPC to start projection fit
0770     mm_match->set_test_windows_printout(false);  // used for tuning search windows only
0771     se->registerSubsystem(mm_match);
0772   }
0773 }
0774 
0775 void alignment(std::string datafilename = "mille_output_data_file",
0776                std::string steeringfilename = "mille_steer")
0777 {
0778   Fun4AllServer* se = Fun4AllServer::instance();
0779   int verbosity = std::max(Enable::VERBOSITY, Enable::TRACKING_VERBOSITY);
0780 
0781   auto mille = new MakeMilleFiles;
0782   mille->Verbosity(verbosity);
0783   mille->set_datafile_name(datafilename + ".bin");
0784   mille->set_steeringfile_name(steeringfilename + ".txt");
0785   se->registerSubsystem(mille);
0786 
0787   auto helical = new HelicalFitter;
0788   helical->Verbosity(verbosity);
0789   helical->set_datafile_name(datafilename + "_helical.bin");
0790   helical->set_steeringfile_name(steeringfilename + "_helical.txt");
0791   se->registerSubsystem(helical);
0792 }
0793 
0794 void Tracking_Reco()
0795 {
0796   /*
0797    * just a wrapper around track seeding and track fitting methods,
0798    * to minimize disruption to existing steering macros
0799    */
0800   if (G4TRACKING::use_alignment)
0801   {
0802     Tracking_Reco_CommissioningTrackSeed();
0803   }
0804   else
0805   {
0806     Tracking_Reco_TrackSeed();
0807   }
0808 
0809   if (G4TRACKING::convert_seeds_to_svtxtracks)
0810   {
0811     convert_seeds();
0812     vertexing();
0813   }
0814   else
0815   {
0816     Tracking_Reco_TrackFit();
0817   }
0818 
0819   if (G4TRACKING::iterative_seeding)
0820   {
0821     Tracking_Reco_TrackSeed_pass1();
0822 
0823     if (G4TRACKING::convert_seeds_to_svtxtracks)
0824     {
0825       convert_seeds();
0826       vertexing();
0827     }
0828   }
0829 
0830   if (G4TRACKING::use_alignment)
0831   {
0832     alignment();
0833   }
0834 }
0835 
0836 void Filter_Conversion_Electrons(std::string ntuple_outfile)
0837 {
0838   Fun4AllServer* se = Fun4AllServer::instance();
0839   int verbosity = std::max(Enable::VERBOSITY, Enable::TRACKING_VERBOSITY);
0840 
0841   SecondaryVertexFinder* secvert = new SecondaryVertexFinder;
0842   secvert->Verbosity(verbosity);
0843   //  secvert->set_write_electrons_node(true);  // writes copy of filtered electron tracks to node tree
0844   //  secvert->set_write_ntuple(false);  // writes ntuple for tuning cuts
0845   secvert->setDecayParticleMass(0.000511);  // for electrons
0846   secvert->setOutfileName(ntuple_outfile);
0847   se->registerSubsystem(secvert);
0848 }
0849 
0850 void Reject_Laser_Events()
0851 {
0852   if (G4TPC::REJECT_LASER_EVENTS)
0853   {
0854     auto se = Fun4AllServer::instance();
0855     int verbosity = std::max(Enable::VERBOSITY, Enable::TRACKING_VERBOSITY);
0856 
0857     LaserEventRejecter *rejecter = new LaserEventRejecter();
0858     se->registerSubsystem(rejecter);
0859   }
0860 }
0861 
0862 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!