cepgen is hosted by Hepforge, IPPP Durham
CepGen 1.2.5
Central exclusive processes event generator
Loading...
Searching...
No Matches
CLAS.cpp
Go to the documentation of this file.
1/*
2 * CepGen: a central exclusive processes event generator
3 * Copyright (C) 2017-2024 Laurent Forthomme
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19#include <array>
20#include <cmath>
21
24#include "CepGen/Physics/PDG.h"
27
28namespace cepgen {
29 namespace strfun {
32 class CLAS : public Parameterisation {
33 public:
34 explicit CLAS(const ParametersList& params) : Parameterisation(params), mpi0_(PDG::get().mass(PDG::piZero)) {
35 const auto& model = steer<std::string>("model");
36 if (model == "proton")
37 mod_params_ = Parameters::standard_proton();
38 else if (model == "neutron")
39 mod_params_ = Parameters::standard_neutron();
40 else if (model == "deuteron")
41 mod_params_ = Parameters::standard_deuteron();
42 else
43 throw CG_FATAL("CLAS") << "Invalid modelling selected: " << model << "!";
44 }
45
48 desc.setDescription("CLAS");
49 desc.add<std::string>("model", "proton")
50 .setDescription("Nucleon modelling ('proton', 'deuteron', or 'neutron' handled)")
51 .allow("proton", "parton-from-proton")
52 .allow("deuteron", "parton-from-deuteron")
53 .allow("neutron", "parton-from-neutron");
54 return desc;
55 }
56
58 struct Parameters {
62
64 struct Resonance {
65 double amplitude{0.}, mass{0.}, width{0.};
67 };
68
69 enum { neutron = 0, proton = 1, deuteron = 2 } mode{proton};
70 std::array<double, 7> c_slac{};
71 // CLAS parameterisation
72 double alpha{0.}, beta{0.}, mu{0.}, mup{0.};
73 std::array<double, 3> x{};
74 std::array<double, 4> b{};
75 std::vector<Resonance> resonances;
76 };
77
78 void eval() override {
79 const double w2 = utils::mX2(args_.xbj, args_.q2, mp2_), w = sqrt(w2);
80 if (w < mx_min_)
81 return;
82 const auto rb = resbkg(args_.q2, w);
83 setF2(f2slac(args_.xbj, args_.q2) * (rb.first + rb.second));
84 }
85
86 private:
90 std::pair<double, double> resbkg(double q2, double w) const;
96 double f2slac(double xbj, double q2) const;
97 Parameters mod_params_;
98 static constexpr double COEFF = 6.08974;
99 double mpi0_;
100 };
101
103 Parameters params;
104 // SLAC fit parameters
105 params.c_slac = {{0.25615, 2.1785, 0.89784, -6.7162, 3.7557, 1.6421, 0.37636}};
106 // CLAS parameterisation
107 params.x = {{-0.599937, 4.76158, 0.411676}};
108 params.b = {{0.755311, 3.35065, 3.51024, 1.74470}};
109 params.alpha = -0.174985;
110 params.beta = 0.00967019;
111 params.mu = -0.0352567;
112 params.mup = 3.51852;
113
115 r0.amplitude = 1.04;
116 r0.mass = 1.22991;
117 r0.width = 0.106254;
118 r0.angular_momentum = 1;
119 params.resonances.emplace_back(r0);
120
122 r1.amplitude = 0.481327;
123 r1.mass = 1.51015;
124 r1.width = 0.0816620;
125 r1.angular_momentum = 2;
126 params.resonances.emplace_back(r1);
127
129 r2.amplitude = 0.655872;
130 r2.mass = 1.71762;
131 r2.width = 0.125520;
132 r2.angular_momentum = 3;
133 params.resonances.emplace_back(r2);
134
136 r3.amplitude = 0.747338;
137 r3.mass = 1.95381;
138 r3.width = 0.198915;
139 r3.angular_momentum = 2;
140 params.resonances.emplace_back(r3);
141
142 return params;
143 }
144
146 Parameters params = standard_proton();
147 params.mode = Parameters::neutron;
148 params.c_slac = {{0.0640, 0.2250, 4.1060, -7.0790, 3.0550, 1.6421, 0.37636}};
149 return params;
150 }
151
153 Parameters params = standard_proton();
154 params.mode = Parameters::deuteron;
155 params.c_slac = {{0.47709, 2.1602, 3.6274, -10.470, 4.9272, 1.5121, 0.35115}};
156 params.x = {{-0.21262, 6.9690, 0.40314}};
157 params.b = {{0.76111, 4.1470, 3.7119, 1.4218}};
158 params.alpha = -0.24480;
159 params.beta = 0.014503;
160
161 params.resonances.clear();
162
164 r0.amplitude = 0.74847;
165 r0.mass = 1.2400;
166 r0.width = 0.12115;
167 r0.angular_momentum = 1;
168 params.resonances.emplace_back(r0);
169
171 r1.amplitude = 0.011500;
172 r1.mass = 1.4772;
173 r1.width = 0.0069580;
174 r1.angular_momentum = 2;
175 params.resonances.emplace_back(r1);
176
178 r2.amplitude = 0.12662;
179 r2.mass = 1.5233;
180 r2.width = 0.084095;
181 r2.angular_momentum = 3;
182 params.resonances.emplace_back(r2);
183
185 r3.amplitude = 0.747338;
186 r3.mass = 1.95381;
187 r3.width = 0.198915;
188 r3.angular_momentum = 2;
189 params.resonances.emplace_back(r3);
190
191 return params;
192 }
193
194 double CLAS::f2slac(double xbj, double q2) const {
195 if (xbj >= 1.)
196 return 0.;
197
198 const double xsxb = (q2 + mod_params_.c_slac[6]) / (q2 + mod_params_.c_slac[5] * xbj);
199 const double xs = xbj * xsxb;
200
201 double f2 = 0.;
202 for (unsigned short i = 0; i < 5; ++i)
203 f2 += mod_params_.c_slac[i] * pow(1. - xs, i);
204
205 if (mod_params_.mode == Parameters::deuteron && xbj > 0.)
206 f2 /= (1. - exp(-7.70 * (1. / xbj - 1. + mp2_ / q2)));
207
208 return f2 * pow(1. - xs, 3) / xsxb;
209 }
210
211 std::pair<double, double> CLAS::resbkg(double q2, double w) const {
212 const double mpi02 = mpi0_ * mpi0_;
213
214 if (w < mx_min_)
215 return std::make_pair(0., 0.);
216 if (w > 4.)
217 return std::make_pair(1., 0.);
218
219 const double w2 = w * w;
220
221 double qs = pow(w2 + mp2_ - mpi02, 2) - 4. * mp2_ * w2;
222 if (qs <= 0.)
223 return std::make_pair(1., 0.);
224 qs = 0.5 * sqrt(qs) / w;
225
226 const double omega = 0.5 * (w2 + q2 - mp2_) / mp_;
227 const double xn = 0.5 * q2 / (mp_ * omega);
228
229 const double bkg2 =
230 (w > mod_params_.b[3]) ? exp(-mod_params_.b[2] * (w2 - mod_params_.b[3] * mod_params_.b[3])) : 1.;
231
232 double f2bkg =
233 (mod_params_.b[0]) * (1. - exp(-mod_params_.b[1] * (w - mx_min_))) + (1. - mod_params_.b[0]) * (1. - bkg2);
234 f2bkg *= (1. + (1. - f2bkg) * (mod_params_.x[0] + mod_params_.x[1] * pow(xn - mod_params_.x[2], 2)));
235
236 double etab = 1., etad = 1.;
237 if (mod_params_.mode != Parameters::deuteron && q2 <= 2. && w <= 2.5) {
238 etab = 1. - 2.5 * q2 * exp(-12.5 * q2 * q2 - 50. * (w - 1.325) * (w - 1.325));
239 etad = 1. + 2.5 * q2 * exp(-12.5 * q2 * q2);
240 }
241 f2bkg *= etab;
242
243 double f2resn = 0.;
244
245 unsigned short i = 0;
246 for (const auto& res : mod_params_.resonances) {
247 const double ai = (i == 0)
248 ? etad * (res.amplitude + q2 * std::min(0., mod_params_.alpha + mod_params_.beta * q2))
249 : res.amplitude;
250 const double dmi = (i == 2) ? res.mass * (1. + mod_params_.mu / (1. + mod_params_.mup * q2)) : res.mass;
251 double qs0 = pow(dmi * dmi + mp2_ - mpi02, 2) - 4. * mp2_ * dmi * dmi;
252 if (qs0 <= 0.)
253 break;
254 qs0 = 0.5 * sqrt(qs0) / dmi;
255 int ji = 2 * res.angular_momentum;
256 const double dg =
257 0.5 * res.width * pow(qs / qs0, ji + 1) * (1. + pow(COEFF * qs0, ji)) / (1. + pow(COEFF * qs, ji));
258 f2resn += ai * dg / ((w - dmi) * (w - dmi) + dg * dg);
259 ++i;
260 }
261 f2resn *= 0.5 * (1. - mod_params_.b[0]) * bkg2 / mp_ * M_1_PI;
262
263 return std::make_pair(f2bkg, f2resn);
264 }
265 } // namespace strfun
266} // namespace cepgen
268REGISTER_STRFUN("CLAS", 103, CLAS);
#define CG_FATAL(mod)
Definition Exception.h:61
#define REGISTER_STRFUN(name, id, obj)
Add a structure functions definition to the list of handled parameterisation.
A singleton holding all physics constants associated to particles.
Definition PDG.h:28
A description object for parameters collection.
CLAS parameterisation for nucleon data at > 0.5 GeV and > 0.15.
Definition CLAS.cpp:32
CLAS(const ParametersList &params)
Definition CLAS.cpp:34
static ParametersDescription description()
Definition CLAS.cpp:46
void eval() override
Local structure functions evaluation method.
Definition CLAS.cpp:78
Base object for the parameterisation of nucleon structure functions.
Arguments args_
Last couple computed.
Parameterisation & setF2(double f2)
const double mx_min_
Minimum diffractive mass, in GeV/c^2.
const double mp2_
Squared proton mass, in GeV^2/c^4.
static ParametersDescription description()
Generic description for the structure functions.
const double mp_
Proton mass, in GeV/c^2.
double mX2(double xbj, double q2, double mp2)
Compute the diffractive mass from virtuality/Bjorken x.
Definition Utils.cpp:29
Common namespace for this Monte Carlo generator.
List of steering parameters for a physics case.
Definition CLAS.cpp:58
std::vector< Resonance > resonances
Definition CLAS.cpp:75
enum cepgen::strfun::CLAS::Parameters::@0 proton
Nucleon type.
static Parameters standard_neutron()
Parton-from-neutron emission.
Definition CLAS.cpp:145
static Parameters standard_proton()
Parton-from-proton emission.
Definition CLAS.cpp:102
std::array< double, 4 > b
Definition CLAS.cpp:74
std::array< double, 3 > x
Definition CLAS.cpp:73
static Parameters standard_deuteron()
Parton-from-deuteron emission.
Definition CLAS.cpp:152
std::array< double, 7 > c_slac
SLAC fit parameters.
Definition CLAS.cpp:70
Physical properties associated to a resonance.
Definition CLAS.cpp:64