VertexFitterResult.cc

Go to the documentation of this file.

#include "RecoVertex/KalmanVertexFit/interface/VertexFitterResult.h"
#include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
#include "TrackingTools/TransientTrack/interface/TrackTransientTrack.h"
#include "TrackingTools/PatternTools/interface/trackingParametersAtClosestApproachToBeamSpot.h"

using namespace reco;
using namespace std;

VertexFitterResult::VertexFitterResult(const int maxTracks, const MagneticField* magField)
  : theMagField(magField)
{
  theMaxTracks = maxTracks;
  if (theMagField==nullptr) theMaxTracks=0;
  for ( int i=0; i<5; i++ ) {
    if ( maxTracks>0 ) {
      simPars[i] = new float[maxTracks];
      recPars[i] = new float[maxTracks];
      refPars[i] = new float[maxTracks];
      recErrs[i] = new float[maxTracks];
      refErrs[i] = new float[maxTracks];
    }
    else {
      simPars[i] = nullptr;
      recPars[i] = nullptr;
      refPars[i] = nullptr;
      recErrs[i] = nullptr;
      refErrs[i] = nullptr;
    }
  }
  trackWeight = new float[maxTracks];
  simIndex = new int[maxTracks];
  recIndex = new int[maxTracks];
  numberOfRecTracks=theMaxTracks;
  numberOfSimTracks=theMaxTracks;
  reset();
}

VertexFitterResult::~VertexFitterResult()
{
    //
    // delete arrays
    //
    for ( int i=0; i<5; i++ ) {
      delete [] simPars[i];
      delete [] recPars[i];
      delete [] refPars[i];
      delete [] recErrs[i];
      delete [] refErrs[i];
    }
    delete trackWeight;
    delete simIndex;
    delete recIndex;
}

void VertexFitterResult::fill(const TransientVertex & recVertex,
    const TrackingVertex * simv, reco::RecoToSimCollection *recSimColl,
    const float &time) 
{
  TTrackCont recTrackV;
  if (recVertex.isValid()) recTrackV = recVertex.originalTracks();
  fill(recVertex, recTrackV, simv, recSimColl, time);
}

void VertexFitterResult::fill(const TransientVertex & recVertex, 
    const TTrackCont & recTrackV, const TrackingVertex * simv,
    reco::RecoToSimCollection *recSimColl, const float &time)
{
  TrackingParticleRefVector simTrackV;

  Basic3DVector<double> vert;
  if (recVertex.isValid()) {
    recPos[0] = recVertex.position().x();
    recPos[1] = recVertex.position().y();
    recPos[2] = recVertex.position().z();

    recErr[0] = sqrt(recVertex.positionError().cxx());
    recErr[1] = sqrt(recVertex.positionError().cyy());
    recErr[2] = sqrt(recVertex.positionError().czz());
    vert = (Basic3DVector<double>) recVertex.position();

    chi[0] = recVertex.totalChiSquared();
    chi[1] = recVertex.degreesOfFreedom();
    chi[2] = ChiSquaredProbability(recVertex.totalChiSquared(), 
                       recVertex.degreesOfFreedom());
    vertex = 2;
    fitTime = time;
    tracks[1] = recVertex.originalTracks().size();
  }

  if (simv!=nullptr) {
    simPos[0] = simv->position().x();
    simPos[1] = simv->position().y();
    simPos[2] = simv->position().z();

    simTrackV = simv->daughterTracks();
    vertex += 1;
    for(TrackingVertex::tp_iterator simTrack = simv->daughterTracks_begin();
                 (simTrack != simv->daughterTracks_end() && (numberOfSimTracks<theMaxTracks));
         simTrack++) {
      
      Basic3DVector<double> momAtVtx((**simTrack).momentum());

      std::pair<bool, reco::TrackBase::ParameterVector> paramPair =
    reco::trackingParametersAtClosestApproachToBeamSpot(vert, momAtVtx, 
                                                            (float) (**simTrack).charge(), 
                                                            *theMagField, 
                                                            recTrackV.front().stateAtBeamLine().beamSpot());
        if (paramPair.first) {
      fillParameters(paramPair.second, simPars, numberOfSimTracks);
      simIndex[numberOfSimTracks] = -1;
      ++numberOfSimTracks;
        }
    }
    tracks[0] = numberOfSimTracks;
  }


  // now store all the recTrack...  

  for(TTrackCont::const_iterator recTrack =recTrackV.begin();
               (recTrack != recTrackV.end() 
        && (numberOfRecTracks<theMaxTracks));
           recTrack++) {
    //    std::cout << "Input; 1/Pt " << 1./(*recTrack).momentumAtVertex().transverse() << std::endl;

    //looking for sim tracks corresponding to our reconstructed tracks:  
    simIndex[numberOfRecTracks] = -1;

    std::vector<std::pair<TrackingParticleRef, double> > simFound;
    try {
      const TrackTransientTrack* ttt = dynamic_cast<const TrackTransientTrack*>(recTrack->basicTransientTrack());
      if ((ttt!=nullptr) && (recSimColl!=nullptr)) simFound = (*recSimColl)[ttt->trackBaseRef()];
//       if (recSimColl!=0) simFound = (*recSimColl)[recTrack->persistentTrackRef()];
//      if (recSimColl!=0) simFound = (*recSimColl)[recTrack];

    } catch (cms::Exception const& e) {
//       LogDebug("TrackValidator") << "reco::Track #" << rT << " with pt=" << track->pt() 
//               << " NOT associated to any TrackingParticle" << "\n";
//       edm::LogError("TrackValidator") << e.what() << "\n";
    }

    if(!simFound.empty()) {
      //OK, it was associated, so get the state on the same surface as the 'SimState'
      TrackingParticleRefVector::const_iterator simTrackI = 
    find(simTrackV.begin(), simTrackV.end(), simFound[0].first);
      if (simTrackI!=simTrackV.end()) ++tracks[2];
      int simTrackIndex = simTrackI-simTrackV.begin();
      if (simTrackIndex<numberOfSimTracks) {
        simIndex[numberOfRecTracks] = simTrackIndex;
        recIndex[simTrackIndex] = numberOfRecTracks;
    //  cout << "Assoc; 1/Pt " << 1./(*recTrack).momentumAtVertex().transverse() << std::endl;
      }
    }

    TrajectoryStateClosestToPoint tscp = recTrack->trajectoryStateClosestToPoint(recVertex.position());
    fillParameters(recTrack->track().parameters(), recPars, numberOfRecTracks);
    fillErrors(tscp.perigeeError(), recErrs, numberOfRecTracks);
//     trackWeight[numberOfRecTracks] = recVertex.trackWeight(*recTrack);
// 
//     if ((recVertex.isValid())&&(recVertex.hasRefittedTracks())) {
//       //looking for corresponding refitted tracks:
//       TrajectoryStateOnSurface refip;
//       RecTrackCont::iterator refTrackI = 
//              find_if(refTracks.begin(), refTracks.end(), RecTrackMatch(*recTrack));
//       if (refTrackI!=refTracks.end()) {
//         // If it was not found, it would mean that it was not used in the fit,
//  // or with a low weight such that the track was then discarded.
//  if(simFound.size() != 0) {
//    refip = refTrackI->stateOnSurface(simFound[0]->impactPointStateOnSurface().surface());
//  } else {
//           refip = refTrackI->innermostState();
//  }
// 
//  fillParameters(refip, refPars, numberOfRecTracks);
//  fillErrors(refip, refErrs, numberOfRecTracks);
//       }
//     }
// 
    ++numberOfRecTracks;
  }
  
}

void VertexFitterResult::fillParameters (const reco::TrackBase::ParameterVector& perigee,
    float* params[5], int trackNumber)
{
  params[0][trackNumber] = perigee[0];
  params[1][trackNumber] = perigee[1];
  params[2][trackNumber] = perigee[2];
  params[3][trackNumber] = perigee[3];
  params[4][trackNumber] = perigee[4];
}

void VertexFitterResult::fillParameters (const PerigeeTrajectoryParameters & ptp,
    float* params[5], int trackNumber)
{
  const AlgebraicVector5 & perigee = ptp.vector();
  params[0][trackNumber] = perigee[0];
  params[1][trackNumber] = perigee[1];
  params[2][trackNumber] = perigee[2];
  params[3][trackNumber] = perigee[3];
  params[4][trackNumber] = perigee[4];
}

void VertexFitterResult::fillErrors (const PerigeeTrajectoryError & pte,
    float* errors[5], int trackNumber)
{
  errors[0][trackNumber] = pte.transverseCurvatureError(); 
  errors[1][trackNumber] = pte.thetaError(); 
  errors[2][trackNumber] = pte.phiError(); 
  errors[3][trackNumber] = pte.transverseImpactParameterError(); 
  errors[4][trackNumber] = pte.longitudinalImpactParameterError(); 
}

void VertexFitterResult::reset()
{
  for ( int i=0; i<3; ++i ) {
    simPos[i] = 0.;
    recPos[i] = 0.;
    recErr[i] = 0.;
    chi[i] = 0.;
    tracks[i] = 0;
  }
  vertex =0;
  fitTime = 0;

  for ( int j=0; j<numberOfRecTracks; ++j ) {
    for ( int i=0; i<5; ++i ) {
       recPars[i][j] = 0;
       refPars[i][j] = 0;
       recErrs[i][j] = 0;
       refErrs[i][j] = 0;
    }
    trackWeight[j] = 0;
    simIndex[j] = -1;
  }
  for ( int j=0; j<numberOfSimTracks; ++j ) {
    for ( int i=0; i<5; ++i ) {
       simPars[i][j] = 0;
    }
    recIndex[j] = -1;
  }

  numberOfRecTracks=0;
  numberOfSimTracks=0;
}