ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/dhidas/OSUAnalysis/Tools/src/TopPairEventCandidate.cpp
Revision: 1.2
Committed: Tue Dec 6 13:56:15 2011 UTC (13 years, 5 months ago) by clseitz
Branch: MAIN
Changes since 1.1: +6 -0 lines
Log Message:
*** empty log message ***

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