Physics/src/EwkMuLumiMonitorDQM.cc

0001 #include "DQM/Physics/src/EwkMuLumiMonitorDQM.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/Common/interface/View.h"
0013 #include "DataFormats/Common/interface/Handle.h"
0014
0015 #include "DataFormats/RecoCandidate/interface/IsoDepositVetos.h"
0016 #include "DataFormats/RecoCandidate/interface/IsoDepositDirection.h"
0017
0018 #include "FWCore/Common/interface/TriggerNames.h"
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "DataFormats/Common/interface/TriggerResults.h"
0021 #include "DataFormats/HLTReco/interface/TriggerObject.h"
0022 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0023 #include "DataFormats/Math/interface/deltaR.h"
0024
0025 #include "DataFormats/METReco/interface/MET.h"
0026 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0027
0028 #include "DataFormats/Math/interface/LorentzVector.h"
0029
0030 using namespace edm;
0031 using namespace std;
0032 using namespace reco;
0033 using namespace isodeposit;
0034
0035 EwkMuLumiMonitorDQM::EwkMuLumiMonitorDQM(const ParameterSet& cfg)
0036     :  // Input collections
0037       trigTag_(cfg.getUntrackedParameter<edm::InputTag>("TrigTag", edm::InputTag("TriggerResults::HLT"))),
0038       trigToken_(consumes<edm::TriggerResults>(trigTag_)),
0039       trigEvToken_(consumes<trigger::TriggerEvent>(cfg.getUntrackedParameter<edm::InputTag>("triggerEvent"))),
0040       beamSpotToken_(consumes<reco::BeamSpot>(
0041           cfg.getUntrackedParameter<edm::InputTag>("offlineBeamSpot", edm::InputTag("offlineBeamSpot")))),
0042       muonToken_(consumes<edm::View<reco::Muon> >(cfg.getUntrackedParameter<edm::InputTag>("muons"))),
0043       trackToken_(consumes<reco::TrackCollection>(cfg.getUntrackedParameter<edm::InputTag>("tracks"))),
0044       caloTowerToken_(consumes<CaloTowerCollection>(cfg.getUntrackedParameter<edm::InputTag>("calotower"))),
0045       metToken_(consumes<edm::View<reco::MET> >(cfg.getUntrackedParameter<edm::InputTag>("metTag"))),
0046       metIncludesMuons_(cfg.getUntrackedParameter<bool>("METIncludesMuons")),
0047       // Main cuts
0048       // massMin_(cfg.getUntrackedParameter<double>("MtMin", 20.)),
0049       //   massMax_(cfg.getUntrackedParameter<double>("MtMax", 2000.))
0050       //  hltPath_(cfg.getUntrackedParameter<std::string> ("hltPath")) ,
0051       //  L3FilterName_(cfg.getUntrackedParameter<std::string>
0052       // ("L3FilterName")),
0053       ptMuCut_(cfg.getUntrackedParameter<double>("ptMuCut")),
0054       etaMuCut_(cfg.getUntrackedParameter<double>("etaMuCut")),
0055       isRelativeIso_(cfg.getUntrackedParameter<bool>("IsRelativeIso")),
0056       isCombinedIso_(cfg.getUntrackedParameter<bool>("IsCombinedIso")),
0057       isoCut03_(cfg.getUntrackedParameter<double>("IsoCut03")),
0058       //  deltaRTrk_(cfg.getUntrackedParameter<double>("deltaRTrk")),
0059       ptThreshold_(cfg.getUntrackedParameter<double>("ptThreshold")),
0060       // deltaRVetoTrk_(cfg.getUntrackedParameter<double>("deltaRVetoTrk")),
0061       maxDPtRel_(cfg.getUntrackedParameter<double>("maxDPtRel")),
0062       maxDeltaR_(cfg.getUntrackedParameter<double>("maxDeltaR")),
0063       mtMin_(cfg.getUntrackedParameter<double>("mtMin")),
0064       mtMax_(cfg.getUntrackedParameter<double>("mtMax")),
0065       acopCut_(cfg.getUntrackedParameter<double>("acopCut")),
0066       dxyCut_(cfg.getUntrackedParameter<double>("DxyCut")) {
0067   // just to initialize
0068   isValidHltConfig_ = false;
0069 }
0070
0071 void EwkMuLumiMonitorDQM::dqmBeginRun(const Run& r, const EventSetup& iSetup) {
0072   nall = 0;
0073   nEvWithHighPtMu = 0;
0074   nInKinRange = 0;
0075   nsel = 0;
0076   niso = 0;
0077   nhlt = 0;
0078   n1hlt = 0;
0079   n2hlt = 0;
0080   nNotIso = 0;
0081   nGlbSta = 0;
0082   nGlbTrk = 0;
0083   nTMass = 0;
0084   nW = 0;
0085
0086   // passed as parameter to HLTConfigProvider::init(), not yet used
0087   bool isConfigChanged = false;
0088
0089   // isValidHltConfig_ used to short-circuit analyze() in case of problems
0090   isValidHltConfig_ = hltConfigProvider_.init(r, iSetup, trigTag_.process(), isConfigChanged);
0091   // std::cout << "hlt config trigger is valid??" << isValidHltConfig_ <<
0092   // std::endl;
0093 }
0094
0095 void EwkMuLumiMonitorDQM::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const&, edm::EventSetup const&) {
0096   ibooker.setCurrentFolder("Physics/EwkMuLumiMonitorDQM");
0097
0098   mass2HLT_ = ibooker.book1D("Z_2HLT_MASS", "Z mass [GeV/c^{2}]", 200, 0., 200.);
0099   mass1HLT_ = ibooker.book1D("Z_1HLT_MASS", "Z mass [GeV/c^{2}]", 200, 0., 200.);
0100   massNotIso_ = ibooker.book1D("Z_NOTISO_MASS", "Z mass [GeV/c^{2}]", 200, 0., 200.);
0101   massGlbSta_ = ibooker.book1D("Z_GLBSTA_MASS", "Z mass [GeV/c^{2}]", 200, 0., 200.);
0102   massGlbTrk_ = ibooker.book1D("Z_GLBTRK_MASS", "Z mass [GeV/c^{2}]", 200, 0., 200.);
0103   massIsBothGlbTrkThanW_ = ibooker.book1D("Z_ISBOTHGLBTRKTHANW_MASS", "Z mass [GeV/c^{2}]", 200, 0., 200.);
0104
0105   highMass2HLT_ = ibooker.book1D("Z_2HLT_HIGHMASS", "Z high mass [GeV/c^{2}]", 2000, 0., 2000.);
0106   highMass1HLT_ = ibooker.book1D("Z_1HLT_HIGHMASS", "Z high mass [GeV/c^{2}]", 2000, 0., 2000.);
0107   highMassNotIso_ = ibooker.book1D("Z_NOTISO_HIGHMASS", "Z high mass [GeV/c^{2}]", 2000, 0., 2000.);
0108   highMassGlbSta_ = ibooker.book1D("Z_GLBSTA_HIGHMASS", "Z high mass [GeV/c^{2}]", 2000, 0., 2000.);
0109   highMassGlbTrk_ = ibooker.book1D("Z_GLBTRK_HIGHMASS", "Z high mass [GeV/c^{2}]", 2000, 0., 2000.);
0110   highMassIsBothGlbTrkThanW_ =
0111       ibooker.book1D("Z_ISBOTHGLBTRKTHANW_HIGHMASS", "Z high mass [GeV/c^{2}]", 2000, 0., 2000.);
0112
0113   TMass_ = ibooker.book1D("TMASS", "Transverse mass [GeV]", 300, 0., 300.);
0114 }
0115
0116 double EwkMuLumiMonitorDQM::muIso(const reco::Muon& mu) {
0117   double isovar = mu.isolationR03().sumPt;
0118   if (isCombinedIso_) {
0119     isovar += mu.isolationR03().emEt;
0120     isovar += mu.isolationR03().hadEt;
0121   }
0122   if (isRelativeIso_)
0123     isovar /= mu.pt();
0124   return isovar;
0125 }
0126
0127 double EwkMuLumiMonitorDQM::tkIso(const reco::Track& tk,
0128                                   Handle<TrackCollection> tracks,
0129                                   Handle<CaloTowerCollection> calotower) {
0130   double ptSum = 0;
0131   for (size_t i = 0; i < tracks->size(); ++i) {
0132     const reco::Track& elem = tracks->at(i);
0133     double elemPt = elem.pt();
0134     // same parameter used for muIsolation: dR [0.01, IsoCut03_], |dZ|<0.2,
0135     // |d_r(xy)|<0.1
0136     double elemVx = elem.vx();
0137     double elemVy = elem.vy();
0138     double elemD0 = sqrt(elemVx * elemVx + elemVy * elemVy);
0139     if (elemD0 > 0.2)
0140       continue;
0141     double dz = fabs(elem.vz() - tk.vz());
0142     if (dz > 0.1)
0143       continue;
0144     // evaluate only for tracks with pt>ptTreshold
0145     if (elemPt < ptThreshold_)
0146       continue;
0147     double dR = deltaR(elem.eta(), elem.phi(), tk.eta(), tk.phi());
0148     // isolation in a cone with dR=0.3, and vetoing the track itself
0149     if ((dR < 0.01) || (dR > 0.3))
0150       continue;
0151     ptSum += elemPt;
0152   }
0153   if (isCombinedIso_) {
0154     // loop on clusters....
0155     for (CaloTowerCollection::const_iterator it = calotower->begin(); it != calotower->end(); ++it) {
0156       double dR = deltaR(it->eta(), it->phi(), tk.outerEta(), tk.outerPhi());
0157       // veto value is 0.1 for towers....
0158       if ((dR < 0.1) || (dR > 0.3))
0159         continue;
0160       ptSum += it->emEnergy();
0161       ptSum += it->hadEnergy();
0162     }
0163   }
0164   if (isRelativeIso_)
0165     ptSum /= tk.pt();
0166   return ptSum;
0167 }
0168
0169 bool EwkMuLumiMonitorDQM::IsMuMatchedToHLTMu(const reco::Muon& mu,
0170                                              const std::vector<reco::Particle>& HLTMu,
0171                                              double DR,
0172                                              double DPtRel) {
0173   size_t dim = HLTMu.size();
0174   size_t nPass = 0;
0175   if (dim == 0)
0176     return false;
0177   for (size_t k = 0; k < dim; k++) {
0178     if ((deltaR(HLTMu[k], mu) < DR) && (fabs(HLTMu[k].pt() - mu.pt()) / HLTMu[k].pt() < DPtRel)) {
0179       nPass++;
0180     }
0181   }
0182   return (nPass > 0);
0183 }
0184
0185 void EwkMuLumiMonitorDQM::analyze(const Event& ev, const EventSetup&) {
0186   nall++;
0187   bool hlt_sel = false;
0188   double iso1 = -1;
0189   double iso2 = -1;
0190   bool isMu1Iso = false;
0191   bool isMu2Iso = false;
0192   bool singleTrigFlag1 = false;
0193   bool singleTrigFlag2 = false;
0194   isZGolden1HLT_ = false;
0195   isZGolden2HLT_ = false;
0196   isZGoldenNoIso_ = false;
0197   isZGlbSta_ = false;
0198   isZGlbTrk_ = false;
0199   isW_ = false;
0200   // Trigger
0201   bool trigger_fired = false;
0202
0203   Handle<TriggerResults> triggerResults;
0204   if (!ev.getByToken(trigToken_, triggerResults)) {
0205     // LogWarning("") << ">>> TRIGGER collection does not exist !!!";
0206     return;
0207   }
0208
0209   ev.getByToken(trigToken_, triggerResults);
0210   /*
0211     const edm::TriggerNames & trigNames = ev.triggerNames(*triggerResults);
0212
0213
0214   for (size_t i=0; i<triggerResults->size(); i++) {
0215   std::string trigName = trigNames.triggerName(i);
0216   //std::cout << " trigName == " << trigName << std::endl;
0217     if ( trigName == hltPath_ && triggerResults->accept(i)) {
0218     trigger_fired = true;
0219     hlt_sel=true;
0220     nhlt++;
0221     }
0222     }
0223   */
0224   // see the trigger single muon which are present
0225   string lowestMuonUnprescaledTrig = "";
0226   bool lowestMuonUnprescaledTrigFound = false;
0227   const std::vector<std::string>& triggerNames = hltConfigProvider_.triggerNames();
0228   for (size_t ts = 0; ts < triggerNames.size(); ts++) {
0229     string trig = triggerNames[ts];
0230     size_t f = trig.find("HLT_Mu");
0231     if ((f != std::string::npos)) {
0232       // std::cout << "single muon trigger present: " << trig << std::endl;
0233       // See if the trigger is prescaled;
0234       /// number of prescale sets available
0235       bool prescaled = false;
0236       const unsigned int prescaleSize = hltConfigProvider_.prescaleSize();
0237       for (unsigned int ps = 0; ps < prescaleSize; ps++) {
0238         const unsigned int prescaleValue = hltConfigProvider_.prescaleValue(ps, trig);
0239         if (prescaleValue != 1)
0240           prescaled = true;
0241         //  std::cout<< " prescaleValue[" << ps << "] =" << prescaleValue
0242         //<<std::endl;
0243       }
0244       if (!prescaled) {
0245         // looking now for the lowest hlt path not prescaled, with name of the
0246         // form HLT_MuX or HLTMuX_vY
0247         for (unsigned int n = 9; n < 100; n++) {
0248           string lowestTrig = "HLT_Mu";
0249           string lowestTrigv0 = "copy";
0250           std::stringstream out;
0251           out << n;
0252           std::string s = out.str();
0253           lowestTrig.append(s);
0254           lowestTrigv0 = lowestTrig;
0255           for (unsigned int v = 1; v < 10; v++) {
0256             lowestTrig.append("_v");
0257             std::stringstream oout;
0258             oout << v;
0259             std::string ss = oout.str();
0260             lowestTrig.append(ss);
0261             if (trig == lowestTrig)
0262               lowestMuonUnprescaledTrig = trig;
0263             if (trig == lowestTrig)
0264               lowestMuonUnprescaledTrigFound = true;
0265             if (trig == lowestTrig)
0266               break;
0267           }
0268           if (lowestMuonUnprescaledTrigFound)
0269             break;
0270
0271           lowestTrig = lowestTrigv0;
0272           if (trig == lowestTrig)
0273             lowestMuonUnprescaledTrig = trig;
0274           //      if (trig==lowestTrig) {std::cout << " before break, lowestTrig
0275           // lowest single muon trigger present unprescaled: " << lowestTrig <<
0276           // std::endl; }
0277           if (trig == lowestTrig)
0278             lowestMuonUnprescaledTrigFound = true;
0279           if (trig == lowestTrig)
0280             break;
0281         }
0282         if (lowestMuonUnprescaledTrigFound)
0283           break;
0284       }
0285     }
0286   }
0287   //  std::cout << "after break, lowest single muon trigger present unprescaled:
0288   // " << lowestMuonUnprescaledTrig << std::endl;
0289   unsigned int triggerIndex;  // index of trigger path
0290
0291   // See if event passed trigger paths
0292   std::string hltPath_ = lowestMuonUnprescaledTrig;
0293
0294   triggerIndex = hltConfigProvider_.triggerIndex(hltPath_);
0295   if (triggerIndex < triggerResults->size())
0296     trigger_fired = triggerResults->accept(triggerIndex);
0297   std::string L3FilterName_ = "";
0298   if (trigger_fired) {
0299     const std::vector<std::string>& moduleLabs = hltConfigProvider_.moduleLabels(hltPath_);
0300     /*for (size_t k =0; k < moduleLabs.size()-1 ; k++){
0301       std::cout << "moduleLabs[" << k << "] == " << moduleLabs[k] << std::endl;
0302     }
0303     */
0304     // the l3 filter name is just the last module....
0305     size_t moduleLabsSizeMinus2 = moduleLabs.size() - 2;
0306     //  std::cout<<"moduleLabs[" << moduleLabsSizeMinus2 << "]== "<<
0307     // moduleLabs[moduleLabsSizeMinus2] << std::endl;
0308
0309     L3FilterName_ = moduleLabs[moduleLabsSizeMinus2];
0310   }
0311
0312   edm::Handle<trigger::TriggerEvent> handleTriggerEvent;
0313   LogTrace("") << ">>> Trigger bit: " << trigger_fired << " (" << hltPath_ << ")";
0314   if (!ev.getByToken(trigEvToken_, handleTriggerEvent)) {
0315     // LogWarning( "errorTriggerEventValid" ) << "trigger::TriggerEvent product
0316     // with InputTag " << trigEv_.encode() << " not in event";
0317     return;
0318   }
0319   ev.getByToken(trigEvToken_, handleTriggerEvent);
0320   const trigger::TriggerObjectCollection& toc(handleTriggerEvent->getObjects());
0321   size_t nMuHLT = 0;
0322   std::vector<reco::Particle> HLTMuMatched;
0323   for (size_t ia = 0; ia < handleTriggerEvent->sizeFilters(); ++ia) {
0324     std::string fullname = handleTriggerEvent->filterTag(ia).encode();
0325     // std::cout<< "fullname::== " << fullname<< std::endl;
0326     std::string name;
0327     size_t p = fullname.find_first_of(':');
0328     if (p != std::string::npos) {
0329       name = fullname.substr(0, p);
0330     } else {
0331       name = fullname;
0332     }
0333     if (!toc.empty()) {
0334       const trigger::Keys& k = handleTriggerEvent->filterKeys(ia);
0335       for (trigger::Keys::const_iterator ki = k.begin(); ki != k.end(); ++ki) {
0336         // looking at all the single muon l3 trigger present, for example
0337         // hltSingleMu15L3Filtered15.....
0338         if (name == L3FilterName_) {
0339           // trigger_fired = true;
0340           hlt_sel = true;
0341           nhlt++;
0342           HLTMuMatched.push_back(toc[*ki].particle());
0343           nMuHLT++;
0344         }
0345       }
0346     }
0347   }
0348
0349   // Beam spot
0350   Handle<reco::BeamSpot> beamSpotHandle;
0351   if (!ev.getByToken(beamSpotToken_, beamSpotHandle)) {
0352     // LogWarning("") << ">>> No beam spot found !!!";
0353     return;
0354   }
0355
0356   //  looping on muon....
0357   Handle<View<Muon> > muons;
0358   if (!ev.getByToken(muonToken_, muons)) {
0359     // LogError("") << ">>> muon collection does not exist !!!";
0360     return;
0361   }
0362
0363   ev.getByToken(muonToken_, muons);
0364   // saving only muons with pt> ptMuCut and eta<etaMuCut, and dxy<dxyCut
0365   std::vector<reco::Muon> highPtGlbMuons;
0366   std::vector<reco::Muon> highPtStaMuons;
0367
0368   for (size_t i = 0; i < muons->size(); i++) {
0369     const reco::Muon& mu = muons->at(i);
0370     double pt = mu.pt();
0371     double eta = mu.eta();
0372     if (pt > ptMuCut_ && fabs(eta) < etaMuCut_) {
0373       if (mu.isGlobalMuon()) {
0374         // check the dxy....
0375         double dxy = mu.innerTrack()->dxy(beamSpotHandle->position());
0376         if (fabs(dxy) > dxyCut_)
0377           continue;
0378         highPtGlbMuons.push_back(mu);
0379       }
0380       if (mu.isGlobalMuon())
0381         continue;
0382       // if is not, look if it is a standalone....
0383       if (mu.isStandAloneMuon())
0384         highPtStaMuons.push_back(mu);
0385       nEvWithHighPtMu++;
0386     }
0387   }
0388   size_t nHighPtGlbMu = highPtGlbMuons.size();
0389   size_t nHighPtStaMu = highPtStaMuons.size();
0390   if (hlt_sel && (nHighPtGlbMu > 0)) {
0391     // loop on high pt muons if there's at least two with opposite charge build
0392     // a Z, more then one z candidate is foreseen.........
0393     // stop the loop after 10 cicles....
0394     (nHighPtGlbMu > 10) ? nHighPtGlbMu = 10 : 1;
0395     // Z selection
0396     if (nHighPtGlbMu > 1) {
0397       for (unsigned int i = 0; i < nHighPtGlbMu; i++) {
0398         reco::Muon muon1 = highPtGlbMuons[i];
0399         const math::XYZTLorentzVector& mu1(muon1.p4());
0400         // double pt1= muon1.pt();
0401         for (unsigned int j = i + 1; j < nHighPtGlbMu; ++j) {
0402           reco::Muon muon2 = highPtGlbMuons[j];
0403           const math::XYZTLorentzVector& mu2(muon2.p4());
0404           // double pt2= muon2.pt();
0405           if (muon1.charge() == muon2.charge())
0406             continue;
0407           math::XYZTLorentzVector pair = mu1 + mu2;
0408           double mass = pair.M();
0409           // checking isolation and hlt maching
0410           iso1 = muIso(muon1);
0411           iso2 = muIso(muon2);
0412           isMu1Iso = (iso1 < isoCut03_);
0413           isMu2Iso = (iso2 < isoCut03_);
0414           singleTrigFlag1 = IsMuMatchedToHLTMu(muon1, HLTMuMatched, maxDeltaR_, maxDPtRel_);
0415           singleTrigFlag2 = IsMuMatchedToHLTMu(muon2, HLTMuMatched, maxDeltaR_, maxDPtRel_);
0416           if (singleTrigFlag1 && singleTrigFlag2)
0417             isZGolden2HLT_ = true;
0418           if ((singleTrigFlag1 && !singleTrigFlag2) || (!singleTrigFlag1 && singleTrigFlag2))
0419             isZGolden1HLT_ = true;
0420           // Z Golden passing all criteria, with 2 or 1 muon matched to an HLT
0421           // muon. Using the two cathegories as a control sample for the HLT
0422           // matching efficiency
0423           if (isZGolden2HLT_ || isZGolden1HLT_) {
0424             // a Z golden has been found, let's remove the two muons from the
0425             // list, dome for avoiding resolution effect enter muons in the
0426             // standalone and tracker collections.........
0427             highPtGlbMuons.erase(highPtGlbMuons.begin() + i);
0428             highPtGlbMuons.erase(highPtGlbMuons.begin() + j);
0429             // and updating the number of high pt muons....
0430             nHighPtGlbMu = highPtGlbMuons.size();
0431
0432             if (isMu1Iso && isMu2Iso) {
0433               niso++;
0434               if (isZGolden2HLT_) {
0435                 n2hlt++;
0436                 mass2HLT_->Fill(mass);
0437                 highMass2HLT_->Fill(mass);
0438                 /*
0439                   if (pt1 > pt2) {
0440                   highest_mupt2HLT_ -> Fill (pt1);
0441                   lowest_mupt2HLT_ -> Fill (pt2);
0442                   } else {
0443                   highest_mupt2HLT_ -> Fill (pt2);
0444                   lowest_mupt2HLT_ -> Fill (pt1);
0445                   }
0446                 */
0447               }
0448               if (isZGolden1HLT_) {
0449                 n1hlt++;
0450                 mass1HLT_->Fill(mass);
0451                 highMass1HLT_->Fill(mass);
0452                 /*
0453                   if (pt1 >pt2) {
0454                   highest_mupt1HLT_ -> Fill (pt1);
0455                   lowest_mupt1HLT_ -> Fill (pt2);
0456                 } else {
0457                   highest_mupt1HLT_ -> Fill (pt2);
0458                   lowest_mupt1HLT_ -> Fill (pt1);
0459                 }
0460                 */
0461               }
0462             } else {
0463               // ZGlbGlb when at least one of the two muon is not isolated and
0464               // at least one HLT matched, used as control sample for the
0465               // isolation efficiency
0466               // filling events with one muon not isolated and both hlt mathced
0467               if (((isMu1Iso && !isMu2Iso) || (!isMu1Iso && isMu2Iso)) && (isZGolden2HLT_ && isZGolden1HLT_)) {
0468                 isZGoldenNoIso_ = true;
0469                 nNotIso++;
0470                 massNotIso_->Fill(mass);
0471                 highMassNotIso_->Fill(mass);
0472               }
0473               /*
0474                 if (pt1 > pt2) {
0475                 highest_muptNotIso_ -> Fill (pt1);
0476                 lowest_muptNotIso_ -> Fill (pt2);
0477                 } else {
0478                 highest_muptNotIso_ -> Fill (pt2);
0479                 lowest_muptNotIso_ -> Fill (pt1);
0480                 }
0481               */
0482             }
0483           }
0484         }
0485       }
0486     }
0487     // W->MuNu selection (if at least one muon has been selected)
0488     // looking for a W if a Z is not found.... let's think if we prefer to
0489     // exclude zMuMuNotIso or zMuSta....
0490     if (!(isZGolden2HLT_ || isZGolden1HLT_)) {
0491       Handle<View<MET> > metCollection;
0492       if (!ev.getByToken(metToken_, metCollection)) {
0493         // LogError("") << ">>> MET collection does not exist !!!";
0494         return;
0495       }
0496       const MET& met = metCollection->at(0);
0497       nW = 0;
0498       for (unsigned int i = 0; i < nHighPtGlbMu; i++) {
0499         reco::Muon muon1 = highPtGlbMuons[i];
0500         /*
0501                quality cut not implemented
0502             Quality Cuts           double dxy =
0503 gm->dxy(beamSpotHandle->position());
0504             Cut in 0.2           double trackerHits =
0505 gm->hitPattern().numberOfValidTrackerHits();
0506             Cut in 11           bool quality = fabs(dxy)<dxyCut_  &&
0507 muon::isGoodMuon(mu,muon::GlobalMuonPromptTight) && trackerHits>=trackerHitsCut_
0508 && mu.isTrackerMuon() ;
0509 if (!quality) continue;
0510         */
0511         // isolation cut and hlt maching
0512         iso1 = muIso(muon1);
0513         isMu1Iso = (iso1 < isoCut03_);
0514         if (!isMu1Iso)
0515           continue;
0516         // checking if muon is matched to any HLT muon....
0517         singleTrigFlag1 = IsMuMatchedToHLTMu(muon1, HLTMuMatched, maxDeltaR_, maxDPtRel_);
0518         if (!singleTrigFlag1)
0519           continue;
0520         // std::cout << " is GlobMu macthecd to HLT: " << singleTrigFlag1 <<
0521         // std::endl;
0522         // MT cuts
0523         double met_px = met.px();
0524         double met_py = met.py();
0525
0526         if (!metIncludesMuons_) {
0527           for (unsigned int i = 0; i < nHighPtGlbMu; i++) {
0528             reco::Muon muon1 = highPtGlbMuons[i];
0529             met_px -= muon1.px();
0530             met_py -= muon1.py();
0531           }
0532         }
0533         double met_et = met.pt();  // sqrt(met_px*met_px+met_py*met_py);
0534         LogTrace("") << ">>> MET, MET_px, MET_py: " << met_et << ", " << met_px << ", " << met_py << " [GeV]";
0535         double w_et = met_et + muon1.pt();
0536         double w_px = met_px + muon1.px();
0537         double w_py = met_py + muon1.py();
0538         double massT = w_et * w_et - w_px * w_px - w_py * w_py;
0539         massT = (massT > 0) ? sqrt(massT) : 0;
0540         // Acoplanarity cuts
0541         Geom::Phi<double> deltaphi(muon1.phi() - atan2(met_py, met_px));
0542         double acop = M_PI - fabs(deltaphi.value());
0543         // std::cout << " is acp of W candidate less then cut: " << (acop<
0544         // acopCut_) << std::endl;
0545         if (acop > acopCut_)
0546           continue;  // Cut in 2.0
0547         TMass_->Fill(massT);
0548         if (massT < mtMin_ || massT > mtMax_)
0549           continue;  // Cut in (50,200) GeV
0550         // we give the number of W only in the tMass selected but we have a look
0551         // at mass tails to check the QCD background
0552         isW_ = true;
0553         nW++;
0554         nTMass++;
0555       }
0556     }
0557     // if a zGlobal is not selected, look at the dimuonGlobalOneStandAlone and
0558     // dimuonGlobalOneTrack...., used as a control sample for the tracking
0559     // efficiency  and muon system efficiency
0560     if (!(isZGolden2HLT_ || isZGolden1HLT_ || isZGoldenNoIso_)) {
0561       for (unsigned int i = 0; i < nHighPtGlbMu; ++i) {
0562         reco::Muon glbMuon = highPtGlbMuons[i];
0563         const math::XYZTLorentzVector& mu1(glbMuon.p4());
0564         // double pt1= glbMuon.pt();
0565         // checking that the global muon is hlt matched otherwise skip the event
0566         singleTrigFlag1 = IsMuMatchedToHLTMu(glbMuon, HLTMuMatched, maxDeltaR_, maxDPtRel_);
0567         if (!singleTrigFlag1)
0568           continue;
0569         // checking that the global muon is isolated matched otherwise skip the
0570         // event
0571         iso1 = muIso(glbMuon);
0572         isMu1Iso = (iso1 < isoCut03_);
0573         if (!isMu1Iso)
0574           continue;
0575         // look at the standalone muon ....
0576         // stop the loop after 10 cicles....
0577         (nHighPtStaMu > 10) ? nHighPtStaMu = 10 : 1;
0578         for (unsigned int j = 0; j < nHighPtStaMu; ++j) {
0579           reco::Muon staMuon = highPtStaMuons[j];
0580           const math::XYZTLorentzVector& mu2(staMuon.p4());
0581           // double pt2= staMuon.pt();
0582           if (glbMuon.charge() == staMuon.charge())
0583             continue;
0584           math::XYZTLorentzVector pair = mu1 + mu2;
0585           double mass = pair.M();
0586           iso2 = muIso(staMuon);
0587           isMu2Iso = (iso2 < isoCut03_);
0588           LogTrace("") << "\t... isolation value" << iso1 << ", isolated? " << isMu1Iso;
0589           LogTrace("") << "\t... isolation value" << iso2 << ", isolated? " << isMu2Iso;
0590           // requiring theat the global mu is mathed to the HLT  and both are
0591           // isolated
0592           if (isMu2Iso) {
0593             isZGlbSta_ = true;
0594             nGlbSta++;
0595             massGlbSta_->Fill(mass);
0596             highMassGlbSta_->Fill(mass);
0597             /*
0598               if (pt1 > pt2) {
0599               highest_muptGlbSta_ -> Fill (pt1);
0600               lowest_muptGlbSta_ -> Fill (pt2);
0601             } else {
0602               highest_muptGlbSta_ -> Fill (pt2);
0603               lowest_muptGlbSta_ -> Fill (pt1);
0604             }
0605             */
0606           }
0607         }
0608         // look at the tracks....
0609         Handle<TrackCollection> tracks;
0610         if (!ev.getByToken(trackToken_, tracks)) {
0611           // LogError("") << ">>> track collection does not exist !!!";
0612           return;
0613         }
0614         ev.getByToken(trackToken_, tracks);
0615         Handle<CaloTowerCollection> calotower;
0616         if (!ev.getByToken(caloTowerToken_, calotower)) {
0617           // LogError("") << ">>> calotower collection does not exist !!!";
0618           return;
0619         }
0620         ev.getByToken(caloTowerToken_, calotower);
0621         // avoid to loop on more than 5000 trks
0622         size_t nTrk = tracks->size();
0623         (nTrk > 5000) ? nTrk = 5000 : 1;
0624         for (unsigned int j = 0; j < nTrk; j++) {
0625           const reco::Track& tk = tracks->at(j);
0626           if (glbMuon.charge() == tk.charge())
0627             continue;
0628           double pt2 = tk.pt();
0629           double eta = tk.eta();
0630           double dxy = tk.dxy(beamSpotHandle->position());
0631           if (pt2 < ptMuCut_ || fabs(eta) > etaMuCut_ || fabs(dxy) > dxyCut_)
0632             continue;
0633           // assuming that the track is a mu....
0634           math::XYZTLorentzVector mu2(tk.px(), tk.py(), tk.pz(), sqrt((tk.p() * tk.p()) + (0.10566 * 0.10566)));
0635           math::XYZTLorentzVector pair = mu1 + mu2;
0636           double mass = pair.M();
0637           // now requiring that the tracks is isolated.......
0638           iso2 = tkIso(tk, tracks, calotower);
0639           isMu2Iso = (iso2 < isoCut03_);
0640           //          std::cout << "found a track with rel/comb iso: " << iso2
0641           //<< std::endl;
0642           if (isMu2Iso) {
0643             isZGlbTrk_ = true;
0644             nGlbTrk++;
0645             massGlbTrk_->Fill(mass);
0646             highMassGlbTrk_->Fill(mass);
0647             if (!isW_)
0648               continue;
0649             massIsBothGlbTrkThanW_->Fill(mass);
0650             highMassIsBothGlbTrkThanW_->Fill(mass);
0651             /*
0652               if (pt1 > pt2) {
0653               highest_muptGlbTrk_ -> Fill (pt1);
0654               lowest_muptGlbTrk_ -> Fill (pt2);
0655               } else {
0656               highest_muptGlbTrk_ -> Fill (pt2);
0657               lowest_muptGlbTrk_ -> Fill (pt1);
0658             }
0659             */
0660           }
0661         }
0662       }
0663     }
0664   }
0665   if ((hlt_sel || isZGolden2HLT_ || isZGolden1HLT_ || isZGoldenNoIso_ || isZGlbSta_ || isZGlbTrk_ || isW_)) {
0666     nsel++;
0667     LogTrace("") << ">>>> Event ACCEPTED";
0668   } else {
0669     LogTrace("") << ">>>> Event REJECTED";
0670   }
0671   return;
0672 }
0673
0674 // Local Variables:
0675 // show-trailing-whitespace: t
0676 // truncate-lines: t
0677 // End: