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