TrackerValidationVariables Class Reference

#include <TrackerValidationVariables.h>

Classes
struct	AVHitStruct
struct	AVTrackStruct

Public Member Functions
void	fillHitQuantities (const Trajectory *trajectory, std::vector< AVHitStruct > &v_avhitout)
void	fillHitQuantities (const edm::Event &, std::vector< AVHitStruct > &v_avhitout)
void	fillTrackQuantities (const edm::Event &, const edm::EventSetup &, std::vector< AVTrackStruct > &v_avtrackout)
	TrackerValidationVariables ()
	TrackerValidationVariables (const edm::ParameterSet &config, edm::ConsumesCollector &&iC)
	~TrackerValidationVariables ()

Private Attributes
edm::EDGetTokenT< std::vector < Trajectory > >	trajCollectionToken_
edm::EDGetTokenT < TrajTrackAssociationCollection >	trajTracksToken_

Detailed Description

Definition at line 20 of file TrackerValidationVariables.h.

Constructor & Destructor Documentation

TrackerValidationVariables::TrackerValidationVariables

(

)

Definition at line 49 of file TrackerValidationVariables.cc.

{

}

TrackerValidationVariables::TrackerValidationVariables	(	const edm::ParameterSet &	config,
		edm::ConsumesCollector &&	iC
	)

Definition at line 54 of file TrackerValidationVariables.cc.

References edm::ParameterSet::getParameter(), HLT_25ns14e33_v1_cff::InputTag, AlCaHLTBitMon_QueryRunRegistry::string, trajCollectionToken_, and trajTracksToken_.

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

TrackerValidationVariables::~TrackerValidationVariables

(

)

Definition at line 61 of file TrackerValidationVariables.cc.

{

}

Member Function Documentation

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

Definition at line 67 of file TrackerValidationVariables.cc.

References RadialStripTopology::angularWidth(), funct::cos(), SiPixelRawToDigiRegional_cfi::deltaPhi, RadialStripTopology::detHeight(), TrackerValidationVariables::AVHitStruct::eta, PV3DBase< T, PVType, FrameType >::eta(), Exception, F(), TrajectoryStateOnSurface::globalDirection(), RectangularPlaneBounds::inside(), TrackerValidationVariables::AVHitStruct::inside, RectangularPlaneBounds::length(), TrackerValidationVariables::AVHitStruct::localAlpha, TrackerValidationVariables::AVHitStruct::localBeta, TrajectoryStateOnSurface::localDirection(), TrajectoryStateOnSurface::localError(), RadialStripTopology::localPosition(), TrajectoryStateOnSurface::localPosition(), RadialStripTopology::localStripLength(), TrackerValidationVariables::AVHitStruct::localX, TrackerValidationVariables::AVHitStruct::localXnorm, TrackerValidationVariables::AVHitStruct::localY, TrackerValidationVariables::AVHitStruct::localYnorm, RadialStripTopology::measurementError(), RadialStripTopology::measurementPosition(), Trajectory::measurements(), NULL, RadialStripTopology::originToIntersection(), PV3DBase< T, PVType, FrameType >::perp(), TrackerValidationVariables::AVHitStruct::phi, PV3DBase< T, PVType, FrameType >::phi(), PixelSubdetector::PixelBarrel, PixelSubdetector::PixelEndcap, LocalTrajectoryError::positionError(), funct::pow(), TrackerValidationVariables::AVHitStruct::rawDetId, DetId::rawId(), TrackerValidationVariables::AVHitStruct::resErrX, TrackerValidationVariables::AVHitStruct::resErrY, TrackerValidationVariables::AVHitStruct::resX, TrackerValidationVariables::AVHitStruct::resXatTrkY, TrackerValidationVariables::AVHitStruct::resXprime, TrackerValidationVariables::AVHitStruct::resXprimeErr, TrackerValidationVariables::AVHitStruct::resY, TrackerValidationVariables::AVHitStruct::resYprime, TrackerValidationVariables::AVHitStruct::resYprimeErr, funct::sin(), mathSSE::sqrt(), RadialStripTopology::stripAngle(), DetId::subdetId(), StripSubdetector::TEC, StripSubdetector::TIB, StripSubdetector::TID, StripSubdetector::TOB, Surface::toGlobal(), GeomDetType::topology(), DetId::Tracker, GeomDet::type(), MeasurementError::uu(), MeasurementError::vv(), RectangularPlaneBounds::width(), create_public_lumi_plots::width, PV2DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::x(), LocalError::xx(), PV2DBase< T, PVType, FrameType >::y(), PV3DBase< T, PVType, FrameType >::y(), LocalError::yy(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by MonitorTrackResiduals::analyze(), fillHitQuantities(), and fillTrackQuantities().

{
  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::fillHitQuantities	(	const edm::Event &	event,
		std::vector< AVHitStruct > &	v_avhitout
	)

Definition at line 359 of file TrackerValidationVariables.cc.

References fillHitQuantities(), LogDebug, and trajCollectionToken_.

{
  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);
  }
}

void TrackerValidationVariables::fillTrackQuantities	(	const edm::Event &	event,
		const edm::EventSetup &	eventSetup,
		std::vector< AVTrackStruct > &	v_avtrackout
	)

Definition at line 309 of file TrackerValidationVariables.cc.

References TrackerValidationVariables::AVTrackStruct::charge, reco::TrackBase::charge(), TrackerValidationVariables::AVTrackStruct::chi2, reco::TrackBase::chi2(), TrackerValidationVariables::AVTrackStruct::chi2Prob, edm::AssociationMap< edm::OneToOne< std::vector< Trajectory >, reco::TrackCollection, unsigned short > >::const_iterator, TrackerValidationVariables::AVTrackStruct::d0, reco::TrackBase::d0(), TrackerValidationVariables::AVTrackStruct::dz, reco::TrackBase::dz(), TrackerValidationVariables::AVTrackStruct::eta, reco::TrackBase::eta(), fillHitQuantities(), edm::EventSetup::get(), TrackerValidationVariables::AVTrackStruct::hits, TrackerValidationVariables::AVTrackStruct::kappa, LogDebug, HLT_25ns14e33_v1_cff::magneticField, reco::TrackBase::ndof(), reco::TrackBase::normalizedChi2(), TrackerValidationVariables::AVTrackStruct::normchi2, TrackerValidationVariables::AVTrackStruct::numberOfLostHits, reco::TrackBase::numberOfLostHits(), TrackerValidationVariables::AVTrackStruct::numberOfValidHits, reco::TrackBase::numberOfValidHits(), TrackerValidationVariables::AVTrackStruct::p, reco::TrackBase::p(), TrackerValidationVariables::AVTrackStruct::phi, reco::TrackBase::phi(), TrackerValidationVariables::AVTrackStruct::pt, reco::TrackBase::pt(), TrackerValidationVariables::AVTrackStruct::ptError, reco::TrackBase::ptError(), TrackerValidationVariables::AVTrackStruct::px, reco::TrackBase::px(), TrackerValidationVariables::AVTrackStruct::py, reco::TrackBase::py(), TrackerValidationVariables::AVTrackStruct::pz, reco::TrackBase::pz(), trajTracksToken_, reco::TrackBase::vx(), reco::TrackBase::vy(), and reco::TrackBase::vz().

Referenced by TrackerOfflineValidation::analyze().

{
  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);
  }
}

Member Data Documentation

edm::EDGetTokenT<std::vector<Trajectory> > TrackerValidationVariables::trajCollectionToken_

private

Definition at line 90 of file TrackerValidationVariables.h.

Referenced by fillHitQuantities(), and TrackerValidationVariables().

edm::EDGetTokenT<TrajTrackAssociationCollection> TrackerValidationVariables::trajTracksToken_

private

Definition at line 91 of file TrackerValidationVariables.h.

Referenced by fillTrackQuantities(), and TrackerValidationVariables().