Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-09-21 02:12:24 

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