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.5 by peiffer, Mon Jul 9 08:37:32 2012 UTC vs.
Revision 1.10 by rkogler, Mon Feb 4 12:40:46 2013 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();
4	>	// global function to define a tagged jet
5
6	<	LorentzVector allsubjets(0,0,0,0);
7	<
8	<	for(int j=0; j<topjet.numberOfDaughters(); ++j){
9	<	allsubjets += topjet.subjets()[j].v4();
10	<	}
11	<	if(!allsubjets.isTimelike()){
12	<	mjet=0;
13	<	mmin=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();
20	<
21	<	if(nsubjets>=3) {
22	<
23	<	std::vector<Particle> subjets = topjet.subjets();
24	<	sort(subjets.begin(), subjets.end(), HigherPt());
25	<
26	<	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();
35	<
36	<	//minimum pairwise mass
37	<	mmin = std::min(m01,std::min(m02,m12));
38	<	}
18	>	//variable HEP Tagger from Rebekka
19
20	<	//at least 3 sub-jets
21	<	if(nsubjets<3) return false;
22	<	//minimum pairwise mass > 50 GeV/c^2
23	<	if(mmin<50) return false;
24	<	//jet mass between 140 and 250 GeV/c^2
25	<	if(mjet<140 || mjet>250) return false;
26	<
27	<	return true;
28	<	}
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	<	Jet* nextJet(const Particle *p, std::vector<Jet> *jets){
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	<	double deltarmin = double_infinity();
47	<	Jet* nextjet=0;
48	<	for(unsigned int i=0; i<jets->size();++i){
49	<	if(jets->at(i).deltaR(*p) < deltarmin){
50	<	deltarmin = jets->at(i).deltaR(*p);
51	<	nextjet = &jets->at(i);
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		}
59	–	}
65
66	<	return nextjet;
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	<	bool WTag(TopJet prunedjet,  double& mjet, int &nsubjets, double& massdrop){
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	<	nsubjets=prunedjet.numberOfDaughters();
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	<	mjet = 0;
146	<	if(prunedjet.v4().isTimelike())
147	<	mjet = prunedjet.v4().M();
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	<	//calculate mass drop for first sub-jet ordered in pt
157	<	massdrop = 0;
74	<	if(nsubjets>=1 && mjet>0){
156	>	return true;
157	>	}
158
76	–	std::vector< Particle > subjets = prunedjet.subjets();
77	–	sort(subjets.begin(), subjets.end(), HigherPt());
159
160	<	double m1 = 0;
161	<	if(subjets[0].v4().isTimelike())
162	<	m1 = subjets[0].v4().M();
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	<	massdrop = m1/mjet;
166	<	}
167	<
168	<	//at least 2 sub-jets
87	<	if(nsubjets<2) return false;
88	<	//60 GeV < pruned jet mass < 100 GeV
89	<	if(mjet <= 60 || mjet >= 100) return false;
90	<	//mass drop < 0.4
91	<	if(massdrop>=0.4) return false;
165	>	}
166	>
167	>	Jet* nextJet(const Particle *p, std::vector<Jet> *jets)
168	>	{
169
170	<	return true;
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	<	double pTrel(const Particle *p, std::vector<Jet> *jets){
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	<	double ptrel=0;
193	>	//calculate mass drop for first sub-jet ordered in pt
194	>	massdrop = 0;
195	>	if(nsubjets>=1 && mjet>0) {
196
197	<	Jet* nextjet =  nextJet(p,jets);
197	>	std::vector< Particle > subjets = prunedjet.subjets();
198	>	sort(subjets.begin(), subjets.end(), HigherPt());
199
200	<	TVector3 p3(p->v4().Px(),p->v4().Py(),p->v4().Pz());
201	<	TVector3 jet3(nextjet->v4().Px(),nextjet->v4().Py(),nextjet->v4().Pz());
200	>	double m1 = 0;
201	>	if(subjets[0].v4().isTimelike())
202	>	m1 = subjets[0].v4().M();
203
204	<	if(p3.Mag()!=0 && jet3.Mag()!=0){
205	<	double sin_alpha = (p3.Cross(jet3)).Mag()/p3.Mag()/jet3.Mag();
206	<	ptrel = p3.Mag()*sin_alpha;
207	<	}
208	<	else{
209	<	std::cout << "something strange happend in the ptrel calculation: either lepton or jet momentum is 0" <<std::endl;
210	<	}
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
114	–	return ptrel;
216		}
217
218	<	double deltaRmin(const Particle *p, std::vector<Jet> *jets){
219	<	return nextJet(p,jets)->deltaR(*p);
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	<	TVector3 toVector(LorentzVector v4){
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;
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){
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;
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){
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;
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){
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;
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){
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	<	Particle p1;
299	<	p1.set_v4(v1);
300	<	Particle p2;
164	<	p2.set_v4(v2);
165	<	return p1.deltaR(p2);
298	>	double double_infinity()
299	>	{
300	>	return std::numeric_limits<double>::infinity() ;
301		}
302
303	<	double double_infinity(){
304	<	return std::numeric_limits<double>::infinity() ;
303	>	int int_infinity()
304	>	{
305	>	return std::numeric_limits<int>::max() ;
306		}
307
308	<	int int_infinity(){
309	<	return std::numeric_limits<int>::max() ;
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 myPow(int x, unsigned int p) {
316	<	int i = 1;
317	<	for (unsigned int j = 1; j <= p; j++)  i *= x;
318	<	return i;
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		}

Diff Legend

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