Rivet analyses referenceMC_SUSYValidate generic SUSY events, including various lepton invariant massExperiment: () Status: VALIDATED Authors:
Beams: * * Beam energies: ANY Run details:
Analysis of generic SUSY events at the LHC, based on Atlas Herwig++ validation analysis contents. Plotted are $\eta$, $\phi$ and $p_\perp$ observables for charged tracks, photons, isolated photons, electrons, muons, and jets, as well as various dilepton mass `edge' plots for different event selection criteria. Source code: MC_SUSY.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/ChargedFinalState.hh"
5#include "Rivet/Projections/FastJets.hh"
6#include "Rivet/Projections/IdentifiedFinalState.hh"
7#include "Rivet/Projections/MissingMomentum.hh"
8#include "Rivet/Projections/LeadingParticlesFinalState.hh"
9
10namespace Rivet {
11
12
13
14
15 /// @brief MC validation analysis for SUSY events
16 /// @author Andy Buckley
17 class MC_SUSY : public Analysis {
18 public:
19
20 /// Constructor
21 MC_SUSY()
22 : Analysis("MC_SUSY")
23 { }
24
25
26 /// @name Analysis methods
27 /// @{
28
29 // Book histograms
30 void init() {
31 // Basic final state
32 const FinalState fs((Cuts::etaIn(-4.0, 4.0) && Cuts::pT >= 10*GeV));
33
34 // Tracks and jets
35 declare(ChargedFinalState(fs), "Tracks");
36 declare(FastJets(fs, JetAlg::ANTIKT, 0.7), "Jets");
37
38 IdentifiedFinalState photonfs(fs);
39 photonfs.acceptId(PID::PHOTON);
40 declare(photonfs, "AllPhotons");
41
42 IdentifiedFinalState efs(fs);
43 efs.acceptIdPair(PID::ELECTRON);
44 declare(efs, "Electrons");
45
46 IdentifiedFinalState mufs(fs);
47 mufs.acceptIdPair(PID::MUON);
48 declare(mufs, "Muons");
49
50 MissingMomentum missing(fs);
51 declare(missing, "MET");
52
53 LeadingParticlesFinalState lpfs(fs);
54 lpfs.addParticleIdPair(PID::ELECTRON);
55 lpfs.addParticleIdPair(PID::MUON);
56 declare(lpfs, "LeadingParticles");
57
58 book(_hist_n_trk ,"n-trk", 50, 0.5, 300.5);
59 book(_hist_phi_trk ,"phi-trk", 50, -PI, PI);
60 book(_hist_eta_trk ,"eta-trk", 50, -4, 4);
61 book(_hist_pt_trk ,"pt-trk", 100, 0.0, 1500);
62
63 book(_hist_n_jet ,"n-jet", 21, -0.5, 20.5);
64 book(_hist_phi_jet ,"phi-jet", 50, -PI, PI);
65 book(_hist_eta_jet ,"eta-jet", 50, -4, 4);
66 book(_hist_pt_jet ,"pt-jet", 100, 0.0, 1500);
67
68 book(_hist_n_e ,"n-e", 11, -0.5, 10.5);
69 book(_hist_phi_e ,"phi-e", 50, -PI, PI);
70 book(_hist_eta_e ,"eta-e", 50, -4, 4);
71 book(_hist_pt_e ,"pt-e", 100, 0.0, 500);
72
73 book(_hist_n_mu ,"n-mu", 11, -0.5, 10.5);
74 book(_hist_phi_mu ,"phi-mu", 50, -PI, PI);
75 book(_hist_eta_mu ,"eta-mu", 50, -4, 4);
76 book(_hist_pt_mu ,"pt-mu", 100, 0.0, 500);
77
78 book(_hist_n_gamma ,"n-gamma", 11, -0.5, 10.5);
79 book(_hist_phi_gamma ,"phi-gamma", 50, -PI, PI);
80 book(_hist_eta_gamma ,"eta-gamma", 50, -4, 4);
81 book(_hist_pt_gamma ,"pt-gamma", 100, 0.0, 500);
82
83 book(_hist_n_gammaiso ,"n-gamma-iso", 11, -0.5, 10.5);
84 book(_hist_phi_gammaiso ,"phi-gamma-iso", 50, -PI, PI);
85 book(_hist_eta_gammaiso ,"eta-gamma-iso", 50, -4, 4);
86 book(_hist_pt_gammaiso ,"pt-gamma-iso", 100, 0.0, 500);
87
88 book(_hist_met ,"Etmiss", 100, 0.0, 1500);
89
90 book(_hist_mll_ossf_ee ,"mll-ossf-ee", 50, 0.0, 500);
91 book(_hist_mll_ossf_mumu ,"mll-ossf-mumu", 50, 0.0, 500);
92 book(_hist_mll_osof_emu ,"mll-osof-emu", 50, 0.0, 500);
93
94 book(_hist_mll_all_ossf_ee ,"mll-all-ossf-ee", 50, 0.0, 500);
95 book(_hist_mll_all_ossf_mumu ,"mll-all-ossf-mumu", 50, 0.0, 500);
96 book(_hist_mll_all_osof_emu ,"mll-all-osof-emu", 50, 0.0, 500);
97
98 book(_hist_mll_2_ossf_ee ,"mll-2-ossf-ee", 50, 0.0, 500);
99 book(_hist_mll_2_ossf_mumu ,"mll-2-ossf-mumu", 50, 0.0, 500);
100 book(_hist_mll_2_osof_emu ,"mll-2-osof-emu", 50, 0.0, 500);
101
102 /// @todo LSP eta, pT, phi, mass: no reliable cross-scenario LSP PID but
103 /// maybe plot for all of chi^0_1, gravitino, sneutrino, gluino, ... or
104 /// identify the LSP as any PID::isSUSY (?) particle with status = 1?
105 }
106
107
108 // Do the analysis
109 void analyze(const Event& evt) {
110 const FinalState& tracks = apply<FinalState>(evt, "Tracks");
111 if (tracks.particles().empty()) {
112 MSG_DEBUG("Failed multiplicity cut");
113 vetoEvent;
114 }
115
116 // Fill track histos
117 _hist_n_trk->fill(tracks.size());
118 for (const Particle& t : tracks.particles()) {
119 const FourMomentum& p = t.momentum();
120 _hist_phi_trk->fill(mapAngleMPiToPi(p.phi()));
121 _hist_eta_trk->fill(p.eta());
122 _hist_pt_trk->fill(p.pT()/GeV);
123 }
124
125 // Get jets and fill jet histos
126 const FastJets& jetpro = apply<FastJets>(evt, "Jets");
127 const Jets jets = jetpro.jetsByPt();
128 MSG_DEBUG("Jet multiplicity = " << jets.size());
129 _hist_n_jet->fill(jets.size());
130 for (const Jet& j : jets) {
131 const FourMomentum& pj = j.momentum();
132 _hist_phi_jet->fill(mapAngleMPiToPi(pj.phi()));
133 _hist_eta_jet->fill(pj.eta());
134 _hist_pt_jet->fill(pj.pT()/GeV);
135 }
136
137 /// @todo Resum photons around electrons
138
139 // Fill final state electron/positron histos
140 const FinalState& efs = apply<FinalState>(evt, "Electrons");
141 _hist_n_e->fill(efs.size());
142 vector<FourMomentum> epluses, eminuses;
143 for (const Particle& e : efs.particles()) {
144 const FourMomentum& p = e.momentum();
145 _hist_phi_e->fill(mapAngleMPiToPi(p.phi()));
146 _hist_eta_e->fill(p.eta());
147 _hist_pt_e->fill(p.pT()/GeV);
148 // Add sufficiently hard leptons to collections for m_ll histo
149 if (p.pT()/GeV > 20) {
150 if (PID::charge3(e.pid()) > 0) epluses += p; else eminuses += p;
151 }
152 }
153
154 /// @todo Resum photons around muons
155
156 // Fill final state muon/antimuon histos
157 const FinalState& mufs = apply<FinalState>(evt, "Muons");
158 _hist_n_mu->fill(mufs.size());
159 vector<FourMomentum> mupluses, muminuses;
160 for (const Particle& mu : mufs.particles()) {
161 const FourMomentum& p = mu.momentum();
162 _hist_phi_mu->fill(mapAngleMPiToPi(p.phi()));
163 _hist_eta_mu->fill(p.eta());
164 _hist_pt_mu->fill(p.pT()/GeV);
165 // Add sufficiently hard leptons to collections for m_ll histo
166 if (p.pT()/GeV > 20) {
167 if (PID::charge3(mu.pid()) > 0) mupluses += p; else muminuses += p;
168 }
169 }
170
171 // Fill final state non-isolated photon histos
172 const FinalState& allphotonfs = apply<FinalState>(evt, "AllPhotons");
173 _hist_n_gamma->fill(allphotonfs.size());
174 Particles isolatedphotons;
175 for (const Particle& ph : allphotonfs.particles()) {
176 const FourMomentum& p = ph.momentum();
177 _hist_phi_gamma->fill(mapAngleMPiToPi(p.phi()));
178 _hist_eta_gamma->fill(p.eta());
179 _hist_pt_gamma->fill(p.pT()/GeV);
180 // Select isolated photons
181 bool isolated = true;
182 for (const Jet& j : jets) {
183 if (deltaR(j.momentum(), p) < 0.2) {
184 isolated = false;
185 break;
186 }
187 }
188 if (isolated) isolatedphotons += ph;
189 }
190
191
192 // Fill final state isolated photon histos
193 _hist_n_gammaiso->fill(isolatedphotons.size());
194 for (const Particle& ph_iso : isolatedphotons) {
195 const FourMomentum& p = ph_iso.momentum();
196 _hist_phi_gammaiso->fill(mapAngleMPiToPi(p.phi()));
197 _hist_eta_gammaiso->fill(p.eta());
198 _hist_pt_gammaiso->fill(p.pT()/GeV);
199 }
200
201 // Calculate and fill missing Et histos
202 const MissingMomentum& met = apply<MissingMomentum>(evt, "MET");
203 _hist_met->fill(met.vectorEt().mod()/GeV);
204
205 // Choose highest-pT leptons of each sign and flavour for dilepton mass edges
206 const FinalState& lpfs = apply<FinalState>(evt, "LeadingParticles");
207 bool eplus_ok(false), eminus_ok(false), muplus_ok(false), muminus_ok(false);
208 FourMomentum peplus, peminus, pmuplus, pmuminus;
209 for (const Particle& p : lpfs.particles()) {
210 // Only use leptons above 20 GeV
211 if (p.pT()/GeV < 20) continue;
212 // Identify the PID
213 const PdgId pid = p.pid();
214 if (pid == PID::ELECTRON) {
215 eminus_ok = true;
216 peminus = p.momentum();
217 } else if (pid == PID::POSITRON) {
218 eplus_ok = true;
219 peplus = p.momentum();
220 } else if (pid == PID::MUON) {
221 muminus_ok = true;
222 pmuminus = p.momentum();
223 } else if (pid == PID::ANTIMUON) {
224 muplus_ok = true;
225 pmuplus = p.momentum();
226 } else {
227 throw Error("Unexpected particle type in leading particles FS!");
228 }
229 }
230 // m_ee
231 if (eminus_ok && eplus_ok) {
232 const double m_ee = FourMomentum(peplus + peminus).mass();
233 _hist_mll_ossf_ee->fill(m_ee/GeV);
234 if (epluses.size() == 1 && eminuses.size() == 1)
235 _hist_mll_2_ossf_ee->fill(m_ee/GeV);
236 }
237 // m_mumu
238 if (muminus_ok && muplus_ok) {
239 const double m_mumu = FourMomentum(pmuplus + pmuminus).mass();
240 _hist_mll_ossf_mumu->fill(m_mumu/GeV);
241 if (mupluses.size() == 1 && muminuses.size() == 1)
242 _hist_mll_2_ossf_mumu->fill(m_mumu/GeV);
243 }
244 // m_emu (both configurations)
245 if (eminus_ok && muplus_ok) {
246 const double m_emu = FourMomentum(pmuplus + peminus).mass();
247 _hist_mll_osof_emu->fill(m_emu/GeV);
248 if (mupluses.size() == 1 && eminuses.size() == 1)
249 _hist_mll_2_osof_emu->fill(m_emu/GeV);
250
251 }
252 if (muminus_ok && eplus_ok) {
253 const double m_mue = FourMomentum(peplus + pmuminus).mass();
254 _hist_mll_osof_emu->fill(m_mue/GeV);
255 if (epluses.size() == 1 && muminuses.size() == 1)
256 _hist_mll_2_osof_emu->fill(m_mue/GeV);
257 }
258
259
260 // m_ll plots using *all* electrons, positrons, muons and antimuons
261 // m_ee
262 for (const FourMomentum& peplus : epluses) {
263 for (const FourMomentum& peminus : eminuses) {
264 const double m_ee = FourMomentum(peplus + peminus).mass();
265 _hist_mll_all_ossf_ee->fill(m_ee/GeV);
266 }
267 }
268 // m_mumu
269 for (const FourMomentum& pmuplus : mupluses) {
270 for (const FourMomentum& pmuminus : muminuses) {
271 const double m_mumu = FourMomentum(pmuplus + pmuminus).mass();
272 _hist_mll_all_ossf_mumu->fill(m_mumu/GeV);
273 }
274 }
275 // m_emu (both configurations)
276 for (const FourMomentum& pmuplus : mupluses) {
277 for (const FourMomentum& peminus : eminuses) {
278 const double m_emu = FourMomentum(pmuplus + peminus).mass();
279 _hist_mll_all_osof_emu->fill(m_emu/GeV);
280 }
281 }
282 for (const FourMomentum& peplus : epluses) {
283 for (const FourMomentum& pmuminus : muminuses) {
284 const double m_mue = FourMomentum(peplus + pmuminus).mass();
285 _hist_mll_all_osof_emu->fill(m_mue/GeV);
286 }
287 }
288
289 }
290
291
292 void finalize() {
293 /// @todo Normalisations
294 }
295
296 /// @}
297
298
299 private:
300
301 Histo1DPtr _hist_n_trk, _hist_phi_trk, _hist_eta_trk, _hist_pt_trk;
302 Histo1DPtr _hist_n_jet, _hist_phi_jet, _hist_eta_jet, _hist_pt_jet;
303 Histo1DPtr _hist_n_e, _hist_phi_e, _hist_eta_e, _hist_pt_e;
304 Histo1DPtr _hist_n_mu, _hist_phi_mu, _hist_eta_mu, _hist_pt_mu;
305 Histo1DPtr _hist_n_gamma, _hist_phi_gamma, _hist_eta_gamma, _hist_pt_gamma;
306 Histo1DPtr _hist_n_gammaiso, _hist_phi_gammaiso, _hist_eta_gammaiso, _hist_pt_gammaiso;
307 Histo1DPtr _hist_met;
308 Histo1DPtr _hist_mll_2_ossf_ee, _hist_mll_2_ossf_mumu, _hist_mll_2_osof_emu;
309 Histo1DPtr _hist_mll_ossf_ee, _hist_mll_ossf_mumu, _hist_mll_osof_emu;
310 Histo1DPtr _hist_mll_all_ossf_ee, _hist_mll_all_ossf_mumu, _hist_mll_all_osof_emu;
311 };
312
313
314
315 RIVET_DECLARE_PLUGIN(MC_SUSY);
316
317}
|