| 1 |
/*
|
| 2 |
* TopPairEventCandidate.cpp
|
| 3 |
*
|
| 4 |
* Created on: 9 Jul 2010
|
| 5 |
* Author: kreczko
|
| 6 |
*/
|
| 7 |
|
| 8 |
#include "../interface/TopPairEventCandidate.h"
|
| 9 |
#include <iostream>
|
| 10 |
#include <iomanip>
|
| 11 |
using namespace std;
|
| 12 |
|
| 13 |
namespace BAT {
|
| 14 |
|
| 15 |
const double TopPairEventCandidate::pdgWmass = 80.399;
|
| 16 |
const double TopPairEventCandidate::wMassWidth = 8.0; // Derived from 10% jet resolution on old jet resolution twiki page
|
| 17 |
const double TopPairEventCandidate::pdgTopMass = 172.9;
|
| 18 |
const double TopPairEventCandidate::topMassWidth = 17.0; // See above
|
| 19 |
// Above values from PDG 2011
|
| 20 |
|
| 21 |
// The following are obsolete and are only retained in
|
| 22 |
// case the old ttbar code might be useful during
|
| 23 |
// testing.
|
| 24 |
const double TopPairEventCandidate::W_mass = TopPairEventCandidate::pdgWmass;
|
| 25 |
const double TopPairEventCandidate::matched_hadronic_W_mass = TopPairEventCandidate::pdgWmass;
|
| 26 |
const double TopPairEventCandidate::matched_hadronic_W_mass_sigma = TopPairEventCandidate::wMassWidth;
|
| 27 |
const double TopPairEventCandidate::matched_leptonic_top_mass = TopPairEventCandidate::pdgTopMass;
|
| 28 |
const double TopPairEventCandidate::matched_leptonic_top_mass_sigma = TopPairEventCandidate::topMassWidth;
|
| 29 |
const double TopPairEventCandidate::matched_hadronic_top_mass = TopPairEventCandidate::pdgTopMass;
|
| 30 |
const double TopPairEventCandidate::matched_hadronic_top_mass_sigma = TopPairEventCandidate::topMassWidth;
|
| 31 |
const double TopPairEventCandidate::matched_ptratio = 0.19;
|
| 32 |
const double TopPairEventCandidate::matched_ptratio_sigma = 0.4;
|
| 33 |
const double TopPairEventCandidate::matched_pt_ttbarSystem = 0.;
|
| 34 |
const double TopPairEventCandidate::matched_pt_ttbarSystem_sigma = 50.;
|
| 35 |
const double TopPairEventCandidate::matched_HTSystem = 1;
|
| 36 |
const double TopPairEventCandidate::matched_HTSystem_sigma = 0.1;
|
| 37 |
const double TopPairEventCandidate::matched_angle = 0.95;
|
| 38 |
const double TopPairEventCandidate::matched_angle_sigma = 0.31;
|
| 39 |
|
| 40 |
NeutrinoSelectionCriterion::value TopPairEventCandidate::usedNeutrinoSelection; // = NeutrinoSelectionCriterion::chi2;
|
| 41 |
TTbarReconstructionCriterion::value TopPairEventCandidate::usedTTbarReconstruction; // = TTbarReconstructionCriterion::TopMassDifference;
|
| 42 |
|
| 43 |
|
| 44 |
TopPairEventCandidate::TopPairEventCandidate() :
|
| 45 |
Event(),
|
| 46 |
leptonFromW(),
|
| 47 |
leptonicBJet(),
|
| 48 |
hadronicBJet(),
|
| 49 |
jet1FromW(),
|
| 50 |
jet2FromW(),
|
| 51 |
neutrino1(),
|
| 52 |
neutrino2(),
|
| 53 |
leptonicW1(),
|
| 54 |
leptonicW2(),
|
| 55 |
hadronicW(),
|
| 56 |
leptonicTop1(),
|
| 57 |
leptonicTop2(),
|
| 58 |
hadronicTop(),
|
| 59 |
selectedNeutrino(0),
|
| 60 |
currentSelectedNeutrino(0),
|
| 61 |
hadronicBIndex(0),
|
| 62 |
leptonicBIndex(0),
|
| 63 |
jet1FromWIndex(0),
|
| 64 |
jet2FromWIndex(0),
|
| 65 |
doneReconstruction(false),
|
| 66 |
conversionTagger(new ConversionTagger()),
|
| 67 |
doneConversionTagging(false),
|
| 68 |
solutions(),
|
| 69 |
compareSolutions(){
|
| 70 |
}
|
| 71 |
|
| 72 |
TopPairEventCandidate::TopPairEventCandidate(const Event& event) :
|
| 73 |
Event(event),
|
| 74 |
leptonFromW(),
|
| 75 |
leptonicBJet(),
|
| 76 |
hadronicBJet(),
|
| 77 |
jet1FromW(),
|
| 78 |
jet2FromW(),
|
| 79 |
neutrino1(),
|
| 80 |
neutrino2(),
|
| 81 |
leptonicW1(),
|
| 82 |
leptonicW2(),
|
| 83 |
hadronicW(),
|
| 84 |
leptonicTop1(),
|
| 85 |
leptonicTop2(),
|
| 86 |
hadronicTop(),
|
| 87 |
selectedNeutrino(0),
|
| 88 |
currentSelectedNeutrino(0),
|
| 89 |
hadronicBIndex(0),
|
| 90 |
leptonicBIndex(0),
|
| 91 |
jet1FromWIndex(0),
|
| 92 |
jet2FromWIndex(0),
|
| 93 |
doneReconstruction(false),
|
| 94 |
conversionTagger(new ConversionTagger()),
|
| 95 |
doneConversionTagging(false) {
|
| 96 |
|
| 97 |
}
|
| 98 |
|
| 99 |
TopPairEventCandidate::~TopPairEventCandidate() {
|
| 100 |
}
|
| 101 |
|
| 102 |
bool TopPairEventCandidate::passesScrapingFilter() const {
|
| 103 |
if (tracks.size() > 10) {
|
| 104 |
if (numberOfHighPurityTracks / (1.0 * tracks.size()) > 0.25)
|
| 105 |
return true;
|
| 106 |
else
|
| 107 |
return false;
|
| 108 |
} else
|
| 109 |
return isBeamScraping == false;
|
| 110 |
}
|
| 111 |
|
| 112 |
bool TopPairEventCandidate::passesHighLevelTrigger() const {
|
| 113 |
if (isRealData()) {
|
| 114 |
if (runNumber < 140041)
|
| 115 |
return HLT(HLTriggers::HLT_Ele10_LW_L1R);
|
| 116 |
else if (runNumber >= 140041 && runNumber <= 143962)
|
| 117 |
return HLT(HLTriggers::HLT_Ele15_SW_L1R);
|
| 118 |
else if (runNumber > 143962 && runNumber <= 146427)
|
| 119 |
return HLT(HLTriggers::HLT_Ele15_SW_CaloEleId_L1R);
|
| 120 |
else if (runNumber > 146427 && runNumber <= 147116)
|
| 121 |
return HLT(HLTriggers::HLT_Ele17_SW_CaloEleId_L1R);
|
| 122 |
else if (runNumber > 147116 && runNumber <= 148818)
|
| 123 |
return HLT(HLTriggers::HLT_Ele17_SW_TightEleId_L1R);
|
| 124 |
else if (runNumber >= 148819 && runNumber < 149181)
|
| 125 |
return HLT(HLTriggers::HLT_Ele22_SW_TighterEleId_L1R);
|
| 126 |
//return HLT(HLTriggers::HLT_Ele22_SW_TighterEleId_L1R_v2);
|
| 127 |
else if(runNumber >= 149181 && runNumber < 160000)
|
| 128 |
return HLT(HLTriggers::HLT_Ele22_SW_TighterEleId_L1R);
|
| 129 |
//return HLT(HLTriggers::HLT_Ele22_SW_TighterEleId_L1R_v3);
|
| 130 |
else if(runNumber > 160000)
|
| 131 |
return (HLT(HLTriggers::HLT_Ele27_CaloIdVT_CaloIsoT_TrkIdT_TrkIsoT) ||
|
| 132 |
HLT(HLTriggers::HLT_Ele25_CaloIdVT_CaloIsoT_TrkIdT_TrkIsoT_TriCentralJet30)
|
| 133 |
|| HLT(HLTriggers::HLT_Ele25_CaloIdVT_CaloIsoT_TrkIdT_TrkIsoT_TriCentralPFJet30)
|
| 134 |
|| HLT(HLTriggers::HLT_Mu30) || HLT(HLTriggers::HLT_IsoMu20)
|
| 135 |
|| HLT(HLTriggers::HLT_IsoMu17) || HLT(HLTriggers::HLT_IsoMu17_eta2p1)
|
| 136 |
|| HLT(HLTriggers::HLT_IsoMu24) || HLT(HLTriggers::HLT_IsoMu24_eta2p1)
|
| 137 |
|| HLT(HLTriggers::HLT_IsoMu30_eta2p1));
|
| 138 |
else return false;
|
| 139 |
} else return true; // do not use HLT for MC
|
| 140 |
}
|
| 141 |
|
| 142 |
bool TopPairEventCandidate::hasOneGoodPrimaryVertex() const {
|
| 143 |
return goodVertices.size() >= 1;
|
| 144 |
}
|
| 145 |
|
| 146 |
bool TopPairEventCandidate::hasOnlyOneGoodIsolatedElectron() const {
|
| 147 |
if(Event::usePFIsolation)
|
| 148 |
return goodPFIsolatedElectrons.size() == 1;
|
| 149 |
else
|
| 150 |
return goodIsolatedElectrons.size() == 1;
|
| 151 |
}
|
| 152 |
|
| 153 |
bool TopPairEventCandidate::isolatedElectronDoesNotComeFromConversion() const {
|
| 154 |
bool passConversion = false;
|
| 155 |
if (Event::usePFIsolation) {
|
| 156 |
if (goodPFIsolatedElectrons.size() > 0)
|
| 157 |
passConversion = goodPFIsolatedElectrons.front()->isFromConversion() == false;
|
| 158 |
} else {
|
| 159 |
if (goodIsolatedElectrons.size() > 0)
|
| 160 |
passConversion = goodIsolatedElectrons.front()->isFromConversion() == false;
|
| 161 |
}
|
| 162 |
|
| 163 |
return passConversion;
|
| 164 |
}
|
| 165 |
|
| 166 |
bool TopPairEventCandidate::isolatedElectronNotTaggedAsFromConversion() const {
|
| 167 |
bool passConversion = false;
|
| 168 |
ElectronPointer electron;
|
| 169 |
if (Event::usePFIsolation) {
|
| 170 |
if (goodPFIsolatedElectrons.size() > 0)
|
| 171 |
electron = goodPFIsolatedElectrons.front();
|
| 172 |
} else {
|
| 173 |
if (goodIsolatedElectrons.size() > 0)
|
| 174 |
electron = goodIsolatedElectrons.front();
|
| 175 |
}
|
| 176 |
if (electron != 0) {
|
| 177 |
if (useCustomConversionTagger) {
|
| 178 |
conversionTagger->calculateConversionVariables(electron, tracks, 3.8, 0.45);
|
| 179 |
passConversion = conversionTagger->isFromConversion(0.02, 0.02) == false;
|
| 180 |
} else {
|
| 181 |
passConversion = electron->isTaggedAsConversion(0.02, 0.02) == false;
|
| 182 |
}
|
| 183 |
}
|
| 184 |
|
| 185 |
return passConversion;
|
| 186 |
}
|
| 187 |
|
| 188 |
bool TopPairEventCandidate::hasNoIsolatedMuon() const {
|
| 189 |
return looseMuons.size() == 0;
|
| 190 |
}
|
| 191 |
|
| 192 |
bool TopPairEventCandidate::hasAtLeastOneGoodJet() const {
|
| 193 |
return goodJets.size() >= 1;
|
| 194 |
}
|
| 195 |
|
| 196 |
bool TopPairEventCandidate::hasAtLeastTwoGoodJets() const {
|
| 197 |
return goodJets.size() >= 2;
|
| 198 |
}
|
| 199 |
|
| 200 |
bool TopPairEventCandidate::hasAtLeastThreeGoodJets() const {
|
| 201 |
return goodJets.size() >= 3;
|
| 202 |
}
|
| 203 |
|
| 204 |
bool TopPairEventCandidate::hasAtLeastFourGoodJets() const {
|
| 205 |
return goodJets.size() >= 4;
|
| 206 |
}
|
| 207 |
|
| 208 |
bool TopPairEventCandidate::isNotAZBosonEvent() const {
|
| 209 |
float invariantMass = 0;
|
| 210 |
bool isZEvent = false;
|
| 211 |
ElectronPointer isoElectron;
|
| 212 |
|
| 213 |
if (Event::usePFIsolation && goodPFIsolatedElectrons.size() > 0)
|
| 214 |
isoElectron = goodPFIsolatedElectrons.front();
|
| 215 |
else if (goodIsolatedElectrons.size() > 0)
|
| 216 |
isoElectron = goodIsolatedElectrons.front();
|
| 217 |
|
| 218 |
if (isoElectron != NULL && allElectrons.size() > 1) {
|
| 219 |
for (unsigned int index = 0; index < allElectrons.size(); ++index) {
|
| 220 |
const ElectronPointer looseElectron = allElectrons.at(index);
|
| 221 |
bool passLooseIso = false;
|
| 222 |
|
| 223 |
if (Event::usePFIsolation)
|
| 224 |
passLooseIso = looseElectron->isLoose() && looseElectron->pfIsolation() < 1.;
|
| 225 |
else
|
| 226 |
passLooseIso = looseElectron->isLoose() && looseElectron->relativeIsolation() < 1.;
|
| 227 |
|
| 228 |
if (passLooseIso)
|
| 229 |
invariantMass = isoElectron->invariantMass(looseElectron);
|
| 230 |
else
|
| 231 |
invariantMass = 0;
|
| 232 |
|
| 233 |
bool passesLowerLimit = invariantMass > 76;
|
| 234 |
bool passesUpperLimit = invariantMass < 106;
|
| 235 |
if (passesLowerLimit && passesUpperLimit)
|
| 236 |
isZEvent = true;
|
| 237 |
}
|
| 238 |
|
| 239 |
}
|
| 240 |
|
| 241 |
|
| 242 |
return isZEvent == false;
|
| 243 |
}
|
| 244 |
|
| 245 |
bool TopPairEventCandidate::passesMETCut() const{
|
| 246 |
return met->et() > 20.;
|
| 247 |
}
|
| 248 |
|
| 249 |
bool TopPairEventCandidate::passesAsymmetricJetCuts() const {
|
| 250 |
if(goodJets.size() < 2)// good jets have a cut of 30 GeV!
|
| 251 |
return false;
|
| 252 |
JetPointer leadingJet = goodJets.front();
|
| 253 |
JetPointer secondLeadingJet = goodJets.at(1);
|
| 254 |
return leadingJet->pt() > 70 && secondLeadingJet->pt() > 50;
|
| 255 |
|
| 256 |
}
|
| 257 |
|
| 258 |
bool TopPairEventCandidate::hasAtLeastOneBtag() const {
|
| 259 |
return goodBJets.size() >= 1;
|
| 260 |
}
|
| 261 |
|
| 262 |
bool TopPairEventCandidate::hasAtLeastTwoBtags() const{
|
| 263 |
return goodBJets.size() >= 2;
|
| 264 |
}
|
| 265 |
|
| 266 |
bool TopPairEventCandidate::passesFullTTbarEPlusJetSelection() const {
|
| 267 |
// unsigned int newstep = (int) TTbarEPlusJetsSelection::NUMBER_OF_SELECTION_STEPS - 1;
|
| 268 |
///////////////////////////// from Lucaz
|
| 269 |
return passesSelectionStepUpTo(TTbarEPlusJetsSelection::AtLeastFourGoodJets);//(TTbarEPlusJetsSelection::Step) newstep);
|
| 270 |
//
|
| 271 |
}
|
| 272 |
|
| 273 |
bool TopPairEventCandidate::passesSelectionStepUpTo(enum TTbarEPlusJetsSelection::Step step) const {
|
| 274 |
if (step == TTbarEPlusJetsSelection::FilterOutScraping)
|
| 275 |
return passesSelectionStep(step);
|
| 276 |
else {
|
| 277 |
unsigned int newstep = (int) step - 1;
|
| 278 |
return passesSelectionStep(step) && passesSelectionStepUpTo((TTbarEPlusJetsSelection::Step) newstep);
|
| 279 |
}
|
| 280 |
}
|
| 281 |
|
| 282 |
bool TopPairEventCandidate::passesNMinus1(enum TTbarEPlusJetsSelection::Step omitted) const {
|
| 283 |
bool passes(true);
|
| 284 |
|
| 285 |
for (unsigned int cut = 0; cut < TTbarEPlusJetsSelection::NUMBER_OF_SELECTION_STEPS; ++cut) {
|
| 286 |
if (cut == (unsigned int) omitted)
|
| 287 |
continue;
|
| 288 |
passes = passes && passesSelectionStep((TTbarEPlusJetsSelection::Step) cut);
|
| 289 |
}
|
| 290 |
return passes;
|
| 291 |
}
|
| 292 |
|
| 293 |
bool TopPairEventCandidate::passesSelectionStep(enum TTbarEPlusJetsSelection::Step step) const {
|
| 294 |
switch (step) {
|
| 295 |
case TTbarEPlusJetsSelection::FilterOutScraping:
|
| 296 |
return passesScrapingFilter();
|
| 297 |
case TTbarEPlusJetsSelection::HighLevelTrigger:
|
| 298 |
return passesHighLevelTrigger();
|
| 299 |
case TTbarEPlusJetsSelection::GoodPrimaryvertex:
|
| 300 |
return hasOneGoodPrimaryVertex();
|
| 301 |
case TTbarEPlusJetsSelection::OneIsolatedElectron:
|
| 302 |
return hasOnlyOneGoodIsolatedElectron();
|
| 303 |
case TTbarEPlusJetsSelection::ConversionRejection:
|
| 304 |
return isolatedElectronDoesNotComeFromConversion();
|
| 305 |
case TTbarEPlusJetsSelection::ConversionFinder:
|
| 306 |
return isolatedElectronNotTaggedAsFromConversion();
|
| 307 |
case TTbarEPlusJetsSelection::LooseMuonVeto:
|
| 308 |
return hasNoIsolatedMuon();
|
| 309 |
case TTbarEPlusJetsSelection::AtLeastOneGoodJets:
|
| 310 |
return hasAtLeastOneGoodJet();
|
| 311 |
case TTbarEPlusJetsSelection::AtLeastTwoGoodJets:
|
| 312 |
return hasAtLeastTwoGoodJets();
|
| 313 |
case TTbarEPlusJetsSelection::AtLeastThreeGoodJets:
|
| 314 |
return hasAtLeastThreeGoodJets();
|
| 315 |
case TTbarEPlusJetsSelection::AtLeastFourGoodJets:
|
| 316 |
return hasAtLeastFourGoodJets();
|
| 317 |
case TTbarEPlusJetsSelection::Zveto:
|
| 318 |
return isNotAZBosonEvent();
|
| 319 |
case TTbarEPlusJetsSelection::MissingTransverseEnergy:
|
| 320 |
return passesMETCut();
|
| 321 |
case TTbarEPlusJetsSelection::AsymmetricJetCuts:
|
| 322 |
return passesAsymmetricJetCuts();
|
| 323 |
case TTbarEPlusJetsSelection::AtLeastOneBtag:
|
| 324 |
return hasAtLeastOneBtag();
|
| 325 |
case TTbarEPlusJetsSelection::AtLeastTwoBtags:
|
| 326 |
return hasAtLeastTwoBtags();
|
| 327 |
default:
|
| 328 |
return false;
|
| 329 |
}
|
| 330 |
}
|
| 331 |
|
| 332 |
bool TopPairEventCandidate::passesRelIsoSelection() const{
|
| 333 |
bool passesFirst3 = passesSelectionStepUpTo(TTbarEPlusJetsSelection::GoodPrimaryvertex);
|
| 334 |
bool passGoodElectrons = goodElectrons.size() > 0 && goodIsolatedElectrons.size() < 2;
|
| 335 |
bool passesBothIsolationvetos = false;
|
| 336 |
if (passGoodElectrons) {
|
| 337 |
const ElectronPointer electron = MostIsolatedElectron();
|
| 338 |
if (electron->isGood()) {
|
| 339 |
if (useCustomConversionTagger) {
|
| 340 |
conversionTagger->calculateConversionVariables(electron, tracks, 3.8, 0.45);
|
| 341 |
passesBothIsolationvetos = electron->isFromConversion() == false && conversionTagger->isFromConversion(
|
| 342 |
0.02, 0.02) == false;
|
| 343 |
}
|
| 344 |
else{
|
| 345 |
passesBothIsolationvetos = electron->isFromConversion() == false && electron->isTaggedAsConversion(
|
| 346 |
0.02, 0.02) == false;
|
| 347 |
}
|
| 348 |
}
|
| 349 |
|
| 350 |
}
|
| 351 |
bool muonVeto = hasNoIsolatedMuon();
|
| 352 |
bool Zveto = isNotAZBosonEvent();
|
| 353 |
return passesFirst3 && passGoodElectrons && passesBothIsolationvetos && muonVeto && Zveto;
|
| 354 |
}
|
| 355 |
|
| 356 |
bool TopPairEventCandidate::passesRelIsoControlSelection() const{
|
| 357 |
bool passesFirst3 = passesSelectionStepUpTo(TTbarEPlusJetsSelection::GoodPrimaryvertex);
|
| 358 |
bool passGoodElectrons = allElectrons.size() > 0 && goodIsolatedElectrons.size() < 2;
|
| 359 |
bool passesBothIsolationvetos = false;
|
| 360 |
if (passGoodElectrons) {
|
| 361 |
const ElectronPointer electron = MostIsolatedElectron();
|
| 362 |
if (electron->isQCDElectron()) {
|
| 363 |
if (useCustomConversionTagger) {
|
| 364 |
conversionTagger->calculateConversionVariables(electron, tracks, 3.8, 0.45);
|
| 365 |
passesBothIsolationvetos = electron->isFromConversion() == false && conversionTagger->isFromConversion(
|
| 366 |
0.02, 0.02) == false;
|
| 367 |
} else {
|
| 368 |
passesBothIsolationvetos = electron->isFromConversion() == false && electron->isTaggedAsConversion(
|
| 369 |
0.02, 0.02) == false;
|
| 370 |
}
|
| 371 |
}
|
| 372 |
|
| 373 |
}
|
| 374 |
bool muonVeto = hasNoIsolatedMuon();
|
| 375 |
bool Zveto = isNotAZBosonEvent();
|
| 376 |
return passesFirst3 && passGoodElectrons && passesBothIsolationvetos && muonVeto && Zveto;
|
| 377 |
}
|
| 378 |
|
| 379 |
bool TopPairEventCandidate::passesPFIsoSelection() const{
|
| 380 |
bool passesFirst3 = passesSelectionStepUpTo(TTbarEPlusJetsSelection::GoodPrimaryvertex);
|
| 381 |
bool passGoodElectrons = goodElectrons.size() > 0 && goodPFIsolatedElectrons.size() < 2;
|
| 382 |
bool passesBothIsolationvetos = false;
|
| 383 |
if (passGoodElectrons) {
|
| 384 |
const ElectronPointer electron = MostPFIsolatedElectron();
|
| 385 |
if (electron->isGood()) {
|
| 386 |
if (useCustomConversionTagger) {
|
| 387 |
conversionTagger->calculateConversionVariables(electron, tracks, 3.8, 0.45);
|
| 388 |
passesBothIsolationvetos = electron->isFromConversion() == false && conversionTagger->isFromConversion(
|
| 389 |
0.02, 0.02) == false;
|
| 390 |
}
|
| 391 |
else{
|
| 392 |
passesBothIsolationvetos = electron->isFromConversion() == false && electron->isTaggedAsConversion(
|
| 393 |
0.02, 0.02) == false;
|
| 394 |
}
|
| 395 |
}
|
| 396 |
|
| 397 |
}
|
| 398 |
bool muonVeto = hasNoIsolatedMuon();
|
| 399 |
bool Zveto = isNotAZBosonEvent();
|
| 400 |
return passesFirst3 && passGoodElectrons && passesBothIsolationvetos && muonVeto && Zveto;
|
| 401 |
}
|
| 402 |
|
| 403 |
bool TopPairEventCandidate::passesPFIsoControlSelection() const{
|
| 404 |
bool passesFirst3 = passesSelectionStepUpTo(TTbarEPlusJetsSelection::GoodPrimaryvertex);
|
| 405 |
bool passGoodElectrons = allElectrons.size() > 0 && goodPFIsolatedElectrons.size() < 2;
|
| 406 |
bool passesBothIsolationvetos = false;
|
| 407 |
if (passGoodElectrons) {
|
| 408 |
const ElectronPointer electron = MostPFIsolatedElectron();
|
| 409 |
if (electron->isQCDElectron()) {
|
| 410 |
if (useCustomConversionTagger) {
|
| 411 |
conversionTagger->calculateConversionVariables(electron, tracks, 3.8, 0.45);
|
| 412 |
passesBothIsolationvetos = electron->isFromConversion() == false && conversionTagger->isFromConversion(
|
| 413 |
0.02, 0.02) == false;
|
| 414 |
} else {
|
| 415 |
passesBothIsolationvetos = electron->isFromConversion() == false && electron->isTaggedAsConversion(
|
| 416 |
0.02, 0.02) == false;
|
| 417 |
}
|
| 418 |
}
|
| 419 |
|
| 420 |
}
|
| 421 |
bool muonVeto = hasNoIsolatedMuon();
|
| 422 |
bool Zveto = isNotAZBosonEvent();
|
| 423 |
return passesFirst3 && passGoodElectrons && passesBothIsolationvetos && muonVeto && Zveto;
|
| 424 |
}
|
| 425 |
|
| 426 |
bool TopPairEventCandidate::passesConversionSelection() const {
|
| 427 |
bool passesFirst6 = passesSelectionStepUpTo(TTbarEPlusJetsSelection::Zveto);
|
| 428 |
bool isConversion1 = isolatedElectronDoesNotComeFromConversion() == false;
|
| 429 |
bool isConversion2 = isolatedElectronNotTaggedAsFromConversion() == false;
|
| 430 |
bool atLeast4Jets = hasAtLeastFourGoodJets();
|
| 431 |
return passesFirst6 && (isConversion1 || isConversion2) && atLeast4Jets;
|
| 432 |
}
|
| 433 |
|
| 434 |
bool TopPairEventCandidate::passesAntiIsolationSelection() const {
|
| 435 |
//require at least one good electron and no isolated good electrons
|
| 436 |
if (!(goodElectrons.size() > 0 && goodIsolatedElectrons.size() == 0))
|
| 437 |
return false;
|
| 438 |
|
| 439 |
bool passesFirst3 = passesSelectionStep(TTbarEPlusJetsSelection::GoodPrimaryvertex);
|
| 440 |
|
| 441 |
|
| 442 |
bool muonVeto = passesSelectionStep(TTbarEPlusJetsSelection::LooseMuonVeto);
|
| 443 |
bool zveto = passesSelectionStep(TTbarEPlusJetsSelection::Zveto);
|
| 444 |
bool conversionVeto = (goodElectrons.front()->isFromConversion() || goodElectrons.front()->isTaggedAsConversion(
|
| 445 |
0.2, 0.2)) == false;
|
| 446 |
bool jets = hasAtLeastFourGoodJets();
|
| 447 |
return passesFirst3 && muonVeto && zveto && conversionVeto && jets;
|
| 448 |
}
|
| 449 |
|
| 450 |
void TopPairEventCandidate::reconstructTTbar(LeptonPointer lepton) {
|
| 451 |
if (goodJets.size() < 4)
|
| 452 |
throw ReconstructionException("Not enough jets available to reconstruct top event using Mass Equality method.");
|
| 453 |
leptonFromW = lepton;
|
| 454 |
selectedNeutrino = 0;
|
| 455 |
currentSelectedNeutrino = 0;
|
| 456 |
reconstructNeutrinos();
|
| 457 |
double chosen_TopMassDifference(9999999.);
|
| 458 |
double chosen_Chi2Total(9999999.);
|
| 459 |
|
| 460 |
for (unsigned short hadBindex = 0; hadBindex < goodJets.size(); ++hadBindex) {
|
| 461 |
for (unsigned short lepBindex = 0; lepBindex < goodJets.size(); ++lepBindex) {
|
| 462 |
if (lepBindex == hadBindex)
|
| 463 |
continue;
|
| 464 |
for (unsigned short jet1Index = 0; jet1Index < goodJets.size(); ++jet1Index) {
|
| 465 |
if (jet1Index == lepBindex || jet1Index == hadBindex)
|
| 466 |
continue;
|
| 467 |
for (unsigned short jet2Index = 0; jet2Index < goodJets.size(); ++jet2Index) {
|
| 468 |
if (jet2Index == jet1Index || jet2Index == lepBindex || jet2Index == hadBindex)
|
| 469 |
continue;
|
| 470 |
hadronicBJet = goodJets.at(hadBindex);
|
| 471 |
leptonicBJet = goodJets.at(lepBindex);
|
| 472 |
jet1FromW = goodJets.at(jet1Index);
|
| 473 |
jet2FromW = goodJets.at(jet2Index);
|
| 474 |
|
| 475 |
leptonicW1 = ParticlePointer(new Particle(*neutrino1 + *leptonFromW));
|
| 476 |
leptonicW2 = ParticlePointer(new Particle(*neutrino2 + *leptonFromW));
|
| 477 |
hadronicW = ParticlePointer(new Particle(*jet1FromW + *jet2FromW));
|
| 478 |
leptonicTop1 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW1));
|
| 479 |
leptonicTop2 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW2));
|
| 480 |
hadronicTop = ParticlePointer(new Particle(*hadronicBJet + *hadronicW));
|
| 481 |
fillHypotheses();
|
| 482 |
selectNeutrinoSolution();
|
| 483 |
double TopMassDifference = calculateTopMassDifference(currentSelectedNeutrino);
|
| 484 |
double chi2 = getTotalChi2(currentSelectedNeutrino);
|
| 485 |
switch (usedTTbarReconstruction) {
|
| 486 |
case TTbarReconstructionCriterion::TopMassDifference:
|
| 487 |
if (TopMassDifference < chosen_TopMassDifference) {
|
| 488 |
hadronicBIndex = hadBindex;
|
| 489 |
leptonicBIndex = lepBindex;
|
| 490 |
jet1FromWIndex = jet1Index;
|
| 491 |
jet2FromWIndex = jet2Index;
|
| 492 |
chosen_TopMassDifference = TopMassDifference;
|
| 493 |
selectedNeutrino = currentSelectedNeutrino;
|
| 494 |
}
|
| 495 |
break;
|
| 496 |
|
| 497 |
case TTbarReconstructionCriterion::chi2:
|
| 498 |
if (chi2 < chosen_Chi2Total) {
|
| 499 |
hadronicBIndex = hadBindex;
|
| 500 |
leptonicBIndex = lepBindex;
|
| 501 |
jet1FromWIndex = jet1Index;
|
| 502 |
jet2FromWIndex = jet2Index;
|
| 503 |
chosen_Chi2Total = chi2;
|
| 504 |
selectedNeutrino = currentSelectedNeutrino;
|
| 505 |
}
|
| 506 |
break;
|
| 507 |
}
|
| 508 |
}
|
| 509 |
}
|
| 510 |
}
|
| 511 |
}
|
| 512 |
std::sort(solutions.begin(), solutions.end(), compareSolutions);
|
| 513 |
hadronicBJet = goodJets.at(hadronicBIndex);
|
| 514 |
leptonicBJet = goodJets.at(leptonicBIndex);
|
| 515 |
jet1FromW = goodJets.at(jet1FromWIndex);
|
| 516 |
jet2FromW = goodJets.at(jet2FromWIndex);
|
| 517 |
leptonicW1 = ParticlePointer(new Particle(*neutrino1 + *leptonFromW));
|
| 518 |
leptonicW2 = ParticlePointer(new Particle(*neutrino2 + *leptonFromW));
|
| 519 |
hadronicW = ParticlePointer(new Particle(*jet1FromW + *jet2FromW));
|
| 520 |
leptonicTop1 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW1));
|
| 521 |
leptonicTop2 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW2));
|
| 522 |
hadronicTop = ParticlePointer(new Particle(*hadronicBJet + *hadronicW));
|
| 523 |
if (selectedNeutrino == 1)
|
| 524 |
ttbarResonance = ParticlePointer(new Particle(*leptonicTop1 + *hadronicTop));
|
| 525 |
else
|
| 526 |
ttbarResonance = ParticlePointer(new Particle(*leptonicTop2 + *hadronicTop));
|
| 527 |
doneReconstruction = true;
|
| 528 |
}
|
| 529 |
|
| 530 |
|
| 531 |
void TopPairEventCandidate::reconstructNeutrinos() {
|
| 532 |
boost::array<double, 2> neutrinoPzs = computeNeutrinoPz();
|
| 533 |
double energy1 = sqrt(met->et() * met->et() + neutrinoPzs.at(0) * neutrinoPzs.at(0));
|
| 534 |
double energy2 = sqrt(met->et() * met->et() + neutrinoPzs.at(1) * neutrinoPzs.at(1));
|
| 535 |
neutrino1 = ParticlePointer(new Particle(energy1, met->px(), met->py(), neutrinoPzs.at(0)));
|
| 536 |
neutrino2 = ParticlePointer(new Particle(energy2, met->px(), met->py(), neutrinoPzs.at(1)));
|
| 537 |
|
| 538 |
if (isnan(neutrino1->energy()) && isnan(neutrino2->energy()))
|
| 539 |
throw ReconstructionException("No physical neutrino solution found");
|
| 540 |
else if (isnan(neutrino1->energy()))
|
| 541 |
neutrino1 = neutrino2;
|
| 542 |
else if (isnan(neutrino2->energy()))
|
| 543 |
neutrino2 = neutrino1;
|
| 544 |
}
|
| 545 |
|
| 546 |
const boost::array<double, 2> TopPairEventCandidate::computeNeutrinoPz() {
|
| 547 |
if (leptonFromW == 0)
|
| 548 |
throw ReconstructionException("Could not reconstruct neutrinos: no isolated leptons found");
|
| 549 |
if (met->energy() == 0)
|
| 550 |
throw ReconstructionException("Could not reconstruct neutrinos: no MET found");
|
| 551 |
boost::array<double, 2> neutrinoPzs;
|
| 552 |
|
| 553 |
// Two unknowns, neutrino E and pz. But since neutrino mass = 0, can
|
| 554 |
// solve for pz. Result is long and messy, as follows.
|
| 555 |
|
| 556 |
double pz1(0), pz2(0);
|
| 557 |
double ee = leptonFromW->energy();
|
| 558 |
double pxe = leptonFromW->px();
|
| 559 |
double pye = leptonFromW->py();
|
| 560 |
double pze = leptonFromW->pz();
|
| 561 |
double pxnu = met->px();
|
| 562 |
double pynu = met->py();
|
| 563 |
|
| 564 |
double a = W_mass * W_mass + 2.0 * pxe * pxnu + 2.0 * pye * pynu;
|
| 565 |
double A = 4.0 * (ee * ee - pze * pze);
|
| 566 |
double B = -4.0 * a * pze;
|
| 567 |
double C = 4.0 * ee * ee * (pxnu * pxnu + pynu * pynu) - a * a;
|
| 568 |
|
| 569 |
double tmproot = B * B - 4.0 * A * C;
|
| 570 |
if (tmproot < 0) {
|
| 571 |
pz1 = pz2 = -B / (2 * A); // Take real part
|
| 572 |
} else {
|
| 573 |
pz1 = (-B + TMath::Sqrt(tmproot)) / (2.0 * A);
|
| 574 |
pz2 = (-B - TMath::Sqrt(tmproot)) / (2.0 * A);
|
| 575 |
|
| 576 |
}
|
| 577 |
neutrinoPzs[0] = pz1;
|
| 578 |
neutrinoPzs[1] = pz2;
|
| 579 |
return neutrinoPzs;
|
| 580 |
}
|
| 581 |
|
| 582 |
|
| 583 |
void TopPairEventCandidate::fillHypotheses() {
|
| 584 |
TtbarHypothesisPointer hypothesis1(fillHypothesis(1));
|
| 585 |
TtbarHypothesisPointer hypothesis2(fillHypothesis(2));
|
| 586 |
solutions.push_back(hypothesis1);
|
| 587 |
solutions.push_back(hypothesis2);
|
| 588 |
|
| 589 |
}
|
| 590 |
|
| 591 |
const TtbarHypothesisPointer TopPairEventCandidate::fillHypothesis(unsigned short int neutrinoSolution) {
|
| 592 |
TtbarHypothesisPointer hypothesis(new TtbarHypothesis());
|
| 593 |
hypothesis->leptonFromW = leptonFromW;
|
| 594 |
hypothesis->leptonicBjet = leptonicBJet;
|
| 595 |
hypothesis->hadronicBJet = hadronicBJet;
|
| 596 |
hypothesis->jet1FromW = jet1FromW;
|
| 597 |
hypothesis->jet2FromW = jet2FromW;
|
| 598 |
hypothesis->hadronicW = hadronicW;
|
| 599 |
hypothesis->hadronicTop = hadronicTop;
|
| 600 |
hypothesis->hadronicChi2 = getHadronicChi2();
|
| 601 |
if(neutrinoSolution == 1) {
|
| 602 |
hypothesis->neutrinoFromW = neutrino1;
|
| 603 |
hypothesis->leptonicW = leptonicW1;
|
| 604 |
hypothesis->leptonicTop = leptonicTop1;
|
| 605 |
}
|
| 606 |
else {
|
| 607 |
hypothesis->neutrinoFromW = neutrino2;
|
| 608 |
hypothesis->leptonicW = leptonicW2;
|
| 609 |
hypothesis->leptonicTop = leptonicTop2;
|
| 610 |
}
|
| 611 |
|
| 612 |
hypothesis->totalChi2 = getTotalChi2(neutrinoSolution);
|
| 613 |
hypothesis->globalChi2 = getGlobalChi2(neutrinoSolution);
|
| 614 |
hypothesis->leptonicChi2 = getLeptonicChi2(neutrinoSolution);
|
| 615 |
ParticlePointer resonance(new Particle(*hypothesis->leptonicTop + *hypothesis->hadronicTop));
|
| 616 |
hypothesis->resonance = resonance;
|
| 617 |
return hypothesis;
|
| 618 |
}
|
| 619 |
|
| 620 |
void TopPairEventCandidate::selectNeutrinoSolution() {
|
| 621 |
|
| 622 |
if (leptonicTop1->mass() < 0 && leptonicTop2->mass() < 0) {
|
| 623 |
inspectReconstructedEvent();
|
| 624 |
throw ReconstructionException("No valid neutrino solution found");
|
| 625 |
} else if (leptonicTop1->mass() < 0 && leptonicTop2->mass() > 0) {
|
| 626 |
currentSelectedNeutrino = 2;
|
| 627 |
} else if (leptonicTop1->mass() > 0 && leptonicTop2->mass() < 0) {
|
| 628 |
currentSelectedNeutrino = 1;
|
| 629 |
} else {// both solutions give positive mass
|
| 630 |
switch (usedNeutrinoSelection) {
|
| 631 |
case NeutrinoSelectionCriterion::TopMassDifference:
|
| 632 |
fabs(leptonicTop1->mass()-hadronicTop->mass()) < fabs(leptonicTop2->mass()-hadronicTop->mass()) ?
|
| 633 |
currentSelectedNeutrino = 1 : currentSelectedNeutrino = 2;
|
| 634 |
break;
|
| 635 |
case NeutrinoSelectionCriterion::chi2:
|
| 636 |
getTotalChi2(1) < getTotalChi2(2) ? currentSelectedNeutrino = 1 : currentSelectedNeutrino = 2;
|
| 637 |
break;
|
| 638 |
|
| 639 |
case NeutrinoSelectionCriterion::pzClosestToLepton:
|
| 640 |
fabs(neutrino1->pz() - leptonFromW->pz()) < fabs(neutrino2->pz()
|
| 641 |
- leptonFromW->pz()) ? currentSelectedNeutrino = 1 : currentSelectedNeutrino = 2;
|
| 642 |
break;
|
| 643 |
|
| 644 |
case NeutrinoSelectionCriterion::mostCentral:
|
| 645 |
fabs(neutrino1->pz()) < fabs(neutrino2->pz()) ? currentSelectedNeutrino = 1 : currentSelectedNeutrino = 2;
|
| 646 |
break;
|
| 647 |
|
| 648 |
case NeutrinoSelectionCriterion::pzClosestToLeptonOrMostcentralIfAbove300:
|
| 649 |
fabs(neutrino1->pz() - leptonFromW->pz()) < fabs(neutrino2->pz()
|
| 650 |
- leptonFromW->pz()) ? currentSelectedNeutrino = 1 : currentSelectedNeutrino = 2;
|
| 651 |
if (fabs(neutrino1->pz()) > 300 || fabs(neutrino2->pz()) > 300)
|
| 652 |
fabs(neutrino1->pz()) < fabs(neutrino2->pz()) ? currentSelectedNeutrino = 1 : currentSelectedNeutrino
|
| 653 |
= 2;
|
| 654 |
break;
|
| 655 |
|
| 656 |
case NeutrinoSelectionCriterion::largestValueOfCosine:
|
| 657 |
TVector3 p3W, p3e;
|
| 658 |
//TODO clean up
|
| 659 |
p3W = leptonicW1->getFourVector().Vect();
|
| 660 |
p3e = leptonFromW->getFourVector().Vect();
|
| 661 |
|
| 662 |
double sinthcm1 = 2. * (p3e.Perp(p3W)) / W_mass;
|
| 663 |
p3W = leptonicW2->getFourVector().Vect();
|
| 664 |
double sinthcm2 = 2. * (p3e.Perp(p3W)) / W_mass;
|
| 665 |
|
| 666 |
double costhcm1 = TMath::Sqrt(1. - sinthcm1 * sinthcm1);
|
| 667 |
double costhcm2 = TMath::Sqrt(1. - sinthcm2 * sinthcm2);
|
| 668 |
costhcm1 > costhcm2 ? currentSelectedNeutrino = 1 : currentSelectedNeutrino = 2;
|
| 669 |
break;
|
| 670 |
|
| 671 |
}
|
| 672 |
}
|
| 673 |
|
| 674 |
}
|
| 675 |
|
| 676 |
double TopPairEventCandidate::calculateTopMassDifference(unsigned short int neutrinoSolution) const {
|
| 677 |
|
| 678 |
double LeptonicTop1MassDifference = fabs(leptonicTop1->mass()-hadronicTop->mass());
|
| 679 |
double LeptonicTop2MassDifference = fabs(leptonicTop2->mass()-hadronicTop->mass());
|
| 680 |
|
| 681 |
if (neutrinoSolution == 1)
|
| 682 |
return LeptonicTop1MassDifference;
|
| 683 |
else
|
| 684 |
return LeptonicTop2MassDifference;
|
| 685 |
|
| 686 |
}
|
| 687 |
|
| 688 |
double TopPairEventCandidate::getTotalChi2() {
|
| 689 |
double totalChi2(9999999);
|
| 690 |
double firstTotalChi2 = getTotalChi2(1);
|
| 691 |
double secondTotalChi2 = getTotalChi2(2);
|
| 692 |
selectedNeutrino == 1 ? totalChi2 = firstTotalChi2 : totalChi2 = secondTotalChi2;
|
| 693 |
// if (firstTotalChi2 < secondTotalChi2) {
|
| 694 |
// selectedNeutrino = 1;
|
| 695 |
// totalChi2 = firstTotalChi2;
|
| 696 |
// } else {
|
| 697 |
// selectedNeutrino = 2;
|
| 698 |
// totalChi2 = secondTotalChi2;
|
| 699 |
// }
|
| 700 |
return totalChi2;
|
| 701 |
}
|
| 702 |
|
| 703 |
double TopPairEventCandidate::getTotalChi2(unsigned short int neutrinoSolution) const {
|
| 704 |
return getLeptonicChi2(neutrinoSolution) + getHadronicChi2() + getGlobalChi2(neutrinoSolution);
|
| 705 |
}
|
| 706 |
|
| 707 |
double TopPairEventCandidate::getLeptonicChi2(unsigned short int neutrinoSolution) const {
|
| 708 |
double topMass(0);
|
| 709 |
double angle = leptonicBJet->angle(leptonFromW);
|
| 710 |
if (neutrinoSolution == 1)
|
| 711 |
topMass = leptonicTop1->mass();
|
| 712 |
else
|
| 713 |
topMass = leptonicTop2->mass();
|
| 714 |
|
| 715 |
return getLeptonicChi2(topMass, angle);
|
| 716 |
}
|
| 717 |
|
| 718 |
double TopPairEventCandidate::getLeptonicChi2(double topMass, double angle) const {
|
| 719 |
cout << "Using old TTBar getLeptonicChi2 \n";
|
| 720 |
double massDifference = TMath::Power(topMass - matched_leptonic_top_mass, 2);
|
| 721 |
double massError = 2 * matched_leptonic_top_mass_sigma * matched_leptonic_top_mass_sigma;
|
| 722 |
double massTerm = massDifference / massError;
|
| 723 |
|
| 724 |
double angleDifference = TMath::Power(angle - matched_angle, 2);
|
| 725 |
double angleError = 2 * matched_angle_sigma * matched_angle_sigma;
|
| 726 |
double angleTerm = angleDifference / angleError;
|
| 727 |
return 1.0 / sqrt(2.0) * (angleTerm + massTerm);
|
| 728 |
}
|
| 729 |
|
| 730 |
double TopPairEventCandidate::getHadronicChi2() const {
|
| 731 |
double ptRatioDifference = TMath::Power(PtRatio() - matched_ptratio, 2);
|
| 732 |
double ptRatioError = 2 * matched_ptratio_sigma * matched_ptratio_sigma;
|
| 733 |
double ptRatioTerm = ptRatioDifference / ptRatioError;
|
| 734 |
|
| 735 |
double WmassDifference = TMath::Power(hadronicW->mass() - matched_hadronic_W_mass, 2);
|
| 736 |
double WmassError = 2 * matched_hadronic_W_mass_sigma * matched_hadronic_W_mass_sigma;
|
| 737 |
double WmassTerm = WmassDifference / WmassError;
|
| 738 |
|
| 739 |
double topMassDifference = TMath::Power(hadronicTop->mass() - matched_hadronic_top_mass, 2);
|
| 740 |
double topMassError = 2 * matched_hadronic_top_mass_sigma * matched_hadronic_top_mass_sigma;
|
| 741 |
double topMassTerm = topMassDifference / topMassError;
|
| 742 |
return 1 / sqrt(3) * (topMassTerm + WmassTerm + ptRatioTerm);
|
| 743 |
|
| 744 |
///////////////////////////// Latest Bristol code may have dropped ptRatioTerm
|
| 745 |
// return 1 / sqrt(2) * (topMassTerm + WmassTerm);// + ptRatioTerm);
|
| 746 |
//
|
| 747 |
}
|
| 748 |
|
| 749 |
double TopPairEventCandidate::PtRatio() const {
|
| 750 |
return TMath::Log(hadronicTop->pt() / hadronicW->pt());
|
| 751 |
}
|
| 752 |
|
| 753 |
double TopPairEventCandidate::getGlobalChi2(unsigned short neutrinoSolution) const {
|
| 754 |
double pttbar = PtTtbarSystem(neutrinoSolution);
|
| 755 |
double pttbarDifference = TMath::Power(pttbar - matched_pt_ttbarSystem, 2);
|
| 756 |
double pttbarError = (2 * matched_pt_ttbarSystem_sigma * matched_pt_ttbarSystem_sigma);
|
| 757 |
double pttbarTerm = pttbarDifference / pttbarError;
|
| 758 |
|
| 759 |
double htSystemDifference = TMath::Power(HTSystem() - matched_HTSystem, 2);
|
| 760 |
double htSystemError = matched_HTSystem_sigma * matched_HTSystem_sigma * 2;
|
| 761 |
double htSystemTerm = htSystemDifference / htSystemError;
|
| 762 |
return 1 / sqrt(2) * (pttbarTerm + htSystemTerm);
|
| 763 |
}
|
| 764 |
|
| 765 |
double TopPairEventCandidate::PtTtbarSystem(unsigned short neutrinoSolution) const {
|
| 766 |
ParticlePointer combined;
|
| 767 |
if (neutrinoSolution == 1)
|
| 768 |
combined = ParticlePointer(new Particle(*leptonicTop1 + *hadronicTop));
|
| 769 |
else
|
| 770 |
combined = ParticlePointer(new Particle(*leptonicTop2 + *hadronicTop));
|
| 771 |
return combined->pt();
|
| 772 |
}
|
| 773 |
|
| 774 |
double TopPairEventCandidate::HT(unsigned short jetLimit) const {
|
| 775 |
double HT(0);
|
| 776 |
unsigned short limit = goodJets.size();
|
| 777 |
if (limit > jetLimit + 1)
|
| 778 |
limit = jetLimit + 1;
|
| 779 |
|
| 780 |
for (unsigned short index = 0; index < limit; ++index)
|
| 781 |
HT += goodJets.at(index)->pt();
|
| 782 |
|
| 783 |
return HT;
|
| 784 |
}
|
| 785 |
|
| 786 |
double TopPairEventCandidate::HTSystem() const {
|
| 787 |
return sumPt() / HT(8);
|
| 788 |
}
|
| 789 |
|
| 790 |
double TopPairEventCandidate::sumPt() const {
|
| 791 |
return leptonicBJet->pt() + hadronicBJet->pt() + jet1FromW->pt() + jet2FromW->pt();
|
| 792 |
}
|
| 793 |
|
| 794 |
void TopPairEventCandidate::throwExpeptionIfNotReconstructed(TString location) const {
|
| 795 |
if (doneReconstruction == false)
|
| 796 |
throw ReconstructionException(TString("Can't access reconstructed particles before reconstruction: ") + location);
|
| 797 |
}
|
| 798 |
|
| 799 |
const LeptonPointer TopPairEventCandidate::getLeptonFromWDecay() const {
|
| 800 |
return leptonFromW;
|
| 801 |
}
|
| 802 |
|
| 803 |
const ParticlePointer TopPairEventCandidate::getNeutrinoFromWDecay() const {
|
| 804 |
throwExpeptionIfNotReconstructed("TopPairEventCandidate::getNeutrinoFromWDecay");
|
| 805 |
if (selectedNeutrino == 1)
|
| 806 |
return neutrino1;
|
| 807 |
else
|
| 808 |
return neutrino2;
|
| 809 |
}
|
| 810 |
|
| 811 |
const JetPointer TopPairEventCandidate::getHadronicBJet() const {
|
| 812 |
throwExpeptionIfNotReconstructed("TopPairEventCandidate::getHadronicBJet");
|
| 813 |
return hadronicBJet;
|
| 814 |
}
|
| 815 |
|
| 816 |
const JetPointer TopPairEventCandidate::getLeptonicBJet() const {
|
| 817 |
throwExpeptionIfNotReconstructed("TopPairEventCandidate::getLeptonicBJet");
|
| 818 |
return leptonicBJet;
|
| 819 |
}
|
| 820 |
|
| 821 |
|
| 822 |
const ParticlePointer TopPairEventCandidate::getLeptonicW() const
|
| 823 |
{
|
| 824 |
throwExpeptionIfNotReconstructed("TopPairEventCandidate::getLeptonicW");
|
| 825 |
if (selectedNeutrino == 1)
|
| 826 |
return leptonicW1;
|
| 827 |
else
|
| 828 |
return leptonicW2;
|
| 829 |
}
|
| 830 |
|
| 831 |
|
| 832 |
const JetPointer TopPairEventCandidate::getJet1FromHadronicW() const {
|
| 833 |
throwExpeptionIfNotReconstructed("TopPairEventCandidate::getJet1FromHadronicW");
|
| 834 |
return jet1FromW;
|
| 835 |
}
|
| 836 |
|
| 837 |
const JetPointer TopPairEventCandidate::getJet2FromHadronicW() const {
|
| 838 |
throwExpeptionIfNotReconstructed("TopPairEventCandidate::getJet2FromHadronicW");
|
| 839 |
return jet2FromW;
|
| 840 |
}
|
| 841 |
|
| 842 |
const ParticlePointer TopPairEventCandidate::getHadronicW() const {
|
| 843 |
throwExpeptionIfNotReconstructed("TopPairEventCandidate::getHadronicW");
|
| 844 |
return hadronicW;
|
| 845 |
}
|
| 846 |
|
| 847 |
const ParticlePointer TopPairEventCandidate::getLeptonicTop() const {
|
| 848 |
throwExpeptionIfNotReconstructed("TopPairEventCandidate::getLeptonicTop");
|
| 849 |
if (selectedNeutrino == 1)
|
| 850 |
return leptonicTop1;
|
| 851 |
else
|
| 852 |
return leptonicTop2;
|
| 853 |
}
|
| 854 |
|
| 855 |
const ParticlePointer TopPairEventCandidate::getHadronicTop() const {
|
| 856 |
throwExpeptionIfNotReconstructed("TopPairEventCandidate::getHadronicTop");
|
| 857 |
return hadronicTop;
|
| 858 |
}
|
| 859 |
|
| 860 |
const ParticlePointer TopPairEventCandidate::getResonance() const {
|
| 861 |
throwExpeptionIfNotReconstructed("TopPairEventCandidate::getResonance");
|
| 862 |
return ttbarResonance;
|
| 863 |
}
|
| 864 |
|
| 865 |
double TopPairEventCandidate::M3() const {
|
| 866 |
double m3(0), max_pt(0);
|
| 867 |
if (goodJets.size() >= 3) {
|
| 868 |
for (unsigned int index1 = 0; index1 < goodJets.size() - 2; ++index1) {
|
| 869 |
for (unsigned int index2 = index1 + 1; index2 < goodJets.size() - 1; ++index2) {
|
| 870 |
for (unsigned int index3 = index2 + 1; index3 < goodJets.size(); ++index3) {
|
| 871 |
FourVector m3Vector(goodJets.at(index1)->getFourVector() + goodJets.at(index2)->getFourVector()
|
| 872 |
+ goodJets.at(index3)->getFourVector());
|
| 873 |
double currentPt = m3Vector.Pt();
|
| 874 |
if (currentPt > max_pt) {
|
| 875 |
max_pt = currentPt;
|
| 876 |
m3 = m3Vector.M();
|
| 877 |
}
|
| 878 |
}
|
| 879 |
}
|
| 880 |
}
|
| 881 |
}
|
| 882 |
|
| 883 |
return m3;
|
| 884 |
}
|
| 885 |
|
| 886 |
|
| 887 |
double TtbarHypothesis::M3() const {
|
| 888 |
double m3(0), max_et(0);
|
| 889 |
JetCollection mcJets;
|
| 890 |
mcJets.clear();
|
| 891 |
mcJets.push_back(jet1FromW);
|
| 892 |
mcJets.push_back(jet2FromW);
|
| 893 |
mcJets.push_back(leptonicBjet);
|
| 894 |
mcJets.push_back(hadronicBJet);
|
| 895 |
if (mcJets.size() >= 3) {
|
| 896 |
for (unsigned int index1 = 0; index1 < mcJets.size() - 2; ++index1) {
|
| 897 |
for (unsigned int index2 = index1 + 1; index2 < mcJets.size() - 1; ++index2) {
|
| 898 |
for (unsigned int index3 = index2 + 1; index3 < mcJets.size(); ++index3) {
|
| 899 |
FourVector m3Vector(mcJets.at(index1)->getFourVector() + mcJets.at(index2)->getFourVector()
|
| 900 |
+ mcJets.at(index3)->getFourVector());
|
| 901 |
double currentEt = m3Vector.Et();
|
| 902 |
if (currentEt > max_et) {
|
| 903 |
max_et = currentEt;
|
| 904 |
m3 = m3Vector.M();
|
| 905 |
}
|
| 906 |
}
|
| 907 |
}
|
| 908 |
}
|
| 909 |
}
|
| 910 |
|
| 911 |
return m3;
|
| 912 |
}
|
| 913 |
|
| 914 |
double TopPairEventCandidate::mttbar() const {
|
| 915 |
return getResonance()->mass();
|
| 916 |
}
|
| 917 |
|
| 918 |
const std::vector<TtbarHypothesisPointer>& TopPairEventCandidate::Solutions() const{
|
| 919 |
return solutions;
|
| 920 |
}
|
| 921 |
|
| 922 |
void TopPairEventCandidate::inspectReconstructedEvent() const {
|
| 923 |
cout << "run " << runNumber << ", event " << eventNumber << endl;
|
| 924 |
cout << "leptonic b jet, goodJet index " << leptonicBIndex << endl;
|
| 925 |
EventContentPrinter::printJet(leptonicBJet);
|
| 926 |
|
| 927 |
cout << "lepton from W" << endl;
|
| 928 |
// EventContentPrinter::printElectron(electronFromW);
|
| 929 |
|
| 930 |
cout << "MET" << endl;
|
| 931 |
EventContentPrinter::printParticle(met);
|
| 932 |
cout << endl;
|
| 933 |
|
| 934 |
cout << "reconstructed neutrino 1(selected: " << selectedNeutrino << ")" << endl;
|
| 935 |
EventContentPrinter::printParticle(neutrino1);
|
| 936 |
cout << endl;
|
| 937 |
|
| 938 |
cout << "reconstructed neutrino 2(selected: " << selectedNeutrino << ")" << endl;
|
| 939 |
EventContentPrinter::printParticle(neutrino2);
|
| 940 |
cout << endl;
|
| 941 |
|
| 942 |
cout << "leptonic W 1 (selected: " << selectedNeutrino << ")" << endl;
|
| 943 |
EventContentPrinter::printParticle(leptonicW1);
|
| 944 |
cout << endl;
|
| 945 |
|
| 946 |
cout << "leptonic W 2 (selected: " << selectedNeutrino << ")" << endl;
|
| 947 |
EventContentPrinter::printParticle(leptonicW2);
|
| 948 |
cout << endl;
|
| 949 |
|
| 950 |
cout << "leptonic top 1 (selected: " << selectedNeutrino << ")" << endl;
|
| 951 |
EventContentPrinter::printParticle(leptonicTop1);
|
| 952 |
cout << endl;
|
| 953 |
|
| 954 |
cout << "leptonic top 2 (selected: " << selectedNeutrino << ")" << endl;
|
| 955 |
EventContentPrinter::printParticle(leptonicTop2);
|
| 956 |
cout << endl;
|
| 957 |
|
| 958 |
cout << "hadronic b jet, goodJet index " << hadronicBIndex << endl;
|
| 959 |
EventContentPrinter::printJet(hadronicBJet);
|
| 960 |
|
| 961 |
cout << "jet1 from W, goodJet index " << jet1FromWIndex << endl;
|
| 962 |
EventContentPrinter::printJet(jet1FromW);
|
| 963 |
|
| 964 |
cout << "jet 2 from W, goodJet index " << jet2FromWIndex << endl;
|
| 965 |
EventContentPrinter::printJet(jet2FromW);
|
| 966 |
|
| 967 |
cout << "hadronic W" << endl;
|
| 968 |
EventContentPrinter::printParticle(hadronicW);
|
| 969 |
cout << endl;
|
| 970 |
|
| 971 |
cout << "hadronic top" << endl;
|
| 972 |
EventContentPrinter::printParticle(hadronicTop);
|
| 973 |
cout << endl;
|
| 974 |
}
|
| 975 |
|
| 976 |
double TopPairEventCandidate::fullHT() const {
|
| 977 |
double ht(met->pt());
|
| 978 |
|
| 979 |
for (unsigned int index = 0; index < goodPFIsolatedElectrons.size(); ++index) {
|
| 980 |
ht += goodPFIsolatedElectrons.at(index)->pt();
|
| 981 |
}
|
| 982 |
|
| 983 |
for (unsigned int index = 0; index < goodIsolatedMuons.size(); ++index) {
|
| 984 |
ht += goodIsolatedMuons.at(index)->pt();
|
| 985 |
}
|
| 986 |
|
| 987 |
for (unsigned int index = 0; index < goodJets.size(); ++index) {
|
| 988 |
ht += goodJets.at(index)->pt();
|
| 989 |
}
|
| 990 |
return ht;
|
| 991 |
}
|
| 992 |
|
| 993 |
double TopPairEventCandidate::transverseWmass(const LeptonPointer lepton) const {
|
| 994 |
double energySquared = pow(lepton->et() + met->et(), 2);
|
| 995 |
double momentumSquared = pow(lepton->px() + met->px(), 2) + pow(lepton->py() + met->py(), 2);
|
| 996 |
double tMassSquared = energySquared - momentumSquared;
|
| 997 |
|
| 998 |
if (tMassSquared > 0)
|
| 999 |
return sqrt(tMassSquared);
|
| 1000 |
else
|
| 1001 |
return -1;
|
| 1002 |
}
|
| 1003 |
|
| 1004 |
void TopPairEventCandidate::reconstruct(Rule rule, const LeptonPointer lepton) {
|
| 1005 |
//TODO: change this into if(!Rule::meetsInitialConditaion())
|
| 1006 |
if (goodJets.size() < 4 || lepton == 0)
|
| 1007 |
throw ReconstructionException("Not enough jets available to reconstruct top event using Mass Equality method.");
|
| 1008 |
leptonFromW = lepton;
|
| 1009 |
selectedNeutrino = 0;
|
| 1010 |
currentSelectedNeutrino = 0;
|
| 1011 |
|
| 1012 |
reconstructNeutrinos();
|
| 1013 |
double chosen_TopMassDifference(9999999.);
|
| 1014 |
double chosen_Chi2Total(9999999.);
|
| 1015 |
for (unsigned short hadBindex = 0; hadBindex < goodJets.size(); ++hadBindex) {
|
| 1016 |
if (rule->hardHadronicBJetCondition(goodJets.at(hadBindex)))
|
| 1017 |
continue;
|
| 1018 |
|
| 1019 |
for (unsigned short lepBindex = 0; lepBindex < goodJets.size(); ++lepBindex) {
|
| 1020 |
if (lepBindex == hadBindex || rule->hardLeptonicBJetCondition(goodJets.at(lepBindex)))
|
| 1021 |
continue;
|
| 1022 |
|
| 1023 |
for (unsigned short jet1Index = 0; jet1Index < goodJets.size(); ++jet1Index) {
|
| 1024 |
if (jet1Index == lepBindex || jet1Index == hadBindex || rule->hardHadronicJetFromWCondition(
|
| 1025 |
goodJets.at(jet1Index)))
|
| 1026 |
continue;
|
| 1027 |
|
| 1028 |
for (unsigned short jet2Index = 0; jet2Index < goodJets.size(); ++jet2Index) {
|
| 1029 |
if (jet2Index == jet1Index || jet2Index == lepBindex || jet2Index == hadBindex
|
| 1030 |
|| rule->hardHadronicJetFromWCondition(goodJets.at(jet2Index)))
|
| 1031 |
continue;
|
| 1032 |
|
| 1033 |
hadronicBJet = goodJets.at(hadBindex);
|
| 1034 |
leptonicBJet = goodJets.at(lepBindex);
|
| 1035 |
jet1FromW = goodJets.at(jet1Index);
|
| 1036 |
jet2FromW = goodJets.at(jet2Index);
|
| 1037 |
|
| 1038 |
leptonicW1 = ParticlePointer(new Particle(*neutrino1 + *leptonFromW));
|
| 1039 |
leptonicW2 = ParticlePointer(new Particle(*neutrino2 + *leptonFromW));
|
| 1040 |
hadronicW = ParticlePointer(new Particle(*jet1FromW + *jet2FromW));
|
| 1041 |
leptonicTop1 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW1));
|
| 1042 |
leptonicTop2 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW2));
|
| 1043 |
hadronicTop = ParticlePointer(new Particle(*hadronicBJet + *hadronicW));
|
| 1044 |
fillHypotheses();
|
| 1045 |
selectNeutrinoSolution();
|
| 1046 |
double TopMassDifference = calculateTopMassDifference(currentSelectedNeutrino);
|
| 1047 |
double chi2 = getTotalChi2(currentSelectedNeutrino);
|
| 1048 |
switch (usedTTbarReconstruction) {
|
| 1049 |
case TTbarReconstructionCriterion::TopMassDifference:
|
| 1050 |
if (TopMassDifference < chosen_TopMassDifference) {
|
| 1051 |
hadronicBIndex = hadBindex;
|
| 1052 |
leptonicBIndex = lepBindex;
|
| 1053 |
jet1FromWIndex = jet1Index;
|
| 1054 |
jet2FromWIndex = jet2Index;
|
| 1055 |
chosen_TopMassDifference = TopMassDifference;
|
| 1056 |
selectedNeutrino = currentSelectedNeutrino;
|
| 1057 |
}
|
| 1058 |
break;
|
| 1059 |
|
| 1060 |
case TTbarReconstructionCriterion::chi2:
|
| 1061 |
if (chi2 < chosen_Chi2Total) {
|
| 1062 |
hadronicBIndex = hadBindex;
|
| 1063 |
leptonicBIndex = lepBindex;
|
| 1064 |
jet1FromWIndex = jet1Index;
|
| 1065 |
jet2FromWIndex = jet2Index;
|
| 1066 |
chosen_Chi2Total = chi2;
|
| 1067 |
selectedNeutrino = currentSelectedNeutrino;
|
| 1068 |
}
|
| 1069 |
break;
|
| 1070 |
}
|
| 1071 |
}
|
| 1072 |
}
|
| 1073 |
}
|
| 1074 |
}
|
| 1075 |
compare_disc sorting;
|
| 1076 |
std::sort(solutions.begin(), solutions.end(), sorting);
|
| 1077 |
hadronicBJet = solutions.front()->hadronicBJet;
|
| 1078 |
leptonicBJet = solutions.front()->leptonicBjet;
|
| 1079 |
jet1FromW = solutions.front()->jet1FromW;
|
| 1080 |
jet2FromW = solutions.front()->jet2FromW;
|
| 1081 |
leptonicW1 = solutions.front()->leptonicW;
|
| 1082 |
leptonicW2 = solutions.front()->leptonicW;
|
| 1083 |
hadronicW = solutions.front()->hadronicW;
|
| 1084 |
// leptonicW1 = ParticlePointer(new Particle(*neutrino1 + *leptonFromW));
|
| 1085 |
// leptonicW2 = ParticlePointer(new Particle(*neutrino2 + *leptonFromW));
|
| 1086 |
// hadronicW = ParticlePointer(new Particle(*jet1FromW + *jet2FromW));
|
| 1087 |
leptonicTop1 = solutions.front()->leptonicTop;
|
| 1088 |
leptonicTop2 = solutions.front()->leptonicTop;
|
| 1089 |
hadronicTop = solutions.front()->hadronicTop;
|
| 1090 |
ttbarResonance = solutions.front()->resonance;
|
| 1091 |
if (selectedNeutrino == 1)
|
| 1092 |
ttbarResonance = ParticlePointer(new Particle(*leptonicTop1 + *hadronicTop));
|
| 1093 |
else
|
| 1094 |
ttbarResonance = ParticlePointer(new Particle(*leptonicTop2 + *hadronicTop));
|
| 1095 |
doneReconstruction = true;
|
| 1096 |
}
|
| 1097 |
|
| 1098 |
}
|
