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_S9212353 : public Analysis {
 15  public:
 16
 17    /// Constructor
 18    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2011_S9212353);
 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, FastJets::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      const double weight = 1.0;
100      Particles veto_e
101        = apply<IdentifiedFinalState>(event, "veto_elecs").particles();
102      if ( ! veto_e.empty() ) {
103        MSG_DEBUG("electrons in veto region");
104        vetoEvent;
105      }
106
107      Jets cand_jets;
108      for ( const Jet& jet :
109          apply<FastJets>(event, "AntiKtJets04").jetsByPt(20.0*GeV) ) {
110        if ( fabs( jet.eta() ) < 2.8 ) {
111          cand_jets.push_back(jet);
112        }
113      }
114
115      Particles candtemp_e =
116        apply<IdentifiedFinalState>(event, "elecs").particlesByPt();
117      Particles candtemp_mu =
118        apply<IdentifiedFinalState>(event,"muons").particlesByPt();
119      Particles chg_tracks =
120        apply<ChargedFinalState>(event, "cfs").particles();
121      Particles cand_mu;
122      Particles cand_e;
123
124
125      // pTcone around muon track
126      for ( const Particle & mu : candtemp_mu ) {
127        double pTinCone = -mu.pT();
128        for ( const Particle & track : chg_tracks ) {
129          if ( deltaR(mu.momentum(),track.momentum()) < 0.2 )
130            pTinCone += track.pT();
131        }
132        if ( pTinCone < 1.8*GeV )
133          cand_mu.push_back(mu);
134      }
135
136      // pTcone around electron
137      for ( const Particle & e : candtemp_e ) {
138        double pTinCone = -e.pT();
139        for ( const Particle & track : chg_tracks ) {
140          if ( deltaR(e.momentum(),track.momentum()) < 0.2 )
141            pTinCone += track.pT();
142        }
143        if ( pTinCone < 0.1 * e.pT() )
144          cand_e.push_back(e);
145      }
146
147      // discard jets that overlap with electrons
148      Jets recon_jets;
149      for ( const Jet& jet : cand_jets ) {
150          bool away_from_e = true;
151          for ( const Particle & e : cand_e ) {
152            if ( deltaR(e.momentum(),jet.momentum()) < 0.2 ) {
153              away_from_e = false;
154              break;
155            }
156          }
157          if ( away_from_e )
158            recon_jets.push_back( jet );
159      }
160
161      // only consider leptons far from jet
162      Particles recon_e, recon_mu;
163      for ( const Particle & e : cand_e ) {
164        bool e_near_jet = false;
165        for ( const Jet& jet : recon_jets ) {
166          if ( deltaR(e.momentum(),jet.momentum()) < 0.4 &&
167               deltaR(e.momentum(),jet.momentum()) > 0.2 )
168            e_near_jet = true;
169        }
170        if ( ! e_near_jet )
171          recon_e.push_back( e );
172      }
173
174      for ( const Particle & mu : cand_mu ) {
175        bool mu_near_jet = false;
176        for ( const Jet& jet : recon_jets ) {
177          if ( deltaR(mu.momentum(),jet.momentum()) < 0.4 )
178            mu_near_jet = true;
179        }
180        if ( ! mu_near_jet )
181          recon_mu.push_back( mu );
182      }
183
184      // pTmiss
185      Particles vfs_particles
186        = apply<VisibleFinalState>(event, "vfs").particles();
187      FourMomentum pTmiss;
188      for ( const Particle & p : vfs_particles ) {
189        pTmiss -= p.momentum();
190      }
191      double eTmiss = pTmiss.pT();
192
193
194      // ==================== observables ====================
195
196
197      // Njets
198      int Njets = 0;
199      double pTmiss_phi = pTmiss.phi();
200      for ( const Jet& jet : recon_jets ) {
201        if ( jet.abseta() < 2.8 )
202          Njets+=1;
203      }
204      if ( Njets < 3 ) {
205        MSG_DEBUG("Only " << Njets << " jets w/ eta<2.8 left");
206        vetoEvent;
207      }
208
209      Particles lepton;
210      if ( recon_mu.empty() && recon_e.empty() ) {
211        MSG_DEBUG("No leptons");
212        vetoEvent;
213      }
214      else {
215        for ( const Particle & mu : recon_mu )
216            lepton.push_back(mu);
217        for ( const Particle & e : recon_e )
218            lepton.push_back(e);
219      }
220
221
222      std::sort(lepton.begin(), lepton.end(), cmpMomByPt);
223
224      double e_id = 11;
225      double mu_id = 13;
226
227      // one hard leading lepton cut
228      if ( lepton[0].abspid() == e_id &&
229           lepton[0].pT() <= 25*GeV ) {
230        vetoEvent;
231      }
232      else if ( lepton[0].abspid() == mu_id &&
233                lepton[0].pT() <= 20*GeV ) {
234        vetoEvent;
235      }
236
237      // exactly one hard leading lepton cut
238      if(lepton.size()>1) {
239        if ( lepton[1].abspid() == e_id &&
240             lepton[1].pT() > 20*GeV ) {
241          vetoEvent;
242        }
243        else if ( lepton[1].abspid() == mu_id &&
244                  lepton[1].pT() > 10*GeV ) {
245          vetoEvent;
246        }
247      }
248
249      // 3JL
250      if ( recon_jets[0].pT() > 60.0*GeV &&
251           recon_jets[1].pT() > 25.0*GeV &&
252           recon_jets[2].pT() > 25.0*GeV &&
253           deltaPhi( pTmiss_phi, recon_jets[0].phi() ) > 0.2 &&
254           deltaPhi( pTmiss_phi, recon_jets[1].phi() ) > 0.2 &&
255           deltaPhi( pTmiss_phi, recon_jets[2].phi() ) > 0.2 ) {
256
257        FourMomentum pT_l = lepton[0].momentum();
258        double dPhi = deltaPhi( pT_l.phi(), pTmiss_phi);
259        double mT = sqrt( 2 * pT_l.pT() * eTmiss * (1 - cos(dPhi)) );
260        double m_eff = eTmiss + pT_l.pT()
261          + recon_jets[0].pT()
262          + recon_jets[1].pT()
263          + recon_jets[2].pT();
264
265        if (  lepton[0].abspid() == e_id ) {
266          _3j_hist_mT_e->fill(mT, weight);
267          _3j_hist_eTmiss_e->fill(eTmiss, weight);
268          _3j_hist_m_eff_e->fill(m_eff, weight);
269          if ( mT > 100*GeV && eTmiss > 125*GeV ) {
270            _3jl_hist_m_eff_e_final->fill(m_eff, weight);
271            if ( m_eff > 500*GeV && eTmiss > 0.25*m_eff ) {
272              _3jl_count_e_channel->fill(0.5,weight);
273            }
274          }
275        }
276
277        else if (  lepton[0].abspid() == mu_id ) {
278          _3j_hist_mT_mu->fill(mT, weight);
279          _3j_hist_eTmiss_mu->fill(eTmiss, weight);
280          _3j_hist_m_eff_mu->fill(m_eff, weight);
281          if ( mT > 100*GeV && eTmiss > 125*GeV ) {
282            _3jl_hist_m_eff_mu_final->fill(m_eff, weight);
283            if ( m_eff > 500*GeV && eTmiss > 0.25*m_eff ) {
284              _3jl_count_mu_channel->fill(0.5,weight);
285            }
286          }
287        }
288
289      }
290
291      // 3JT
292      if ( recon_jets[0].pT() > 80.0*GeV &&
293           recon_jets[1].pT() > 25.0*GeV &&
294           recon_jets[2].pT() > 25.0*GeV &&
295           deltaPhi( pTmiss_phi, recon_jets[0].phi() ) > 0.2 &&
296           deltaPhi( pTmiss_phi, recon_jets[1].phi() ) > 0.2 &&
297           deltaPhi( pTmiss_phi, recon_jets[2].phi() ) > 0.2 ) {
298
299        FourMomentum pT_l = lepton[0].momentum();
300        double dPhi = deltaPhi( pT_l.phi(), pTmiss_phi);
301        double mT = sqrt( 2 * pT_l.pT() * eTmiss * (1 - cos(dPhi)) );
302        double m_eff = eTmiss + pT_l.pT()
303          + recon_jets[0].pT()
304          + recon_jets[1].pT()
305          + recon_jets[2].pT();
306
307
308        if (  lepton[0].abspid() == e_id ) {
309          if ( mT > 100*GeV && eTmiss > 240*GeV ) {
310            _3jt_hist_m_eff_e_final->fill(m_eff, weight);
311            if ( m_eff > 600*GeV && eTmiss > 0.15*m_eff ) {
312              _3jt_count_e_channel->fill(0.5,weight);
313            }
314          }
315        }
316
317        else if (  lepton[0].abspid() == mu_id ) {
318          if ( mT > 100*GeV && eTmiss > 240*GeV ) {
319            _3jt_hist_m_eff_mu_final->fill(m_eff, weight);
320            if ( m_eff > 600*GeV && eTmiss > 0.15*m_eff ) {
321              _3jt_count_mu_channel->fill(0.5,weight);
322            }
323          }
324        }
325
326      }
327
328      if ( Njets < 4 ) {
329        MSG_DEBUG("Only " << Njets << " jets w/ eta<2.8 left");
330        vetoEvent;
331      }
332
333
334
335      // 4JL
336      if ( recon_jets[0].pT() > 60.0*GeV &&
337           recon_jets[1].pT() > 25.0*GeV &&
338           recon_jets[2].pT() > 25.0*GeV &&
339           recon_jets[3].pT() > 25.0*GeV &&
340           deltaPhi( pTmiss_phi, recon_jets[0].phi() ) > 0.2 &&
341           deltaPhi( pTmiss_phi, recon_jets[1].phi() ) > 0.2 &&
342           deltaPhi( pTmiss_phi, recon_jets[2].phi() ) > 0.2 &&
343           deltaPhi( pTmiss_phi, recon_jets[3].phi() ) > 0.2 ) {
344
345        FourMomentum pT_l = lepton[0].momentum();
346        double dPhi = deltaPhi( pT_l.phi(), pTmiss_phi);
347        double mT = sqrt( 2 * pT_l.pT() * eTmiss * (1 - cos(dPhi)) );
348        double m_eff = eTmiss + pT_l.pT()
349          + recon_jets[0].pT()
350          + recon_jets[1].pT()
351          + recon_jets[2].pT()
352          + recon_jets[3].pT();
353
354
355        if (  lepton[0].abspid() == e_id ) {
356          _4j_hist_mT_e->fill(mT, weight);
357          _4j_hist_eTmiss_e->fill(eTmiss, weight);
358          _4j_hist_m_eff_e->fill(m_eff, weight);
359          if ( mT > 100*GeV && eTmiss > 140*GeV ) {
360            _4jl_hist_m_eff_e_final->fill(m_eff, weight);
361            if ( m_eff > 300*GeV && eTmiss > 0.3*m_eff ) {
362              _4jl_count_e_channel->fill(0.5,weight);
363            }
364          }
365        }
366
367        // Muon channel signal region
368        else if (  lepton[0].abspid() == mu_id ) {
369          _4j_hist_mT_mu->fill(mT, weight);
370          _4j_hist_eTmiss_mu->fill(eTmiss, weight);
371          _4j_hist_m_eff_mu->fill(m_eff, weight);
372          if ( mT > 100*GeV && eTmiss > 140*GeV ) {
373            _4jl_hist_m_eff_mu_final->fill(m_eff, weight);
374            if ( m_eff > 300*GeV && eTmiss > 0.3*m_eff ) {
375              _4jl_count_mu_channel->fill(0.5,weight);
376            }
377          }
378        }
379
380      }
381
382      // 4JT
383      if ( recon_jets[0].pT() > 60.0*GeV &&
384           recon_jets[1].pT() > 40.0*GeV &&
385           recon_jets[2].pT() > 40.0*GeV &&
386           recon_jets[3].pT() > 40.0*GeV &&
387           deltaPhi( pTmiss_phi, recon_jets[0].phi() ) > 0.2 &&
388           deltaPhi( pTmiss_phi, recon_jets[1].phi() ) > 0.2 &&
389           deltaPhi( pTmiss_phi, recon_jets[2].phi() ) > 0.2 &&
390           deltaPhi( pTmiss_phi, recon_jets[3].phi() ) > 0.2 ) {
391
392        FourMomentum pT_l = lepton[0].momentum();
393
394        double m_eff = eTmiss + pT_l.pT()
395          + recon_jets[0].pT()
396          + recon_jets[1].pT()
397          + recon_jets[2].pT()
398          + recon_jets[3].pT();
399
400
401        if (  lepton[0].abspid() == e_id ) {
402          if ( eTmiss > 200*GeV ) {
403            _4jt_hist_m_eff_e_final->fill(m_eff, weight);
404            if ( m_eff > 500*GeV && eTmiss > 0.15*m_eff ) {
405              _4jt_count_e_channel->fill(0.5,weight);
406            }
407          }
408        }
409
410        // Muon channel signal region
411        else if (  lepton[0].abspid() == mu_id ) {
412          if ( eTmiss > 200*GeV ) {
413            _4jt_hist_m_eff_mu_final->fill(m_eff, weight);
414            if ( m_eff > 500*GeV && eTmiss > 0.15*m_eff ) {
415              _4jt_count_mu_channel->fill(0.5,weight);
416            }
417          }
418        }
419
420      }
421     }
422
423    //@}
424
425
426    void finalize() {
427
428      scale( _3j_hist_eTmiss_e, 10. * 1.04e3 * crossSection()/sumOfWeights() );
429      scale( _3j_hist_eTmiss_mu, 10. * 1.04e3 * crossSection()/sumOfWeights() );
430      scale( _3j_hist_m_eff_e, 50. * 1.04e3 * crossSection()/sumOfWeights() );
431      scale( _3j_hist_m_eff_mu, 50. * 1.04e3 * crossSection()/sumOfWeights() );
432      scale( _3j_hist_mT_e, 10. * 1.04e3 * crossSection()/sumOfWeights() );
433      scale( _3j_hist_mT_mu, 10. * 1.04e3 * crossSection()/sumOfWeights() );
434      scale( _3jl_hist_m_eff_e_final, 100. * 1.04e3 * crossSection()/sumOfWeights() );
435      scale( _3jl_hist_m_eff_mu_final, 100. * 1.04e3 * crossSection()/sumOfWeights() );
436      scale( _3jt_hist_m_eff_e_final, 100. * 1.04e3 * crossSection()/sumOfWeights() );
437      scale( _3jt_hist_m_eff_mu_final, 100. * 1.04e3 * crossSection()/sumOfWeights() );
438
439      scale( _4j_hist_eTmiss_e, 10. * 1.04e3 * crossSection()/sumOfWeights() );
440      scale( _4j_hist_eTmiss_mu, 10. * 1.04e3 * crossSection()/sumOfWeights() );
441      scale( _4j_hist_m_eff_e, 50. * 1.04e3 * crossSection()/sumOfWeights() );
442      scale( _4j_hist_m_eff_mu, 50. * 1.04e3 * crossSection()/sumOfWeights() );
443      scale( _4j_hist_mT_e, 10. * 1.04e3 * crossSection()/sumOfWeights() );
444      scale( _4j_hist_mT_mu, 10. * 1.04e3 * crossSection()/sumOfWeights() );
445      scale( _4jl_hist_m_eff_e_final, 100. * 1.04e3 * crossSection()/sumOfWeights() );
446      scale( _4jl_hist_m_eff_mu_final, 100. * 1.04e3 * crossSection()/sumOfWeights() );
447      scale( _4jt_hist_m_eff_e_final, 100. * 1.04e3 * crossSection()/sumOfWeights() );
448      scale( _4jt_hist_m_eff_mu_final, 100. * 1.04e3 * crossSection()/sumOfWeights() );
449
450
451    }
452
453  private:
454
455    /// @name Histograms
456    //@{
457    Histo1DPtr _3jl_count_e_channel;
458    Histo1DPtr _3jl_count_mu_channel;
459    Histo1DPtr _3jt_count_e_channel;
460    Histo1DPtr _3jt_count_mu_channel;
461    Histo1DPtr _3j_hist_eTmiss_e;
462    Histo1DPtr _3j_hist_eTmiss_mu;
463    Histo1DPtr _3j_hist_m_eff_e;
464    Histo1DPtr _3j_hist_m_eff_mu;
465    Histo1DPtr _3j_hist_mT_e;
466    Histo1DPtr _3j_hist_mT_mu;
467    Histo1DPtr _3jl_hist_m_eff_e_final;
468    Histo1DPtr _3jl_hist_m_eff_mu_final;
469    Histo1DPtr _3jt_hist_m_eff_e_final;
470    Histo1DPtr _3jt_hist_m_eff_mu_final;
471
472
473
474    Histo1DPtr _4jl_count_e_channel;
475    Histo1DPtr _4jl_count_mu_channel;
476    Histo1DPtr _4jt_count_e_channel;
477    Histo1DPtr _4jt_count_mu_channel;
478    Histo1DPtr _4j_hist_eTmiss_e;
479    Histo1DPtr _4j_hist_eTmiss_mu;
480    Histo1DPtr _4j_hist_m_eff_e;
481    Histo1DPtr _4j_hist_m_eff_mu;
482    Histo1DPtr _4j_hist_mT_e;
483    Histo1DPtr _4j_hist_mT_mu;
484    Histo1DPtr _4jl_hist_m_eff_e_final;
485    Histo1DPtr _4jl_hist_m_eff_mu_final;
486    Histo1DPtr _4jt_hist_m_eff_e_final;
487    Histo1DPtr _4jt_hist_m_eff_mu_final;
488    //@}
489
490  };
491
492
493
494  RIVET_DECLARE_ALIASED_PLUGIN(ATLAS_2011_S9212353, ATLAS_2011_I930005);
495
496}