Rivet analyses
title: ZEUS_2021_I1869927
Azimuthal correlations in photoproduction and DIS at ZEUS
Experiment: ZEUS (HERA)
Inspire ID: 1869927
Status: VALIDATED
Authors: - Ilkka Helenius
References: - arXiv: 2106.12377 - JHEP 12 (2021) 102
Beams: p+ e-, p+ e+, e+ p+, e- p+
Beam energies: (920.0, 27.5); (920.0, 27.5); (27.5, 920.0); (27.5, 920.0)GeV
Run details: - Particle multiplicities, $p_T$ spectrum and correlations for particles with $0.1 < p_T < 5.0$ GeV in electron-proton collisions at HERA, both DIS ($Q^2>5~\mathrm{GeV^2}$) and photoproduction ($Q^2<1~\mathrm{GeV^2}$).
Collective behaviour of final-state hadrons, and multiparton interactions are studied in high-multiplicity ep scattering at a centre-of-mass energy $\sqrt{s} = 318~\mathrm{GeV}$ with the ZEUS detector at HERA. Two- and four-particle azimuthal correlations, as well as multiplicity, transverse momentum, and pseudorapidity distributions for charged-particle multiplicities $N_{\mathrm{ch}} \geq 20$ are measured.
Source code:ZEUS_2021_I1869927.cc
```c++ // -- C++ --
include "Rivet/Analysis.hh"
include "Rivet/Projections/ChargedFinalState.hh"
include "Rivet/Projections/DISKinematics.hh"
include "Rivet/Projections/FinalState.hh"
include "Rivet/Projections/PrimaryParticles.hh"
namespace Rivet {
/// @brief Azimuthal correlations in photoproduction and DIS at ZEUS class ZEUS_2021_I1869927 : public Analysis { public:
/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(ZEUS_2021_I1869927);
/// @name Analysis methods
///@{
// Remove particles if they were decay products of these particles
bool unwantedDecay(const Particle& p) const {
for (PdgId decayID : parentDecayIDs) {
if (p.hasParentWith(Cuts::pid == decayID)) return true;
}
return false;
}
/// Book histograms and initialise projections before the run
void init() {
// Particles to be used for analysis
// 0.1 < pT < 5.0 GeV, -1.5 < eta < 2.0
// Charged pions, kaons, protons, Xi, Sigma, Omega
const ChargedFinalState cfs(
Cuts::pT > 0.1 * GeV && Cuts::pT < 5.0 * GeV && Cuts::eta > -2.0 && Cuts::eta < 2.0
&& (Cuts::abspid == 211 || Cuts::abspid == 321 || Cuts::abspid == 2212 || Cuts::abspid == 3312
|| Cuts::abspid == 3222 || Cuts::abspid == 3112 || Cuts::abspid == 3334));
// Final-state and kinematics projections
declare(cfs, "SelectedParticles");
declare(DISKinematics(), "Kinematics");
// Counter for calculating the sum of weights within the acceptance
book(_c["sow_c1"], "_sow_c1");
book(_c["sow_c2"], "_sow_c2");
book(_c["sow"], "_sow");
// Virtuality dependent c_1{2}
book(_r["c12_0"], 1, 1, 1);
book(_r["c12_1"], 2, 1, 1);
book(_r["c12_2"], 3, 1, 1);
book(_r["c12_3"], 4, 1, 1);
// Virtuality dependent c_2{2}
book(_r["c22_0"], 5, 1, 1);
book(_r["c22_1"], 6, 1, 1);
book(_r["c22_2"], 7, 1, 1);
book(_r["c22_3"], 8, 1, 1);
// Histograms used for bin-by-bin average calculation for the above estimates
book(_h["c12_0_num"], "_c12_0_num", _r["c12_0"].binning().edges<0>());
book(_h["c12_1_num"], "_c12_1_num", _r["c12_0"].binning().edges<0>());
book(_h["c12_2_num"], "_c12_2_num", _r["c12_0"].binning().edges<0>());
book(_h["c12_3_num"], "_c12_3_num", _r["c12_0"].binning().edges<0>());
book(_h["c22_0_num"], "_c22_0_num", _r["c22_0"].binning().edges<0>());
book(_h["c22_1_num"], "_c22_1_num", _r["c22_0"].binning().edges<0>());
book(_h["c22_2_num"], "_c22_2_num", _r["c22_0"].binning().edges<0>());
book(_h["c22_3_num"], "_c22_3_num", _r["c22_0"].binning().edges<0>());
book(_h["c_0_den"], "_c_0_den", _r["c22_0"].binning().edges<0>());
book(_h["c_1_den"], "_c_1_den", _r["c22_0"].binning().edges<0>());
book(_h["c_2_den"], "_c_2_den", _r["c22_0"].binning().edges<0>());
book(_h["c_3_den"], "_c_3_den", _r["c22_0"].binning().edges<0>());
// Normalized single-particle distributions
book(_h_Nch, 9, 1, 1);
book(_h_pT, 10, 1, 1);
book(_h_eta, 11, 1, 1);
// Further two-particle correlations
book(_r["c2eta0"], 12, 1, 1);
book(_r["c2eta1"], 13, 1, 1);
book(_r["c2pT0"], 14, 1, 1);
book(_r["c2pT1"], 15, 1, 1);
// Histograms to allow bin-by-bin averages with non-unit event weights
book(_h["c2eta0_num"], "_c2eta0_num", _r["c2eta0"].binning().edges<0>());
book(_h["c2eta1_num"], "_c2eta1_num", _r["c2eta1"].binning().edges<0>());
book(_h["c2eta_den"], "_c2eta_den", _r["c2eta0"].binning().edges<0>());
book(_h["c2pT0_num"], "_c2pT0_num", _r["c2pT0"].binning().edges<0>());
book(_h["c2pT1_num"], "_c2pT1_num", _r["c2pT1"].binning().edges<0>());
book(_h["c2pT_den"], "_c2pT_den", _r["c2pT0"].binning().edges<0>());
// Four-particle correlations
book(_r["c14pT1"], 16, 1, 1);
book(_r["c24pT1"], 17, 1, 1);
// Additional histograms are need for c4 calculations, use the same binning
book(_r["c12pT1"], "_c12pT1", _h_pT.binning().edges<0>());
book(_r["c22pT1"], "_c22pT1", _h_pT.binning().edges<0>());
book(_h["c12pT1_num"], "_c12pT1_num", _h_pT.binning().edges<0>());
book(_h["c22pT1_num"], "_c22pT1_num", _h_pT.binning().edges<0>());
book(_h["c12pT1_den"], "_c12pT1_den", _h_pT.binning().edges<0>());
book(_h["c14pT1_num"], "_c14pT1_num", _h_pT.binning().edges<0>());
book(_h["c24pT1_num"], "_c24pT1_num", _h_pT.binning().edges<0>());
book(_h["c14pT1_den"], "_c14pT1_den", _h_pT.binning().edges<0>());
}
/// Perform the per-event analysis
void analyze(const Event& event) {
// Use scattered lepton momentum to derive orientation and kinematics
const DISKinematics& kin = apply<DISKinematics>(event, "Kinematics");
int orientation = kin.orientation();
double Q2 = kin.Q2() * GeV2;
// Cut out the intermediate Q2 range between photoproduction and DIS excluded from the analysis
if (inRange(Q2, 1.0, 5.0)) vetoEvent;
// Final-state particles
const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "SelectedParticles");
// Find particles from considered decays within acceptance and calculate multiplicity
// Acceptance -1.5 < eta < 2.0 with orientation == 1, otherwise reversed
int Nch = 0;
vector<Particle> particles;
for (const Particle& p : cfs.particles()) {
if (unwantedDecay(p)) continue;
if ((orientation * p.eta() < -1.5)) continue;
particles.push_back(p);
++Nch;
}
// Cut on Nch, fill discrete histogram for photoproduction
if (Nch < 20) vetoEvent;
if (inRange(Q2, 0. * GeV2, 1. * GeV2)) _h_Nch->fill(Nch);
// Calculate the correlations
// 1st particle loop
for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
const Particle p1 = particles[ip1];
// Fill in single-particle quantities
if (inRange(Q2, 0. * GeV2, 1. * GeV2)) {
_h_pT->fill(p1.pT() / GeV);
_h_eta->fill(orientation * p1.eta());
}
// 2nd particle loop for two-particle correlations, now the unwanted decays already removed.
for (size_t ip2 = 0; ip2 < particles.size(); ++ip2) {
if (ip2 == ip1) continue;
// Calculate two-particle angular correlations
const Particle p2 = particles[ip2];
double dEta = fabs(p1.eta() - p2.eta());
double pTave = 0.5 * (p1.pT() + p2.pT());
// Calculate cosines of delta phi.
double dPhi = p1.phi() - p2.phi();
double c1 = cos(dPhi);
double c2 = cos(2 * dPhi);
// Virtuality dependence without pT or gap conditions
_h["c12_0_num"]->fill(Q2 / GeV2, c1);
_h["c22_0_num"]->fill(Q2 / GeV2, c2);
_h["c_0_den"]->fill(Q2 / GeV2);
// Check the event classification
bool largeGap = (dEta > 2.0) ? true : false;
bool highpT = ((p1.pT() > 0.5 * GeV) && (p2.pT() > 0.5 * GeV)) ? true : false;
// Fill in the Q^2-dependent histograms
if (largeGap) {
_h["c12_1_num"]->fill(Q2 / GeV2, c1);
_h["c22_1_num"]->fill(Q2 / GeV2, c2);
_h["c_1_den"]->fill(Q2 / GeV2);
}
if (highpT) {
_h["c12_2_num"]->fill(Q2 / GeV2, c1);
_h["c22_2_num"]->fill(Q2 / GeV2, c2);
_h["c_2_den"]->fill(Q2 / GeV2);
}
if (largeGap && highpT) {
_h["c12_3_num"]->fill(Q2 / GeV2, c1);
_h["c22_3_num"]->fill(Q2 / GeV2, c2);
_h["c_3_den"]->fill(Q2 / GeV2);
}
// Only photoproduction from this on
if (Q2 > 1.0 * GeV2) continue;
// Fill two-particle correlation histograms
_h["c2eta0_num"]->fill(dEta, c1);
_h["c2eta1_num"]->fill(dEta, c2);
_h["c2eta_den"]->fill(dEta);
_h["c2pT0_num"]->fill(pTave, c1);
_h["c2pT1_num"]->fill(pTave, c2);
_h["c2pT_den"]->fill(pTave);
// Additional histograms used for c_4 calculations
_h["c12pT1_num"]->fill(p1.pT() / GeV, c1);
_h["c22pT1_num"]->fill(p1.pT() / GeV, c2);
_h["c12pT1_den"]->fill(p1.pT() / GeV);
_c["sow_c1"]->fill(c1);
_c["sow_c2"]->fill(c2);
_c["sow"]->fill(1.);
// 3rd particle loop for two-particle correlations.
for (size_t ip3 = 0; ip3 < particles.size(); ++ip3) {
if (ip3 == ip1) continue;
if (ip3 == ip2) continue;
const Particle p3 = particles[ip3];
// 4th particle loop for two-particle correlations.
for (size_t ip4 = 0; ip4 < particles.size(); ++ip4) {
if (ip4 == ip1) continue;
if (ip4 == ip2) continue;
if (ip4 == ip3) continue;
const Particle p4 = particles[ip4];
// Calculate delta-phi and cosines
double dPhi4 = p1.phi() + p2.phi() - p3.phi() - p4.phi();
double c14 = cos(dPhi4);
double c24 = cos(2 * dPhi4);
// Fill final histograms
_h["c14pT1_num"]->fill(p1.pT() / GeV, c14);
_h["c24pT1_num"]->fill(p1.pT() / GeV, c24);
_h["c14pT1_den"]->fill(p1.pT() / GeV);
}
}
}
}
}
/// Normalise histograms etc., after the run
void finalize() {
// Normalized single-particle quantities, scale pT and eta with bin widhts
normalize(_h_Nch);
normalize({_h_pT, _h_eta});
for (auto& item : {_h_pT, _h_eta}) {
for (auto& b : item->bins()) {
b.scaleW(b.xWidth());
}
}
// Calculate bin-averaged correlators
divide(_h["c12_0_num"], _h["c_0_den"], _r["c12_0"]);
divide(_h["c12_1_num"], _h["c_1_den"], _r["c12_1"]);
divide(_h["c12_2_num"], _h["c_2_den"], _r["c12_2"]);
divide(_h["c12_3_num"], _h["c_3_den"], _r["c12_3"]);
divide(_h["c22_0_num"], _h["c_0_den"], _r["c22_0"]);
divide(_h["c22_1_num"], _h["c_1_den"], _r["c22_1"]);
divide(_h["c22_2_num"], _h["c_2_den"], _r["c22_2"]);
divide(_h["c22_3_num"], _h["c_3_den"], _r["c22_3"]);
divide(_h["c2eta0_num"], _h["c2eta_den"], _r["c2eta0"]);
divide(_h["c2eta1_num"], _h["c2eta_den"], _r["c2eta1"]);
divide(_h["c2pT0_num"], _h["c2pT_den"], _r["c2pT0"]);
divide(_h["c2pT1_num"], _h["c2pT_den"], _r["c2pT1"]);
// Calculate c_n{4}(pT) as C_n{4}(pT) - 2*c_n{2}(pT)*c_n{2}
divide(_h["c12pT1_num"], _h["c12pT1_den"], _r["c12pT1"]);
divide(_h["c22pT1_num"], _h["c12pT1_den"], _r["c22pT1"]);
divide(_h["c14pT1_num"], _h["c14pT1_den"], _r["c14pT1"]);
divide(_h["c24pT1_num"], _h["c14pT1_den"], _r["c24pT1"]);
_r["c12pT1"]->scale(2.0 * _c["sow_c1"]->sumW() / _c["sow"]->sumW());
_r["c22pT1"]->scale(2.0 * _c["sow_c2"]->sumW() / _c["sow"]->sumW());
*_r["c14pT1"] -= *_r["c12pT1"];
*_r["c24pT1"] -= *_r["c22pT1"];
}
/// @}
/// @name Histograms, estimates and counters
/// @{
// Single-particle quantities as simple histograms
Histo1DPtr _h_pT, _h_eta;
BinnedHistoPtr<int> _h_Nch;
// Estimates to store final result, histogram and counters to account bin-by-bin averaging of correlations
map<string, Estimate1DPtr> _r;
map<string, Histo1DPtr> _h;
map<string, CounterPtr> _c;
// Do not consider decay products of these particles.
const vector<PdgId> parentDecayIDs = {PID::PROTON, PID::PHOTON, PID::K0, PID::ELECTRON,
PID::NEUTRON, PID::MUON, PID::K0L, PID::PIPLUS,
PID::KPLUS, PID::XI0, PID::LAMBDA, PID::XIMINUS,
PID::SIGMAMINUS, PID::K0S, PID::OMEGAMINUS, PID::SIGMAPLUS};
/// @}
};
RIVET_DECLARE_PLUGIN(ZEUS_2021_I1869927);
} ```