1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/ChargedFinalState.hh"
  5#include "Rivet/Projections/VisibleFinalState.hh"
  6#include "Rivet/Projections/VetoedFinalState.hh"
  7#include "Rivet/Projections/IdentifiedFinalState.hh"
  8#include "Rivet/Projections/FastJets.hh"
  9#include "Rivet/Tools/RivetMT2.hh"
 10
 11namespace Rivet {
 12
 13
 14  class ATLAS_2012_I1117704 : public Analysis {
 15  public:
 16
 17    /// @name Constructors etc.
 18    /// @{
 19
 20    /// Constructor
 21    ATLAS_2012_I1117704()
 22      : Analysis("ATLAS_2012_I1117704")
 23    {    }
 24
 25    /// @}
 26
 27
 28  public:
 29
 30    /// @name Analysis methods
 31    /// @{
 32
 33    /// Book histograms and initialise projections before the run
 34    void init() {
 35
 36      // projection to find the electrons
 37      IdentifiedFinalState elecs(Cuts::abseta < 2.47 && Cuts::pT > 20*GeV);
 38      elecs.acceptIdPair(PID::ELECTRON);
 39      declare(elecs, "elecs");
 40
 41      // projection to find the muons
 42      IdentifiedFinalState muons(Cuts::abseta < 2.4 && Cuts::pT > 10*GeV);
 43      muons.acceptIdPair(PID::MUON);
 44      declare(muons, "muons");
 45
 46      // for pTmiss
 47      declare(VisibleFinalState(Cuts::abseta < 4.9), "vfs");
 48
 49      VetoedFinalState vfs;
 50      vfs.addVetoPairId(PID::MUON);
 51
 52      /// Jet finder
 53      declare(FastJets(vfs, JetAlg::ANTIKT, 0.4), "AntiKtJets04");
 54
 55      // all tracks (to do deltaR with leptons)
 56      declare(ChargedFinalState(Cuts::abseta < 3),"cfs");
 57
 58      /// Book histograms
 59      book(_etmiss_HT_7j55 ,"etmiss_HT_7j55", 8, 0., 16.);
 60      book(_etmiss_HT_8j55 ,"etmiss_HT_8j55", 8, 0., 16.);
 61      book(_etmiss_HT_9j55 ,"etmiss_HT_9j55", 8, 0., 16.);
 62      book(_etmiss_HT_6j80 ,"etmiss_HT_6j80", 8, 0., 16.);
 63      book(_etmiss_HT_7j80 ,"etmiss_HT_7j80", 8, 0., 16.);
 64      book(_etmiss_HT_8j80 ,"etmiss_HT_8j80", 8, 0., 16.);
 65
 66      book(_hist_njet55 ,"hist_njet55", 11, 2.5, 13.5);
 67      book(_hist_njet80 ,"hist_njet80", 11, 2.5, 13.5);
 68
 69      book(_count_7j55 ,"count_7j55", 1, 0., 1.);
 70      book(_count_8j55 ,"count_8j55", 1, 0., 1.);
 71      book(_count_9j55 ,"count_9j55", 1, 0., 1.);
 72      book(_count_6j80 ,"count_6j80", 1, 0., 1.);
 73      book(_count_7j80 ,"count_7j80", 1, 0., 1.);
 74      book(_count_8j80 ,"count_8j80", 1, 0., 1.);
 75
 76    }
 77
 78
 79    /// Perform the per-event analysis
 80    void analyze(const Event& event) {
 81
 82      // get the jet candidates
 83      Jets cand_jets = apply<FastJets>(event, "AntiKtJets04").jetsByPt(Cuts::pT > 20.0*GeV && Cuts::abseta < 2.8);
 84
 85      // candidate muons
 86      Particles cand_mu;
 87      Particles chg_tracks =
 88        apply<ChargedFinalState>(event, "cfs").particles();
 89      for ( const Particle& mu :
 90                apply<IdentifiedFinalState>(event, "muons").particlesByPt() ) {
 91        double pTinCone = -mu.pT();
 92        for ( const Particle& track : chg_tracks ) {
 93          if ( deltaR(mu.momentum(),track.momentum()) <= 0.2 )
 94            pTinCone += track.pT();
 95        }
 96        if ( pTinCone < 1.8*GeV )
 97          cand_mu.push_back(mu);
 98      }
 99
100      // candidate electrons
101      Particles cand_e  =
102        apply<IdentifiedFinalState>(event, "elecs").particlesByPt();
103
104      // resolve jet/lepton ambiguity
105      Jets recon_jets;
106      for ( const Jet& jet : cand_jets ) {
107        // candidates after |eta| < 2.8
108        if ( fabs( jet.eta() ) >= 2.8 ) continue;
109        bool away_from_e = true;
110        for ( const Particle& e : cand_e ) {
111          if ( deltaR(e.momentum(),jet.momentum()) <= 0.2 ) {
112            away_from_e = false;
113            break;
114          }
115        }
116        if ( away_from_e ) recon_jets.push_back( jet );
117      }
118
119      // only keep electrons more than R=0.4 from jets
120      Particles recon_e;
121      for ( const Particle& e : cand_e ) {
122        bool away = true;
123        for ( const Jet& jet : recon_jets ) {
124          if ( deltaR(e.momentum(),jet.momentum()) < 0.4 ) {
125            away = false;
126            break;
127          }
128        }
129        if ( away )
130          recon_e.push_back( e );
131      }
132
133      // only keep muons more than R=0.4 from jets
134      Particles recon_mu;
135      for ( const Particle& mu : cand_mu ) {
136        bool away = true;
137        for ( const Jet& jet : recon_jets ) {
138          if ( deltaR(mu.momentum(),jet.momentum()) < 0.4 ) {
139            away = false;
140            break;
141          }
142        }
143        if ( away )
144          recon_mu.push_back( mu );
145      }
146
147      // pTmiss
148      Particles vfs_particles =
149        apply<VisibleFinalState>(event, "vfs").particles();
150      FourMomentum pTmiss;
151      for ( const Particle& p  : vfs_particles ) {
152        pTmiss -= p.momentum();
153      }
154      double eTmiss = pTmiss.pT();
155
156      // now only use recon_jets, recon_mu, recon_e
157
158      // reject events with electrons and muons
159      if ( ! ( recon_mu.empty() && recon_e.empty() ) ) {
160        MSG_DEBUG("Charged leptons left after selection");
161        vetoEvent;
162      }
163
164      // calculate H_T
165      double HT=0;
166      for ( const Jet& jet : recon_jets ) {
167        if ( jet.pT() > 40 * GeV )
168          HT += jet.pT() ;
169      }
170
171      // number of jets
172      unsigned int njet55=0, njet80=0;
173      for (unsigned int ix=0;ix<recon_jets.size();++ix) {
174        if(recon_jets[ix].pT()>80.*GeV) ++njet80;
175        if(recon_jets[ix].pT()>55.*GeV) ++njet55;
176      }
177
178      if(njet55==0) vetoEvent;
179
180      double ratio = eTmiss/sqrt(HT);
181
182      if(ratio>4.) {
183        _hist_njet55->fill(njet55);
184        _hist_njet80->fill(njet80);
185        // 7j55
186        if(njet55>=7)
187          _count_7j55->fill( 0.5);
188        // 8j55
189        if(njet55>=8)
190          _count_8j55->fill( 0.5) ;
191        // 8j55
192        if(njet55>=9)
193          _count_9j55->fill( 0.5) ;
194        // 6j80
195        if(njet80>=6)
196          _count_6j80->fill( 0.5) ;
197        // 7j80
198        if(njet80>=7)
199          _count_7j80->fill( 0.5) ;
200        // 8j80
201        if(njet80>=8)
202          _count_8j80->fill( 0.5) ;
203      }
204
205      if(njet55>=7)
206        _etmiss_HT_7j55->fill( ratio);
207      // 8j55
208      if(njet55>=8)
209        _etmiss_HT_8j55->fill( ratio) ;
210      // 8j55
211      if(njet55>=9)
212        _etmiss_HT_9j55->fill( ratio) ;
213      // 6j80
214      if(njet80>=6)
215        _etmiss_HT_6j80->fill( ratio) ;
216      // 7j80
217      if(njet80>=7)
218        _etmiss_HT_7j80->fill( ratio) ;
219      // 8j80
220      if(njet80>=8)
221        _etmiss_HT_8j80->fill( ratio) ;
222
223    }
224
225    /// @}
226
227    void finalize() {
228      double norm = crossSection()/femtobarn*4.7/sumOfWeights();
229
230      scale(_etmiss_HT_7j55,2.*norm);
231      scale(_etmiss_HT_8j55,2.*norm);
232      scale(_etmiss_HT_9j55,2.*norm);
233      scale(_etmiss_HT_6j80,2.*norm);
234      scale(_etmiss_HT_7j80,2.*norm);
235      scale(_etmiss_HT_8j80,2.*norm);
236
237      scale(_hist_njet55,norm);
238      scale(_hist_njet80,norm);
239
240      scale(_count_7j55,norm);
241      scale(_count_8j55,norm);
242      scale(_count_9j55,norm);
243      scale(_count_6j80,norm);
244      scale(_count_7j80,norm);
245      scale(_count_8j80,norm);
246    }
247
248  private:
249
250    /// @name Histograms
251    /// @{
252    Histo1DPtr _etmiss_HT_7j55;
253    Histo1DPtr _etmiss_HT_8j55;
254    Histo1DPtr _etmiss_HT_9j55;
255    Histo1DPtr _etmiss_HT_6j80;
256    Histo1DPtr _etmiss_HT_7j80;
257    Histo1DPtr _etmiss_HT_8j80;
258
259    Histo1DPtr _hist_njet55;
260    Histo1DPtr _hist_njet80;
261
262    Histo1DPtr _count_7j55;
263    Histo1DPtr _count_8j55;
264    Histo1DPtr _count_9j55;
265    Histo1DPtr _count_6j80;
266    Histo1DPtr _count_7j80;
267    Histo1DPtr _count_8j80;
268    /// @}
269
270  };
271
272  RIVET_DECLARE_PLUGIN(ATLAS_2012_I1117704);
273
274}