Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:13:51

0001 #include "DQMServices/Examples/interface/DQMExample_Step1.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 
0004 // Geometry
0005 #include "DataFormats/Math/interface/LorentzVector.h"
0006 #include "DataFormats/Math/interface/deltaPhi.h"
0007 #include "DataFormats/Math/interface/deltaR.h"
0008 #include "TLorentzVector.h"
0009 
0010 #include <cmath>
0011 #include <cstdio>
0012 #include <iomanip>
0013 #include <iostream>
0014 #include <sstream>
0015 #include <string>
0016 
0017 //
0018 // -------------------------------------- Constructor
0019 // --------------------------------------------
0020 //
0021 DQMExample_Step1::DQMExample_Step1(const edm::ParameterSet &ps) {
0022   edm::LogInfo("DQMExample_Step1") << "Constructor  DQMExample_Step1::DQMExample_Step1 " << std::endl;
0023 
0024   // Get parameters from configuration file
0025   theElectronCollection_ = consumes<reco::GsfElectronCollection>(ps.getParameter<edm::InputTag>("electronCollection"));
0026   theCaloJetCollection_ = consumes<reco::CaloJetCollection>(ps.getParameter<edm::InputTag>("caloJetCollection"));
0027   thePfMETCollection_ = consumes<reco::PFMETCollection>(ps.getParameter<edm::InputTag>("pfMETCollection"));
0028   theConversionCollection_ =
0029       consumes<reco::ConversionCollection>(ps.getParameter<edm::InputTag>("conversionsCollection"));
0030   thePVCollection_ = consumes<reco::VertexCollection>(ps.getParameter<edm::InputTag>("PVCollection"));
0031   theBSCollection_ = consumes<reco::BeamSpot>(ps.getParameter<edm::InputTag>("beamSpotCollection"));
0032   triggerEvent_ = consumes<trigger::TriggerEvent>(ps.getParameter<edm::InputTag>("TriggerEvent"));
0033   triggerResults_ = consumes<edm::TriggerResults>(ps.getParameter<edm::InputTag>("TriggerResults"));
0034   triggerFilter_ = ps.getParameter<edm::InputTag>("TriggerFilter");
0035   triggerPath_ = ps.getParameter<std::string>("TriggerPath");
0036 
0037   // cuts:
0038   ptThrL1_ = ps.getUntrackedParameter<double>("PtThrL1");
0039   ptThrL2_ = ps.getUntrackedParameter<double>("PtThrL2");
0040   ptThrJet_ = ps.getUntrackedParameter<double>("PtThrJet");
0041   ptThrMet_ = ps.getUntrackedParameter<double>("PtThrMet");
0042 }
0043 
0044 //
0045 // -- Destructor
0046 //
0047 DQMExample_Step1::~DQMExample_Step1() {
0048   edm::LogInfo("DQMExample_Step1") << "Destructor DQMExample_Step1::~DQMExample_Step1 " << std::endl;
0049 }
0050 
0051 //
0052 // -------------------------------------- beginRun
0053 // --------------------------------------------
0054 //
0055 void DQMExample_Step1::dqmBeginRun(edm::Run const &, edm::EventSetup const &) {
0056   edm::LogInfo("DQMExample_Step1") << "DQMExample_Step1::beginRun" << std::endl;
0057 }
0058 //
0059 // -------------------------------------- bookHistos
0060 // --------------------------------------------
0061 //
0062 void DQMExample_Step1::bookHistograms(DQMStore::IBooker &ibooker_, edm::Run const &, edm::EventSetup const &) {
0063   edm::LogInfo("DQMExample_Step1") << "DQMExample_Step1::bookHistograms" << std::endl;
0064 
0065   // book at beginRun
0066   bookHistos(ibooker_);
0067 }
0068 //
0069 // -------------------------------------- Analyze
0070 // --------------------------------------------
0071 //
0072 void DQMExample_Step1::analyze(edm::Event const &e, edm::EventSetup const &eSetup) {
0073   edm::LogInfo("DQMExample_Step1") << "DQMExample_Step1::analyze" << std::endl;
0074 
0075   //-------------------------------
0076   //--- Vertex Info
0077   //-------------------------------
0078   edm::Handle<reco::VertexCollection> vertexHandle;
0079   e.getByToken(thePVCollection_, vertexHandle);
0080   if (!vertexHandle.isValid()) {
0081     edm::LogError("DQMClientExample") << "invalid collection: vertex"
0082                                       << "\n";
0083     return;
0084   }
0085 
0086   int vertex_number = vertexHandle->size();
0087   reco::VertexCollection::const_iterator v = vertexHandle->begin();
0088 
0089   math::XYZPoint PVPoint(-999, -999, -999);
0090   if (vertex_number != 0)
0091     PVPoint = math::XYZPoint(v->position().x(), v->position().y(), v->position().z());
0092 
0093   PVPoint_ = PVPoint;
0094 
0095   //-------------------------------
0096   //--- MET
0097   //-------------------------------
0098   edm::Handle<reco::PFMETCollection> pfMETCollection;
0099   e.getByToken(thePfMETCollection_, pfMETCollection);
0100   if (!pfMETCollection.isValid()) {
0101     edm::LogError("DQMClientExample") << "invalid collection: MET"
0102                                       << "\n";
0103     return;
0104   }
0105   //-------------------------------
0106   //--- Electrons
0107   //-------------------------------
0108   edm::Handle<reco::GsfElectronCollection> electronCollection;
0109   e.getByToken(theElectronCollection_, electronCollection);
0110   if (!electronCollection.isValid()) {
0111     edm::LogError("DQMClientExample") << "invalid collection: electrons"
0112                                       << "\n";
0113     return;
0114   }
0115 
0116   float nEle = 0;
0117   int posEle = 0, negEle = 0;
0118   const reco::GsfElectron *ele1 = nullptr;
0119   const reco::GsfElectron *ele2 = nullptr;
0120   for (reco::GsfElectronCollection::const_iterator recoElectron = electronCollection->begin();
0121        recoElectron != electronCollection->end();
0122        ++recoElectron) {
0123     // decreasing pT
0124     if (MediumEle(e, eSetup, *recoElectron)) {
0125       if (!ele1 && recoElectron->pt() > ptThrL1_)
0126         ele1 = &(*recoElectron);
0127 
0128       else if (!ele2 && recoElectron->pt() > ptThrL2_)
0129         ele2 = &(*recoElectron);
0130     }
0131 
0132     if (recoElectron->charge() == 1)
0133       posEle++;
0134     else if (recoElectron->charge() == -1)
0135       negEle++;
0136 
0137   }  // end of loop over electrons
0138 
0139   nEle = posEle + negEle;
0140 
0141   //-------------------------------
0142   //--- Jets
0143   //-------------------------------
0144   edm::Handle<reco::CaloJetCollection> caloJetCollection;
0145   e.getByToken(theCaloJetCollection_, caloJetCollection);
0146   if (!caloJetCollection.isValid()) {
0147     edm::LogError("DQMClientExample") << "invalid collection: jets"
0148                                       << "\n";
0149     return;
0150   }
0151 
0152   int nJet = 0;
0153   const reco::CaloJet *jet1 = nullptr;
0154   const reco::CaloJet *jet2 = nullptr;
0155 
0156   for (reco::CaloJetCollection::const_iterator i_calojet = caloJetCollection->begin();
0157        i_calojet != caloJetCollection->end();
0158        ++i_calojet) {
0159     // remove jet-ele matching
0160     if (ele1)
0161       if (Distance(*i_calojet, *ele1) < 0.3)
0162         continue;
0163 
0164     if (ele2)
0165       if (Distance(*i_calojet, *ele2) < 0.3)
0166         continue;
0167 
0168     if (i_calojet->pt() < ptThrJet_)
0169       continue;
0170 
0171     nJet++;
0172 
0173     if (!jet1)
0174       jet1 = &(*i_calojet);
0175 
0176     else if (!jet2)
0177       jet2 = &(*i_calojet);
0178   }
0179 
0180   // ---------------------------
0181   // ---- Analyze Trigger Event
0182   // ---------------------------
0183 
0184   // check what is in the menu
0185   edm::Handle<edm::TriggerResults> hltresults;
0186   e.getByToken(triggerResults_, hltresults);
0187 
0188   if (!hltresults.isValid()) {
0189     edm::LogError("DQMClientExample") << "invalid collection: TriggerResults"
0190                                       << "\n";
0191     return;
0192   }
0193 
0194   bool hasFired = false;
0195   const edm::TriggerNames &trigNames = e.triggerNames(*hltresults);
0196   unsigned int numTriggers = trigNames.size();
0197 
0198   for (unsigned int hltIndex = 0; hltIndex < numTriggers; ++hltIndex) {
0199     if (trigNames.triggerName(hltIndex) == triggerPath_ && hltresults->wasrun(hltIndex) && hltresults->accept(hltIndex))
0200       hasFired = true;
0201   }
0202 
0203   // access the trigger event
0204   edm::Handle<trigger::TriggerEvent> triggerEvent;
0205   e.getByToken(triggerEvent_, triggerEvent);
0206   if (triggerEvent.failedToGet()) {
0207     edm::LogError("DQMClientExample") << "invalid collection: TriggerEvent"
0208                                       << "\n";
0209     return;
0210   }
0211 
0212   reco::Particle *ele1_HLT = nullptr;
0213   int nEle_HLT = 0;
0214 
0215   size_t filterIndex = triggerEvent->filterIndex(triggerFilter_);
0216   trigger::TriggerObjectCollection triggerObjects = triggerEvent->getObjects();
0217   if (!(filterIndex >= triggerEvent->sizeFilters())) {
0218     const trigger::Keys &keys = triggerEvent->filterKeys(filterIndex);
0219     std::vector<reco::Particle> triggeredEle;
0220 
0221     for (size_t j = 0; j < keys.size(); ++j) {
0222       trigger::TriggerObject foundObject = triggerObjects[keys[j]];
0223       if (abs(foundObject.particle().pdgId()) != 11)
0224         continue;  // make sure that it is an electron
0225 
0226       triggeredEle.push_back(foundObject.particle());
0227       ++nEle_HLT;
0228     }
0229 
0230     if (!triggeredEle.empty())
0231       ele1_HLT = &(triggeredEle.at(0));
0232   }
0233 
0234   //-------------------------------
0235   //--- Fill the histos
0236   //-------------------------------
0237 
0238   // vertex
0239   h_vertex_number->Fill(vertex_number);
0240 
0241   // met
0242   h_pfMet->Fill(pfMETCollection->begin()->et());
0243 
0244   // multiplicities
0245   h_eMultiplicity->Fill(nEle);
0246   h_jMultiplicity->Fill(nJet);
0247   h_eMultiplicity_HLT->Fill(nEle_HLT);
0248 
0249   // leading not matched
0250   if (ele1) {
0251     h_ePt_leading->Fill(ele1->pt());
0252     h_eEta_leading->Fill(ele1->eta());
0253     h_ePhi_leading->Fill(ele1->phi());
0254   }
0255   if (ele1_HLT) {
0256     h_ePt_leading_HLT->Fill(ele1_HLT->pt());
0257     h_eEta_leading_HLT->Fill(ele1_HLT->eta());
0258     h_ePhi_leading_HLT->Fill(ele1_HLT->phi());
0259   }
0260   // leading Jet
0261   if (jet1) {
0262     h_jPt_leading->Fill(jet1->pt());
0263     h_jEta_leading->Fill(jet1->eta());
0264     h_jPhi_leading->Fill(jet1->phi());
0265   }
0266 
0267   // fill only when the trigger candidate mathes with the reco one
0268   if (ele1 && ele1_HLT && deltaR(*ele1_HLT, *ele1) < 0.3 && hasFired == true) {
0269     h_ePt_leading_matched->Fill(ele1->pt());
0270     h_eEta_leading_matched->Fill(ele1->eta());
0271     h_ePhi_leading_matched->Fill(ele1->phi());
0272 
0273     h_ePt_leading_HLT_matched->Fill(ele1_HLT->pt());
0274     h_eEta_leading_HLT_matched->Fill(ele1_HLT->eta());
0275     h_ePhi_leading_HLT_matched->Fill(ele1_HLT->phi());
0276 
0277     h_ePt_diff->Fill(ele1->pt() - ele1_HLT->pt());
0278   }
0279 }
0280 
0281 //
0282 // -------------------------------------- endRun
0283 // --------------------------------------------
0284 //
0285 //
0286 // -------------------------------------- book histograms
0287 // --------------------------------------------
0288 //
0289 void DQMExample_Step1::bookHistos(DQMStore::IBooker &ibooker_) {
0290   ibooker_.cd();
0291   ibooker_.setCurrentFolder("Physics/TopTest");
0292 
0293   h_vertex_number = ibooker_.book1D("Vertex_number", "Number of event vertices in collection", 40, -0.5, 39.5);
0294 
0295   h_pfMet = ibooker_.book1D("pfMet", "Pf Missing E_{T}; GeV", 20, 0.0, 100);
0296 
0297   h_eMultiplicity = ibooker_.book1D("NElectrons", "# of electrons per event", 10, 0., 10.);
0298   h_ePt_leading_matched = ibooker_.book1D("ElePt_leading_matched", "Pt of leading electron", 50, 0., 100.);
0299   h_eEta_leading_matched = ibooker_.book1D("EleEta_leading_matched", "Eta of leading electron", 50, -5., 5.);
0300   h_ePhi_leading_matched = ibooker_.book1D("ElePhi_leading_matched", "Phi of leading electron", 50, -3.5, 3.5);
0301 
0302   h_ePt_leading = ibooker_.book1D("ElePt_leading", "Pt of leading electron", 50, 0., 100.);
0303   h_eEta_leading = ibooker_.book1D("EleEta_leading", "Eta of leading electron", 50, -5., 5.);
0304   h_ePhi_leading = ibooker_.book1D("ElePhi_leading", "Phi of leading electron", 50, -3.5, 3.5);
0305 
0306   h_jMultiplicity = ibooker_.book1D("NJets", "# of electrons per event", 10, 0., 10.);
0307   h_jPt_leading = ibooker_.book1D("JetPt_leading", "Pt of leading Jet", 150, 0., 300.);
0308   h_jEta_leading = ibooker_.book1D("JetEta_leading", "Eta of leading Jet", 50, -5., 5.);
0309   h_jPhi_leading = ibooker_.book1D("JetPhi_leading", "Phi of leading Jet", 50, -3.5, 3.5);
0310 
0311   h_eMultiplicity_HLT = ibooker_.book1D("NElectrons_HLT", "# of electrons per event @HLT", 10, 0., 10.);
0312   h_ePt_leading_HLT = ibooker_.book1D("ElePt_leading_HLT", "Pt of leading electron @HLT", 50, 0., 100.);
0313   h_eEta_leading_HLT = ibooker_.book1D("EleEta_leading_HLT", "Eta of leading electron @HLT", 50, -5., 5.);
0314   h_ePhi_leading_HLT = ibooker_.book1D("ElePhi_leading_HLT", "Phi of leading electron @HLT", 50, -3.5, 3.5);
0315 
0316   h_ePt_leading_HLT_matched = ibooker_.book1D("ElePt_leading_HLT_matched", "Pt of leading electron @HLT", 50, 0., 100.);
0317   h_eEta_leading_HLT_matched =
0318       ibooker_.book1D("EleEta_leading_HLT_matched", "Eta of leading electron @HLT", 50, -5., 5.);
0319   h_ePhi_leading_HLT_matched =
0320       ibooker_.book1D("ElePhi_leading_HLT_matched", "Phi of leading electron @HLT", 50, -3.5, 3.5);
0321 
0322   h_ePt_diff = ibooker_.book1D("ElePt_diff_matched", "pT(RECO) - pT(HLT) for mathed candidates", 100, -10, 10.);
0323 
0324   ibooker_.cd();
0325 }
0326 
0327 //
0328 // -------------------------------------- functions
0329 // --------------------------------------------
0330 //
0331 double DQMExample_Step1::Distance(const reco::Candidate &c1, const reco::Candidate &c2) { return deltaR(c1, c2); }
0332 
0333 double DQMExample_Step1::DistancePhi(const reco::Candidate &c1, const reco::Candidate &c2) {
0334   return deltaPhi(c1.p4().phi(), c2.p4().phi());
0335 }
0336 
0337 // This always returns only a positive deltaPhi
0338 double DQMExample_Step1::calcDeltaPhi(double phi1, double phi2) {
0339   double deltaPhi = phi1 - phi2;
0340   if (deltaPhi < 0)
0341     deltaPhi = -deltaPhi;
0342   if (deltaPhi > 3.1415926) {
0343     deltaPhi = 2 * 3.1415926 - deltaPhi;
0344   }
0345   return deltaPhi;
0346 }
0347 
0348 //
0349 // -------------------------------------- electronID
0350 // --------------------------------------------
0351 //
0352 bool DQMExample_Step1::MediumEle(const edm::Event &iEvent,
0353                                  const edm::EventSetup &iESetup,
0354                                  const reco::GsfElectron &electron) {
0355   //********* CONVERSION TOOLS
0356   edm::Handle<reco::ConversionCollection> conversions_h;
0357   iEvent.getByToken(theConversionCollection_, conversions_h);
0358 
0359   bool isMediumEle = false;
0360 
0361   float pt = electron.pt();
0362   float eta = electron.eta();
0363 
0364   int isEB = electron.isEB();
0365   float sigmaIetaIeta = electron.sigmaIetaIeta();
0366   float DetaIn = electron.deltaEtaSuperClusterTrackAtVtx();
0367   float DphiIn = electron.deltaPhiSuperClusterTrackAtVtx();
0368   float HOverE = electron.hadronicOverEm();
0369   float ooemoop = (1.0 / electron.ecalEnergy() - electron.eSuperClusterOverP() / electron.ecalEnergy());
0370 
0371   int mishits = electron.gsfTrack()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS);
0372   int nAmbiguousGsfTracks = electron.ambiguousGsfTracksSize();
0373 
0374   reco::GsfTrackRef eleTrack = electron.gsfTrack();
0375   float dxy = eleTrack->dxy(PVPoint_);
0376   float dz = eleTrack->dz(PVPoint_);
0377 
0378   edm::Handle<reco::BeamSpot> BSHandle;
0379   iEvent.getByToken(theBSCollection_, BSHandle);
0380   const reco::BeamSpot BS = *BSHandle;
0381 
0382   bool isConverted = ConversionTools::hasMatchedConversion(electron, *conversions_h, BS.position());
0383 
0384   // default
0385   if ((pt > 12.) && (fabs(eta) < 2.5) &&
0386       (((isEB == 1) && (fabs(DetaIn) < 0.004)) || ((isEB == 0) && (fabs(DetaIn) < 0.007))) &&
0387       (((isEB == 1) && (fabs(DphiIn) < 0.060)) || ((isEB == 0) && (fabs(DphiIn) < 0.030))) &&
0388       (((isEB == 1) && (sigmaIetaIeta < 0.010)) || ((isEB == 0) && (sigmaIetaIeta < 0.030))) &&
0389       (((isEB == 1) && (HOverE < 0.120)) || ((isEB == 0) && (HOverE < 0.100))) &&
0390       (((isEB == 1) && (fabs(ooemoop) < 0.050)) || ((isEB == 0) && (fabs(ooemoop) < 0.050))) &&
0391       (((isEB == 1) && (fabs(dxy) < 0.020)) || ((isEB == 0) && (fabs(dxy) < 0.020))) &&
0392       (((isEB == 1) && (fabs(dz) < 0.100)) || ((isEB == 0) && (fabs(dz) < 0.100))) &&
0393       (((isEB == 1) && (!isConverted)) || ((isEB == 0) && (!isConverted))) && (mishits == 0) &&
0394       (nAmbiguousGsfTracks == 0))
0395     isMediumEle = true;
0396 
0397   return isMediumEle;
0398 }