ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/dhidas/OSUAnalysis/Tools/src/TopPairEventCandidate.cpp
Revision: 1.1
Committed: Thu Dec 1 16:28:48 2011 UTC (13 years, 5 months ago) by dhidas
Branch: MAIN
Branch point for: dhidas
Log Message:
Initial revision

File Contents

# User Rev Content
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     }