packages/trackreco/ActsEvaluator.cc

0001 #include "ActsEvaluator.h"
0002
0003 /// General fun4all and subsysreco includes
0004 #include <fun4all/Fun4AllReturnCodes.h>
0005 #include <fun4all/Fun4AllServer.h>
0006 #include <g4main/PHG4Hit.h>
0007 #include <g4main/PHG4Particle.h>
0008 #include <g4main/PHG4TruthInfoContainer.h>
0009 #include <phool/PHCompositeNode.h>
0010 #include <phool/getClass.h>
0011
0012 #include <trackbase/TrkrCluster.h>
0013 #include <trackbase/TrkrClusterContainer.h>
0014
0015 #include <trackbase_historic/SvtxTrack.h>
0016 #include <trackbase_historic/SvtxTrackMap.h>
0017 #include <trackbase_historic/TrackSeed.h>
0018 #include <trackbase_historic/TrackSeedContainer.h>
0019 #include <trackbase_historic/TrackSeedHelper.h>
0020
0021 #include <trackbase/InttDefs.h>
0022 #include <trackbase/MvtxDefs.h>
0023 #include <trackbase/TpcDefs.h>
0024
0025 #include <g4eval/SvtxClusterEval.h>
0026 #include <g4eval/SvtxEvalStack.h>
0027 #include <g4eval/SvtxTrackEval.h>
0028
0029 #include <globalvertex/SvtxVertex.h>
0030 #include <globalvertex/SvtxVertexMap.h>
0031
0032 #include <g4main/PHG4VtxPoint.h>
0033
0034 #include <Acts/Definitions/Algebra.hpp>
0035 #include <Acts/Definitions/Units.hpp>
0036 #include <Acts/EventData/MultiTrajectoryHelpers.hpp>
0037
0038 #include <TFile.h>
0039 #include <TTree.h>
0040
0041 #include <cmath>
0042
0043
0044 ActsEvaluator::ActsEvaluator(const std::string& name)
0045   : m_filename(name)
0046 {
0047 }
0048
0049 void ActsEvaluator::Init(PHCompositeNode* topNode)
0050 {
0051   if (m_verbosity > 1)
0052   {
0053     std::cout << "Starting ActsEvaluator::Init" << std::endl;
0054   }
0055
0056   if (getNodes(topNode) == Fun4AllReturnCodes::ABORTEVENT)
0057   {
0058     std::cout << "Error: Nodes not available in ActsEvaluator::Init"
0059               << std::endl;
0060     return;
0061   }
0062
0063   initializeTree();
0064
0065   if (m_verbosity > 1)
0066   {
0067     std::cout << "Finished ActsEvaluator::Init" << std::endl;
0068   }
0069 }
0070 void ActsEvaluator::next_event(PHCompositeNode* topNode)
0071 {
0072   m_eventNr++;
0073   if(m_isData)
0074   {
0075     return;
0076   }
0077   if (!m_svtxEvalStack)
0078   {
0079     m_svtxEvalStack = new SvtxEvalStack(topNode);
0080   }
0081
0082   m_svtxEvalStack->next_event(topNode);
0083
0084 }
0085 void ActsEvaluator::process_track(const ActsTrackFittingAlgorithm::TrackContainer& tracks,
0086                   std::vector<Acts::MultiTrajectoryTraits::IndexType>& trackTips,
0087                   Trajectory::IndexedParameters& paramsMap,
0088                                   SvtxTrack* track,
0089                                   const TrackSeed* seed,
0090                                   const ActsTrackFittingAlgorithm::MeasurementContainer& measurements)
0091 {
0092   if (m_verbosity > 1)
0093   {
0094     std::cout << "Starting ActsEvaluator at event " << m_eventNr
0095               << std::endl;
0096   }
0097
0098   evaluateTrackFit(tracks,trackTips,paramsMap, track, seed, measurements);
0099
0100
0101   if (m_verbosity > 1)
0102   {
0103     std::cout << "Finished ActsEvaluator::process_event" << std::endl;
0104   }
0105 }
0106
0107 void ActsEvaluator::evaluateTrackFit(const ActsTrackFittingAlgorithm::TrackContainer& tracks,
0108                      std::vector<Acts::MultiTrajectoryTraits::IndexType>& trackTips,
0109                      Trajectory::IndexedParameters& paramsMap,
0110                                      SvtxTrack* track,
0111                                      const TrackSeed* seed,
0112                                      const ActsTrackFittingAlgorithm::MeasurementContainer& measurements)
0113 {
0114   if (m_verbosity > 5)
0115     {
0116       std::cout << "Evaluating Acts track fits" << std::endl;
0117     }
0118   /// Skip failed fits
0119   if (trackTips.empty())
0120   {
0121     if (m_verbosity > 1)
0122     {
0123       std::cout << "TrackTips empty in ActsEvaluator" << std::endl;
0124     }
0125     return;
0126   }
0127   SvtxTrackEval* trackeval = nullptr;
0128   if (m_svtxEvalStack)
0129   {
0130     trackeval = m_svtxEvalStack->get_track_eval();
0131   }
0132   int iTrack = track->get_id();
0133   int iTraj = iTrack;
0134   if (m_verbosity > 2)
0135   {
0136     std::cout << "Starting trajectory with trackKey " << track->get_id()
0137               << " and corresponding to tpc track seed "
0138               << seed->get_tpc_seed_index() << std::endl;
0139   }
0140
0141   const auto& mj = tracks.trackStateContainer();
0142   const auto& trackTip = trackTips.front();
0143   m_trajNr = iTraj;
0144
0145   if (m_verbosity > 2)
0146   {
0147     std::cout << "beginning trackTip " << trackTip
0148               << " corresponding to track " << iTrack
0149               << " in trajectroy " << iTraj << std::endl;
0150   }
0151
0152   if (m_verbosity > 2)
0153   {
0154     std::cout << "Evaluating track key " << iTrack
0155               << " for track tip " << trackTip << std::endl;
0156   }
0157   PHG4Particle* g4particle = nullptr;
0158   if (trackeval)
0159   {
0160     g4particle = trackeval->max_truth_particle_by_nclusters(track);
0161   }
0162   if (m_verbosity > 1)
0163   {
0164     std::cout << "Analyzing SvtxTrack " << iTrack << std::endl;
0165     if(g4particle){
0166     std::cout << "TruthParticle : " << g4particle->get_px()
0167               << ", " << g4particle->get_py() << ", "
0168               << g4particle->get_pz() << ", " << g4particle->get_e()
0169               << std::endl;
0170     }
0171   }
0172
0173   m_trackNr = iTrack;
0174
0175   auto trajState =
0176       Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip);
0177
0178   const auto& params = paramsMap.find(trackTip)->second;
0179
0180   if (m_verbosity > 1)
0181   {
0182
0183       std::cout << "Fitted params : "
0184                 << params.position(m_tGeometry->geometry().getGeoContext())
0185                 << std::endl
0186                 << params.momentum()
0187                 << std::endl;
0188       std::cout << "Track has " << trajState.nMeasurements
0189                 << " measurements and " << trajState.nHoles
0190                 << " holes and " << trajState.nOutliers
0191                 << " outliers and " << trajState.nStates
0192                 << " states " << std::endl;
0193
0194   }
0195
0196   m_nMeasurements = trajState.nMeasurements;
0197   m_nStates = trajState.nStates;
0198   m_nOutliers = trajState.nOutliers;
0199   m_nSharedHits = trajState.nSharedHits;
0200   m_nHoles = trajState.nHoles;
0201   m_chi2_fit = trajState.chi2Sum;
0202   m_ndf_fit = trajState.NDF;
0203   m_quality = track->get_quality();
0204
0205 if(g4particle){
0206   fillG4Particle(g4particle);
0207   }
0208   fillProtoTrack(seed);
0209   fillFittedTrackParams(paramsMap, trackTip);
0210   visitTrackStates(mj, trackTip, measurements);
0211
0212   m_trackTree->Fill();
0213
0214   /// Start fresh for the next track
0215   clearTrackVariables();
0216   if (m_verbosity > 1)
0217   {
0218     std::cout << "Finished track " << iTrack << std::endl;
0219   }
0220
0221   if (m_verbosity > 1)
0222   {
0223     std::cout << "Analyzed " << iTrack << " tracks in trajectory number "
0224               << iTraj << std::endl;
0225   }
0226
0227   if (m_verbosity > 5)
0228   {
0229     std::cout << "Finished evaluating track fits" << std::endl;
0230   }
0231   return;
0232 }
0233
0234 void ActsEvaluator::End()
0235 {
0236   m_trackFile->cd();
0237   m_trackTree->Write();
0238   m_trackFile->Close();
0239 }
0240
0241 void ActsEvaluator::visitTrackStates(const Acts::ConstVectorMultiTrajectory& traj,
0242                                      const size_t& trackTip,
0243                                      const ActsTrackFittingAlgorithm::MeasurementContainer& measurements)
0244 {
0245   if (m_verbosity > 2)
0246   {
0247     std::cout << "Begin visit track states" << std::endl;
0248   }
0249
0250   traj.visitBackwards(trackTip, [&](const auto& state)
0251                       {
0252     /// Only fill the track states with non-outlier measurement
0253     auto typeFlags = state.typeFlags();
0254     if (not typeFlags.test(Acts::TrackStateFlag::MeasurementFlag))
0255     {
0256       return true;
0257     }
0258
0259     const auto& surface = state.referenceSurface();
0260
0261     /// Get the geometry ID
0262     auto geoID = surface.geometryId();
0263     m_volumeID.push_back(geoID.volume());
0264     m_layerID.push_back(geoID.layer());
0265     m_moduleID.push_back(geoID.sensitive());
0266     if(m_verbosity > 3)
0267       {
0268       std::cout << "Cluster volume : layer : sensitive " << geoID.volume()
0269         << " : " << geoID.layer() << " : "
0270         << geoID.sensitive() << std::endl;
0271       }
0272     auto sourceLink = state.getUncalibratedSourceLink().template get<ActsSourceLink>();
0273     const auto& cluskey = sourceLink.cluskey();
0274
0275     m_sphenixlayer.push_back(TrkrDefs::getLayer(cluskey));
0276     Acts::Vector2 local = Acts::Vector2::Zero();
0277
0278     /// get the local measurement that acts used
0279     std::visit([&](const auto& meas) {
0280     local(0) = meas.parameters()[0];
0281     local(1) = meas.parameters()[1];
0282       }, measurements[sourceLink.index()]);
0283
0284     /// Get global position
0285     /// This is an arbitrary vector. Doesn't matter in coordinate transformation
0286     /// in Acts code
0287     Acts::Vector3 mom(1, 1, 1);
0288     Acts::Vector3 global = surface.localToGlobal(m_tGeometry->geometry().getGeoContext(),
0289                          local, mom);
0290
0291     /// Get measurement covariance
0292     /// These are the same
0293     //auto slcov = sourceLink.covariance();
0294     auto cov = state.effectiveCalibratedCovariance();
0295
0296     m_lx_hit.push_back(local.x());
0297     m_ly_hit.push_back(local.y());
0298     m_x_hit.push_back(global.x());
0299     m_y_hit.push_back(global.y());
0300     m_z_hit.push_back(global.z());
0301
0302     /// Get the truth hit corresponding to this trackState
0303     /// We go backwards from hitID -> TrkrDefs::cluskey to g4hit with
0304     /// the map created in PHActsSourceLinks
0305     float gt = -9999;
0306     Acts::Vector3 globalTruthPos = getGlobalTruthHit(cluskey, gt);
0307     float truthLOC0 = std::numeric_limits<float>::quiet_NaN();
0308     float truthLOC1 = std::numeric_limits<float>::quiet_NaN();
0309     float truthPHI = std::numeric_limits<float>::quiet_NaN();
0310     float truthTHETA = std::numeric_limits<float>::quiet_NaN();
0311     float truthQOP = std::numeric_limits<float>::quiet_NaN();
0312     float truthTIME = std::numeric_limits<float>::quiet_NaN();
0313     float momentum = std::numeric_limits<float>::quiet_NaN();
0314     if(!std::isnan(globalTruthPos(0)))
0315     {
0316     float gx = globalTruthPos(0);
0317     float gy = globalTruthPos(1);
0318     float gz = globalTruthPos(2);
0319
0320     /// Get local truth position
0321     const float r = std::sqrt(gx * gx + gy * gy + gz * gz);
0322     Acts::Vector3 globalTruthUnitDir(gx / r, gy / r, gz / r);
0323
0324     auto vecResult = surface.globalToLocal(
0325         m_tGeometry->geometry().getGeoContext(),
0326         globalTruthPos,
0327         globalTruthUnitDir);
0328
0329     /// Push the truth hit info
0330     m_t_x.push_back(gx);
0331     m_t_y.push_back(gy);
0332     m_t_z.push_back(gz);
0333     m_t_r.push_back(std::sqrt(gx * gx + gy * gy));
0334     m_t_dx.push_back(gx / r);
0335     m_t_dy.push_back(gy / r);
0336     m_t_dz.push_back(gz / r);
0337
0338     /// Get the truth track parameter at this track State
0339
0340     if(!std::isnan(m_t_px))
0341     {
0342         momentum = std::sqrt(m_t_px * m_t_px +
0343            m_t_py * m_t_py +
0344            m_t_pz * m_t_pz);
0345     }
0346     if(vecResult.ok())
0347       {
0348     Acts::Vector2 truthLocVec = vecResult.value();
0349     truthLOC0 = truthLocVec.x();
0350     truthLOC1 = truthLocVec.y();
0351       }
0352     else
0353       {
0354     truthLOC0 = -9999.;
0355     truthLOC1 = -9999.;
0356       }
0357
0358     truthPHI = phi(globalTruthUnitDir);
0359     truthTHETA = theta(globalTruthUnitDir);
0360     truthQOP =
0361         m_t_charge / momentum;
0362     truthTIME = gt;
0363
0364     m_t_eLOC0.push_back(truthLOC0);
0365     m_t_eLOC1.push_back(truthLOC1);
0366     m_t_ePHI.push_back(truthPHI);
0367     m_t_eTHETA.push_back(truthTHETA);
0368     m_t_eQOP.push_back(truthQOP);
0369     m_t_eT.push_back(truthTIME);
0370     }
0371     else
0372     {
0373       m_t_x.push_back(std::numeric_limits<float>::quiet_NaN());
0374       m_t_y.push_back(std::numeric_limits<float>::quiet_NaN());
0375       m_t_z.push_back(std::numeric_limits<float>::quiet_NaN());
0376       m_t_r.push_back(std::numeric_limits<float>::quiet_NaN());
0377       m_t_dx.push_back(std::numeric_limits<float>::quiet_NaN());
0378       m_t_dy.push_back(std::numeric_limits<float>::quiet_NaN());
0379       m_t_dz.push_back(std::numeric_limits<float>::quiet_NaN());
0380
0381       m_t_eLOC0.push_back(std::numeric_limits<float>::quiet_NaN());
0382       m_t_eLOC1.push_back(std::numeric_limits<float>::quiet_NaN());
0383       m_t_ePHI.push_back(std::numeric_limits<float>::quiet_NaN());
0384       m_t_eTHETA.push_back(std::numeric_limits<float>::quiet_NaN());
0385       m_t_eQOP.push_back(std::numeric_limits<float>::quiet_NaN());
0386       m_t_eT.push_back(std::numeric_limits<float>::quiet_NaN());
0387     }
0388     /// Get the predicted parameter for this state
0389     bool predicted = false;
0390     if (state.hasPredicted())
0391     {
0392       predicted = true;
0393       m_nPredicted++;
0394
0395       auto parameters = state.predicted();
0396       auto covariance = state.predictedCovariance();
0397
0398       /// Local hit residual info
0399       auto H = state.effectiveProjector();
0400       auto resCov = cov + H * covariance * H.transpose();
0401       auto residual = state.effectiveCalibrated() - H * parameters;
0402       m_res_x_hit.push_back(residual(Acts::eBoundLoc0));
0403       m_res_y_hit.push_back(residual(Acts::eBoundLoc1));
0404       m_err_x_hit.push_back(
0405           std::sqrt(resCov(Acts::eBoundLoc0, Acts::eBoundLoc0)));
0406       m_err_y_hit.push_back(
0407           std::sqrt(resCov(Acts::eBoundLoc1, Acts::eBoundLoc1)));
0408       m_pull_x_hit.push_back(
0409           residual(Acts::eBoundLoc0) /
0410           std::sqrt(resCov(Acts::eBoundLoc0, Acts::eBoundLoc0)));
0411       m_pull_y_hit.push_back(
0412           residual(Acts::eBoundLoc1) /
0413           std::sqrt(resCov(Acts::eBoundLoc1, Acts::eBoundLoc1)));
0414       m_dim_hit.push_back(state.calibratedSize());
0415
0416       /// Predicted parameter
0417       m_eLOC0_prt.push_back(parameters[Acts::eBoundLoc0]);
0418       m_eLOC1_prt.push_back(parameters[Acts::eBoundLoc1]);
0419       m_ePHI_prt.push_back(parameters[Acts::eBoundPhi]);
0420       m_eTHETA_prt.push_back(parameters[Acts::eBoundTheta]);
0421       m_eQOP_prt.push_back(parameters[Acts::eBoundQOverP]);
0422       m_eT_prt.push_back(parameters[Acts::eBoundTime]);
0423
0424       /// Predicted residual
0425       m_res_eLOC0_prt.push_back(parameters[Acts::eBoundLoc0] -
0426                                 truthLOC0);
0427       m_res_eLOC1_prt.push_back(parameters[Acts::eBoundLoc1] -
0428                                 truthLOC1);
0429       m_res_ePHI_prt.push_back(parameters[Acts::eBoundPhi] -
0430                                truthPHI);
0431       m_res_eTHETA_prt.push_back(
0432           parameters[Acts::eBoundTheta] - truthTHETA);
0433       m_res_eQOP_prt.push_back(parameters[Acts::eBoundQOverP] -
0434                                truthQOP);
0435       m_res_eT_prt.push_back(parameters[Acts::eBoundTime] -
0436                              truthTIME);
0437
0438       /// Predicted parameter Uncertainties
0439       m_err_eLOC0_prt.push_back(
0440           std::sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
0441       m_err_eLOC1_prt.push_back(
0442           std::sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
0443       m_err_ePHI_prt.push_back(
0444           std::sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
0445       m_err_eTHETA_prt.push_back(
0446           std::sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
0447       m_err_eQOP_prt.push_back(
0448           std::sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
0449       m_err_eT_prt.push_back(
0450           std::sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
0451
0452       /// Predicted parameter pulls
0453       m_pull_eLOC0_prt.push_back(
0454           (parameters[Acts::eBoundLoc0] - truthLOC0) /
0455           std::sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
0456       m_pull_eLOC1_prt.push_back(
0457           (parameters[Acts::eBoundLoc1] - truthLOC1) /
0458           std::sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
0459       m_pull_ePHI_prt.push_back(
0460           (parameters[Acts::eBoundPhi] - truthPHI) /
0461           std::sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
0462       m_pull_eTHETA_prt.push_back(
0463           (parameters[Acts::eBoundTheta] - truthTHETA) /
0464           std::sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
0465       m_pull_eQOP_prt.push_back(
0466           (parameters[Acts::eBoundQOverP] - truthQOP) /
0467           std::sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
0468       m_pull_eT_prt.push_back(
0469           (parameters[Acts::eBoundTime] - truthTIME) /
0470           std::sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
0471
0472       Acts::FreeVector freeParams =
0473     Acts::transformBoundToFreeParameters(state.referenceSurface(),
0474                              m_tGeometry->geometry().getGeoContext(),
0475                              parameters);
0476
0477       m_x_prt.push_back(freeParams[Acts::eFreePos0]);
0478       m_y_prt.push_back(freeParams[Acts::eFreePos1]);
0479       m_z_prt.push_back(freeParams[Acts::eFreePos2]);
0480       auto p = std::abs(1 / freeParams[Acts::eFreeQOverP]);
0481       m_px_prt.push_back(p * freeParams[Acts::eFreeDir0]);
0482       m_py_prt.push_back(p * freeParams[Acts::eFreeDir1]);
0483       m_pz_prt.push_back(p * freeParams[Acts::eFreeDir2]);
0484       m_pT_prt.push_back(p * std::hypot(freeParams[Acts::eFreeDir0],
0485                     freeParams[Acts::eFreeDir1]));
0486       m_eta_prt.push_back(
0487          Acts::VectorHelpers::eta(freeParams.segment<3>(Acts::eFreeDir0)));
0488     }
0489     else
0490     {
0491       /// Push bad values if no predicted parameter
0492       m_res_x_hit.push_back(-9999);
0493       m_res_y_hit.push_back(-9999);
0494       m_err_x_hit.push_back(-9999);
0495       m_err_y_hit.push_back(-9999);
0496       m_pull_x_hit.push_back(-9999);
0497       m_pull_y_hit.push_back(-9999);
0498       m_dim_hit.push_back(-9999);
0499       m_eLOC0_prt.push_back(-9999);
0500       m_eLOC1_prt.push_back(-9999);
0501       m_ePHI_prt.push_back(-9999);
0502       m_eTHETA_prt.push_back(-9999);
0503       m_eQOP_prt.push_back(-9999);
0504       m_eT_prt.push_back(-9999);
0505       m_res_eLOC0_prt.push_back(-9999);
0506       m_res_eLOC1_prt.push_back(-9999);
0507       m_res_ePHI_prt.push_back(-9999);
0508       m_res_eTHETA_prt.push_back(-9999);
0509       m_res_eQOP_prt.push_back(-9999);
0510       m_res_eT_prt.push_back(-9999);
0511       m_err_eLOC0_prt.push_back(-9999);
0512       m_err_eLOC1_prt.push_back(-9999);
0513       m_err_ePHI_prt.push_back(-9999);
0514       m_err_eTHETA_prt.push_back(-9999);
0515       m_err_eQOP_prt.push_back(-9999);
0516       m_err_eT_prt.push_back(-9999);
0517       m_pull_eLOC0_prt.push_back(-9999);
0518       m_pull_eLOC1_prt.push_back(-9999);
0519       m_pull_ePHI_prt.push_back(-9999);
0520       m_pull_eTHETA_prt.push_back(-9999);
0521       m_pull_eQOP_prt.push_back(-9999);
0522       m_pull_eT_prt.push_back(-9999);
0523       m_x_prt.push_back(-9999);
0524       m_y_prt.push_back(-9999);
0525       m_z_prt.push_back(-9999);
0526       m_px_prt.push_back(-9999);
0527       m_py_prt.push_back(-9999);
0528       m_pz_prt.push_back(-9999);
0529       m_pT_prt.push_back(-9999);
0530       m_eta_prt.push_back(-9999);
0531     }
0532
0533     bool filtered = false;
0534     if (state.hasFiltered())
0535     {
0536       filtered = true;
0537       m_nFiltered++;
0538
0539       auto parameter = state.filtered();
0540       auto covariance = state.filteredCovariance();
0541
0542       m_eLOC0_flt.push_back(parameter[Acts::eBoundLoc0]);
0543       m_eLOC1_flt.push_back(parameter[Acts::eBoundLoc1]);
0544       m_ePHI_flt.push_back(parameter[Acts::eBoundPhi]);
0545       m_eTHETA_flt.push_back(parameter[Acts::eBoundTheta]);
0546       m_eQOP_flt.push_back(parameter[Acts::eBoundQOverP]);
0547       m_eT_flt.push_back(parameter[Acts::eBoundTime]);
0548
0549       m_res_eLOC0_flt.push_back(parameter[Acts::eBoundLoc0] -
0550                                 truthLOC0);
0551       m_res_eLOC1_flt.push_back(parameter[Acts::eBoundLoc1] -
0552                                 truthLOC1);
0553       m_res_ePHI_flt.push_back(parameter[Acts::eBoundPhi] -
0554                                truthPHI);
0555       m_res_eTHETA_flt.push_back(parameter[Acts::eBoundTheta] -
0556                                  truthTHETA);
0557       m_res_eQOP_flt.push_back(parameter[Acts::eBoundQOverP] -
0558                                truthQOP);
0559       m_res_eT_flt.push_back(parameter[Acts::eBoundTime] -
0560                              truthTIME);
0561
0562       /// Filtered parameter uncertainties
0563       m_err_eLOC0_flt.push_back(
0564           std::sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
0565       m_err_eLOC1_flt.push_back(
0566           std::sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
0567       m_err_ePHI_flt.push_back(
0568           std::sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
0569       m_err_eTHETA_flt.push_back(
0570           std::sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
0571       m_err_eQOP_flt.push_back(
0572           std::sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
0573       m_err_eT_flt.push_back(
0574           std::sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
0575
0576       /// Filtered parameter pulls
0577       m_pull_eLOC0_flt.push_back(
0578           (parameter[Acts::eBoundLoc0] - truthLOC0) /
0579           std::sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
0580       m_pull_eLOC1_flt.push_back(
0581           (parameter[Acts::eBoundLoc1] - truthLOC1) /
0582           std::sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
0583       m_pull_ePHI_flt.push_back(
0584           (parameter[Acts::eBoundPhi] - truthPHI) /
0585           std::sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
0586       m_pull_eTHETA_flt.push_back(
0587           (parameter[Acts::eBoundTheta] - truthTHETA) /
0588           std::sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
0589       m_pull_eQOP_flt.push_back(
0590           (parameter[Acts::eBoundQOverP] - truthQOP) /
0591           std::sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
0592       m_pull_eT_flt.push_back(
0593           (parameter[Acts::eBoundTime] - truthTIME) /
0594           std::sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
0595
0596       Acts::FreeVector freeparams = Acts::transformBoundToFreeParameters(surface, m_tGeometry->geometry().getGeoContext(), parameter);
0597
0598       /// Other filtered parameter info
0599       m_x_flt.push_back(freeparams[Acts::eFreePos0]);
0600       m_y_flt.push_back(freeparams[Acts::eFreePos1]);
0601       m_z_flt.push_back(freeparams[Acts::eFreePos2]);
0602       auto p = std::abs(1/freeparams[Acts::eFreeQOverP]);
0603       m_px_flt.push_back(p * freeparams[Acts::eFreeDir0]);
0604       m_py_flt.push_back(p * freeparams[Acts::eFreeDir1]);
0605       m_pz_flt.push_back(p * freeparams[Acts::eFreeDir2]);
0606       m_pT_flt.push_back(std::sqrt(p * std::hypot(freeparams[Acts::eFreeDir0],
0607                          freeparams[Acts::eFreeDir1])));
0608       m_eta_flt.push_back(eta(freeparams.segment<3>(Acts::eFreeDir0)));
0609       m_chi2.push_back(state.chi2());
0610
0611     }
0612     else
0613     {
0614       /// Push bad values if no filtered parameter
0615       m_eLOC0_flt.push_back(-9999);
0616       m_eLOC1_flt.push_back(-9999);
0617       m_ePHI_flt.push_back(-9999);
0618       m_eTHETA_flt.push_back(-9999);
0619       m_eQOP_flt.push_back(-9999);
0620       m_eT_flt.push_back(-9999);
0621       m_res_eLOC0_flt.push_back(-9999);
0622       m_res_eLOC1_flt.push_back(-9999);
0623       m_res_ePHI_flt.push_back(-9999);
0624       m_res_eTHETA_flt.push_back(-9999);
0625       m_res_eQOP_flt.push_back(-9999);
0626       m_res_eT_flt.push_back(-9999);
0627       m_err_eLOC0_flt.push_back(-9999);
0628       m_err_eLOC1_flt.push_back(-9999);
0629       m_err_ePHI_flt.push_back(-9999);
0630       m_err_eTHETA_flt.push_back(-9999);
0631       m_err_eQOP_flt.push_back(-9999);
0632       m_err_eT_flt.push_back(-9999);
0633       m_pull_eLOC0_flt.push_back(-9999);
0634       m_pull_eLOC1_flt.push_back(-9999);
0635       m_pull_ePHI_flt.push_back(-9999);
0636       m_pull_eTHETA_flt.push_back(-9999);
0637       m_pull_eQOP_flt.push_back(-9999);
0638       m_pull_eT_flt.push_back(-9999);
0639       m_x_flt.push_back(-9999);
0640       m_y_flt.push_back(-9999);
0641       m_z_flt.push_back(-9999);
0642       m_py_flt.push_back(-9999);
0643       m_pz_flt.push_back(-9999);
0644       m_pT_flt.push_back(-9999);
0645       m_eta_flt.push_back(-9999);
0646       m_chi2.push_back(-9999);
0647     }
0648
0649     bool smoothed = false;
0650     if (state.hasSmoothed())
0651     {
0652       smoothed = true;
0653       m_nSmoothed++;
0654
0655       auto parameter = state.smoothed();
0656       auto covariance = state.smoothedCovariance();
0657
0658       m_eLOC0_smt.push_back(parameter[Acts::eBoundLoc0]);
0659       m_eLOC1_smt.push_back(parameter[Acts::eBoundLoc1]);
0660       m_ePHI_smt.push_back(parameter[Acts::eBoundPhi]);
0661       m_eTHETA_smt.push_back(parameter[Acts::eBoundTheta]);
0662       m_eQOP_smt.push_back(parameter[Acts::eBoundQOverP]);
0663       m_eT_smt.push_back(parameter[Acts::eBoundTime]);
0664
0665       m_res_eLOC0_smt.push_back(parameter[Acts::eBoundLoc0] -
0666                                 truthLOC0);
0667       m_res_eLOC1_smt.push_back(parameter[Acts::eBoundLoc1] -
0668                                 truthLOC1);
0669       m_res_ePHI_smt.push_back(parameter[Acts::eBoundPhi] -
0670                                truthPHI);
0671       m_res_eTHETA_smt.push_back(parameter[Acts::eBoundTheta] -
0672                                  truthTHETA);
0673       m_res_eQOP_smt.push_back(parameter[Acts::eBoundQOverP] -
0674                                truthQOP);
0675       m_res_eT_smt.push_back(parameter[Acts::eBoundTime] -
0676                              truthTIME);
0677
0678       /// Smoothed parameter uncertainties
0679       m_err_eLOC0_smt.push_back(
0680           std::sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
0681       m_err_eLOC1_smt.push_back(
0682           std::sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
0683       m_err_ePHI_smt.push_back(
0684           std::sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
0685       m_err_eTHETA_smt.push_back(
0686           std::sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
0687       m_err_eQOP_smt.push_back(
0688           std::sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
0689       m_err_eT_smt.push_back(
0690           std::sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
0691
0692       /// Smoothed parameter pulls
0693       m_pull_eLOC0_smt.push_back(
0694           (parameter[Acts::eBoundLoc0] - truthLOC0) /
0695           std::sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
0696       m_pull_eLOC1_smt.push_back(
0697           (parameter[Acts::eBoundLoc1] - truthLOC1) /
0698           std::sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
0699       m_pull_ePHI_smt.push_back(
0700           (parameter[Acts::eBoundPhi] - truthPHI) /
0701           std::sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
0702       m_pull_eTHETA_smt.push_back(
0703           (parameter[Acts::eBoundTheta] - truthTHETA) /
0704           std::sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
0705       m_pull_eQOP_smt.push_back(
0706           (parameter[Acts::eBoundQOverP] - truthQOP) /
0707           std::sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
0708       m_pull_eT_smt.push_back(
0709           (parameter[Acts::eBoundTime] - truthTIME) /
0710           std::sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
0711
0712       Acts::FreeVector freeparams = Acts::transformBoundToFreeParameters(surface, m_tGeometry->geometry().getGeoContext(), parameter);
0713
0714       /// Other smoothed parameter info
0715       m_x_smt.push_back(freeparams[Acts::eFreePos0]);
0716       m_y_smt.push_back(freeparams[Acts::eFreePos1]);
0717       m_z_smt.push_back(freeparams[Acts::eFreePos2]);
0718       auto p = std::abs(1/freeparams[Acts::eFreeQOverP]);
0719       m_px_smt.push_back(p * freeparams[Acts::eFreeDir0]);
0720       m_py_smt.push_back(p * freeparams[Acts::eFreeDir1]);
0721       m_pz_smt.push_back(p * freeparams[Acts::eFreeDir2]);
0722       m_pT_smt.push_back(std::sqrt(p * std::hypot(freeparams[Acts::eFreeDir0],
0723                          freeparams[Acts::eFreeDir1])));
0724       m_eta_smt.push_back(eta(freeparams.segment<3>(Acts::eFreeDir0)));
0725
0726     }
0727     else
0728     {
0729       /// Push bad values if no smoothed parameter
0730       m_eLOC0_smt.push_back(-9999);
0731       m_eLOC1_smt.push_back(-9999);
0732       m_ePHI_smt.push_back(-9999);
0733       m_eTHETA_smt.push_back(-9999);
0734       m_eQOP_smt.push_back(-9999);
0735       m_eT_smt.push_back(-9999);
0736       m_res_eLOC0_smt.push_back(-9999);
0737       m_res_eLOC1_smt.push_back(-9999);
0738       m_res_ePHI_smt.push_back(-9999);
0739       m_res_eTHETA_smt.push_back(-9999);
0740       m_res_eQOP_smt.push_back(-9999);
0741       m_res_eT_smt.push_back(-9999);
0742       m_err_eLOC0_smt.push_back(-9999);
0743       m_err_eLOC1_smt.push_back(-9999);
0744       m_err_ePHI_smt.push_back(-9999);
0745       m_err_eTHETA_smt.push_back(-9999);
0746       m_err_eQOP_smt.push_back(-9999);
0747       m_err_eT_smt.push_back(-9999);
0748       m_pull_eLOC0_smt.push_back(-9999);
0749       m_pull_eLOC1_smt.push_back(-9999);
0750       m_pull_ePHI_smt.push_back(-9999);
0751       m_pull_eTHETA_smt.push_back(-9999);
0752       m_pull_eQOP_smt.push_back(-9999);
0753       m_pull_eT_smt.push_back(-9999);
0754       m_x_smt.push_back(-9999);
0755       m_y_smt.push_back(-9999);
0756       m_z_smt.push_back(-9999);
0757       m_px_smt.push_back(-9999);
0758       m_py_smt.push_back(-9999);
0759       m_pz_smt.push_back(-9999);
0760       m_pT_smt.push_back(-9999);
0761       m_eta_smt.push_back(-9999);
0762     }
0763
0764     /// Save whether or not states had various KF steps
0765     m_prt.push_back(predicted);
0766     m_flt.push_back(filtered);
0767     m_smt.push_back(smoothed);
0768
0769     return true; }  /// Finish lambda function
0770   );                       /// Finish multi trajectory visitBackwards call
0771
0772   if (m_verbosity > 2)
0773     std::cout << "Finished track states" << std::endl;
0774
0775   return;
0776 }
0777
0778 Acts::Vector3 ActsEvaluator::getGlobalTruthHit(TrkrDefs::cluskey cluskey,
0779                                                float& _gt)
0780 {
0781   Acts::Vector3 ret(std::numeric_limits<float>::quiet_NaN(),
0782                     std::numeric_limits<float>::quiet_NaN(),
0783                     std::numeric_limits<float>::quiet_NaN());
0784   if(m_svtxEvalStack == nullptr)
0785                     {
0786                       return ret;
0787                     }
0788   SvtxClusterEval* clustereval = m_svtxEvalStack->get_cluster_eval();
0789
0790   const auto [truth_ckey, truth_cluster] = clustereval->max_truth_cluster_by_energy(cluskey);
0791
0792   float gx = -9999;
0793   float gy = -9999;
0794   float gz = -9999;
0795   float gt = -9999;
0796
0797   if (truth_cluster)
0798   {
0799     gx = truth_cluster->getX();
0800     gy = truth_cluster->getY();
0801     gz = truth_cluster->getZ();
0802   }
0803
0804   /// Convert to acts units of mm
0805   gx *= Acts::UnitConstants::cm;
0806   gy *= Acts::UnitConstants::cm;
0807   gz *= Acts::UnitConstants::cm;
0808   ret(0) = gx;
0809   ret(1) = gy;
0810   ret(2) = gz;
0811   _gt = gt;
0812   return ret;
0813 }
0814
0815 //___________________________________________________________________________________
0816 Surface ActsEvaluator::getSurface(TrkrDefs::cluskey cluskey, TrkrCluster* cluster)
0817 {
0818   return m_tGeometry->maps().getSurface(cluskey, cluster);
0819 }
0820
0821 void ActsEvaluator::fillProtoTrack(const TrackSeed* seed)
0822 {
0823   if (m_verbosity > 2)
0824   {
0825     std::cout << "Filling proto track seed quantities" << std::endl;
0826   }
0827
0828   unsigned int tpcid = seed->get_tpc_seed_index();
0829   unsigned int siid = seed->get_silicon_seed_index();
0830
0831   auto siseed = m_siliconSeeds->get(siid);
0832   auto tpcseed = m_tpcSeeds->get(tpcid);
0833   if(!tpcseed) { return; }
0834
0835   const Acts::Vector3 position = TrackSeedHelper::get_xyz( siseed?siseed:tpcseed )*Acts::UnitConstants::cm;
0836   const Acts::Vector3 momentum(tpcseed->get_px(),tpcseed->get_py(),tpcseed->get_pz());
0837
0838   m_protoTrackPx = momentum(0);
0839   m_protoTrackPy = momentum(1);
0840   m_protoTrackPz = momentum(2);
0841   m_protoTrackX = position(0);
0842   m_protoTrackY = position(1);
0843   m_protoTrackZ = position(2);
0844
0845   for (auto svtxseed : {siseed, tpcseed})
0846   {
0847     //! protect against tpc only tracks
0848     if (!svtxseed)
0849     {
0850       continue;
0851     }
0852     for (auto clusIter = svtxseed->begin_cluster_keys();
0853          clusIter != svtxseed->end_cluster_keys();
0854          ++clusIter)
0855     {
0856       auto key = *clusIter;
0857       auto cluster = m_clusterContainer->findCluster(key);
0858
0859       /// Get source link global position
0860       Acts::Vector2 loc = m_tGeometry->getLocalCoords(key, cluster);
0861
0862       Acts::Vector3 mom(0, 0, 0);
0863       Acts::Vector3 globalPos = m_tGeometry->getGlobalPosition(key, cluster) * Acts::UnitConstants::cm;
0864
0865       m_SLx.push_back(globalPos(0));
0866       m_SLy.push_back(globalPos(1));
0867       m_SLz.push_back(globalPos(2));
0868       m_SL_lx.push_back(loc(0));
0869       m_SL_ly.push_back(loc(1));
0870
0871       /// Get corresponding truth hit position
0872       float gt = -9999;
0873
0874       Acts::Vector3 globalTruthPos = getGlobalTruthHit(key, gt);
0875       if(std::isnan(globalTruthPos(0)))
0876       {
0877         m_t_SL_lx.push_back(std::numeric_limits<float>::quiet_NaN());
0878         m_t_SL_ly.push_back(std::numeric_limits<float>::quiet_NaN());
0879         m_t_SL_gx.push_back(std::numeric_limits<float>::quiet_NaN());
0880         m_t_SL_gy.push_back(std::numeric_limits<float>::quiet_NaN());
0881         m_t_SL_gz.push_back(std::numeric_limits<float>::quiet_NaN());
0882         return;
0883       }
0884       float gx = globalTruthPos(0);
0885       float gy = globalTruthPos(1);
0886       float gz = globalTruthPos(2);
0887
0888       /// Get local truth position
0889       const float r = std::sqrt(gx * gx + gy * gy + gz * gz);
0890       Acts::Vector3 globalTruthUnitDir(gx / r, gy / r, gz / r);
0891
0892       auto surf = getSurface(key, cluster);
0893
0894       auto truthLocal = (*surf).globalToLocal(m_tGeometry->geometry().getGeoContext(),
0895                                               globalTruthPos,
0896                                               globalTruthUnitDir);
0897
0898       if (truthLocal.ok())
0899       {
0900         Acts::Vector2 truthLocalVec = truthLocal.value();
0901
0902         m_t_SL_lx.push_back(truthLocalVec(0));
0903         m_t_SL_ly.push_back(truthLocalVec(1));
0904         m_t_SL_gx.push_back(gx);
0905         m_t_SL_gy.push_back(gy);
0906         m_t_SL_gz.push_back(gz);
0907       }
0908       else
0909       {
0910         Acts::Vector3 loct = (*surf).transform(m_tGeometry->geometry().getGeoContext()).inverse() * globalTruthPos;
0911
0912         m_t_SL_lx.push_back(loct(0));
0913         m_t_SL_ly.push_back(loct(1));
0914         m_t_SL_gx.push_back(gx);
0915         m_t_SL_gy.push_back(gy);
0916         m_t_SL_gz.push_back(gz);
0917       }
0918     }
0919   }
0920   if (m_verbosity > 2)
0921   {
0922     std::cout << "Filled proto track" << std::endl;
0923   }
0924 }
0925
0926 void ActsEvaluator::fillFittedTrackParams(const Trajectory::IndexedParameters& paramsMap,
0927                                           const size_t& trackTip)
0928 {
0929   m_hasFittedParams = false;
0930
0931   if (m_verbosity > 2)
0932   {
0933     std::cout << "Filling fitted track parameters" << std::endl;
0934   }
0935
0936   /// If it has track parameters, fill the values
0937   if (paramsMap.find(trackTip) != paramsMap.end())
0938   {
0939     m_hasFittedParams = true;
0940     const auto& boundParam = paramsMap.find(trackTip)->second;
0941     const auto& parameter = boundParam.parameters();
0942     const auto& covariance = *boundParam.covariance();
0943     m_charge_fit = boundParam.charge();
0944     m_eLOC0_fit = parameter[Acts::eBoundLoc0];
0945     m_eLOC1_fit = parameter[Acts::eBoundLoc1];
0946     m_ePHI_fit = parameter[Acts::eBoundPhi];
0947     m_eTHETA_fit = parameter[Acts::eBoundTheta];
0948     m_eQOP_fit = parameter[Acts::eBoundQOverP];
0949     m_eT_fit = parameter[Acts::eBoundTime];
0950     m_err_eLOC0_fit =
0951         std::sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0));
0952     m_err_eLOC1_fit =
0953         std::sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1));
0954     m_err_ePHI_fit = std::sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi));
0955     m_err_eTHETA_fit =
0956         std::sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta));
0957     m_err_eQOP_fit = std::sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP));
0958     m_err_eT_fit = std::sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime));
0959
0960     m_px_fit = boundParam.momentum()(0);
0961     m_py_fit = boundParam.momentum()(1);
0962     m_pz_fit = boundParam.momentum()(2);
0963     m_x_fit = boundParam.position(m_tGeometry->geometry().getGeoContext())(0);
0964     m_y_fit = boundParam.position(m_tGeometry->geometry().getGeoContext())(1);
0965     m_z_fit = boundParam.position(m_tGeometry->geometry().getGeoContext())(2);
0966
0967     return;
0968   }
0969
0970   /// Otherwise mark it as a bad fit
0971   m_eLOC0_fit = -9999;
0972   m_eLOC1_fit = -9999;
0973   m_ePHI_fit = -9999;
0974   m_eTHETA_fit = -9999;
0975   m_eQOP_fit = -9999;
0976   m_eT_fit = -9999;
0977   m_charge_fit = -9999;
0978   m_err_eLOC0_fit = -9999;
0979   m_err_eLOC1_fit = -9999;
0980   m_err_ePHI_fit = -9999;
0981   m_err_eTHETA_fit = -9999;
0982   m_err_eQOP_fit = -9999;
0983   m_err_eT_fit = -9999;
0984   m_px_fit = -9999;
0985   m_py_fit = -9999;
0986   m_pz_fit = -9999;
0987   m_x_fit = -9999;
0988   m_y_fit = -9999;
0989   m_z_fit = -9999;
0990
0991   if (m_verbosity > 2)
0992   {
0993     std::cout << "Finished fitted track params" << std::endl;
0994   }
0995   return;
0996 }
0997
0998 void ActsEvaluator::fillG4Particle(PHG4Particle* part)
0999 {
1000   SvtxTruthEval* trutheval = m_svtxEvalStack->get_truth_eval();
1001
1002   if (part)
1003   {
1004     m_t_barcode = part->get_track_id();
1005     const auto pid = part->get_pid();
1006     m_t_charge = pid < 0 ? -1 : 1;
1007     const auto vtx = trutheval->get_vertex(part);
1008     m_t_vx = vtx->get_x() * Acts::UnitConstants::cm;
1009     m_t_vy = vtx->get_y() * Acts::UnitConstants::cm;
1010     m_t_vz = vtx->get_z() * Acts::UnitConstants::cm;
1011     if (m_verbosity > 1)
1012       std::cout << "truth vertex : (" << m_t_vx << ", " << m_t_vy
1013                 << ", " << m_t_vz << ")" << std::endl;
1014     m_t_px = part->get_px();
1015     m_t_py = part->get_py();
1016     m_t_pz = part->get_pz();
1017     const double p = std::sqrt(m_t_px * m_t_px + m_t_py * m_t_py + m_t_pz * m_t_pz);
1018     m_t_theta = std::acos(m_t_pz / p);
1019     m_t_phi = std::atan(m_t_py / m_t_px);
1020     m_t_pT = std::sqrt(m_t_px * m_t_px + m_t_py * m_t_py);
1021     m_t_eta = std::atanh(m_t_pz / p);
1022
1023     return;
1024   }
1025
1026   /// If particle doesn't exist, just fill with -999
1027   m_t_barcode = -9999;
1028   m_t_charge = -9999;
1029   m_t_vx = -9999;
1030   m_t_vy = -9999;
1031   m_t_vz = -9999;
1032   m_t_px = -9999;
1033   m_t_py = -9999;
1034   m_t_pz = -9999;
1035   m_t_theta = -9999;
1036   m_t_phi = -9999;
1037   m_t_pT = -9999;
1038   m_t_eta = -9999;
1039
1040   return;
1041 }
1042
1043 int ActsEvaluator::getNodes(PHCompositeNode* topNode)
1044 {
1045   m_tpcSeeds = findNode::getClass<TrackSeedContainer>(topNode, "TpcTrackSeedContainer");
1046   m_siliconSeeds = findNode::getClass<TrackSeedContainer>(topNode, "SiliconTrackSeedContainer");
1047
1048   if (!m_tpcSeeds or !m_siliconSeeds)
1049   {
1050     std::cout << PHWHERE << "Seed containers not found, cannot continue!"
1051               << std::endl;
1052     return Fun4AllReturnCodes::ABORTEVENT;
1053   }
1054
1055
1056
1057   m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
1058
1059   if (!m_tGeometry)
1060   {
1061     std::cout << PHWHERE << "No Acts Tracking geometry on node tree. Bailing"
1062               << std::endl;
1063
1064     return Fun4AllReturnCodes::ABORTEVENT;
1065   }
1066
1067   m_trackMap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
1068
1069   if (!m_trackMap)
1070   {
1071     std::cout << PHWHERE << "No SvtxTrackMap on node tree. Bailing."
1072               << std::endl;
1073
1074     return Fun4AllReturnCodes::ABORTEVENT;
1075   }
1076
1077   m_clusterContainer = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
1078   if (!m_clusterContainer)
1079   {
1080     std::cout << PHWHERE << "No clusters, bailing"
1081               << std::endl;
1082     return Fun4AllReturnCodes::ABORTEVENT;
1083   }
1084
1085   m_truthInfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
1086
1087   if (!m_truthInfo)
1088   {
1089     std::cout << PHWHERE << "PHG4TruthInfoContainer not found! If you are not running the Acts Evaluator on data, this will crash"
1090               << std::endl;
1091
1092   }
1093
1094   return Fun4AllReturnCodes::EVENT_OK;
1095 }
1096
1097 void ActsEvaluator::clearTrackVariables()
1098 {
1099   m_t_x.clear();
1100   m_t_y.clear();
1101   m_t_z.clear();
1102   m_t_r.clear();
1103   m_t_dx.clear();
1104   m_t_dy.clear();
1105   m_t_dz.clear();
1106   m_t_eLOC0.clear();
1107   m_t_eLOC1.clear();
1108   m_t_ePHI.clear();
1109   m_t_eTHETA.clear();
1110   m_t_eQOP.clear();
1111   m_t_eT.clear();
1112
1113   m_volumeID.clear();
1114   m_layerID.clear();
1115   m_moduleID.clear();
1116   m_sphenixlayer.clear();
1117   m_lx_hit.clear();
1118   m_ly_hit.clear();
1119   m_x_hit.clear();
1120   m_y_hit.clear();
1121   m_z_hit.clear();
1122   m_res_x_hit.clear();
1123   m_res_y_hit.clear();
1124   m_err_x_hit.clear();
1125   m_err_y_hit.clear();
1126   m_pull_x_hit.clear();
1127   m_pull_y_hit.clear();
1128   m_dim_hit.clear();
1129
1130   m_prt.clear();
1131   m_eLOC0_prt.clear();
1132   m_eLOC1_prt.clear();
1133   m_ePHI_prt.clear();
1134   m_eTHETA_prt.clear();
1135   m_eQOP_prt.clear();
1136   m_eT_prt.clear();
1137   m_res_eLOC0_prt.clear();
1138   m_res_eLOC1_prt.clear();
1139   m_res_ePHI_prt.clear();
1140   m_res_eTHETA_prt.clear();
1141   m_res_eQOP_prt.clear();
1142   m_res_eT_prt.clear();
1143   m_err_eLOC0_prt.clear();
1144   m_err_eLOC1_prt.clear();
1145   m_err_ePHI_prt.clear();
1146   m_err_eTHETA_prt.clear();
1147   m_err_eQOP_prt.clear();
1148   m_err_eT_prt.clear();
1149   m_pull_eLOC0_prt.clear();
1150   m_pull_eLOC1_prt.clear();
1151   m_pull_ePHI_prt.clear();
1152   m_pull_eTHETA_prt.clear();
1153   m_pull_eQOP_prt.clear();
1154   m_pull_eT_prt.clear();
1155   m_x_prt.clear();
1156   m_y_prt.clear();
1157   m_z_prt.clear();
1158   m_px_prt.clear();
1159   m_py_prt.clear();
1160   m_pz_prt.clear();
1161   m_eta_prt.clear();
1162   m_pT_prt.clear();
1163
1164   m_flt.clear();
1165   m_eLOC0_flt.clear();
1166   m_eLOC1_flt.clear();
1167   m_ePHI_flt.clear();
1168   m_eTHETA_flt.clear();
1169   m_eQOP_flt.clear();
1170   m_eT_flt.clear();
1171   m_res_eLOC0_flt.clear();
1172   m_res_eLOC1_flt.clear();
1173   m_res_ePHI_flt.clear();
1174   m_res_eTHETA_flt.clear();
1175   m_res_eQOP_flt.clear();
1176   m_res_eT_flt.clear();
1177   m_err_eLOC0_flt.clear();
1178   m_err_eLOC1_flt.clear();
1179   m_err_ePHI_flt.clear();
1180   m_err_eTHETA_flt.clear();
1181   m_err_eQOP_flt.clear();
1182   m_err_eT_flt.clear();
1183   m_pull_eLOC0_flt.clear();
1184   m_pull_eLOC1_flt.clear();
1185   m_pull_ePHI_flt.clear();
1186   m_pull_eTHETA_flt.clear();
1187   m_pull_eQOP_flt.clear();
1188   m_pull_eT_flt.clear();
1189   m_x_flt.clear();
1190   m_y_flt.clear();
1191   m_z_flt.clear();
1192   m_px_flt.clear();
1193   m_py_flt.clear();
1194   m_pz_flt.clear();
1195   m_eta_flt.clear();
1196   m_pT_flt.clear();
1197   m_chi2.clear();
1198
1199   m_smt.clear();
1200   m_eLOC0_smt.clear();
1201   m_eLOC1_smt.clear();
1202   m_ePHI_smt.clear();
1203   m_eTHETA_smt.clear();
1204   m_eQOP_smt.clear();
1205   m_eT_smt.clear();
1206   m_res_eLOC0_smt.clear();
1207   m_res_eLOC1_smt.clear();
1208   m_res_ePHI_smt.clear();
1209   m_res_eTHETA_smt.clear();
1210   m_res_eQOP_smt.clear();
1211   m_res_eT_smt.clear();
1212   m_err_eLOC0_smt.clear();
1213   m_err_eLOC1_smt.clear();
1214   m_err_ePHI_smt.clear();
1215   m_err_eTHETA_smt.clear();
1216   m_err_eQOP_smt.clear();
1217   m_err_eT_smt.clear();
1218   m_pull_eLOC0_smt.clear();
1219   m_pull_eLOC1_smt.clear();
1220   m_pull_ePHI_smt.clear();
1221   m_pull_eTHETA_smt.clear();
1222   m_pull_eQOP_smt.clear();
1223   m_pull_eT_smt.clear();
1224   m_x_smt.clear();
1225   m_y_smt.clear();
1226   m_z_smt.clear();
1227   m_px_smt.clear();
1228   m_py_smt.clear();
1229   m_pz_smt.clear();
1230   m_eta_smt.clear();
1231   m_pT_smt.clear();
1232
1233   m_SLx.clear();
1234   m_SLy.clear();
1235   m_SLz.clear();
1236   m_SL_lx.clear();
1237   m_SL_ly.clear();
1238   m_t_SL_lx.clear();
1239   m_t_SL_ly.clear();
1240   m_t_SL_gx.clear();
1241   m_t_SL_gy.clear();
1242   m_t_SL_gz.clear();
1243
1244   m_protoTrackPx = -9999.;
1245   m_protoTrackPy = -9999.;
1246   m_protoTrackPz = -9999.;
1247   m_protoTrackX = -9999.;
1248   m_protoTrackY = -9999.;
1249   m_protoTrackZ = -9999.;
1250   m_protoD0Cov = -9999.;
1251   m_protoZ0Cov = -9999.;
1252   m_protoPhiCov = -9999.;
1253   m_protoThetaCov = -9999.;
1254   m_protoQopCov = -9999.;
1255
1256   return;
1257 }
1258
1259 void ActsEvaluator::initializeTree()
1260 {
1261
1262   m_trackFile = new TFile(m_filename.c_str(), "RECREATE");
1263
1264   m_trackTree = new TTree("tracktree", "A tree with Acts KF track information");
1265
1266   m_trackTree->Branch("event_nr", &m_eventNr);
1267   m_trackTree->Branch("traj_nr", &m_trajNr);
1268   m_trackTree->Branch("track_nr", &m_trackNr);
1269   m_trackTree->Branch("t_barcode", &m_t_barcode, "t_barcode/l");
1270   m_trackTree->Branch("t_charge", &m_t_charge);
1271   m_trackTree->Branch("t_time", &m_t_time);
1272   m_trackTree->Branch("t_vx", &m_t_vx);
1273   m_trackTree->Branch("t_vy", &m_t_vy);
1274   m_trackTree->Branch("t_vz", &m_t_vz);
1275   m_trackTree->Branch("t_px", &m_t_px);
1276   m_trackTree->Branch("t_py", &m_t_py);
1277   m_trackTree->Branch("t_pz", &m_t_pz);
1278   m_trackTree->Branch("t_theta", &m_t_theta);
1279   m_trackTree->Branch("t_phi", &m_t_phi);
1280   m_trackTree->Branch("t_eta", &m_t_eta);
1281   m_trackTree->Branch("t_pT", &m_t_pT);
1282
1283   m_trackTree->Branch("t_x", &m_t_x);
1284   m_trackTree->Branch("t_y", &m_t_y);
1285   m_trackTree->Branch("t_z", &m_t_z);
1286   m_trackTree->Branch("t_r", &m_t_r);
1287   m_trackTree->Branch("t_dx", &m_t_dx);
1288   m_trackTree->Branch("t_dy", &m_t_dy);
1289   m_trackTree->Branch("t_dz", &m_t_dz);
1290   m_trackTree->Branch("t_eLOC0", &m_t_eLOC0);
1291   m_trackTree->Branch("t_eLOC1", &m_t_eLOC1);
1292   m_trackTree->Branch("t_ePHI", &m_t_ePHI);
1293   m_trackTree->Branch("t_eTHETA", &m_t_eTHETA);
1294   m_trackTree->Branch("t_eQOP", &m_t_eQOP);
1295   m_trackTree->Branch("t_eT", &m_t_eT);
1296
1297   m_trackTree->Branch("hasFittedParams", &m_hasFittedParams);
1298   m_trackTree->Branch("chi2_fit", &m_chi2_fit);
1299   m_trackTree->Branch("quality", &m_quality);
1300   m_trackTree->Branch("ndf_fit", &m_ndf_fit);
1301   m_trackTree->Branch("eLOC0_fit", &m_eLOC0_fit);
1302   m_trackTree->Branch("eLOC1_fit", &m_eLOC1_fit);
1303   m_trackTree->Branch("ePHI_fit", &m_ePHI_fit);
1304   m_trackTree->Branch("eTHETA_fit", &m_eTHETA_fit);
1305   m_trackTree->Branch("eQOP_fit", &m_eQOP_fit);
1306   m_trackTree->Branch("eT_fit", &m_eT_fit);
1307   m_trackTree->Branch("charge_fit", &m_charge_fit);
1308   m_trackTree->Branch("err_eLOC0_fit", &m_err_eLOC0_fit);
1309   m_trackTree->Branch("err_eLOC1_fit", &m_err_eLOC1_fit);
1310   m_trackTree->Branch("err_ePHI_fit", &m_err_ePHI_fit);
1311   m_trackTree->Branch("err_eTHETA_fit", &m_err_eTHETA_fit);
1312   m_trackTree->Branch("err_eQOP_fit", &m_err_eQOP_fit);
1313   m_trackTree->Branch("err_eT_fit", &m_err_eT_fit);
1314   m_trackTree->Branch("g_px_fit", &m_px_fit);
1315   m_trackTree->Branch("g_py_fit", &m_py_fit);
1316   m_trackTree->Branch("g_pz_fit", &m_pz_fit);
1317   m_trackTree->Branch("g_x_fit", &m_x_fit);
1318   m_trackTree->Branch("g_y_fit", &m_y_fit);
1319   m_trackTree->Branch("g_z_fit", &m_z_fit);
1320   m_trackTree->Branch("g_dca3Dxy_fit", &m_dca3Dxy);
1321   m_trackTree->Branch("g_dca3Dz_fit", &m_dca3Dz);
1322   m_trackTree->Branch("g_dca3Dxy_cov", &m_dca3DxyCov);
1323   m_trackTree->Branch("g_dca3Dz_cov", &m_dca3DzCov);
1324
1325   m_trackTree->Branch("g_protoTrackPx", &m_protoTrackPx);
1326   m_trackTree->Branch("g_protoTrackPy", &m_protoTrackPy);
1327   m_trackTree->Branch("g_protoTrackPz", &m_protoTrackPz);
1328   m_trackTree->Branch("g_protoTrackX", &m_protoTrackX);
1329   m_trackTree->Branch("g_protoTrackY", &m_protoTrackY);
1330   m_trackTree->Branch("g_protoTrackZ", &m_protoTrackZ);
1331   m_trackTree->Branch("g_protoTrackD0Cov", &m_protoD0Cov);
1332   m_trackTree->Branch("g_protoTrackZ0Cov", &m_protoZ0Cov);
1333   m_trackTree->Branch("g_protoTrackPhiCov", &m_protoPhiCov);
1334   m_trackTree->Branch("g_protoTrackThetaCov", &m_protoThetaCov);
1335   m_trackTree->Branch("g_protoTrackQopCov", &m_protoQopCov);
1336   m_trackTree->Branch("g_SLx", &m_SLx);
1337   m_trackTree->Branch("g_SLy", &m_SLy);
1338   m_trackTree->Branch("g_SLz", &m_SLz);
1339   m_trackTree->Branch("l_SLx", &m_SL_lx);
1340   m_trackTree->Branch("l_SLy", &m_SL_ly);
1341   m_trackTree->Branch("t_SL_lx", &m_t_SL_lx);
1342   m_trackTree->Branch("t_SL_ly", &m_t_SL_ly);
1343   m_trackTree->Branch("t_SL_gx", &m_t_SL_gx);
1344   m_trackTree->Branch("t_SL_gy", &m_t_SL_gy);
1345   m_trackTree->Branch("t_SL_gz", &m_t_SL_gz);
1346
1347   m_trackTree->Branch("nSharedHits", &m_nSharedHits);
1348   m_trackTree->Branch("nHoles", &m_nHoles);
1349   m_trackTree->Branch("nOutliers", &m_nOutliers);
1350   m_trackTree->Branch("nStates", &m_nStates);
1351   m_trackTree->Branch("nMeasurements", &m_nMeasurements);
1352   m_trackTree->Branch("volume_id", &m_volumeID);
1353   m_trackTree->Branch("layer_id", &m_layerID);
1354   m_trackTree->Branch("module_id", &m_moduleID);
1355   m_trackTree->Branch("sphenixlayer",&m_sphenixlayer);
1356   m_trackTree->Branch("l_x_hit", &m_lx_hit);
1357   m_trackTree->Branch("l_y_hit", &m_ly_hit);
1358   m_trackTree->Branch("g_x_hit", &m_x_hit);
1359   m_trackTree->Branch("g_y_hit", &m_y_hit);
1360   m_trackTree->Branch("g_z_hit", &m_z_hit);
1361   m_trackTree->Branch("res_x_hit", &m_res_x_hit);
1362   m_trackTree->Branch("res_y_hit", &m_res_y_hit);
1363   m_trackTree->Branch("err_x_hit", &m_err_x_hit);
1364   m_trackTree->Branch("err_y_hit", &m_err_y_hit);
1365   m_trackTree->Branch("pull_x_hit", &m_pull_x_hit);
1366   m_trackTree->Branch("pull_y_hit", &m_pull_y_hit);
1367   m_trackTree->Branch("dim_hit", &m_dim_hit);
1368
1369   m_trackTree->Branch("nPredicted", &m_nPredicted);
1370   m_trackTree->Branch("predicted", &m_prt);
1371   m_trackTree->Branch("eLOC0_prt", &m_eLOC0_prt);
1372   m_trackTree->Branch("eLOC1_prt", &m_eLOC1_prt);
1373   m_trackTree->Branch("ePHI_prt", &m_ePHI_prt);
1374   m_trackTree->Branch("eTHETA_prt", &m_eTHETA_prt);
1375   m_trackTree->Branch("eQOP_prt", &m_eQOP_prt);
1376   m_trackTree->Branch("eT_prt", &m_eT_prt);
1377   m_trackTree->Branch("res_eLOC0_prt", &m_res_eLOC0_prt);
1378   m_trackTree->Branch("res_eLOC1_prt", &m_res_eLOC1_prt);
1379   m_trackTree->Branch("res_ePHI_prt", &m_res_ePHI_prt);
1380   m_trackTree->Branch("res_eTHETA_prt", &m_res_eTHETA_prt);
1381   m_trackTree->Branch("res_eQOP_prt", &m_res_eQOP_prt);
1382   m_trackTree->Branch("res_eT_prt", &m_res_eT_prt);
1383   m_trackTree->Branch("err_eLOC0_prt", &m_err_eLOC0_prt);
1384   m_trackTree->Branch("err_eLOC1_prt", &m_err_eLOC1_prt);
1385   m_trackTree->Branch("err_ePHI_prt", &m_err_ePHI_prt);
1386   m_trackTree->Branch("err_eTHETA_prt", &m_err_eTHETA_prt);
1387   m_trackTree->Branch("err_eQOP_prt", &m_err_eQOP_prt);
1388   m_trackTree->Branch("err_eT_prt", &m_err_eT_prt);
1389   m_trackTree->Branch("pull_eLOC0_prt", &m_pull_eLOC0_prt);
1390   m_trackTree->Branch("pull_eLOC1_prt", &m_pull_eLOC1_prt);
1391   m_trackTree->Branch("pull_ePHI_prt", &m_pull_ePHI_prt);
1392   m_trackTree->Branch("pull_eTHETA_prt", &m_pull_eTHETA_prt);
1393   m_trackTree->Branch("pull_eQOP_prt", &m_pull_eQOP_prt);
1394   m_trackTree->Branch("pull_eT_prt", &m_pull_eT_prt);
1395   m_trackTree->Branch("g_x_prt", &m_x_prt);
1396   m_trackTree->Branch("g_y_prt", &m_y_prt);
1397   m_trackTree->Branch("g_z_prt", &m_z_prt);
1398   m_trackTree->Branch("px_prt", &m_px_prt);
1399   m_trackTree->Branch("py_prt", &m_py_prt);
1400   m_trackTree->Branch("pz_prt", &m_pz_prt);
1401   m_trackTree->Branch("eta_prt", &m_eta_prt);
1402   m_trackTree->Branch("pT_prt", &m_pT_prt);
1403
1404   m_trackTree->Branch("nFiltered", &m_nFiltered);
1405   m_trackTree->Branch("filtered", &m_flt);
1406   m_trackTree->Branch("eLOC0_flt", &m_eLOC0_flt);
1407   m_trackTree->Branch("eLOC1_flt", &m_eLOC1_flt);
1408   m_trackTree->Branch("ePHI_flt", &m_ePHI_flt);
1409   m_trackTree->Branch("eTHETA_flt", &m_eTHETA_flt);
1410   m_trackTree->Branch("eQOP_flt", &m_eQOP_flt);
1411   m_trackTree->Branch("eT_flt", &m_eT_flt);
1412   m_trackTree->Branch("res_eLOC0_flt", &m_res_eLOC0_flt);
1413   m_trackTree->Branch("res_eLOC1_flt", &m_res_eLOC1_flt);
1414   m_trackTree->Branch("res_ePHI_flt", &m_res_ePHI_flt);
1415   m_trackTree->Branch("res_eTHETA_flt", &m_res_eTHETA_flt);
1416   m_trackTree->Branch("res_eQOP_flt", &m_res_eQOP_flt);
1417   m_trackTree->Branch("res_eT_flt", &m_res_eT_flt);
1418   m_trackTree->Branch("err_eLOC0_flt", &m_err_eLOC0_flt);
1419   m_trackTree->Branch("err_eLOC1_flt", &m_err_eLOC1_flt);
1420   m_trackTree->Branch("err_ePHI_flt", &m_err_ePHI_flt);
1421   m_trackTree->Branch("err_eTHETA_flt", &m_err_eTHETA_flt);
1422   m_trackTree->Branch("err_eQOP_flt", &m_err_eQOP_flt);
1423   m_trackTree->Branch("err_eT_flt", &m_err_eT_flt);
1424   m_trackTree->Branch("pull_eLOC0_flt", &m_pull_eLOC0_flt);
1425   m_trackTree->Branch("pull_eLOC1_flt", &m_pull_eLOC1_flt);
1426   m_trackTree->Branch("pull_ePHI_flt", &m_pull_ePHI_flt);
1427   m_trackTree->Branch("pull_eTHETA_flt", &m_pull_eTHETA_flt);
1428   m_trackTree->Branch("pull_eQOP_flt", &m_pull_eQOP_flt);
1429   m_trackTree->Branch("pull_eT_flt", &m_pull_eT_flt);
1430   m_trackTree->Branch("g_x_flt", &m_x_flt);
1431   m_trackTree->Branch("g_y_flt", &m_y_flt);
1432   m_trackTree->Branch("g_z_flt", &m_z_flt);
1433   m_trackTree->Branch("px_flt", &m_px_flt);
1434   m_trackTree->Branch("py_flt", &m_py_flt);
1435   m_trackTree->Branch("pz_flt", &m_pz_flt);
1436   m_trackTree->Branch("eta_flt", &m_eta_flt);
1437   m_trackTree->Branch("pT_flt", &m_pT_flt);
1438   m_trackTree->Branch("chi2", &m_chi2);
1439
1440   m_trackTree->Branch("nSmoothed", &m_nSmoothed);
1441   m_trackTree->Branch("smoothed", &m_smt);
1442   m_trackTree->Branch("eLOC0_smt", &m_eLOC0_smt);
1443   m_trackTree->Branch("eLOC1_smt", &m_eLOC1_smt);
1444   m_trackTree->Branch("ePHI_smt", &m_ePHI_smt);
1445   m_trackTree->Branch("eTHETA_smt", &m_eTHETA_smt);
1446   m_trackTree->Branch("eQOP_smt", &m_eQOP_smt);
1447   m_trackTree->Branch("eT_smt", &m_eT_smt);
1448   m_trackTree->Branch("res_eLOC0_smt", &m_res_eLOC0_smt);
1449   m_trackTree->Branch("res_eLOC1_smt", &m_res_eLOC1_smt);
1450   m_trackTree->Branch("res_ePHI_smt", &m_res_ePHI_smt);
1451   m_trackTree->Branch("res_eTHETA_smt", &m_res_eTHETA_smt);
1452   m_trackTree->Branch("res_eQOP_smt", &m_res_eQOP_smt);
1453   m_trackTree->Branch("res_eT_smt", &m_res_eT_smt);
1454   m_trackTree->Branch("err_eLOC0_smt", &m_err_eLOC0_smt);
1455   m_trackTree->Branch("err_eLOC1_smt", &m_err_eLOC1_smt);
1456   m_trackTree->Branch("err_ePHI_smt", &m_err_ePHI_smt);
1457   m_trackTree->Branch("err_eTHETA_smt", &m_err_eTHETA_smt);
1458   m_trackTree->Branch("err_eQOP_smt", &m_err_eQOP_smt);
1459   m_trackTree->Branch("err_eT_smt", &m_err_eT_smt);
1460   m_trackTree->Branch("pull_eLOC0_smt", &m_pull_eLOC0_smt);
1461   m_trackTree->Branch("pull_eLOC1_smt", &m_pull_eLOC1_smt);
1462   m_trackTree->Branch("pull_ePHI_smt", &m_pull_ePHI_smt);
1463   m_trackTree->Branch("pull_eTHETA_smt", &m_pull_eTHETA_smt);
1464   m_trackTree->Branch("pull_eQOP_smt", &m_pull_eQOP_smt);
1465   m_trackTree->Branch("pull_eT_smt", &m_pull_eT_smt);
1466   m_trackTree->Branch("g_x_smt", &m_x_smt);
1467   m_trackTree->Branch("g_y_smt", &m_y_smt);
1468   m_trackTree->Branch("g_z_smt", &m_z_smt);
1469   m_trackTree->Branch("px_smt", &m_px_smt);
1470   m_trackTree->Branch("py_smt", &m_py_smt);
1471   m_trackTree->Branch("pz_smt", &m_pz_smt);
1472   m_trackTree->Branch("eta_smt", &m_eta_smt);
1473   m_trackTree->Branch("pT_smt", &m_pT_smt);
1474 }