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