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