HbbAnalysis/src/HistosMuons.cc

#include <iostream>
#include <fstream>

#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "DataFormats/MuonReco/interface/Muon.h"
#include "DataFormats/MuonReco/interface/MuonIsolation.h"


#include "UserCode/HbbAnalysis/interface/HistosMuons.hh"

namespace HbbAnalysis {//namespace

  void HistosMuons::FillEventHistograms(const edm::Handle<std::vector<pat::Muon> > & aCol){

    if (doGenMatched_) {
      unsigned int nMatched = 0;
      for (std::vector<pat::Muon>::const_iterator iMuon = aCol->begin(); 
           iMuon != aCol->end(); 
           iMuon++) 
        {
          if (MatchesGenMuon(*iMuon)) nMatched++;
        }
      p_nMuons->Fill(nMatched);
    }
    else p_nMuons->Fill(aCol->size());

  }

  void HistosMuons::FillHistograms(const pat::Muon & aMuon, const edm::Handle<std::vector<reco::Vertex> > & aRecoVertices, bool isLead){//FillHistograms

    const reco::Muon* recoMuon = dynamic_cast<const reco::Muon*>(aMuon.originalObject());
    bool lIsGenMatched = !doGenMatched_ || (doGenMatched_ && MatchesGenMuon(aMuon));
    if (lIsGenMatched) {
      if (isLead) {
        FillBaseHistograms(aMuon.pt(),aMuon.eta(),aMuon.phi(),aMuon.charge());
        p_caloCompat->Fill(recoMuon->caloCompatibility());
        p_segCompat->Fill(muon::segmentCompatibility(*recoMuon));
        p_nChambers->Fill(recoMuon->numberOfChambers());
        p_nMatchesLoose->Fill(recoMuon->numberOfMatches(reco::Muon::NoArbitration));
        p_nMatchesMedium->Fill(recoMuon->numberOfMatches(reco::Muon::SegmentArbitration));
        p_nMatchesTight->Fill(recoMuon->numberOfMatches(reco::Muon::SegmentAndTrackArbitration));
        p_type->Fill(recoMuon->type());
        
        if (recoMuon->isGlobalMuon()) {
          p_muonType->Fill(1);
          p_muonTypevsPt->Fill(aMuon.pt(),1);
          p_muonTypevsEta->Fill(aMuon.eta(),1);
        }
        else if (recoMuon->isTrackerMuon()) {
          p_muonType->Fill(2);
          p_muonTypevsPt->Fill(aMuon.pt(),2);
          p_muonTypevsEta->Fill(aMuon.eta(),2);
        }
        else if (recoMuon->isStandAloneMuon()) {
          p_muonType->Fill(3);
          p_muonTypevsPt->Fill(aMuon.pt(),3);
          p_muonTypevsEta->Fill(aMuon.eta(),3);
        }
        else if (recoMuon->isCaloMuon()) {
          p_muonType->Fill(4);
          p_muonTypevsPt->Fill(aMuon.pt(),4);
          p_muonTypevsEta->Fill(aMuon.eta(),4);
        }
        else {
          p_muonType->Fill(0);
          p_muonTypevsPt->Fill(aMuon.pt(),0);
          p_muonTypevsEta->Fill(aMuon.eta(),0);
        }

        if ( aMuon.innerTrack().isAvailable() && aMuon.innerTrack().isNonnull() ) {
          if ( aRecoVertices->size() >= 1 ) {
            const reco::Vertex& thePrimaryEventVertex = (*aRecoVertices->begin());
            p_trkIPd0->Fill(-aMuon.innerTrack()->dxy(thePrimaryEventVertex.position()));
            p_trkIPdz->Fill(aMuon.innerTrack()->dz(thePrimaryEventVertex.position()));
          }
          else {
            p_trkIPd0->Fill(0);
            p_trkIPdz->Fill(0);
          }
          p_trknHits->Fill(aMuon.innerTrack()->numberOfValidHits());
        }
        else {
          p_trkIPd0->Fill(0);
          p_trkIPdz->Fill(0);
          p_trknHits->Fill(0);
        }

        for (unsigned int id(0); id<idEff_.numberOfBins(); id++){
          idEff_.incrementTotal(id);
          for (unsigned ieta(0); ieta<idEffEta_[id].numberOfBins(); ieta++){
            double lCutMin = idEffEta_[id].xMin()+idEffEta_[id].stepSize()*ieta;
            double lCutMax = idEffEta_[id].xMin()+idEffEta_[id].stepSize()*(ieta+1);
            
            if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
              idEffEta_[id].incrementTotal(ieta);
            }
          }
        }

        for (unsigned ieta(0); ieta<idEffEta_[0].numberOfBins(); ieta++){
          double lCutMin = idEffEta_[0].xMin()+idEffEta_[0].stepSize()*ieta;
          double lCutMax = idEffEta_[0].xMin()+idEffEta_[0].stepSize()*(ieta+1);
            

          if (muon::isGoodMuon(*recoMuon,muon::AllGlobalMuons)){
            if (ieta == 0) {
              p_muonID->Fill(1);
              p_muonIDvsPt->Fill(aMuon.pt(),1);
              p_muonIDvsEta->Fill(aMuon.eta(),1);
              idEff_.incrementPass(1);
            }
            if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
              idEffEta_[1].incrementPass(ieta);
            }
          }
          if (muon::isGoodMuon(*recoMuon,muon::AllStandAloneMuons)){
            if (ieta == 0) {
              p_muonID->Fill(2);
              p_muonIDvsPt->Fill(aMuon.pt(),2);
              p_muonIDvsEta->Fill(aMuon.eta(),2);
              idEff_.incrementPass(2);
            }
            if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
              idEffEta_[2].incrementPass(ieta);
            }
          }
          if (muon::isGoodMuon(*recoMuon,muon::AllTrackerMuons)){
            if (ieta == 0) {
              p_muonID->Fill(3);
              p_muonIDvsPt->Fill(aMuon.pt(),3);
              p_muonIDvsEta->Fill(aMuon.eta(),3);
              idEff_.incrementPass(3);
            }
            if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
              idEffEta_[3].incrementPass(ieta);
            }
          }
          if (muon::isGoodMuon(*recoMuon,muon::TrackerMuonArbitrated)){
            if (ieta == 0) {
              p_muonID->Fill(4);
              p_muonIDvsPt->Fill(aMuon.pt(),4);
              p_muonIDvsEta->Fill(aMuon.eta(),4);
              idEff_.incrementPass(4);
            }
            if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
              idEffEta_[4].incrementPass(ieta);
            }
          }
          if (muon::isGoodMuon(*recoMuon,muon::AllArbitrated)){
            if (ieta == 0) {
              p_muonID->Fill(5);
              p_muonIDvsPt->Fill(aMuon.pt(),5);
              p_muonIDvsEta->Fill(aMuon.eta(),5);
              idEff_.incrementPass(5);
            }
            if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
              idEffEta_[5].incrementPass(ieta);
            }
          }
          if (muon::isGoodMuon(*recoMuon,muon::GlobalMuonPromptTight)){
            if (ieta == 0) {
              p_muonID->Fill(6);
              p_muonIDvsPt->Fill(aMuon.pt(),6);
              p_muonIDvsEta->Fill(aMuon.eta(),6);
              idEff_.incrementPass(6);
            }
            if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
              idEffEta_[6].incrementPass(ieta);
            }
          }
          if (muon::isGoodMuon(*recoMuon,muon::TMLastStationLoose)){
            if (ieta == 0) {
              p_muonID->Fill(7);
              p_muonIDvsPt->Fill(aMuon.pt(),7);
              p_muonIDvsEta->Fill(aMuon.eta(),7);
              idEff_.incrementPass(7);
            }
            if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
              idEffEta_[7].incrementPass(ieta);
            }
          }
          if (muon::isGoodMuon(*recoMuon,muon::TMLastStationTight)){
            if (ieta == 0) {
              p_muonID->Fill(8);
              p_muonIDvsPt->Fill(aMuon.pt(),8);
              p_muonIDvsEta->Fill(aMuon.eta(),8);
              idEff_.incrementPass(8);
            }
            if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
              idEffEta_[8].incrementPass(ieta);
            }
          }
          if (muon::isGoodMuon(*recoMuon,muon::TM2DCompatibilityLoose)){
            if (ieta == 0) {
              p_muonID->Fill(9);
              p_muonIDvsPt->Fill(aMuon.pt(),9);
              p_muonIDvsEta->Fill(aMuon.eta(),9);
              idEff_.incrementPass(9);
            }
            if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
              idEffEta_[9].incrementPass(ieta);
            }
          }
          if (muon::isGoodMuon(*recoMuon,muon::TM2DCompatibilityTight)){
            if (ieta == 0) {
              p_muonID->Fill(10);
              p_muonIDvsPt->Fill(aMuon.pt(),10);
              p_muonIDvsEta->Fill(aMuon.eta(),10);
              idEff_.incrementPass(10);
            }
            if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
              idEffEta_[10].incrementPass(ieta);
            }
          }
          if (muon::isGoodMuon(*recoMuon,muon::TMOneStationLoose)){
            if (ieta == 0) {
              p_muonID->Fill(11);
              p_muonIDvsPt->Fill(aMuon.pt(),11);
              p_muonIDvsEta->Fill(aMuon.eta(),11);
              idEff_.incrementPass(11);
            }
            if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
              idEffEta_[11].incrementPass(ieta);
            }
          }
          if (muon::isGoodMuon(*recoMuon,muon::TMOneStationTight)){
            if (ieta == 0) {
              p_muonID->Fill(12);
              p_muonIDvsPt->Fill(aMuon.pt(),12);
              p_muonIDvsEta->Fill(aMuon.eta(),12);
              idEff_.incrementPass(12);
            }
            if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
              idEffEta_[12].incrementPass(ieta);
            }
          }
          if (muon::isGoodMuon(*recoMuon,muon::TMLastStationOptimizedLowPtLoose)){
            if (ieta == 0) {
              p_muonID->Fill(13);
              p_muonIDvsPt->Fill(aMuon.pt(),13);
              p_muonIDvsEta->Fill(aMuon.eta(),13);
              idEff_.incrementPass(13);
            }
            if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
              idEffEta_[13].incrementPass(ieta);
            }
          }
          if (muon::isGoodMuon(*recoMuon,muon::TMLastStationOptimizedLowPtTight)){
            if (ieta == 0) {
              p_muonID->Fill(14);
              p_muonIDvsPt->Fill(aMuon.pt(),14);
              p_muonIDvsEta->Fill(aMuon.eta(),14);
              idEff_.incrementPass(14);
            }
            if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
              idEffEta_[14].incrementPass(ieta);
            }
          }

        }

        p_caloCompatvsPt->Fill(aMuon.pt(),recoMuon->caloCompatibility());
        p_segCompatvsPt->Fill(aMuon.pt(),muon::segmentCompatibility(*recoMuon));
        p_nChambersvsPt->Fill(aMuon.pt(),recoMuon->numberOfChambers());
        p_nMatchesvsPt->Fill(aMuon.pt(),recoMuon->numberOfMatches());

        p_caloCompatvsEta->Fill(aMuon.eta(),recoMuon->caloCompatibility());
        p_segCompatvsEta->Fill(aMuon.eta(),muon::segmentCompatibility(*recoMuon));
        p_nChambersvsEta->Fill(aMuon.eta(),recoMuon->numberOfChambers());
        p_nMatchesvsEta->Fill(aMuon.eta(),recoMuon->numberOfMatches());

      }

      p_isoR03_emEt->Fill(recoMuon->isolationR03().emEt);
      p_isoR05_emEt->Fill(recoMuon->isolationR05().emEt);
      p_isoR03_hadEt->Fill(recoMuon->isolationR03().hadEt);
      p_isoR05_hadEt->Fill(recoMuon->isolationR05().hadEt);
      p_isoR03_nTracks->Fill(recoMuon->isolationR03().nTracks);
      p_isoR05_nTracks->Fill(recoMuon->isolationR05().nTracks);
      p_isoR03_nJets->Fill(recoMuon->isolationR03().nJets);
      p_isoR05_nJets->Fill(recoMuon->isolationR05().nJets);

      double lIsoR03Var[4] = {
        recoMuon->isolationR03().sumPt,
        recoMuon->isolationR03().sumPt/aMuon.pt(),
        recoMuon->isolationR03().sumPt+recoMuon->isolationR03().emEt,
        (recoMuon->isolationR03().sumPt+recoMuon->isolationR03().emEt)/aMuon.pt() 
      };
      double lIsoR05Var[4] = {
        recoMuon->isolationR05().sumPt,
        recoMuon->isolationR05().sumPt/aMuon.pt(),
        recoMuon->isolationR05().sumPt+recoMuon->isolationR05().emEt,
        (recoMuon->isolationR05().sumPt+recoMuon->isolationR05().emEt)/aMuon.pt() 
      };

      double lIsoR03Cut[4] = {3.,0.09,6.,0.18};

      for (unsigned int i(0); i<getNVar(); i++){//loop on iso variables
        p_isoR03[i]->Fill(lIsoR03Var[i]);
        p_isoR05[i]->Fill(lIsoR05Var[i]);
 
        //R03
        for (unsigned int lBin(0); lBin<isoR03Eff_[i].numberOfBins(); lBin++){
          double lCut = isoR03Eff_[i].xMin()+isoR03Eff_[i].stepSize()*lBin;

          //std::cout << "******** var R03 #" << i << ", bin " << lBin << " isolation = " << lIsoR03Var[i] << ", cut = " << lCut << std::endl;

          isoR03Eff_[i].incrementTotal(lBin);
          if (lIsoR03Var[i] <= lCut){
            isoR03Eff_[i].incrementPass(lBin);
          }
        }
  
        //R05
        for (unsigned int lBin(0); lBin<isoR05Eff_[i].numberOfBins(); lBin++){
          double lCut = isoR05Eff_[i].xMin()+isoR05Eff_[i].stepSize()*lBin;

          //std::cout << "******** var R05 #" << i << ", bin " << lBin << " isolation = " << lIsoR03Var[i] << ", cut = " << lCut << std::endl;

          isoR05Eff_[i].incrementTotal(lBin);
          if (lIsoR05Var[i] <= lCut){
            isoR05Eff_[i].incrementPass(lBin);
          }
        }
  
        //eff vs pT
        for (unsigned int lBin(0); lBin<isoEffPt_[i].numberOfBins(); lBin++){
          double lCutMin = isoEffPt_[i].xMin()+isoEffPt_[i].stepSize()*lBin;
          double lCutMax = isoEffPt_[i].xMin()+isoEffPt_[i].stepSize()*(lBin+1);

          if (aMuon.pt() >= lCutMin && aMuon.pt() < lCutMax){
            isoEffPt_[i].incrementTotal(lBin);

            if (lIsoR03Var[i] <= lIsoR03Cut[i]){
              isoEffPt_[i].incrementPass(lBin);
            }
          }
        }

        //eff vs eta
        for (unsigned int lBin(0); lBin<isoEffEta_[i].numberOfBins(); lBin++){
          double lCutMin = isoEffEta_[i].xMin()+isoEffEta_[i].stepSize()*lBin;
          double lCutMax = isoEffEta_[i].xMin()+isoEffEta_[i].stepSize()*(lBin+1);

          if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
            isoEffEta_[i].incrementTotal(lBin);

            if (lIsoR03Var[i] <= lIsoR03Cut[i]){
              isoEffEta_[i].incrementPass(lBin);
            }
          }
        }

      }//loop on iso variables
    }//lIsGenMatched
         
  }//FillHistograms

  bool HistosMuons::MatchesGenMuon(const pat::Muon& aPatMuon)
  {
    
    bool isGenMatched = false;

    //if (aPatMuon.genParticlesSize() == 0) {
    //  edm::LogWarning("MatchesGenMuon") << "No ref to gen particles for pat::Muon ! ";
    //  return false;
    // }

    const std::vector<reco::GenParticleRef> & lVec = aPatMuon.genParticleRefs();
    //for ( std::vector<reco::GenParticleRef>::const_iterator it = aPatMuon.genParticleRefs().begin(); 
    //    it != aPatMuon.genParticleRefs().end(); ++it ) {
    for ( std::vector<reco::GenParticleRef>::const_iterator it = lVec.begin(); 
          it != lVec.end(); ++it ) {
      if ( it->ref().isNonnull() && it->ref().isValid() ) {
        const reco::GenParticleRef & genParticle = (*it);
        //assert (genParticle->isNonnull() && genParticle->isValid());
        if ( genParticle->pdgId() == -13 || genParticle->pdgId() == +13 ) isGenMatched = true;
      } else {
        edm::LogWarning("MatchesGenMuon") << " edm::Ref of genParticle associated to pat::Muon is invalid !!";
      }
    }
    return isGenMatched;
  }

}//namespace


Revision:	1.7
Committed:	Tue Feb 9 14:52:25 2010 UTC (15 years, 3 months ago) by amagnan
Content type:	text/plain
Branch:	MAIN
CVS Tags:	v00-05-00, HbbAnaFor35X, v00-04-02, v00-04-01, v00-04-00
Changes since 1.6:	+17 -17 lines
Log Message:	export code to CMSSW_3_4_X
#	User	Rev	Content
1	amagnan	1.1	#include <iostream>
2			#include <fstream>
3
4			#include "FWCore/MessageLogger/interface/MessageLogger.h"
5
6			#include "DataFormats/MuonReco/interface/Muon.h"
7			#include "DataFormats/MuonReco/interface/MuonIsolation.h"
8
9
10			#include "UserCode/HbbAnalysis/interface/HistosMuons.hh"
11
12			namespace HbbAnalysis {//namespace
13
14	amagnan	1.3	void HistosMuons::FillEventHistograms(const edm::Handle<std::vector<pat::Muon> > & aCol){
15	amagnan	1.1
16			if (doGenMatched_) {
17	amagnan	1.2	unsigned int nMatched = 0;
18			for (std::vector<pat::Muon>::const_iterator iMuon = aCol->begin();
19			iMuon != aCol->end();
20	amagnan	1.1	iMuon++)
21			{
22	amagnan	1.2	if (MatchesGenMuon(*iMuon)) nMatched++;
23	amagnan	1.1	}
24	amagnan	1.2	p_nMuons->Fill(nMatched);
25	amagnan	1.1	}
26	amagnan	1.2	else p_nMuons->Fill(aCol->size());
27	amagnan	1.1
28			}
29
30	amagnan	1.6	void HistosMuons::FillHistograms(const pat::Muon & aMuon, const edm::Handle<std::vector<reco::Vertex> > & aRecoVertices, bool isLead){//FillHistograms
31	amagnan	1.1
32			const reco::Muon* recoMuon = dynamic_cast<const reco::Muon*>(aMuon.originalObject());
33			bool lIsGenMatched = !doGenMatched_ \|\| (doGenMatched_ && MatchesGenMuon(aMuon));
34			if (lIsGenMatched) {
35			if (isLead) {
36	amagnan	1.4	FillBaseHistograms(aMuon.pt(),aMuon.eta(),aMuon.phi(),aMuon.charge());
37	amagnan	1.1	p_caloCompat->Fill(recoMuon->caloCompatibility());
38	amagnan	1.7	p_segCompat->Fill(muon::segmentCompatibility(*recoMuon));
39	amagnan	1.1	p_nChambers->Fill(recoMuon->numberOfChambers());
40			p_nMatchesLoose->Fill(recoMuon->numberOfMatches(reco::Muon::NoArbitration));
41			p_nMatchesMedium->Fill(recoMuon->numberOfMatches(reco::Muon::SegmentArbitration));
42			p_nMatchesTight->Fill(recoMuon->numberOfMatches(reco::Muon::SegmentAndTrackArbitration));
43			p_type->Fill(recoMuon->type());
44
45			if (recoMuon->isGlobalMuon()) {
46			p_muonType->Fill(1);
47			p_muonTypevsPt->Fill(aMuon.pt(),1);
48			p_muonTypevsEta->Fill(aMuon.eta(),1);
49			}
50			else if (recoMuon->isTrackerMuon()) {
51			p_muonType->Fill(2);
52			p_muonTypevsPt->Fill(aMuon.pt(),2);
53			p_muonTypevsEta->Fill(aMuon.eta(),2);
54			}
55			else if (recoMuon->isStandAloneMuon()) {
56			p_muonType->Fill(3);
57			p_muonTypevsPt->Fill(aMuon.pt(),3);
58			p_muonTypevsEta->Fill(aMuon.eta(),3);
59			}
60			else if (recoMuon->isCaloMuon()) {
61			p_muonType->Fill(4);
62			p_muonTypevsPt->Fill(aMuon.pt(),4);
63			p_muonTypevsEta->Fill(aMuon.eta(),4);
64			}
65			else {
66			p_muonType->Fill(0);
67			p_muonTypevsPt->Fill(aMuon.pt(),0);
68			p_muonTypevsEta->Fill(aMuon.eta(),0);
69			}
70
71	amagnan	1.6	if ( aMuon.innerTrack().isAvailable() && aMuon.innerTrack().isNonnull() ) {
72			if ( aRecoVertices->size() >= 1 ) {
73			const reco::Vertex& thePrimaryEventVertex = (*aRecoVertices->begin());
74			p_trkIPd0->Fill(-aMuon.innerTrack()->dxy(thePrimaryEventVertex.position()));
75			p_trkIPdz->Fill(aMuon.innerTrack()->dz(thePrimaryEventVertex.position()));
76			}
77			else {
78			p_trkIPd0->Fill(0);
79			p_trkIPdz->Fill(0);
80			}
81			p_trknHits->Fill(aMuon.innerTrack()->numberOfValidHits());
82			}
83			else {
84			p_trkIPd0->Fill(0);
85			p_trkIPdz->Fill(0);
86			p_trknHits->Fill(0);
87			}
88
89	amagnan	1.1	for (unsigned int id(0); id<idEff_.numberOfBins(); id++){
90			idEff_.incrementTotal(id);
91			for (unsigned ieta(0); ieta<idEffEta_[id].numberOfBins(); ieta++){
92			double lCutMin = idEffEta_[id].xMin()+idEffEta_[id].stepSize()*ieta;
93			double lCutMax = idEffEta_[id].xMin()+idEffEta_[id].stepSize()*(ieta+1);
94
95			if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
96			idEffEta_[id].incrementTotal(ieta);
97			}
98			}
99			}
100
101			for (unsigned ieta(0); ieta<idEffEta_[0].numberOfBins(); ieta++){
102			double lCutMin = idEffEta_[0].xMin()+idEffEta_[0].stepSize()*ieta;
103			double lCutMax = idEffEta_[0].xMin()+idEffEta_[0].stepSize()*(ieta+1);
104
105
106	amagnan	1.7	if (muon::isGoodMuon(*recoMuon,muon::AllGlobalMuons)){
107	amagnan	1.1	if (ieta == 0) {
108			p_muonID->Fill(1);
109			p_muonIDvsPt->Fill(aMuon.pt(),1);
110			p_muonIDvsEta->Fill(aMuon.eta(),1);
111			idEff_.incrementPass(1);
112			}
113			if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
114			idEffEta_[1].incrementPass(ieta);
115			}
116			}
117	amagnan	1.7	if (muon::isGoodMuon(*recoMuon,muon::AllStandAloneMuons)){
118	amagnan	1.1	if (ieta == 0) {
119			p_muonID->Fill(2);
120			p_muonIDvsPt->Fill(aMuon.pt(),2);
121			p_muonIDvsEta->Fill(aMuon.eta(),2);
122			idEff_.incrementPass(2);
123			}
124			if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
125			idEffEta_[2].incrementPass(ieta);
126			}
127			}
128	amagnan	1.7	if (muon::isGoodMuon(*recoMuon,muon::AllTrackerMuons)){
129	amagnan	1.1	if (ieta == 0) {
130			p_muonID->Fill(3);
131			p_muonIDvsPt->Fill(aMuon.pt(),3);
132			p_muonIDvsEta->Fill(aMuon.eta(),3);
133			idEff_.incrementPass(3);
134			}
135			if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
136			idEffEta_[3].incrementPass(ieta);
137			}
138			}
139	amagnan	1.7	if (muon::isGoodMuon(*recoMuon,muon::TrackerMuonArbitrated)){
140	amagnan	1.1	if (ieta == 0) {
141			p_muonID->Fill(4);
142			p_muonIDvsPt->Fill(aMuon.pt(),4);
143			p_muonIDvsEta->Fill(aMuon.eta(),4);
144			idEff_.incrementPass(4);
145			}
146			if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
147			idEffEta_[4].incrementPass(ieta);
148			}
149			}
150	amagnan	1.7	if (muon::isGoodMuon(*recoMuon,muon::AllArbitrated)){
151	amagnan	1.1	if (ieta == 0) {
152			p_muonID->Fill(5);
153			p_muonIDvsPt->Fill(aMuon.pt(),5);
154			p_muonIDvsEta->Fill(aMuon.eta(),5);
155			idEff_.incrementPass(5);
156			}
157			if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
158			idEffEta_[5].incrementPass(ieta);
159			}
160			}
161	amagnan	1.7	if (muon::isGoodMuon(*recoMuon,muon::GlobalMuonPromptTight)){
162	amagnan	1.1	if (ieta == 0) {
163			p_muonID->Fill(6);
164			p_muonIDvsPt->Fill(aMuon.pt(),6);
165			p_muonIDvsEta->Fill(aMuon.eta(),6);
166			idEff_.incrementPass(6);
167			}
168			if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
169			idEffEta_[6].incrementPass(ieta);
170			}
171			}
172	amagnan	1.7	if (muon::isGoodMuon(*recoMuon,muon::TMLastStationLoose)){
173	amagnan	1.1	if (ieta == 0) {
174			p_muonID->Fill(7);
175			p_muonIDvsPt->Fill(aMuon.pt(),7);
176			p_muonIDvsEta->Fill(aMuon.eta(),7);
177			idEff_.incrementPass(7);
178			}
179			if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
180			idEffEta_[7].incrementPass(ieta);
181			}
182			}
183	amagnan	1.7	if (muon::isGoodMuon(*recoMuon,muon::TMLastStationTight)){
184	amagnan	1.1	if (ieta == 0) {
185			p_muonID->Fill(8);
186			p_muonIDvsPt->Fill(aMuon.pt(),8);
187			p_muonIDvsEta->Fill(aMuon.eta(),8);
188			idEff_.incrementPass(8);
189			}
190			if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
191			idEffEta_[8].incrementPass(ieta);
192			}
193			}
194	amagnan	1.7	if (muon::isGoodMuon(*recoMuon,muon::TM2DCompatibilityLoose)){
195	amagnan	1.1	if (ieta == 0) {
196			p_muonID->Fill(9);
197			p_muonIDvsPt->Fill(aMuon.pt(),9);
198			p_muonIDvsEta->Fill(aMuon.eta(),9);
199			idEff_.incrementPass(9);
200			}
201			if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
202			idEffEta_[9].incrementPass(ieta);
203			}
204			}
205	amagnan	1.7	if (muon::isGoodMuon(*recoMuon,muon::TM2DCompatibilityTight)){
206	amagnan	1.1	if (ieta == 0) {
207			p_muonID->Fill(10);
208			p_muonIDvsPt->Fill(aMuon.pt(),10);
209			p_muonIDvsEta->Fill(aMuon.eta(),10);
210			idEff_.incrementPass(10);
211			}
212			if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
213			idEffEta_[10].incrementPass(ieta);
214			}
215			}
216	amagnan	1.7	if (muon::isGoodMuon(*recoMuon,muon::TMOneStationLoose)){
217	amagnan	1.1	if (ieta == 0) {
218			p_muonID->Fill(11);
219			p_muonIDvsPt->Fill(aMuon.pt(),11);
220			p_muonIDvsEta->Fill(aMuon.eta(),11);
221			idEff_.incrementPass(11);
222			}
223			if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
224			idEffEta_[11].incrementPass(ieta);
225			}
226			}
227	amagnan	1.7	if (muon::isGoodMuon(*recoMuon,muon::TMOneStationTight)){
228	amagnan	1.1	if (ieta == 0) {
229			p_muonID->Fill(12);
230			p_muonIDvsPt->Fill(aMuon.pt(),12);
231			p_muonIDvsEta->Fill(aMuon.eta(),12);
232			idEff_.incrementPass(12);
233			}
234			if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
235			idEffEta_[12].incrementPass(ieta);
236			}
237			}
238	amagnan	1.7	if (muon::isGoodMuon(*recoMuon,muon::TMLastStationOptimizedLowPtLoose)){
239	amagnan	1.1	if (ieta == 0) {
240			p_muonID->Fill(13);
241			p_muonIDvsPt->Fill(aMuon.pt(),13);
242			p_muonIDvsEta->Fill(aMuon.eta(),13);
243			idEff_.incrementPass(13);
244			}
245			if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
246			idEffEta_[13].incrementPass(ieta);
247			}
248			}
249	amagnan	1.7	if (muon::isGoodMuon(*recoMuon,muon::TMLastStationOptimizedLowPtTight)){
250	amagnan	1.1	if (ieta == 0) {
251			p_muonID->Fill(14);
252			p_muonIDvsPt->Fill(aMuon.pt(),14);
253			p_muonIDvsEta->Fill(aMuon.eta(),14);
254			idEff_.incrementPass(14);
255			}
256			if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
257			idEffEta_[14].incrementPass(ieta);
258			}
259			}
260
261			}
262
263			p_caloCompatvsPt->Fill(aMuon.pt(),recoMuon->caloCompatibility());
264	amagnan	1.7	p_segCompatvsPt->Fill(aMuon.pt(),muon::segmentCompatibility(*recoMuon));
265	amagnan	1.1	p_nChambersvsPt->Fill(aMuon.pt(),recoMuon->numberOfChambers());
266			p_nMatchesvsPt->Fill(aMuon.pt(),recoMuon->numberOfMatches());
267
268			p_caloCompatvsEta->Fill(aMuon.eta(),recoMuon->caloCompatibility());
269	amagnan	1.7	p_segCompatvsEta->Fill(aMuon.eta(),muon::segmentCompatibility(*recoMuon));
270	amagnan	1.1	p_nChambersvsEta->Fill(aMuon.eta(),recoMuon->numberOfChambers());
271			p_nMatchesvsEta->Fill(aMuon.eta(),recoMuon->numberOfMatches());
272
273			}
274
275			p_isoR03_emEt->Fill(recoMuon->isolationR03().emEt);
276			p_isoR05_emEt->Fill(recoMuon->isolationR05().emEt);
277			p_isoR03_hadEt->Fill(recoMuon->isolationR03().hadEt);
278			p_isoR05_hadEt->Fill(recoMuon->isolationR05().hadEt);
279			p_isoR03_nTracks->Fill(recoMuon->isolationR03().nTracks);
280			p_isoR05_nTracks->Fill(recoMuon->isolationR05().nTracks);
281			p_isoR03_nJets->Fill(recoMuon->isolationR03().nJets);
282			p_isoR05_nJets->Fill(recoMuon->isolationR05().nJets);
283
284			double lIsoR03Var[4] = {
285			recoMuon->isolationR03().sumPt,
286			recoMuon->isolationR03().sumPt/aMuon.pt(),
287			recoMuon->isolationR03().sumPt+recoMuon->isolationR03().emEt,
288			(recoMuon->isolationR03().sumPt+recoMuon->isolationR03().emEt)/aMuon.pt()
289			};
290			double lIsoR05Var[4] = {
291			recoMuon->isolationR05().sumPt,
292			recoMuon->isolationR05().sumPt/aMuon.pt(),
293			recoMuon->isolationR05().sumPt+recoMuon->isolationR05().emEt,
294			(recoMuon->isolationR05().sumPt+recoMuon->isolationR05().emEt)/aMuon.pt()
295			};
296
297			double lIsoR03Cut[4] = {3.,0.09,6.,0.18};
298
299			for (unsigned int i(0); i<getNVar(); i++){//loop on iso variables
300			p_isoR03[i]->Fill(lIsoR03Var[i]);
301			p_isoR05[i]->Fill(lIsoR05Var[i]);
302
303			//R03
304			for (unsigned int lBin(0); lBin<isoR03Eff_[i].numberOfBins(); lBin++){
305			double lCut = isoR03Eff_[i].xMin()+isoR03Eff_[i].stepSize()*lBin;
306
307			//std::cout << "******** var R03 #" << i << ", bin " << lBin << " isolation = " << lIsoR03Var[i] << ", cut = " << lCut << std::endl;
308
309			isoR03Eff_[i].incrementTotal(lBin);
310			if (lIsoR03Var[i] <= lCut){
311			isoR03Eff_[i].incrementPass(lBin);
312			}
313			}
314
315			//R05
316			for (unsigned int lBin(0); lBin<isoR05Eff_[i].numberOfBins(); lBin++){
317			double lCut = isoR05Eff_[i].xMin()+isoR05Eff_[i].stepSize()*lBin;
318
319			//std::cout << "******** var R05 #" << i << ", bin " << lBin << " isolation = " << lIsoR03Var[i] << ", cut = " << lCut << std::endl;
320
321			isoR05Eff_[i].incrementTotal(lBin);
322			if (lIsoR05Var[i] <= lCut){
323			isoR05Eff_[i].incrementPass(lBin);
324			}
325			}
326
327			//eff vs pT
328			for (unsigned int lBin(0); lBin<isoEffPt_[i].numberOfBins(); lBin++){
329			double lCutMin = isoEffPt_[i].xMin()+isoEffPt_[i].stepSize()*lBin;
330			double lCutMax = isoEffPt_[i].xMin()+isoEffPt_[i].stepSize()*(lBin+1);
331
332			if (aMuon.pt() >= lCutMin && aMuon.pt() < lCutMax){
333			isoEffPt_[i].incrementTotal(lBin);
334
335			if (lIsoR03Var[i] <= lIsoR03Cut[i]){
336			isoEffPt_[i].incrementPass(lBin);
337			}
338			}
339			}
340
341			//eff vs eta
342			for (unsigned int lBin(0); lBin<isoEffEta_[i].numberOfBins(); lBin++){
343			double lCutMin = isoEffEta_[i].xMin()+isoEffEta_[i].stepSize()*lBin;
344			double lCutMax = isoEffEta_[i].xMin()+isoEffEta_[i].stepSize()*(lBin+1);
345
346			if (aMuon.eta() >= lCutMin && aMuon.eta() < lCutMax){
347			isoEffEta_[i].incrementTotal(lBin);
348
349			if (lIsoR03Var[i] <= lIsoR03Cut[i]){
350			isoEffEta_[i].incrementPass(lBin);
351			}
352			}
353			}
354
355			}//loop on iso variables
356			}//lIsGenMatched
357
358			}//FillHistograms
359
360			bool HistosMuons::MatchesGenMuon(const pat::Muon& aPatMuon)
361			{
362
363	amagnan	1.2	bool isGenMatched = false;
364	amagnan	1.1
365			//if (aPatMuon.genParticlesSize() == 0) {
366			// edm::LogWarning("MatchesGenMuon") << "No ref to gen particles for pat::Muon ! ";
367			// return false;
368			// }
369
370			const std::vector<reco::GenParticleRef> & lVec = aPatMuon.genParticleRefs();
371			//for ( std::vector<reco::GenParticleRef>::const_iterator it = aPatMuon.genParticleRefs().begin();
372			// it != aPatMuon.genParticleRefs().end(); ++it ) {
373			for ( std::vector<reco::GenParticleRef>::const_iterator it = lVec.begin();
374			it != lVec.end(); ++it ) {
375			if ( it->ref().isNonnull() && it->ref().isValid() ) {
376			const reco::GenParticleRef & genParticle = (*it);
377			//assert (genParticle->isNonnull() && genParticle->isValid());
378	amagnan	1.2	if ( genParticle->pdgId() == -13 \|\| genParticle->pdgId() == +13 ) isGenMatched = true;
379	amagnan	1.1	} else {
380			edm::LogWarning("MatchesGenMuon") << " edm::Ref of genParticle associated to pat::Muon is invalid !!";
381			}
382			}
383	amagnan	1.2	return isGenMatched;
384	amagnan	1.1	}
385
386			}//namespace
387
388