ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/UHHAnalysis/SFrameTools/src/Utils.cxx

Comparing UserCode/UHHAnalysis/SFrameTools/src/Utils.cxx (file contents):
Revision 1.2 by peiffer, Thu May 31 08:56:47 2012 UTC vs.
Revision 1.8 by bazterra, Wed Dec 19 00:02:35 2012 UTC

#	Line 1 | Line 1
1		#include "../include/Utils.h"
2
3
4	<	bool TopTag(TopJet topjet,  double &mjet, int &nsubjets, double &mmin){
5	<
6	<	nsubjets=topjet.numberOfDaughters();
7	<
8	<	LorentzVector allsubjets(0,0,0,0);
9	<
10	<	for(int j=0; j<topjet.numberOfDaughters(); ++j){
11	<	allsubjets += topjet.subjets()[j].v4();
12	<	}
13	<	if(!allsubjets.isTimelike()){
14	<	mjet=0;
15	<	mmin=0;
16	<	return false;
17	<	}
4	>	// global function to define a tagged jet
5
6	<	mjet = allsubjets.M();
7	<
8	<	if(nsubjets>=3) {
9	<
10	<	std::vector<Particle> subjets = topjet.subjets();
11	<	sort(subjets.begin(), subjets.end(), HigherPt());
12	<
13	<	double m01 = 0;
27	<	if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
28	<	m01=(subjets[0].v4()+subjets[1].v4()).M();
29	<	double m02 = 0;
30	<	if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
31	<	m02=(subjets[0].v4()+subjets[2].v4()).M();
32	<	double m12 = 0;
33	<	if( (subjets[1].v4()+subjets[2].v4()).isTimelike()  )
34	<	m12 = (subjets[1].v4()+subjets[2].v4()).M();
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	+
100	+
101	+
102	+
103	+	//HEP Tagger from Ivan
104	+
105	+	bool HepTopTag(TopJet topjet)
106	+	{
107	+
108	+	double mjet;
109	+	double ptjet;
110	+	int nsubjets;
111	+
112	+	double topmass=172.3;
113	+	double wmass=80.4;
114	+
115	+	nsubjets=topjet.numberOfDaughters();
116	+
117	+	LorentzVector allsubjets(0,0,0,0);
118	+
119	+	for(int j=0; j<topjet.numberOfDaughters(); ++j) {
120	+	allsubjets += topjet.subjets()[j].v4();
121	+	}
122	+	if(!allsubjets.isTimelike()) {
123	+	mjet=0;
124	+	return false;
125	+	}
126	+
127	+	mjet = allsubjets.M();
128	+	ptjet= allsubjets.Pt();
129	+
130	+	double m12, m13, m23;
131	+
132	+	//The subjetcs have to be three
133	+	if(nsubjets==3) {
134	+
135	+	std::vector<Particle> subjets = topjet.subjets();
136	+	sort(subjets.begin(), subjets.end(), HigherPt());
137	+
138	+	m12 = 0;
139	+	if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
140	+	m12=(subjets[0].v4()+subjets[1].v4()).M();
141	+	m13 = 0;
142	+	if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
143	+	m13=(subjets[0].v4()+subjets[2].v4()).M();
144	+	m23 = 0;
145	+	if( (subjets[1].v4()+subjets[2].v4()).isTimelike()  )
146	+	m23 = (subjets[1].v4()+subjets[2].v4()).M();
147	+
148	+	} else {
149	+	return false;
150	+	}
151	+
152	+	double rmin=0.85*wmass/topmass;
153	+	double rmax=1.15*wmass/topmass;
154	+
155	+	int keep=0;
156	+
157	+	//Conditions on the subjects: at least one has to be true
158	+	//1 condition
159	+	if(atan(m13/m12)>0.2 && atan(m13/m12)<1.3 && m23/mjet>rmin && m23/mjet<rmax) keep=1;
160	+
161	+	//2 condition
162	+	double cond2left=pow(rmin,2)*(1+pow((m13/m12),2));
163	+	double cond2cent=1-pow(m23/mjet,2);
164	+	double cond2right=pow(rmax,2)*(1+pow(m13/m12,2));
165	+
166	+	if(cond2left<cond2cent && cond2cent<cond2right && m23/mjet>0.35) keep=1;
167	+
168	+	//3 condition
169	+	double cond3left=pow(rmin,2)*(1+pow((m12/m13),2));
170	+	double cond3cent=1-pow(m23/mjet,2);
171	+	double cond3right=pow(rmax,2)*(1+pow(m12/m13,2));
172	+
173	+	if(cond3left<cond3cent && cond3cent<cond3right && m23/mjet>0.35) keep=1;
174	+
175	+	//Final requirement: at least one of the three subjets conditions and total pt
176	+	if(keep==1 && ptjet>200) {
177	+	return true;
178	+	} else {
179	+	return false;
180	+	}
181	+
182	+	}
183	+
184	+
185	+	//default values (mminLower=50., mjetLower=140, mjetUpper=250.) defined in Utils.h
186	+	bool variableTopTag(TopJet topjet, double &mjet, int &nsubjets, double &mmin, double mminLower, double mjetLower, double mjetUpper)
187	+	{
188	+
189	+	nsubjets=topjet.numberOfDaughters();
190	+
191	+	LorentzVector allsubjets(0,0,0,0);
192	+
193	+	for(int j=0; j<topjet.numberOfDaughters(); ++j) {
194	+	allsubjets += topjet.subjets()[j].v4();
195	+	}
196	+	if(!allsubjets.isTimelike()) {
197	+	mjet=0;
198	+	mmin=0;
199	+	//  mminLower=50;
200	+	//   mjetLower=140;
201	+	//   mjetUpper=250;
202	+	return false;
203	+	}
204	+
205	+	mjet = allsubjets.M();
206	+
207	+	if(nsubjets>=3) {
208	+
209	+	std::vector<Particle> subjets = topjet.subjets();
210	+	sort(subjets.begin(), subjets.end(), HigherPt());
211	+
212	+	double m01 = 0;
213	+	if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
214	+	m01=(subjets[0].v4()+subjets[1].v4()).M();
215	+	double m02 = 0;
216	+	if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
217	+	m02=(subjets[0].v4()+subjets[2].v4()).M();
218	+	double m12 = 0;
219	+	if( (subjets[1].v4()+subjets[2].v4()).isTimelike() )
220	+	m12 = (subjets[1].v4()+subjets[2].v4()).M();
221	+
222	+	//minimum pairwise mass
223	+	mmin = std::min(m01,std::min(m02,m12));
224	+	}
225	+
226	+	//at least 3 sub-jets
227	+	if(nsubjets<3) return false;
228		//minimum pairwise mass > 50 GeV/c^2
229	<	mmin = std::min(m01,std::min(m02,m12));
230	<	}
229	>	if(mmin<mminLower) return false;
230	>	//jet mass between 140 and 250 GeV/c^2
231	>	if(mjet<mjetLower || mjet>mjetUpper) return false;
232
233	<	//at least 3 sub-jets
41	<	if(nsubjets<3) return false;
42	<	//minimum pairwise mass > 50 GeV/c^2
43	<	if(mmin<50) return false;
44	<	//jet mass between 140 and 250 GeV/c^2
45	<	if(mjet<140 || mjet>250) return false;
46	<
47	<	return true;
233	>	return true;
234		}
235
50	–	double HTlep(const BaseCycleContainer *bcc){
236
237	<	double htlep=0;
237	>	bool TopTag(TopJet topjet,  double &mjet, int &nsubjets, double &mmin)
238	>	{
239	>
240	>	nsubjets=topjet.numberOfDaughters();
241
242	<	if(bcc->electrons){
243	<	for(unsigned int i=0; i<bcc->electrons->size(); ++i){
244	<	htlep += bcc->electrons->at(i).pt();
242	>	LorentzVector allsubjets(0,0,0,0);
243	>
244	>	for(int j=0; j<topjet.numberOfDaughters(); ++j) {
245	>	allsubjets += topjet.subjets()[j].v4();
246	>	}
247	>	if(!allsubjets.isTimelike()) {
248	>	mjet=0;
249	>	mmin=0;
250	>	return false;
251		}
252	<	}
253	<	if(bcc->muons){
254	<	for(unsigned int i=0; i<bcc->muons->size(); ++i){
255	<	htlep += bcc->muons->at(i).pt();
252	>
253	>	mjet = allsubjets.M();
254	>
255	>	if(nsubjets>=3) {
256	>
257	>	std::vector<Particle> subjets = topjet.subjets();
258	>	sort(subjets.begin(), subjets.end(), HigherPt());
259	>
260	>	double m01 = 0;
261	>	if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
262	>	m01=(subjets[0].v4()+subjets[1].v4()).M();
263	>	double m02 = 0;
264	>	if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
265	>	m02=(subjets[0].v4()+subjets[2].v4()).M();
266	>	double m12 = 0;
267	>	if( (subjets[1].v4()+subjets[2].v4()).isTimelike()  )
268	>	m12 = (subjets[1].v4()+subjets[2].v4()).M();
269	>
270	>	//minimum pairwise mass
271	>	mmin = std::min(m01,std::min(m02,m12));
272		}
63	–	}
64	–	if(bcc->met) htlep += bcc->met->pt();
273
274	<	return htlep;
274	>	//at least 3 sub-jets
275	>	if(nsubjets<3) return false;
276	>	//minimum pairwise mass > 50 GeV/c^2
277	>	if(mmin<50) return false;
278	>	//jet mass between 140 and 250 GeV/c^2
279	>	if(mjet<140 || mjet>250) return false;
280
281	+	return true;
282		}
283
284	<	Jet* nextJet(const Particle *p, std::vector<Jet> *jets){
284	>	Jet* nextJet(const Particle *p, std::vector<Jet> *jets)
285	>	{
286
287	<	double deltarmin = double_infinity();
288	<	Jet* nextjet=0;
289	<	for(unsigned int i=0; i<jets->size();++i){
290	<	if(jets->at(i).deltaR(*p) < deltarmin){
291	<	deltarmin = jets->at(i).deltaR(*p);
292	<	nextjet = &jets->at(i);
287	>	double deltarmin = double_infinity();
288	>	Jet* nextjet=0;
289	>	for(unsigned int i=0; i<jets->size(); ++i) {
290	>	if(jets->at(i).deltaR(*p) < deltarmin) {
291	>	deltarmin = jets->at(i).deltaR(*p);
292	>	nextjet = &jets->at(i);
293	>	}
294		}
79	–	}
295
296	<	return nextjet;
296	>	return nextjet;
297		}
298
299	<	double pTrel(const Particle *p, std::vector<Jet> *jets){
299	>	bool WTag(TopJet prunedjet,  double& mjet, int &nsubjets, double& massdrop)
300	>	{
301	>
302	>	nsubjets=prunedjet.numberOfDaughters();
303	>
304	>	mjet = 0;
305	>	if(prunedjet.v4().isTimelike())
306	>	mjet = prunedjet.v4().M();
307	>
308	>	//calculate mass drop for first sub-jet ordered in pt
309	>	massdrop = 0;
310	>	if(nsubjets>=1 && mjet>0) {
311	>
312	>	std::vector< Particle > subjets = prunedjet.subjets();
313	>	sort(subjets.begin(), subjets.end(), HigherPt());
314
315	<	double ptrel=0;
315	>	double m1 = 0;
316	>	if(subjets[0].v4().isTimelike())
317	>	m1 = subjets[0].v4().M();
318
319	<	Jet* nextjet =  nextJet(p,jets);
319	>	massdrop = m1/mjet;
320	>	}
321
322	<	TVector3 p3(p->v4().Px(),p->v4().Py(),p->v4().Pz());
323	<	TVector3 jet3(nextjet->v4().Px(),nextjet->v4().Py(),nextjet->v4().Pz());
322	>	//at least 2 sub-jets
323	>	if(nsubjets<2) return false;
324	>	//60 GeV < pruned jet mass < 100 GeV
325	>	if(mjet <= 60 || mjet >= 100) return false;
326	>	//mass drop < 0.4
327	>	if(massdrop>=0.4) return false;
328
329	<	if(p3.Mag()!=0 && jet3.Mag()!=0){
94	<	double sin_alpha = (p3.Cross(jet3)).Mag()/p3.Mag()/jet3.Mag();
95	<	ptrel = p3.Mag()*sin_alpha;
96	<	}
97	<	else{
98	<	std::cout << "something strange happend in the ptrel calculation: either lepton or jet momentum is 0" <<std::endl;
99	<	}
329	>	return true;
330
101	–	return ptrel;
331		}
332
333	<	double deltaRmin(const Particle *p, std::vector<Jet> *jets){
334	<	return nextJet(p,jets)->deltaR(*p);
333	>	double pTrel(const Particle *p, std::vector<Jet> *jets)
334	>	{
335	>
336	>	double ptrel=0;
337	>
338	>	Jet* nextjet =  nextJet(p,jets);
339	>
340	>	TVector3 p3(p->v4().Px(),p->v4().Py(),p->v4().Pz());
341	>	TVector3 jet3(nextjet->v4().Px(),nextjet->v4().Py(),nextjet->v4().Pz());
342	>
343	>	if(p3.Mag()!=0 && jet3.Mag()!=0) {
344	>	double sin_alpha = (p3.Cross(jet3)).Mag()/p3.Mag()/jet3.Mag();
345	>	ptrel = p3.Mag()*sin_alpha;
346	>	} else {
347	>	std::cout << "something strange happend in the ptrel calculation: either lepton or jet momentum is 0" <<std::endl;
348	>	}
349	>
350	>	return ptrel;
351		}
352
353	+	double deltaRmin(const Particle *p, std::vector<Jet> *jets)
354	+	{
355	+	return nextJet(p,jets)->deltaR(*p);
356	+	}
357	+
358	+	TVector3 toVector(LorentzVector v4)
359	+	{
360	+
361	+	TVector3 v3(0,0,0);
362	+	v3.SetX(v4.X());
363	+	v3.SetY(v4.Y());
364	+	v3.SetZ(v4.Z());
365	+	return v3;
366	+	}
367	+
368	+	TVector3 toVector(LorentzVectorXYZE v4)
369	+	{
370	+
371	+	TVector3 v3(0,0,0);
372	+	v3.SetX(v4.X());
373	+	v3.SetY(v4.Y());
374	+	v3.SetZ(v4.Z());
375	+	return v3;
376	+	}
377	+
378	+	LorentzVectorXYZE toXYZ(LorentzVector v4)
379	+	{
380	+
381	+	LorentzVectorXYZE v4_new(0,0,0,0);
382	+	v4_new.SetPx(v4.Px());
383	+	v4_new.SetPy(v4.Py());
384	+	v4_new.SetPz(v4.Pz());
385	+	v4_new.SetE(v4.E());
386	+	return v4_new;
387	+	}
388	+
389	+	LorentzVector toPtEtaPhi(LorentzVectorXYZE v4)
390	+	{
391	+
392	+	LorentzVector v4_new(0,0,0,0);
393	+	v4_new.SetPt(v4.Pt());
394	+	v4_new.SetEta(v4.Eta());
395	+	v4_new.SetPhi(v4.Phi());
396	+	v4_new.SetE(v4.E());
397	+	return v4_new;
398	+	}
399	+
400	+	double deltaR(LorentzVector v1, LorentzVector v2)
401	+	{
402
403	<	double double_infinity(){
404	<	return std::numeric_limits<double>::infinity() ;
403	>	Particle p1;
404	>	p1.set_v4(v1);
405	>	Particle p2;
406	>	p2.set_v4(v2);
407	>	return p1.deltaR(p2);
408		}
409
410	<	int int_infinity(){
411	<	return std::numeric_limits<int>::max() ;
410	>	double double_infinity()
411	>	{
412	>	return std::numeric_limits<double>::infinity() ;
413	>	}
414	>
415	>	int int_infinity()
416	>	{
417	>	return std::numeric_limits<int>::max() ;
418	>	}
419	>
420	>	int myPow(int x, unsigned int p)
421	>	{
422	>	int i = 1;
423	>	for (unsigned int j = 1; j <= p; j++)  i *= x;
424	>	return i;
425	>	}
426	>
427	>	int JetFlavor(Jet *jet)
428	>	{
429	>
430	>	EventCalc* calc = EventCalc::Instance();
431	>
432	>	std::vector< GenParticle >* genparticles = calc->GetGenParticles();
433	>	if(genparticles) {
434	>
435	>	//fill pdg IDs of all matched GenParticles
436	>	std::vector<int> matched_genparticle_ids;
437	>	for(unsigned int i=0; i<genparticles->size(); ++i) {
438	>
439	>	GenParticle genp = genparticles->at(i);
440	>
441	>	//only take status 3 particles into account
442	>	if(genp.status()!=3) continue;
443	>
444	>	if(jet->deltaR(genp)<0.5)
445	>	matched_genparticle_ids.push_back(genp.pdgId());
446	>	}
447	>
448	>	//search for b quarks first
449	>	for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
450	>	if(abs(matched_genparticle_ids[i])==5) return matched_genparticle_ids[i];
451	>	}
452	>	//no b quark -> search for c quarks
453	>	for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
454	>	if(abs(matched_genparticle_ids[i])==4) return matched_genparticle_ids[i];
455	>	}
456	>	for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
457	>	if(abs(matched_genparticle_ids[i])==3) return matched_genparticle_ids[i];
458	>	}
459	>	for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
460	>	if(abs(matched_genparticle_ids[i])==2) return matched_genparticle_ids[i];
461	>	}
462	>	for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
463	>	if(abs(matched_genparticle_ids[i])==1) return matched_genparticle_ids[i];
464	>	}
465	>	for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
466	>	if(abs(matched_genparticle_ids[i])==21) return matched_genparticle_ids[i];
467	>	}
468	>
469	>	}
470	>
471	>	//no matched GenParticle -> return default value 0
472	>	return 0;
473	>
474		}

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines