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 #include "DataFormats/VertexReco/interface/Vertex.h"
00016 
00017 #include "DataFormats/MuonReco/interface/Muon.h"
00018 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
00019 #include "DataFormats/METReco/interface/MET.h"
00020 #include "DataFormats/JetReco/interface/Jet.h"
00021 #include "DataFormats/EgammaCandidates/interface/Photon.h" 
00022 
00023 #include "DataFormats/GeometryVector/interface/Phi.h"
00024 
00025 #include "FWCore/Common/interface/TriggerNames.h"
00026 #include "FWCore/Framework/interface/Event.h"
00027 #include "DataFormats/Common/interface/TriggerResults.h"
00028 
00029 #include "DataFormats/Common/interface/View.h"
00030  
00031 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
00032  
00033 using namespace edm;
00034 using namespace std;
00035 using namespace reco;
00036 
00037 EwkMuDQM::EwkMuDQM( const ParameterSet & cfg ) :
00038       // Input collections
00039       trigTag_      (cfg.getUntrackedParameter<edm::InputTag> ("TrigTag", edm::InputTag("TriggerResults::HLT"))),
00040       muonTag_      (cfg.getUntrackedParameter<edm::InputTag> ("MuonTag", edm::InputTag("muons"))),
00041       metTag_       (cfg.getUntrackedParameter<edm::InputTag> ("METTag", edm::InputTag("pfmet"))),
00042       jetTag_       (cfg.getUntrackedParameter<edm::InputTag> ("JetTag", edm::InputTag("ak5PFJets"))),
00043       phoTag_       (cfg.getUntrackedParameter<edm::InputTag> ("phoTag", edm::InputTag("photons"))),
00044       pfPhoTag_     (cfg.getUntrackedParameter<edm::InputTag> ("pfPhoTag", edm::InputTag("pfPhotonTranslator","pfPhot"))),
00045       vertexTag_    (cfg.getUntrackedParameter<edm::InputTag> ("VertexTag", edm::InputTag("offlinePrimaryVertices"))),
00046       trigPathNames_(cfg.getUntrackedParameter<std::vector <std::string> >("TrigPathNames")),           
00047 
00048       // Muon quality cuts
00049       isAlsoTrackerMuon_(cfg.getUntrackedParameter<bool>("IsAlsoTrackerMuon", true)),  // Glb muon also tracker muon 
00050       dxyCut_           (cfg.getUntrackedParameter<double>("DxyCut", 0.2)),            // dxy < 0.2 cm 
00051       normalizedChi2Cut_(cfg.getUntrackedParameter<double>("NormalizedChi2Cut", 10.)), // chi2/ndof (of global fit) <10.0
00052       trackerHitsCut_   (cfg.getUntrackedParameter<int>("TrackerHitsCut", 11)),        // Tracker Hits >10 
00053       pixelHitsCut_     (cfg.getUntrackedParameter<int>("PixelHitsCut", 1)),           // Pixel Hits >0
00054       muonHitsCut_      (cfg.getUntrackedParameter<int>("MuonHitsCut", 1)),            // Valid Muon Hits >0 
00055       nMatchesCut_      (cfg.getUntrackedParameter<int>("NMatchesCut", 2)),            // At least 2 Chambers with matches 
00056 
00057       // W-boson cuts 
00058       isRelativeIso_(cfg.getUntrackedParameter<bool>("IsRelativeIso", true)),
00059       isCombinedIso_(cfg.getUntrackedParameter<bool>("IsCombinedIso", false)),
00060       isoCut03_     (cfg.getUntrackedParameter<double>("IsoCut03", 0.1)),
00061       acopCut_      (cfg.getUntrackedParameter<double>("AcopCut", 999.)),
00062       metMin_       (cfg.getUntrackedParameter<double>("MetMin", -999999.)),
00063       metMax_       (cfg.getUntrackedParameter<double>("MetMax", 999999.)),
00064       mtMin_        (cfg.getUntrackedParameter<double>("MtMin", 50.)),
00065       mtMax_        (cfg.getUntrackedParameter<double>("MtMax", 200.)),
00066       ptCut_        (cfg.getUntrackedParameter<double>("PtCut", 20.)),
00067       etaCut_       (cfg.getUntrackedParameter<double>("EtaCut", 2.4)),
00068 
00069       // Z rejection
00070       ptThrForZ1_  (cfg.getUntrackedParameter<double>("PtThrForZ1", 20.)),
00071       ptThrForZ2_  (cfg.getUntrackedParameter<double>("PtThrForZ2", 10.)),
00072 
00073       // Z selection
00074       dimuonMassMin_(cfg.getUntrackedParameter<double>("dimuonMassMin", 80.)),
00075       dimuonMassMax_(cfg.getUntrackedParameter<double>("dimuonMassMax", 120.)), 
00076 
00077       // Top rejection
00078       eJetMin_     (cfg.getUntrackedParameter<double>("EJetMin", 999999.)),
00079       nJetMax_     (cfg.getUntrackedParameter<int>("NJetMax", 999999)), 
00080 
00081       // Photon cuts 
00082       ptThrForPhoton_(cfg.getUntrackedParameter<double>("ptThrForPhoton",5.)),
00083       nPhoMax_(cfg.getUntrackedParameter<int>("nPhoMax", 999999)) 
00084 {
00085   isValidHltConfig_ = false;
00086 
00087   theDbe = Service<DQMStore>().operator->();
00088   theDbe->setCurrentFolder("Physics/EwkMuDQM");
00089   init_histograms();
00090 
00091 }
00092 
00093 void EwkMuDQM::beginRun(const Run& iRun, const EventSetup& iSet) {
00094       nall = 0;
00095       nsel = 0;
00096       nz   = 0;
00097 
00098       nrec = 0; 
00099       niso = 0; 
00100       nhlt = 0; 
00101       nmet = 0;
00102 
00103      // passed as parameter to HLTConfigProvider::init(), not yet used
00104      bool isConfigChanged = false;
00105      // isValidHltConfig_ used to short-circuit analyze() in case of problems
00106      isValidHltConfig_ = hltConfigProvider_.init( iRun, iSet, "HLT", isConfigChanged );
00107 
00108 }
00109 
00110 void EwkMuDQM::beginJob() {
00111 
00112 }
00113 
00114 void EwkMuDQM::init_histograms() {
00115 
00116   char chtitle[256] = "";
00117 
00118   pt_before_ = theDbe->book1D("PT_BEFORECUTS","Muon transverse momentum (global muon) [GeV]",100,0.,100.);
00119   pt_after_ = theDbe->book1D("PT_AFTERWCUTS","Muon transverse momentum (global muon) [GeV]",100,0.,100.);
00120 
00121   eta_before_ = theDbe->book1D("ETA_BEFORECUTS","Muon pseudo-rapidity",50,-2.5,2.5);
00122   eta_after_ = theDbe->book1D("ETA_AFTERWCUTS","Muon pseudo-rapidity",50,-2.5,2.5);
00123 
00124   dxy_before_ = theDbe->book1D("DXY_BEFORECUTS","Muon transverse distance to beam spot [cm]",1000,-0.5,0.5);
00125   dxy_after_ = theDbe->book1D("DXY_AFTERWCUTS","Muon transverse distance to beam spot [cm]",1000,-0.5,0.5);
00126 
00127   goodewkmuon_before_ = theDbe->book1D("GOODEWKMUON_BEFORECUTS","Quality-muon flag",2,-0.5,1.5);
00128   goodewkmuon_after_ = theDbe->book1D("GOODEWKMUON_AFTERWCUTS","Quality-muon flag",2,-0.5,1.5);
00129 
00130   if (isRelativeIso_) {
00131     if (isCombinedIso_) {
00132       iso_before_ = theDbe->book1D("ISO_BEFORECUTS","Relative (combined) isolation variable",100, 0., 1.);
00133       iso_after_ = theDbe->book1D("ISO_AFTERWCUTS","Relative (combined) isolation variable",100, 0., 1.);
00134     } else {
00135       iso_before_ = theDbe->book1D("ISO_BEFORECUTS","Relative (tracker) isolation variable",100, 0., 1.);
00136       iso_after_ = theDbe->book1D("ISO_AFTERWCUTS","Relative (tracker) isolation variable",100, 0., 1.);
00137     }
00138   } else {
00139     if (isCombinedIso_) {
00140       iso_before_ = theDbe->book1D("ISO_BEFORECUTS","Absolute (combined) isolation variable [GeV]",100, 0., 20.);
00141       iso_after_ = theDbe->book1D("ISO_AFTERWCUTS","Absolute (combined) isolation variable [GeV]",100, 0., 20.);
00142     } else {
00143       iso_before_ = theDbe->book1D("ISO_BEFORECUTS","Absolute (tracker) isolation variable [GeV]",100, 0., 20.);
00144       iso_after_ = theDbe->book1D("ISO_AFTERWCUTS","Absolute (tracker) isolation variable [GeV]",100, 0., 20.);
00145     }
00146   }
00147 
00148   trig_before_ = theDbe->book1D("TRIG_BEFORECUTS","Trigger response (boolean of muon triggers)",2,-0.5,1.5);
00149   trig_after_ = theDbe->book1D("TRIG_AFTERWCUTS","Trigger response (boolean of muon triggers)",2,-0.5,1.5);
00150 
00151   snprintf(chtitle, 255, "Transverse mass (%s) [GeV]", metTag_.label().data());
00152   mt_before_ = theDbe->book1D("MT_BEFORECUTS",chtitle,150,0.,300.);
00153   mt_after_ = theDbe->book1D("MT_AFTERWCUTS",chtitle,150,0.,300.);
00154 
00155   snprintf(chtitle, 255, "Missing transverse energy (%s) [GeV]", metTag_.label().data());
00156   met_before_ = theDbe->book1D("MET_BEFORECUTS",chtitle,100,0.,200.);
00157   met_after_ = theDbe->book1D("MET_AFTERWCUTS",chtitle,100,0.,200.);
00158   met_afterZ_ = theDbe->book1D("MET_AFTERZCUTS",chtitle,100,0.,200.);
00159 
00160   snprintf(chtitle, 255, "MU-MET (%s) acoplanarity", metTag_.label().data());
00161   acop_before_ = theDbe->book1D("ACOP_BEFORECUTS",chtitle,50,0.,M_PI);
00162   acop_after_ = theDbe->book1D("ACOP_AFTERWCUTS",chtitle,50,0.,M_PI);
00163 
00164   snprintf(chtitle, 255, "Z selection: muons above %.2f GeV", ptThrForZ1_);
00165   n_zselPt1thr_ = theDbe->book1D("NZSELPT1THR",chtitle,10,-0.5,9.5);
00166   snprintf(chtitle, 255, "Z selection: muons above %.2f GeV", ptThrForZ2_);
00167   n_zselPt2thr_ = theDbe->book1D("NZSELPT2THR",chtitle,10,-0.5,9.5);
00168 
00169   snprintf(chtitle, 255, "Number of jets (%s) above %.2f GeV", jetTag_.label().data(), eJetMin_);
00170   njets_before_ = theDbe->book1D("NJETS_BEFORECUTS",chtitle,16,-0.5,15.5);
00171   njets_after_ = theDbe->book1D("NJETS_AFTERWCUTS",chtitle,16,-0.5,15.5);
00172   njets_afterZ_ = theDbe->book1D("NJETS_AFTERZCUTS",chtitle,16,-0.5,15.5);
00173 
00174   leadingjet_pt_before_ = theDbe->book1D("LEADINGJET_PT_BEFORECUTS","Leading Jet transverse momentum",300,0.,300.);
00175   leadingjet_pt_after_ = theDbe->book1D("LEADINGJET_PT_AFTERWCUTS","Leading Jet transverse momentum",300,0.,300.);
00176   leadingjet_pt_afterZ_ = theDbe->book1D("LEADINGJET_PT_AFTERZCUTS","Leading Jet transverse momentum",300,0.,300.);
00177 
00178   leadingjet_eta_before_ = theDbe->book1D("LEADINGJET_ETA_BEFORECUTS","Leading Jet pseudo-rapidity",50,-2.5,2.5);
00179   leadingjet_eta_after_ = theDbe->book1D("LEADINGJET_ETA_AFTERWCUTS","Leading Jet pseudo-rapidity",50,-2.5,2.5);
00180   leadingjet_eta_afterZ_ = theDbe->book1D("LEADINGJET_ETA_AFTERZCUTS","Leading Jet pseudo-rapidity",50,-2.5,2.5);
00181 
00184   //ptPlus_before_ = theDbe->book1D("PTPLUS_BEFORE_CUTS","Muon+ transverse momentum before cuts [GeV]",100,0.,100.);
00185   //ptMinus_before_ = theDbe->book1D("PTMINUS_BEFORE_CUTS","Muon- transverse momentum before cuts [GeV]",100,0.,100.);
00186   //ptPlus_afterW_ = theDbe->book1D("PTPLUS_AFTERW_CUTS","Muon+ transverse momentum after W cuts [GeV]",100,0.,100.);
00187   //ptMinus_afterW_ = theDbe->book1D("PTMINUS_AFTERW_CUTS","Muon- transverse momentum after W cuts [GeV]",100,0.,100.);
00188   //ptPlus_afterZ_ = theDbe->book1D("PTPLUS_AFTERZ_CUTS","Muon+ transverse momentum after Z cuts [GeV]",100,0.,100.);
00189   //ptMinus_afterZ_ = theDbe->book1D("PTMINUS_AFTERZ_CUTS","Muon- transverse momentum after Z cuts [GeV]",100,0.,100.);
00190   ptDiffPM_before_ = theDbe->book1D("PTDIFFPM_BEFORE_CUTS","pt(Muon+)-pt(Muon-) after Z cuts [GeV]",200,-100.,100.);
00191   ptDiffPM_afterZ_ = theDbe->book1D("PTDIFFPM_AFTERZ_CUTS","pt(Muon+)-pt(Muon-) after Z cuts [GeV]",200,-100.,100.);
00192 
00195   pt1_afterZ_ = theDbe->book1D("PT1_AFTERZCUTS","Muon transverse momentum (global muon) [GeV]",100,0.,100.);
00196   eta1_afterZ_ = theDbe->book1D("ETA1_AFTERZCUTS","Muon pseudo-rapidity",50,-2.5,2.5);
00197   dxy1_afterZ_ = theDbe->book1D("DXY1_AFTERZCUTS","Muon transverse distance to beam spot [cm]",1000,-0.5,0.5);
00198   goodewkmuon1_afterZ_ = theDbe->book1D("GOODEWKMUON1_AFTERZCUTS","Quality-muon flag",2,-0.5,1.5);
00199 
00200   if (isRelativeIso_) {
00201     if (isCombinedIso_) {
00202       iso1_afterZ_ = theDbe->book1D("ISO1_AFTERZCUTS","Relative (combined) isolation variable",100, 0., 1.);
00203       iso2_afterZ_ = theDbe->book1D("ISO2_AFTERZCUTS","Relative (combined) isolation variable",100, 0., 1.);
00204     } else {
00205       iso1_afterZ_ = theDbe->book1D("ISO1_AFTERZCUTS","Relative (tracker) isolation variable",100, 0., 1.);
00206       iso2_afterZ_ = theDbe->book1D("ISO2_AFTERZCUTS","Relative (tracker) isolation variable",100, 0., 1.);
00207     }
00208   } else {
00209     if (isCombinedIso_) {
00210       iso1_afterZ_ = theDbe->book1D("ISO1_AFTERZCUTS","Absolute (combined) isolation variable [GeV]",100, 0., 20.);
00211       iso2_afterZ_ = theDbe->book1D("ISO2_AFTERZCUTS","Absolute (combined) isolation variable [GeV]",100, 0., 20.);
00212     } else {
00213       iso1_afterZ_ = theDbe->book1D("ISO1_AFTERZCUTS","Absolute (tracker) isolation variable [GeV]",100, 0., 20.);
00214       iso2_afterZ_ = theDbe->book1D("ISO2_AFTERZCUTS","Absolute (tracker) isolation variable [GeV]",100, 0., 20.);
00215     }
00216   }
00217 
00218   pt2_afterZ_ = theDbe->book1D("PT2_AFTERZCUTS","Muon transverse momentum (global muon) [GeV]",100,0.,100.);
00219   eta2_afterZ_ = theDbe->book1D("ETA2_AFTERZCUTS","Muon pseudo-rapidity",50,-2.5,2.5);
00220   dxy2_afterZ_ = theDbe->book1D("DXY2_AFTERZCUTS","Muon transverse distance to beam spot [cm]",1000,-0.5,0.5);
00221   goodewkmuon2_afterZ_ = theDbe->book1D("GOODEWKMUON2_AFTERZCUTS","Quality-muon flag",2,-0.5,1.5);
00222   ztrig_afterZ_ = theDbe->book1D("ZTRIG_AFTERZCUTS","Trigger response (boolean of muon triggers)",2,-0.5,1.5); 
00223   dimuonmass_before_= theDbe->book1D("DIMUONMASS_BEFORECUTS","DiMuonMass (2 globals)",100,0,200);
00224   dimuonmass_afterZ_= theDbe->book1D("DIMUONMASS_AFTERZCUTS","DiMuonMass (2 globals)",100,0,200);
00225   npvs_before_ = theDbe->book1D("NPVs_BEFORECUTS","Number of Valid Primary Vertices",51,-0.5,50.5);
00226   npvs_after_ = theDbe->book1D("NPVs_AFTERWCUTS","Number of Valid Primary Vertices",51,-0.5,50.5);
00227   npvs_afterZ_ = theDbe->book1D("NPVs_AFTERZCUTS","Number of Valid Primary Vertices",51,-0.5,50.5);
00228   muoncharge_before_ = theDbe->book1D("MUONCHARGE_BEFORECUTS","Muon Charge",3,-1.5,1.5);
00229   muoncharge_after_ = theDbe->book1D("MUONCHARGE_AFTERWCUTS","Muon Charge",3,-1.5,1.5);
00230   muoncharge_afterZ_ = theDbe->book1D("MUONCHARGE_AFTERZCUTS","Muon Charge",3,-1.5,1.5);
00231 
00232   // Adding these to replace the NZ ones (more useful, since they are more general?)     
00233   nmuons_ = theDbe->book1D("NMuons","Number of muons in the event",10,-0.5,9.5);
00234   ngoodmuons_ = theDbe->book1D("NGoodMuons","Number of muons passing the quality criteria",10,-0.5,9.5);
00235 
00236   nph_ = theDbe->book1D("nph","Number of photons in the event",20,0.,20.); 
00237   //npfph_ = theDbe->book1D("npfph","Number of PF photons in the event",20,0.,20.); 
00238   phPt_ = theDbe->book1D("phPt","Photon transverse momentum [GeV]",1000,0.,1000.);
00239   //pfphPt_ = theDbe->book1D("pfphPt","PF Photon transverse momentum [GeV]",1000,0.,1000.); 
00240   snprintf(chtitle, 255, "Photon pseudorapidity (pT>%4.1f)",ptThrForPhoton_);
00241   phEta_ = theDbe->book1D("phEta",chtitle,100,-2.5,2.5); 
00242   //pfphEta_ = theDbe->book1D("pfphEta","PF Photon pseudorapidity",100,-2.5,2.5); 
00243 
00244 }
00245 
00246 
00247 void EwkMuDQM::endJob() {
00248 }
00249 
00250 void EwkMuDQM::endRun(const Run& r, const EventSetup& iSet) {
00251 
00252 }
00253 
00254 void EwkMuDQM::analyze (const Event & ev, const EventSetup & iSet) {
00255       
00256       // Muon collection
00257       Handle<View<Muon> > muonCollection;
00258       if (!ev.getByLabel(muonTag_, muonCollection)) {
00259         //LogWarning("") << ">>> Muon collection does not exist !!!";
00260         return;
00261       }
00262       unsigned int muonCollectionSize = muonCollection->size();
00263 
00264       // Beam spot
00265       Handle<reco::BeamSpot> beamSpotHandle;
00266       if (!ev.getByLabel(InputTag("offlineBeamSpot"), beamSpotHandle)) {
00267         //LogWarning("") << ">>> No beam spot found !!!";
00268         return;
00269       }
00270 
00271 
00272       // Loop to reject/control Z->mumu is done separately
00273       unsigned int nmuonsForZ1 = 0;
00274       unsigned int nmuonsForZ2 = 0;
00275       bool cosmic = false;
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             double dxy = mu.innerTrack()->dxy(beamSpotHandle->position());
00281 
00282             if (fabs(dxy)>1) { cosmic=true; break;} 
00283 
00284             if (pt>ptThrForZ1_) nmuonsForZ1++;
00285             if (pt>ptThrForZ2_) nmuonsForZ2++;
00286 
00287             for (unsigned int j=i+1; j<muonCollectionSize; j++) {
00288                  const Muon& mu2 = muonCollection->at(j);
00289                  if (mu2.isGlobalMuon() && (mu.charge()*mu2.charge()==-1) ){
00290                          const math::XYZTLorentzVector ZRecoGlb (mu.px()+mu2.px(), mu.py()+mu2.py() , mu.pz()+mu2.pz(), mu.p()+mu2.p());
00291                          dimuonmass_before_->Fill(ZRecoGlb.mass());
00292                          if (mu.charge()>0) {
00293                            ptDiffPM_before_->Fill(mu.pt()-mu2.pt());
00294                          }
00295                          else {
00296                            ptDiffPM_before_->Fill(mu2.pt()-mu.pt());
00297                          }
00298                  }
00299             }
00300       }
00301      if(cosmic) return;
00302 
00303       LogTrace("") << "> Z rejection: muons above " << ptThrForZ1_ << " [GeV]: " << nmuonsForZ1;
00304       LogTrace("") << "> Z rejection: muons above " << ptThrForZ2_ << " [GeV]: " << nmuonsForZ2;
00305 
00306       // MET
00307       Handle<View<MET> > metCollection;
00308       if (!ev.getByLabel(metTag_, metCollection)) {
00309         //LogWarning("") << ">>> MET collection does not exist !!!";
00310         return;
00311       }
00312       const MET& met = metCollection->at(0);
00313       double met_et = met.pt();
00314       LogTrace("") << ">>> MET, MET_px, MET_py: " << met_et << ", " << met.px() << ", " << met.py() << " [GeV]";
00315       met_before_->Fill(met_et);
00316 
00317       // Vertices in the event
00318       Handle<View<reco::Vertex> > vertexCollection;
00319            if (!ev.getByLabel(vertexTag_, vertexCollection)) {
00320                  LogError("") << ">>> Vertex collection does not exist !!!";
00321                  return;
00322             }
00323       unsigned int vertexCollectionSize = vertexCollection->size();
00324 
00325       
00326 
00327       int nvvertex = 0;
00328       for (unsigned int i=0; i<vertexCollectionSize; i++) {
00329             const Vertex& vertex = vertexCollection->at(i);
00330             if (vertex.isValid()) nvvertex++;
00331       }
00332 
00333       npvs_before_->Fill(nvvertex);
00334 
00335       bool trigger_fired = false;
00336       Handle<TriggerResults> triggerResults;
00337       if (!ev.getByLabel(trigTag_, triggerResults)) {
00338         //LogWarning("") << ">>> TRIGGER collection does not exist !!!";
00339         return;
00340       }
00341       const edm::TriggerNames & trigNames = ev.triggerNames(*triggerResults);
00342       //  LogWarning("")<<"Loop over triggers";
00343 
00344 
00345       for (unsigned int i=0; i<triggerResults->size(); i++)
00346       {
00347               const std::string trigName = trigNames.triggerName(i);
00348 
00349               bool found=false; 
00350               for(unsigned int index=0; index<trigPathNames_.size() && found==false; index++) {
00351                    size_t trigPath = trigName.find(trigPathNames_[index]); // 0 if found, pos if not
00352                    if (trigPath==0) found=true;
00353               }
00354               if(!found) {continue;}
00355               
00356               bool prescaled=false;    
00357               for (unsigned int ps= 0; ps<  hltConfigProvider_.prescaleSize(); ps++){
00358                   const unsigned int prescaleValue = hltConfigProvider_.prescaleValue(ps, trigName) ;
00359                   if (prescaleValue != 1) prescaled =true;
00360               }
00361             
00362               if( triggerResults->accept(i) && !prescaled){   trigger_fired=true;}
00363                         // LogWarning("")<<"TrigNo: "<<i<<"  "<<found<<"  "<<trigName<<" ---> FIRED";}
00364       }     
00365       trig_before_->Fill(trigger_fired);
00366 
00367       // Jet collection
00368       Handle<View<Jet> > jetCollection;
00369       if (!ev.getByLabel(jetTag_, jetCollection)) {
00370         //LogError("") << ">>> JET collection does not exist !!!";
00371         return;
00372       }
00373       unsigned int jetCollectionSize = jetCollection->size();
00374       int njets = 0; int LEADJET=-1;  double max_pt=0;
00375       for (unsigned int i=0; i<jetCollectionSize; i++) {
00376             const Jet& jet = jetCollection->at(i);
00377                   double minDistance=99999; // This is in order to use PFJets
00378                   for (unsigned int j=0; j<muonCollectionSize; j++) {
00379                         const Muon& mu = muonCollection->at(j);
00380                         double distance = sqrt( (mu.eta()-jet.eta())*(mu.eta()-jet.eta()) +(mu.phi()-jet.phi())*(mu.phi()-jet.phi()) );      
00381                         if (minDistance>distance) minDistance=distance;
00382                   }
00383                   if (minDistance<0.3) continue; // 0.3 is the isolation cone around the muon
00384             if(jet.et()>max_pt) { LEADJET=i; max_pt=jet.et();}
00385             if (jet.et()>eJetMin_) {njets++;}
00386       }
00387 
00388 
00389       LogTrace("") << ">>> Total number of jets: " << jetCollectionSize;
00390       LogTrace("") << ">>> Number of jets above " << eJetMin_ << " [GeV]: " << njets;
00391       njets_before_->Fill(njets);
00392       double lead_jet_pt=-1;
00393       double lead_jet_eta=-100;
00394       if(LEADJET!=-1){
00395       const Jet& leadJet = jetCollection->at(LEADJET);
00396       leadingjet_pt_before_->Fill(leadJet.pt());
00397       leadingjet_eta_before_->Fill(leadJet.eta());
00398       lead_jet_pt=leadJet.pt();
00399       lead_jet_eta=leadJet.eta();
00400       }
00401       //Photon Collection
00402       Handle<View<Photon> > photonCollection;
00403       if(!ev.getByLabel(phoTag_,photonCollection)){
00404       //LogError("")
00405       return;
00406       }
00407       unsigned int ngam=0;
00408       
00409       for (unsigned int i=0; i<photonCollection->size(); i++){
00410         const Photon &ph = photonCollection->at(i);
00411         double photonPt = ph.pt();
00412         if (photonPt> ptThrForPhoton_) {
00413           ngam++;
00414           phEta_->Fill(ph.eta());
00415         }
00416         phPt_->Fill(photonPt); 
00417         }
00418       nph_->Fill(ngam); 
00419       LogTrace("") << " >>> N photons " << ngam << std::endl;
00420 
00421       nmuons_->Fill(muonCollectionSize);
00422 
00424       //Handle<View<Photon> > pfPhotonCollection;
00426       //if(!ev.getByLabel(pfPhoTag_,pfPhotonCollection)){
00428       //return;
00429       //}
00430       //unsigned int npfgam=0;
00431       //
00432       //for (unsigned int i=0; i<pfPhotonCollection->size(); i++){
00433       //        const Photon &ph = pfPhotonCollection->at(i);
00434       //        double photonPt = ph.pt();
00435       //        if (photonPt> ptThrForPhoton_) npfgam++;
00436       //  pfphPt_->Fill(photonPt); 
00437       //        }
00438       //npfph_->Fill(npfgam); 
00439       //LogTrace("") << " >>> N PF photons " << npfgam << std::endl;
00440 
00441       // Start counting
00442       nall++;
00443 
00444       // Histograms per event should be done only once, so keep track of them
00445       bool hlt_hist_done = false;
00446       bool zhlt_hist_done = false;
00447       bool zjets_hist_done = false;
00448       bool zfullsel_hist_done = false;
00449       bool met_hist_done = false;
00450       bool njets_hist_done = false;
00451       bool wfullsel_hist_done = false;
00452 
00453       // Central W->mu nu selection criteria
00454       const int NFLAGS = 11;
00455       bool muon_sel[NFLAGS];
00456       const int NFLAGSZ = 13;
00457       bool zmuon_sel[NFLAGSZ];
00458       bool muon4Z=false;
00459 
00460       double number_of_muons=0;
00461       double number_of_goodMuons=0;
00462 
00463 
00464       for (unsigned int i=0; i<muonCollectionSize; i++) {
00465             for (int j=0; j<NFLAGS; ++j) {
00466                   muon_sel[j] = false;
00467             }
00468 
00469             number_of_muons++;
00470 
00471             const Muon& mu = muonCollection->at(i);
00472             if (!mu.isGlobalMuon()) continue;
00473             if (mu.globalTrack().isNull()) continue;
00474             if (mu.innerTrack().isNull()) continue;
00475 
00476             LogTrace("") << "> Wsel: processing muon number " << i << "...";
00477             reco::TrackRef gm = mu.globalTrack();
00478             reco::TrackRef tk = mu.innerTrack();
00479 
00480             // Pt,eta cuts
00481             double pt = mu.pt();
00482             double eta = mu.eta();
00483             LogTrace("") << "\t... pt, eta: " << pt << " [GeV], " << eta;;
00484             if (pt>ptCut_) muon_sel[0] = true; 
00485             if (fabs(eta)<etaCut_) muon_sel[1] = true; 
00486 
00487             double charge=mu.charge();
00488 
00489             // d0, chi2, nhits quality cuts
00490             double dxy = gm->dxy(beamSpotHandle->position());
00491             double normalizedChi2 = gm->normalizedChi2();
00492             double trackerHits = tk->hitPattern().numberOfValidTrackerHits();
00493             int pixelHits = tk->hitPattern().numberOfValidPixelHits();
00494             int muonHits = gm->hitPattern().numberOfValidMuonHits();
00495             int nMatches = mu.numberOfMatches();
00496 
00497             LogTrace("") << "\t... dxy, normalizedChi2, trackerHits, isTrackerMuon?: " << dxy << " [cm], " << normalizedChi2 << ", " << trackerHits << ", " << mu.isTrackerMuon();
00498             if (fabs(dxy)<dxyCut_) muon_sel[2] = true; 
00499 
00500             bool quality=true;
00501             
00502             if (normalizedChi2>normalizedChi2Cut_) quality =false; 
00503             if (trackerHits<trackerHitsCut_) quality =false;
00504             if (pixelHits<pixelHitsCut_) quality =false;
00505             if (muonHits<muonHitsCut_) quality=false;;
00506             if (!mu.isTrackerMuon()) quality=false;
00507             if (nMatches<nMatchesCut_) quality=false;
00508             muon_sel[3]=quality;
00509             if(quality) number_of_goodMuons++;
00510 
00511             pt_before_->Fill(pt);
00512             eta_before_->Fill(eta);
00513             dxy_before_->Fill(dxy);
00514             muoncharge_before_->Fill(charge);
00515             goodewkmuon_before_->Fill(quality);
00516 
00517             // Charge asymmetry
00518             //if (quality) {
00519             //  if (charge>0) ptPlus_before_->Fill(pt);
00520             //  if (charge<0) ptMinus_before_->Fill(pt);
00521             //}
00522 
00523 
00524             // Isolation cuts
00525             double isovar = mu.isolationR03().sumPt;
00526             if (isCombinedIso_) {
00527                   isovar += mu.isolationR03().emEt;
00528                   isovar += mu.isolationR03().hadEt;
00529             }
00530             if (isRelativeIso_) isovar /= pt;
00531             if (isovar<isoCut03_) muon_sel[4] = true; 
00532 
00533             LogTrace("") << "\t... isolation value" << isovar <<", isolated? " << muon_sel[6];
00534             iso_before_->Fill(isovar);
00535 
00536 
00537             // HLT (not mtched to muon for the time being)
00538             if (trigger_fired) muon_sel[5] = true; 
00539 
00540             // For Z:
00541             if (pt>ptThrForZ1_ && fabs(eta)<etaCut_ && fabs(dxy)<dxyCut_ && quality && trigger_fired && isovar<isoCut03_) { muon4Z = true;}
00542 
00543 
00544             // MET/MT cuts
00545             double w_et = met_et+mu.pt();
00546             double w_px = met.px()+mu.px();
00547             double w_py = met.py()+mu.py();
00548             
00549             double massT = w_et*w_et - w_px*w_px - w_py*w_py;
00550             massT = (massT>0) ? sqrt(massT) : 0;
00551 
00552             LogTrace("") << "\t... W mass, W_et, W_px, W_py: " << massT << ", " << w_et << ", " << w_px << ", " << w_py << " [GeV]";
00553             if (massT>mtMin_ && massT<mtMax_) muon_sel[6] = true; 
00554             mt_before_->Fill(massT);
00555             if (met_et>metMin_ && met_et<metMax_) muon_sel[7] = true; 
00556 
00557             // Acoplanarity cuts
00558             Geom::Phi<double> deltaphi(mu.phi()-atan2(met.py(),met.px()));
00559             double acop = deltaphi.value();
00560             if (acop<0) acop = - acop;
00561             acop = M_PI - acop;
00562             LogTrace("") << "\t... acoplanarity: " << acop;
00563             if (acop<acopCut_) muon_sel[8] = true; 
00564             acop_before_->Fill(acop);
00565 
00566             // Remaining flags (from global event information)
00567             if (nmuonsForZ1<1 || nmuonsForZ2<2) muon_sel[9] = true; 
00568             if (njets<=nJetMax_) muon_sel[10] = true; 
00569 
00570             // Collect necessary flags "per muon"
00571             int flags_passed = 0;
00572             for (int j=0; j<NFLAGS; ++j) {
00573               if (muon_sel[j]) flags_passed += 1;
00574             }
00575 
00576             // Do N-1 histograms now (and only once for global event quantities)
00577             if (flags_passed >= (NFLAGS-1)) {
00578                   if (!muon_sel[0] || flags_passed==NFLAGS) 
00579                         pt_after_->Fill(pt);
00580                   if (!muon_sel[1] || flags_passed==NFLAGS) 
00581                         eta_after_->Fill(eta);
00582                   if (!muon_sel[2] || flags_passed==NFLAGS) 
00583                         dxy_after_->Fill(dxy);
00584                   if (!muon_sel[3] || flags_passed==NFLAGS)
00585                         goodewkmuon_after_->Fill(quality);
00586                   if (!muon_sel[4] || flags_passed==NFLAGS) 
00587                         iso_after_->Fill(isovar);
00588                   if (!muon_sel[5] || flags_passed==NFLAGS) 
00589                         if (!hlt_hist_done) trig_after_->Fill(trigger_fired);
00590                         hlt_hist_done = true;
00591                   if (!muon_sel[6] || flags_passed==NFLAGS) 
00592                         mt_after_->Fill(massT);
00593                   if (!muon_sel[7] || flags_passed==NFLAGS) 
00594                         if (!met_hist_done) met_after_->Fill(met_et);
00595                         met_hist_done = true;
00596                   if (!muon_sel[8] || flags_passed==NFLAGS) 
00597                         acop_after_->Fill(acop);
00598                   // no action here for muon_sel[9]
00599                   if (!muon_sel[10] || flags_passed==NFLAGS) { 
00600                         if (!njets_hist_done) {
00601                                     njets_after_->Fill(njets);
00602                                     leadingjet_pt_after_->Fill(lead_jet_pt);
00603                                     leadingjet_eta_after_->Fill(lead_jet_eta);
00604                         }
00605                         njets_hist_done = true;
00606                   }
00607                   if( flags_passed==NFLAGS ) {
00608                     if (!wfullsel_hist_done){
00609                         npvs_after_->Fill(nvvertex);
00610                         muoncharge_after_->Fill(charge);
00611                         //if (charge>0) ptPlus_afterW_->Fill(pt);
00612                         //if (charge<0) ptMinus_afterW_->Fill(pt);
00613                     }
00614                     wfullsel_hist_done=true;    
00615                   } 
00616             }
00617             
00618 
00619             // The cases in which the event is rejected as a Z are considered independently:
00620             if ( muon4Z &&  !muon_sel[9]){
00621 
00622                    // Plots for 2 muons       
00623                    for (unsigned int j=i+1; j<muonCollectionSize; j++) {
00624 
00625                          for (int ij=0; ij<NFLAGSZ; ++ij) {
00626                                zmuon_sel[ij] = false;
00627                          }
00628 
00629                          for (int ji=0; ji<5; ++ji ) {
00630                            zmuon_sel[ji] = muon_sel[ji];
00631                          }
00632 
00633                          const Muon& mu2 = muonCollection->at(j);
00634                               if (!mu2.isGlobalMuon()) continue;
00635                               if (mu2.charge() * charge != -1 ) continue;
00636                                     reco::TrackRef gm2 = mu2.globalTrack();
00637                                     reco::TrackRef tk2 = mu2.innerTrack();
00638                                     double pt2 = mu2.pt(); if (pt2>ptThrForZ2_) zmuon_sel[5] = true;
00639                                     double eta2=mu2.eta(); if (fabs(eta2)<etaCut_) zmuon_sel[6] = true;
00640                                     double dxy2 = gm2->dxy(beamSpotHandle->position()); if (fabs(dxy2)<dxyCut_) zmuon_sel[7] = true;
00641                                     double normalizedChi22 = gm2->normalizedChi2();
00642                                     double trackerHits2    = tk2->hitPattern().numberOfValidTrackerHits();
00643                                     int    pixelHits2      = tk2->hitPattern().numberOfValidPixelHits();
00644                                     int    muonHits2       = gm2->hitPattern().numberOfValidMuonHits();
00645                                     int    nMatches2       = mu2.numberOfMatches();
00646                                     bool quality2=true;
00647                                     if (normalizedChi22>normalizedChi2Cut_) quality2 = false; 
00648                                     if (trackerHits2<trackerHitsCut_)       quality2 = false;
00649                                     if (pixelHits2<pixelHitsCut_)           quality2 = false;
00650                                     if (muonHits2<muonHitsCut_)             quality2 = false;
00651                                     if (!mu2.isTrackerMuon())               quality2 = false;
00652                                     if (nMatches2<nMatchesCut_)             quality2 = false;
00653                                     zmuon_sel[8]=quality2;
00654                                     double isovar2 = mu2.isolationR03().sumPt; 
00655                                     if (isCombinedIso_) {
00656                                           isovar2 += mu2.isolationR03().emEt;
00657                                           isovar2 += mu2.isolationR03().hadEt;
00658                                     }
00659                                     if (isRelativeIso_) isovar2 /= pt2;
00660                                     if (isovar2<isoCut03_) zmuon_sel[9] = true;
00661                                     if (trigger_fired) zmuon_sel[10] = true; 
00662                                const math::XYZTLorentzVector ZRecoGlb (mu.px()+mu2.px(), mu.py()+mu2.py() , mu.pz()+mu2.pz(), mu.p()+mu2.p());
00663                                if (ZRecoGlb.mass()>dimuonMassMin_ && ZRecoGlb.mass()<dimuonMassMax_) zmuon_sel[11] = true; 
00664 
00665                                         //jet flag
00666                                         if (njets <=nJetMax_) zmuon_sel[12] = true; 
00667 
00668                                // start filling histos: N-1 plots
00669                                int  flags_passed_z = 0;
00670 
00671                                for (int jj=0; jj<NFLAGSZ; ++jj) {
00672                                  if (zmuon_sel[jj]) ++flags_passed_z ;
00673                                }
00674 
00675                                if (flags_passed_z >= (NFLAGSZ-1)) {
00676                                        if (!zmuon_sel[0]  || flags_passed_z==NFLAGSZ) {pt1_afterZ_->Fill(pt); } 
00677                                        if (!zmuon_sel[1]  || flags_passed_z==NFLAGSZ) {eta1_afterZ_->Fill(eta); }  
00678                                        if (!zmuon_sel[2]  || flags_passed_z==NFLAGSZ) {dxy1_afterZ_->Fill(dxy); }  
00679                                        if (!zmuon_sel[3]  || flags_passed_z==NFLAGSZ) {goodewkmuon1_afterZ_->Fill(quality); }  
00680                                        if (!zmuon_sel[4]  || flags_passed_z==NFLAGSZ) {iso1_afterZ_->Fill(isovar); }  
00681                                        if (!zmuon_sel[5]  || flags_passed_z==NFLAGSZ) { pt2_afterZ_->Fill(pt2); }  
00682                                        if (!zmuon_sel[6]  || flags_passed_z==NFLAGSZ) { eta2_afterZ_->Fill(eta2); }  
00683                                        if (!zmuon_sel[7]  || flags_passed_z==NFLAGSZ) {dxy2_afterZ_->Fill(dxy2); }  
00684                                        if (!zmuon_sel[8]  || flags_passed_z==NFLAGSZ) {goodewkmuon2_afterZ_->Fill(quality2); }  
00685                                        if (!zmuon_sel[9]  || flags_passed_z==NFLAGSZ) {iso2_afterZ_->Fill(isovar2); }   
00686                                        if (!zmuon_sel[10] || flags_passed_z==NFLAGSZ) { 
00687                                          if (!zhlt_hist_done) ztrig_afterZ_->Fill(trigger_fired); 
00688                                           zhlt_hist_done = true; 
00689                                        }  
00690                                        if (!zmuon_sel[11] || flags_passed_z==NFLAGSZ) {dimuonmass_afterZ_->Fill(ZRecoGlb.mass()); } 
00691                                        if (!zmuon_sel[12] || flags_passed_z==NFLAGSZ ){
00692                                          if(!zjets_hist_done){
00693                                            njets_afterZ_->Fill(njets);
00694                                            leadingjet_pt_afterZ_->Fill(lead_jet_pt);
00695                                            leadingjet_eta_afterZ_->Fill(lead_jet_eta);
00696                                          }
00697                                          zjets_hist_done=true;
00698                                        }
00699                                        if(flags_passed_z==NFLAGSZ) {
00700                                          met_afterZ_->Fill(met_et);
00701                                          if(!zfullsel_hist_done){
00702                                            npvs_afterZ_->Fill(nvvertex);
00703                                            muoncharge_afterZ_->Fill(charge);
00704                                            if (charge>0) {
00705                                              //ptPlus_afterZ_->Fill(mu.pt());
00706                                              //ptMinus_afterZ_->Fill(mu2.pt());
00707                                              ptDiffPM_afterZ_->Fill(mu.pt()-mu2.pt());
00708                                            }
00709                                            else {
00710                                              //ptPlus_afterZ_->Fill(mu2.pt());
00711                                              //ptMinus_afterZ_->Fill(mu.pt());
00712                                              ptDiffPM_afterZ_->Fill(mu2.pt()-mu.pt());
00713                                            }
00714                                          }
00715                                          zfullsel_hist_done=true;
00716                                        }
00717                                }
00718                    }
00719             }
00720       }
00721 
00722       if (zfullsel_hist_done) {
00723         // here was a Z candidate
00724         n_zselPt1thr_->Fill(nmuonsForZ1);
00725         n_zselPt2thr_->Fill(nmuonsForZ2);
00726       }
00727 
00728       //nmuons_->Fill(number_of_muons);
00729       //nmuons_->Fill(muonCollectionSize);
00730       ngoodmuons_->Fill(number_of_goodMuons);
00731 
00732       return;
00733 
00734 }
00735
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/DQM/Physics/src/EwkMuDQM.cc
