Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:08:03 

0001 #include "DQM/Physics/src/EwkMuDQM.h"
0002 
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005 #include "DataFormats/Common/interface/Handle.h"
0006 
0007 #include "FWCore/ServiceRegistry/interface/Service.h"
0008 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0009 
0010 #include "DQMServices/Core/interface/DQMStore.h"
0011 
0012 #include "DataFormats/TrackReco/interface/Track.h"
0013 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0014 #include "DataFormats/VertexReco/interface/Vertex.h"
0015 
0016 #include "DataFormats/MuonReco/interface/Muon.h"
0017 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0018 #include "DataFormats/METReco/interface/MET.h"
0019 #include "DataFormats/JetReco/interface/Jet.h"
0020 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0021 
0022 #include "DataFormats/GeometryVector/interface/Phi.h"
0023 
0024 #include "FWCore/Common/interface/TriggerNames.h"
0025 #include "FWCore/Framework/interface/Event.h"
0026 #include "DataFormats/Common/interface/TriggerResults.h"
0027 
0028 #include "DataFormats/Common/interface/View.h"
0029 
0030 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0031 
0032 using namespace edm;
0033 using namespace std;
0034 using namespace reco;
0035 
0036 EwkMuDQM::EwkMuDQM(const ParameterSet& cfg)
0037     :  // Input collections
0038       metTag_(cfg.getUntrackedParameter<edm::InputTag>("METTag", edm::InputTag("pfmet"))),
0039       jetTag_(cfg.getUntrackedParameter<edm::InputTag>("JetTag", edm::InputTag("ak4PFJets"))),
0040       //trigTag_(consumes<edm::TriggerResults>(
0041       // cfg.getUntrackedParameter<edm::InputTag>(
0042       // "TrigTag", edm::InputTag("TriggerResults::HLT")))),
0043       muonTag_(consumes<edm::View<reco::Muon> >(
0044           cfg.getUntrackedParameter<edm::InputTag>("MuonTag", edm::InputTag("muons")))),
0045       metToken_(
0046           consumes<edm::View<reco::MET> >(cfg.getUntrackedParameter<edm::InputTag>("METTag", edm::InputTag("pfmet")))),
0047       jetToken_(consumes<edm::View<reco::Jet> >(
0048           cfg.getUntrackedParameter<edm::InputTag>("JetTag", edm::InputTag("ak4PFJets")))),
0049       phoTag_(consumes<edm::View<reco::Photon> >(
0050           cfg.getUntrackedParameter<edm::InputTag>("phoTag", edm::InputTag("photons")))),
0051       vertexTag_(consumes<edm::View<reco::Vertex> >(
0052           cfg.getUntrackedParameter<edm::InputTag>("VertexTag", edm::InputTag("offlinePrimaryVertices")))),
0053       beamSpotTag_(consumes<reco::BeamSpot>(
0054           cfg.getUntrackedParameter<edm::InputTag>("beamSpotTag", edm::InputTag("offlineBeamSpot")))),
0055       // trigPathNames_(cfg.getUntrackedParameter<std::vector<std::string> >(
0056       //   "TrigPathNames")),
0057 
0058       // Muon quality cuts
0059       isAlsoTrackerMuon_(cfg.getUntrackedParameter<bool>("IsAlsoTrackerMuon", true)),  // Glb muon also tracker muon
0060       dxyCut_(cfg.getUntrackedParameter<double>("DxyCut", 0.2)),                       // dxy < 0.2 cm
0061       normalizedChi2Cut_(
0062           cfg.getUntrackedParameter<double>("NormalizedChi2Cut", 10.)),  // chi2/ndof (of global fit) <10.0
0063       trackerHitsCut_(cfg.getUntrackedParameter<int>("TrackerHitsCut",
0064                                                      11)),               // Tracker Hits >10
0065       pixelHitsCut_(cfg.getUntrackedParameter<int>("PixelHitsCut", 1)),  // Pixel Hits >0
0066       muonHitsCut_(cfg.getUntrackedParameter<int>("MuonHitsCut",
0067                                                   1)),                 // Valid Muon Hits >0
0068       nMatchesCut_(cfg.getUntrackedParameter<int>("NMatchesCut", 2)),  // At least 2 Chambers with matches
0069 
0070       // W-boson cuts
0071       isRelativeIso_(cfg.getUntrackedParameter<bool>("IsRelativeIso", true)),
0072       isCombinedIso_(cfg.getUntrackedParameter<bool>("IsCombinedIso", false)),
0073       isoCut03_(cfg.getUntrackedParameter<double>("IsoCut03", 0.1)),
0074       acopCut_(cfg.getUntrackedParameter<double>("AcopCut", 999.)),
0075       metMin_(cfg.getUntrackedParameter<double>("MetMin", -999999.)),
0076       metMax_(cfg.getUntrackedParameter<double>("MetMax", 999999.)),
0077       mtMin_(cfg.getUntrackedParameter<double>("MtMin", 50.)),
0078       mtMax_(cfg.getUntrackedParameter<double>("MtMax", 200.)),
0079       ptCut_(cfg.getUntrackedParameter<double>("PtCut", 20.)),
0080       etaCut_(cfg.getUntrackedParameter<double>("EtaCut", 2.4)),
0081 
0082       // Z rejection
0083       ptThrForZ1_(cfg.getUntrackedParameter<double>("PtThrForZ1", 20.)),
0084       ptThrForZ2_(cfg.getUntrackedParameter<double>("PtThrForZ2", 10.)),
0085 
0086       // Z selection
0087       dimuonMassMin_(cfg.getUntrackedParameter<double>("dimuonMassMin", 80.)),
0088       dimuonMassMax_(cfg.getUntrackedParameter<double>("dimuonMassMax", 120.)),
0089 
0090       // Top rejection
0091       eJetMin_(cfg.getUntrackedParameter<double>("EJetMin", 999999.)),
0092       nJetMax_(cfg.getUntrackedParameter<int>("NJetMax", 999999)),
0093 
0094       // Photon cuts
0095       ptThrForPhoton_(cfg.getUntrackedParameter<double>("ptThrForPhoton", 5.)),
0096       nPhoMax_(cfg.getUntrackedParameter<int>("nPhoMax", 999999)),
0097       hltPrescaleProvider_(cfg, consumesCollector(), *this) {
0098   isValidHltConfig_ = false;
0099 }
0100 
0101 void EwkMuDQM::dqmBeginRun(const Run& iRun, const EventSetup& iSet) {
0102   nall = 0;
0103   nsel = 0;
0104   nz = 0;
0105 
0106   nrec = 0;
0107   niso = 0;
0108   nhlt = 0;
0109   nmet = 0;
0110 
0111   // passed as parameter to HLTConfigProvider::init(), not yet used
0112   bool isConfigChanged = false;
0113   // isValidHltConfig_ used to short-circuit analyze() in case of problems
0114   isValidHltConfig_ = hltPrescaleProvider_.init(iRun, iSet, "HLT", isConfigChanged);
0115 }
0116 
0117 void EwkMuDQM::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const&, edm::EventSetup const&) {
0118   ibooker.setCurrentFolder("Physics/EwkMuDQM");
0119 
0120   char chtitle[256] = "";
0121 
0122   pt_before_ = ibooker.book1D("PT_BEFORECUTS", "Muon transverse momentum (global muon) [GeV]", 100, 0., 100.);
0123   pt_after_ = ibooker.book1D("PT_AFTERWCUTS", "Muon transverse momentum (global muon) [GeV]", 100, 0., 100.);
0124 
0125   eta_before_ = ibooker.book1D("ETA_BEFORECUTS", "Muon pseudo-rapidity", 50, -2.5, 2.5);
0126   eta_after_ = ibooker.book1D("ETA_AFTERWCUTS", "Muon pseudo-rapidity", 50, -2.5, 2.5);
0127 
0128   dxy_before_ = ibooker.book1D("DXY_BEFORECUTS", "Muon transverse distance to beam spot [cm]", 100, -0.5, 0.5);
0129   dxy_after_ = ibooker.book1D("DXY_AFTERWCUTS", "Muon transverse distance to beam spot [cm]", 100, -0.5, 0.5);
0130 
0131   goodewkmuon_before_ = ibooker.book1D("GOODEWKMUON_BEFORECUTS", "Quality-muon flag", 2, -0.5, 1.5);
0132   goodewkmuon_after_ = ibooker.book1D("GOODEWKMUON_AFTERWCUTS", "Quality-muon flag", 2, -0.5, 1.5);
0133 
0134   if (isRelativeIso_) {
0135     if (isCombinedIso_) {
0136       iso_before_ = ibooker.book1D("ISO_BEFORECUTS", "Relative (combined) isolation variable", 100, 0., 1.);
0137       iso_after_ = ibooker.book1D("ISO_AFTERWCUTS", "Relative (combined) isolation variable", 100, 0., 1.);
0138     } else {
0139       iso_before_ = ibooker.book1D("ISO_BEFORECUTS", "Relative (tracker) isolation variable", 100, 0., 1.);
0140       iso_after_ = ibooker.book1D("ISO_AFTERWCUTS", "Relative (tracker) isolation variable", 100, 0., 1.);
0141     }
0142   } else {
0143     if (isCombinedIso_) {
0144       iso_before_ = ibooker.book1D("ISO_BEFORECUTS", "Absolute (combined) isolation variable [GeV]", 100, 0., 20.);
0145       iso_after_ = ibooker.book1D("ISO_AFTERWCUTS", "Absolute (combined) isolation variable [GeV]", 100, 0., 20.);
0146     } else {
0147       iso_before_ = ibooker.book1D("ISO_BEFORECUTS", "Absolute (tracker) isolation variable [GeV]", 100, 0., 20.);
0148       iso_after_ = ibooker.book1D("ISO_AFTERWCUTS", "Absolute (tracker) isolation variable [GeV]", 100, 0., 20.);
0149     }
0150   }
0151 
0152   /*  trig_before_ = ibooker.book1D("TRIG_BEFORECUTS",
0153       "Trigger response (boolean of muon triggers)", 2, -0.5, 1.5);
0154   trig_after_ = ibooker.book1D("TRIG_AFTERWCUTS",
0155       "Trigger response (boolean of muon triggers)", 2, -0.5, 1.5);
0156 */
0157   snprintf(chtitle, 255, "Transverse mass (%s) [GeV]", metTag_.label().data());
0158   mt_before_ = ibooker.book1D("MT_BEFORECUTS", chtitle, 150, 0., 300.);
0159   mt_after_ = ibooker.book1D("MT_AFTERWCUTS", chtitle, 150, 0., 300.);
0160 
0161   snprintf(chtitle, 255, "Missing transverse energy (%s) [GeV]", metTag_.label().data());
0162   met_before_ = ibooker.book1D("MET_BEFORECUTS", chtitle, 100, 0., 200.);
0163   met_after_ = ibooker.book1D("MET_AFTERWCUTS", chtitle, 100, 0., 200.);
0164   met_afterZ_ = ibooker.book1D("MET_AFTERZCUTS", chtitle, 100, 0., 200.);
0165 
0166   snprintf(chtitle, 255, "MU-MET (%s) acoplanarity", metTag_.label().data());
0167   acop_before_ = ibooker.book1D("ACOP_BEFORECUTS", chtitle, 50, 0., M_PI);
0168   acop_after_ = ibooker.book1D("ACOP_AFTERWCUTS", chtitle, 50, 0., M_PI);
0169 
0170   snprintf(chtitle, 255, "Z selection: muons above %.2f GeV", ptThrForZ1_);
0171   n_zselPt1thr_ = ibooker.book1D("NZSELPT1THR", chtitle, 10, -0.5, 9.5);
0172   snprintf(chtitle, 255, "Z selection: muons above %.2f GeV", ptThrForZ2_);
0173   n_zselPt2thr_ = ibooker.book1D("NZSELPT2THR", chtitle, 10, -0.5, 9.5);
0174 
0175   snprintf(chtitle, 255, "Number of jets (%s) above %.2f GeV", jetTag_.label().data(), eJetMin_);
0176   njets_before_ = ibooker.book1D("NJETS_BEFORECUTS", chtitle, 16, -0.5, 15.5);
0177   njets_after_ = ibooker.book1D("NJETS_AFTERWCUTS", chtitle, 16, -0.5, 15.5);
0178   njets_afterZ_ = ibooker.book1D("NJETS_AFTERZCUTS", chtitle, 16, -0.5, 15.5);
0179 
0180   leadingjet_pt_before_ = ibooker.book1D("LEADINGJET_PT_BEFORECUTS", "Leading Jet transverse momentum", 300, 0., 300.);
0181   leadingjet_pt_after_ = ibooker.book1D("LEADINGJET_PT_AFTERWCUTS", "Leading Jet transverse momentum", 300, 0., 300.);
0182   leadingjet_pt_afterZ_ = ibooker.book1D("LEADINGJET_PT_AFTERZCUTS", "Leading Jet transverse momentum", 300, 0., 300.);
0183 
0184   leadingjet_eta_before_ = ibooker.book1D("LEADINGJET_ETA_BEFORECUTS", "Leading Jet pseudo-rapidity", 50, -2.5, 2.5);
0185   leadingjet_eta_after_ = ibooker.book1D("LEADINGJET_ETA_AFTERWCUTS", "Leading Jet pseudo-rapidity", 50, -2.5, 2.5);
0186   leadingjet_eta_afterZ_ = ibooker.book1D("LEADINGJET_ETA_AFTERZCUTS", "Leading Jet pseudo-rapidity", 50, -2.5, 2.5);
0187 
0188   ptDiffPM_before_ = ibooker.book1D("PTDIFFPM_BEFORE_CUTS", "pt(Muon+)-pt(Muon-) after Z cuts [GeV]", 200, -100., 100.);
0189   ptDiffPM_afterZ_ = ibooker.book1D("PTDIFFPM_AFTERZ_CUTS", "pt(Muon+)-pt(Muon-) after Z cuts [GeV]", 200, -100., 100.);
0190 
0191   /**\ For Z-boson events  */
0192 
0193   pt1_afterZ_ = ibooker.book1D("PT1_AFTERZCUTS", "Muon transverse momentum (global muon) [GeV]", 100, 0., 100.);
0194   eta1_afterZ_ = ibooker.book1D("ETA1_AFTERZCUTS", "Muon pseudo-rapidity", 50, -2.5, 2.5);
0195   dxy1_afterZ_ = ibooker.book1D("DXY1_AFTERZCUTS", "Muon transverse distance to beam spot [cm]", 100, -0.5, 0.5);
0196   goodewkmuon1_afterZ_ = ibooker.book1D("GOODEWKMUON1_AFTERZCUTS", "Quality-muon flag", 2, -0.5, 1.5);
0197 
0198   if (isRelativeIso_) {
0199     if (isCombinedIso_) {
0200       iso1_afterZ_ = ibooker.book1D("ISO1_AFTERZCUTS", "Relative (combined) isolation variable", 100, 0., 1.);
0201       iso2_afterZ_ = ibooker.book1D("ISO2_AFTERZCUTS", "Relative (combined) isolation variable", 100, 0., 1.);
0202     } else {
0203       iso1_afterZ_ = ibooker.book1D("ISO1_AFTERZCUTS", "Relative (tracker) isolation variable", 100, 0., 1.);
0204       iso2_afterZ_ = ibooker.book1D("ISO2_AFTERZCUTS", "Relative (tracker) isolation variable", 100, 0., 1.);
0205     }
0206   } else {
0207     if (isCombinedIso_) {
0208       iso1_afterZ_ = ibooker.book1D("ISO1_AFTERZCUTS", "Absolute (combined) isolation variable [GeV]", 100, 0., 20.);
0209       iso2_afterZ_ = ibooker.book1D("ISO2_AFTERZCUTS", "Absolute (combined) isolation variable [GeV]", 100, 0., 20.);
0210     } else {
0211       iso1_afterZ_ = ibooker.book1D("ISO1_AFTERZCUTS", "Absolute (tracker) isolation variable [GeV]", 100, 0., 20.);
0212       iso2_afterZ_ = ibooker.book1D("ISO2_AFTERZCUTS", "Absolute (tracker) isolation variable [GeV]", 100, 0., 20.);
0213     }
0214   }
0215 
0216   pt2_afterZ_ = ibooker.book1D("PT2_AFTERZCUTS", "Muon transverse momentum (global muon) [GeV]", 100, 0., 100.);
0217   eta2_afterZ_ = ibooker.book1D("ETA2_AFTERZCUTS", "Muon pseudo-rapidity", 50, -2.5, 2.5);
0218   dxy2_afterZ_ = ibooker.book1D("DXY2_AFTERZCUTS", "Muon transverse distance to beam spot [cm]", 100, -0.5, 0.5);
0219   goodewkmuon2_afterZ_ = ibooker.book1D("GOODEWKMUON2_AFTERZCUTS", "Quality-muon flag", 2, -0.5, 1.5);
0220   /*  ztrig_afterZ_ = ibooker.book1D("ZTRIG_AFTERZCUTS",
0221       "Trigger response (boolean of muon triggers)", 2, -0.5, 1.5);
0222  */
0223   dimuonmass_before_ = ibooker.book1D("DIMUONMASS_BEFORECUTS", "DiMuonMass (2 globals)", 100, 0, 200);
0224   dimuonmass_afterZ_ = ibooker.book1D("DIMUONMASS_AFTERZCUTS", "DiMuonMass (2 globals)", 100, 0, 200);
0225   npvs_before_ = ibooker.book1D("NPVs_BEFORECUTS", "Number of Valid Primary Vertices", 51, -0.5, 50.5);
0226   npvs_after_ = ibooker.book1D("NPVs_AFTERWCUTS", "Number of Valid Primary Vertices", 51, -0.5, 50.5);
0227   npvs_afterZ_ = ibooker.book1D("NPVs_AFTERZCUTS", "Number of Valid Primary Vertices", 51, -0.5, 50.5);
0228   muoncharge_before_ = ibooker.book1D("MUONCHARGE_BEFORECUTS", "Muon Charge", 3, -1.5, 1.5);
0229   muoncharge_after_ = ibooker.book1D("MUONCHARGE_AFTERWCUTS", "Muon Charge", 3, -1.5, 1.5);
0230   muoncharge_afterZ_ = ibooker.book1D("MUONCHARGE_AFTERZCUTS", "Muon Charge", 3, -1.5, 1.5);
0231 
0232   // Adding these to replace the NZ ones (more useful, since they are more
0233   // general?)
0234   nmuons_ = ibooker.book1D("NMuons", "Number of muons in the event", 10, -0.5, 9.5);
0235   ngoodmuons_ = ibooker.book1D("NGoodMuons", "Number of muons passing the quality criteria", 10, -0.5, 9.5);
0236 
0237   nph_ = ibooker.book1D("nph", "Number of photons in the event", 20, 0., 20.);
0238   phPt_ = ibooker.book1D("phPt", "Photon transverse momentum [GeV]", 100, 0., 1000.);
0239   snprintf(chtitle, 255, "Photon pseudorapidity (pT>%4.1f)", ptThrForPhoton_);
0240   phEta_ = ibooker.book1D("phEta", chtitle, 100, -2.5, 2.5);
0241 }
0242 
0243 void EwkMuDQM::analyze(const Event& ev, const EventSetup& iSet) {
0244   // Muon collection
0245   Handle<View<Muon> > muonCollection;
0246   if (!ev.getByToken(muonTag_, muonCollection)) {
0247     // LogWarning("") << ">>> Muon collection does not exist !!!";
0248     return;
0249   }
0250   unsigned int muonCollectionSize = muonCollection->size();
0251 
0252   // Beam spot
0253   Handle<reco::BeamSpot> beamSpotHandle;
0254   if (!ev.getByToken(beamSpotTag_, beamSpotHandle)) {
0255     // LogWarning("") << ">>> No beam spot found !!!";
0256     return;
0257   }
0258 
0259   // Loop to reject/control Z->mumu is done separately
0260   unsigned int nmuonsForZ1 = 0;
0261   unsigned int nmuonsForZ2 = 0;
0262   bool cosmic = false;
0263   for (unsigned int i = 0; i < muonCollectionSize; i++) {
0264     const Muon& mu = muonCollection->at(i);
0265     if (!mu.isGlobalMuon())
0266       continue;
0267     double pt = mu.pt();
0268     double dxy = mu.innerTrack()->dxy(beamSpotHandle->position());
0269 
0270     if (fabs(dxy) > 1) {
0271       cosmic = true;
0272       break;
0273     }
0274 
0275     if (pt > ptThrForZ1_)
0276       nmuonsForZ1++;
0277     if (pt > ptThrForZ2_)
0278       nmuonsForZ2++;
0279 
0280     for (unsigned int j = i + 1; j < muonCollectionSize; j++) {
0281       const Muon& mu2 = muonCollection->at(j);
0282       if (mu2.isGlobalMuon() && (mu.charge() * mu2.charge() == -1)) {
0283         const math::XYZTLorentzVector ZRecoGlb(
0284             mu.px() + mu2.px(), mu.py() + mu2.py(), mu.pz() + mu2.pz(), mu.p() + mu2.p());
0285         dimuonmass_before_->Fill(ZRecoGlb.mass());
0286         if (mu.charge() > 0) {
0287           ptDiffPM_before_->Fill(mu.pt() - mu2.pt());
0288         } else {
0289           ptDiffPM_before_->Fill(mu2.pt() - mu.pt());
0290         }
0291       }
0292     }
0293   }
0294   if (cosmic)
0295     return;
0296 
0297   LogTrace("") << "> Z rejection: muons above " << ptThrForZ1_ << " [GeV]: " << nmuonsForZ1;
0298   LogTrace("") << "> Z rejection: muons above " << ptThrForZ2_ << " [GeV]: " << nmuonsForZ2;
0299 
0300   // MET
0301   Handle<View<MET> > metCollection;
0302   if (!ev.getByToken(metToken_, metCollection)) {
0303     // LogWarning("") << ">>> MET collection does not exist !!!";
0304     return;
0305   }
0306   const MET& met = metCollection->at(0);
0307   double met_et = met.pt();
0308   LogTrace("") << ">>> MET, MET_px, MET_py: " << met_et << ", " << met.px() << ", " << met.py() << " [GeV]";
0309   met_before_->Fill(met_et);
0310 
0311   // Vertices in the event
0312   Handle<View<reco::Vertex> > vertexCollection;
0313   if (!ev.getByToken(vertexTag_, vertexCollection)) {
0314     LogError("") << ">>> Vertex collection does not exist !!!";
0315     return;
0316   }
0317   unsigned int vertexCollectionSize = vertexCollection->size();
0318 
0319   int nvvertex = 0;
0320   for (unsigned int i = 0; i < vertexCollectionSize; i++) {
0321     const Vertex& vertex = vertexCollection->at(i);
0322     if (vertex.isValid())
0323       nvvertex++;
0324   }
0325 
0326   npvs_before_->Fill(nvvertex);
0327 
0328   // bool trigger_fired = false;
0329   //Handle<TriggerResults> triggerResults;
0330   // if (!ev.getByToken(trigTag_, triggerResults)) {
0331   // LogWarning("") << ">>> TRIGGER collection does not exist !!!";
0332   // return;
0333   // }
0334   // const edm::TriggerNames& trigNames = ev.triggerNames(*triggerResults);
0335   //  LogWarning("")<<"Loop over triggers";
0336 
0337   //HLTConfigProvider const&  hltConfigProvider = hltPrescaleProvider_.hltConfigProvider();
0338 
0339   /*  change faulty logic of triggering
0340   for (unsigned int i=0; i<triggerResults->size(); i++)
0341   {
0342           const std::string trigName = trigNames.triggerName(i);
0343 
0344           bool found=false;
0345           for(unsigned int index=0; index<trigPathNames_.size() && found==false;
0346   index++) {
0347                size_t trigPath = trigName.find(trigPathNames_[index]); // 0 if
0348   found, pos if not
0349                if (trigPath==0) found=true;
0350           }
0351           if(!found) {continue;}
0352 
0353           bool prescaled=false;
0354           for (unsigned int ps= 0; ps<  hltConfigProvider.prescaleSize();
0355   ps++){
0356               const unsigned int prescaleValue =
0357   hltConfigProvider.prescaleValue(ps, trigName) ;
0358               if (prescaleValue != 1) prescaled =true;
0359           }
0360 
0361           if( triggerResults->accept(i) && !prescaled){   trigger_fired=true;}
0362                     // LogWarning("")<<"TrigNo: "<<i<<"  "<<found<<"
0363   "<<trigName<<" ---> FIRED";}
0364   }
0365   */
0366 
0367   // get the prescale set for this event
0368   const int prescaleSet = hltPrescaleProvider_.prescaleSet(ev, iSet);
0369   if (prescaleSet == -1) {
0370     LogTrace("") << "Failed to determine prescaleSet\n";
0371     // std::cout << "Failed to determine prescaleSet. Check the GlobalTag in
0372     // cfg\n";
0373     return;
0374   }
0375 
0376   // for (unsigned int i = 0;
0377   //    (i < triggerResults->size()) && (trigger_fired == false); i++) {
0378   // skip trigger, if it did not fire
0379   //if (!triggerResults->accept(i)) continue;
0380 
0381   // skip trigger, if it is not on our list
0382   //bool found = false;
0383   //const std::string trigName = trigNames.triggerName(i);
0384   //for (unsigned int index = 0;
0385   //   index < trigPathNames_.size() && found == false; index++) {
0386   // if (trigName.find(trigPathNames_.at(index)) == 0) found = true;
0387   // }
0388   // if (!found) continue;
0389 
0390   // skip trigger, if it is prescaled
0391   /* if (prescaleSet != -1) {
0392       if (hltConfigProvider.prescaleValue(prescaleSet, trigName) != 1)
0393         continue;
0394     } else {
0395       // prescaleSet is not known.
0396       // This branch is not needed, if prescaleSet=-1 forces to skip event
0397       int prescaled = 0;
0398       for (unsigned int ps = 0;
0399            !prescaled && (ps < hltConfigProvider.prescaleSize()); ++ps) {
0400         if (hltConfigProvider.prescaleValue(ps, trigName) != 1) {
0401           prescaled = 1;
0402         }
0403       }
0404       if (prescaled) {
0405         // std::cout << "trigger prescaled\n";
0406         continue;
0407       }
0408     }
0409 */
0410   // std::cout << "found unprescaled trigger that fired: " << trigName <<
0411   // "\n";
0412   // trigger_fired = true;
0413   // }
0414   // if (trigger_fired) std::cout << "\n\tGot Trigger\n";
0415 
0416   // trig_before_->Fill(trigger_fired);
0417 
0418   // Jet collection
0419   Handle<View<Jet> > jetCollection;
0420   if (!ev.getByToken(jetToken_, jetCollection)) {
0421     // LogError("") << ">>> JET collection does not exist !!!";
0422     return;
0423   }
0424   unsigned int jetCollectionSize = jetCollection->size();
0425   int njets = 0;
0426   int LEADJET = -1;
0427   double max_pt = 0;
0428   for (unsigned int i = 0; i < jetCollectionSize; i++) {
0429     const Jet& jet = jetCollection->at(i);
0430     double minDistance = 99999;  // This is in order to use PFJets
0431     for (unsigned int j = 0; j < muonCollectionSize; j++) {
0432       const Muon& mu = muonCollection->at(j);
0433       double distance =
0434           sqrt((mu.eta() - jet.eta()) * (mu.eta() - jet.eta()) + (mu.phi() - jet.phi()) * (mu.phi() - jet.phi()));
0435       if (minDistance > distance)
0436         minDistance = distance;
0437     }
0438     if (minDistance < 0.3)
0439       continue;  // 0.3 is the isolation cone around the muon
0440     if (jet.et() > max_pt) {
0441       LEADJET = i;
0442       max_pt = jet.et();
0443     }
0444     if (jet.et() > eJetMin_) {
0445       njets++;
0446     }
0447   }
0448 
0449   LogTrace("") << ">>> Total number of jets: " << jetCollectionSize;
0450   LogTrace("") << ">>> Number of jets above " << eJetMin_ << " [GeV]: " << njets;
0451   njets_before_->Fill(njets);
0452   double lead_jet_pt = -1;
0453   double lead_jet_eta = -100;
0454   if (LEADJET != -1) {
0455     const Jet& leadJet = jetCollection->at(LEADJET);
0456     leadingjet_pt_before_->Fill(leadJet.pt());
0457     leadingjet_eta_before_->Fill(leadJet.eta());
0458     lead_jet_pt = leadJet.pt();
0459     lead_jet_eta = leadJet.eta();
0460   }
0461   // Photon Collection
0462   Handle<View<Photon> > photonCollection;
0463   if (!ev.getByToken(phoTag_, photonCollection)) {
0464     // LogError("")
0465     return;
0466   }
0467   unsigned int ngam = 0;
0468 
0469   for (unsigned int i = 0; i < photonCollection->size(); i++) {
0470     const Photon& ph = photonCollection->at(i);
0471     double photonPt = ph.pt();
0472     if (photonPt > ptThrForPhoton_) {
0473       ngam++;
0474       phEta_->Fill(ph.eta());
0475     }
0476     phPt_->Fill(photonPt);
0477   }
0478   nph_->Fill(ngam);
0479   LogTrace("") << " >>> N photons " << ngam << std::endl;
0480 
0481   nmuons_->Fill(muonCollectionSize);
0482 
0483   // Start counting
0484   nall++;
0485 
0486   // Histograms per event should be done only once, so keep track of them
0487   // bool hlt_hist_done = false;
0488   // bool zhlt_hist_done = false;
0489   bool zjets_hist_done = false;
0490   bool zfullsel_hist_done = false;
0491   bool met_hist_done = false;
0492   bool njets_hist_done = false;
0493   bool wfullsel_hist_done = false;
0494 
0495   // Central W->mu nu selection criteria
0496   const int NFLAGS = 10;
0497   bool muon_sel[NFLAGS];
0498   const int NFLAGSZ = 12;
0499   bool zmuon_sel[NFLAGSZ];
0500   bool muon4Z = false;
0501 
0502   double number_of_goodMuons = 0;
0503 
0504   for (unsigned int i = 0; i < muonCollectionSize; i++) {
0505     for (int j = 0; j < NFLAGS; ++j) {
0506       muon_sel[j] = false;
0507     }
0508 
0509     const Muon& mu = muonCollection->at(i);
0510     if (!mu.isGlobalMuon())
0511       continue;
0512     if (mu.globalTrack().isNull())
0513       continue;
0514     if (mu.innerTrack().isNull())
0515       continue;
0516 
0517     LogTrace("") << "> Wsel: processing muon number " << i << "...";
0518     reco::TrackRef gm = mu.globalTrack();
0519     reco::TrackRef tk = mu.innerTrack();
0520 
0521     // Pt,eta cuts
0522     double pt = mu.pt();
0523     double eta = mu.eta();
0524     LogTrace("") << "\t... pt, eta: " << pt << " [GeV], " << eta;
0525     ;
0526     if (pt > ptCut_)
0527       muon_sel[0] = true;
0528     if (fabs(eta) < etaCut_)
0529       muon_sel[1] = true;
0530 
0531     double charge = mu.charge();
0532 
0533     // d0, chi2, nhits quality cuts
0534     double dxy = gm->dxy(beamSpotHandle->position());
0535     double normalizedChi2 = gm->normalizedChi2();
0536     double trackerHits = tk->hitPattern().numberOfValidTrackerHits();
0537     int pixelHits = tk->hitPattern().numberOfValidPixelHits();
0538     int muonHits = gm->hitPattern().numberOfValidMuonHits();
0539     int nMatches = mu.numberOfMatches();
0540 
0541     LogTrace("") << "\t... dxy, normalizedChi2, trackerHits, isTrackerMuon?: " << dxy << " [cm], " << normalizedChi2
0542                  << ", " << trackerHits << ", " << mu.isTrackerMuon();
0543     if (fabs(dxy) < dxyCut_)
0544       muon_sel[2] = true;
0545 
0546     bool quality = true;
0547 
0548     if (normalizedChi2 > normalizedChi2Cut_)
0549       quality = false;
0550     if (trackerHits < trackerHitsCut_)
0551       quality = false;
0552     if (pixelHits < pixelHitsCut_)
0553       quality = false;
0554     if (muonHits < muonHitsCut_)
0555       quality = false;
0556     ;
0557     if (!mu.isTrackerMuon())
0558       quality = false;
0559     if (nMatches < nMatchesCut_)
0560       quality = false;
0561     muon_sel[3] = quality;
0562     if (quality)
0563       number_of_goodMuons++;
0564 
0565     pt_before_->Fill(pt);
0566     eta_before_->Fill(eta);
0567     dxy_before_->Fill(dxy);
0568     muoncharge_before_->Fill(charge);
0569     goodewkmuon_before_->Fill(quality);
0570 
0571     // Charge asymmetry
0572     // if (quality) {
0573     //  if (charge>0) ptPlus_before_->Fill(pt);
0574     //  if (charge<0) ptMinus_before_->Fill(pt);
0575     //}
0576 
0577     // Isolation cuts
0578     double isovar = mu.isolationR03().sumPt;
0579     if (isCombinedIso_) {
0580       isovar += mu.isolationR03().emEt;
0581       isovar += mu.isolationR03().hadEt;
0582     }
0583     if (isRelativeIso_)
0584       isovar /= pt;
0585     if (isovar < isoCut03_)
0586       muon_sel[4] = true;
0587 
0588     LogTrace("") << "\t... isolation value" << isovar << ", isolated? " << muon_sel[6];
0589     iso_before_->Fill(isovar);
0590 
0591     // HLT (not mtched to muon for the time being)
0592     // if (trigger_fired) muon_sel[5] = true;
0593 
0594     // For Z:
0595     if (pt > ptThrForZ1_ && fabs(eta) < etaCut_ && fabs(dxy) < dxyCut_ && quality && isovar < isoCut03_) {
0596       muon4Z = true;
0597     }
0598 
0599     // MET/MT cuts
0600     double w_et = met_et + mu.pt();
0601     double w_px = met.px() + mu.px();
0602     double w_py = met.py() + mu.py();
0603 
0604     double massT = w_et * w_et - w_px * w_px - w_py * w_py;
0605     massT = (massT > 0) ? sqrt(massT) : 0;
0606 
0607     LogTrace("") << "\t... W mass, W_et, W_px, W_py: " << massT << ", " << w_et << ", " << w_px << ", " << w_py
0608                  << " [GeV]";
0609     if (massT > mtMin_ && massT < mtMax_)
0610       muon_sel[5] = true;
0611     mt_before_->Fill(massT);
0612     if (met_et > metMin_ && met_et < metMax_)
0613       muon_sel[6] = true;
0614 
0615     // Acoplanarity cuts
0616     Geom::Phi<double> deltaphi(mu.phi() - atan2(met.py(), met.px()));
0617     double acop = deltaphi.value();
0618     if (acop < 0)
0619       acop = -acop;
0620     acop = M_PI - acop;
0621     LogTrace("") << "\t... acoplanarity: " << acop;
0622     if (acop < acopCut_)
0623       muon_sel[7] = true;
0624     acop_before_->Fill(acop);
0625 
0626     // Remaining flags (from global event information)
0627     if (nmuonsForZ1 < 1 || nmuonsForZ2 < 2)
0628       muon_sel[8] = true;
0629     if (njets <= nJetMax_)
0630       muon_sel[9] = true;
0631 
0632     // Collect necessary flags "per muon"
0633     int flags_passed = 0;
0634     for (int j = 0; j < NFLAGS; ++j) {
0635       if (muon_sel[j])
0636         flags_passed += 1;
0637     }
0638 
0639     // Do N-1 histograms now (and only once for global event quantities)
0640     if (flags_passed >= (NFLAGS - 1)) {
0641       if (!muon_sel[0] || flags_passed == NFLAGS)
0642         pt_after_->Fill(pt);
0643       if (!muon_sel[1] || flags_passed == NFLAGS)
0644         eta_after_->Fill(eta);
0645       if (!muon_sel[2] || flags_passed == NFLAGS)
0646         dxy_after_->Fill(dxy);
0647       if (!muon_sel[3] || flags_passed == NFLAGS)
0648         goodewkmuon_after_->Fill(quality);
0649       if (!muon_sel[4] || flags_passed == NFLAGS)
0650         iso_after_->Fill(isovar);
0651       // if (!muon_sel[5] || flags_passed == NFLAGS)
0652       //   if (!hlt_hist_done) trig_after_->Fill(trigger_fired);
0653       //hlt_hist_done = true;
0654       if (!muon_sel[5] || flags_passed == NFLAGS)
0655         mt_after_->Fill(massT);
0656       if (!muon_sel[6] || flags_passed == NFLAGS)
0657         if (!met_hist_done)
0658           met_after_->Fill(met_et);
0659       met_hist_done = true;
0660       if (!muon_sel[7] || flags_passed == NFLAGS)
0661         acop_after_->Fill(acop);
0662       // no action here for muon_sel[8]
0663       if (!muon_sel[9] || flags_passed == NFLAGS) {
0664         if (!njets_hist_done) {
0665           njets_after_->Fill(njets);
0666           leadingjet_pt_after_->Fill(lead_jet_pt);
0667           leadingjet_eta_after_->Fill(lead_jet_eta);
0668         }
0669         njets_hist_done = true;
0670       }
0671       if (flags_passed == NFLAGS) {
0672         if (!wfullsel_hist_done) {
0673           npvs_after_->Fill(nvvertex);
0674           muoncharge_after_->Fill(charge);
0675           // if (charge>0) ptPlus_afterW_->Fill(pt);
0676           // if (charge<0) ptMinus_afterW_->Fill(pt);
0677         }
0678         wfullsel_hist_done = true;
0679       }
0680     }
0681 
0682     // The cases in which the event is rejected as a Z are considered
0683     // independently:
0684     if (muon4Z && !muon_sel[8]) {
0685       // Plots for 2 muons
0686       for (unsigned int j = i + 1; j < muonCollectionSize; j++) {
0687         for (int ij = 0; ij < NFLAGSZ; ++ij) {
0688           zmuon_sel[ij] = false;
0689         }
0690 
0691         for (int ji = 0; ji < 5; ++ji) {
0692           zmuon_sel[ji] = muon_sel[ji];
0693         }
0694 
0695         const Muon& mu2 = muonCollection->at(j);
0696         if (!mu2.isGlobalMuon())
0697           continue;
0698         if (mu2.charge() * charge != -1)
0699           continue;
0700         reco::TrackRef gm2 = mu2.globalTrack();
0701         reco::TrackRef tk2 = mu2.innerTrack();
0702         double pt2 = mu2.pt();
0703         if (pt2 > ptThrForZ2_)
0704           zmuon_sel[5] = true;
0705         double eta2 = mu2.eta();
0706         if (fabs(eta2) < etaCut_)
0707           zmuon_sel[6] = true;
0708         double dxy2 = gm2->dxy(beamSpotHandle->position());
0709         if (fabs(dxy2) < dxyCut_)
0710           zmuon_sel[7] = true;
0711         double normalizedChi22 = gm2->normalizedChi2();
0712         double trackerHits2 = tk2->hitPattern().numberOfValidTrackerHits();
0713         int pixelHits2 = tk2->hitPattern().numberOfValidPixelHits();
0714         int muonHits2 = gm2->hitPattern().numberOfValidMuonHits();
0715         int nMatches2 = mu2.numberOfMatches();
0716         bool quality2 = true;
0717         if (normalizedChi22 > normalizedChi2Cut_)
0718           quality2 = false;
0719         if (trackerHits2 < trackerHitsCut_)
0720           quality2 = false;
0721         if (pixelHits2 < pixelHitsCut_)
0722           quality2 = false;
0723         if (muonHits2 < muonHitsCut_)
0724           quality2 = false;
0725         if (!mu2.isTrackerMuon())
0726           quality2 = false;
0727         if (nMatches2 < nMatchesCut_)
0728           quality2 = false;
0729         zmuon_sel[8] = quality2;
0730         double isovar2 = mu2.isolationR03().sumPt;
0731         if (isCombinedIso_) {
0732           isovar2 += mu2.isolationR03().emEt;
0733           isovar2 += mu2.isolationR03().hadEt;
0734         }
0735         if (isRelativeIso_)
0736           isovar2 /= pt2;
0737         if (isovar2 < isoCut03_)
0738           zmuon_sel[9] = true;
0739         //  if (trigger_fired) zmuon_sel[10] = true;
0740         const math::XYZTLorentzVector ZRecoGlb(
0741             mu.px() + mu2.px(), mu.py() + mu2.py(), mu.pz() + mu2.pz(), mu.p() + mu2.p());
0742         if (ZRecoGlb.mass() > dimuonMassMin_ && ZRecoGlb.mass() < dimuonMassMax_)
0743           zmuon_sel[10] = true;
0744 
0745         // jet flag
0746         if (njets <= nJetMax_)
0747           zmuon_sel[11] = true;
0748 
0749         // start filling histos: N-1 plots
0750         int flags_passed_z = 0;
0751 
0752         for (int jj = 0; jj < NFLAGSZ; ++jj) {
0753           if (zmuon_sel[jj])
0754             ++flags_passed_z;
0755         }
0756 
0757         if (flags_passed_z >= (NFLAGSZ - 1)) {
0758           if (!zmuon_sel[0] || flags_passed_z == NFLAGSZ) {
0759             pt1_afterZ_->Fill(pt);
0760           }
0761           if (!zmuon_sel[1] || flags_passed_z == NFLAGSZ) {
0762             eta1_afterZ_->Fill(eta);
0763           }
0764           if (!zmuon_sel[2] || flags_passed_z == NFLAGSZ) {
0765             dxy1_afterZ_->Fill(dxy);
0766           }
0767           if (!zmuon_sel[3] || flags_passed_z == NFLAGSZ) {
0768             goodewkmuon1_afterZ_->Fill(quality);
0769           }
0770           if (!zmuon_sel[4] || flags_passed_z == NFLAGSZ) {
0771             iso1_afterZ_->Fill(isovar);
0772           }
0773           if (!zmuon_sel[5] || flags_passed_z == NFLAGSZ) {
0774             pt2_afterZ_->Fill(pt2);
0775           }
0776           if (!zmuon_sel[6] || flags_passed_z == NFLAGSZ) {
0777             eta2_afterZ_->Fill(eta2);
0778           }
0779           if (!zmuon_sel[7] || flags_passed_z == NFLAGSZ) {
0780             dxy2_afterZ_->Fill(dxy2);
0781           }
0782           if (!zmuon_sel[8] || flags_passed_z == NFLAGSZ) {
0783             goodewkmuon2_afterZ_->Fill(quality2);
0784           }
0785           if (!zmuon_sel[9] || flags_passed_z == NFLAGSZ) {
0786             iso2_afterZ_->Fill(isovar2);
0787           }
0788           //     if (!zmuon_sel[10] || flags_passed_z == NFLAGSZ) {
0789           //     if (!zhlt_hist_done) ztrig_afterZ_->Fill(trigger_fired);
0790           //   zhlt_hist_done = true;
0791           //  }
0792           if (!zmuon_sel[10] || flags_passed_z == NFLAGSZ) {
0793             dimuonmass_afterZ_->Fill(ZRecoGlb.mass());
0794           }
0795           if (!zmuon_sel[11] || flags_passed_z == NFLAGSZ) {
0796             if (!zjets_hist_done) {
0797               njets_afterZ_->Fill(njets);
0798               leadingjet_pt_afterZ_->Fill(lead_jet_pt);
0799               leadingjet_eta_afterZ_->Fill(lead_jet_eta);
0800             }
0801             zjets_hist_done = true;
0802           }
0803           if (flags_passed_z == NFLAGSZ) {
0804             met_afterZ_->Fill(met_et);
0805             if (!zfullsel_hist_done) {
0806               npvs_afterZ_->Fill(nvvertex);
0807               muoncharge_afterZ_->Fill(charge);
0808               if (charge > 0) {
0809                 // ptPlus_afterZ_->Fill(mu.pt());
0810                 // ptMinus_afterZ_->Fill(mu2.pt());
0811                 ptDiffPM_afterZ_->Fill(mu.pt() - mu2.pt());
0812               } else {
0813                 // ptPlus_afterZ_->Fill(mu2.pt());
0814                 // ptMinus_afterZ_->Fill(mu.pt());
0815                 ptDiffPM_afterZ_->Fill(mu2.pt() - mu.pt());
0816               }
0817             }
0818             zfullsel_hist_done = true;
0819           }
0820         }
0821       }
0822     }
0823   }
0824 
0825   if (zfullsel_hist_done) {
0826     // here was a Z candidate
0827     n_zselPt1thr_->Fill(nmuonsForZ1);
0828     n_zselPt2thr_->Fill(nmuonsForZ2);
0829   }
0830 
0831   // nmuons_->Fill(muonCollectionSize);
0832   ngoodmuons_->Fill(number_of_goodMuons);
0833 
0834   return;
0835 }
0836 
0837 // Local Variables:
0838 // show-trailing-whitespace: t
0839 // truncate-lines: t
0840 // End:
[ 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!