File indexing completed on 2023-10-25 09:43:59
0001
0002
0003
0004
0005
0006 #include "DQMOffline/Alignment/interface/TkAlCaRecoMonitor.h"
0007 #include "DQMServices/Core/interface/DQMStore.h"
0008 #include "DataFormats/Math/interface/deltaR.h"
0009 #include "DataFormats/TrackReco/interface/Track.h"
0010 #include "DataFormats/TrackReco/interface/TrackBase.h"
0011 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013 #include "FWCore/Utilities/interface/InputTag.h"
0014 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0015
0016 #include <string>
0017 #include "TLorentzVector.h"
0018
0019 TkAlCaRecoMonitor::TkAlCaRecoMonitor(const edm::ParameterSet &iConfig)
0020 : tkGeomToken_(esConsumes()),
0021 mfToken_(esConsumes()),
0022 trackProducer_(consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("TrackProducer"))),
0023 referenceTrackProducer_(
0024 consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("ReferenceTrackProducer"))),
0025 jetCollection_(mayConsume<reco::CaloJetCollection>(iConfig.getParameter<edm::InputTag>("CaloJetCollection"))),
0026 daughterMass_(iConfig.getParameter<double>("daughterMass")),
0027 maxJetPt_(iConfig.getParameter<double>("maxJetPt")),
0028 fillInvariantMass_(iConfig.getParameter<bool>("fillInvariantMass")),
0029 fillRawIdMap_(iConfig.getParameter<bool>("fillRawIdMap")),
0030 runsOnReco_(iConfig.getParameter<bool>("runsOnReco")),
0031 useSignedR_(iConfig.getParameter<bool>("useSignedR")) {
0032
0033 conf_ = iConfig;
0034 }
0035
0036 void TkAlCaRecoMonitor::bookHistograms(DQMStore::IBooker &iBooker, edm::Run const &, edm::EventSetup const &) {
0037 std::string histname;
0038
0039 std::string AlgoName = conf_.getParameter<std::string>("AlgoName");
0040 std::string MEFolderName = conf_.getParameter<std::string>("FolderName");
0041 iBooker.setCurrentFolder(MEFolderName + "/TkAlignmentSpecific");
0042
0043
0044 unsigned int MassBin = conf_.getParameter<unsigned int>("MassBin");
0045 double MassMin = conf_.getParameter<double>("MassMin");
0046 double MassMax = conf_.getParameter<double>("MassMax");
0047
0048 if (fillInvariantMass_) {
0049 histname = "InvariantMass_";
0050 invariantMass_ = iBooker.book1D(histname + AlgoName, histname + AlgoName, MassBin, MassMin, MassMax);
0051 invariantMass_->setAxisTitle("invariant Mass / GeV");
0052 } else {
0053 invariantMass_ = nullptr;
0054 }
0055
0056 unsigned int TrackPtPositiveBin = conf_.getParameter<unsigned int>("TrackPtBin");
0057 double TrackPtPositiveMin = conf_.getParameter<double>("TrackPtMin");
0058 double TrackPtPositiveMax = conf_.getParameter<double>("TrackPtMax");
0059
0060 histname = "TrackPtPositive_";
0061 TrackPtPositive_ = iBooker.book1D(
0062 histname + AlgoName, histname + AlgoName, TrackPtPositiveBin, TrackPtPositiveMin, TrackPtPositiveMax);
0063 TrackPtPositive_->setAxisTitle("p_{T} of tracks charge > 0");
0064
0065 unsigned int TrackPtNegativeBin = conf_.getParameter<unsigned int>("TrackPtBin");
0066 double TrackPtNegativeMin = conf_.getParameter<double>("TrackPtMin");
0067 double TrackPtNegativeMax = conf_.getParameter<double>("TrackPtMax");
0068
0069 histname = "TrackPtNegative_";
0070 TrackPtNegative_ = iBooker.book1D(
0071 histname + AlgoName, histname + AlgoName, TrackPtNegativeBin, TrackPtNegativeMin, TrackPtNegativeMax);
0072 TrackPtNegative_->setAxisTitle("p_{T} of tracks charge < 0");
0073
0074 histname = "TrackQuality_";
0075 TrackQuality_ = iBooker.book1D(
0076 histname + AlgoName, histname + AlgoName, reco::TrackBase::qualitySize, -0.5, reco::TrackBase::qualitySize - 0.5);
0077 TrackQuality_->setAxisTitle("quality");
0078 for (int i = 0; i < reco::TrackBase::qualitySize; ++i) {
0079 TrackQuality_->getTH1()->GetXaxis()->SetBinLabel(
0080 i + 1, reco::TrackBase::qualityName(reco::TrackBase::TrackQuality(i)).c_str());
0081 }
0082
0083 unsigned int SumChargeBin = conf_.getParameter<unsigned int>("SumChargeBin");
0084 double SumChargeMin = conf_.getParameter<double>("SumChargeMin");
0085 double SumChargeMax = conf_.getParameter<double>("SumChargeMax");
0086
0087 histname = "SumCharge_";
0088 sumCharge_ = iBooker.book1D(histname + AlgoName, histname + AlgoName, SumChargeBin, SumChargeMin, SumChargeMax);
0089 sumCharge_->setAxisTitle("#SigmaCharge");
0090
0091 unsigned int TrackCurvatureBin = conf_.getParameter<unsigned int>("TrackCurvatureBin");
0092 double TrackCurvatureMin = conf_.getParameter<double>("TrackCurvatureMin");
0093 double TrackCurvatureMax = conf_.getParameter<double>("TrackCurvatureMax");
0094
0095 histname = "TrackCurvature_";
0096 TrackCurvature_ =
0097 iBooker.book1D(histname + AlgoName, histname + AlgoName, TrackCurvatureBin, TrackCurvatureMin, TrackCurvatureMax);
0098 TrackCurvature_->setAxisTitle("#kappa track");
0099
0100 if (runsOnReco_) {
0101 unsigned int JetPtBin = conf_.getParameter<unsigned int>("JetPtBin");
0102 double JetPtMin = conf_.getParameter<double>("JetPtMin");
0103 double JetPtMax = conf_.getParameter<double>("JetPtMax");
0104
0105 histname = "JetPt_";
0106 jetPt_ = iBooker.book1D(histname + AlgoName, histname + AlgoName, JetPtBin, JetPtMin, JetPtMax);
0107 jetPt_->setAxisTitle("jet p_{T} / GeV");
0108
0109 unsigned int MinJetDeltaRBin = conf_.getParameter<unsigned int>("MinJetDeltaRBin");
0110 double MinJetDeltaRMin = conf_.getParameter<double>("MinJetDeltaRMin");
0111 double MinJetDeltaRMax = conf_.getParameter<double>("MinJetDeltaRMax");
0112
0113 histname = "MinJetDeltaR_";
0114 minJetDeltaR_ =
0115 iBooker.book1D(histname + AlgoName, histname + AlgoName, MinJetDeltaRBin, MinJetDeltaRMin, MinJetDeltaRMax);
0116 minJetDeltaR_->setAxisTitle("minimal Jet #DeltaR / rad");
0117 } else {
0118 jetPt_ = nullptr;
0119 minJetDeltaR_ = nullptr;
0120 }
0121
0122 unsigned int MinTrackDeltaRBin = conf_.getParameter<unsigned int>("MinTrackDeltaRBin");
0123 double MinTrackDeltaRMin = conf_.getParameter<double>("MinTrackDeltaRMin");
0124 double MinTrackDeltaRMax = conf_.getParameter<double>("MinTrackDeltaRMax");
0125
0126 histname = "MinTrackDeltaR_";
0127 minTrackDeltaR_ =
0128 iBooker.book1D(histname + AlgoName, histname + AlgoName, MinTrackDeltaRBin, MinTrackDeltaRMin, MinTrackDeltaRMax);
0129 minTrackDeltaR_->setAxisTitle("minimal Track #DeltaR / rad");
0130
0131 unsigned int TrackEfficiencyBin = conf_.getParameter<unsigned int>("TrackEfficiencyBin");
0132 double TrackEfficiencyMin = conf_.getParameter<double>("TrackEfficiencyMin");
0133 double TrackEfficiencyMax = conf_.getParameter<double>("TrackEfficiencyMax");
0134
0135 histname = "AlCaRecoTrackEfficiency_";
0136 AlCaRecoTrackEfficiency_ = iBooker.book1D(
0137 histname + AlgoName, histname + AlgoName, TrackEfficiencyBin, TrackEfficiencyMin, TrackEfficiencyMax);
0138 Labels l_tp, l_rtp;
0139 labelsForToken(referenceTrackProducer_, l_rtp);
0140 labelsForToken(trackProducer_, l_tp);
0141 AlCaRecoTrackEfficiency_->setAxisTitle("n(" + std::string(l_tp.module) + ") / n(" + std::string(l_rtp.module) + ")");
0142
0143 int zBin = conf_.getParameter<unsigned int>("HitMapsZBin");
0144 double zMax = conf_.getParameter<double>("HitMapZMax");
0145
0146 int rBin = conf_.getParameter<unsigned int>("HitMapsRBin");
0147 double rMax = conf_.getParameter<double>("HitMapRMax");
0148
0149 histname = "Hits_ZvsR_";
0150 double rMin = 0.0;
0151 if (useSignedR_)
0152 rMin = -rMax;
0153
0154 Hits_ZvsR_ = iBooker.book2D(histname + AlgoName, histname + AlgoName, zBin, -zMax, zMax, rBin, rMin, rMax);
0155
0156 histname = "Hits_XvsY_";
0157 Hits_XvsY_ = iBooker.book2D(histname + AlgoName, histname + AlgoName, rBin, -rMax, rMax, rBin, -rMax, rMax);
0158
0159 if (fillRawIdMap_) {
0160 histname = "Hits_perDetId_";
0161
0162
0163
0164
0165
0166
0167 Hits_perDetId_ = iBooker.book1D(histname + AlgoName, histname + AlgoName, 16601, -0.5, 16600.5);
0168 Hits_perDetId_->setAxisTitle("rawId Bins");
0169
0170
0171
0172
0173
0174
0175
0176
0177
0178
0179 }
0180 }
0181
0182
0183
0184
0185 void TkAlCaRecoMonitor::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup)
0186
0187 {
0188 edm::Handle<reco::TrackCollection> trackCollection;
0189 iEvent.getByToken(trackProducer_, trackCollection);
0190 if (!trackCollection.isValid()) {
0191 edm::LogError("Alignment") << "invalid trackcollection encountered!";
0192 return;
0193 }
0194
0195 edm::Handle<reco::TrackCollection> referenceTrackCollection;
0196 iEvent.getByToken(referenceTrackProducer_, referenceTrackCollection);
0197 if (!trackCollection.isValid()) {
0198 edm::LogError("Alignment") << "invalid reference track-collection encountered!";
0199 return;
0200 }
0201
0202 const auto &geometry = iSetup.getHandle(tkGeomToken_);
0203 if (!geometry.isValid()) {
0204 edm::LogError("Alignment") << "invalid geometry found in event setup!";
0205 }
0206
0207 const auto &magneticField = iSetup.getHandle(mfToken_);
0208 if (!magneticField.isValid()) {
0209 edm::LogError("Alignment") << "invalid magnetic field configuration encountered!";
0210 return;
0211 }
0212
0213 edm::Handle<reco::CaloJetCollection> jets;
0214 if (runsOnReco_) {
0215 iEvent.getByToken(jetCollection_, jets);
0216 if (!jets.isValid()) {
0217 edm::LogError("Alignment") << "no jet collection found in event!";
0218 }
0219 }
0220
0221 if (fillRawIdMap_ && binByRawId_.empty())
0222 this->fillRawIdMap(*geometry);
0223
0224 AlCaRecoTrackEfficiency_->Fill(static_cast<double>((*trackCollection).size()) / (*referenceTrackCollection).size());
0225
0226 double sumOfCharges = 0;
0227 for (const auto &track : *trackCollection) {
0228 double dR2 = 0;
0229 if (runsOnReco_) {
0230 double minJetDeltaR2 = 10 * 10;
0231 for (const auto &itJet : *jets) {
0232 jetPt_->Fill(itJet.pt());
0233 dR2 = deltaR2(track, itJet);
0234 if (itJet.pt() > maxJetPt_ && dR2 < minJetDeltaR2)
0235 minJetDeltaR2 = dR2;
0236
0237
0238
0239 }
0240 minJetDeltaR_->Fill(std::sqrt(minJetDeltaR2));
0241 }
0242
0243 double minTrackDeltaR2 = 10 * 10;
0244 for (const auto &track2 : *trackCollection) {
0245 dR2 = deltaR2(track, track2);
0246 if (dR2 < minTrackDeltaR2 && dR2 > 1e-12)
0247 minTrackDeltaR2 = dR2;
0248 }
0249
0250 for (int i = 0; i < reco::TrackBase::qualitySize; ++i) {
0251 if (track.quality(reco::TrackBase::TrackQuality(i))) {
0252 TrackQuality_->Fill(i);
0253 }
0254 }
0255
0256 GlobalPoint gPoint(track.vx(), track.vy(), track.vz());
0257 double B = magneticField->inTesla(gPoint).z();
0258 double curv = -track.charge() * 0.002998 * B / track.pt();
0259 TrackCurvature_->Fill(curv);
0260
0261 track.charge() > 0 ? TrackPtPositive_->Fill(track.pt()) : TrackPtNegative_->Fill(track.pt());
0262
0263 minTrackDeltaR_->Fill(std::sqrt(minTrackDeltaR2));
0264 fillHitmaps(track, *geometry);
0265 sumOfCharges += track.charge();
0266 }
0267
0268 sumCharge_->Fill(sumOfCharges);
0269
0270 if (fillInvariantMass_) {
0271 if ((*trackCollection).size() == 2) {
0272 TLorentzVector track0(
0273 (*trackCollection).at(0).px(),
0274 (*trackCollection).at(0).py(),
0275 (*trackCollection).at(0).pz(),
0276 sqrt(((*trackCollection).at(0).p() * (*trackCollection).at(0).p()) + daughterMass_ * daughterMass_));
0277 TLorentzVector track1(
0278 (*trackCollection).at(1).px(),
0279 (*trackCollection).at(1).py(),
0280 (*trackCollection).at(1).pz(),
0281 sqrt(((*trackCollection).at(1).p() * (*trackCollection).at(1).p()) + daughterMass_ * daughterMass_));
0282 TLorentzVector mother = track0 + track1;
0283
0284 invariantMass_->Fill(mother.M());
0285 } else {
0286 edm::LogInfo("Alignment") << "wrong number of tracks trackcollection encountered: " << (*trackCollection).size();
0287 }
0288 }
0289 }
0290
0291
0292 void TkAlCaRecoMonitor::fillHitmaps(const reco::Track &track, const TrackerGeometry &geometry)
0293
0294 {
0295 for (auto const &iHit : track.recHits()) {
0296 if (iHit->isValid()) {
0297 const DetId geoId(iHit->geographicalId());
0298 const GeomDet *gd = geometry.idToDet(geoId);
0299
0300
0301
0302
0303 const GlobalPoint globP(gd->toGlobal(Local3DPoint(0., 0., 0.)));
0304 double r = sqrt(globP.x() * globP.x() + globP.y() * globP.y());
0305 if (useSignedR_)
0306 r *= globP.y() / fabs(globP.y());
0307
0308 Hits_ZvsR_->Fill(globP.z(), r);
0309 Hits_XvsY_->Fill(globP.x(), globP.y());
0310
0311 if (fillRawIdMap_)
0312 Hits_perDetId_->Fill(binByRawId_[geoId.rawId()]);
0313 }
0314 }
0315 }
0316
0317
0318 void TkAlCaRecoMonitor::fillRawIdMap(const TrackerGeometry &geometry)
0319
0320 {
0321 std::vector<int> sortedRawIds;
0322 for (const auto &iDetId : geometry.detUnitIds()) {
0323 sortedRawIds.push_back(iDetId.rawId());
0324 }
0325 std::sort(sortedRawIds.begin(), sortedRawIds.end());
0326
0327 int i = 0;
0328 for (const auto &iRawId : sortedRawIds) {
0329 binByRawId_[iRawId] = i;
0330 ++i;
0331 }
0332 }
0333
0334
0335 void TkAlCaRecoMonitor::fillDescriptions(edm::ConfigurationDescriptions &descriptions)
0336
0337 {
0338 edm::ParameterSetDescription desc;
0339 desc.setComment("Generic track analyzer to check ALCARECO Tracker Alignment specific sample quantities");
0340 desc.add<edm::InputTag>("TrackProducer", edm::InputTag("generalTracks"));
0341 desc.add<edm::InputTag>("ReferenceTrackProducer", edm::InputTag("generalTrakcs"));
0342 desc.add<edm::InputTag>("CaloJetCollection", edm::InputTag("ak4CaloJets"));
0343 desc.add<std::string>("AlgoName", "testTkAlCaReco");
0344 desc.add<std::string>("FolderName", "TkAlCaRecoMonitor");
0345 desc.add<double>("daughterMass", kMuonMass_)->setComment("GeV");
0346 desc.add<double>("maxJetPt", 10.)->setComment("GeV");
0347 desc.add<bool>("fillInvariantMass", false);
0348 desc.add<bool>("runsOnReco", false);
0349 desc.add<bool>("useSignedR", false);
0350 desc.add<bool>("fillRawIdMap", false);
0351 desc.add<unsigned int>("MassBin", 100);
0352 desc.add<double>("MassMin", 0.0);
0353 desc.add<double>("MassMax", 100.0);
0354 desc.add<unsigned int>("TrackPtBin", 110);
0355 desc.add<double>("TrackPtMin", 0.0);
0356 desc.add<double>("TrackPtMax", 110.);
0357 desc.add<unsigned int>("TrackCurvatureBin", 2000);
0358 desc.add<double>("TrackCurvatureMin", -0.01)->setComment("1/GeV");
0359 desc.add<double>("TrackCurvatureMax", 0.01)->setComment("1/GeV");
0360 desc.add<unsigned int>("SumChargeBin", 11);
0361 desc.add<double>("SumChargeMin", -5.5);
0362 desc.add<double>("SumChargeMax", 5.5);
0363 desc.add<unsigned int>("JetPtBin", 100);
0364 desc.add<double>("JetPtMin", 0.0);
0365 desc.add<double>("JetPtMax", 50.0);
0366 desc.add<unsigned int>("MinJetDeltaRBin", 100);
0367 desc.add<double>("MinJetDeltaRMin", 0);
0368 desc.add<double>("MinJetDeltaRMax", 10);
0369 desc.add<unsigned int>("MinTrackDeltaRBin", 100);
0370 desc.add<double>("MinTrackDeltaRMin", 0);
0371 desc.add<double>("MinTrackDeltaRMax", 3.2);
0372 desc.add<unsigned int>("TrackEfficiencyBin", 102);
0373 desc.add<double>("TrackEfficiencyMin", -0.01);
0374 desc.add<double>("TrackEfficiencyMax", 1.01);
0375 desc.add<unsigned int>("HitMapsZBin", 300);
0376 desc.add<double>("HitMapZMax", 300.)->setComment("cm");
0377 desc.add<unsigned int>("HitMapsRBin", 120);
0378 desc.add<double>("HitMapRMax", 120.)->setComment("cm");
0379 descriptions.addWithDefaultLabel(desc);
0380 }
0381
0382 DEFINE_FWK_MODULE(TkAlCaRecoMonitor);