ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/dhidas/OSUAnalysis/Tools/src/ToplikeCandidate.cpp
Revision: 1.3
Committed: Sun Dec 18 16:37:29 2011 UTC (13 years, 4 months ago) by clseitz
Branch: MAIN
CVS Tags: HEAD
Changes since 1.2: +1 -1 lines
Log Message: 
general update

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 clseitz 1.3 return fullHT() > 700.0;
877 dhidas 1.1 }
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 clseitz 1.2 //CS
968     bool TopPlusXCandidates::passesFullEPlusJetSelectionMicro() const {
969     return (passesSelectionStepUpTo(toplikeElectronSelMicroSize - 1, toplikeElectronSelectionMicro));
970     }
971    
972    
973     bool TopPlusXCandidates::passesFullMuPlusJetSelectionMicro() const {
974     return (passesSelectionStepUpTo(toplikeMuonSelMicroSize - 1, toplikeMuonSelectionMicro));
975     }
976 dhidas 1.1
977    
978     bool TopPlusXCandidates::passesFullEPlusJetSelectionKK() const {
979     return (passesSelectionStepUpTo(toplikeElectronSelSize - 1, toplikeElectronSelection));
980     }
981    
982    
983     double TopPlusXCandidates::jetsHTNoTop(unsigned int numJets) const {
984     double ht = 0.0;
985     JetCollection topJets;
986    
987     for (unsigned int index = 0; index < goodJets.size(); ++index) {
988     if (goodJets.at(index) != leptonicBJet) {
989     if (topJets.size() <= 0 || topJets.back()->pt() > goodJets.at(index)->pt())
990     topJets.push_back(goodJets.at(index));
991     else {
992     unsigned int srchind = 0;
993     while (topJets.at(srchind)->pt() > goodJets.at(index)->pt())
994     ++srchind;
995     topJets.insert(topJets.begin() + srchind, goodJets.at(index));
996     }
997     }
998     }
999     for (unsigned int index = 0; index < topJets.size() && index < numJets; ++index) {
1000     ht += topJets.at(index)->pt();
1001     }
1002     return ht;
1003     }
1004    
1005    
1006     double TopPlusXCandidates::mass2jets(unsigned int jet1, unsigned int jet2) const
1007     {
1008     int hijets[3] = {-1, -1, -1};
1009     unsigned int hijindex = 0;
1010    
1011     for (unsigned int index = 0; index < goodJets.size(); ++index) {
1012     if (goodJets.at(index) != leptonicBJet && hijets[hijindex] == -1) {
1013     hijets[hijindex] = index;
1014     if (hijindex < 2)
1015     hijindex++;
1016     }
1017     }
1018     if (hijets[jet1] > -1 && hijets[jet2] > -1) {
1019     ParticlePointer candidate =
1020     ParticlePointer(new Particle(*goodJets.at(hijets[jet1]) +
1021     *goodJets.at(hijets[jet2])));
1022     return (candidate->mass());
1023     }
1024     return (0.0);
1025     }
1026    
1027    
1028     double TopPlusXCandidates::mass4jets() const
1029     {
1030     unsigned int cnt = 0;
1031    
1032     ParticlePointer candidate;
1033     for (unsigned int index = 0; index < goodJets.size(); ++index) {
1034     if (goodJets.at(index) != leptonicBJet) {
1035     if (cnt == 0) {
1036     candidate = ParticlePointer(new Particle(*goodJets.at(index)));
1037     } else if (cnt < 3) {
1038     candidate = ParticlePointer(new Particle(*goodJets.at(index) + *candidate));
1039     }
1040     ++cnt;
1041     }
1042     }
1043     if (cnt > 0)
1044     return (candidate->mass());
1045     return (0.0);
1046     }
1047    
1048    
1049     JetPointer TopPlusXCandidates::leadnontopjet() const {
1050     double hiJetPt = 0;
1051     unsigned int hiInd = 0;
1052    
1053     for (unsigned int index = 0; index < goodJets.size(); ++index) {
1054     if (goodJets.at(index) != leptonicBJet && goodJets.at(index)->pt() > hiJetPt) {
1055     hiJetPt = goodJets.at(index)->pt();
1056     hiInd = index;
1057     }
1058     }
1059     return (goodJets.at(hiInd));
1060     }
1061    
1062    
1063     // ST is the sum of jet and lepton pT minus the pT of the leading lepton and jet.
1064     double TopPlusXCandidates::ST() const {
1065     double st = fullHT() - met->pt();
1066    
1067     double hiLepPt = 0;
1068     for (unsigned int index = 0; index < goodPFIsolatedElectrons.size(); ++index) {
1069     if (goodPFIsolatedElectrons.at(index)->pt() > hiLepPt)
1070     hiLepPt = goodPFIsolatedElectrons.at(index)->pt();
1071     }
1072     for (unsigned int index = 0; index < goodIsolatedMuons.size(); ++index) {
1073     if (goodIsolatedMuons.at(index)->pt() > hiLepPt)
1074     hiLepPt = goodIsolatedMuons.at(index)->pt();
1075     }
1076     double hiJetPt = 0;
1077     for (unsigned int index = 0; index < goodJets.size(); ++index) {
1078     if (goodJets.at(index)->pt() > hiJetPt)
1079     hiJetPt = goodJets.at(index)->pt();
1080     }
1081     st -= hiLepPt + hiJetPt;
1082     return st;
1083     }
1084    
1085     double TopPlusXCandidates::getTopMetAngle() const
1086     {
1087     return (getLeptonicTop()->getFourVector().Theta() - met->getFourVector().Theta());
1088     }
1089    
1090    
1091     void TopPlusXCandidates::recoBestLeptoTop(LeptonPointer lepton,
1092     const JetPointer hadTopJet, const JetPointer wJet1, const JetPointer wJet2)
1093     {
1094     if (goodBJets.size() < 1)
1095     throw ReconstructionException("Not enough jets available to reconstruct top using Chi2 method.");
1096     double bestLepChi2(9999999.);
1097     // JetCollection *topJets = &goodBJets;
1098     JetCollection *topJets = &goodJets;
1099     if (topJets->size() <= 1 && hadTopJet != 0)
1100     topJets = &goodJets;
1101     leptonFromW = lepton;
1102     selectedNeutrino = 0;
1103     currentSelectedNeutrino = 0;
1104     reconstructNeutrinos();
1105    
1106     for (unsigned short jet1Index = 0; jet1Index < topJets->size(); ++jet1Index) {
1107     if (hadTopJet == topJets->at(jet1Index) ||
1108     dJetFromWp == topJets->at(jet1Index) ||
1109     wJet1 == topJets->at(jet1Index) ||
1110     wJet2 == topJets->at(jet1Index))
1111     continue;
1112     leptonicBJet = topJets->at(jet1Index);
1113     leptonicW1 = ParticlePointer(new Particle(*neutrino1 + *leptonFromW));
1114     leptonicW2 = ParticlePointer(new Particle(*neutrino2 + *leptonFromW));
1115     leptonicTop1 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW1));
1116     leptonicTop2 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW2));
1117     selectNeutrinoSolution();
1118     double ltopchi2 = TopPairEventCandidate::getLeptonicChi2(currentSelectedNeutrino);
1119     if (ltopchi2 < bestLepChi2) {
1120     selectedNeutrino = currentSelectedNeutrino;
1121     leptonicBIndex = jet1Index;
1122     bestLepChi2 = ltopchi2;
1123     }
1124     }
1125     leptonicBJet = topJets->at(leptonicBIndex);
1126     leptonicW1 = ParticlePointer(new Particle(*neutrino1 + *leptonFromW));
1127     leptonicW2 = ParticlePointer(new Particle(*neutrino2 + *leptonFromW));
1128     leptonicTop1 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW1));
1129     leptonicTop2 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW2));
1130     doneReconstructiontop = true;
1131     doneReconstruction = true;
1132     }
1133    
1134    
1135     void TopPlusXCandidates::findDJet(const JetPointer topJet)
1136     {
1137     dJetFromWpIndex = -1;
1138     // dJetFromWp = 0; // unnecessary
1139     unsigned short cnt = 0;
1140     do {
1141     for (unsigned short jet1Index = 0; jet1Index < goodJets.size(); ++jet1Index)
1142     {
1143     if ((goodJets.at(jet1Index)->isBJet() == false || cnt >= 1) &&
1144     goodJets.at(jet1Index) != topJet &&
1145     jet1FromW != goodJets.at(jet1Index) &&
1146     jet2FromW != goodJets.at(jet1Index))
1147     {
1148     dJetFromWpIndex = jet1Index;
1149     dJetFromWp = goodJets.at(jet1Index);
1150     break;
1151     }
1152     }
1153     ++cnt;
1154     } while (dJetFromWpIndex == -1 && cnt < 2);
1155     if (dJetFromWpIndex == -1)
1156     throw ReconstructionException("No d jet found -- all jets used for tops.");
1157     }
1158    
1159    
1160     void TopPlusXCandidates::recoHadronicTop(LeptonPointer lepton,
1161     const JetPointer lepTopJet)
1162     {
1163     if (goodJets.size() < 4)
1164     throw ReconstructionException("Not enough jets available to reconstruct top using Chi2 method.");
1165     double bestHadChi2(9999999.);
1166     // JetCollection *topJets = &goodBJets;
1167     JetCollection *topJets = &goodJets;
1168     if (topJets->size() <= 1 && lepTopJet != 0)
1169     topJets = &goodJets;
1170    
1171     for (unsigned short hadBindex = 0; hadBindex < topJets->size(); ++hadBindex) {
1172     if (lepTopJet == topJets->at(hadBindex) ||
1173     dJetFromWp == topJets->at(hadBindex))
1174     continue;
1175     for (unsigned short jet1Index = 0; jet1Index < goodJets.size(); ++jet1Index) {
1176     if (lepTopJet == goodJets.at(jet1Index) ||
1177     goodJets.at(jet1Index) == topJets->at(hadBindex) ||
1178     jet1Index == dJetFromWpIndex)
1179     continue;
1180     for (unsigned short jet2Index = 0; jet2Index < goodJets.size();
1181     ++jet2Index) {
1182     if (jet2Index == jet1Index ||
1183     goodJets.at(jet2Index) == topJets->at(hadBindex) ||
1184     jet2Index == dJetFromWpIndex ||
1185     lepTopJet == goodJets.at(jet2Index))
1186     continue;
1187     hadronicBJet = topJets->at(hadBindex);
1188     jet1FromW = goodJets.at(jet1Index);
1189     jet2FromW = goodJets.at(jet2Index);
1190     hadronicW = ParticlePointer(new Particle(*jet1FromW + *jet2FromW));
1191     hadronicTop = ParticlePointer(new Particle(*hadronicBJet + *hadronicW));
1192     double htopchi2 = getHadronicChi2();
1193     if (htopchi2 < bestHadChi2) {
1194     hadronicBIndex = hadBindex;
1195     jet1FromWIndex = jet1Index;
1196     jet2FromWIndex = jet2Index;
1197     bestHadChi2 = htopchi2;
1198     }
1199     }
1200     }
1201     }
1202     hadronicBJet = topJets->at(hadronicBIndex);
1203     jet1FromW = goodJets.at(jet1FromWIndex);
1204     jet2FromW = goodJets.at(jet2FromWIndex);
1205     hadronicW = ParticlePointer(new Particle(*jet1FromW + *jet2FromW));
1206     hadronicTop = ParticlePointer(new Particle(*hadronicBJet + *hadronicW));
1207     doneReconstructiontop = true;
1208     doneReconstruction = true;
1209     }
1210    
1211    
1212     TwoNonResTops::TwoNonResTops() :
1213     wpGenParticlesSet(false)
1214     {
1215     }
1216    
1217    
1218     TwoNonResTops::TwoNonResTops(const Event& event) :
1219     TopPlusXCandidates(event),
1220     wpGenParticlesSet(false)
1221     {
1222     }
1223    
1224    
1225     TwoNonResTops::~TwoNonResTops()
1226     {
1227     }
1228    
1229     double TwoNonResTops::getGlobalChi2() const { // No global chi2 for this class
1230     return (0.0);
1231     }
1232    
1233    
1234     void TwoNonResTops::recoNonResTops(LeptonPointer lepton) {
1235     if (goodJets.size() < 4)
1236     throw ReconstructionException("Not enough jets available to reconstruct two tops");
1237     // const JetCollection &topJets = goodBJets;
1238     const JetCollection &topJets = goodJets;
1239     leptonFromW = lepton;
1240     selectedNeutrino = 0;
1241     currentSelectedNeutrino = 0;
1242     reconstructNeutrinos();
1243     double chosen_TopMassDifference(9999999.);
1244     double chosen_Chi2Total(9999999.);
1245    
1246     for (unsigned short hadBindex = 0; hadBindex < topJets.size(); ++hadBindex) {
1247     for (unsigned short lepBindex = 0; lepBindex < topJets.size(); ++lepBindex) {
1248     if (lepBindex == hadBindex)
1249     continue;
1250     for (unsigned short jet1Index = 0; jet1Index < goodJets.size(); ++jet1Index) {
1251     if (goodJets.at(jet1Index) == topJets.at(lepBindex) ||
1252     goodJets.at(jet1Index) == topJets.at(hadBindex))
1253     continue;
1254     for (unsigned short jet2Index = 0; jet2Index < goodJets.size(); ++jet2Index) {
1255     if (jet2Index == jet1Index ||
1256     goodJets.at(jet2Index) == topJets.at(lepBindex) ||
1257     goodJets.at(jet2Index) == topJets.at(hadBindex))
1258     continue;
1259     hadronicBJet = topJets.at(hadBindex);
1260     leptonicBJet = topJets.at(lepBindex);
1261     jet1FromW = goodJets.at(jet1Index);
1262     jet2FromW = goodJets.at(jet2Index);
1263    
1264     leptonicW1 = ParticlePointer(new Particle(*neutrino1 + *leptonFromW));
1265     leptonicW2 = ParticlePointer(new Particle(*neutrino2 + *leptonFromW));
1266     hadronicW = ParticlePointer(new Particle(*jet1FromW + *jet2FromW));
1267     leptonicTop1 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW1));
1268     leptonicTop2 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW2));
1269     hadronicTop = ParticlePointer(new Particle(*hadronicBJet + *hadronicW));
1270     fillHypotheses();
1271     selectNeutrinoSolution();
1272     double TopMassDifference = calculateTopMassDifference(currentSelectedNeutrino);
1273     double chi2 = getTotalChi2(currentSelectedNeutrino);
1274     switch (usedTTbarReconstruction) {
1275     case TTbarReconstructionCriterion::TopMassDifference:
1276     if (TopMassDifference < chosen_TopMassDifference) {
1277     hadronicBIndex = hadBindex;
1278     leptonicBIndex = lepBindex;
1279     jet1FromWIndex = jet1Index;
1280     jet2FromWIndex = jet2Index;
1281     chosen_TopMassDifference = TopMassDifference;
1282     selectedNeutrino = currentSelectedNeutrino;
1283     }
1284     break;
1285     case TTbarReconstructionCriterion::chi2:
1286     if (chi2 < chosen_Chi2Total) {
1287     hadronicBIndex = hadBindex;
1288     leptonicBIndex = lepBindex;
1289     jet1FromWIndex = jet1Index;
1290     jet2FromWIndex = jet2Index;
1291     chosen_Chi2Total = chi2;
1292     selectedNeutrino = currentSelectedNeutrino;
1293     }
1294     break;
1295     }
1296     }
1297     }
1298     }
1299     }
1300     std::sort(solutions.begin(), solutions.end(), compareSolutions);
1301     hadronicBJet = topJets.at(hadronicBIndex);
1302     leptonicBJet = topJets.at(leptonicBIndex);
1303     jet1FromW = goodJets.at(jet1FromWIndex);
1304     jet2FromW = goodJets.at(jet2FromWIndex);
1305     leptonicW1 = ParticlePointer(new Particle(*neutrino1 + *leptonFromW));
1306     leptonicW2 = ParticlePointer(new Particle(*neutrino2 + *leptonFromW));
1307     hadronicW = ParticlePointer(new Particle(*jet1FromW + *jet2FromW));
1308     leptonicTop1 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW1));
1309     leptonicTop2 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW2));
1310     hadronicTop = ParticlePointer(new Particle(*hadronicBJet + *hadronicW));
1311     if (selectedNeutrino == 1) // Resonance may be real of fake -- can exclude ttbar?
1312     ttbarResonance = ParticlePointer(new Particle(*leptonicTop1 + *hadronicTop));
1313     else
1314     ttbarResonance = ParticlePointer(new Particle(*leptonicTop2 + *hadronicTop));
1315     doneReconstruction = true;
1316     }
1317    
1318    
1319     double TwoNonResTops::getTotalChi2DJ(unsigned short neutrinoSolution) const
1320     {
1321     double chi2 = TopPairEventCandidate::getTotalChi2(neutrinoSolution);
1322     double pt = dJetFromWp->pt();
1323     if (pt < jet1FromW->pt())
1324     chi2 += (0.5/sqrt(2.0));
1325     if (pt < jet2FromW->pt())
1326     chi2 += (0.5/sqrt(2.0));
1327     double deltaphi = dJetFromWp->deltaPhi(hadronicBJet);
1328     double deltaphi2 = dJetFromWp->deltaPhi(leptonicBJet);
1329     if (deltaphi2 > deltaphi)
1330     deltaphi = deltaphi2;
1331     double dphichi = deltaphi - pi;
1332     chi2 += (1.0/sqrt(2.0)) * (dphichi * dphichi)/(pi * pi/2.0);
1333     if (dJetFromWp->isBJet())
1334     chi2 += (0.5/sqrt(2.0));
1335     return (chi2);
1336     }
1337    
1338    
1339     void TwoNonResTops::recoNonResTopsWdjet(LeptonPointer lepton) {
1340     if (goodJets.size() < 5)
1341     throw ReconstructionException("Not enough jets available to reconstruct two tops");
1342     const JetCollection &topJets = goodJets;
1343     leptonFromW = lepton;
1344     selectedNeutrino = 0;
1345     currentSelectedNeutrino = 0;
1346     reconstructNeutrinos();
1347     double chosen_TopMassDifference(9999999.);
1348     double chosen_Chi2Total(9999999.);
1349    
1350     for (unsigned short hadBindex = 0; hadBindex < topJets.size(); ++hadBindex) {
1351     for (unsigned short lepBindex = 0; lepBindex < topJets.size(); ++lepBindex) {
1352     if (lepBindex == hadBindex)
1353     continue;
1354     for (unsigned short jet1Index = 0; jet1Index < goodJets.size(); ++jet1Index) {
1355     if (goodJets.at(jet1Index) == topJets.at(lepBindex) ||
1356     goodJets.at(jet1Index) == topJets.at(hadBindex))
1357     continue;
1358     for (unsigned short jet2Index = 0; jet2Index < goodJets.size(); ++jet2Index) {
1359     if (jet2Index == jet1Index ||
1360     goodJets.at(jet2Index) == topJets.at(lepBindex) ||
1361     goodJets.at(jet2Index) == topJets.at(hadBindex))
1362     continue;
1363     for (unsigned short djetIndex = 0; djetIndex < goodJets.size(); ++djetIndex) {
1364     if (djetIndex == jet1Index || jet2Index == djetIndex ||
1365     goodJets.at(djetIndex) == topJets.at(lepBindex) ||
1366     goodJets.at(djetIndex) == topJets.at(hadBindex))
1367     continue;
1368     hadronicBJet = topJets.at(hadBindex);
1369     leptonicBJet = topJets.at(lepBindex);
1370     jet1FromW = goodJets.at(jet1Index);
1371     jet2FromW = goodJets.at(jet2Index);
1372     dJetFromWp = goodJets.at(djetIndex);
1373    
1374     leptonicW1 = ParticlePointer(new Particle(*neutrino1 + *leptonFromW));
1375     leptonicW2 = ParticlePointer(new Particle(*neutrino2 + *leptonFromW));
1376     hadronicW = ParticlePointer(new Particle(*jet1FromW + *jet2FromW));
1377     leptonicTop1 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW1));
1378     leptonicTop2 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW2));
1379     hadronicTop = ParticlePointer(new Particle(*hadronicBJet + *hadronicW));
1380     fillHypotheses();
1381     selectNeutrinoSolution();
1382     double TopMassDifference = calculateTopMassDifference(currentSelectedNeutrino);
1383     double chi2 = getTotalChi2DJ(currentSelectedNeutrino);
1384     switch (usedTTbarReconstruction) {
1385     case TTbarReconstructionCriterion::TopMassDifference:
1386     if (TopMassDifference < chosen_TopMassDifference) {
1387     hadronicBIndex = hadBindex;
1388     leptonicBIndex = lepBindex;
1389     jet1FromWIndex = jet1Index;
1390     jet2FromWIndex = jet2Index;
1391     chosen_TopMassDifference = TopMassDifference;
1392     selectedNeutrino = currentSelectedNeutrino;
1393     }
1394     break;
1395     case TTbarReconstructionCriterion::chi2:
1396     if (chi2 < chosen_Chi2Total) {
1397     hadronicBIndex = hadBindex;
1398     leptonicBIndex = lepBindex;
1399     jet1FromWIndex = jet1Index;
1400     jet2FromWIndex = jet2Index;
1401     dJetFromWpIndex = djetIndex;
1402     chosen_Chi2Total = chi2;
1403     selectedNeutrino = currentSelectedNeutrino;
1404     }
1405     break;
1406     }
1407     }
1408     }
1409     }
1410     }
1411     }
1412     std::sort(solutions.begin(), solutions.end(), compareSolutions);
1413     hadronicBJet = topJets.at(hadronicBIndex);
1414     leptonicBJet = topJets.at(leptonicBIndex);
1415     jet1FromW = goodJets.at(jet1FromWIndex);
1416     jet2FromW = goodJets.at(jet2FromWIndex);
1417     dJetFromWp = goodJets.at(dJetFromWpIndex);
1418     leptonicW1 = ParticlePointer(new Particle(*neutrino1 + *leptonFromW));
1419     leptonicW2 = ParticlePointer(new Particle(*neutrino2 + *leptonFromW));
1420     hadronicW = ParticlePointer(new Particle(*jet1FromW + *jet2FromW));
1421     leptonicTop1 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW1));
1422     leptonicTop2 = ParticlePointer(new Particle(*leptonicBJet + *leptonicW2));
1423     hadronicTop = ParticlePointer(new Particle(*hadronicBJet + *hadronicW));
1424     if (selectedNeutrino == 1) // Resonance may be real of fake -- can exclude ttbar?
1425     ttbarResonance = ParticlePointer(new Particle(*leptonicTop1 + *hadronicTop));
1426     else
1427     ttbarResonance = ParticlePointer(new Particle(*leptonicTop2 + *hadronicTop));
1428     doneReconstruction = true;
1429     }
1430    
1431    
1432     const JetPointer TwoNonResTops::getLtNonTopJet() const {
1433     throwExpeptionIfNotReconstructed("TwoNonResTops::getLtNonTopJet");
1434     if (goodJets.size() < 5)
1435     throw ReconstructionException("Not enough jets available to reconstruct Wprime");
1436     int ltJetInd = -1, bJetInd = -1;
1437     unsigned int index = 0;
1438     while (index < goodJets.size() && ltJetInd < 0) {
1439     JetPointer currJet = goodJets.at(index);
1440     if (currJet->isBJet() == false && currJet != leptonicBJet &&
1441     currJet != hadronicBJet && currJet != jet1FromW && currJet != jet2FromW) {
1442     ltJetInd = index;
1443     } else if (bJetInd < 0 && currJet != leptonicBJet && currJet != hadronicBJet &&
1444     currJet != jet1FromW && currJet != jet2FromW) {
1445     bJetInd = index;
1446     }
1447     ++index;
1448     }
1449     if (ltJetInd < 0 && bJetInd >= 0) // Take a B jet if no others found.
1450     ltJetInd = bJetInd;
1451     if (ltJetInd >= 0)
1452     return (goodJets.at(ltJetInd));
1453     throw ReconstructionException(TString("getLtNonTopJet:Fifth jet missing"));
1454     }
1455    
1456    
1457     const ParticlePointer TwoNonResTops::getWprime() const {
1458     throwExpeptionIfNotReconstructed("TwoNonResTops::getWprime");
1459     ParticlePointer ltJet = getLtNonTopJet();
1460     ParticlePointer hiPtTop = getLeptonicTop();
1461    
1462     // Generator study indicates lower Pt top is usually spectator
1463     if (hiPtTop->pt() < hadronicTop->pt())
1464     hiPtTop = hadronicTop;
1465     return (ParticlePointer(new Particle(*hiPtTop + *ltJet)));
1466     }
1467    
1468    
1469     const ParticlePointer TwoNonResTops::getNegWprime() const {
1470     throwExpeptionIfNotReconstructed("TwoNonResTops::getNegWprime");
1471     ParticlePointer ltJet = getLtNonTopJet();
1472     ParticlePointer negTop = hadronicTop;
1473     if (leptonFromW->charge() < 0)
1474     negTop = getLeptonicTop();
1475    
1476     return (ParticlePointer(new Particle(*negTop + *ltJet)));
1477     // return (ParticlePointer(new Particle(*negTop + *dJetFromWp)));
1478     }
1479    
1480    
1481     // Deliberately reco the incorrect W' to see if ttbar still looks the same
1482     const ParticlePointer TwoNonResTops::getNegWprimeBad() const {
1483     throwExpeptionIfNotReconstructed("TwoNonResTops::getNegWprimeBad");
1484     ParticlePointer ltJet = getLtNonTopJet();
1485     ParticlePointer negTop = hadronicTop;
1486     if (leptonFromW->charge() > 0)
1487     negTop = getLeptonicTop(); // Use wrong top
1488    
1489     return (ParticlePointer(new Particle(*negTop + *ltJet)));
1490     // return (ParticlePointer(new Particle(*negTop + *dJetFromWp)));
1491     }
1492    
1493    
1494     const ParticlePointer TwoNonResTops::getBtBWprime() const {
1495     throwExpeptionIfNotReconstructed("TwoNonResTops::getBtBWprim");
1496     ParticlePointer ltJet = getLtNonTopJet();
1497     ParticlePointer hiAnglTop = getLeptonicTop();
1498     if (getHadTopDAngle() > hiAnglTop->angle(ltJet))
1499     hiAnglTop = hadronicTop;
1500     return (ParticlePointer(new Particle(*hiAnglTop + *ltJet)));
1501     }
1502    
1503    
1504     double TwoNonResTops::getHadTopDAngle() const
1505     {
1506     ParticlePointer jet = getLtNonTopJet();
1507     return (hadronicTop->angle(jet));
1508     }
1509    
1510    
1511     void TwoNonResTops::recoWpFrom1Top(LeptonPointer lepton) {
1512     if (goodJets.size() < 5)
1513     throw
1514     ReconstructionException("Not enough jets available to reconstruct W'");
1515     ParticlePointer negTop, posTop;
1516     if (lepton->charge() < 0) {
1517     // findDJet();
1518     recoBestLeptoTop(lepton);
1519     // findDJet(leptonicBJet);
1520     recoHadronicTop(lepton);
1521     // recoHadronicTop(lepton, leptonicBJet);
1522     negTop = getLeptonicTop();
1523     posTop = hadronicTop;
1524     } else {
1525     // findDJet();
1526     recoHadronicTop(lepton);
1527     // findDJet(hadronicBJet);
1528     recoBestLeptoTop(lepton);
1529     // recoBestLeptoTop(lepton, hadronicBJet, jet1FromW, jet2FromW);
1530     negTop = hadronicTop;
1531     posTop = getLeptonicTop();
1532     }
1533     JetPointer ltJet = getLtNonTopJet();
1534     wpFrom1Top = ParticlePointer(new Particle(*negTop + *ltJet));
1535     wpBadFrom1Top = ParticlePointer(new Particle(*posTop + *ltJet));
1536     topFromWp = negTop;
1537     dJetFromWp = ltJet;
1538     // wpFrom1Top = ParticlePointer(new Particle(*negTop + *dJetFromWp));
1539     // wpBadFrom1Top = ParticlePointer(new Particle(*posTop + *dJetFromWp));
1540     }
1541    
1542    
1543     void TwoNonResTops::getTrueWpChain()
1544     {
1545     const MCParticleCollection &mcColl = genParticles;
1546     const int nullInd = -99;
1547     int wpInd = nullInd, wpTopInd = nullInd, specTopInd = nullInd, dQInd = nullInd;
1548     int wpTopBInd = nullInd, specTopBInd = nullInd;
1549     if (mcColl.size() > 0) {
1550     for (unsigned int i = 0; i < mcColl.size(); ++i) {
1551     if (fabs(mcColl.at(i)->pdgId()) == WPRIME_ID) {
1552     if (wpInd == nullInd)
1553     wpInd = i;
1554     else if (mcColl.at(i)->motherIndex() != wpInd)
1555     cout << "Stray W': " << i << " mother " << mcColl.at(i)->motherIndex() << endl;
1556     }
1557     }
1558     for (unsigned int i = 0; i < mcColl.size(); ++i) {
1559     if (fabs(mcColl.at(i)->pdgId()) == TOP_ID) {
1560     if (mcColl.at(i)->motherIndex() == wpInd) {
1561     if (wpTopInd == nullInd)
1562     wpTopInd = i;
1563     else cout << "Additional W' top: " << i << endl;
1564     } else if (specTopInd == nullInd)
1565     specTopInd = i;
1566     else cout << "Additional spec top: " << i << endl;
1567     } else if (fabs(mcColl.at(i)->pdgId()) == DQ_ID) {
1568     if (mcColl.at(i)->motherIndex() == wpInd) {
1569     if (dQInd == nullInd)
1570     dQInd = i;
1571     else cout << "Additional d daughter: " << i << endl;
1572     }
1573     }
1574     }
1575     if (wpTopInd >= 0 && specTopInd >= 0) {
1576     int wptprevPdgId = 0;
1577     int sptprevPdgId = 0;
1578     for (unsigned int i = 0; i < mcColl.size(); ++i) {
1579     int pdgId = fabs(mcColl.at(i)->pdgId());
1580     if (mcColl.at(i)->motherIndex() == wpTopInd ||
1581     mcColl.at(i)->motherIndex() == wpInd) {
1582     if (mcColl.at(i)->isQuark() && pdgId % 2 == 1) {
1583     if (wpTopBInd == nullInd || pdgId > wptprevPdgId) {
1584     wpTopBInd = i;
1585     wptprevPdgId = pdgId;
1586     } else if (pdgId == BQ_ID)
1587     cout << "Additional W' top b quark: " << i << " mother " << mcColl.at(i)->motherIndex() << endl;
1588     }
1589     } else if (mcColl.at(i)->motherIndex() == specTopInd ||
1590     mcColl.at(i)->motherIndex() == mcColl.at(specTopInd)->motherIndex()) {
1591     if (mcColl.at(i)->isQuark() && pdgId % 2 == 1) {
1592     if (specTopBInd == nullInd || pdgId > sptprevPdgId) {
1593     specTopBInd = i;
1594     sptprevPdgId = pdgId;
1595     } else if (pdgId == BQ_ID)
1596     cout << "Additional spec top b quark: " << i << " mother " << mcColl.at(i)->motherIndex() << endl;
1597     }
1598     } else if (pdgId == BQ_ID && fabs(mcColl.at(mcColl.at(i)->motherIndex())->pdgId()) != BQ_ID) {
1599     // cout << "Additional b quark: " << i << " mother " << mcColl.at(i)->motherIndex();
1600     // cout << " ID " << mcColl.at(mcColl.at(i)->motherIndex())->pdgId() << endl;
1601     }
1602     }
1603     }
1604     if (wpInd >= 0 && wpTopInd >= 0 && specTopInd >= 0 && dQInd >= 0 &&
1605     wpTopBInd >= 0 && specTopBInd >= 0) {
1606     wpTru = mcColl.at(wpInd);
1607     wpTopTru = mcColl.at(wpTopInd);
1608     specTopTru = mcColl.at(specTopInd);
1609     wpTopBTru = mcColl.at(wpTopBInd);
1610     specTopBTru = mcColl.at(specTopBInd);
1611     wpDDauTru = mcColl.at(dQInd);
1612     wpGenParticlesSet = true;
1613     } else {
1614     // cout << "#11 : " << mcColl.at(11)->pdgId() << " mother " << mcColl.at(11)->motherIndex();
1615     // cout << " ID " << mcColl.at(mcColl.at(11)->motherIndex())->pdgId() << endl;
1616     cout << "wp wpTop specTop dq wpTopB specTopB " << wpInd << " " << wpTopInd << " " << specTopInd << " " << dQInd;
1617     cout << " " << wpTopBInd << " " << specTopBInd << endl;
1618     }
1619     }
1620     }
1621    
1622    
1623     MCParticlePointer TwoNonResTops::getWpTru() const
1624     {
1625     if (wpGenParticlesSet == false)
1626     throw ReconstructionException("Need to run TwoNonResTops::getTrueWpChain first");
1627     return (wpTru);
1628     }
1629    
1630     MCParticlePointer TwoNonResTops::getWpTopTru() const
1631     {
1632     if (wpGenParticlesSet == false)
1633     throw ReconstructionException("Need to run TwoNonResTops::getTrueWpChain first");
1634     return (wpTopTru);
1635     }
1636    
1637     MCParticlePointer TwoNonResTops::getSpecTop() const
1638     {
1639     if (wpGenParticlesSet == false)
1640     throw ReconstructionException("Need to run TwoNonResTops::getTrueWpChain first");
1641     return (specTopTru);
1642     }
1643    
1644     MCParticlePointer TwoNonResTops::getWpTopBTru() const
1645     {
1646     if (wpGenParticlesSet == false)
1647     throw ReconstructionException("Need to run TwoNonResTops::getTrueWpChain first");
1648     return (wpTopBTru);
1649     }
1650    
1651     MCParticlePointer TwoNonResTops::getSpecTopB() const
1652     {
1653     if (wpGenParticlesSet == false)
1654     throw ReconstructionException("Need to run TwoNonResTops::getTrueWpChain first");
1655     return (specTopBTru);
1656     }
1657    
1658     MCParticlePointer TwoNonResTops::getWpDDauTru() const
1659     {
1660     if (wpGenParticlesSet == false)
1661     throw ReconstructionException("Need to run TwoNonResTops::getTrueWpChain first");
1662     return (wpDDauTru);
1663     }
1664    
1665    
1666     TopNoMassConstraint::TopNoMassConstraint()
1667     {
1668     }
1669    
1670    
1671     TopNoMassConstraint::TopNoMassConstraint(const Event& event) :
1672     ToplikeCandidate(event)
1673     {
1674     }
1675    
1676    
1677     TopNoMassConstraint::~TopNoMassConstraint()
1678     {
1679     }
1680    
1681    
1682     double TopNoMassConstraint::getLoneHadChi2() const
1683     {
1684     double ptRatioDifference = TMath::Power(PtRatio() - matched_ptratio, 2);
1685     double ptRatioError = 2 * matched_ptratio_sigma * matched_ptratio_sigma;
1686     double ptRatioTerm = ptRatioDifference / ptRatioError;
1687    
1688     double WmassDifference = TMath::Power(hadronicW->mass() - pdgWmass, 2);
1689     double WmassError = 2 * wMassWidth * wMassWidth;
1690     double WmassTerm = WmassDifference / WmassError;
1691     return 1.0 / sqrt(2.0) * (WmassTerm + ptRatioTerm);
1692     }
1693    
1694     double TopNoMassConstraint::getLepBAngle() const
1695     {
1696     return (leptonicBJet->angle(leptonFromW));
1697     }
1698    
1699    
1700     LeptoTopNoMassConstraint::LeptoTopNoMassConstraint()
1701     {
1702     }
1703    
1704    
1705     LeptoTopNoMassConstraint::LeptoTopNoMassConstraint(const Event& event) :
1706     TopNoMassConstraint(event)
1707     {
1708     }
1709    
1710    
1711     LeptoTopNoMassConstraint::~LeptoTopNoMassConstraint()
1712     {
1713     }
1714    
1715    
1716     /*
1717     Tried for neutrino selection, but getGlobalChi2 still uses hadronic top quantities.
1718     double LeptoTopNoMassConstraint::getTotalChi2(unsigned short int neutrinoSolution) const {
1719     return TopPairEventCandidate::getLeptonicChi2(neutrinoSolution) + TopPairEventCandidate::getGlobalChi2(neutrinoSolution);
1720     }
1721     */
1722    
1723     void LeptoTopNoMassConstraint::inspectReconstructedEvent() const {
1724     cout << "run " << runNumber << ", event " << eventNumber << endl;
1725     cout << "leptonic b jet, goodJet index " << leptonicBIndex << endl;
1726     EventContentPrinter::printJet(leptonicBJet);
1727    
1728     cout << "lepton from W" << endl;
1729     // EventContentPrinter::printElectron(electronFromW);
1730    
1731     cout << "MET" << endl;
1732     EventContentPrinter::printParticle(met);
1733     cout << endl;
1734    
1735     cout << "reconstructed neutrino 1(selected: " << selectedNeutrino << ")" << endl;
1736     EventContentPrinter::printParticle(neutrino1);
1737     cout << endl;
1738    
1739     cout << "reconstructed neutrino 2(selected: " << selectedNeutrino << ")" << endl;
1740     EventContentPrinter::printParticle(neutrino2);
1741     cout << endl;
1742    
1743     cout << "leptonic W 1 (selected: " << selectedNeutrino << ")" << endl;
1744     EventContentPrinter::printParticle(leptonicW1);
1745     cout << endl;
1746    
1747     cout << "leptonic W 2 (selected: " << selectedNeutrino << ")" << endl;
1748     EventContentPrinter::printParticle(leptonicW2);
1749     cout << endl;
1750    
1751     cout << "leptonic top 1 (selected: " << selectedNeutrino << ")" << endl;
1752     EventContentPrinter::printParticle(leptonicTop1);
1753     cout << endl;
1754    
1755     cout << "leptonic top 2 (selected: " << selectedNeutrino << ")" << endl;
1756     EventContentPrinter::printParticle(leptonicTop2);
1757     cout << endl;
1758    
1759     cout << endl;
1760     }
1761    
1762    
1763     }