TrackerValidationVariables.cc

Go to the documentation of this file.

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Framework/interface/ConsumesCollector.h"
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Utilities/interface/InputTag.h"

#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
#include "Geometry/CommonDetUnit/interface/GeomDetType.h"
#include "Geometry/CommonTopologies/interface/RadialStripTopology.h"

#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
#include "MagneticField/Engine/interface/MagneticField.h"

#include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
#include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
#include "TrackingTools/PatternTools/interface/Trajectory.h"
#include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
#include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"

#include "CondFormats/Alignment/interface/Definitions.h"

#include "DataFormats/Math/interface/deltaPhi.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
#include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
#include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
#include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementPoint.h"
#include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementError.h"
#include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
#include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"

#include "Alignment/TrackerAlignment/interface/TrackerAlignableId.h"

#include "Alignment/OfflineValidation/interface/TrackerValidationVariables.h"

#include "TMath.h"

#include <string>

TrackerValidationVariables::TrackerValidationVariables()
{

}

TrackerValidationVariables::TrackerValidationVariables(const edm::ParameterSet& config,
                                                       edm::ConsumesCollector && iC)
{
  trajCollectionToken_ = iC.consumes<std::vector<Trajectory> >(edm::InputTag(config.getParameter<std::string>("trajectoryInput")));
  trajTracksToken_ = iC.consumes<TrajTrackAssociationCollection>(config.getParameter<edm::InputTag>("Tracks"));
}

TrackerValidationVariables::~TrackerValidationVariables()
{

}

void
TrackerValidationVariables::fillHitQuantities(const Trajectory* trajectory, std::vector<AVHitStruct> & v_avhitout)
{
  TrajectoryStateCombiner tsoscomb;

  const std::vector<TrajectoryMeasurement> &tmColl = trajectory->measurements();
  for(std::vector<TrajectoryMeasurement>::const_iterator itTraj = tmColl.begin();
      itTraj != tmColl.end();
      ++itTraj) {
    
    if (!itTraj->updatedState().isValid()) continue;
    
    TrajectoryStateOnSurface tsos = tsoscomb( itTraj->forwardPredictedState(), itTraj->backwardPredictedState() );
    TransientTrackingRecHit::ConstRecHitPointer hit = itTraj->recHit();
    
    if(!hit->isValid() || hit->geographicalId().det() != DetId::Tracker) continue;
    
    AVHitStruct hitStruct;
    const DetId& hit_detId = hit->geographicalId();
    unsigned int IntRawDetID = (hit_detId.rawId()); 
    unsigned int IntSubDetID = (hit_detId.subdetId());
    
    if(IntSubDetID == 0) continue;
    
    //first calculate residuals in cartesian coordinates in the local module coordinate system
    
    LocalPoint lPHit = hit->localPosition();
    LocalPoint lPTrk = tsos.localPosition();
    LocalVector lVTrk = tsos.localDirection();

    hitStruct.localAlpha = atan2(lVTrk.x(), lVTrk.z()); // wrt. normal tg(alpha)=x/z
    hitStruct.localBeta  = atan2(lVTrk.y(), lVTrk.z()); // wrt. normal tg(beta)= y/z

    LocalError errHit = hit->localPositionError();
    // no need to add  APE to hitError anymore
    // AlgebraicROOTObject<2>::SymMatrix mat = asSMatrix<2>(hit->parametersError());
    // LocalError errHit = LocalError( mat(0,0),mat(0,1),mat(1,1) );
    LocalError errTrk = tsos.localError().positionError();
    
    //check for negative error values: track error can have negative value, if matrix inversion fails (very rare case)
    //hit error should always give positive values
    if(errHit.xx()<0. || errHit.yy()<0. || errTrk.xx()<0. || errTrk.yy()<0.){
      edm::LogError("TrackerValidationVariables") << "@SUB=TrackerValidationVariables::fillHitQuantities"
                          << "One of the squared error methods gives negative result"
                          << "\n\terrHit.xx()\terrHit.yy()\terrTrk.xx()\terrTrk.yy()"
                          << "\n\t" << errHit.xx()
                          << "\t" << errHit.yy()
                          << "\t" << errTrk.xx()
                          << "\t" << errTrk.yy();
      continue;
    }
    
    align::LocalVector res = lPTrk - lPHit;
    
    float resXErr = std::sqrt( errHit.xx() + errTrk.xx() );
    float resYErr = std::sqrt( errHit.yy() + errTrk.yy() );
    
    hitStruct.resX = res.x();
    hitStruct.resY = res.y();
    hitStruct.resErrX = resXErr;
    hitStruct.resErrY = resYErr;

    // hitStruct.localX = lPhit.x();
    // hitStruct.localY = lPhit.y();
    // EM: use predictions for local coordinates
    hitStruct.localX = lPTrk.x();
    hitStruct.localY = lPTrk.y();

    // now calculate residuals taking global orientation of modules and radial topology in TID/TEC into account
    float resXprime(999.F), resYprime(999.F);
    float resXatTrkY(999.F);
    float resXprimeErr(999.F), resYprimeErr(999.F);
    
    if(hit->detUnit()){ // is it a single physical module?
      const GeomDetUnit& detUnit = *(hit->detUnit());
      float uOrientation(-999.F), vOrientation(-999.F);
      float resXTopol(999.F), resYTopol(999.F);
      float resXatTrkYTopol(999.F);       

      const Surface& surface = hit->detUnit()->surface();
      const BoundPlane& boundplane = hit->detUnit()->surface();
      const Bounds& bound = boundplane.bounds();
      
      float length = 0;
      float width = 0;

      LocalPoint lPModule(0.,0.,0.), lUDirection(1.,0.,0.), lVDirection(0.,1.,0.);
      GlobalPoint gPModule    = surface.toGlobal(lPModule),
              gUDirection = surface.toGlobal(lUDirection),
              gVDirection = surface.toGlobal(lVDirection);
      
      if (IntSubDetID == PixelSubdetector::PixelBarrel ||
      IntSubDetID == StripSubdetector::TIB || 
      IntSubDetID == StripSubdetector::TOB) {
    
    uOrientation = deltaPhi(gUDirection.phi(),gPModule.phi()) >= 0. ? +1.F : -1.F;
    vOrientation = gVDirection.z() - gPModule.z() >= 0 ? +1.F : -1.F;
    resXTopol = res.x();
    resXatTrkYTopol = res.x();
    resYTopol = res.y();
    resXprimeErr = resXErr;
    resYprimeErr = resYErr;

    const RectangularPlaneBounds *rectangularBound = dynamic_cast<const RectangularPlaneBounds*>(&bound);
    if (rectangularBound!=NULL) {
      hitStruct.inside = rectangularBound->inside(lPTrk);
      length = rectangularBound->length();
      width = rectangularBound->width();
      hitStruct.localXnorm = 2*hitStruct.localX/width;
      hitStruct.localYnorm = 2*hitStruct.localY/length;
    } else {
      throw cms::Exception("Geometry Error")
        << "[TrackerValidationVariables] Cannot cast bounds to RectangularPlaneBounds as expected for TPB, TIB and TOB";
    }

      } else if (IntSubDetID == PixelSubdetector::PixelEndcap) {
    
    uOrientation = gUDirection.perp() - gPModule.perp() >= 0 ? +1.F : -1.F;
    vOrientation = deltaPhi(gVDirection.phi(),gPModule.phi()) >= 0. ? +1.F : -1.F;
    resXTopol = res.x();
    resXatTrkYTopol = res.x();
    resYTopol = res.y();
    resXprimeErr = resXErr;
    resYprimeErr = resYErr;
    
    const RectangularPlaneBounds *rectangularBound = dynamic_cast<const RectangularPlaneBounds*>(&bound);
    if (rectangularBound!=NULL) {
      hitStruct.inside = rectangularBound->inside(lPTrk);
      length = rectangularBound->length();
      width = rectangularBound->width();
      hitStruct.localXnorm = 2*hitStruct.localX/width;
      hitStruct.localYnorm = 2*hitStruct.localY/length;
    } else {
      throw cms::Exception("Geometry Error")
        << "[TrackerValidationVariables] Cannot cast bounds to RectangularPlaneBounds as expected for TPE";
    }

      } else if (IntSubDetID == StripSubdetector::TID ||
         IntSubDetID == StripSubdetector::TEC) {
    
    uOrientation = deltaPhi(gUDirection.phi(),gPModule.phi()) >= 0. ? +1.F : -1.F;
    vOrientation = gVDirection.perp() - gPModule.perp() >= 0. ? +1.F : -1.F;
    
    if (!dynamic_cast<const RadialStripTopology*>(&detUnit.type().topology()))continue;
    const RadialStripTopology& topol = dynamic_cast<const RadialStripTopology&>(detUnit.type().topology());
    
    MeasurementPoint measHitPos = topol.measurementPosition(lPHit);
    MeasurementPoint measTrkPos = topol.measurementPosition(lPTrk);
    
    MeasurementError measHitErr = topol.measurementError(lPHit,errHit);
    MeasurementError measTrkErr = topol.measurementError(lPTrk,errTrk);
    
    if (measHitErr.uu()<0. ||
        measHitErr.vv()<0. ||
        measTrkErr.uu()<0. ||
        measTrkErr.vv()<0.){
      edm::LogError("TrackerValidationVariables") << "@SUB=TrackerValidationVariables::fillHitQuantities"
                              << "One of the squared error methods gives negative result"
                              << "\n\tmeasHitErr.uu()\tmeasHitErr.vv()\tmeasTrkErr.uu()\tmeasTrkErr.vv()"
                              << "\n\t" << measHitErr.uu()
                              << "\t" << measHitErr.vv()
                              << "\t" << measTrkErr.uu()
                              << "\t" << measTrkErr.vv();
      continue;
    }
      
    float localStripLengthHit = topol.localStripLength(lPHit);
    float localStripLengthTrk = topol.localStripLength(lPTrk);
    float phiHit = topol.stripAngle(measHitPos.x());
    float phiTrk = topol.stripAngle(measTrkPos.x());
    float r_0 = topol.originToIntersection();
    
    resXTopol = (phiTrk-phiHit)*r_0;
//      resXTopol = (tan(phiTrk)-tan(phiHit))*r_0;
        
    LocalPoint LocalHitPosCor= topol.localPosition(MeasurementPoint(measHitPos.x(), measTrkPos.y()));                       
       resXatTrkYTopol = lPTrk.x() - LocalHitPosCor.x();  

    
    //resYTopol = measTrkPos.y()*localStripLengthTrk - measHitPos.y()*localStripLengthHit;
    float cosPhiHit(cos(phiHit)), cosPhiTrk(cos(phiTrk)),
          sinPhiHit(sin(phiHit)), sinPhiTrk(sin(phiTrk));
    float l_0 = r_0 - topol.detHeight()/2;
    resYTopol = measTrkPos.y()*localStripLengthTrk - measHitPos.y()*localStripLengthHit + l_0*(1/cosPhiTrk - 1/cosPhiHit);
    
    resXprimeErr = std::sqrt(measHitErr.uu()+measTrkErr.uu())*topol.angularWidth()*r_0;
    //resYprimeErr = std::sqrt(measHitErr.vv()*localStripLengthHit*localStripLengthHit + measTrkErr.vv()*localStripLengthTrk*localStripLengthTrk);
    float helpSummand = l_0*l_0*topol.angularWidth()*topol.angularWidth()*(sinPhiHit*sinPhiHit/pow(cosPhiHit,4)*measHitErr.uu()
                                           + sinPhiTrk*sinPhiTrk/pow(cosPhiTrk,4)*measTrkErr.uu() );
    resYprimeErr = std::sqrt(measHitErr.vv()*localStripLengthHit*localStripLengthHit
                 + measTrkErr.vv()*localStripLengthTrk*localStripLengthTrk + helpSummand );  


    const TrapezoidalPlaneBounds *trapezoidalBound = dynamic_cast < const TrapezoidalPlaneBounds * >(& bound);
    if (trapezoidalBound!=NULL) {
      hitStruct.inside = trapezoidalBound->inside(lPTrk);
      length = trapezoidalBound->length();
      width  = trapezoidalBound->width();
      //float widthAtHalfLength = trapezoidalBound->widthAtHalfLength();

      //    int yAxisOrientation=trapezoidalBound->yAxisOrientation(); 
      // for trapezoidal shape modules, scale with as function of local y coordinate 
      //    float widthAtlocalY=width-(1-yAxisOrientation*2*lPTrk.y()/length)*(width-widthAtHalfLength); 
      //    hitStruct.localXnorm = 2*hitStruct.localX/widthAtlocalY;  
      hitStruct.localXnorm = 2*hitStruct.localX/width; 
      hitStruct.localYnorm = 2*hitStruct.localY/length;
    } else {
      throw cms::Exception("Geometry Error")
        << "[TrackerValidationVariables] Cannot cast bounds to TrapezoidalPlaneBounds as expected for TID and TEC";
    }

      } else {
    edm::LogWarning("TrackerValidationVariables") << "@SUB=TrackerValidationVariables::fillHitQuantities" 
                              << "No valid tracker subdetector " << IntSubDetID;
    continue;
      }  
      
      resXprime = resXTopol*uOrientation;
      resXatTrkY = resXatTrkYTopol;
      resYprime = resYTopol*vOrientation;
      
    } else { // not a detUnit, so must be a virtual 2D-Module
      // FIXME: at present only for det units residuals are calculated and filled in the hitStruct
      // But in principle this method should also be useable for the gluedDets (2D modules in TIB, TID, TOB, TEC)
      // In this case, only orientation should be taken into account for primeResiduals, but not the radial topology
      // At present, default values (999.F) are given out
    }
    
    hitStruct.resXprime = resXprime;
    hitStruct.resXatTrkY = resXatTrkY;
    hitStruct.resYprime = resYprime;
    hitStruct.resXprimeErr = resXprimeErr;
    hitStruct.resYprimeErr = resYprimeErr;
    
    hitStruct.rawDetId = IntRawDetID;
    hitStruct.phi = tsos.globalDirection().phi();
    hitStruct.eta = tsos.globalDirection().eta();
    
    v_avhitout.push_back(hitStruct);
  } 
}

void
TrackerValidationVariables::fillTrackQuantities(const edm::Event& event,
                                                const edm::EventSetup& eventSetup,
                                                std::vector<AVTrackStruct> & v_avtrackout)
{
  edm::ESHandle<MagneticField> magneticField;
  eventSetup.get<IdealMagneticFieldRecord>().get(magneticField);

  edm::Handle<TrajTrackAssociationCollection> TrajTracksMap;
  event.getByToken(trajTracksToken_, TrajTracksMap);
  LogDebug("TrackerValidationVariables") << "TrajTrack collection size " << TrajTracksMap->size();
  
  const Trajectory* trajectory;
  const reco::Track* track;
  
  for ( TrajTrackAssociationCollection::const_iterator iPair = TrajTracksMap->begin();
    iPair != TrajTracksMap->end();
    ++iPair) {
    
    trajectory = &(*(*iPair).key);
    track = &(*(*iPair).val);
    
    AVTrackStruct trackStruct;
    
    trackStruct.p = track->p();
    trackStruct.pt = track->pt();
    trackStruct.ptError = track->ptError();
    trackStruct.px = track->px();
    trackStruct.py = track->py();
    trackStruct.pz = track->pz();
    trackStruct.eta = track->eta();
    trackStruct.phi = track->phi();
    trackStruct.chi2 = track->chi2();
    trackStruct.chi2Prob= TMath::Prob(track->chi2(),track->ndof());
    trackStruct.normchi2 = track->normalizedChi2();
    GlobalPoint gPoint(track->vx(), track->vy(), track->vz());
    double theLocalMagFieldInInverseGeV = magneticField->inInverseGeV(gPoint).z();
    trackStruct.kappa = -track->charge()*theLocalMagFieldInInverseGeV/track->pt();
    trackStruct.charge = track->charge();
    trackStruct.d0 = track->d0();
    trackStruct.dz = track->dz();
    trackStruct.numberOfValidHits = track->numberOfValidHits();
    trackStruct.numberOfLostHits = track->numberOfLostHits();
    
    fillHitQuantities(trajectory, trackStruct.hits);
    
    v_avtrackout.push_back(trackStruct);
  }
}

void
TrackerValidationVariables::fillHitQuantities(const edm::Event& event, std::vector<AVHitStruct> & v_avhitout)
{
  edm::Handle<std::vector<Trajectory> > trajCollectionHandle;
  event.getByToken(trajCollectionToken_, trajCollectionHandle);
  
  LogDebug("TrackerValidationVariables") << "trajColl->size(): " << trajCollectionHandle->size() ;

  for (std::vector<Trajectory>::const_iterator it = trajCollectionHandle->begin(), itEnd = trajCollectionHandle->end(); 
       it!=itEnd;
       ++it) {
    
    fillHitQuantities(&(*it), v_avhitout);
  }
}