1#include "Rivet/Analysis.hh"
  2#include "Rivet/Projections/FinalState.hh"
  3#include "Rivet/Projections/VetoedFinalState.hh"
  4#include "Rivet/Projections/IdentifiedFinalState.hh"
  5#include "Rivet/Projections/PromptFinalState.hh"
  6#include "Rivet/Projections/LeptonFinder.hh"
  7#include "Rivet/Projections/FastJets.hh"
  8#include "Rivet/Projections/MissingMomentum.hh"
  9
 10namespace Rivet {
 11
 12
 13  /// ttbar l+jets cross sections at 13 TeV
 14  class ATLAS_2018_I1656578 : public Analysis {
 15  public:
 16
 17    /// Constructor
 18    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2018_I1656578);
 19
 20
 21    /// Book cuts and projections
 22    void init() {
 23      // Eta ranges
 24      Cut eta_full = (Cuts::abseta < 5.0);
 25      Cut lep_cuts = (Cuts::abseta < 2.5) && (Cuts::pT > 25*GeV);
 26
 27      // All final state particles
 28      FinalState fs(eta_full);
 29
 30      // Get photons to dress leptons
 31      IdentifiedFinalState all_photons(fs);
 32      all_photons.acceptIdPair(PID::PHOTON);
 33
 34      PromptFinalState photons(Cuts::abspid == PID::PHOTON, TauDecaysAs::PROMPT);
 35      declare(photons, "photons");
 36
 37      // Projection to find the electrons
 38      PromptFinalState electrons(Cuts::abspid == PID::ELECTRON, TauDecaysAs::PROMPT);
 39
 40      LeptonFinder dressedelectrons(electrons, photons, 0.1, lep_cuts);
 41      declare(dressedelectrons, "elecs");
 42
 43      LeptonFinder ewdressedelectrons(electrons, all_photons, 0.1, eta_full);
 44
 45      // Projection to find the muons
 46      PromptFinalState muons(Cuts::abspid == PID::MUON, TauDecaysAs::PROMPT);
 47
 48      LeptonFinder dressedmuons(muons, photons, 0.1, lep_cuts);
 49      declare(dressedmuons, "muons");
 50
 51      LeptonFinder ewdressedmuons(muons, all_photons, 0.1, eta_full);
 52
 53      // Projection to find MET
 54      declare(MissingMomentum(fs), "MET");
 55
 56      // Jet clustering.
 57      VetoedFinalState vfs(fs);
 58      vfs.addVetoOnThisFinalState(ewdressedelectrons);
 59      vfs.addVetoOnThisFinalState(ewdressedmuons);
 60      FastJets jets(vfs, JetAlg::ANTIKT, 0.4, JetMuons::ALL, JetInvisibles::DECAY);
 61      declare(jets, "jets");
 62
 63
 64      book(_h["absPout_inc"],                 114, 1, 1);
 65      book(_h["absPout_inc_norm"],            115, 1, 1);
 66      book(_h["ptpseudotophadron_r1"],        98, 1, 1);
 67      book(_h["ptpseudotophadron_r1_norm"],   99, 1, 1);
 68      book(_h["ptttbar_r1"],                  100, 1, 1);
 69      book(_h["ptttbar_r1_norm"],             101, 1, 1);
 70      book(_h["absPout_r1"],                  96, 1, 1);
 71      book(_h["absPout_r1_norm"],             97, 1, 1);
 72      book(_h["ptpseudotophadron_r2"],        110, 1, 1);
 73      book(_h["ptpseudotophadron_r2_norm"],   111, 1, 1);
 74      book(_h["ptttbar_r2"],                  112, 1, 1);
 75      book(_h["ptttbar_r2_norm"],             113, 1, 1);
 76      book(_h["absPout_r2"],                  108, 1, 1);
 77      book(_h["absPout_r2_norm"],             109, 1, 1);
 78      book(_h["ptpseudotophadron_r3"],        104, 1, 1);
 79      book(_h["ptpseudotophadron_r3_norm"],   105, 1, 1);
 80      book(_h["ptttbar_r3"],                  106, 1, 1);
 81      book(_h["ptttbar_r3_norm"],             107, 1, 1);
 82      book(_h["absPout_r3"],                  102, 1, 1);
 83      book(_h["absPout_r3_norm"],             103, 1, 1);
 84
 85    }
 86
 87
 88    void analyze(const Event& event) {
 89
 90      // Get the selected objects, using the projections.
 91      DressedLeptons electrons = apply<LeptonFinder>(event, "elecs").dressedLeptons();
 92      DressedLeptons muons     = apply<LeptonFinder>(event, "muons").dressedLeptons();
 93      const Jets& jets = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5);
 94
 95      const Vector3 met = apply<MissingMomentum>(event, "MET").vectorMPT();
 96
 97      Jets bjets, lightjets;
 98      for (Jet jet : jets) {
 99        bool b_tagged = jet.bTagged(Cuts::pT > 5*GeV);
100        if ( b_tagged && bjets.size() < 2) bjets += jet;
101        else lightjets += jet;
102      }
103
104      bool single_electron = (electrons.size() == 1) && (muons.empty());
105      bool single_muon = (muons.size() == 1) && (electrons.empty());
106
107
108      DressedLepton *lepton = NULL;
109      if (single_electron)   lepton = &electrons[0];
110      else if (single_muon)  lepton = &muons[0];
111
112      if (!single_electron && !single_muon)  vetoEvent;
113
114      bool num_b_tagged_jets = (bjets.size() == 2);
115      if (!num_b_tagged_jets)  vetoEvent;
116
117      if (jets.size() < 4)  vetoEvent;
118
119      bool reg_4jex2bin = (jets.size() == 4);
120      bool reg_5jex2bin = (jets.size() == 5);
121      bool reg_6jin2bin = (jets.size() >= 6);
122
123
124      FourMomentum pbjet1; //Momentum of bjet1
125      FourMomentum pbjet2; //Momentum of bjet
126
127      if ( deltaR(bjets[0], *lepton) <= deltaR(bjets[1], *lepton) ) {
128        pbjet1 = bjets[0].momentum();
129        pbjet2 = bjets[1].momentum();
130      } else {
131        pbjet1 = bjets[1].momentum();
132        pbjet2 = bjets[0].momentum();
133      }
134
135      double bestWmass = 1000.0*TeV;
136      double mWPDG = 80.399*GeV;
137      int Wj1index = -1, Wj2index = -1;
138      for (unsigned int i = 0; i < (lightjets.size() - 1); ++i) {
139        for (unsigned int j = i + 1; j < lightjets.size(); ++j) {
140          double wmass = (lightjets[i].momentum() + lightjets[j].momentum()).mass();
141          if (fabs(wmass - mWPDG) < fabs(bestWmass - mWPDG)) {
142            bestWmass = wmass;
143            Wj1index = i;
144            Wj2index = j;
145          }
146        }
147      }
148
149      FourMomentum pjet1 = lightjets[Wj1index].momentum();
150      FourMomentum pjet2 = lightjets[Wj2index].momentum();
151
152      // compute hadronic W boson
153      FourMomentum pWhadron = pjet1 + pjet2;
154      double pz = computeneutrinoz(lepton->momentum(), met);
155      FourMomentum ppseudoneutrino( sqrt(sqr(met.x()) + sqr(met.y()) + sqr(pz)), met.x(), met.y(), pz);
156
157      //compute leptonic, hadronic, combined pseudo-top
158      FourMomentum ppseudotoplepton = lepton->momentum() + ppseudoneutrino + pbjet1;
159      FourMomentum ppseudotophadron = pbjet2 + pWhadron;
160      FourMomentum pttbar = ppseudotoplepton + ppseudotophadron;
161
162      Vector3 z_versor(0,0,1);
163      Vector3 vpseudotophadron = ppseudotophadron.vector3();
164      Vector3 vpseudotoplepton = ppseudotoplepton.vector3();
165      // Variables
166      double absPout = fabs(vpseudotophadron.dot((vpseudotoplepton.cross(z_versor))/(vpseudotoplepton.cross(z_versor).mod())));
167
168      //pseudotop hadrons and leptons fill histogram
169      if (reg_4jex2bin) {
170        _h["ptpseudotophadron_r1"]->fill(ppseudotophadron.pt()); //pT of pseudo top hadron
171        _h["ptpseudotophadron_r1_norm"]->fill(ppseudotophadron.pt()); //pT of pseudo top hadron
172        _h["ptttbar_r1"]->fill(pttbar.pt()); //fill pT of ttbar in combined channel
173        _h["ptttbar_r1_norm"]->fill(pttbar.pt()); //fill pT of ttbar in combined channel
174        _h["absPout_r1"]->fill(absPout);
175        _h["absPout_r1_norm"]->fill(absPout);
176      }
177      if (reg_5jex2bin) {
178        _h["ptpseudotophadron_r2"]->fill(ppseudotophadron.pt()); //pT of pseudo top hadron
179        _h["ptpseudotophadron_r2_norm"]->fill(ppseudotophadron.pt()); //pT of pseudo top hadron
180        _h["ptttbar_r2"]->fill(pttbar.pt()); //fill pT of ttbar in combined channel
181        _h["ptttbar_r2_norm"]->fill(pttbar.pt()); //fill pT of ttbar in combined channel
182        _h["absPout_r2"]->fill(absPout);
183        _h["absPout_r2_norm"]->fill(absPout);
184      }
185      if (reg_6jin2bin) {
186        _h["ptpseudotophadron_r3"]->fill(ppseudotophadron.pt()); //pT of pseudo top hadron
187        _h["ptpseudotophadron_r3_norm"]->fill(ppseudotophadron.pt()); //pT of pseudo top hadron
188        _h["ptttbar_r3"]->fill(pttbar.pt()); //fill pT of ttbar in combined channel
189        _h["ptttbar_r3_norm"]->fill(pttbar.pt()); //fill pT of ttbar in combined channel
190        _h["absPout_r3"]->fill(absPout);
191        _h["absPout_r3_norm"]->fill(absPout);
192      }
193      _h["absPout_inc"]->fill(absPout);
194      _h["absPout_inc_norm"]->fill(absPout);
195    }
196
197
198    void finalize() {
199      // Normalize to cross-section
200      const double sf = (crossSection()/picobarn / sumOfWeights());
201      for (auto hist : _h) {
202        scale(hist.second, sf);
203        // Normalized distributions
204        if (hist.first.find("_norm") != string::npos)  normalize(hist.second);
205      }
206    }
207
208
209    double computeneutrinoz(const FourMomentum& lepton, const Vector3 &met) const {
210      // computing z component of neutrino momentum given lepton and met
211      double pzneutrino;
212      double m_W = 80.399; // in GeV, given in the paper
213      double k = (( sqr( m_W ) - sqr( lepton.mass() ) ) / 2 ) + (lepton.px() * met.x() + lepton.py() * met.y());
214      double a = sqr ( lepton.E() )- sqr ( lepton.pz() );
215      double b = -2*k*lepton.pz();
216      double c = sqr( lepton.E() ) * sqr( met.mod() ) - sqr( k );
217      double discriminant = sqr(b) - 4 * a * c;
218      double quad[2] = { (- b - sqrt(discriminant)) / (2 * a), (- b + sqrt(discriminant)) / (2 * a) }; //two possible quadratic solns
219      if (discriminant < 0)  pzneutrino = - b / (2 * a); //if the discriminant is negative
220      else { //if the discriminant is greater than or equal to zero, take the soln with smallest absolute value
221        double absquad[2];
222        for (int n=0; n<2; ++n)  absquad[n] = fabs(quad[n]);
223        if (absquad[0] < absquad[1])  pzneutrino = quad[0];
224        else                          pzneutrino = quad[1];
225      }
226      return pzneutrino;
227    }
228
229
230  private:
231
232    /// @name Objects that are used by the event selection decisions
233    map<string, Histo1DPtr> _h;
234
235  };
236
237
238  RIVET_DECLARE_PLUGIN(ATLAS_2018_I1656578);
239
240
241}