Rivet analyses
H1 energy flow in DIS
Experiment: H1 (HERA)
Inspire ID: 503947
Status: VALIDATED
Authors: - Peter Richardson
References: - Eur.Phys.J.C12:595-607,2000 - DOI: 10.1007/s100520000287 - arXiv: hep-ex/9907027
Beams: p+ e+
Beam energies: (820.0, 27.5)GeV
Run details: - e+p deep inelastic scattering with p at 820~GeV, e+ at 27.5 GeV $\sqrt{s} = 300~\GeV$
Measurements of transverse energy flow for neutral current deep- inelastic scattering events produced in positron-proton collisions at HERA. The kinematic range covers squared momentum transfers Q2 from 3.2 to 2200 GeV2; the Bjorken scaling variable x from 8 × 10−5 to 0.11 and the hadronic mass W from 66 to 233 GeV. The transverse energy flow is measured in the hadronic centre of mass frame and is studied as a function of Q2, x, W and pseudorapidity. The behaviour of the mean transverse energy in the central pseudorapidity region and an interval corresponding to the photon fragmentation region are analysed as a function of Q2 and W. This analysis is useful for exploring the effect of photon PDFs and for tuning models of parton evolution and treatment of fragmentation and the proton remnant in DIS.
Source
code:H1_2000_I503947.cc
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Math/Constants.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/DISKinematics.hh"
namespace Rivet {
/// @brief H1 energy flow and charged particle spectra
///
/// @author Peter Richardson
///
/// Based on the HZTOOL analysis HZ99091
class H1_2000_I503947 : public Analysis {
public:
/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(H1_2000_I503947);
/// @name Analysis methods
/// @{
/// Initialise projections and histograms
void init() {
// Projections
const DISLepton dl;
declare(dl, "Lepton");
declare(DISKinematics(), "Kinematics");
declare(dl.remainingFinalState(), "FS");
// Histograms and weight vectors for low Q^2 a
_histETLowQa.resize(17);
for (size_t ix = 0; ix < 17; ++ix) {
book(_histETLowQa[ix], ix+1, 1, 1);
book(_weightETLowQa[ix], "TMP/ETLowQa" + to_string(ix));
}
// Histograms and weight vectors for high Q^2 a
_histETHighQa.resize(7);
for (size_t ix = 0; ix < 7; ++ix) {
book(_histETHighQa[ix], ix+18, 1, 1);
book(_weightETHighQa[ix], "TMP/ETHighQa" + to_string(ix));
}
// Histograms and weight vectors for low Q^2 b
_histETLowQb.resize(5);
for (size_t ix = 0; ix < 5; ++ix) {
book(_histETLowQb[ix], ix+25, 1, 1);
book(_weightETLowQb[ix], "TMP/ETLowQb" + to_string(ix));
}
// Histograms and weight vectors for high Q^2 b
_histETHighQb.resize(5);
for (size_t ix = 0; ix < 3; ++ix) {
book(_histETHighQb[ix], 30+ix, 1, 1);
book(_weightETHighQb[ix], "TMP/ETHighQb" + to_string(ix));
}
// Histograms for the averages
book(_histAverETCentral, 33, 1, 1);
book(_histAverETFrag, 34, 1, 1);
}
/// Analyze each event
void analyze(const Event& event) {
// DIS kinematics
const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
if ( dk.failed() ) vetoEvent;
double q2 = dk.Q2();
double x = dk.x();
double y = dk.y();
double w2 = dk.W2();
// Kinematics of the scattered lepton
const DISLepton& dl = apply<DISLepton>(event,"Lepton");
if ( dl.failed() ) vetoEvent;
const FourMomentum leptonMom = dl.out();
const double enel = leptonMom.E();
const double thel = 180 - leptonMom.angle(dl.in().mom())/degree;
// Extract the particles other than the lepton
Particles particles= apply<FinalState>(event, "FS").particles();
// Cut on the forward energy
double efwd = 0.;
for (const Particle& p : particles) {
const double th = 180 - p.angle(dl.in())/degree;
if (inRange(th, 4.4, 15.0)) efwd += p.E();
}
// There are four possible selections for events
bool evcut[4];
// Low Q2 selection a
evcut[0] = enel/GeV > 12. && w2 >= 4400.*GeV2 && efwd/GeV > 0.5 && inRange(thel,157.,176.);
// Low Q2 selection b
evcut[1] = enel/GeV > 12. && inRange(y,0.3,0.5);
// High Q2 selection a
evcut[2] = inRange(thel,12.,150.) && inRange(y,0.05,0.6) && w2 >= 4400.*GeV2 && efwd > 0.5;
// High Q2 selection b
evcut[3] = inRange(thel,12.,150.) && inRange(y,0.05,0.6) && inRange(w2,27110.*GeV2,45182.*GeV2);
// Veto if fails all cuts
/// @todo Can we use all()?
if (! (evcut[0] || evcut[1] || evcut[2] || evcut[3]) ) vetoEvent;
// Find the bins
int bin[4] = {-1,-1,-1,-1};
// For the low Q2 selection a)
if (q2 > 2.5*GeV2 && q2 <= 5.*GeV2) {
if (x > 0.00005 && x <= 0.0001 ) bin[0] = 0;
if (x > 0.0001 && x <= 0.0002 ) bin[0] = 1;
if (x > 0.0002 && x <= 0.00035) bin[0] = 2;
if (x > 0.00035 && x <= 0.0010 ) bin[0] = 3;
}
else if (q2 > 5.*GeV2 && q2 <= 10.*GeV2) {
if (x > 0.0001 && x <= 0.0002 ) bin[0] = 4;
if (x > 0.0002 && x <= 0.00035) bin[0] = 5;
if (x > 0.00035 && x <= 0.0007 ) bin[0] = 6;
if (x > 0.0007 && x <= 0.0020 ) bin[0] = 7;
}
else if (q2 > 10.*GeV2 && q2 <= 20.*GeV2) {
if (x > 0.0002 && x <= 0.0005) bin[0] = 8;
if (x > 0.0005 && x <= 0.0008) bin[0] = 9;
if (x > 0.0008 && x <= 0.0015) bin[0] = 10;
if (x > 0.0015 && x <= 0.040 ) bin[0] = 11;
}
else if (q2 > 20.*GeV2 && q2 <= 50.*GeV2) {
if (x > 0.0005 && x <= 0.0014) bin[0] = 12;
if (x > 0.0014 && x <= 0.0030) bin[0] = 13;
if (x > 0.0030 && x <= 0.0100) bin[0] = 14;
}
else if (q2 > 50.*GeV2 && q2 <= 100.*GeV2) {
if (x >0.0008 && x <= 0.0030) bin[0] = 15;
if (x >0.0030 && x <= 0.0200) bin[0] = 16;
}
// check in one of the bins
evcut[0] &= bin[0] >= 0;
// For the low Q2 selection b)
if (q2 > 2.5*GeV2 && q2 <= 5. *GeV2) bin[1] = 0;
if (q2 > 5. *GeV2 && q2 <= 10. *GeV2) bin[1] = 1;
if (q2 > 10.*GeV2 && q2 <= 20. *GeV2) bin[1] = 2;
if (q2 > 20.*GeV2 && q2 <= 50. *GeV2) bin[1] = 3;
if (q2 > 50.*GeV2 && q2 <= 100.*GeV2) bin[1] = 4;
// check in one of the bins
evcut[1] &= bin[1] >= 0;
// for the high Q2 selection a)
if (q2 > 100.*GeV2 && q2 <= 400.*GeV2) {
if (x > 0.00251 && x <= 0.00631) bin[2] = 0;
if (x > 0.00631 && x <= 0.0158 ) bin[2] = 1;
if (x > 0.0158 && x <= 0.0398 ) bin[2] = 2;
}
else if (q2 > 400.*GeV2 && q2 <= 1100.*GeV2) {
if (x > 0.00631 && x <= 0.0158 ) bin[2] = 3;
if (x > 0.0158 && x <= 0.0398 ) bin[2] = 4;
if (x > 0.0398 && x <= 1. ) bin[2] = 5;
}
else if (q2 > 1100.*GeV2 && q2 <= 100000.*GeV2) {
if (x > 0. && x <= 1.) bin[2] = 6;
}
// check in one of the bins
evcut[2] &= bin[2] >= 0;
// for the high Q2 selection b)
if (q2 > 100.*GeV2 && q2 <= 220.*GeV2) bin[3] = 0;
else if (q2 > 220.*GeV2 && q2 <= 400.*GeV2) bin[3] = 1;
else if (q2 > 400. ) bin[3] = 2;
// check in one of*GeV the bins
evcut[3] &= bin[3] >= 0;
// Veto if fails all cuts after bin selection
/// @todo Can we use all()?
if (! (evcut[0] || evcut[1] || evcut[2] || evcut[3])) vetoEvent;
// Increment the count for normalisation
if (evcut[0]) _weightETLowQa [bin[0]]->fill();
if (evcut[1]) _weightETLowQb [bin[1]]->fill();
if (evcut[2]) _weightETHighQa[bin[2]]->fill();
if (evcut[3]) _weightETHighQb[bin[3]]->fill();
// Boost to hadronic CoM
const LorentzTransform hcmboost = dk.boostHCM();
// Loop over the particles
double etcent = 0;
double etfrag = 0;
for (const Particle& p : particles) {
// Boost momentum to CMS
const FourMomentum hcmMom = hcmboost.transform(p.momentum());
double et = fabs(hcmMom.Et());
double eta = hcmMom.eta();
// Averages in central and forward region
if (fabs(eta) < 0.5 ) etcent += et;
if (eta > 2 && eta <= 3.) etfrag += et;
// Histograms of Et flow
if (evcut[0]) _histETLowQa [bin[0]]->fill(eta, et);
if (evcut[1]) _histETLowQb [bin[1]]->fill(eta, et);
if (evcut[2]) _histETHighQa[bin[2]]->fill(eta, et);
if (evcut[3]) _histETHighQb[bin[3]]->fill(eta, et);
}
// Fill histograms for the average quantities
if (evcut[1] || evcut[3]) {
_histAverETCentral->fill(q2, etcent);
_histAverETFrag ->fill(q2, etfrag);
}
}
// Finalize
void finalize() {
// Normalization of the Et distributions to unit cross-section
// These histograms are filled once per particles though,
// so no unit area to be expected
/// @todo Simplify by using normalize() instead? Are all these being normalized to area=1?
for (size_t ix = 0; ix < _weightETLowQa.size(); ++ix) {
if (_weightETLowQa[ix]->val() ) scale(_histETLowQa[ix], 1/ *_weightETLowQa[ix]);
}
for (size_t ix = 0; ix < _weightETHighQa.size(); ++ix) {
if (_weightETHighQa[ix]->val()) scale(_histETHighQa[ix], 1/ *_weightETHighQa[ix]);
}
for (size_t ix = 0; ix < _weightETLowQb.size(); ++ix) {
if (_weightETLowQb[ix]->val() ) scale(_histETLowQb[ix], 1/ *_weightETLowQb[ix]);
}
for (size_t ix = 0; ix < _weightETHighQb.size(); ++ix) {
if (_weightETHighQb[ix]->val()) scale(_histETHighQb[ix], 1/ *_weightETHighQb[ix]);
}
}
/// @}
private:
/// @name Histograms
/// @{
vector<Histo1DPtr> _histETLowQa;
vector<Histo1DPtr> _histETHighQa;
vector<Histo1DPtr> _histETLowQb;
vector<Histo1DPtr> _histETHighQb;
Profile1DPtr _histAverETCentral;
Profile1DPtr _histAverETFrag;
/// @}
/// @name Storage of weights for normalisation
/// @{
array<CounterPtr,17> _weightETLowQa;
array<CounterPtr, 7> _weightETHighQa;
array<CounterPtr, 5> _weightETLowQb;
array<CounterPtr, 3> _weightETHighQb;
/// @}
};
RIVET_DECLARE_ALIASED_PLUGIN(H1_2000_I503947, H1_2000_S4129130);
}Aliases: - H1_2000_S4129130