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