rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

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: 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
  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}