File indexing completed on 2022-01-24 01:11:31
0001
0002
0003
0004
0005
0006
0007 #include "DQM/TrackingMonitor/interface/TrackSplittingMonitor.h"
0008 #include "DQMServices/Core/interface/DQMStore.h"
0009 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0010 #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
0011 #include "DataFormats/TrackReco/interface/Track.h"
0012 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0013 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 #include "FWCore/Utilities/interface/InputTag.h"
0016 #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
0017 #include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
0018 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0019 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0020 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0021 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0022 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0023
0024 #include <string>
0025
0026 TrackSplittingMonitor::TrackSplittingMonitor(const edm::ParameterSet& iConfig)
0027 : dqmStore_(edm::Service<DQMStore>().operator->()),
0028 conf_(iConfig),
0029 mfToken_(esConsumes()),
0030 tkGeomToken_(esConsumes()),
0031 dtGeomToken_(esConsumes()),
0032 cscGeomToken_(esConsumes()),
0033 rpcGeomToken_(esConsumes()) {
0034 splitTracks_ = conf_.getParameter<edm::InputTag>("splitTrackCollection");
0035 splitMuons_ = conf_.getParameter<edm::InputTag>("splitMuonCollection");
0036 splitTracksToken_ = consumes<std::vector<reco::Track> >(splitTracks_);
0037 splitMuonsToken_ = mayConsume<std::vector<reco::Muon> >(splitMuons_);
0038
0039 plotMuons_ = conf_.getParameter<bool>("ifPlotMuons");
0040
0041
0042 pixelHitsPerLeg_ = conf_.getParameter<int>("pixelHitsPerLeg");
0043 totalHitsPerLeg_ = conf_.getParameter<int>("totalHitsPerLeg");
0044 d0Cut_ = conf_.getParameter<double>("d0Cut");
0045 dzCut_ = conf_.getParameter<double>("dzCut");
0046 ptCut_ = conf_.getParameter<double>("ptCut");
0047 norchiCut_ = conf_.getParameter<double>("norchiCut");
0048 }
0049
0050 TrackSplittingMonitor::~TrackSplittingMonitor() = default;
0051
0052 void TrackSplittingMonitor::bookHistograms(DQMStore::IBooker& ibooker,
0053 edm::Run const& ,
0054 edm::EventSetup const& )
0055
0056 {
0057 std::string MEFolderName = conf_.getParameter<std::string>("FolderName");
0058 ibooker.setCurrentFolder(MEFolderName);
0059
0060
0061 int ddxyBin = conf_.getParameter<int>("ddxyBin");
0062 double ddxyMin = conf_.getParameter<double>("ddxyMin");
0063 double ddxyMax = conf_.getParameter<double>("ddxyMax");
0064
0065 int ddzBin = conf_.getParameter<int>("ddzBin");
0066 double ddzMin = conf_.getParameter<double>("ddzMin");
0067 double ddzMax = conf_.getParameter<double>("ddzMax");
0068
0069 int dphiBin = conf_.getParameter<int>("dphiBin");
0070 double dphiMin = conf_.getParameter<double>("dphiMin");
0071 double dphiMax = conf_.getParameter<double>("dphiMax");
0072
0073 int dthetaBin = conf_.getParameter<int>("dthetaBin");
0074 double dthetaMin = conf_.getParameter<double>("dthetaMin");
0075 double dthetaMax = conf_.getParameter<double>("dthetaMax");
0076
0077 int dptBin = conf_.getParameter<int>("dptBin");
0078 double dptMin = conf_.getParameter<double>("dptMin");
0079 double dptMax = conf_.getParameter<double>("dptMax");
0080
0081 int dcurvBin = conf_.getParameter<int>("dcurvBin");
0082 double dcurvMin = conf_.getParameter<double>("dcurvMin");
0083 double dcurvMax = conf_.getParameter<double>("dcurvMax");
0084
0085 int normBin = conf_.getParameter<int>("normBin");
0086 double normMin = conf_.getParameter<double>("normMin");
0087 double normMax = conf_.getParameter<double>("normMax");
0088
0089
0090 ddxyAbsoluteResiduals_tracker_ =
0091 ibooker.book1D("ddxyAbsoluteResiduals_tracker", "ddxyAbsoluteResiduals_tracker", ddxyBin, ddxyMin, ddxyMax);
0092 ddzAbsoluteResiduals_tracker_ =
0093 ibooker.book1D("ddzAbsoluteResiduals_tracker", "ddzAbsoluteResiduals_tracker", ddzBin, ddzMin, ddzMax);
0094 dphiAbsoluteResiduals_tracker_ =
0095 ibooker.book1D("dphiAbsoluteResiduals_tracker", "dphiAbsoluteResiduals_tracker", dphiBin, dphiMin, dphiMax);
0096 dthetaAbsoluteResiduals_tracker_ = ibooker.book1D(
0097 "dthetaAbsoluteResiduals_tracker", "dthetaAbsoluteResiduals_tracker", dthetaBin, dthetaMin, dthetaMax);
0098 dptAbsoluteResiduals_tracker_ =
0099 ibooker.book1D("dptAbsoluteResiduals_tracker", "dptAbsoluteResiduals_tracker", dptBin, dptMin, dptMax);
0100 dcurvAbsoluteResiduals_tracker_ =
0101 ibooker.book1D("dcurvAbsoluteResiduals_tracker", "dcurvAbsoluteResiduals_tracker", dcurvBin, dcurvMin, dcurvMax);
0102
0103 ddxyNormalizedResiduals_tracker_ =
0104 ibooker.book1D("ddxyNormalizedResiduals_tracker", "ddxyNormalizedResiduals_tracker", normBin, normMin, normMax);
0105 ddzNormalizedResiduals_tracker_ =
0106 ibooker.book1D("ddzNormalizedResiduals_tracker", "ddzNormalizedResiduals_tracker", normBin, normMin, normMax);
0107 dphiNormalizedResiduals_tracker_ =
0108 ibooker.book1D("dphiNormalizedResiduals_tracker", "dphiNormalizedResiduals_tracker", normBin, normMin, normMax);
0109 dthetaNormalizedResiduals_tracker_ = ibooker.book1D(
0110 "dthetaNormalizedResiduals_tracker", "dthetaNormalizedResiduals_tracker", normBin, normMin, normMax);
0111 dptNormalizedResiduals_tracker_ =
0112 ibooker.book1D("dptNormalizedResiduals_tracker", "dptNormalizedResiduals_tracker", normBin, normMin, normMax);
0113 dcurvNormalizedResiduals_tracker_ =
0114 ibooker.book1D("dcurvNormalizedResiduals_tracker", "dcurvNormalizedResiduals_tracker", normBin, normMin, normMax);
0115
0116 if (plotMuons_) {
0117 ddxyAbsoluteResiduals_global_ =
0118 ibooker.book1D("ddxyAbsoluteResiduals_global", "ddxyAbsoluteResiduals_global", ddxyBin, ddxyMin, ddxyMax);
0119 ddzAbsoluteResiduals_global_ =
0120 ibooker.book1D("ddzAbsoluteResiduals_global", "ddzAbsoluteResiduals_global", ddzBin, ddzMin, ddzMax);
0121 dphiAbsoluteResiduals_global_ =
0122 ibooker.book1D("dphiAbsoluteResiduals_global", "dphiAbsoluteResiduals_global", dphiBin, dphiMin, dphiMax);
0123 dthetaAbsoluteResiduals_global_ = ibooker.book1D(
0124 "dthetaAbsoluteResiduals_global", "dthetaAbsoluteResiduals_global", dthetaBin, dthetaMin, dthetaMax);
0125 dptAbsoluteResiduals_global_ =
0126 ibooker.book1D("dptAbsoluteResiduals_global", "dptAbsoluteResiduals_global", dptBin, dptMin, dptMax);
0127 dcurvAbsoluteResiduals_global_ =
0128 ibooker.book1D("dcurvAbsoluteResiduals_global", "dcurvAbsoluteResiduals_global", dcurvBin, dcurvMin, dcurvMax);
0129
0130 ddxyNormalizedResiduals_global_ =
0131 ibooker.book1D("ddxyNormalizedResiduals_global", "ddxyNormalizedResiduals_global", normBin, normMin, normMax);
0132 ddzNormalizedResiduals_global_ =
0133 ibooker.book1D("ddzNormalizedResiduals_global", "ddzNormalizedResiduals_global", normBin, normMin, normMax);
0134 dphiNormalizedResiduals_global_ =
0135 ibooker.book1D("dphiNormalizedResiduals_global", "dphiNormalizedResiduals_global", normBin, normMin, normMax);
0136 dthetaNormalizedResiduals_global_ = ibooker.book1D(
0137 "dthetaNormalizedResiduals_global", "dthetaNormalizedResiduals_global", normBin, normMin, normMax);
0138 dptNormalizedResiduals_global_ =
0139 ibooker.book1D("dptNormalizedResiduals_global", "dptNormalizedResiduals_global", normBin, normMin, normMax);
0140 dcurvNormalizedResiduals_global_ =
0141 ibooker.book1D("dcurvNormalizedResiduals_global", "dcurvNormalizedResiduals_global", normBin, normMin, normMax);
0142 }
0143
0144 ddxyAbsoluteResiduals_tracker_->setAxisTitle("(#delta d_{xy})/#sqrt{2} [#mum]");
0145 ddxyAbsoluteResiduals_tracker_->setAxisTitle("(#delta d_{z})/#sqrt{2} [#mum]");
0146 ddxyAbsoluteResiduals_tracker_->setAxisTitle("(#delta #phi)/#sqrt{2} [mrad]");
0147 ddxyAbsoluteResiduals_tracker_->setAxisTitle("(#delta #theta)/#sqrt{2} [mrad]");
0148 ddxyAbsoluteResiduals_tracker_->setAxisTitle("(#delta pT)/#sqrt{2} [GeV]");
0149 ddxyAbsoluteResiduals_tracker_->setAxisTitle("(#delta (1/pT))/#sqrt{2} [GeV^{-1}]");
0150
0151 ddxyNormalizedResiduals_tracker_->setAxisTitle("#delta d_{xy}/#sigma(d_{xy}");
0152 ddxyNormalizedResiduals_tracker_->setAxisTitle("#delta d_{z}/#sigma(d_{z})");
0153 ddxyNormalizedResiduals_tracker_->setAxisTitle("#delta #phi/#sigma(d_{#phi})");
0154 ddxyNormalizedResiduals_tracker_->setAxisTitle("#delta #theta/#sigma(d_{#theta})");
0155 ddxyNormalizedResiduals_tracker_->setAxisTitle("#delta p_{T}/#sigma(p_{T})");
0156 ddxyNormalizedResiduals_tracker_->setAxisTitle("#delta 1/p_{T}/#sigma(1/p_{T})");
0157
0158 if (plotMuons_) {
0159 ddxyAbsoluteResiduals_global_->setAxisTitle("(#delta d_{xy})/#sqrt{2} [#mum]");
0160 ddxyAbsoluteResiduals_global_->setAxisTitle("(#delta d_{z})/#sqrt{2} [#mum]");
0161 ddxyAbsoluteResiduals_global_->setAxisTitle("(#delta #phi)/#sqrt{2} [mrad]");
0162 ddxyAbsoluteResiduals_global_->setAxisTitle("(#delta #theta)/#sqrt{2} [mrad]");
0163 ddxyAbsoluteResiduals_global_->setAxisTitle("(#delta pT)/#sqrt{2} [GeV]");
0164 ddxyAbsoluteResiduals_global_->setAxisTitle("(#delta (1/pT))/#sqrt{2} [GeV^{-1}]");
0165
0166 ddxyNormalizedResiduals_global_->setAxisTitle("#delta d_{xy}/#sigma(d_{xy}");
0167 ddxyNormalizedResiduals_global_->setAxisTitle("#delta d_{z}/#sigma(d_{z})");
0168 ddxyNormalizedResiduals_global_->setAxisTitle("#delta #phi/#sigma(d_{#phi})");
0169 ddxyNormalizedResiduals_global_->setAxisTitle("#delta #theta/#sigma(d_{#theta})");
0170 ddxyNormalizedResiduals_global_->setAxisTitle("#delta p_{T}/#sigma(p_{T})");
0171 ddxyNormalizedResiduals_global_->setAxisTitle("#delta 1/p_{T}/#sigma(1/p_{T})");
0172 }
0173 }
0174
0175
0176
0177
0178 void TrackSplittingMonitor::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0179 theMagField = &iSetup.getData(mfToken_);
0180 theGeometry = &iSetup.getData(tkGeomToken_);
0181 dtGeometry = &iSetup.getData(dtGeomToken_);
0182 cscGeometry = &iSetup.getData(cscGeomToken_);
0183 rpcGeometry = &iSetup.getData(rpcGeomToken_);
0184
0185 edm::Handle<std::vector<reco::Track> > splitTracks = iEvent.getHandle(splitTracksToken_);
0186 if (!splitTracks.isValid())
0187 return;
0188
0189 edm::Handle<std::vector<reco::Muon> > splitMuons;
0190 if (plotMuons_) {
0191 splitMuons = iEvent.getHandle(splitMuonsToken_);
0192 }
0193
0194 if (splitTracks->size() == 2) {
0195
0196 edm::LogInfo("TrackSplittingMonitor") << "Split Track size: " << splitTracks->size();
0197
0198
0199 reco::Track track1 = splitTracks->at(0);
0200 reco::Track track2 = splitTracks->at(1);
0201
0202
0203
0204
0205
0206 double nRechits1 = 0;
0207 double nRechitinBPIX1 = 0;
0208 for (trackingRecHit_iterator iHit = track1.recHitsBegin(); iHit != track1.recHitsEnd(); ++iHit) {
0209 if ((*iHit)->isValid()) {
0210 nRechits1++;
0211 int type = (*iHit)->geographicalId().subdetId();
0212 if (type == int(PixelSubdetector::PixelBarrel)) {
0213 ++nRechitinBPIX1;
0214 }
0215 }
0216 }
0217
0218 double nRechits2 = 0;
0219 double nRechitinBPIX2 = 0;
0220 for (trackingRecHit_iterator iHit = track2.recHitsBegin(); iHit != track2.recHitsEnd(); ++iHit) {
0221 if ((*iHit)->isValid()) {
0222 nRechits2++;
0223 int type = (*iHit)->geographicalId().subdetId();
0224 if (type == int(PixelSubdetector::PixelBarrel)) {
0225 ++nRechitinBPIX2;
0226 }
0227 }
0228 }
0229
0230
0231 double d01 = track1.d0();
0232 double dz1 = track1.dz();
0233 double d02 = track2.d0();
0234 double dz2 = track2.dz();
0235
0236
0237 double pt1 = track1.pt();
0238 double pt2 = track2.pt();
0239
0240
0241 double norchi1 = track1.normalizedChi2();
0242 double norchi2 = track2.normalizedChi2();
0243
0244
0245
0246 if ((nRechitinBPIX1 >= pixelHitsPerLeg_) && (nRechitinBPIX1 >= pixelHitsPerLeg_) &&
0247 (nRechits1 >= totalHitsPerLeg_) && (nRechits2 >= totalHitsPerLeg_)) {
0248
0249 if (((fabs(d01) < d0Cut_)) && (fabs(d02) < d0Cut_) && (fabs(dz2) < dzCut_) && (fabs(dz2) < dzCut_)) {
0250
0251 if ((pt1 + pt2) / 2 < ptCut_) {
0252
0253 if ((norchi1 < norchiCut_) && (norchi2 < norchiCut_)) {
0254
0255 edm::LogInfo("TrackSplittingMonitor") << " Setected after all cuts ?";
0256
0257 double ddxyVal = d01 - d02;
0258 double ddzVal = dz1 - dz2;
0259 double dphiVal = track1.phi() - track2.phi();
0260 double dthetaVal = track1.theta() - track2.theta();
0261 double dptVal = pt1 - pt2;
0262 double dcurvVal = (1 / pt1) - (1 / pt2);
0263
0264 double d01ErrVal = track1.d0Error();
0265 double d02ErrVal = track2.d0Error();
0266 double dz1ErrVal = track1.dzError();
0267 double dz2ErrVal = track2.dzError();
0268 double phi1ErrVal = track1.phiError();
0269 double phi2ErrVal = track2.phiError();
0270 double theta1ErrVal = track1.thetaError();
0271 double theta2ErrVal = track2.thetaError();
0272 double pt1ErrVal = track1.ptError();
0273 double pt2ErrVal = track2.ptError();
0274
0275 ddxyAbsoluteResiduals_tracker_->Fill(10000.0 * ddxyVal / sqrt(2.0));
0276 ddxyAbsoluteResiduals_tracker_->Fill(10000.0 * ddzVal / sqrt(2.0));
0277 ddxyAbsoluteResiduals_tracker_->Fill(1000.0 * dphiVal / sqrt(2.0));
0278 ddxyAbsoluteResiduals_tracker_->Fill(1000.0 * dthetaVal / sqrt(2.0));
0279 ddxyAbsoluteResiduals_tracker_->Fill(dptVal / sqrt(2.0));
0280 ddxyAbsoluteResiduals_tracker_->Fill(dcurvVal / sqrt(2.0));
0281
0282 ddxyNormalizedResiduals_tracker_->Fill(ddxyVal / sqrt(d01ErrVal * d01ErrVal + d02ErrVal * d02ErrVal));
0283 ddxyNormalizedResiduals_tracker_->Fill(ddzVal / sqrt(dz1ErrVal * dz1ErrVal + dz2ErrVal * dz2ErrVal));
0284 ddxyNormalizedResiduals_tracker_->Fill(dphiVal / sqrt(phi1ErrVal * phi1ErrVal + phi2ErrVal * phi2ErrVal));
0285 ddxyNormalizedResiduals_tracker_->Fill(dthetaVal /
0286 sqrt(theta1ErrVal * theta1ErrVal + theta2ErrVal * theta2ErrVal));
0287 ddxyNormalizedResiduals_tracker_->Fill(dptVal / sqrt(pt1ErrVal * pt1ErrVal + pt2ErrVal * pt2ErrVal));
0288 ddxyNormalizedResiduals_tracker_->Fill(
0289 dcurvVal / sqrt(pow(pt1ErrVal, 2) / pow(pt1, 4) + pow(pt2ErrVal, 2) / pow(pt2, 4)));
0290
0291
0292 if (plotMuons_ && splitMuons.isValid()) {
0293 int gmCtr = 0;
0294 bool topGlobalMuonFlag = false;
0295 bool bottomGlobalMuonFlag = false;
0296 int topGlobalMuon = -1;
0297 int bottomGlobalMuon = -1;
0298 double topGlobalMuonNorchi2 = 1e10;
0299 double bottomGlobalMuonNorchi2 = 1e10;
0300
0301
0302 for (std::vector<reco::Muon>::const_iterator gmI = splitMuons->begin(); gmI != splitMuons->end(); gmI++) {
0303 if (gmI->isTrackerMuon() && gmI->isStandAloneMuon() && gmI->isGlobalMuon()) {
0304 reco::TrackRef trackerTrackRef1(splitTracks, 0);
0305 reco::TrackRef trackerTrackRef2(splitTracks, 1);
0306
0307 if (gmI->innerTrack() == trackerTrackRef1) {
0308 if (gmI->globalTrack()->normalizedChi2() < topGlobalMuonNorchi2) {
0309 topGlobalMuonFlag = true;
0310 topGlobalMuonNorchi2 = gmI->globalTrack()->normalizedChi2();
0311 topGlobalMuon = gmCtr;
0312 }
0313 }
0314 if (gmI->innerTrack() == trackerTrackRef2) {
0315 if (gmI->globalTrack()->normalizedChi2() < bottomGlobalMuonNorchi2) {
0316 bottomGlobalMuonFlag = true;
0317 bottomGlobalMuonNorchi2 = gmI->globalTrack()->normalizedChi2();
0318 bottomGlobalMuon = gmCtr;
0319 }
0320 }
0321 }
0322 gmCtr++;
0323 }
0324
0325 if (bottomGlobalMuonFlag && topGlobalMuonFlag) {
0326 reco::Muon muonTop = splitMuons->at(topGlobalMuon);
0327 reco::Muon muonBottom = splitMuons->at(bottomGlobalMuon);
0328
0329 reco::TrackRef glb1 = muonTop.globalTrack();
0330 reco::TrackRef glb2 = muonBottom.globalTrack();
0331
0332 double ddxyValGlb = glb1->d0() - glb2->d0();
0333 double ddzValGlb = glb1->dz() - glb2->dz();
0334 double dphiValGlb = glb1->phi() - glb2->phi();
0335 double dthetaValGlb = glb1->theta() - glb2->theta();
0336 double dptValGlb = glb1->pt() - glb2->pt();
0337 double dcurvValGlb = (1 / glb1->pt()) - (1 / glb2->pt());
0338
0339 double d01ErrValGlb = glb1->d0Error();
0340 double d02ErrValGlb = glb2->d0Error();
0341 double dz1ErrValGlb = glb1->dzError();
0342 double dz2ErrValGlb = glb2->dzError();
0343 double phi1ErrValGlb = glb1->phiError();
0344 double phi2ErrValGlb = glb2->phiError();
0345 double theta1ErrValGlb = glb1->thetaError();
0346 double theta2ErrValGlb = glb2->thetaError();
0347 double pt1ErrValGlb = glb1->ptError();
0348 double pt2ErrValGlb = glb2->ptError();
0349
0350 ddxyAbsoluteResiduals_global_->Fill(10000.0 * ddxyValGlb / sqrt(2.0));
0351 ddxyAbsoluteResiduals_global_->Fill(10000.0 * ddzValGlb / sqrt(2.0));
0352 ddxyAbsoluteResiduals_global_->Fill(1000.0 * dphiValGlb / sqrt(2.0));
0353 ddxyAbsoluteResiduals_global_->Fill(1000.0 * dthetaValGlb / sqrt(2.0));
0354 ddxyAbsoluteResiduals_global_->Fill(dptValGlb / sqrt(2.0));
0355 ddxyAbsoluteResiduals_global_->Fill(dcurvValGlb / sqrt(2.0));
0356
0357 ddxyNormalizedResiduals_global_->Fill(ddxyValGlb /
0358 sqrt(d01ErrValGlb * d01ErrValGlb + d02ErrValGlb * d02ErrValGlb));
0359 ddxyNormalizedResiduals_global_->Fill(ddzValGlb /
0360 sqrt(dz1ErrValGlb * dz1ErrValGlb + dz2ErrValGlb * dz2ErrValGlb));
0361 ddxyNormalizedResiduals_global_->Fill(
0362 dphiValGlb / sqrt(phi1ErrValGlb * phi1ErrValGlb + phi2ErrValGlb * phi2ErrValGlb));
0363 ddxyNormalizedResiduals_global_->Fill(
0364 dthetaValGlb / sqrt(theta1ErrValGlb * theta1ErrValGlb + theta2ErrValGlb * theta2ErrValGlb));
0365 ddxyNormalizedResiduals_global_->Fill(dptValGlb /
0366 sqrt(pt1ErrValGlb * pt1ErrValGlb + pt2ErrValGlb * pt2ErrValGlb));
0367 ddxyNormalizedResiduals_global_->Fill(
0368 dcurvValGlb / sqrt(pow(pt1ErrValGlb, 2) / pow(pt1, 4) + pow(pt2ErrValGlb, 2) / pow(pt2, 4)));
0369 }
0370
0371 }
0372 }
0373 }
0374 }
0375 }
0376 }
0377 }
0378
0379 DEFINE_FWK_MODULE(TrackSplittingMonitor);