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/VetoedFinalState.hh"
  8#include "Rivet/Projections/FastJets.hh"
  9
 10namespace Rivet {
 11
 12
 13  /// Single-lepton search for supersymmetry
 14  class ATLAS_2011_I930005 : public Analysis {
 15  public:
 16
 17    /// Constructor
 18    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2011_I930005);
 19
 20
 21    /// @name Analysis methods
 22    /// @{
 23
 24    /// Book histograms and initialize projections before the run
 25    void init() {
 26
 27      // projection to find the electrons
 28      IdentifiedFinalState elecs(Cuts::abseta < 2.47 && Cuts::pT > 20*GeV);
 29      elecs.acceptIdPair(PID::ELECTRON);
 30      declare(elecs, "elecs");
 31
 32
 33      // veto region electrons (from 2010 arXiv:1102.2357v2)
 34      Cut vetocut = Cuts::absetaIn(1.37, 1.52);
 35      IdentifiedFinalState veto_elecs(vetocut && Cuts::pT > 10*GeV);
 36      veto_elecs.acceptIdPair(PID::ELECTRON);
 37      declare(veto_elecs, "veto_elecs");
 38
 39
 40      // projection to find the muons
 41      IdentifiedFinalState muons(Cuts::abseta < 2.4 && Cuts::pT > 10*GeV);
 42      muons.acceptIdPair(PID::MUON);
 43      declare(muons, "muons");
 44
 45
 46      // Jet finder
 47      VetoedFinalState vfs;
 48      vfs.addVetoPairId(PID::MUON);
 49      declare(FastJets(vfs, JetAlg::ANTIKT, 0.4), "AntiKtJets04");
 50
 51
 52      // all tracks (to do deltaR with leptons)
 53      declare(ChargedFinalState(Cuts::abseta < 3 && Cuts::pT > 0.5*GeV), "cfs");
 54
 55
 56      // for pTmiss
 57      declare(VisibleFinalState(Cuts::abseta < 4.5),"vfs");
 58
 59
 60      /// Book histograms
 61      book(_3jl_count_mu_channel ,"3jl_count_muon_channel", 1, 0., 1.);
 62      book(_3jl_count_e_channel ,"3jl_count_electron_channel", 1, 0., 1.);
 63      book(_3jt_count_mu_channel ,"3jt_count_muon_channel", 1, 0., 1.);
 64      book(_3jt_count_e_channel ,"3jt_count_electron_channel", 1, 0., 1.);
 65      book(_3j_hist_eTmiss_e ,"3j_Et_miss_e", 65, 0., 650.);
 66      book(_3j_hist_eTmiss_mu ,"3j_Et_miss_mu", 65, 0., 650.);
 67      book(_3j_hist_mT_e ,"3j_mT_e", 58, 0., 580.);
 68      book(_3j_hist_mT_mu ,"3j_mT_mu", 58, 0., 580.);
 69      book(_3j_hist_m_eff_e ,"3j_m_eff_e", 46, 0., 2300.);
 70      book(_3j_hist_m_eff_mu ,"3j_m_eff_mu", 46, 0., 2300.);
 71      book(_3jl_hist_m_eff_e_final ,"3jl_m_eff_e_final", 15, 0., 1500.);
 72      book(_3jl_hist_m_eff_mu_final ,"3jl_m_eff_mu_final", 15, 0., 1500.);
 73      book(_3jt_hist_m_eff_e_final ,"3jt_m_eff_e_final", 15, 0., 1500.);
 74      book(_3jt_hist_m_eff_mu_final ,"3jt_m_eff_mu_final", 15, 0., 1500.);
 75
 76
 77      book(_4jl_count_mu_channel ,"4jl_count_muon_channel", 1, 0., 1.);
 78      book(_4jl_count_e_channel ,"4jl_count_electron_channel", 1, 0., 1.);
 79      book(_4jt_count_mu_channel ,"4jt_count_muon_channel", 1, 0., 1.);
 80      book(_4jt_count_e_channel ,"4jt_count_electron_channel", 1, 0., 1.);
 81      book(_4j_hist_eTmiss_e ,"4j_Et_miss_e", 65, 0., 650.);
 82      book(_4j_hist_eTmiss_mu ,"4j_Et_miss_mu", 65, 0., 650.);
 83      book(_4j_hist_mT_e ,"4j_mT_e", 58, 0., 580.);
 84      book(_4j_hist_mT_mu ,"4j_mT_mu", 58, 0., 580.);
 85      book(_4j_hist_m_eff_e ,"4j_m_eff_e", 46, 0., 2300.);
 86      book(_4j_hist_m_eff_mu ,"4j_m_eff_mu", 46, 0., 2300.);
 87      book(_4jl_hist_m_eff_e_final ,"4jl_m_eff_e_final", 15, 0., 1500.);
 88      book(_4jl_hist_m_eff_mu_final ,"4jl_m_eff_mu_final", 15, 0., 1500.);
 89      book(_4jt_hist_m_eff_e_final ,"4jt_m_eff_e_final", 15, 0., 1500.);
 90      book(_4jt_hist_m_eff_mu_final ,"4jt_m_eff_mu_final", 15, 0., 1500.);
 91
 92
 93    }
 94
 95
 96
 97    /// Perform the per-event analysis
 98    void analyze(const Event& event) {
 99      Particles veto_e
100        = apply<IdentifiedFinalState>(event, "veto_elecs").particles();
101      if ( ! veto_e.empty() ) {
102        MSG_DEBUG("electrons in veto region");
103        vetoEvent;
104      }
105
106      Jets cand_jets = apply<FastJets>(event, "AntiKtJets04").jetsByPt(Cuts::pT > 20*GeV && Cuts::abseta < 2.8);
107
108      Particles candtemp_e =
109        apply<IdentifiedFinalState>(event, "elecs").particlesByPt();
110      Particles candtemp_mu =
111        apply<IdentifiedFinalState>(event,"muons").particlesByPt();
112      Particles chg_tracks =
113        apply<ChargedFinalState>(event, "cfs").particles();
114      Particles cand_mu;
115      Particles cand_e;
116
117
118      // pTcone around muon track
119      for ( const Particle & mu : candtemp_mu ) {
120        double pTinCone = -mu.pT();
121        for ( const Particle & track : chg_tracks ) {
122          if ( deltaR(mu.momentum(),track.momentum()) < 0.2 )
123            pTinCone += track.pT();
124        }
125        if ( pTinCone < 1.8*GeV )
126          cand_mu.push_back(mu);
127      }
128
129      // pTcone around electron
130      for ( const Particle & e : candtemp_e ) {
131        double pTinCone = -e.pT();
132        for ( const Particle & track : chg_tracks ) {
133          if ( deltaR(e.momentum(),track.momentum()) < 0.2 )
134            pTinCone += track.pT();
135        }
136        if ( pTinCone < 0.1 * e.pT() )
137          cand_e.push_back(e);
138      }
139
140      // discard jets that overlap with electrons
141      Jets recon_jets;
142      for ( const Jet& jet : cand_jets ) {
143          bool away_from_e = true;
144          for ( const Particle & e : cand_e ) {
145            if ( deltaR(e.momentum(),jet.momentum()) < 0.2 ) {
146              away_from_e = false;
147              break;
148            }
149          }
150          if ( away_from_e )
151            recon_jets.push_back( jet );
152      }
153
154      // only consider leptons far from jet
155      Particles recon_e, recon_mu;
156      for ( const Particle & e : cand_e ) {
157        bool e_near_jet = false;
158        for ( const Jet& jet : recon_jets ) {
159          if ( deltaR(e.momentum(),jet.momentum()) < 0.4 &&
160               deltaR(e.momentum(),jet.momentum()) > 0.2 )
161            e_near_jet = true;
162        }
163        if ( ! e_near_jet )
164          recon_e.push_back( e );
165      }
166
167      for ( const Particle & mu : cand_mu ) {
168        bool mu_near_jet = false;
169        for ( const Jet& jet : recon_jets ) {
170          if ( deltaR(mu.momentum(),jet.momentum()) < 0.4 )
171            mu_near_jet = true;
172        }
173        if ( ! mu_near_jet )
174          recon_mu.push_back( mu );
175      }
176
177      // pTmiss
178      Particles vfs_particles
179        = apply<VisibleFinalState>(event, "vfs").particles();
180      FourMomentum pTmiss;
181      for ( const Particle & p : vfs_particles ) {
182        pTmiss -= p.momentum();
183      }
184      double eTmiss = pTmiss.pT();
185
186
187      // ==================== observables ====================
188
189
190      // Njets
191      int Njets = 0;
192      double pTmiss_phi = pTmiss.phi();
193      for ( const Jet& jet : recon_jets ) {
194        if ( jet.abseta() < 2.8 )
195          Njets+=1;
196      }
197      if ( Njets < 3 ) {
198        MSG_DEBUG("Only " << Njets << " jets w/ eta<2.8 left");
199        vetoEvent;
200      }
201
202      Particles lepton;
203      if ( recon_mu.empty() && recon_e.empty() ) {
204        MSG_DEBUG("No leptons");
205        vetoEvent;
206      }
207      else {
208        for ( const Particle & mu : recon_mu )
209            lepton.push_back(mu);
210        for ( const Particle & e : recon_e )
211            lepton.push_back(e);
212      }
213
214
215      std::sort(lepton.begin(), lepton.end(), cmpMomByPt);
216
217      double e_id = 11;
218      double mu_id = 13;
219
220      // one hard leading lepton cut
221      if ( lepton[0].abspid() == e_id &&
222           lepton[0].pT() <= 25*GeV ) {
223        vetoEvent;
224      }
225      else if ( lepton[0].abspid() == mu_id &&
226                lepton[0].pT() <= 20*GeV ) {
227        vetoEvent;
228      }
229
230      // exactly one hard leading lepton cut
231      if(lepton.size()>1) {
232        if ( lepton[1].abspid() == e_id &&
233             lepton[1].pT() > 20*GeV ) {
234          vetoEvent;
235        }
236        else if ( lepton[1].abspid() == mu_id &&
237                  lepton[1].pT() > 10*GeV ) {
238          vetoEvent;
239        }
240      }
241
242      // 3JL
243      if ( recon_jets[0].pT() > 60.0*GeV &&
244           recon_jets[1].pT() > 25.0*GeV &&
245           recon_jets[2].pT() > 25.0*GeV &&
246           deltaPhi( pTmiss_phi, recon_jets[0].phi() ) > 0.2 &&
247           deltaPhi( pTmiss_phi, recon_jets[1].phi() ) > 0.2 &&
248           deltaPhi( pTmiss_phi, recon_jets[2].phi() ) > 0.2 ) {
249
250        FourMomentum pT_l = lepton[0].momentum();
251        double dPhi = deltaPhi( pT_l.phi(), pTmiss_phi);
252        double mT = sqrt( 2 * pT_l.pT() * eTmiss * (1 - cos(dPhi)) );
253        double m_eff = eTmiss + pT_l.pT()
254          + recon_jets[0].pT()
255          + recon_jets[1].pT()
256          + recon_jets[2].pT();
257
258        if (  lepton[0].abspid() == e_id ) {
259          _3j_hist_mT_e->fill(mT);
260          _3j_hist_eTmiss_e->fill(eTmiss);
261          _3j_hist_m_eff_e->fill(m_eff);
262          if ( mT > 100*GeV && eTmiss > 125*GeV ) {
263            _3jl_hist_m_eff_e_final->fill(m_eff);
264            if ( m_eff > 500*GeV && eTmiss > 0.25*m_eff ) {
265              _3jl_count_e_channel->fill(0.5);
266            }
267          }
268        }
269
270        else if (  lepton[0].abspid() == mu_id ) {
271          _3j_hist_mT_mu->fill(mT);
272          _3j_hist_eTmiss_mu->fill(eTmiss);
273          _3j_hist_m_eff_mu->fill(m_eff);
274          if ( mT > 100*GeV && eTmiss > 125*GeV ) {
275            _3jl_hist_m_eff_mu_final->fill(m_eff);
276            if ( m_eff > 500*GeV && eTmiss > 0.25*m_eff ) {
277              _3jl_count_mu_channel->fill(0.5);
278            }
279          }
280        }
281
282      }
283
284      // 3JT
285      if ( recon_jets[0].pT() > 80.0*GeV &&
286           recon_jets[1].pT() > 25.0*GeV &&
287           recon_jets[2].pT() > 25.0*GeV &&
288           deltaPhi( pTmiss_phi, recon_jets[0].phi() ) > 0.2 &&
289           deltaPhi( pTmiss_phi, recon_jets[1].phi() ) > 0.2 &&
290           deltaPhi( pTmiss_phi, recon_jets[2].phi() ) > 0.2 ) {
291
292        FourMomentum pT_l = lepton[0].momentum();
293        double dPhi = deltaPhi( pT_l.phi(), pTmiss_phi);
294        double mT = sqrt( 2 * pT_l.pT() * eTmiss * (1 - cos(dPhi)) );
295        double m_eff = eTmiss + pT_l.pT()
296          + recon_jets[0].pT()
297          + recon_jets[1].pT()
298          + recon_jets[2].pT();
299
300
301        if (  lepton[0].abspid() == e_id ) {
302          if ( mT > 100*GeV && eTmiss > 240*GeV ) {
303            _3jt_hist_m_eff_e_final->fill(m_eff);
304            if ( m_eff > 600*GeV && eTmiss > 0.15*m_eff ) {
305              _3jt_count_e_channel->fill(0.5);
306            }
307          }
308        }
309
310        else if (  lepton[0].abspid() == mu_id ) {
311          if ( mT > 100*GeV && eTmiss > 240*GeV ) {
312            _3jt_hist_m_eff_mu_final->fill(m_eff);
313            if ( m_eff > 600*GeV && eTmiss > 0.15*m_eff ) {
314              _3jt_count_mu_channel->fill(0.5);
315            }
316          }
317        }
318
319      }
320
321      if ( Njets < 4 ) {
322        MSG_DEBUG("Only " << Njets << " jets w/ eta<2.8 left");
323        vetoEvent;
324      }
325
326
327
328      // 4JL
329      if ( recon_jets[0].pT() > 60.0*GeV &&
330           recon_jets[1].pT() > 25.0*GeV &&
331           recon_jets[2].pT() > 25.0*GeV &&
332           recon_jets[3].pT() > 25.0*GeV &&
333           deltaPhi( pTmiss_phi, recon_jets[0].phi() ) > 0.2 &&
334           deltaPhi( pTmiss_phi, recon_jets[1].phi() ) > 0.2 &&
335           deltaPhi( pTmiss_phi, recon_jets[2].phi() ) > 0.2 &&
336           deltaPhi( pTmiss_phi, recon_jets[3].phi() ) > 0.2 ) {
337
338        FourMomentum pT_l = lepton[0].momentum();
339        double dPhi = deltaPhi( pT_l.phi(), pTmiss_phi);
340        double mT = sqrt( 2 * pT_l.pT() * eTmiss * (1 - cos(dPhi)) );
341        double m_eff = eTmiss + pT_l.pT()
342          + recon_jets[0].pT()
343          + recon_jets[1].pT()
344          + recon_jets[2].pT()
345          + recon_jets[3].pT();
346
347
348        if (  lepton[0].abspid() == e_id ) {
349          _4j_hist_mT_e->fill(mT);
350          _4j_hist_eTmiss_e->fill(eTmiss);
351          _4j_hist_m_eff_e->fill(m_eff);
352          if ( mT > 100*GeV && eTmiss > 140*GeV ) {
353            _4jl_hist_m_eff_e_final->fill(m_eff);
354            if ( m_eff > 300*GeV && eTmiss > 0.3*m_eff ) {
355              _4jl_count_e_channel->fill(0.5);
356            }
357          }
358        }
359
360        // Muon channel signal region
361        else if (  lepton[0].abspid() == mu_id ) {
362          _4j_hist_mT_mu->fill(mT);
363          _4j_hist_eTmiss_mu->fill(eTmiss);
364          _4j_hist_m_eff_mu->fill(m_eff);
365          if ( mT > 100*GeV && eTmiss > 140*GeV ) {
366            _4jl_hist_m_eff_mu_final->fill(m_eff);
367            if ( m_eff > 300*GeV && eTmiss > 0.3*m_eff ) {
368              _4jl_count_mu_channel->fill(0.5);
369            }
370          }
371        }
372
373      }
374
375      // 4JT
376      if ( recon_jets[0].pT() > 60.0*GeV &&
377           recon_jets[1].pT() > 40.0*GeV &&
378           recon_jets[2].pT() > 40.0*GeV &&
379           recon_jets[3].pT() > 40.0*GeV &&
380           deltaPhi( pTmiss_phi, recon_jets[0].phi() ) > 0.2 &&
381           deltaPhi( pTmiss_phi, recon_jets[1].phi() ) > 0.2 &&
382           deltaPhi( pTmiss_phi, recon_jets[2].phi() ) > 0.2 &&
383           deltaPhi( pTmiss_phi, recon_jets[3].phi() ) > 0.2 ) {
384
385        FourMomentum pT_l = lepton[0].momentum();
386
387        double m_eff = eTmiss + pT_l.pT()
388          + recon_jets[0].pT()
389          + recon_jets[1].pT()
390          + recon_jets[2].pT()
391          + recon_jets[3].pT();
392
393
394        if (  lepton[0].abspid() == e_id ) {
395          if ( eTmiss > 200*GeV ) {
396            _4jt_hist_m_eff_e_final->fill(m_eff);
397            if ( m_eff > 500*GeV && eTmiss > 0.15*m_eff ) {
398              _4jt_count_e_channel->fill(0.5);
399            }
400          }
401        }
402
403        // Muon channel signal region
404        else if (  lepton[0].abspid() == mu_id ) {
405          if ( eTmiss > 200*GeV ) {
406            _4jt_hist_m_eff_mu_final->fill(m_eff);
407            if ( m_eff > 500*GeV && eTmiss > 0.15*m_eff ) {
408              _4jt_count_mu_channel->fill(0.5);
409            }
410          }
411        }
412
413      }
414     }
415
416    /// @}
417
418
419    void finalize() {
420
421      scale( _3j_hist_eTmiss_e, 10. * 1.04e3 * crossSection()/picobarn/sumOfWeights() );
422      scale( _3j_hist_eTmiss_mu, 10. * 1.04e3 * crossSection()/picobarn/sumOfWeights() );
423      scale( _3j_hist_m_eff_e, 50. * 1.04e3 * crossSection()/picobarn/sumOfWeights() );
424      scale( _3j_hist_m_eff_mu, 50. * 1.04e3 * crossSection()/picobarn/sumOfWeights() );
425      scale( _3j_hist_mT_e, 10. * 1.04e3 * crossSection()/picobarn/sumOfWeights() );
426      scale( _3j_hist_mT_mu, 10. * 1.04e3 * crossSection()/picobarn/sumOfWeights() );
427      scale( _3jl_hist_m_eff_e_final, 100. * 1.04e3 * crossSection()/picobarn/sumOfWeights() );
428      scale( _3jl_hist_m_eff_mu_final, 100. * 1.04e3 * crossSection()/picobarn/sumOfWeights() );
429      scale( _3jt_hist_m_eff_e_final, 100. * 1.04e3 * crossSection()/picobarn/sumOfWeights() );
430      scale( _3jt_hist_m_eff_mu_final, 100. * 1.04e3 * crossSection()/picobarn/sumOfWeights() );
431
432      scale( _4j_hist_eTmiss_e, 10. * 1.04e3 * crossSection()/picobarn/sumOfWeights() );
433      scale( _4j_hist_eTmiss_mu, 10. * 1.04e3 * crossSection()/picobarn/sumOfWeights() );
434      scale( _4j_hist_m_eff_e, 50. * 1.04e3 * crossSection()/picobarn/sumOfWeights() );
435      scale( _4j_hist_m_eff_mu, 50. * 1.04e3 * crossSection()/picobarn/sumOfWeights() );
436      scale( _4j_hist_mT_e, 10. * 1.04e3 * crossSection()/picobarn/sumOfWeights() );
437      scale( _4j_hist_mT_mu, 10. * 1.04e3 * crossSection()/picobarn/sumOfWeights() );
438      scale( _4jl_hist_m_eff_e_final, 100. * 1.04e3 * crossSection()/picobarn/sumOfWeights() );
439      scale( _4jl_hist_m_eff_mu_final, 100. * 1.04e3 * crossSection()/picobarn/sumOfWeights() );
440      scale( _4jt_hist_m_eff_e_final, 100. * 1.04e3 * crossSection()/picobarn/sumOfWeights() );
441      scale( _4jt_hist_m_eff_mu_final, 100. * 1.04e3 * crossSection()/picobarn/sumOfWeights() );
442
443
444    }
445
446  private:
447
448    /// @name Histograms
449    /// @{
450    Histo1DPtr _3jl_count_e_channel;
451    Histo1DPtr _3jl_count_mu_channel;
452    Histo1DPtr _3jt_count_e_channel;
453    Histo1DPtr _3jt_count_mu_channel;
454    Histo1DPtr _3j_hist_eTmiss_e;
455    Histo1DPtr _3j_hist_eTmiss_mu;
456    Histo1DPtr _3j_hist_m_eff_e;
457    Histo1DPtr _3j_hist_m_eff_mu;
458    Histo1DPtr _3j_hist_mT_e;
459    Histo1DPtr _3j_hist_mT_mu;
460    Histo1DPtr _3jl_hist_m_eff_e_final;
461    Histo1DPtr _3jl_hist_m_eff_mu_final;
462    Histo1DPtr _3jt_hist_m_eff_e_final;
463    Histo1DPtr _3jt_hist_m_eff_mu_final;
464
465
466
467    Histo1DPtr _4jl_count_e_channel;
468    Histo1DPtr _4jl_count_mu_channel;
469    Histo1DPtr _4jt_count_e_channel;
470    Histo1DPtr _4jt_count_mu_channel;
471    Histo1DPtr _4j_hist_eTmiss_e;
472    Histo1DPtr _4j_hist_eTmiss_mu;
473    Histo1DPtr _4j_hist_m_eff_e;
474    Histo1DPtr _4j_hist_m_eff_mu;
475    Histo1DPtr _4j_hist_mT_e;
476    Histo1DPtr _4j_hist_mT_mu;
477    Histo1DPtr _4jl_hist_m_eff_e_final;
478    Histo1DPtr _4jl_hist_m_eff_mu_final;
479    Histo1DPtr _4jt_hist_m_eff_e_final;
480    Histo1DPtr _4jt_hist_m_eff_mu_final;
481    /// @}
482
483  };
484
485
486
487  RIVET_DECLARE_ALIASED_PLUGIN(ATLAS_2011_I930005, ATLAS_2011_S9212353);
488
489}