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/IdentifiedFinalState.hh"
  7#include "Rivet/Projections/FastJets.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// @brief 0-lepton squark and gluino search
 13  ///
 14  /// @author Chris Wymant
 15  class ATLAS_2011_I930002 : public Analysis {
 16  public:
 17
 18    /// Constructor
 19    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2011_I930002);
 20
 21
 22    /// @name Analysis methods
 23    /// @{
 24
 25    /// Book histograms and initialise projections before the run
 26    void init() {
 27
 28      // Projection to find the electrons
 29      IdentifiedFinalState elecs(Cuts::abseta < 2.47 && Cuts::pT > 20*GeV);
 30      elecs.acceptIdPair(PID::ELECTRON);
 31      declare(elecs, "elecs");
 32
 33      // Projection to find the muons
 34      IdentifiedFinalState muons(Cuts::abseta < 2.4 && Cuts::pT > 10*GeV);
 35      muons.acceptIdPair(PID::MUON);
 36      declare(muons, "muons");
 37
 38      // Jet finder
 39      declare(FastJets(FinalState(), JetAlg::ANTIKT, 0.4), "AntiKtJets04");
 40
 41      // All tracks (to do deltaR with leptons)
 42      declare(ChargedFinalState(Cuts::abseta < 3.0), "cfs");
 43
 44      // Used for pTmiss (N.B. the real 'vfs' extends beyond 4.5 to |eta| = 4.9)
 45      declare(VisibleFinalState(Cuts::abseta < 4.5), "vfs");
 46
 47
 48      // Book histograms
 49      book(_count_2j   ,"count_2j",   1, 0., 1.);
 50      book(_count_3j   ,"count_3j",   1, 0., 1.);
 51      book(_count_4j5  ,"count_4j5",  1, 0., 1.);
 52      book(_count_4j10 ,"count_4j10", 1, 0., 1.);
 53      book(_count_HM   ,"count_HM",   1, 0., 1.);
 54
 55      book(_hist_meff_2j  ,1, 1, 1);
 56      book(_hist_meff_3j  ,2, 1, 1);
 57      book(_hist_meff_4j  ,3, 1, 1);
 58      book(_hist_meff_HM  ,4, 1, 1);
 59
 60      book(_hist_eTmiss  ,"Et_miss", 20, 0., 1000.);
 61    }
 62
 63
 64    /// Perform the per-event analysis
 65    void analyze(const Event& event) {
 66
 67      const Jets& cand_jets = apply<FastJets>(event, "AntiKtJets04").jetsByPt(Cuts::pT > 20*GeV && Cuts::abseta < 4.9);
 68
 69      const Particles cand_e  = apply<IdentifiedFinalState>(event, "elecs").particlesByPt();
 70
 71      // Muon isolation not mentioned in hep-exp 1109.6572, unlike in 1102.5290,
 72      // but assumed to still be applicable
 73      Particles cand_mu;
 74      const Particles chg_tracks = apply<ChargedFinalState>(event, "cfs").particles();
 75      const Particles muons = apply<IdentifiedFinalState>(event, "muons").particlesByPt();
 76      for (const Particle& mu : muons) {
 77        double pTinCone = -mu.pT();
 78        for (const Particle& track : chg_tracks) {
 79          if ( deltaR(mu.momentum(),track.momentum()) <= 0.2 ) {
 80            pTinCone += track.pT();
 81          }
 82        }
 83        if ( pTinCone < 1.8*GeV ) cand_mu.push_back(mu);
 84      }
 85
 86      // Resolve jet-lepton overlap for jets with |eta| < 2.8
 87      Jets cand_jets_2;
 88      for ( const Jet& jet : cand_jets ) {
 89        if ( fabs( jet.eta() ) >= 2.8 ) {
 90          cand_jets_2.push_back( jet );
 91        } else {
 92          bool away_from_e = true;
 93          for ( const Particle & e : cand_e ) {
 94            if ( deltaR(e.momentum(),jet.momentum()) <= 0.2 ) {
 95              away_from_e = false;
 96              break;
 97            }
 98          }
 99          if ( away_from_e ) cand_jets_2.push_back( jet );
100        }
101      }
102
103
104      Particles recon_e, recon_mu;
105
106      for ( const Particle & e : cand_e ) {
107        bool away = true;
108        for ( const Jet& jet : cand_jets_2 ) {
109          if ( deltaR(e.momentum(),jet.momentum()) < 0.4 ) {
110            away = false;
111            break;
112          }
113        }
114        if ( away ) recon_e.push_back( e );
115      }
116
117      for ( const Particle & mu : cand_mu ) {
118        bool away = true;
119        for ( const Jet& jet : cand_jets_2 ) {
120          if ( deltaR(mu.momentum(),jet.momentum()) < 0.4 ) {
121            away = false;
122            break;
123          }
124        }
125        if ( away ) recon_mu.push_back( mu );
126      }
127
128
129      // pTmiss
130      // Based on all candidate electrons, muons and jets, plus everything else with |eta| < 4.5
131      // i.e. everything in our projection "vfs" plus the jets with |eta| > 4.5
132      Particles vfs_particles = apply<VisibleFinalState>(event, "vfs").particles();
133      FourMomentum pTmiss;
134      for ( const Particle & p : vfs_particles ) {
135        pTmiss -= p.momentum();
136      }
137      for ( const Jet& jet : cand_jets_2 ) {
138        if ( fabs( jet.eta() ) > 4.5 ) pTmiss -= jet.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 ( fabs( jet.eta() ) <= 2.8 ) recon_jets.push_back( jet );
147      }
148      // NB. It seems that jets with |eta| > 2.8 could have been thrown away at
149      // the start; we don't do so, in order to follow both the structure of
150      // the paper and the similar Rivet analysis ATLAS_2011_S8983313
151
152      // 'candidate' muons needed only 10 GeV, to cause a veto they need 20 GeV
153      Particles veto_mu;
154      for ( const Particle & mu : cand_mu ) {
155        if ( mu.pT() >= 20.0*GeV ) veto_mu.push_back(mu);
156      }
157
158      if ( ! ( veto_mu.empty() && recon_e.empty() ) ) {
159        MSG_DEBUG("Charged leptons left after selection");
160        vetoEvent;
161      }
162
163      if ( eTmiss <= 130 * GeV ) {
164        MSG_DEBUG("Not enough eTmiss: " << eTmiss << " < 130");
165        vetoEvent;
166      }
167
168      if ( recon_jets.empty() || recon_jets[0].pT() <= 130.0 * GeV ) {
169        MSG_DEBUG("No hard leading jet in " << recon_jets.size() << " jets");
170        vetoEvent;
171      }
172
173      // ==================== observables ====================
174
175      int Njets = 0;
176      double min_dPhi = 999.999;
177      double pTmiss_phi = pTmiss.phi();
178      for ( const Jet& jet : recon_jets ) {
179        if ( jet.pT() > 40 * GeV ) {
180          if ( Njets < 3 ) {
181            min_dPhi = min( min_dPhi, deltaPhi( pTmiss_phi, jet.phi() ) );
182          }
183          ++Njets;
184        }
185      }
186
187      int NjetsHighMass = 0;
188      for ( const Jet& jet : recon_jets ) {
189        if ( jet.pT() > 80.0 * GeV ) {
190          ++NjetsHighMass;
191        }
192      }
193
194      if ( Njets < 2 ) {
195        MSG_DEBUG("Only " << Njets << " >40 GeV jets left");
196        vetoEvent;
197      }
198
199      if ( min_dPhi <= 0.4 ) {
200        MSG_DEBUG("dPhi too small");
201        vetoEvent;
202      }
203
204      // m_eff
205      double m_eff_2j = eTmiss + recon_jets[0].pT() + recon_jets[1].pT();
206      double m_eff_3j = recon_jets.size() < 3 ? -999.0 : m_eff_2j + recon_jets[2].pT();
207      double m_eff_4j = recon_jets.size() < 4 ? -999.0 : m_eff_3j + recon_jets[3].pT();
208      double m_eff_HM = eTmiss;
209      for ( const Jet& jet : recon_jets ) {
210        if ( jet.pT() > 40.0 * GeV ) m_eff_HM += jet.pT();
211      }
212
213      double et_meff_2j = eTmiss / m_eff_2j;
214      double et_meff_3j = eTmiss / m_eff_3j;
215      double et_meff_4j = eTmiss / m_eff_4j;
216      double et_meff_HM = eTmiss / m_eff_HM;
217
218
219      // ==================== FILL ====================
220
221      MSG_DEBUG( "Trying to fill "
222                 << Njets << ' '
223                 << m_eff_2j << ' '
224                 << et_meff_2j << ' '
225                 << m_eff_3j << ' '
226                 << et_meff_3j << ' '
227                 << m_eff_4j << ' '
228                 << et_meff_4j << ' '
229                 << m_eff_HM << ' '
230                 << et_meff_HM );
231
232
233      _hist_eTmiss->fill(eTmiss);
234
235
236      // 2j region
237      if ( et_meff_2j > 0.3 ) {
238        _hist_meff_2j->fill(m_eff_2j);
239        if ( m_eff_2j > 1000 * GeV ) {
240          MSG_DEBUG("Hits 2j");
241          _count_2j->fill(0.5);
242        }
243      }
244
245
246      // 3j region
247      if ( Njets >= 3 && et_meff_3j > 0.25 ) {
248        _hist_meff_3j->fill(m_eff_3j);
249        if ( m_eff_3j > 1000 * GeV ) {
250          MSG_DEBUG("Hits 3j");
251          _count_3j->fill(0.5);
252        }
253      }
254
255
256      // 4j5 & 4j10 regions
257      if ( Njets >= 4 && et_meff_4j > 0.25 ) {
258        _hist_meff_4j->fill(m_eff_4j);
259        if ( m_eff_4j > 500 * GeV ) {
260          MSG_DEBUG("Hits 4j5");
261          _count_4j5->fill(0.5);
262        }
263        if ( m_eff_4j > 1000 * GeV ) {
264          MSG_DEBUG("Hits 4j10");
265          _count_4j10->fill(0.5);
266        }
267      }
268
269
270      // High mass region
271      if ( NjetsHighMass >= 4 && et_meff_HM > 0.2 ) {
272        _hist_meff_HM->fill(m_eff_HM);
273        if ( m_eff_HM > 1100 * GeV ) {
274          MSG_DEBUG("Hits HM");
275          _count_HM->fill(0.5);
276        }
277      }
278
279    }
280
281
282    void finalize() {
283      // Two, three and four jet channels have bin width = 100 (GeV)
284      // High mass channel has bin width = 150 (GeV)
285      // Integrated luminosity = 1040 (pb)
286      scale( _hist_meff_2j, 100. * 1040 * crossSection()/picobarn/sumOfWeights() );
287      scale( _hist_meff_3j, 100. * 1040 * crossSection()/picobarn/sumOfWeights() );
288      scale( _hist_meff_4j, 100. * 1040 * crossSection()/picobarn/sumOfWeights() );
289      scale( _hist_meff_HM, 150. * 1040 * crossSection()/picobarn/sumOfWeights() );
290    }
291
292    /// @}
293
294
295  private:
296
297    Histo1DPtr _count_2j;
298    Histo1DPtr _count_3j;
299    Histo1DPtr _count_4j5;
300    Histo1DPtr _count_4j10;
301    Histo1DPtr _count_HM;
302    Histo1DPtr _hist_meff_2j;
303    Histo1DPtr _hist_meff_3j;
304    Histo1DPtr _hist_meff_4j;
305    Histo1DPtr _hist_meff_HM;
306    Histo1DPtr _hist_eTmiss;
307
308  };
309
310
311  RIVET_DECLARE_ALIASED_PLUGIN(ATLAS_2011_I930002, ATLAS_2011_S9212183);
312
313}