Rivet analyses referenceATLAS_2014_I1306615Higgs diphoton events at 8 TeV in ATLASExperiment: ATLAS (LHC) Inspire ID: 1306615 Status: VALIDATED Authors:
Beam energies: (4000.0, 4000.0) GeV Run details:
Measurements of fiducial and differential cross sections are presented for Higgs boson production in proton-proton collisions at a centre-of-mass energy of 8TeV. The analysis is performed in the Higgs diphoton decay channel using 20.3/fb of data recorded by the ATLAS experiment at the CERN Large Hadron Collider. The signal yields are corrected for the effects of detector inefficiency and resolution. Differential cross sections are also presented, as a function of variables related to the diphoton kinematics and the jet activity produced in the Higgs boson events. The observed spectra are statistically limited but broadly in line with the theoretical expectations. Note on reinterpretation: Since the background is subtracted from a fit, any BSM signal would presumably also have beed subtracted. Therefore should not be used to constrain BSM signals unless they comes from a genuine excess of SM Higgs->gammagamma. Source code: ATLAS_2014_I1306615.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/PromptFinalState.hh"
5#include "Rivet/Projections/DressedLeptons.hh"
6#include "Rivet/Projections/VetoedFinalState.hh"
7#include "Rivet/Projections/FastJets.hh"
8
9namespace Rivet {
10
11
12 /// @brief ATLAS H->yy differential cross-sections measurement
13 ///
14 /// @author Michaela Queitsch-Maitland <michaela.queitsch-maitland@cern.ch>
15 //
16 // arXiv: http://arxiv.org/abs/ARXIV:1407.4222
17 // HepData: http://hepdata.cedar.ac.uk/view/ins1306615
18 class ATLAS_2014_I1306615 : public Analysis {
19 public:
20
21 // Constructor
22 ATLAS_2014_I1306615()
23 : Analysis("ATLAS_2014_I1306615")
24 { }
25
26
27 // Book histograms and initialise projections before the run
28 void init() {
29
30 // Final state
31 // All particles within |eta| < 5.0
32 const FinalState FS(Cuts::abseta<5.0);
33 declare(FS,"FS");
34
35 // Project photons with pT > 25 GeV and |eta| < 2.37
36 PromptFinalState ph_FS(Cuts::abseta<2.37 && Cuts::pT>25*GeV && Cuts::pid == PID::PHOTON);
37 declare(ph_FS, "PH_FS");
38
39 // Project photons for dressing
40 FinalState ph_dressing_FS(Cuts::abspid == PID::PHOTON);
41
42 // Project bare electrons
43 PromptFinalState el_bare_FS(Cuts::abseta < 5.0 && Cuts::abspid == PID::ELECTRON);
44
45 // Project dressed electrons with pT > 15 GeV and |eta| < 2.47
46 DressedLeptons el_dressed_FS(ph_dressing_FS, el_bare_FS, 0.1, Cuts::abseta < 2.47 && Cuts::pT > 15*GeV);
47 declare(el_dressed_FS,"EL_DRESSED_FS");
48
49 // Project bare muons
50 PromptFinalState mu_bare_FS(Cuts::abseta < 5.0 && Cuts::abspid == PID::MUON);
51
52 // Project dressed muons with pT > 15 GeV and |eta| < 2.47
53 //DressedLeptons mu_dressed_FS(ph_dressing_FS, mu_bare_FS, 0.1, true, -2.47, 2.47, 15.0*GeV, false);
54 DressedLeptons mu_dressed_FS(ph_dressing_FS, mu_bare_FS, 0.1, Cuts::abseta < 2.47 && Cuts::pT > 15*GeV);
55 declare(mu_dressed_FS,"MU_DRESSED_FS");
56
57 // Final state excluding muons and neutrinos (for jet building and photon isolation)
58 VetoedFinalState veto_mu_nu_FS(FS);
59 veto_mu_nu_FS.vetoNeutrinos();
60 veto_mu_nu_FS.addVetoPairId(PID::MUON);
61 declare(veto_mu_nu_FS, "VETO_MU_NU_FS");
62
63 // Build the anti-kT R=0.4 jets, using FinalState particles (vetoing muons and neutrinos)
64 FastJets jets(veto_mu_nu_FS, FastJets::ANTIKT, 0.4);
65 declare(jets, "JETS");
66
67 // Book histograms
68 // 1D distributions
69 book(_h_pT_yy ,1,1,1);
70 book(_h_y_yy ,2,1,1);
71 book(_h_Njets30 ,3,1,1);
72 book(_h_Njets50 ,4,1,1);
73 book(_h_pT_j1 ,5,1,1);
74 book(_h_y_j1 ,6,1,1);
75 book(_h_HT ,7,1,1);
76 book(_h_pT_j2 ,8,1,1);
77 book(_h_Dy_jj ,9,1,1);
78 book(_h_Dphi_yy_jj ,10,1,1);
79 book(_h_cosTS_CS ,11,1,1);
80 book(_h_cosTS_CS_5bin ,12,1,1);
81 book(_h_Dphi_jj ,13,1,1);
82 book(_h_pTt_yy ,14,1,1);
83 book(_h_Dy_yy ,15,1,1);
84 book(_h_tau_jet ,16,1,1);
85 book(_h_sum_tau_jet ,17,1,1);
86 book(_h_y_j2 ,18,1,1);
87 book(_h_pT_j3 ,19,1,1);
88 book(_h_m_jj ,20,1,1);
89 book(_h_pT_yy_jj ,21,1,1);
90
91 // 2D distributions of cosTS_CS x pT_yy
92 book(_h_cosTS_pTyy_low, 22,1,1);
93 book(_h_cosTS_pTyy_high, 22,1,2);
94 book(_h_cosTS_pTyy_rest, 22,1,3);
95
96 // 2D distributions of Njets x pT_yy
97 book(_h_pTyy_Njets0, 23,1,1);
98 book(_h_pTyy_Njets1, 23,1,2);
99 book(_h_pTyy_Njets2, 23,1,3);
100
101 book(_h_pTj1_excl, 24,1,1);
102
103 // Fiducial regions
104 book(_h_fidXSecs, 30,1,1);
105 }
106
107 // Perform the per-event analysis
108 void analyze(const Event& event) {
109
110 // Get final state particles
111 const Particles& FS_ptcls = apply<FinalState>(event, "FS").particles();
112 const Particles& ptcls_veto_mu_nu = apply<VetoedFinalState>(event, "VETO_MU_NU_FS").particles();
113 const Particles& photons = apply<PromptFinalState>(event, "PH_FS").particlesByPt();
114 vector<DressedLepton> good_el = apply<DressedLeptons>(event, "EL_DRESSED_FS").dressedLeptons();
115 vector<DressedLepton> good_mu = apply<DressedLeptons>(event, "MU_DRESSED_FS").dressedLeptons();
116
117 // For isolation calculation
118 float dR_iso = 0.4;
119 float ETcut_iso = 14.0;
120 FourMomentum ET_iso;
121
122 // Fiducial selection: pT > 25 GeV, |eta| < 2.37 and isolation (in cone deltaR = 0.4) is < 14 GeV
123 Particles fid_photons;
124 for (const Particle& ph : photons) {
125
126 // Calculate isolation
127 ET_iso = - ph.momentum();
128 // Loop over fs truth particles (excluding muons and neutrinos)
129 for (const Particle& p : ptcls_veto_mu_nu) {
130 // Check if the truth particle is in a cone of 0.4
131 if ( deltaR(ph, p) < dR_iso ) ET_iso += p.momentum();
132 }
133
134 // Check isolation
135 if ( ET_iso.Et() > ETcut_iso ) continue;
136
137 // Fill vector of photons passing fiducial selection
138 fid_photons.push_back(ph);
139 }
140
141 if (fid_photons.size() < 2) vetoEvent;
142
143 const FourMomentum y1 = fid_photons[0].momentum();
144 const FourMomentum y2 = fid_photons[1].momentum();
145
146 double m_yy = (y1 + y2).mass();
147
148 // Relative pT cuts
149 if ( y1.pT() < 0.35 * m_yy || y2.pT() < 0.25 * m_yy ) vetoEvent;
150
151 // Mass window cut
152 if ( m_yy < 105 || m_yy > 160 ) vetoEvent;
153
154 // -------------------------------------------- //
155 // Passed diphoton baseline fiducial selection! //
156 // -------------------------------------------- //
157
158 // Muon and Electron selection
159 ifilter_discard(good_mu, [&](const DressedLepton& lep) { return deltaR(lep, y1) < 0.4 || deltaR(lep, y2) < 0.4; });
160 ifilter_discard(good_el, [&](const DressedLepton& lep) { return deltaR(lep, y1) < 0.4 || deltaR(lep, y2) < 0.4; });
161
162 // Find prompt, invisible particles for missing ET calculation
163 // Based on VisibleFinalState projection
164 FourMomentum invisible(0,0,0,0);
165 for (const Particle& p : FS_ptcls) {
166
167 // Veto non-prompt particles (from hadron or tau decay)
168 if ( !p.isPrompt() ) continue;
169 // Charged particles are visible
170 if ( PID::charge3( p.pid() ) != 0 ) continue;
171 // Neutral hadrons are visible
172 if ( PID::isHadron( p.pid() ) ) continue;
173 // Photons are visible
174 if ( p.pid() == PID::PHOTON ) continue;
175 // Gluons are visible (for parton level analyses)
176 if ( p.pid() == PID::GLUON ) continue;
177 // Everything else is invisible
178 invisible += p.momentum();
179 }
180 double MET = invisible.Et();
181
182 // Jet selection
183 // Get jets with pT > 25 GeV and |rapidity| < 4.4
184 //const Jets& jets = apply<FastJets>(event, "JETS").jetsByPt(25.0*GeV, MAXDOUBLE, -4.4, 4.4, RAPIDITY);
185 const Jets& jets = apply<FastJets>(event, "JETS").jetsByPt(Cuts::pT>25*GeV && Cuts::absrap <4.4);
186
187 Jets jets_25, jets_30, jets_50;
188
189 for (const Jet& jet : jets) {
190
191 bool passOverlap = true;
192 // Overlap with leading photons
193 if ( deltaR(y1, jet.momentum()) < 0.4 ) passOverlap = false;
194 if ( deltaR(y2, jet.momentum()) < 0.4 ) passOverlap = false;
195
196 // Overlap with good electrons
197 for (const auto& el : good_el) {
198 if ( deltaR(el, jet) < 0.2 ) passOverlap = false;
199 }
200
201 if ( ! passOverlap ) continue;
202
203 if ( jet.abseta() < 2.4 || ( jet.abseta() > 2.4 && jet.pT() > 30*GeV) ) jets_25 += jet;
204 if ( jet.pT() > 30*GeV ) jets_30 += jet;
205 if ( jet.pT() > 50*GeV ) jets_50 += jet;
206 }
207
208 // Fiducial regions
209 _h_fidXSecs->fill(1);
210 if ( jets_30.size() >= 1 ) _h_fidXSecs->fill(2);
211 if ( jets_30.size() >= 2 ) _h_fidXSecs->fill(3);
212 if ( jets_30.size() >= 3 ) _h_fidXSecs->fill(4);
213 if ( jets_30.size() >= 2 && passVBFCuts(y1 + y2, jets_30.at(0).momentum(), jets_30.at(1).momentum()) ) _h_fidXSecs->fill(5);
214 if ( (good_el.size() + good_mu.size()) > 0 ) _h_fidXSecs->fill(6);
215 if ( MET > 80 ) _h_fidXSecs->fill(7);
216
217 // Fill histograms
218 // Inclusive variables
219 _pT_yy = (y1 + y2).pT();
220 _y_yy = (y1 + y2).absrap();
221 _cosTS_CS = cosTS_CS(y1, y2);
222 _pTt_yy = pTt(y1, y2);
223 _Dy_yy = fabs( deltaRap(y1, y2) );
224
225 _Njets30 = jets_30.size() > 3 ? 3 : jets_30.size();
226 _Njets50 = jets_50.size() > 3 ? 3 : jets_50.size();
227 _h_Njets30->fill(_Njets30);
228 _h_Njets50->fill(_Njets50);
229
230 _pT_j1 = jets_30.size() > 0 ? jets_30[0].momentum().pT() : 0.;
231 _pT_j2 = jets_30.size() > 1 ? jets_30[1].momentum().pT() : 0.;
232 _pT_j3 = jets_30.size() > 2 ? jets_30[2].momentum().pT() : 0.;
233
234 _HT = 0.0;
235 for (const Jet& jet : jets_30) { _HT += jet.pT(); }
236
237 _tau_jet = tau_jet_max(y1 + y2, jets_25);
238 _sum_tau_jet = sum_tau_jet(y1 + y2, jets_25);
239
240 _h_pT_yy ->fill(_pT_yy);
241 _h_y_yy ->fill(_y_yy);
242 _h_pT_j1 ->fill(_pT_j1);
243 _h_cosTS_CS ->fill(_cosTS_CS);
244 _h_cosTS_CS_5bin->fill(_cosTS_CS);
245 _h_HT ->fill(_HT);
246 _h_pTt_yy ->fill(_pTt_yy);
247 _h_Dy_yy ->fill(_Dy_yy);
248 _h_tau_jet ->fill(_tau_jet);
249 _h_sum_tau_jet ->fill(_sum_tau_jet);
250
251 // >=1 jet variables
252 if ( jets_30.size() >= 1 ) {
253 FourMomentum j1 = jets_30[0].momentum();
254 _y_j1 = j1.absrap();
255
256 _h_pT_j2->fill(_pT_j2);
257 _h_y_j1 ->fill(_y_j1);
258 }
259
260 // >=2 jet variables
261 if ( jets_30.size() >= 2 ) {
262 FourMomentum j1 = jets_30[0].momentum();
263 FourMomentum j2 = jets_30[1].momentum();
264
265 _Dy_jj = fabs( deltaRap(j1, j2) );
266 _Dphi_jj = fabs( deltaPhi(j1, j2) );
267 _Dphi_yy_jj = fabs( deltaPhi(y1 + y2, j1 + j2) );
268 _m_jj = (j1 + j2).mass();
269 _pT_yy_jj = (y1 + y2 + j1 + j2).pT();
270 _y_j2 = j2.absrap();
271
272 _h_Dy_jj ->fill(_Dy_jj);
273 _h_Dphi_jj ->fill(_Dphi_jj);
274 _h_Dphi_yy_jj ->fill(_Dphi_yy_jj);
275 _h_m_jj ->fill(_m_jj);
276 _h_pT_yy_jj ->fill(_pT_yy_jj);
277 _h_pT_j3 ->fill(_pT_j3);
278 _h_y_j2 ->fill(_y_j2);
279 }
280
281 // 2D distributions of cosTS_CS x pT_yy
282 if ( _pT_yy < 80 )
283 _h_cosTS_pTyy_low->fill(_cosTS_CS);
284 else if ( _pT_yy > 80 && _pT_yy < 200 )
285 _h_cosTS_pTyy_high->fill(_cosTS_CS);
286 else if ( _pT_yy > 200 )
287 _h_cosTS_pTyy_rest->fill(_cosTS_CS);
288
289 // 2D distributions of pT_yy x Njets
290 if ( _Njets30 == 0 )
291 _h_pTyy_Njets0->fill(_pT_yy);
292 else if ( _Njets30 == 1 )
293 _h_pTyy_Njets1->fill(_pT_yy);
294 else if ( _Njets30 >= 2 )
295 _h_pTyy_Njets2->fill(_pT_yy);
296
297 if ( _Njets30 == 1 ) _h_pTj1_excl->fill(_pT_j1);
298
299 }
300
301 // Normalise histograms after the run
302 void finalize() {
303
304 const double xs = crossSectionPerEvent()/femtobarn;
305
306 scale(_h_pT_yy, xs);
307 scale(_h_y_yy, xs);
308 scale(_h_pT_j1, xs);
309 scale(_h_y_j1, xs);
310 scale(_h_HT, xs);
311 scale(_h_pT_j2, xs);
312 scale(_h_Dy_jj, xs);
313 scale(_h_Dphi_yy_jj, xs);
314 scale(_h_cosTS_CS, xs);
315 scale(_h_cosTS_CS_5bin, xs);
316 scale(_h_Dphi_jj, xs);
317 scale(_h_pTt_yy, xs);
318 scale(_h_Dy_yy, xs);
319 scale(_h_tau_jet, xs);
320 scale(_h_sum_tau_jet, xs);
321 scale(_h_y_j2, xs);
322 scale(_h_pT_j3, xs);
323 scale(_h_m_jj, xs);
324 scale(_h_pT_yy_jj, xs);
325 scale(_h_cosTS_pTyy_low, xs);
326 scale(_h_cosTS_pTyy_high, xs);
327 scale(_h_cosTS_pTyy_rest, xs);
328 scale(_h_pTyy_Njets0, xs);
329 scale(_h_pTyy_Njets1, xs);
330 scale(_h_pTyy_Njets2, xs);
331 scale(_h_pTj1_excl, xs);
332 scale(_h_Njets30, xs);
333 scale(_h_Njets50, xs);
334 scale(_h_fidXSecs, xs);
335 }
336
337 // VBF-enhanced dijet topology selection cuts
338 bool passVBFCuts(const FourMomentum &H, const FourMomentum &j1, const FourMomentum &j2) {
339 return ( fabs(deltaRap(j1, j2)) > 2.8 && (j1 + j2).mass() > 400 && fabs(deltaPhi(H, j1 + j2)) > 2.6 );
340 }
341
342 // Cosine of the decay angle in the Collins-Soper frame
343 double cosTS_CS(const FourMomentum &y1, const FourMomentum &y2) {
344 return fabs( ( (y1.E() + y1.pz())* (y2.E() - y2.pz()) - (y1.E() - y1.pz()) * (y2.E() + y2.pz()) )
345 / ((y1 + y2).mass() * sqrt(pow((y1 + y2).mass(), 2) + pow((y1 + y2).pt(), 2)) ) );
346 }
347
348 // Diphoton pT along thrust axis
349 double pTt(const FourMomentum &y1, const FourMomentum &y2) {
350 return fabs(y1.px() * y2.py() - y2.px() * y1.py()) / (y1 - y2).pT()*2;
351 }
352
353 // Tau of jet (see paper for description)
354 // tau_jet = mT/(2*cosh(y*)), where mT = pT (+) m, and y* = rapidty in Higgs rest frame
355 double tau_jet( const FourMomentum &H, const FourMomentum &jet ) {
356 return sqrt( pow(jet.pT(),2) + pow(jet.mass(),2) ) / (2.0 * cosh( jet.rapidity() - H.rapidity() ) );
357 }
358
359 // Maximal (leading) tau_jet (see paper for description)
360 double tau_jet_max(const FourMomentum &H, const Jets& jets, double tau_jet_cut = 8.) {
361 double max_tj = 0;
362 for (const auto& jet : jets) {
363 FourMomentum j = jet.momentum();
364 if (tau_jet(H, j) > tau_jet_cut) max_tj = max(tau_jet(H, j), max_tj);
365 }
366 return max_tj;
367 }
368
369 // Scalar sum of tau for all jets (see paper for description)
370 double sum_tau_jet(const FourMomentum &H, const Jets& jets, double tau_jet_cut = 8.) {
371 double sum_tj = 0;
372 for (const auto& jet : jets) {
373 FourMomentum j = jet.momentum();
374 if (tau_jet(H, j) > tau_jet_cut) sum_tj += tau_jet(H, j);
375 }
376 return sum_tj;
377 }
378
379 private:
380
381 Histo1DPtr _h_pT_yy;
382 Histo1DPtr _h_y_yy;
383 Histo1DPtr _h_Njets30;
384 Histo1DPtr _h_Njets50;
385 Histo1DPtr _h_pT_j1;
386 Histo1DPtr _h_y_j1;
387 Histo1DPtr _h_HT;
388 Histo1DPtr _h_pT_j2;
389 Histo1DPtr _h_Dy_jj;
390 Histo1DPtr _h_Dphi_yy_jj;
391 Histo1DPtr _h_cosTS_CS;
392 Histo1DPtr _h_cosTS_CS_5bin;
393 Histo1DPtr _h_Dphi_jj;
394 Histo1DPtr _h_pTt_yy;
395 Histo1DPtr _h_Dy_yy;
396 Histo1DPtr _h_tau_jet;
397 Histo1DPtr _h_sum_tau_jet;
398 Histo1DPtr _h_y_j2;
399 Histo1DPtr _h_pT_j3;
400 Histo1DPtr _h_m_jj;
401 Histo1DPtr _h_pT_yy_jj;
402 Histo1DPtr _h_cosTS_pTyy_low;
403 Histo1DPtr _h_cosTS_pTyy_high;
404 Histo1DPtr _h_cosTS_pTyy_rest;
405 Histo1DPtr _h_pTyy_Njets0;
406 Histo1DPtr _h_pTyy_Njets1;
407 Histo1DPtr _h_pTyy_Njets2;
408 Histo1DPtr _h_pTj1_excl;
409 Histo1DPtr _h_fidXSecs;
410
411 int _Njets30;
412 int _Njets50;
413 double _pT_yy;
414 double _y_yy;
415 double _cosTS_CS;
416 double _pT_j1;
417 double _m_jj;
418 double _y_j1;
419 double _HT;
420 double _pT_j2;
421 double _y_j2;
422 double _Dphi_yy_jj;
423 double _pT_yy_jj;
424 double _Dphi_jj;
425 double _Dy_jj;
426 double _pT_j3;
427 double _pTt_yy;
428 double _Dy_yy;
429 double _tau_jet;
430 double _sum_tau_jet;
431 };
432
433
434 RIVET_DECLARE_PLUGIN(ATLAS_2014_I1306615);
435
436}
|