ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/dhidas/OSUAnalysis/Tools/src/ToplikeCandidate.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 #include <iostream>
2     #include <iomanip>
3    
4     #include "../interface/ToplikeCandidate.h"
5    
6    
7     using namespace std;
8     using namespace BAT;
9    
10     namespace BAT {
11    
12     static const double pi = 3.1415927;
13     static const int DQ_ID = 1, BQ_ID = 5, TOP_ID = 6, W_ID = 24, ELECTRON_ID = 11;
14     static const int TAU_ID = 15, NU_E_ID = 12, MU_ID = 13, NU_MU_ID = 14;
15     static const int WPRIME_ID = 4000023;
16     static const unsigned int MCLIMIT = 90; // Number of MC particles to check
17     static const float DELTA_R_LIMIT = 0.5; // Biggest deltaR to consider in range
18     static const double pt_tprime_sys_mean = 48.49;
19     static const double pt_tprime_sys_sigma = 27.87;
20    
21     static map <ToplikeSelectionSteps::Step, TTbarEPlusJetsSelection::Step> selSteps;
22    
23    
24     ToplikeCandidate::ToplikeCandidate() :
25     jet3FromWIndex(0),
26     jet4FromWIndex(0),
27     doneReconstructiontop(false)
28     {
29     }
30    
31    
32     ToplikeCandidate::ToplikeCandidate(const Event& event) :
33     TopPairEventCandidate(event),
34     jet3FromWIndex(0),
35     jet4FromWIndex(0),
36     doneReconstructiontop(false)
37     {
38     static bool mapset = false;
39     if (mapset == false) {
40     for(unsigned int ind = 0; ind < sizeof(selStepArr)/sizeof(selPair); ++ind) {
41     selSteps[selStepArr[ind].key] = selStepArr[ind].value;
42     }
43     mapset= true;
44     }
45     }
46    
47    
48     ToplikeCandidate::~ToplikeCandidate()
49     {
50     }
51    
52    
53     bool ToplikeCandidate::hasAtLeastFiveGoodJets() const {
54     return goodJets.size() >= 5;
55     }
56    
57     /*
58     const ParticlePointer ToplikeCandidate::getResonance() const {
59     throwExpeptionIfNotReconstructed("ToplikeCandidate::getResonance");
60     return tPrime;
61     }
62     */
63    
64    
65     double ToplikeCandidate::tpmass() const {
66     /*
67     ParticlePointer tmp = ParticlePointer(new Particle(*leptonicTop1 + *hadronicW));
68     cout << "tmp px = " << tmp->px() << " tp py = " << tmp->py() << " tp pz = " << tmp->pz() << endl;
69     cout << "energy = " << tmp->energy() << endl;
70     return tmp->mass();
71     cout << "tp px = " << tPrime->px() << " tp py = " << tPrime->py() << " tp pz = " << tPrime->pz() << endl;
72     cout << "energy = " << tPrime->energy() << endl;
73     if (tPrime->px() <= 0)
74     tPrime->setMass(0.0);
75     */
76     return tPrime->mass();
77     }
78    
79    
80     double ToplikeCandidate::sumPt() const {
81     return (jet1FromW->pt() + jet2FromW->pt() + hadronicBJet->pt());
82     }
83    
84    
85     double ToplikeCandidate::getGlobalChi2() const {
86     double pttbar = PtTtbarSystem();
87     double pttbarTerm = TMath::Power(pttbar - pt_tprime_sys_mean, 2) / (2 * pt_tprime_sys_sigma * pt_tprime_sys_sigma);
88    
89     double htSystemDifference = TMath::Power(HTSystem() - matched_HTSystem, 2);
90     double htSystemError = matched_HTSystem_sigma * matched_HTSystem_sigma * 2;
91     double htSystemTerm = htSystemDifference / htSystemError;
92     return 1 / sqrt(2) * (pttbarTerm + htSystemTerm);
93     }
94    
95    
96     double ToplikeCandidate::PtTtbarSystem() const {
97     ParticlePointer combined;
98     combined = ParticlePointer(new Particle(*hadronicWtPrime + *hadronicTop));
99     return combined->pt() / HTSystem();
100     }
101    
102    
103     double ToplikeCandidate::TPrimeHTSystem() const {
104     return htSystem;
105     }
106    
107    
108     double ToplikeCandidate::PtTPrimeSystem() const {
109     return ptTprimeSystem;
110     }
111    
112    
113     double ToplikeCandidate::PtTtbarSystem(unsigned short neutrinoSolution) const {
114     ParticlePointer combined;
115     if (neutrinoSolution == 1)
116     combined = ParticlePointer(new Particle(*leptonicTop1 + *hadronicW));
117     else
118     combined = ParticlePointer(new Particle(*leptonicTop2 + *hadronicW));
119     return combined->pt() / HTSystem();
120     }
121    
122    
123     double ToplikeCandidate::getWChi2() const {
124     double WmassDifference = TMath::Power(hadronicWtPrime->mass() - pdgWmass, 2);
125     double WmassError = 2 * wMassWidth * wMassWidth;
126     double WmassTerm = WmassDifference / WmassError;
127     return (WmassTerm);
128     }
129    
130    
131     double ToplikeCandidate::getLoneHadChi2() const
132     {
133     return (TopPairEventCandidate::getHadronicChi2());
134     }
135    
136    
137     /*
138     double ToplikeCandidate::getTotalChi2(unsigned short int neutrinoSolution) const {
139     return getLeptonicChi2(neutrinoSolution) + getHadronicChi2() + getGlobalChi2(neutrinoSolution);
140     }
141     */
142    
143    
144     void ToplikeCandidate::recoTprimeUsingChi2(ElectronPointer electron) {
145     // cout << "Starting tp reco " << endl;
146     if (goodJets.size() < 3)
147     throw ReconstructionException("Not enough jets available to reconstruct Tprime using Chi2 method.");
148     double chosen_Chi2Total(9999999.);
149     leptonFromW = electron;
150     selectedNeutrino = 0;
151     currentSelectedNeutrino = 0;
152     reconstructNeutrinos();
153    
154     // cout << "Starting jet loop" << endl;
155     for (unsigned short hadBindex = 0; hadBindex < goodJets.size(); ++hadBindex) {
156     for (unsigned short jet1Index = 0; jet1Index < goodJets.size(); ++jet1Index) {
157     if (jet1Index == hadBindex)
158     continue;
159     for (unsigned short jet2Index = 0; jet2Index < goodJets.size(); ++jet2Index) {
160     if (jet2Index == jet1Index || jet2Index == hadBindex)
161     continue;
162     hadronicBJet = goodJets.at(hadBindex);
163     leptonicBJet = hadronicBJet; // hack
164     jet1FromW = goodJets.at(jet1Index);
165     jet2FromW = goodJets.at(jet2Index);
166     leptonicW1 = ParticlePointer(new Particle(*neutrino1 + *leptonFromW));
167     leptonicW2 = ParticlePointer(new Particle(*neutrino2 + *leptonFromW));
168     leptonicTop1 = ParticlePointer(new Particle(*hadronicBJet + *leptonicW1));
169     leptonicTop2 = ParticlePointer(new Particle(*hadronicBJet + *leptonicW2));
170     hadronicW = ParticlePointer(new Particle(*jet1FromW + *jet2FromW));
171     // fillHypotheses();
172     selectNeutrinoSolution();
173     double chi2 = TopPairEventCandidate::getTotalChi2(currentSelectedNeutrino);
174     if (chi2 < chosen_Chi2Total) {
175     hadronicBIndex = hadBindex;
176     jet1FromWIndex = jet1Index;
177     jet2FromWIndex = jet2Index;
178     chosen_Chi2Total = chi2;
179     selectedNeutrino = currentSelectedNeutrino;
180     }
181     }
182     }
183     }
184     // cout << "Done with loop" << endl;
185     hadronicBJet = goodJets.at(hadronicBIndex);
186     leptonicBJet = hadronicBJet; // hack
187     jet1FromW = goodJets.at(jet1FromWIndex);
188     jet2FromW = goodJets.at(jet2FromWIndex);
189     hadronicW = ParticlePointer(new Particle(*jet1FromW + *jet2FromW));
190     leptonicW1 = ParticlePointer(new Particle(*neutrino1 + *leptonFromW));
191     leptonicW2 = ParticlePointer(new Particle(*neutrino2 + *leptonFromW));
192     leptonicTop1 = ParticlePointer(new Particle(*hadronicBJet + *leptonicW1));
193     leptonicTop2 = ParticlePointer(new Particle(*hadronicBJet + *leptonicW2));
194     ptTprimeSystem = PtTtbarSystem(selectedNeutrino);
195     htSystem = HTSystem();
196     if (selectedNeutrino == 1)
197     tPrime = ParticlePointer(new Particle(*leptonicTop1 + *hadronicW));
198     else
199     tPrime = ParticlePointer(new Particle(*leptonicTop2 + *hadronicW));
200     // cout << "tPrime mass = " << tPrime->mass() << endl;
201     doneReconstruction = true;
202     }
203    
204    
205     // Omit lepton-jet angle from criteria.
206     double ToplikeCandidate::getLeptonicChi2(double topMass, double unused)
207     const
208     {
209     double massDifference = TMath::Power(topMass - pdgTopMass, 2);
210     double massError = 2 * topMassWidth * topMassWidth;
211     double massTerm = massDifference / massError;
212     return (massTerm);
213     }
214    
215    
216     double ToplikeCandidate::getHadronicChi2() const {
217     double WmassDifference = TMath::Power(hadronicW->mass() - pdgWmass, 2);
218     double WmassError = 2 * wMassWidth * wMassWidth;
219     double WmassTerm = WmassDifference / WmassError;
220    
221     double topMassDifference = TMath::Power(hadronicTop->mass() - pdgTopMass, 2);
222     double topMassError = 2 * topMassWidth * topMassWidth;
223     double topMassTerm = topMassDifference / topMassError;
224     return ((1.0 / sqrt(2.0)) * (topMassTerm + WmassTerm));
225     }
226    
227    
228     void ToplikeCandidate::recoBestSingleTop(LeptonPointer lepton) {
229     if (goodJets.size() < 3)
230     throw ReconstructionException("Not enough jets available to reconstruct top using Chi2 method.");
231     double bestLepChi2(9999999.), bestHadChi2(9999999.);
232     leptonFromW = lepton;
233     selectedNeutrino = 0;
234     currentSelectedNeutrino = 0;
235     reconstructNeutrinos();
236    
237     for (unsigned short hadBindex = 0; hadBindex < goodJets.size(); ++hadBindex) {
238     for (unsigned short jet1Index = 0; jet1Index < goodJets.size(); ++jet1Index) {
239     if (hadBindex == jet1Index)
240     continue;
241     for (unsigned short jet2Index = 0; jet2Index < goodJets.size(); ++jet2Index) {
242     if (jet2Index == jet1Index || jet2Index == hadBindex)
243     continue;
244     hadronicBJet = goodJets.at(hadBindex);
245     leptonicBJet = hadronicBJet; // hack
246     jet1FromW = goodJets.at(jet1Index);
247     jet2FromW = goodJets.at(jet2Index);
248     leptonicW1 = ParticlePointer(new Particle(*neutrino1 + *leptonFromW));
249     leptonicW2 = ParticlePointer(new Particle(*neutrino2 + *leptonFromW));
250     leptonicTop1 = ParticlePointer(new Particle(*hadronicBJet + *leptonicW1));
251     leptonicTop2 = ParticlePointer(new Particle(*hadronicBJet + *leptonicW2));
252     hadronicW = ParticlePointer(new Particle(*jet1FromW + *jet2FromW));
253     hadronicTop = ParticlePointer(new Particle(*hadronicBJet + *hadronicW));
254     // fillHypotheses();
255     selectNeutrinoSolution();
256     double ltopchi2 = TopPairEventCandidate::getLeptonicChi2(currentSelectedNeutrino);
257     double htopchi2 = getHadronicChi2();
258     if (ltopchi2 < bestLepChi2) {
259     selectedNeutrino = currentSelectedNeutrino;
260     leptonicBIndex = hadBindex;
261     bestLepChi2 = ltopchi2;
262     }
263     if (htopchi2 < bestHadChi2) {
264     hadronicBIndex = hadBindex;
265     jet1FromWIndex = jet1Index;
266     jet2FromWIndex = jet2Index;
267     bestHadChi2 = htopchi2;
268     }
269     }
270     }
271     }
272     hadronicBJet = goodJets.at(hadronicBIndex);
273     leptonicBJet = goodJets.at(leptonicBIndex);
274     jet1FromW = goodJets.at(jet1FromWIndex);
275     jet2FromW = goodJets.at(jet2FromWIndex);
276     hadronicW = ParticlePointer(new Particle(*jet1FromW + *jet2FromW));
277     leptonicW1 = ParticlePointer(new Particle(*neutrino1 + *leptonFromW));
278     leptonicW2 = ParticlePointer(new Particle(*neutrino2 + *leptonFromW));
279     hadronicTop = ParticlePointer(new Particle(*hadronicBJet + *hadronicW));
280     leptonicTop1 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW1));
281     leptonicTop2 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW2));
282     // cout << "Found lone tops " << endl;
283     doneReconstructiontop = true;
284     doneReconstruction = true;
285     }
286    
287    
288    
289     static bool cmpMcEntries(const mcObj &firstArg, const mcObj &secondArg)
290     {
291     return (firstArg.deltaR < secondArg.deltaR);
292     }
293    
294    
295     mcList ToplikeCandidate::getMCList(const Particle *const recoObject) const
296     {
297     const MCParticleCollection &mcColl = genParticles;
298     mcList mcObjs;
299     if (mcColl.size() > 0) {
300     for (unsigned int i = 0; i < mcColl.size(); ++i) {
301     mcObj mcEntry;
302     mcEntry.deltaR = recoObject->deltaR(mcColl.at(i));
303     mcEntry.ptr = mcColl.at(i);
304     mcEntry.gpInd = i;
305     mcObjs.push_back(mcEntry);
306     }
307     sort(mcObjs.begin(), mcObjs.end(), cmpMcEntries);
308     }
309     return (mcObjs);
310     }
311    
312    
313     mcList ToplikeCandidate::getMCListForJets(const Jet *const recoObject)
314     const
315     {
316     const JetCollection &mcColl = genJets;
317     mcList mcObjs;
318     if (mcColl.size() > 0) {
319     for (unsigned int i = 0; i < mcColl.size(); ++i) {
320     mcObj mcEntry;
321     mcEntry.deltaR = recoObject->deltaR(mcColl.at(i));
322     mcEntry.gpInd = i;
323     mcObjs.push_back(mcEntry);
324     }
325     sort(mcObjs.begin(), mcObjs.end(), cmpMcEntries);
326     }
327     return (mcObjs);
328     }
329    
330    
331     int ToplikeCandidate::missRecoJetsForGenJets() const
332     {
333     const JetCollection &mcColl = genJets;
334     int numGenJets = mcColl.size();
335     int matchJets = 0;
336     for (int i = 0; i < numGenJets; ++i) {
337     bool match = false;
338     for (unsigned int recoInd = 0;
339     recoInd < GoodJets().size() && match == false; ++recoInd)
340     if (GoodJets().at(recoInd)->deltaR(mcColl.at(i)) < 0.3) {
341     ++matchJets;
342     match = true;
343     }
344     }
345     return (numGenJets - matchJets);
346     }
347    
348    
349    
350     topTruth ToplikeCandidate::getMCMatches(const ParticlePointer &recoObject)
351     const
352     {
353     const float nullDeltaR = 3;
354     bool deltaRSet = false;
355     topTruth result = {nullDeltaR, false, false};
356     mcList mcObjs = getMCList(recoObject.get());
357     if (mcObjs.size() <= 0)
358     return (result);
359     /*
360     if (mcObjs.begin()->first > 0.9) // Nothing close enough
361     return ("");
362     TString idList;
363     idList += mcObjs.at(0).first;
364     idList += " ";
365     for (unsigned int i = 0; i < mcObjs.size() && i < 12; ++i) {
366     idList += PdgStrs::pdgStr(mcObjs.at(i).second);
367     idList += " ";
368     }
369     */
370     float topDeltaR = nullDeltaR;
371     bool wFound = false, lepFound = false;
372     for (unsigned int i = 0; i < mcObjs.size() && i < MCLIMIT &&
373     (deltaRSet == false || lepFound == false || wFound == false); ++i) {
374     if (deltaRSet == false && fabs(mcObjs.at(i).ptr->pdgId()) == TOP_ID) {
375     topDeltaR = mcObjs.at(i).deltaR;
376     deltaRSet = true;
377     if (topDeltaR > nullDeltaR)
378     topDeltaR = nullDeltaR;
379     }
380     if (fabs(mcObjs.at(i).ptr->pdgId()) == ELECTRON_ID ||
381     fabs(mcObjs.at(i).ptr->pdgId()) == MU_ID)
382     lepFound = true;
383     if (fabs(mcObjs.at(i).ptr->pdgId()) == W_ID)
384     wFound = true;
385     }
386     if (deltaRSet == false)
387     cout << " No top found in GPs! Increase MCLIMIT.\n" ;
388     result.deltaR = topDeltaR;
389     result.leptonic = wFound && lepFound;
390     result.hadronic = wFound;
391     return (result);
392     }
393    
394    
395     static int findMCParticle(const mcList &mcObjs, int pdgId, int antiId = 0,
396     unsigned int startInd = 0)
397     {
398     if (antiId == 0)
399     antiId = pdgId;
400     for (unsigned int i = startInd; i < mcObjs.size(); ++i) {
401     if (pdgId == 0 || mcObjs.at(i).ptr->pdgId() == pdgId ||
402     mcObjs.at(i).ptr->pdgId() == antiId) {
403     if (mcObjs.at(i).deltaR < DELTA_R_LIMIT)
404     return (i);
405     else return (-1);
406     }
407     }
408     return (-1);
409     }
410    
411    
412     // Scans up the decay chain when mother & daughter are the same.
413    
414     bool ToplikeCandidate::chkIfChild(int gpInd, int motherGpInd) const
415     {
416     if (motherGpInd < 0)
417     return (false);
418     const MCParticleCollection &mcColl = genParticles;
419     if (mcColl.size() <= 0)
420     return (false);
421     int mother = mcColl.at(gpInd)->motherIndex();
422     if (mother == motherGpInd)
423     return (true);
424     while (mother > 0 &&
425     mcColl.at(gpInd)->pdgId() == mcColl.at(mother)->pdgId())
426     {
427     gpInd = mother;
428     mother = mcColl.at(gpInd)->motherIndex();
429     if (mother == motherGpInd)
430     return (true);
431     }
432     int gMother = mcColl.at(motherGpInd)->motherIndex();
433    
434     // W' sample b's have W' as mother rather than top.
435     if (gMother > 0 &&
436     fabs(mcColl.at(gMother)->pdgId()) == WPRIME_ID && mother == gMother)
437     return (true);
438     while (gMother > 0 &&
439     mcColl.at(motherGpInd)->pdgId() == mcColl.at(gMother)->pdgId())
440     {
441     motherGpInd = gMother;
442     gMother = mcColl.at(motherGpInd)->motherIndex();
443     if (mother == motherGpInd)
444     return (true);
445     }
446     /*
447     if (mcColl.at(mother)->motherIndex() == mcColl.at(motherGpInd)->motherIndex())
448     {
449     cout << "gm match: ind cpdg mi mpdg gm pdg omi ompdg " << gpInd << " ";
450     cout << mcColl.at(gpInd)->pdgId() << " " << mother << " ";
451     cout << mcColl.at(mother)->pdgId() << " ";
452     int gm = mcColl.at(mother)->motherIndex();
453     cout << gm << " " << mcColl.at(gm)->pdgId() << " ";
454     cout << motherGpInd << " " << mcColl.at(motherGpInd)->pdgId() << endl;
455     return (true);
456     }
457     */
458     return (false);
459     }
460    
461    
462     int ToplikeCandidate::chkMCMatch(const Particle *const recoObject,
463     int motherInd, const int pdgId, int antiId) const
464     {
465     mcList truList = getMCList(recoObject);
466     if (truList.size() > 0) {
467     int index = -1;
468     while ((index = findMCParticle(truList, pdgId, antiId, index + 1)) >= 0)
469     if (motherInd == -1 || chkIfChild(truList.at(index).gpInd, motherInd))
470     return (truList.at(index).gpInd);
471     }
472     return (-1);
473     }
474    
475    
476     // Look for a non-top quark whose mother is specified by motherInd and which
477     // doesn't match badInd. If motherInd == -1, no mother is checked for.
478    
479     int ToplikeCandidate::chkForQuark(const JetPointer recoObject,
480     const int motherInd, const int badInd) const
481     {
482     mcList truList = getMCList(recoObject.get());
483     if (truList.size() > 0) {
484     for (int i = 0; i < (int) truList.size() &&
485     truList.at(i).deltaR < DELTA_R_LIMIT; ++i) {
486     if (truList.at(i).ptr->isQuark() && fabs(truList.at(i).ptr->pdgId()) !=
487     TOP_ID && (motherInd == -1 || chkIfChild(truList.at(i).gpInd, motherInd)) && truList.at(i).gpInd != badInd)
488     return (truList.at(i).gpInd);
489     // Checking for q-qbar pair not necessary since both daughters of mother
490     }
491     }
492     return (-1);
493     }
494    
495    
496     // Checks if ind1 is a b quark, other two are daughters of a W, and
497     // b & W are siblings.
498     bool ToplikeCandidate::chkBW(const int ind1, const int ind2, const int ind3)
499     const
500     {
501     const MCParticleCollection &mcColl = genParticles;
502     if (mcColl.size() <= 0)
503     return (false);
504     if (fabs(mcColl.at(ind1)->pdgId()) != BQ_ID)
505     return (false);
506     int jetMother = mcColl.at(ind2)->motherIndex();
507     if (jetMother > 0 && fabs(mcColl.at(jetMother)->pdgId()) == W_ID &&
508     mcColl.at(ind3)->motherIndex() == jetMother &&
509     mcColl.at(ind1)->motherIndex() == mcColl.at(jetMother)->motherIndex())
510     return (true);
511     return (false);
512     }
513    
514    
515     int ToplikeCandidate::goodHTopMatching() const
516     {
517     int numFound = 0, mostFound = 0;
518     mcList bTruList = getMCList(hadronicBJet.get());
519     mcList j1TruList = getMCList(jet1FromW.get());
520     mcList j2TruList = getMCList(jet2FromW.get());
521     int bInd = -1;
522     while (bTruList.size() > 0 &&
523     (bInd = findMCParticle(bTruList, 0, 0, bInd + 1)) >= 0)
524     {
525     if (bTruList.at(bInd).ptr->isQuark())
526     {
527     numFound = 1;
528     int bGPInd = bTruList.at(bInd).gpInd;
529     int j1Ind = -1;
530     while (j1TruList.size() > 0 &&
531     (j1Ind = findMCParticle(j1TruList, 0, 0, j1Ind + 1)) >= 0)
532     {
533     int j1GPInd = j1TruList.at(j1Ind).gpInd;
534     if (j1GPInd != bGPInd && j1TruList.at(j1Ind).ptr->isQuark())
535     {
536     numFound = 2;
537     int j2Ind = -1;
538     while (j2TruList.size() > 0 &&
539     (j2Ind = findMCParticle(j2TruList, 0, 0, j2Ind + 1)) >= 0)
540     {
541     int j2GPInd = j2TruList.at(j2Ind).gpInd;
542     if (j2GPInd != bGPInd && j2GPInd != j1GPInd &&
543     j2TruList.at(j2Ind).ptr->isQuark())
544     {
545     numFound = 3;
546     if (chkBW(bGPInd, j1GPInd, j2GPInd) ||
547     chkBW(j1GPInd, bGPInd, j2GPInd) ||
548     chkBW(j2GPInd, j1GPInd, bGPInd))
549     return (4);
550     }
551     }
552     if (numFound > mostFound)
553     mostFound = numFound;
554     }
555     }
556     }
557     }
558     if (numFound > mostFound)
559     mostFound = numFound;
560     return (mostFound);
561     }
562    
563    
564    
565     int ToplikeCandidate::strictHTopMatching(const mcObj &topTru) const
566     {
567     int numFound = 1; // Start with finding the top.
568     short antiMult = 1;
569     if (topTru.ptr->pdgId() < 0)
570     antiMult = -1;
571     // cout << " hti " << topTru.gpInd << " ";
572     if (chkMCMatch(hadronicBJet.get(), topTru.gpInd, BQ_ID * antiMult) >= 0)
573     ++numFound;
574     int wInd = chkMCMatch(hadronicW.get(), topTru.gpInd, W_ID * antiMult);
575     if (wInd >= 0)
576     ++numFound;
577     else return (numFound);
578    
579     // Try to find 2 quarks from W. If both fail, trying again will still fail
580     // but is harmless. If one fails, try again to see if can get both if we
581     // let the other jet grab the 1st quark. If that doesn't work, one quark
582     // will still be found.
583     bool pairFnd = false, last1 = false;
584     int takenInd2 = -1;
585     int takenInd = chkForQuark(jet1FromW, wInd);
586     if (takenInd >= 0 && chkForQuark(jet2FromW, wInd, takenInd) >= 0) {
587     numFound += 2;
588     pairFnd = true;
589     } else {
590     takenInd2 = chkForQuark(jet2FromW, wInd);
591     if (takenInd2 >= 0)
592     ++numFound;
593     if (chkForQuark(jet1FromW, wInd, takenInd2) >= 0) {
594     ++numFound;
595     last1 = true;
596     }
597     }
598     /*
599     if (numFound == 5 || numFound == 4) {
600     cout << "num pair wInd takenInd taken2 last1 " << numFound << " ";
601     cout << pairFnd << " " << wInd << " ";
602     cout << takenInd << " " << takenInd2 << " " << last1 << endl;
603     }
604     */
605     return (numFound);
606     }
607    
608    
609     int ToplikeCandidate::strictLTopMatching(const mcObj &topTru) const
610     {
611     int numFound = 1;
612     short antiMult = 1;
613     if (topTru.ptr->pdgId() < 0)
614     antiMult = -1;
615     if (chkMCMatch(getLeptonicBJet().get(), topTru.gpInd, BQ_ID * antiMult) >= 0)
616     ++numFound;
617     int wInd = chkMCMatch(getLeptonicW().get(), topTru.gpInd, W_ID * antiMult);
618     if (wInd >= 0)
619     ++numFound;
620     else return (numFound);
621     // Lepton is opposite from top to preserve lepton number.
622     if (chkMCMatch(getLeptonFromWDecay().get(), wInd, -ELECTRON_ID * antiMult)
623     >= 0 ||
624     chkMCMatch(getLeptonFromWDecay().get(), wInd, -MU_ID * antiMult) >= 0)
625     ++numFound;
626     if (chkMCMatch(getNeutrinoFromWDecay().get(), wInd, NU_E_ID * antiMult) >= 0
627     || chkMCMatch(getNeutrinoFromWDecay().get(), wInd, NU_MU_ID * antiMult)
628     >= 0)
629     ++numFound;
630     return (numFound);
631     }
632    
633    
634     int ToplikeCandidate::goodLTopMatching() const
635     {
636     const MCParticleCollection &mcColl = genParticles;
637     int numParticles = mcColl.size();
638     if (numParticles <= 0)
639     return (0);
640     int numFound = 0, mostFound = 0;
641     mcList lepTruList = getMCList(getLeptonFromWDecay().get());
642     mcList bTruList = getMCList(getLeptonicBJet().get());
643     for (int lepTypes = 0; lepTypes < 2; ++lepTypes) {
644     int lepInd = -1;
645     int lepPdgId = ELECTRON_ID;
646     if (lepTypes == 1)
647     lepPdgId = MU_ID;
648     while (lepTruList.size() > 0 && (lepInd =
649     findMCParticle(lepTruList, lepPdgId, -lepPdgId, lepInd + 1)) >= 0)
650     {
651     numFound = 1;
652     int bInd = -1;
653     int lepGPInd = lepTruList.at(lepInd).gpInd;
654     if (lepGPInd >= 0 && lepGPInd < numParticles) {
655     int lepMothInd = mcColl.at(lepGPInd)->motherIndex();
656     if (lepMothInd >= 0 && lepMothInd < numParticles &&
657     fabs(mcColl.at(lepMothInd)->pdgId()) == W_ID) {
658     numFound = 2;
659     while (bTruList.size() > 0 &&
660     (bInd = findMCParticle(bTruList, BQ_ID, -BQ_ID, bInd + 1)) >= 0)
661     {
662     numFound = 3;
663     int bGPInd = bTruList.at(bInd).gpInd;
664     if ((bGPInd >= 0 && bGPInd < numParticles) &&
665     (chkIfChild(bGPInd, mcColl.at(lepMothInd)->motherIndex()) ||
666     chkIfChild(lepMothInd, mcColl.at(bGPInd)->motherIndex())))
667     return(4);
668     }
669     }
670     }
671     if (numFound > mostFound)
672     mostFound = numFound;
673     }
674     }
675     return (mostFound);
676     }
677    
678    
679    
680     /*
681     int ToplikeCandidate::getNumMCPermutsForTop(const ParticlePointer &recoTop,
682     bool hadronic) const
683     {
684     mcList topTru = getMCList(recoTop.get());
685     if (topTru.size() <= 0)
686     return (0);
687     int highestNumFound = 0;
688     int topInd = -1;
689     while ((topInd = findMCParticle(topTru, TOP_ID, -TOP_ID, topInd + 1)) >= 0) {
690     int numFound = 0;
691     if (hadronic)
692     numFound = hTopDauMatching(topTru.at(topInd));
693     else numFound = lTopDauMatching(topTru.at(topInd));
694     if (numFound > highestNumFound)
695     highestNumFound = numFound;
696     }
697     return (highestNumFound);
698     }
699     */
700    
701    
702     int ToplikeCandidate::getNumMCMatchesForTop(const ParticlePointer &recoTop,
703     bool hadronic)
704     const
705     {
706     mcList topTru = getMCList(recoTop.get());
707     if (topTru.size() <= 0)
708     return (0);
709     int highestNumFound = 0;
710     int topInd = -1;
711     while ((topInd = findMCParticle(topTru, TOP_ID, -TOP_ID, topInd + 1)) >= 0) {
712     int numFound = 0;
713     if (hadronic)
714     numFound = strictHTopMatching(topTru.at(topInd));
715     else numFound = strictLTopMatching(topTru.at(topInd));
716     if (numFound > highestNumFound)
717     highestNumFound = numFound;
718     }
719     return (highestNumFound);
720     }
721    
722    
723     int ToplikeCandidate::getNumMCMatchesHTop() const
724     {
725     return (getNumMCMatchesForTop(getHadronicTop(), true));
726     }
727    
728     /*
729     int ToplikeCandidate::getNumDauMatchesHTop() const
730     {
731     return (getNumMCPermutsForTop(getHadronicTop(), true));
732     }
733    
734    
735     int ToplikeCandidate::getNumDauMatchesLTop() const
736     {
737     return (getNumMCPermutsForTop(getLeptonicTop(), false));
738     }
739     */
740    
741    
742     int ToplikeCandidate::getNumCorrectIDHTop() const
743     {
744     mcList topTru = getMCList(getHadronicTop().get());
745     int numFound = 0;
746     bool topFnd = false, bFnd = false, jet1Fnd = false, jet2Fnd = false;
747     if (topTru.size() > 0 && findMCParticle(topTru, TOP_ID, -TOP_ID) >= 0) {
748     topFnd = true;
749     ++numFound;
750     }
751    
752     // Take either particle or anti-particle.
753     if (chkMCMatch(hadronicBJet.get(), -1, BQ_ID, -BQ_ID) >= 0) {
754     ++numFound;
755     bFnd = true;
756     }
757     int wInd = chkMCMatch(hadronicW.get(), -1, W_ID, -W_ID);
758     if (wInd >= 0)
759     ++numFound;
760    
761     // Look for any quark near a jet. Even allow the same quark for both jets.
762     if (chkForQuark(jet1FromW, -1) >= 0) {
763     ++numFound;
764     jet1Fnd = true;
765     }
766     if (chkForQuark(jet2FromW, -1) >= 0) {
767     jet2Fnd = true;
768     ++numFound;
769     }
770     /*
771     if (numFound == 4 || numFound == 3) {
772     cout << "had num t b w j1 j2 " << numFound << " " << topFnd << " " << bFnd << " " << wInd;
773     cout << " " << jet1Fnd << " " << jet2Fnd << endl;
774     }
775     */
776     return (numFound);
777     }
778    
779    
780     int ToplikeCandidate::getNumMCMatchesLTop() const
781     {
782     return (getNumMCMatchesForTop(getLeptonicTop(), false));
783     }
784    
785    
786     int ToplikeCandidate::getNumCorrectIDLTop() const
787     {
788     mcList topTru = getMCList(getLeptonicTop().get());
789     int numFound = 0;
790     if (topTru.size() > 0 && findMCParticle(topTru, TOP_ID, -TOP_ID) >= 0)
791     ++numFound;
792     if (chkMCMatch(getLeptonicBJet().get(), -1, BQ_ID, -BQ_ID) >= 0)
793     ++numFound;
794     if (chkMCMatch(getLeptonicW().get(), -1, W_ID, -W_ID) >= 0)
795     ++numFound;
796     if (chkMCMatch(getLeptonFromWDecay().get(), -1, -ELECTRON_ID,
797     ELECTRON_ID) >= 0 ||
798     chkMCMatch(getLeptonFromWDecay().get(), -1, -MU_ID, MU_ID) >= 0)
799     ++numFound;
800     if (chkMCMatch(getNeutrinoFromWDecay().get(), -1, NU_E_ID, -NU_E_ID) >= 0 ||
801     chkMCMatch(getNeutrinoFromWDecay().get(), -1, NU_MU_ID, -NU_MU_ID) >= 0)
802     ++numFound;
803     return (numFound);
804     }
805    
806    
807     int ToplikeCandidate::getTausInJets() const
808     {
809     int numTaus = 0;
810     for (unsigned int index = 0; index < GoodJets().size(); ++index) {
811     mcList truList = getMCList(GoodJets().at(index).get());
812     int tauIndex = -1;
813     if (truList.size() > 0 &&
814     (tauIndex = findMCParticle(truList, TAU_ID, -TAU_ID)) >= 0)
815     ++numTaus;
816     }
817     // cout << " Found taus: " << numTaus;
818     return (numTaus);
819     }
820    
821    
822     int ToplikeCandidate::numTaus() const
823     {
824     const MCParticleCollection &mcColl = genParticles;
825     int numTaus = 0;
826     if (mcColl.size() > 0) {
827     for (unsigned int i = 0; i < mcColl.size(); ++i) {
828     if (fabs(mcColl.at(i)->pdgId()) == TAU_ID)
829     ++numTaus;
830     }
831     }
832     return (numTaus);
833     }
834    
835    
836     int ToplikeCandidate::missGenJetsForRecoJets() const
837     {
838     int matchJets = 0, numJets = GoodJets().size();
839     for (int index = 0; index < numJets; ++index) {
840     mcList mcObjs = getMCListForJets(GoodJets().at(index).get());
841     if (mcObjs.size() > 0 && mcObjs.at(0).deltaR < 0.3)
842     ++matchJets;
843     }
844     return (numJets - matchJets);
845     }
846    
847    
848    
849     TopPlusXCandidates::TopPlusXCandidates() :
850     dJetFromWpIndex(-1)
851     {
852     }
853    
854    
855     TopPlusXCandidates::TopPlusXCandidates(const Event& event) :
856     ToplikeCandidate(event),
857     dJetFromWpIndex(-1)
858     {
859     }
860    
861    
862     TopPlusXCandidates::~TopPlusXCandidates()
863     {
864     }
865    
866    
867     const JetPointer TopPlusXCandidates::nulljet;
868    
869    
870     bool TopPlusXCandidates::passesMETCut() const {
871     return met->pt() > 20.;
872     }
873    
874    
875     bool TopPlusXCandidates::passesHTCut() const {
876     return fullHT() > 700.;
877     }
878    
879    
880     bool TopPlusXCandidates::hasAtLeastOneGoodJet() const {
881     return goodJets.size() >= 1 && goodJets.front()->pt() > 180.0;
882     }
883    
884    
885     bool TopPlusXCandidates::hasAtLeastTwoGoodJets() const {
886     return goodJets.size() >= 2 && goodJets.at(1)->pt() > 90.0;
887     }
888    
889    
890     bool TopPlusXCandidates::hasOnlyOneGoodIsolatedMuon() const {
891     return goodIsolatedMuons.size() == 1;
892     }
893    
894    
895     bool TopPlusXCandidates::hasNoGoodIsolatedMuon() const {
896     return goodIsolatedMuons.size() <= 0;
897     }
898    
899    
900     bool TopPlusXCandidates::hasNoGoodIsolatedElectron() const {
901     if(Event::usePFIsolation)
902     return goodPFIsolatedElectrons.size() <= 0;
903     else
904     return goodIsolatedElectrons.size() <= 0;
905     }
906    
907    
908     bool TopPlusXCandidates::hasNoIsolatedElectron() const {
909     return looseElectrons.size() <= 0;
910     }
911    
912     bool TopPlusXCandidates::passesExtraElectronCut() const{
913     ElectronPointer electron;
914     if(Event::usePFIsolation) {
915     if (goodPFIsolatedElectrons.size() <= 0)
916     return (false);
917     electron = GoodPFIsolatedElectrons().front();
918     } else if (goodIsolatedElectrons.size() > 0) {
919     electron = GoodIsolatedElectrons().front();
920     } else return (false);
921     return (electron->pt() > 45.);
922     }
923    
924    
925     bool TopPlusXCandidates::passesSelectionStep(enum
926     ToplikeSelectionSteps::Step step) const {
927     switch (step) {
928     case ToplikeSelectionSteps::AtLeastFiveGoodJets:
929     return (hasAtLeastFiveGoodJets());
930     case ToplikeSelectionSteps::FortyfiveGeVElectron:
931     return (passesExtraElectronCut());
932     case ToplikeSelectionSteps::OneIsolatedMuon:
933     return (hasOnlyOneGoodIsolatedMuon());
934     case ToplikeSelectionSteps::LooseElectronVeto:
935     return (hasNoIsolatedElectron());
936     case ToplikeSelectionSteps::TightElectronVeto:
937     return (hasNoGoodIsolatedElectron());
938     case ToplikeSelectionSteps::TightMuonVeto:
939     return (hasNoGoodIsolatedMuon());
940     case ToplikeSelectionSteps::HT700GeV:
941     return (passesHTCut());
942     default:
943     return (TopPairEventCandidate::passesSelectionStep(selSteps[step]));
944     }
945     }
946    
947    
948     bool TopPlusXCandidates::passesSelectionStepUpTo(unsigned int step,
949     const ToplikeSelectionSteps::Step steps[]) const {
950     if (step <= 0)
951     return (passesSelectionStep(steps[0]));
952     else {
953     return (passesSelectionStep(steps[step]) &&
954     passesSelectionStepUpTo(step - 1, steps));
955     }
956     }
957    
958    
959     bool TopPlusXCandidates::passesFullEPlusJetSelection() const {
960     return (passesSelectionStepUpTo(toplikeElectronSelSize - 1, toplikeElectronSelection));
961     }
962    
963    
964     bool TopPlusXCandidates::passesFullMuPlusJetSelection() const {
965     return (passesSelectionStepUpTo(toplikeMuonSelSize - 1, toplikeMuonSelection));
966     }
967    
968    
969     bool TopPlusXCandidates::passesFullEPlusJetSelectionKK() const {
970     return (passesSelectionStepUpTo(toplikeElectronSelSize - 1, toplikeElectronSelection));
971     }
972    
973    
974     double TopPlusXCandidates::jetsHTNoTop(unsigned int numJets) const {
975     double ht = 0.0;
976     JetCollection topJets;
977    
978     for (unsigned int index = 0; index < goodJets.size(); ++index) {
979     if (goodJets.at(index) != leptonicBJet) {
980     if (topJets.size() <= 0 || topJets.back()->pt() > goodJets.at(index)->pt())
981     topJets.push_back(goodJets.at(index));
982     else {
983     unsigned int srchind = 0;
984     while (topJets.at(srchind)->pt() > goodJets.at(index)->pt())
985     ++srchind;
986     topJets.insert(topJets.begin() + srchind, goodJets.at(index));
987     }
988     }
989     }
990     for (unsigned int index = 0; index < topJets.size() && index < numJets; ++index) {
991     ht += topJets.at(index)->pt();
992     }
993     return ht;
994     }
995    
996    
997     double TopPlusXCandidates::mass2jets(unsigned int jet1, unsigned int jet2) const
998     {
999     int hijets[3] = {-1, -1, -1};
1000     unsigned int hijindex = 0;
1001    
1002     for (unsigned int index = 0; index < goodJets.size(); ++index) {
1003     if (goodJets.at(index) != leptonicBJet && hijets[hijindex] == -1) {
1004     hijets[hijindex] = index;
1005     if (hijindex < 2)
1006     hijindex++;
1007     }
1008     }
1009     if (hijets[jet1] > -1 && hijets[jet2] > -1) {
1010     ParticlePointer candidate =
1011     ParticlePointer(new Particle(*goodJets.at(hijets[jet1]) +
1012     *goodJets.at(hijets[jet2])));
1013     return (candidate->mass());
1014     }
1015     return (0.0);
1016     }
1017    
1018    
1019     double TopPlusXCandidates::mass4jets() const
1020     {
1021     unsigned int cnt = 0;
1022    
1023     ParticlePointer candidate;
1024     for (unsigned int index = 0; index < goodJets.size(); ++index) {
1025     if (goodJets.at(index) != leptonicBJet) {
1026     if (cnt == 0) {
1027     candidate = ParticlePointer(new Particle(*goodJets.at(index)));
1028     } else if (cnt < 3) {
1029     candidate = ParticlePointer(new Particle(*goodJets.at(index) + *candidate));
1030     }
1031     ++cnt;
1032     }
1033     }
1034     if (cnt > 0)
1035     return (candidate->mass());
1036     return (0.0);
1037     }
1038    
1039    
1040     JetPointer TopPlusXCandidates::leadnontopjet() const {
1041     double hiJetPt = 0;
1042     unsigned int hiInd = 0;
1043    
1044     for (unsigned int index = 0; index < goodJets.size(); ++index) {
1045     if (goodJets.at(index) != leptonicBJet && goodJets.at(index)->pt() > hiJetPt) {
1046     hiJetPt = goodJets.at(index)->pt();
1047     hiInd = index;
1048     }
1049     }
1050     return (goodJets.at(hiInd));
1051     }
1052    
1053    
1054     // ST is the sum of jet and lepton pT minus the pT of the leading lepton and jet.
1055     double TopPlusXCandidates::ST() const {
1056     double st = fullHT() - met->pt();
1057    
1058     double hiLepPt = 0;
1059     for (unsigned int index = 0; index < goodPFIsolatedElectrons.size(); ++index) {
1060     if (goodPFIsolatedElectrons.at(index)->pt() > hiLepPt)
1061     hiLepPt = goodPFIsolatedElectrons.at(index)->pt();
1062     }
1063     for (unsigned int index = 0; index < goodIsolatedMuons.size(); ++index) {
1064     if (goodIsolatedMuons.at(index)->pt() > hiLepPt)
1065     hiLepPt = goodIsolatedMuons.at(index)->pt();
1066     }
1067     double hiJetPt = 0;
1068     for (unsigned int index = 0; index < goodJets.size(); ++index) {
1069     if (goodJets.at(index)->pt() > hiJetPt)
1070     hiJetPt = goodJets.at(index)->pt();
1071     }
1072     st -= hiLepPt + hiJetPt;
1073     return st;
1074     }
1075    
1076     double TopPlusXCandidates::getTopMetAngle() const
1077     {
1078     return (getLeptonicTop()->getFourVector().Theta() - met->getFourVector().Theta());
1079     }
1080    
1081    
1082     void TopPlusXCandidates::recoBestLeptoTop(LeptonPointer lepton,
1083     const JetPointer hadTopJet, const JetPointer wJet1, const JetPointer wJet2)
1084     {
1085     if (goodBJets.size() < 1)
1086     throw ReconstructionException("Not enough jets available to reconstruct top using Chi2 method.");
1087     double bestLepChi2(9999999.);
1088     // JetCollection *topJets = &goodBJets;
1089     JetCollection *topJets = &goodJets;
1090     if (topJets->size() <= 1 && hadTopJet != 0)
1091     topJets = &goodJets;
1092     leptonFromW = lepton;
1093     selectedNeutrino = 0;
1094     currentSelectedNeutrino = 0;
1095     reconstructNeutrinos();
1096    
1097     for (unsigned short jet1Index = 0; jet1Index < topJets->size(); ++jet1Index) {
1098     if (hadTopJet == topJets->at(jet1Index) ||
1099     dJetFromWp == topJets->at(jet1Index) ||
1100     wJet1 == topJets->at(jet1Index) ||
1101     wJet2 == topJets->at(jet1Index))
1102     continue;
1103     leptonicBJet = topJets->at(jet1Index);
1104     leptonicW1 = ParticlePointer(new Particle(*neutrino1 + *leptonFromW));
1105     leptonicW2 = ParticlePointer(new Particle(*neutrino2 + *leptonFromW));
1106     leptonicTop1 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW1));
1107     leptonicTop2 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW2));
1108     selectNeutrinoSolution();
1109     double ltopchi2 = TopPairEventCandidate::getLeptonicChi2(currentSelectedNeutrino);
1110     if (ltopchi2 < bestLepChi2) {
1111     selectedNeutrino = currentSelectedNeutrino;
1112     leptonicBIndex = jet1Index;
1113     bestLepChi2 = ltopchi2;
1114     }
1115     }
1116     leptonicBJet = topJets->at(leptonicBIndex);
1117     leptonicW1 = ParticlePointer(new Particle(*neutrino1 + *leptonFromW));
1118     leptonicW2 = ParticlePointer(new Particle(*neutrino2 + *leptonFromW));
1119     leptonicTop1 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW1));
1120     leptonicTop2 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW2));
1121     doneReconstructiontop = true;
1122     doneReconstruction = true;
1123     }
1124    
1125    
1126     void TopPlusXCandidates::findDJet(const JetPointer topJet)
1127     {
1128     dJetFromWpIndex = -1;
1129     // dJetFromWp = 0; // unnecessary
1130     unsigned short cnt = 0;
1131     do {
1132     for (unsigned short jet1Index = 0; jet1Index < goodJets.size(); ++jet1Index)
1133     {
1134     if ((goodJets.at(jet1Index)->isBJet() == false || cnt >= 1) &&
1135     goodJets.at(jet1Index) != topJet &&
1136     jet1FromW != goodJets.at(jet1Index) &&
1137     jet2FromW != goodJets.at(jet1Index))
1138     {
1139     dJetFromWpIndex = jet1Index;
1140     dJetFromWp = goodJets.at(jet1Index);
1141     break;
1142     }
1143     }
1144     ++cnt;
1145     } while (dJetFromWpIndex == -1 && cnt < 2);
1146     if (dJetFromWpIndex == -1)
1147     throw ReconstructionException("No d jet found -- all jets used for tops.");
1148     }
1149    
1150    
1151     void TopPlusXCandidates::recoHadronicTop(LeptonPointer lepton,
1152     const JetPointer lepTopJet)
1153     {
1154     if (goodJets.size() < 4)
1155     throw ReconstructionException("Not enough jets available to reconstruct top using Chi2 method.");
1156     double bestHadChi2(9999999.);
1157     // JetCollection *topJets = &goodBJets;
1158     JetCollection *topJets = &goodJets;
1159     if (topJets->size() <= 1 && lepTopJet != 0)
1160     topJets = &goodJets;
1161    
1162     for (unsigned short hadBindex = 0; hadBindex < topJets->size(); ++hadBindex) {
1163     if (lepTopJet == topJets->at(hadBindex) ||
1164     dJetFromWp == topJets->at(hadBindex))
1165     continue;
1166     for (unsigned short jet1Index = 0; jet1Index < goodJets.size(); ++jet1Index) {
1167     if (lepTopJet == goodJets.at(jet1Index) ||
1168     goodJets.at(jet1Index) == topJets->at(hadBindex) ||
1169     jet1Index == dJetFromWpIndex)
1170     continue;
1171     for (unsigned short jet2Index = 0; jet2Index < goodJets.size();
1172     ++jet2Index) {
1173     if (jet2Index == jet1Index ||
1174     goodJets.at(jet2Index) == topJets->at(hadBindex) ||
1175     jet2Index == dJetFromWpIndex ||
1176     lepTopJet == goodJets.at(jet2Index))
1177     continue;
1178     hadronicBJet = topJets->at(hadBindex);
1179     jet1FromW = goodJets.at(jet1Index);
1180     jet2FromW = goodJets.at(jet2Index);
1181     hadronicW = ParticlePointer(new Particle(*jet1FromW + *jet2FromW));
1182     hadronicTop = ParticlePointer(new Particle(*hadronicBJet + *hadronicW));
1183     double htopchi2 = getHadronicChi2();
1184     if (htopchi2 < bestHadChi2) {
1185     hadronicBIndex = hadBindex;
1186     jet1FromWIndex = jet1Index;
1187     jet2FromWIndex = jet2Index;
1188     bestHadChi2 = htopchi2;
1189     }
1190     }
1191     }
1192     }
1193     hadronicBJet = topJets->at(hadronicBIndex);
1194     jet1FromW = goodJets.at(jet1FromWIndex);
1195     jet2FromW = goodJets.at(jet2FromWIndex);
1196     hadronicW = ParticlePointer(new Particle(*jet1FromW + *jet2FromW));
1197     hadronicTop = ParticlePointer(new Particle(*hadronicBJet + *hadronicW));
1198     doneReconstructiontop = true;
1199     doneReconstruction = true;
1200     }
1201    
1202    
1203     TwoNonResTops::TwoNonResTops() :
1204     wpGenParticlesSet(false)
1205     {
1206     }
1207    
1208    
1209     TwoNonResTops::TwoNonResTops(const Event& event) :
1210     TopPlusXCandidates(event),
1211     wpGenParticlesSet(false)
1212     {
1213     }
1214    
1215    
1216     TwoNonResTops::~TwoNonResTops()
1217     {
1218     }
1219    
1220     double TwoNonResTops::getGlobalChi2() const { // No global chi2 for this class
1221     return (0.0);
1222     }
1223    
1224    
1225     void TwoNonResTops::recoNonResTops(LeptonPointer lepton) {
1226     if (goodJets.size() < 4)
1227     throw ReconstructionException("Not enough jets available to reconstruct two tops");
1228     // const JetCollection &topJets = goodBJets;
1229     const JetCollection &topJets = goodJets;
1230     leptonFromW = lepton;
1231     selectedNeutrino = 0;
1232     currentSelectedNeutrino = 0;
1233     reconstructNeutrinos();
1234     double chosen_TopMassDifference(9999999.);
1235     double chosen_Chi2Total(9999999.);
1236    
1237     for (unsigned short hadBindex = 0; hadBindex < topJets.size(); ++hadBindex) {
1238     for (unsigned short lepBindex = 0; lepBindex < topJets.size(); ++lepBindex) {
1239     if (lepBindex == hadBindex)
1240     continue;
1241     for (unsigned short jet1Index = 0; jet1Index < goodJets.size(); ++jet1Index) {
1242     if (goodJets.at(jet1Index) == topJets.at(lepBindex) ||
1243     goodJets.at(jet1Index) == topJets.at(hadBindex))
1244     continue;
1245     for (unsigned short jet2Index = 0; jet2Index < goodJets.size(); ++jet2Index) {
1246     if (jet2Index == jet1Index ||
1247     goodJets.at(jet2Index) == topJets.at(lepBindex) ||
1248     goodJets.at(jet2Index) == topJets.at(hadBindex))
1249     continue;
1250     hadronicBJet = topJets.at(hadBindex);
1251     leptonicBJet = topJets.at(lepBindex);
1252     jet1FromW = goodJets.at(jet1Index);
1253     jet2FromW = goodJets.at(jet2Index);
1254    
1255     leptonicW1 = ParticlePointer(new Particle(*neutrino1 + *leptonFromW));
1256     leptonicW2 = ParticlePointer(new Particle(*neutrino2 + *leptonFromW));
1257     hadronicW = ParticlePointer(new Particle(*jet1FromW + *jet2FromW));
1258     leptonicTop1 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW1));
1259     leptonicTop2 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW2));
1260     hadronicTop = ParticlePointer(new Particle(*hadronicBJet + *hadronicW));
1261     fillHypotheses();
1262     selectNeutrinoSolution();
1263     double TopMassDifference = calculateTopMassDifference(currentSelectedNeutrino);
1264     double chi2 = getTotalChi2(currentSelectedNeutrino);
1265     switch (usedTTbarReconstruction) {
1266     case TTbarReconstructionCriterion::TopMassDifference:
1267     if (TopMassDifference < chosen_TopMassDifference) {
1268     hadronicBIndex = hadBindex;
1269     leptonicBIndex = lepBindex;
1270     jet1FromWIndex = jet1Index;
1271     jet2FromWIndex = jet2Index;
1272     chosen_TopMassDifference = TopMassDifference;
1273     selectedNeutrino = currentSelectedNeutrino;
1274     }
1275     break;
1276     case TTbarReconstructionCriterion::chi2:
1277     if (chi2 < chosen_Chi2Total) {
1278     hadronicBIndex = hadBindex;
1279     leptonicBIndex = lepBindex;
1280     jet1FromWIndex = jet1Index;
1281     jet2FromWIndex = jet2Index;
1282     chosen_Chi2Total = chi2;
1283     selectedNeutrino = currentSelectedNeutrino;
1284     }
1285     break;
1286     }
1287     }
1288     }
1289     }
1290     }
1291     std::sort(solutions.begin(), solutions.end(), compareSolutions);
1292     hadronicBJet = topJets.at(hadronicBIndex);
1293     leptonicBJet = topJets.at(leptonicBIndex);
1294     jet1FromW = goodJets.at(jet1FromWIndex);
1295     jet2FromW = goodJets.at(jet2FromWIndex);
1296     leptonicW1 = ParticlePointer(new Particle(*neutrino1 + *leptonFromW));
1297     leptonicW2 = ParticlePointer(new Particle(*neutrino2 + *leptonFromW));
1298     hadronicW = ParticlePointer(new Particle(*jet1FromW + *jet2FromW));
1299     leptonicTop1 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW1));
1300     leptonicTop2 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW2));
1301     hadronicTop = ParticlePointer(new Particle(*hadronicBJet + *hadronicW));
1302     if (selectedNeutrino == 1) // Resonance may be real of fake -- can exclude ttbar?
1303     ttbarResonance = ParticlePointer(new Particle(*leptonicTop1 + *hadronicTop));
1304     else
1305     ttbarResonance = ParticlePointer(new Particle(*leptonicTop2 + *hadronicTop));
1306     doneReconstruction = true;
1307     }
1308    
1309    
1310     double TwoNonResTops::getTotalChi2DJ(unsigned short neutrinoSolution) const
1311     {
1312     double chi2 = TopPairEventCandidate::getTotalChi2(neutrinoSolution);
1313     double pt = dJetFromWp->pt();
1314     if (pt < jet1FromW->pt())
1315     chi2 += (0.5/sqrt(2.0));
1316     if (pt < jet2FromW->pt())
1317     chi2 += (0.5/sqrt(2.0));
1318     double deltaphi = dJetFromWp->deltaPhi(hadronicBJet);
1319     double deltaphi2 = dJetFromWp->deltaPhi(leptonicBJet);
1320     if (deltaphi2 > deltaphi)
1321     deltaphi = deltaphi2;
1322     double dphichi = deltaphi - pi;
1323     chi2 += (1.0/sqrt(2.0)) * (dphichi * dphichi)/(pi * pi/2.0);
1324     if (dJetFromWp->isBJet())
1325     chi2 += (0.5/sqrt(2.0));
1326     return (chi2);
1327     }
1328    
1329    
1330     void TwoNonResTops::recoNonResTopsWdjet(LeptonPointer lepton) {
1331     if (goodJets.size() < 5)
1332     throw ReconstructionException("Not enough jets available to reconstruct two tops");
1333     const JetCollection &topJets = goodJets;
1334     leptonFromW = lepton;
1335     selectedNeutrino = 0;
1336     currentSelectedNeutrino = 0;
1337     reconstructNeutrinos();
1338     double chosen_TopMassDifference(9999999.);
1339     double chosen_Chi2Total(9999999.);
1340    
1341     for (unsigned short hadBindex = 0; hadBindex < topJets.size(); ++hadBindex) {
1342     for (unsigned short lepBindex = 0; lepBindex < topJets.size(); ++lepBindex) {
1343     if (lepBindex == hadBindex)
1344     continue;
1345     for (unsigned short jet1Index = 0; jet1Index < goodJets.size(); ++jet1Index) {
1346     if (goodJets.at(jet1Index) == topJets.at(lepBindex) ||
1347     goodJets.at(jet1Index) == topJets.at(hadBindex))
1348     continue;
1349     for (unsigned short jet2Index = 0; jet2Index < goodJets.size(); ++jet2Index) {
1350     if (jet2Index == jet1Index ||
1351     goodJets.at(jet2Index) == topJets.at(lepBindex) ||
1352     goodJets.at(jet2Index) == topJets.at(hadBindex))
1353     continue;
1354     for (unsigned short djetIndex = 0; djetIndex < goodJets.size(); ++djetIndex) {
1355     if (djetIndex == jet1Index || jet2Index == djetIndex ||
1356     goodJets.at(djetIndex) == topJets.at(lepBindex) ||
1357     goodJets.at(djetIndex) == topJets.at(hadBindex))
1358     continue;
1359     hadronicBJet = topJets.at(hadBindex);
1360     leptonicBJet = topJets.at(lepBindex);
1361     jet1FromW = goodJets.at(jet1Index);
1362     jet2FromW = goodJets.at(jet2Index);
1363     dJetFromWp = goodJets.at(djetIndex);
1364    
1365     leptonicW1 = ParticlePointer(new Particle(*neutrino1 + *leptonFromW));
1366     leptonicW2 = ParticlePointer(new Particle(*neutrino2 + *leptonFromW));
1367     hadronicW = ParticlePointer(new Particle(*jet1FromW + *jet2FromW));
1368     leptonicTop1 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW1));
1369     leptonicTop2 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW2));
1370     hadronicTop = ParticlePointer(new Particle(*hadronicBJet + *hadronicW));
1371     fillHypotheses();
1372     selectNeutrinoSolution();
1373     double TopMassDifference = calculateTopMassDifference(currentSelectedNeutrino);
1374     double chi2 = getTotalChi2DJ(currentSelectedNeutrino);
1375     switch (usedTTbarReconstruction) {
1376     case TTbarReconstructionCriterion::TopMassDifference:
1377     if (TopMassDifference < chosen_TopMassDifference) {
1378     hadronicBIndex = hadBindex;
1379     leptonicBIndex = lepBindex;
1380     jet1FromWIndex = jet1Index;
1381     jet2FromWIndex = jet2Index;
1382     chosen_TopMassDifference = TopMassDifference;
1383     selectedNeutrino = currentSelectedNeutrino;
1384     }
1385     break;
1386     case TTbarReconstructionCriterion::chi2:
1387     if (chi2 < chosen_Chi2Total) {
1388     hadronicBIndex = hadBindex;
1389     leptonicBIndex = lepBindex;
1390     jet1FromWIndex = jet1Index;
1391     jet2FromWIndex = jet2Index;
1392     dJetFromWpIndex = djetIndex;
1393     chosen_Chi2Total = chi2;
1394     selectedNeutrino = currentSelectedNeutrino;
1395     }
1396     break;
1397     }
1398     }
1399     }
1400     }
1401     }
1402     }
1403     std::sort(solutions.begin(), solutions.end(), compareSolutions);
1404     hadronicBJet = topJets.at(hadronicBIndex);
1405     leptonicBJet = topJets.at(leptonicBIndex);
1406     jet1FromW = goodJets.at(jet1FromWIndex);
1407     jet2FromW = goodJets.at(jet2FromWIndex);
1408     dJetFromWp = goodJets.at(dJetFromWpIndex);
1409     leptonicW1 = ParticlePointer(new Particle(*neutrino1 + *leptonFromW));
1410     leptonicW2 = ParticlePointer(new Particle(*neutrino2 + *leptonFromW));
1411     hadronicW = ParticlePointer(new Particle(*jet1FromW + *jet2FromW));
1412     leptonicTop1 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW1));
1413     leptonicTop2 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW2));
1414     hadronicTop = ParticlePointer(new Particle(*hadronicBJet + *hadronicW));
1415     if (selectedNeutrino == 1) // Resonance may be real of fake -- can exclude ttbar?
1416     ttbarResonance = ParticlePointer(new Particle(*leptonicTop1 + *hadronicTop));
1417     else
1418     ttbarResonance = ParticlePointer(new Particle(*leptonicTop2 + *hadronicTop));
1419     doneReconstruction = true;
1420     }
1421    
1422    
1423     const JetPointer TwoNonResTops::getLtNonTopJet() const {
1424     throwExpeptionIfNotReconstructed("TwoNonResTops::getLtNonTopJet");
1425     if (goodJets.size() < 5)
1426     throw ReconstructionException("Not enough jets available to reconstruct Wprime");
1427     int ltJetInd = -1, bJetInd = -1;
1428     unsigned int index = 0;
1429     while (index < goodJets.size() && ltJetInd < 0) {
1430     JetPointer currJet = goodJets.at(index);
1431     if (currJet->isBJet() == false && currJet != leptonicBJet &&
1432     currJet != hadronicBJet && currJet != jet1FromW && currJet != jet2FromW) {
1433     ltJetInd = index;
1434     } else if (bJetInd < 0 && currJet != leptonicBJet && currJet != hadronicBJet &&
1435     currJet != jet1FromW && currJet != jet2FromW) {
1436     bJetInd = index;
1437     }
1438     ++index;
1439     }
1440     if (ltJetInd < 0 && bJetInd >= 0) // Take a B jet if no others found.
1441     ltJetInd = bJetInd;
1442     if (ltJetInd >= 0)
1443     return (goodJets.at(ltJetInd));
1444     throw ReconstructionException(TString("getLtNonTopJet:Fifth jet missing"));
1445     }
1446    
1447    
1448     const ParticlePointer TwoNonResTops::getWprime() const {
1449     throwExpeptionIfNotReconstructed("TwoNonResTops::getWprime");
1450     ParticlePointer ltJet = getLtNonTopJet();
1451     ParticlePointer hiPtTop = getLeptonicTop();
1452    
1453     // Generator study indicates lower Pt top is usually spectator
1454     if (hiPtTop->pt() < hadronicTop->pt())
1455     hiPtTop = hadronicTop;
1456     return (ParticlePointer(new Particle(*hiPtTop + *ltJet)));
1457     }
1458    
1459    
1460     const ParticlePointer TwoNonResTops::getNegWprime() const {
1461     throwExpeptionIfNotReconstructed("TwoNonResTops::getNegWprime");
1462     ParticlePointer ltJet = getLtNonTopJet();
1463     ParticlePointer negTop = hadronicTop;
1464     if (leptonFromW->charge() < 0)
1465     negTop = getLeptonicTop();
1466    
1467     return (ParticlePointer(new Particle(*negTop + *ltJet)));
1468     // return (ParticlePointer(new Particle(*negTop + *dJetFromWp)));
1469     }
1470    
1471    
1472     // Deliberately reco the incorrect W' to see if ttbar still looks the same
1473     const ParticlePointer TwoNonResTops::getNegWprimeBad() const {
1474     throwExpeptionIfNotReconstructed("TwoNonResTops::getNegWprimeBad");
1475     ParticlePointer ltJet = getLtNonTopJet();
1476     ParticlePointer negTop = hadronicTop;
1477     if (leptonFromW->charge() > 0)
1478     negTop = getLeptonicTop(); // Use wrong top
1479    
1480     return (ParticlePointer(new Particle(*negTop + *ltJet)));
1481     // return (ParticlePointer(new Particle(*negTop + *dJetFromWp)));
1482     }
1483    
1484    
1485     const ParticlePointer TwoNonResTops::getBtBWprime() const {
1486     throwExpeptionIfNotReconstructed("TwoNonResTops::getBtBWprim");
1487     ParticlePointer ltJet = getLtNonTopJet();
1488     ParticlePointer hiAnglTop = getLeptonicTop();
1489     if (getHadTopDAngle() > hiAnglTop->angle(ltJet))
1490     hiAnglTop = hadronicTop;
1491     return (ParticlePointer(new Particle(*hiAnglTop + *ltJet)));
1492     }
1493    
1494    
1495     double TwoNonResTops::getHadTopDAngle() const
1496     {
1497     ParticlePointer jet = getLtNonTopJet();
1498     return (hadronicTop->angle(jet));
1499     }
1500    
1501    
1502     void TwoNonResTops::recoWpFrom1Top(LeptonPointer lepton) {
1503     if (goodJets.size() < 5)
1504     throw
1505     ReconstructionException("Not enough jets available to reconstruct W'");
1506     ParticlePointer negTop, posTop;
1507     if (lepton->charge() < 0) {
1508     // findDJet();
1509     recoBestLeptoTop(lepton);
1510     // findDJet(leptonicBJet);
1511     recoHadronicTop(lepton);
1512     // recoHadronicTop(lepton, leptonicBJet);
1513     negTop = getLeptonicTop();
1514     posTop = hadronicTop;
1515     } else {
1516     // findDJet();
1517     recoHadronicTop(lepton);
1518     // findDJet(hadronicBJet);
1519     recoBestLeptoTop(lepton);
1520     // recoBestLeptoTop(lepton, hadronicBJet, jet1FromW, jet2FromW);
1521     negTop = hadronicTop;
1522     posTop = getLeptonicTop();
1523     }
1524     JetPointer ltJet = getLtNonTopJet();
1525     wpFrom1Top = ParticlePointer(new Particle(*negTop + *ltJet));
1526     wpBadFrom1Top = ParticlePointer(new Particle(*posTop + *ltJet));
1527     topFromWp = negTop;
1528     dJetFromWp = ltJet;
1529     // wpFrom1Top = ParticlePointer(new Particle(*negTop + *dJetFromWp));
1530     // wpBadFrom1Top = ParticlePointer(new Particle(*posTop + *dJetFromWp));
1531     }
1532    
1533    
1534     void TwoNonResTops::getTrueWpChain()
1535     {
1536     const MCParticleCollection &mcColl = genParticles;
1537     const int nullInd = -99;
1538     int wpInd = nullInd, wpTopInd = nullInd, specTopInd = nullInd, dQInd = nullInd;
1539     int wpTopBInd = nullInd, specTopBInd = nullInd;
1540     if (mcColl.size() > 0) {
1541     for (unsigned int i = 0; i < mcColl.size(); ++i) {
1542     if (fabs(mcColl.at(i)->pdgId()) == WPRIME_ID) {
1543     if (wpInd == nullInd)
1544     wpInd = i;
1545     else if (mcColl.at(i)->motherIndex() != wpInd)
1546     cout << "Stray W': " << i << " mother " << mcColl.at(i)->motherIndex() << endl;
1547     }
1548     }
1549     for (unsigned int i = 0; i < mcColl.size(); ++i) {
1550     if (fabs(mcColl.at(i)->pdgId()) == TOP_ID) {
1551     if (mcColl.at(i)->motherIndex() == wpInd) {
1552     if (wpTopInd == nullInd)
1553     wpTopInd = i;
1554     else cout << "Additional W' top: " << i << endl;
1555     } else if (specTopInd == nullInd)
1556     specTopInd = i;
1557     else cout << "Additional spec top: " << i << endl;
1558     } else if (fabs(mcColl.at(i)->pdgId()) == DQ_ID) {
1559     if (mcColl.at(i)->motherIndex() == wpInd) {
1560     if (dQInd == nullInd)
1561     dQInd = i;
1562     else cout << "Additional d daughter: " << i << endl;
1563     }
1564     }
1565     }
1566     if (wpTopInd >= 0 && specTopInd >= 0) {
1567     int wptprevPdgId = 0;
1568     int sptprevPdgId = 0;
1569     for (unsigned int i = 0; i < mcColl.size(); ++i) {
1570     int pdgId = fabs(mcColl.at(i)->pdgId());
1571     if (mcColl.at(i)->motherIndex() == wpTopInd ||
1572     mcColl.at(i)->motherIndex() == wpInd) {
1573     if (mcColl.at(i)->isQuark() && pdgId % 2 == 1) {
1574     if (wpTopBInd == nullInd || pdgId > wptprevPdgId) {
1575     wpTopBInd = i;
1576     wptprevPdgId = pdgId;
1577     } else if (pdgId == BQ_ID)
1578     cout << "Additional W' top b quark: " << i << " mother " << mcColl.at(i)->motherIndex() << endl;
1579     }
1580     } else if (mcColl.at(i)->motherIndex() == specTopInd ||
1581     mcColl.at(i)->motherIndex() == mcColl.at(specTopInd)->motherIndex()) {
1582     if (mcColl.at(i)->isQuark() && pdgId % 2 == 1) {
1583     if (specTopBInd == nullInd || pdgId > sptprevPdgId) {
1584     specTopBInd = i;
1585     sptprevPdgId = pdgId;
1586     } else if (pdgId == BQ_ID)
1587     cout << "Additional spec top b quark: " << i << " mother " << mcColl.at(i)->motherIndex() << endl;
1588     }
1589     } else if (pdgId == BQ_ID && fabs(mcColl.at(mcColl.at(i)->motherIndex())->pdgId()) != BQ_ID) {
1590     // cout << "Additional b quark: " << i << " mother " << mcColl.at(i)->motherIndex();
1591     // cout << " ID " << mcColl.at(mcColl.at(i)->motherIndex())->pdgId() << endl;
1592     }
1593     }
1594     }
1595     if (wpInd >= 0 && wpTopInd >= 0 && specTopInd >= 0 && dQInd >= 0 &&
1596     wpTopBInd >= 0 && specTopBInd >= 0) {
1597     wpTru = mcColl.at(wpInd);
1598     wpTopTru = mcColl.at(wpTopInd);
1599     specTopTru = mcColl.at(specTopInd);
1600     wpTopBTru = mcColl.at(wpTopBInd);
1601     specTopBTru = mcColl.at(specTopBInd);
1602     wpDDauTru = mcColl.at(dQInd);
1603     wpGenParticlesSet = true;
1604     } else {
1605     // cout << "#11 : " << mcColl.at(11)->pdgId() << " mother " << mcColl.at(11)->motherIndex();
1606     // cout << " ID " << mcColl.at(mcColl.at(11)->motherIndex())->pdgId() << endl;
1607     cout << "wp wpTop specTop dq wpTopB specTopB " << wpInd << " " << wpTopInd << " " << specTopInd << " " << dQInd;
1608     cout << " " << wpTopBInd << " " << specTopBInd << endl;
1609     }
1610     }
1611     }
1612    
1613    
1614     MCParticlePointer TwoNonResTops::getWpTru() const
1615     {
1616     if (wpGenParticlesSet == false)
1617     throw ReconstructionException("Need to run TwoNonResTops::getTrueWpChain first");
1618     return (wpTru);
1619     }
1620    
1621     MCParticlePointer TwoNonResTops::getWpTopTru() const
1622     {
1623     if (wpGenParticlesSet == false)
1624     throw ReconstructionException("Need to run TwoNonResTops::getTrueWpChain first");
1625     return (wpTopTru);
1626     }
1627    
1628     MCParticlePointer TwoNonResTops::getSpecTop() const
1629     {
1630     if (wpGenParticlesSet == false)
1631     throw ReconstructionException("Need to run TwoNonResTops::getTrueWpChain first");
1632     return (specTopTru);
1633     }
1634    
1635     MCParticlePointer TwoNonResTops::getWpTopBTru() const
1636     {
1637     if (wpGenParticlesSet == false)
1638     throw ReconstructionException("Need to run TwoNonResTops::getTrueWpChain first");
1639     return (wpTopBTru);
1640     }
1641    
1642     MCParticlePointer TwoNonResTops::getSpecTopB() const
1643     {
1644     if (wpGenParticlesSet == false)
1645     throw ReconstructionException("Need to run TwoNonResTops::getTrueWpChain first");
1646     return (specTopBTru);
1647     }
1648    
1649     MCParticlePointer TwoNonResTops::getWpDDauTru() const
1650     {
1651     if (wpGenParticlesSet == false)
1652     throw ReconstructionException("Need to run TwoNonResTops::getTrueWpChain first");
1653     return (wpDDauTru);
1654     }
1655    
1656    
1657     TopNoMassConstraint::TopNoMassConstraint()
1658     {
1659     }
1660    
1661    
1662     TopNoMassConstraint::TopNoMassConstraint(const Event& event) :
1663     ToplikeCandidate(event)
1664     {
1665     }
1666    
1667    
1668     TopNoMassConstraint::~TopNoMassConstraint()
1669     {
1670     }
1671    
1672    
1673     double TopNoMassConstraint::getLoneHadChi2() const
1674     {
1675     double ptRatioDifference = TMath::Power(PtRatio() - matched_ptratio, 2);
1676     double ptRatioError = 2 * matched_ptratio_sigma * matched_ptratio_sigma;
1677     double ptRatioTerm = ptRatioDifference / ptRatioError;
1678    
1679     double WmassDifference = TMath::Power(hadronicW->mass() - pdgWmass, 2);
1680     double WmassError = 2 * wMassWidth * wMassWidth;
1681     double WmassTerm = WmassDifference / WmassError;
1682     return 1.0 / sqrt(2.0) * (WmassTerm + ptRatioTerm);
1683     }
1684    
1685     double TopNoMassConstraint::getLepBAngle() const
1686     {
1687     return (leptonicBJet->angle(leptonFromW));
1688     }
1689    
1690    
1691     LeptoTopNoMassConstraint::LeptoTopNoMassConstraint()
1692     {
1693     }
1694    
1695    
1696     LeptoTopNoMassConstraint::LeptoTopNoMassConstraint(const Event& event) :
1697     TopNoMassConstraint(event)
1698     {
1699     }
1700    
1701    
1702     LeptoTopNoMassConstraint::~LeptoTopNoMassConstraint()
1703     {
1704     }
1705    
1706    
1707     /*
1708     Tried for neutrino selection, but getGlobalChi2 still uses hadronic top quantities.
1709     double LeptoTopNoMassConstraint::getTotalChi2(unsigned short int neutrinoSolution) const {
1710     return TopPairEventCandidate::getLeptonicChi2(neutrinoSolution) + TopPairEventCandidate::getGlobalChi2(neutrinoSolution);
1711     }
1712     */
1713    
1714     void LeptoTopNoMassConstraint::inspectReconstructedEvent() const {
1715     cout << "run " << runNumber << ", event " << eventNumber << endl;
1716     cout << "leptonic b jet, goodJet index " << leptonicBIndex << endl;
1717     EventContentPrinter::printJet(leptonicBJet);
1718    
1719     cout << "lepton from W" << endl;
1720     // EventContentPrinter::printElectron(electronFromW);
1721    
1722     cout << "MET" << endl;
1723     EventContentPrinter::printParticle(met);
1724     cout << endl;
1725    
1726     cout << "reconstructed neutrino 1(selected: " << selectedNeutrino << ")" << endl;
1727     EventContentPrinter::printParticle(neutrino1);
1728     cout << endl;
1729    
1730     cout << "reconstructed neutrino 2(selected: " << selectedNeutrino << ")" << endl;
1731     EventContentPrinter::printParticle(neutrino2);
1732     cout << endl;
1733    
1734     cout << "leptonic W 1 (selected: " << selectedNeutrino << ")" << endl;
1735     EventContentPrinter::printParticle(leptonicW1);
1736     cout << endl;
1737    
1738     cout << "leptonic W 2 (selected: " << selectedNeutrino << ")" << endl;
1739     EventContentPrinter::printParticle(leptonicW2);
1740     cout << endl;
1741    
1742     cout << "leptonic top 1 (selected: " << selectedNeutrino << ")" << endl;
1743     EventContentPrinter::printParticle(leptonicTop1);
1744     cout << endl;
1745    
1746     cout << "leptonic top 2 (selected: " << selectedNeutrino << ")" << endl;
1747     EventContentPrinter::printParticle(leptonicTop2);
1748     cout << endl;
1749    
1750     cout << endl;
1751     }
1752    
1753    
1754     }