Rivet analyses referenceH1_1996_I424463Transverse momentum spectra of charged particles in DIS (H1 1996)Experiment: H1 (HERA) Inspire ID: 424463 Status: VALIDATED Authors:
Beam energies: (27.5, 820.0); (820.0, 27.5) GeV
Transverse momentum spectra of charged particles produced in deep inelastic scattering are measured as a function of the kinematic variables x and Q2 using the H1 detector at the epcollider HERA. The data are compared to different patton emission models, either with or without ordering of the emissions in transverse momentum. The data provide evidence for a relatively large amount of patton radiation between the current and the remnant systems. Source code: H1_1996_I424463.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 Transverse momentum spectra of charged particles in DIS
13 class H1_1996_I424463 : public Analysis {
14 public:
15
16 /// Constructor
17 RIVET_DEFAULT_ANALYSIS_CTOR(H1_1996_I424463);
18
19
20 /// @name Analysis methods
21 /// @{
22
23 /// Book histograms and initialise projections before the run
24 void init() {
25
26 // Book projections
27 declare(DISLepton(), "Lepton");
28 declare(DISKinematics(), "Kinematics");
29 declare(ChargedFinalState(), "CFS");
30 declare(FinalState(), "FS");
31
32 // Book histograms
33 book(_NevAll, "TMP/Nev_all");
34 int ixx = 0 ;
35 for (size_t ix = 0; ix < 10; ++ix) {
36 book(_Nevt_after_cuts[ix], "TMP/Nevt_after_cuts" + to_string(ix));
37 for(unsigned int ih=0;ih<2;++ih) {
38 book(_h_dndpt_eta_bin[ih][ix], ih*10 + ix+1 , 1, 1);
39 book(_hdndeta_bin [ih][ix], ih*10 + ix+29, 1, 1);
40 }
41 if (ix != 6 && ix != 9) {
42 book(_hdndptmax_low_eta_bin[ixx], ixx+21, 1, 1);
43 ixx=ixx+1;
44 }
45 }
46 }
47
48
49 /// Perform the per-event analysis
50 void analyze(const Event& event) {
51 if (_edgespT.empty()) {
52 _edgespT = _h_dndpt_eta_bin [0][0]->xEdges();
53 _edgesEta = _hdndeta_bin [0][0]->xEdges();
54 _edgespTMax = _hdndptmax_low_eta_bin[0]->xEdges();
55 }
56 const FinalState& fs = apply<FinalState>(event, "FS");
57 const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
58 const DISLepton& dl = apply<DISLepton>(event,"Lepton");
59
60 // Get the DIS kinematics
61 double x = dk.x();
62 double y = dk.y();
63 double Q2 = dk.Q2()/GeV;
64 double W2 = dk.W2()/GeV;
65
66 // Momentum of the scattered lepton
67 FourMomentum leptonMom = dl.out().momentum();
68 double enel = leptonMom.E()/GeV;
69 double thel = 180.-leptonMom.angle(dl.in().momentum())/degree;
70
71 _NevAll -> fill() ;
72
73 bool cut = y > 0.05 && Q2 > 5. && Q2 < 100.&& enel > 12. && W2 > 4400. && thel > 157. && thel < 173.;
74 if (!cut) vetoEvent;
75
76 int ibin[10] ;
77 for (int i = 0; i < 10; i++) ibin[i] = 0;
78 if ( 5. < Q2 && Q2 < 50. && x > 0.0001 && x < 0.0010) ibin[0] = 1;
79 if ( 5. < Q2 && Q2 < 10. && x > 0.0001 && x < 0.0002) ibin[1] = 1;
80 if ( 6. < Q2 && Q2 < 10. && x > 0.0002 && x < 0.0005) ibin[2] = 1;
81 if (10. < Q2 && Q2 < 20. && x > 0.0002 && x < 0.0005) ibin[3] = 1;
82 if (10. < Q2 && Q2 < 20. && x > 0.0005 && x < 0.0008) ibin[4] = 1;
83 if (10. < Q2 && Q2 < 20. && x > 0.0008 && x < 0.0015) ibin[5] = 1;
84 if (10. < Q2 && Q2 < 20. && x > 0.0015 && x < 0.0040) ibin[6] = 1;
85 if (20. < Q2 && Q2 < 50. && x > 0.0005 && x < 0.0014) ibin[7] = 1;
86 if (20. < Q2 && Q2 < 50. && x > 0.0014 && x < 0.0030) ibin[8] = 1;
87 if (20. < Q2 && Q2 < 50. && x > 0.0030 && x < 0.0100) ibin[9] = 1;
88 for (int i = 0; i < 10; i++) {
89 if (ibin[i]==1) _Nevt_after_cuts[i]->fill();
90 }
91
92 // Extract the particles other than the lepton
93 /// @todo Improve to avoid HepMC digging
94 Particles particles;
95 particles.reserve(fs.particles().size());
96 ConstGenParticlePtr dislepGP = dl.out().genParticle();
97 for (const Particle& p : fs.particles()) {
98 ConstGenParticlePtr loopGP = p.genParticle();
99 if (loopGP == dislepGP) continue;
100 particles.push_back(p);
101 }
102
103 // Boost to hadronic CM
104 const LorentzTransform hcmboost = dk.boostHCM();
105
106 // Loop over the particles
107 int mult = 0 ;
108 double ptmax_high[10], ptmax_low[10];
109 for (int i = 0; i < 10; i++) {
110 ptmax_high[i] = 0.; ptmax_low[i] = 0.;
111 }
112 double EtSum = 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 // Boost to hcm
118 const FourMomentum hcmMom = hcmboost.transform(p.momentum());
119
120 // Apply safety cuts
121 if (eta > -5 && eta < 10.) {
122 mult = mult + 1;
123 double pThcm =hcmMom.pT() ;
124 double etahcm = hcmMom.pseudorapidity();
125 if (etahcm > 0. && etahcm < 2.0) {
126 EtSum = EtSum + hcmMom.Et();
127 }
128 if (PID::charge(p.pid()) != 0) {
129 if (etahcm > 0.5 && etahcm < 1.5) {
130 for ( int i=0; i< 10; i++) {
131 if (ibin[i]==1) {
132 fillpT(1,i,pThcm);
133 if (pThcm > ptmax_low[i] ) ptmax_low[i] = pThcm;
134 }
135 }
136
137 }
138 if (etahcm > 1.5 && etahcm < 2.5){
139 for ( int i=0; i< 10; i++) {
140 if (ibin[i]==1) {
141 fillpT(0,i,pThcm);
142 if (pThcm > ptmax_high[i] ) ptmax_high[i] = pThcm;
143 }
144 }
145 }
146 for ( int i=0; i< 10; i++) {
147 if (ibin[i]==1) fillEta(1,i,etahcm);
148 if (ibin[i]==1 && pThcm > 1.) fillEta(0,i,etahcm);
149 }
150 }
151 } // end of loop over the particles
152 }
153 int ii=0;
154 for ( int i=0; i< 10; i++) {
155 if (i != 6 && i != 9 ) {
156 if ( ibin[i]==1 && EtSum > 6. ) {
157 fillpTMax(ii,ptmax_low[i]);
158 }
159 ii=ii+1;
160 }
161 }
162 }
163
164
165 /// Normalise histograms etc., after the run
166 void finalize() {
167 MSG_DEBUG("All events: " << _NevAll->val() << " after cuts: "<< _Nevt_after_cuts[0]->val());
168 MSG_DEBUG("Cut1 events: " << _NevAll->val() << " after cuts: "<< _Nevt_after_cuts[1]->val());
169 int ii = 0;
170 for ( int i=0; i< 10; i++) {
171 if (_Nevt_after_cuts[i]->val() != 0) {
172 for(unsigned int ih=0;ih<2;++ih) {
173 scale(_h_dndpt_eta_bin[ih][i], 1./ *_Nevt_after_cuts[i]);
174 size_t ioff=0;
175 if( (ih==0 && (i==6 || i==9)) ||
176 (ih==1 && (i<=4 || i==7))) ioff=1;
177 else if (ih==1 && (i==5 || i==8)) ioff=2;
178 for(auto & b : _h_dndpt_eta_bin[ih][i]->bins()) {
179 const size_t idx = b.index()+ioff;
180 b.scaleW(1./_axispT.width(idx));
181 }
182 scale(_hdndeta_bin [ih][i], 1./ *_Nevt_after_cuts[i]);
183 ioff=0;
184 if( (ih==0 && (i==5 || i==6 || i==8 || i==9)) ||
185 (ih==1 && (i==0 || i==4 || i==5 || i==8))) ioff=1;
186 else if (ih==1 && (i==6 || i==9)) ioff=2;
187 for(auto & b : _hdndeta_bin[ih][i]->bins()) {
188 const size_t idx = b.index()+ioff;
189 b.scaleW(1./_axisEta.width(idx));
190 }
191 }
192 }
193 if ( i !=6 && i!=9) {
194 if (_Nevt_after_cuts[i]->val() != 0) {
195 normalize(_hdndptmax_low_eta_bin[ii]);
196 for(auto & b : _hdndptmax_low_eta_bin[ii]->bins()) {
197 const size_t idx = b.index();
198 b.scaleW(1./_axispTMax.width(idx));
199 }
200 }
201 ii=ii+1;
202 }
203 }
204
205 }
206 /// @}
207
208 void fillEta(const unsigned int ix, const unsigned int iy, const double value) {
209 string edge = "OTHER";
210 const size_t idx = _axisEta.index(value);
211 if (idx && idx <= _edgesEta.size()) {
212 if ( ( (ix==0 && iy>=1&& iy<=3) || (ix==0 && (iy==1 || iy==3)) ) && idx == _edgesEta.size())
213 ;
214 else if (idx==1 && ( (ix==0 && (iy==5 || iy==6 || iy==8 || iy==9)) ||
215 (ix==1 && (iy==0 || iy==4 || iy==5 || iy==8))))
216 ;
217 else if (idx<=2 && (ix==1 && (iy==6 || iy==9)) )
218 ;
219 else {
220 edge = _edgesEta[idx-1];
221 }
222 }
223 _hdndeta_bin[ix][iy]->fill(edge);
224 }
225
226 void fillpT(const unsigned int ix, const unsigned int iy, const double value) {
227 string edge = "OTHER";
228 const size_t idx = _axispT.index(value);
229 if (idx && idx <= _edgespT.size()) {
230 if(((ix==0 && iy==8) || (ix==1 && (iy==2 || iy==4 || iy==6))) && idx == _edgespT.size())
231 ;
232 else if (idx==1 && ( (ix==0 && (iy==6 || iy==9)) ||
233 (ix==1 && (iy<=4 || iy==7))))
234 ;
235 else if (idx<=2 && (ix==1 && (iy==5 || iy==8)))
236 ;
237 else
238 edge = _edgespT[idx-1];
239 }
240 _h_dndpt_eta_bin[ix][iy]->fill(edge);
241 }
242
243 void fillpTMax(const unsigned int ix, const double value) {
244 string edge = "OTHER";
245 const size_t idx = _axispTMax.index(value);
246 if (idx && idx <= _edgespTMax.size()) {
247 if(idx==_edgespTMax.size() && (ix==1 || ix==2 || ix==4))
248 ;
249 else if(idx>=_edgespTMax.size()-1 && ix==5)
250 ;
251 else
252 edge = _edgespTMax[idx-1];
253 }
254 _hdndptmax_low_eta_bin[ix]->fill(edge);
255 }
256
257 /// @name Histograms
258 /// @{
259 array<CounterPtr,10> _Nevt_after_cuts;
260 BinnedHistoPtr<string> _h_dndpt_eta_bin[2][10], _hdndeta_bin[2][10], _hdndptmax_low_eta_bin[8];
261 CounterPtr _NevAll ;
262
263 vector<string> _edgespT,_edgesEta,_edgespTMax;
264 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};
265 YODA::Axis<double> _axisEta = YODA::Axis<double>{ 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.00};
266 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};
267 /// @}
268
269 };
270
271
272 RIVET_DECLARE_ALIASED_PLUGIN(H1_1996_I424463, H1_1997_I424463);
273
274}
|