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         if (hltConfigProvider_.prescaleValue<double>(ps, trig) != 1)
0239           prescaled = true;
0240       }
0241       if (!prescaled) {
0242         // looking now for the lowest hlt path not prescaled, with name of the
0243         // form HLT_MuX or HLTMuX_vY
0244         for (unsigned int n = 9; n < 100; n++) {
0245           string lowestTrig = "HLT_Mu";
0246           string lowestTrigv0 = "copy";
0247           std::stringstream out;
0248           out << n;
0249           std::string s = out.str();
0250           lowestTrig.append(s);
0251           lowestTrigv0 = lowestTrig;
0252           for (unsigned int v = 1; v < 10; v++) {
0253             lowestTrig.append("_v");
0254             std::stringstream oout;
0255             oout << v;
0256             std::string ss = oout.str();
0257             lowestTrig.append(ss);
0258             if (trig == lowestTrig)
0259               lowestMuonUnprescaledTrig = trig;
0260             if (trig == lowestTrig)
0261               lowestMuonUnprescaledTrigFound = true;
0262             if (trig == lowestTrig)
0263               break;
0264           }
0265           if (lowestMuonUnprescaledTrigFound)
0266             break;
0267
0268           lowestTrig = lowestTrigv0;
0269           if (trig == lowestTrig)
0270             lowestMuonUnprescaledTrig = trig;
0271           //      if (trig==lowestTrig) {std::cout << " before break, lowestTrig
0272           // lowest single muon trigger present unprescaled: " << lowestTrig <<
0273           // std::endl; }
0274           if (trig == lowestTrig)
0275             lowestMuonUnprescaledTrigFound = true;
0276           if (trig == lowestTrig)
0277             break;
0278         }
0279         if (lowestMuonUnprescaledTrigFound)
0280           break;
0281       }
0282     }
0283   }
0284   //  std::cout << "after break, lowest single muon trigger present unprescaled:
0285   // " << lowestMuonUnprescaledTrig << std::endl;
0286   unsigned int triggerIndex;  // index of trigger path
0287
0288   // See if event passed trigger paths
0289   std::string hltPath_ = lowestMuonUnprescaledTrig;
0290
0291   triggerIndex = hltConfigProvider_.triggerIndex(hltPath_);
0292   if (triggerIndex < triggerResults->size())
0293     trigger_fired = triggerResults->accept(triggerIndex);
0294   std::string L3FilterName_ = "";
0295   if (trigger_fired) {
0296     const std::vector<std::string>& moduleLabs = hltConfigProvider_.moduleLabels(hltPath_);
0297     /*for (size_t k =0; k < moduleLabs.size()-1 ; k++){
0298       std::cout << "moduleLabs[" << k << "] == " << moduleLabs[k] << std::endl;
0299     }
0300     */
0301     // the l3 filter name is just the last module....
0302     size_t moduleLabsSizeMinus2 = moduleLabs.size() - 2;
0303     //  std::cout<<"moduleLabs[" << moduleLabsSizeMinus2 << "]== "<<
0304     // moduleLabs[moduleLabsSizeMinus2] << std::endl;
0305
0306     L3FilterName_ = moduleLabs[moduleLabsSizeMinus2];
0307   }
0308
0309   edm::Handle<trigger::TriggerEvent> handleTriggerEvent;
0310   LogTrace("") << ">>> Trigger bit: " << trigger_fired << " (" << hltPath_ << ")";
0311   if (!ev.getByToken(trigEvToken_, handleTriggerEvent)) {
0312     // LogWarning( "errorTriggerEventValid" ) << "trigger::TriggerEvent product
0313     // with InputTag " << trigEv_.encode() << " not in event";
0314     return;
0315   }
0316   ev.getByToken(trigEvToken_, handleTriggerEvent);
0317   const trigger::TriggerObjectCollection& toc(handleTriggerEvent->getObjects());
0318   std::vector<reco::Particle> HLTMuMatched;
0319   for (size_t ia = 0; ia < handleTriggerEvent->sizeFilters(); ++ia) {
0320     std::string fullname = handleTriggerEvent->filterTag(ia).encode();
0321     // std::cout<< "fullname::== " << fullname<< std::endl;
0322     std::string name;
0323     size_t p = fullname.find_first_of(':');
0324     if (p != std::string::npos) {
0325       name = fullname.substr(0, p);
0326     } else {
0327       name = fullname;
0328     }
0329     if (!toc.empty()) {
0330       const trigger::Keys& k = handleTriggerEvent->filterKeys(ia);
0331       for (trigger::Keys::const_iterator ki = k.begin(); ki != k.end(); ++ki) {
0332         // looking at all the single muon l3 trigger present, for example
0333         // hltSingleMu15L3Filtered15.....
0334         if (name == L3FilterName_) {
0335           // trigger_fired = true;
0336           hlt_sel = true;
0337           nhlt++;
0338           HLTMuMatched.push_back(toc[*ki].particle());
0339         }
0340       }
0341     }
0342   }
0343
0344   // Beam spot
0345   Handle<reco::BeamSpot> beamSpotHandle;
0346   if (!ev.getByToken(beamSpotToken_, beamSpotHandle)) {
0347     // LogWarning("") << ">>> No beam spot found !!!";
0348     return;
0349   }
0350
0351   //  looping on muon....
0352   Handle<View<Muon> > muons;
0353   if (!ev.getByToken(muonToken_, muons)) {
0354     // LogError("") << ">>> muon collection does not exist !!!";
0355     return;
0356   }
0357
0358   ev.getByToken(muonToken_, muons);
0359   // saving only muons with pt> ptMuCut and eta<etaMuCut, and dxy<dxyCut
0360   std::vector<reco::Muon> highPtGlbMuons;
0361   std::vector<reco::Muon> highPtStaMuons;
0362
0363   for (size_t i = 0; i < muons->size(); i++) {
0364     const reco::Muon& mu = muons->at(i);
0365     double pt = mu.pt();
0366     double eta = mu.eta();
0367     if (pt > ptMuCut_ && fabs(eta) < etaMuCut_) {
0368       if (mu.isGlobalMuon()) {
0369         // check the dxy....
0370         double dxy = mu.innerTrack()->dxy(beamSpotHandle->position());
0371         if (fabs(dxy) > dxyCut_)
0372           continue;
0373         highPtGlbMuons.push_back(mu);
0374       }
0375       if (mu.isGlobalMuon())
0376         continue;
0377       // if is not, look if it is a standalone....
0378       if (mu.isStandAloneMuon())
0379         highPtStaMuons.push_back(mu);
0380       nEvWithHighPtMu++;
0381     }
0382   }
0383   size_t nHighPtGlbMu = highPtGlbMuons.size();
0384   size_t nHighPtStaMu = highPtStaMuons.size();
0385   if (hlt_sel && (nHighPtGlbMu > 0)) {
0386     // loop on high pt muons if there's at least two with opposite charge build
0387     // a Z, more then one z candidate is foreseen.........
0388     // stop the loop after 10 cicles....
0389     (nHighPtGlbMu > 10) ? nHighPtGlbMu = 10 : 1;
0390     // Z selection
0391     if (nHighPtGlbMu > 1) {
0392       for (unsigned int i = 0; i < nHighPtGlbMu; i++) {
0393         reco::Muon muon1 = highPtGlbMuons[i];
0394         const math::XYZTLorentzVector& mu1(muon1.p4());
0395         // double pt1= muon1.pt();
0396         for (unsigned int j = i + 1; j < nHighPtGlbMu; ++j) {
0397           reco::Muon muon2 = highPtGlbMuons[j];
0398           const math::XYZTLorentzVector& mu2(muon2.p4());
0399           // double pt2= muon2.pt();
0400           if (muon1.charge() == muon2.charge())
0401             continue;
0402           math::XYZTLorentzVector pair = mu1 + mu2;
0403           double mass = pair.M();
0404           // checking isolation and hlt maching
0405           iso1 = muIso(muon1);
0406           iso2 = muIso(muon2);
0407           isMu1Iso = (iso1 < isoCut03_);
0408           isMu2Iso = (iso2 < isoCut03_);
0409           singleTrigFlag1 = IsMuMatchedToHLTMu(muon1, HLTMuMatched, maxDeltaR_, maxDPtRel_);
0410           singleTrigFlag2 = IsMuMatchedToHLTMu(muon2, HLTMuMatched, maxDeltaR_, maxDPtRel_);
0411           if (singleTrigFlag1 && singleTrigFlag2)
0412             isZGolden2HLT_ = true;
0413           if ((singleTrigFlag1 && !singleTrigFlag2) || (!singleTrigFlag1 && singleTrigFlag2))
0414             isZGolden1HLT_ = true;
0415           // Z Golden passing all criteria, with 2 or 1 muon matched to an HLT
0416           // muon. Using the two cathegories as a control sample for the HLT
0417           // matching efficiency
0418           if (isZGolden2HLT_ || isZGolden1HLT_) {
0419             // a Z golden has been found, let's remove the two muons from the
0420             // list, dome for avoiding resolution effect enter muons in the
0421             // standalone and tracker collections.........
0422             highPtGlbMuons.erase(highPtGlbMuons.begin() + i);
0423             highPtGlbMuons.erase(highPtGlbMuons.begin() + j);
0424             // and updating the number of high pt muons....
0425             nHighPtGlbMu = highPtGlbMuons.size();
0426
0427             if (isMu1Iso && isMu2Iso) {
0428               niso++;
0429               if (isZGolden2HLT_) {
0430                 n2hlt++;
0431                 mass2HLT_->Fill(mass);
0432                 highMass2HLT_->Fill(mass);
0433                 /*
0434                   if (pt1 > pt2) {
0435                   highest_mupt2HLT_ -> Fill (pt1);
0436                   lowest_mupt2HLT_ -> Fill (pt2);
0437                   } else {
0438                   highest_mupt2HLT_ -> Fill (pt2);
0439                   lowest_mupt2HLT_ -> Fill (pt1);
0440                   }
0441                 */
0442               }
0443               if (isZGolden1HLT_) {
0444                 n1hlt++;
0445                 mass1HLT_->Fill(mass);
0446                 highMass1HLT_->Fill(mass);
0447                 /*
0448                   if (pt1 >pt2) {
0449                   highest_mupt1HLT_ -> Fill (pt1);
0450                   lowest_mupt1HLT_ -> Fill (pt2);
0451                 } else {
0452                   highest_mupt1HLT_ -> Fill (pt2);
0453                   lowest_mupt1HLT_ -> Fill (pt1);
0454                 }
0455                 */
0456               }
0457             } else {
0458               // ZGlbGlb when at least one of the two muon is not isolated and
0459               // at least one HLT matched, used as control sample for the
0460               // isolation efficiency
0461               // filling events with one muon not isolated and both hlt mathced
0462               if (((isMu1Iso && !isMu2Iso) || (!isMu1Iso && isMu2Iso)) && (isZGolden2HLT_ && isZGolden1HLT_)) {
0463                 isZGoldenNoIso_ = true;
0464                 nNotIso++;
0465                 massNotIso_->Fill(mass);
0466                 highMassNotIso_->Fill(mass);
0467               }
0468               /*
0469                 if (pt1 > pt2) {
0470                 highest_muptNotIso_ -> Fill (pt1);
0471                 lowest_muptNotIso_ -> Fill (pt2);
0472                 } else {
0473                 highest_muptNotIso_ -> Fill (pt2);
0474                 lowest_muptNotIso_ -> Fill (pt1);
0475                 }
0476               */
0477             }
0478           }
0479         }
0480       }
0481     }
0482     // W->MuNu selection (if at least one muon has been selected)
0483     // looking for a W if a Z is not found.... let's think if we prefer to
0484     // exclude zMuMuNotIso or zMuSta....
0485     if (!(isZGolden2HLT_ || isZGolden1HLT_)) {
0486       Handle<View<MET> > metCollection;
0487       if (!ev.getByToken(metToken_, metCollection)) {
0488         // LogError("") << ">>> MET collection does not exist !!!";
0489         return;
0490       }
0491       const MET& met = metCollection->at(0);
0492       nW = 0;
0493       for (unsigned int i = 0; i < nHighPtGlbMu; i++) {
0494         reco::Muon muon1 = highPtGlbMuons[i];
0495         /*
0496                quality cut not implemented
0497             Quality Cuts           double dxy =
0498 gm->dxy(beamSpotHandle->position());
0499             Cut in 0.2           double trackerHits =
0500 gm->hitPattern().numberOfValidTrackerHits();
0501             Cut in 11           bool quality = fabs(dxy)<dxyCut_  &&
0502 muon::isGoodMuon(mu,muon::GlobalMuonPromptTight) && trackerHits>=trackerHitsCut_
0503 && mu.isTrackerMuon() ;
0504 if (!quality) continue;
0505         */
0506         // isolation cut and hlt maching
0507         iso1 = muIso(muon1);
0508         isMu1Iso = (iso1 < isoCut03_);
0509         if (!isMu1Iso)
0510           continue;
0511         // checking if muon is matched to any HLT muon....
0512         singleTrigFlag1 = IsMuMatchedToHLTMu(muon1, HLTMuMatched, maxDeltaR_, maxDPtRel_);
0513         if (!singleTrigFlag1)
0514           continue;
0515         // std::cout << " is GlobMu macthecd to HLT: " << singleTrigFlag1 <<
0516         // std::endl;
0517         // MT cuts
0518         double met_px = met.px();
0519         double met_py = met.py();
0520
0521         if (!metIncludesMuons_) {
0522           for (unsigned int i = 0; i < nHighPtGlbMu; i++) {
0523             reco::Muon muon1 = highPtGlbMuons[i];
0524             met_px -= muon1.px();
0525             met_py -= muon1.py();
0526           }
0527         }
0528         double met_et = met.pt();  // sqrt(met_px*met_px+met_py*met_py);
0529         LogTrace("") << ">>> MET, MET_px, MET_py: " << met_et << ", " << met_px << ", " << met_py << " [GeV]";
0530         double w_et = met_et + muon1.pt();
0531         double w_px = met_px + muon1.px();
0532         double w_py = met_py + muon1.py();
0533         double massT = w_et * w_et - w_px * w_px - w_py * w_py;
0534         massT = (massT > 0) ? sqrt(massT) : 0;
0535         // Acoplanarity cuts
0536         Geom::Phi<double> deltaphi(muon1.phi() - atan2(met_py, met_px));
0537         double acop = M_PI - fabs(deltaphi.value());
0538         // std::cout << " is acp of W candidate less then cut: " << (acop<
0539         // acopCut_) << std::endl;
0540         if (acop > acopCut_)
0541           continue;  // Cut in 2.0
0542         TMass_->Fill(massT);
0543         if (massT < mtMin_ || massT > mtMax_)
0544           continue;  // Cut in (50,200) GeV
0545         // we give the number of W only in the tMass selected but we have a look
0546         // at mass tails to check the QCD background
0547         isW_ = true;
0548         nW++;
0549         nTMass++;
0550       }
0551     }
0552     // if a zGlobal is not selected, look at the dimuonGlobalOneStandAlone and
0553     // dimuonGlobalOneTrack...., used as a control sample for the tracking
0554     // efficiency  and muon system efficiency
0555     if (!(isZGolden2HLT_ || isZGolden1HLT_ || isZGoldenNoIso_)) {
0556       for (unsigned int i = 0; i < nHighPtGlbMu; ++i) {
0557         reco::Muon glbMuon = highPtGlbMuons[i];
0558         const math::XYZTLorentzVector& mu1(glbMuon.p4());
0559         // double pt1= glbMuon.pt();
0560         // checking that the global muon is hlt matched otherwise skip the event
0561         singleTrigFlag1 = IsMuMatchedToHLTMu(glbMuon, HLTMuMatched, maxDeltaR_, maxDPtRel_);
0562         if (!singleTrigFlag1)
0563           continue;
0564         // checking that the global muon is isolated matched otherwise skip the
0565         // event
0566         iso1 = muIso(glbMuon);
0567         isMu1Iso = (iso1 < isoCut03_);
0568         if (!isMu1Iso)
0569           continue;
0570         // look at the standalone muon ....
0571         // stop the loop after 10 cicles....
0572         (nHighPtStaMu > 10) ? nHighPtStaMu = 10 : 1;
0573         for (unsigned int j = 0; j < nHighPtStaMu; ++j) {
0574           reco::Muon staMuon = highPtStaMuons[j];
0575           const math::XYZTLorentzVector& mu2(staMuon.p4());
0576           // double pt2= staMuon.pt();
0577           if (glbMuon.charge() == staMuon.charge())
0578             continue;
0579           math::XYZTLorentzVector pair = mu1 + mu2;
0580           double mass = pair.M();
0581           iso2 = muIso(staMuon);
0582           isMu2Iso = (iso2 < isoCut03_);
0583           LogTrace("") << "\t... isolation value" << iso1 << ", isolated? " << isMu1Iso;
0584           LogTrace("") << "\t... isolation value" << iso2 << ", isolated? " << isMu2Iso;
0585           // requiring theat the global mu is mathed to the HLT  and both are
0586           // isolated
0587           if (isMu2Iso) {
0588             isZGlbSta_ = true;
0589             nGlbSta++;
0590             massGlbSta_->Fill(mass);
0591             highMassGlbSta_->Fill(mass);
0592             /*
0593               if (pt1 > pt2) {
0594               highest_muptGlbSta_ -> Fill (pt1);
0595               lowest_muptGlbSta_ -> Fill (pt2);
0596             } else {
0597               highest_muptGlbSta_ -> Fill (pt2);
0598               lowest_muptGlbSta_ -> Fill (pt1);
0599             }
0600             */
0601           }
0602         }
0603         // look at the tracks....
0604         Handle<TrackCollection> tracks;
0605         if (!ev.getByToken(trackToken_, tracks)) {
0606           // LogError("") << ">>> track collection does not exist !!!";
0607           return;
0608         }
0609         ev.getByToken(trackToken_, tracks);
0610         Handle<CaloTowerCollection> calotower;
0611         if (!ev.getByToken(caloTowerToken_, calotower)) {
0612           // LogError("") << ">>> calotower collection does not exist !!!";
0613           return;
0614         }
0615         ev.getByToken(caloTowerToken_, calotower);
0616         // avoid to loop on more than 5000 trks
0617         size_t nTrk = tracks->size();
0618         (nTrk > 5000) ? nTrk = 5000 : 1;
0619         for (unsigned int j = 0; j < nTrk; j++) {
0620           const reco::Track& tk = tracks->at(j);
0621           if (glbMuon.charge() == tk.charge())
0622             continue;
0623           double pt2 = tk.pt();
0624           double eta = tk.eta();
0625           double dxy = tk.dxy(beamSpotHandle->position());
0626           if (pt2 < ptMuCut_ || fabs(eta) > etaMuCut_ || fabs(dxy) > dxyCut_)
0627             continue;
0628           // assuming that the track is a mu....
0629           math::XYZTLorentzVector mu2(tk.px(), tk.py(), tk.pz(), sqrt((tk.p() * tk.p()) + (0.10566 * 0.10566)));
0630           math::XYZTLorentzVector pair = mu1 + mu2;
0631           double mass = pair.M();
0632           // now requiring that the tracks is isolated.......
0633           iso2 = tkIso(tk, tracks, calotower);
0634           isMu2Iso = (iso2 < isoCut03_);
0635           //          std::cout << "found a track with rel/comb iso: " << iso2
0636           //<< std::endl;
0637           if (isMu2Iso) {
0638             isZGlbTrk_ = true;
0639             nGlbTrk++;
0640             massGlbTrk_->Fill(mass);
0641             highMassGlbTrk_->Fill(mass);
0642             if (!isW_)
0643               continue;
0644             massIsBothGlbTrkThanW_->Fill(mass);
0645             highMassIsBothGlbTrkThanW_->Fill(mass);
0646             /*
0647               if (pt1 > pt2) {
0648               highest_muptGlbTrk_ -> Fill (pt1);
0649               lowest_muptGlbTrk_ -> Fill (pt2);
0650               } else {
0651               highest_muptGlbTrk_ -> Fill (pt2);
0652               lowest_muptGlbTrk_ -> Fill (pt1);
0653             }
0654             */
0655           }
0656         }
0657       }
0658     }
0659   }
0660   if ((hlt_sel || isZGolden2HLT_ || isZGolden1HLT_ || isZGoldenNoIso_ || isZGlbSta_ || isZGlbTrk_ || isW_)) {
0661     nsel++;
0662     LogTrace("") << ">>>> Event ACCEPTED";
0663   } else {
0664     LogTrace("") << ">>>> Event REJECTED";
0665   }
0666   return;
0667 }
0668
0669 // Local Variables:
0670 // show-trailing-whitespace: t
0671 // truncate-lines: t
0672 // End: