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/VetoedFinalState.hh"
  7#include "Rivet/Projections/IdentifiedFinalState.hh"
  8#include "Rivet/Projections/FastJets.hh"
  9#include "Rivet/Tools/RivetMT2.hh"
 10
 11namespace Rivet {
 12
 13
 14  /// @author Peter Richardson
 15  class ATLAS_2012_I1190891 : public Analysis {
 16  public:
 17
 18    /// Constructor
 19    ATLAS_2012_I1190891()
 20      : Analysis("ATLAS_2012_I1190891")
 21    {    }
 22
 23
 24    /// @name Analysis methods
 25    /// @{
 26
 27    /// Book histograms and initialise projections before the run
 28    void init() {
 29
 30      // projection to find the electrons
 31      IdentifiedFinalState elecs(Cuts::abseta < 2.47 && Cuts::pT > 10*GeV);
 32      elecs.acceptIdPair(PID::ELECTRON);
 33      declare(elecs, "elecs");
 34
 35      // projection to find the muons
 36      IdentifiedFinalState muons(Cuts::abseta < 2.4 && Cuts::pT > 10*GeV);
 37      muons.acceptIdPair(PID::MUON);
 38      declare(muons, "muons");
 39
 40      // for pTmiss
 41      declare(VisibleFinalState(Cuts::abseta < 4.9),"vfs");
 42
 43      VetoedFinalState vfs;
 44      vfs.addVetoPairId(PID::MUON);
 45
 46      /// Jet finder
 47      declare(FastJets(vfs, JetAlg::ANTIKT, 0.4), "AntiKtJets04");
 48
 49      // all tracks (to do deltaR with leptons)
 50      declare(ChargedFinalState((Cuts::etaIn(-3.0,3.0))),"cfs");
 51
 52      // Book histograms
 53      book(_hist_etmiss ,"hist_etmiss",10,0.,500.);
 54      book(_hist_meff   ,"hist_m_eff",7,0.,1050.);
 55      book(_count_SR1 ,"count_SR1", 1, 0., 1.);
 56      book(_count_SR2 ,"count_SR2", 1, 0., 1.);
 57    }
 58
 59
 60    /// Perform the per-event analysis
 61    void analyze(const Event& event) {
 62
 63      // get the jet candidates
 64      Jets cand_jets = apply<FastJets>(event, "AntiKtJets04").jetsByPt(Cuts::pT > 20*GeV && Cuts::abseta < 2.5);
 65
 66      // candidate muons
 67      Particles cand_mu;
 68      Particles chg_tracks =
 69        apply<ChargedFinalState>(event, "cfs").particles();
 70      for ( const Particle & mu :
 71                apply<IdentifiedFinalState>(event, "muons").particlesByPt() ) {
 72        double pTinCone = -mu.pT();
 73        for ( const Particle & track : chg_tracks ) {
 74          if ( deltaR(mu.momentum(),track.momentum()) <= 0.2 )
 75            pTinCone += track.pT();
 76        }
 77        if ( pTinCone < 1.8*GeV )
 78          cand_mu.push_back(mu);
 79      }
 80
 81      // candidate electrons
 82      Particles cand_e;
 83      for ( const Particle & e :
 84                apply<IdentifiedFinalState>(event, "elecs").particlesByPt() ) {
 85        double pTinCone = -e.perp();
 86        for ( const Particle & track : chg_tracks ) {
 87          if ( deltaR(e.momentum(),track.momentum()) <= 0.2 )
 88            pTinCone += track.pT();
 89        }
 90        if (pTinCone/e.perp()<0.1) {
 91          cand_e.push_back(e);
 92        }
 93      }
 94
 95      // resolve jet/lepton ambiguity
 96      Jets recon_jets;
 97      for ( const Jet& jet : cand_jets ) {
 98        bool away_from_e = true;
 99        for ( const Particle & e : cand_e ) {
100          if ( deltaR(e.momentum(),jet.momentum()) <= 0.2 ) {
101            away_from_e = false;
102            break;
103          }
104        }
105        if ( away_from_e )
106          recon_jets.push_back( jet );
107      }
108
109      // only keep electrons more than R=0.4 from jets
110      Particles cand2_e;
111      for(unsigned int ie=0;ie<cand_e.size();++ie) {
112        const Particle & e = cand_e[ie];
113        // at least 0.4 from any jets
114        bool away = true;
115        for ( const Jet& jet : recon_jets ) {
116          if ( deltaR(e.momentum(),jet.momentum()) < 0.4 ) {
117            away = false;
118            break;
119          }
120        }
121        // and 0.1 from any muons
122        if ( away ) {
123          for ( const Particle & mu : cand_mu ) {
124            if ( deltaR(mu.momentum(),e.momentum()) < 0.1 ) {
125              away = false;
126              break;
127            }
128          }
129        }
130        // and 0.1 from electrons ( keep higher energy)
131        for(unsigned int ie2=0;ie2<cand_e.size();++ie2) {
132          if(ie==ie2) continue;
133          if ( deltaR(e.momentum(),cand_e[ie2].momentum()) < 0.1 &&
134               e.E() < cand_e[ie2].E() ) {
135            away = false;
136            break;
137          }
138        }
139        // if isolated keep it
140        if ( away )
141          cand2_e.push_back( e );
142      }
143
144      // remove e+e- pairs with mass < 20.
145      Particles recon_e;
146      for(unsigned int ie=0;ie<cand2_e.size();++ie) {
147        bool pass = true;
148        for(unsigned int ie2=0;ie2<cand2_e.size();++ie2) {
149          if(cand2_e[ie].pid()*cand2_e[ie2].pid()>0) continue;
150          double mtest = (cand2_e[ie].momentum()+cand2_e[ie2].momentum()).mass();
151          if(mtest<=20.) {
152            pass = false;
153            break;
154          }
155        }
156        if(pass) recon_e.push_back(cand2_e[ie]);
157      }
158
159      // only keep muons more than R=0.4 from jets
160      Particles cand2_mu;
161      for(unsigned int imu=0;imu<cand_mu.size();++imu) {
162        const Particle & mu = cand_mu[imu];
163        bool away = true;
164        // at least 0.4 from any jets
165        for ( const Jet& jet : recon_jets ) {
166          if ( deltaR(mu.momentum(),jet.momentum()) < 0.4 ) {
167            away = false;
168            break;
169          }
170        }
171        // and 0.1 from any electrona
172        if ( away ) {
173          for ( const Particle & e : cand_e ) {
174            if ( deltaR(mu.momentum(),e.momentum()) < 0.1 ) {
175              away = false;
176              break;
177            }
178          }
179        }
180        if ( away )
181          cand2_mu.push_back( mu );
182      }
183
184      // remove mu+mu- pairs with mass < 20.
185      Particles recon_mu;
186      for(unsigned int imu=0;imu<cand2_mu.size();++imu) {
187        bool pass = true;
188        for(unsigned int imu2=0;imu2<cand2_mu.size();++imu2) {
189          if(cand2_mu[imu].pid()*cand2_mu[imu2].pid()>0) continue;
190          double mtest = (cand2_mu[imu].momentum()+cand2_mu[imu2].momentum()).mass();
191          if(mtest<=20.) {
192            pass = false;
193            break;
194          }
195        }
196        if(pass) recon_mu.push_back(cand2_mu[imu]);
197      }
198
199      // pTmiss
200      Particles vfs_particles =
201        apply<VisibleFinalState>(event, "vfs").particles();
202      FourMomentum pTmiss;
203      for ( const Particle & p : vfs_particles ) {
204        pTmiss -= p.momentum();
205      }
206      double eTmiss = pTmiss.pT();
207
208      // now only use recon_jets, recon_mu, recon_e
209
210      // reject events with less than 4 electrons and muons
211      if ( recon_mu.size() + recon_e.size() < 4 ) {
212        MSG_DEBUG("To few charged leptons left after selection");
213        vetoEvent;
214      }
215
216      // check if passes single lepton trigger
217      bool passSingle =
218        ( !recon_e .empty() && recon_e[0] .perp()>25. )||
219        ( !recon_mu.empty() && recon_mu[0].perp()>20.);
220
221      // or two lepton trigger
222      bool passDouble =
223        ( recon_mu.size()>=2 && recon_mu[1].perp()>12.) ||
224        ( recon_e .size()>=2 && recon_e [1].perp()>17.) ||
225        ( !recon_e.empty() && !recon_mu.empty() &&
226          recon_e[0].perp()>15. &&  recon_mu[0].perp()>10.);
227
228      // must pass a trigger
229      if( !passSingle && !passDouble ) {
230        MSG_DEBUG("Hardest lepton fails trigger");
231        vetoEvent;
232      }
233
234      // calculate meff
235      double meff = eTmiss;
236      for ( const Particle & e  : recon_e  )
237        meff += e.perp();
238      for ( const Particle & mu : recon_mu )
239        meff += mu.perp();
240      for ( const Jet & jet : recon_jets ) {
241        double pT = jet.perp();
242        if(pT>40.) meff += pT;
243      }
244
245      // mass of SFOS pairs closest to the Z mass
246      for(unsigned int ix=0;ix<recon_e.size();++ix) {
247        for(unsigned int iy=ix+1;iy<recon_e.size();++iy) {
248          if(recon_e[ix].pid()*recon_e[iy].pid()>0) continue;
249          double mtest = (recon_e[ix].momentum()+recon_e[iy].momentum()).mass();
250          if(mtest>81.2 && mtest<101.2) vetoEvent;
251        }
252      }
253      for(unsigned int ix=0;ix<recon_mu.size();++ix) {
254        for(unsigned int iy=ix+1;iy<recon_mu.size();++iy) {
255          if(recon_mu[ix].pid()*recon_mu[iy].pid()>0) continue;
256          double mtest = (recon_mu[ix].momentum()+recon_mu[iy].momentum()).mass();
257          if(mtest>81.2 && mtest<101.2) vetoEvent;
258        }
259      }
260
261      // make the control plots
262      _hist_etmiss ->fill(eTmiss);
263      _hist_meff   ->fill(meff  );
264      // finally the counts
265      if(eTmiss>50.) _count_SR1->fill(0.5);
266      if(meff >300.) _count_SR2->fill(0.5);
267    }
268
269    /// @}
270
271    void finalize() {
272      double norm = crossSection()/femtobarn*4.7/sumOfWeights();
273      scale(_hist_etmiss,norm* 50.);
274      scale(_hist_meff  ,norm*150.);
275      scale(_count_SR1,norm);
276      scale(_count_SR2,norm);
277    }
278
279  private:
280
281    /// @name Histograms
282    /// @{
283    Histo1DPtr _hist_etmiss;
284    Histo1DPtr _hist_meff;
285    Histo1DPtr _count_SR1;
286    Histo1DPtr _count_SR2;
287    /// @}
288
289  };
290
291  RIVET_DECLARE_PLUGIN(ATLAS_2012_I1190891);
292
293}