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.8 by bazterra, Wed Dec 19 00:02:35 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	+
12	+
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	+	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	+
81	+
82	+	bool HepTopTagFull(TopJet topjet){
83	+
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	+	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	+
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	+	std::vector<fastjet::PseudoJet> SortedJets = sorted_by_pt(tops);
121	+
122	+
123	+	//Run the HEPTopTagger
124	+	external::HEPTopTagger tagger(*JetFinder, SortedJets[0]);
125	+
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	+	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
#	Line 96 | Line 237 | bool variableHepTopTag(TopJet topjet, do
237		}
238
239
99	–
100	–
101	–
102	–
240		//HEP Tagger from Ivan
241
242		bool HepTopTag(TopJet topjet)
243		{
244	<
245	<	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	<	}
244	>	//call variable tagger with default parameters
245	>	return variableHepTopTag(topjet);
246
247		}
248
184	–
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		{
#	Line 236 | Line 300 | bool variableTopTag(TopJet topjet, doubl
300
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
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	–	}
273	–
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;
306		}
307	<
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);
#	Line 334 | Line 360 | double pTrel(const Particle *p, std::vec
360		{
361
362		double ptrel=0;
337	–
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());
#	Line 352 | Line 378 | double pTrel(const Particle *p, std::vec
378
379		double deltaRmin(const Particle *p, std::vector<Jet> *jets)
380		{
381	<	return nextJet(p,jets)->deltaR(*p);
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)
#	Line 439 | Line 468 | int JetFlavor(Jet *jet)
468		GenParticle genp = genparticles->at(i);
469
470		//only take status 3 particles into account
471	<	if(genp.status()!=3) continue;
471	>	//if(genp.status()!=3) continue;
472
473		if(jet->deltaR(genp)<0.5)
474		matched_genparticle_ids.push_back(genp.pdgId());

Diff Legend

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