HbbAnalysis/interface/JetFlavour.hh

#ifndef HbbAnalysis_JetFlavour_hh
#define HbbAnalysis_JetFlavour_hh

#include <vector>
#include <string>
#include <cassert>

#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
#include "DataFormats/PatCandidates/interface/Jet.h"

#include "FWCore/ServiceRegistry/interface/Service.h"
#include "PhysicsTools/UtilAlgos/interface/TFileService.h"

#include "TH1F.h"
#include "TH2F.h"


namespace HbbAnalysis {//namespace

  class JetFlavour {


  public:

    //to reset counters and initialise histograms
    void Initialise(TFileDirectory & aDir, unsigned short aDebug, const unsigned int aFlavour);

    //find partons of given flavour: looking for quarks with 
    //more than one parent, and status 3
    //written for Zbb/Zcc or QCD events. If W is found, return 0.
    //Otherwise, return the number of partons found.
    unsigned int fillPartons(const edm::Handle<reco::GenParticleCollection> & aGenParticles);

    inline const flavourIni(){
      return flavour_;
    };

    inline const std::vector<unsigned int> & partons(){
      return partons_;
    };

    inline const unsigned int nPartons(){
      assert (eid_.size() == partons_.size());
      assert (muid_.size() == partons_.size());
      assert (flav_.size() == partons_.size());
      return partons_.size();
    };

    inline const reco::GenParticle & parton(const edm::Handle<reco::GenParticleCollection> & aGenParticles,unsigned int aIndex){
      assert(aIndex < partons_.size());
      return (*aGenParticles).at(partons_.at(aIndex));
    };

    inline const std::vector<unsigned int> & eids() {
      assert (eid_.size() == partons_.size());
      return eid_;
    }

    inline const unsigned int nElectronDecays(){
      assert (eid_.size() == partons_.size());
      return eid_.size();
    };

    //eid_ filled by 0 if no stable electron/muon was found
    //to ensure same size as partons_ vector...
    inline const bool hasStableElectron(unsigned int aIndex){
      return eid_.at(aIndex);
    };

    inline const bool hasStableMuon(unsigned int aIndex){
      return muid_.at(aIndex);
    };

    inline const reco::GenParticle & stableElectron(const edm::Handle<reco::GenParticleCollection> & aGenParticles,unsigned int aIndex){
      assert (eid_.size() == partons_.size());
      assert(aIndex < eid_.size());
      return (*aGenParticles).at(eid_.at(aIndex));
    };

    inline const std::vector<unsigned int> & muids() {
      assert (muid_.size() == partons_.size());
      return muid_;
    };

    inline const unsigned int nMuonDecays(){
      assert (muid_.size() == partons_.size());
      return muid_.size();
    };

    inline const reco::GenParticle & stableMuon(const edm::Handle<reco::GenParticleCollection> & aGenParticles,unsigned int aIndex){
      assert (muid_.size() == partons_.size());
      assert(aIndex < muid_.size());
      return (*aGenParticles).at(muid_.at(aIndex));
    };

    inline const std::vector<unsigned int> & flavours() {
      assert (flav_.size() == partons_.size());
      return flav_;
    };

    inline const unsigned int nFlavours(){
      assert (flav_.size() == partons_.size());
      return flav_.size();
    };

    inline const unsigned int partonFlavour(unsigned int aIndex){
      assert (flav_.size() == partons_.size());
      assert(aIndex < flav_.size());
      assert (flav_.at(aIndex) < 4);
      return flav_.at(aIndex);
    };


    std::pair<int,unsigned int> partonMatchingGenJet(const pat::Jet & aJet, const edm::Handle<reco::GenParticleCollection> & aGenParticles, double aDR);

    unsigned int leptonMatchingGenJet(const pat::Jet & aJet, const edm::Handle<reco::GenParticleCollection> & aGenParticles, double aDR);

    unsigned int partonMatchingRecoJet(const pat::Jet & aJet, const edm::Handle<reco::GenParticleCollection> & aGenParticles, double aDR);

    unsigned int leptonMatchingRecoJet(const pat::Jet & aJet, const edm::Handle<reco::GenParticleCollection> & aGenParticles, double aDR);

    bool compareParticles(const reco::Candidate * aMum, const reco::GenParticle & aP);

    void printParticle(const reco::GenParticle & aP) const;

    void printSummary();

  private:

    unsigned short debug_;

    //flavour of partons required: 3=u,d,s,g,
    // 4=c, 5=b, 21=g, 45=b and/or c.
    unsigned int flavour_;

    //vector of references to the initial partons of the hard-scattering
    //found by looking for quarks with 2 parents and status 3
    std::vector<unsigned int> partons_;

    //flavour: 
    //for b or c jets: 0=all, 1=had, 2=e,     3=mu
    //for light-jet  : 0=all,        2=u,d,s, 3=gluon
    std::vector<unsigned int> flav_;

    //vector of references to the status 1 final lepton
    std::vector<unsigned int> eid_;
    std::vector<unsigned int> muid_;

    //counters to print event summary
    unsigned int nTot_;
    unsigned int nLight_;
    unsigned int nGluon_;
    unsigned int nCharm_;
    unsigned int nBottom_;
    unsigned int nTotal_;
    unsigned int nHeavy_;
    unsigned int nSemiHeavy_;

    unsigned int nZEvts_;
    unsigned int nWEvts_;

    //control histograms
    TH1F *p_nQuarks;
    TH1F *p_nAntiQuarks;
    TH1F *p_partonFlavour;

    TH1F *p_parentOfStableFlavour;
    TH1F *p_stableFlavour;
    TH1F *p_numberOfParentsOfStable;


  };

}//namespace


#endif //HbbAnalysis_JetFlavour_hh
Revision:	1.3
Committed:	Thu Jul 2 15:38:21 2009 UTC (15 years, 10 months ago) by amagnan
Content type:	text/plain
Branch:	MAIN
Changes since 1.2:	+17 -13 lines
Log Message:	bug in jetFlavour for light jets
#	User	Rev	Content
1	amagnan	1.1	#ifndef HbbAnalysis_JetFlavour_hh
2			#define HbbAnalysis_JetFlavour_hh
3
4			#include <vector>
5			#include <string>
6			#include <cassert>
7
8			#include "DataFormats/Common/interface/Handle.h"
9			#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
10			#include "DataFormats/PatCandidates/interface/Jet.h"
11
12			#include "FWCore/ServiceRegistry/interface/Service.h"
13			#include "PhysicsTools/UtilAlgos/interface/TFileService.h"
14
15			#include "TH1F.h"
16			#include "TH2F.h"
17
18
19			namespace HbbAnalysis {//namespace
20
21			class JetFlavour {
22
23
24			public:
25
26			//to reset counters and initialise histograms
27			void Initialise(TFileDirectory & aDir, unsigned short aDebug, const unsigned int aFlavour);
28
29			//find partons of given flavour: looking for quarks with
30			//more than one parent, and status 3
31			//written for Zbb/Zcc or QCD events. If W is found, return 0.
32			//Otherwise, return the number of partons found.
33			unsigned int fillPartons(const edm::Handle<reco::GenParticleCollection> & aGenParticles);
34
35	amagnan	1.3	inline const flavourIni(){
36			return flavour_;
37			};
38
39	amagnan	1.1	inline const std::vector<unsigned int> & partons(){
40			return partons_;
41	amagnan	1.3	};
42	amagnan	1.1
43			inline const unsigned int nPartons(){
44			assert (eid_.size() == partons_.size());
45			assert (muid_.size() == partons_.size());
46			assert (flav_.size() == partons_.size());
47			return partons_.size();
48	amagnan	1.3	};
49	amagnan	1.1
50			inline const reco::GenParticle & parton(const edm::Handle<reco::GenParticleCollection> & aGenParticles,unsigned int aIndex){
51			assert(aIndex < partons_.size());
52			return (*aGenParticles).at(partons_.at(aIndex));
53	amagnan	1.3	};
54	amagnan	1.1
55			inline const std::vector<unsigned int> & eids() {
56			assert (eid_.size() == partons_.size());
57			return eid_;
58			}
59
60			inline const unsigned int nElectronDecays(){
61			assert (eid_.size() == partons_.size());
62			return eid_.size();
63	amagnan	1.3	};
64	amagnan	1.1
65	amagnan	1.2	//eid_ filled by 0 if no stable electron/muon was found
66			//to ensure same size as partons_ vector...
67			inline const bool hasStableElectron(unsigned int aIndex){
68			return eid_.at(aIndex);
69	amagnan	1.3	};
70	amagnan	1.2
71			inline const bool hasStableMuon(unsigned int aIndex){
72			return muid_.at(aIndex);
73	amagnan	1.3	};
74	amagnan	1.2
75	amagnan	1.1	inline const reco::GenParticle & stableElectron(const edm::Handle<reco::GenParticleCollection> & aGenParticles,unsigned int aIndex){
76			assert (eid_.size() == partons_.size());
77			assert(aIndex < eid_.size());
78			return (*aGenParticles).at(eid_.at(aIndex));
79	amagnan	1.3	};
80	amagnan	1.1
81			inline const std::vector<unsigned int> & muids() {
82			assert (muid_.size() == partons_.size());
83			return muid_;
84	amagnan	1.3	};
85	amagnan	1.1
86			inline const unsigned int nMuonDecays(){
87			assert (muid_.size() == partons_.size());
88			return muid_.size();
89	amagnan	1.3	};
90	amagnan	1.1
91			inline const reco::GenParticle & stableMuon(const edm::Handle<reco::GenParticleCollection> & aGenParticles,unsigned int aIndex){
92			assert (muid_.size() == partons_.size());
93			assert(aIndex < muid_.size());
94			return (*aGenParticles).at(muid_.at(aIndex));
95	amagnan	1.3	};
96	amagnan	1.1
97			inline const std::vector<unsigned int> & flavours() {
98			assert (flav_.size() == partons_.size());
99			return flav_;
100	amagnan	1.3	};
101	amagnan	1.1
102			inline const unsigned int nFlavours(){
103			assert (flav_.size() == partons_.size());
104			return flav_.size();
105	amagnan	1.3	};
106	amagnan	1.1
107			inline const unsigned int partonFlavour(unsigned int aIndex){
108			assert (flav_.size() == partons_.size());
109			assert(aIndex < flav_.size());
110	amagnan	1.2	assert (flav_.at(aIndex) < 4);
111	amagnan	1.1	return flav_.at(aIndex);
112	amagnan	1.3	};
113	amagnan	1.1
114
115			std::pair<int,unsigned int> partonMatchingGenJet(const pat::Jet & aJet, const edm::Handle<reco::GenParticleCollection> & aGenParticles, double aDR);
116
117			unsigned int leptonMatchingGenJet(const pat::Jet & aJet, const edm::Handle<reco::GenParticleCollection> & aGenParticles, double aDR);
118
119			unsigned int partonMatchingRecoJet(const pat::Jet & aJet, const edm::Handle<reco::GenParticleCollection> & aGenParticles, double aDR);
120
121			unsigned int leptonMatchingRecoJet(const pat::Jet & aJet, const edm::Handle<reco::GenParticleCollection> & aGenParticles, double aDR);
122
123			bool compareParticles(const reco::Candidate * aMum, const reco::GenParticle & aP);
124
125			void printParticle(const reco::GenParticle & aP) const;
126
127			void printSummary();
128
129			private:
130
131			unsigned short debug_;
132
133			//flavour of partons required: 3=u,d,s,g,
134			// 4=c, 5=b, 21=g, 45=b and/or c.
135			unsigned int flavour_;
136
137			//vector of references to the initial partons of the hard-scattering
138			//found by looking for quarks with 2 parents and status 3
139			std::vector<unsigned int> partons_;
140
141			//flavour:
142			//for b or c jets: 0=all, 1=had, 2=e, 3=mu
143			//for light-jet : 0=all, 2=u,d,s, 3=gluon
144			std::vector<unsigned int> flav_;
145
146			//vector of references to the status 1 final lepton
147			std::vector<unsigned int> eid_;
148			std::vector<unsigned int> muid_;
149
150			//counters to print event summary
151			unsigned int nTot_;
152			unsigned int nLight_;
153			unsigned int nGluon_;
154			unsigned int nCharm_;
155			unsigned int nBottom_;
156			unsigned int nTotal_;
157			unsigned int nHeavy_;
158			unsigned int nSemiHeavy_;
159
160			unsigned int nZEvts_;
161			unsigned int nWEvts_;
162
163			//control histograms
164			TH1F *p_nQuarks;
165			TH1F *p_nAntiQuarks;
166			TH1F *p_partonFlavour;
167
168			TH1F *p_parentOfStableFlavour;
169			TH1F *p_stableFlavour;
170			TH1F *p_numberOfParentsOfStable;
171
172
173			};
174
175			}//namespace
176
177
178			#endif //HbbAnalysis_JetFlavour_hh