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_2012_CONF_2012_104 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    ATLAS_2012_CONF_2012_104()
 18      : Analysis("ATLAS_2012_CONF_2012_104")
 19    {    }
 20
 21
 22    /// @name Analysis methods
 23    //@{
 24
 25    /// Book histograms and initialize projections before the run
 26    void init() {
 27
 28      // projection to find the electrons
 29      IdentifiedFinalState elecs(Cuts::abseta < 2.47 && Cuts::pT > 10*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      VetoedFinalState vfs;
 40      vfs.addVetoPairId(PID::MUON);
 41      declare(FastJets(vfs, FastJets::ANTIKT, 0.4), "AntiKtJets04");
 42
 43      // all tracks (to do deltaR with leptons)
 44      declare(ChargedFinalState(Cuts::abseta < 3 && Cuts::pT > 0.5*GeV), "cfs");
 45
 46      // for pTmiss
 47      declare(VisibleFinalState(Cuts::abseta < 4.9),"vfs");
 48
 49      // Book histograms
 50      book(_count_e  ,"count_e" , 1, 0., 1.);
 51      book(_count_mu ,"count_mu", 1, 0., 1.);
 52
 53      book(_hist_eTmiss_e  ,"hist_eTmiss_e"  , 25, 0., 1000.);
 54      book(_hist_eTmiss_mu ,"hist_eTmiss_mu" , 25, 0., 1000.);
 55
 56    }
 57
 58    /// Perform the per-event analysis
 59    void analyze(const Event& event) {
 60      const double weight = 1.0;
 61
 62      // get the candiate jets
 63      Jets cand_jets;
 64      for ( const Jet& jet :
 65                apply<FastJets>(event, "AntiKtJets04").jetsByPt(20.0*GeV) ) {
 66        if ( fabs( jet.eta() ) < 2.8 ) {
 67          cand_jets.push_back(jet);
 68        }
 69      }
 70
 71      // get the candidate "medium" leptons without isolation
 72      Particles cand_e;
 73      for( const Particle & e :
 74               apply<IdentifiedFinalState>(event, "elecs").particlesByPt()) {
 75        // remove any leptons within 0.4 of any candidate jets
 76        bool e_near_jet = false;
 77        for ( const Jet& jet : cand_jets ) {
 78          double dR = deltaR(e.momentum(),jet.momentum());
 79          if ( dR < 0.4 && dR > 0.2 ) {
 80            e_near_jet = true;
 81            break;
 82          }
 83        }
 84        if ( ! e_near_jet ) cand_e.push_back(e);
 85      }
 86      Particles cand_mu;
 87      for( const Particle & mu :
 88               apply<IdentifiedFinalState>(event, "muons").particlesByPt()) {
 89        // remove any leptons within 0.4 of any candidate jets
 90        bool mu_near_jet = false;
 91        for ( const Jet& jet : cand_jets ) {
 92          if ( deltaR(mu.momentum(),jet.momentum()) < 0.4 ) {
 93            mu_near_jet = true;
 94            break;
 95          }
 96        }
 97        if ( ! mu_near_jet ) cand_mu.push_back(mu);
 98      }
 99      // apply the isolation
100      Particles chg_tracks =
101        apply<ChargedFinalState>(event, "cfs").particles();
102      // pTcone around muon track (hard)
103      Particles recon_mu;
104      for ( const Particle & mu : cand_mu ) {
105        double pTinCone = -mu.pT();
106        if(-pTinCone<25.) continue;
107        for ( const Particle & track : chg_tracks ) {
108          if ( deltaR(mu.momentum(),track.momentum()) < 0.2 )
109            pTinCone += track.pT();
110        }
111        if ( pTinCone < 1.8*GeV ) recon_mu.push_back(mu);
112      }
113      // pTcone around electron track (hard)
114      Particles recon_e;
115      for ( const Particle & e : cand_e ) {
116        double pTinCone = -e.pT();
117        if(-pTinCone<25.) continue;
118        for ( const Particle & track : chg_tracks ) {
119          if ( deltaR(e.momentum(),track.momentum()) < 0.2 )
120            pTinCone += track.pT();
121        }
122        if ( pTinCone < 0.1 * e.pT() ) recon_e.push_back(e);
123      }
124
125      // discard jets that overlap with electrons
126      Jets recon_jets;
127      for ( const Jet& jet : cand_jets ) {
128        if(jet.abseta()>2.5||
129           jet.perp()<25.) continue;
130        bool away_from_e = true;
131        for ( const Particle & e : cand_e ) {
132          if ( deltaR(e.momentum(),jet.momentum()) < 0.2 ) {
133            away_from_e = false;
134            break;
135          }
136        }
137        if ( away_from_e ) recon_jets.push_back( jet );
138      }
139
140      // pTmiss
141      FourMomentum pTmiss;
142      for ( const Particle & p :
143                apply<VisibleFinalState>(event, "vfs").particles() ) {
144        pTmiss -= p.momentum();
145      }
146      double eTmiss = pTmiss.pT();
147
148      // at least 4 jets with pT>80.
149      if(recon_jets.size()<4 || recon_jets[3].perp()<80.) vetoEvent;
150
151      // only 1 signal lepton
152      if( recon_e.size() + recon_mu.size() != 1 )
153        vetoEvent;
154      if( cand_e .size() + cand_mu .size() != 1 )
155        vetoEvent;
156
157      // start of meff calculation
158      double HT=0.;
159      for( const Jet & jet : recon_jets) {
160        double pT = jet.perp();
161        if(pT>40.) HT += pT;
162      }
163
164      // get the lepton
165      Particle lepton = recon_e.empty() ? recon_mu[0] : recon_e[0];
166
167      // lepton variables
168      double pT = lepton.perp();
169
170      double mT  = 2.*(pT*eTmiss -
171                       lepton.px()*pTmiss.px() -
172                       lepton.py()*pTmiss.py());
173      mT = sqrt(mT);
174      HT += pT;
175      double m_eff_inc  = HT + eTmiss + pT;
176      double m_eff_4 = eTmiss + pT;
177      for(unsigned int ix=0;ix<4;++ix)
178        m_eff_4 +=  recon_jets[ix].perp();
179
180      // four jet selecton
181      if(mT>100.&& eTmiss/m_eff_4>0.2 &&
182         m_eff_inc > 800.) {
183        if( eTmiss > 250. ) {
184          if(lepton.abspid()==PID::ELECTRON)
185            _count_e->fill(0.5,weight);
186          else if(lepton.abspid()==PID::MUON)
187            _count_mu->fill(0.5,weight);
188        }
189        if(lepton.abspid()==PID::ELECTRON)
190          _hist_eTmiss_e ->fill(eTmiss,weight);
191        else if(lepton.abspid()==PID::MUON)
192          _hist_eTmiss_mu->fill(eTmiss,weight);
193      }
194    }
195    //@}
196
197
198    void finalize() {
199
200      double norm = 5.8* crossSection()/sumOfWeights()/femtobarn;
201      scale(_count_e ,norm);
202      scale(_count_mu,norm);
203      scale(_hist_eTmiss_e  ,40.*norm);
204      scale(_hist_eTmiss_mu ,40.*norm);
205
206    }
207
208  private:
209
210    /// @name Histograms
211    //@{
212    Histo1DPtr _count_e ;
213    Histo1DPtr _count_mu;
214
215    Histo1DPtr _hist_eTmiss_e ;
216    Histo1DPtr _hist_eTmiss_mu;
217    //@}
218
219  };
220
221  // The hook for the plugin system
222  RIVET_DECLARE_PLUGIN(ATLAS_2012_CONF_2012_104);
223
224}