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