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_I1186556 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    ATLAS_2012_I1186556()
 18      : Analysis("ATLAS_2012_I1186556")
 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 > 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      VetoedFinalState vfs;
 40      vfs.addVetoPairId(PID::MUON);
 41      declare(FastJets(vfs, JetAlg::ANTIKT, 0.4), "AntiKtJets04");
 42
 43      // all tracks (to do deltaR with leptons)
 44      declare(ChargedFinalState(Cuts::abseta < 3.0 && Cuts::pT > 1*GeV), "cfs");
 45
 46      // for pTmiss
 47      declare(VisibleFinalState(Cuts::abseta < 4.9),"vfs");
 48
 49      // Book histograms
 50      book(_count_SR_SF     ,"count_SR_SF"    , 1, 0., 1.);
 51      book(_count_SR_OF     ,"count_SR_OF"    , 1, 0., 1.);
 52
 53      book(_hist_mT2_SF_exp ,"hist_mT2_SF_exp", 40 , 0., 200. );
 54      book(_hist_mT2_OF_exp ,"hist_mT2_OF_exp", 40 , 0., 200. );
 55      book(_hist_mT2_SF_MC  ,"hist_mT2_SF_MC" , 500, 0., 1000.);
 56      book(_hist_mT2_OF_MC  ,"hist_mT2_OF_MC" , 500, 0., 1000.);
 57
 58    }
 59
 60    /// Perform the per-event analysis
 61    void analyze(const Event& event) {
 62
 63      // get the candiate jets
 64      Jets cand_jets;
 65      for ( const Jet& jet :
 66                apply<FastJets>(event, "AntiKtJets04").jetsByPt(Cuts::pT > 20*GeV && Cuts::abseta < 4.5) ) {
 67        cand_jets.push_back(jet);
 68      }
 69      // charged tracks for isolation
 70      Particles chg_tracks =
 71        apply<ChargedFinalState>(event, "cfs").particles();
 72      // find the electrons
 73      Particles cand_e;
 74      for( const Particle & e :
 75               apply<IdentifiedFinalState>(event, "elecs").particlesByPt()) {
 76        // remove any leptons within 0.4 of any candidate jets
 77        bool e_near_jet = false;
 78        for ( const Jet& jet : cand_jets ) {
 79          double dR = deltaR(e.momentum(),jet.momentum());
 80          if ( dR < 0.4 && dR > 0.2 ) {
 81            e_near_jet = true;
 82            break;
 83          }
 84        }
 85	if ( e_near_jet ) continue;
 86	cand_e.push_back(e);
 87      }
 88      Particles cand_mu;
 89      for( const Particle & mu :
 90               apply<IdentifiedFinalState>(event, "muons").particlesByPt()) {
 91        // remove any leptons within 0.4 of any candidate jets
 92        bool mu_near_jet = false;
 93        for ( const Jet& jet : cand_jets ) {
 94          if ( deltaR(mu.momentum(),jet.momentum()) < 0.4 ) {
 95            mu_near_jet = true;
 96            break;
 97          }
 98        }
 99        if ( mu_near_jet ) continue;
100	cand_mu.push_back(mu);
101      }
102      // pTcone around muon track
103      Particles recon_mu;
104      for ( const Particle & mu : cand_mu ) {
105        double pTinCone = -mu.pT();
106        for ( const Particle & track : chg_tracks ) {
107          if ( deltaR(mu.momentum(),track.momentum()) < 0.2 )
108            pTinCone += track.pT();
109        }
110        if ( pTinCone < 1.8*GeV ) recon_mu.push_back(mu);
111      }
112      // pTcone around electron track
113      Particles recon_e;
114      for ( const Particle & e : cand_e ) {
115        double pTinCone = -e.pT();
116        for ( const Particle & track : chg_tracks ) {
117          if ( deltaR(e.momentum(),track.momentum()) < 0.2 )
118            pTinCone += track.pT();
119        }
120        if ( pTinCone < 0.1 * e.pT() ) recon_e.push_back(e);
121      }
122
123      // pTmiss
124      FourMomentum pTmiss;
125      for ( const Particle & p :
126                apply<VisibleFinalState>(event, "vfs").particles() ) {
127        pTmiss -= p.momentum();
128      }
129
130      // discard jets that overlap with electrons
131      Jets recon_jets;
132      for ( const Jet& jet : cand_jets ) {
133        if(jet.abseta()>2.5||
134           jet.perp()<20.) continue;
135	bool away_from_e = true;
136	for ( const Particle & e : cand_e ) {
137	  if ( deltaR(e.momentum(),jet.momentum()) < 0.2 ) {
138	    away_from_e = false;
139	    break;
140	  }
141	}
142	if ( away_from_e ) recon_jets.push_back( jet );
143      }
144
145      // put leptons into 1 vector and order by pT
146      Particles leptons(recon_e.begin(),recon_e.end());
147      leptons.insert(leptons.begin(),recon_mu.begin(),recon_mu.end());
148      sort(leptons.begin(),leptons.end(),cmpMomByPt);
149
150      // exactly two leptons
151      if(leptons.size() !=2) vetoEvent;
152
153      // hardest lepton pT greater the 25 (20) e(mu)
154      if( (leptons[0].abspid()==PID::ELECTRON && leptons[0].perp()<25.) ||
155	  (leptons[0].abspid()==PID::ELECTRON && leptons[0].perp()<20.))
156	vetoEvent;
157
158      // require opposite sign
159      if(leptons[0].pid()*leptons[1].pid()>0) vetoEvent;
160
161      // and invariant mass > 20
162      double mll = (leptons[0].momentum() + leptons[1].momentum()).mass();
163      if(mll<20.) vetoEvent;
164
165      // two jets 1st pT > 50 and second pT> 25
166      if(recon_jets.size()<2 || recon_jets[0].perp()<50. ||
167	 recon_jets[1].perp()<25.) vetoEvent;
168
169      // calculate mT2
170      double m_T2 = mT2( leptons[0], leptons[1], pTmiss,0.0 ); // zero mass invisibles
171
172      // same flavour region
173      if(leptons[0].pid()==-leptons[1].pid()) {
174	// remove Z region
175	if(mll>71.&&mll<111.) vetoEvent;
176	// require at least 1 b jet
177	unsigned int n_b=0;
178	for(unsigned int ix=0;ix<recon_jets.size();++ix) {
179	   if(recon_jets[ix].bTagged() && rand()/static_cast<double>(RAND_MAX)<=0.60)
180	     ++n_b;
181	}
182	if(n_b==0) vetoEvent;
183	_hist_mT2_SF_exp->fill(m_T2);
184	_hist_mT2_SF_MC ->fill(m_T2);
185	if(m_T2>120.) _count_SR_SF->fill(0.5);
186      }
187      // opposite flavour region
188      else {
189	_hist_mT2_OF_exp->fill(m_T2);
190	_hist_mT2_OF_MC ->fill(m_T2);
191	if(m_T2>120.) _count_SR_OF->fill(0.5);
192      }
193    }
194    /// @}
195
196
197    void finalize() {
198
199      double norm = 4.7* crossSection()/sumOfWeights()/femtobarn;
200      scale(_count_SR_SF    ,   norm);
201      scale(_count_SR_OF    ,   norm);
202      scale(_hist_mT2_SF_exp,5.*norm);
203      scale(_hist_mT2_OF_exp,5.*norm);
204      scale(_hist_mT2_SF_MC ,   norm/4.7);
205      scale(_hist_mT2_OF_MC ,   norm/4.7);
206
207    }
208
209  private:
210
211    /// @name Histograms
212    /// @{
213    Histo1DPtr _count_SR_SF;
214    Histo1DPtr _count_SR_OF;
215
216    Histo1DPtr _hist_mT2_SF_exp;
217    Histo1DPtr _hist_mT2_OF_exp;
218    Histo1DPtr _hist_mT2_SF_MC;
219    Histo1DPtr _hist_mT2_OF_MC;
220    /// @}
221
222  };
223
224  RIVET_DECLARE_PLUGIN(ATLAS_2012_I1186556);
225
226}