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