1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/IdentifiedFinalState.hh"
  4#include "Rivet/Projections/VetoedFinalState.hh"
  5#include "Rivet/Projections/MissingMomentum.hh"
  6#include "Rivet/Projections/FastJets.hh"
  7#include "Rivet/Projections/DressedLeptons.hh"
  8#include "Rivet/Projections/LeadingParticlesFinalState.hh"
  9
 10namespace Rivet {
 11
 12
 13  /// ATLAS W + jets production at 7 TeV
 14  class ATLAS_2012_I1083318 : public Analysis {
 15  public:
 16
 17    /// Constructor
 18    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2012_I1083318);
 19
 20
 21    /// @name Analysis methods
 22    //@{
 23
 24    /// Book histograms and initialise projections before the run
 25    void init() {
 26
 27      FinalState fs;
 28      IdentifiedFinalState allleptons;
 29      allleptons.acceptIdPair(PID::ELECTRON);
 30      allleptons.acceptIdPair(PID::MUON);
 31      Cut cuts = Cuts::abseta < 2.5 && Cuts::pT > 20*GeV;
 32      DressedLeptons leptons(fs, allleptons, 0.1, cuts);
 33      declare(leptons, "leptons");
 34
 35      // Leading neutrinos for Etmiss
 36      LeadingParticlesFinalState neutrinos(fs);
 37      neutrinos.addParticleIdPair(PID::NU_E);
 38      neutrinos.addParticleIdPair(PID::NU_MU);
 39      neutrinos.setLeadingOnly(true);
 40      declare(neutrinos, "neutrinos");
 41
 42      // Input for the jets: "Neutrinos, electrons, and muons from decays of the
 43      // massive W boson were not used"
 44      VetoedFinalState veto;
 45      veto.addVetoOnThisFinalState(leptons);
 46      veto.addVetoOnThisFinalState(neutrinos);
 47      FastJets jets(veto, FastJets::ANTIKT, 0.4, JetAlg::Muons::ALL, JetAlg::Invisibles::DECAY);
 48      declare(jets, "jets");
 49
 50      for (size_t i = 0; i < 2; ++i) {
 51        book(_h_NjetIncl[i] ,1, 1, i+1);
 52        book(_h_RatioNjetIncl[i], 2, 1, i+1);
 53        book(_h_FirstJetPt_1jet[i] ,3, 1, i+1);
 54        book(_h_FirstJetPt_2jet[i] ,4, 1, i+1);
 55        book(_h_FirstJetPt_3jet[i] ,5, 1, i+1);
 56        book(_h_FirstJetPt_4jet[i] ,6, 1, i+1);
 57        book(_h_SecondJetPt_2jet[i] ,7, 1, i+1);
 58        book(_h_SecondJetPt_3jet[i] ,8, 1, i+1);
 59        book(_h_SecondJetPt_4jet[i] ,9, 1, i+1);
 60        book(_h_ThirdJetPt_3jet[i] ,10, 1, i+1);
 61        book(_h_ThirdJetPt_4jet[i] ,11, 1, i+1);
 62        book(_h_FourthJetPt_4jet[i] ,12, 1, i+1);
 63        book(_h_Ht_1jet[i] ,13, 1, i+1);
 64        book(_h_Ht_2jet[i] ,14, 1, i+1);
 65        book(_h_Ht_3jet[i] ,15, 1, i+1);
 66        book(_h_Ht_4jet[i] ,16, 1, i+1);
 67        book(_h_Minv_2jet[i] ,17, 1, i+1);
 68        book(_h_Minv_3jet[i] ,18, 1, i+1);
 69        book(_h_Minv_4jet[i] ,19, 1, i+1);
 70        book(_h_JetRapidity[i] ,20, 1, i+1);
 71        book(_h_DeltaYElecJet[i] ,21, 1, i+1);
 72        book(_h_SumYElecJet[i] ,22, 1, i+1);
 73        book(_h_DeltaR_2jet[i] ,23, 1, i+1);
 74        book(_h_DeltaY_2jet[i] ,24, 1, i+1);
 75        book(_h_DeltaPhi_2jet[i] ,25, 1, i+1);
 76      }
 77    }
 78
 79
 80    /// Perform the per-event analysis
 81    void analyze(const Event& event) {
 82      const vector<DressedLepton>& leptons = apply<DressedLeptons>(event, "leptons").dressedLeptons();
 83      Particles neutrinos = apply<FinalState>(event, "neutrinos").particlesByPt();
 84
 85      if (leptons.size() != 1 || (neutrinos.size() == 0)) vetoEvent;
 86
 87      FourMomentum lepton = leptons[0].momentum();
 88      FourMomentum p_miss = neutrinos[0].momentum();
 89      if (p_miss.Et() < 25.0*GeV) vetoEvent;
 90
 91      double mT = sqrt(2.0 * lepton.pT() * p_miss.Et() * (1.0 - cos( lepton.phi()-p_miss.phi()) ) );
 92      if (mT < 40.0*GeV) vetoEvent;
 93
 94      double jetcuts[] = { 30.0*GeV, 20.0*GeV };
 95      const FastJets& jetpro = apply<FastJets>(event, "jets");
 96
 97      for (size_t i = 0; i < 2; ++i) {
 98        vector<FourMomentum> jets;
 99        double HT = lepton.pT() + p_miss.pT();
100        for (const Jet& jet : jetpro.jetsByPt(jetcuts[i])) {
101          if (jet.absrap() < 4.4 && deltaR(lepton, jet.momentum()) > 0.5) {
102            jets.push_back(jet.momentum());
103            HT += jet.pT();
104          }
105        }
106
107        _h_NjetIncl[i]->fill(0.0);
108
109        // Njet>=1 observables
110        if (jets.size() < 1) continue;
111        _h_NjetIncl[i]->fill(1.0);
112        _h_FirstJetPt_1jet[i]->fill(jets[0].pT());
113        _h_JetRapidity[i]->fill(jets[0].rapidity());
114        _h_Ht_1jet[i]->fill(HT);
115        _h_DeltaYElecJet[i]->fill(lepton.rapidity()-jets[0].rapidity());
116        _h_SumYElecJet[i]->fill(lepton.rapidity()+jets[0].rapidity());
117
118        // Njet>=2 observables
119        if (jets.size() < 2) continue;
120        _h_NjetIncl[i]->fill(2.0);
121        _h_FirstJetPt_2jet[i]->fill(jets[0].pT());
122        _h_SecondJetPt_2jet[i]->fill(jets[1].pT());
123        _h_Ht_2jet[i]->fill(HT);
124        double m2_2jet = FourMomentum(jets[0]+jets[1]).mass2();
125        _h_Minv_2jet[i]->fill(m2_2jet>0.0 ? sqrt(m2_2jet) : 0.0);
126        _h_DeltaR_2jet[i]->fill(deltaR(jets[0], jets[1]));
127        _h_DeltaY_2jet[i]->fill(jets[0].rapidity()-jets[1].rapidity());
128        _h_DeltaPhi_2jet[i]->fill(deltaPhi(jets[0], jets[1]));
129
130        // Njet>=3 observables
131        if (jets.size() < 3) continue;
132        _h_NjetIncl[i]->fill(3.0);
133        _h_FirstJetPt_3jet[i]->fill(jets[0].pT());
134        _h_SecondJetPt_3jet[i]->fill(jets[1].pT());
135        _h_ThirdJetPt_3jet[i]->fill(jets[2].pT());
136        _h_Ht_3jet[i]->fill(HT);
137        double m2_3jet = FourMomentum(jets[0]+jets[1]+jets[2]).mass2();
138        _h_Minv_3jet[i]->fill(m2_3jet>0.0 ? sqrt(m2_3jet) : 0.0);
139
140        // Njet>=4 observables
141        if (jets.size() < 4) continue;
142        _h_NjetIncl[i]->fill(4.0);
143        _h_FirstJetPt_4jet[i]->fill(jets[0].pT());
144        _h_SecondJetPt_4jet[i]->fill(jets[1].pT());
145        _h_ThirdJetPt_4jet[i]->fill(jets[2].pT());
146        _h_FourthJetPt_4jet[i]->fill(jets[3].pT());
147        _h_Ht_4jet[i]->fill(HT);
148        double m2_4jet = FourMomentum(jets[0]+jets[1]+jets[2]+jets[3]).mass2();
149        _h_Minv_4jet[i]->fill(m2_4jet>0.0 ? sqrt(m2_4jet) : 0.0);
150
151        // Njet>=5 observables
152        if (jets.size() < 5) continue;
153        _h_NjetIncl[i]->fill(5.0);
154      }
155    }
156
157
158    /// Normalise histograms etc., after the run
159    void finalize() {
160      for (size_t i = 0; i < 2; ++i) {
161
162        // Construct jet multiplicity ratio
163        for (size_t n = 1; n < _h_NjetIncl[i]->numBins(); ++n) {
164          YODA::HistoBin1D& b0 = _h_NjetIncl[i]->bin(n-1);
165          YODA::HistoBin1D& b1 = _h_NjetIncl[i]->bin(n);
166          double val = 0.0, err= 0.0;
167          if (b0.height() && b1.height()) {
168            val = b1.height() / b0.height();
169            err = b1.height() / b0.height() * (b0.relErr() + b1.relErr());
170          }
171          _h_RatioNjetIncl[i]->addPoint(n, val, 0.5, err);
172        }
173
174        // Scale all histos to the cross section
175        const double factor = crossSection()/sumOfWeights();
176        scale(_h_DeltaPhi_2jet[i], factor);
177        scale(_h_DeltaR_2jet[i], factor);
178        scale(_h_DeltaY_2jet[i], factor);
179        scale(_h_DeltaYElecJet[i], factor);
180        scale(_h_FirstJetPt_1jet[i], factor);
181        scale(_h_FirstJetPt_2jet[i], factor);
182        scale(_h_FirstJetPt_3jet[i], factor);
183        scale(_h_FirstJetPt_4jet[i], factor);
184        scale(_h_FourthJetPt_4jet[i], factor);
185        scale(_h_Ht_1jet[i], factor);
186        scale(_h_Ht_2jet[i], factor);
187        scale(_h_Ht_3jet[i], factor);
188        scale(_h_Ht_4jet[i], factor);
189        scale(_h_JetRapidity[i], factor);
190        scale(_h_Minv_2jet[i], factor);
191        scale(_h_Minv_3jet[i], factor);
192        scale(_h_Minv_4jet[i], factor);
193        scale(_h_NjetIncl[i], factor);
194        scale(_h_SecondJetPt_2jet[i], factor);
195        scale(_h_SecondJetPt_3jet[i], factor);
196        scale(_h_SecondJetPt_4jet[i], factor);
197        scale(_h_SumYElecJet[i], factor);
198        scale(_h_ThirdJetPt_3jet[i], factor);
199        scale(_h_ThirdJetPt_4jet[i], factor);
200      }
201    }
202
203    //@}
204
205
206  private:
207
208    /// @name Histograms
209    //@{
210    Histo1DPtr _h_DeltaPhi_2jet[2];
211    Histo1DPtr _h_DeltaR_2jet[2];
212    Histo1DPtr _h_DeltaY_2jet[2];
213    Histo1DPtr _h_DeltaYElecJet[2];
214    Histo1DPtr _h_FirstJetPt_1jet[2];
215    Histo1DPtr _h_FirstJetPt_2jet[2];
216    Histo1DPtr _h_FirstJetPt_3jet[2];
217    Histo1DPtr _h_FirstJetPt_4jet[2];
218    Histo1DPtr _h_FourthJetPt_4jet[2];
219    Histo1DPtr _h_Ht_1jet[2];
220    Histo1DPtr _h_Ht_2jet[2];
221    Histo1DPtr _h_Ht_3jet[2];
222    Histo1DPtr _h_Ht_4jet[2];
223    Histo1DPtr _h_JetRapidity[2];
224    Histo1DPtr _h_Minv_2jet[2];
225    Histo1DPtr _h_Minv_3jet[2];
226    Histo1DPtr _h_Minv_4jet[2];
227    Histo1DPtr _h_NjetIncl[2];
228    Scatter2DPtr _h_RatioNjetIncl[2];
229    Histo1DPtr _h_SecondJetPt_2jet[2];
230    Histo1DPtr _h_SecondJetPt_3jet[2];
231    Histo1DPtr _h_SecondJetPt_4jet[2];
232    Histo1DPtr _h_SumYElecJet[2];
233    Histo1DPtr _h_ThirdJetPt_3jet[2];
234    Histo1DPtr _h_ThirdJetPt_4jet[2];
235    //@}
236
237
238  };
239
240
241
242  // The hook for the plugin system
243  RIVET_DECLARE_PLUGIN(ATLAS_2012_I1083318);
244
245}