KalmanVertexFit/src/VertexFitterResult.cc

0001 #include "RecoVertex/KalmanVertexFit/interface/VertexFitterResult.h"
0002 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0003 #include "TrackingTools/TransientTrack/interface/TrackTransientTrack.h"
0004 #include "TrackingTools/PatternTools/interface/trackingParametersAtClosestApproachToBeamSpot.h"
0005
0006 using namespace reco;
0007 using namespace std;
0008
0009 VertexFitterResult::VertexFitterResult(const int maxTracks, const MagneticField* magField) : theMagField(magField) {
0010   theMaxTracks = maxTracks;
0011   if (theMagField == nullptr)
0012     theMaxTracks = 0;
0013   for (int i = 0; i < 5; i++) {
0014     if (maxTracks > 0) {
0015       simPars[i] = new float[maxTracks];
0016       recPars[i] = new float[maxTracks];
0017       refPars[i] = new float[maxTracks];
0018       recErrs[i] = new float[maxTracks];
0019       refErrs[i] = new float[maxTracks];
0020     } else {
0021       simPars[i] = nullptr;
0022       recPars[i] = nullptr;
0023       refPars[i] = nullptr;
0024       recErrs[i] = nullptr;
0025       refErrs[i] = nullptr;
0026     }
0027   }
0028   trackWeight = new float[maxTracks];
0029   simIndex = new int[maxTracks];
0030   recIndex = new int[maxTracks];
0031   numberOfRecTracks = theMaxTracks;
0032   numberOfSimTracks = theMaxTracks;
0033   reset();
0034 }
0035
0036 VertexFitterResult::~VertexFitterResult() {
0037   //
0038   // delete arrays
0039   //
0040   for (int i = 0; i < 5; i++) {
0041     delete[] simPars[i];
0042     delete[] recPars[i];
0043     delete[] refPars[i];
0044     delete[] recErrs[i];
0045     delete[] refErrs[i];
0046   }
0047   delete trackWeight;
0048   delete simIndex;
0049   delete recIndex;
0050 }
0051
0052 void VertexFitterResult::fill(const TransientVertex& recVertex,
0053                               const TrackingVertex* simv,
0054                               reco::RecoToSimCollection* recSimColl,
0055                               const float& time) {
0056   TTrackCont recTrackV;
0057   if (recVertex.isValid())
0058     recTrackV = recVertex.originalTracks();
0059   fill(recVertex, recTrackV, simv, recSimColl, time);
0060 }
0061
0062 void VertexFitterResult::fill(const TransientVertex& recVertex,
0063                               const TTrackCont& recTrackV,
0064                               const TrackingVertex* simv,
0065                               reco::RecoToSimCollection* recSimColl,
0066                               const float& time) {
0067   TrackingParticleRefVector simTrackV;
0068
0069   Basic3DVector<double> vert;
0070   if (recVertex.isValid()) {
0071     recPos[0] = recVertex.position().x();
0072     recPos[1] = recVertex.position().y();
0073     recPos[2] = recVertex.position().z();
0074
0075     recErr[0] = sqrt(recVertex.positionError().cxx());
0076     recErr[1] = sqrt(recVertex.positionError().cyy());
0077     recErr[2] = sqrt(recVertex.positionError().czz());
0078     vert = (Basic3DVector<double>)recVertex.position();
0079
0080     chi[0] = recVertex.totalChiSquared();
0081     chi[1] = recVertex.degreesOfFreedom();
0082     chi[2] = ChiSquaredProbability(recVertex.totalChiSquared(), recVertex.degreesOfFreedom());
0083     vertex = 2;
0084     fitTime = time;
0085     tracks[1] = recVertex.originalTracks().size();
0086   }
0087
0088   if (simv != nullptr) {
0089     simPos[0] = simv->position().x();
0090     simPos[1] = simv->position().y();
0091     simPos[2] = simv->position().z();
0092
0093     simTrackV = simv->daughterTracks();
0094     vertex += 1;
0095     for (TrackingVertex::tp_iterator simTrack = simv->daughterTracks_begin();
0096          (simTrack != simv->daughterTracks_end() && (numberOfSimTracks < theMaxTracks));
0097          simTrack++) {
0098       Basic3DVector<double> momAtVtx((**simTrack).momentum());
0099
0100       std::pair<bool, reco::TrackBase::ParameterVector> paramPair = reco::trackingParametersAtClosestApproachToBeamSpot(
0101           vert, momAtVtx, (float)(**simTrack).charge(), *theMagField, recTrackV.front().stateAtBeamLine().beamSpot());
0102       if (paramPair.first) {
0103         fillParameters(paramPair.second, simPars, numberOfSimTracks);
0104         simIndex[numberOfSimTracks] = -1;
0105         ++numberOfSimTracks;
0106       }
0107     }
0108     tracks[0] = numberOfSimTracks;
0109   }
0110
0111   // now store all the recTrack...
0112
0113   for (TTrackCont::const_iterator recTrack = recTrackV.begin();
0114        (recTrack != recTrackV.end() && (numberOfRecTracks < theMaxTracks));
0115        recTrack++) {
0116     //    std::cout << "Input; 1/Pt " << 1./(*recTrack).momentumAtVertex().transverse() << std::endl;
0117
0118     //looking for sim tracks corresponding to our reconstructed tracks:
0119     simIndex[numberOfRecTracks] = -1;
0120
0121     std::vector<std::pair<TrackingParticleRef, double> > simFound;
0122     try {
0123       const TrackTransientTrack* ttt = dynamic_cast<const TrackTransientTrack*>(recTrack->basicTransientTrack());
0124       if ((ttt != nullptr) && (recSimColl != nullptr))
0125         simFound = (*recSimColl)[ttt->trackBaseRef()];
0126       //       if (recSimColl!=0) simFound = (*recSimColl)[recTrack->persistentTrackRef()];
0127       //      if (recSimColl!=0) simFound = (*recSimColl)[recTrack];
0128
0129     } catch (cms::Exception const& e) {
0130       //       LogDebug("TrackValidator") << "reco::Track #" << rT << " with pt=" << track->pt()
0131       //                 << " NOT associated to any TrackingParticle" << "\n";
0132       //       edm::LogError("TrackValidator") << e.what() << "\n";
0133     }
0134
0135     if (!simFound.empty()) {
0136       //OK, it was associated, so get the state on the same surface as the 'SimState'
0137       TrackingParticleRefVector::const_iterator simTrackI = find(simTrackV.begin(), simTrackV.end(), simFound[0].first);
0138       if (simTrackI != simTrackV.end())
0139         ++tracks[2];
0140       int simTrackIndex = simTrackI - simTrackV.begin();
0141       if (simTrackIndex < numberOfSimTracks) {
0142         simIndex[numberOfRecTracks] = simTrackIndex;
0143         recIndex[simTrackIndex] = numberOfRecTracks;
0144         //  cout << "Assoc; 1/Pt " << 1./(*recTrack).momentumAtVertex().transverse() << std::endl;
0145       }
0146     }
0147
0148     TrajectoryStateClosestToPoint tscp = recTrack->trajectoryStateClosestToPoint(recVertex.position());
0149     fillParameters(recTrack->track().parameters(), recPars, numberOfRecTracks);
0150     fillErrors(tscp.perigeeError(), recErrs, numberOfRecTracks);
0151     //     trackWeight[numberOfRecTracks] = recVertex.trackWeight(*recTrack);
0152     //
0153     //     if ((recVertex.isValid())&&(recVertex.hasRefittedTracks())) {
0154     //       //looking for corresponding refitted tracks:
0155     //       TrajectoryStateOnSurface refip;
0156     //       RecTrackCont::iterator refTrackI =
0157     //              find_if(refTracks.begin(), refTracks.end(), RecTrackMatch(*recTrack));
0158     //       if (refTrackI!=refTracks.end()) {
0159     //         // If it was not found, it would mean that it was not used in the fit,
0160     //  // or with a low weight such that the track was then discarded.
0161     //  if(simFound.size() != 0) {
0162     //    refip = refTrackI->stateOnSurface(simFound[0]->impactPointStateOnSurface().surface());
0163     //  } else {
0164     //           refip = refTrackI->innermostState();
0165     //  }
0166     //
0167     //  fillParameters(refip, refPars, numberOfRecTracks);
0168     //  fillErrors(refip, refErrs, numberOfRecTracks);
0169     //       }
0170     //     }
0171     //
0172     ++numberOfRecTracks;
0173   }
0174 }
0175
0176 void VertexFitterResult::fillParameters(const reco::TrackBase::ParameterVector& perigee,
0177                                         float* params[5],
0178                                         int trackNumber) {
0179   params[0][trackNumber] = perigee[0];
0180   params[1][trackNumber] = perigee[1];
0181   params[2][trackNumber] = perigee[2];
0182   params[3][trackNumber] = perigee[3];
0183   params[4][trackNumber] = perigee[4];
0184 }
0185
0186 void VertexFitterResult::fillParameters(const PerigeeTrajectoryParameters& ptp, float* params[5], int trackNumber) {
0187   const AlgebraicVector5& perigee = ptp.vector();
0188   params[0][trackNumber] = perigee[0];
0189   params[1][trackNumber] = perigee[1];
0190   params[2][trackNumber] = perigee[2];
0191   params[3][trackNumber] = perigee[3];
0192   params[4][trackNumber] = perigee[4];
0193 }
0194
0195 void VertexFitterResult::fillErrors(const PerigeeTrajectoryError& pte, float* errors[5], int trackNumber) {
0196   errors[0][trackNumber] = pte.transverseCurvatureError();
0197   errors[1][trackNumber] = pte.thetaError();
0198   errors[2][trackNumber] = pte.phiError();
0199   errors[3][trackNumber] = pte.transverseImpactParameterError();
0200   errors[4][trackNumber] = pte.longitudinalImpactParameterError();
0201 }
0202
0203 void VertexFitterResult::reset() {
0204   for (int i = 0; i < 3; ++i) {
0205     simPos[i] = 0.;
0206     recPos[i] = 0.;
0207     recErr[i] = 0.;
0208     chi[i] = 0.;
0209     tracks[i] = 0;
0210   }
0211   vertex = 0;
0212   fitTime = 0;
0213
0214   for (int j = 0; j < numberOfRecTracks; ++j) {
0215     for (int i = 0; i < 5; ++i) {
0216       recPars[i][j] = 0;
0217       refPars[i][j] = 0;
0218       recErrs[i][j] = 0;
0219       refErrs[i][j] = 0;
0220     }
0221     trackWeight[j] = 0;
0222     simIndex[j] = -1;
0223   }
0224   for (int j = 0; j < numberOfSimTracks; ++j) {
0225     for (int i = 0; i < 5; ++i) {
0226       simPars[i][j] = 0;
0227     }
0228     recIndex[j] = -1;
0229   }
0230
0231   numberOfRecTracks = 0;
0232   numberOfSimTracks = 0;
0233 }