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