1#include "Rivet/Analysis.hh"
  2#include "Rivet/Projections/FinalState.hh"
  3#include "Rivet/Projections/VetoedFinalState.hh"
  4#include "Rivet/Projections/ChargedLeptons.hh"
  5#include "Rivet/Projections/MissingMomentum.hh"
  6#include "Rivet/Projections/FastJets.hh"
  7#include "Rivet/AnalysisLoader.hh"
  8
  9namespace Rivet {
 10
 11  
 12
 13
 14  class MC_TTBAR : public Analysis {
 15  public:
 16
 17    /// Minimal constructor
 18    RIVET_DEFAULT_ANALYSIS_CTOR(MC_TTBAR);
 19
 20
 21    /// @name Analysis methods
 22    //@{
 23
 24    /// Set up projections and book histograms
 25    void init() {
 26
 27      _mode = 1; string pre = "onelep_"; // default is single-lepton decay mode
 28      if ( getOption("TTMODE") == "ALLHAD" ) { _mode = 0; pre = "allhad_"; }
 29      if ( getOption("TTMODE") == "ONELEP" ) { _mode = 1; pre = "onelep_"; }
 30      if ( getOption("TTMODE") == "TWOLEP" ) { _mode = 2; pre = "twolep_"; }
 31      if ( getOption("TTMODE") == "ANYLEP" ) { _mode = 3; pre = "anylep_"; }
 32
 33      // A FinalState is used to select particles within |eta| < 4.2 and with pT
 34      // > 30 GeV, out of which the ChargedLeptons projection picks only the
 35      // electrons and muons, to be accessed later as "LFS".
 36      ChargedLeptons lfs(FinalState(Cuts::abseta < 4.2 && Cuts::pT > 30*GeV));
 37      declare(lfs, "LFS");
 38
 39      // A second FinalState is used to select all particles in |eta| < 4.2,
 40      // with no pT cut. This is used to construct jets and measure missing
 41      // transverse energy.
 42      VetoedFinalState fs(FinalState(Cuts::abseta < 4.2));
 43      fs.addVetoOnThisFinalState(lfs);
 44      declare(FastJets(fs, FastJets::ANTIKT, 0.6), "Jets");
 45      declare(MissingMomentum(fs), "MissingET");
 46
 47      // Booking of histograms
 48      book(_h["njets"], pre + "jet_mult", 11, -0.5, 10.5);
 49      //
 50      book(_h["jet_1_pT"], pre + "jet_1_pT", logspace(50, 20.0, 500.0));
 51      book(_h["jet_2_pT"], pre + "jet_2_pT", logspace(50, 20.0, 400.0));
 52      book(_h["jet_3_pT"], pre + "jet_3_pT", logspace(50, 20.0, 300.0));
 53      book(_h["jet_4_pT"], pre + "jet_4_pT", logspace(50, 20.0, 200.0));
 54      book(_h["jet_HT"],   pre + "jet_HT", logspace(50, 100.0, 2000.0));
 55      //
 56      book(_h["bjet_1_pT"], pre + "jetb_1_pT", logspace(50, 20.0, 400.0));
 57      book(_h["bjet_2_pT"], pre + "jetb_2_pT", logspace(50, 20.0, 300.0));
 58      //
 59      book(_h["ljet_1_pT"], pre + "jetl_1_pT", logspace(50, 20.0, 400.0));
 60      book(_h["ljet_2_pT"], pre + "jetl_2_pT", logspace(50, 20.0, 300.0));
 61      //
 62      if (_mode != 2)  book(_h["tt_mass"], pre + "tt_mass", 200, 300.0, 700.0);
 63      //
 64      if (_mode < 2) { // these rely on a hadronic W being part of the ttbar decay
 65        book(_h["W_mass"], pre + "W_mass", 75, 30, 180);
 66        book(_h["t_mass"], pre + "t_mass", 150, 130, 430);
 67        book(_h["t_mass_W_cut"], pre + "t_mass_W_cut", 150, 130, 430);
 68        book(_h["jetb_1_W_dR"],  pre + "jetb_1_W_dR", 20, 0.0, 7.0);
 69        book(_h["jetb_1_W_deta"], pre + "jetb_1_W_deta", 20, 0.0, 7.0);
 70        book(_h["jetb_1_W_dphi"], pre + "jetb_1_W_dphi", 20, 0.0, M_PI);
 71      }
 72      //
 73      book(_h["jetb_1_jetb_2_dR"],   pre + "jetb_1_jetb_2_dR", 20, 0.0, 7.0);
 74      book(_h["jetb_1_jetb_2_deta"], pre + "jetb_1_jetb_2_deta", 20, 0.0, 7.0);
 75      book(_h["jetb_1_jetb_2_dphi"], pre + "jetb_1_jetb_2_dphi", 20, 0.0, M_PI);
 76      book(_h["jetb_1_jetl_1_dR"],   pre + "jetb_1_jetl_1_dR", 20, 0.0, 7.0);
 77      book(_h["jetb_1_jetl_1_deta"], pre + "jetb_1_jetl_1_deta", 20, 0.0, 7.0);
 78      book(_h["jetb_1_jetl_1_dphi"], pre + "jetb_1_jetl_1_dphi", 20, 0.0, M_PI);
 79      book(_h["jetl_1_jetl_2_dR"],   pre + "jetl_1_jetl_2_dR", 20, 0.0, 7.0);
 80      book(_h["jetl_1_jetl_2_deta"], pre + "jetl_1_jetl_2_deta", 20, 0.0, 7.0);
 81      book(_h["jetl_1_jetl_2_dphi"], pre + "jetl_1_jetl_2_dphi", 20, 0.0, M_PI);
 82      if (_mode > 0) { // these rely on at least one leptonic decay mode
 83        book(_h["jetb_1_l_dR"],   pre + "jetb_1_l_dR", 20, 0.0, 7.0);
 84        book(_h["jetb_1_l_deta"], pre + "jetb_1_l_deta", 20, 0.0, 7.0);
 85        book(_h["jetb_1_l_dphi"], pre + "jetb_1_l_dphi", 20, 0.0, M_PI);
 86        book(_h["jetb_1_l_mass"], pre + "jetb_1_l_mass", 40, 0.0, 500.0);
 87        if (_mode > 1) {
 88          book(_h["jetb_1_l2_dR"],   pre + "jetb_1_l2_dR", 20, 0.0, 7.0);
 89          book(_h["jetb_1_l2_deta"], pre + "jetb_1_l2_deta", 20, 0.0, 7.0);
 90          book(_h["jetb_1_l2_dphi"], pre + "jetb_1_l2_dphi", 20, 0.0, M_PI);
 91          book(_h["jetb_1_l2_mass"], pre + "jetb_1_l2_mass", 40, 0.0, 500.0);
 92        }
 93      }
 94    }
 95
 96
 97    void analyze(const Event& event) {
 98      const double weight = 1.0;
 99
100      // Use the "LFS" projection to require at least one hard charged
101      // lepton. This is an experimental signature for the leptonically decaying
102      // W. This helps to reduce pure QCD backgrounds.
103      const ChargedLeptons& lfs = apply<ChargedLeptons>(event, "LFS");
104      MSG_DEBUG("Charged lepton multiplicity = " << lfs.chargedLeptons().size());
105      for (const Particle& lepton : lfs.chargedLeptons()) {
106        MSG_DEBUG("Lepton pT = " << lepton.pT());
107      }
108
109      size_t nLeps = lfs.chargedLeptons().size();
110      bool leptonMultiFail = _mode == 3 && nLeps == 0; // non-all-hadronic
111      leptonMultiFail |= _mode == 2 && nLeps != 2; // dilepton
112      leptonMultiFail |= _mode == 1 && nLeps != 1; // single lepton
113      leptonMultiFail |= _mode == 0 && nLeps != 0; // all-hadronic
114      if (leptonMultiFail) {
115        MSG_DEBUG("Event failed lepton multiplicity cut");
116        vetoEvent;
117      }
118
119      // Use a missing ET cut to bias toward events with a hard neutrino from
120      // the leptonically decaying W. This helps to reduce pure QCD backgrounds.
121      // not applied in all-hadronic mode
122      const Vector3& met = apply<MissingMomentum>(event, "MissingET").vectorMissingPt();
123      MSG_DEBUG("Vector pT = " << met.mod() << " GeV");
124      if (_mode > 0 && met.mod() < 30*GeV) {
125        MSG_DEBUG("Event failed missing ET cut");
126        vetoEvent;
127      }
128
129      // Use the "Jets" projection to check how many jets with pT > 30 GeV there are
130      // remove jets overlapping with any lepton (dR < 0.3)
131      // cut on jet multiplicity depending on ttbar decay mode
132      const FastJets& jetpro = apply<FastJets>(event, "Jets");
133      const Jets jets = discardIfAnyDeltaRLess(jetpro.jetsByPt(30*GeV), lfs.chargedLeptons(), 0.3);
134
135      if (     _mode == 0 && jets.size() < 6)  vetoEvent; // all-hadronic
136      else if (_mode == 1 && jets.size() < 4)  vetoEvent; // single lepton
137      else if (_mode == 2 && jets.size() < 2)  vetoEvent; // dilepton
138      else if (_mode == 3 && nLeps == 1 && jets.size() < 4)  vetoEvent; // non-allhadronic
139      else if (_mode == 3 && nLeps == 2 && jets.size() < 2)  vetoEvent;
140      MSG_DEBUG("Event failed jet multiplicity cut");
141
142      // Fill histograms for inclusive jet kinematics 
143      _h["njets"]->fill(jets.size(), weight);
144      if (jets.size() > 0)  _h["jet_1_pT"]->fill(jets[0].pT()/GeV, weight);
145      if (jets.size() > 1)  _h["jet_2_pT"]->fill(jets[1].pT()/GeV, weight);
146      if (jets.size() > 2)  _h["jet_3_pT"]->fill(jets[2].pT()/GeV, weight);
147      if (jets.size() > 3)  _h["jet_4_pT"]->fill(jets[3].pT()/GeV, weight);
148      double ht = 0.0;
149      for (const Jet& j : jets) { ht += j.pT(); }
150      _h["jet_HT"]->fill(ht/GeV, weight);
151
152      // Sort the jets into b-jets and light jets. We expect one hard b-jet from
153      // each top decay, so our 4 hardest jets should include two b-jets. The
154      // Jet::bTagged() method is equivalent to perfect experimental
155      // b-tagging, in a generator-independent way.
156      Jets bjets, ljets;
157      for (const Jet& jet : jets) {
158        if (jet.bTagged())  bjets += jet;
159        else                ljets += jet;
160      }
161      MSG_DEBUG("Number of b-jets = " << bjets.size());
162      MSG_DEBUG("Number of l-jets = " << ljets.size());
163      if (bjets.size() != 2) {
164        MSG_DEBUG("Event failed post-lepton-isolation b-tagging cut");
165        vetoEvent;
166      }
167      if (_mode == 0 && ljets.size() < 4)  vetoEvent;
168      else if (_mode == 1 && ljets.size() < 2)  vetoEvent;
169      else if (_mode == 3 && nLeps == 1 && ljets.size() < 2)  vetoEvent;
170
171      // Plot the pTs of the identified jets.
172      _h["bjet_1_pT"]->fill(bjets[0].pT(), weight);
173      _h["bjet_2_pT"]->fill(bjets[1].pT(), weight);
174      // need to check size to cater for dileptonic mode
175      if (ljets.size() > 0)  _h["ljet_1_pT"]->fill(ljets[0].pT(), weight);
176      if (ljets.size() > 1)  _h["ljet_2_pT"]->fill(ljets[1].pT(), weight);
177
178
179      // Try to reconstruct ttbar pair (doesn't really work in the dileptonic mode)
180      FourMomentum ttpair = bjets[0].mom() + bjets[1].mom();
181      if (_mode == 0) {
182        ttpair += ljets[0].mom() + ljets[1].mom() + ljets[2].mom() + ljets[3].mom();
183      }
184      else if (nLeps < 2) {
185        ttpair += ljets[0].mom() + ljets[1].mom();
186        const FourMomentum lep = lfs.chargedLeptons()[0].mom();
187        double pz = findZcomponent(lep, met);
188        FourMomentum neutrino(sqrt(sqr(met.x()) + sqr(met.y()) + sqr(pz)), met.x(), met.y(), pz);
189        ttpair += lep + neutrino;
190      }
191      if (nLeps < 2)  _h["tt_mass"]->fill(ttpair.mass()/GeV, weight);
192
193      if (_mode < 2) {
194        // Construct the hadronically decaying W momentum 4-vector from pairs of
195        // non-b-tagged jets. The pair which best matches the W mass is used. We start
196        // with an always terrible 4-vector estimate which should always be "beaten" by
197        // a real jet pair.
198        FourMomentum W(10*(sqrtS()>0.?sqrtS():14000.), 0, 0, 0);
199        for (size_t i = 0; i < ljets.size()-1; ++i) {
200          for (size_t j = i + 1; j < ljets.size(); ++j) {
201            const FourMomentum Wcand = ljets[i].momentum() + ljets[j].momentum();
202            MSG_TRACE(i << "," << j << ": candidate W mass = " << Wcand.mass()/GeV
203                      << " GeV, vs. incumbent candidate with " << W.mass()/GeV << " GeV");
204            if (fabs(Wcand.mass() - 80.4*GeV) < fabs(W.mass() - 80.4*GeV)) {
205              W = Wcand;
206            }
207          }
208        }
209        MSG_DEBUG("Candidate W mass = " << W.mass() << " GeV");
210
211        // There are two b-jets with which this can be combined to make the
212        // hadronically decaying top, one of which is correct and the other is
213        // not... but we have no way to identify which is which, so we construct
214        // both possible top momenta and fill the histograms with both.
215        const FourMomentum t1 = W + bjets[0].momentum();
216        const FourMomentum t2 = W + bjets[1].momentum();
217        _h["W_mass"]->fill(W.mass(), weight);
218        _h["t_mass"]->fill(t1.mass(), weight);
219        _h["t_mass"]->fill(t2.mass(), weight);
220
221        // Placing a cut on the well-known W mass helps to reduce backgrounds
222        // only done for all-hadronic and semileptonic mode (since W is hadronic)
223        if (!inRange(W.mass()/GeV, 75.0, 85.0))  vetoEvent;
224        MSG_DEBUG("W found with mass " << W.mass()/GeV << " GeV");
225
226        _h["t_mass_W_cut"]->fill(t1.mass(), weight);
227        _h["t_mass_W_cut"]->fill(t2.mass(), weight);
228
229        _h["jetb_1_W_dR"]->fill(deltaR(bjets[0].momentum(), W),weight);
230        _h["jetb_1_W_deta"]->fill(fabs(bjets[0].eta()-W.eta()),weight);
231        _h["jetb_1_W_dphi"]->fill(deltaPhi(bjets[0].momentum(),W),weight);
232      }
233
234      _h["jetb_1_jetb_2_dR"]->fill(deltaR(bjets[0].momentum(), bjets[1].momentum()),weight);
235      _h["jetb_1_jetb_2_deta"]->fill(fabs(bjets[0].eta()-bjets[1].eta()),weight);
236      _h["jetb_1_jetb_2_dphi"]->fill(deltaPhi(bjets[0].momentum(),bjets[1].momentum()),weight);
237
238      if (ljets.size() > 0) {
239        _h["jetb_1_jetl_1_dR"]->fill(deltaR(bjets[0].momentum(), ljets[0].momentum()),weight);
240        _h["jetb_1_jetl_1_deta"]->fill(fabs(bjets[0].eta()-ljets[0].eta()),weight);
241        _h["jetb_1_jetl_1_dphi"]->fill(deltaPhi(bjets[0].momentum(),ljets[0].momentum()),weight);
242        if (ljets.size() > 1) {
243          _h["jetl_1_jetl_2_dR"]->fill(deltaR(ljets[0].momentum(), ljets[1].momentum()),weight);
244          _h["jetl_1_jetl_2_deta"]->fill(fabs(ljets[0].eta()-ljets[1].eta()),weight);
245          _h["jetl_1_jetl_2_dphi"]->fill(deltaPhi(ljets[0].momentum(),ljets[1].momentum()),weight);
246        }
247      }
248
249      // lepton-centric plots
250      if (_mode > 0) {
251        FourMomentum l=lfs.chargedLeptons()[0].momentum();
252        _h["jetb_1_l_dR"]->fill(deltaR(bjets[0].momentum(), l),weight);
253        _h["jetb_1_l_deta"]->fill(fabs(bjets[0].eta()-l.eta()),weight);
254        _h["jetb_1_l_dphi"]->fill(deltaPhi(bjets[0].momentum(),l),weight);
255        _h["jetb_1_l_mass"]->fill(FourMomentum(bjets[0].momentum()+l).mass(), weight);
256
257        if (nLeps > 1) {
258          FourMomentum l=lfs.chargedLeptons()[1].momentum();
259          _h["jetb_1_l2_dR"]->fill(deltaR(bjets[0].momentum(), l),weight);
260          _h["jetb_1_l2_deta"]->fill(fabs(bjets[0].eta()-l.eta()),weight);
261          _h["jetb_1_l2_dphi"]->fill(deltaPhi(bjets[0].momentum(),l),weight);
262          _h["jetb_1_l2_mass"]->fill(FourMomentum(bjets[0].momentum()+l).mass(), weight);
263        }
264      }
265
266    }
267
268    double findZcomponent(const FourMomentum& lepton, const Vector3& met) const {
269      // estimate z-component of momentum given lepton 4-vector and MET 3-vector
270      double pz_estimate;
271      double m_W = 80.399*GeV;
272      double k = (( sqr( m_W ) - sqr( lepton.mass() ) ) / 2 ) + (lepton.px() * met.x() + lepton.py() * met.y());
273      double a = sqr ( lepton.E() )- sqr ( lepton.pz() );
274      double b = -2*k*lepton.pz();
275      double c = sqr( lepton.E() ) * sqr( met.perp() ) - sqr( k );
276      double discriminant = sqr(b) - 4 * a * c;
277      double quad[2] = { (- b - sqrt(discriminant)) / (2 * a), (- b + sqrt(discriminant)) / (2 * a) }; //two possible quadratic solns
278      if (discriminant < 0)  pz_estimate = - b / (2 * a); //if the discriminant is negative
279      else { //if the discriminant is greater than or equal to zero, take the soln with smallest absolute value
280        double absquad[2];
281        for (int n=0; n<2; ++n)  absquad[n] = fabs(quad[n]);
282        if (absquad[0] < absquad[1])  pz_estimate = quad[0];
283        else                          pz_estimate = quad[1];
284      }
285      return pz_estimate;
286    }
287
288    void finalize() {
289      const double sf = crossSection() / sumOfWeights();
290      for (auto hist : _h) { scale(hist.second, sf); }
291    }
292
293    //@}
294
295  protected:
296
297      size_t _mode;
298
299
300  private:
301
302    // @name Histogram data members
303    //@{
304    map<string, Histo1DPtr> _h;
305    //@}
306
307  };
308
309
310
311  // The hook for the plugin system
312  RIVET_DECLARE_PLUGIN(MC_TTBAR);
313
314}