9#include <framework/logging/Logger.h>
10#include <framework/gearbox/Unit.h>
11#include <generators/particlegun/ParticleGun.h>
37 double randomPolyline(
size_t n,
const double* x,
const double* y)
41 std::vector<double> weights(n - 1);
43 for (
size_t i = 0; i < n - 1; ++i) {
44 weights[i] = (x[i + 1] - x[i]) * max(y[i], y[i + 1]);
50 double weight = gRandom->Uniform(0, sumw);
52 for (; segment < n - 1; ++segment) {
53 weight -= weights[segment];
54 if (weight <= 0)
break;
56 const double x1 = x[segment];
57 const double x2 = x[segment + 1];
58 const double xx = gRandom->Uniform(x1, x2);
60 const double y1 = y[segment];
61 const double y2 = y[segment + 1];
62 const double yy = (y2 == y1) ? y1 : y1 + (xx - x1) / (x2 - x1) * (y2 - y1);
64 const double randY = gRandom->Uniform(0, max(y1, y2));
66 if (randY < yy)
return xx;
81 return gRandom->Uniform(params[0], params[1]);
83 return 1. / gRandom->Uniform(1. / params[1], 1. / params[0]);
85 return acos(gRandom->Uniform(cos(params[0]), cos(params[1])));
88 return exp(gRandom->Uniform(log(params[0]), log(params[1])));
91 return gRandom->Gaus(params[0], params[1]);
93 return acos(gRandom->Gaus(params[0], params[1]));
98 for (
size_t i = params.size() / 2; i < params.size(); ++i) rand += params[i];
100 rand = gRandom->Uniform(0, rand);
102 for (
size_t i = params.size() / 2; i < params.size(); ++i) {
105 if (rand <= 0)
return params[i - (params.size() / 2)];
107 B2FATAL(
"Something wrong with picking fixed spectra values");
110 return randomPolyline(params.size() / 2, params.data(), params.data() + params.size() / 2);
112 return acos(randomPolyline(params.size() / 2, params.data(), params.data() + params.size() / 2));
114 B2FATAL(
"Unknown distribution");
130 int nTracks =
static_cast<int>(
m_params.nTracks);
133 }
else if (
m_params.varyNumberOfTracks) {
134 nTracks = gRandom->Poisson(
m_params.nTracks);
138 double firstMomentum;
139 for (
int i = 0; i < nTracks; ++i) {
140 MCParticleGraph::GraphParticle& p = graph.
addParticle();
142 if (
m_params.pdgCodes.size() == 1) {
150 int index =
static_cast<int>(gRandom->Uniform(
m_params.pdgCodes.size()));
154 p.setFirstDaughter(0);
155 p.setLastDaughter(0);
161 if (
m_params.fixedMomentumPerEvent) {
162 i == 0 ? firstMomentum = momentum : momentum = firstMomentum;
167 double pt = momentum * sin(theta);
173 momentum = (sin(theta) > 0) ? (pt / sin(theta)) : numeric_limits<double>::max();
176 double pz = momentum * cos(theta);
177 double px = pt * cos(phi);
178 double py = pt * sin(phi);
179 double m = p.getMass();
180 double e =
sqrt(momentum * momentum + m * m);
182 p.setMomentum(px, py, pz);
184 p.setProductionVertex(vx, vy, vz);
189 p.setDecayTime(numeric_limits<double>::infinity());
192 p.setProductionTime(timeOffset);
217 std::map<EDistribution, std::string> distributionNames = {
235 auto getDist = [&p](
const std::string & dist) ->
EDistribution {
236 if (dist ==
"momentum")
return p.momentumDist;
237 if (dist ==
"xVertex")
return p.xVertexDist;
238 if (dist ==
"yVertex")
return p.yVertexDist;
239 if (dist ==
"zVertex")
return p.zVertexDist;
240 if (dist ==
"theta")
return p.thetaDist;
241 if (dist ==
"phi")
return p.phiDist;
242 if (dist ==
"timeOffset")
return p.timeDist;
243 throw std::runtime_error(
"wrong parameter");
246 auto getPars = [&p](
const std::string & dist) ->
const std::vector<double>& {
247 if (dist ==
"momentum")
return p.momentumParams;
248 if (dist ==
"xVertex")
return p.xVertexParams;
249 if (dist ==
"yVertex")
return p.yVertexParams;
250 if (dist ==
"zVertex")
return p.zVertexParams;
251 if (dist ==
"theta")
return p.thetaParams;
252 if (dist ==
"phi")
return p.phiParams;
253 if (dist ==
"timeOffset")
return p.timeParams;
254 throw std::runtime_error(
"wrong parameter");
258 auto excludeDistribution = [&](
const std::string & dist,
EDistribution excluded) {
259 if (getDist(dist) == excluded) {
260 B2ERROR(distributionNames[excluded] <<
" is not allowed for " << dist <<
" generation");
266 for (
auto dist : {
"momentum",
"xVertex",
"yVertex",
"zVertex"}) {
271 for (
auto dist : {
"xVertex",
"yVertex",
"zVertex",
"phi",
"theta",
"timeOffset"}) {
280 if (p.pdgCodes.empty()) {
281 B2ERROR(
"No pdgCodes specified for generation");
286 for (
auto par : {
"momentum",
"xVertex",
"yVertex",
"zVertex",
"theta",
"phi",
"timeOffset"}) {
290 const size_t hasParams = getPars(par).size();
291 if (hasParams < minParams) {
292 B2ERROR(par <<
" generation: " << distributionNames[dist]
293 <<
" distribution requires at least " << minParams
294 <<
" parameters, " << hasParams <<
" given");
300 if ((hasParams % 2) != 0) {
301 B2ERROR(par <<
" generation: " << distributionNames[dist] <<
" requires an even number of parameters");
306 const std::vector<double>& pvec = getPars(par);
308 ?
"weight" :
"y coordinate";
311 for (
size_t i = 0; i < (hasParams / 2) - 1; ++i) {
312 if (pvec[i] > pvec[i + 1]) {
313 B2ERROR(par <<
" generation: " << distributionNames[dist] <<
" requires x coordinates in ascending order");
319 auto minmax = std::minmax_element(pvec.begin() + hasParams / 2, pvec.end());
320 if (*minmax.first < 0) {
321 B2ERROR(par <<
" generation: " << distributionNames[dist] <<
" requires "
322 << parname <<
"s to be non-negative");
325 if (*minmax.second <= 0) {
326 B2ERROR(par <<
" generation: " << distributionNames[dist] <<
" requires at least one "
327 << parname <<
" to be positive");
336 if (p.momentumParams[0] == 0 || p.momentumParams[1] == 0) {
337 B2ERROR(
"inversePt distribution does not allow zero momentum");
Class to build, validate and sort a particle decay chain.
@ c_PrimaryParticle
bit 0: Particle is primary particle.
@ c_StableInGenerator
bit 1: Particle is stable, i.e., not decaying in the generator.
EDistribution
enum containing all known distributions available for generation of values
@ c_polylineDistribution
Distribution given as list of (x,y) points.
@ c_uniformLogPtDistribution
Uniform distribution of the logarithm of the Pt, parameters are min and max value.
@ c_normalDistribution
Normal distribution, parameters are mean and sigma.
@ c_normalCosDistribution
Normal distribution of the cosinge, parameters are mean and sigma.
@ c_uniformPtDistribution
Uniform distribution of Pt, parameters are min and max value.
@ c_polylinePtDistribution
Same as polylineDistribution but for the transverse momentum.
@ c_discretePtSpectrum
Discrete pt spectrum, parameters are first the values and then the weights (non-negative) for each va...
@ c_uniformDistribution
Uniform distribution, parameters are min and max value.
@ c_discreteSpectrum
Discrete spectrum, parameters are first the values and then the weights (non-negative) for each value...
@ c_fixedValue
Fixed value, no random generation at all, 1 parameter.
@ c_normalPtDistribution
Normal distribution of Pt, parameters are mean and sigma.
@ c_uniformCosDistribution
Uniform distribution of the cosine of the values, parameters are min and max value.
@ c_inversePtDistribution
Distribution uniform in the inverse pt distribution, that is uniform in track curvature.
@ c_polylineCosDistribution
Same as polylineDistribution but for the cosine of the angle.
@ c_uniformLogDistribution
Uniform distribution of the logarithm of the values, parameters are min and max value.
Parameters m_params
All relevant parameters.
bool generateEvent(MCParticleGraph &graph)
Generate the next event and store the result in the given MCParticle graph.
bool setParameters(const Parameters ¶meters)
Set the parameters for generating the Particles.
double generateValue(EDistribution dist, const std::vector< double > ¶ms)
Generate a value according to the given distribution with the given parameters.
double sqrt(double a)
sqrt for double
GraphParticle & addParticle()
Add new particle to the graph.
Abstract base class for different kinds of events.