Rivet analyses referenceH1_1995_I394793A Study of the Fragmentation of Quarks in ep Collisions at HERAExperiment: H1 (HERA) Inspire ID: 394793 Status: VALIDATED HEPDATA Authors:
Beam energies: (26.7, 820.0); (820.0, 26.7) GeV
Deep inelastic scattering (DIS) events, selected from 1993 data taken by the H1 experiment at HERA, are studied in the Breit frame of reference. It is shown that certain aspects of the quarks emerging from within the proton in $ep$ interactions are essentially the same as those of quarks pair-created from the vacuum in $e^+e^-$ annihilation. Source code: H1_1995_I394793.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/PromptFinalState.hh"
4#include "Rivet/Projections/DISKinematics.hh"
5#include "Rivet/Projections/ChargedFinalState.hh"
6#include "Rivet/Projections/DISLepton.hh"
7
8namespace Rivet {
9
10
11 /// @brief A Study of the Fragmentation of Quarks in ep Collisions at HERA (H1)
12 class H1_1995_I394793 : public Analysis {
13 public:
14
15 /// Constructor
16 RIVET_DEFAULT_ANALYSIS_CTOR(H1_1995_I394793);
17
18
19 /// @name Analysis methods
20 ///@{
21
22 /// Book histograms and initialise projections before the run
23 void init() {
24
25 // Initialise and register projections
26 declare(DISLepton(), "Lepton");
27 declare(DISKinematics(), "Kinematics");
28
29 // The basic final-state projection:
30 // all final-state particles within
31 // the given eta acceptance
32 const FinalState fs(Cuts::abseta < 4.9);
33 declare(fs, "FS");
34 const ChargedFinalState cfs;
35 declare(cfs, "CFS");
36
37 // book a counter
38 book(_Nevt_after_cuts, "TMP/Nevt_after_cuts");
39 book(_Nevt_afterfwd_cuts, "TMP/Nevt_afterfwd_cuts");
40 book(_Nevt_afterh_cuts, "TMP/Nevt_afterh_cuts");
41 book(_Nevt_afterhfwd_cuts, "TMP/Nevt_afterhfwd_cuts");
42
43
44 // Book histograms
45 // specify custom binning
46 // take binning from reference data using HEPData ID (digits in "d01-x01-y01" etc.)
47 book(_h["costh_lowQ"], 1, 1, 1);
48 book(_h["costh_highQ"], 1, 1, 2);
49 book(_h["costh_lowQ_noEfwd"], 2, 1, 1);
50 book(_h["costh_highQ_noEfwd"], 2, 1, 2);
51 book(_h["xp_posCharge_lowQ"], 3, 1, 1);
52 book(_h["xp_negCharge_lowQ"], 3, 1, 2);
53 book(_h["xp_posCharge_highQ"], 3, 1, 3);
54 book(_h["xp_negCharge_highQ"], 3, 1, 4);
55 book(_h["ksi_lowQ"], 4, 1, 1);
56 book(_h["ksi_highQ"], 4, 1, 2);
57 book(_s["Mult_vrs_Q2"], 5, 1, 1);
58 book(_s["Mult_vrs_Q2_noEfwd"], 6, 1, 1);
59
60 book(_h["Mult_vrs_Q2_nchrg"],"TMP/Mult_vrs_Q2_nchrg", refData(5,1,1));
61 book(_h["Mult_vrs_Q2_noEfwd_nchrg"],"TMP/Mult_vrs_Q2_noEfwd_nchrg", refData(6,1,1));
62 book(_h["Mult_vrs_Q2_count"],"TMP/Mult_vrs_Q2_count", refData(5,1,1));
63 book(_h["Mult_vrs_Q2_noEfwd_count"],"TMP/Mult_vrs_Q2_noEfwd_count", refData(6,1,1));
64
65 }
66
67
68 /// Perform the per-event analysis
69 void analyze(const Event& event) {
70
71 const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
72
73
74 //DIS kinematics
75 const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
76 if ( dk.failed() ) vetoEvent;
77 double y = dk.y();
78 double w2 = dk.W2();
79 double Q2 = dk.Q2();
80
81
82
83 bool cut ;
84
85 cut = Q2 > 12 && y < 0.6 && w2 > 3000 ;
86
87
88 if ( !cut ) vetoEvent ;
89
90 const DISLepton& dl = applyProjection<DISLepton>(event,"Lepton");
91 if ( dl.failed() ) vetoEvent;
92 /*
93 cout << " scattered lepton angle " << 180.- dl.out().momentum().angle(dl.in().momentum())/degree << endl;
94 cout << " in lepton " << dl.in().momentum() << endl;
95 cout << " out lepton " << dl.out().momentum() << endl;
96 */
97
98 const FinalState& fs = apply<FinalState>(event, "FS");
99 Particles particles; particles.reserve(fs.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
108 double efwd = 0.;
109
110 for (const Particle& p : particles) {
111 const double th = 180. - p.momentum().angle(dl.in().momentum())/degree;
112
113 if (inRange(th, 4.4, 15.0)) {
114 efwd += p.E();
115 //cout << " angle " << th << " pid " << p.pid() << " Efwd = " << efwd << endl;
116 }
117 }
118
119 bool evcut[2];
120 evcut[0] = efwd > 0.5;
121
122 // fill the counter
123 _Nevt_after_cuts -> fill();
124 if (Q2 > 100 ) _Nevt_afterh_cuts -> fill();
125
126 if ( evcut[0] && (Q2 < 80 ) ) _Nevt_afterfwd_cuts -> fill();
127 if ( evcut[0] && (Q2 > 100 ) ) _Nevt_afterhfwd_cuts -> fill();
128
129 double n_charg = 0;
130
131 // Boost to Breit
132 const LorentzTransform breitboost = dk.boostBreit();
133
134 for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) {
135 const Particle& p = particles[ip1];
136 const FourMomentum BreMom = breitboost.transform(p.momentum());
137 // cout << BreMom.pz() << endl;
138 double x = cos(BreMom.theta());
139
140 if (Q2 < 80 ) _h["costh_lowQ"] ->fill(x);
141 if (Q2 > 100 ) _h["costh_highQ"] ->fill(x);
142
143 if (Q2 < 80 && evcut[0]) _h["costh_lowQ_noEfwd"] ->fill(x);
144 if (Q2 > 100 && evcut[0]) _h["costh_highQ_noEfwd"] ->fill(x);
145
146
147 if ( BreMom.pz() > 0. ) continue;
148 double pcal= sqrt(BreMom.px2() + BreMom.py2()+ BreMom.pz2()) ;
149 double xp = 2*pcal/(sqrt(Q2));
150 double xi = log(1/xp);
151
152 double charge = p.charge() ;
153 // cout << " charge " << charge << endl;
154
155 if (charge > 0 ) {
156 if (Q2 < 80 ) _h["xp_posCharge_lowQ"] -> fill(xp);
157 if (Q2 > 100 ) _h["xp_posCharge_highQ"] -> fill(xp);
158 } else {
159 if (Q2 < 80 ) _h["xp_negCharge_lowQ"] -> fill(xp);
160 if (Q2 > 100 ) _h["xp_negCharge_highQ"] -> fill(xp);
161 }
162
163 if (Q2 < 80 ) _h["ksi_lowQ"] -> fill(xi);
164 if (Q2 > 100 ) _h["ksi_highQ"] -> fill (xi);
165
166 n_charg = n_charg + 1;
167
168 }
169 _h["Mult_vrs_Q2_nchrg"] -> fill(Q2,n_charg) ;
170 _h["Mult_vrs_Q2_count"] -> fill(Q2) ;
171 if ( evcut[0]) _h["Mult_vrs_Q2_noEfwd_nchrg"] -> fill(Q2,n_charg) ;
172 if ( evcut[0]) _h["Mult_vrs_Q2_noEfwd_count"] -> fill(Q2) ;
173}
174
175 /// Normalise histograms etc., after the run
176 void finalize() {
177
178
179 normalize(_h["xp_posCharge_lowQ"]);
180 if(dbl(*_Nevt_afterh_cuts)>0) scale(_h["xp_posCharge_highQ"], 1.0/ *_Nevt_afterh_cuts);
181
182 normalize(_h["xp_negCharge_lowQ"]);
183 if(dbl(*_Nevt_afterh_cuts)>0) scale(_h["xp_negCharge_highQ"], 1.0/ *_Nevt_afterh_cuts);
184
185 scale(_h["costh_lowQ"], 1.0/ *_Nevt_after_cuts);
186 if(dbl(*_Nevt_afterh_cuts)>0) scale(_h["costh_highQ"], 1.0/ *_Nevt_afterh_cuts);
187 //cout << " after fwd cuts " << dbl(*_Nevt_afterfwd_cuts) << endl;
188
189 if(dbl(*_Nevt_afterfwd_cuts)>0) scale(_h["costh_lowQ_noEfwd"], 1.0/ *_Nevt_afterfwd_cuts);
190 if(dbl(*_Nevt_afterhfwd_cuts)>0) scale(_h["costh_highQ_noEfwd"], 1.0/ *_Nevt_afterhfwd_cuts);
191
192 scale(_h["ksi_lowQ"], 1.0/ *_Nevt_after_cuts);
193 if(dbl(*_Nevt_afterh_cuts)>0) scale(_h["ksi_highQ"], 1.0/ *_Nevt_afterh_cuts);
194 //cout << " Nevt " << dbl(*_Nevt_after_cuts) << endl;
195
196 divide(_h["Mult_vrs_Q2_nchrg"], _h["Mult_vrs_Q2_count"], _s["Mult_vrs_Q2"]);
197 divide(_h["Mult_vrs_Q2_noEfwd_nchrg"], _h["Mult_vrs_Q2_noEfwd_count"], _s["Mult_vrs_Q2_noEfwd"]);
198 }
199
200 ///@}
201
202
203 /// @name Histograms
204 ///@{
205 map<string, Histo1DPtr> _h;
206 map<string, Profile1DPtr> _p;
207 map<string, CounterPtr> _c;
208 map<string, Scatter2DPtr> _s;
209 CounterPtr _Nevt_after_cuts;
210 CounterPtr _Nevt_afterfwd_cuts;
211 CounterPtr _Nevt_afterh_cuts;
212 CounterPtr _Nevt_afterhfwd_cuts;
213 ///@}
214
215
216 };
217
218
219 RIVET_DECLARE_PLUGIN(H1_1995_I394793);
220
221}
|