Rivet analyses
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, pT spectrum and correlations for particles with 0.1 < pT < 5.0 GeV in electron-proton collisions at HERA, both DIS (Q2 > 5 GeV2) and photoproduction (Q2 < 1 GeV2).
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 Nch ≥ 20 are measured.
Source
code:ZEUS_2021_I1869927.cc
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/ChargedFinalState.hh"
#include "Rivet/Projections/PrimaryParticles.hh"
#include "Rivet/Projections/DISKinematics.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);
}