File indexing completed on 2021-11-22 06:26:42
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);
0075 if (prechit->isOnEdge())
0076 hitStruct.isOnEdgePixel = true;
0077 if (prechit->hasBadPixels())
0078 hitStruct.isOtherBadPixel = true;
0079 }
0080
0081 auto lPTrk = trajParams[h].position();
0082 auto lVTrk = trajParams[h].direction();
0083
0084 auto gtrkdirup = hit->surface()->toGlobal(lVTrk);
0085
0086 hitStruct.rawDetId = IntRawDetID;
0087 hitStruct.phi = gtrkdirup.phi();
0088 hitStruct.eta = gtrkdirup.eta();
0089
0090 hitStruct.localAlpha = std::atan2(lVTrk.x(), lVTrk.z());
0091 hitStruct.localBeta = std::atan2(lVTrk.y(), lVTrk.z());
0092
0093 hitStruct.resX = residuals.residualX(h);
0094 hitStruct.resY = residuals.residualY(h);
0095 hitStruct.resErrX = hitStruct.resX / residuals.pullX(h);
0096 hitStruct.resErrY = hitStruct.resY / residuals.pullY(h);
0097
0098
0099
0100
0101 hitStruct.localX = lPTrk.x();
0102 hitStruct.localY = lPTrk.y();
0103
0104
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()) {
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
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 {
0166
0167
0168
0169
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());
0209 if (prechit->isOnEdge())
0210 hitStruct.isOnEdgePixel = true;
0211 if (prechit->hasBadPixels())
0212 hitStruct.isOtherBadPixel = true;
0213 }
0214
0215
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());
0222 hitStruct.localBeta = atan2(lVTrk.y(), lVTrk.z());
0223
0224 LocalError errHit = hit->localPositionError();
0225
0226
0227
0228 LocalError errTrk = tsos.localError().positionError();
0229
0230
0231
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
0252
0253
0254 hitStruct.localX = lPTrk.x();
0255 hitStruct.localY = lPTrk.y();
0256
0257
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()) {
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
0354
0355 LocalPoint LocalHitPosCor = topol.localPosition(MeasurementPoint(measHitPos.x(), measTrkPos.y()));
0356 resXatTrkYTopol = lPTrk.x() - LocalHitPosCor.x();
0357
0358
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
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
0378
0379
0380
0381
0382
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 {
0401
0402
0403
0404
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 }