#include <TrackerValidationVariables.h>

Classes
struct	AVHitStruct

struct	AVTrackStruct

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

void	fillHitQuantities (reco::Track const &track, std::vector< AVHitStruct > &v_avhitout)

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

void	fillTrackQuantities (const edm::Event &event, const edm::EventSetup &eventSetup, std::function< bool(const reco::Track &)> trackFilter, std::vector< AVTrackStruct > &v_avtrackout)

	TrackerValidationVariables ()

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

	~TrackerValidationVariables ()

Private Attributes
edm::EDGetTokenT< std::vector < reco::Track > >	tracksToken_

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

Detailed Description

Definition at line 20 of file TrackerValidationVariables.h.

Constructor & Destructor Documentation

TrackerValidationVariables::TrackerValidationVariables ( )

Definition at line 48 of file TrackerValidationVariables.cc.

 {
 
 }

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

Definition at line 53 of file TrackerValidationVariables.cc.

References edm::ParameterSet::getParameter(), HLT_FULL_cff::InputTag, AlCaHLTBitMon_QueryRunRegistry::string, tracksToken_, and trajCollectionToken_.

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

TrackerValidationVariables::~TrackerValidationVariables ( )

Definition at line 60 of file TrackerValidationVariables.cc.

 {
 
 }

Member Function Documentation

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

Definition at line 194 of file TrackerValidationVariables.cc.

Referenced by 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() );
     if(!tsos.isValid()) continue;
     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.barePhi(),gPModule.barePhi()) >= 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.barePhi(),gPModule.barePhi()) >= 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.barePhi(),gPModule.barePhi()) >= 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	(	reco::Track const &	track,
		std::vector< AVHitStruct > &	v_avhitout
	)

Definition at line 65 of file TrackerValidationVariables.cc.

                                                                                                                {
 
     auto const & trajParams = track.extra()->trajParams();
     auto const & residuals = track.extra()->residuals();
 
     assert(trajParams.size()==track.recHitsSize());
     auto hb = track.recHitsBegin();
     for(unsigned int h=0;h<track.recHitsSize();h++){
       auto hit = *(hb+h);
       if(!hit->isValid()) continue;
 
       AVHitStruct hitStruct;
       const DetId & hit_detId = hit->geographicalId();
       auto IntRawDetID = hit_detId.rawId();
       auto IntSubDetID = hit_detId.subdetId();
 
       if(IntSubDetID == 0) continue;
 
       auto lPTrk = trajParams[h].position();   // update state
       auto lVTrk = trajParams[h].direction();
 
       auto gtrkdirup= hit->surface()->toGlobal(lVTrk);
 
       hitStruct.rawDetId = IntRawDetID;
       hitStruct.phi = gtrkdirup.phi();  // direction, not position
       hitStruct.eta = gtrkdirup.eta();  // same
 
 
       hitStruct.localAlpha = std::atan2(lVTrk.x(), lVTrk.z()); // wrt. normal tg(alpha)=x/z
       hitStruct.localBeta  = std::atan2(lVTrk.y(), lVTrk.z()); // wrt. normal tg(beta)= y/z
 
       hitStruct.resX = residuals.residualX(h);
       hitStruct.resY = residuals.residualY(h);
       hitStruct.resErrX = hitStruct.resX/residuals.pullX(h);  // for backward compatibility....
       hitStruct.resErrY = hitStruct.resY/residuals.pullY(h);
 
       // 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?
         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 == PixelSubdetector::PixelEndcap ||
           IntSubDetID == StripSubdetector::TIB ||
           IntSubDetID == StripSubdetector::TOB) {
 
           if (IntSubDetID == PixelSubdetector::PixelEndcap) {
             uOrientation = gUDirection.perp() - gPModule.perp() >= 0 ? +1.F : -1.F;
             vOrientation = deltaPhi(gVDirection.barePhi(),gPModule.barePhi()) >= 0. ? +1.F : -1.F;
           } else {
             uOrientation = deltaPhi(gUDirection.barePhi(),gPModule.barePhi()) >= 0. ? +1.F : -1.F;
             vOrientation = gVDirection.z() - gPModule.z() >= 0 ? +1.F : -1.F;  
           }
 
           resXTopol = hitStruct.resX;
           resXatTrkYTopol = hitStruct.resX;
           resYTopol = hitStruct.resY;
           resXprimeErr = hitStruct.resErrX;
           resYprimeErr = hitStruct.resErrY;
 
           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) {
         // not possible to compute precisely as with Trajectory
       } 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;
     
     
     v_avhitout.push_back(hitStruct);
 
    }
 }

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

Definition at line 437 of file TrackerValidationVariables.cc.

Referenced by SiPixelPhase1TrackResiduals::analyze(), and TrackerOfflineValidation::analyze().

 {
   fillTrackQuantities(event, 
                       eventSetup,
                       [](const reco::Track&) -> bool { return true; },
                       v_avtrackout);
 }

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

Definition at line 448 of file TrackerValidationVariables.cc.

 {
   edm::ESHandle<MagneticField> magneticField;
   eventSetup.get<IdealMagneticFieldRecord>().get(magneticField);
 
   edm::Handle<reco::TrackCollection> tracksH;
   event.getByToken(tracksToken_, tracksH);
   if(!tracksH.isValid()) return;
   auto const & tracks = *tracksH;
   auto ntrk = tracks.size();
   LogDebug("TrackerValidationVariables") << "Track collection size " << ntrk;
   
 
   edm::Handle<std::vector<Trajectory>> trajsH;
   event.getByToken(trajCollectionToken_, trajsH);
   bool yesTraj = trajsH.isValid();
   std::vector<Trajectory> const * trajs = nullptr;
   if (yesTraj) trajs = &(*trajsH);
   if (yesTraj) assert (trajs->size()==tracks.size());
 
   Trajectory const * trajectory =  nullptr;
   for (unsigned int i=0; i<ntrk; ++i) {
     auto const & track = tracks[i];  
     if (yesTraj) trajectory = &(*trajs)[i];
     
     if (!trackFilter(track)) continue;
     
     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();
     if (trajectory)
       fillHitQuantities(trajectory, trackStruct.hits);
     else fillHitQuantities(track, trackStruct.hits);
     
     v_avtrackout.push_back(trackStruct);
   }
 }

Member Data Documentation

edm::EDGetTokenT<std::vector<reco::Track> > TrackerValidationVariables::tracksToken_

private

Definition at line 94 of file TrackerValidationVariables.h.

Referenced by fillTrackQuantities(), and TrackerValidationVariables().

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

private

Definition at line 93 of file TrackerValidationVariables.h.

Referenced by fillTrackQuantities(), and TrackerValidationVariables().

Classes

Public Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation