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, FastJets::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      const double weight = 1.0;
 82
 83      // get the jet candidates
 84      Jets cand_jets;
 85      for (const Jet& jet :
 86               apply<FastJets>(event, "AntiKtJets04").jetsByPt(20.0*GeV) ) {
 87        if ( fabs( jet.eta() ) < 2.8 ) {
 88          cand_jets.push_back(jet);
 89        }
 90      }
 91
 92      // candidate muons
 93      Particles cand_mu;
 94      Particles chg_tracks =
 95        apply<ChargedFinalState>(event, "cfs").particles();
 96      for ( const Particle& mu :
 97                apply<IdentifiedFinalState>(event, "muons").particlesByPt() ) {
 98        double pTinCone = -mu.pT();
 99        for ( const Particle& track : chg_tracks ) {
100          if ( deltaR(mu.momentum(),track.momentum()) <= 0.2 )
101            pTinCone += track.pT();
102        }
103        if ( pTinCone < 1.8*GeV )
104          cand_mu.push_back(mu);
105      }
106
107      // candidate electrons
108      Particles cand_e  =
109        apply<IdentifiedFinalState>(event, "elecs").particlesByPt();
110
111      // resolve jet/lepton ambiguity
112      Jets recon_jets;
113      for ( const Jet& jet : cand_jets ) {
114        // candidates after |eta| < 2.8
115        if ( fabs( jet.eta() ) >= 2.8 ) continue;
116        bool away_from_e = true;
117        for ( const Particle& e : cand_e ) {
118          if ( deltaR(e.momentum(),jet.momentum()) <= 0.2 ) {
119            away_from_e = false;
120            break;
121          }
122        }
123        if ( away_from_e ) recon_jets.push_back( jet );
124      }
125
126      // only keep electrons more than R=0.4 from jets
127      Particles recon_e;
128      for ( const Particle& e : cand_e ) {
129        bool away = true;
130        for ( const Jet& jet : recon_jets ) {
131          if ( deltaR(e.momentum(),jet.momentum()) < 0.4 ) {
132            away = false;
133            break;
134          }
135        }
136        if ( away )
137          recon_e.push_back( e );
138      }
139
140      // only keep muons more than R=0.4 from jets
141      Particles recon_mu;
142      for ( const Particle& mu : cand_mu ) {
143        bool away = true;
144        for ( const Jet& jet : recon_jets ) {
145          if ( deltaR(mu.momentum(),jet.momentum()) < 0.4 ) {
146            away = false;
147            break;
148          }
149        }
150        if ( away )
151          recon_mu.push_back( mu );
152      }
153
154      // pTmiss
155      Particles vfs_particles =
156        apply<VisibleFinalState>(event, "vfs").particles();
157      FourMomentum pTmiss;
158      for ( const Particle& p  : vfs_particles ) {
159        pTmiss -= p.momentum();
160      }
161      double eTmiss = pTmiss.pT();
162
163      // now only use recon_jets, recon_mu, recon_e
164
165      // reject events with electrons and muons
166      if ( ! ( recon_mu.empty() && recon_e.empty() ) ) {
167        MSG_DEBUG("Charged leptons left after selection");
168        vetoEvent;
169      }
170
171      // calculate H_T
172      double HT=0;
173      for ( const Jet& jet : recon_jets ) {
174        if ( jet.pT() > 40 * GeV )
175          HT += jet.pT() ;
176      }
177
178      // number of jets
179      unsigned int njet55=0, njet80=0;
180      for (unsigned int ix=0;ix<recon_jets.size();++ix) {
181        if(recon_jets[ix].pT()>80.*GeV) ++njet80;
182        if(recon_jets[ix].pT()>55.*GeV) ++njet55;
183      }
184
185      if(njet55==0) vetoEvent;
186
187      double ratio = eTmiss/sqrt(HT);
188
189      if(ratio>4.) {
190        _hist_njet55->fill(njet55,weight);
191        _hist_njet80->fill(njet80,weight);
192        // 7j55
193        if(njet55>=7)
194          _count_7j55->fill( 0.5, weight);
195        // 8j55
196        if(njet55>=8)
197          _count_8j55->fill( 0.5, weight) ;
198        // 8j55
199        if(njet55>=9)
200          _count_9j55->fill( 0.5, weight) ;
201        // 6j80
202        if(njet80>=6)
203          _count_6j80->fill( 0.5, weight) ;
204        // 7j80
205        if(njet80>=7)
206          _count_7j80->fill( 0.5, weight) ;
207        // 8j80
208        if(njet80>=8)
209          _count_8j80->fill( 0.5, weight) ;
210      }
211
212      if(njet55>=7)
213        _etmiss_HT_7j55->fill( ratio, weight);
214      // 8j55
215      if(njet55>=8)
216        _etmiss_HT_8j55->fill( ratio, weight) ;
217      // 8j55
218      if(njet55>=9)
219        _etmiss_HT_9j55->fill( ratio, weight) ;
220      // 6j80
221      if(njet80>=6)
222        _etmiss_HT_6j80->fill( ratio, weight) ;
223      // 7j80
224      if(njet80>=7)
225        _etmiss_HT_7j80->fill( ratio, weight) ;
226      // 8j80
227      if(njet80>=8)
228        _etmiss_HT_8j80->fill( ratio, weight) ;
229
230    }
231
232    //@}
233
234    void finalize() {
235      double norm = crossSection()/femtobarn*4.7/sumOfWeights();
236
237      scale(_etmiss_HT_7j55,2.*norm);
238      scale(_etmiss_HT_8j55,2.*norm);
239      scale(_etmiss_HT_9j55,2.*norm);
240      scale(_etmiss_HT_6j80,2.*norm);
241      scale(_etmiss_HT_7j80,2.*norm);
242      scale(_etmiss_HT_8j80,2.*norm);
243
244      scale(_hist_njet55,norm);
245      scale(_hist_njet80,norm);
246
247      scale(_count_7j55,norm);
248      scale(_count_8j55,norm);
249      scale(_count_9j55,norm);
250      scale(_count_6j80,norm);
251      scale(_count_7j80,norm);
252      scale(_count_8j80,norm);
253    }
254
255  private:
256
257    /// @name Histograms
258    //@{
259    Histo1DPtr _etmiss_HT_7j55;
260    Histo1DPtr _etmiss_HT_8j55;
261    Histo1DPtr _etmiss_HT_9j55;
262    Histo1DPtr _etmiss_HT_6j80;
263    Histo1DPtr _etmiss_HT_7j80;
264    Histo1DPtr _etmiss_HT_8j80;
265
266    Histo1DPtr _hist_njet55;
267    Histo1DPtr _hist_njet80;
268
269    Histo1DPtr _count_7j55;
270    Histo1DPtr _count_8j55;
271    Histo1DPtr _count_9j55;
272    Histo1DPtr _count_6j80;
273    Histo1DPtr _count_7j80;
274    Histo1DPtr _count_8j80;
275    //@}
276
277  };
278
279  // The hook for the plugin system
280  RIVET_DECLARE_PLUGIN(ATLAS_2012_I1117704);
281
282}