rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2012_CONF_2012_001

4 or more lepton plus missing transverse energy SUSY search
Experiment: ATLAS (LHC)
Status: PRELIMINARY
Authors:
  • Peter Richardson
References:
  • ATLAS-CONF-2012-001
  • ATLAS-CONF-2012-035
Beams: p+ p+
Beam energies: (3500.0, 3500.0) GeV
Run details:
  • BSM signal events at 7000 GeV.

Search for SUSY using events with 4 or more leptons in association with missing transverse energy in proton-proton collisions at a centre-of-mass energy of 7 TeV. The data sample has a total integrated luminosity of 2.06 fb$^{-1}$. There is no reference data and in addition to the control plots from the paper the number of events in the two signal regions, correctly normalized to an integrated luminosity 2.06 fb$^{-1}$, are calculated.

Source code: ATLAS_2012_CONF_2012_001.cc
  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_CONF_2012_001 : public Analysis {
 16  public:
 17
 18    /// Constructor
 19    ATLAS_2012_CONF_2012_001()
 20      : Analysis("ATLAS_2012_CONF_2012_001")
 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, FastJets::ANTIKT, 0.4), "AntiKtJets04");
 48
 49      // all tracks (to do deltaR with leptons)
 50      declare(ChargedFinalState(Cuts::abseta < 3.0),"cfs");
 51
 52      // Book histograms
 53      {Histo1DPtr tmp; _hist_leptonpT.push_back(book(tmp,1,1,1));}
 54      {Histo1DPtr tmp; _hist_leptonpT.push_back(book(tmp,2,1,1));}
 55      {Histo1DPtr tmp; _hist_leptonpT.push_back(book(tmp,3,1,1));}
 56      {Histo1DPtr tmp; _hist_leptonpT.push_back(book(tmp,4,1,1));}
 57      book(_hist_njet   ,5,1,1);
 58      book(_hist_etmiss ,6,1,1);
 59      book(_hist_mSFOS  ,7,1,1);
 60      book(_hist_meff   ,8,1,1);
 61
 62      {Histo1DPtr tmp; _hist_leptonpT_MC.push_back(book(tmp, "hist_lepton_pT_1", 26, 0., 260));}
 63      {Histo1DPtr tmp; _hist_leptonpT_MC.push_back(book(tmp, "hist_lepton_pT_2", 15, 0., 150));}
 64      {Histo1DPtr tmp; _hist_leptonpT_MC.push_back(book(tmp, "hist_lepton_pT_3", 20, 0., 100));}
 65      {Histo1DPtr tmp; _hist_leptonpT_MC.push_back(book(tmp, "hist_lepton_pT_4", 20, 0., 100));}
 66      book(_hist_njet_MC   ,"hist_njet", 7, -0.5, 6.5);
 67      book(_hist_etmiss_MC ,"hist_etmiss",11,0.,220.);
 68      book(_hist_mSFOS_MC  ,"hist_m_SFOS",13,0.,260.);
 69      book(_hist_meff_MC   ,"hist_m_eff",19,0.,950.);
 70
 71      book(_count_SR1 ,"count_SR1", 1, 0., 1.);
 72      book(_count_SR2 ,"count_SR2", 1, 0., 1.);
 73    }
 74
 75
 76    /// Perform the per-event analysis
 77    void analyze(const Event& event) {
 78      const double weight = 1.0;
 79      // get the jet candidates
 80      Jets cand_jets;
 81      for (const Jet& jet :
 82               apply<FastJets>(event, "AntiKtJets04").jetsByPt(20.0*GeV) ) {
 83        if ( fabs( jet.eta() ) < 2.8 ) {
 84          cand_jets.push_back(jet);
 85        }
 86      }
 87
 88      // candidate muons
 89      Particles cand_mu;
 90      Particles chg_tracks =
 91        apply<ChargedFinalState>(event, "cfs").particles();
 92      for ( const Particle & mu :
 93                apply<IdentifiedFinalState>(event, "muons").particlesByPt() ) {
 94        double pTinCone = -mu.pT();
 95        for ( const Particle & track : chg_tracks ) {
 96          if ( deltaR(mu.momentum(),track.momentum()) <= 0.2 )
 97            pTinCone += track.pT();
 98        }
 99        if ( pTinCone < 1.8*GeV )
100          cand_mu.push_back(mu);
101      }
102
103      // candidate electrons
104      Particles cand_e;
105      for ( const Particle & e :
106                apply<IdentifiedFinalState>(event, "elecs").particlesByPt() ) {
107        double eta = e.eta();
108        // remove electrons with pT<15 in old veto region
109        if( fabs(eta)>1.37 && fabs(eta) < 1.52 && e.perp()< 15.*GeV)
110          continue;
111        double pTinCone = -e.perp();
112        for ( const Particle & track : chg_tracks ) {
113          if ( deltaR(e.momentum(),track.momentum()) <= 0.2 )
114            pTinCone += track.pT();
115        }
116        if (pTinCone/e.perp()<0.1) {
117          cand_e.push_back(e);
118        }
119      }
120
121      // resolve jet/lepton ambiguity
122      Jets recon_jets;
123      for ( const Jet& jet : cand_jets ) {
124        bool away_from_e = true;
125        for ( const Particle & e : cand_e ) {
126          if ( deltaR(e.momentum(),jet.momentum()) <= 0.2 ) {
127            away_from_e = false;
128            break;
129          }
130        }
131        if ( away_from_e )
132          recon_jets.push_back( jet );
133      }
134
135      // only keep electrons more than R=0.4 from jets
136      Particles cand2_e;
137      for(unsigned int ie=0;ie<cand_e.size();++ie) {
138        const Particle & e = cand_e[ie];
139        // at least 0.4 from any jets
140        bool away = true;
141        for ( const Jet& jet : recon_jets ) {
142          if ( deltaR(e.momentum(),jet.momentum()) < 0.4 ) {
143            away = false;
144            break;
145          }
146        }
147        // and 0.1 from any muons
148        if ( away ) {
149          for ( const Particle & mu : cand_mu ) {
150            if ( deltaR(mu.momentum(),e.momentum()) < 0.1 ) {
151              away = false;
152              break;
153            }
154          }
155        }
156        // and 0.1 from electrons
157        for(unsigned int ie2=0;ie2<cand_e.size();++ie2) {
158          if(ie==ie2) continue;
159          if ( deltaR(e.momentum(),cand_e[ie2].momentum()) < 0.1 ) {
160            away = false;
161            break;
162          }
163        }
164        // if isolated keep it
165        if ( away ) cand2_e.push_back( e );
166      }
167      // remove e+e- pairs with mass < 20.
168      Particles recon_e;
169      for(unsigned int ie=0;ie<cand2_e.size();++ie) {
170	bool pass = true;
171	for(unsigned int ie2=0;ie2<cand2_e.size();++ie2) {
172	  if(cand2_e[ie].pid()*cand2_e[ie2].pid()>0) continue;
173	  double mtest = (cand2_e[ie].momentum()+cand2_e[ie2].momentum()).mass();
174	  if(mtest<=20.) {
175	    pass = false;
176	    break;
177	  }
178	}
179	if(pass) recon_e.push_back(cand2_e[ie]);
180      }
181
182      // only keep muons more than R=0.4 from jets
183      Particles cand2_mu;
184      for(unsigned int imu=0;imu<cand_mu.size();++imu) {
185        const Particle & mu = cand_mu[imu];
186        bool away = true;
187        // at least 0.4 from any jets
188        for ( const Jet& jet : recon_jets ) {
189          if ( deltaR(mu.momentum(),jet.momentum()) < 0.4 ) {
190            away = false;
191            break;
192          }
193        }
194        // and 0.1 from any electrona
195        if ( away ) {
196          for ( const Particle & e : cand_e ) {
197            if ( deltaR(mu.momentum(),e.momentum()) < 0.1 ) {
198              away = false;
199              break;
200            }
201          }
202        }
203        // and 0.1 from muons
204        for(unsigned int imu2=0;imu2<cand_mu.size();++imu2) {
205          if(imu==imu2) continue;
206          if ( deltaR(mu.momentum(),cand_mu[imu2].momentum()) < 0.1 ) {
207            away = false;
208            break;
209          }
210        }
211        if ( away )
212          cand2_mu.push_back( mu );
213      }
214
215      // remove mu+mu- pairs with mass < 20.
216      Particles recon_mu;
217      for(unsigned int imu=0;imu<cand2_mu.size();++imu) {
218	bool pass = true;
219	for(unsigned int imu2=0;imu2<cand2_mu.size();++imu2) {
220	  if(cand2_mu[imu].pid()*cand2_mu[imu2].pid()>0) continue;
221	  double mtest = (cand2_mu[imu].momentum()+cand2_mu[imu2].momentum()).mass();
222	  if(mtest<=20.) {
223	    pass = false;
224	    break;
225	  }
226	}
227	if(pass) recon_mu.push_back(cand2_mu[imu]);
228      }
229
230      // pTmiss
231      Particles vfs_particles =
232        apply<VisibleFinalState>(event, "vfs").particles();
233      FourMomentum pTmiss;
234      for ( const Particle & p : vfs_particles ) {
235        pTmiss -= p.momentum();
236      }
237      double eTmiss = pTmiss.pT();
238
239      // now only use recon_jets, recon_mu, recon_e
240
241      // reject events with less than 4 electrons and muons
242      if ( recon_mu.size() + recon_e.size() < 4 ) {
243        MSG_DEBUG("To few charged leptons left after selection");
244        vetoEvent;
245      }
246
247      // ATLAS calo problem
248      if(rand()/static_cast<double>(RAND_MAX)<=0.42) {
249        for ( const Particle & e : recon_e ) {
250          double eta = e.eta();
251          double phi = e.azimuthalAngle(MINUSPI_PLUSPI);
252          if(eta>-0.1&&eta<1.5&&phi>-0.9&&phi<-0.5)
253            vetoEvent;
254        }
255        for ( const Jet & jet : recon_jets ) {
256          double eta = jet.rapidity();
257          double phi = jet.azimuthalAngle(MINUSPI_PLUSPI);
258          if(jet.perp()>40 && eta>-0.1&&eta<1.5&&phi>-0.9&&phi<-0.5)
259            vetoEvent;
260        }
261      }
262
263      // check at least one e/mu passing trigger
264      if( !( !recon_e .empty() && recon_e[0] .perp()>25.)  &&
265          !( !recon_mu.empty() && recon_mu[0].perp()>20.) ) {
266        MSG_DEBUG("Hardest lepton fails trigger");
267        vetoEvent;
268      }
269
270      // calculate meff
271      double meff = eTmiss;
272      for ( const Particle & e  : recon_e  )
273        meff += e.perp();
274      for ( const Particle & mu : recon_mu )
275        meff += mu.perp();
276      for ( const Jet & jet : recon_jets ) {
277        double pT = jet.perp();
278        if(pT>40.) meff += pT;
279      }
280
281      double mSFOS=1e30, mdiff=1e30;
282      // mass of SFOS pairs closest to the Z mass
283      for(unsigned int ix=0;ix<recon_e.size();++ix) {
284        for(unsigned int iy=ix+1;iy<recon_e.size();++iy) {
285          if(recon_e[ix].pid()*recon_e[iy].pid()>0) continue;
286          double mtest = (recon_e[ix].momentum()+recon_e[iy].momentum()).mass();
287          if(fabs(mtest-90.)<mdiff) {
288            mSFOS = mtest;
289            mdiff = fabs(mtest-90.);
290          }
291        }
292      }
293      for(unsigned int ix=0;ix<recon_mu.size();++ix) {
294        for(unsigned int iy=ix+1;iy<recon_mu.size();++iy) {
295          if(recon_mu[ix].pid()*recon_mu[iy].pid()>0) continue;
296          double mtest = (recon_mu[ix].momentum()+recon_mu[iy].momentum()).mass();
297          if(fabs(mtest-91.118)<mdiff) {
298            mSFOS = mtest;
299            mdiff = fabs(mtest-91.118);
300          }
301        }
302      }
303
304      // make the control plots
305      // lepton pT
306      unsigned int ie=0,imu=0;
307      for(unsigned int ix=0;ix<4;++ix) {
308        double pTe  = ie <recon_e .size() ?
309          recon_e [ie ].perp() : -1*GeV;
310        double pTmu = imu<recon_mu.size() ?
311          recon_mu[imu].perp() : -1*GeV;
312        if(pTe>pTmu) {
313          _hist_leptonpT   [ix]->fill(pTe ,weight);
314          _hist_leptonpT_MC[ix]->fill(pTe ,weight);
315          ++ie;
316        }
317        else {
318          _hist_leptonpT   [ix]->fill(pTmu,weight);
319          _hist_leptonpT_MC[ix]->fill(pTmu,weight);
320          ++imu;
321        }
322      }
323      // njet
324      _hist_njet   ->fill(recon_jets.size(),weight);
325      _hist_njet_MC->fill(recon_jets.size(),weight);
326      // etmiss
327      _hist_etmiss   ->fill(eTmiss,weight);
328      _hist_etmiss_MC->fill(eTmiss,weight);
329      if(mSFOS<1e30) {
330        _hist_mSFOS   ->fill(mSFOS,weight);
331        _hist_mSFOS_MC->fill(mSFOS,weight);
332      }
333      _hist_meff   ->fill(meff,weight);
334      _hist_meff_MC->fill(meff,weight);
335
336      // finally the counts
337      if(eTmiss>50.) {
338        _count_SR1->fill(0.5,weight);
339        if(mdiff>10.) _count_SR2->fill(0.5,weight);
340      }
341    }
342
343    //@}
344
345    void finalize() {
346      double norm = crossSection()/femtobarn*2.06/sumOfWeights();
347      // these are number of events at 2.06fb^-1 per 10 GeV
348      scale(_hist_leptonpT   [0],norm*10.);
349      scale(_hist_leptonpT   [1],norm*10.);
350      scale(_hist_leptonpT_MC[0],norm*10.);
351      scale(_hist_leptonpT_MC[1],norm*10.);
352      // these are number of events at 2.06fb^-1 per 5 GeV
353      scale(_hist_leptonpT   [2],norm*5.);
354      scale(_hist_leptonpT   [3],norm*5.);
355      scale(_hist_leptonpT_MC[2],norm*5.);
356      scale(_hist_leptonpT_MC[3],norm*5.);
357      // these are number of events at 2.06fb^-1 per 20 GeV
358      scale(_hist_etmiss      ,norm*20.);
359      scale(_hist_mSFOS       ,norm*20.);
360      scale(_hist_etmiss_MC   ,norm*20.);
361      scale(_hist_mSFOS_MC    ,norm*20.);
362      // these are number of events at 2.06fb^-1 per 50 GeV
363      scale(_hist_meff        ,norm*50.);
364      scale(_hist_meff_MC     ,norm*50.);
365      // these are number of events at 2.06fb^-1
366      scale(_hist_njet        ,norm);
367      scale(_hist_njet_MC     ,norm);
368      scale(_count_SR1,norm);
369      scale(_count_SR2,norm);
370    }
371
372  private:
373
374    /// @name Histograms
375    //@{
376    vector<Histo1DPtr> _hist_leptonpT,_hist_leptonpT_MC;
377    Histo1DPtr _hist_njet;
378    Histo1DPtr _hist_njet_MC;
379    Histo1DPtr _hist_etmiss;
380    Histo1DPtr _hist_etmiss_MC;
381    Histo1DPtr _hist_mSFOS;
382    Histo1DPtr _hist_mSFOS_MC;
383    Histo1DPtr _hist_meff;
384    Histo1DPtr _hist_meff_MC;
385    Histo1DPtr _count_SR1;
386    Histo1DPtr _count_SR2;
387    //@}
388
389  };
390
391  // The hook for the plugin system
392  RIVET_DECLARE_PLUGIN(ATLAS_2012_CONF_2012_001);
393
394}