File indexing completed on 2025-08-05 08:18:11
0001 #include "PHG4SimpleEventGenerator.h"
0002
0003 #include "PHG4InEvent.h"
0004 #include "PHG4Particle.h" // for PHG4Particle
0005 #include "PHG4Particlev2.h"
0006 #include "PHG4Utils.h"
0007
0008 #include <fun4all/Fun4AllReturnCodes.h>
0009
0010 #include <phool/PHCompositeNode.h>
0011 #include <phool/PHDataNode.h> // for PHDataNode
0012 #include <phool/PHNode.h> // for PHNode
0013 #include <phool/PHNodeIterator.h> // for PHNodeIterator
0014 #include <phool/PHObject.h> // for PHObject
0015 #include <phool/getClass.h>
0016 #include <phool/phool.h> // for PHWHERE
0017
0018 #include <TSystem.h>
0019
0020 #include <gsl/gsl_randist.h>
0021 #include <gsl/gsl_rng.h> // for gsl_rng_uniform_pos
0022
0023 #include <cassert>
0024 #include <cmath>
0025 #include <cstdlib>
0026 #include <iostream> // for operator<<, endl, basic_ostream
0027 #include <memory> // for allocator_traits<>::value_type
0028
0029 PHG4SimpleEventGenerator::PHG4SimpleEventGenerator(const std::string &name)
0030 : PHG4ParticleGeneratorBase(name)
0031 {
0032 return;
0033 }
0034
0035 void PHG4SimpleEventGenerator::add_particles(const std::string &name, const unsigned int num)
0036 {
0037 _particle_names.emplace_back(name, num);
0038 return;
0039 }
0040
0041 void PHG4SimpleEventGenerator::add_particles(const int pid, const unsigned int num)
0042 {
0043 _particle_codes.emplace_back(pid, num);
0044 return;
0045 }
0046
0047 void PHG4SimpleEventGenerator::set_eta_range(const double min, const double max)
0048 {
0049 if (min > max)
0050 {
0051 std::cout << "not setting eta bc etamin " << min << " > etamax: " << max << std::endl;
0052 gSystem->Exit(1);
0053 }
0054 m_EtaMin = min;
0055 m_EtaMax = max;
0056 m_ThetaMin = NAN;
0057 m_ThetaMax = NAN;
0058 return;
0059 }
0060
0061 void PHG4SimpleEventGenerator::set_theta_range(const double min, const double max)
0062 {
0063 if (min > max)
0064 {
0065 std::cout << __PRETTY_FUNCTION__ << " thetamin " << min << " > thetamax: " << max << std::endl;
0066 gSystem->Exit(1);
0067 }
0068 if (min < 0 || max > M_PI)
0069 {
0070 std::cout << __PRETTY_FUNCTION__ << " min or max outside range (range is 0 to pi) min: " << min << ", max: " << max << std::endl;
0071 gSystem->Exit(1);
0072 }
0073 m_ThetaMin = min;
0074 m_ThetaMax = max;
0075 m_EtaMin = NAN;
0076 m_EtaMax = NAN;
0077 return;
0078 }
0079
0080 void PHG4SimpleEventGenerator::set_phi_range(const double min, const double max)
0081 {
0082 if (min > max)
0083 {
0084 std::cout << __PRETTY_FUNCTION__ << " phimin " << min << " > phimax: " << max << std::endl;
0085 gSystem->Exit(1);
0086 return;
0087 }
0088 if (min < -M_PI || max > M_PI)
0089 {
0090 std::cout << __PRETTY_FUNCTION__ << "min or max outside range (range is -pi to pi), min: " << min << ", max: " << max << std::endl;
0091 gSystem->Exit(1);
0092 }
0093
0094 m_PhiMin = min;
0095 m_PhiMax = max;
0096 return;
0097 }
0098
0099 void PHG4SimpleEventGenerator::set_power_law_n(const double n)
0100 {
0101 m_powerLawN = n;
0102 }
0103
0104 void PHG4SimpleEventGenerator::set_pt_range(const double min, const double max, const double pt_gaus_width)
0105 {
0106 if (min > max)
0107 {
0108 std::cout << __PRETTY_FUNCTION__ << " ptmin " << min << " > ptmax: " << max << std::endl;
0109 gSystem->Exit(1);
0110 }
0111 if (min < 0 || max < 0 || pt_gaus_width < 0)
0112 {
0113 std::cout << __PRETTY_FUNCTION__ << " values need to be >= 0, min: " << min
0114 << ", max: " << max << ", pt_gaus_width: " << pt_gaus_width << std::endl;
0115 gSystem->Exit(1);
0116 }
0117
0118 m_Pt_Min = min;
0119 m_Pt_Max = max;
0120 m_Pt_GausWidth = pt_gaus_width;
0121 m_P_Min = NAN;
0122 m_P_Max = NAN;
0123 m_P_GausWidth = NAN;
0124 return;
0125 }
0126
0127 void PHG4SimpleEventGenerator::set_p_range(const double min, const double max, const double p_gaus_width)
0128 {
0129 if (min > max)
0130 {
0131 std::cout << __PRETTY_FUNCTION__ << " pmin " << min << " > pmax: " << max << std::endl;
0132 gSystem->Exit(1);
0133 }
0134 if (min < 0 || max < 0 || p_gaus_width < 0)
0135 {
0136 std::cout << __PRETTY_FUNCTION__ << " values need to be >= 0, min: " << min
0137 << ", max: " << max << ", p_gaus_width: " << p_gaus_width << std::endl;
0138 gSystem->Exit(1);
0139 }
0140 m_Pt_Min = NAN;
0141 m_Pt_Max = NAN;
0142 m_Pt_GausWidth = NAN;
0143 m_P_Min = min;
0144 m_P_Max = max;
0145 m_P_GausWidth = p_gaus_width;
0146 return;
0147 }
0148
0149 void PHG4SimpleEventGenerator::set_vertex_distribution_function(FUNCTION x, FUNCTION y, FUNCTION z)
0150 {
0151 m_VertexFunc_x = x;
0152 m_VertexFunc_y = y;
0153 m_VertexFunc_z = z;
0154 return;
0155 }
0156
0157 void PHG4SimpleEventGenerator::set_vertex_distribution_mean(const double x, const double y, const double z)
0158 {
0159 m_Vertex_x = x;
0160 m_Vertex_y = y;
0161 m_Vertex_z = z;
0162 return;
0163 }
0164
0165 void PHG4SimpleEventGenerator::set_vertex_distribution_width(const double x, const double y, const double z)
0166 {
0167 m_VertexWidth_x = x;
0168 m_VertexWidth_y = y;
0169 m_VertexWidth_z = z;
0170 return;
0171 }
0172
0173 void PHG4SimpleEventGenerator::set_existing_vertex_offset_vector(const double x, const double y, const double z)
0174 {
0175 m_VertexOffset_x = x;
0176 m_VertexOffset_y = y;
0177 m_VertexOffset_z = z;
0178 return;
0179 }
0180
0181 void PHG4SimpleEventGenerator::set_vertex_size_parameters(const double mean, const double width)
0182 {
0183 m_VertexSizeMean = mean;
0184 m_VertexSizeWidth = width;
0185 return;
0186 }
0187
0188 int PHG4SimpleEventGenerator::InitRun(PHCompositeNode *topNode)
0189 {
0190 if (m_FunctionNames.find(m_VertexFunc_x) == m_FunctionNames.end())
0191 {
0192 std::cout << PHWHERE << "::Error - unknown x vertex distribution function requested" << std::endl;
0193 gSystem->Exit(1);
0194 }
0195 if (m_FunctionNames.find(m_VertexFunc_y) == m_FunctionNames.end())
0196 {
0197 std::cout << PHWHERE << "::Error - unknown y vertex distribution function requested" << std::endl;
0198 gSystem->Exit(1);
0199 }
0200 if (m_FunctionNames.find(m_VertexFunc_z) == m_FunctionNames.end())
0201 {
0202 std::cout << PHWHERE << "::Error - unknown z vertex distribution function requested" << std::endl;
0203 gSystem->Exit(1);
0204 }
0205
0206 m_InEvent = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
0207 if (!m_InEvent)
0208 {
0209 PHNodeIterator iter(topNode);
0210 PHCompositeNode *dstNode;
0211 dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0212
0213 m_InEvent = new PHG4InEvent();
0214 PHDataNode<PHObject> *newNode = new PHDataNode<PHObject>(m_InEvent, "PHG4INEVENT", "PHObject");
0215 dstNode->addNode(newNode);
0216 }
0217
0218 if (Verbosity() > 0)
0219 {
0220 std::cout << "================ PHG4SimpleEventGenerator::InitRun() ======================" << std::endl;
0221 std::cout << " Random seed = " << get_seed() << std::endl;
0222 std::cout << " Particles:" << std::endl;
0223 for (auto &_particle_code : _particle_codes)
0224 {
0225 std::cout << " " << _particle_code.first << ", count = " << _particle_code.second << std::endl;
0226 }
0227 for (auto &_particle_name : _particle_names)
0228 {
0229 std::cout << " " << _particle_name.first << ", count = " << _particle_name.second << std::endl;
0230 }
0231 if (get_reuse_existing_vertex())
0232 {
0233 std::cout << " Vertex Distribution: Set to reuse a previously generated sim vertex" << std::endl;
0234 std::cout << " Vertex offset vector (x,y,z) = (" << m_VertexOffset_x << "," << m_VertexOffset_y << "," << m_VertexOffset_z << ")" << std::endl;
0235 }
0236 else
0237 {
0238 std::cout << " Vertex Distribution Function (x,y,z) = ("
0239 << m_FunctionNames.find(m_VertexFunc_x)->second << ","
0240 << m_FunctionNames.find(m_VertexFunc_y)->second << ","
0241 << m_FunctionNames.find(m_VertexFunc_z)->second << ")"
0242 << std::endl;
0243 std::cout << " Vertex mean (x,y,z) = (" << m_Vertex_x << "," << m_Vertex_y << "," << m_Vertex_z << ")" << std::endl;
0244 std::cout << " Vertex width (x,y,z) = (" << m_VertexWidth_x << "," << m_VertexWidth_y << "," << m_VertexWidth_z << ")" << std::endl;
0245 }
0246 std::cout << " Vertex size function (r) = ("
0247 << m_FunctionNames.find(m_VertexSizeFunc_r)->second << ")"
0248 << std::endl;
0249 std::cout << " Vertex size (mean) = (" << m_VertexSizeMean << ")" << std::endl;
0250 std::cout << " Vertex size (width) = (" << m_VertexSizeWidth << ")" << std::endl;
0251 if (std::isfinite(m_EtaMin) && std::isfinite(m_EtaMax))
0252 {
0253 std::cout << " Eta range = " << m_EtaMin << " - " << m_EtaMax << std::endl;
0254 }
0255 if (std::isfinite(m_ThetaMin) && std::isfinite(m_ThetaMax))
0256 {
0257 std::cout << " Theta range = " << m_ThetaMin << " - " << m_ThetaMax
0258 << ", deg: " << m_ThetaMin / M_PI * 180. << " - " << m_ThetaMax / M_PI * 180. << std::endl;
0259 }
0260 std::cout << " Phi range = " << m_PhiMin << " - " << m_PhiMax
0261 << ", deg: " << m_PhiMin / M_PI * 180. << " - " << m_PhiMax / M_PI * 180. << std::endl;
0262 if (std::isfinite(m_Pt_Min) && std::isfinite(m_Pt_Max))
0263 {
0264 std::cout << " pT range = " << m_Pt_Min << " - " << m_Pt_Max << std::endl;
0265 }
0266 if (std::isfinite(m_P_Min) && std::isfinite(m_P_Max))
0267 {
0268 std::cout << " p range = " << m_P_Min << " - " << m_P_Max << std::endl;
0269 }
0270 std::cout << " t0 = " << get_t0() << std::endl;
0271 std::cout << "===========================================================================" << std::endl;
0272 }
0273
0274
0275 for (auto &_particle_code : _particle_codes)
0276 {
0277 int pdgcode = _particle_code.first;
0278 unsigned int count = _particle_code.second;
0279 std::string pdgname = get_pdgname(pdgcode);
0280 _particle_names.emplace_back(pdgname, count);
0281 }
0282
0283 return Fun4AllReturnCodes::EVENT_OK;
0284 }
0285
0286 int PHG4SimpleEventGenerator::process_event(PHCompositeNode *topNode)
0287 {
0288 if (Verbosity() > 0)
0289 {
0290 std::cout << "====================== PHG4SimpleEventGenerator::process_event() =====================" << std::endl;
0291 std::cout << "PHG4SimpleEventGenerator::process_event - reuse_existing_vertex = " << get_reuse_existing_vertex() << std::endl;
0292 }
0293
0294 if (!ReuseExistingVertex(topNode))
0295 {
0296
0297 set_vtx(smearvtx(m_Vertex_x, m_VertexWidth_x, m_VertexFunc_x),
0298 smearvtx(m_Vertex_y, m_VertexWidth_y, m_VertexFunc_y),
0299 smearvtx(m_Vertex_z, m_VertexWidth_z, m_VertexFunc_z));
0300 }
0301 set_vtx(get_vtx_x() + m_VertexOffset_x,
0302 get_vtx_y() + m_VertexOffset_y,
0303 get_vtx_z() + m_VertexOffset_z);
0304
0305 if (Verbosity() > 0)
0306 {
0307 std::cout << "PHG4SimpleEventGenerator::process_event - vertex center" << get_reuse_existing_vertex()
0308 << get_vtx_x() << ", " << get_vtx_y() << ", " << get_vtx_z() << " cm"
0309 << std::endl;
0310 }
0311
0312 int vtxindex = -1;
0313 int trackid = -1;
0314 for (unsigned int i = 0; i < _particle_names.size(); ++i)
0315 {
0316 std::string pdgname = _particle_names[i].first;
0317 int pdgcode = get_pdgcode(pdgname);
0318 unsigned int nparticles = _particle_names[i].second;
0319
0320 for (unsigned int j = 0; j < nparticles; ++j)
0321 {
0322 if ((m_VertexSizeWidth > 0.0) || (m_VertexSizeMean != 0.0))
0323 {
0324 double r = smearvtx(m_VertexSizeMean, m_VertexSizeWidth, m_VertexSizeFunc_r);
0325
0326 double x = 0.0;
0327 double y = 0.0;
0328 double z = 0.0;
0329 gsl_ran_dir_3d(RandomGenerator(), &x, &y, &z);
0330 x *= r;
0331 y *= r;
0332 z *= r;
0333
0334 vtxindex = m_InEvent->AddVtx(get_vtx_x() + x, get_vtx_y() + y, get_vtx_z() + z, get_t0());
0335 }
0336 else if ((i == 0) && (j == 0))
0337 {
0338 vtxindex = m_InEvent->AddVtx(get_vtx_x(), get_vtx_y(), get_vtx_z(), get_t0());
0339 }
0340
0341 ++trackid;
0342
0343 double eta;
0344 if (!std::isnan(m_EtaMin) && !std::isnan(m_EtaMax))
0345 {
0346 eta = (m_EtaMax - m_EtaMin) * gsl_rng_uniform_pos(RandomGenerator()) + m_EtaMin;
0347 }
0348 else if (!std::isnan(m_ThetaMin) && !std::isnan(m_ThetaMax))
0349 {
0350 double theta = (m_ThetaMax - m_ThetaMin) * gsl_rng_uniform_pos(RandomGenerator()) + m_ThetaMin;
0351 eta = PHG4Utils::get_eta(theta);
0352 }
0353 else
0354 {
0355 std::cout << PHWHERE << "Error: neither eta range or theta range was specified" << std::endl;
0356 std::cout << "That should not happen, please inform the software group howthis happened" << std::endl;
0357 exit(-1);
0358 }
0359
0360 double phi = (m_PhiMax - m_PhiMin) * gsl_rng_uniform_pos(RandomGenerator()) + m_PhiMin;
0361
0362 double pt;
0363
0364 if (!std::isnan(m_P_Min) && !std::isnan(m_P_Max) && !std::isnan(m_P_GausWidth))
0365 {
0366 pt = ((m_P_Max - m_P_Min) * gsl_rng_uniform_pos(RandomGenerator()) + m_P_Min + gsl_ran_gaussian(RandomGenerator(), m_P_GausWidth)) / cosh(eta);
0367 if (!std::isnan(m_powerLawN))
0368 {
0369 double y = gsl_rng_uniform_pos(RandomGenerator());
0370 double x1 = pow(m_Pt_Max, m_powerLawN + 1);
0371 double x0 = pow(m_Pt_Min, m_powerLawN + 1);
0372 pt = pow((x1 - x0) * y + x0, 1. / (m_powerLawN + 1.));
0373 }
0374 }
0375 else if (!std::isnan(m_Pt_Min) && !std::isnan(m_Pt_Max) && !std::isnan(m_Pt_GausWidth))
0376 {
0377 pt = (m_Pt_Max - m_Pt_Min) * gsl_rng_uniform_pos(RandomGenerator()) + m_Pt_Min + gsl_ran_gaussian(RandomGenerator(), m_Pt_GausWidth);
0378 if (!std::isnan(m_powerLawN))
0379 {
0380 double y = gsl_rng_uniform_pos(RandomGenerator());
0381 double x1 = pow(m_Pt_Max, m_powerLawN + 1);
0382 double x0 = pow(m_Pt_Min, m_powerLawN + 1);
0383 pt = pow((x1 - x0) * y + x0, 1. / (m_powerLawN + 1.));
0384 }
0385 }
0386 else
0387 {
0388 std::cout << PHWHERE << "Error: neither a p range or pt range was specified" << std::endl;
0389 exit(-1);
0390 }
0391
0392 double px = pt * cos(phi);
0393 double py = pt * sin(phi);
0394 double pz = pt * sinh(eta);
0395 double m = get_mass(pdgcode);
0396 double e = sqrt(px * px + py * py + pz * pz + m * m);
0397
0398 PHG4Particle *particle = new PHG4Particlev2();
0399 particle->set_track_id(trackid);
0400 particle->set_vtx_id(vtxindex);
0401 particle->set_parent_id(0);
0402 particle->set_name(pdgname);
0403 particle->set_pid(pdgcode);
0404 particle->set_px(px);
0405 particle->set_py(py);
0406 particle->set_pz(pz);
0407 particle->set_e(e);
0408
0409 m_InEvent->AddParticle(vtxindex, particle);
0410 if (EmbedFlag() != 0)
0411 {
0412 m_InEvent->AddEmbeddedParticle(particle, EmbedFlag());
0413 }
0414 }
0415 }
0416
0417 if (Verbosity() > 0)
0418 {
0419 m_InEvent->identify();
0420 std::cout << "======================================================================================" << std::endl;
0421 }
0422
0423 return Fun4AllReturnCodes::EVENT_OK;
0424 }
0425
0426 double
0427 PHG4SimpleEventGenerator::smearvtx(const double position, const double width, FUNCTION dist) const
0428 {
0429 double res = position;
0430 if (dist == Uniform)
0431 {
0432 res = (position - width) + 2 * gsl_rng_uniform_pos(RandomGenerator()) * width;
0433 }
0434 else if (dist == Gaus)
0435 {
0436 res = position + gsl_ran_gaussian(RandomGenerator(), width);
0437 }
0438 else
0439 {
0440 std::cout << __PRETTY_FUNCTION__ << " invalid distribution function " << dist
0441 << " (" << m_FunctionNames.find(dist)->second << ")" << std::endl;
0442 gSystem->Exit(1);
0443 }
0444 return res;
0445 }
