SFrameTools/src/Utils.cxx

#include "../include/Utils.h"


// global function to define a tagged jet

bool IsTagged(Jet & jet, E_BtagType type)
{
    if(type==e_CSVL && jet.btag_combinedSecondaryVertex()>0.244) return true;
    if(type==e_CSVM && jet.btag_combinedSecondaryVertex()>0.679) return true;
    if(type==e_CSVT && jet.btag_combinedSecondaryVertex()>0.898) return true;
    if(type==e_JPL && jet.btag_jetProbability()>0.275) return true;
    if(type==e_JPM && jet.btag_jetProbability()>0.545) return true;
    if(type==e_JPT && jet.btag_jetProbability()>0.790) return true;       
    
    return false;
}

//variable HEP Tagger from Rebekka

bool variableHepTopTag(TopJet topjet, double ptJetMin, double massWindowLower, double massWindowUpper, double cutCondition2, double cutCondition3)
{
    double mjet;
    double ptjet;
    int nsubjets;

    double topmass=172.3;
    double wmass=80.4;

    nsubjets=topjet.numberOfDaughters();

    LorentzVector allsubjets(0,0,0,0);

    for(int j=0; j<topjet.numberOfDaughters(); ++j) {
        allsubjets += topjet.subjets()[j].v4();
    }
    if(!allsubjets.isTimelike()) {
        mjet=0;
        return false;
    }

    mjet = allsubjets.M();
    ptjet= allsubjets.Pt();

    double m12, m13, m23;

    //The subjetcs have to be three
    if(nsubjets==3) {

        std::vector<Particle> subjets = topjet.subjets();
        sort(subjets.begin(), subjets.end(), HigherPt());

        m12 = 0;
        if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
            m12=(subjets[0].v4()+subjets[1].v4()).M();
        m13 = 0;
        if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
            m13=(subjets[0].v4()+subjets[2].v4()).M();
        m23 = 0;
        if( (subjets[1].v4()+subjets[2].v4()).isTimelike()  )
            m23 = (subjets[1].v4()+subjets[2].v4()).M();

    } else {
        return false;
    }

    double rmin=massWindowLower*wmass/topmass;
    double rmax=massWindowUpper*wmass/topmass;

    int keep=0;

    //Conditions on the subjects: at least one has to be true
    //1 condition
    if(atan(m13/m12)>0.2 && atan(m13/m12)<1.3 && m23/mjet>rmin && m23/mjet<rmax) keep=1;

    //2 condition
    double cond2left=pow(rmin,2)*(1+pow((m13/m12),2));
    double cond2cent=1-pow(m23/mjet,2);
    double cond2right=pow(rmax,2)*(1+pow(m13/m12,2));

    if(cond2left<cond2cent && cond2cent<cond2right && m23/mjet>cutCondition2) keep=1;

    //3 condition
    double cond3left=pow(rmin,2)*(1+pow((m12/m13),2));
    double cond3cent=1-pow(m23/mjet,2);
    double cond3right=pow(rmax,2)*(1+pow(m12/m13,2));

    if(cond3left<cond3cent && cond3cent<cond3right && m23/mjet>cutCondition3) keep=1;

    //Final requirement: at least one of the three subjets conditions and total pt
    if(keep==1 && ptjet>ptJetMin) {
        return true;
    } else {
        return false;
    }

}


//HEP Tagger from Ivan

bool HepTopTag(TopJet topjet)
{
  //call variable tagger with default parameters
  return variableHepTopTag(topjet);

}

//default values (mminLower=50., mjetLower=140, mjetUpper=250.) defined in Utils.h
bool variableTopTag(TopJet topjet, double &mjet, int &nsubjets, double &mmin, double mminLower, double mjetLower, double mjetUpper)
{

    nsubjets=topjet.numberOfDaughters();

    LorentzVector allsubjets(0,0,0,0);

    for(int j=0; j<topjet.numberOfDaughters(); ++j) {
        allsubjets += topjet.subjets()[j].v4();
    }
    if(!allsubjets.isTimelike()) {
        mjet=0;
        mmin=0;
//  mminLower=50;
//   mjetLower=140;
//   mjetUpper=250;
        return false;
    }

    mjet = allsubjets.M();

    if(nsubjets>=3) {

        std::vector<Particle> subjets = topjet.subjets();
        sort(subjets.begin(), subjets.end(), HigherPt());

        double m01 = 0;
        if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
            m01=(subjets[0].v4()+subjets[1].v4()).M();
        double m02 = 0;
        if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
            m02=(subjets[0].v4()+subjets[2].v4()).M();
        double m12 = 0;
        if( (subjets[1].v4()+subjets[2].v4()).isTimelike() )
            m12 = (subjets[1].v4()+subjets[2].v4()).M();

        //minimum pairwise mass
        mmin = std::min(m01,std::min(m02,m12));
    }

    //at least 3 sub-jets
    if(nsubjets<3) return false;
    //minimum pairwise mass > 50 GeV/c^2
    if(mmin<mminLower) return false;
    //jet mass between 140 and 250 GeV/c^2
    if(mjet<mjetLower || mjet>mjetUpper) return false;

    return true;
}


bool TopTag(TopJet topjet,  double &mjet, int &nsubjets, double &mmin)
{
  //call variable tagger with default parameters
  return variableTopTag(topjet, mjet, nsubjets, mmin);

}
 
Jet* nextJet(const Particle *p, std::vector<Jet> *jets)
{

    double deltarmin = double_infinity();
    Jet* nextjet=0;
    for(unsigned int i=0; i<jets->size(); ++i) {
        Jet ji = jets->at(i);
        if (fabs(p->pt() - ji.pt())<1e-8) continue; // skip identical particle
        if(jets->at(i).deltaR(*p) < deltarmin) {
            deltarmin = jets->at(i).deltaR(*p);
            nextjet = &jets->at(i);
        }
    }

    return nextjet;
}

bool WTag(TopJet prunedjet,  double& mjet, int &nsubjets, double& massdrop)
{

    nsubjets=prunedjet.numberOfDaughters();

    mjet = 0;
    if(prunedjet.v4().isTimelike())
        mjet = prunedjet.v4().M();

    //calculate mass drop for first sub-jet ordered in pt
    massdrop = 0;
    if(nsubjets>=1 && mjet>0) {

        std::vector< Particle > subjets = prunedjet.subjets();
        sort(subjets.begin(), subjets.end(), HigherPt());

        double m1 = 0;
        if(subjets[0].v4().isTimelike())
            m1 = subjets[0].v4().M();

        massdrop = m1/mjet;
    }

    //at least 2 sub-jets
    if(nsubjets<2) return false;
    //60 GeV < pruned jet mass < 100 GeV
    if(mjet <= 60 || mjet >= 100) return false;
    //mass drop < 0.4
    if(massdrop>=0.4) return false;

    return true;

}

double pTrel(const Particle *p, std::vector<Jet> *jets)
{

    double ptrel=0;
    Jet* nextjet =  nextJet(p,jets);
    if (!nextjet) return ptrel;

    TVector3 p3(p->v4().Px(),p->v4().Py(),p->v4().Pz());
    TVector3 jet3(nextjet->v4().Px(),nextjet->v4().Py(),nextjet->v4().Pz());

    if(p3.Mag()!=0 && jet3.Mag()!=0) {
        double sin_alpha = (p3.Cross(jet3)).Mag()/p3.Mag()/jet3.Mag();
        ptrel = p3.Mag()*sin_alpha;
    } else {
        std::cout << "something strange happend in the ptrel calculation: either lepton or jet momentum is 0" <<std::endl;
    }

    return ptrel;
}

double deltaRmin(const Particle *p, std::vector<Jet> *jets)
{
    Jet* j = nextJet(p,jets);
    double dr = 999.;
    if (j) dr = j->deltaR(*p);
    return dr;
}

TVector3 toVector(LorentzVector v4)
{

    TVector3 v3(0,0,0);
    v3.SetX(v4.X());
    v3.SetY(v4.Y());
    v3.SetZ(v4.Z());
    return v3;
}

TVector3 toVector(LorentzVectorXYZE v4)
{

    TVector3 v3(0,0,0);
    v3.SetX(v4.X());
    v3.SetY(v4.Y());
    v3.SetZ(v4.Z());
    return v3;
}

LorentzVectorXYZE toXYZ(LorentzVector v4)
{

    LorentzVectorXYZE v4_new(0,0,0,0);
    v4_new.SetPx(v4.Px());
    v4_new.SetPy(v4.Py());
    v4_new.SetPz(v4.Pz());
    v4_new.SetE(v4.E());
    return v4_new;
}

LorentzVector toPtEtaPhi(LorentzVectorXYZE v4)
{

    LorentzVector v4_new(0,0,0,0);
    v4_new.SetPt(v4.Pt());
    v4_new.SetEta(v4.Eta());
    v4_new.SetPhi(v4.Phi());
    v4_new.SetE(v4.E());
    return v4_new;
}

double deltaR(LorentzVector v1, LorentzVector v2)
{

    Particle p1;
    p1.set_v4(v1);
    Particle p2;
    p2.set_v4(v2);
    return p1.deltaR(p2);
}

double double_infinity()
{
    return std::numeric_limits<double>::infinity() ;
}

int int_infinity()
{
    return std::numeric_limits<int>::max() ;
}

int myPow(int x, unsigned int p)
{
    int i = 1;
    for (unsigned int j = 1; j <= p; j++)  i *= x;
    return i;
}

int JetFlavor(Jet *jet)
{

    EventCalc* calc = EventCalc::Instance();

    std::vector< GenParticle >* genparticles = calc->GetGenParticles();
    if(genparticles) {

        //fill pdg IDs of all matched GenParticles
        std::vector<int> matched_genparticle_ids;
        for(unsigned int i=0; i<genparticles->size(); ++i) {

            GenParticle genp = genparticles->at(i);

            //only take status 3 particles into account
            //if(genp.status()!=3) continue;

            if(jet->deltaR(genp)<0.5)
                matched_genparticle_ids.push_back(genp.pdgId());
        }

        //search for b quarks first
        for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
            if(abs(matched_genparticle_ids[i])==5) return matched_genparticle_ids[i];
        }
        //no b quark -> search for c quarks
        for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
            if(abs(matched_genparticle_ids[i])==4) return matched_genparticle_ids[i];
        }
        for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
            if(abs(matched_genparticle_ids[i])==3) return matched_genparticle_ids[i];
        }
        for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
            if(abs(matched_genparticle_ids[i])==2) return matched_genparticle_ids[i];
        }
        for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
            if(abs(matched_genparticle_ids[i])==1) return matched_genparticle_ids[i];
        }
        for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
            if(abs(matched_genparticle_ids[i])==21) return matched_genparticle_ids[i];
        }

    }

    //no matched GenParticle -> return default value 0
    return 0;

}
Revision:	1.11
Committed:	Wed Feb 13 16:30:24 2013 UTC (12 years, 2 months ago) by peiffer
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Feb-15-2013-v1, Feb-14-2013
Changes since 1.10:	+1 -1 lines
Log Message:	remove status=3 in jet flavor
#	User	Rev	Content
1	peiffer	1.1	#include "../include/Utils.h"
2
3
4	bazterra	1.8	// global function to define a tagged jet
5
6			bool IsTagged(Jet & jet, E_BtagType type)
7			{
8			if(type==e_CSVL && jet.btag_combinedSecondaryVertex()>0.244) return true;
9			if(type==e_CSVM && jet.btag_combinedSecondaryVertex()>0.679) return true;
10			if(type==e_CSVT && jet.btag_combinedSecondaryVertex()>0.898) return true;
11			if(type==e_JPL && jet.btag_jetProbability()>0.275) return true;
12			if(type==e_JPM && jet.btag_jetProbability()>0.545) return true;
13			if(type==e_JPT && jet.btag_jetProbability()>0.790) return true;
14
15			return false;
16			}
17
18	hoeing	1.6	//variable HEP Tagger from Rebekka
19
20	bazterra	1.8	bool variableHepTopTag(TopJet topjet, double ptJetMin, double massWindowLower, double massWindowUpper, double cutCondition2, double cutCondition3)
21			{
22			double mjet;
23			double ptjet;
24			int nsubjets;
25
26			double topmass=172.3;
27			double wmass=80.4;
28
29			nsubjets=topjet.numberOfDaughters();
30	hoeing	1.6
31	bazterra	1.8	LorentzVector allsubjets(0,0,0,0);
32	hoeing	1.6
33	bazterra	1.8	for(int j=0; j<topjet.numberOfDaughters(); ++j) {
34			allsubjets += topjet.subjets()[j].v4();
35			}
36			if(!allsubjets.isTimelike()) {
37			mjet=0;
38			return false;
39			}
40	hoeing	1.6
41	bazterra	1.8	mjet = allsubjets.M();
42			ptjet= allsubjets.Pt();
43	hoeing	1.6
44	bazterra	1.8	double m12, m13, m23;
45	hoeing	1.6
46	bazterra	1.8	//The subjetcs have to be three
47			if(nsubjets==3) {
48	hoeing	1.6
49	bazterra	1.8	std::vector<Particle> subjets = topjet.subjets();
50			sort(subjets.begin(), subjets.end(), HigherPt());
51
52			m12 = 0;
53			if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
54			m12=(subjets[0].v4()+subjets[1].v4()).M();
55			m13 = 0;
56			if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
57			m13=(subjets[0].v4()+subjets[2].v4()).M();
58			m23 = 0;
59			if( (subjets[1].v4()+subjets[2].v4()).isTimelike() )
60			m23 = (subjets[1].v4()+subjets[2].v4()).M();
61	hoeing	1.6
62	bazterra	1.8	} else {
63			return false;
64			}
65	hoeing	1.6
66	bazterra	1.8	double rmin=massWindowLower*wmass/topmass;
67			double rmax=massWindowUpper*wmass/topmass;
68	hoeing	1.6
69	bazterra	1.8	int keep=0;
70	hoeing	1.6
71	bazterra	1.8	//Conditions on the subjects: at least one has to be true
72			//1 condition
73			if(atan(m13/m12)>0.2 && atan(m13/m12)<1.3 && m23/mjet>rmin && m23/mjet<rmax) keep=1;
74
75			//2 condition
76			double cond2left=pow(rmin,2)*(1+pow((m13/m12),2));
77			double cond2cent=1-pow(m23/mjet,2);
78			double cond2right=pow(rmax,2)*(1+pow(m13/m12,2));
79
80			if(cond2left<cond2cent && cond2cent<cond2right && m23/mjet>cutCondition2) keep=1;
81
82			//3 condition
83			double cond3left=pow(rmin,2)*(1+pow((m12/m13),2));
84			double cond3cent=1-pow(m23/mjet,2);
85			double cond3right=pow(rmax,2)*(1+pow(m12/m13,2));
86
87			if(cond3left<cond3cent && cond3cent<cond3right && m23/mjet>cutCondition3) keep=1;
88
89			//Final requirement: at least one of the three subjets conditions and total pt
90			if(keep==1 && ptjet>ptJetMin) {
91			return true;
92			} else {
93			return false;
94			}
95	hoeing	1.6
96			}
97
98
99			//HEP Tagger from Ivan
100
101	bazterra	1.8	bool HepTopTag(TopJet topjet)
102			{
103	peiffer	1.9	//call variable tagger with default parameters
104			return variableHepTopTag(topjet);
105	hoeing	1.6
106			}
107
108			//default values (mminLower=50., mjetLower=140, mjetUpper=250.) defined in Utils.h
109	bazterra	1.8	bool variableTopTag(TopJet topjet, double &mjet, int &nsubjets, double &mmin, double mminLower, double mjetLower, double mjetUpper)
110			{
111
112			nsubjets=topjet.numberOfDaughters();
113
114			LorentzVector allsubjets(0,0,0,0);
115
116			for(int j=0; j<topjet.numberOfDaughters(); ++j) {
117			allsubjets += topjet.subjets()[j].v4();
118			}
119			if(!allsubjets.isTimelike()) {
120			mjet=0;
121			mmin=0;
122			// mminLower=50;
123	hoeing	1.6	// mjetLower=140;
124			// mjetUpper=250;
125	bazterra	1.8	return false;
126			}
127
128			mjet = allsubjets.M();
129
130			if(nsubjets>=3) {
131
132			std::vector<Particle> subjets = topjet.subjets();
133			sort(subjets.begin(), subjets.end(), HigherPt());
134
135			double m01 = 0;
136			if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
137			m01=(subjets[0].v4()+subjets[1].v4()).M();
138			double m02 = 0;
139			if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
140			m02=(subjets[0].v4()+subjets[2].v4()).M();
141			double m12 = 0;
142			if( (subjets[1].v4()+subjets[2].v4()).isTimelike() )
143			m12 = (subjets[1].v4()+subjets[2].v4()).M();
144
145			//minimum pairwise mass
146			mmin = std::min(m01,std::min(m02,m12));
147			}
148
149			//at least 3 sub-jets
150			if(nsubjets<3) return false;
151			//minimum pairwise mass > 50 GeV/c^2
152			if(mmin<mminLower) return false;
153			//jet mass between 140 and 250 GeV/c^2
154			if(mjet<mjetLower \|\| mjet>mjetUpper) return false;
155
156			return true;
157			}
158
159
160			bool TopTag(TopJet topjet, double &mjet, int &nsubjets, double &mmin)
161			{
162	peiffer	1.9	//call variable tagger with default parameters
163			return variableTopTag(topjet, mjet, nsubjets, mmin);
164	bazterra	1.8
165	peiffer	1.1	}
166	peiffer	1.9
167	bazterra	1.8	Jet* nextJet(const Particle p, std::vector<Jet> jets)
168			{
169	peiffer	1.2
170	bazterra	1.8	double deltarmin = double_infinity();
171			Jet* nextjet=0;
172			for(unsigned int i=0; i<jets->size(); ++i) {
173	rkogler	1.10	Jet ji = jets->at(i);
174			if (fabs(p->pt() - ji.pt())<1e-8) continue; // skip identical particle
175	bazterra	1.8	if(jets->at(i).deltaR(*p) < deltarmin) {
176			deltarmin = jets->at(i).deltaR(*p);
177			nextjet = &jets->at(i);
178			}
179	peiffer	1.2	}
180
181	bazterra	1.8	return nextjet;
182	peiffer	1.2	}
183
184	bazterra	1.8	bool WTag(TopJet prunedjet, double& mjet, int &nsubjets, double& massdrop)
185			{
186	peiffer	1.5
187	bazterra	1.8	nsubjets=prunedjet.numberOfDaughters();
188	peiffer	1.5
189	bazterra	1.8	mjet = 0;
190			if(prunedjet.v4().isTimelike())
191			mjet = prunedjet.v4().M();
192	peiffer	1.5
193	bazterra	1.8	//calculate mass drop for first sub-jet ordered in pt
194			massdrop = 0;
195			if(nsubjets>=1 && mjet>0) {
196	peiffer	1.5
197	bazterra	1.8	std::vector< Particle > subjets = prunedjet.subjets();
198			sort(subjets.begin(), subjets.end(), HigherPt());
199	peiffer	1.5
200	bazterra	1.8	double m1 = 0;
201			if(subjets[0].v4().isTimelike())
202			m1 = subjets[0].v4().M();
203	peiffer	1.5
204	bazterra	1.8	massdrop = m1/mjet;
205			}
206	peiffer	1.5
207	bazterra	1.8	//at least 2 sub-jets
208			if(nsubjets<2) return false;
209			//60 GeV < pruned jet mass < 100 GeV
210			if(mjet <= 60 \|\| mjet >= 100) return false;
211			//mass drop < 0.4
212			if(massdrop>=0.4) return false;
213
214			return true;
215	peiffer	1.5
216			}
217
218	bazterra	1.8	double pTrel(const Particle p, std::vector<Jet> jets)
219			{
220	peiffer	1.2
221	bazterra	1.8	double ptrel=0;
222			Jet* nextjet = nextJet(p,jets);
223	rkogler	1.10	if (!nextjet) return ptrel;
224	peiffer	1.2
225	bazterra	1.8	TVector3 p3(p->v4().Px(),p->v4().Py(),p->v4().Pz());
226			TVector3 jet3(nextjet->v4().Px(),nextjet->v4().Py(),nextjet->v4().Pz());
227	peiffer	1.2
228	bazterra	1.8	if(p3.Mag()!=0 && jet3.Mag()!=0) {
229			double sin_alpha = (p3.Cross(jet3)).Mag()/p3.Mag()/jet3.Mag();
230			ptrel = p3.Mag()*sin_alpha;
231			} else {
232			std::cout << "something strange happend in the ptrel calculation: either lepton or jet momentum is 0" <<std::endl;
233			}
234	peiffer	1.2
235	bazterra	1.8	return ptrel;
236	peiffer	1.2	}
237
238	bazterra	1.8	double deltaRmin(const Particle p, std::vector<Jet> jets)
239			{
240	rkogler	1.10	Jet* j = nextJet(p,jets);
241			double dr = 999.;
242			if (j) dr = j->deltaR(*p);
243			return dr;
244	peiffer	1.2	}
245
246	bazterra	1.8	TVector3 toVector(LorentzVector v4)
247			{
248	peiffer	1.4
249	bazterra	1.8	TVector3 v3(0,0,0);
250			v3.SetX(v4.X());
251			v3.SetY(v4.Y());
252			v3.SetZ(v4.Z());
253			return v3;
254	peiffer	1.4	}
255
256	bazterra	1.8	TVector3 toVector(LorentzVectorXYZE v4)
257			{
258	peiffer	1.4
259	bazterra	1.8	TVector3 v3(0,0,0);
260			v3.SetX(v4.X());
261			v3.SetY(v4.Y());
262			v3.SetZ(v4.Z());
263			return v3;
264	peiffer	1.4	}
265
266	bazterra	1.8	LorentzVectorXYZE toXYZ(LorentzVector v4)
267			{
268	peiffer	1.4
269	bazterra	1.8	LorentzVectorXYZE v4_new(0,0,0,0);
270			v4_new.SetPx(v4.Px());
271			v4_new.SetPy(v4.Py());
272			v4_new.SetPz(v4.Pz());
273			v4_new.SetE(v4.E());
274			return v4_new;
275	peiffer	1.4	}
276
277	bazterra	1.8	LorentzVector toPtEtaPhi(LorentzVectorXYZE v4)
278			{
279	peiffer	1.4
280	bazterra	1.8	LorentzVector v4_new(0,0,0,0);
281			v4_new.SetPt(v4.Pt());
282			v4_new.SetEta(v4.Eta());
283			v4_new.SetPhi(v4.Phi());
284			v4_new.SetE(v4.E());
285			return v4_new;
286	peiffer	1.4	}
287
288	bazterra	1.8	double deltaR(LorentzVector v1, LorentzVector v2)
289			{
290	peiffer	1.4
291	bazterra	1.8	Particle p1;
292			p1.set_v4(v1);
293			Particle p2;
294			p2.set_v4(v2);
295			return p1.deltaR(p2);
296	peiffer	1.4	}
297	peiffer	1.2
298	bazterra	1.8	double double_infinity()
299			{
300			return std::numeric_limits<double>::infinity() ;
301	peiffer	1.2	}
302
303	bazterra	1.8	int int_infinity()
304			{
305			return std::numeric_limits<int>::max() ;
306	peiffer	1.2	}
307	peiffer	1.4
308	bazterra	1.8	int myPow(int x, unsigned int p)
309			{
310			int i = 1;
311			for (unsigned int j = 1; j <= p; j++) i *= x;
312			return i;
313	peiffer	1.4	}
314	peiffer	1.7
315	bazterra	1.8	int JetFlavor(Jet *jet)
316			{
317
318			EventCalc* calc = EventCalc::Instance();
319	peiffer	1.7
320	bazterra	1.8	std::vector< GenParticle >* genparticles = calc->GetGenParticles();
321			if(genparticles) {
322	peiffer	1.7
323	bazterra	1.8	//fill pdg IDs of all matched GenParticles
324			std::vector<int> matched_genparticle_ids;
325			for(unsigned int i=0; i<genparticles->size(); ++i) {
326	peiffer	1.7
327	bazterra	1.8	GenParticle genp = genparticles->at(i);
328
329			//only take status 3 particles into account
330	peiffer	1.11	//if(genp.status()!=3) continue;
331	peiffer	1.7
332	bazterra	1.8	if(jet->deltaR(genp)<0.5)
333			matched_genparticle_ids.push_back(genp.pdgId());
334			}
335	peiffer	1.7
336	bazterra	1.8	//search for b quarks first
337			for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
338			if(abs(matched_genparticle_ids[i])==5) return matched_genparticle_ids[i];
339			}
340			//no b quark -> search for c quarks
341			for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
342			if(abs(matched_genparticle_ids[i])==4) return matched_genparticle_ids[i];
343			}
344			for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
345			if(abs(matched_genparticle_ids[i])==3) return matched_genparticle_ids[i];
346			}
347			for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
348			if(abs(matched_genparticle_ids[i])==2) return matched_genparticle_ids[i];
349			}
350			for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
351			if(abs(matched_genparticle_ids[i])==1) return matched_genparticle_ids[i];
352			}
353			for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
354			if(abs(matched_genparticle_ids[i])==21) return matched_genparticle_ids[i];
355			}
356	peiffer	1.7
357			}
358
359	bazterra	1.8	//no matched GenParticle -> return default value 0
360			return 0;
361	peiffer	1.7
362			}