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      /// @name Standard constructors and destructors.
 18      //@{
 19
 20      /// The default constructor. May specify the minimum and maximum
 21      /// pseudorapidity \f$ \eta \f$ and the min \f$ p_T \f$ (in GeV).
 22      PseudoTop(double lepR = 0.1, double lepMinPt = 20, double lepMaxEta = 2.4,
 23                double jetR = 0.4, double jetMinPt = 30, double jetMaxEta = 4.7)
 24        : FinalState(),
 25          _lepR(lepR), _lepMinPt(lepMinPt), _lepMaxEta(lepMaxEta),
 26          _jetR(jetR), _jetMinPt(jetMinPt), _jetMaxEta(jetMaxEta)
 27      {
 28        setName("PseudoTop");
 29      }
 30
 31      enum TTbarMode {CH_NONE=-1, CH_FULLHADRON = 0, CH_SEMILEPTON, CH_FULLLEPTON};
 32      enum DecayMode {CH_HADRON = 0, CH_MUON, CH_ELECTRON};
 33
 34      TTbarMode mode() const {
 35        if (!_isValid) return CH_NONE;
 36        if (_mode1 == CH_HADRON && _mode2 == CH_HADRON) return CH_FULLHADRON;
 37        else if ( _mode1 != CH_HADRON && _mode2 != CH_HADRON) return CH_FULLLEPTON;
 38        else return CH_SEMILEPTON;
 39      }
 40      virtual const DecayMode& mode1() const {return _mode1;}
 41      virtual const DecayMode& mode2() const {return _mode2;}
 42
 43      /// Clone on the heap.
 44      DEFAULT_RIVET_PROJ_CLONE(PseudoTop);
 45
 46      //@}
 47
 48    public:
 49      virtual const Particle& t1() const {return _t1;}
 50      virtual const Particle& t2() const {return _t2;}
 51      virtual const Particle& b1() const {return _b1;}
 52      virtual const Particle& b2() const {return _b2;}
 53      virtual const Particles& wDecays1() const {return _wDecays1;}
 54      virtual const Particles& wDecays2() const {return _wDecays2;}
 55      virtual const Jets& jets() const {return _jets;}
 56      virtual const Jets& bjets() const {return _bjets;}
 57      virtual const Jets& ljets() const {return _ljets;}
 58
 59      // Apply the projection to the event
 60      void project(const Event& e); // override; ///< @todo Re-enable when C++11 allowed
 61      void cleanup(std::map<double, std::pair<size_t, size_t> >& v, const bool doCrossCleanup=false) const;
 62      CmpState compare(const Projection& p) const;
 63
 64
 65    private:
 66      const double _lepR, _lepMinPt, _lepMaxEta;
 67      const double _jetR, _jetMinPt, _jetMaxEta;
 68
 69      //constexpr ///< @todo Re-enable when C++11 allowed
 70      static double _tMass; // = 172.5*GeV; ///< @todo Re-enable when C++11 allowed
 71      //constexpr ///< @todo Re-enable when C++11 allowed
 72      static double _wMass; // = 80.4*GeV; ///< @todo Re-enable when C++11 allowed
 73
 74    private:
 75      bool _isValid;
 76      DecayMode _mode1, _mode2;
 77
 78      Particle _t1, _t2;
 79      Particle _b1, _b2;
 80      Particles _wDecays1, _wDecays2;
 81      Jets _jets, _bjets, _ljets;
 82
 83    };
 84
 85    // More implementation below the analysis code
 86
 87  }
 88
 89
 90
 91  /// Pseudo-top analysis from CMS
 92  class CMS_2015_I1370682 : public Analysis {
 93  public:
 94
 95    CMS_2015_I1370682()
 96      : Analysis("CMS_2015_I1370682"),
 97        _applyCorrection(true),
 98        _doShapeOnly(true)
 99    {    }
100
101
102    void init() {
103      declare(PseudoTop(0.1, 20, 2.4, 0.5, 30, 2.4), "ttbar");
104
105      // Lepton + Jet channel
106      book(_hSL_topPt         ,"d15-x01-y01"); // 1/sigma dsigma/dpt(top)
107      book(_hSL_topPtTtbarSys ,"d16-x01-y01"); // 1/sigma dsigma/dpt*(top)
108      book(_hSL_topY          ,"d17-x01-y01"); // 1/sigma dsigma/dy(top)
109      book(_hSL_ttbarDelPhi   ,"d18-x01-y01"); // 1/sigma dsigma/ddeltaphi(t,tbar)
110      book(_hSL_topPtLead     ,"d19-x01-y01"); // 1/sigma dsigma/dpt(t1)
111      book(_hSL_topPtSubLead  ,"d20-x01-y01"); // 1/sigma dsigma/dpt(t2)
112      book(_hSL_ttbarPt       ,"d21-x01-y01"); // 1/sigma dsigma/dpt(ttbar)
113      book(_hSL_ttbarY        ,"d22-x01-y01"); // 1/sigma dsigma/dy(ttbar)
114      book(_hSL_ttbarMass     ,"d23-x01-y01"); // 1/sigma dsigma/dm(ttbar)
115
116      // Dilepton channel
117      book(_hDL_topPt         ,"d24-x01-y01"); // 1/sigma dsigma/dpt(top)
118      book(_hDL_topPtTtbarSys ,"d25-x01-y01"); // 1/sigma dsigma/dpt*(top)
119      book(_hDL_topY          ,"d26-x01-y01"); // 1/sigma dsigma/dy(top)
120      book(_hDL_ttbarDelPhi   ,"d27-x01-y01"); // 1/sigma dsigma/ddeltaphi(t,tbar)
121      book(_hDL_topPtLead     ,"d28-x01-y01"); // 1/sigma dsigma/dpt(t1)
122      book(_hDL_topPtSubLead  ,"d29-x01-y01"); // 1/sigma dsigma/dpt(t2)
123      book(_hDL_ttbarPt       ,"d30-x01-y01"); // 1/sigma dsigma/dpt(ttbar)
124      book(_hDL_ttbarY        ,"d31-x01-y01"); // 1/sigma dsigma/dy(ttbar)
125      book(_hDL_ttbarMass     ,"d32-x01-y01"); // 1/sigma dsigma/dm(ttbar)
126
127    }
128
129
130    void analyze(const Event& event) {
131
132      // Get the ttbar candidate
133      const PseudoTop& ttbar = apply<PseudoTop>(event, "ttbar");
134      if ( ttbar.mode() == PseudoTop::CH_NONE ) vetoEvent;
135
136      const FourMomentum& t1P4 = ttbar.t1().momentum();
137      const FourMomentum& t2P4 = ttbar.t2().momentum();
138      const double pt1 = std::max(t1P4.pT(), t2P4.pT());
139      const double pt2 = std::min(t1P4.pT(), t2P4.pT());
140      const double dPhi = deltaPhi(t1P4, t2P4);
141      const FourMomentum ttP4 = t1P4 + t2P4;
142      const FourMomentum t1P4AtCM = LorentzTransform::mkFrameTransformFromBeta(ttP4.betaVec()).transform(t1P4);
143
144      const double weight = 1.0;
145
146      if ( ttbar.mode() == PseudoTop::CH_SEMILEPTON ) {
147        const Particle lCand1 = ttbar.wDecays1()[0]; // w1 dau0 is the lepton in the PseudoTop
148        if (lCand1.pT() < 33*GeV || lCand1.abseta() > 2.1) vetoEvent;
149        _hSL_topPt->fill(t1P4.pT(), weight);
150        _hSL_topPt->fill(t2P4.pT(), weight);
151        _hSL_topPtTtbarSys->fill(t1P4AtCM.pT(), weight);
152        _hSL_topY->fill(t1P4.rapidity(), weight);
153        _hSL_topY->fill(t2P4.rapidity(), weight);
154        _hSL_ttbarDelPhi->fill(dPhi, weight);
155        _hSL_topPtLead->fill(pt1, weight);
156        _hSL_topPtSubLead->fill(pt2, weight);
157        _hSL_ttbarPt->fill(ttP4.pT(), weight);
158        _hSL_ttbarY->fill(ttP4.rapidity(), weight);
159        _hSL_ttbarMass->fill(ttP4.mass(), weight);
160      }
161      else if ( ttbar.mode() == PseudoTop::CH_FULLLEPTON ) {
162        const Particle lCand1 = ttbar.wDecays1()[0]; // dau0 are the lepton in the PseudoTop
163        const Particle lCand2 = ttbar.wDecays2()[0]; // dau0 are the lepton in the PseudoTop
164        if (lCand1.pT() < 20*GeV || lCand1.abseta() > 2.4) vetoEvent;
165        if (lCand2.pT() < 20*GeV || lCand2.abseta() > 2.4) vetoEvent;
166        _hDL_topPt->fill(t1P4.pT(), weight);
167        _hDL_topPt->fill(t2P4.pT(), weight);
168        _hDL_topPtTtbarSys->fill(t1P4AtCM.pT(), weight);
169        _hDL_topY->fill(t1P4.rapidity(), weight);
170        _hDL_topY->fill(t2P4.rapidity(), weight);
171        _hDL_ttbarDelPhi->fill(dPhi, weight);
172        _hDL_topPtLead->fill(pt1, weight);
173        _hDL_topPtSubLead->fill(pt2, weight);
174        _hDL_ttbarPt->fill(ttP4.pT(), weight);
175        _hDL_ttbarY->fill(ttP4.rapidity(), weight);
176        _hDL_ttbarMass->fill(ttP4.mass(), weight);
177      }
178
179    }
180
181
182    void finalize() {
183      if ( _applyCorrection ) {
184        // Correction functions for TOP-12-028 paper, (parton bin height)/(pseudotop bin height)
185        const double ch15[] = { 5.473609, 4.941048, 4.173346, 3.391191, 2.785644, 2.371346, 2.194161, 2.197167, };
186        const double ch16[] = { 5.470905, 4.948201, 4.081982, 3.225532, 2.617519, 2.239217, 2.127878, 2.185918, };
187        const double ch17[] = { 10.003667, 4.546519, 3.828115, 3.601018, 3.522194, 3.524694, 3.600951, 3.808553, 4.531891, 9.995370, };
188        const double ch18[] = { 4.406683, 4.054041, 3.885393, 4.213646, };
189        const double ch19[] = { 6.182537, 5.257703, 4.422280, 3.568402, 2.889408, 2.415878, 2.189974, 2.173210, };
190        const double ch20[] = { 5.199874, 4.693318, 3.902882, 3.143785, 2.607877, 2.280189, 2.204124, 2.260829, };
191        const double ch21[] = { 6.053523, 3.777506, 3.562251, 3.601356, 3.569347, 3.410472, };
192        const double ch22[] = { 11.932351, 4.803773, 3.782709, 3.390775, 3.226806, 3.218982, 3.382678, 3.773653, 4.788191, 11.905338, };
193        const double ch23[] = { 7.145255, 5.637595, 4.049882, 3.025917, 2.326430, 1.773824, 1.235329, };
194
195        const double ch24[] = { 2.268193, 2.372063, 2.323975, 2.034655, 1.736793, };
196        const double ch25[] = { 2.231852, 2.383086, 2.341894, 2.031318, 1.729672, 1.486993, };
197        const double ch26[] = { 3.993526, 2.308249, 2.075136, 2.038297, 2.036302, 2.078270, 2.295817, 4.017713, };
198        const double ch27[] = { 2.205978, 2.175010, 2.215376, 2.473144, };
199        const double ch28[] = { 2.321077, 2.371895, 2.338871, 2.057821, 1.755382, };
200        const double ch29[] = { 2.222707, 2.372591, 2.301688, 1.991162, 1.695343, };
201        const double ch30[] = { 2.599677, 2.026855, 2.138620, 2.229553, };
202        const double ch31[] = { 5.791779, 2.636219, 2.103642, 1.967198, 1.962168, 2.096514, 2.641189, 5.780828, };
203        const double ch32[] = { 2.006685, 2.545525, 2.477745, 2.335747, 2.194226, 2.076500, };
204
205        applyCorrection(_hSL_topPt, ch15);
206        applyCorrection(_hSL_topPtTtbarSys, ch16);
207        applyCorrection(_hSL_topY, ch17);
208        applyCorrection(_hSL_ttbarDelPhi, ch18);
209        applyCorrection(_hSL_topPtLead, ch19);
210        applyCorrection(_hSL_topPtSubLead, ch20);
211        applyCorrection(_hSL_ttbarPt, ch21);
212        applyCorrection(_hSL_ttbarY, ch22);
213        applyCorrection(_hSL_ttbarMass, ch23);
214
215        applyCorrection(_hDL_topPt, ch24);
216        applyCorrection(_hDL_topPtTtbarSys, ch25);
217        applyCorrection(_hDL_topY, ch26);
218        applyCorrection(_hDL_ttbarDelPhi, ch27);
219        applyCorrection(_hDL_topPtLead, ch28);
220        applyCorrection(_hDL_topPtSubLead, ch29);
221        applyCorrection(_hDL_ttbarPt, ch30);
222        applyCorrection(_hDL_ttbarY, ch31);
223        applyCorrection(_hDL_ttbarMass, ch32);
224      }
225
226      if ( _doShapeOnly ) {
227        normalize(_hSL_topPt        );
228        normalize(_hSL_topPtTtbarSys);
229        normalize(_hSL_topY         );
230        normalize(_hSL_ttbarDelPhi  );
231        normalize(_hSL_topPtLead    );
232        normalize(_hSL_topPtSubLead );
233        normalize(_hSL_ttbarPt      );
234        normalize(_hSL_ttbarY       );
235        normalize(_hSL_ttbarMass    );
236
237        normalize(_hDL_topPt        );
238        normalize(_hDL_topPtTtbarSys);
239        normalize(_hDL_topY         );
240        normalize(_hDL_ttbarDelPhi  );
241        normalize(_hDL_topPtLead    );
242        normalize(_hDL_topPtSubLead );
243        normalize(_hDL_ttbarPt      );
244        normalize(_hDL_ttbarY       );
245        normalize(_hDL_ttbarMass    );
246      }
247      else {
248        const double s = 1./sumOfWeights();
249        scale(_hSL_topPt        , s);
250        scale(_hSL_topPtTtbarSys, s);
251        scale(_hSL_topY         , s);
252        scale(_hSL_ttbarDelPhi  , s);
253        scale(_hSL_topPtLead    , s);
254        scale(_hSL_topPtSubLead , s);
255        scale(_hSL_ttbarPt      , s);
256        scale(_hSL_ttbarY       , s);
257        scale(_hSL_ttbarMass    , s);
258        scale(_hDL_topPt        , s);
259        scale(_hDL_topPtTtbarSys, s);
260        scale(_hDL_topY         , s);
261        scale(_hDL_ttbarDelPhi  , s);
262        scale(_hDL_topPtLead    , s);
263        scale(_hDL_topPtSubLead , s);
264        scale(_hDL_ttbarPt      , s);
265        scale(_hDL_ttbarY       , s);
266        scale(_hDL_ttbarMass    , s);
267      }
268
269    }
270
271
272    void applyCorrection(Histo1DPtr h, const double* cf) {
273      vector<YODA::HistoBin1D>& bins = h->bins();
274      for (size_t i=0, n=bins.size(); i<n; ++i ) {
275        const double s = cf[i];
276        YODA::HistoBin1D& bin = bins[i];
277        bin.scaleW(s);
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      sortByPt(neutrinos);
437
438      // Proceed to lepton dressing
439      FastJets fjLep(FinalState(), FastJets::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(_lepMinPt)) {
446        if (lep.abseta() > _lepMaxEta) continue;
447        double leadingPt = -1;
448        int leptonId = 0;
449        for (const Particle& p : lep.particles()) {
450          /// @warning Barcodes aren't future-proof in HepMC
451          dressedIdxs.insert(HepMCUtils::uniqueId(p.genParticle()));
452          if (p.isLepton() && p.pT() > leadingPt) {
453            leadingPt = p.pT();
454            leptonId = p.pid();
455          }
456        }
457        if (leptonId == 0) continue;
458        leptons.push_back(lep);
459        leptonsId.push_back(leptonId);
460      }
461
462      // Re-use particles not used in lepton dressing
463      for (const Particle& rp : pForLep) {
464        const int barcode = HepMCUtils::uniqueId(rp.genParticle());
465        // Skip if the particle is used in dressing
466        if (dressedIdxs.find(barcode) != dressedIdxs.end()) continue;
467        // Put back to be used in jet clustering
468        pForJet.push_back(rp);
469      }
470
471      // Then do the jet clustering
472      FastJets fjJet(FinalState(), FastJets::ANTIKT, _jetR);
473      //fjJet.useInvisibles(); // NOTE: CMS proposal to remove neutrinos (AB: wouldn't work anyway, since they were excluded from clustering inputs)
474      fjJet.calc(pForJet);
475      for (const Jet& jet : fjJet.jetsByPt(_jetMinPt)) {
476        if (jet.abseta() > _jetMaxEta) continue;
477        _jets.push_back(jet);
478        bool isBJet = false;
479        for (const Particle& rp : jet.particles()) {
480          if (PID::hasBottom(rp.pid())) {
481            isBJet = true;
482            break;
483          }
484        }
485        if ( isBJet ) _bjets.push_back(jet);
486        else _ljets.push_back(jet);
487      }
488
489      // Every building blocks are ready. Continue to pseudo-W and pseudo-top combination
490
491      if (_bjets.size() < 2) return; // Ignore single top for now
492      map<double, pair<size_t, size_t> > wLepCandIdxs;
493      map<double, pair<size_t, size_t> > wHadCandIdxs;
494
495      // Collect leptonic-decaying W's
496      for (size_t iLep = 0, nLep = leptons.size(); iLep < nLep; ++iLep) {
497        const Jet& lep = leptons.at(iLep);
498        for (size_t iNu = 0, nNu = neutrinos.size(); iNu < nNu; ++iNu) {
499          const Particle& nu = neutrinos.at(iNu);
500          const double m = (lep.momentum()+nu.momentum()).mass();
501          const double dm = std::abs(m-_wMass);
502          wLepCandIdxs[dm] = make_pair(iLep, iNu);
503        }
504      }
505
506      // Continue to hadronic decaying W's
507      for (size_t i = 0, nLjet = _ljets.size(); i < nLjet; ++i) {
508        const Jet& ljet1 = _ljets[i];
509        for (size_t j = i+1; j < nLjet; ++j) {
510          const Jet& ljet2 = _ljets[j];
511          const double m = (ljet1.momentum()+ljet2.momentum()).mass();
512          const double dm = std::abs(m-_wMass);
513          wHadCandIdxs[dm] = make_pair(i, j);
514        }
515      }
516
517      // Cleanup W candidate, choose pairs with minimum dm if they share decay products
518      cleanup(wLepCandIdxs);
519      cleanup(wHadCandIdxs, true);
520      const size_t nWLepCand = wLepCandIdxs.size();
521      const size_t nWHadCand = wHadCandIdxs.size();
522
523      if (nWLepCand + nWHadCand < 2) return; // We skip single top
524
525      int w1Q = 1, w2Q = -1;
526      int w1dau1Id = 1, w2dau1Id = -1;
527      FourMomentum w1dau1LVec, w1dau2LVec;
528      FourMomentum w2dau1LVec, w2dau2LVec;
529      if (nWLepCand == 0) { // Full hadronic case
530        const pair<size_t, size_t>& idPair1 = wHadCandIdxs.begin()->second;
531        const pair<size_t, size_t>& idPair2 = (++wHadCandIdxs.begin())->second;  ///< @todo Reinstate std::next
532        const Jet& w1dau1 = _ljets[idPair1.first];
533        const Jet& w1dau2 = _ljets[idPair1.second];
534        const Jet& w2dau1 = _ljets[idPair2.first];
535        const Jet& w2dau2 = _ljets[idPair2.second];
536        w1dau1LVec = w1dau1.momentum();
537        w1dau2LVec = w1dau2.momentum();
538        w2dau1LVec = w2dau1.momentum();
539        w2dau2LVec = w2dau2.momentum();
540      } else if (nWLepCand == 1) { // Semi-leptonic case
541        const pair<size_t, size_t>& idPair1 = wLepCandIdxs.begin()->second;
542        const pair<size_t, size_t>& idPair2 = wHadCandIdxs.begin()->second;
543        const Jet& w1dau1 = leptons[idPair1.first];
544        const Particle& w1dau2 = neutrinos[idPair1.second];
545        const Jet& w2dau1 = _ljets[idPair2.first];
546        const Jet& w2dau2 = _ljets[idPair2.second];
547        w1dau1LVec = w1dau1.momentum();
548        w1dau2LVec = w1dau2.momentum();
549        w2dau1LVec = w2dau1.momentum();
550        w2dau2LVec = w2dau2.momentum();
551        w1dau1Id = leptonsId[idPair1.first];
552        w1Q = w1dau1Id > 0 ? -1 : 1;
553        w2Q = -w1Q;
554        switch (w1dau1Id) {
555        case 13: case -13: _mode1 = CH_MUON; break;
556        case 11: case -11: _mode1 = CH_ELECTRON; break;
557        }
558      } else { // Full leptonic case
559        const pair<size_t, size_t>& idPair1 = wLepCandIdxs.begin()->second;
560        const pair<size_t, size_t>& idPair2 = (++wLepCandIdxs.begin())->second;  ///< @todo Reinstate std::next
561        const Jet& w1dau1 = leptons[idPair1.first];
562        const Particle& w1dau2 = neutrinos[idPair1.second];
563        const Jet& w2dau1 = leptons[idPair2.first];
564        const Particle& w2dau2 = neutrinos[idPair2.second];
565        w1dau1LVec = w1dau1.momentum();
566        w1dau2LVec = w1dau2.momentum();
567        w2dau1LVec = w2dau1.momentum();
568        w2dau2LVec = w2dau2.momentum();
569        w1dau1Id = leptonsId[idPair1.first];
570        w2dau1Id = leptonsId[idPair2.first];
571        w1Q = w1dau1Id > 0 ? -1 : 1;
572        w2Q = w2dau1Id > 0 ? -1 : 1;
573        switch (w1dau1Id) {
574        case 13: case -13: _mode1 = CH_MUON; break;
575        case 11: case -11: _mode1 = CH_ELECTRON; break;
576        }
577        switch (w2dau1Id) {
578        case 13: case -13: _mode2 = CH_MUON; break;
579        case 11: case -11: _mode2 = CH_ELECTRON; break;
580        }
581      }
582      const FourMomentum w1LVec = w1dau1LVec+w1dau2LVec;
583      const FourMomentum w2LVec = w2dau1LVec+w2dau2LVec;
584
585      // Combine b jets
586      double sumDm = 1e9;
587      FourMomentum b1LVec, b2LVec;
588      for (size_t i = 0, n = _bjets.size(); i < n; ++i) {
589        const Jet& bjet1 = _bjets[i];
590        const double mtop1 = (w1LVec+bjet1.momentum()).mass();
591        const double dmtop1 = std::abs(mtop1-_tMass);
592        for (size_t j=0; j<n; ++j) {
593          if (i == j) continue;
594          const Jet& bjet2 = _bjets[j];
595          const double mtop2 = (w2LVec+bjet2.momentum()).mass();
596          const double dmtop2 = std::abs(mtop2-_tMass);
597
598          if (sumDm <= dmtop1+dmtop2) continue;
599
600          sumDm = dmtop1+dmtop2;
601          b1LVec = bjet1.momentum();
602          b2LVec = bjet2.momentum();
603        }
604      }
605      if (sumDm >= 1e9) return; // Failed to make top, but this should not happen.
606
607      const FourMomentum t1LVec = w1LVec + b1LVec;
608      const FourMomentum t2LVec = w2LVec + b2LVec;
609
610      // Put all of them into candidate collection
611      _t1 = Particle(w1Q*6, t1LVec);
612      _b1 = Particle(w1Q*5, b1LVec);
613      _wDecays1.push_back(Particle(w1dau1Id, w1dau1LVec));
614      _wDecays1.push_back(Particle(-w1dau1Id+w1Q, w1dau2LVec));
615
616      _t2 = Particle(w2Q*6, t2LVec);
617      _b2 = Particle(w2Q*5, b2LVec);
618      _wDecays2.push_back(Particle(w2dau1Id, w2dau1LVec));
619      _wDecays2.push_back(Particle(-w2dau1Id+w2Q, w2dau2LVec));
620
621      _isValid = true;
622    }
623
624  }
625
626
627}