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 "CommonTools/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 unsigned int flavourIni() const{
      return flavour_;
    };

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

    inline const unsigned int nPartons() const{
      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) const{
      assert(aIndex < partons_.size());
      return (*aGenParticles).at(partons_.at(aIndex));
    };

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

    inline const unsigned int nElectronDecays() const{
      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) const{
      return eid_.at(aIndex);
    };

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

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

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

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

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

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

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

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


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

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

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

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

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

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

    void printSummary() const;

  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.6
Committed:	Thu Mar 4 15:17:41 2010 UTC (15 years, 2 months ago) by amagnan
Content type:	text/plain
Branch:	MAIN
CVS Tags:	v01-00-00, beforeMETHacks, v00-05-03, v00-05-02, v00-05-01, v00-05-00, HbbAnaFor35X, v00-04-02, v00-04-01, v00-04-00, HEAD
Changes since 1.5:	+1 -1 lines
Log Message:	update beginJob and TFileService from CommonTools
#	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	amagnan	1.6	#include "CommonTools/UtilAlgos/interface/TFileService.h"
14	amagnan	1.1
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.5	inline const unsigned int flavourIni() const{
36	amagnan	1.3	return flavour_;
37			};
38
39	amagnan	1.5	inline const std::vector<unsigned int> & partons() const{
40	amagnan	1.1	return partons_;
41	amagnan	1.3	};
42	amagnan	1.1
43	amagnan	1.5	inline const unsigned int nPartons() const{
44	amagnan	1.1	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	amagnan	1.5	inline const reco::GenParticle & parton(const edm::Handle<reco::GenParticleCollection> & aGenParticles,unsigned int aIndex) const{
51	amagnan	1.1	assert(aIndex < partons_.size());
52			return (*aGenParticles).at(partons_.at(aIndex));
53	amagnan	1.3	};
54	amagnan	1.1
55	amagnan	1.5	inline const std::vector<unsigned int> & eids() const{
56	amagnan	1.1	assert (eid_.size() == partons_.size());
57			return eid_;
58			}
59
60	amagnan	1.5	inline const unsigned int nElectronDecays() const{
61	amagnan	1.1	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	amagnan	1.5	inline const bool hasStableElectron(unsigned int aIndex) const{
68	amagnan	1.2	return eid_.at(aIndex);
69	amagnan	1.3	};
70	amagnan	1.2
71	amagnan	1.5	inline const bool hasStableMuon(unsigned int aIndex) const{
72	amagnan	1.2	return muid_.at(aIndex);
73	amagnan	1.3	};
74	amagnan	1.2
75	amagnan	1.5	inline const reco::GenParticle & stableElectron(const edm::Handle<reco::GenParticleCollection> & aGenParticles,unsigned int aIndex) const{
76	amagnan	1.1	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	amagnan	1.5	inline const std::vector<unsigned int> & muids() const{
82	amagnan	1.1	assert (muid_.size() == partons_.size());
83			return muid_;
84	amagnan	1.3	};
85	amagnan	1.1
86	amagnan	1.5	inline const unsigned int nMuonDecays() const{
87	amagnan	1.1	assert (muid_.size() == partons_.size());
88			return muid_.size();
89	amagnan	1.3	};
90	amagnan	1.1
91	amagnan	1.5	inline const reco::GenParticle & stableMuon(const edm::Handle<reco::GenParticleCollection> & aGenParticles,unsigned int aIndex) const{
92	amagnan	1.1	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	amagnan	1.5	inline const std::vector<unsigned int> & flavours() const{
98	amagnan	1.1	assert (flav_.size() == partons_.size());
99			return flav_;
100	amagnan	1.3	};
101	amagnan	1.1
102	amagnan	1.5	inline const unsigned int nFlavours() const{
103	amagnan	1.1	assert (flav_.size() == partons_.size());
104			return flav_.size();
105	amagnan	1.3	};
106	amagnan	1.1
107	amagnan	1.5	inline const unsigned int partonFlavour(unsigned int aIndex) const{
108	amagnan	1.1	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	amagnan	1.5	const std::pair<int,unsigned int> partonMatchingGenJet(const pat::Jet & aJet, const edm::Handle<reco::GenParticleCollection> & aGenParticles, double aDR) const ;
116	amagnan	1.1
117	amagnan	1.5	const unsigned int leptonMatchingGenJet(const pat::Jet & aJet, const edm::Handle<reco::GenParticleCollection> & aGenParticles, double aDR) const ;
118	amagnan	1.1
119	amagnan	1.5	const unsigned int partonMatchingRecoJet(const pat::Jet & aJet, const edm::Handle<reco::GenParticleCollection> & aGenParticles, double aDR) const ;
120	amagnan	1.1
121	amagnan	1.5	const unsigned int leptonMatchingRecoJet(const pat::Jet & aJet, const edm::Handle<reco::GenParticleCollection> & aGenParticles, double aDR) const ;
122	amagnan	1.1
123			bool compareParticles(const reco::Candidate * aMum, const reco::GenParticle & aP);
124
125			void printParticle(const reco::GenParticle & aP) const;
126
127	amagnan	1.5	void printSummary() const;
128	amagnan	1.1
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