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
#	Content
1	#include "../include/Utils.h"
2
3
4	// 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	//variable HEP Tagger from Rebekka
19
20	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
31	LorentzVector allsubjets(0,0,0,0);
32
33	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
41	mjet = allsubjets.M();
42	ptjet= allsubjets.Pt();
43
44	double m12, m13, m23;
45
46	//The subjetcs have to be three
47	if(nsubjets==3) {
48
49	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
62	} else {
63	return false;
64	}
65
66	double rmin=massWindowLower*wmass/topmass;
67	double rmax=massWindowUpper*wmass/topmass;
68
69	int keep=0;
70
71	//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
96	}
97
98
99	//HEP Tagger from Ivan
100
101	bool HepTopTag(TopJet topjet)
102	{
103	//call variable tagger with default parameters
104	return variableHepTopTag(topjet);
105
106	}
107
108	//default values (mminLower=50., mjetLower=140, mjetUpper=250.) defined in Utils.h
109	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	// mjetLower=140;
124	// mjetUpper=250;
125	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	//call variable tagger with default parameters
163	return variableTopTag(topjet, mjet, nsubjets, mmin);
164
165	}
166
167	Jet* nextJet(const Particle p, std::vector<Jet> jets)
168	{
169
170	double deltarmin = double_infinity();
171	Jet* nextjet=0;
172	for(unsigned int i=0; i<jets->size(); ++i) {
173	Jet ji = jets->at(i);
174	if (fabs(p->pt() - ji.pt())<1e-8) continue; // skip identical particle
175	if(jets->at(i).deltaR(*p) < deltarmin) {
176	deltarmin = jets->at(i).deltaR(*p);
177	nextjet = &jets->at(i);
178	}
179	}
180
181	return nextjet;
182	}
183
184	bool WTag(TopJet prunedjet, double& mjet, int &nsubjets, double& massdrop)
185	{
186
187	nsubjets=prunedjet.numberOfDaughters();
188
189	mjet = 0;
190	if(prunedjet.v4().isTimelike())
191	mjet = prunedjet.v4().M();
192
193	//calculate mass drop for first sub-jet ordered in pt
194	massdrop = 0;
195	if(nsubjets>=1 && mjet>0) {
196
197	std::vector< Particle > subjets = prunedjet.subjets();
198	sort(subjets.begin(), subjets.end(), HigherPt());
199
200	double m1 = 0;
201	if(subjets[0].v4().isTimelike())
202	m1 = subjets[0].v4().M();
203
204	massdrop = m1/mjet;
205	}
206
207	//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
216	}
217
218	double pTrel(const Particle p, std::vector<Jet> jets)
219	{
220
221	double ptrel=0;
222	Jet* nextjet = nextJet(p,jets);
223	if (!nextjet) return ptrel;
224
225	TVector3 p3(p->v4().Px(),p->v4().Py(),p->v4().Pz());
226	TVector3 jet3(nextjet->v4().Px(),nextjet->v4().Py(),nextjet->v4().Pz());
227
228	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
235	return ptrel;
236	}
237
238	double deltaRmin(const Particle p, std::vector<Jet> jets)
239	{
240	Jet* j = nextJet(p,jets);
241	double dr = 999.;
242	if (j) dr = j->deltaR(*p);
243	return dr;
244	}
245
246	TVector3 toVector(LorentzVector v4)
247	{
248
249	TVector3 v3(0,0,0);
250	v3.SetX(v4.X());
251	v3.SetY(v4.Y());
252	v3.SetZ(v4.Z());
253	return v3;
254	}
255
256	TVector3 toVector(LorentzVectorXYZE v4)
257	{
258
259	TVector3 v3(0,0,0);
260	v3.SetX(v4.X());
261	v3.SetY(v4.Y());
262	v3.SetZ(v4.Z());
263	return v3;
264	}
265
266	LorentzVectorXYZE toXYZ(LorentzVector v4)
267	{
268
269	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	}
276
277	LorentzVector toPtEtaPhi(LorentzVectorXYZE v4)
278	{
279
280	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	}
287
288	double deltaR(LorentzVector v1, LorentzVector v2)
289	{
290
291	Particle p1;
292	p1.set_v4(v1);
293	Particle p2;
294	p2.set_v4(v2);
295	return p1.deltaR(p2);
296	}
297
298	double double_infinity()
299	{
300	return std::numeric_limits<double>::infinity() ;
301	}
302
303	int int_infinity()
304	{
305	return std::numeric_limits<int>::max() ;
306	}
307
308	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	}
314
315	int JetFlavor(Jet *jet)
316	{
317
318	EventCalc* calc = EventCalc::Instance();
319
320	std::vector< GenParticle >* genparticles = calc->GetGenParticles();
321	if(genparticles) {
322
323	//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
327	GenParticle genp = genparticles->at(i);
328
329	//only take status 3 particles into account
330	//if(genp.status()!=3) continue;
331
332	if(jet->deltaR(genp)<0.5)
333	matched_genparticle_ids.push_back(genp.pdgId());
334	}
335
336	//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
357	}
358
359	//no matched GenParticle -> return default value 0
360	return 0;
361
362	}