Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 09:43:59 

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(
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");  // 300
0144   double zMax = conf_.getParameter<double>("HitMapZMax");      // 300.0; //cm
0145 
0146   int rBin = conf_.getParameter<unsigned int>("HitMapsRBin");  // 120;
0147   double rMax = conf_.getParameter<double>("HitMapRMax");      // 120.0; //cm
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     // leads to differences in axsis between samples??
0163     // int nModules = binByRawId_.size();
0164     // Hits_perDetId_ = iBooker.book1D(histname+AlgoName, histname+AlgoName,
0165     // nModules, static_cast<double>(nModules) -0.5,
0166     // static_cast<double>(nModules) -0.5);
0167     Hits_perDetId_ = iBooker.book1D(histname + AlgoName, histname + AlgoName, 16601, -0.5, 16600.5);
0168     Hits_perDetId_->setAxisTitle("rawId Bins");
0169 
0170     //// impossible takes too much memory :(
0171     //  std::stringstream binLabel;
0172     //  for( std::map<int,int>::iterator it = binByRawId_.begin(); it !=
0173     //  binByRawId_.end(); ++it ){
0174     //    binLabel.str() = "";
0175     //    binLabel << (*it).first;
0176     //    Hits_perDetId_->getTH1()->GetXaxis()->SetBinLabel( (*it).second +1,
0177     //    binLabel.str().c_str());
0178     //  }
0179   }
0180 }
0181 //
0182 // -- Analyse
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   // fill only once - not yet in beginJob since no access to geometry
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;  // some number > 2pi
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         // edm::LogInfo("Alignment") <<">  isolated: "<< isolated << " jetPt "<<
0238         // (*itJet).pt() <<" deltaR: "<< deltaR(*(*it),(*itJet)) ;
0239       }
0240       minJetDeltaR_->Fill(std::sqrt(minJetDeltaR2));
0241     }
0242 
0243     double minTrackDeltaR2 = 10 * 10;  // some number > 2pi
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       // since 2_1_X local hit positions are transient. taking center of the hit
0300       // module for now. The alternative would be the coarse estimation or a
0301       // refit.
0302       // const GlobalPoint globP( gd->toGlobal( hit->localPosition() ) );
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);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!