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  class ATLAS_2011_CONF_2011_090 : public Analysis {
 14  public:
 15
 16    /// @name Constructors etc.
 17    //@{
 18
 19    /// Constructor
 20
 21    ATLAS_2011_CONF_2011_090()
 22      : Analysis("ATLAS_2011_CONF_2011_090")
 23    {    }
 24
 25    //@}
 26
 27
 28  public:
 29
 30    /// @name Analysis methods
 31    //@{
 32
 33    /// Book histograms and initialize projections before the run
 34    void init() {
 35
 36      // projection to find the electrons
 37      IdentifiedFinalState elecs(Cuts::abseta < 2.47 && Cuts::pT >= 20*GeV);
 38      elecs.acceptIdPair(PID::ELECTRON);
 39      declare(elecs, "elecs");
 40
 41      // veto region electrons (from 2010 arXiv:1102.2357v2)
 42      Cut vetocut = Cuts::absetaIn(1.37, 1.52);
 43      IdentifiedFinalState veto_elecs(vetocut && Cuts::pT > 10*GeV);
 44      veto_elecs.acceptIdPair(PID::ELECTRON);
 45      declare(veto_elecs, "veto_elecs");
 46
 47      // projection to find the muons
 48      IdentifiedFinalState muons(Cuts::abseta < 2.4 && Cuts::pT > 10*GeV);
 49      muons.acceptIdPair(PID::MUON);
 50      declare(muons, "muons");
 51
 52      // Jet finder
 53      VetoedFinalState vfs;
 54      vfs.addVetoPairId(PID::MUON);
 55      declare(FastJets(vfs, FastJets::ANTIKT, 0.4), "AntiKtJets04");
 56
 57      // all tracks (to do deltaR with leptons)
 58      declare(ChargedFinalState(Cuts::abseta < 3.0 && Cuts::pT > 0.5*GeV), "cfs");
 59
 60      // for pTmiss
 61      declare(VisibleFinalState(Cuts::abseta < 4.9), "vfs");
 62
 63      /// Book histograms
 64      book(_count_mu_channel ,"count_muon_channel", 1, 0., 1.);
 65      book(_count_e_channel ,"count_electron_channel", 1, 0., 1.);
 66      book(_hist_eTmiss_e ,"Et_miss_e", 50, 0., 500.);
 67      book(_hist_eTmiss_mu ,"Et_miss_mu", 50, 0., 500.);
 68      book(_hist_m_eff_e ,"m_eff_e", 60, 0., 1500.);
 69      book(_hist_m_eff_mu ,"m_eff_mu", 60, 0., 1500.);
 70      book(_hist_m_eff_e_final ,"m_eff_e_final", 15, 0., 1500.);
 71      book(_hist_m_eff_mu_final ,"m_eff_mu_final", 15, 0., 1500.);
 72    }
 73
 74
 75
 76    /// Perform the per-event analysis
 77    void analyze(const Event& event) {
 78
 79      const double weight = 1.0;
 80
 81      Particles veto_e = apply<IdentifiedFinalState>(event, "veto_elecs").particles();
 82      if ( ! veto_e.empty() ) {
 83        MSG_DEBUG("electrons in veto region");
 84        vetoEvent;
 85      }
 86
 87      Jets cand_jets = apply<FastJets>(event, "AntiKtJets04").jetsByPt(Cuts::pT > 20*GeV && Cuts::abseta < 2.8);
 88
 89      Particles candtemp_e = apply<IdentifiedFinalState>(event, "elecs").particlesByPt();
 90      Particles candtemp_mu = apply<IdentifiedFinalState>(event,"muons").particlesByPt();
 91      Particles chg_tracks = apply<ChargedFinalState>(event, "cfs").particles();
 92      Particles cand_mu;
 93      Particles cand_e;
 94
 95      // pTcone around muon track
 96      for ( const Particle& mu : candtemp_mu ) {
 97        double pTinCone = -mu.pT();
 98        for ( const Particle& track : chg_tracks ) {
 99          if ( deltaR(mu, track) < 0.2 )
100            pTinCone += track.pT();
101        }
102        if ( pTinCone < 1.8*GeV )
103          cand_mu.push_back(mu);
104      }
105
106
107      // pTcone around electron
108      for ( const Particle& e : candtemp_e ) {
109        double pTinCone = -e.pT();
110        for ( const Particle& track : chg_tracks ) {
111          if ( deltaR(e, track) < 0.2 )
112            pTinCone += track.pT();
113        }
114        if ( pTinCone < 0.10 * e.pT() )
115          cand_e.push_back(e);
116      }
117
118
119      // discard jets that overlap with electrons
120      Jets cand_jets_2;
121      for ( const Jet& jet : cand_jets ) {
122        bool away_from_e = true;
123        for ( const Particle & e : cand_e ) {
124          if ( deltaR(e, jet) <= 0.2 ) {
125            away_from_e = false;
126            break;
127          }
128        }
129        if ( away_from_e )
130          cand_jets_2.push_back( jet );
131      }
132
133      // only consider leptons far from jet
134      Particles recon_e, recon_mu;
135      for ( const Particle & e : cand_e ) {
136        bool e_near_jet = false;
137        for ( const Jet& jet : cand_jets_2 ) {
138          if (inRange(deltaR(e, jet), 0.2, 0.4)) e_near_jet = true;
139        }
140        if ( ! e_near_jet ) recon_e.push_back( e );
141      }
142
143      for ( const Particle & mu : cand_mu ) {
144        bool mu_near_jet = false;
145        for ( const Jet& jet : cand_jets_2 ) {
146          if ( deltaR(mu, jet) < 0.4 ) mu_near_jet = true;
147        }
148        if ( ! mu_near_jet ) recon_mu.push_back( mu );
149      }
150
151      // pTmiss
152      Particles vfs_particles
153        = apply<VisibleFinalState>(event, "vfs").particles();
154      FourMomentum pTmiss;
155      for ( const Particle & p : vfs_particles ) {
156        pTmiss -= p.momentum();
157      }
158      double eTmiss = pTmiss.pT();
159
160
161      // final jet filter
162      Jets recon_jets;
163      for ( const Jet& jet : cand_jets_2 ) {
164        recon_jets.push_back( jet );
165      }
166
167
168
169
170      // ==================== observables ====================
171
172
173      // Njets
174
175      int Njets = 0;
176      double pTmiss_phi = pTmiss.phi();
177      for ( const Jet& jet : recon_jets ) {
178        if ( jet.abseta() < 2.8 )
179          Njets+=1;
180      }
181      if ( Njets < 3 ) {
182        MSG_DEBUG("Only " << Njets << " jets w/ eta<2.8 left");
183        vetoEvent;
184      }
185
186      if ( recon_jets[0].pT() <= 60.0 * GeV ) {
187        MSG_DEBUG("No hard leading jet in " << recon_jets.size() << " jets");
188        vetoEvent;
189      }
190      for ( int i = 1; i < 3; ++i ) {
191        if ( recon_jets[i].pT() <= 25*GeV ) {
192          vetoEvent;
193        }
194      }
195
196      for ( int i = 0; i < 3; ++i ) {
197        double dPhi = deltaPhi( pTmiss_phi, recon_jets[i].phi() );
198        if ( dPhi <= 0.2 ) {
199          MSG_DEBUG("dPhi too small");
200          vetoEvent;
201          break;
202        }
203      }
204
205
206      Particles lepton;
207      if ( recon_mu.empty() && recon_e.empty() ) {
208        MSG_DEBUG("No leptons");
209        vetoEvent;
210      }
211      else {
212        for ( const Particle & mu : recon_mu )
213          lepton.push_back(mu);
214        for ( const Particle & e : recon_e )
215          lepton.push_back(e);
216      }
217
218      std::sort(lepton.begin(), lepton.end(), cmpMomByPt);
219
220      double e_id = 11;
221      double mu_id = 13;
222
223      // one hard leading lepton cut
224      if ( lepton[0].abspid() == e_id &&
225           lepton[0].pT() <= 25*GeV ) {
226        vetoEvent;
227      }
228      else if ( lepton[0].abspid() == mu_id &&
229                lepton[0].pT() <= 20*GeV ) {
230        vetoEvent;
231      }
232
233      // exactly one hard leading lepton cut
234      if(lepton.size()>1) {
235        if ( lepton[1].abspid() == e_id &&
236             lepton[1].pT() > 20*GeV ) {
237          vetoEvent;
238        }
239        else if ( lepton[1].abspid() == mu_id &&
240                  lepton[1].pT() > 10*GeV ) {
241          vetoEvent;
242        }
243      }
244
245
246      // ==================== FILL ====================
247
248
249      FourMomentum pT_l = lepton[0].momentum();
250
251
252      double dPhi = deltaPhi( pT_l.phi(), pTmiss_phi);
253      double mT = sqrt( 2 * pT_l.pT() * eTmiss * (1 - cos(dPhi)) );
254
255
256      // effective mass
257      double m_eff = eTmiss + pT_l.pT()
258        + recon_jets[0].pT()
259        + recon_jets[1].pT()
260        + recon_jets[2].pT();
261
262
263      // Electron channel signal region
264
265      if (  lepton[0].abspid() == e_id ) {
266
267        _hist_eTmiss_e->fill(eTmiss, weight);
268        _hist_m_eff_e->fill(m_eff, weight);
269
270        if ( mT > 100*GeV && eTmiss > 125*GeV ) {
271          _hist_m_eff_e_final->fill(m_eff, weight);
272          if ( m_eff > 500*GeV && eTmiss > 0.25*m_eff ) {
273            _count_e_channel->fill(0.5,weight);
274          }
275        }
276      }
277
278      // Muon channel signal region
279
280      else if (  lepton[0].abspid() == mu_id ) {
281
282        _hist_eTmiss_mu->fill(eTmiss, weight);
283        _hist_m_eff_mu->fill(m_eff, weight);
284
285        if ( mT > 100*GeV && eTmiss > 125*GeV ) {
286          _hist_m_eff_mu_final->fill(m_eff, weight);
287          if ( m_eff > 500*GeV && eTmiss > 0.25*m_eff ) {
288            _count_mu_channel->fill(0.5,weight);
289          }
290        }
291
292      }
293
294
295    }
296
297    //@}
298
299
300    void finalize() {
301      scale( _hist_eTmiss_e      , 10.  * 165. * crossSection()/picobarn/sumOfWeights() );
302      scale( _hist_eTmiss_mu     , 10.  * 165. * crossSection()/picobarn/sumOfWeights() );
303      scale( _hist_m_eff_e       , 25.  * 165. * crossSection()/picobarn/sumOfWeights() );
304      scale( _hist_m_eff_mu      , 25.  * 165. * crossSection()/picobarn/sumOfWeights() );
305      scale( _hist_m_eff_e_final , 100. * 165. * crossSection()/picobarn/sumOfWeights() );
306      scale( _hist_m_eff_mu_final, 100. * 165. * crossSection()/picobarn/sumOfWeights() );
307    }
308
309  private:
310
311    /// @name Histograms
312    //@{
313    Histo1DPtr _count_e_channel;
314    Histo1DPtr _count_mu_channel;
315
316    Histo1DPtr _hist_eTmiss_e;
317    Histo1DPtr _hist_eTmiss_mu;
318
319    Histo1DPtr _hist_m_eff_e;
320    Histo1DPtr _hist_m_eff_mu;
321    Histo1DPtr _hist_m_eff_e_final;
322    Histo1DPtr _hist_m_eff_mu_final;
323
324
325    //@}
326
327
328  };
329
330
331
332  // The hook for the plugin system
333  RIVET_DECLARE_PLUGIN(ATLAS_2011_CONF_2011_090);
334
335}