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
 10namespace Rivet {
 11
 12
 13  class ATLAS_2011_S8983313 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2011_S8983313);
 18
 19
 20    /// @name Analysis methods
 21    //@{
 22
 23    /// Book histograms and initialise projections before the run
 24    void init() {
 25
 26      // projection to find the electrons
 27      Cut pt10 = Cuts::pT > 10.0*GeV;
 28      IdentifiedFinalState elecs(Cuts::abseta < 2.47 && pt10);
 29      elecs.acceptIdPair(PID::ELECTRON);
 30      declare(elecs, "elecs");
 31
 32      // veto region electrons
 33      Cut vetocut = Cuts::absetaIn(1.37, 1.52);
 34      IdentifiedFinalState veto_elecs(vetocut && pt10);
 35      veto_elecs.acceptIdPair(PID::ELECTRON);
 36      declare(veto_elecs, "veto_elecs");
 37
 38      // projection to find the muons
 39      IdentifiedFinalState muons(Cuts::abseta < 2.4 && pt10 );
 40      muons.acceptIdPair(PID::MUON);
 41      declare(muons, "muons");
 42
 43      VetoedFinalState vfs;
 44      vfs.addVetoPairId(PID::MUON);
 45
 46      /// Jet finder
 47      declare(FastJets(vfs, FastJets::ANTIKT, 0.4), "AntiKtJets04");
 48
 49      // all tracks (to do deltaR with leptons)
 50      declare(ChargedFinalState(Cuts::abseta < 3.0),"cfs");
 51
 52      // for pTmiss
 53      declare(VisibleFinalState(Cuts::abseta < 4.9),"vfs");
 54
 55      /// Book histograms
 56      book(_count_A ,"count_A", 1, 0., 1.);
 57      book(_count_B ,"count_B", 1, 0., 1.);
 58      book(_count_C ,"count_C", 1, 0., 1.);
 59      book(_count_D ,"count_D", 1, 0., 1.);
 60
 61      book(_hist_meff_A  ,"m_eff_A", 30, 0., 3000.);
 62      book(_hist_mT2_B   ,"m_T2", 25, 0., 1000.);
 63      book(_hist_meff_CD ,"m_eff_C_D", 30, 0., 3000.);
 64      book(_hist_eTmiss  ,"Et_miss", 20, 0., 1000.);
 65    }
 66
 67
 68    /// Perform the per-event analysis
 69    void analyze(const Event& event) {
 70      const double weight = 1.0;
 71
 72      Particles veto_e = apply<IdentifiedFinalState>(event, "veto_elecs").particles();
 73      if ( ! veto_e.empty() ) {
 74        MSG_DEBUG("electrons in veto region");
 75        vetoEvent;
 76      }
 77
 78
 79      Jets cand_jets = apply<FastJets>(event, "AntiKtJets04").jetsByPt(Cuts::pT > 20*GeV && Cuts::abseta < 4.9);
 80      Particles cand_e  = apply<IdentifiedFinalState>(event, "elecs").particlesByPt();
 81
 82      Particles cand_mu;
 83      Particles chg_tracks = apply<ChargedFinalState>(event, "cfs").particles();
 84      for ( const Particle & mu : apply<IdentifiedFinalState>(event, "muons").particlesByPt() ) {
 85        double pTinCone = -mu.pT();
 86        for ( const Particle & track : chg_tracks ) {
 87          if ( deltaR(mu, track) <= 0.2 ) pTinCone += track.pT();
 88        }
 89        if ( pTinCone < 1.8*GeV ) cand_mu.push_back(mu);
 90      }
 91
 92      Jets cand_jets_2;
 93      for (const Jet& jet : cand_jets) {
 94        if (jet.abseta() >= 2.5)
 95          cand_jets_2.push_back(jet);
 96        else {
 97          bool away_from_e = true;
 98          for (const Particle& e : cand_e) {
 99            if (deltaR(e, jet) <= 0.2 ) {
100              away_from_e = false;
101              break;
102            }
103          }
104          if ( away_from_e )
105            cand_jets_2.push_back( jet );
106        }
107      }
108
109      Particles recon_e, recon_mu;
110      for ( const Particle & e : cand_e ) {
111        bool away = true;
112        for ( const Jet& jet : cand_jets_2 ) {
113          if ( deltaR(e, jet) < 0.4 ) {
114            away = false;
115            break;
116          }
117        }
118        if ( away )
119          recon_e.push_back( e );
120      }
121
122      for ( const Particle & mu : cand_mu ) {
123        bool away = true;
124        for ( const Jet& jet : cand_jets_2 ) {
125          if ( deltaR(mu, jet) < 0.4 ) {
126            away = false;
127            break;
128          }
129        }
130        if ( away )
131          recon_mu.push_back( mu );
132      }
133
134
135      // pTmiss
136      Particles vfs_particles = apply<VisibleFinalState>(event, "vfs").particles();
137      FourMomentum pTmiss;
138      for ( const Particle & p : vfs_particles ) {
139        pTmiss -= p.momentum();
140      }
141      double eTmiss = pTmiss.pT();
142
143
144      // final jet filter
145      Jets recon_jets;
146      for ( const Jet& jet : cand_jets_2 ) {
147        if ( jet.abseta() <= 2.5 ) recon_jets.push_back( jet );
148      }
149
150
151      // now only use recon_jets, recon_mu, recon_e
152
153      if ( ! ( recon_mu.empty() && recon_e.empty() ) ) {
154        MSG_DEBUG("Charged leptons left after selection");
155        vetoEvent;
156      }
157
158      if ( eTmiss <= 100*GeV ) {
159        MSG_DEBUG("Not enough eTmiss: " << eTmiss << " < 100");
160        vetoEvent;
161      }
162
163
164      if ( recon_jets.empty() || recon_jets[0].pT() <= 120.0*GeV ) {
165        MSG_DEBUG("No hard leading jet in " << recon_jets.size() << " jets");
166        vetoEvent;
167      }
168
169      // ==================== observables ====================
170
171      // Njets, min_dPhi
172
173      int Njets = 0;
174      double min_dPhi = 999.999;
175      double pTmiss_phi = pTmiss.phi();
176      for ( const Jet& jet : recon_jets ) {
177        if ( jet.pT() > 40 * GeV ) {
178          if ( Njets < 3 )
179            min_dPhi = min( min_dPhi, deltaPhi( pTmiss_phi, jet.phi() ) );
180          ++Njets;
181        }
182      }
183
184      if ( Njets < 2 ) {
185        MSG_DEBUG("Only " << Njets << " >40 GeV jets left");
186        vetoEvent;
187      }
188
189      if ( min_dPhi <= 0.4 ) {
190        MSG_DEBUG("dPhi too small");
191        vetoEvent;
192      }
193
194      // m_eff
195
196      double m_eff_2j = eTmiss
197        + recon_jets[0].pT()
198        + recon_jets[1].pT();
199
200      double m_eff_3j = recon_jets.size() < 3 ? -999.0 : m_eff_2j + recon_jets[2].pT();
201
202      // etmiss / m_eff
203
204      double et_meff_2j = eTmiss / m_eff_2j;
205      double et_meff_3j = eTmiss / m_eff_3j;
206
207      FourMomentum a = recon_jets[0].momentum();
208      FourMomentum b = recon_jets[1].momentum();
209
210      double m_T2 = mT2( a, b, pTmiss, 0.0 ); // zero mass invisibles
211
212
213      // ==================== FILL ====================
214
215      MSG_DEBUG( "Trying to fill "
216                 << Njets << ' '
217                 << m_eff_2j << ' '
218                 << et_meff_2j << ' '
219                 << m_eff_3j << ' '
220                 << et_meff_3j << ' '
221                 << m_T2 );
222
223      _hist_eTmiss->fill(eTmiss, weight);
224
225      // AAAAAAAAAA
226      if ( et_meff_2j > 0.3 ) {
227        _hist_meff_A->fill(m_eff_2j, weight);
228        if ( m_eff_2j > 500 * GeV ) {
229          MSG_DEBUG("Hits A");
230          _count_A->fill(0.5, weight);
231        }
232      }
233
234      // BBBBBBBBBB
235      _hist_mT2_B->fill(m_T2, weight);
236      if ( m_T2 > 300 * GeV ) {
237        MSG_DEBUG("Hits B");
238        _count_B->fill(0.5, weight);
239      }
240
241      // need 3 jets for C and D
242      if ( Njets >= 3 && et_meff_3j > 0.25 ) {
243
244        _hist_meff_CD->fill(m_eff_3j, weight);
245
246        // CCCCCCCCCC
247        if ( m_eff_3j > 500 * GeV ) {
248          MSG_DEBUG("Hits C");
249          _count_C->fill(0.5, weight);
250        }
251
252        // DDDDDDDDDD
253        if ( m_eff_3j > 1000 * GeV ) {
254          MSG_DEBUG("Hits D");
255          _count_D->fill(0.5, weight);
256        }
257      }
258
259    }
260
261    //@}
262
263    void finalize() {
264
265      double norm = crossSection()/picobarn*35.0/sumOfWeights();
266      scale(_hist_meff_A ,100.*norm);
267      scale(_hist_mT2_B  ,100.*norm);
268      scale(_hist_meff_CD, 40.*norm);
269      scale(_hist_eTmiss , 50.*norm);
270    }
271
272
273  private:
274
275    /// @name Histograms
276    //@{
277    Histo1DPtr _count_A;
278    Histo1DPtr _count_B;
279    Histo1DPtr _count_C;
280    Histo1DPtr _count_D;
281    Histo1DPtr _hist_meff_A;
282    Histo1DPtr _hist_mT2_B;
283    Histo1DPtr _hist_meff_CD;
284    Histo1DPtr _hist_eTmiss;
285    //@}
286
287  };
288
289
290
291  RIVET_DECLARE_ALIASED_PLUGIN(ATLAS_2011_S8983313, ATLAS_2011_I890749);
292
293}