cepgen is hosted by Hepforge, IPPP Durham
CepGen 1.2.5
Central exclusive processes event generator
Loading...
Searching...
No Matches
Shamov.cpp
Go to the documentation of this file.
1/*
2 * CepGen: a central exclusive processes event generator
3 * Copyright (C) 2013-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 "CepGen/Modules/StructureFunctionsFactory.h"
20#include "CepGen/Physics/Constants.h"
21#include "CepGen/Physics/Utils.h"
22#include "CepGen/StructureFunctions/Parameterisation.h"
23#include "CepGen/Utils/GridHandler.h"
24
25namespace cepgen {
26 namespace strfun {
27 class Shamov final : public Parameterisation {
28 public:
29 explicit Shamov(const ParametersList&);
30
31 static ParametersDescription description() {
32 auto desc = Parameterisation::description();
33 desc.setDescription("Shamov (hybrid, soft)");
34 desc.addAs<int, Mode>("mode", Mode::RealResAndNonRes).setDescription("sub-structure functions choice");
35 desc.add<int>("fitModel", 2);
36 desc.add<double>("q20", 0.65 /* 0.36 */)
37 .setDescription("first parameter for non-resonant gamma-p cross section q^2 dependence");
38 desc.add<double>("rPower", 0.71 /* 0.52 */)
39 .setDescription("second parameter for non-resonant gamma-p cross section q^2 dependence");
40 desc.add<double>("gm0", 1.).setDescription("scaling factor for the magnetic form factor template fit");
41 desc.add<double>("gmb", 0.984)
42 .setDescription("exponential parameter for the magnetic form factor template fit");
43 desc.add<double>("lowQ2", 1.e-7);
44 return desc;
45 }
46
47 void eval() override;
48
49 private:
50 static constexpr double prefac_ =
51 2. * M_PI * M_PI * constants::ALPHA_EM * constants::GEVM2_TO_PB * 1e-9; // pb/GeV -> mb/GeV
52 static constexpr std::array<float, 19> gmq_ = {{0.000,
53 0.200,
54 0.300,
55 0.400,
56 0.470,
57 0.480,
58 0.500,
59 0.600,
60 0.630,
61 0.631,
62 0.770,
63 0.780,
64 0.790,
65 0.970,
66 0.980,
67 1.150,
68 1.340,
69 1.570,
70 2.340}};
71 static constexpr std::array<float, 19> gmv_ = {{3.000,
72 1.770,
73 1.380,
74 1.170,
75 0.978,
76 0.961,
77 0.964,
78 0.766,
79 0.735,
80 0.719,
81 0.570,
82 0.572,
83 0.553,
84 0.460,
85 0.446,
86 0.326,
87 0.269,
88 0.209,
89 0.102}};
90 static constexpr std::array<float, 291> gp_en_ = {
91 {1.070, 1.100, 1.130, 1.150, 1.174, 1.194, 1.213, 1.232, 1.251, 1.270, 1.288, 1.306, 1.324, 1.342, 1.359,
92 1.376, 1.393, 1.410, 1.426, 1.443, 1.456, 1.459, 1.475, 1.491, 1.506, 1.522, 1.537, 1.543, 1.552, 1.567,
93 1.582, 1.597, 1.612, 1.617, 1.626, 1.640, 1.655, 1.669, 1.683, 1.697, 1.710, 1.724, 1.738, 1.743, 1.751,
94 1.764, 1.778, 1.791, 1.804, 1.817, 1.830, 1.831, 1.843, 1.855, 1.868, 1.880, 1.893, 1.898, 1.905, 1.917,
95 1.922, 1.930, 1.942, 1.954, 1.966, 1.978, 1.989, 2.001, 2.013, 2.025, 2.036, 2.041, 2.048, 2.059, 2.070,
96 2.082, 2.093, 2.104, 2.115, 2.117, 2.126, 2.137, 2.148, 2.153, 2.154, 2.159, 2.170, 2.174, 2.181, 2.191,
97 2.202, 2.213, 2.223, 2.234, 2.238, 2.244, 2.255, 2.265, 2.275, 2.286, 2.296, 2.300, 2.306, 2.316, 2.320,
98 2.326, 2.336, 2.346, 2.356, 2.360, 2.366, 2.376, 2.386, 2.396, 2.400, 2.406, 2.415, 2.419, 2.425, 2.435,
99 2.444, 2.454, 2.464, 2.473, 2.477, 2.483, 2.492, 2.501, 2.511, 2.513, 2.520, 2.529, 2.533, 2.539, 2.548,
100 2.551, 2.557, 2.566, 2.575, 2.584, 2.593, 2.602, 2.611, 2.620, 2.624, 2.629, 2.638, 2.642, 2.647, 2.649,
101 2.656, 2.665, 2.674, 2.682, 2.691, 2.695, 2.700, 2.708, 2.717, 2.726, 2.734, 2.743, 2.746, 2.751, 2.760,
102 2.763, 2.768, 2.777, 2.785, 2.794, 2.797, 2.802, 2.810, 2.819, 2.827, 2.830, 2.835, 2.844, 2.847, 2.852,
103 2.860, 2.868, 2.876, 2.877, 2.885, 2.893, 2.896, 2.897, 2.901, 2.909, 2.917, 2.919, 2.925, 2.933, 2.941,
104 2.944, 2.949, 2.957, 2.958, 2.965, 3.032, 3.035, 3.038, 3.114, 3.130, 3.147, 3.207, 3.218, 3.253, 3.296,
105 3.305, 3.383, 3.389, 3.416, 3.471, 3.479, 3.540, 3.551, 3.582, 3.634, 3.683, 3.784, 3.867, 3.868, 3.911,
106 3.937, 4.029, 4.087, 4.144, 4.204, 4.257, 4.282, 4.327, 4.334, 4.379, 4.396, 4.458, 4.514, 4.580, 4.581,
107 4.645, 4.715, 4.773, 4.860, 4.912, 4.931, 5.001, 5.065, 5.140, 5.141, 5.213, 5.356, 5.374, 5.541, 5.623,
108 5.704, 5.862, 5.935, 6.665, 7.271, 7.672, 7.733, 8.066, 8.329, 8.485, 9.125, 9.186, 9.576, 9.979, 15.10,
109 20.00, 30.00, 40.00, 50.00, 60.00, 70.00, 80.00, 90.00, 100.0, 120.0, 150.0, 170.0, 200.0, 240.0, 270.0,
110 310.0, 350.0, 400.0, 500.0, 7000., 14000.}};
111 static constexpr std::array<float, 291> gp_cs_ = {
112 {0.000, 0.052, 0.108, 0.175, 0.424, 0.487, 0.527, 0.478, 0.407, 0.334, 0.244, 0.225, 0.200, 0.178, 0.177,
113 0.187, 0.194, 0.212, 0.223, 0.233, 0.211, 0.240, 0.265, 0.279, 0.276, 0.261, 0.245, 0.201, 0.221, 0.206,
114 0.214, 0.209, 0.202, 0.193, 0.205, 0.201, 0.212, 0.218, 0.215, 0.192, 0.191, 0.175, 0.165, 0.182, 0.159,
115 0.162, 0.150, 0.149, 0.144, 0.156, 0.150, 0.147, 0.154, 0.154, 0.154, 0.147, 0.154, 0.154, 0.144, 0.152,
116 0.151, 0.156, 0.154, 0.146, 0.139, 0.156, 0.150, 0.150, 0.145, 0.139, 0.145, 0.146, 0.142, 0.141, 0.142,
117 0.143, 0.149, 0.154, 0.135, 0.115, 0.148, 0.144, 0.144, 0.143, 0.142, 0.149, 0.144, 0.148, 0.138, 0.132,
118 0.145, 0.138, 0.145, 0.136, 0.138, 0.138, 0.139, 0.136, 0.129, 0.136, 0.140, 0.144, 0.133, 0.139, 0.139,
119 0.143, 0.140, 0.140, 0.139, 0.134, 0.141, 0.130, 0.136, 0.124, 0.132, 0.128, 0.130, 0.142, 0.132, 0.134,
120 0.139, 0.133, 0.144, 0.133, 0.135, 0.136, 0.130, 0.134, 0.134, 0.111, 0.130, 0.131, 0.130, 0.129, 0.140,
121 0.134, 0.138, 0.129, 0.144, 0.128, 0.133, 0.132, 0.127, 0.128, 0.127, 0.124, 0.124, 0.133, 0.127, 0.127,
122 0.121, 0.134, 0.129, 0.134, 0.123, 0.129, 0.132, 0.121, 0.137, 0.123, 0.130, 0.135, 0.127, 0.129, 0.128,
123 0.126, 0.123, 0.122, 0.120, 0.119, 0.128, 0.134, 0.132, 0.125, 0.122, 0.126, 0.114, 0.127, 0.135, 0.122,
124 0.140, 0.125, 0.136, 0.116, 0.136, 0.124, 0.131, 0.127, 0.146, 0.136, 0.120, 0.131, 0.142, 0.132, 0.129,
125 0.132, 0.134, 0.133, 0.133, 0.127, 0.128, 0.116, 0.125, 0.127, 0.124, 0.122, 0.134, 0.122, 0.128, 0.122,
126 0.122, 0.130, 0.118, 0.125, 0.124, 0.121, 0.116, 0.122, 0.129, 0.122, 0.120, 0.118, 0.126, 0.122, 0.124,
127 0.124, 0.122, 0.120, 0.119, 0.123, 0.115, 0.124, 0.114, 0.123, 0.117, 0.120, 0.122, 0.128, 0.124, 0.112,
128 0.128, 0.121, 0.122, 0.115, 0.115, 0.118, 0.116, 0.111, 0.115, 0.118, 0.119, 0.113, 0.113, 0.115, 0.114,
129 0.113, 0.115, 0.117, 0.115, 0.114, 0.114, 0.114, 0.115, 0.112, 0.114, 0.115, 0.115, 0.115, 0.114, 0.113,
130 0.114, 0.115, 0.115, 0.114, 0.114, 0.117, 0.116, 0.116, 0.116, 0.114, 0.117, 0.116, 0.114, 0.116, 0.120,
131 0.116, 0.118, 0.143, 0.143, 0.143, 0.143}};
132 static constexpr std::array<float, 291> gp_nr_ = {
133 {0.683, 0.539, 0.374, 0.267, 0.378, 0.143, 0.085, 0.137, 0.233, 0.297, 0.249, 0.333, 0.364, 0.365, 0.408,
134 0.460, 0.476, 0.494, 0.472, 0.431, 0.303, 0.374, 0.376, 0.384, 0.394, 0.412, 0.444, 0.357, 0.458, 0.486,
135 0.553, 0.576, 0.577, 0.559, 0.580, 0.545, 0.530, 0.531, 0.568, 0.590, 0.655, 0.677, 0.699, 0.740, 0.720,
136 0.750, 0.749, 0.765, 0.771, 0.800, 0.803, 0.799, 0.817, 0.825, 0.831, 0.830, 0.844, 0.846, 0.839, 0.852,
137 0.853, 0.860, 0.863, 0.859, 0.856, 0.876, 0.874, 0.877, 0.876, 0.873, 0.881, 0.883, 0.881, 0.883, 0.886,
138 0.889, 0.895, 0.901, 0.889, 0.870, 0.900, 0.900, 0.901, 0.901, 0.900, 0.906, 0.904, 0.908, 0.902, 0.899,
139 0.909, 0.906, 0.912, 0.908, 0.910, 0.910, 0.912, 0.911, 0.908, 0.914, 0.917, 0.920, 0.914, 0.919, 0.919,
140 0.922, 0.921, 0.922, 0.923, 0.920, 0.925, 0.919, 0.924, 0.917, 0.923, 0.921, 0.923, 0.929, 0.925, 0.927,
141 0.930, 0.927, 0.934, 0.929, 0.930, 0.931, 0.929, 0.931, 0.932, 0.918, 0.931, 0.932, 0.931, 0.931, 0.937,
142 0.935, 0.937, 0.933, 0.940, 0.934, 0.937, 0.937, 0.935, 0.936, 0.936, 0.935, 0.935, 0.939, 0.937, 0.937,
143 0.935, 0.941, 0.939, 0.942, 0.937, 0.941, 0.942, 0.937, 0.945, 0.939, 0.943, 0.945, 0.942, 0.943, 0.943,
144 0.942, 0.941, 0.941, 0.941, 0.941, 0.945, 0.948, 0.947, 0.944, 0.943, 0.946, 0.940, 0.946, 0.950, 0.945,
145 0.952, 0.947, 0.951, 0.943, 0.952, 0.947, 0.950, 0.948, 0.955, 0.952, 0.946, 0.951, 0.955, 0.952, 0.951,
146 0.952, 0.953, 0.953, 0.953, 0.951, 0.954, 0.949, 0.953, 0.956, 0.955, 0.955, 0.960, 0.957, 0.960, 0.959,
147 0.959, 0.963, 0.960, 0.962, 0.963, 0.962, 0.962, 0.964, 0.966, 0.966, 0.966, 0.967, 0.970, 0.969, 0.970,
148 0.971, 0.971, 0.971, 0.972, 0.973, 0.972, 0.974, 0.973, 0.975, 0.974, 0.975, 0.976, 0.977, 0.977, 0.975,
149 0.978, 0.978, 0.978, 0.978, 0.978, 0.979, 0.979, 0.978, 0.979, 0.980, 0.981, 0.980, 0.980, 0.982, 0.982,
150 0.982, 0.983, 0.984, 0.986, 0.987, 0.988, 0.988, 0.989, 0.989, 0.990, 0.991, 0.991, 0.991, 0.992, 0.992,
151 0.992, 0.993, 0.993, 0.993, 0.993, 0.994, 0.994, 0.994, 0.994, 0.995, 0.995, 0.995, 0.995, 0.995, 0.996,
152 0.996, 0.996, 1.000, 1.000, 1.000, 1.000}};
153
154 const enum class Mode {
155 SuriYennie = 0,
158 RealRes = 1,
162 RealResAndNonRes = 2,
165 RealAndSuriYennieNonRes = 3,
169 RealAndFitNonRes = 4
170 } mode_;
171 size_t fit_model_;
172 const double gm0_;
173 const double gmb_;
174 const double q20_;
175 const double r_power_;
176 const double lowq2_;
177
178 GridHandler<1, 2> sigma_grid_{GridType::linear};
179 GridHandler<1, 1> gm_grid_{GridType::linear};
180 std::unique_ptr<Parameterisation> sy_sf_;
181 bool non_resonant_{true};
182 };
183
184 constexpr std::array<float, 19> Shamov::gmq_, Shamov::gmv_;
185 constexpr std::array<float, 291> Shamov::gp_en_, Shamov::gp_cs_, Shamov::gp_nr_;
186
187 Shamov::Shamov(const ParametersList& params)
188 : Parameterisation(params),
189 mode_(steerAs<int, Mode>("mode")),
190 fit_model_(steer<int>("fitModel")),
191 gm0_(steer<double>("gm0")),
192 gmb_(steer<double>("gmb")),
193 q20_(steer<double>("q20")),
194 r_power_(steer<double>("rPower")),
195 lowq2_(steer<double>("lowQ2")),
196 sy_sf_(StructureFunctionsFactory::get().build("SuriYennie", steer<ParametersList>("syParams"))) {
197 //----- initialise the interpolation grids
198
199 //--- grid E -> (cross section, norm)
200 for (size_t i = 0; i < gp_en_.size(); ++i)
201 sigma_grid_.insert({gp_en_.at(i)}, {gp_cs_.at(i), gp_nr_.at(i)});
202 sigma_grid_.initialise();
203
204 //--- grid Q -> gamma_v
205 for (size_t i = 0; i < gmv_.size(); ++i)
206 gm_grid_.insert({gmq_.at(i)}, {gmv_.at(i)});
207 gm_grid_.initialise();
208 if (mode_ == Mode::SuriYennie || mode_ == Mode::RealAndSuriYennieNonRes || mode_ == Mode::RealResAndNonRes ||
209 mode_ == Mode::RealAndFitNonRes)
210 non_resonant_ = true;
211 }
212
213 void Shamov::eval() {
214 //--- Suri & Yennie structure functions
215 const double mx2 = utils::mX2(args_.xbj, args_.q2, mp2_), mx = sqrt(mx2);
216 if (mode_ == Mode::SuriYennie) { // trivial composite case
217 setW1(sy_sf_->W1(args_.xbj, args_.q2));
218 setW2(sy_sf_->W2(args_.xbj, args_.q2));
219 } else {
220 const auto sigma = sigma_grid_.eval({mx});
221 double sgp = sigma[0]; // cross section value at MX
222
223 if (mode_ == Mode::RealAndFitNonRes && mx > 1.5)
224 non_resonant_ = true;
225 //if (mx > sigma_grid_.max()[0] || rand() * 1. / RAND_MAX < sigma[1])
226 // non_resonant_ = true;
227
228 //--- q^2 dependence of sigma
229 double Gm;
230 if (non_resonant_) {
231 if (mode_ == Mode::RealAndFitNonRes)
234 Gm = std::pow(1. + std::pow(args_.q2 / q20_, 2), -r_power_);
235 else
236 Gm = sy_sf_->W1(args_.xbj, args_.q2) / sy_sf_->W1(args_.xbj, lowq2_);
237 } else { // resonant
238 if (args_.q2 >= gm_grid_.max()[0]) // above grid range
239 Gm = gm0_ * exp(-gmb_ * args_.q2);
240 else
241 Gm = gm_grid_.eval({args_.q2})[0];
242 Gm /= 3.; // due to data normalization
243 }
244 sgp *= Gm; // cross section with some q^2 dependence
245
246 //--- for W -> cross section
247 const double s1 = prefac_ * 4. * mp_ / (mx2 - mp2_); // mb/GeV
248 const double s2 = prefac_ * (std::pow(mx2 - mp2_, 2) + 2. * (mx2 + mp2_) * args_.q2 + args_.q2 * args_.q2) /
249 mp_ / args_.q2 / (mx2 - mp2_); // mb/GeV
250 //--- ratio (\sigma_T+\sigma_L)/\sigma_T according to S & Y
251 const double ratio = (s2 * sy_sf_->W2(args_.xbj, args_.q2)) / (s1 * sy_sf_->W1(args_.xbj, args_.q2));
252
253 if (mode_ == Mode::RealAndFitNonRes) { // transverse sigma only
254 setW1(sgp / s1);
255 setW2(sgp * ratio / s2);
256 } else { // longitudinal + transverse component for sigma
257 setW1(sgp / ratio / s1);
258 setW2(sgp / s2);
259 }
260 }
261 const double nu = 0.5 * (args_.q2 + mx2 - mp2_) / mp_;
262 setF2(W2(args_.xbj, args_.q2) * nu / mp_);
263 }
264 } // namespace strfun
265} // namespace cepgen
266using cepgen::strfun::Shamov;
267REGISTER_STRFUN("Shamov", 302, Shamov);
REGISTER_STRFUN
#define REGISTER_STRFUN(name, id, obj)
Add a structure functions definition to the list of handled parameterisation.
Definition StructureFunctionsFactory.h:25
cepgen::GridHandler::initialise
void initialise()
Initialise the grid and all useful interpolators/accelerators.
Definition GridHandler.cpp:150
cepgen::GridHandler::insert
void insert(coord_t coord, values_t value)
Insert a new value in the grid.
Definition GridHandler.cpp:129
cepgen::GridHandler::max
std::array< double, D > max() const
Highest bound of the grid coordinates.
Definition GridHandler.cpp:288
cepgen::GridHandler::eval
values_t eval(coord_t in_coords) const
Interpolate a point to a given coordinate.
Definition GridHandler.cpp:37
cepgen::ParametersDescription
A description object for parameters collection.
Definition ParametersDescription.h:26
cepgen::strfun::Parameterisation
Base object for the parameterisation of nucleon structure functions.
Definition Parameterisation.h:30
cepgen::strfun::Parameterisation::args_
Arguments args_
Last couple computed.
Definition Parameterisation.h:109
cepgen::strfun::Parameterisation::setF2
Parameterisation & setF2(double f2)
Definition Parameterisation.cpp:83
cepgen::strfun::Parameterisation::setW2
Parameterisation & setW2(double w2)
Definition Parameterisation.cpp:99
cepgen::strfun::Parameterisation::mp2_
const double mp2_
Squared proton mass, in GeV^2/c^4.
Definition Parameterisation.h:106
cepgen::strfun::Parameterisation::W2
double W2(double xbj, double q2)
Definition Parameterisation.cpp:67
cepgen::strfun::Parameterisation::description
static ParametersDescription description()
Generic description for the structure functions.
Definition Parameterisation.cpp:148
cepgen::strfun::Parameterisation::setW1
Parameterisation & setW1(double w1)
Definition Parameterisation.cpp:94
cepgen::strfun::Parameterisation::mp_
const double mp_
Proton mass, in GeV/c^2.
Definition Parameterisation.h:105
cepgen::strfun::Shamov::description
static ParametersDescription description()
Definition Shamov.cpp:31
cepgen::strfun::Shamov::Shamov
Shamov(const ParametersList &)
Definition Shamov.cpp:187
cepgen::strfun::Shamov::eval
void eval() override
Local structure functions evaluation method.
Definition Shamov.cpp:213
cepgen::constants::GEVM2_TO_PB
constexpr double GEVM2_TO_PB
Conversion factor between GeV and barn i.e. in GeV .
Definition Constants.h:37
cepgen::constants::ALPHA_EM
constexpr double ALPHA_EM
Electromagnetic coupling constant .
Definition Constants.h:28
cepgen::utils::mX2
double mX2(double xbj, double q2, double mp2)
Compute the diffractive mass from virtuality/Bjorken x.
Definition Utils.cpp:29
cepgen
Common namespace for this Monte Carlo generator.
Definition CommandLineHandler.cpp:36