#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 &, std::vector< AVTrackStruct > &v_avtrackout)

	TrackerValidationVariables ()

	TrackerValidationVariables (const edm::EventSetup &, const edm::ParameterSet &)

	~TrackerValidationVariables ()

Private Attributes
const edm::ParameterSet	conf_

edm::ESHandle< MagneticField >	magneticField_

edm::ESHandle< TrackerGeometry >	tkGeom_

Detailed Description

Definition at line 16 of file TrackerValidationVariables.h.

Constructor & Destructor Documentation

TrackerValidationVariables::TrackerValidationVariables ( )

Definition at line 45 of file TrackerValidationVariables.cc.

 {
 
 }

TrackerValidationVariables::TrackerValidationVariables	(	const edm::EventSetup &	es,
		const edm::ParameterSet &	iSetup
	)

Definition at line 50 of file TrackerValidationVariables.cc.

References edm::EventSetup::get(), magneticField_, and tkGeom_.

   : conf_(iSetup)
 {
   es.get<TrackerDigiGeometryRecord>().get(tkGeom_);
   es.get<IdealMagneticFieldRecord>().get(magneticField_);
 }

TrackerValidationVariables::~TrackerValidationVariables ( )

Definition at line 57 of file TrackerValidationVariables.cc.

 {
 
 }

Member Function Documentation

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

Definition at line 63 of file TrackerValidationVariables.cc.

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();
     // adding APE to hitError
     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 351 of file TrackerValidationVariables.cc.

References conf_, fillHitQuantities(), edm::ParameterSet::getParameter(), and LogDebug.

 {
   edm::Handle<std::vector<Trajectory> > trajCollectionHandle;
   event.getByLabel(conf_.getParameter<std::string>("trajectoryInput"), 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,
		std::vector< AVTrackStruct > &	v_avtrackout
	)

Definition at line 305 of file TrackerValidationVariables.cc.

Referenced by TrackerOfflineValidation::analyze().

 {
   edm::InputTag TrjTrackTag = conf_.getParameter<edm::InputTag>("Tracks");
   edm::Handle<TrajTrackAssociationCollection> TrajTracksMap;
   event.getByLabel(TrjTrackTag, 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);
   }
 }