OfflineValidation/src/TrackerValidationVariables.cc

0001 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0002 #include "FWCore/Framework/interface/ConsumesCollector.h"
0003 #include "FWCore/Framework/interface/Frameworkfwd.h"
0004 #include "FWCore/Framework/interface/MakerMacros.h"
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/Framework/interface/EventSetup.h"
0007 #include "FWCore/Utilities/interface/InputTag.h"
0008
0009 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0010 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0011 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0012 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0013 #include "Geometry/CommonTopologies/interface/RadialStripTopology.h"
0014
0015 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
0016 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0017 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0018 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0019 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0020 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0021 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0022
0023 #include "CondFormats/Alignment/interface/Definitions.h"
0024
0025 #include "DataFormats/Math/interface/deltaPhi.h"
0026 #include "DataFormats/TrackReco/interface/Track.h"
0027 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0028 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0029 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0030 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0031 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementPoint.h"
0032 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementError.h"
0033 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
0034 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
0035
0036 #include "Alignment/TrackerAlignment/interface/TrackerAlignableId.h"
0037
0038 #include "Alignment/OfflineValidation/interface/TrackerValidationVariables.h"
0039
0040 #include "TMath.h"
0041
0042 #include <string>
0043
0044 TrackerValidationVariables::TrackerValidationVariables(const edm::ParameterSet& config, edm::ConsumesCollector&& iC)
0045     : magneticFieldToken_{iC.esConsumes<MagneticField, IdealMagneticFieldRecord>()} {
0046   trajCollectionToken_ =
0047       iC.consumes<std::vector<Trajectory>>(edm::InputTag(config.getParameter<std::string>("trajectoryInput")));
0048   tracksToken_ = iC.consumes<reco::TrackCollection>(config.getParameter<edm::InputTag>("Tracks"));
0049 }
0050
0051 TrackerValidationVariables::~TrackerValidationVariables() {}
0052
0053 void TrackerValidationVariables::fillHitQuantities(reco::Track const& track, std::vector<AVHitStruct>& v_avhitout) {
0054   auto const& trajParams = track.extra()->trajParams();
0055   auto const& residuals = track.extra()->residuals();
0056
0057   assert(trajParams.size() == track.recHitsSize());
0058   auto hb = track.recHitsBegin();
0059   for (unsigned int h = 0; h < track.recHitsSize(); h++) {
0060     auto hit = *(hb + h);
0061     if (!hit->isValid())
0062       continue;
0063
0064     AVHitStruct hitStruct;
0065     const DetId& hit_detId = hit->geographicalId();
0066     auto IntRawDetID = hit_detId.rawId();
0067     auto IntSubDetID = hit_detId.subdetId();
0068
0069     if (IntSubDetID == 0)
0070       continue;
0071
0072     if (IntSubDetID == PixelSubdetector::PixelBarrel || IntSubDetID == PixelSubdetector::PixelEndcap) {
0073       const SiPixelRecHit* prechit = dynamic_cast<const SiPixelRecHit*>(
0074           hit);  //to be used to get the associated cluster and the cluster probability
0075       if (prechit->isOnEdge())
0076         hitStruct.isOnEdgePixel = true;
0077       if (prechit->hasBadPixels())
0078         hitStruct.isOtherBadPixel = true;
0079     }
0080
0081     auto lPTrk = trajParams[h].position();  // update state
0082     auto lVTrk = trajParams[h].direction();
0083
0084     auto gtrkdirup = hit->surface()->toGlobal(lVTrk);
0085
0086     hitStruct.rawDetId = IntRawDetID;
0087     hitStruct.phi = gtrkdirup.phi();  // direction, not position
0088     hitStruct.eta = gtrkdirup.eta();  // same
0089
0090     hitStruct.localAlpha = std::atan2(lVTrk.x(), lVTrk.z());  // wrt. normal tg(alpha)=x/z
0091     hitStruct.localBeta = std::atan2(lVTrk.y(), lVTrk.z());   // wrt. normal tg(beta)= y/z
0092
0093     hitStruct.resX = residuals.residualX(h);
0094     hitStruct.resY = residuals.residualY(h);
0095     hitStruct.resErrX = hitStruct.resX / residuals.pullX(h);  // for backward compatibility....
0096     hitStruct.resErrY = hitStruct.resY / residuals.pullY(h);
0097
0098     // hitStruct.localX = lPhit.x();
0099     // hitStruct.localY = lPhit.y();
0100     // EM: use predictions for local coordinates
0101     hitStruct.localX = lPTrk.x();
0102     hitStruct.localY = lPTrk.y();
0103
0104     // now calculate residuals taking global orientation of modules and radial topology in TID/TEC into account
0105     float resXprime(999.F), resYprime(999.F);
0106     float resXatTrkY(999.F);
0107     float resXprimeErr(999.F), resYprimeErr(999.F);
0108
0109     if (hit->detUnit()) {  // is it a single physical module?
0110       float uOrientation(-999.F), vOrientation(-999.F);
0111       float resXTopol(999.F), resYTopol(999.F);
0112       float resXatTrkYTopol(999.F);
0113
0114       const Surface& surface = hit->detUnit()->surface();
0115       const BoundPlane& boundplane = hit->detUnit()->surface();
0116       const Bounds& bound = boundplane.bounds();
0117
0118       float length = 0;
0119       float width = 0;
0120
0121       LocalPoint lPModule(0., 0., 0.), lUDirection(1., 0., 0.), lVDirection(0., 1., 0.);
0122       GlobalPoint gPModule = surface.toGlobal(lPModule), gUDirection = surface.toGlobal(lUDirection),
0123                   gVDirection = surface.toGlobal(lVDirection);
0124
0125       if (IntSubDetID == PixelSubdetector::PixelBarrel || IntSubDetID == PixelSubdetector::PixelEndcap ||
0126           IntSubDetID == StripSubdetector::TIB || IntSubDetID == StripSubdetector::TOB) {
0127         if (IntSubDetID == PixelSubdetector::PixelEndcap) {
0128           uOrientation = gUDirection.perp() - gPModule.perp() >= 0 ? +1.F : -1.F;
0129           vOrientation = deltaPhi(gVDirection.barePhi(), gPModule.barePhi()) >= 0. ? +1.F : -1.F;
0130         } else {
0131           uOrientation = deltaPhi(gUDirection.barePhi(), gPModule.barePhi()) >= 0. ? +1.F : -1.F;
0132           vOrientation = gVDirection.z() - gPModule.z() >= 0 ? +1.F : -1.F;
0133         }
0134
0135         resXTopol = hitStruct.resX;
0136         resXatTrkYTopol = hitStruct.resX;
0137         resYTopol = hitStruct.resY;
0138         resXprimeErr = hitStruct.resErrX;
0139         resYprimeErr = hitStruct.resErrY;
0140
0141         const RectangularPlaneBounds* rectangularBound = dynamic_cast<const RectangularPlaneBounds*>(&bound);
0142         if (rectangularBound != nullptr) {
0143           hitStruct.inside = rectangularBound->inside(lPTrk);
0144           length = rectangularBound->length();
0145           width = rectangularBound->width();
0146           hitStruct.localXnorm = 2 * hitStruct.localX / width;
0147           hitStruct.localYnorm = 2 * hitStruct.localY / length;
0148         } else {
0149           throw cms::Exception("Geometry Error")
0150               << "[TrackerValidationVariables] Cannot cast bounds to RectangularPlaneBounds as expected for TPE";
0151         }
0152
0153       } else if (IntSubDetID == StripSubdetector::TID || IntSubDetID == StripSubdetector::TEC) {
0154         // not possible to compute precisely as with Trajectory
0155       } else {
0156         edm::LogWarning("TrackerValidationVariables") << "@SUB=TrackerValidationVariables::fillHitQuantities"
0157                                                       << "No valid tracker subdetector " << IntSubDetID;
0158         continue;
0159       }
0160
0161       resXprime = resXTopol * uOrientation;
0162       resXatTrkY = resXatTrkYTopol;
0163       resYprime = resYTopol * vOrientation;
0164
0165     } else {  // not a detUnit, so must be a virtual 2D-Module
0166       // FIXME: at present only for det units residuals are calculated and filled in the hitStruct
0167       // But in principle this method should also be useable for the gluedDets (2D modules in TIB, TID, TOB, TEC)
0168       // In this case, only orientation should be taken into account for primeResiduals, but not the radial topology
0169       // At present, default values (999.F) are given out
0170     }
0171
0172     hitStruct.resXprime = resXprime;
0173     hitStruct.resXatTrkY = resXatTrkY;
0174     hitStruct.resYprime = resYprime;
0175     hitStruct.resXprimeErr = resXprimeErr;
0176     hitStruct.resYprimeErr = resYprimeErr;
0177
0178     v_avhitout.push_back(hitStruct);
0179   }
0180 }
0181
0182 void TrackerValidationVariables::fillHitQuantities(const Trajectory* trajectory, std::vector<AVHitStruct>& v_avhitout) {
0183   TrajectoryStateCombiner tsoscomb;
0184
0185   const std::vector<TrajectoryMeasurement>& tmColl = trajectory->measurements();
0186   for (std::vector<TrajectoryMeasurement>::const_iterator itTraj = tmColl.begin(); itTraj != tmColl.end(); ++itTraj) {
0187     if (!itTraj->updatedState().isValid())
0188       continue;
0189
0190     TrajectoryStateOnSurface tsos = tsoscomb(itTraj->forwardPredictedState(), itTraj->backwardPredictedState());
0191     if (!tsos.isValid())
0192       continue;
0193     TransientTrackingRecHit::ConstRecHitPointer hit = itTraj->recHit();
0194
0195     if (!hit->isValid() || hit->geographicalId().det() != DetId::Tracker)
0196       continue;
0197
0198     AVHitStruct hitStruct;
0199     const DetId& hit_detId = hit->geographicalId();
0200     unsigned int IntRawDetID = (hit_detId.rawId());
0201     unsigned int IntSubDetID = (hit_detId.subdetId());
0202
0203     if (IntSubDetID == 0)
0204       continue;
0205
0206     if (IntSubDetID == PixelSubdetector::PixelBarrel || IntSubDetID == PixelSubdetector::PixelEndcap) {
0207       const SiPixelRecHit* prechit = dynamic_cast<const SiPixelRecHit*>(
0208           hit.get());  //to be used to get the associated cluster and the cluster probability
0209       if (prechit->isOnEdge())
0210         hitStruct.isOnEdgePixel = true;
0211       if (prechit->hasBadPixels())
0212         hitStruct.isOtherBadPixel = true;
0213     }
0214
0215     //first calculate residuals in cartesian coordinates in the local module coordinate system
0216
0217     LocalPoint lPHit = hit->localPosition();
0218     LocalPoint lPTrk = tsos.localPosition();
0219     LocalVector lVTrk = tsos.localDirection();
0220
0221     hitStruct.localAlpha = atan2(lVTrk.x(), lVTrk.z());  // wrt. normal tg(alpha)=x/z
0222     hitStruct.localBeta = atan2(lVTrk.y(), lVTrk.z());   // wrt. normal tg(beta)= y/z
0223
0224     LocalError errHit = hit->localPositionError();
0225     // no need to add  APE to hitError anymore
0226     // AlgebraicROOTObject<2>::SymMatrix mat = asSMatrix<2>(hit->parametersError());
0227     // LocalError errHit = LocalError( mat(0,0),mat(0,1),mat(1,1) );
0228     LocalError errTrk = tsos.localError().positionError();
0229
0230     //check for negative error values: track error can have negative value, if matrix inversion fails (very rare case)
0231     //hit error should always give positive values
0232     if (errHit.xx() < 0. || errHit.yy() < 0. || errTrk.xx() < 0. || errTrk.yy() < 0.) {
0233       edm::LogError("TrackerValidationVariables")
0234           << "@SUB=TrackerValidationVariables::fillHitQuantities"
0235           << "One of the squared error methods gives negative result"
0236           << "\n\terrHit.xx()\terrHit.yy()\terrTrk.xx()\terrTrk.yy()"
0237           << "\n\t" << errHit.xx() << "\t" << errHit.yy() << "\t" << errTrk.xx() << "\t" << errTrk.yy();
0238       continue;
0239     }
0240
0241     align::LocalVector res = lPTrk - lPHit;
0242
0243     float resXErr = std::sqrt(errHit.xx() + errTrk.xx());
0244     float resYErr = std::sqrt(errHit.yy() + errTrk.yy());
0245
0246     hitStruct.resX = res.x();
0247     hitStruct.resY = res.y();
0248     hitStruct.resErrX = resXErr;
0249     hitStruct.resErrY = resYErr;
0250
0251     // hitStruct.localX = lPhit.x();
0252     // hitStruct.localY = lPhit.y();
0253     // EM: use predictions for local coordinates
0254     hitStruct.localX = lPTrk.x();
0255     hitStruct.localY = lPTrk.y();
0256
0257     // now calculate residuals taking global orientation of modules and radial topology in TID/TEC into account
0258     float resXprime(999.F), resYprime(999.F);
0259     float resXatTrkY(999.F);
0260     float resXprimeErr(999.F), resYprimeErr(999.F);
0261
0262     if (hit->detUnit()) {  // is it a single physical module?
0263       const GeomDetUnit& detUnit = *(hit->detUnit());
0264       float uOrientation(-999.F), vOrientation(-999.F);
0265       float resXTopol(999.F), resYTopol(999.F);
0266       float resXatTrkYTopol(999.F);
0267
0268       const Surface& surface = hit->detUnit()->surface();
0269       const BoundPlane& boundplane = hit->detUnit()->surface();
0270       const Bounds& bound = boundplane.bounds();
0271
0272       float length = 0;
0273       float width = 0;
0274
0275       LocalPoint lPModule(0., 0., 0.), lUDirection(1., 0., 0.), lVDirection(0., 1., 0.);
0276       GlobalPoint gPModule = surface.toGlobal(lPModule), gUDirection = surface.toGlobal(lUDirection),
0277                   gVDirection = surface.toGlobal(lVDirection);
0278
0279       if (IntSubDetID == PixelSubdetector::PixelBarrel || IntSubDetID == StripSubdetector::TIB ||
0280           IntSubDetID == StripSubdetector::TOB) {
0281         uOrientation = deltaPhi(gUDirection.barePhi(), gPModule.barePhi()) >= 0. ? +1.F : -1.F;
0282         vOrientation = gVDirection.z() - gPModule.z() >= 0 ? +1.F : -1.F;
0283         resXTopol = res.x();
0284         resXatTrkYTopol = res.x();
0285         resYTopol = res.y();
0286         resXprimeErr = resXErr;
0287         resYprimeErr = resYErr;
0288
0289         const RectangularPlaneBounds* rectangularBound = dynamic_cast<const RectangularPlaneBounds*>(&bound);
0290         if (rectangularBound != nullptr) {
0291           hitStruct.inside = rectangularBound->inside(lPTrk);
0292           length = rectangularBound->length();
0293           width = rectangularBound->width();
0294           hitStruct.localXnorm = 2 * hitStruct.localX / width;
0295           hitStruct.localYnorm = 2 * hitStruct.localY / length;
0296         } else {
0297           throw cms::Exception("Geometry Error") << "[TrackerValidationVariables] Cannot cast bounds to "
0298                                                     "RectangularPlaneBounds as expected for TPB, TIB and TOB";
0299         }
0300
0301       } else if (IntSubDetID == PixelSubdetector::PixelEndcap) {
0302         uOrientation = gUDirection.perp() - gPModule.perp() >= 0 ? +1.F : -1.F;
0303         vOrientation = deltaPhi(gVDirection.barePhi(), gPModule.barePhi()) >= 0. ? +1.F : -1.F;
0304         resXTopol = res.x();
0305         resXatTrkYTopol = res.x();
0306         resYTopol = res.y();
0307         resXprimeErr = resXErr;
0308         resYprimeErr = resYErr;
0309
0310         const RectangularPlaneBounds* rectangularBound = dynamic_cast<const RectangularPlaneBounds*>(&bound);
0311         if (rectangularBound != nullptr) {
0312           hitStruct.inside = rectangularBound->inside(lPTrk);
0313           length = rectangularBound->length();
0314           width = rectangularBound->width();
0315           hitStruct.localXnorm = 2 * hitStruct.localX / width;
0316           hitStruct.localYnorm = 2 * hitStruct.localY / length;
0317         } else {
0318           throw cms::Exception("Geometry Error")
0319               << "[TrackerValidationVariables] Cannot cast bounds to RectangularPlaneBounds as expected for TPE";
0320         }
0321
0322       } else if (IntSubDetID == StripSubdetector::TID || IntSubDetID == StripSubdetector::TEC) {
0323         uOrientation = deltaPhi(gUDirection.barePhi(), gPModule.barePhi()) >= 0. ? +1.F : -1.F;
0324         vOrientation = gVDirection.perp() - gPModule.perp() >= 0. ? +1.F : -1.F;
0325
0326         if (!dynamic_cast<const RadialStripTopology*>(&detUnit.type().topology()))
0327           continue;
0328         const RadialStripTopology& topol = dynamic_cast<const RadialStripTopology&>(detUnit.type().topology());
0329
0330         MeasurementPoint measHitPos = topol.measurementPosition(lPHit);
0331         MeasurementPoint measTrkPos = topol.measurementPosition(lPTrk);
0332
0333         MeasurementError measHitErr = topol.measurementError(lPHit, errHit);
0334         MeasurementError measTrkErr = topol.measurementError(lPTrk, errTrk);
0335
0336         if (measHitErr.uu() < 0. || measHitErr.vv() < 0. || measTrkErr.uu() < 0. || measTrkErr.vv() < 0.) {
0337           edm::LogError("TrackerValidationVariables")
0338               << "@SUB=TrackerValidationVariables::fillHitQuantities"
0339               << "One of the squared error methods gives negative result"
0340               << "\n\tmeasHitErr.uu()\tmeasHitErr.vv()\tmeasTrkErr.uu()\tmeasTrkErr.vv()"
0341               << "\n\t" << measHitErr.uu() << "\t" << measHitErr.vv() << "\t" << measTrkErr.uu() << "\t"
0342               << measTrkErr.vv();
0343           continue;
0344         }
0345
0346         float localStripLengthHit = topol.localStripLength(lPHit);
0347         float localStripLengthTrk = topol.localStripLength(lPTrk);
0348         float phiHit = topol.stripAngle(measHitPos.x());
0349         float phiTrk = topol.stripAngle(measTrkPos.x());
0350         float r_0 = topol.originToIntersection();
0351
0352         resXTopol = (phiTrk - phiHit) * r_0;
0353         //      resXTopol = (tan(phiTrk)-tan(phiHit))*r_0;
0354
0355         LocalPoint LocalHitPosCor = topol.localPosition(MeasurementPoint(measHitPos.x(), measTrkPos.y()));
0356         resXatTrkYTopol = lPTrk.x() - LocalHitPosCor.x();
0357
0358         //resYTopol = measTrkPos.y()*localStripLengthTrk - measHitPos.y()*localStripLengthHit;
0359         float cosPhiHit(cos(phiHit)), cosPhiTrk(cos(phiTrk)), sinPhiHit(sin(phiHit)), sinPhiTrk(sin(phiTrk));
0360         float l_0 = r_0 - topol.detHeight() / 2;
0361         resYTopol = measTrkPos.y() * localStripLengthTrk - measHitPos.y() * localStripLengthHit +
0362                     l_0 * (1 / cosPhiTrk - 1 / cosPhiHit);
0363
0364         resXprimeErr = std::sqrt(measHitErr.uu() + measTrkErr.uu()) * topol.angularWidth() * r_0;
0365         //resYprimeErr = std::sqrt(measHitErr.vv()*localStripLengthHit*localStripLengthHit + measTrkErr.vv()*localStripLengthTrk*localStripLengthTrk);
0366         float helpSummand = l_0 * l_0 * topol.angularWidth() * topol.angularWidth() *
0367                             (sinPhiHit * sinPhiHit / pow(cosPhiHit, 4) * measHitErr.uu() +
0368                              sinPhiTrk * sinPhiTrk / pow(cosPhiTrk, 4) * measTrkErr.uu());
0369         resYprimeErr = std::sqrt(measHitErr.vv() * localStripLengthHit * localStripLengthHit +
0370                                  measTrkErr.vv() * localStripLengthTrk * localStripLengthTrk + helpSummand);
0371
0372         const TrapezoidalPlaneBounds* trapezoidalBound = dynamic_cast<const TrapezoidalPlaneBounds*>(&bound);
0373         if (trapezoidalBound != nullptr) {
0374           hitStruct.inside = trapezoidalBound->inside(lPTrk);
0375           length = trapezoidalBound->length();
0376           width = trapezoidalBound->width();
0377           //float widthAtHalfLength = trapezoidalBound->widthAtHalfLength();
0378
0379           //    int yAxisOrientation=trapezoidalBound->yAxisOrientation();
0380           // for trapezoidal shape modules, scale with as function of local y coordinate
0381           //    float widthAtlocalY=width-(1-yAxisOrientation*2*lPTrk.y()/length)*(width-widthAtHalfLength);
0382           //    hitStruct.localXnorm = 2*hitStruct.localX/widthAtlocalY;
0383           hitStruct.localXnorm = 2 * hitStruct.localX / width;
0384           hitStruct.localYnorm = 2 * hitStruct.localY / length;
0385         } else {
0386           throw cms::Exception("Geometry Error") << "[TrackerValidationVariables] Cannot cast bounds to "
0387                                                     "TrapezoidalPlaneBounds as expected for TID and TEC";
0388         }
0389
0390       } else {
0391         edm::LogWarning("TrackerValidationVariables") << "@SUB=TrackerValidationVariables::fillHitQuantities"
0392                                                       << "No valid tracker subdetector " << IntSubDetID;
0393         continue;
0394       }
0395
0396       resXprime = resXTopol * uOrientation;
0397       resXatTrkY = resXatTrkYTopol;
0398       resYprime = resYTopol * vOrientation;
0399
0400     } else {  // not a detUnit, so must be a virtual 2D-Module
0401       // FIXME: at present only for det units residuals are calculated and filled in the hitStruct
0402       // But in principle this method should also be useable for the gluedDets (2D modules in TIB, TID, TOB, TEC)
0403       // In this case, only orientation should be taken into account for primeResiduals, but not the radial topology
0404       // At present, default values (999.F) are given out
0405     }
0406
0407     hitStruct.resXprime = resXprime;
0408     hitStruct.resXatTrkY = resXatTrkY;
0409     hitStruct.resYprime = resYprime;
0410     hitStruct.resXprimeErr = resXprimeErr;
0411     hitStruct.resYprimeErr = resYprimeErr;
0412
0413     hitStruct.rawDetId = IntRawDetID;
0414     hitStruct.phi = tsos.globalDirection().phi();
0415     hitStruct.eta = tsos.globalDirection().eta();
0416
0417     v_avhitout.push_back(hitStruct);
0418   }
0419 }
0420
0421 void TrackerValidationVariables::fillTrackQuantities(const edm::Event& event,
0422                                                      const edm::EventSetup& eventSetup,
0423                                                      std::vector<AVTrackStruct>& v_avtrackout) {
0424   fillTrackQuantities(
0425       event, eventSetup, [](const reco::Track&) -> bool { return true; }, v_avtrackout);
0426 }
0427
0428 void TrackerValidationVariables::fillTrackQuantities(const edm::Event& event,
0429                                                      const edm::EventSetup& eventSetup,
0430                                                      std::function<bool(const reco::Track&)> trackFilter,
0431                                                      std::vector<AVTrackStruct>& v_avtrackout) {
0432   const MagneticField& magneticField = eventSetup.getData(magneticFieldToken_);
0433
0434   edm::Handle<reco::TrackCollection> tracksH;
0435   event.getByToken(tracksToken_, tracksH);
0436   if (!tracksH.isValid())
0437     return;
0438   auto const& tracks = *tracksH;
0439   auto ntrk = tracks.size();
0440   LogDebug("TrackerValidationVariables") << "Track collection size " << ntrk;
0441
0442   edm::Handle<std::vector<Trajectory>> trajsH;
0443   event.getByToken(trajCollectionToken_, trajsH);
0444   bool yesTraj = trajsH.isValid();
0445   std::vector<Trajectory> const* trajs = nullptr;
0446   if (yesTraj)
0447     trajs = &(*trajsH);
0448   if (yesTraj)
0449     assert(trajs->size() == tracks.size());
0450
0451   Trajectory const* trajectory = nullptr;
0452   for (unsigned int i = 0; i < ntrk; ++i) {
0453     auto const& track = tracks[i];
0454     if (yesTraj)
0455       trajectory = &(*trajs)[i];
0456
0457     if (!trackFilter(track))
0458       continue;
0459
0460     AVTrackStruct trackStruct;
0461
0462     trackStruct.p = track.p();
0463     trackStruct.pt = track.pt();
0464     trackStruct.ptError = track.ptError();
0465     trackStruct.px = track.px();
0466     trackStruct.py = track.py();
0467     trackStruct.pz = track.pz();
0468     trackStruct.eta = track.eta();
0469     trackStruct.phi = track.phi();
0470     trackStruct.chi2 = track.chi2();
0471     trackStruct.chi2Prob = TMath::Prob(track.chi2(), track.ndof());
0472     trackStruct.normchi2 = track.normalizedChi2();
0473     GlobalPoint gPoint(track.vx(), track.vy(), track.vz());
0474     double theLocalMagFieldInInverseGeV = magneticField.inInverseGeV(gPoint).z();
0475     trackStruct.kappa = -track.charge() * theLocalMagFieldInInverseGeV / track.pt();
0476     trackStruct.charge = track.charge();
0477     trackStruct.d0 = track.d0();
0478     trackStruct.dz = track.dz();
0479     trackStruct.numberOfValidHits = track.numberOfValidHits();
0480     trackStruct.numberOfLostHits = track.numberOfLostHits();
0481     if (trajectory)
0482       fillHitQuantities(trajectory, trackStruct.hits);
0483     else
0484       fillHitQuantities(track, trackStruct.hits);
0485
0486     v_avtrackout.push_back(trackStruct);
0487   }
0488 }