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.7 by peiffer, Fri Dec 7 10:56:32 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	<	//variable HEP Tagger from Rebekka
4	>	// global function to define a tagged jet
5
6	<	bool variableHepTopTag(TopJet topjet, double ptJetMin, double massWindowLower, double massWindowUpper, double cutCondition2, double cutCondition3){
7	<
8	<	double mjet;
9	<	double ptjet;
10	<	int nsubjets;
11	<
12	<	double topmass=172.3;
13	<	double wmass=80.4;
14	<
15	<	nsubjets=topjet.numberOfDaughters();
16	<
17	<	LorentzVector allsubjets(0,0,0,0);
18	<
19	<	for(int j=0; j<topjet.numberOfDaughters(); ++j){
20	<	allsubjets += topjet.subjets()[j].v4();
21	<	}
22	<	if(!allsubjets.isTimelike()){
23	<	mjet=0;
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	<	}
16	>	}
17
18	<	mjet = allsubjets.M();
28	<	ptjet= allsubjets.Pt();
18	>	//variable HEP Tagger from Rebekka
19
20	<	double m12, m13, m23;
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	<	//The subjetcs have to be three
27	<	if(nsubjets==3) {
34	<
35	<	std::vector<Particle> subjets = topjet.subjets();
36	<	sort(subjets.begin(), subjets.end(), HigherPt());
37	<
38	<	m12 = 0;
39	<	if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
40	<	m12=(subjets[0].v4()+subjets[1].v4()).M();
41	<	m13 = 0;
42	<	if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
43	<	m13=(subjets[0].v4()+subjets[2].v4()).M();
44	<	m23 = 0;
45	<	if( (subjets[1].v4()+subjets[2].v4()).isTimelike()  )
46	<	m23 = (subjets[1].v4()+subjets[2].v4()).M();
47	<
48	<	}
49	<	else{
50	<	return false;
51	<	}
26	>	double topmass=172.3;
27	>	double wmass=80.4;
28
29	<	double rmin=massWindowLower*wmass/topmass;
54	<	double rmax=massWindowUpper*wmass/topmass;
29	>	nsubjets=topjet.numberOfDaughters();
30
31	<	int keep=0;
31	>	LorentzVector allsubjets(0,0,0,0);
32
33	<	//Conditions on the subjects: at least one has to be true
34	<	//1 condition
35	<	if(atan(m13/m12)>0.2 && atan(m13/m12)<1.3 && m23/mjet>rmin && m23/mjet<rmax) keep=1;
36	<
37	<	//2 condition
38	<	double cond2left=pow(rmin,2)*(1+pow((m13/m12),2));
39	<	double cond2cent=1-pow(m23/mjet,2);
65	<	double cond2right=pow(rmax,2)*(1+pow(m13/m12,2));
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	<	if(cond2left<cond2cent && cond2cent<cond2right && m23/mjet>cutCondition2) keep=1;
41	>	mjet = allsubjets.M();
42	>	ptjet= allsubjets.Pt();
43
44	<	//3 condition
70	<	double cond3left=pow(rmin,2)*(1+pow((m12/m13),2));
71	<	double cond3cent=1-pow(m23/mjet,2);
72	<	double cond3right=pow(rmax,2)*(1+pow(m12/m13,2));
44	>	double m12, m13, m23;
45
46	<	if(cond3left<cond3cent && cond3cent<cond3right && m23/mjet>cutCondition3) keep=1;
46	>	//The subjetcs have to be three
47	>	if(nsubjets==3) {
48
49	<	//Final requirement: at least one of the three subjets conditions and total pt
50	<	if(keep==1 && ptjet>ptJetMin){
51	<	return true;
52	<	}
53	<	else{
54	<	return false;
55	<	}
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
#	Line 90 | Line 102 | bool variableHepTopTag(TopJet topjet, do
102
103		//HEP Tagger from Ivan
104
105	<	bool HepTopTag(TopJet topjet){
106	<
95	<	double mjet;
96	<	double ptjet;
97	<	int nsubjets;
98	<
99	<	double topmass=172.3;
100	<	double wmass=80.4;
101	<
102	<	nsubjets=topjet.numberOfDaughters();
103	<
104	<	LorentzVector allsubjets(0,0,0,0);
105	<
106	<	for(int j=0; j<topjet.numberOfDaughters(); ++j){
107	<	allsubjets += topjet.subjets()[j].v4();
108	<	}
109	<	if(!allsubjets.isTimelike()){
110	<	mjet=0;
111	<	return false;
112	<	}
105	>	bool HepTopTag(TopJet topjet)
106	>	{
107
108	<	mjet = allsubjets.M();
109	<	ptjet= allsubjets.Pt();
108	>	double mjet;
109	>	double ptjet;
110	>	int nsubjets;
111
112	<	double m12, m13, m23;
112	>	double topmass=172.3;
113	>	double wmass=80.4;
114
115	<	//The subjetcs have to be three
120	<	if(nsubjets==3) {
121	<
122	<	std::vector<Particle> subjets = topjet.subjets();
123	<	sort(subjets.begin(), subjets.end(), HigherPt());
124	<
125	<	m12 = 0;
126	<	if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
127	<	m12=(subjets[0].v4()+subjets[1].v4()).M();
128	<	m13 = 0;
129	<	if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
130	<	m13=(subjets[0].v4()+subjets[2].v4()).M();
131	<	m23 = 0;
132	<	if( (subjets[1].v4()+subjets[2].v4()).isTimelike()  )
133	<	m23 = (subjets[1].v4()+subjets[2].v4()).M();
134	<
135	<	}
136	<	else{
137	<	return false;
138	<	}
115	>	nsubjets=topjet.numberOfDaughters();
116
117	<	double rmin=0.85*wmass/topmass;
141	<	double rmax=1.15*wmass/topmass;
117	>	LorentzVector allsubjets(0,0,0,0);
118
119	<	int keep=0;
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	<	//Conditions on the subjects: at least one has to be true
128	<	//1 condition
147	<	if(atan(m13/m12)>0.2 && atan(m13/m12)<1.3 && m23/mjet>rmin && m23/mjet<rmax) keep=1;
148	<
149	<	//2 condition
150	<	double cond2left=pow(rmin,2)*(1+pow((m13/m12),2));
151	<	double cond2cent=1-pow(m23/mjet,2);
152	<	double cond2right=pow(rmax,2)*(1+pow(m13/m12,2));
127	>	mjet = allsubjets.M();
128	>	ptjet= allsubjets.Pt();
129
130	<	if(cond2left<cond2cent && cond2cent<cond2right && m23/mjet>0.35) keep=1;
130	>	double m12, m13, m23;
131
132	<	//3 condition
133	<	double cond3left=pow(rmin,2)*(1+pow((m12/m13),2));
158	<	double cond3cent=1-pow(m23/mjet,2);
159	<	double cond3right=pow(rmax,2)*(1+pow(m12/m13,2));
132	>	//The subjetcs have to be three
133	>	if(nsubjets==3) {
134
135	<	if(cond3left<cond3cent && cond3cent<cond3right && m23/mjet>0.35) keep=1;
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	<	//Final requirement: at least one of the three subjets conditions and total pt
149	<	if(keep==1 && ptjet>200){
150	<	return true;
151	<	}
152	<	else{
153	<	return false;
154	<	}
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	<	nsubjets=topjet.numberOfDaughters();
189	<
190	<	LorentzVector allsubjets(0,0,0,0);
191	<
192	<	for(int j=0; j<topjet.numberOfDaughters(); ++j){
193	<	allsubjets += topjet.subjets()[j].v4();
194	<	}
195	<	if(!allsubjets.isTimelike()){
196	<	mjet=0;
197	<	mmin=0;
198	<	//  mminLower=50;
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	<	}
192	<
193	<	mjet = allsubjets.M();
194	<
195	<	if(nsubjets>=3) {
196	<
197	<	std::vector<Particle> subjets = topjet.subjets();
198	<	sort(subjets.begin(), subjets.end(), HigherPt());
199	<
200	<	double m01 = 0;
201	<	if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
202	<	m01=(subjets[0].v4()+subjets[1].v4()).M();
203	<	double m02 = 0;
204	<	if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
205	<	m02=(subjets[0].v4()+subjets[2].v4()).M();
206	<	double m12 = 0;
207	<	if( (subjets[1].v4()+subjets[2].v4()).isTimelike() )
208	<	m12 = (subjets[1].v4()+subjets[2].v4()).M();
209	<
210	<	//minimum pairwise mass
211	<	mmin = std::min(m01,std::min(m02,m12));
212	<	}
213	<
214	<	//at least 3 sub-jets
215	<	if(nsubjets<3) return false;
216	<	//minimum pairwise mass > 50 GeV/c^2
217	<	if(mmin<mminLower) return false;
218	<	//jet mass between 140 and 250 GeV/c^2
219	<	if(mjet<mjetLower || mjet>mjetUpper) return false;
220	<
221	<	return true;
222	<	}
223	<
224	<
225	<	bool TopTag(TopJet topjet,  double &mjet, int &nsubjets, double &mmin){
226	<
227	<	nsubjets=topjet.numberOfDaughters();
228	<
229	<	LorentzVector allsubjets(0,0,0,0);
230	<
231	<	for(int j=0; j<topjet.numberOfDaughters(); ++j){
232	<	allsubjets += topjet.subjets()[j].v4();
233	<	}
234	<	if(!allsubjets.isTimelike()){
235	<	mjet=0;
236	<	mmin=0;
237	<	return false;
238	<	}
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());
246	<
247	<	double m01 = 0;
248	<	if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
249	<	m01=(subjets[0].v4()+subjets[1].v4()).M();
250	<	double m02 = 0;
251	<	if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
252	<	m02=(subjets[0].v4()+subjets[2].v4()).M();
253	<	double m12 = 0;
254	<	if( (subjets[1].v4()+subjets[2].v4()).isTimelike()  )
255	<	m12 = (subjets[1].v4()+subjets[2].v4()).M();
256	<
257	<	//minimum pairwise mass
258	<	mmin = std::min(m01,std::min(m02,m12));
259	<	}
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	<	//at least 3 sub-jets
213	<	if(nsubjets<3) return false;
214	<	//minimum pairwise mass > 50 GeV/c^2
215	<	if(mmin<50) return false;
216	<	//jet mass between 140 and 250 GeV/c^2
217	<	if(mjet<140 || mjet>250) return false;
218	<
219	<	return true;
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	>	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	>	return true;
234		}
235
271	–	Jet* nextJet(const Particle *p, std::vector<Jet> *jets){
236
237	<	double deltarmin = double_infinity();
238	<	Jet* nextjet=0;
239	<	for(unsigned int i=0; i<jets->size();++i){
240	<	if(jets->at(i).deltaR(*p) < deltarmin){
241	<	deltarmin = jets->at(i).deltaR(*p);
242	<	nextjet = &jets->at(i);
237	>	bool TopTag(TopJet topjet,  double &mjet, int &nsubjets, double &mmin)
238	>	{
239	>
240	>	nsubjets=topjet.numberOfDaughters();
241	>
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	>	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		}
280	–	}
273
274	<	return nextjet;
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	<	bool WTag(TopJet prunedjet,  double& mjet, int &nsubjets, double& massdrop){
284	>	Jet* nextJet(const Particle *p, std::vector<Jet> *jets)
285	>	{
286
287	<	nsubjets=prunedjet.numberOfDaughters();
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	>	}
295
296	<	mjet = 0;
297	<	if(prunedjet.v4().isTimelike())
298	<	mjet = prunedjet.v4().M();
296	>	return nextjet;
297	>	}
298	>
299	>	bool WTag(TopJet prunedjet,  double& mjet, int &nsubjets, double& massdrop)
300	>	{
301	>
302	>	nsubjets=prunedjet.numberOfDaughters();
303
304	<	//calculate mass drop for first sub-jet ordered in pt
305	<	massdrop = 0;
306	<	if(nsubjets>=1 && mjet>0){
304	>	mjet = 0;
305	>	if(prunedjet.v4().isTimelike())
306	>	mjet = prunedjet.v4().M();
307
308	<	std::vector< Particle > subjets = prunedjet.subjets();
309	<	sort(subjets.begin(), subjets.end(), HigherPt());
308	>	//calculate mass drop for first sub-jet ordered in pt
309	>	massdrop = 0;
310	>	if(nsubjets>=1 && mjet>0) {
311
312	<	double m1 = 0;
313	<	if(subjets[0].v4().isTimelike())
302	<	m1 = subjets[0].v4().M();
312	>	std::vector< Particle > subjets = prunedjet.subjets();
313	>	sort(subjets.begin(), subjets.end(), HigherPt());
314
315	<	massdrop = m1/mjet;
316	<	}
317	<
318	<	//at least 2 sub-jets
319	<	if(nsubjets<2) return false;
320	<	//60 GeV < pruned jet mass < 100 GeV
310	<	if(mjet <= 60 || mjet >= 100) return false;
311	<	//mass drop < 0.4
312	<	if(massdrop>=0.4) return false;
315	>	double m1 = 0;
316	>	if(subjets[0].v4().isTimelike())
317	>	m1 = subjets[0].v4().M();
318	>
319	>	massdrop = m1/mjet;
320	>	}
321
322	<	return true;
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	>	return true;
330
331		}
332
333	<	double pTrel(const Particle *p, std::vector<Jet> *jets){
333	>	double pTrel(const Particle *p, std::vector<Jet> *jets)
334	>	{
335
336	<	double ptrel=0;
336	>	double ptrel=0;
337
338	<	Jet* nextjet =  nextJet(p,jets);
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());
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	<	}
347	<	else{
348	<	std::cout << "something strange happend in the ptrel calculation: either lepton or jet momentum is 0" <<std::endl;
333	<	}
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;
350	>	return ptrel;
351		}
352
353	<	double deltaRmin(const Particle *p, std::vector<Jet> *jets){
354	<	return nextJet(p,jets)->deltaR(*p);
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){
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;
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){
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;
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){
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;
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){
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;
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){
400	>	double deltaR(LorentzVector v1, LorentzVector v2)
401	>	{
402
403	<	Particle p1;
404	<	p1.set_v4(v1);
405	<	Particle p2;
406	<	p2.set_v4(v2);
407	<	return p1.deltaR(p2);
403	>	Particle p1;
404	>	p1.set_v4(v1);
405	>	Particle p2;
406	>	p2.set_v4(v2);
407	>	return p1.deltaR(p2);
408		}
409
410	<	double double_infinity(){
411	<	return std::numeric_limits<double>::infinity() ;
410	>	double double_infinity()
411	>	{
412	>	return std::numeric_limits<double>::infinity() ;
413		}
414
415	<	int int_infinity(){
416	<	return std::numeric_limits<int>::max() ;
415	>	int int_infinity()
416	>	{
417	>	return std::numeric_limits<int>::max() ;
418		}
419
420	<	int myPow(int x, unsigned int p) {
421	<	int i = 1;
422	<	for (unsigned int j = 1; j <= p; j++)  i *= x;
423	<	return i;
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){
427	>	int JetFlavor(Jet *jet)
428	>	{
429
430	<	EventCalc* calc = EventCalc::Instance();
430	>	EventCalc* calc = EventCalc::Instance();
431
432	<	std::vector< GenParticle >* genparticles = calc->GetGenParticles();
433	<	if(genparticles){
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){
413	<
414	<	GenParticle genp = genparticles->at(i);
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	<	//only take status 3 particles into account
417	<	if(genp.status()!=3) continue;
439	>	GenParticle genp = genparticles->at(i);
440
441	<	if(jet->deltaR(genp)<0.5)
442	<	matched_genparticle_ids.push_back(genp.pdgId());
443	<	}
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
423	–	//search for b quarks first
424	–	for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i){
425	–	if(abs(matched_genparticle_ids[i])==5) return matched_genparticle_ids[i];
426	–	}
427	–	//no b quark -> search for c quarks
428	–	for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i){
429	–	if(abs(matched_genparticle_ids[i])==4) return matched_genparticle_ids[i];
430	–	}
431	–	for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i){
432	–	if(abs(matched_genparticle_ids[i])==3) return matched_genparticle_ids[i];
433	–	}
434	–	for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i){
435	–	if(abs(matched_genparticle_ids[i])==2) return matched_genparticle_ids[i];
436	–	}
437	–	for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i){
438	–	if(abs(matched_genparticle_ids[i])==1) return matched_genparticle_ids[i];
439	–	}
440	–	for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i){
441	–	if(abs(matched_genparticle_ids[i])==21) return matched_genparticle_ids[i];
469		}
443	–
444	–	}
470
471	<	//no matched GenParticle -> return default value 0
472	<	return 0;
471	>	//no matched GenParticle -> return default value 0
472	>	return 0;
473
474		}

Diff Legend

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