Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-06-20 00:21:34

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()), mfToken_(esConsumes()) {
0021   conf_ = iConfig;
0022   trackProducer_ = consumes<reco::TrackCollection>(conf_.getParameter<edm::InputTag>("TrackProducer"));
0023   referenceTrackProducer_ =
0024       consumes<reco::TrackCollection>(conf_.getParameter<edm::InputTag>("ReferenceTrackProducer"));
0025   jetCollection_ = mayConsume<reco::CaloJetCollection>(conf_.getParameter<edm::InputTag>("CaloJetCollection"));
0026 }
0027 
0028 void TkAlCaRecoMonitor::bookHistograms(DQMStore::IBooker &iBooker, edm::Run const &, edm::EventSetup const &) {
0029   std::string histname;  // for naming the histograms according to algorithm used
0030 
0031   std::string AlgoName = conf_.getParameter<std::string>("AlgoName");
0032   std::string MEFolderName = conf_.getParameter<std::string>("FolderName");
0033 
0034   daughterMass_ = conf_.getParameter<double>("daughterMass");
0035 
0036   maxJetPt_ = conf_.getParameter<double>("maxJetPt");
0037 
0038   iBooker.setCurrentFolder(MEFolderName + "/TkAlignmentSpecific");
0039   fillInvariantMass_ = conf_.getParameter<bool>("fillInvariantMass");
0040   runsOnReco_ = conf_.getParameter<bool>("runsOnReco");
0041   useSignedR_ = conf_.getParameter<bool>("useSignedR");
0042   fillRawIdMap_ = conf_.getParameter<bool>("fillRawIdMap");
0043 
0044   //
0045   unsigned int MassBin = conf_.getParameter<unsigned int>("MassBin");
0046   double MassMin = conf_.getParameter<double>("MassMin");
0047   double MassMax = conf_.getParameter<double>("MassMax");
0048 
0049   if (fillInvariantMass_) {
0050     histname = "InvariantMass_";
0051     invariantMass_ = iBooker.book1D(histname + AlgoName, histname + AlgoName, MassBin, MassMin, MassMax);
0052     invariantMass_->setAxisTitle("invariant Mass / GeV");
0053   } else {
0054     invariantMass_ = nullptr;
0055   }
0056 
0057   unsigned int TrackPtPositiveBin = conf_.getParameter<unsigned int>("TrackPtBin");
0058   double TrackPtPositiveMin = conf_.getParameter<double>("TrackPtMin");
0059   double TrackPtPositiveMax = conf_.getParameter<double>("TrackPtMax");
0060 
0061   histname = "TrackPtPositive_";
0062   TrackPtPositive_ = iBooker.book1D(
0063       histname + AlgoName, histname + AlgoName, TrackPtPositiveBin, TrackPtPositiveMin, TrackPtPositiveMax);
0064   TrackPtPositive_->setAxisTitle("p_{T} of tracks charge > 0");
0065 
0066   unsigned int TrackPtNegativeBin = conf_.getParameter<unsigned int>("TrackPtBin");
0067   double TrackPtNegativeMin = conf_.getParameter<double>("TrackPtMin");
0068   double TrackPtNegativeMax = conf_.getParameter<double>("TrackPtMax");
0069 
0070   histname = "TrackPtNegative_";
0071   TrackPtNegative_ = iBooker.book1D(
0072       histname + AlgoName, histname + AlgoName, TrackPtNegativeBin, TrackPtNegativeMin, TrackPtNegativeMax);
0073   TrackPtNegative_->setAxisTitle("p_{T} of tracks charge < 0");
0074 
0075   histname = "TrackQuality_";
0076   TrackQuality_ = iBooker.book1D(
0077       histname + AlgoName, histname + AlgoName, reco::TrackBase::qualitySize, -0.5, reco::TrackBase::qualitySize - 0.5);
0078   TrackQuality_->setAxisTitle("quality");
0079   for (int i = 0; i < reco::TrackBase::qualitySize; ++i) {
0080     TrackQuality_->getTH1()->GetXaxis()->SetBinLabel(
0081         i + 1, reco::TrackBase::qualityName(reco::TrackBase::TrackQuality(i)).c_str());
0082   }
0083 
0084   unsigned int SumChargeBin = conf_.getParameter<unsigned int>("SumChargeBin");
0085   double SumChargeMin = conf_.getParameter<double>("SumChargeMin");
0086   double SumChargeMax = conf_.getParameter<double>("SumChargeMax");
0087 
0088   histname = "SumCharge_";
0089   sumCharge_ = iBooker.book1D(histname + AlgoName, histname + AlgoName, SumChargeBin, SumChargeMin, SumChargeMax);
0090   sumCharge_->setAxisTitle("#SigmaCharge");
0091 
0092   unsigned int TrackCurvatureBin = conf_.getParameter<unsigned int>("TrackCurvatureBin");
0093   double TrackCurvatureMin = conf_.getParameter<double>("TrackCurvatureMin");
0094   double TrackCurvatureMax = conf_.getParameter<double>("TrackCurvatureMax");
0095 
0096   histname = "TrackCurvature_";
0097   TrackCurvature_ =
0098       iBooker.book1D(histname + AlgoName, histname + AlgoName, TrackCurvatureBin, TrackCurvatureMin, TrackCurvatureMax);
0099   TrackCurvature_->setAxisTitle("#kappa track");
0100 
0101   if (runsOnReco_) {
0102     unsigned int JetPtBin = conf_.getParameter<unsigned int>("JetPtBin");
0103     double JetPtMin = conf_.getParameter<double>("JetPtMin");
0104     double JetPtMax = conf_.getParameter<double>("JetPtMax");
0105 
0106     histname = "JetPt_";
0107     jetPt_ = iBooker.book1D(histname + AlgoName, histname + AlgoName, JetPtBin, JetPtMin, JetPtMax);
0108     jetPt_->setAxisTitle("jet p_{T} / GeV");
0109 
0110     unsigned int MinJetDeltaRBin = conf_.getParameter<unsigned int>("MinJetDeltaRBin");
0111     double MinJetDeltaRMin = conf_.getParameter<double>("MinJetDeltaRMin");
0112     double MinJetDeltaRMax = conf_.getParameter<double>("MinJetDeltaRMax");
0113 
0114     histname = "MinJetDeltaR_";
0115     minJetDeltaR_ =
0116         iBooker.book1D(histname + AlgoName, histname + AlgoName, MinJetDeltaRBin, MinJetDeltaRMin, MinJetDeltaRMax);
0117     minJetDeltaR_->setAxisTitle("minimal Jet #DeltaR / rad");
0118   } else {
0119     jetPt_ = nullptr;
0120     minJetDeltaR_ = nullptr;
0121   }
0122 
0123   unsigned int MinTrackDeltaRBin = conf_.getParameter<unsigned int>("MinTrackDeltaRBin");
0124   double MinTrackDeltaRMin = conf_.getParameter<double>("MinTrackDeltaRMin");
0125   double MinTrackDeltaRMax = conf_.getParameter<double>("MinTrackDeltaRMax");
0126 
0127   histname = "MinTrackDeltaR_";
0128   minTrackDeltaR_ =
0129       iBooker.book1D(histname + AlgoName, histname + AlgoName, MinTrackDeltaRBin, MinTrackDeltaRMin, MinTrackDeltaRMax);
0130   minTrackDeltaR_->setAxisTitle("minimal Track #DeltaR / rad");
0131 
0132   unsigned int TrackEfficiencyBin = conf_.getParameter<unsigned int>("TrackEfficiencyBin");
0133   double TrackEfficiencyMin = conf_.getParameter<double>("TrackEfficiencyMin");
0134   double TrackEfficiencyMax = conf_.getParameter<double>("TrackEfficiencyMax");
0135 
0136   histname = "AlCaRecoTrackEfficiency_";
0137   AlCaRecoTrackEfficiency_ = iBooker.book1D(
0138       histname + AlgoName, histname + AlgoName, TrackEfficiencyBin, TrackEfficiencyMin, TrackEfficiencyMax);
0139   Labels l_tp, l_rtp;
0140   labelsForToken(referenceTrackProducer_, l_rtp);
0141   labelsForToken(trackProducer_, l_tp);
0142   AlCaRecoTrackEfficiency_->setAxisTitle("n(" + std::string(l_tp.module) + ") / n(" + std::string(l_rtp.module) + ")");
0143 
0144   int zBin = conf_.getParameter<unsigned int>("HitMapsZBin");  // 300
0145   double zMax = conf_.getParameter<double>("HitMapZMax");      // 300.0; //cm
0146 
0147   int rBin = conf_.getParameter<unsigned int>("HitMapsRBin");  // 120;
0148   double rMax = conf_.getParameter<double>("HitMapRMax");      // 120.0; //cm
0149 
0150   histname = "Hits_ZvsR_";
0151   double rMin = 0.0;
0152   if (useSignedR_)
0153     rMin = -rMax;
0154 
0155   Hits_ZvsR_ = iBooker.book2D(histname + AlgoName, histname + AlgoName, zBin, -zMax, zMax, rBin, rMin, rMax);
0156 
0157   histname = "Hits_XvsY_";
0158   Hits_XvsY_ = iBooker.book2D(histname + AlgoName, histname + AlgoName, rBin, -rMax, rMax, rBin, -rMax, rMax);
0159 
0160   if (fillRawIdMap_) {
0161     histname = "Hits_perDetId_";
0162 
0163     // leads to differences in axsis between samples??
0164     // int nModules = binByRawId_.size();
0165     // Hits_perDetId_ = iBooker.book1D(histname+AlgoName, histname+AlgoName,
0166     // nModules, static_cast<double>(nModules) -0.5,
0167     // static_cast<double>(nModules) -0.5);
0168     Hits_perDetId_ = iBooker.book1D(histname + AlgoName, histname + AlgoName, 16601, -0.5, 16600.5);
0169     Hits_perDetId_->setAxisTitle("rawId Bins");
0170 
0171     //// impossible takes too much memory :(
0172     //  std::stringstream binLabel;
0173     //  for( std::map<int,int>::iterator it = binByRawId_.begin(); it !=
0174     //  binByRawId_.end(); ++it ){
0175     //    binLabel.str() = "";
0176     //    binLabel << (*it).first;
0177     //    Hits_perDetId_->getTH1()->GetXaxis()->SetBinLabel( (*it).second +1,
0178     //    binLabel.str().c_str());
0179     //  }
0180   }
0181 }
0182 //
0183 // -- Analyse
0184 //
0185 void TkAlCaRecoMonitor::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0186   edm::Handle<reco::TrackCollection> trackCollection;
0187   iEvent.getByToken(trackProducer_, trackCollection);
0188   if (!trackCollection.isValid()) {
0189     edm::LogError("Alignment") << "invalid trackcollection encountered!";
0190     return;
0191   }
0192 
0193   edm::Handle<reco::TrackCollection> referenceTrackCollection;
0194   iEvent.getByToken(referenceTrackProducer_, referenceTrackCollection);
0195   if (!trackCollection.isValid()) {
0196     edm::LogError("Alignment") << "invalid reference track-collection encountered!";
0197     return;
0198   }
0199 
0200   const auto &geometry = iSetup.getHandle(tkGeomToken_);
0201   if (!geometry.isValid()) {
0202     edm::LogError("Alignment") << "invalid geometry found in event setup!";
0203   }
0204 
0205   const auto &magneticField = iSetup.getHandle(mfToken_);
0206   if (!magneticField.isValid()) {
0207     edm::LogError("Alignment") << "invalid magnetic field configuration encountered!";
0208     return;
0209   }
0210 
0211   edm::Handle<reco::CaloJetCollection> jets;
0212   if (runsOnReco_) {
0213     iEvent.getByToken(jetCollection_, jets);
0214     if (!jets.isValid()) {
0215       edm::LogError("Alignment") << "no jet collection found in event!";
0216     }
0217   }
0218   // fill only once - not yet in beginJob since no access to geometry
0219   if (fillRawIdMap_ && binByRawId_.empty())
0220     this->fillRawIdMap(*geometry);
0221 
0222   AlCaRecoTrackEfficiency_->Fill(static_cast<double>((*trackCollection).size()) / (*referenceTrackCollection).size());
0223 
0224   double sumOfCharges = 0;
0225   for (reco::TrackCollection::const_iterator track = (*trackCollection).begin(); track < (*trackCollection).end();
0226        ++track) {
0227     double dR = 0;
0228     if (runsOnReco_) {
0229       double minJetDeltaR = 10;  // some number > 2pi
0230       for (reco::CaloJetCollection::const_iterator itJet = jets->begin(); itJet != jets->end(); ++itJet) {
0231         jetPt_->Fill((*itJet).pt());
0232         dR = deltaR((*track), (*itJet));
0233         if ((*itJet).pt() > maxJetPt_ && dR < minJetDeltaR)
0234           minJetDeltaR = dR;
0235 
0236         // edm::LogInfo("Alignment") <<">  isolated: "<< isolated << " jetPt "<<
0237         // (*itJet).pt() <<" deltaR: "<< deltaR(*(*it),(*itJet)) ;
0238       }
0239       minJetDeltaR_->Fill(minJetDeltaR);
0240     }
0241 
0242     double minTrackDeltaR = 10;  // some number > 2pi
0243     for (reco::TrackCollection::const_iterator track2 = (*trackCollection).begin(); track2 < (*trackCollection).end();
0244          ++track2) {
0245       dR = deltaR((*track), (*track2));
0246       if (dR < minTrackDeltaR && dR > 1e-6)
0247         minTrackDeltaR = dR;
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     if ((*track).charge() > 0)
0262       TrackPtPositive_->Fill((*track).pt());
0263     if ((*track).charge() < 0)
0264       TrackPtNegative_->Fill((*track).pt());
0265 
0266     minTrackDeltaR_->Fill(minTrackDeltaR);
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 void TkAlCaRecoMonitor::fillHitmaps(const reco::Track &track, const TrackerGeometry &geometry) {
0295   for (trackingRecHit_iterator iHit = track.recHitsBegin(); iHit != track.recHitsEnd(); ++iHit) {
0296     if ((*iHit)->isValid()) {
0297       const TrackingRecHit *hit = (*iHit);
0298       const DetId geoId(hit->geographicalId());
0299       const GeomDet *gd = geometry.idToDet(geoId);
0300       // since 2_1_X local hit positions are transient. taking center of the hit
0301       // module for now. The alternative would be the coarse estimation or a
0302       // refit.
0303       // const GlobalPoint globP( gd->toGlobal( hit->localPosition() ) );
0304       const GlobalPoint globP(gd->toGlobal(Local3DPoint(0., 0., 0.)));
0305       double r = sqrt(globP.x() * globP.x() + globP.y() * globP.y());
0306       if (useSignedR_)
0307         r *= globP.y() / fabs(globP.y());
0308 
0309       Hits_ZvsR_->Fill(globP.z(), r);
0310       Hits_XvsY_->Fill(globP.x(), globP.y());
0311 
0312       if (fillRawIdMap_)
0313         Hits_perDetId_->Fill(binByRawId_[geoId.rawId()]);
0314     }
0315   }
0316 }
0317 
0318 void TkAlCaRecoMonitor::fillRawIdMap(const TrackerGeometry &geometry) {
0319   std::vector<int> sortedRawIds;
0320   for (std::vector<DetId>::const_iterator iDetId = geometry.detUnitIds().begin(); iDetId != geometry.detUnitIds().end();
0321        ++iDetId) {
0322     sortedRawIds.push_back((*iDetId).rawId());
0323   }
0324   std::sort(sortedRawIds.begin(), sortedRawIds.end());
0325 
0326   int i = 0;
0327   for (std::vector<int>::iterator iRawId = sortedRawIds.begin(); iRawId != sortedRawIds.end(); ++iRawId) {
0328     binByRawId_[*iRawId] = i;
0329     ++i;
0330   }
0331 }
0332 
0333 DEFINE_FWK_MODULE(TkAlCaRecoMonitor);