Rivet analyses
title: H1_2013_I1217865
Charged particle production in deep-inelastic ep scattering at H1
Experiment: H1 (HERA)
Inspire ID: 1217865
Status: VALIDATED
Authors: - Anastasia Grebenyuk - Hannes Jung
References: - Eur.Phys.J. C73 (2013) 2406, - arXiv: 1302.1321
Beams: e+ p+, p+ e+, e- p+, p+ e-
Beam energies: ANY
Run details: - Inclusive DIS
Charged particle production in deep-inelastic ep scattering is measured with the H1 detector at HERA. The kinematic range of the analysis covers low photon virtualities, 5 < Q2 < 100 GeV2 and small values of Bjorken-x, 10-2 < x < 10-2. The analysis is performed in the hadronic centre-of-mass system. The charged particle densities are measured as a function of pseudorapidity eta and transverse momentum pT in the range 0 <eta <5 and 0 < pT < 10 GeV in bins of x and Q2.
Source code:H1_2013_I1217865.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 Charged particle production in deep-inelastic ep scattering at H1 class H1_2013_I1217865 : public Analysis { public:
/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(H1_2013_I1217865);
/// @name Analysis methods
/// @{
/// Book histograms and initialise projections before the run
void init() {
// Initialise and register projections
//declare(FinalState(Cuts::abseta < 5 && Cuts::pT > 100*MeV), "FS");
// Book histograms
declare(DISLepton(), "Lepton");
declare(DISKinematics(), "Kinematics");
declare(ChargedFinalState(), "CFS");
declare(FinalState(), "FS");
_h_dn_dpT_cen.resize(9);
_h_dn_dpT_curr.resize(9);
_h_dn_deta_soft.resize(9);
_h_dn_deta_hard.resize(9);
book(_h_dn_dpT_cen[0], 19, 1, 1);
book(_h_dn_dpT_curr[0], 20, 1, 1);
book(_h_dn_deta_soft[0], 1, 1, 1);
book(_h_dn_deta_hard[0], 2, 1, 1);
for (size_t ix = 0; ix < 9; ++ix) {
book(_Nevt_after_cuts[ix], "TMP/Nevt_after_cuts" + to_string(ix));
if (ix > 0) {
book(_h_dn_dpT_cen[ix], ix + 20, 1, 1);
book(_h_dn_dpT_curr[ix], ix + 28, 1, 1);
book(_h_dn_deta_soft[ix], ix + 2, 1, 1);
book(_h_dn_deta_hard[ix], ix + 10, 1, 1);
}
}
}
/// Perform the per-event analysis
void analyze(const Event& event) {
const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
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;
// Momentum of the scattered lepton
FourMomentum leptonMom = dl.out().momentum();
double enel = leptonMom.E();
double thel = 180. - leptonMom.angle(dl.in().momentum()) / degree;
getLog() << Log::DEBUG << "enel/GeV = " << enel / GeV << ", thel = " << thel << ", y = " << y
<< ", x = " << x << std::endl;
bool cut = y > 0.05 && y < 0.6 && Q2 > 5. && Q2 < 100.;
if (!cut) vetoEvent;
int ibin[10];
for (int i = 0; i < 9; i++) {
ibin[i] = 0;
}
ibin[0] = 1;
if (5. < Q2 && Q2 < 10. && x > 0.0001 && x < 0.00024) ibin[1] = 1;
if (5. < Q2 && Q2 < 10. && x > 0.00024 && x < 0.0005) ibin[2] = 1;
if (5. < Q2 && Q2 < 10. && x > 0.0005 && x < 0.002) ibin[3] = 1;
if (10. < Q2 && Q2 < 20. && x > 0.0002 && x < 0.00052) ibin[4] = 1;
if (10. < Q2 && Q2 < 20. && x > 0.00052 && x < 0.0011) ibin[5] = 1;
if (10. < Q2 && Q2 < 20. && x > 0.0011 && x < 0.0037) ibin[6] = 1;
if (20. < Q2 && Q2 < 100. && x > 0.0004 && x < 0.0017) ibin[7] = 1;
if (20. < Q2 && Q2 < 100. && x > 0.0017 && x < 0.01) ibin[8] = 1;
for (int i = 0; i < 9; i++) {
if (ibin[i] == 1) _Nevt_after_cuts[i]->fill();
}
// Extract the particles other than the lepton
Particles particles;
particles.reserve(cfs.particles().size());
ConstGenParticlePtr dislepGP = dl.out().genParticle();
for (const Particle& p : cfs.particles()) {
ConstGenParticlePtr loopGP = p.genParticle();
if (loopGP == dislepGP) continue;
particles.push_back(p);
}
// Boost to hadronic CM
const LorentzTransform hcmboost = dk.boostHCM();
int mult = 0;
// Loop over the particles
// long ncharged(0);
for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
const Particle& p = particles[ip1];
double eta = p.momentum().pseudorapidity();
double pT = p.momentum().pT() / GeV;
// Boost to hcm
const FourMomentum hcmMom = hcmboost.transform(p.momentum());
if (pT > 0.15 && eta > -2. && eta < 2.5) {
mult = mult + 1;
double pThcm = hcmMom.pT();
double etahcm = hcmMom.pseudorapidity();
if (etahcm > 0. && etahcm < 1.5) {
_h_dn_dpT_cen[0]->fill(pThcm);
for (int i = 1; i < 9; i++) {
if (ibin[i] == 1) {
_h_dn_dpT_cen[i]->fill(pThcm);
}
}
}
if (etahcm > 1.5 && etahcm < 5.) {
_h_dn_dpT_curr[0]->fill(pThcm);
for (int i = 1; i < 9; i++) {
if (ibin[i] == 1) {
_h_dn_dpT_curr[i]->fill(pThcm);
}
}
}
if (pThcm < 1.) {
_h_dn_deta_soft[0]->fill(etahcm);
for (int i = 1; i < 9; i++) {
if (ibin[i] == 1) {
_h_dn_deta_soft[i]->fill(etahcm);
}
}
}
if (pThcm > 1. && pThcm < 10.) {
_h_dn_deta_hard[0]->fill(etahcm);
for (int i = 1; i < 9; i++) {
if (ibin[i] == 1) {
_h_dn_deta_hard[i]->fill(etahcm);
}
}
}
} // if (etahcm > 0. && etahcm < 1.5){
} // end of loop over the particles
}
/// Normalise histograms etc., after the run
void finalize() {
// normalize(_h_dn_dpT_cen);
if (_Nevt_after_cuts[0]->val() != 0) scale(_h_dn_dpT_cen[0], 1. / *_Nevt_after_cuts[0]);
if (_Nevt_after_cuts[0]->val() != 0) scale(_h_dn_dpT_curr[0], 1. / *_Nevt_after_cuts[0]);
if (_Nevt_after_cuts[0]->val() != 0) scale(_h_dn_deta_soft[0], 1. / *_Nevt_after_cuts[0]);
if (_Nevt_after_cuts[0]->val() != 0) scale(_h_dn_deta_hard[0], 1. / *_Nevt_after_cuts[0]);
for (int i = 1; i < 9; i++) {
if (_Nevt_after_cuts[i]->val() != 0) {
scale(_h_dn_dpT_cen[i], 1. / *_Nevt_after_cuts[i]);
scale(_h_dn_dpT_curr[i], 1. / *_Nevt_after_cuts[i]);
scale(_h_dn_deta_soft[i], 1. / *_Nevt_after_cuts[i]);
scale(_h_dn_deta_hard[i], 1. / *_Nevt_after_cuts[i]);
}
}
}
private:
/**
* Polar angle with right direction of the beam
*/
inline double beamAngle(const FourVector& v, const bool& order) {
double thel = v.polarAngle() / degree;
if (thel < 0.) thel += 180.;
if (!order) thel = 180. - thel;
return thel;
}
/// @}
/// @name Histograms
/// @{
Histo1DPtr _h_dn_dpT_2r;
Histo1DPtr _h_dn_dpT_2l;
vector<Histo1DPtr> _h_dn_dpT_cen;
vector<Histo1DPtr> _h_dn_dpT_curr;
vector<Histo1DPtr> _h_dn_deta_soft;
vector<Histo1DPtr> _h_dn_deta_hard;
array<CounterPtr, 9> _Nevt_after_cuts;
/// @}
};
RIVET_DECLARE_PLUGIN(H1_2013_I1217865);
} ```