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.6 by hoeing, Mon Oct 15 15:35:51 2012 UTC vs.
Revision 1.12 by hoeing, Tue Mar 26 13:15:59 2013 UTC

#	Line 1 | Line 1
1		#include "../include/Utils.h"
2
3	+	#include <fastjet/JetDefinition.hh>
4	+	#include <fastjet/PseudoJet.hh>
5	+	#include <fastjet/ClusterSequence.hh>
6	+	#include <fastjet/ClusterSequenceArea.hh>
7	+	#include <fastjet/GhostedAreaSpec.hh>
8	+	namespace external {
9	+	#include "../include/HEPTopTagger.h"
10	+	}
11
4	–	//variable HEP Tagger from Rebekka
12
13	<	bool variableHepTopTag(TopJet topjet, double ptJetMin, double massWindowLower, double massWindowUpper, double cutCondition2, double cutCondition3){
14	<
15	<	double mjet;
16	<	double ptjet;
17	<	int nsubjets;
18	<
19	<	double topmass=172.3;
20	<	double wmass=80.4;
21	<
22	<	nsubjets=topjet.numberOfDaughters();
23	<
24	<	LorentzVector allsubjets(0,0,0,0);
25	<
26	<	for(int j=0; j<topjet.numberOfDaughters(); ++j){
27	<	allsubjets += topjet.subjets()[j].v4();
13	>
14	>	//subjet b-tagger, returns number of b-tagged subjets
15	>
16	>	int subJetBTag(TopJet topjet, E_BtagType type){
17	>	int nBTagsSub = 0;
18	>	float discriminator_cut;
19	>	if(type==e_CSVL) discriminator_cut = 0.244;
20	>	if(type==e_CSVM) discriminator_cut = 0.679;
21	>	if(type==e_CSVT) discriminator_cut = 0.898;
22	>
23	>	//Create a vector of subjets
24	>	std::vector<Particle> subjets_top;
25	>	//Create a float vector of the subjets discriminators
26	>	std::vector<float> btagsub_combinedSecondaryVertex_top;
27	>
28	>	//Fill the vector of subjets with the subjets of a topjet
29	>	subjets_top=topjet.subjets();
30	>	//Fill the vector of discriminators with the discriminators of the subjets of a certain topjet
31	>	btagsub_combinedSecondaryVertex_top=topjet.btagsub_combinedSecondaryVertex();
32	>
33	>	//Looping over subjets and checking if they are b-tagged
34	>	for(unsigned int i=0; i < btagsub_combinedSecondaryVertex_top.size(); ++i){
35	>	float test=btagsub_combinedSecondaryVertex_top[i];
36	>	if(test>discriminator_cut){
37	>	nBTagsSub += 1;
38	>	//This means it is b-tagged
39	>	}
40		}
41	<	if(!allsubjets.isTimelike()){
42	<	mjet=0;
43	<	return false;
41	>	return nBTagsSub;
42	>	}
43	>
44	>
45	>	bool HiggsTag(TopJet topjet, E_BtagType type1, E_BtagType type2){
46	>	int nBTagsSub1 = 0;
47	>	int nBTagsSub2 = 0;
48	>	float discriminator_cut1;
49	>	float discriminator_cut2;
50	>	if(type1==e_CSVL) discriminator_cut1 = 0.244;
51	>	if(type1==e_CSVM) discriminator_cut1 = 0.679;
52	>	if(type1==e_CSVT) discriminator_cut1 = 0.898;
53	>	if(type2==e_CSVL) discriminator_cut2 = 0.244;
54	>	if(type2==e_CSVM) discriminator_cut2 = 0.679;
55	>	if(type2==e_CSVT) discriminator_cut2 = 0.898;
56	>	// std::cout << "discriminator_cut1: " <<  discriminator_cut1 << " discriminator_cut2: "<< discriminator_cut1 << std::endl;
57	>	//Create a vector of subjets
58	>	std::vector<Particle> subjets_top;
59	>	//Create a float vector of the subjets discriminators
60	>	std::vector<float> btagsub_combinedSecondaryVertex_top;
61	>
62	>	//Fill the vector of subjets with the subjets of a topjet
63	>	subjets_top=topjet.subjets();
64	>	//Fill the vector of discriminators with the discriminators of the subjets of a certain topjet
65	>	btagsub_combinedSecondaryVertex_top=topjet.btagsub_combinedSecondaryVertex();
66	>
67	>	//Looping over subjets and checking if they are b-tagged
68	>	for(unsigned int i=0; i < btagsub_combinedSecondaryVertex_top.size(); ++i){
69	>	float test=btagsub_combinedSecondaryVertex_top[i];
70	>	if (nBTagsSub1 != 0 && test>discriminator_cut2) nBTagsSub2 =+ 1;
71	>	if(test>discriminator_cut1){
72	>	if(test>discriminator_cut2 && nBTagsSub2==0) nBTagsSub2+=1;
73	>	else nBTagsSub1 += 1;      //This means it is b-tagged
74	>
75	>	}
76		}
77	+	if (nBTagsSub1!=0 && nBTagsSub2!=0) return true;
78	+	else return false;
79	+	}
80
27	–	mjet = allsubjets.M();
28	–	ptjet= allsubjets.Pt();
81
82	<	double m12, m13, m23;
82	>	bool HepTopTagFull(TopJet topjet){
83
84	<	//The subjetcs have to be three
85	<	if(nsubjets==3) {
84	>	//Transform the SFrame TopJet object in a fastjet::PseudoJet
85	>
86	>	std::vector<fastjet::PseudoJet> jetpart;
87	>	std::vector<Particle> pfconstituents_jet;
88	>
89	>	fastjet::ClusterSequence* JetFinder;
90	>	fastjet::JetDefinition* JetDef ;
91
92	<	std::vector<Particle> subjets = topjet.subjets();
93	<	sort(subjets.begin(), subjets.end(), HigherPt());
94	<
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();
92	>	pfconstituents_jet=topjet.pfconstituents();
93	>
94	>	for(unsigned int ic=0; ic<pfconstituents_jet.size(); ++ic){
95
96	+	jetpart.push_back( fastjet::PseudoJet(
97	+	pfconstituents_jet[ic].pt()*cos(pfconstituents_jet[ic].phi()),
98	+	pfconstituents_jet[ic].pt()*sin(pfconstituents_jet[ic].phi()),
99	+	pfconstituents_jet[ic].pt()*sinh(pfconstituents_jet[ic].eta()),
100	+	pfconstituents_jet[ic].energy() ) );
101	+
102		}
103	<	else{
103	>
104	>	//Clustering definition
105	>	double conesize=3;
106	>	JetDef = new
107	>	fastjet::JetDefinition(fastjet::cambridge_algorithm,conesize);
108	>
109	>	JetFinder = new fastjet::ClusterSequence(jetpart, *JetDef);
110	>
111	>	std::vector<fastjet::PseudoJet> tops = JetFinder->inclusive_jets(10.);
112	>
113	>	if (tops.size() != 1){
114	>	std::cout << "Problem! Doesn't give exactly one jet!!!!!!!!!!!!!!Gives " << tops.size() << " jets" << std::endl;
115	>	delete JetFinder;
116	>	delete JetDef;
117		return false;
118		}
119
120	<	double rmin=massWindowLower*wmass/topmass;
54	<	double rmax=massWindowUpper*wmass/topmass;
120	>	std::vector<fastjet::PseudoJet> SortedJets = sorted_by_pt(tops);
121
56	–	int keep=0;
122
123	<	//Conditions on the subjects: at least one has to be true
124	<	//1 condition
60	<	if(atan(m13/m12)>0.2 && atan(m13/m12)<1.3 && m23/mjet>rmin && m23/mjet<rmax) keep=1;
61	<
62	<	//2 condition
63	<	double cond2left=pow(rmin,2)*(1+pow((m13/m12),2));
64	<	double cond2cent=1-pow(m23/mjet,2);
65	<	double cond2right=pow(rmax,2)*(1+pow(m13/m12,2));
66	<
67	<	if(cond2left<cond2cent && cond2cent<cond2right && m23/mjet>cutCondition2) keep=1;
68	<
69	<	//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));
123	>	//Run the HEPTopTagger
124	>	external::HEPTopTagger tagger(*JetFinder, SortedJets[0]);
125
126	<	if(cond3left<cond3cent && cond3cent<cond3right && m23/mjet>cutCondition3) keep=1;
126	>	//Mass window to be applied in a second step
127	>	tagger.set_top_range(0.0, 10000.0);
128	>	tagger.set_mass_drop_threshold(0.8);
129	>	tagger.set_max_subjet_mass(30);
130
131	<	//Final requirement: at least one of the three subjets conditions and total pt
132	<	if(keep==1 && ptjet>ptJetMin){
133	<	return true;
134	<	}
135	<	else{
136	<	return false;
137	<	}
131	>	tagger.run_tagger();
132	>
133	>	delete JetFinder;
134	>	delete JetDef;
135	>
136	>	if (tagger.is_masscut_passed()) return true;
137	>	else return false;
138
139	+	return true;
140	+
141		}
142
143
144
145	+	// global function to define a tagged jet
146
147	+	bool IsTagged(Jet & jet, E_BtagType type)
148	+	{
149	+	if(type==e_CSVL && jet.btag_combinedSecondaryVertex()>0.244) return true;
150	+	if(type==e_CSVM && jet.btag_combinedSecondaryVertex()>0.679) return true;
151	+	if(type==e_CSVT && jet.btag_combinedSecondaryVertex()>0.898) return true;
152	+	if(type==e_JPL && jet.btag_jetProbability()>0.275) return true;
153	+	if(type==e_JPM && jet.btag_jetProbability()>0.545) return true;
154	+	if(type==e_JPT && jet.btag_jetProbability()>0.790) return true;
155	+
156	+	return false;
157	+	}
158
159	+	//variable HEP Tagger from Rebekka
160
161	<	//HEP Tagger from Ivan
161	>	bool variableHepTopTag(TopJet topjet, double ptJetMin, double massWindowLower, double massWindowUpper, double cutCondition2, double cutCondition3)
162	>	{
163	>	double mjet;
164	>	double ptjet;
165	>	int nsubjets;
166
167	<	bool HepTopTag(TopJet topjet){
168	<
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	<	}
167	>	double topmass=172.3;
168	>	double wmass=80.4;
169
170	<	mjet = allsubjets.M();
115	<	ptjet= allsubjets.Pt();
170	>	nsubjets=topjet.numberOfDaughters();
171
172	<	double m12, m13, m23;
172	>	LorentzVector allsubjets(0,0,0,0);
173
174	<	//The subjetcs have to be three
175	<	if(nsubjets==3) {
174	>	for(int j=0; j<topjet.numberOfDaughters(); ++j) {
175	>	allsubjets += topjet.subjets()[j].v4();
176	>	}
177	>	if(!allsubjets.isTimelike()) {
178	>	mjet=0;
179	>	return false;
180	>	}
181
182	<	std::vector<Particle> subjets = topjet.subjets();
183	<	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	<	}
182	>	mjet = allsubjets.M();
183	>	ptjet= allsubjets.Pt();
184
185	<	double rmin=0.85*wmass/topmass;
141	<	double rmax=1.15*wmass/topmass;
185	>	double m12, m13, m23;
186
187	<	int keep=0;
187	>	//The subjetcs have to be three
188	>	if(nsubjets==3) {
189
190	<	//Conditions on the subjects: at least one has to be true
191	<	//1 condition
192	<	if(atan(m13/m12)>0.2 && atan(m13/m12)<1.3 && m23/mjet>rmin && m23/mjet<rmax) keep=1;
193	<
194	<	//2 condition
195	<	double cond2left=pow(rmin,2)*(1+pow((m13/m12),2));
196	<	double cond2cent=1-pow(m23/mjet,2);
197	<	double cond2right=pow(rmax,2)*(1+pow(m13/m12,2));
198	<
199	<	if(cond2left<cond2cent && cond2cent<cond2right && m23/mjet>0.35) keep=1;
200	<
201	<	//3 condition
157	<	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));
190	>	std::vector<Particle> subjets = topjet.subjets();
191	>	sort(subjets.begin(), subjets.end(), HigherPt());
192	>
193	>	m12 = 0;
194	>	if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
195	>	m12=(subjets[0].v4()+subjets[1].v4()).M();
196	>	m13 = 0;
197	>	if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
198	>	m13=(subjets[0].v4()+subjets[2].v4()).M();
199	>	m23 = 0;
200	>	if( (subjets[1].v4()+subjets[2].v4()).isTimelike()  )
201	>	m23 = (subjets[1].v4()+subjets[2].v4()).M();
202
203	<	if(cond3left<cond3cent && cond3cent<cond3right && m23/mjet>0.35) keep=1;
203	>	} else {
204	>	return false;
205	>	}
206
207	<	//Final requirement: at least one of the three subjets conditions and total pt
208	<	if(keep==1 && ptjet>200){
209	<	return true;
210	<	}
211	<	else{
212	<	return false;
213	<	}
207	>	double rmin=massWindowLower*wmass/topmass;
208	>	double rmax=massWindowUpper*wmass/topmass;
209	>
210	>	int keep=0;
211	>
212	>	//Conditions on the subjects: at least one has to be true
213	>	//1 condition
214	>	if(atan(m13/m12)>0.2 && atan(m13/m12)<1.3 && m23/mjet>rmin && m23/mjet<rmax) keep=1;
215	>
216	>	//2 condition
217	>	double cond2left=pow(rmin,2)*(1+pow((m13/m12),2));
218	>	double cond2cent=1-pow(m23/mjet,2);
219	>	double cond2right=pow(rmax,2)*(1+pow(m13/m12,2));
220	>
221	>	if(cond2left<cond2cent && cond2cent<cond2right && m23/mjet>cutCondition2) keep=1;
222	>
223	>	//3 condition
224	>	double cond3left=pow(rmin,2)*(1+pow((m12/m13),2));
225	>	double cond3cent=1-pow(m23/mjet,2);
226	>	double cond3right=pow(rmax,2)*(1+pow(m12/m13,2));
227	>
228	>	if(cond3left<cond3cent && cond3cent<cond3right && m23/mjet>cutCondition3) keep=1;
229	>
230	>	//Final requirement: at least one of the three subjets conditions and total pt
231	>	if(keep==1 && ptjet>ptJetMin) {
232	>	return true;
233	>	} else {
234	>	return false;
235	>	}
236
237		}
238
239
240	+	//HEP Tagger from Ivan
241	+
242	+	bool HepTopTag(TopJet topjet)
243	+	{
244	+	//call variable tagger with default parameters
245	+	return variableHepTopTag(topjet);
246	+
247	+	}
248	+
249		//default values (mminLower=50., mjetLower=140, mjetUpper=250.) defined in Utils.h
250	<	bool variableTopTag(TopJet topjet, double &mjet, int &nsubjets, double &mmin, double mminLower, double mjetLower, double mjetUpper){
251	<
252	<	nsubjets=topjet.numberOfDaughters();
253	<
254	<	LorentzVector allsubjets(0,0,0,0);
255	<
256	<	for(int j=0; j<topjet.numberOfDaughters(); ++j){
257	<	allsubjets += topjet.subjets()[j].v4();
258	<	}
259	<	if(!allsubjets.isTimelike()){
260	<	mjet=0;
261	<	mmin=0;
262	<	//  mminLower=50;
250	>	bool variableTopTag(TopJet topjet, double &mjet, int &nsubjets, double &mmin, double mminLower, double mjetLower, double mjetUpper)
251	>	{
252	>
253	>	nsubjets=topjet.numberOfDaughters();
254	>
255	>	LorentzVector allsubjets(0,0,0,0);
256	>
257	>	for(int j=0; j<topjet.numberOfDaughters(); ++j) {
258	>	allsubjets += topjet.subjets()[j].v4();
259	>	}
260	>	if(!allsubjets.isTimelike()) {
261	>	mjet=0;
262	>	mmin=0;
263	>	//  mminLower=50;
264		//   mjetLower=140;
265		//   mjetUpper=250;
266	<	return false;
267	<	}
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	<	}
266	>	return false;
267	>	}
268
269	+	mjet = allsubjets.M();
270
271	<	bool TopTag(TopJet topjet,  double &mjet, int &nsubjets, double &mmin){
226	<
227	<	nsubjets=topjet.numberOfDaughters();
271	>	if(nsubjets>=3) {
272
273	<	LorentzVector allsubjets(0,0,0,0);
274	<
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	<	}
273	>	std::vector<Particle> subjets = topjet.subjets();
274	>	sort(subjets.begin(), subjets.end(), HigherPt());
275
276	<	mjet = allsubjets.M();
277	<
278	<	if(nsubjets>=3) {
279	<
280	<	std::vector<Particle> subjets = topjet.subjets();
281	<	sort(subjets.begin(), subjets.end(), HigherPt());
282	<
283	<	double m01 = 0;
284	<	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	<	}
276	>	double m01 = 0;
277	>	if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
278	>	m01=(subjets[0].v4()+subjets[1].v4()).M();
279	>	double m02 = 0;
280	>	if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
281	>	m02=(subjets[0].v4()+subjets[2].v4()).M();
282	>	double m12 = 0;
283	>	if( (subjets[1].v4()+subjets[2].v4()).isTimelike() )
284	>	m12 = (subjets[1].v4()+subjets[2].v4()).M();
285
286	<	//at least 3 sub-jets
287	<	if(nsubjets<3) return false;
288	<	//minimum pairwise mass > 50 GeV/c^2
289	<	if(mmin<50) return false;
290	<	//jet mass between 140 and 250 GeV/c^2
291	<	if(mjet<140 || mjet>250) return false;
292	<
293	<	return true;
286	>	//minimum pairwise mass
287	>	mmin = std::min(m01,std::min(m02,m12));
288	>	}
289	>
290	>	//at least 3 sub-jets
291	>	if(nsubjets<3) return false;
292	>	//minimum pairwise mass > 50 GeV/c^2
293	>	if(mmin<mminLower) return false;
294	>	//jet mass between 140 and 250 GeV/c^2
295	>	if(mjet<mjetLower || mjet>mjetUpper) return false;
296	>
297	>	return true;
298		}
299
271	–	Jet* nextJet(const Particle *p, std::vector<Jet> *jets){
300
301	<	double deltarmin = double_infinity();
302	<	Jet* nextjet=0;
303	<	for(unsigned int i=0; i<jets->size();++i){
304	<	if(jets->at(i).deltaR(*p) < deltarmin){
277	<	deltarmin = jets->at(i).deltaR(*p);
278	<	nextjet = &jets->at(i);
279	<	}
280	<	}
301	>	bool TopTag(TopJet topjet,  double &mjet, int &nsubjets, double &mmin)
302	>	{
303	>	//call variable tagger with default parameters
304	>	return variableTopTag(topjet, mjet, nsubjets, mmin);
305
282	–	return nextjet;
306		}
307	+
308	+	Jet* nextJet(const Particle *p, std::vector<Jet> *jets)
309	+	{
310	+
311	+	double deltarmin = double_infinity();
312	+	Jet* nextjet=0;
313	+	for(unsigned int i=0; i<jets->size(); ++i) {
314	+	Jet ji = jets->at(i);
315	+	if (fabs(p->pt() - ji.pt())<1e-8) continue; // skip identical particle
316	+	if(jets->at(i).deltaR(*p) < deltarmin) {
317	+	deltarmin = jets->at(i).deltaR(*p);
318	+	nextjet = &jets->at(i);
319	+	}
320	+	}
321
322	<	bool WTag(TopJet prunedjet,  double& mjet, int &nsubjets, double& massdrop){
323	<
287	<	nsubjets=prunedjet.numberOfDaughters();
322	>	return nextjet;
323	>	}
324
325	<	mjet = 0;
326	<	if(prunedjet.v4().isTimelike())
291	<	mjet = prunedjet.v4().M();
325	>	bool WTag(TopJet prunedjet,  double& mjet, int &nsubjets, double& massdrop)
326	>	{
327
328	<	//calculate mass drop for first sub-jet ordered in pt
294	<	massdrop = 0;
295	<	if(nsubjets>=1 && mjet>0){
328	>	nsubjets=prunedjet.numberOfDaughters();
329
330	<	std::vector< Particle > subjets = prunedjet.subjets();
331	<	sort(subjets.begin(), subjets.end(), HigherPt());
330	>	mjet = 0;
331	>	if(prunedjet.v4().isTimelike())
332	>	mjet = prunedjet.v4().M();
333
334	<	double m1 = 0;
335	<	if(subjets[0].v4().isTimelike())
336	<	m1 = subjets[0].v4().M();
334	>	//calculate mass drop for first sub-jet ordered in pt
335	>	massdrop = 0;
336	>	if(nsubjets>=1 && mjet>0) {
337
338	<	massdrop = m1/mjet;
339	<	}
306	<
307	<	//at least 2 sub-jets
308	<	if(nsubjets<2) return false;
309	<	//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;
338	>	std::vector< Particle > subjets = prunedjet.subjets();
339	>	sort(subjets.begin(), subjets.end(), HigherPt());
340
341	<	return true;
341	>	double m1 = 0;
342	>	if(subjets[0].v4().isTimelike())
343	>	m1 = subjets[0].v4().M();
344
345	<	}
345	>	massdrop = m1/mjet;
346	>	}
347
348	<	double pTrel(const Particle *p, std::vector<Jet> *jets){
348	>	//at least 2 sub-jets
349	>	if(nsubjets<2) return false;
350	>	//60 GeV < pruned jet mass < 100 GeV
351	>	if(mjet <= 60 || mjet >= 100) return false;
352	>	//mass drop < 0.4
353	>	if(massdrop>=0.4) return false;
354
355	<	double ptrel=0;
355	>	return true;
356
357	<	Jet* nextjet =  nextJet(p,jets);
357	>	}
358
359	<	TVector3 p3(p->v4().Px(),p->v4().Py(),p->v4().Pz());
360	<	TVector3 jet3(nextjet->v4().Px(),nextjet->v4().Py(),nextjet->v4().Pz());
359	>	double pTrel(const Particle *p, std::vector<Jet> *jets)
360	>	{
361
362	<	if(p3.Mag()!=0 && jet3.Mag()!=0){
363	<	double sin_alpha = (p3.Cross(jet3)).Mag()/p3.Mag()/jet3.Mag();
364	<	ptrel = p3.Mag()*sin_alpha;
365	<	}
366	<	else{
367	<	std::cout << "something strange happend in the ptrel calculation: either lepton or jet momentum is 0" <<std::endl;
368	<	}
362	>	double ptrel=0;
363	>	Jet* nextjet =  nextJet(p,jets);
364	>	if (!nextjet) return ptrel;
365	>
366	>	TVector3 p3(p->v4().Px(),p->v4().Py(),p->v4().Pz());
367	>	TVector3 jet3(nextjet->v4().Px(),nextjet->v4().Py(),nextjet->v4().Pz());
368	>
369	>	if(p3.Mag()!=0 && jet3.Mag()!=0) {
370	>	double sin_alpha = (p3.Cross(jet3)).Mag()/p3.Mag()/jet3.Mag();
371	>	ptrel = p3.Mag()*sin_alpha;
372	>	} else {
373	>	std::cout << "something strange happend in the ptrel calculation: either lepton or jet momentum is 0" <<std::endl;
374	>	}
375
376	<	return ptrel;
376	>	return ptrel;
377		}
378
379	<	double deltaRmin(const Particle *p, std::vector<Jet> *jets){
380	<	return nextJet(p,jets)->deltaR(*p);
379	>	double deltaRmin(const Particle *p, std::vector<Jet> *jets)
380	>	{
381	>	Jet* j = nextJet(p,jets);
382	>	double dr = 999.;
383	>	if (j) dr = j->deltaR(*p);
384	>	return dr;
385		}
386
387	<	TVector3 toVector(LorentzVector v4){
387	>	TVector3 toVector(LorentzVector v4)
388	>	{
389
390	<	TVector3 v3(0,0,0);
391	<	v3.SetX(v4.X());
392	<	v3.SetY(v4.Y());
393	<	v3.SetZ(v4.Z());
394	<	return v3;
390	>	TVector3 v3(0,0,0);
391	>	v3.SetX(v4.X());
392	>	v3.SetY(v4.Y());
393	>	v3.SetZ(v4.Z());
394	>	return v3;
395		}
396
397	<	TVector3 toVector(LorentzVectorXYZE v4){
397	>	TVector3 toVector(LorentzVectorXYZE v4)
398	>	{
399
400	<	TVector3 v3(0,0,0);
401	<	v3.SetX(v4.X());
402	<	v3.SetY(v4.Y());
403	<	v3.SetZ(v4.Z());
404	<	return v3;
400	>	TVector3 v3(0,0,0);
401	>	v3.SetX(v4.X());
402	>	v3.SetY(v4.Y());
403	>	v3.SetZ(v4.Z());
404	>	return v3;
405		}
406
407	<	LorentzVectorXYZE toXYZ(LorentzVector v4){
407	>	LorentzVectorXYZE toXYZ(LorentzVector v4)
408	>	{
409
410	<	LorentzVectorXYZE v4_new(0,0,0,0);
411	<	v4_new.SetPx(v4.Px());
412	<	v4_new.SetPy(v4.Py());
413	<	v4_new.SetPz(v4.Pz());
414	<	v4_new.SetE(v4.E());
415	<	return v4_new;
410	>	LorentzVectorXYZE v4_new(0,0,0,0);
411	>	v4_new.SetPx(v4.Px());
412	>	v4_new.SetPy(v4.Py());
413	>	v4_new.SetPz(v4.Pz());
414	>	v4_new.SetE(v4.E());
415	>	return v4_new;
416		}
417
418	<	LorentzVector toPtEtaPhi(LorentzVectorXYZE v4){
418	>	LorentzVector toPtEtaPhi(LorentzVectorXYZE v4)
419	>	{
420
421	<	LorentzVector v4_new(0,0,0,0);
422	<	v4_new.SetPt(v4.Pt());
423	<	v4_new.SetEta(v4.Eta());
424	<	v4_new.SetPhi(v4.Phi());
425	<	v4_new.SetE(v4.E());
426	<	return v4_new;
421	>	LorentzVector v4_new(0,0,0,0);
422	>	v4_new.SetPt(v4.Pt());
423	>	v4_new.SetEta(v4.Eta());
424	>	v4_new.SetPhi(v4.Phi());
425	>	v4_new.SetE(v4.E());
426	>	return v4_new;
427		}
428
429	<	double deltaR(LorentzVector v1, LorentzVector v2){
429	>	double deltaR(LorentzVector v1, LorentzVector v2)
430	>	{
431	>
432	>	Particle p1;
433	>	p1.set_v4(v1);
434	>	Particle p2;
435	>	p2.set_v4(v2);
436	>	return p1.deltaR(p2);
437	>	}
438
439	<	Particle p1;
440	<	p1.set_v4(v1);
441	<	Particle p2;
385	<	p2.set_v4(v2);
386	<	return p1.deltaR(p2);
439	>	double double_infinity()
440	>	{
441	>	return std::numeric_limits<double>::infinity() ;
442		}
443
444	<	double double_infinity(){
445	<	return std::numeric_limits<double>::infinity() ;
444	>	int int_infinity()
445	>	{
446	>	return std::numeric_limits<int>::max() ;
447		}
448
449	<	int int_infinity(){
450	<	return std::numeric_limits<int>::max() ;
449	>	int myPow(int x, unsigned int p)
450	>	{
451	>	int i = 1;
452	>	for (unsigned int j = 1; j <= p; j++)  i *= x;
453	>	return i;
454		}
455
456	<	int myPow(int x, unsigned int p) {
457	<	int i = 1;
458	<	for (unsigned int j = 1; j <= p; j++)  i *= x;
459	<	return i;
456	>	int JetFlavor(Jet *jet)
457	>	{
458	>
459	>	EventCalc* calc = EventCalc::Instance();
460	>
461	>	std::vector< GenParticle >* genparticles = calc->GetGenParticles();
462	>	if(genparticles) {
463	>
464	>	//fill pdg IDs of all matched GenParticles
465	>	std::vector<int> matched_genparticle_ids;
466	>	for(unsigned int i=0; i<genparticles->size(); ++i) {
467	>
468	>	GenParticle genp = genparticles->at(i);
469	>
470	>	//only take status 3 particles into account
471	>	//if(genp.status()!=3) continue;
472	>
473	>	if(jet->deltaR(genp)<0.5)
474	>	matched_genparticle_ids.push_back(genp.pdgId());
475	>	}
476	>
477	>	//search for b quarks first
478	>	for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
479	>	if(abs(matched_genparticle_ids[i])==5) return matched_genparticle_ids[i];
480	>	}
481	>	//no b quark -> search for c quarks
482	>	for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
483	>	if(abs(matched_genparticle_ids[i])==4) return matched_genparticle_ids[i];
484	>	}
485	>	for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
486	>	if(abs(matched_genparticle_ids[i])==3) return matched_genparticle_ids[i];
487	>	}
488	>	for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
489	>	if(abs(matched_genparticle_ids[i])==2) return matched_genparticle_ids[i];
490	>	}
491	>	for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
492	>	if(abs(matched_genparticle_ids[i])==1) return matched_genparticle_ids[i];
493	>	}
494	>	for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
495	>	if(abs(matched_genparticle_ids[i])==21) return matched_genparticle_ids[i];
496	>	}
497	>
498	>	}
499	>
500	>	//no matched GenParticle -> return default value 0
501	>	return 0;
502	>
503		}

Diff Legend

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