Rivet analyses
title: H1_1996_I424463
Transverse momentum spectra of charged particles in DIS
Experiment: H1 (HERA)
Inspire ID: 424463
Status: VALIDATED
Authors: - Suraj Kumar Singh - Hannes Jung - Andrii Verbytskyi
References: - Nucl. Phys. B485 (1997) 3 - DOI: 10.1016/S0550-3213(96)00675-X - arXiv: hep-ex/9610006 - DESY-96-215
Beams: e+ p+, p+ e+
Beam energies: (27.5, 820.0); (820.0, 27.5)GeV
Run details: none listed
Transverse momentum spectra of charged particles produced in deep inelastic scattering are measured as a function of the kinematic variables $x$ and $Q^2$ using the H1 detector at the ep collider HERA. The data are compared to different parton-emission models, either with or without ordering of the emissions in transverse momentum. The data provide evidence for a relatively large amount of parton radiation between the current and the remnant systems.
Source code:H1_1996_I424463.cc
```c++ // -- C++ --
include "Rivet/Analysis.hh"
include "Rivet/Projections/ChargedFinalState.hh"
include "Rivet/Projections/DISKinematics.hh"
include "Rivet/Projections/FinalState.hh"
include "Rivet/Tools/ParticleIdUtils.hh"
namespace Rivet {
/// @brief Transverse momentum spectra of charged particles in DIS class H1_1996_I424463 : public Analysis { public:
/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(H1_1996_I424463);
/// @name Analysis methods
/// @{
/// Book histograms and initialise projections before the run
void init() {
// Book projections
declare(DISLepton(), "Lepton");
declare(DISKinematics(), "Kinematics");
declare(ChargedFinalState(), "CFS");
declare(FinalState(), "FS");
// Book histograms
book(_NevAll, "TMP/Nev_all");
int ixx = 0;
for (size_t ix = 0; ix < 10; ++ix) {
book(_Nevt_after_cuts[ix], "TMP/Nevt_after_cuts" + to_string(ix));
for (unsigned int ih = 0; ih < 2; ++ih) {
book(_h_dndpt_eta_bin[ih][ix], ih * 10 + ix + 1, 1, 1);
book(_hdndeta_bin[ih][ix], ih * 10 + ix + 29, 1, 1);
}
if (ix != 6 && ix != 9) {
book(_hdndptmax_low_eta_bin[ixx], ixx + 21, 1, 1);
ixx = ixx + 1;
}
}
}
/// Perform the per-event analysis
void analyze(const Event& event) {
if (_edgespT.empty()) {
_edgespT = _h_dndpt_eta_bin[0][0]->xEdges();
_edgesEta = _hdndeta_bin[0][0]->xEdges();
_edgespTMax = _hdndptmax_low_eta_bin[0]->xEdges();
}
const FinalState& fs = apply<FinalState>(event, "FS");
const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
const DISLepton& dl = apply<DISLepton>(event, "Lepton");
// Get the DIS kinematics
double x = dk.x();
double y = dk.y();
double Q2 = dk.Q2() / GeV;
double W2 = dk.W2() / GeV;
// Momentum of the scattered lepton
FourMomentum leptonMom = dl.out().momentum();
double enel = leptonMom.E() / GeV;
double thel = 180. - leptonMom.angle(dl.in().momentum()) / degree;
_NevAll->fill();
bool cut = y > 0.05 && Q2 > 5. && Q2 < 100. && enel > 12. && W2 > 4400. && thel > 157. && thel < 173.;
if (!cut) vetoEvent;
int ibin[10];
for (int i = 0; i < 10; i++) ibin[i] = 0;
if (5. < Q2 && Q2 < 50. && x > 0.0001 && x < 0.0010) ibin[0] = 1;
if (5. < Q2 && Q2 < 10. && x > 0.0001 && x < 0.0002) ibin[1] = 1;
if (6. < Q2 && Q2 < 10. && x > 0.0002 && x < 0.0005) ibin[2] = 1;
if (10. < Q2 && Q2 < 20. && x > 0.0002 && x < 0.0005) ibin[3] = 1;
if (10. < Q2 && Q2 < 20. && x > 0.0005 && x < 0.0008) ibin[4] = 1;
if (10. < Q2 && Q2 < 20. && x > 0.0008 && x < 0.0015) ibin[5] = 1;
if (10. < Q2 && Q2 < 20. && x > 0.0015 && x < 0.0040) ibin[6] = 1;
if (20. < Q2 && Q2 < 50. && x > 0.0005 && x < 0.0014) ibin[7] = 1;
if (20. < Q2 && Q2 < 50. && x > 0.0014 && x < 0.0030) ibin[8] = 1;
if (20. < Q2 && Q2 < 50. && x > 0.0030 && x < 0.0100) ibin[9] = 1;
for (int i = 0; i < 10; i++) {
if (ibin[i] == 1) _Nevt_after_cuts[i]->fill();
}
// Extract the particles other than the lepton
/// @todo Improve to avoid HepMC digging
Particles particles;
particles.reserve(fs.particles().size());
ConstGenParticlePtr dislepGP = dl.out().genParticle();
for (const Particle& p : fs.particles()) {
ConstGenParticlePtr loopGP = p.genParticle();
if (loopGP == dislepGP) continue;
particles.push_back(p);
}
// Boost to hadronic CM
const LorentzTransform hcmboost = dk.boostHCM();
// Loop over the particles
int mult = 0;
double ptmax_high[10], ptmax_low[10];
for (int i = 0; i < 10; i++) {
ptmax_high[i] = 0.;
ptmax_low[i] = 0.;
}
double EtSum = 0;
for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
const Particle& p = particles[ip1];
double eta = p.momentum().pseudorapidity();
// Boost to hcm
const FourMomentum hcmMom = hcmboost.transform(p.momentum());
// Apply safety cuts
if (eta > -5 && eta < 10.) {
mult = mult + 1;
double pThcm = hcmMom.pT();
double etahcm = hcmMom.pseudorapidity();
if (etahcm > 0. && etahcm < 2.0) {
EtSum = EtSum + hcmMom.Et();
}
if (PID::charge(p.pid()) != 0) {
if (etahcm > 0.5 && etahcm < 1.5) {
for (int i = 0; i < 10; i++) {
if (ibin[i] == 1) {
fillpT(1, i, pThcm);
if (pThcm > ptmax_low[i]) ptmax_low[i] = pThcm;
}
}
}
if (etahcm > 1.5 && etahcm < 2.5) {
for (int i = 0; i < 10; i++) {
if (ibin[i] == 1) {
fillpT(0, i, pThcm);
if (pThcm > ptmax_high[i]) ptmax_high[i] = pThcm;
}
}
}
for (int i = 0; i < 10; i++) {
if (ibin[i] == 1) fillEta(1, i, etahcm);
if (ibin[i] == 1 && pThcm > 1.) fillEta(0, i, etahcm);
}
}
} // end of loop over the particles
}
int ii = 0;
for (int i = 0; i < 10; i++) {
if (i != 6 && i != 9) {
if (ibin[i] == 1 && EtSum > 6.) {
fillpTMax(ii, ptmax_low[i]);
}
ii = ii + 1;
}
}
}
/// Normalise histograms etc., after the run
void finalize() {
MSG_DEBUG("All events: " << _NevAll->val() << " after cuts: " << _Nevt_after_cuts[0]->val());
MSG_DEBUG("Cut1 events: " << _NevAll->val() << " after cuts: " << _Nevt_after_cuts[1]->val());
int ii = 0;
for (int i = 0; i < 10; i++) {
if (_Nevt_after_cuts[i]->val() != 0) {
for (unsigned int ih = 0; ih < 2; ++ih) {
scale(_h_dndpt_eta_bin[ih][i], 1. / *_Nevt_after_cuts[i]);
size_t ioff = 0;
if ((ih == 0 && (i == 6 || i == 9)) || (ih == 1 && (i <= 4 || i == 7)))
ioff = 1;
else if (ih == 1 && (i == 5 || i == 8))
ioff = 2;
for (auto& b : _h_dndpt_eta_bin[ih][i]->bins()) {
const size_t idx = b.index() + ioff;
b.scaleW(1. / _axispT.width(idx));
}
scale(_hdndeta_bin[ih][i], 1. / *_Nevt_after_cuts[i]);
ioff = 0;
if ((ih == 0 && (i == 5 || i == 6 || i == 8 || i == 9))
|| (ih == 1 && (i == 0 || i == 4 || i == 5 || i == 8)))
ioff = 1;
else if (ih == 1 && (i == 6 || i == 9))
ioff = 2;
for (auto& b : _hdndeta_bin[ih][i]->bins()) {
const size_t idx = b.index() + ioff;
b.scaleW(1. / _axisEta.width(idx));
}
}
}
if (i != 6 && i != 9) {
if (_Nevt_after_cuts[i]->val() != 0) {
normalize(_hdndptmax_low_eta_bin[ii]);
for (auto& b : _hdndptmax_low_eta_bin[ii]->bins()) {
const size_t idx = b.index();
b.scaleW(1. / _axispTMax.width(idx));
}
}
ii = ii + 1;
}
}
}
/// @}
void fillEta(const unsigned int ix, const unsigned int iy, const double value) {
string edge = "OTHER";
const size_t idx = _axisEta.index(value);
if (idx && idx <= _edgesEta.size()) {
if (((ix == 0 && iy >= 1 && iy <= 3) || (ix == 0 && (iy == 1 || iy == 3))) && idx == _edgesEta.size())
;
else if (idx == 1
&& ((ix == 0 && (iy == 5 || iy == 6 || iy == 8 || iy == 9))
|| (ix == 1 && (iy == 0 || iy == 4 || iy == 5 || iy == 8))))
;
else if (idx <= 2 && (ix == 1 && (iy == 6 || iy == 9)))
;
else {
edge = _edgesEta[idx - 1];
}
}
_hdndeta_bin[ix][iy]->fill(edge);
}
void fillpT(const unsigned int ix, const unsigned int iy, const double value) {
string edge = "OTHER";
const size_t idx = _axispT.index(value);
if (idx && idx <= _edgespT.size()) {
if (((ix == 0 && iy == 8) || (ix == 1 && (iy == 2 || iy == 4 || iy == 6))) && idx == _edgespT.size())
;
else if (idx == 1 && ((ix == 0 && (iy == 6 || iy == 9)) || (ix == 1 && (iy <= 4 || iy == 7))))
;
else if (idx <= 2 && (ix == 1 && (iy == 5 || iy == 8)))
;
else
edge = _edgespT[idx - 1];
}
_h_dndpt_eta_bin[ix][iy]->fill(edge);
}
void fillpTMax(const unsigned int ix, const double value) {
string edge = "OTHER";
const size_t idx = _axispTMax.index(value);
if (idx && idx <= _edgespTMax.size()) {
if (idx == _edgespTMax.size() && (ix == 1 || ix == 2 || ix == 4))
;
else if (idx >= _edgespTMax.size() - 1 && ix == 5)
;
else
edge = _edgespTMax[idx - 1];
}
_hdndptmax_low_eta_bin[ix]->fill(edge);
}
/// @name Histograms
/// @{
array<CounterPtr, 10> _Nevt_after_cuts;
BinnedHistoPtr<string> _h_dndpt_eta_bin[2][10], _hdndeta_bin[2][10], _hdndptmax_low_eta_bin[8];
CounterPtr _NevAll;
vector<string> _edgespT, _edgesEta, _edgespTMax;
YODA::Axis<double> _axispT = YODA::Axis<double>{0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4,
1.6, 1.825, 2.125, 2.525, 3.125, 4.0, 5.0};
YODA::Axis<double> _axisEta = YODA::Axis<double>{0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.00};
YODA::Axis<double> _axispTMax = YODA::Axis<double>{0.005, 0.255, 0.505, 0.755, 1.005, 1.255, 1.505,
1.755, 2.065, 2.5, 3.125, 4.0, 5.0};
/// @}
};
RIVET_DECLARE_ALIASED_PLUGIN(H1_1996_I424463, H1_1997_I424463);
} ```
Aliases: - H1_1997_I424463