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