1 |
dhidas |
1.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 |
|
|
}
|