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