Reference Manual

Go to the documentation of this file.
00001 #include "DQM/Physics/src/EwkMuDQM.h"
00002 
00003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00005 #include "DataFormats/Common/interface/Handle.h"
00006 
00007 #include "FWCore/ServiceRegistry/interface/Service.h"
00008 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00009 
00010 #include "DQMServices/Core/interface/DQMStore.h"
00011 #include "DQMServices/Core/interface/MonitorElement.h"
00012 
00013 #include "DataFormats/TrackReco/interface/Track.h"
00014 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00015 
00016 #include "DataFormats/MuonReco/interface/Muon.h"
00017 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
00018 #include "DataFormats/METReco/interface/MET.h"
00019 #include "DataFormats/JetReco/interface/Jet.h"
00020 
00021 #include "DataFormats/GeometryVector/interface/Phi.h"
00022 
00023 #include "FWCore/Common/interface/TriggerNames.h"
00024 #include "FWCore/Framework/interface/Event.h"
00025 #include "DataFormats/Common/interface/TriggerResults.h"
00026 
00027 #include "DataFormats/Common/interface/View.h"
00028   
00029 using namespace edm;
00030 using namespace std;
00031 using namespace reco;
00032 
00033 EwkMuDQM::EwkMuDQM( const ParameterSet & cfg ) :
00034       // Input collections
00035       trigTag_(cfg.getUntrackedParameter<edm::InputTag> ("TrigTag", edm::InputTag("TriggerResults::HLT"))),
00036       muonTag_(cfg.getUntrackedParameter<edm::InputTag> ("MuonTag", edm::InputTag("muons"))),
00037       metTag_(cfg.getUntrackedParameter<edm::InputTag> ("METTag", edm::InputTag("met"))),
00038       metIncludesMuons_(cfg.getUntrackedParameter<bool> ("METIncludesMuons", false)),
00039       jetTag_(cfg.getUntrackedParameter<edm::InputTag> ("JetTag", edm::InputTag("sisCone5CaloJets"))),
00040 
00041       // Main cuts 
00042       muonTrig_(cfg.getUntrackedParameter<std::string> ("MuonTrig", "HLT_Mu9")),
00043       ptCut_(cfg.getUntrackedParameter<double>("PtCut", 25.)),
00044       etaCut_(cfg.getUntrackedParameter<double>("EtaCut", 2.1)),
00045       isRelativeIso_(cfg.getUntrackedParameter<bool>("IsRelativeIso", true)),
00046       isCombinedIso_(cfg.getUntrackedParameter<bool>("IsCombinedIso", false)),
00047       isoCut03_(cfg.getUntrackedParameter<double>("IsoCut03", 0.1)),
00048       mtMin_(cfg.getUntrackedParameter<double>("MtMin", 50.)),
00049       mtMax_(cfg.getUntrackedParameter<double>("MtMax", 200.)),
00050       metMin_(cfg.getUntrackedParameter<double>("MetMin", -999999.)),
00051       metMax_(cfg.getUntrackedParameter<double>("MetMax", 999999.)),
00052       acopCut_(cfg.getUntrackedParameter<double>("AcopCut", 2.)),
00053 
00054       // Muon quality cuts
00055       dxyCut_(cfg.getUntrackedParameter<double>("DxyCut", 0.2)),
00056       normalizedChi2Cut_(cfg.getUntrackedParameter<double>("NormalizedChi2Cut", 10.)),
00057       trackerHitsCut_(cfg.getUntrackedParameter<int>("TrackerHitsCut", 11)),
00058       isAlsoTrackerMuon_(cfg.getUntrackedParameter<bool>("IsAlsoTrackerMuon", true)),
00059 
00060       // Z rejection
00061       ptThrForZ1_(cfg.getUntrackedParameter<double>("PtThrForZ1", 20.)),
00062       ptThrForZ2_(cfg.getUntrackedParameter<double>("PtThrForZ2", 10.)),
00063 
00064       // Top rejection
00065       eJetMin_(cfg.getUntrackedParameter<double>("EJetMin", 999999.)),
00066       nJetMax_(cfg.getUntrackedParameter<int>("NJetMax", 999999))
00067 {
00068 }
00069 
00070 void EwkMuDQM::beginRun(const Run& r, const EventSetup&) {
00071       nall = 0;
00072       nsel = 0;
00073 
00074       nrec = 0; 
00075       niso = 0; 
00076       nhlt = 0; 
00077       nmet = 0;
00078 }
00079 
00080 
00081 void EwkMuDQM::beginJob() {
00082       theDbe = Service<DQMStore>().operator->();
00083       theDbe->setCurrentFolder("Physics/EwkMuDQM");
00084 
00085       init_histograms();
00086 }
00087 
00088 void EwkMuDQM::init_histograms() {
00089 
00090       char chtitle[256] = "";
00091       for (int i=0; i<2; ++i) {
00092             snprintf(chtitle, 255, "Muon transverse momentum (global muon) [GeV]");
00093             pt_before_ = theDbe->book1D("PT_BEFORECUTS",chtitle,100,0.,100.);
00094             pt_after_ = theDbe->book1D("PT_LASTCUT",chtitle,100,0.,100.);
00095 
00096             snprintf(chtitle, 255, "Muon pseudo-rapidity");
00097             eta_before_ = theDbe->book1D("ETA_BEFORECUTS",chtitle,50,-2.5,2.5);
00098             eta_after_ = theDbe->book1D("ETA_LASTCUT",chtitle,50,-2.5,2.5);
00099 
00100             snprintf(chtitle, 255, "Muon transverse distance to beam spot [cm]");
00101             dxy_before_ = theDbe->book1D("DXY_BEFORECUTS",chtitle,100,-0.5,0.5);
00102             dxy_after_ = theDbe->book1D("DXY_LASTCUT",chtitle,100,-0.5,0.5);
00103 
00104             snprintf(chtitle, 255, "Normalized Chi2, inner track fit");
00105             chi2_before_ = theDbe->book1D("CHI2_BEFORECUTS",chtitle,100,0.,100.);
00106             chi2_after_ = theDbe->book1D("CHI2_LASTCUT",chtitle,100,0.,100.);
00107 
00108             snprintf(chtitle, 255, "Number of hits, inner track");
00109             nhits_before_ = theDbe->book1D("NHITS_BEFORECUTS",chtitle,40,-0.5,39.5);
00110             nhits_after_ = theDbe->book1D("NHITS_LASTCUT",chtitle,40,-0.5,39.5);
00111 
00112             snprintf(chtitle, 255, "number Of Valid Muon Hits");
00113             muonhits_before_= theDbe->book1D("MUONHITS_BEFORECUTS",chtitle,40,-0.5,39.5);
00114             muonhits_after_= theDbe->book1D("MUONHITS_LASTCUT",chtitle,40,-0.5,39.5);
00115 
00116             snprintf(chtitle, 255, "Tracker-muon flag (for global muons)");
00117             tkmu_before_ = theDbe->book1D("TKMU_BEFORECUTS",chtitle,2,-0.5,1.5);
00118             tkmu_after_ = theDbe->book1D("TKMU_LASTCUT",chtitle,2,-0.5,1.5);
00119 
00120             snprintf(chtitle, 255, "Quality-muon flag");
00121             goodewkmuon_before_ = theDbe->book1D("GOODEWKMUON_BEFORECUTS",chtitle,2,-0.5,1.5);
00122             goodewkmuon_after_ = theDbe->book1D("GOODEWKMUON_LASTCUT",chtitle,2,-0.5,1.5);
00123 
00124             if (isRelativeIso_) {
00125                   if (isCombinedIso_) {
00126                         snprintf(chtitle, 255, "Relative (combined) isolation variable");
00127                   } else {
00128                         snprintf(chtitle, 255, "Relative (tracker) isolation variable");
00129                   }
00130                   iso_before_ = theDbe->book1D("ISO_BEFORECUTS",chtitle,100, 0., 1.);
00131                   iso_after_ = theDbe->book1D("ISO_LASTCUT",chtitle,100, 0., 1.);
00132             } else {
00133                   if (isCombinedIso_) {
00134                         snprintf(chtitle, 255, "Absolute (combined) isolation variable [GeV]");
00135                   } else {
00136                         snprintf(chtitle, 255, "Absolute (tracker) isolation variable [GeV]");
00137                   }
00138                   iso_before_ = theDbe->book1D("ISO_BEFORECUTS",chtitle,100, 0., 20.);
00139                   iso_after_ = theDbe->book1D("ISO_LASTCUT",chtitle,100, 0., 20.);
00140             }
00141 
00142             snprintf(chtitle, 255, "Trigger response (bit %s)", muonTrig_.data());
00143             trig_before_ = theDbe->book1D("TRIG_BEFORECUTS",chtitle,2,-0.5,1.5);
00144             trig_after_ = theDbe->book1D("TRIG_LASTCUT",chtitle,2,-0.5,1.5);
00145 
00146             snprintf(chtitle, 255, "Transverse mass (%s) [GeV]", metTag_.label().data());
00147             mt_before_ = theDbe->book1D("MT_BEFORECUTS",chtitle,150,0.,300.);
00148             mt_after_ = theDbe->book1D("MT_LASTCUT",chtitle,150,0.,300.);
00149 
00150             snprintf(chtitle, 255, "Missing transverse energy (%s) [GeV]", metTag_.label().data());
00151             met_before_ = theDbe->book1D("MET_BEFORECUTS",chtitle,100,0.,200.);
00152             met_after_ = theDbe->book1D("MET_LASTCUT",chtitle,100,0.,200.);
00153 
00154             snprintf(chtitle, 255, "MU-MET (%s) acoplanarity", metTag_.label().data());
00155             acop_before_ = theDbe->book1D("ACOP_BEFORECUTS",chtitle,50,0.,M_PI);
00156             acop_after_ = theDbe->book1D("ACOP_LASTCUT",chtitle,50,0.,M_PI);
00157 
00158             snprintf(chtitle, 255, "Z rejection: number of muons above %.2f GeV", ptThrForZ1_);
00159             nz1_before_ = theDbe->book1D("NZ1_BEFORECUTS",chtitle,10,-0.5,9.5);
00160             nz1_after_ = theDbe->book1D("NZ1_LASTCUT",chtitle,10,-0.5,9.5);
00161 
00162             snprintf(chtitle, 255, "Z rejection: number of muons above %.2f GeV", ptThrForZ2_);
00163             nz2_before_ = theDbe->book1D("NZ2_BEFORECUTS",chtitle,10,-0.5,9.5);
00164             nz2_after_ = theDbe->book1D("NZ2_LASTCUT",chtitle,10,-0.5,9.5);
00165 
00166             snprintf(chtitle, 255, "Number of jets (%s) above %.2f GeV", jetTag_.label().data(), eJetMin_);
00167             njets_before_ = theDbe->book1D("NJETS_BEFORECUTS",chtitle,10,-0.5,9.5);
00168             njets_after_ = theDbe->book1D("NJETS_LASTCUT",chtitle,10,-0.5,9.5);
00169 
00170             snprintf(chtitle, 255, "DiMuonMass (2 globals)");
00171             dimuonmass_before_= theDbe->book1D("DIMUONMASS_BEFORECUTS",chtitle,100,0,200);
00172             dimuonmass_after_= theDbe->book1D("DIMUONMASS_AFTERZCUTS",chtitle,100,0,200);
00173 
00174             snprintf(chtitle, 255, "DiMuon Mass (global pt + StandAlone pt");
00175             dimuonSAmass_before_= theDbe->book1D("DIMUONSTAMASS_BEFORECUTS",chtitle,100,0,200);
00176             dimuonSAmass_after_= theDbe->book1D("DIMUONSTAMASS_AFTERZCUTS",chtitle,100,0,200);
00177 
00178             snprintf(chtitle, 255, "DiMuon Mass (StandAlone pt + StandAlone pt");
00179             dimuonSASAmass_before_= theDbe->book1D("DIMUONSTASTAMASS_BEFORECUTS",chtitle,100,0,200); 
00180             dimuonSASAmass_after_= theDbe->book1D("DIMUONSTASTAMASS_AFTERZCUTS",chtitle,100,0,200);
00181             
00182             snprintf(chtitle, 255, "Global pt for Muons in Z");
00183             ptmuonZ_after_= theDbe->book1D("PT_AFTERZCUT",chtitle,100,0.,100.);
00184       }
00185 }
00186 
00187 
00188 void EwkMuDQM::endJob() {
00189 }
00190 
00191 void EwkMuDQM::endRun(const Run& r, const EventSetup&) {
00192 
00193 }
00194 
00195 void EwkMuDQM::analyze (const Event & ev, const EventSetup &) {
00196       
00197       // Reset global event selection flags
00198       bool rec_sel = false;
00199       bool iso_sel = false;
00200       bool hlt_sel = false;
00201       bool met_sel = false;
00202       bool all_sel = false;
00203 
00204       // Muon collection
00205       Handle<View<Muon> > muonCollection;
00206       if (!ev.getByLabel(muonTag_, muonCollection)) {
00207         //LogWarning("") << ">>> Muon collection does not exist !!!";
00208         return;
00209       }
00210       unsigned int muonCollectionSize = muonCollection->size();
00211 
00212       // Beam spot
00213       Handle<reco::BeamSpot> beamSpotHandle;
00214       if (!ev.getByLabel(InputTag("offlineBeamSpot"), beamSpotHandle)) {
00215         //LogWarning("") << ">>> No beam spot found !!!";
00216         return;
00217       }
00218   
00219       // MET
00220       double met_px = 0.;
00221       double met_py = 0.;
00222       Handle<View<MET> > metCollection;
00223       if (!ev.getByLabel(metTag_, metCollection)) {
00224         //LogWarning("") << ">>> MET collection does not exist !!!";
00225         return;
00226       }
00227       const MET& met = metCollection->at(0);
00228       met_px = met.px();
00229       met_py = met.py();
00230       if (!metIncludesMuons_) {
00231             for (unsigned int i=0; i<muonCollectionSize; i++) {
00232                   const Muon& mu = muonCollection->at(i);
00233                   if (!mu.isGlobalMuon()) continue;
00234                   met_px -= mu.px();
00235                   met_py -= mu.py();
00236             }
00237       }
00238       double met_et = sqrt(met_px*met_px+met_py*met_py);
00239       LogTrace("") << ">>> MET, MET_px, MET_py: " << met_et << ", " << met_px << ", " << met_py << " [GeV]";
00240       met_before_->Fill(met_et);
00241 
00242       // Trigger
00243       Handle<TriggerResults> triggerResults;
00244       if (!ev.getByLabel(trigTag_, triggerResults)) {
00245         //LogWarning("") << ">>> TRIGGER collection does not exist !!!";
00246         return;
00247       }
00248       const edm::TriggerNames & trigNames = ev.triggerNames(*triggerResults);
00249       bool trigger_fired = false;
00250       /*
00251       for (unsigned int i=0; i<triggerResults->size(); i++) {
00252             if (triggerResults->accept(i)) {
00253                   LogTrace("") << "Accept by: " << i << ", Trigger: " << trigNames.triggerName(i);
00254             }
00255       }
00256       */
00257 
00258       // the following gives error on CRAFT08 data where itrig1=19 (vector index out of range)
00259       /*
00260       int itrig1 = trigNames.triggerIndex(muonTrig_);
00261       if (triggerResults->accept(itrig1)) trigger_fired = true;
00262       */
00263       //suggested replacement: lm250909
00264       for (unsigned int i=0; i<triggerResults->size(); i++) {
00265         std::string trigName = trigNames.triggerName(i);
00266         if ( trigName == muonTrig_ && triggerResults->accept(i)) trigger_fired = true;
00267       }
00268 
00269 
00270       LogTrace("") << ">>> Trigger bit: " << trigger_fired << " (" << muonTrig_ << ")";
00271       trig_before_->Fill(trigger_fired);
00272 
00273       // Loop to reject/control Z->mumu is done separately
00274       unsigned int nmuonsForZ1 = 0;
00275       unsigned int nmuonsForZ2 = 0;
00276       for (unsigned int i=0; i<muonCollectionSize; i++) {
00277             const Muon& mu = muonCollection->at(i);
00278             if (!mu.isGlobalMuon()) continue;
00279             double pt = mu.pt();
00280             if (pt>ptThrForZ1_) nmuonsForZ1++;
00281             if (pt>ptThrForZ2_) nmuonsForZ2++;
00282 
00283             for (unsigned int j=0; j<muonCollectionSize; j++) {
00284                   if (i==j) continue;
00285                   const Muon& mu2 = muonCollection->at(j);
00286                  // Glb + Glb  
00287                  if (mu2.isGlobalMuon() && j>i ){
00288                          const math::XYZTLorentzVector ZRecoGlb (mu.px()+mu2.px(), mu.py()+mu2.py() , mu.pz()+mu2.pz(), mu.p()+mu2.p());
00289                          dimuonmass_before_->Fill(ZRecoGlb.mass());
00290                  }
00291                   // Glb + Standalone 
00292                  if (mu2.isStandAloneMuon()){
00293                          const math::XYZTLorentzVector ZRecoSta (mu2.outerTrack()->px()+mu.px(), mu.py()+mu.outerTrack()->py() , mu.pz()+mu2.outerTrack()->pz(), mu.p()+mu2.outerTrack()->p());
00294                          dimuonSAmass_before_->Fill(ZRecoSta.mass());
00295                  }
00296                   // Standalone + Standalone 
00297                  if (mu2.isStandAloneMuon() && j>i){
00298                          const math::XYZTLorentzVector ZRecoStaSta (mu2.outerTrack()->px()+mu.outerTrack()->px(), mu.outerTrack()->py()+mu.outerTrack()->py() , mu.outerTrack()->pz()+mu2.outerTrack()->pz(), mu.outerTrack()->p()+mu2.outerTrack()->p());
00299                          dimuonSASAmass_before_->Fill(ZRecoStaSta.mass());
00300                  }
00301             }
00302       
00303       }
00304       LogTrace("") << "> Z rejection: muons above " << ptThrForZ1_ << " [GeV]: " << nmuonsForZ1;
00305       LogTrace("") << "> Z rejection: muons above " << ptThrForZ2_ << " [GeV]: " << nmuonsForZ2;
00306       nz1_before_->Fill(nmuonsForZ1);
00307       nz2_before_->Fill(nmuonsForZ2);
00308       
00309       // Jet collection
00310       Handle<View<Jet> > jetCollection;
00311       if (!ev.getByLabel(jetTag_, jetCollection)) {
00312         //LogError("") << ">>> JET collection does not exist !!!";
00313         return;
00314       }
00315       unsigned int jetCollectionSize = jetCollection->size();
00316       int njets = 0;
00317       for (unsigned int i=0; i<jetCollectionSize; i++) {
00318             const Jet& jet = jetCollection->at(i);
00319             if (jet.et()>eJetMin_) njets++;
00320       }
00321       LogTrace("") << ">>> Total number of jets: " << jetCollectionSize;
00322       LogTrace("") << ">>> Number of jets above " << eJetMin_ << " [GeV]: " << njets;
00323       njets_before_->Fill(njets);
00324 
00325       // Start counting
00326       nall++;
00327 
00328       // Histograms per event shouldbe done only once, so keep track of them
00329       bool hlt_hist_done = false;
00330       bool met_hist_done = false;
00331       bool nz1_hist_done = false;
00332       bool nz2_hist_done = false;
00333       bool njets_hist_done = false;
00334 
00335       // Central W->mu nu selection criteria
00336       const int NFLAGS = 13;
00337       bool muon_sel[NFLAGS];
00338       bool muon4Z=true;
00339 
00340       for (unsigned int i=0; i<muonCollectionSize; i++) {
00341             for (int j=0; j<NFLAGS; ++j) {
00342                   muon_sel[j] = false;
00343             }
00344 
00345             const Muon& mu = muonCollection->at(i);
00346             if (!mu.isGlobalMuon()) continue;
00347             if (mu.globalTrack().isNull()) continue;
00348             if (mu.innerTrack().isNull()) continue;
00349 
00350             LogTrace("") << "> Wsel: processing muon number " << i << "...";
00351             reco::TrackRef gm = mu.globalTrack();
00352             reco::TrackRef tk = mu.innerTrack();
00353 
00354             // Pt,eta cuts
00355             double pt = mu.pt();
00356             double eta = mu.eta();
00357             LogTrace("") << "\t... pt, eta: " << pt << " [GeV], " << eta;;
00358             if (pt>ptCut_) muon_sel[0] = true; 
00359             if (fabs(eta)<etaCut_) muon_sel[1] = true; 
00360             if (pt<ptThrForZ1_) { muon4Z = false;}
00361 
00362             // d0, chi2, nhits quality cuts
00363             double dxy = tk->dxy(beamSpotHandle->position());
00364             double normalizedChi2 = gm->normalizedChi2();
00365             double trackerHits = tk->numberOfValidHits();
00366             double validmuonhits=gm->hitPattern().numberOfValidMuonHits();
00367             LogTrace("") << "\t... dxy, normalizedChi2, trackerHits, isTrackerMuon?: " << dxy << " [cm], " << normalizedChi2 << ", " << trackerHits << ", " << mu.isTrackerMuon();
00368             if (fabs(dxy)<dxyCut_) muon_sel[2] = true; 
00369 //            if (normalizedChi2<normalizedChi2Cut_) muon_sel[3] = true; 
00370             if (muon::isGoodMuon(mu,muon::GlobalMuonPromptTight)) muon_sel[3] = true;
00371             if (trackerHits>=trackerHitsCut_) muon_sel[4] = true; 
00372             if (mu.isTrackerMuon()) muon_sel[5] = true; 
00373 
00374             pt_before_->Fill(pt);
00375             eta_before_->Fill(eta);
00376             dxy_before_->Fill(dxy);
00377             chi2_before_->Fill(normalizedChi2);
00378             nhits_before_->Fill(trackerHits);
00379             muonhits_before_->Fill(validmuonhits);
00380             tkmu_before_->Fill(mu.isTrackerMuon());
00381 
00382             bool quality = muon_sel[4]*muon_sel[2]* muon_sel[3]* muon_sel[5];
00383             goodewkmuon_before_->Fill(quality);
00384 
00385             // Isolation cuts
00386             double isovar = mu.isolationR03().sumPt;
00387             if (isCombinedIso_) {
00388                   isovar += mu.isolationR03().emEt;
00389                   isovar += mu.isolationR03().hadEt;
00390             }
00391             if (isRelativeIso_) isovar /= pt;
00392             if (isovar<isoCut03_) muon_sel[6] = true; 
00393             if (isovar>=isoCut03_) { muon4Z = false;}
00394 
00395             LogTrace("") << "\t... isolation value" << isovar <<", isolated? " << muon_sel[6];
00396             iso_before_->Fill(isovar);
00397 
00398 
00399             // HLT (not mtched to muon for the time being)
00400             if (trigger_fired) muon_sel[7] = true; 
00401             else { muon4Z = false;}
00402 
00403             // MET/MT cuts
00404             double w_et = met_et+mu.pt();
00405             double w_px = met_px+mu.px();
00406             double w_py = met_py+mu.py();
00407             
00408             double massT = w_et*w_et - w_px*w_px - w_py*w_py;
00409             massT = (massT>0) ? sqrt(massT) : 0;
00410 
00411             LogTrace("") << "\t... W mass, W_et, W_px, W_py: " << massT << ", " << w_et << ", " << w_px << ", " << w_py << " [GeV]";
00412             if (massT>mtMin_ && massT<mtMax_) muon_sel[8] = true; 
00413             mt_before_->Fill(massT);
00414             if (met_et>metMin_ && met_et<metMax_) muon_sel[9] = true; 
00415 
00416             // Acoplanarity cuts
00417             Geom::Phi<double> deltaphi(mu.phi()-atan2(met_py,met_px));
00418             double acop = deltaphi.value();
00419             if (acop<0) acop = - acop;
00420             acop = M_PI - acop;
00421             LogTrace("") << "\t... acoplanarity: " << acop;
00422             if (acop<acopCut_) muon_sel[10] = true; 
00423             acop_before_->Fill(acop);
00424 
00425             // Remaining flags (from global event information)
00426             if (nmuonsForZ1<1 || nmuonsForZ2<2) muon_sel[11] = true; 
00427             if (njets<=nJetMax_) muon_sel[12] = true; 
00428 
00429             // Collect necessary flags "per muon"
00430             int flags_passed = 0;
00431             bool rec_sel_this = true;
00432             bool iso_sel_this = true;
00433             bool hlt_sel_this = true;
00434             bool met_sel_this = true;
00435             bool all_sel_this = true;
00436             for (int j=0; j<NFLAGS; ++j) {
00437                   if (muon_sel[j]) flags_passed += 1;
00438                   if (j<6 && !muon_sel[j]) rec_sel_this = false;
00439                   if (j<7 && !muon_sel[j]) iso_sel_this = false;
00440                   if (j<8 && !muon_sel[j]) hlt_sel_this = false;
00441                   if (j<11 && !muon_sel[j]) met_sel_this = false;
00442                   if (!muon_sel[j]) all_sel_this = false;
00443             }
00444 
00445             // "rec" => pt,eta and quality cuts are satisfied
00446             if (rec_sel_this) rec_sel = true;
00447             // "iso" => "rec" AND "muon is isolated"
00448             if (iso_sel_this) iso_sel = true;
00449             // "hlt" => "iso" AND "event is triggered"
00450             if (hlt_sel_this) hlt_sel = true;
00451             // "met" => "hlt" AND "MET/MT and acoplanarity cuts"
00452             if (met_sel_this) met_sel = true;
00453             // "all" => "met" AND "Z/top rejection cuts"
00454             if (all_sel_this) all_sel = true;
00455 
00456             // Do N-1 histograms now (and only once for global event quantities)
00457             if (flags_passed >= (NFLAGS-1)) {
00458                   if (!muon_sel[0] || flags_passed==NFLAGS) 
00459                         pt_after_->Fill(pt);
00460                   if (!muon_sel[1] || flags_passed==NFLAGS) 
00461                         eta_after_->Fill(eta);
00462                   if (!muon_sel[2] || flags_passed==NFLAGS) 
00463                         dxy_after_->Fill(dxy);
00464                   if (!muon_sel[3] || flags_passed==NFLAGS){ 
00465                         chi2_after_->Fill(normalizedChi2);
00466                         muonhits_after_->Fill(validmuonhits);
00467                   }
00468                   if (!muon_sel[4] || flags_passed==NFLAGS) 
00469                         nhits_after_->Fill(trackerHits);
00470                   if (!muon_sel[5] || flags_passed==NFLAGS) 
00471                         tkmu_after_->Fill(mu.isTrackerMuon());
00472                   if (!muon_sel[6] || flags_passed==NFLAGS) 
00473                         iso_after_->Fill(isovar);
00474                   if (!muon_sel[2]||!muon_sel[3] || !muon_sel[4] || !muon_sel[5] || flags_passed==NFLAGS) 
00475                         goodewkmuon_after_->Fill(quality);
00476                   if (!muon_sel[7] || flags_passed==NFLAGS) 
00477                         if (!hlt_hist_done) trig_after_->Fill(trigger_fired);
00478                         hlt_hist_done = true;
00479                   if (!muon_sel[8] || flags_passed==NFLAGS) 
00480                         mt_after_->Fill(massT);
00481                   if (!muon_sel[9] || flags_passed==NFLAGS) 
00482                         if (!met_hist_done) met_after_->Fill(met_et);
00483                         met_hist_done = true;
00484                   if (!muon_sel[10] || flags_passed==NFLAGS) 
00485                         acop_after_->Fill(acop);
00486                   if (!muon_sel[11] || flags_passed==NFLAGS) 
00487                         if (!nz1_hist_done) nz1_after_->Fill(nmuonsForZ1);
00488                         nz1_hist_done = true;
00489                   if (!muon_sel[11] || flags_passed==NFLAGS) 
00490                         if (!nz2_hist_done) nz2_after_->Fill(nmuonsForZ2);
00491                         nz2_hist_done = true;
00492                   if (!muon_sel[12] || flags_passed==NFLAGS) 
00493                         if (!njets_hist_done) njets_after_->Fill(njets);
00494                         njets_hist_done = true;
00495             }
00496 
00497 
00498             // The cases in which the event is rejected as a Z are considered independently:
00499             if ( muon4Z &&  !muon_sel[11]){
00500                    // Plots for 2 muons       
00501                    bool usedMuon=false;
00502                    for (unsigned int j=0; j<muonCollectionSize; j++) {
00503                          if (i==j) continue;
00504                          const Muon& mu2 = muonCollection->at(j);
00505                                     double pt2 = mu2.pt();
00506                                     double isovar2 = mu2.isolationR03().sumPt;
00507                                     if (isCombinedIso_) {
00508                                           isovar2 += mu2.isolationR03().emEt;
00509                                           isovar2 += mu2.isolationR03().hadEt;
00510                                     }
00511                                     if (isRelativeIso_) isovar2 /= pt2;
00512 
00513                           if (pt2<=ptThrForZ1_ || isovar2>=isoCut03_) continue;
00514                   
00515                   // Glb + Glb  
00516                              if (mu2.isGlobalMuon() && j>i ){
00517                                const math::XYZTLorentzVector ZRecoGlb (mu.px()+mu2.px(), mu.py()+mu2.py() , mu.pz()+mu2.pz(), mu.p()+mu2.p());
00518                                dimuonmass_after_->Fill(ZRecoGlb.mass());
00519                                if(!usedMuon){ptmuonZ_after_->Fill(mu.pt()); usedMuon=true;}
00520                              }
00521                   // Glb + Standalone 
00522                              if (mu2.isStandAloneMuon()){
00523                               const math::XYZTLorentzVector ZRecoSta (mu2.outerTrack()->px()+mu.px(), mu.py()+mu.outerTrack()->py() , mu.pz()+mu2.outerTrack()->pz(), mu.p()+mu2.outerTrack()->p());
00524                               dimuonSAmass_after_->Fill(ZRecoSta.mass());
00525                              }
00526                   // Standalone + Standalone 
00527                              if (mu2.isStandAloneMuon() && j>i){
00528                               const math::XYZTLorentzVector ZRecoStaSta (mu2.outerTrack()->px()+mu.outerTrack()->px(), mu.outerTrack()->py()+mu.outerTrack()->py() , mu.outerTrack()->pz()+mu2.outerTrack()->pz(), mu.outerTrack()->p()+mu2.outerTrack()->p());
00529                               dimuonSASAmass_after_->Fill(ZRecoStaSta.mass());
00530                              }
00531             }
00532 
00533 
00534 
00535       }
00536 
00537 
00538 
00539 
00540       }
00541 
00542       return;
00543 
00544 }
00545
/data/refman/pasoursint/CMSSW_4_1_8_patch13/src/DQM/Physics/src/EwkMuDQM.cc
