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_I890749 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2011_I890749);
 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, JetAlg::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
 71      Particles veto_e = apply<IdentifiedFinalState>(event, "veto_elecs").particles();
 72      if ( ! veto_e.empty() ) {
 73        MSG_DEBUG("electrons in veto region");
 74        vetoEvent;
 75      }
 76
 77
 78      Jets cand_jets = apply<FastJets>(event, "AntiKtJets04").jetsByPt(Cuts::pT > 20*GeV && Cuts::abseta < 4.9);
 79      Particles cand_e  = apply<IdentifiedFinalState>(event, "elecs").particlesByPt();
 80
 81      Particles cand_mu;
 82      Particles chg_tracks = apply<ChargedFinalState>(event, "cfs").particles();
 83      for ( const Particle & mu : apply<IdentifiedFinalState>(event, "muons").particlesByPt() ) {
 84        double pTinCone = -mu.pT();
 85        for ( const Particle & track : chg_tracks ) {
 86          if ( deltaR(mu, track) <= 0.2 ) pTinCone += track.pT();
 87        }
 88        if ( pTinCone < 1.8*GeV ) cand_mu.push_back(mu);
 89      }
 90
 91      Jets cand_jets_2;
 92      for (const Jet& jet : cand_jets) {
 93        if (jet.abseta() >= 2.5)
 94          cand_jets_2.push_back(jet);
 95        else {
 96          bool away_from_e = true;
 97          for (const Particle& e : cand_e) {
 98            if (deltaR(e, jet) <= 0.2 ) {
 99              away_from_e = false;
100              break;
101            }
102          }
103          if ( away_from_e )
104            cand_jets_2.push_back( jet );
105        }
106      }
107
108      Particles recon_e, recon_mu;
109      for ( const Particle & e : cand_e ) {
110        bool away = true;
111        for ( const Jet& jet : cand_jets_2 ) {
112          if ( deltaR(e, jet) < 0.4 ) {
113            away = false;
114            break;
115          }
116        }
117        if ( away )
118          recon_e.push_back( e );
119      }
120
121      for ( const Particle & mu : cand_mu ) {
122        bool away = true;
123        for ( const Jet& jet : cand_jets_2 ) {
124          if ( deltaR(mu, jet) < 0.4 ) {
125            away = false;
126            break;
127          }
128        }
129        if ( away )
130          recon_mu.push_back( mu );
131      }
132
133
134      // pTmiss
135      Particles vfs_particles = 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
143      // final jet filter
144      Jets recon_jets;
145      for ( const Jet& jet : cand_jets_2 ) {
146        if ( jet.abseta() <= 2.5 ) recon_jets.push_back( jet );
147      }
148
149
150      // now only use recon_jets, recon_mu, recon_e
151
152      if ( ! ( recon_mu.empty() && recon_e.empty() ) ) {
153        MSG_DEBUG("Charged leptons left after selection");
154        vetoEvent;
155      }
156
157      if ( eTmiss <= 100*GeV ) {
158        MSG_DEBUG("Not enough eTmiss: " << eTmiss << " < 100");
159        vetoEvent;
160      }
161
162
163      if ( recon_jets.empty() || recon_jets[0].pT() <= 120.0*GeV ) {
164        MSG_DEBUG("No hard leading jet in " << recon_jets.size() << " jets");
165        vetoEvent;
166      }
167
168      // ==================== observables ====================
169
170      // Njets, min_dPhi
171
172      int Njets = 0;
173      double min_dPhi = 999.999;
174      double pTmiss_phi = pTmiss.phi();
175      for ( const Jet& jet : recon_jets ) {
176        if ( jet.pT() > 40 * GeV ) {
177          if ( Njets < 3 )
178            min_dPhi = min( min_dPhi, deltaPhi( pTmiss_phi, jet.phi() ) );
179          ++Njets;
180        }
181      }
182
183      if ( Njets < 2 ) {
184        MSG_DEBUG("Only " << Njets << " >40 GeV jets left");
185        vetoEvent;
186      }
187
188      if ( min_dPhi <= 0.4 ) {
189        MSG_DEBUG("dPhi too small");
190        vetoEvent;
191      }
192
193      // m_eff
194
195      double m_eff_2j = eTmiss
196        + recon_jets[0].pT()
197        + recon_jets[1].pT();
198
199      double m_eff_3j = recon_jets.size() < 3 ? -999.0 : m_eff_2j + recon_jets[2].pT();
200
201      // etmiss / m_eff
202
203      double et_meff_2j = eTmiss / m_eff_2j;
204      double et_meff_3j = eTmiss / m_eff_3j;
205
206      FourMomentum a = recon_jets[0].momentum();
207      FourMomentum b = recon_jets[1].momentum();
208
209      double m_T2 = mT2( a, b, pTmiss, 0.0 ); // zero mass invisibles
210
211
212      // ==================== FILL ====================
213
214      MSG_DEBUG( "Trying to fill "
215                 << Njets << ' '
216                 << m_eff_2j << ' '
217                 << et_meff_2j << ' '
218                 << m_eff_3j << ' '
219                 << et_meff_3j << ' '
220                 << m_T2 );
221
222      _hist_eTmiss->fill(eTmiss);
223
224      // AAAAAAAAAA
225      if ( et_meff_2j > 0.3 ) {
226        _hist_meff_A->fill(m_eff_2j);
227        if ( m_eff_2j > 500 * GeV ) {
228          MSG_DEBUG("Hits A");
229          _count_A->fill(0.5);
230        }
231      }
232
233      // BBBBBBBBBB
234      _hist_mT2_B->fill(m_T2);
235      if ( m_T2 > 300 * GeV ) {
236        MSG_DEBUG("Hits B");
237        _count_B->fill(0.5);
238      }
239
240      // need 3 jets for C and D
241      if ( Njets >= 3 && et_meff_3j > 0.25 ) {
242
243        _hist_meff_CD->fill(m_eff_3j);
244
245        // CCCCCCCCCC
246        if ( m_eff_3j > 500 * GeV ) {
247          MSG_DEBUG("Hits C");
248          _count_C->fill(0.5);
249        }
250
251        // DDDDDDDDDD
252        if ( m_eff_3j > 1000 * GeV ) {
253          MSG_DEBUG("Hits D");
254          _count_D->fill(0.5);
255        }
256      }
257
258    }
259
260    /// @}
261
262    void finalize() {
263
264      double norm = crossSection()/picobarn*35.0/sumOfWeights();
265      scale(_hist_meff_A ,100.*norm);
266      scale(_hist_mT2_B  ,100.*norm);
267      scale(_hist_meff_CD, 40.*norm);
268      scale(_hist_eTmiss , 50.*norm);
269    }
270
271
272  private:
273
274    /// @name Histograms
275    /// @{
276    Histo1DPtr _count_A;
277    Histo1DPtr _count_B;
278    Histo1DPtr _count_C;
279    Histo1DPtr _count_D;
280    Histo1DPtr _hist_meff_A;
281    Histo1DPtr _hist_mT2_B;
282    Histo1DPtr _hist_meff_CD;
283    Histo1DPtr _hist_eTmiss;
284    /// @}
285
286  };
287
288
289
290  RIVET_DECLARE_ALIASED_PLUGIN(ATLAS_2011_I890749, ATLAS_2011_S8983313);
291
292}