Rivet analyses referenceATLAS_2020_I1790439Higgs to 4-lepton at 13 TeVExperiment: ATLAS (ATLAS) Inspire ID: 1790439 Status: VALIDATED Authors:
Beam energies: (6500.0, 6500.0) GeV Run details:
Inclusive and differential fiducial cross sections of the Higgs boson are measured in the $H \to ZZ^* \to 4\ell$ ($\ell = e, \mu$) decay channel. The results are based on proton-proton collision data produced at the Large Hadron Collider at a centre-of-mass energy of 13 TeV and recorded by the ATLAS detector from 2015 to 2018, equivalent to an integrated luminosity of 139 fb$^{-1}$. The inclusive fiducial cross section for the $H \to ZZ^* \to 4\ell$ process is measured to be $\sigma_\text{fid} = 3.28 \pm 0.32$ fb, in agreement with the Standard Model prediction of $\sigma_\text{fid, SM} = 3.41 \pm 0.18$ fb. Differential fiducial cross sections are measured for a variety of observables which are sensitive to the production and decay of the Higgs boson. All measurements are in agreement with the Standard Model predictions. The results are used to constrain anomalous Higgs boson interactions with Standard Model particles. For fiducial cross sections, different bin represents cross section measurement in different channels. For opposite-flavor cross sections, the branching ratio of Higgs to $4\ell$ of 1.18e$^{-4}$ is used while for same-flavor cross sections, 1.3e$^{-4}$ is used, which takes interference into consideration. For all other cross sections, a combined branching ratio of 1.24e$^{-4}$ is used. In the routine, a Matrix Element calculation based on LO MG is used to assist the pairing of leptons. The Higgs to $ZZ$ branching ratio of 0.02641 is multiplied by the $ZZ$ to four lepton branching ratio 0.004736842 and used to scale all the histograms except for the $\mathtt{xs\_flavor}$ histograms where bins are scaled by the Higgs to four same flavour leptons branching ratio 0.00013 or the Higgs to $ee\mu\mu$ branching ratio 0.000118. Source code: ATLAS_2020_I1790439.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FastJets.hh"
4#include "Rivet/Projections/FinalState.hh"
5#include "Rivet/Projections/PromptFinalState.hh"
6#include "Rivet/Projections/VetoedFinalState.hh"
7#include "Rivet/Projections/LeptonFinder.hh"
8
9namespace Rivet {
10
11
12 /// @brief H->ZZ->4l 13 TeV analysis
13 class ATLAS_2020_I1790439 : public Analysis {
14 public:
15
16 /// Default constructor
17 ATLAS_2020_I1790439()
18 : Analysis("ATLAS_2020_I1790439"),
19 MGME()
20 { }
21
22
23 void init() {
24
25 /// Dressed leptons
26 Cut cut_lep = (Cuts::abseta < 2.7) && (Cuts::pT > 5*GeV);
27 PromptFinalState prompt_photons(Cuts::abspid == PID::PHOTON);
28 PromptFinalState prompt_leptons(Cuts::abspid == PID::MUON || Cuts::abspid == PID::ELECTRON);
29 LeptonFinder dLeptons(prompt_leptons, prompt_photons, 0.1, cut_lep);
30 declare(dLeptons, "AllLeptons");
31
32 /// Jet inputs
33 FinalState fs_jet(Cuts::abseta < 5.0);
34 VetoedFinalState jet_input(fs_jet);
35
36 // reject all leptons dressed with only prompt photons from jet input
37 FinalState leptons(Cuts::abspid == PID::ELECTRON || Cuts::abspid == PID::MUON);
38 LeptonFinder reject_leptons(leptons, prompt_photons, 0.1);
39 jet_input.addVetoOnThisFinalState(reject_leptons);
40
41 // reject prompt invisibles, including from tau decays
42 VetoedFinalState invis_fs_jet(fs_jet);
43 invis_fs_jet.addVetoOnThisFinalState(VisibleFinalState(fs_jet));
44 PromptFinalState invis_pfs_jet = PromptFinalState(invis_fs_jet, TauDecaysAs::PROMPT);
45 jet_input.addVetoOnThisFinalState(invis_pfs_jet);
46
47 // declare jets
48 FastJets jets(jet_input, JetAlg::ANTIKT, 0.4, JetMuons::NONE, JetInvisibles::DECAY);
49 declare(jets, "Jets");
50
51 // Book histograms with continuous (float) binning
52 book(_h["H4l_pt"], 5, 1, 1);
53 book(_h["Z1_m"], 7, 1, 1);
54 book(_h["Z2_m"], 9, 1, 1);
55 book(_h["abshiggs_y"], 11, 1, 1);
56 book(_h["abscthstr"], 13, 1, 1);
57 book(_h["cth1"], 15, 1, 1);
58 book(_h["cth2"], 17, 1, 1);
59 book(_h["phi"], 19, 1, 1);
60 book(_h["phi1"], 21, 1, 1);
61 book(_h["pt4lvy4l_0_0p5"], 51, 1, 1);
62 book(_h["pt4lvy4l_0p5_1"], 51, 1, 3);
63 book(_h["pt4lvy4l_1_1p5"], 51, 1, 5);
64 book(_h["pt4lvy4l_1p5_2p5"], 51, 1, 7);
65 book(_h["pt4lvnjet_0"], 53, 1, 1);
66 book(_h["pt4lvnjet_1"], 53, 1, 3);
67 book(_h["pt4lvnjet_2"], 53, 1, 5);
68 book(_h["pt4lvnjet_3"], 53, 1, 7);
69 book(_h["Z1_m_4l"], 65, 1, 1);
70 book(_h["Z1_m_2l2l"], 66, 1, 1);
71 book(_h["Z2_m_4l"], 68, 1, 1);
72 book(_h["Z2_m_2l2l"], 69, 1, 1);
73 book(_h["phi_4l"], 71, 1, 1);
74 book(_h["phi_2l2l"], 72, 1, 1);
75
76 // Book histograms with discrete (string) binning
77 book(_s["xs_flavour"], 3, 1, 1);
78 book(_s["n_jets"], 23, 1, 1);
79 book(_s["n_jets_incl"], 25, 1, 1);
80 book(_s["n_bjets"], 26, 1, 1);
81 book(_s["jet_pt_leading"], 28, 1, 1);
82 book(_s["jet_pt_subleading"], 30, 1, 1);
83 book(_s["dijet_m"], 32, 1, 1);
84 book(_s["dijet_deltaeta"], 34, 1, 1);
85 book(_s["dijet_deltaphi"], 36, 1, 1);
86 book(_s["pt4lj"], 38, 1, 1);
87 book(_s["pt4ljj"], 40, 1, 1);
88 book(_s["m4lj"], 42, 1, 1);
89 book(_s["m4ljj"], 44, 1, 1);
90 book(_s["m12vsm34"], 46, 1, 1);
91 book(_s["m12vsm34_2l2m"], 48, 1, 1);
92 book(_s["m12vsm34_2l2e"], 49, 1, 1);
93 book(_s["pt4lvpt4lj"], 55, 1, 1);
94 book(_s["pt4ljvm4lj"], 57, 1, 1);
95 book(_s["pt4lvptj0"], 59, 1, 1);
96 book(_s["ptj0vyj0"], 61, 1, 1);
97 book(_s["ptj0vptj1"], 63, 1, 1);
98 book(_s["m12vsm34_4l"], 74, 1, 1);
99 book(_s["m12vsm34_2l2l"], 75, 1, 1);
100
101 // helper axes to map onto discrete binning labels
102 _axisMap["jet_pt_leading"] = YODA::Axis<double>({30., 60., 120., 350.});
103 _axisMap["pt4lj"] = YODA::Axis<double>({0., 60., 120., 350.});
104 _axisMap["m4lj"] = YODA::Axis<double>({120., 180., 220., 300., 400., 600., 2000.});
105 _axisMap["jet_pt_subleading"] = YODA::Axis<double>({30., 60., 120., 120., 350. });
106 _axisMap["dijet_m"] = YODA::Axis<double>({0., 120., 450., 3000.});
107 _axisMap["dijet_deltaeta"] = YODA::Axis<double>({0., 1., 2.5, 9.});
108 _axisMap["dijet_deltaphi"] = YODA::Axis<double>({0.0, 1.571, 3.142, 4.712, 6.283});
109 _axisMap["pt4ljj"] = YODA::Axis<double>({0., 60., 120.});
110 _axisMap["m4ljj"] = YODA::Axis<double>({180., 320., 450., 600., 1000., 2500.});
111 }
112
113 /// Do the per-event analysis
114 void analyze(const Event& e) {
115
116 if (edges.empty()) {
117 for (const string& label : vector<string>{ "xs_flavour", "n_jets", "n_jets_incl", "n_bjets",
118 "jet_pt_leading", "jet_pt_subleading", "dijet_m",
119 "dijet_deltaeta", "dijet_deltaphi", "pt4lj",
120 "pt4ljj", "m4lj", "m4ljj", "m12vsm34", "m12vsm34_2l2m",
121 "m12vsm34_2l2e", "pt4lvpt4lj", "pt4ljvm4lj", "pt4lvptj0",
122 "ptj0vyj0", "ptj0vptj1", "m12vsm34_4l","m12vsm34_2l2l" }) {
123
124 edges[label] = _s[label]->xEdges();
125 }
126 }
127 _s["xs_flavour"]->fill(edges["xs_flavour"][8]);
128
129 const DressedLeptons& all_leps = apply<LeptonFinder>(e, "AllLeptons").dressedLeptons();
130 unsigned int n_parts = all_leps.size();
131 unsigned int n_OSSF_pairs = 0;
132 std::vector<Zstate> dileptons;
133 Particles leptons;
134
135 // Form Z candidate (opposite-sign same-flavour) lepton pairs
136 for (unsigned int i = 0; i < n_parts; i++) {
137 for (unsigned int j = i + 1; j < n_parts; j++) {
138 if (isOSSF(all_leps[i], all_leps[j])){
139 n_OSSF_pairs += 1;
140 // Set positive charge first for later ME calculation
141 if (all_leps[i].charge() > 0) {
142 dileptons.push_back(Zstate( ParticlePair(all_leps[i], all_leps[j]) ));
143 } else {
144 dileptons.push_back(Zstate( ParticlePair(all_leps[j], all_leps[i]) ));
145 }
146 }
147 }
148 }
149
150 // At least two pairs required to select ZZ->llll final state
151 if (n_OSSF_pairs < 2) vetoEvent;
152
153
154 // Form the quadruplet of two lepon pairs passing kinematics cuts
155 std::vector<Quadruplet> quadruplets;
156 for (unsigned int i = 0; i < dileptons.size(); i++) {
157 for (unsigned int j = i+1; j < dileptons.size(); j++) {
158 // Only use unique leptons
159 if (isSame( dileptons[i].first , dileptons[j].first )) continue;
160 if (isSame( dileptons[i].first , dileptons[j].second )) continue;
161 if (isSame( dileptons[i].second, dileptons[j].first )) continue;
162 if (isSame( dileptons[i].second, dileptons[j].second )) continue;
163
164 leptons.clear();
165 leptons.push_back( dileptons[i].first );
166 leptons.push_back( dileptons[i].second );
167 leptons.push_back( dileptons[j].first );
168 leptons.push_back( dileptons[j].second );
169 leptons = sortByPt(leptons);
170
171 // Apply kinematic cuts
172 if ( leptons[0].pt() < 20*GeV) continue;
173 if ( leptons[1].pt() < 15*GeV) continue;
174 if ( leptons[2].pt() < 10*GeV) continue;
175
176 // Form the quad with pair closest to Z pole first
177 if (dileptons[i].Zdist() < dileptons[j].Zdist()) {
178 quadruplets.push_back(Quadruplet(dileptons[i], dileptons[j]));
179 } else {
180 quadruplets.push_back(Quadruplet(dileptons[j], dileptons[i]));
181 }
182 }
183 }
184
185 // Veto if no quad passes kinematic selection
186 unsigned int n_quads = quadruplets.size();
187 if(n_quads == 0) vetoEvent;
188
189 // To resolve ambiguities in lepton pairing order quads by channel priority first, then m12 - mz and m34 - mz
190 // The first in every channel is considered nominal
191 std::sort(quadruplets.begin(), quadruplets.end(),
192 [](const Quadruplet & q1, const Quadruplet & q2) {
193 if (q1.ch_priority == q2.ch_priority) {
194 // if rarely, Z1 the same distance from the Z pole, compare Z2
195 if (fabs( q1.Z1().Zdist() - q2.Z1().Zdist() ) < 1.e-5)
196 return q1.Z2().Zdist() < q2.Z2().Zdist();
197 else
198 return q1.Z1().Zdist() < q2.Z1().Zdist();
199 } else
200 return q1.ch_priority < q2.ch_priority;
201 });
202
203
204 // Select the best quad
205 Particles leptons_sel4l;
206 Quadruplet quadSel;
207 float MEHZZ_max = -999.;
208 int prevQuadType = -999.;
209 bool atleastonequadpassed = false;
210 bool isNominalQuad = false;
211 bool extraLep = false;
212
213 for (unsigned int iquad = 0; iquad < n_quads; iquad++) {
214
215 // Veto event if nominal quad was not selected in 4 lepton case
216 if (n_parts == 4 && iquad > 0) vetoEvent;
217
218 int quadType = (int) quadruplets[iquad].type();
219
220 if (quadType != prevQuadType) {
221 isNominalQuad = true;
222 } else {
223 isNominalQuad = false;
224 }
225 prevQuadType = quadType;
226
227 Quadruplet & quad = quadruplets[iquad];
228
229 // Z invariant mass requirements
230 if (!(inRange(quad.Z1().mom().mass(), 50*GeV, 106*GeV))) continue;
231 if (!(inRange(quad.Z2().mom().mass(), 12*GeV, 115*GeV))) continue;
232
233 // Lepton separation and J/Psi veto
234 bool b_pass_leptonseparation = true;
235 bool b_pass_jpsi = true;
236 leptons_sel4l.clear();
237 leptons_sel4l.push_back(quad.Z1().first);
238 leptons_sel4l.push_back(quad.Z1().second);
239 leptons_sel4l.push_back(quad.Z2().first);
240 leptons_sel4l.push_back(quad.Z2().second);
241
242 for (unsigned int i = 0; i < 4; i++) {
243 for (unsigned int j = i+1; j < 4; j++) {
244 if ( deltaR( leptons_sel4l[i], leptons_sel4l[j]) < 0.1) b_pass_leptonseparation = false;
245 if ( isOSSF(leptons_sel4l[i], leptons_sel4l[j])
246 && (leptons_sel4l[i].mom() + leptons_sel4l[j].mom()).mass() <= 5.*GeV) b_pass_jpsi = false;
247 }
248 }
249 if(b_pass_leptonseparation == false || b_pass_jpsi == false) continue;
250
251 // Only consider the event if at least one nominal quadruplet passes cuts
252 if (isNominalQuad) atleastonequadpassed = true;
253
254 // Direct selection for case with only 4 prompt leptons
255 if (n_parts == 4 ) {
256 quadSel = quad;
257 break;
258 }
259
260 // In cases with extra leptons meeting further requirements use max ME to select quad
261 float MEHZZ;
262 if (quadType == 0 || quadType == 1) MEHZZ = MGME.Compute(quadruplets[iquad]) / 2;
263 else MEHZZ = MGME.Compute(quadruplets[iquad]);
264
265 // Check leptons other than the ones from the channel's nominal quadruplet
266 // and don't need to recheck extraLep for a second channel
267 if(isNominalQuad && !extraLep) {
268 for (const Particle &lep: all_leps) {
269 bool lep_in_quad = false;
270 bool lep_too_close = false;
271
272 // pT requirement
273 if (lep.pt() < 12*GeV) continue;
274
275 // skip if lepton included in quad or not isolated from quad leptons
276 for (unsigned int i = 0; i < 4; i++) {
277 if (isSame(lep, leptons_sel4l[i])) lep_in_quad = true ;
278 if (deltaR(lep, leptons_sel4l[i]) < 0.1) lep_too_close = true ;
279 }
280 if (lep_in_quad || lep_too_close) continue;
281
282 extraLep = true;
283 break;
284 }
285
286 if(!extraLep) {
287 // In case of no suitable extra leptons select the quad directly
288 quadSel = quad;
289 break;
290 }
291 }
292
293 // Use ME to select the quad when there are extra leptons
294 if (MEHZZ > MEHZZ_max) {
295 quadSel = quad;
296 MEHZZ_max = MEHZZ;
297 }
298 }
299
300 if(!atleastonequadpassed) vetoEvent;
301
302 // Veto if quad not in Higgs mass window
303 FourMomentum Higgs = quadSel.mom();
304 double H4l_mass = Higgs.mass();
305 if (!(inRange(H4l_mass, 105.*GeV, 160.*GeV))) vetoEvent;
306
307 // Higgs observables
308 double H4l_pt = Higgs.pt();
309 double H4l_rapidity = Higgs.absrapidity();
310 LorentzTransform HRF_boost = LorentzTransform::mkFrameTransformFromBeta(Higgs.betaVec());
311
312 FourMomentum Z1_in_HRF = HRF_boost.transform( quadSel.Z1().mom() );
313 FourMomentum Z2_in_HRF = HRF_boost.transform( quadSel.Z2().mom() );
314 double H4l_costheta = fabs(cos( Z1_in_HRF.theta()));
315 double H4l_m12 = quadSel.Z1().mom().mass();
316 double H4l_m34 = quadSel.Z2().mom().mass();
317
318 FourMomentum v11_HRF = HRF_boost.transform( quadSel.Z1().second.mom() );
319 FourMomentum v12_HRF = HRF_boost.transform( quadSel.Z1().first.mom() );
320 FourMomentum v21_HRF = HRF_boost.transform( quadSel.Z2().second.mom() );
321 FourMomentum v22_HRF = HRF_boost.transform( quadSel.Z2().first.mom() );
322
323 Vector3 v11 = v11_HRF.p3();
324 Vector3 v12 = v12_HRF.p3();
325 Vector3 v21 = v21_HRF.p3();
326 Vector3 v22 = v22_HRF.p3();
327 Vector3 nz(0, 0, 1);
328 Vector3 qz1 = Z1_in_HRF.p3();
329 Vector3 n1p = v11.cross(v12).unit();
330 Vector3 n2p = v21.cross(v22).unit();
331 Vector3 nscp = nz.cross(qz1).unit();
332
333 double H4l_Phi = qz1.dot(n1p.cross(n2p)) / fabs(qz1.dot(n1p.cross(n2p))) * acos(- n1p.dot(n2p));
334 double H4l_Phi1 = qz1.dot(n1p.cross(nscp)) / fabs(qz1.dot(n1p.cross(nscp))) * acos( n1p.dot(nscp));
335
336 LorentzTransform Z1RF_boost = LorentzTransform::mkFrameTransformFromBeta(Z1_in_HRF.betaVec());
337 LorentzTransform Z2RF_boost = LorentzTransform::mkFrameTransformFromBeta(Z2_in_HRF.betaVec());
338
339 FourMomentum Z1_in_Z2RF = Z2RF_boost.transform( Z1_in_HRF );
340 FourMomentum Z2_in_Z1RF = Z1RF_boost.transform( Z2_in_HRF );
341
342 Vector3 Z1_p3 = Z1_in_Z2RF.p3();
343 Vector3 Z2_p3 = Z2_in_Z1RF.p3();
344
345 FourMomentum n_Z1 = Z1RF_boost.transform( v11_HRF );
346 FourMomentum n_Z2 = Z2RF_boost.transform( v21_HRF );
347
348 // angle is negative with its Z
349 double H4l_cth1 = - cos( n_Z1.p3().angle(Z2_p3));
350 double H4l_cth2 = - cos( n_Z2.p3().angle(Z1_p3));
351
352
353 // Jet observables
354 Jets jets = apply<FastJets>(e, "Jets").jetsByPt(Cuts::pT > 30*GeV && Cuts::absrap < 4.4);
355
356 // discard jets which overlap leptons
357 for (const Particle& l: all_leps) idiscard(jets, deltaRLess(l, 0.1));
358
359 // collect b-tagged jets
360 const Jets b_jets_btagged = select(jets, hasBTag(Cuts::pT>5*GeV));
361
362 unsigned int n_bjets = b_jets_btagged.size();
363 unsigned int n_jets = jets.size();
364 double leading_jet_pt = (n_jets > 0 ? jets[0].pT() : 0);
365 double leading_jet_y = (n_jets > 0 ? jets[0].absrapidity() : 0);
366 double subleading_jet_pt = (n_jets > 1 ? jets[1].pT() : 0);
367
368 FourMomentum Dijets = {0,0,0,0};
369 if(n_jets > 1) Dijets = jets[0].mom() + jets[1].mom();
370 double mjj = (n_jets > 1 ? Dijets.mass() : -999);
371
372 double dphijj = -1;
373 if(n_jets > 1){
374 if(jets[0].eta() > jets[1].eta()) dphijj = jets[0].phi() - jets[1].phi();
375 else dphijj = jets[1].phi() - jets[0].phi();
376 if(dphijj < 0) dphijj = dphijj + TWOPI;
377 }
378
379 double detajj = (n_jets > 1 ? fabs(jets[0].eta() - jets[1].eta()): -1);
380
381 FourMomentum m4lj, m4ljj;
382 if (n_jets > 0) m4lj = Higgs + jets[0];
383 if (n_jets > 1) m4ljj = Higgs + Dijets;
384
385 double H4l_m4lj = m4lj.mass();
386 double H4l_m4ljj = m4ljj.mass();
387 double H4l_pt4lj = m4lj.pt();
388 double H4l_pt4ljj = m4ljj.pt();
389
390 // Fill histograms
391 // Branching ratios for 4l and 2l2l channels
392 double BR_SF = 0.00013;
393 double BR_OF = 0.000118;
394 if (inRange(H4l_mass, 115.*GeV, 130.*GeV)){
395 if (quadSel.type() == Quadruplet::FlavCombi::mm
396 || quadSel.type() == Quadruplet::FlavCombi::ee ) {
397 _s["xs_flavour"]->fill(edges["xs_flavour"][(int)quadSel.type()], BR_SF);
398 _s["xs_flavour"]->fill(edges["xs_flavour"][4], BR_SF);
399 }
400 else if (quadSel.type() == Quadruplet::FlavCombi::em
401 || quadSel.type() == Quadruplet::FlavCombi::me ) {
402 _s["xs_flavour"]->fill(edges["xs_flavour"][(int)quadSel.type()], BR_OF);
403 _s["xs_flavour"]->fill(edges["xs_flavour"][5], BR_OF);
404 }
405 _s["xs_flavour"]->fill(edges["xs_flavour"][6], Br);
406 _s["xs_flavour"]->fill(edges["xs_flavour"][7], Br);
407 }
408
409 // Higgs variables
410 for(const auto & p: std::map<std::string, double>{
411 {"H4l_pt", H4l_pt},
412 {"Z1_m", H4l_m12},
413 {"Z2_m", H4l_m34},
414 {"abshiggs_y", H4l_rapidity},
415 {"abscthstr", H4l_costheta},
416 {"cth1", H4l_cth1},
417 {"cth2", H4l_cth2},
418 {"phi", H4l_Phi},
419 {"phi1", H4l_Phi1}}) {
420 _h[ p.first ]->fill(p.second);
421 }
422
423 // Jet variables
424 discreteFill("n_jets", min((int)n_jets, 3));
425
426 discreteFill("n_jets_incl", 0);
427 if (n_jets >= 1) discreteFill("n_jets_incl", 1);
428 if (n_jets >= 2) discreteFill("n_jets_incl", 2);
429 if (n_jets >= 3) discreteFill("n_jets_incl", 3);
430
431 if (n_jets == 0) discreteFill("n_bjets", 0);
432 else if (n_bjets == 0) discreteFill("n_bjets", 1);
433 else if (n_bjets >= 1) discreteFill("n_bjets", 2);
434
435 if (n_jets == 0) {
436 discreteFill("jet_pt_leading", 0);
437 discreteFill("pt4lj", 0);
438 discreteFill("m4lj", 0);
439 }
440 else if(n_jets >= 1){
441 discreteFill("jet_pt_leading", leading_jet_pt);
442 discreteFill("pt4lj", H4l_pt4lj);
443 discreteFill("m4lj", H4l_m4lj);
444 }
445
446 if (n_jets < 2) {
447 discreteFill("jet_pt_subleading", 0);
448 discreteFill("dijet_m", 0);
449 discreteFill("dijet_deltaeta", 0);
450 discreteFill("dijet_deltaphi", 0);
451 discreteFill("pt4ljj", 0);
452 discreteFill("m4ljj", 0);
453 }
454 else if (n_jets >= 2) {
455 discreteFill("jet_pt_subleading", subleading_jet_pt);
456 discreteFill("dijet_m", mjj);
457 discreteFill("dijet_deltaeta", detajj);
458 discreteFill("dijet_deltaphi", dphijj);
459 discreteFill("pt4ljj", H4l_pt4ljj);
460 discreteFill("m4ljj", H4l_m4ljj);
461 }
462
463 // m12 vs m34 (all channels)
464 if(H4l_m12 < 82 && H4l_m34 < 32) discreteFill("m12vsm34", 0);
465 else if(H4l_m12 < 74 && H4l_m34 > 32) discreteFill("m12vsm34", 1);
466 else if(H4l_m12 > 74 && H4l_m34 > 32) discreteFill("m12vsm34", 2);
467 else if(H4l_m12 > 82 && H4l_m34 < 32 && H4l_m34 > 24) discreteFill("m12vsm34", 3);
468 else if(H4l_m12 > 82 && H4l_m34 < 24) discreteFill("m12vsm34", 4);
469
470 if (quadSel.type() == Quadruplet::FlavCombi::em
471 || quadSel.type() == Quadruplet::FlavCombi::mm ){
472 if(H4l_m12 < 82 && H4l_m34 < 32) discreteFill("m12vsm34_2l2m", 0);
473 else if(H4l_m12 < 74 && H4l_m34 > 32) discreteFill("m12vsm34_2l2m", 1);
474 else if(H4l_m12 > 74 && H4l_m34 > 32) discreteFill("m12vsm34_2l2m", 2);
475 else if(H4l_m12 > 82 && H4l_m34 < 32 && H4l_m34 > 24) discreteFill("m12vsm34_2l2m", 3);
476 else if(H4l_m12 > 82 && H4l_m34 < 24) discreteFill("m12vsm34_2l2m", 4);
477 }
478 else if (quadSel.type() == Quadruplet::FlavCombi::me
479 || quadSel.type() == Quadruplet::FlavCombi::ee ){
480 if(H4l_m12 < 82 && H4l_m34 < 32) discreteFill("m12vsm34_2l2e", 0);
481 else if(H4l_m12 < 74 && H4l_m34 > 32) discreteFill("m12vsm34_2l2e", 1);
482 else if(H4l_m12 > 74 && H4l_m34 > 32) discreteFill("m12vsm34_2l2e", 2);
483 else if(H4l_m12 > 82 && H4l_m34 < 32 && H4l_m34 > 24) discreteFill("m12vsm34_2l2e", 3);
484 else if(H4l_m12 > 82 && H4l_m34 < 24) discreteFill("m12vsm34_2l2e", 4);
485 }
486
487 // m12 vs m34 (4l channels only)
488 if (quadSel.type() == Quadruplet::FlavCombi::ee
489 || quadSel.type() == Quadruplet::FlavCombi::mm ){
490 if(H4l_m12 < 82 && H4l_m34 < 32) discreteFill("m12vsm34_4l", 0);
491 else if(H4l_m12 < 74 && H4l_m34 > 32) discreteFill("m12vsm34_4l", 1);
492 else if(H4l_m12 > 74 && H4l_m34 > 32) discreteFill("m12vsm34_4l", 2);
493 else if(H4l_m12 > 82 && H4l_m34 < 32 && H4l_m34 > 24) discreteFill("m12vsm34_4l", 3);
494 else if(H4l_m12 > 82 && H4l_m34 < 24) discreteFill("m12vsm34_4l", 4);
495 _h["Z1_m_4l"]->fill(H4l_m12);
496 _h["Z2_m_4l"]->fill(H4l_m34);
497 _h["phi_4l"]->fill(H4l_Phi);
498 }
499
500 // m12 vs m34 (2l2l channels only)
501 if (quadSel.type() == Quadruplet::FlavCombi::me
502 || quadSel.type() == Quadruplet::FlavCombi::em ){
503 if(H4l_m12 < 82 && H4l_m34 < 32) discreteFill("m12vsm34_2l2l", 0);
504 else if(H4l_m12 < 74 && H4l_m34 > 32) discreteFill("m12vsm34_2l2l", 1);
505 else if(H4l_m12 > 74 && H4l_m34 > 32) discreteFill("m12vsm34_2l2l", 2);
506 else if(H4l_m12 > 82 && H4l_m34 < 32 && H4l_m34 > 24) discreteFill("m12vsm34_2l2l", 3);
507 else if(H4l_m12 > 82 && H4l_m34 < 24) discreteFill("m12vsm34_2l2l", 4);
508 _h["Z1_m_2l2l"]->fill(H4l_m12);
509 _h["Z2_m_2l2l"]->fill(H4l_m34);
510 _h["phi_2l2l"]->fill(H4l_Phi);
511 }
512
513 // 2d differential variables
514 if (0 < H4l_rapidity && H4l_rapidity < 0.5) _h["pt4lvy4l_0_0p5"]->fill(H4l_pt);
515 else if (0.5 < H4l_rapidity && H4l_rapidity < 1) _h["pt4lvy4l_0p5_1"]->fill(H4l_pt);
516 else if (1 < H4l_rapidity && H4l_rapidity < 1.5) _h["pt4lvy4l_1_1p5"]->fill(H4l_pt);
517 else if (1.5 < H4l_rapidity && H4l_rapidity < 2.5) _h["pt4lvy4l_1p5_2p5"]->fill(H4l_pt);
518
519 if (n_jets == 0) _h["pt4lvnjet_0"]->fill(H4l_pt);
520 else if (n_jets == 1) _h["pt4lvnjet_1"]->fill(H4l_pt);
521 else if (n_jets == 2) _h["pt4lvnjet_2"]->fill(H4l_pt);
522 else if (n_jets > 2) _h["pt4lvnjet_3"]->fill(H4l_pt);
523
524 if (n_jets == 0) {
525 discreteFill("pt4lvptj0", 0);
526 discreteFill("pt4lvpt4lj", 0);
527 discreteFill("pt4ljvm4lj", 0);
528 discreteFill("ptj0vptj1", 0);
529 discreteFill("ptj0vyj0", 0);
530 }
531 else {
532
533 if ( 0 < H4l_pt4lj && H4l_pt4lj < 60 && 0 < H4l_pt && H4l_pt < 120) discreteFill("pt4lvpt4lj", 1);
534 else if (0 < H4l_pt4lj && H4l_pt4lj < 60 && 120 < H4l_pt && H4l_pt < 350) discreteFill("pt4lvpt4lj", 2);
535 else if (60 < H4l_pt4lj && H4l_pt4lj < 350 && 0 < H4l_pt && H4l_pt < 120) discreteFill("pt4lvpt4lj", 3);
536 else if (60 < H4l_pt4lj && H4l_pt4lj < 350 && 120 < H4l_pt && H4l_pt < 350) discreteFill("pt4lvpt4lj", 4);
537
538 if ( 120 < H4l_m4lj && H4l_m4lj < 220 && 0 < H4l_pt4lj && H4l_pt4lj < 350) discreteFill("pt4ljvm4lj", 1);
539 else if (220 < H4l_m4lj && H4l_m4lj < 350 && 0 < H4l_pt4lj && H4l_pt4lj < 60) discreteFill("pt4ljvm4lj", 2);
540 else if (220 < H4l_m4lj && H4l_m4lj < 350 && 60 < H4l_pt4lj && H4l_pt4lj < 350) discreteFill("pt4ljvm4lj", 3);
541 else if (350 < H4l_m4lj && H4l_m4lj < 2000 && 0 < H4l_pt4lj && H4l_pt4lj < 350) discreteFill("pt4ljvm4lj", 4);
542
543 if ( 30 < leading_jet_pt && leading_jet_pt < 60 && 0 < H4l_pt && H4l_pt < 80) discreteFill("pt4lvptj0", 1);
544 else if (30 < leading_jet_pt && leading_jet_pt < 60 && 80 < H4l_pt && H4l_pt < 350) discreteFill("pt4lvptj0", 2);
545 else if (60 < leading_jet_pt && leading_jet_pt < 120 && 0 < H4l_pt && H4l_pt < 120) discreteFill("pt4lvptj0", 3);
546 else if (60 < leading_jet_pt && leading_jet_pt < 120 && 120 < H4l_pt && H4l_pt < 350) discreteFill("pt4lvptj0", 4);
547 else if (120 < leading_jet_pt && leading_jet_pt < 350 && 0 < H4l_pt && H4l_pt < 120) discreteFill("pt4lvptj0", 5);
548 else if (120 < leading_jet_pt && leading_jet_pt < 350 && 120 < H4l_pt && H4l_pt < 350) discreteFill("pt4lvptj0", 6);
549
550 if (30 < leading_jet_pt && leading_jet_pt < 120 && 0 < leading_jet_y && leading_jet_y < 0.8) {
551 discreteFill("ptj0vyj0", 1);
552 }
553 else if (30 < leading_jet_pt && leading_jet_pt < 120 && 0.8 < leading_jet_y && leading_jet_y < 1.7) {
554 discreteFill("ptj0vyj0", 2);
555 }
556 else if (30 < leading_jet_pt && leading_jet_pt < 120 && 1.7 < leading_jet_y) {
557 discreteFill("ptj0vyj0", 3);
558 }
559 else if (120 < leading_jet_pt && leading_jet_pt < 350 && 0 < leading_jet_y && leading_jet_y < 1.7) {
560 discreteFill("ptj0vyj0", 4);
561 }
562 else if (120 < leading_jet_pt && leading_jet_pt < 350 && 1.7 < leading_jet_y) {
563 discreteFill("ptj0vyj0", 5);
564 }
565
566 if ( n_jets == 1 && 30 < leading_jet_pt && leading_jet_pt < 60) discreteFill("ptj0vptj1", 1);
567 else if (n_jets == 1 && 60 < leading_jet_pt && leading_jet_pt < 350) discreteFill("ptj0vptj1", 2);
568 else if (30 < leading_jet_pt && leading_jet_pt < 60 && 30 < subleading_jet_pt && subleading_jet_pt < 60) {
569 discreteFill("ptj0vptj1", 3);
570 }
571 else if (60 < leading_jet_pt && leading_jet_pt < 350 && 30 < subleading_jet_pt && subleading_jet_pt < 60) {
572 discreteFill("ptj0vptj1", 4);
573 }
574 else if (60 < leading_jet_pt && leading_jet_pt < 350 && 60 < subleading_jet_pt && subleading_jet_pt < 350) {
575 discreteFill("ptj0vptj1", 5);
576 }
577 }
578 }
579
580 template<typename T>
581 void discreteFill(const string& label, const T coord) {
582 if constexpr( std::is_floating_point<T>::value) {
583 discreteFill(label, _axisMap[label].index(coord));
584 }
585 else {
586 _s[label]->fill(edges[label][coord]);
587 }
588 }
589
590 void finalize() {
591
592 const double sf = crossSection() / femtobarn / sumOfWeights();
593 scale(_h, sf * Br);
594 for (auto& hist : _s) {
595 if (hist.first == "xs_flavour") {
596 scale(hist.second, sf);
597 }
598 else {
599 scale(hist.second, sf * Br);
600 }
601 }
602 }
603
604 private:
605
606 // Br(H-->ZZ) * BR(ZZ-->4l)
607 const double Br = 0.02641 * 0.004736842;
608
609 map<string, Histo1DPtr> _h;
610 map<string, BinnedHistoPtr<string>> _s;
611 map<string, vector<string>> edges;
612 map<string, YODA::Axis<double>> _axisMap;
613
614 /// Generic Z candidate
615 struct Zstate : public ParticlePair {
616 Zstate() { }
617 Zstate(ParticlePair _particlepair) : ParticlePair(_particlepair) { }
618 FourMomentum mom() const { return first.momentum() + second.momentum(); }
619 double Zdist() const { return fabs(mom().mass() - 91.1876*GeV); }
620 int flavour() const { return first.abspid(); }
621 };
622
623 /// Generic quadruplet
624 struct Quadruplet {
625
626 // find out which type it is: 4mu = 0, 4e = 1, 2mu2e = 2, 2e2mu = 3 (mm, ee, me, em)
627 // channel priority is 4m, 2e2m, 2m2e, 4e
628 enum class FlavCombi { mm=0, ee, me, em, undefined };
629
630 Quadruplet() { }
631 Quadruplet(Zstate z1, Zstate z2) : _z1(z1), _z2(z2) {
632 if ( _z1.flavour() == 13 && _z2.flavour() == 13) { _type = FlavCombi::mm; ch_priority = 0;}
633 else if (_z1.flavour() == 11 && _z2.flavour() == 11) { _type = FlavCombi::ee; ch_priority = 3;}
634 else if (_z1.flavour() == 13 && _z2.flavour() == 11) { _type = FlavCombi::me; ch_priority = 2;}
635 else if (_z1.flavour() == 11 && _z2.flavour() == 13) { _type = FlavCombi::em; ch_priority = 1;}
636 else {_type = FlavCombi::undefined;}
637 }
638
639 Quadruplet(Quadruplet const & quad) :
640 _z1(quad._z1),
641 _z2(quad._z2),
642 _type(quad._type),
643 ch_priority(quad.ch_priority) {}
644
645 Zstate _z1, _z2;
646 FlavCombi _type;
647 int ch_priority;
648
649 const Zstate& Z1() const { return _z1; }
650 const Zstate& Z2() const { return _z2; }
651 FourMomentum mom() const { return _z1.mom() + _z2.mom(); }
652 FlavCombi type() const {return _type; }
653 };
654
655
656 // save and calculate parameters
657 class Parameters_heft {
658
659 public:
660
661 // Model parameters independent of aS
662 double mdl_WH, mdl_WZ,
663 aS, mdl_Gf, aEWM1, mdl_MH, mdl_MZ, mdl_MTA, mdl_MT, mdl_MB,
664 mdl_MP, mdl_conjg__CKM3x3, mdl_CKM3x3, mdl_MZ__exp__2, mdl_MZ__exp__4,
665 mdl_MH__exp__4, mdl_MT__exp__4, mdl_MH__exp__2,
666 mdl_MT__exp__2,
667 mdl_MH__exp__6, mdl_MT__exp__6, mdl_aEW, mdl_MW, mdl_ee,
668 mdl_MW__exp__2, mdl_sw2, mdl_cw,
669 mdl_sw,
670 mdl_v, mdl_ee__exp__2;
671 std::complex<double> mdl_complexi;
672 // Model parameters dependent on aS
673 double mdl_GH ;
674 // Model couplings independent of aS
675 std::complex<double> GC_40, GC_54, GC_73;
676 // Model couplings dependent on aS
677 std::complex<double> GC_13;
678
679 // Set parameters and couplings that are unchanged during the run
680 Parameters_heft() {
681 mdl_WH = 6.382339e-03;
682 mdl_WZ = 2.441404e+00;
683 aS = 1.180000e-01;
684 mdl_Gf = 1.166390e-05;
685 aEWM1 = 1.325070e+02;
686 mdl_MH = 1.250000e+02;
687 mdl_MZ = 9.118800e+01;
688 mdl_MT = 1.730000e+02;
689 mdl_complexi = std::complex<double> (0., 1.);
690 mdl_MZ__exp__2 = pow(mdl_MZ, 2.);
691 mdl_MZ__exp__4 = pow(mdl_MZ, 4.);
692 mdl_MH__exp__2 = pow(mdl_MH, 2.);
693 mdl_MT__exp__4 = pow(mdl_MT, 4.);
694 mdl_MH__exp__4 = pow(mdl_MH, 4.);
695 mdl_MT__exp__2 = pow(mdl_MT, 2.);
696 mdl_MH__exp__6 = pow(mdl_MH, 6.);
697 mdl_MT__exp__6 = pow(mdl_MT, 6.);
698 mdl_aEW = 1./aEWM1;
699 mdl_MW = sqrt(mdl_MZ__exp__2/2. + sqrt(mdl_MZ__exp__4/4. - (mdl_aEW * M_PI * mdl_MZ__exp__2)/(mdl_Gf * sqrt(2.)))); mdl_ee = 2. * sqrt(mdl_aEW) * sqrt(M_PI);
700 mdl_MW__exp__2 = pow(mdl_MW, 2.);
701 mdl_sw2 = 1. - mdl_MW__exp__2/mdl_MZ__exp__2;
702 mdl_cw = sqrt(1. - mdl_sw2);
703 mdl_sw = sqrt(mdl_sw2);
704 mdl_v = (2. * mdl_MW * mdl_sw)/mdl_ee;
705 mdl_ee__exp__2 = pow(mdl_ee, 2.);
706 GC_40 = -(mdl_ee * mdl_complexi * mdl_cw)/(2. * mdl_sw);
707 GC_54 = (mdl_ee * mdl_complexi * mdl_sw)/(2. * mdl_cw);
708 GC_73 = mdl_ee__exp__2 * mdl_complexi * mdl_v + ((1. - mdl_sw2) * mdl_ee__exp__2 * mdl_complexi * mdl_v)/(2. * mdl_sw2) +
709 (mdl_ee__exp__2 * mdl_complexi * mdl_sw2 * mdl_v)/(2. * (1. - mdl_sw2));
710 }
711
712 // Set Mass
713 void set4lepMass(double m_m4l){
714 mdl_MH = m_m4l;
715 mdl_WH = setWidth(m_m4l);
716 mdl_MH__exp__2 = pow(mdl_MH, 2.);
717 mdl_MH__exp__4 = pow(mdl_MH, 4.);
718 mdl_MH__exp__6 = pow(mdl_MH, 6.);
719 mdl_GH = -(4 * aS * M_PI * (1. + (13. * mdl_MH__exp__6)/(16800. * mdl_MT__exp__6) + mdl_MH__exp__4/(168. * mdl_MT__exp__4) +
720 (7. * mdl_MH__exp__2)/(120. * mdl_MT__exp__2)))/(12. * pow(M_PI, 2.) * mdl_v);
721 GC_13 = -(mdl_complexi * mdl_GH);
722 }
723
724 private:
725
726 // Set Width
727 long double setWidth(double m_m4l){
728 long double Higgs_width_Poly_Fit_Zone1_coeff0 = -1.450308902710193E+03;
729 long double Higgs_width_Poly_Fit_Zone1_coeff1 = 1.129291251156317E+02;
730 long double Higgs_width_Poly_Fit_Zone1_coeff2 = -3.893063071316150E+00;
731 long double Higgs_width_Poly_Fit_Zone1_coeff3 = 7.798666884832531E-02;
732 long double Higgs_width_Poly_Fit_Zone1_coeff4 = -1.000455877406390E-03;
733 long double Higgs_width_Poly_Fit_Zone1_coeff5 = 8.523735379647125E-06;
734 long double Higgs_width_Poly_Fit_Zone1_coeff6 = -4.823164754652171E-08;
735 long double Higgs_width_Poly_Fit_Zone1_coeff7 = 1.747954506786346E-10;
736 long double Higgs_width_Poly_Fit_Zone1_coeff8 = -3.681723572169337E-13;
737 long double Higgs_width_Poly_Fit_Zone1_coeff9 = 3.434207075968898E-16;
738
739 long double Higgs_width_Poly_Fit_Zone2_coeff0 = 2.563291882845993E+02;
740 long double Higgs_width_Poly_Fit_Zone2_coeff1 = -1.037082025855304E+01;
741 long double Higgs_width_Poly_Fit_Zone2_coeff2 = 1.780260502696301E-01;
742 long double Higgs_width_Poly_Fit_Zone2_coeff3 = -1.720311784419889E-03;
743 long double Higgs_width_Poly_Fit_Zone2_coeff4 = 1.038418605369741E-05;
744 long double Higgs_width_Poly_Fit_Zone2_coeff5 = -4.092496883922424E-08;
745 long double Higgs_width_Poly_Fit_Zone2_coeff6 = 1.067667966800388E-10;
746 long double Higgs_width_Poly_Fit_Zone2_coeff7 = -1.823343280081685E-13;
747 long double Higgs_width_Poly_Fit_Zone2_coeff8 = 1.955637395597351E-16;
748 long double Higgs_width_Poly_Fit_Zone2_coeff9 = -1.193287048560413E-19;
749 long double Higgs_width_Poly_Fit_Zone2_coeff10 = 3.156196649452213E-23;
750
751 long double Higgs_width_Poly_Fit_Zone3_coeff0 = -5.255605465437446E+02;
752 long double Higgs_width_Poly_Fit_Zone3_coeff1 = 1.036972988796150E+01;
753 long double Higgs_width_Poly_Fit_Zone3_coeff2 = -6.817022987365029E-02;
754 long double Higgs_width_Poly_Fit_Zone3_coeff3 = 1.493275723660056E-04;
755
756 long double m_m4l__2 = m_m4l * m_m4l;
757 long double m_m4l__3 = m_m4l__2 * m_m4l;
758
759 if( m_m4l < 156.5 ) return ( Higgs_width_Poly_Fit_Zone1_coeff0
760 + Higgs_width_Poly_Fit_Zone1_coeff1 * m_m4l
761 + Higgs_width_Poly_Fit_Zone1_coeff2 * m_m4l__2
762 + Higgs_width_Poly_Fit_Zone1_coeff3 * m_m4l__3
763 + Higgs_width_Poly_Fit_Zone1_coeff4 * m_m4l__2 * m_m4l__2
764 + Higgs_width_Poly_Fit_Zone1_coeff5 * m_m4l__2 * m_m4l__3
765 + Higgs_width_Poly_Fit_Zone1_coeff6 * m_m4l__2 * m_m4l__2 * m_m4l__2
766 + Higgs_width_Poly_Fit_Zone1_coeff7 * m_m4l__2 * m_m4l__2 * m_m4l__3
767 + Higgs_width_Poly_Fit_Zone1_coeff8 * m_m4l__2 * m_m4l__2 * m_m4l__2 * m_m4l__2
768 + Higgs_width_Poly_Fit_Zone1_coeff9 * m_m4l__2 * m_m4l__2 * m_m4l__2 * m_m4l__3 );
769 else if( m_m4l >= 156.5 && m_m4l <= 162 ) return ( Higgs_width_Poly_Fit_Zone3_coeff0
770 + Higgs_width_Poly_Fit_Zone3_coeff1 * m_m4l
771 + Higgs_width_Poly_Fit_Zone3_coeff2 * m_m4l__2
772 + Higgs_width_Poly_Fit_Zone3_coeff3 * m_m4l__3 );
773 else return ( Higgs_width_Poly_Fit_Zone2_coeff0
774 + Higgs_width_Poly_Fit_Zone2_coeff1 * m_m4l
775 + Higgs_width_Poly_Fit_Zone2_coeff2 * m_m4l__2
776 + Higgs_width_Poly_Fit_Zone2_coeff3 * m_m4l__3
777 + Higgs_width_Poly_Fit_Zone2_coeff4 * m_m4l__2 * m_m4l__2
778 + Higgs_width_Poly_Fit_Zone2_coeff5 * m_m4l__2 * m_m4l__3
779 + Higgs_width_Poly_Fit_Zone2_coeff6 * m_m4l__2 * m_m4l__2 * m_m4l__2
780 + Higgs_width_Poly_Fit_Zone2_coeff7 * m_m4l__2 * m_m4l__2 * m_m4l__3
781 + Higgs_width_Poly_Fit_Zone2_coeff8 * m_m4l__2 * m_m4l__2 * m_m4l__2 * m_m4l__2
782 + Higgs_width_Poly_Fit_Zone2_coeff9 * m_m4l__2 * m_m4l__2 * m_m4l__2 * m_m4l__3
783 + Higgs_width_Poly_Fit_Zone2_coeff10 * m_m4l__2 * m_m4l__2 * m_m4l__2 * m_m4l__2 * m_m4l__2 );
784 }
785
786 }; // class Parameters_heft
787
788 // calculate LO ME
789 class CPPProcess_P0_Sigma_heft_pp_H_ZZ_4l_heft_gg_epemmupmum {
790
791 public:
792
793 // Constructor.
794 CPPProcess_P0_Sigma_heft_pp_H_ZZ_4l_heft_gg_epemmupmum() :
795 pars(),
796 pout(4, vector<double>(4, 0.)) { }
797
798 // Destructor.
799 virtual ~CPPProcess_P0_Sigma_heft_pp_H_ZZ_4l_heft_gg_epemmupmum() {}
800
801 float Compute(const Quadruplet& quad) {
802 FourMomentum cms = quad.mom();
803 // use the cms mass as a value for MH
804 pars.set4lepMass(cms.mass());
805
806 const FourMomentum* fermionsMom[4] = {
807 &quad.Z1().first.momentum(),
808 &quad.Z1().second.momentum(),
809 &quad.Z2().first.momentum(),
810 &quad.Z2().second.momentum()
811 };
812
813 // boost to center-of-mass frame
814 LorentzTransform HRF_boost = LorentzTransform::mkFrameTransformFromBeta(cms.betaVec());
815 for(std::size_t i = 0; i < 4; ++i) {
816 FourMomentum tmpMom = HRF_boost.transform(*fermionsMom[i]);
817 pout[i][0] = tmpMom.E();
818 pout[i][1] = tmpMom.px();
819 pout[i][2] = tmpMom.py();
820 pout[i][3] = tmpMom.pz();
821 }
822
823 // Evaluate matrix element
824 return sigmaKin();
825 }
826
827 private:
828
829 // Calculate flavour-independent parts of cross section. Evaluate |M|^2, part independent of
830 // incoming flavour. Return matrix element
831
832 double sigmaKin() {
833 // Local variables and constants
834 static const int ncomb = 16;
835 // Helicities for the process
836 static const int helicities[ncomb][4] = {
837 {-1, -1, -1, -1}, {-1, -1, -1, 1}, {-1, -1, 1, -1}, {-1, -1, 1, 1},
838 {-1, 1, -1, -1}, {-1, 1, -1, 1}, {-1, 1, 1, -1}, {-1, 1, 1, 1},
839 { 1, -1, -1, -1}, { 1, -1, -1, 1}, { 1, -1, 1, -1}, { 1, -1, 1, 1},
840 { 1, 1, -1, -1}, { 1, 1, -1, 1}, { 1, 1, 1, -1}, { 1, 1, 1, 1}
841 };
842
843 // Reset the matrix elements
844 double matrix_element = 0.;
845
846 // Calculate the matrix element for all helicities
847 for(int ihel = 0; ihel < ncomb; ihel++) {
848 matrix_element += matrix_gg_h_h_zz_z_epem_z_mupmum(helicities[ihel]);
849 }
850
851 // Denominators: spins, colors and identical particles
852 return matrix_element /= 128.;
853 }
854
855 // Calculate wavefunctions and matrix elements for all subprocesses
856 double matrix_gg_h_h_zz_z_epem_z_mupmum(const int hel[]){
857 // Calculate all wavefunctions
858 ixx(pout[0], hel[0], w[2], true);
859 ixx(pout[1], hel[1], w[3], false);
860 FFV2_4_3(w[2], w[3], pars.GC_40, pars.GC_54,
861 pars.mdl_MZ, pars.mdl_WZ, w[4]);
862 ixx(pout[2], hel[2], w[5], true);
863 ixx(pout[3], hel[3], w[6], false);
864 FFV2_4_3(w[5], w[6], pars.GC_40, pars.GC_54,
865 pars.mdl_MZ, pars.mdl_WZ, w[7]);
866 VVS2_3(w[4], w[7], pars.GC_73, pars.mdl_MH, pars.mdl_WH, w[8]);
867
868 // Calculate all amplitudes
869 // Amplitude(s) for diagram number 0
870 std::complex<double> amp = VVS3_0(pars.mdl_MH, w[8], pars.GC_13);
871
872 // Calculate color flows
873 // Sum and square the color flows to get the matrix element
874 return real(8. * amp * conj(amp));
875 }
876
877 // wave function,
878 // ixx true takes anti-particle, false takes particle
879 void ixx(std::vector<double> p, int nhel, std::complex<double> fi[6], bool isixx) {
880 std::complex<double> chi[2];
881 double sqp0p3;
882 fi[0] = std::complex<double> (p[0], p[3]);
883 fi[1] = std::complex<double> (p[1], p[2]);
884 if (p[1] == 0.0 and p[2] == 0.0 and p[3] < 0.0) sqp0p3 = 0.0;
885 else sqp0p3 = pow(max(p[0] + p[3], 0.0), 0.5);
886 if (isixx) chi[0] = std::complex<double> (-sqp0p3, 0.0);
887 else chi[0] = std::complex<double> (sqp0p3, 0.0);
888 if (sqp0p3 == 0.0) chi[1] = std::complex<double> (-nhel * pow(2.0 * p[0], 0.5), 0.0);
889 else chi[1] = std::complex<double> (nhel * p[1], -p[2])/sqp0p3;
890 if (isixx) {
891 if (nhel == 1) {
892 fi[2] = chi[1];
893 fi[3] = chi[0];
894 fi[4] = std::complex<double> (0.0, 0.0);
895 fi[5] = std::complex<double> (0.0, 0.0);
896 } else {
897 fi[2] = std::complex<double> (0.0, 0.0);
898 fi[3] = std::complex<double> (0.0, 0.0);
899 fi[4] = chi[0];
900 fi[5] = chi[1];
901 }
902 } else {
903 if(nhel == 1){
904 fi[2] = chi[0];
905 fi[3] = chi[1];
906 fi[4] = std::complex<double> (0.00, 0.00);
907 fi[5] = std::complex<double> (0.00, 0.00);
908 } else {
909 fi[2] = std::complex<double> (0.00, 0.00);
910 fi[3] = std::complex<double> (0.00, 0.00);
911 fi[4] = chi[1];
912 fi[5] = chi[0];
913 }
914 }
915 return;
916 }
917
918 // vertices
919 void VVS2_3(std::complex<double> V1[], std::complex<double> V2[],
920 std::complex<double> COUP,
921 double M3, double W3, std::complex<double> S3[]) {
922 std::complex<double> cI = std::complex<double> (0., 1.);
923 std::complex<double> TMP1;
924 double P3[4];
925 std::complex<double> denom;
926 S3[0] = +V1[0] + V2[0];
927 S3[1] = +V1[1] + V2[1];
928 P3[0] = -S3[0].real();
929 P3[1] = -S3[1].real();
930 P3[2] = -S3[1].imag();
931 P3[3] = -S3[0].imag();
932 TMP1 = (V2[2] * V1[2] - V2[3] * V1[3] - V2[4] * V1[4] - V2[5] * V1[5]);
933 denom = COUP/(pow(P3[0], 2) - pow(P3[1], 2) - pow(P3[2], 2) - pow(P3[3], 2) -
934 M3 * (M3 - cI * W3));
935 S3[2] = denom * cI * TMP1;
936 }
937
938 void FFV2_4_3(std::complex<double> F1[], std::complex<double> F2[],
939 std::complex<double> COUP1, std::complex<double> COUP2,
940 double M3, double W3, std::complex<double> V3[]) {
941
942 std::complex<double> cI = std::complex<double> (0., 1.);
943 std::complex<double> denom;
944 std::complex<double> TMP11;
945 double P3[4];
946 std::complex<double> TMP14;
947 double OM3 = 1./pow(M3, 2);
948 V3[0] = +F1[0] + F2[0];
949 V3[1] = +F1[1] + F2[1];
950 P3[0] = -V3[0].real();
951 P3[1] = -V3[1].real();
952 P3[2] = -V3[1].imag();
953 P3[3] = -V3[0].imag();
954 TMP14 = (F1[4] * (F2[2] * (P3[0] - P3[3]) - F2[3] * (P3[1] + cI * (P3[2]))) +
955 F1[5] * (F2[2] * (+cI * (P3[2]) - P3[1]) + F2[3] * (P3[0] + P3[3])));
956 TMP11 = (F1[2] * (F2[4] * (P3[0] + P3[3]) + F2[5] * (P3[1] + cI * (P3[2]))) +
957 F1[3] * (F2[4] * (P3[1] - cI * (P3[2])) + F2[5] * (P3[0] - P3[3])));
958 denom = 1. / (pow(P3[0], 2) - pow(P3[1], 2) - pow(P3[2], 2) - pow(P3[3], 2) -
959 M3 * (M3 - cI * W3));
960 V3[2] = COUP2 * denom * - 2. * cI * (OM3 * - 1./2. * P3[0] * (TMP11 + 2. * (TMP14)) +
961 (+1./2. * (F2[4] * F1[2] + F2[5] * F1[3]) + F2[2] * F1[4] + F2[3] * F1[5]));
962 V3[3] = COUP2 * denom * - 2. * cI * (OM3 * - 1./2. * P3[1] * (TMP11 + 2. * (TMP14)) +
963 (-1./2. * (F2[5] * F1[2] + F2[4] * F1[3]) + F2[3] * F1[4] + F2[2] * F1[5]));
964 V3[4] = COUP2 * denom * 2. * cI * (OM3 * 1./2. * P3[2] * (TMP11 + 2. * (TMP14)) +
965 (+1./2. * cI * (F2[5] * F1[2]) - 1./2. * cI * (F2[4] * F1[3]) - cI *
966 (F2[3] * F1[4]) + cI * (F2[2] * F1[5])));
967 V3[5] = COUP2 * denom * 2. * cI * (OM3 * 1./2. * P3[3] * (TMP11 + 2. * (TMP14)) +
968 (+1./2. * (F2[4] * F1[2]) - 1./2. * (F2[5] * F1[3]) - F2[2] * F1[4] + F2[3] * F1[5]));
969 V3[2] += COUP1 * denom * - cI * (F2[4] * F1[2] + F2[5] * F1[3] - P3[0] * OM3 * TMP11);
970 V3[3] += COUP1 * denom * - cI * (-F2[5] * F1[2] - F2[4] * F1[3] - P3[1] * OM3 * TMP11);
971 V3[4] += COUP1 * denom * - cI * (-cI * (F2[5] * F1[2]) + cI * (F2[4] * F1[3]) - P3[2] * OM3 * TMP11);
972 V3[5] += COUP1 * denom * - cI * (F2[5] * F1[3] - F2[4] * F1[2] - P3[3] * OM3 * TMP11);
973 }
974
975 std::complex<double> VVS3_0(double mass,
976 std::complex<double> S3[],
977 std::complex<double> COUP) {
978 std::complex<double> TMP15;
979 TMP15 = std::complex<double> (0., pow(mass, 2) / 2.);
980 return COUP * S3[2] * (TMP15);
981 }
982
983 static const int nwavefuncs = 9;
984 std::complex<double> w[nwavefuncs][18];
985
986 // Pointer to the model parameters
987 Parameters_heft pars;
988
989 // vector with momenta (to be changed each event)
990 std::vector < std::vector<double> > pout;
991 };
992
993 CPPProcess_P0_Sigma_heft_pp_H_ZZ_4l_heft_gg_epemmupmum MGME;
994
995 };
996
997
998
999 RIVET_DECLARE_PLUGIN(ATLAS_2020_I1790439);
1000
1001}
|