Rivet analyses
title: ATLAS_2020_I1790439
Higgs to 4-lepton at 13 TeV
Experiment: ATLAS (ATLAS)
Inspire ID: 1790439
Status: VALIDATED
Authors: - Rongkun Wang - Neil Warrack - Roger Naranjo
References: - Expt page: ATLAS-HIGG-2018-29 - arXiv: 2004.03969 - submitted to EPJC
Beams: p+ p+
Beam energies: (6500.0, 6500.0)GeV
Run details: - routine expects Higgs to 4l processes including ggH, VBF, VH and ttH. The total production cross section including tau decays should be passed and the routine applies branching ratios.
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
```c++ // -- C++ --
include "Rivet/Analysis.hh"
include "Rivet/Projections/FastJets.hh"
include "Rivet/Projections/FinalState.hh"
include "Rivet/Projections/LeptonFinder.hh"
include "Rivet/Projections/PromptFinalState.hh"
include "Rivet/Projections/VetoedFinalState.hh"
namespace Rivet {
/// @brief H->ZZ->4l 13 TeV analysis class ATLAS_2020_I1790439 : public Analysis { public:
/// Default constructor
ATLAS_2020_I1790439()
: Analysis("ATLAS_2020_I1790439"), MGME() { }
void init() {
/// Dressed leptons
Cut cut_lep = (Cuts::abseta < 2.7) && (Cuts::pT > 5 * GeV);
PromptFinalState prompt_photons(Cuts::abspid == PID::PHOTON);
PromptFinalState prompt_leptons(Cuts::abspid == PID::MUON || Cuts::abspid == PID::ELECTRON);
LeptonFinder dLeptons(prompt_leptons, prompt_photons, 0.1, cut_lep);
declare(dLeptons, "AllLeptons");
/// Jet inputs
FinalState fs_jet(Cuts::abseta < 5.0);
VetoedFinalState jet_input(fs_jet);
// reject all leptons dressed with only prompt photons from jet input
FinalState leptons(Cuts::abspid == PID::ELECTRON || Cuts::abspid == PID::MUON);
LeptonFinder reject_leptons(leptons, prompt_photons, 0.1);
jet_input.addVetoOnThisFinalState(reject_leptons);
// reject prompt invisibles, including from tau decays
VetoedFinalState invis_fs_jet(fs_jet);
invis_fs_jet.addVetoOnThisFinalState(VisibleFinalState(fs_jet));
PromptFinalState invis_pfs_jet = PromptFinalState(invis_fs_jet, TauDecaysAs::PROMPT);
jet_input.addVetoOnThisFinalState(invis_pfs_jet);
// declare jets
FastJets jets(jet_input, JetAlg::ANTIKT, 0.4, JetMuons::NONE, JetInvisibles::DECAY);
declare(jets, "Jets");
// Book histograms with continuous (float) binning
book(_h["H4l_pt"], 5, 1, 1);
book(_h["Z1_m"], 7, 1, 1);
book(_h["Z2_m"], 9, 1, 1);
book(_h["abshiggs_y"], 11, 1, 1);
book(_h["abscthstr"], 13, 1, 1);
book(_h["cth1"], 15, 1, 1);
book(_h["cth2"], 17, 1, 1);
book(_h["phi"], 19, 1, 1);
book(_h["phi1"], 21, 1, 1);
book(_h["pt4lvy4l_0_0p5"], 51, 1, 1);
book(_h["pt4lvy4l_0p5_1"], 51, 1, 3);
book(_h["pt4lvy4l_1_1p5"], 51, 1, 5);
book(_h["pt4lvy4l_1p5_2p5"], 51, 1, 7);
book(_h["pt4lvnjet_0"], 53, 1, 1);
book(_h["pt4lvnjet_1"], 53, 1, 3);
book(_h["pt4lvnjet_2"], 53, 1, 5);
book(_h["pt4lvnjet_3"], 53, 1, 7);
book(_h["Z1_m_4l"], 65, 1, 1);
book(_h["Z1_m_2l2l"], 66, 1, 1);
book(_h["Z2_m_4l"], 68, 1, 1);
book(_h["Z2_m_2l2l"], 69, 1, 1);
book(_h["phi_4l"], 71, 1, 1);
book(_h["phi_2l2l"], 72, 1, 1);
// Book histograms with discrete (string) binning
book(_s["xs_flavour"], 3, 1, 1);
book(_s["n_jets"], 23, 1, 1);
book(_s["n_jets_incl"], 25, 1, 1);
book(_s["n_bjets"], 26, 1, 1);
book(_s["jet_pt_leading"], 28, 1, 1);
book(_s["jet_pt_subleading"], 30, 1, 1);
book(_s["dijet_m"], 32, 1, 1);
book(_s["dijet_deltaeta"], 34, 1, 1);
book(_s["dijet_deltaphi"], 36, 1, 1);
book(_s["pt4lj"], 38, 1, 1);
book(_s["pt4ljj"], 40, 1, 1);
book(_s["m4lj"], 42, 1, 1);
book(_s["m4ljj"], 44, 1, 1);
book(_s["m12vsm34"], 46, 1, 1);
book(_s["m12vsm34_2l2m"], 48, 1, 1);
book(_s["m12vsm34_2l2e"], 49, 1, 1);
book(_s["pt4lvpt4lj"], 55, 1, 1);
book(_s["pt4ljvm4lj"], 57, 1, 1);
book(_s["pt4lvptj0"], 59, 1, 1);
book(_s["ptj0vyj0"], 61, 1, 1);
book(_s["ptj0vptj1"], 63, 1, 1);
book(_s["m12vsm34_4l"], 74, 1, 1);
book(_s["m12vsm34_2l2l"], 75, 1, 1);
// helper axes to map onto discrete binning labels
_axisMap["jet_pt_leading"] = YODA::Axis<double>({30., 60., 120., 350.});
_axisMap["pt4lj"] = YODA::Axis<double>({0., 60., 120., 350.});
_axisMap["m4lj"] = YODA::Axis<double>({120., 180., 220., 300., 400., 600., 2000.});
_axisMap["jet_pt_subleading"] = YODA::Axis<double>({30., 60., 120., 120., 350.});
_axisMap["dijet_m"] = YODA::Axis<double>({0., 120., 450., 3000.});
_axisMap["dijet_deltaeta"] = YODA::Axis<double>({0., 1., 2.5, 9.});
_axisMap["dijet_deltaphi"] = YODA::Axis<double>({0.0, 1.571, 3.142, 4.712, 6.283});
_axisMap["pt4ljj"] = YODA::Axis<double>({0., 60., 120.});
_axisMap["m4ljj"] = YODA::Axis<double>({180., 320., 450., 600., 1000., 2500.});
}
/// Do the per-event analysis
void analyze(const Event& e) {
if (edges.empty()) {
for (const string& label : vector<string>{"xs_flavour", "n_jets", "n_jets_incl",
"n_bjets", "jet_pt_leading", "jet_pt_subleading",
"dijet_m", "dijet_deltaeta", "dijet_deltaphi",
"pt4lj", "pt4ljj", "m4lj",
"m4ljj", "m12vsm34", "m12vsm34_2l2m",
"m12vsm34_2l2e", "pt4lvpt4lj", "pt4ljvm4lj",
"pt4lvptj0", "ptj0vyj0", "ptj0vptj1",
"m12vsm34_4l", "m12vsm34_2l2l"}) {
edges[label] = _s[label]->xEdges();
}
}
_s["xs_flavour"]->fill(edges["xs_flavour"][8]);
const DressedLeptons& all_leps = apply<LeptonFinder>(e, "AllLeptons").dressedLeptons();
unsigned int n_parts = all_leps.size();
unsigned int n_OSSF_pairs = 0;
std::vector<Zstate> dileptons;
Particles leptons;
// Form Z candidate (opposite-sign same-flavour) lepton pairs
for (unsigned int i = 0; i < n_parts; i++) {
for (unsigned int j = i + 1; j < n_parts; j++) {
if (isOSSF(all_leps[i], all_leps[j])) {
n_OSSF_pairs += 1;
// Set positive charge first for later ME calculation
if (all_leps[i].charge() > 0) {
dileptons.push_back(Zstate(ParticlePair(all_leps[i], all_leps[j])));
}
else {
dileptons.push_back(Zstate(ParticlePair(all_leps[j], all_leps[i])));
}
}
}
}
// At least two pairs required to select ZZ->llll final state
if (n_OSSF_pairs < 2) vetoEvent;
// Form the quadruplet of two lepon pairs passing kinematics cuts
std::vector<Quadruplet> quadruplets;
for (unsigned int i = 0; i < dileptons.size(); i++) {
for (unsigned int j = i + 1; j < dileptons.size(); j++) {
// Only use unique leptons
if (isSame(dileptons[i].first, dileptons[j].first)) continue;
if (isSame(dileptons[i].first, dileptons[j].second)) continue;
if (isSame(dileptons[i].second, dileptons[j].first)) continue;
if (isSame(dileptons[i].second, dileptons[j].second)) continue;
leptons.clear();
leptons.push_back(dileptons[i].first);
leptons.push_back(dileptons[i].second);
leptons.push_back(dileptons[j].first);
leptons.push_back(dileptons[j].second);
leptons = sortByPt(leptons);
// Apply kinematic cuts
if (leptons[0].pt() < 20 * GeV) continue;
if (leptons[1].pt() < 15 * GeV) continue;
if (leptons[2].pt() < 10 * GeV) continue;
// Form the quad with pair closest to Z pole first
if (dileptons[i].Zdist() < dileptons[j].Zdist()) {
quadruplets.push_back(Quadruplet(dileptons[i], dileptons[j]));
}
else {
quadruplets.push_back(Quadruplet(dileptons[j], dileptons[i]));
}
}
}
// Veto if no quad passes kinematic selection
unsigned int n_quads = quadruplets.size();
if (n_quads == 0) vetoEvent;
// To resolve ambiguities in lepton pairing order quads by channel priority first, then m12 - mz and m34 - mz
// The first in every channel is considered nominal
std::sort(quadruplets.begin(), quadruplets.end(), [](const Quadruplet& q1, const Quadruplet& q2) {
if (q1.ch_priority == q2.ch_priority) {
// if rarely, Z1 the same distance from the Z pole, compare Z2
if (fabs(q1.Z1().Zdist() - q2.Z1().Zdist()) < 1.e-5)
return q1.Z2().Zdist() < q2.Z2().Zdist();
else
return q1.Z1().Zdist() < q2.Z1().Zdist();
}
else
return q1.ch_priority < q2.ch_priority;
});
// Select the best quad
Particles leptons_sel4l;
Quadruplet quadSel;
float MEHZZ_max = -999.;
int prevQuadType = -999.;
bool atleastonequadpassed = false;
bool isNominalQuad = false;
bool extraLep = false;
for (unsigned int iquad = 0; iquad < n_quads; iquad++) {
// Veto event if nominal quad was not selected in 4 lepton case
if (n_parts == 4 && iquad > 0) vetoEvent;
int quadType = (int)quadruplets[iquad].type();
if (quadType != prevQuadType) {
isNominalQuad = true;
}
else {
isNominalQuad = false;
}
prevQuadType = quadType;
Quadruplet& quad = quadruplets[iquad];
// Z invariant mass requirements
if (!(inRange(quad.Z1().mom().mass(), 50 * GeV, 106 * GeV))) continue;
if (!(inRange(quad.Z2().mom().mass(), 12 * GeV, 115 * GeV))) continue;
// Lepton separation and J/Psi veto
bool b_pass_leptonseparation = true;
bool b_pass_jpsi = true;
leptons_sel4l.clear();
leptons_sel4l.push_back(quad.Z1().first);
leptons_sel4l.push_back(quad.Z1().second);
leptons_sel4l.push_back(quad.Z2().first);
leptons_sel4l.push_back(quad.Z2().second);
for (unsigned int i = 0; i < 4; i++) {
for (unsigned int j = i + 1; j < 4; j++) {
if (deltaR(leptons_sel4l[i], leptons_sel4l[j]) < 0.1) b_pass_leptonseparation = false;
if (isOSSF(leptons_sel4l[i], leptons_sel4l[j])
&& (leptons_sel4l[i].mom() + leptons_sel4l[j].mom()).mass() <= 5. * GeV)
b_pass_jpsi = false;
}
}
if (b_pass_leptonseparation == false || b_pass_jpsi == false) continue;
// Only consider the event if at least one nominal quadruplet passes cuts
if (isNominalQuad) atleastonequadpassed = true;
// Direct selection for case with only 4 prompt leptons
if (n_parts == 4) {
quadSel = quad;
break;
}
// In cases with extra leptons meeting further requirements use max ME to select quad
float MEHZZ;
if (quadType == 0 || quadType == 1)
MEHZZ = MGME.Compute(quadruplets[iquad]) / 2;
else
MEHZZ = MGME.Compute(quadruplets[iquad]);
// Check leptons other than the ones from the channel's nominal quadruplet
// and don't need to recheck extraLep for a second channel
if (isNominalQuad && !extraLep) {
for (const Particle& lep : all_leps) {
bool lep_in_quad = false;
bool lep_too_close = false;
// pT requirement
if (lep.pt() < 12 * GeV) continue;
// skip if lepton included in quad or not isolated from quad leptons
for (unsigned int i = 0; i < 4; i++) {
if (isSame(lep, leptons_sel4l[i])) lep_in_quad = true;
if (deltaR(lep, leptons_sel4l[i]) < 0.1) lep_too_close = true;
}
if (lep_in_quad || lep_too_close) continue;
extraLep = true;
break;
}
if (!extraLep) {
// In case of no suitable extra leptons select the quad directly
quadSel = quad;
break;
}
}
// Use ME to select the quad when there are extra leptons
if (MEHZZ > MEHZZ_max) {
quadSel = quad;
MEHZZ_max = MEHZZ;
}
}
if (!atleastonequadpassed) vetoEvent;
// Veto if quad not in Higgs mass window
FourMomentum Higgs = quadSel.mom();
double H4l_mass = Higgs.mass();
if (!(inRange(H4l_mass, 105. * GeV, 160. * GeV))) vetoEvent;
// Higgs observables
double H4l_pt = Higgs.pt();
double H4l_rapidity = Higgs.absrapidity();
LorentzTransform HRF_boost = LorentzTransform::mkFrameTransformFromBeta(Higgs.betaVec());
FourMomentum Z1_in_HRF = HRF_boost.transform(quadSel.Z1().mom());
FourMomentum Z2_in_HRF = HRF_boost.transform(quadSel.Z2().mom());
double H4l_costheta = fabs(cos(Z1_in_HRF.theta()));
double H4l_m12 = quadSel.Z1().mom().mass();
double H4l_m34 = quadSel.Z2().mom().mass();
FourMomentum v11_HRF = HRF_boost.transform(quadSel.Z1().second.mom());
FourMomentum v12_HRF = HRF_boost.transform(quadSel.Z1().first.mom());
FourMomentum v21_HRF = HRF_boost.transform(quadSel.Z2().second.mom());
FourMomentum v22_HRF = HRF_boost.transform(quadSel.Z2().first.mom());
Vector3 v11 = v11_HRF.p3();
Vector3 v12 = v12_HRF.p3();
Vector3 v21 = v21_HRF.p3();
Vector3 v22 = v22_HRF.p3();
Vector3 nz(0, 0, 1);
Vector3 qz1 = Z1_in_HRF.p3();
Vector3 n1p = v11.cross(v12).unit();
Vector3 n2p = v21.cross(v22).unit();
Vector3 nscp = nz.cross(qz1).unit();
double H4l_Phi = qz1.dot(n1p.cross(n2p)) / fabs(qz1.dot(n1p.cross(n2p))) * acos(-n1p.dot(n2p));
double H4l_Phi1 = qz1.dot(n1p.cross(nscp)) / fabs(qz1.dot(n1p.cross(nscp))) * acos(n1p.dot(nscp));
LorentzTransform Z1RF_boost = LorentzTransform::mkFrameTransformFromBeta(Z1_in_HRF.betaVec());
LorentzTransform Z2RF_boost = LorentzTransform::mkFrameTransformFromBeta(Z2_in_HRF.betaVec());
FourMomentum Z1_in_Z2RF = Z2RF_boost.transform(Z1_in_HRF);
FourMomentum Z2_in_Z1RF = Z1RF_boost.transform(Z2_in_HRF);
Vector3 Z1_p3 = Z1_in_Z2RF.p3();
Vector3 Z2_p3 = Z2_in_Z1RF.p3();
FourMomentum n_Z1 = Z1RF_boost.transform(v11_HRF);
FourMomentum n_Z2 = Z2RF_boost.transform(v21_HRF);
// angle is negative with its Z
double H4l_cth1 = -cos(n_Z1.p3().angle(Z2_p3));
double H4l_cth2 = -cos(n_Z2.p3().angle(Z1_p3));
// Jet observables
Jets jets = apply<FastJets>(e, "Jets").jetsByPt(Cuts::pT > 30 * GeV && Cuts::absrap < 4.4);
// discard jets which overlap leptons
for (const Particle& l : all_leps) idiscard(jets, deltaRLess(l, 0.1));
// collect b-tagged jets
const Jets b_jets_btagged = select(jets, hasBTag(Cuts::pT > 5 * GeV));
unsigned int n_bjets = b_jets_btagged.size();
unsigned int n_jets = jets.size();
double leading_jet_pt = (n_jets > 0 ? jets[0].pT() : 0);
double leading_jet_y = (n_jets > 0 ? jets[0].absrapidity() : 0);
double subleading_jet_pt = (n_jets > 1 ? jets[1].pT() : 0);
FourMomentum Dijets = {0, 0, 0, 0};
if (n_jets > 1) Dijets = jets[0].mom() + jets[1].mom();
double mjj = (n_jets > 1 ? Dijets.mass() : -999);
double dphijj = -1;
if (n_jets > 1) {
if (jets[0].eta() > jets[1].eta())
dphijj = jets[0].phi() - jets[1].phi();
else
dphijj = jets[1].phi() - jets[0].phi();
if (dphijj < 0) dphijj = dphijj + TWOPI;
}
double detajj = (n_jets > 1 ? fabs(jets[0].eta() - jets[1].eta()) : -1);
FourMomentum m4lj, m4ljj;
if (n_jets > 0) m4lj = Higgs + jets[0];
if (n_jets > 1) m4ljj = Higgs + Dijets;
double H4l_m4lj = m4lj.mass();
double H4l_m4ljj = m4ljj.mass();
double H4l_pt4lj = m4lj.pt();
double H4l_pt4ljj = m4ljj.pt();
// Fill histograms
// Branching ratios for 4l and 2l2l channels
double BR_SF = 0.00013;
double BR_OF = 0.000118;
if (inRange(H4l_mass, 115. * GeV, 130. * GeV)) {
if (quadSel.type() == Quadruplet::FlavCombi::mm || quadSel.type() == Quadruplet::FlavCombi::ee) {
_s["xs_flavour"]->fill(edges["xs_flavour"][(int)quadSel.type()], BR_SF);
_s["xs_flavour"]->fill(edges["xs_flavour"][4], BR_SF);
}
else if (quadSel.type() == Quadruplet::FlavCombi::em || quadSel.type() == Quadruplet::FlavCombi::me) {
_s["xs_flavour"]->fill(edges["xs_flavour"][(int)quadSel.type()], BR_OF);
_s["xs_flavour"]->fill(edges["xs_flavour"][5], BR_OF);
}
_s["xs_flavour"]->fill(edges["xs_flavour"][6], Br);
_s["xs_flavour"]->fill(edges["xs_flavour"][7], Br);
}
// Higgs variables
for (const auto& p : std::map<std::string, double>{{"H4l_pt", H4l_pt},
{"Z1_m", H4l_m12},
{"Z2_m", H4l_m34},
{"abshiggs_y", H4l_rapidity},
{"abscthstr", H4l_costheta},
{"cth1", H4l_cth1},
{"cth2", H4l_cth2},
{"phi", H4l_Phi},
{"phi1", H4l_Phi1}}) {
_h[p.first]->fill(p.second);
}
// Jet variables
discreteFill("n_jets", min((int)n_jets, 3));
discreteFill("n_jets_incl", 0);
if (n_jets >= 1) discreteFill("n_jets_incl", 1);
if (n_jets >= 2) discreteFill("n_jets_incl", 2);
if (n_jets >= 3) discreteFill("n_jets_incl", 3);
if (n_jets == 0)
discreteFill("n_bjets", 0);
else if (n_bjets == 0)
discreteFill("n_bjets", 1);
else if (n_bjets >= 1)
discreteFill("n_bjets", 2);
if (n_jets == 0) {
discreteFill("jet_pt_leading", 0);
discreteFill("pt4lj", 0);
discreteFill("m4lj", 0);
}
else if (n_jets >= 1) {
discreteFill("jet_pt_leading", leading_jet_pt);
discreteFill("pt4lj", H4l_pt4lj);
discreteFill("m4lj", H4l_m4lj);
}
if (n_jets < 2) {
discreteFill("jet_pt_subleading", 0);
discreteFill("dijet_m", 0);
discreteFill("dijet_deltaeta", 0);
discreteFill("dijet_deltaphi", 0);
discreteFill("pt4ljj", 0);
discreteFill("m4ljj", 0);
}
else if (n_jets >= 2) {
discreteFill("jet_pt_subleading", subleading_jet_pt);
discreteFill("dijet_m", mjj);
discreteFill("dijet_deltaeta", detajj);
discreteFill("dijet_deltaphi", dphijj);
discreteFill("pt4ljj", H4l_pt4ljj);
discreteFill("m4ljj", H4l_m4ljj);
}
// m12 vs m34 (all channels)
if (H4l_m12 < 82 && H4l_m34 < 32)
discreteFill("m12vsm34", 0);
else if (H4l_m12 < 74 && H4l_m34 > 32)
discreteFill("m12vsm34", 1);
else if (H4l_m12 > 74 && H4l_m34 > 32)
discreteFill("m12vsm34", 2);
else if (H4l_m12 > 82 && H4l_m34 < 32 && H4l_m34 > 24)
discreteFill("m12vsm34", 3);
else if (H4l_m12 > 82 && H4l_m34 < 24)
discreteFill("m12vsm34", 4);
if (quadSel.type() == Quadruplet::FlavCombi::em || quadSel.type() == Quadruplet::FlavCombi::mm) {
if (H4l_m12 < 82 && H4l_m34 < 32)
discreteFill("m12vsm34_2l2m", 0);
else if (H4l_m12 < 74 && H4l_m34 > 32)
discreteFill("m12vsm34_2l2m", 1);
else if (H4l_m12 > 74 && H4l_m34 > 32)
discreteFill("m12vsm34_2l2m", 2);
else if (H4l_m12 > 82 && H4l_m34 < 32 && H4l_m34 > 24)
discreteFill("m12vsm34_2l2m", 3);
else if (H4l_m12 > 82 && H4l_m34 < 24)
discreteFill("m12vsm34_2l2m", 4);
}
else if (quadSel.type() == Quadruplet::FlavCombi::me || quadSel.type() == Quadruplet::FlavCombi::ee) {
if (H4l_m12 < 82 && H4l_m34 < 32)
discreteFill("m12vsm34_2l2e", 0);
else if (H4l_m12 < 74 && H4l_m34 > 32)
discreteFill("m12vsm34_2l2e", 1);
else if (H4l_m12 > 74 && H4l_m34 > 32)
discreteFill("m12vsm34_2l2e", 2);
else if (H4l_m12 > 82 && H4l_m34 < 32 && H4l_m34 > 24)
discreteFill("m12vsm34_2l2e", 3);
else if (H4l_m12 > 82 && H4l_m34 < 24)
discreteFill("m12vsm34_2l2e", 4);
}
// m12 vs m34 (4l channels only)
if (quadSel.type() == Quadruplet::FlavCombi::ee || quadSel.type() == Quadruplet::FlavCombi::mm) {
if (H4l_m12 < 82 && H4l_m34 < 32)
discreteFill("m12vsm34_4l", 0);
else if (H4l_m12 < 74 && H4l_m34 > 32)
discreteFill("m12vsm34_4l", 1);
else if (H4l_m12 > 74 && H4l_m34 > 32)
discreteFill("m12vsm34_4l", 2);
else if (H4l_m12 > 82 && H4l_m34 < 32 && H4l_m34 > 24)
discreteFill("m12vsm34_4l", 3);
else if (H4l_m12 > 82 && H4l_m34 < 24)
discreteFill("m12vsm34_4l", 4);
_h["Z1_m_4l"]->fill(H4l_m12);
_h["Z2_m_4l"]->fill(H4l_m34);
_h["phi_4l"]->fill(H4l_Phi);
}
// m12 vs m34 (2l2l channels only)
if (quadSel.type() == Quadruplet::FlavCombi::me || quadSel.type() == Quadruplet::FlavCombi::em) {
if (H4l_m12 < 82 && H4l_m34 < 32)
discreteFill("m12vsm34_2l2l", 0);
else if (H4l_m12 < 74 && H4l_m34 > 32)
discreteFill("m12vsm34_2l2l", 1);
else if (H4l_m12 > 74 && H4l_m34 > 32)
discreteFill("m12vsm34_2l2l", 2);
else if (H4l_m12 > 82 && H4l_m34 < 32 && H4l_m34 > 24)
discreteFill("m12vsm34_2l2l", 3);
else if (H4l_m12 > 82 && H4l_m34 < 24)
discreteFill("m12vsm34_2l2l", 4);
_h["Z1_m_2l2l"]->fill(H4l_m12);
_h["Z2_m_2l2l"]->fill(H4l_m34);
_h["phi_2l2l"]->fill(H4l_Phi);
}
// 2d differential variables
if (0 < H4l_rapidity && H4l_rapidity < 0.5)
_h["pt4lvy4l_0_0p5"]->fill(H4l_pt);
else if (0.5 < H4l_rapidity && H4l_rapidity < 1)
_h["pt4lvy4l_0p5_1"]->fill(H4l_pt);
else if (1 < H4l_rapidity && H4l_rapidity < 1.5)
_h["pt4lvy4l_1_1p5"]->fill(H4l_pt);
else if (1.5 < H4l_rapidity && H4l_rapidity < 2.5)
_h["pt4lvy4l_1p5_2p5"]->fill(H4l_pt);
if (n_jets == 0)
_h["pt4lvnjet_0"]->fill(H4l_pt);
else if (n_jets == 1)
_h["pt4lvnjet_1"]->fill(H4l_pt);
else if (n_jets == 2)
_h["pt4lvnjet_2"]->fill(H4l_pt);
else if (n_jets > 2)
_h["pt4lvnjet_3"]->fill(H4l_pt);
if (n_jets == 0) {
discreteFill("pt4lvptj0", 0);
discreteFill("pt4lvpt4lj", 0);
discreteFill("pt4ljvm4lj", 0);
discreteFill("ptj0vptj1", 0);
discreteFill("ptj0vyj0", 0);
}
else {
if (0 < H4l_pt4lj && H4l_pt4lj < 60 && 0 < H4l_pt && H4l_pt < 120)
discreteFill("pt4lvpt4lj", 1);
else if (0 < H4l_pt4lj && H4l_pt4lj < 60 && 120 < H4l_pt && H4l_pt < 350)
discreteFill("pt4lvpt4lj", 2);
else if (60 < H4l_pt4lj && H4l_pt4lj < 350 && 0 < H4l_pt && H4l_pt < 120)
discreteFill("pt4lvpt4lj", 3);
else if (60 < H4l_pt4lj && H4l_pt4lj < 350 && 120 < H4l_pt && H4l_pt < 350)
discreteFill("pt4lvpt4lj", 4);
if (120 < H4l_m4lj && H4l_m4lj < 220 && 0 < H4l_pt4lj && H4l_pt4lj < 350)
discreteFill("pt4ljvm4lj", 1);
else if (220 < H4l_m4lj && H4l_m4lj < 350 && 0 < H4l_pt4lj && H4l_pt4lj < 60)
discreteFill("pt4ljvm4lj", 2);
else if (220 < H4l_m4lj && H4l_m4lj < 350 && 60 < H4l_pt4lj && H4l_pt4lj < 350)
discreteFill("pt4ljvm4lj", 3);
else if (350 < H4l_m4lj && H4l_m4lj < 2000 && 0 < H4l_pt4lj && H4l_pt4lj < 350)
discreteFill("pt4ljvm4lj", 4);
if (30 < leading_jet_pt && leading_jet_pt < 60 && 0 < H4l_pt && H4l_pt < 80)
discreteFill("pt4lvptj0", 1);
else if (30 < leading_jet_pt && leading_jet_pt < 60 && 80 < H4l_pt && H4l_pt < 350)
discreteFill("pt4lvptj0", 2);
else if (60 < leading_jet_pt && leading_jet_pt < 120 && 0 < H4l_pt && H4l_pt < 120)
discreteFill("pt4lvptj0", 3);
else if (60 < leading_jet_pt && leading_jet_pt < 120 && 120 < H4l_pt && H4l_pt < 350)
discreteFill("pt4lvptj0", 4);
else if (120 < leading_jet_pt && leading_jet_pt < 350 && 0 < H4l_pt && H4l_pt < 120)
discreteFill("pt4lvptj0", 5);
else if (120 < leading_jet_pt && leading_jet_pt < 350 && 120 < H4l_pt && H4l_pt < 350)
discreteFill("pt4lvptj0", 6);
if (30 < leading_jet_pt && leading_jet_pt < 120 && 0 < leading_jet_y && leading_jet_y < 0.8) {
discreteFill("ptj0vyj0", 1);
}
else if (30 < leading_jet_pt && leading_jet_pt < 120 && 0.8 < leading_jet_y && leading_jet_y < 1.7) {
discreteFill("ptj0vyj0", 2);
}
else if (30 < leading_jet_pt && leading_jet_pt < 120 && 1.7 < leading_jet_y) {
discreteFill("ptj0vyj0", 3);
}
else if (120 < leading_jet_pt && leading_jet_pt < 350 && 0 < leading_jet_y && leading_jet_y < 1.7) {
discreteFill("ptj0vyj0", 4);
}
else if (120 < leading_jet_pt && leading_jet_pt < 350 && 1.7 < leading_jet_y) {
discreteFill("ptj0vyj0", 5);
}
if (n_jets == 1 && 30 < leading_jet_pt && leading_jet_pt < 60)
discreteFill("ptj0vptj1", 1);
else if (n_jets == 1 && 60 < leading_jet_pt && leading_jet_pt < 350)
discreteFill("ptj0vptj1", 2);
else if (30 < leading_jet_pt && leading_jet_pt < 60 && 30 < subleading_jet_pt
&& subleading_jet_pt < 60) {
discreteFill("ptj0vptj1", 3);
}
else if (60 < leading_jet_pt && leading_jet_pt < 350 && 30 < subleading_jet_pt
&& subleading_jet_pt < 60) {
discreteFill("ptj0vptj1", 4);
}
else if (60 < leading_jet_pt && leading_jet_pt < 350 && 60 < subleading_jet_pt
&& subleading_jet_pt < 350) {
discreteFill("ptj0vptj1", 5);
}
}
}
template <typename T>
void discreteFill(const string& label, const T coord) {
if constexpr (std::is_floating_point<T>::value) {
discreteFill(label, _axisMap[label].index(coord));
}
else {
_s[label]->fill(edges[label][coord]);
}
}
void finalize() {
const double sf = crossSection() / femtobarn / sumOfWeights();
scale(_h, sf * Br);
for (auto& hist : _s) {
if (hist.first == "xs_flavour") {
scale(hist.second, sf);
}
else {
scale(hist.second, sf * Br);
}
}
}
private:
// Br(H-->ZZ) * BR(ZZ-->4l)
const double Br = 0.02641 * 0.004736842;
map<string, Histo1DPtr> _h;
map<string, BinnedHistoPtr<string>> _s;
map<string, vector<string>> edges;
map<string, YODA::Axis<double>> _axisMap;
/// Generic Z candidate
struct Zstate : public ParticlePair {
Zstate() { }
Zstate(ParticlePair _particlepair)
: ParticlePair(_particlepair) { }
FourMomentum mom() const {
return first.momentum() + second.momentum();
}
double Zdist() const {
return fabs(mom().mass() - 91.1876 * GeV);
}
int flavour() const {
return first.abspid();
}
};
/// Generic quadruplet
struct Quadruplet {
// find out which type it is: 4mu = 0, 4e = 1, 2mu2e = 2, 2e2mu = 3 (mm, ee, me, em)
// channel priority is 4m, 2e2m, 2m2e, 4e
enum class FlavCombi { mm = 0, ee, me, em, undefined };
Quadruplet() { }
Quadruplet(Zstate z1, Zstate z2)
: _z1(z1), _z2(z2) {
if (_z1.flavour() == 13 && _z2.flavour() == 13) {
_type = FlavCombi::mm;
ch_priority = 0;
}
else if (_z1.flavour() == 11 && _z2.flavour() == 11) {
_type = FlavCombi::ee;
ch_priority = 3;
}
else if (_z1.flavour() == 13 && _z2.flavour() == 11) {
_type = FlavCombi::me;
ch_priority = 2;
}
else if (_z1.flavour() == 11 && _z2.flavour() == 13) {
_type = FlavCombi::em;
ch_priority = 1;
}
else {
_type = FlavCombi::undefined;
}
}
Quadruplet(const Quadruplet& quad)
: _z1(quad._z1), _z2(quad._z2), _type(quad._type), ch_priority(quad.ch_priority) { }
Zstate _z1, _z2;
FlavCombi _type;
int ch_priority;
const Zstate& Z1() const {
return _z1;
}
const Zstate& Z2() const {
return _z2;
}
FourMomentum mom() const {
return _z1.mom() + _z2.mom();
}
FlavCombi type() const {
return _type;
}
};
// save and calculate parameters
class Parameters_heft {
public:
// Model parameters independent of aS
double mdl_WH, mdl_WZ, aS, mdl_Gf, aEWM1, mdl_MH, mdl_MZ, mdl_MTA, mdl_MT, mdl_MB, mdl_MP,
mdl_conjg__CKM3x3, mdl_CKM3x3, mdl_MZ__exp__2, mdl_MZ__exp__4, mdl_MH__exp__4, mdl_MT__exp__4,
mdl_MH__exp__2, mdl_MT__exp__2, mdl_MH__exp__6, mdl_MT__exp__6, mdl_aEW, mdl_MW, mdl_ee,
mdl_MW__exp__2, mdl_sw2, mdl_cw, mdl_sw, mdl_v, mdl_ee__exp__2;
std::complex<double> mdl_complexi;
// Model parameters dependent on aS
double mdl_GH;
// Model couplings independent of aS
std::complex<double> GC_40, GC_54, GC_73;
// Model couplings dependent on aS
std::complex<double> GC_13;
// Set parameters and couplings that are unchanged during the run
Parameters_heft() {
mdl_WH = 6.382339e-03;
mdl_WZ = 2.441404e+00;
aS = 1.180000e-01;
mdl_Gf = 1.166390e-05;
aEWM1 = 1.325070e+02;
mdl_MH = 1.250000e+02;
mdl_MZ = 9.118800e+01;
mdl_MT = 1.730000e+02;
mdl_complexi = std::complex<double>(0., 1.);
mdl_MZ__exp__2 = pow(mdl_MZ, 2.);
mdl_MZ__exp__4 = pow(mdl_MZ, 4.);
mdl_MH__exp__2 = pow(mdl_MH, 2.);
mdl_MT__exp__4 = pow(mdl_MT, 4.);
mdl_MH__exp__4 = pow(mdl_MH, 4.);
mdl_MT__exp__2 = pow(mdl_MT, 2.);
mdl_MH__exp__6 = pow(mdl_MH, 6.);
mdl_MT__exp__6 = pow(mdl_MT, 6.);
mdl_aEW = 1. / aEWM1;
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);
mdl_MW__exp__2 = pow(mdl_MW, 2.);
mdl_sw2 = 1. - mdl_MW__exp__2 / mdl_MZ__exp__2;
mdl_cw = sqrt(1. - mdl_sw2);
mdl_sw = sqrt(mdl_sw2);
mdl_v = (2. * mdl_MW * mdl_sw) / mdl_ee;
mdl_ee__exp__2 = pow(mdl_ee, 2.);
GC_40 = -(mdl_ee * mdl_complexi * mdl_cw) / (2. * mdl_sw);
GC_54 = (mdl_ee * mdl_complexi * mdl_sw) / (2. * mdl_cw);
GC_73 = mdl_ee__exp__2 * mdl_complexi * mdl_v
+ ((1. - mdl_sw2) * mdl_ee__exp__2 * mdl_complexi * mdl_v) / (2. * mdl_sw2)
+ (mdl_ee__exp__2 * mdl_complexi * mdl_sw2 * mdl_v) / (2. * (1. - mdl_sw2));
}
// Set Mass
void set4lepMass(double m_m4l) {
mdl_MH = m_m4l;
mdl_WH = setWidth(m_m4l);
mdl_MH__exp__2 = pow(mdl_MH, 2.);
mdl_MH__exp__4 = pow(mdl_MH, 4.);
mdl_MH__exp__6 = pow(mdl_MH, 6.);
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)
+ (7. * mdl_MH__exp__2) / (120. * mdl_MT__exp__2)))
/ (12. * pow(M_PI, 2.) * mdl_v);
GC_13 = -(mdl_complexi * mdl_GH);
}
private:
// Set Width
long double setWidth(double m_m4l) {
long double Higgs_width_Poly_Fit_Zone1_coeff0 = -1.450308902710193E+03;
long double Higgs_width_Poly_Fit_Zone1_coeff1 = 1.129291251156317E+02;
long double Higgs_width_Poly_Fit_Zone1_coeff2 = -3.893063071316150E+00;
long double Higgs_width_Poly_Fit_Zone1_coeff3 = 7.798666884832531E-02;
long double Higgs_width_Poly_Fit_Zone1_coeff4 = -1.000455877406390E-03;
long double Higgs_width_Poly_Fit_Zone1_coeff5 = 8.523735379647125E-06;
long double Higgs_width_Poly_Fit_Zone1_coeff6 = -4.823164754652171E-08;
long double Higgs_width_Poly_Fit_Zone1_coeff7 = 1.747954506786346E-10;
long double Higgs_width_Poly_Fit_Zone1_coeff8 = -3.681723572169337E-13;
long double Higgs_width_Poly_Fit_Zone1_coeff9 = 3.434207075968898E-16;
long double Higgs_width_Poly_Fit_Zone2_coeff0 = 2.563291882845993E+02;
long double Higgs_width_Poly_Fit_Zone2_coeff1 = -1.037082025855304E+01;
long double Higgs_width_Poly_Fit_Zone2_coeff2 = 1.780260502696301E-01;
long double Higgs_width_Poly_Fit_Zone2_coeff3 = -1.720311784419889E-03;
long double Higgs_width_Poly_Fit_Zone2_coeff4 = 1.038418605369741E-05;
long double Higgs_width_Poly_Fit_Zone2_coeff5 = -4.092496883922424E-08;
long double Higgs_width_Poly_Fit_Zone2_coeff6 = 1.067667966800388E-10;
long double Higgs_width_Poly_Fit_Zone2_coeff7 = -1.823343280081685E-13;
long double Higgs_width_Poly_Fit_Zone2_coeff8 = 1.955637395597351E-16;
long double Higgs_width_Poly_Fit_Zone2_coeff9 = -1.193287048560413E-19;
long double Higgs_width_Poly_Fit_Zone2_coeff10 = 3.156196649452213E-23;
long double Higgs_width_Poly_Fit_Zone3_coeff0 = -5.255605465437446E+02;
long double Higgs_width_Poly_Fit_Zone3_coeff1 = 1.036972988796150E+01;
long double Higgs_width_Poly_Fit_Zone3_coeff2 = -6.817022987365029E-02;
long double Higgs_width_Poly_Fit_Zone3_coeff3 = 1.493275723660056E-04;
long double m_m4l__2 = m_m4l * m_m4l;
long double m_m4l__3 = m_m4l__2 * m_m4l;
if (m_m4l < 156.5)
return (Higgs_width_Poly_Fit_Zone1_coeff0 + Higgs_width_Poly_Fit_Zone1_coeff1 * m_m4l
+ Higgs_width_Poly_Fit_Zone1_coeff2 * m_m4l__2
+ Higgs_width_Poly_Fit_Zone1_coeff3 * m_m4l__3
+ Higgs_width_Poly_Fit_Zone1_coeff4 * m_m4l__2 * m_m4l__2
+ Higgs_width_Poly_Fit_Zone1_coeff5 * m_m4l__2 * m_m4l__3
+ Higgs_width_Poly_Fit_Zone1_coeff6 * m_m4l__2 * m_m4l__2 * m_m4l__2
+ Higgs_width_Poly_Fit_Zone1_coeff7 * m_m4l__2 * m_m4l__2 * m_m4l__3
+ Higgs_width_Poly_Fit_Zone1_coeff8 * m_m4l__2 * m_m4l__2 * m_m4l__2 * m_m4l__2
+ Higgs_width_Poly_Fit_Zone1_coeff9 * m_m4l__2 * m_m4l__2 * m_m4l__2 * m_m4l__3);
else if (m_m4l >= 156.5 && m_m4l <= 162)
return (Higgs_width_Poly_Fit_Zone3_coeff0 + Higgs_width_Poly_Fit_Zone3_coeff1 * m_m4l
+ Higgs_width_Poly_Fit_Zone3_coeff2 * m_m4l__2
+ Higgs_width_Poly_Fit_Zone3_coeff3 * m_m4l__3);
else
return (
Higgs_width_Poly_Fit_Zone2_coeff0 + Higgs_width_Poly_Fit_Zone2_coeff1 * m_m4l
+ Higgs_width_Poly_Fit_Zone2_coeff2 * m_m4l__2 + Higgs_width_Poly_Fit_Zone2_coeff3 * m_m4l__3
+ Higgs_width_Poly_Fit_Zone2_coeff4 * m_m4l__2 * m_m4l__2
+ Higgs_width_Poly_Fit_Zone2_coeff5 * m_m4l__2 * m_m4l__3
+ Higgs_width_Poly_Fit_Zone2_coeff6 * m_m4l__2 * m_m4l__2 * m_m4l__2
+ Higgs_width_Poly_Fit_Zone2_coeff7 * m_m4l__2 * m_m4l__2 * m_m4l__3
+ Higgs_width_Poly_Fit_Zone2_coeff8 * m_m4l__2 * m_m4l__2 * m_m4l__2 * m_m4l__2
+ Higgs_width_Poly_Fit_Zone2_coeff9 * m_m4l__2 * m_m4l__2 * m_m4l__2 * m_m4l__3
+ Higgs_width_Poly_Fit_Zone2_coeff10 * m_m4l__2 * m_m4l__2 * m_m4l__2 * m_m4l__2 * m_m4l__2);
}
}; // class Parameters_heft
// calculate LO ME
class CPPProcess_P0_Sigma_heft_pp_H_ZZ_4l_heft_gg_epemmupmum {
public:
// Constructor.
CPPProcess_P0_Sigma_heft_pp_H_ZZ_4l_heft_gg_epemmupmum()
: pars(), pout(4, vector<double>(4, 0.)) { }
// Destructor.
virtual ~CPPProcess_P0_Sigma_heft_pp_H_ZZ_4l_heft_gg_epemmupmum() { }
float Compute(const Quadruplet& quad) {
FourMomentum cms = quad.mom();
// use the cms mass as a value for MH
pars.set4lepMass(cms.mass());
const FourMomentum* fermionsMom[4] = {&quad.Z1().first.momentum(), &quad.Z1().second.momentum(),
&quad.Z2().first.momentum(), &quad.Z2().second.momentum()};
// boost to center-of-mass frame
LorentzTransform HRF_boost = LorentzTransform::mkFrameTransformFromBeta(cms.betaVec());
for (std::size_t i = 0; i < 4; ++i) {
FourMomentum tmpMom = HRF_boost.transform(*fermionsMom[i]);
pout[i][0] = tmpMom.E();
pout[i][1] = tmpMom.px();
pout[i][2] = tmpMom.py();
pout[i][3] = tmpMom.pz();
}
// Evaluate matrix element
return sigmaKin();
}
private:
// Calculate flavour-independent parts of cross section. Evaluate |M|^2, part independent of
// incoming flavour. Return matrix element
double sigmaKin() {
// Local variables and constants
static const int ncomb = 16;
// Helicities for the process
static const int helicities[ncomb][4] = {{-1, -1, -1, -1}, {-1, -1, -1, 1}, {-1, -1, 1, -1},
{-1, -1, 1, 1}, {-1, 1, -1, -1}, {-1, 1, -1, 1},
{-1, 1, 1, -1}, {-1, 1, 1, 1}, {1, -1, -1, -1},
{1, -1, -1, 1}, {1, -1, 1, -1}, {1, -1, 1, 1},
{1, 1, -1, -1}, {1, 1, -1, 1}, {1, 1, 1, -1},
{1, 1, 1, 1}};
// Reset the matrix elements
double matrix_element = 0.;
// Calculate the matrix element for all helicities
for (int ihel = 0; ihel < ncomb; ihel++) {
matrix_element += matrix_gg_h_h_zz_z_epem_z_mupmum(helicities[ihel]);
}
// Denominators: spins, colors and identical particles
return matrix_element /= 128.;
}
// Calculate wavefunctions and matrix elements for all subprocesses
double matrix_gg_h_h_zz_z_epem_z_mupmum(const int hel[]) {
// Calculate all wavefunctions
ixx(pout[0], hel[0], w[2], true);
ixx(pout[1], hel[1], w[3], false);
FFV2_4_3(w[2], w[3], pars.GC_40, pars.GC_54, pars.mdl_MZ, pars.mdl_WZ, w[4]);
ixx(pout[2], hel[2], w[5], true);
ixx(pout[3], hel[3], w[6], false);
FFV2_4_3(w[5], w[6], pars.GC_40, pars.GC_54, pars.mdl_MZ, pars.mdl_WZ, w[7]);
VVS2_3(w[4], w[7], pars.GC_73, pars.mdl_MH, pars.mdl_WH, w[8]);
// Calculate all amplitudes
// Amplitude(s) for diagram number 0
std::complex<double> amp = VVS3_0(pars.mdl_MH, w[8], pars.GC_13);
// Calculate color flows
// Sum and square the color flows to get the matrix element
return real(8. * amp * conj(amp));
}
// wave function,
// ixx true takes anti-particle, false takes particle
void ixx(std::vector<double> p, int nhel, std::complex<double> fi[6], bool isixx) {
std::complex<double> chi[2];
double sqp0p3;
fi[0] = std::complex<double>(p[0], p[3]);
fi[1] = std::complex<double>(p[1], p[2]);
if (p[1] == 0.0 and p[2] == 0.0 and p[3] < 0.0)
sqp0p3 = 0.0;
else
sqp0p3 = pow(max(p[0] + p[3], 0.0), 0.5);
if (isixx)
chi[0] = std::complex<double>(-sqp0p3, 0.0);
else
chi[0] = std::complex<double>(sqp0p3, 0.0);
if (sqp0p3 == 0.0)
chi[1] = std::complex<double>(-nhel * pow(2.0 * p[0], 0.5), 0.0);
else
chi[1] = std::complex<double>(nhel * p[1], -p[2]) / sqp0p3;
if (isixx) {
if (nhel == 1) {
fi[2] = chi[1];
fi[3] = chi[0];
fi[4] = std::complex<double>(0.0, 0.0);
fi[5] = std::complex<double>(0.0, 0.0);
}
else {
fi[2] = std::complex<double>(0.0, 0.0);
fi[3] = std::complex<double>(0.0, 0.0);
fi[4] = chi[0];
fi[5] = chi[1];
}
}
else {
if (nhel == 1) {
fi[2] = chi[0];
fi[3] = chi[1];
fi[4] = std::complex<double>(0.00, 0.00);
fi[5] = std::complex<double>(0.00, 0.00);
}
else {
fi[2] = std::complex<double>(0.00, 0.00);
fi[3] = std::complex<double>(0.00, 0.00);
fi[4] = chi[1];
fi[5] = chi[0];
}
}
return;
}
// vertices
void VVS2_3(std::complex<double> V1[],
std::complex<double> V2[],
std::complex<double> COUP,
double M3,
double W3,
std::complex<double> S3[]) {
std::complex<double> cI = std::complex<double>(0., 1.);
std::complex<double> TMP1;
double P3[4];
std::complex<double> denom;
S3[0] = +V1[0] + V2[0];
S3[1] = +V1[1] + V2[1];
P3[0] = -S3[0].real();
P3[1] = -S3[1].real();
P3[2] = -S3[1].imag();
P3[3] = -S3[0].imag();
TMP1 = (V2[2] * V1[2] - V2[3] * V1[3] - V2[4] * V1[4] - V2[5] * V1[5]);
denom = COUP / (pow(P3[0], 2) - pow(P3[1], 2) - pow(P3[2], 2) - pow(P3[3], 2) - M3 * (M3 - cI * W3));
S3[2] = denom * cI * TMP1;
}
void FFV2_4_3(std::complex<double> F1[],
std::complex<double> F2[],
std::complex<double> COUP1,
std::complex<double> COUP2,
double M3,
double W3,
std::complex<double> V3[]) {
std::complex<double> cI = std::complex<double>(0., 1.);
std::complex<double> denom;
std::complex<double> TMP11;
double P3[4];
std::complex<double> TMP14;
double OM3 = 1. / pow(M3, 2);
V3[0] = +F1[0] + F2[0];
V3[1] = +F1[1] + F2[1];
P3[0] = -V3[0].real();
P3[1] = -V3[1].real();
P3[2] = -V3[1].imag();
P3[3] = -V3[0].imag();
TMP14 = (F1[4] * (F2[2] * (P3[0] - P3[3]) - F2[3] * (P3[1] + cI * (P3[2])))
+ F1[5] * (F2[2] * (+cI * (P3[2]) - P3[1]) + F2[3] * (P3[0] + P3[3])));
TMP11 = (F1[2] * (F2[4] * (P3[0] + P3[3]) + F2[5] * (P3[1] + cI * (P3[2])))
+ F1[3] * (F2[4] * (P3[1] - cI * (P3[2])) + F2[5] * (P3[0] - P3[3])));
denom = 1. / (pow(P3[0], 2) - pow(P3[1], 2) - pow(P3[2], 2) - pow(P3[3], 2) - M3 * (M3 - cI * W3));
V3[2] = COUP2 * denom * -2. * cI
* (OM3 * -1. / 2. * P3[0] * (TMP11 + 2. * (TMP14))
+ (+1. / 2. * (F2[4] * F1[2] + F2[5] * F1[3]) + F2[2] * F1[4] + F2[3] * F1[5]));
V3[3] = COUP2 * denom * -2. * cI
* (OM3 * -1. / 2. * P3[1] * (TMP11 + 2. * (TMP14))
+ (-1. / 2. * (F2[5] * F1[2] + F2[4] * F1[3]) + F2[3] * F1[4] + F2[2] * F1[5]));
V3[4] = COUP2 * denom * 2. * cI
* (OM3 * 1. / 2. * P3[2] * (TMP11 + 2. * (TMP14))
+ (+1. / 2. * cI * (F2[5] * F1[2]) - 1. / 2. * cI * (F2[4] * F1[3]) - cI * (F2[3] * F1[4])
+ cI * (F2[2] * F1[5])));
V3[5] = COUP2 * denom * 2. * cI
* (OM3 * 1. / 2. * P3[3] * (TMP11 + 2. * (TMP14))
+ (+1. / 2. * (F2[4] * F1[2]) - 1. / 2. * (F2[5] * F1[3]) - F2[2] * F1[4] + F2[3] * F1[5]));
V3[2] += COUP1 * denom * -cI * (F2[4] * F1[2] + F2[5] * F1[3] - P3[0] * OM3 * TMP11);
V3[3] += COUP1 * denom * -cI * (-F2[5] * F1[2] - F2[4] * F1[3] - P3[1] * OM3 * TMP11);
V3[4] += COUP1 * denom * -cI * (-cI * (F2[5] * F1[2]) + cI * (F2[4] * F1[3]) - P3[2] * OM3 * TMP11);
V3[5] += COUP1 * denom * -cI * (F2[5] * F1[3] - F2[4] * F1[2] - P3[3] * OM3 * TMP11);
}
std::complex<double> VVS3_0(double mass, std::complex<double> S3[], std::complex<double> COUP) {
std::complex<double> TMP15;
TMP15 = std::complex<double>(0., pow(mass, 2) / 2.);
return COUP * S3[2] * (TMP15);
}
static const int nwavefuncs = 9;
std::complex<double> w[nwavefuncs][18];
// Pointer to the model parameters
Parameters_heft pars;
// vector with momenta (to be changed each event)
std::vector<std::vector<double>> pout;
};
CPPProcess_P0_Sigma_heft_pp_H_ZZ_4l_heft_gg_epemmupmum MGME;
};
RIVET_DECLARE_PLUGIN(ATLAS_2020_I1790439);
} ```