1#include "Rivet/Analysis.hh"
  2#include "Rivet/Math/LorentzTrans.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5
  6namespace Rivet {
  7
  8
  9  namespace { //< only visible in this compilation unit
 10
 11    /// @brief Pseudo top finder
 12    ///
 13    /// Find top quark in the particle level.
 14    /// The definition is based on the agreement at the LHC working group.
 15    class PseudoTop : public FinalState {
 16    public:
 17
 18      /// @name Standard constructors and destructors.
 19      /// @{
 20
 21      /// The default constructor. May specify the minimum and maximum
 22      /// pseudorapidity \f$ \eta \f$ and the min \f$ p_T \f$ (in GeV).
 23      PseudoTop(double lepR = 0.1, double lepMinPt = 20, double lepMaxEta = 2.4,
 24                double jetR = 0.4, double jetMinPt = 30, double jetMaxEta = 4.7)
 25        : FinalState(),
 26          _lepR(lepR), _lepMinPt(lepMinPt), _lepMaxEta(lepMaxEta),
 27          _jetR(jetR), _jetMinPt(jetMinPt), _jetMaxEta(jetMaxEta)
 28      {
 29        setName("PseudoTop");
 30      }
 31
 32      enum TTbarMode {CH_NONE=-1, CH_FULLHADRON = 0, CH_SEMILEPTON, CH_FULLLEPTON};
 33      enum DecayMode {CH_HADRON = 0, CH_MUON, CH_ELECTRON};
 34
 35      TTbarMode mode() const {
 36        if (!_isValid) return CH_NONE;
 37        if (_mode1 == CH_HADRON && _mode2 == CH_HADRON) return CH_FULLHADRON;
 38        else if ( _mode1 != CH_HADRON && _mode2 != CH_HADRON) return CH_FULLLEPTON;
 39        else return CH_SEMILEPTON;
 40      }
 41      virtual const DecayMode& mode1() const {return _mode1;}
 42      virtual const DecayMode& mode2() const {return _mode2;}
 43
 44      /// Clone on the heap.
 45      RIVET_DEFAULT_PROJ_CLONE(PseudoTop);
 46
 47      /// @}
 48
 49      /// Import to avoid warnings about overload-hiding
 50      using Projection::operator =;
 51
 52    public:
 53
 54      virtual const Particle& t1() const {return _t1;}
 55      virtual const Particle& t2() const {return _t2;}
 56      virtual const Particle& b1() const {return _b1;}
 57      virtual const Particle& b2() const {return _b2;}
 58      virtual const Particles& wDecays1() const {return _wDecays1;}
 59      virtual const Particles& wDecays2() const {return _wDecays2;}
 60      virtual const Jets& jets() const {return _jets;}
 61      virtual const Jets& bjets() const {return _bjets;}
 62      virtual const Jets& ljets() const {return _ljets;}
 63
 64      // Apply the projection to the event
 65      void project(const Event& e); // override; ///< @todo Re-enable when C++11 allowed
 66      void cleanup(std::map<double, std::pair<size_t, size_t> >& v, const bool doCrossCleanup=false) const;
 67      CmpState compare(const Projection& p) const;
 68
 69
 70    private:
 71      const double _lepR, _lepMinPt, _lepMaxEta;
 72      const double _jetR, _jetMinPt, _jetMaxEta;
 73
 74      //constexpr ///< @todo Re-enable when C++11 allowed
 75      static double _tMass; // = 172.5*GeV; ///< @todo Re-enable when C++11 allowed
 76      //constexpr ///< @todo Re-enable when C++11 allowed
 77      static double _wMass; // = 80.4*GeV; ///< @todo Re-enable when C++11 allowed
 78
 79    private:
 80      bool _isValid;
 81      DecayMode _mode1, _mode2;
 82
 83      Particle _t1, _t2;
 84      Particle _b1, _b2;
 85      Particles _wDecays1, _wDecays2;
 86      Jets _jets, _bjets, _ljets;
 87
 88    };
 89
 90    // More implementation below the analysis code
 91
 92  }
 93
 94
 95
 96  /// Pseudo-top analysis from CMS
 97  class CMS_2015_I1370682 : public Analysis {
 98  public:
 99
100    CMS_2015_I1370682()
101      : Analysis("CMS_2015_I1370682"),
102        _applyCorrection(true),
103        _doShapeOnly(true)
104    {    }
105
106
107    void init() {
108      declare(PseudoTop(0.1, 20, 2.4, 0.5, 30, 2.4), "ttbar");
109
110      // Lepton + Jet channel
111      book(_hSL_topPt         ,"d15-x01-y01"); // 1/sigma dsigma/dpt(top)
112      book(_hSL_topPtTtbarSys ,"d16-x01-y01"); // 1/sigma dsigma/dpt*(top)
113      book(_hSL_topY          ,"d17-x01-y01"); // 1/sigma dsigma/dy(top)
114      book(_hSL_ttbarDelPhi   ,"d18-x01-y01"); // 1/sigma dsigma/ddeltaphi(t,tbar)
115      book(_hSL_topPtLead     ,"d19-x01-y01"); // 1/sigma dsigma/dpt(t1)
116      book(_hSL_topPtSubLead  ,"d20-x01-y01"); // 1/sigma dsigma/dpt(t2)
117      book(_hSL_ttbarPt       ,"d21-x01-y01"); // 1/sigma dsigma/dpt(ttbar)
118      book(_hSL_ttbarY        ,"d22-x01-y01"); // 1/sigma dsigma/dy(ttbar)
119      book(_hSL_ttbarMass     ,"d23-x01-y01"); // 1/sigma dsigma/dm(ttbar)
120
121      // Dilepton channel
122      book(_hDL_topPt         ,"d24-x01-y01"); // 1/sigma dsigma/dpt(top)
123      book(_hDL_topPtTtbarSys ,"d25-x01-y01"); // 1/sigma dsigma/dpt*(top)
124      book(_hDL_topY          ,"d26-x01-y01"); // 1/sigma dsigma/dy(top)
125      book(_hDL_ttbarDelPhi   ,"d27-x01-y01"); // 1/sigma dsigma/ddeltaphi(t,tbar)
126      book(_hDL_topPtLead     ,"d28-x01-y01"); // 1/sigma dsigma/dpt(t1)
127      book(_hDL_topPtSubLead  ,"d29-x01-y01"); // 1/sigma dsigma/dpt(t2)
128      book(_hDL_ttbarPt       ,"d30-x01-y01"); // 1/sigma dsigma/dpt(ttbar)
129      book(_hDL_ttbarY        ,"d31-x01-y01"); // 1/sigma dsigma/dy(ttbar)
130      book(_hDL_ttbarMass     ,"d32-x01-y01"); // 1/sigma dsigma/dm(ttbar)
131
132    }
133
134
135    void analyze(const Event& event) {
136
137      // Get the ttbar candidate
138      const PseudoTop& ttbar = apply<PseudoTop>(event, "ttbar");
139      if ( ttbar.mode() == PseudoTop::CH_NONE ) vetoEvent;
140
141      const FourMomentum& t1P4 = ttbar.t1().momentum();
142      const FourMomentum& t2P4 = ttbar.t2().momentum();
143      const double pt1 = std::max(t1P4.pT(), t2P4.pT());
144      const double pt2 = std::min(t1P4.pT(), t2P4.pT());
145      const double dPhi = deltaPhi(t1P4, t2P4);
146      const FourMomentum ttP4 = t1P4 + t2P4;
147      const FourMomentum t1P4AtCM = LorentzTransform::mkFrameTransformFromBeta(ttP4.betaVec()).transform(t1P4);
148
149      if ( ttbar.mode() == PseudoTop::CH_SEMILEPTON ) {
150        const Particle lCand1 = ttbar.wDecays1()[0]; // w1 dau0 is the lepton in the PseudoTop
151        if (lCand1.pT() < 33*GeV || lCand1.abseta() > 2.1) vetoEvent;
152        _hSL_topPt->fill(t1P4.pT());
153        _hSL_topPt->fill(t2P4.pT());
154        _hSL_topPtTtbarSys->fill(t1P4AtCM.pT());
155        _hSL_topY->fill(t1P4.rapidity());
156        _hSL_topY->fill(t2P4.rapidity());
157        _hSL_ttbarDelPhi->fill(dPhi);
158        _hSL_topPtLead->fill(pt1);
159        _hSL_topPtSubLead->fill(pt2);
160        _hSL_ttbarPt->fill(ttP4.pT());
161        _hSL_ttbarY->fill(ttP4.rapidity());
162        _hSL_ttbarMass->fill(ttP4.mass());
163      }
164      else if ( ttbar.mode() == PseudoTop::CH_FULLLEPTON ) {
165        const Particle lCand1 = ttbar.wDecays1()[0]; // dau0 are the lepton in the PseudoTop
166        const Particle lCand2 = ttbar.wDecays2()[0]; // dau0 are the lepton in the PseudoTop
167        if (lCand1.pT() < 20*GeV || lCand1.abseta() > 2.4) vetoEvent;
168        if (lCand2.pT() < 20*GeV || lCand2.abseta() > 2.4) vetoEvent;
169        _hDL_topPt->fill(t1P4.pT());
170        _hDL_topPt->fill(t2P4.pT());
171        _hDL_topPtTtbarSys->fill(t1P4AtCM.pT());
172        _hDL_topY->fill(t1P4.rapidity());
173        _hDL_topY->fill(t2P4.rapidity());
174        _hDL_ttbarDelPhi->fill(dPhi);
175        _hDL_topPtLead->fill(pt1);
176        _hDL_topPtSubLead->fill(pt2);
177        _hDL_ttbarPt->fill(ttP4.pT());
178        _hDL_ttbarY->fill(ttP4.rapidity());
179        _hDL_ttbarMass->fill(ttP4.mass());
180      }
181
182    }
183
184
185    void finalize() {
186      if ( _applyCorrection ) {
187        // Correction functions for TOP-12-028 paper, (parton bin height)/(pseudotop bin height)
188        const double ch15[] = { 5.473609, 4.941048, 4.173346, 3.391191, 2.785644, 2.371346, 2.194161, 2.197167, };
189        const double ch16[] = { 5.470905, 4.948201, 4.081982, 3.225532, 2.617519, 2.239217, 2.127878, 2.185918, };
190        const double ch17[] = { 10.003667, 4.546519, 3.828115, 3.601018, 3.522194, 3.524694, 3.600951, 3.808553, 4.531891, 9.995370, };
191        const double ch18[] = { 4.406683, 4.054041, 3.885393, 4.213646, };
192        const double ch19[] = { 6.182537, 5.257703, 4.422280, 3.568402, 2.889408, 2.415878, 2.189974, 2.173210, };
193        const double ch20[] = { 5.199874, 4.693318, 3.902882, 3.143785, 2.607877, 2.280189, 2.204124, 2.260829, };
194        const double ch21[] = { 6.053523, 3.777506, 3.562251, 3.601356, 3.569347, 3.410472, };
195        const double ch22[] = { 11.932351, 4.803773, 3.782709, 3.390775, 3.226806, 3.218982, 3.382678, 3.773653, 4.788191, 11.905338, };
196        const double ch23[] = { 7.145255, 5.637595, 4.049882, 3.025917, 2.326430, 1.773824, 1.235329, };
197
198        const double ch24[] = { 2.268193, 2.372063, 2.323975, 2.034655, 1.736793, };
199        const double ch25[] = { 2.231852, 2.383086, 2.341894, 2.031318, 1.729672, 1.486993, };
200        const double ch26[] = { 3.993526, 2.308249, 2.075136, 2.038297, 2.036302, 2.078270, 2.295817, 4.017713, };
201        const double ch27[] = { 2.205978, 2.175010, 2.215376, 2.473144, };
202        const double ch28[] = { 2.321077, 2.371895, 2.338871, 2.057821, 1.755382, };
203        const double ch29[] = { 2.222707, 2.372591, 2.301688, 1.991162, 1.695343, };
204        const double ch30[] = { 2.599677, 2.026855, 2.138620, 2.229553, };
205        const double ch31[] = { 5.791779, 2.636219, 2.103642, 1.967198, 1.962168, 2.096514, 2.641189, 5.780828, };
206        const double ch32[] = { 2.006685, 2.545525, 2.477745, 2.335747, 2.194226, 2.076500, };
207
208        applyCorrection(_hSL_topPt, ch15);
209        applyCorrection(_hSL_topPtTtbarSys, ch16);
210        applyCorrection(_hSL_topY, ch17);
211        applyCorrection(_hSL_ttbarDelPhi, ch18);
212        applyCorrection(_hSL_topPtLead, ch19);
213        applyCorrection(_hSL_topPtSubLead, ch20);
214        applyCorrection(_hSL_ttbarPt, ch21);
215        applyCorrection(_hSL_ttbarY, ch22);
216        applyCorrection(_hSL_ttbarMass, ch23);
217
218        applyCorrection(_hDL_topPt, ch24);
219        applyCorrection(_hDL_topPtTtbarSys, ch25);
220        applyCorrection(_hDL_topY, ch26);
221        applyCorrection(_hDL_ttbarDelPhi, ch27);
222        applyCorrection(_hDL_topPtLead, ch28);
223        applyCorrection(_hDL_topPtSubLead, ch29);
224        applyCorrection(_hDL_ttbarPt, ch30);
225        applyCorrection(_hDL_ttbarY, ch31);
226        applyCorrection(_hDL_ttbarMass, ch32);
227      }
228
229      if ( _doShapeOnly ) {
230        normalize(_hSL_topPt        );
231        normalize(_hSL_topPtTtbarSys);
232        normalize(_hSL_topY         );
233        normalize(_hSL_ttbarDelPhi  );
234        normalize(_hSL_topPtLead    );
235        normalize(_hSL_topPtSubLead );
236        normalize(_hSL_ttbarPt      );
237        normalize(_hSL_ttbarY       );
238        normalize(_hSL_ttbarMass    );
239
240        normalize(_hDL_topPt        );
241        normalize(_hDL_topPtTtbarSys);
242        normalize(_hDL_topY         );
243        normalize(_hDL_ttbarDelPhi  );
244        normalize(_hDL_topPtLead    );
245        normalize(_hDL_topPtSubLead );
246        normalize(_hDL_ttbarPt      );
247        normalize(_hDL_ttbarY       );
248        normalize(_hDL_ttbarMass    );
249      }
250      else {
251        const double s = 1./sumOfWeights();
252        scale(_hSL_topPt        , s);
253        scale(_hSL_topPtTtbarSys, s);
254        scale(_hSL_topY         , s);
255        scale(_hSL_ttbarDelPhi  , s);
256        scale(_hSL_topPtLead    , s);
257        scale(_hSL_topPtSubLead , s);
258        scale(_hSL_ttbarPt      , s);
259        scale(_hSL_ttbarY       , s);
260        scale(_hSL_ttbarMass    , s);
261        scale(_hDL_topPt        , s);
262        scale(_hDL_topPtTtbarSys, s);
263        scale(_hDL_topY         , s);
264        scale(_hDL_ttbarDelPhi  , s);
265        scale(_hDL_topPtLead    , s);
266        scale(_hDL_topPtSubLead , s);
267        scale(_hDL_ttbarPt      , s);
268        scale(_hDL_ttbarY       , s);
269        scale(_hDL_ttbarMass    , s);
270      }
271
272    }
273
274
275    void applyCorrection(Histo1DPtr& h, const double* cf) {
276      for (auto& bin : h->bins()) {
277        bin.scaleW( cf[bin.index()-1] );
278      }
279    }
280
281
282  private:
283
284    const bool _applyCorrection, _doShapeOnly;
285    Histo1DPtr _hSL_topPt, _hSL_topPtTtbarSys, _hSL_topY, _hSL_ttbarDelPhi, _hSL_topPtLead,
286      _hSL_topPtSubLead, _hSL_ttbarPt, _hSL_ttbarY, _hSL_ttbarMass;
287    Histo1DPtr _hDL_topPt, _hDL_topPtTtbarSys, _hDL_topY, _hDL_ttbarDelPhi, _hDL_topPtLead,
288      _hDL_topPtSubLead, _hDL_ttbarPt, _hDL_ttbarY, _hDL_ttbarMass;
289
290  };
291
292
293
294  RIVET_DECLARE_PLUGIN(CMS_2015_I1370682);
295
296
297  ///////////////
298
299  // More PseudoTop implementation
300  namespace {
301
302
303    double PseudoTop::_tMass = 172.5*GeV;
304    double PseudoTop::_wMass = 80.4*GeV;
305
306    CmpState PseudoTop::compare(const Projection& p) const {
307      const PCmp fscmp = mkNamedPCmp(p, "FS");
308      if (fscmp != CmpState::EQ) return fscmp;
309
310      const PseudoTop& other = dynamic_cast<const PseudoTop&>(p);
311      CmpState cs_lepR = cmp(_lepR, other._lepR);
312      if (cs_lepR != CmpState::EQ) return cs_lepR;
313
314      CmpState cs_jetR = cmp(_jetR, other._jetR);
315      if (cs_jetR != CmpState::EQ) return cs_jetR;
316
317      CmpState cs_lepMinPt = cmp(_lepMinPt, other._lepMinPt);
318      if (cs_lepMinPt != CmpState::EQ) return cs_lepMinPt;
319
320      CmpState cs_jetMinPt = cmp(_jetMinPt, other._jetMinPt);
321      if (cs_jetMinPt != CmpState::EQ) return cs_jetMinPt;
322
323      CmpState cs_lepMaxEta = cmp(_lepMaxEta, other._lepMaxEta);
324      if (cs_lepMaxEta != CmpState::EQ) return cs_lepMaxEta;
325
326      CmpState cs_jetMaxEta = cmp(_jetMaxEta, other._jetMaxEta);
327      return cs_jetMaxEta;
328    }
329
330    void PseudoTop::cleanup(map<double, pair<size_t, size_t> >& v, const bool doCrossCleanup) const {
331      vector<map<double, pair<size_t, size_t> >::iterator> toErase;
332      set<size_t> usedLeg1, usedLeg2;
333      if ( !doCrossCleanup ) {
334        /// @todo Reinstate when C++11 allowed: for (auto key = v.begin(); key != v.end(); ++key) {
335        for (map<double, pair<size_t, size_t> >::iterator key = v.begin(); key != v.end(); ++key) {
336          const size_t leg1 = key->second.first;
337          const size_t leg2 = key->second.second;
338          if (usedLeg1.find(leg1) == usedLeg1.end() and
339              usedLeg2.find(leg2) == usedLeg2.end()) {
340            usedLeg1.insert(leg1);
341            usedLeg2.insert(leg2);
342          } else {
343            toErase.push_back(key);
344          }
345        }
346      }
347      else {
348        /// @todo Reinstate when C++11 allowed: for (auto key = v.begin(); key != v.end(); ++key) {
349        for (map<double, pair<size_t, size_t> >::iterator key = v.begin(); key != v.end(); ++key) {
350          const size_t leg1 = key->second.first;
351          const size_t leg2 = key->second.second;
352          if (usedLeg1.find(leg1) == usedLeg1.end() and
353              usedLeg1.find(leg2) == usedLeg1.end()) {
354            usedLeg1.insert(leg1);
355            usedLeg1.insert(leg2);
356          } else {
357            toErase.push_back(key);
358          }
359        }
360      }
361      /// @todo Reinstate when C++11 allowed:  for (auto& key : toErase) v.erase(key);
362      for (size_t i = 0; i < toErase.size(); ++i) v.erase(toErase[i]);
363    }
364
365
366    void PseudoTop::project(const Event& e) {
367      // Leptons : do the lepton clustering anti-kt R=0.1 using stable photons and leptons not from hadron decay
368      // Neutrinos : neutrinos not from hadron decay
369      // MET : vector sum of all invisible particles in x-y plane
370      // Jets : anti-kt R=0.4 using all particles excluding neutrinos and particles used in lepton clustering
371      //        add ghost B hadrons during the jet clustering to identify B jets.
372
373      // W->lv : dressed lepton and neutrino pairs
374      // W->jj : light flavored dijet
375      // W candidate : select lv or jj pairs which minimise |mW1-80.4|+|mW2-80.4|
376      //               lepton-neutrino pair will be selected with higher priority
377
378      // t->Wb : W candidate + b jet
379      // t candidate : select Wb pairs which minimise |mtop1-172.5|+|mtop2-172.5|
380
381      _isValid = false;
382      _wDecays1.clear();
383      _wDecays2.clear();
384      _jets.clear();
385      _bjets.clear();
386      _ljets.clear();
387      _mode1 = _mode2 = CH_HADRON;
388
389      // Collect final state particles
390      Particles pForLep, pForJet;
391      Particles neutrinos; // Prompt neutrinos
392      /// @todo Avoid this unsafe jump into HepMC -- all this can be done properly via VisibleFS and HeavyHadrons projections
393      for (ConstGenParticlePtr p : HepMCUtils::particles(e.genEvent())) {//
394        const int status = p->status();
395        const int pid = p->pdg_id();
396        if (status == 1) {
397          Particle rp = *p;
398          if (!PID::isHadron(pid) && !rp.fromHadron()) {
399            // Collect particles not from hadron decay
400            if (rp.isNeutrino()) {
401              // Prompt neutrinos are kept in separate collection
402              neutrinos.push_back(rp);
403            } else if (pid == PID::PHOTON || rp.isLepton()) {
404              // Leptons and photons for the dressing
405              pForLep.push_back(rp);
406            }
407          } else if (!rp.isNeutrino()) {
408            // Use all particles from hadron decay
409            pForJet.push_back(rp);
410          }
411        } else if (PID::isHadron(pid) && PID::hasBottom(pid)) {
412          // NOTE: Consider B hadrons with pT > 5GeV - not in CMS proposal
413          //if ( p->momentum().perp() < 5 ) continue;
414
415          // Do unstable particles, to be used in the ghost B clustering
416          // Use last B hadrons only
417          bool isLast = true;
418          for (ConstGenParticlePtr pp : HepMCUtils::particles(p->end_vertex(), Relatives::CHILDREN)) {
419            if (PID::hasBottom(pp->pdg_id())) {
420              isLast = false;
421              break;
422            }
423          }
424          if (!isLast) continue;
425
426          // Rescale momentum by 10^-20
427          /// @todo Why the factor of 1/rho() as well?
428          //Particle ghost(pdgId, FourMomentum(p->momentum())*1e-20/p->momentum().rho());
429          Particle ghost(pid, FourMomentum(p->momentum()));
430          ghost.setMomentum(ghost.momentum()*1.e-20 / ghost.momentum().rho());
431          pForJet.push_back(ghost);
432        }
433      }
434
435      // Start object building from trivial thing - prompt neutrinos
436      isortByPt(neutrinos);
437
438      // Proceed to lepton dressing
439      FastJets fjLep(FinalState(), JetAlg::ANTIKT, _lepR);
440      fjLep.calc(pForLep);
441
442      Jets leptons;
443      vector<int> leptonsId;
444      set<int> dressedIdxs;
445      for (const Jet& lep : fjLep.jetsByPt(Cuts::pT > _lepMinPt && Cuts::abseta < _lepMaxEta)) {
446        double leadingPt = -1;
447        int leptonId = 0;
448        for (const Particle& p : lep.particles()) {
449          /// @warning Barcodes aren't future-proof in HepMC
450          dressedIdxs.insert(HepMCUtils::uniqueId(p.genParticle()));
451          if (p.isLepton() && p.pT() > leadingPt) {
452            leadingPt = p.pT();
453            leptonId = p.pid();
454          }
455        }
456        if (leptonId == 0) continue;
457        leptons.push_back(lep);
458        leptonsId.push_back(leptonId);
459      }
460
461      // Re-use particles not used in lepton dressing
462      for (const Particle& rp : pForLep) {
463        const int barcode = HepMCUtils::uniqueId(rp.genParticle());
464        // Skip if the particle is used in dressing
465        if (dressedIdxs.find(barcode) != dressedIdxs.end()) continue;
466        // Put back to be used in jet clustering
467        pForJet.push_back(rp);
468      }
469
470      // Then do the jet clustering
471      FastJets fjJet(FinalState(), JetAlg::ANTIKT, _jetR);
472      //fjJet.useInvisibles(); // NOTE: CMS proposal to remove neutrinos (AB: wouldn't work anyway, since they were excluded from clustering inputs)
473      fjJet.calc(pForJet);
474      for (const Jet& jet : fjJet.jetsByPt(Cuts::pT > _jetMinPt && Cuts::abseta < _jetMaxEta)) {
475        _jets.push_back(jet);
476        bool isBJet = false;
477        for (const Particle& rp : jet.particles()) {
478          if (PID::hasBottom(rp.pid())) {
479            isBJet = true;
480            break;
481          }
482        }
483        if ( isBJet ) _bjets.push_back(jet);
484        else _ljets.push_back(jet);
485      }
486
487      // Every building blocks are ready. Continue to pseudo-W and pseudo-top combination
488
489      if (_bjets.size() < 2) return; // Ignore single top for now
490      map<double, pair<size_t, size_t> > wLepCandIdxs;
491      map<double, pair<size_t, size_t> > wHadCandIdxs;
492
493      // Collect leptonic-decaying W's
494      for (size_t iLep = 0, nLep = leptons.size(); iLep < nLep; ++iLep) {
495        const Jet& lep = leptons.at(iLep);
496        for (size_t iNu = 0, nNu = neutrinos.size(); iNu < nNu; ++iNu) {
497          const Particle& nu = neutrinos.at(iNu);
498          const double m = (lep.momentum()+nu.momentum()).mass();
499          const double dm = std::abs(m-_wMass);
500          wLepCandIdxs[dm] = make_pair(iLep, iNu);
501        }
502      }
503
504      // Continue to hadronic decaying W's
505      for (size_t i = 0, nLjet = _ljets.size(); i < nLjet; ++i) {
506        const Jet& ljet1 = _ljets[i];
507        for (size_t j = i+1; j < nLjet; ++j) {
508          const Jet& ljet2 = _ljets[j];
509          const double m = (ljet1.momentum()+ljet2.momentum()).mass();
510          const double dm = std::abs(m-_wMass);
511          wHadCandIdxs[dm] = make_pair(i, j);
512        }
513      }
514
515      // Cleanup W candidate, choose pairs with minimum dm if they share decay products
516      cleanup(wLepCandIdxs);
517      cleanup(wHadCandIdxs, true);
518      const size_t nWLepCand = wLepCandIdxs.size();
519      const size_t nWHadCand = wHadCandIdxs.size();
520
521      if (nWLepCand + nWHadCand < 2) return; // We skip single top
522
523      int w1Q = 1, w2Q = -1;
524      int w1dau1Id = 1, w2dau1Id = -1;
525      FourMomentum w1dau1LVec, w1dau2LVec;
526      FourMomentum w2dau1LVec, w2dau2LVec;
527      if (nWLepCand == 0) { // Full hadronic case
528        const pair<size_t, size_t>& idPair1 = wHadCandIdxs.begin()->second;
529        const pair<size_t, size_t>& idPair2 = (++wHadCandIdxs.begin())->second;  ///< @todo Reinstate std::next
530        const Jet& w1dau1 = _ljets[idPair1.first];
531        const Jet& w1dau2 = _ljets[idPair1.second];
532        const Jet& w2dau1 = _ljets[idPair2.first];
533        const Jet& w2dau2 = _ljets[idPair2.second];
534        w1dau1LVec = w1dau1.momentum();
535        w1dau2LVec = w1dau2.momentum();
536        w2dau1LVec = w2dau1.momentum();
537        w2dau2LVec = w2dau2.momentum();
538      } else if (nWLepCand == 1) { // Semi-leptonic case
539        const pair<size_t, size_t>& idPair1 = wLepCandIdxs.begin()->second;
540        const pair<size_t, size_t>& idPair2 = wHadCandIdxs.begin()->second;
541        const Jet& w1dau1 = leptons[idPair1.first];
542        const Particle& w1dau2 = neutrinos[idPair1.second];
543        const Jet& w2dau1 = _ljets[idPair2.first];
544        const Jet& w2dau2 = _ljets[idPair2.second];
545        w1dau1LVec = w1dau1.momentum();
546        w1dau2LVec = w1dau2.momentum();
547        w2dau1LVec = w2dau1.momentum();
548        w2dau2LVec = w2dau2.momentum();
549        w1dau1Id = leptonsId[idPair1.first];
550        w1Q = w1dau1Id > 0 ? -1 : 1;
551        w2Q = -w1Q;
552        switch (w1dau1Id) {
553        case 13: case -13: _mode1 = CH_MUON; break;
554        case 11: case -11: _mode1 = CH_ELECTRON; break;
555        }
556      } else { // Full leptonic case
557        const pair<size_t, size_t>& idPair1 = wLepCandIdxs.begin()->second;
558        const pair<size_t, size_t>& idPair2 = (++wLepCandIdxs.begin())->second;  ///< @todo Reinstate std::next
559        const Jet& w1dau1 = leptons[idPair1.first];
560        const Particle& w1dau2 = neutrinos[idPair1.second];
561        const Jet& w2dau1 = leptons[idPair2.first];
562        const Particle& w2dau2 = neutrinos[idPair2.second];
563        w1dau1LVec = w1dau1.momentum();
564        w1dau2LVec = w1dau2.momentum();
565        w2dau1LVec = w2dau1.momentum();
566        w2dau2LVec = w2dau2.momentum();
567        w1dau1Id = leptonsId[idPair1.first];
568        w2dau1Id = leptonsId[idPair2.first];
569        w1Q = w1dau1Id > 0 ? -1 : 1;
570        w2Q = w2dau1Id > 0 ? -1 : 1;
571        switch (w1dau1Id) {
572        case 13: case -13: _mode1 = CH_MUON; break;
573        case 11: case -11: _mode1 = CH_ELECTRON; break;
574        }
575        switch (w2dau1Id) {
576        case 13: case -13: _mode2 = CH_MUON; break;
577        case 11: case -11: _mode2 = CH_ELECTRON; break;
578        }
579      }
580      const FourMomentum w1LVec = w1dau1LVec+w1dau2LVec;
581      const FourMomentum w2LVec = w2dau1LVec+w2dau2LVec;
582
583      // Combine b jets
584      double sumDm = 1e9;
585      FourMomentum b1LVec, b2LVec;
586      for (size_t i = 0, n = _bjets.size(); i < n; ++i) {
587        const Jet& bjet1 = _bjets[i];
588        const double mtop1 = (w1LVec+bjet1.momentum()).mass();
589        const double dmtop1 = std::abs(mtop1-_tMass);
590        for (size_t j=0; j<n; ++j) {
591          if (i == j) continue;
592          const Jet& bjet2 = _bjets[j];
593          const double mtop2 = (w2LVec+bjet2.momentum()).mass();
594          const double dmtop2 = std::abs(mtop2-_tMass);
595
596          if (sumDm <= dmtop1+dmtop2) continue;
597
598          sumDm = dmtop1+dmtop2;
599          b1LVec = bjet1.momentum();
600          b2LVec = bjet2.momentum();
601        }
602      }
603      if (sumDm >= 1e9) return; // Failed to make top, but this should not happen.
604
605      const FourMomentum t1LVec = w1LVec + b1LVec;
606      const FourMomentum t2LVec = w2LVec + b2LVec;
607
608      // Put all of them into candidate collection
609      _t1 = Particle(w1Q*6, t1LVec);
610      _b1 = Particle(w1Q*5, b1LVec);
611      _wDecays1.push_back(Particle(w1dau1Id, w1dau1LVec));
612      _wDecays1.push_back(Particle(-w1dau1Id+w1Q, w1dau2LVec));
613
614      _t2 = Particle(w2Q*6, t2LVec);
615      _b2 = Particle(w2Q*5, b2LVec);
616      _wDecays2.push_back(Particle(w2dau1Id, w2dau1LVec));
617      _wDecays2.push_back(Particle(-w2dau1Id+w2Q, w2dau2LVec));
618
619      _isValid = true;
620    }
621
622  }
623
624
625}