96 |
|
} |
97 |
|
|
98 |
|
|
99 |
– |
|
100 |
– |
|
101 |
– |
|
102 |
– |
|
99 |
|
//HEP Tagger from Ivan |
100 |
|
|
101 |
|
bool HepTopTag(TopJet topjet) |
102 |
|
{ |
103 |
< |
|
104 |
< |
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 |
< |
} |
103 |
> |
//call variable tagger with default parameters |
104 |
> |
return variableHepTopTag(topjet); |
105 |
|
|
106 |
|
} |
107 |
|
|
184 |
– |
|
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 |
|
{ |
159 |
|
|
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 |
|
|
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; |
165 |
|
} |
166 |
< |
|
166 |
> |
|
167 |
|
Jet* nextJet(const Particle *p, std::vector<Jet> *jets) |
168 |
|
{ |
169 |
|
|
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); |
219 |
|
{ |
220 |
|
|
221 |
|
double ptrel=0; |
337 |
– |
|
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()); |
237 |
|
|
238 |
|
double deltaRmin(const Particle *p, std::vector<Jet> *jets) |
239 |
|
{ |
240 |
< |
return nextJet(p,jets)->deltaR(*p); |
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) |
327 |
|
GenParticle genp = genparticles->at(i); |
328 |
|
|
329 |
|
//only take status 3 particles into account |
330 |
< |
if(genp.status()!=3) continue; |
330 |
> |
//if(genp.status()!=3) continue; |
331 |
|
|
332 |
|
if(jet->deltaR(genp)<0.5) |
333 |
|
matched_genparticle_ids.push_back(genp.pdgId()); |