Alignment/src/TkAlCaRecoMonitor.cc

0001 /*
0002  *  See header file for a description of this class.
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   // copy configuration object to use it in bookHistograms
0033   conf_ = iConfig;
0034 }
0035
0036 void TkAlCaRecoMonitor::bookHistograms(DQMStore::IBooker &iBooker, edm::Run const &, edm::EventSetup const &) {
0037   std::string histname;  // for naming the histograms according to algorithm used
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(histname + AlgoName,
0076                                  histname + AlgoName,
0077                                  reco::TrackBase::qualitySize,
0078                                  -0.5,
0079                                  static_cast<double>(reco::TrackBase::qualitySize) - 0.5);
0080   TrackQuality_->setAxisTitle("quality");
0081   for (int i = 0; i < reco::TrackBase::qualitySize; ++i) {
0082     TrackQuality_->getTH1()->GetXaxis()->SetBinLabel(
0083         i + 1, reco::TrackBase::qualityName(reco::TrackBase::TrackQuality(i)).c_str());
0084   }
0085
0086   unsigned int SumChargeBin = conf_.getParameter<unsigned int>("SumChargeBin");
0087   double SumChargeMin = conf_.getParameter<double>("SumChargeMin");
0088   double SumChargeMax = conf_.getParameter<double>("SumChargeMax");
0089
0090   histname = "SumCharge_";
0091   sumCharge_ = iBooker.book1D(histname + AlgoName, histname + AlgoName, SumChargeBin, SumChargeMin, SumChargeMax);
0092   sumCharge_->setAxisTitle("#SigmaCharge");
0093
0094   unsigned int TrackCurvatureBin = conf_.getParameter<unsigned int>("TrackCurvatureBin");
0095   double TrackCurvatureMin = conf_.getParameter<double>("TrackCurvatureMin");
0096   double TrackCurvatureMax = conf_.getParameter<double>("TrackCurvatureMax");
0097
0098   histname = "TrackCurvature_";
0099   TrackCurvature_ =
0100       iBooker.book1D(histname + AlgoName, histname + AlgoName, TrackCurvatureBin, TrackCurvatureMin, TrackCurvatureMax);
0101   TrackCurvature_->setAxisTitle("#kappa track");
0102
0103   if (runsOnReco_) {
0104     unsigned int JetPtBin = conf_.getParameter<unsigned int>("JetPtBin");
0105     double JetPtMin = conf_.getParameter<double>("JetPtMin");
0106     double JetPtMax = conf_.getParameter<double>("JetPtMax");
0107
0108     histname = "JetPt_";
0109     jetPt_ = iBooker.book1D(histname + AlgoName, histname + AlgoName, JetPtBin, JetPtMin, JetPtMax);
0110     jetPt_->setAxisTitle("jet p_{T} / GeV");
0111
0112     unsigned int MinJetDeltaRBin = conf_.getParameter<unsigned int>("MinJetDeltaRBin");
0113     double MinJetDeltaRMin = conf_.getParameter<double>("MinJetDeltaRMin");
0114     double MinJetDeltaRMax = conf_.getParameter<double>("MinJetDeltaRMax");
0115
0116     histname = "MinJetDeltaR_";
0117     minJetDeltaR_ =
0118         iBooker.book1D(histname + AlgoName, histname + AlgoName, MinJetDeltaRBin, MinJetDeltaRMin, MinJetDeltaRMax);
0119     minJetDeltaR_->setAxisTitle("minimal Jet #DeltaR / rad");
0120   } else {
0121     jetPt_ = nullptr;
0122     minJetDeltaR_ = nullptr;
0123   }
0124
0125   unsigned int MinTrackDeltaRBin = conf_.getParameter<unsigned int>("MinTrackDeltaRBin");
0126   double MinTrackDeltaRMin = conf_.getParameter<double>("MinTrackDeltaRMin");
0127   double MinTrackDeltaRMax = conf_.getParameter<double>("MinTrackDeltaRMax");
0128
0129   histname = "MinTrackDeltaR_";
0130   minTrackDeltaR_ =
0131       iBooker.book1D(histname + AlgoName, histname + AlgoName, MinTrackDeltaRBin, MinTrackDeltaRMin, MinTrackDeltaRMax);
0132   minTrackDeltaR_->setAxisTitle("minimal Track #DeltaR / rad");
0133
0134   unsigned int TrackEfficiencyBin = conf_.getParameter<unsigned int>("TrackEfficiencyBin");
0135   double TrackEfficiencyMin = conf_.getParameter<double>("TrackEfficiencyMin");
0136   double TrackEfficiencyMax = conf_.getParameter<double>("TrackEfficiencyMax");
0137
0138   histname = "AlCaRecoTrackEfficiency_";
0139   AlCaRecoTrackEfficiency_ = iBooker.book1D(
0140       histname + AlgoName, histname + AlgoName, TrackEfficiencyBin, TrackEfficiencyMin, TrackEfficiencyMax);
0141   Labels l_tp, l_rtp;
0142   labelsForToken(referenceTrackProducer_, l_rtp);
0143   labelsForToken(trackProducer_, l_tp);
0144   AlCaRecoTrackEfficiency_->setAxisTitle("n(" + std::string(l_tp.module) + ") / n(" + std::string(l_rtp.module) + ")");
0145
0146   int zBin = conf_.getParameter<unsigned int>("HitMapsZBin");  // 300
0147   double zMax = conf_.getParameter<double>("HitMapZMax");      // 300.0; //cm
0148
0149   int rBin = conf_.getParameter<unsigned int>("HitMapsRBin");  // 120;
0150   double rMax = conf_.getParameter<double>("HitMapRMax");      // 120.0; //cm
0151
0152   histname = "Hits_ZvsR_";
0153   double rMin = 0.0;
0154   if (useSignedR_)
0155     rMin = -rMax;
0156
0157   Hits_ZvsR_ = iBooker.book2D(histname + AlgoName, histname + AlgoName, zBin, -zMax, zMax, rBin, rMin, rMax);
0158
0159   histname = "Hits_XvsY_";
0160   Hits_XvsY_ = iBooker.book2D(histname + AlgoName, histname + AlgoName, rBin, -rMax, rMax, rBin, -rMax, rMax);
0161
0162   if (fillRawIdMap_) {
0163     histname = "Hits_perDetId_";
0164
0165     // leads to differences in axsis between samples??
0166     // int nModules = binByRawId_.size();
0167     // Hits_perDetId_ = iBooker.book1D(histname+AlgoName, histname+AlgoName,
0168     // nModules, static_cast<double>(nModules) -0.5,
0169     // static_cast<double>(nModules) -0.5);
0170     Hits_perDetId_ = iBooker.book1D(histname + AlgoName, histname + AlgoName, 16601, -0.5, 16600.5);
0171     Hits_perDetId_->setAxisTitle("rawId Bins");
0172
0173     //// impossible takes too much memory :(
0174     //  std::stringstream binLabel;
0175     //  for( std::map<int,int>::iterator it = binByRawId_.begin(); it !=
0176     //  binByRawId_.end(); ++it ){
0177     //    binLabel.str() = "";
0178     //    binLabel << (*it).first;
0179     //    Hits_perDetId_->getTH1()->GetXaxis()->SetBinLabel( (*it).second +1,
0180     //    binLabel.str().c_str());
0181     //  }
0182   }
0183 }
0184 //
0185 // -- Analyse
0186 //
0187 //*************************************************************
0188 void TkAlCaRecoMonitor::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup)
0189 //*************************************************************
0190 {
0191   edm::Handle<reco::TrackCollection> trackCollection;
0192   iEvent.getByToken(trackProducer_, trackCollection);
0193   if (!trackCollection.isValid()) {
0194     edm::LogError("Alignment") << "invalid trackcollection encountered!";
0195     return;
0196   }
0197
0198   edm::Handle<reco::TrackCollection> referenceTrackCollection;
0199   iEvent.getByToken(referenceTrackProducer_, referenceTrackCollection);
0200   if (!trackCollection.isValid()) {
0201     edm::LogError("Alignment") << "invalid reference track-collection encountered!";
0202     return;
0203   }
0204
0205   const auto &geometry = iSetup.getHandle(tkGeomToken_);
0206   if (!geometry.isValid()) {
0207     edm::LogError("Alignment") << "invalid geometry found in event setup!";
0208   }
0209
0210   const auto &magneticField = iSetup.getHandle(mfToken_);
0211   if (!magneticField.isValid()) {
0212     edm::LogError("Alignment") << "invalid magnetic field configuration encountered!";
0213     return;
0214   }
0215
0216   edm::Handle<reco::CaloJetCollection> jets;
0217   if (runsOnReco_) {
0218     iEvent.getByToken(jetCollection_, jets);
0219     if (!jets.isValid()) {
0220       edm::LogError("Alignment") << "no jet collection found in event!";
0221     }
0222   }
0223   // fill only once - not yet in beginJob since no access to geometry
0224   if (fillRawIdMap_ && binByRawId_.empty())
0225     this->fillRawIdMap(*geometry);
0226
0227   AlCaRecoTrackEfficiency_->Fill(static_cast<double>((*trackCollection).size()) / (*referenceTrackCollection).size());
0228
0229   double sumOfCharges = 0;
0230   for (const auto &track : *trackCollection) {
0231     double dR2 = 0;
0232     if (runsOnReco_) {
0233       double minJetDeltaR2 = 10 * 10;  // some number > 2pi
0234       for (const auto &itJet : *jets) {
0235         jetPt_->Fill(itJet.pt());
0236         dR2 = deltaR2(track, itJet);
0237         if (itJet.pt() > maxJetPt_ && dR2 < minJetDeltaR2)
0238           minJetDeltaR2 = dR2;
0239
0240         // edm::LogInfo("Alignment") <<">  isolated: "<< isolated << " jetPt "<<
0241         // (*itJet).pt() <<" deltaR: "<< deltaR(*(*it),(*itJet)) ;
0242       }
0243       minJetDeltaR_->Fill(std::sqrt(minJetDeltaR2));
0244     }
0245
0246     double minTrackDeltaR2 = 10 * 10;  // some number > 2pi
0247     for (const auto &track2 : *trackCollection) {
0248       dR2 = deltaR2(track, track2);
0249       if (dR2 < minTrackDeltaR2 && dR2 > 1e-12)
0250         minTrackDeltaR2 = dR2;
0251     }
0252
0253     for (int i = 0; i < reco::TrackBase::qualitySize; ++i) {
0254       if (track.quality(reco::TrackBase::TrackQuality(i))) {
0255         TrackQuality_->Fill(i);
0256       }
0257     }
0258
0259     GlobalPoint gPoint(track.vx(), track.vy(), track.vz());
0260     double B = magneticField->inTesla(gPoint).z();
0261     double curv = -track.charge() * 0.002998 * B / track.pt();
0262     TrackCurvature_->Fill(curv);
0263
0264     track.charge() > 0 ? TrackPtPositive_->Fill(track.pt()) : TrackPtNegative_->Fill(track.pt());
0265
0266     minTrackDeltaR_->Fill(std::sqrt(minTrackDeltaR2));
0267     fillHitmaps(track, *geometry);
0268     sumOfCharges += track.charge();
0269   }
0270
0271   sumCharge_->Fill(sumOfCharges);
0272
0273   if (fillInvariantMass_) {
0274     if ((*trackCollection).size() == 2) {
0275       TLorentzVector track0(
0276           (*trackCollection).at(0).px(),
0277           (*trackCollection).at(0).py(),
0278           (*trackCollection).at(0).pz(),
0279           sqrt(((*trackCollection).at(0).p() * (*trackCollection).at(0).p()) + daughterMass_ * daughterMass_));
0280       TLorentzVector track1(
0281           (*trackCollection).at(1).px(),
0282           (*trackCollection).at(1).py(),
0283           (*trackCollection).at(1).pz(),
0284           sqrt(((*trackCollection).at(1).p() * (*trackCollection).at(1).p()) + daughterMass_ * daughterMass_));
0285       TLorentzVector mother = track0 + track1;
0286
0287       invariantMass_->Fill(mother.M());
0288     } else {
0289       edm::LogInfo("Alignment") << "wrong number of tracks trackcollection encountered: " << (*trackCollection).size();
0290     }
0291   }
0292 }
0293
0294 //*************************************************************
0295 void TkAlCaRecoMonitor::fillHitmaps(const reco::Track &track, const TrackerGeometry &geometry)
0296 //*************************************************************
0297 {
0298   for (auto const &iHit : track.recHits()) {
0299     if (iHit->isValid()) {
0300       const DetId geoId(iHit->geographicalId());
0301       const GeomDet *gd = geometry.idToDet(geoId);
0302       // since 2_1_X local hit positions are transient. taking center of the hit
0303       // module for now. The alternative would be the coarse estimation or a
0304       // refit.
0305       // const GlobalPoint globP( gd->toGlobal( hit->localPosition() ) );
0306       const GlobalPoint globP(gd->toGlobal(Local3DPoint(0., 0., 0.)));
0307       double r = sqrt(globP.x() * globP.x() + globP.y() * globP.y());
0308       if (useSignedR_)
0309         r *= globP.y() / fabs(globP.y());
0310
0311       Hits_ZvsR_->Fill(globP.z(), r);
0312       Hits_XvsY_->Fill(globP.x(), globP.y());
0313
0314       if (fillRawIdMap_)
0315         Hits_perDetId_->Fill(binByRawId_[geoId.rawId()]);
0316     }
0317   }
0318 }
0319
0320 //*************************************************************
0321 void TkAlCaRecoMonitor::fillRawIdMap(const TrackerGeometry &geometry)
0322 //*************************************************************
0323 {
0324   std::vector<int> sortedRawIds;
0325   for (const auto &iDetId : geometry.detUnitIds()) {
0326     sortedRawIds.push_back(iDetId.rawId());
0327   }
0328   std::sort(sortedRawIds.begin(), sortedRawIds.end());
0329
0330   int i = 0;
0331   for (const auto &iRawId : sortedRawIds) {
0332     binByRawId_[iRawId] = i;
0333     ++i;
0334   }
0335 }
0336
0337 //*************************************************************
0338 void TkAlCaRecoMonitor::fillDescriptions(edm::ConfigurationDescriptions &descriptions)
0339 //*************************************************************
0340 {
0341   edm::ParameterSetDescription desc;
0342   desc.setComment("Generic track analyzer to check ALCARECO Tracker Alignment specific sample quantities");
0343   desc.add<edm::InputTag>("TrackProducer", edm::InputTag("generalTracks"));
0344   desc.add<edm::InputTag>("ReferenceTrackProducer", edm::InputTag("generalTrakcs"));
0345   desc.add<edm::InputTag>("CaloJetCollection", edm::InputTag("ak4CaloJets"));
0346   desc.add<std::string>("AlgoName", "testTkAlCaReco");
0347   desc.add<std::string>("FolderName", "TkAlCaRecoMonitor");
0348   desc.add<double>("daughterMass", kMuonMass_)->setComment("GeV");
0349   desc.add<double>("maxJetPt", 10.)->setComment("GeV");
0350   desc.add<bool>("fillInvariantMass", false);
0351   desc.add<bool>("runsOnReco", false);
0352   desc.add<bool>("useSignedR", false);
0353   desc.add<bool>("fillRawIdMap", false);
0354   desc.add<unsigned int>("MassBin", 100);
0355   desc.add<double>("MassMin", 0.0);
0356   desc.add<double>("MassMax", 100.0);
0357   desc.add<unsigned int>("TrackPtBin", 110);
0358   desc.add<double>("TrackPtMin", 0.0);
0359   desc.add<double>("TrackPtMax", 110.);
0360   desc.add<unsigned int>("TrackCurvatureBin", 2000);
0361   desc.add<double>("TrackCurvatureMin", -0.01)->setComment("1/GeV");
0362   desc.add<double>("TrackCurvatureMax", 0.01)->setComment("1/GeV");
0363   desc.add<unsigned int>("SumChargeBin", 11);
0364   desc.add<double>("SumChargeMin", -5.5);
0365   desc.add<double>("SumChargeMax", 5.5);
0366   desc.add<unsigned int>("JetPtBin", 100);
0367   desc.add<double>("JetPtMin", 0.0);
0368   desc.add<double>("JetPtMax", 50.0);
0369   desc.add<unsigned int>("MinJetDeltaRBin", 100);
0370   desc.add<double>("MinJetDeltaRMin", 0);
0371   desc.add<double>("MinJetDeltaRMax", 10);
0372   desc.add<unsigned int>("MinTrackDeltaRBin", 100);
0373   desc.add<double>("MinTrackDeltaRMin", 0);
0374   desc.add<double>("MinTrackDeltaRMax", 3.2);
0375   desc.add<unsigned int>("TrackEfficiencyBin", 102);
0376   desc.add<double>("TrackEfficiencyMin", -0.01);
0377   desc.add<double>("TrackEfficiencyMax", 1.01);
0378   desc.add<unsigned int>("HitMapsZBin", 300);
0379   desc.add<double>("HitMapZMax", 300.)->setComment("cm");
0380   desc.add<unsigned int>("HitMapsRBin", 120);
0381   desc.add<double>("HitMapRMax", 120.)->setComment("cm");
0382   descriptions.addWithDefaultLabel(desc);
0383 }
0384
0385 DEFINE_FWK_MODULE(TkAlCaRecoMonitor);