Rivet analyses referenceH1_2013_I1217865Charged particle production in deep-inelastic ep scattering at H1Experiment: H1 (HERA) Inspire ID: 1217865 Status: VALIDATED Authors:
Beam energies: ANY Run details:
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 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Tools/ParticleIdUtils.hh"
4#include "Rivet/Projections/FinalState.hh"
5#include "Rivet/Projections/DISKinematics.hh"
6#include "Rivet/Projections/ChargedFinalState.hh"
7#include "Rivet/Projections/DISKinematics.hh"
8
9namespace Rivet {
10
11
12 /// @brief Charged particle production in deep-inelastic ep scattering at H1
13 class H1_2013_I1217865 : public Analysis {
14 public:
15
16 /// Constructor
17
18 RIVET_DEFAULT_ANALYSIS_CTOR(H1_2013_I1217865);
19
20
21 /// @name Analysis methods
22 /// @{
23
24 /// Book histograms and initialise projections before the run
25 void init() {
26
27 // Initialise and register projections
28 //declare(FinalState(Cuts::abseta < 5 && Cuts::pT > 100*MeV), "FS");
29
30 // Book histograms
31 declare(DISLepton(), "Lepton");
32 declare(DISKinematics(), "Kinematics");
33 declare(ChargedFinalState(), "CFS");
34 declare(FinalState(), "FS");
35 _h_dn_dpT_cen.resize(9);
36 _h_dn_dpT_curr.resize(9);
37 _h_dn_deta_soft.resize(9);
38 _h_dn_deta_hard.resize(9);
39
40 book(_h_dn_dpT_cen[0],19,1,1);
41 book(_h_dn_dpT_curr[0] ,20,1,1);
42 book(_h_dn_deta_soft[0],1,1,1);
43 book(_h_dn_deta_hard[0],2,1,1);
44 for (size_t ix = 0; ix < 9; ++ix) {
45 book(_Nevt_after_cuts[ix], "TMP/Nevt_after_cuts" + to_string(ix));
46 if (ix > 0 ) {
47 book(_h_dn_dpT_cen[ix], ix+20, 1, 1);
48 book(_h_dn_dpT_curr[ix],ix+28,1,1);
49 book(_h_dn_deta_soft[ix],ix+2,1,1);
50 book(_h_dn_deta_hard[ix],ix+10,1,1);
51 }
52 }
53
54
55 }
56
57
58 /// Perform the per-event analysis
59 void analyze(const Event& event) {
60 const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
61 const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
62 const DISLepton& dl = apply<DISLepton>(event,"Lepton");
63
64 // Get the DIS kinematics
65 double x = dk.x();
66 double y = dk.y();
67 double Q2 = dk.Q2()/GeV;
68
69 // Momentum of the scattered lepton
70 FourMomentum leptonMom = dl.out().momentum();
71 double enel = leptonMom.E();
72 double thel = 180.-leptonMom.angle(dl.in().momentum())/degree;
73
74
75 getLog()<<Log::DEBUG<<"enel/GeV = "<<enel/GeV<<", thel = "<<thel<<", y = "<<y<<", x = "<<x<<std::endl;
76 bool cut = y > 0.05 && y < 0.6 && Q2 > 5. && Q2 < 100.;
77 if (!cut) vetoEvent;
78
79
80 int ibin[10] ;
81 for ( int i=0; i< 9; i++) {
82 ibin[i] = 0; }
83
84 ibin[0] = 1 ;
85 if(5.<Q2&&Q2<10.&&x>0.0001&&x<0.00024) ibin[1]=1;
86 if(5.<Q2&&Q2<10.&&x>0.00024&&x<0.0005) ibin[2]=1;
87 if(5.<Q2&&Q2<10.&&x>0.0005&&x<0.002) ibin[3]=1;
88 if(10.<Q2&&Q2<20.&&x>0.0002&&x<0.00052) ibin[4]=1;
89 if(10.<Q2&&Q2<20.&&x>0.00052&&x<0.0011) ibin[5]=1;
90 if(10.<Q2&&Q2<20.&&x>0.0011&&x<0.0037) ibin[6]=1;
91 if(20.<Q2&&Q2<100.&&x>0.0004&&x<0.0017) ibin[7]=1;
92 if(20.<Q2&&Q2<100.&&x>0.0017&&x<0.01) ibin[8]=1;
93
94 for ( int i=0; i< 9; i++) { if(ibin[i]==1) _Nevt_after_cuts[i] ->fill(); }
95
96
97 // Extract the particles other than the lepton
98 Particles particles;
99 particles.reserve(cfs.particles().size());
100 ConstGenParticlePtr dislepGP = dl.out().genParticle();
101 for (const Particle& p : cfs.particles()) {
102 ConstGenParticlePtr loopGP = p.genParticle();
103 if (loopGP == dislepGP) continue;
104 particles.push_back(p);
105 }
106
107 // Boost to hadronic CM
108 const LorentzTransform hcmboost = dk.boostHCM();
109
110 int mult = 0 ;
111 // Loop over the particles
112 // long ncharged(0);
113 for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
114 const Particle& p = particles[ip1];
115
116 double eta = p.momentum().pseudorapidity();
117 double pT = p.momentum().pT()/GeV;
118
119 // Boost to hcm
120 const FourMomentum hcmMom = hcmboost.transform(p.momentum());
121
122 if (pT > 0.15 && eta > -2. && eta < 2.5){
123
124 mult = mult + 1;
125
126 double pThcm =hcmMom.pT() ;
127 double etahcm = hcmMom.pseudorapidity();
128
129
130 if (etahcm > 0. && etahcm < 1.5){
131
132 _h_dn_dpT_cen[0]->fill(pThcm );
133 for ( int i=1; i< 9; i++) {
134 if(ibin[i]==1) { _h_dn_dpT_cen[i]->fill(pThcm );}}
135 }
136
137 if (etahcm > 1.5 && etahcm < 5.){
138 _h_dn_dpT_curr[0]->fill(pThcm );
139 for ( int i=1; i< 9; i++) {
140 if(ibin[i]==1) { _h_dn_dpT_curr[i]->fill(pThcm );}}
141 }
142
143 if(pThcm < 1.){
144 _h_dn_deta_soft[0]->fill(etahcm );
145 for ( int i=1; i< 9; i++) {
146 if(ibin[i]==1) {_h_dn_deta_soft[i]->fill(etahcm );}}
147 }
148
149 if(pThcm > 1. && pThcm < 10.){
150 _h_dn_deta_hard[0]->fill(etahcm );
151 for ( int i=1; i< 9; i++) {
152 if(ibin[i]==1) {_h_dn_deta_hard[i]->fill(etahcm );}}
153
154 }
155 } // if (etahcm > 0. && etahcm < 1.5){
156 } // end of loop over the particles
157
158
159
160 }
161
162
163 /// Normalise histograms etc., after the run
164 void finalize() {
165 // normalize(_h_dn_dpT_cen);
166
167 if (_Nevt_after_cuts[0]->val() != 0) scale(_h_dn_dpT_cen[0], 1./ *_Nevt_after_cuts[0]);
168 if (_Nevt_after_cuts[0]->val() != 0) scale(_h_dn_dpT_curr[0], 1./ *_Nevt_after_cuts[0]);
169 if (_Nevt_after_cuts[0]->val() != 0) scale(_h_dn_deta_soft[0], 1./ *_Nevt_after_cuts[0]);
170 if (_Nevt_after_cuts[0]->val() != 0) scale(_h_dn_deta_hard[0], 1./ *_Nevt_after_cuts[0]);
171
172
173 for ( int i=1; i< 9; i++) {
174 if (_Nevt_after_cuts[i]->val() != 0) {
175 scale(_h_dn_dpT_cen[i], 1./ *_Nevt_after_cuts[i]);
176 scale(_h_dn_dpT_curr[i], 1./ *_Nevt_after_cuts[i]);
177 scale(_h_dn_deta_soft[i],1./ *_Nevt_after_cuts[i]);
178 scale(_h_dn_deta_hard[i],1./ *_Nevt_after_cuts[i]);
179 }
180 }
181
182
183
184 }
185 private:
186
187 /**
188 * Polar angle with right direction of the beam
189 */
190 inline double beamAngle(const FourVector& v, const bool & order) {
191 double thel = v.polarAngle()/degree;
192 if(thel<0.) thel+=180.;
193 if(!order) thel = 180.-thel;
194 return thel;
195 }
196
197 /// @}
198
199
200 /// @name Histograms
201 /// @{
202 Histo1DPtr _h_dn_dpT_2r;
203 Histo1DPtr _h_dn_dpT_2l;
204
205 vector<Histo1DPtr> _h_dn_dpT_cen;
206 vector<Histo1DPtr> _h_dn_dpT_curr;
207 vector<Histo1DPtr> _h_dn_deta_soft;
208 vector<Histo1DPtr> _h_dn_deta_hard;
209 array<CounterPtr,9> _Nevt_after_cuts;
210
211
212 /// @}
213
214
215 };
216
217
218 RIVET_DECLARE_PLUGIN(H1_2013_I1217865);
219
220
221}
|