Thermal-FIST 1.5
Package for hadron resonance gas model applications
Loading...
Searching...
No Matches
EventGeneratorBase.cpp
Go to the documentation of this file.
1/*
2 * Thermal-FIST package
3 *
4 * Copyright (c) 2015-2019 Volodymyr Vovchenko
5 *
6 * GNU General Public License (GPLv3 or later)
7 */
8#include "HRGEventGenerator/EventGeneratorBase.h"
9
10#include <functional>
11#include <algorithm>
12#include <stdexcept>
13
14#include "HRGBase/xMath.h"
15#include "HRGBase/ThermalModelBase.h"
16#include "HRGBase/ThermalModelIdeal.h"
17#include "HRGBase/ThermalModelCanonicalStrangeness.h"
18#include "HRGBase/ThermalModelCanonical.h"
19#include "HRGEV/ThermalModelEVDiagonal.h"
20#include "HRGEV/ThermalModelEVCrossterms.h"
21#include "HRGEV/ExcludedVolumeHelper.h"
22#include "HRGVDW/ThermalModelVDW.h"
23#include "HRGRealGas/ThermalModelRealGas.h"
24#include "HRGEventGenerator/ParticleDecaysMC.h"
25#include "HRGBase/ExtraParticles.h"
26
27namespace thermalfist {
28
29 int EventGeneratorBase::fCEAccepted, EventGeneratorBase::fCETotal;
30 double EventGeneratorBase::m_LastWeight = 1., EventGeneratorBase::m_LastLogWeight = 0., EventGeneratorBase::m_LastNormWeight = 1.;
31
32 std::vector<double> LorentzBoost(const std::vector<double>& fourvector, double vx, double vy, double vz)
33 {
34 std::vector<double> ret(4, 0);
35 double v2 = vx * vx + vy * vy + vz * vz;
36 if (v2 == 0.0)
37 return fourvector;
38 double gamma = 1. / sqrt(1. - v2);
39
40 const double& r0 = fourvector[0];
41 const double& rx = fourvector[1];
42 const double& ry = fourvector[2];
43 const double& rz = fourvector[3];
44
45 ret[0] = gamma * r0 - gamma * (vx * rx + vy * ry + vz * rz);
46 ret[1] = -gamma * vx * r0 + (1. + (gamma - 1.) * vx * vx / v2) * rx +
47 (gamma - 1.) * vx * vy / v2 * ry + (gamma - 1.) * vx * vz / v2 * rz;
48 ret[2] = -gamma * vy * r0 + (1. + (gamma - 1.) * vy * vy / v2) * ry +
49 (gamma - 1.) * vy * vx / v2 * rx + (gamma - 1.) * vy * vz / v2 * rz;
50 ret[3] = -gamma * vz * r0 + (1. + (gamma - 1.) * vz * vz / v2) * rz +
51 (gamma - 1.) * vz * vx / v2 * rx + (gamma - 1.) * vz * vy / v2 * ry;
52
53 return ret;
54 }
55
56 EventGeneratorBase::EventGeneratorBase() :
57 m_THM(NULL), m_ParametersSet(false),
58 m_MeanB(0.), m_MeanAB(0.),
59 m_MeanSM(0.), m_MeanASM(0.),
60 m_MeanCM(0.), m_MeanACM(0.),
61 m_MeanCHRM(0.), m_MeanACHRM(0.),
62 m_MeanCHRMM(0.), m_MeanACHRMM(0.)
63 {
64 fCEAccepted = fCETotal = 0;
65 }
66
71
72 void EventGeneratorBase::ClearMomentumGenerators()
73 {
74 for (size_t i = 0; i < m_MomentumGens.size(); ++i) {
75 if (m_MomentumGens[i] != NULL)
76 delete m_MomentumGens[i];
77 }
78 m_MomentumGens.resize(0);
79
80 for (size_t i = 0; i < m_BWGens.size(); ++i) {
81 if (m_BWGens[i] != NULL)
82 delete m_BWGens[i];
83 }
84 m_BWGens.resize(0);
85 }
86
87 //void EventGeneratorBase::SetConfiguration(const ThermalModelParameters& params, EventGeneratorConfiguration::Ensemble ensemble, EventGeneratorConfiguration::ModelType modeltype, ThermalParticleSystem *TPS, ThermalModelBase *THMEVVDW)
88 void EventGeneratorBase::SetConfiguration(ThermalParticleSystem *TPS,
89 const EventGeneratorConfiguration& config)
90 {
91 m_Config = config;
92 SetEVUseSPR(config.fUseEVUseSPRApproximation);
93
94 if (TPS == NULL)
95 return;
96
97 if (m_Config.fEnsemble == EventGeneratorConfiguration::CCE) {
98 m_Config.CFOParameters.muC = 0.;
99 }
100 else if (m_Config.fEnsemble == EventGeneratorConfiguration::SCE) {
101 m_Config.CFOParameters.muS = 0.;
102 m_Config.CFOParameters.muC = 0.;
103 }
104 else if (0 && m_Config.fEnsemble == EventGeneratorConfiguration::CE) {
105 if (m_Config.CanonicalB)
106 m_Config.CFOParameters.muB = 0.;
107 if (m_Config.CanonicalS)
108 m_Config.CFOParameters.muS = 0.;
109 if (m_Config.CanonicalQ)
110 m_Config.CFOParameters.muQ = 0.;
111 if (m_Config.CanonicalC)
112 m_Config.CFOParameters.muC = 0.;
113 }
114
115 if (m_Config.fModelType == EventGeneratorConfiguration::PointParticle) {
116 m_THM = new ThermalModelIdeal(TPS, m_Config.CFOParameters);
117 }
118 else if (m_Config.fModelType == EventGeneratorConfiguration::DiagonalEV) {
119 m_THM = new ThermalModelEVDiagonal(TPS, m_Config.CFOParameters);
120 }
121 else if (m_Config.fModelType == EventGeneratorConfiguration::CrosstermsEV) {
122 m_THM = new ThermalModelEVCrossterms(TPS, m_Config.CFOParameters);
123 }
124 else if (m_Config.fModelType == EventGeneratorConfiguration::QvdW) {
125 m_THM = new ThermalModelVDWFull(TPS, m_Config.CFOParameters);
126 }
127 else if (m_Config.fModelType == EventGeneratorConfiguration::RealGas) {
128 m_THM = new ThermalModelRealGas(TPS, m_Config.CFOParameters);
129 }
130
131 if (m_Config.fEnsemble == EventGeneratorConfiguration::CE && m_Config.fUseGCEConservedCharges) {
132 if (!m_Config.fUsePCE)
133 m_THM->FillChemicalPotentials();
134 else
135 m_THM->SetChemicalPotentials(m_Config.fPCEChems);
136
137 m_THM->CalculatePrimordialDensities();
138
139 // Round to nearest integer
140 m_Config.B = lround(m_THM->BaryonDensity() * m_THM->Volume());
141 m_Config.Q = lround(m_THM->ElectricChargeDensity() * m_THM->Volume());
142 m_Config.S = lround(m_THM->StrangenessDensity() * m_THM->Volume());
143 m_Config.C = lround(m_THM->CharmDensity() * m_THM->Volume());
144 }
145
146 m_THM->SetUseWidth(TPS->ResonanceWidthIntegrationType());
147
148 m_THM->ConstrainMuB(false);
149 m_THM->ConstrainMuQ(false);
150 m_THM->ConstrainMuS(false);
151 m_THM->ConstrainMuC(false);
152
153 if (!m_Config.fUsePCE)
154 m_THM->FillChemicalPotentials();
155 else
156 m_THM->SetChemicalPotentials(m_Config.fPCEChems);
157
158 if (m_Config.fModelType == EventGeneratorConfiguration::DiagonalEV ||
159 m_Config.fModelType == EventGeneratorConfiguration::CrosstermsEV ||
160 m_Config.fModelType == EventGeneratorConfiguration::QvdW) {
161 for (size_t i = 0; i < m_THM->Densities().size(); ++i) {
162 for (size_t j = 0; j < m_THM->Densities().size(); ++j) {
163 if (m_Config.bij.size() == m_THM->Densities().size()
164 && m_Config.bij[i].size() == m_THM->Densities().size())
165 m_THM->SetVirial(i, j, m_Config.bij[i][j]);
166
167 if (m_Config.aij.size() == m_THM->Densities().size()
168 && m_Config.aij[i].size() == m_THM->Densities().size())
169 m_THM->SetAttraction(i, j, m_Config.aij[i][j]);
170 }
171 }
172 }
173
174 if (m_Config.fModelType == EventGeneratorConfiguration::RealGas) {
175 ExcludedVolumeModelBase* evmod;
176 if (config.RealGasExcludedVolumePrescription == 0) {
177 evmod = new ExcludedVolumeModelVDW();
178 }
179 else if (config.RealGasExcludedVolumePrescription == 1) {
180 evmod = new ExcludedVolumeModelCS();
181 }
182 else if (config.RealGasExcludedVolumePrescription == 2) {
183 evmod = new ExcludedVolumeModelVirial();
184 }
185 else if (config.RealGasExcludedVolumePrescription == 3) {
186 evmod = new ExcludedVolumeModelTVM();
187 }
188 else {
189 evmod = new ExcludedVolumeModelVDW();
190 }
191 static_cast<ThermalModelRealGas*>(m_THM)->SetExcludedVolumeModel(new ExcludedVolumeModelCrosstermsGeneralized(evmod, m_Config.bij));
192 static_cast<ThermalModelRealGas*>(m_THM)->SetMeanFieldModel(new MeanFieldModelMultiVDW(m_Config.aij));
193 }
194
195 if (m_Config.fEnsemble == EventGeneratorConfiguration::CE) {
196 // The following procedure currently not used,
197 // but can be considered if SolveChemicalPotentials routine fails
198 // Note to take into account PCE possiblity if this option is considered
199 // TODO: Properly for mixed-canonical ensembles and charm
200 if (0 && !(m_Config.B == 0 && m_Config.Q == 0 && m_Config.S == 0) && !m_Config.fUsePCE) {
201 if (m_Config.S == 0 && !(m_Config.Q == 0 || m_Config.B == 0)) {
202 double QBrat = (double)(m_Config.Q) / m_Config.B;
203
204 double left = 0.000, right = 0.900, center;
205
206 m_THM->SetBaryonChemicalPotential(left);
207 m_THM->SetQoverB(QBrat);
208 m_THM->FixParameters();
209 m_THM->CalculatePrimordialDensities();
210 double valleft = m_THM->CalculateBaryonDensity() * m_THM->Volume() - m_Config.B;
211
212 m_THM->SetBaryonChemicalPotential(right);
213 m_THM->SetQoverB(QBrat);
214 m_THM->FixParameters();
215 m_THM->CalculatePrimordialDensities();
216 double valright = m_THM->CalculateBaryonDensity() * m_THM->Volume() - m_Config.B;
217
218 double valcenter;
219
220 while ((right - left) > 0.00001) {
221 center = (left + right) / 2.;
222 m_THM->SetBaryonChemicalPotential(center);
223 m_THM->SetQoverB(QBrat);
224 m_THM->FixParameters();
225 m_THM->CalculatePrimordialDensities();
226 valcenter = m_THM->CalculateBaryonDensity() * m_THM->Volume() - m_Config.B;
227
228 if (valleft*valcenter > 0.) {
229 left = center;
230 valleft = valcenter;
231 }
232 else {
233 right = center;
234 valright = valcenter;
235 }
236 }
237
238 m_THM->SetBaryonChemicalPotential(center);
239 m_THM->SetQoverB(QBrat);
240 m_THM->FixParameters();
241
242 m_Config.CFOParameters.muB = m_THM->Parameters().muB;
243 m_Config.CFOParameters.muS = m_THM->Parameters().muS;
244 m_Config.CFOParameters.muQ = m_THM->Parameters().muQ;
245
246 m_THM->FillChemicalPotentials();
247 }
248 }
249
250 if (!m_Config.fUsePCE)
251 {
252 bool solve_chems = m_THM->SolveChemicalPotentials(m_Config.B, m_Config.Q, m_Config.S, m_Config.C,
253 m_THM->Parameters().muB, m_THM->Parameters().muQ, m_THM->Parameters().muS, m_THM->Parameters().muC,
254 m_Config.CanonicalB, m_Config.CanonicalQ, m_Config.CanonicalS, m_Config.CanonicalC);
255 if (m_THM->Parameters().muB != m_THM->Parameters().muB ||
256 m_THM->Parameters().muQ != m_THM->Parameters().muQ ||
257 m_THM->Parameters().muS != m_THM->Parameters().muS ||
258 m_THM->Parameters().muC != m_THM->Parameters().muC ||
259 !solve_chems)
260 {
261 printf("**WARNING** Could not constrain chemical potentials. Setting all for exactly conserved charges to zero...\n");
262 if (m_Config.CanonicalB)
263 m_THM->SetBaryonChemicalPotential(0.);
264 else
265 m_THM->SetBaryonChemicalPotential(m_Config.CFOParameters.muB);
266
267 if (m_Config.CanonicalQ)
268 m_THM->SetElectricChemicalPotential(0.);
269 else
270 m_THM->SetElectricChemicalPotential(m_Config.CFOParameters.muQ);
271
272 if (m_Config.CanonicalS)
273 m_THM->SetStrangenessChemicalPotential(0.);
274 else
275 m_THM->SetStrangenessChemicalPotential(m_Config.CFOParameters.muS);
276
277 if (m_Config.CanonicalC)
278 m_THM->SetCharmChemicalPotential(0.);
279 else
280 m_THM->SetCharmChemicalPotential(m_Config.CFOParameters.muC);
281
282 m_THM->FillChemicalPotentials();
283 }
284 m_Config.CFOParameters.muB = m_THM->Parameters().muB;
285 m_Config.CFOParameters.muS = m_THM->Parameters().muS;
286 m_Config.CFOParameters.muQ = m_THM->Parameters().muQ;
287 m_Config.CFOParameters.muC = m_THM->Parameters().muC;
288
289 //std::cout << "Chemical potentials constrained!\n"
290 // << "muB = " << m_Config.CFOParameters.muB
291 // << " muQ = " << m_Config.CFOParameters.muQ
292 // << " muS = " << m_Config.CFOParameters.muS
293 // << " muC = " << m_Config.CFOParameters.muC << "\n";
294
295 //std::cout << "B = " << m_THM->CalculateBaryonDensity() * m_THM->Parameters().V
296 // << " Q = " << m_THM->CalculateChargeDensity() * m_THM->Parameters().V
297 // << " S = " << m_THM->CalculateStrangenessDensity() * m_THM->Parameters().V
298 // << " C = " << m_THM->CalculateCharmDensity() * m_THM->Parameters().V
299 // << "\n";
300 }
301 }
302
303 if (m_Config.fEnsemble == EventGeneratorConfiguration::SCE && !m_Config.fUsePCE) {
304 m_THM->ConstrainMuS(true);
305 m_THM->ConstrainMuC(true);
306 m_THM->ConstrainChemicalPotentials();
307 m_THM->FillChemicalPotentials();
308 m_Config.CFOParameters.muS = m_THM->Parameters().muS;
309 m_Config.CFOParameters.muC = m_THM->Parameters().muC;
310 }
311
312 if (m_Config.fEnsemble == EventGeneratorConfiguration::CCE && !m_Config.fUsePCE) {
313 m_THM->ConstrainMuC(true);
314 m_THM->ConstrainChemicalPotentials();
315 m_THM->FillChemicalPotentials();
316 m_Config.CFOParameters.muC = m_THM->Parameters().muC;
317 }
318
319 //m_THM->CalculateDensitiesGCE();
320 m_THM->CalculatePrimordialDensities();
321 m_DensitiesIdeal = m_THM->GetIdealGasDensities();
322
323 if (m_Config.fEnsemble != EventGeneratorConfiguration::GCE)
324 PrepareMultinomials();
325
326 m_Radii = ComputeEVRadii();
327 }
328
329
330 void EventGeneratorBase::PrepareMultinomials() {
331 if (!m_THM->IsCalculated()) {
332 m_THM->CalculatePrimordialDensities();
333 }
334
335 m_Baryons.resize(0);
336 m_AntiBaryons.resize(0);
337 m_StrangeMesons.resize(0);
338 m_AntiStrangeMesons.resize(0);
339 m_ChargeMesons.resize(0);
340 m_AntiChargeMesons.resize(0);
341 m_CharmMesons.resize(0);
342 m_AntiCharmMesons.resize(0);
343 m_CharmAll.resize(0);
344 m_AntiCharmAll.resize(0);
345 m_MeanB = 0.;
346 m_MeanAB = 0.;
347 m_MeanSM = 0.;
348 m_MeanASM = 0.;
349 m_MeanCM = 0.;
350 m_MeanACM = 0.;
351 m_MeanCHRMM = 0.;
352 m_MeanACHRMM = 0.;
353 m_MeanCHRM = 0.;
354 m_MeanACHRM = 0.;
355
356 std::vector<double> yields = m_THM->Densities();
357 for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i)
358 yields[i] *= m_THM->Volume();
359
360 for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i) {
361 if (m_THM->TPS()->Particles()[i].BaryonCharge() == 1) {
362 m_Baryons.push_back(std::make_pair(yields[i], i));
363 m_MeanB += yields[i];
364 }
365 else if (m_THM->TPS()->Particles()[i].BaryonCharge() == -1) {
366 m_AntiBaryons.push_back(std::make_pair(yields[i], i));
367 m_MeanAB += yields[i];
368 }
369 else if (m_THM->TPS()->Particles()[i].BaryonCharge() == 0 && m_THM->TPS()->Particles()[i].Strangeness() == 1) {
370 m_StrangeMesons.push_back(std::make_pair(yields[i], i));
371 m_MeanSM += yields[i];
372 }
373 else if (m_THM->TPS()->Particles()[i].BaryonCharge() == 0 && m_THM->TPS()->Particles()[i].Strangeness() == -1) {
374 m_AntiStrangeMesons.push_back(std::make_pair(yields[i], i));
375 m_MeanASM += yields[i];
376 }
377 else if (m_THM->TPS()->Particles()[i].BaryonCharge() == 0 && m_THM->TPS()->Particles()[i].Strangeness() == 0 && m_THM->TPS()->Particles()[i].ElectricCharge() == 1) {
378 m_ChargeMesons.push_back(std::make_pair(yields[i], i));
379 m_MeanCM += yields[i];
380 }
381 else if (m_THM->TPS()->Particles()[i].BaryonCharge() == 0 && m_THM->TPS()->Particles()[i].Strangeness() == 0 && m_THM->TPS()->Particles()[i].ElectricCharge() == -1) {
382 m_AntiChargeMesons.push_back(std::make_pair(yields[i], i));
383 m_MeanACM += yields[i];
384 }
385 else if (m_THM->TPS()->Particles()[i].BaryonCharge() == 0 && m_THM->TPS()->Particles()[i].Strangeness() == 0 && m_THM->TPS()->Particles()[i].ElectricCharge() == 0 && m_THM->TPS()->Particles()[i].Charm() == 1) {
386 m_CharmMesons.push_back(std::make_pair(yields[i], i));
387 m_MeanCHRMM += yields[i];
388 }
389 else if (m_THM->TPS()->Particles()[i].BaryonCharge() == 0 && m_THM->TPS()->Particles()[i].Strangeness() == 0 && m_THM->TPS()->Particles()[i].ElectricCharge() == 0 && m_THM->TPS()->Particles()[i].Charm() == -1) {
390 m_AntiCharmMesons.push_back(std::make_pair(yields[i], i));
391 m_MeanACHRMM += yields[i];
392 }
393
394 if (m_THM->TPS()->Particles()[i].Charm() == 1) {
395 m_CharmAll.push_back(std::make_pair(yields[i], i));
396 m_MeanCHRM += yields[i];
397 }
398 else if (m_THM->TPS()->Particles()[i].Charm() == -1) {
399 m_AntiCharmAll.push_back(std::make_pair(yields[i], i));
400 m_MeanACHRM += yields[i];
401 }
402 }
403
404 // sort in descending order and convert to prefix sums
405 std::sort(m_Baryons.begin(), m_Baryons.end(), std::greater< std::pair<double, int> >());
406 std::sort(m_AntiBaryons.begin(), m_AntiBaryons.end(), std::greater< std::pair<double, int> >());
407 std::sort(m_StrangeMesons.begin(), m_StrangeMesons.end(), std::greater< std::pair<double, int> >());
408 std::sort(m_AntiStrangeMesons.begin(), m_AntiStrangeMesons.end(), std::greater< std::pair<double, int> >());
409 std::sort(m_ChargeMesons.begin(), m_ChargeMesons.end(), std::greater< std::pair<double, int> >());
410 std::sort(m_AntiChargeMesons.begin(), m_AntiChargeMesons.end(), std::greater< std::pair<double, int> >());
411 std::sort(m_CharmMesons.begin(), m_CharmMesons.end(), std::greater< std::pair<double, int> >());
412 std::sort(m_AntiCharmMesons.begin(), m_AntiCharmMesons.end(), std::greater< std::pair<double, int> >());
413 std::sort(m_CharmAll.begin(), m_CharmAll.end(), std::greater< std::pair<double, int> >());
414 std::sort(m_AntiCharmAll.begin(), m_AntiCharmAll.end(), std::greater< std::pair<double, int> >());
415
416 for (int i = 1; i < static_cast<int>(m_Baryons.size()); ++i) m_Baryons[i].first += m_Baryons[i - 1].first;
417 for (int i = 1; i < static_cast<int>(m_AntiBaryons.size()); ++i) m_AntiBaryons[i].first += m_AntiBaryons[i - 1].first;
418 for (int i = 1; i < static_cast<int>(m_StrangeMesons.size()); ++i) m_StrangeMesons[i].first += m_StrangeMesons[i - 1].first;
419 for (int i = 1; i < static_cast<int>(m_AntiStrangeMesons.size()); ++i) m_AntiStrangeMesons[i].first += m_AntiStrangeMesons[i - 1].first;
420 for (int i = 1; i < static_cast<int>(m_ChargeMesons.size()); ++i) m_ChargeMesons[i].first += m_ChargeMesons[i - 1].first;
421 for (int i = 1; i < static_cast<int>(m_AntiChargeMesons.size()); ++i) m_AntiChargeMesons[i].first += m_AntiChargeMesons[i - 1].first;
422 for (int i = 1; i < static_cast<int>(m_CharmMesons.size()); ++i) m_CharmMesons[i].first += m_CharmMesons[i - 1].first;
423 for (int i = 1; i < static_cast<int>(m_AntiCharmMesons.size()); ++i) m_AntiCharmMesons[i].first += m_AntiCharmMesons[i - 1].first;
424 for (int i = 1; i < static_cast<int>(m_CharmAll.size()); ++i) m_CharmAll[i].first += m_CharmAll[i - 1].first;
425 for (int i = 1; i < static_cast<int>(m_AntiCharmAll.size()); ++i) m_AntiCharmAll[i].first += m_AntiCharmAll[i - 1].first;
426 }
427
428 std::vector<int> EventGeneratorBase::GenerateTotals() const {
429 const_cast<EventGeneratorBase*>(this)->CheckSetParameters();
430
431 if (!m_THM->IsCalculated())
432 m_THM->CalculatePrimordialDensities();
433
434 std::vector<int> totals(m_THM->TPS()->Particles().size());
435
436 m_LastWeight = 1.;
437 m_LastNormWeight = 1.;
438
439 while (true) {
440 // First generate a configuration which satisfies conservation laws
441 if (m_Config.fEnsemble == EventGeneratorConfiguration::CCE)
442 totals = GenerateTotalsCCE();
443 else if (m_Config.fEnsemble == EventGeneratorConfiguration::SCE)
444 totals = GenerateTotalsSCE();
445 else if (m_Config.fEnsemble == EventGeneratorConfiguration::CE)
446 totals = GenerateTotalsCE();
447 else
448 totals = GenerateTotalsGCE();
449
450 double weight = ComputeWeightNew(totals);
451 //std::cout << ComputeWeight(totals) << " " << ComputeWeightNew(totals) << "\n";
452 if (weight < 0.)
453 continue;
454
455 break;
456 }
457
458 fCEAccepted++;
459
460 return totals;
461 }
462
463 std::vector<int> EventGeneratorBase::GenerateTotalsGCE() const
464 {
465 fCETotal++;
466
467 //if (!m_THM->IsGCECalculated())
468 // m_THM->CalculateDensitiesGCE();
469 if (!m_THM->IsCalculated())
470 m_THM->CalculatePrimordialDensities();
471 std::vector<int> totals(m_THM->TPS()->Particles().size(), 0);
472
473 const std::vector<double>& densities = m_THM->Densities();
474 for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i) {
475 double mean = densities[i] * m_THM->Volume();
476 int total = RandomGenerators::RandomPoisson(mean);
477 totals[i] = total;
478 }
479
480 return totals;
481 }
482
483 std::vector<int> EventGeneratorBase::GenerateTotalsSCE() const
484 {
485 if (!m_THM->IsCalculated())
486 m_THM->CalculatePrimordialDensities();
487
488 std::vector<int> totals(m_THM->TPS()->Particles().size(), 0);
489
490 // Generate SCE configuration depending on whether Vc > V, Vc = V, or Vc < V
491 while (true) {
492 const std::vector<double>& densities = m_THM->Densities();
493
494 // If Vc = V then just generate yields from a single ensemble
495 if (m_THM->Volume() == m_THM->CanonicalVolume()) {
496 totals = GenerateTotalsSCESubVolume(m_THM->Volume());
497 }
498 // If Vc > V then generate yields from a single ensemble with volume Vc
499 // particle from this ensemble is within a smaller volume V
500 // with probability V/Vc
501 else if (m_THM->Volume() < m_THM->CanonicalVolume()) {
502 std::vector<int> totalsaux = GenerateTotalsSCESubVolume(m_THM->CanonicalVolume());
503 double prob = m_THM->Volume() / m_THM->CanonicalVolume();
504 for (size_t i = 0; i < totalsaux.size(); ++i) {
505 for (int j = 0; j < totalsaux[i]; ++j) {
506 if (RandomGenerators::randgenMT.rand() < prob)
507 totals[i]++;
508 }
509 }
510 }
511 // If V > Vc then generate yields from (int)(V/Vc) SCE ensembles
512 // plus special treatment of one more ensemble if (V mod Vc) != 0
513 else {
514 // Number of normal SCE sub-ensembles
515 int multiples = static_cast<int>(m_THM->Volume() / m_THM->CanonicalVolume());
516
517 for (int iter = 0; iter < multiples; ++iter) {
518 std::vector<int> totalsaux = GenerateTotalsSCESubVolume(m_THM->CanonicalVolume());
519 for (size_t i = 0; i < totalsaux.size(); ++i)
520 totals[i] += totalsaux[i];
521 }
522
523 // For (V mod Vc) != 0 there is one more subvolume Vsub < Vc
524 double fraction = (m_THM->Volume() - multiples * m_THM->CanonicalVolume()) / m_THM->CanonicalVolume();
525
526 if (fraction > 0.0) {
527
528 bool successflag = false;
529
530 while (!successflag) {
531
532 std::vector<int> totalsaux = GenerateTotalsSCESubVolume(m_THM->CanonicalVolume());
533 std::vector<int> totalsaux2(m_THM->TPS()->Particles().size(), 0);
534 int netS = 0;
535 for (size_t i = 0; i < totalsaux.size(); ++i) {
536 if (m_THM->TPS()->Particles()[i].Strangeness() > 0) {
537 for (int j = 0; j < totalsaux[i]; ++j) {
538 if (RandomGenerators::randgenMT.rand() < fraction) {
539 totalsaux2[i]++;
540 netS += m_THM->TPS()->Particles()[i].Strangeness();
541 }
542 }
543 }
544 }
545
546 // Check if possible to match S+ with S-
547 // If S = 1 then must be at least one S = -1 hadron
548 if (netS == 1) {
549 bool fl = false;
550 for (size_t i = 0; i < totalsaux.size(); ++i) {
551 if (m_THM->TPS()->Particles()[i].Strangeness() < 0 && totalsaux[i] > 0) {
552 fl = true;
553 break;
554 }
555 }
556 if (!fl) {
557 printf("**WARNING** SCE Event generator: Cannot match S- with S+ = 1, discarding subconfiguration...\n");
558 continue;
559 }
560 }
561
562 int repeatmax = 10000;
563 int repeat = 0;
564 while (true) {
565 int netS2 = netS;
566 for (size_t i = 0; i < totalsaux.size(); ++i) {
567 if (m_THM->TPS()->Particles()[i].Strangeness() < 0) {
568 totalsaux2[i] = 0;
569 for (int j = 0; j < totalsaux[i]; ++j) {
570 if (RandomGenerators::randgenMT.rand() < fraction) {
571 totalsaux2[i]++;
572 netS2 += m_THM->TPS()->Particles()[i].Strangeness();
573 }
574 }
575 }
576 }
577 if (netS2 == 0) {
578 for (size_t i = 0; i < totalsaux2.size(); ++i)
579 totals[i] += totalsaux2[i];
580 successflag = true;
581 break;
582 }
583 repeat++;
584 if (repeat >= repeatmax) {
585 printf("**WARNING** SCE event generator: Cannot match S- with S+, too many tries, discarding configuration...\n");
586 break;
587 }
588 }
589 }
590 }
591 }
592
593 for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i) {
594 if (m_THM->TPS()->Particles()[i].Strangeness() == 0) {
595 double mean = densities[i] * m_THM->Volume();
596 int total = RandomGenerators::RandomPoisson(mean);
597 totals[i] = total;
598 }
599 }
600
601 return totals;
602 }
603
604 return totals;
605 }
606
607 std::vector<int> EventGeneratorBase::GenerateTotalsSCESubVolume(double VolumeSC) const
608 {
609 if (!m_THM->IsCalculated())
610 m_THM->CalculatePrimordialDensities();
611 std::vector<int> totals(m_THM->TPS()->Particles().size(), 0);
612
613 std::vector< std::pair<double, int> > fStrangeMesonsc = m_StrangeMesons;
614 std::vector< std::pair<double, int> > fAntiStrangeMesonsc = m_AntiStrangeMesons;
615
616 for (size_t i = 0; i < fStrangeMesonsc.size(); ++i)
617 fStrangeMesonsc[i].first *= VolumeSC / m_THM->Volume();
618
619 for (size_t i = 0; i < fAntiStrangeMesonsc.size(); ++i)
620 fAntiStrangeMesonsc[i].first *= VolumeSC / m_THM->Volume();
621
622 double fMeanSMc = m_MeanSM * VolumeSC / m_THM->Volume();
623 double fMeanASMc = m_MeanASM * VolumeSC / m_THM->Volume();
624
625 while (1) {
626 fCETotal++;
627 const std::vector<double>& densities = m_THM->Densities();
628
629 for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i) totals[i] = 0;
630 int netS = 0;
631 for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i) {
632 if (m_THM->TPS()->Particles()[i].BaryonCharge() != 0 && m_THM->TPS()->Particles()[i].Strangeness() != 0) {
633 double mean = densities[i] * VolumeSC;
634 int total = RandomGenerators::RandomPoisson(mean);
635 totals[i] = total;
636 netS += totals[i] * m_THM->TPS()->Particles()[i].Strangeness();
637 }
638 }
639 int tSM = RandomGenerators::RandomPoisson(fMeanSMc);
640 int tASM = RandomGenerators::RandomPoisson(fMeanASMc);
641
642 if (netS != tASM - tSM) continue;
643
644 for (int i = 0; i < tSM; ++i) {
645 std::vector< std::pair<double, int> >::iterator it = lower_bound(fStrangeMesonsc.begin(), fStrangeMesonsc.end(), std::make_pair(fMeanSMc*RandomGenerators::randgenMT.rand(), 0));
646 int tind = std::distance(fStrangeMesonsc.begin(), it);
647 if (tind < 0) tind = 0;
648 if (tind >= static_cast<int>(fStrangeMesonsc.size())) tind = fStrangeMesonsc.size() - 1;
649 totals[fStrangeMesonsc[tind].second]++;
650 }
651 for (int i = 0; i < tASM; ++i) {
652 std::vector< std::pair<double, int> >::iterator it = lower_bound(fAntiStrangeMesonsc.begin(), fAntiStrangeMesonsc.end(), std::make_pair(fMeanASMc*RandomGenerators::randgenMT.rand(), 0));
653 int tind = std::distance(fAntiStrangeMesonsc.begin(), it);
654 if (tind < 0) tind = 0;
655 if (tind >= static_cast<int>(fAntiStrangeMesonsc.size())) tind = fAntiStrangeMesonsc.size() - 1;
656 totals[fAntiStrangeMesonsc[tind].second]++;
657 }
658
659 // Cross-check that all resulting strangeness is zero
660 int finS = 0;
661 for (size_t i = 0; i < totals.size(); ++i) {
662 finS += totals[i] * m_THM->TPS()->Particles()[i].Strangeness();
663 }
664
665 if (finS != 0) {
666 throw std::runtime_error("**ERROR** EventGeneratorBase::GenerateTotalsSCESubVolume(): Generated strangeness is non-zero!");
667 }
668
669 return totals;
670 }
671 return totals;
672 }
673
674 std::vector<int> EventGeneratorBase::GenerateTotalsCCE() const
675 {
676 if (!m_THM->IsCalculated())
677 m_THM->CalculatePrimordialDensities();
678
679 // Check there are no multi-charmed particles, otherwise error
680 for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i) {
681 if (abs(m_THM->TPS()->Particles()[i].Charm()) > 1) {
682 throw std::runtime_error("CCE Event generator does not support lists with multi-charmed particles");
683 }
684 }
685
686 std::vector<int> totals(m_THM->TPS()->Particles().size(), 0);
687
688 // Generate CCE configuration depending on whether Vc > V, Vc = V, or Vc < V
689 const std::vector<double>& densities = m_THM->Densities();
690
691 // If Vc = V then just generate yields from a single ensemble
692 if (m_THM->Volume() == m_THM->CanonicalVolume()) {
693 totals = GenerateTotalsCCESubVolume(m_THM->Volume());
694 }
695 // If Vc > V then generate yields from a single ensemble with volume Vc
696 // particle from this ensemble is within a smaller volume V
697 // with probability V/Vc
698 else if (m_THM->Volume() < m_THM->CanonicalVolume()) {
699 std::vector<int> totalsaux = GenerateTotalsCCESubVolume(m_THM->CanonicalVolume());
700 double prob = m_THM->Volume() / m_THM->CanonicalVolume();
701 for (size_t i = 0; i < totalsaux.size(); ++i) {
702 for (int j = 0; j < totalsaux[i]; ++j) {
703 if (RandomGenerators::randgenMT.rand() < prob)
704 totals[i]++;
705 }
706 }
707 }
708 // If V > Vc then generate yields from (int)(V/Vc) SCE ensembles
709 // plus special treatment of one more ensemble if (V mod Vc) != 0
710 else {
711 // Number of normal CCE sub-ensembles
712 int multiples = static_cast<int>(m_THM->Volume() / m_THM->CanonicalVolume());
713
714 for (int iter = 0; iter < multiples; ++iter) {
715 std::vector<int> totalsaux = GenerateTotalsCCESubVolume(m_THM->CanonicalVolume());
716 for (size_t i = 0; i < totalsaux.size(); ++i)
717 totals[i] += totalsaux[i];
718 }
719
720 // For (V mod Vc) != 0 there is one more subvolume Vsub < Vc
721 double fraction = (m_THM->Volume() - multiples * m_THM->CanonicalVolume()) / m_THM->CanonicalVolume();
722
723 if (fraction > 0.0) {
724
725 bool successflag = false;
726
727 while (!successflag) {
728
729 std::vector<int> totalsaux = GenerateTotalsCCESubVolume(m_THM->CanonicalVolume());
730 std::vector<int> totalsaux2(m_THM->TPS()->Particles().size(), 0);
731 int netC = 0;
732 for (size_t i = 0; i < totalsaux.size(); ++i) {
733 if (m_THM->TPS()->Particles()[i].Charm() > 0) {
734 for (int j = 0; j < totalsaux[i]; ++j) {
735 if (RandomGenerators::randgenMT.rand() < fraction) {
736 totalsaux2[i]++;
737 netC += m_THM->TPS()->Particles()[i].Charm();
738 }
739 }
740 }
741 }
742
743 // Check if possible to match C+ with C-
744 // If C = 1 then must be at least one C = -1 hadron
745 if (netC == 1) {
746 bool fl = false;
747 for (size_t i = 0; i < totalsaux.size(); ++i) {
748 if (m_THM->TPS()->Particles()[i].Charm() < 0 && totalsaux[i] > 0) {
749 fl = true;
750 break;
751 }
752 }
753 if (!fl) {
754 printf("**WARNING** CCE Event generator: Cannot match C- with C+ = 1, discarding...\n");
755 continue;
756 }
757 }
758
759 int repeatmax = 10000;
760 int repeat = 0;
761 while (true) {
762 int netC2 = netC;
763 for (size_t i = 0; i < totalsaux.size(); ++i) {
764 if (m_THM->TPS()->Particles()[i].Charm() < 0) {
765 totalsaux2[i] = 0;
766 for (int j = 0; j < totalsaux[i]; ++j) {
767 if (RandomGenerators::randgenMT.rand() < fraction) {
768 totalsaux2[i]++;
769 netC2 += m_THM->TPS()->Particles()[i].Charm();
770 }
771 }
772 }
773 }
774 if (netC2 == 0) {
775 for (size_t i = 0; i < totalsaux2.size(); ++i)
776 totals[i] += totalsaux2[i];
777 successflag = true;
778 break;
779 }
780 repeat++;
781 if (repeat >= repeatmax) {
782 printf("**WARNING** CCE event generator: Cannot match S- with S+, too many tries, discarding...\n");
783 break;
784 }
785 }
786 }
787 }
788 }
789
790 for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i) {
791 if (m_THM->TPS()->Particles()[i].Charm() == 0) {
792 double mean = densities[i] * m_THM->Volume();
793 int total = RandomGenerators::RandomPoisson(mean);
794 totals[i] = total;
795 }
796 }
797
798 return totals;
799 }
800
801 std::vector<int> EventGeneratorBase::GenerateTotalsCCESubVolume(double VolumeSC) const
802 {
803 if (!m_THM->IsCalculated())
804 m_THM->CalculatePrimordialDensities();
805
806 std::vector<int> totals(m_THM->TPS()->Particles().size(), 0);
807
808 std::vector< std::pair<double, int> > fCharmAllc = m_CharmAll;
809 std::vector< std::pair<double, int> > fAntiCharmAllc = m_AntiCharmAll;
810
811 for (size_t i = 0; i < fCharmAllc.size(); ++i)
812 fCharmAllc[i].first *= VolumeSC / m_THM->Volume();
813
814 for (size_t i = 0; i < fAntiCharmAllc.size(); ++i)
815 fAntiCharmAllc[i].first *= VolumeSC / m_THM->Volume();
816
817 // Assuming no multi-charmed particles
818 double fMeanCharmc = m_MeanCHRM * VolumeSC / m_THM->Volume();
819 double fMeanAntiCharmc = m_MeanACHRM * VolumeSC / m_THM->Volume();
820
821 fCETotal++;
822
823 int netC = 0;
824 int tC = RandomGenerators::RandomPoisson(fMeanCharmc);
825 int tAC = RandomGenerators::RandomPoisson(fMeanAntiCharmc);
826 while (tC - tAC != m_Config.C - netC) {
827 fCETotal++;
828 tC = RandomGenerators::RandomPoisson(fMeanCharmc);
829 tAC = RandomGenerators::RandomPoisson(fMeanAntiCharmc);
830 }
831
832 for (int i = 0; i < tC; ++i) {
833 std::vector< std::pair<double, int> >::iterator it = lower_bound(fCharmAllc.begin(), fCharmAllc.end(), std::make_pair(fMeanCharmc*RandomGenerators::randgenMT.rand(), 0));
834 int tind = std::distance(fCharmAllc.begin(), it);
835 if (tind < 0) tind = 0;
836 if (tind >= static_cast<int>(fCharmAllc.size())) tind = fCharmAllc.size() - 1;
837 totals[fCharmAllc[tind].second]++;
838 }
839 for (int i = 0; i < tAC; ++i) {
840 std::vector< std::pair<double, int> >::iterator it = lower_bound(fAntiCharmAllc.begin(), fAntiCharmAllc.end(), std::make_pair(fMeanAntiCharmc*RandomGenerators::randgenMT.rand(), 0));
841 int tind = std::distance(fAntiCharmAllc.begin(), it);
842 if (tind < 0) tind = 0;
843 if (tind >= static_cast<int>(fAntiCharmAllc.size())) tind = fAntiCharmAllc.size() - 1;
844 totals[fAntiCharmAllc[tind].second]++;
845 }
846
847 // Cross-check that total resulting net charm is zero
848 int finC = 0;
849 for (size_t i = 0; i < totals.size(); ++i) {
850 finC += totals[i] * m_THM->TPS()->Particles()[i].Charm();
851 }
852
853 if (finC != m_Config.C) {
854 throw std::runtime_error("**ERROR** EventGeneratorBase::GenerateTotalsCCESubVolume(): Generated charm is non-zero!");
855 }
856
857 return totals;
858 }
859
865
866 std::vector<std::vector<double>> EventGeneratorBase::ComputeEVRadii() const
867 {
868 int Nspecies = m_THM->TPS()->ComponentsNumber();
869 std::vector<std::vector<double>> radii = std::vector<std::vector<double>>(Nspecies,
870 std::vector<double>(Nspecies, 0.0));
871
872 for (int i = 0; i < Nspecies; ++i) {
873 for (int j = 0; j < Nspecies; ++j) {
874 double b = 0.0;
875 if (i < m_Config.bij.size() && j < m_Config.bij[i].size())
876 b = 0.5 * (m_Config.bij[i][j] + m_Config.bij[j][i]);
877 radii[i][j] = CuteHRGHelper::rv(b);
878 }
879 }
880
881 return radii;
882 }
883
884 bool EventGeneratorBase::CheckEVOverlap(
885 const std::vector<SimpleParticle>& particles,
886 const SimpleParticle& cand,
887 const std::vector<int>& ids,
888 const std::vector<std::vector<double>>& radii)
889 const
890 {
891 if (m_Config.fModelType != EventGeneratorConfiguration::DiagonalEV
892 && m_Config.fModelType != EventGeneratorConfiguration::CrosstermsEV
893 && m_Config.fModelType != EventGeneratorConfiguration::QvdW)
894 return false;
895 int idcand = m_THM->TPS()->PdgToId(cand.PDGID);
896 for (int ipart = 0; ipart < particles.size(); ++ipart) {
897 const SimpleParticle& part = particles[ipart];
898 int idpart = 0;
899 if (ids.size())
900 idpart = ids[ipart];
901 else
902 idpart = m_THM->TPS()->PdgToId(part.PDGID);
903
904 double r = 0.0;
905 if (radii.size()) {
906 r = radii[idpart][idcand];
907 }
908 else {
909 double b = 0.5 * (m_Config.bij[idpart][idcand] + m_Config.bij[idcand][idpart]);
910 r = CuteHRGHelper::rv(b);
911 }
912
913 if (r == 0.0)
914 continue;
915
916 double dist2 = ParticleDecaysMC::ParticleDistanceSquared(part, cand);
917 if (dist2 <= 4. * r * r)
918 return true;
919 }
920 return false;
921 }
922
923
924 std::vector<int> EventGeneratorBase::GenerateTotalsCE() const {
925 if (!m_THM->IsCalculated())
926 m_THM->CalculatePrimordialDensities();
927 std::vector<int> totals(m_THM->TPS()->Particles().size(), 0);
928
929 const std::vector< std::pair<double, int> >& fBaryonsc = m_Baryons;
930 const std::vector< std::pair<double, int> >& fAntiBaryonsc = m_AntiBaryons;
931 const std::vector< std::pair<double, int> >& fStrangeMesonsc = m_StrangeMesons;
932 const std::vector< std::pair<double, int> >& fAntiStrangeMesonsc = m_AntiStrangeMesons;
933 const std::vector< std::pair<double, int> >& fChargeMesonsc = m_ChargeMesons;
934 const std::vector< std::pair<double, int> >& fAntiChargeMesonsc = m_AntiChargeMesons;
935 const std::vector< std::pair<double, int> >& fCharmMesonsc = m_CharmMesons;
936 const std::vector< std::pair<double, int> >& fAntiCharmMesonsc = m_AntiCharmMesons;
937
938 // Primitive rejection sampling (not used, but can be explored for comparisons)
939 while (0) {
940 fCETotal++;
941 int netB = 0, netS = 0, netQ = 0, netC = 0;
942 for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i) {
943 double mean = m_THM->Densities()[i] * m_THM->Volume();
944 int total = RandomGenerators::RandomPoisson(mean);
945 totals[i] = total;
946 netB += totals[i] * m_THM->TPS()->Particles()[i].BaryonCharge();
947 netS += totals[i] * m_THM->TPS()->Particles()[i].Strangeness();
948 netQ += totals[i] * m_THM->TPS()->Particles()[i].ElectricCharge();
949 netC += totals[i] * m_THM->TPS()->Particles()[i].Charm();
950 }
951 //if ((!m_Config.CanonicalB || netB == m_THM->Parameters().B)
952 // && (!m_Config.CanonicalS || netS == m_THM->Parameters().S)
953 // && (!m_Config.CanonicalQ || netQ == m_THM->Parameters().Q)
954 // && (!m_Config.CanonicalC || netC == m_THM->Parameters().C)) {
955 // fCEAccepted++;
956 // return totals;
957 //}
958 if ((!m_Config.CanonicalB || netB == m_Config.B)
959 && (!m_Config.CanonicalS || netS == m_Config.S)
960 && (!m_Config.CanonicalQ || netQ == m_Config.Q)
961 && (!m_Config.CanonicalC || netC == m_Config.C)) {
962 fCEAccepted++;
963 return totals;
964 }
965 }
966
967 // Multi-step procedure as described in F. Becattini, L. Ferroni, hep-ph/0307061
968 while (1) {
969 fCETotal++;
970
971 const std::vector<double>& densities = m_THM->Densities();
972
973 for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i) totals[i] = 0;
974 int netB = 0, netS = 0, netQ = 0, netC = 0;
975
976
977 // Light nuclei first
978 bool flNuclei = false; // Whether light nuclei appear at all
979
980 for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i) {
981 if (abs(m_THM->TPS()->Particles()[i].BaryonCharge()) > 1) {
982 flNuclei = true;
983 double mean = densities[i] * m_THM->Volume();
984 int total = RandomGenerators::RandomPoisson(mean);
985 totals[i] = total;
986 netB += totals[i] * m_THM->TPS()->Particles()[i].BaryonCharge();
987 netS += totals[i] * m_THM->TPS()->Particles()[i].Strangeness();
988 netQ += totals[i] * m_THM->TPS()->Particles()[i].ElectricCharge();
989 netC += totals[i] * m_THM->TPS()->Particles()[i].Charm();
990 }
991 }
992
993 // Then all hadrons
994
995 int tB = 0, tAB = 0;
996 // First total baryons and antibaryons from the Poisson distribution
997 if (flNuclei || !m_Config.CanonicalB) {
998 tB = RandomGenerators::RandomPoisson(m_MeanB);
999 tAB = RandomGenerators::RandomPoisson(m_MeanAB);
1000 if (m_Config.CanonicalB && tB - tAB != m_Config.B - netB) continue;
1001 //if (RandomGenerators::randgenMT.rand() > RandomGenerators::SkellamProbability(m_THM->Parameters().B - netB, m_MeanB, m_MeanAB))
1002 // continue;
1003 }
1004 else
1005 // Generate from the Bessel distribution, using Devroye's method, if no light nuclei
1006 {
1007 int nu = m_Config.B - netB;
1008 if (nu < 0) nu = -nu;
1009 double a = 2. * sqrt(m_MeanB * m_MeanAB);
1010 //int BessN = RandomGenerators::BesselDistributionGenerator::RandomBesselDevroye3(a, nu);
1011 //int BessN = RandomGenerators::BesselDistributionGenerator::RandomBesselPoisson(a, nu);
1012 //int BessN = RandomGenerators::BesselDistributionGenerator::RandomBesselCombined(a, nu);
1013 int BessN = RandomGenerators::BesselDistributionGenerator::RandomBesselDevroye1(a, nu);
1014 if (m_Config.B - netB < 0) {
1015 tB = BessN;
1016 tAB = nu + tB;
1017 }
1018 else {
1019 tAB = BessN;
1020 tB = nu + tAB;
1021 }
1022 }
1023
1024 // Then individual baryons and antibaryons from the multinomial distribution
1025 for (int i = 0; i < tB; ++i) {
1026 std::vector< std::pair<double, int> >::const_iterator it = lower_bound(fBaryonsc.begin(), fBaryonsc.end(), std::make_pair(m_MeanB*RandomGenerators::randgenMT.rand(), 0));
1027 int tind = std::distance(fBaryonsc.begin(), it);
1028 if (tind < 0) tind = 0;
1029 if (tind >= static_cast<int>(fBaryonsc.size())) tind = fBaryonsc.size() - 1;
1030 totals[fBaryonsc[tind].second]++;
1031 netS += m_THM->TPS()->Particles()[fBaryonsc[tind].second].Strangeness();
1032 netQ += m_THM->TPS()->Particles()[fBaryonsc[tind].second].ElectricCharge();
1033 netC += m_THM->TPS()->Particles()[fBaryonsc[tind].second].Charm();
1034 }
1035 for (int i = 0; i < tAB; ++i) {
1036 std::vector< std::pair<double, int> >::const_iterator it = lower_bound(fAntiBaryonsc.begin(), fAntiBaryonsc.end(), std::make_pair(m_MeanAB*RandomGenerators::randgenMT.rand(), 0));
1037 int tind = std::distance(fAntiBaryonsc.begin(), it);
1038 if (tind < 0) tind = 0;
1039 if (tind >= static_cast<int>(fAntiBaryonsc.size())) tind = fAntiBaryonsc.size() - 1;
1040 totals[fAntiBaryonsc[tind].second]++;
1041 netS += m_THM->TPS()->Particles()[fAntiBaryonsc[tind].second].Strangeness();
1042 netQ += m_THM->TPS()->Particles()[fAntiBaryonsc[tind].second].ElectricCharge();
1043 netC += m_THM->TPS()->Particles()[fAntiBaryonsc[tind].second].Charm();
1044 }
1045
1046 // Total numbers of (anti)strange mesons
1047
1048 int tSM = RandomGenerators::RandomPoisson(m_MeanSM);
1049 int tASM = RandomGenerators::RandomPoisson(m_MeanASM);
1050 if (m_Config.CanonicalS && netS != tASM - tSM + m_Config.S) continue;
1051
1052
1053 // Multinomial distribution for individual numbers of (anti)strange mesons
1054 for (int i = 0; i < tSM; ++i) {
1055 std::vector< std::pair<double, int> >::const_iterator it = lower_bound(fStrangeMesonsc.begin(), fStrangeMesonsc.end(), std::make_pair(m_MeanSM*RandomGenerators::randgenMT.rand(), 0));
1056 int tind = std::distance(fStrangeMesonsc.begin(), it);
1057 if (tind < 0) tind = 0;
1058 if (tind >= static_cast<int>(fStrangeMesonsc.size())) tind = fStrangeMesonsc.size() - 1;
1059 totals[fStrangeMesonsc[tind].second]++;
1060 netQ += m_THM->TPS()->Particles()[fStrangeMesonsc[tind].second].ElectricCharge();
1061 netC += m_THM->TPS()->Particles()[fStrangeMesonsc[tind].second].Charm();
1062 }
1063 for (int i = 0; i < tASM; ++i) {
1064 std::vector< std::pair<double, int> >::const_iterator it = lower_bound(fAntiStrangeMesonsc.begin(), fAntiStrangeMesonsc.end(), std::make_pair(m_MeanASM*RandomGenerators::randgenMT.rand(), 0));
1065 int tind = std::distance(fAntiStrangeMesonsc.begin(), it);
1066 if (tind < 0) tind = 0;
1067 if (tind >= static_cast<int>(fAntiStrangeMesonsc.size())) tind = fAntiStrangeMesonsc.size() - 1;
1068 totals[fAntiStrangeMesonsc[tind].second]++;
1069 netQ += m_THM->TPS()->Particles()[fAntiStrangeMesonsc[tind].second].ElectricCharge();
1070 netC += m_THM->TPS()->Particles()[fAntiStrangeMesonsc[tind].second].Charm();
1071 }
1072
1073 // Total numbers of remaining electrically charged mesons
1074 int tCM = RandomGenerators::RandomPoisson(m_MeanCM);
1075 int tACM = RandomGenerators::RandomPoisson(m_MeanACM);
1076 if (m_Config.CanonicalQ && netQ != tACM - tCM + m_Config.Q) continue;
1077
1078 // Multinomial distribution for individual numbers of remaining electrically charged mesons
1079 for (int i = 0; i < tCM; ++i) {
1080 std::vector< std::pair<double, int> >::const_iterator it = lower_bound(fChargeMesonsc.begin(), fChargeMesonsc.end(), std::make_pair(m_MeanCM*RandomGenerators::randgenMT.rand(), 0));
1081 int tind = std::distance(fChargeMesonsc.begin(), it);
1082 if (tind < 0) tind = 0;
1083 if (tind >= static_cast<int>(fChargeMesonsc.size())) tind = fChargeMesonsc.size() - 1;
1084 totals[fChargeMesonsc[tind].second]++;
1085 netC += m_THM->TPS()->Particles()[fChargeMesonsc[tind].second].Charm();
1086 }
1087 for (int i = 0; i < tACM; ++i) {
1088 std::vector< std::pair<double, int> >::const_iterator it = lower_bound(fAntiChargeMesonsc.begin(), fAntiChargeMesonsc.end(), std::make_pair(m_MeanACM*RandomGenerators::randgenMT.rand(), 0));
1089 int tind = std::distance(fAntiChargeMesonsc.begin(), it);
1090 if (tind < 0) tind = 0;
1091 if (tind >= static_cast<int>(fAntiChargeMesonsc.size())) tind = fAntiChargeMesonsc.size() - 1;
1092 totals[fAntiChargeMesonsc[tind].second]++;
1093 netC += m_THM->TPS()->Particles()[fAntiChargeMesonsc[tind].second].Charm();
1094 }
1095
1096 // Total numbers of remaining charmed mesons
1097 int tCHRMM = RandomGenerators::RandomPoisson(m_MeanCHRMM);
1098 int tACHRNMM = RandomGenerators::RandomPoisson(m_MeanACHRMM);
1099
1100 if (m_Config.CanonicalC && netC != tACHRNMM - tCHRMM + m_Config.C) continue;
1101
1102 // Multinomial distribution for individual numbers of the remaining charmed mesons
1103 for (int i = 0; i < tCHRMM; ++i) {
1104 std::vector< std::pair<double, int> >::const_iterator it = lower_bound(fCharmMesonsc.begin(), fCharmMesonsc.end(), std::make_pair(m_MeanCHRMM*RandomGenerators::randgenMT.rand(), 0));
1105 int tind = std::distance(fCharmMesonsc.begin(), it);
1106 if (tind < 0) tind = 0;
1107 if (tind >= static_cast<int>(fCharmMesonsc.size())) tind = fCharmMesonsc.size() - 1;
1108 totals[fCharmMesonsc[tind].second]++;
1109 }
1110 for (int i = 0; i < tACHRNMM; ++i) {
1111 std::vector< std::pair<double, int> >::const_iterator it = lower_bound(fAntiCharmMesonsc.begin(), fAntiCharmMesonsc.end(), std::make_pair(m_MeanACHRMM*RandomGenerators::randgenMT.rand(), 0));
1112 int tind = std::distance(fAntiCharmMesonsc.begin(), it);
1113 if (tind < 0) tind = 0;
1114 if (tind >= static_cast<int>(fAntiCharmMesonsc.size())) tind = fAntiCharmMesonsc.size() - 1;
1115 totals[fAntiCharmMesonsc[tind].second]++;
1116 }
1117
1118 // Poisson distribution for all neutral particles
1119 for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i) {
1120 if (m_THM->TPS()->Particles()[i].BaryonCharge() == 0
1121 && m_THM->TPS()->Particles()[i].Strangeness() == 0
1122 && m_THM->TPS()->Particles()[i].ElectricCharge() == 0
1123 && m_THM->TPS()->Particles()[i].Charm() == 0) {
1124 double mean = densities[i] * m_THM->Volume();
1125 int total = RandomGenerators::RandomPoisson(mean);
1126 totals[i] = total;
1127 }
1128 }
1129
1130 // Cross-check that all resulting charges are OK
1131 int finB = 0, finQ = 0, finS = 0, finC = 0;
1132 for (size_t i = 0; i < totals.size(); ++i) {
1133 finB += totals[i] * m_THM->TPS()->Particles()[i].BaryonCharge();
1134 finQ += totals[i] * m_THM->TPS()->Particles()[i].ElectricCharge();
1135 finS += totals[i] * m_THM->TPS()->Particles()[i].Strangeness();
1136 finC += totals[i] * m_THM->TPS()->Particles()[i].Charm();
1137 }
1138
1139 if (m_Config.CanonicalB && finB != m_Config.B) {
1140 throw std::runtime_error("**ERROR** EventGeneratorBase::GenerateTotalsCE(): Generated baryon number does not match the input!");
1141 }
1142
1143 if (m_Config.CanonicalQ && finQ != m_Config.Q) {
1144 throw std::runtime_error("**ERROR** EventGeneratorBase::GenerateTotalsCE(): Generated electric charge does not match the input!");
1145 }
1146
1147 if (m_Config.CanonicalS && finS != m_Config.S) {
1148 throw std::runtime_error("**ERROR** EventGeneratorBase::GenerateTotalsCE(): Generated strangeness does not match the input!");
1149 }
1150
1151 if (m_Config.CanonicalC && finC != m_Config.C) {
1152 throw std::runtime_error("**ERROR** EventGeneratorBase::GenerateTotalsCE(): Generated charm does not match the input!");
1153 }
1154
1155 return totals;
1156 }
1157 return totals;
1158 }
1159
1160 std::pair<std::vector<int>, double> EventGeneratorBase::SampleYields() const
1161 {
1162 std::vector<int> totals = GenerateTotals();
1163 return make_pair(totals, m_LastNormWeight);
1164 //std::vector<int> totals = GenerateTotalsGCE();
1165 //return make_pair(totals, 1.);
1166 }
1167
1168 SimpleParticle EventGeneratorBase::SampleParticle(int id) const
1169 {
1170 const_cast<EventGeneratorBase*>(this)->CheckSetParameters();
1171 const ThermalParticle& species = m_THM->TPS()->Particles()[id];
1172 double tmass = species.Mass();
1173 if (m_THM->UseWidth() && !species.ZeroWidthEnforced() && !(species.GetResonanceWidthIntegrationType() == ThermalParticle::ZeroWidth))
1174 tmass = m_BWGens[id]->GetRandom();
1175
1176 // Check for Bose-Einstein condensation
1177 // Force m = mu if the sampled mass is too small
1178 double tmu = m_THM->FullIdealChemicalPotential(id);
1179 if (species.Statistics() == -1 && tmu > tmass) {
1180 tmass = tmu;
1181 }
1182
1183 std::vector<double> momentum = m_MomentumGens[id]->GetMomentum(tmass);
1184
1185 return SimpleParticle(momentum[0], momentum[1], momentum[2], tmass, species.PdgId(), 0,
1186 momentum[3], momentum[4], momentum[5], momentum[6]);
1187 }
1188
1189 SimpleParticle EventGeneratorBase::SampleParticleByPdg(long long pdgid) const
1190 {
1191 int id = m_THM->TPS()->PdgToId(pdgid);
1192 if (id == -1) {
1193 throw std::invalid_argument("EventGeneratorBase::SampleParticleByPdg(): The input pdg code does not exist in the particle list!");
1194 }
1195 return SampleParticle(id);
1196 }
1197
1198 SimpleEvent EventGeneratorBase::SampleMomenta(const std::vector<int>& yields) const
1199 {
1200 return SampleMomentaWithShuffle(yields);
1201 //const_cast<EventGeneratorBase*>(this)->CheckSetParameters();
1202
1203 //SimpleEvent ret;
1204
1205 //std::vector<int> ids;
1206
1207 //for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i) {
1208 // const ThermalParticle& species = m_THM->TPS()->Particles()[i];
1209 // //primParticles[i].resize(0);
1210 // int total = yields[i];
1211 // for (int part = 0; part < total; ++part) {
1212 // SimpleParticle cand = SampleParticle(i);
1213 // if (!CheckEVOverlap(ret.Particles, cand, ids, m_Radii)) {
1214 // ret.Particles.push_back(cand);
1215 // ids.push_back(i);
1216 // }
1217 // else {
1218 // part--;
1219 // }
1220 // }
1221 //}
1222
1223 //ret.AllParticles = ret.Particles;
1224
1225 //ret.DecayMap.resize(ret.Particles.size());
1226 //fill(ret.DecayMap.begin(), ret.DecayMap.end(), -1);
1227
1228 //ret.DecayMapFinal.resize(ret.Particles.size());
1229 //for (int i = 0; i < ret.DecayMapFinal.size(); ++i)
1230 // ret.DecayMapFinal[i] = i;
1231
1232 //return ret;
1233 }
1234
1235 SimpleEvent EventGeneratorBase::SampleMomentaWithShuffle(const std::vector<int>& yields) const
1236 {
1237 const_cast<EventGeneratorBase*>(this)->CheckSetParameters();
1238
1239 SimpleEvent ret;
1240
1241 std::vector<int> ids;
1242 for (int i = 0; i < m_THM->TPS()->Particles().size(); ++i)
1243 for (int part = 0; part < yields[i]; ++part)
1244 ids.push_back(i);
1245 //std::random_shuffle(ids.begin(), ids.end()); // Removed in C++17
1246 std::shuffle(ids.begin(), ids.end(), RandomGenerators::rng_std);
1247
1248 ret.Particles.resize(ids.size());
1249
1250
1251 bool flOverlap = true;
1252 while (flOverlap) {
1253 int sampled = 0;
1254 while (sampled < ids.size()) {
1255 flOverlap = false;
1256 int i = ids[sampled];
1257 const ThermalParticle& species = m_THM->TPS()->Particles()[i];
1258 SimpleParticle cand = SampleParticle(i);
1259
1260 if (m_Config.fUseEVRejectionCoordinates &&
1261 (m_Config.fModelType == EventGeneratorConfiguration::DiagonalEV
1262 || m_Config.fModelType == EventGeneratorConfiguration::CrosstermsEV
1263 || m_Config.fModelType == EventGeneratorConfiguration::QvdW)) {
1264 for (int i = 0; i < sampled; ++i) {
1265 double r = m_Radii[ids[i]][ids[sampled]];
1266 if (r != 0.0) {
1267 flOverlap |= (ParticleDecaysMC::ParticleDistanceSquared(ret.Particles[i], cand) <= 4. * r * r);
1268 }
1269 }
1270
1271 if (flOverlap) {
1272 if (EVUseSPR()) {
1273 continue;
1274 }
1275 else {
1276 break;
1277 }
1278 }
1279 }
1280
1281 ret.Particles[sampled] = cand;
1282 sampled++;
1283 }
1284 }
1285
1286 ret.AllParticles = ret.Particles;
1287
1288 ret.DecayMap.resize(ret.Particles.size());
1289 fill(ret.DecayMap.begin(), ret.DecayMap.end(), -1);
1290
1291 ret.DecayMapFinal.resize(ret.Particles.size());
1292 for (int i = 0; i < ret.DecayMapFinal.size(); ++i)
1293 ret.DecayMapFinal[i] = i;
1294
1295 return ret;
1296 }
1297
1298 SimpleEvent EventGeneratorBase::GetEvent(bool DoDecays) const
1299 {
1300 const_cast<EventGeneratorBase*>(this)->CheckSetParameters();
1301
1302 if (!m_THM->IsCalculated())
1303 m_THM->CalculatePrimordialDensities();
1304
1305 std::vector<int> totals = GenerateTotals();
1306 if ((m_Config.fModelType == EventGeneratorConfiguration::DiagonalEV
1307 || m_Config.fModelType == EventGeneratorConfiguration::CrosstermsEV)
1308 && m_Config.fUseEVRejectionMultiplicity) {
1309 while (RandomGenerators::randgenMT.rand() > m_LastNormWeight) {
1310 if (m_LastNormWeight > 1.) {
1311 printf("**WARNING** Event weight %lf > 1 in Monte Carlo rejection sampling!", m_LastNormWeight);
1312 }
1313 totals = GenerateTotals();
1314 }
1315 m_LastWeight = 1.0;
1316 m_LastLogWeight = 0.;
1317 m_LastNormWeight = 1.0;
1318 }
1319
1320 SimpleEvent ret = SampleMomenta(totals);
1321 ret.weight = m_LastWeight;
1322 ret.logweight = m_LastLogWeight;
1323 ret.weight = m_LastNormWeight;
1324
1325 if (DoDecays)
1326 return PerformDecays(ret, m_THM->TPS());
1327 else
1328 return ret;
1329 }
1330
1331 // SimpleEvent EventGeneratorBase::PerformDecaysAlternativeWay(const SimpleEvent& evtin, ThermalParticleSystem* TPS)
1332 // {
1333 // SimpleEvent ret;
1334 // ret.weight = evtin.weight;
1335 // ret.logweight = evtin.logweight;
1336 // ret.AllParticles = evtin.Particles;
1337 // reverse(ret.AllParticles.begin(), ret.AllParticles.end());
1338
1339 // for (auto&& part : ret.AllParticles)
1340 // part.processed = false;
1341
1342
1343 // bool flag_repeat = true;
1344 // while (flag_repeat) {
1345 // flag_repeat = false;
1346
1347 // bool current_stable_flag = false;
1348 // long long current_pdgcode = 0;
1349 // long long current_tid = -1;
1350 // for (int i = ret.AllParticles.size() - 1; i >= 0; --i) {
1351 // SimpleParticle& particle = ret.AllParticles[i];
1352
1353 // if (particle.processed)
1354 // continue;
1355
1356 // long long tpdgcode = particle.PDGID;
1357 // if (!(tpdgcode == current_pdgcode))
1358 // {
1359 // current_tid = TPS->PdgToId(tpdgcode);
1360 // if (current_tid != -1)
1361 // current_stable_flag = TPS->Particle(current_tid).IsStable();
1362 // else
1363 // current_stable_flag = true;
1364 // current_pdgcode = tpdgcode;
1365 // }
1366
1367
1368 // if (current_stable_flag) {
1369 // //SimpleParticle prt = primParticles[i][j];
1370 // //double tpt = prt.GetPt();
1371 // //double ty = prt.GetY();
1372 // //if (static_cast<int>(m_acc.size()) < i || !m_acc[i].init || m_acc[i].getAcceptance(ty + m_ycm, tpt) > RandomGenerators::randgenMT.rand())
1373 // // ret.Particles.push_back(prt);
1374 // //primParticles[i][j].processed = true;
1375 // ret.Particles.push_back(particle);
1376 // particle.processed = true;
1377 // }
1378 // else {
1379 // flag_repeat = true;
1380 // double DecParam = RandomGenerators::randgenMT.rand(), tsum = 0.;
1381
1382 // std::vector<double> Bratios;
1383 // if (particle.MotherPDGID != 0 ||
1384 // TPS->ResonanceWidthIntegrationType() != ThermalParticle::eBW) {
1385 // Bratios = TPS->Particles()[current_tid].BranchingRatiosM(particle.m, false);
1386 // }
1387 // else {
1388 // Bratios = TPS->Particles()[current_tid].BranchingRatiosM(particle.m, true);
1389 // }
1390
1391 // int DecayIndex = 0;
1392 // for (DecayIndex = 0; DecayIndex < static_cast<int>(Bratios.size()); ++DecayIndex) {
1393 // tsum += Bratios[DecayIndex];
1394 // if (tsum > DecParam) break;
1395 // }
1396 // if (DecayIndex < static_cast<int>(TPS->Particles()[current_tid].Decays().size())) {
1397 // std::vector<double> masses(0);
1398 // std::vector<long long> pdgids(0);
1399 // for (size_t di = 0; di < TPS->Particles()[current_tid].Decays()[DecayIndex].mDaughters.size(); di++) {
1400 // long long dpdg = TPS->Particles()[current_tid].Decays()[DecayIndex].mDaughters[di];
1401 // if (TPS->PdgToId(dpdg) == -1) {
1402 // continue;
1403 // }
1404 // masses.push_back(TPS->ParticleByPDG(dpdg).Mass());
1405 // pdgids.push_back(dpdg);
1406 // }
1407 // std::vector<SimpleParticle> decres = ParticleDecaysMC::ManyBodyDecay(particle, masses, pdgids);
1408 // for (size_t ind = 0; ind < decres.size(); ind++) {
1409 // decres[ind].processed = false;
1410 // if (TPS->PdgToId(decres[ind].PDGID) != -1)
1411 // ret.AllParticles.push_back(decres[ind]);
1412 // }
1413 // ret.AllParticles[i].processed = true;
1414 // }
1415 // else {
1416 // // Decay through unknown branching ratio, presumably radiative, no hadrons, just ignore the decay products
1417 // ret.AllParticles[i].processed = true;
1418 // }
1419 // }
1420 // }
1421 // }
1422
1423 // return ret;
1424 // }
1425
1426 SimpleEvent EventGeneratorBase::PerformDecays(const SimpleEvent& evtin, const ThermalParticleSystem* TPS, const DecayerFlags& decayerFlags)
1427 {
1428 SimpleEvent ret;
1429 ret.weight = evtin.weight;
1430 ret.logweight = evtin.logweight;
1431
1432 ret.AllParticles = evtin.AllParticles;
1433 ret.DecayMap = evtin.DecayMap;
1434
1435 std::vector< std::vector<SimpleParticle> > primParticles(TPS->Particles().size(), std::vector<SimpleParticle>(0));
1436 std::vector< std::vector<int> > AllParticlesMap(TPS->Particles().size(), std::vector<int>(0));
1437 for(int i = 0; i < evtin.Particles.size(); ++i) {
1438 const SimpleParticle& particle = evtin.Particles[i];
1439 long long tid = TPS->PdgToId(particle.PDGID);
1440 if (tid != -1) {
1441 primParticles[tid].push_back(particle);
1442 AllParticlesMap[tid].push_back(i);
1443 }
1444 }
1445
1446 bool flag_repeat = true;
1447 while (flag_repeat) {
1448 flag_repeat = false;
1449 for (int i = static_cast<int>(primParticles.size()) - 1; i >= 0; --i) {
1450 if (TPS->Particles()[i].IsStable()) {
1451 for (size_t j = 0; j < primParticles[i].size(); ++j) {
1452 if (!primParticles[i][j].processed) {
1453 SimpleParticle prt = primParticles[i][j];
1454 ret.Particles.push_back(prt);
1455 primParticles[i][j].processed = true;
1456 int tid = AllParticlesMap[i][j];
1457 while (tid >= 0 && tid < ret.DecayMap.size() && ret.DecayMap[tid] != -1)
1458 tid = ret.DecayMap[tid];
1459 ret.DecayMapFinal.push_back(tid);
1460 }
1461 }
1462 }
1463 else {
1464 for (size_t j = 0; j < primParticles[i].size(); ++j) {
1465 if (!primParticles[i][j].processed) {
1466 flag_repeat = true;
1467 double DecParam = RandomGenerators::randgenMT.rand(), tsum = 0.;
1468
1469 std::vector<double> Bratios;
1470 double Width = 0.;
1471 if (primParticles[i][j].MotherPDGID != 0 ||
1472 TPS->ResonanceWidthIntegrationType() != ThermalParticle::eBW) {
1473 Bratios = TPS->Particles()[i].BranchingRatiosM(primParticles[i][j].m, false);
1474 Width = TPS->Particles()[i].ResonanceWidth();
1475 }
1476 else {
1477 Bratios = TPS->Particles()[i].BranchingRatiosM(primParticles[i][j].m, true);
1478 Width = TPS->Particles()[i].TotalWidtheBW(primParticles[i][j].m);
1479 }
1480
1481 int DecayIndex = 0;
1482 for (DecayIndex = 0; DecayIndex < static_cast<int>(Bratios.size()); ++DecayIndex) {
1483 tsum += Bratios[DecayIndex];
1484 if (tsum > DecParam) break;
1485 }
1486 if (DecayIndex < static_cast<int>(TPS->Particles()[i].Decays().size())) {
1487 std::vector<double> masses(0);
1488 std::vector<long long> pdgids(0);
1489 for (size_t di = 0; di < TPS->Particles()[i].Decays()[DecayIndex].mDaughters.size(); di++) {
1490 long long dpdg = TPS->Particles()[i].Decays()[DecayIndex].mDaughters[di];
1491 if (TPS->PdgToId(dpdg) == -1) {
1492 // Try to see if the daughter particle is a photon/lepton
1493 if (ExtraParticles::PdgToId(dpdg) == -1)
1494 continue;
1495 else
1496 masses.push_back(ExtraParticles::ParticleByPdg(dpdg).Mass());
1497 }
1498 else {
1499 masses.push_back(TPS->ParticleByPDG(dpdg).Mass());
1500 }
1501 pdgids.push_back(dpdg);
1502 }
1503
1504 // Propagate along straight line before decay
1505 double prop_dx = 0., prop_dy = 0., prop_dz = 0., prop_dt = 0.;
1506 if (decayerFlags.propagateParticles) {
1507 double ct = 0.;
1508 if (Width != 0.)
1509 ct = 1. / Width / xMath::GeVtoifm();
1510 else
1511 ct = DecayLifetimes::GetLifetime(primParticles[i][j].PDGID);
1512
1513 if (ct == 0.) {
1514 if (abs(primParticles[i][j].PDGID) != 311) {
1515 // K0s "decaying" into K0S and K0L is an exception
1516 printf("**WARNING** Could not find the lifetime for decaying particle %lld. Setting to zero...\n", primParticles[i][j].PDGID);
1517 }
1518 }
1519
1520 // Sample lifetime in rest frame
1521 ct = -ct * log(1. - RandomGenerators::randgenMT.randDblExc());
1522
1523 // Lorentz time delay
1524 double vx = primParticles[i][j].px / primParticles[i][j].p0;
1525 double vy = primParticles[i][j].py / primParticles[i][j].p0;
1526 double vz = primParticles[i][j].pz / primParticles[i][j].p0;
1527 double gamma = 1. / sqrt(1. - vx*vx - vy*vy - vz*vz);
1528 ct *= gamma;
1529
1530 // Propagate
1531 prop_dx += vx * ct;
1532 prop_dy += vy * ct;
1533 prop_dz += vz * ct;
1534 prop_dt += ct;
1535 }
1536
1537 std::vector<SimpleParticle> decres = ParticleDecaysMC::ManyBodyDecay(primParticles[i][j], masses, pdgids);
1538 for (size_t ind = 0; ind < decres.size(); ind++) {
1539 decres[ind].processed = false;
1540
1541 // Propagate before decay
1542 if (decayerFlags.propagateParticles && prop_dt != 0.0) {
1543 decres[ind].r0 += prop_dt;
1544 decres[ind].rx += prop_dx;
1545 decres[ind].ry += prop_dy;
1546 decres[ind].rz += prop_dz;
1547 }
1548
1549 if (TPS->PdgToId(decres[ind].PDGID) != -1) {
1550 int tid = TPS->PdgToId(decres[ind].PDGID);
1551 SimpleParticle& dprt = decres[ind];
1552 primParticles[tid].push_back(dprt);
1553 ret.AllParticles.push_back(dprt);
1554 AllParticlesMap[tid].push_back(static_cast<int>(ret.AllParticles.size()) - 1);
1555 ret.DecayMap.push_back(AllParticlesMap[i][j]);
1556 }
1557 else if (ExtraParticles::PdgToId(decres[ind].PDGID) != -1) {
1558 SimpleParticle& dprt = decres[ind];
1559 ret.AllParticles.push_back(dprt);
1560 ret.DecayMap.push_back(AllParticlesMap[i][j]);
1561 ret.PhotonsLeptons.push_back(dprt);
1562 }
1563 }
1564 primParticles[i][j].processed = true;
1565 }
1566 else {
1567 // Decay through unknown branching ratio, presumably radiative, no hadrons, just ignore decay products
1568 primParticles[i][j].processed = true;
1569 }
1570 }
1571 }
1572 }
1573 }
1574 }
1575
1576 return ret;
1577 }
1578
1579 std::vector<double> EventGeneratorBase::GCEMeanYields() const
1580 {
1581 std::vector<double> ret = m_THM->Densities();
1582 for(size_t i = 0; i < ret.size(); ++i) {
1583 ret[i] *= m_THM->Volume();
1584 }
1585 return ret;
1586 }
1587
1588 void EventGeneratorBase::SetVolume(double V)
1589 {
1590 if (m_Config.fEnsemble != EventGeneratorConfiguration::GCE)
1591 RescaleCEMeans(V / m_THM->Volume());
1592 m_THM->SetVolume(V);
1593 m_Config.CFOParameters.V = V;
1594 m_THM->SetCanonicalVolume(V);
1595 m_Config.CFOParameters.SVc = V;
1596 }
1597
1598 void EventGeneratorBase::RescaleCEMeans(double Vmod)
1599 {
1600 m_MeanB *= Vmod;
1601 m_MeanAB *= Vmod;
1602 m_MeanSM *= Vmod;
1603 m_MeanASM *= Vmod;
1604 m_MeanCM *= Vmod;
1605 m_MeanACM *= Vmod;
1606 m_MeanCHRMM *= Vmod;
1607 m_MeanACHRMM *= Vmod;
1608 m_MeanCHRM *= Vmod;
1609 m_MeanACHRM *= Vmod;
1610
1611 for (int i = 0; i < static_cast<int>(m_Baryons.size()); ++i) m_Baryons[i].first *= Vmod;
1612 for (int i = 0; i < static_cast<int>(m_AntiBaryons.size()); ++i) m_AntiBaryons[i].first *= Vmod;
1613 for (int i = 0; i < static_cast<int>(m_StrangeMesons.size()); ++i) m_StrangeMesons[i].first *= Vmod;
1614 for (int i = 0; i < static_cast<int>(m_AntiStrangeMesons.size()); ++i) m_AntiStrangeMesons[i].first *= Vmod;
1615 for (int i = 0; i < static_cast<int>(m_ChargeMesons.size()); ++i) m_ChargeMesons[i].first *= Vmod;
1616 for (int i = 0; i < static_cast<int>(m_AntiChargeMesons.size()); ++i) m_AntiChargeMesons[i].first *= Vmod;
1617 for (int i = 0; i < static_cast<int>(m_CharmMesons.size()); ++i) m_CharmMesons[i].first *= Vmod;
1618 for (int i = 0; i < static_cast<int>(m_AntiCharmMesons.size()); ++i) m_AntiCharmMesons[i].first *= Vmod;
1619 for (int i = 0; i < static_cast<int>(m_CharmAll.size()); ++i) m_CharmAll[i].first *= Vmod;
1620 for (int i = 0; i < static_cast<int>(m_AntiCharmAll.size()); ++i) m_AntiCharmAll[i].first *= Vmod;
1621 }
1622
1623 double EventGeneratorBase::ComputeWeight(const std::vector<int>& totals) const
1624 {
1625 // Compute the normalized weight factor due to EV/vdW interactions
1626 // If V - bN < 0, returns -1.
1627 std::vector<double>* densitiesid = NULL;
1628 std::vector<double> tmpdens;
1629 const std::vector<double>& densities = m_THM->Densities();
1630 if (m_THM->InteractionModel() != ThermalModelBase::Ideal) {
1631 tmpdens = m_DensitiesIdeal;
1632 densitiesid = &tmpdens;
1633 }
1634
1635 double ret = 1.;
1636
1637 if (m_THM->InteractionModel() == ThermalModelBase::DiagonalEV) {
1638 ThermalModelEVDiagonal* model = static_cast<ThermalModelEVDiagonal*>(m_THM);
1639 double V = m_THM->Volume();
1640 double VVN = m_THM->Volume();
1641
1642 for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i)
1643 VVN -= model->ExcludedVolume(i) * totals[i];
1644
1645 if (VVN < 0.)
1646 return -1.;
1647
1648 double weight = 1.;
1649 double logweight = 0.;
1650
1651 double normweight = 1.;
1652 double weightev = 1.;
1653 double VVNev = m_THM->Volume();
1654 for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i)
1655 VVNev -= model->ExcludedVolume(i) * densities[i] * m_THM->Volume();
1656
1657 for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i) {
1658 weight *= pow(VVN / m_THM->Volume() * densitiesid->operator[](i) / densities[i], totals[i]);
1659 if (densitiesid->operator[](i) > 0. && densities[i] > 0.)
1660 logweight += totals[i] * log(VVN / m_THM->Volume() * densitiesid->operator[](i) / densities[i]);
1661
1662 weightev *= pow(VVNev / V * densitiesid->operator[](i) / densities[i], densities[i] * V);
1663
1664 if (densitiesid->operator[](i) > 0. && densities[i] > 0.)
1665 //normweight *= pow(VVN / V, totals[i]) / pow(VVNev / V, densities[i] * V) * pow(densitiesid->operator[](i) / densities[i], totals[i] - (densities[i] * V));
1666 normweight *= pow(VVN / VVNev, totals[i]) * pow(VVNev / V, totals[i] - densities[i] * V) * pow(densitiesid->operator[](i) / densities[i], totals[i] - (densities[i] * V));
1667 }
1668
1669 m_LastWeight = weight;
1670 m_LastLogWeight = logweight;
1671 m_LastNormWeight = normweight;
1672
1673 ret = normweight;
1674 }
1675
1676
1677 if (m_THM->InteractionModel() == ThermalModelBase::CrosstermsEV) {
1678 ThermalModelEVCrossterms* model = static_cast<ThermalModelEVCrossterms*>(m_THM);
1679 double V = m_THM->Volume();
1680
1681 double weight = 1.;
1682 double logweight = 0.;
1683 double normweight = 1.;
1684 double weightev = 1.;
1685 bool fl = 1;
1686 int Nspecies = m_THM->TPS()->Particles().size();
1687 for (size_t i = 0; i < Nspecies; ++i) {
1688 double VVN = m_THM->Volume();
1689
1690 for (size_t j = 0; j < Nspecies; ++j)
1691 VVN -= model->VirialCoefficient(j, i) * totals[j];
1692
1693 if (VVN < 0.) { fl = false; break; }
1694
1695 double VVNev = m_THM->Volume();
1696 for (size_t j = 0; j < Nspecies; ++j)
1697 VVNev -= model->VirialCoefficient(j, i) * densities[j] * V;
1698
1699 weight *= pow(VVN / m_THM->Volume() * densitiesid->operator[](i) / densities[i], totals[i]);
1700 if (densitiesid->operator[](i) > 0. && densities[i] > 0.)
1701 logweight += totals[i] * log(VVN / m_THM->Volume() * densitiesid->operator[](i) / densities[i]);
1702
1703 weightev *= pow(VVNev / m_THM->Volume() * densitiesid->operator[](i) / densities[i], densities[i] * V);
1704 if (densitiesid->operator[](i) > 0. && densities[i] > 0.)
1705 normweight *= pow(VVN / VVNev, totals[i]) * pow(VVNev / V, totals[i] - densities[i] * V) * pow(densitiesid->operator[](i) / densities[i], totals[i] - (densities[i] * V));
1706 }
1707
1708 if (!fl)
1709 return -1.;
1710
1711 m_LastWeight = weight;
1712 m_LastLogWeight = logweight;
1713 m_LastNormWeight = normweight;
1714
1715 ret = normweight;
1716 }
1717
1718 if (m_THM->InteractionModel() == ThermalModelBase::QvdW) {
1719 ThermalModelVDW* model = static_cast<ThermalModelVDW*>(m_THM);
1720 double V = m_THM->Volume();
1721
1722 double weight = 1.;
1723 double logweight = 0.;
1724 double normweight = 1.;
1725 double weightvdw = 1.;
1726 bool fl = 1;
1727 for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i) {
1728 double VVN = m_THM->Volume();
1729
1730 for (size_t j = 0; j < m_THM->TPS()->Particles().size(); ++j)
1731 VVN -= model->VirialCoefficient(j, i) * totals[j];
1732
1733 if (VVN < 0.) { fl = false; break; }
1734
1735 double VVNev = m_THM->Volume();
1736 for (size_t j = 0; j < m_THM->TPS()->Particles().size(); ++j)
1737 VVNev -= model->VirialCoefficient(j, i) * densities[j] * V;
1738
1739 weight *= pow(VVN / m_THM->Volume() * densitiesid->operator[](i) / densities[i], totals[i]);
1740 if (densitiesid->operator[](i) > 0. && densities[i] > 0.)
1741 logweight += totals[i] * log(VVN / m_THM->Volume() * densitiesid->operator[](i) / densities[i]);
1742
1743 for (size_t j = 0; j < m_THM->TPS()->Particles().size(); ++j) {
1744 double aij = model->AttractionCoefficient(i, j);
1745 weight *= exp(aij * totals[j] / m_THM->Parameters().T / m_THM->Volume() * totals[i]);
1746 logweight += totals[i] * aij * totals[j] / m_THM->Parameters().T / m_THM->Volume();
1747 }
1748
1749 weightvdw *= pow(VVNev / m_THM->Volume() * densitiesid->operator[](i) / densities[i], densities[i] * V);
1750 if (densitiesid->operator[](i) > 0. && densities[i] > 0.)
1751 normweight *= pow(VVN / VVNev, totals[i]) * pow(VVNev / V, totals[i] - densities[i] * V) * pow(densitiesid->operator[](i) / densities[i], totals[i] - (densities[i] * V));
1752
1753 for (size_t j = 0; j < m_THM->TPS()->Particles().size(); ++j) {
1754 double aij = model->AttractionCoefficient(i, j);
1755 weightvdw *= exp(aij * densities[j] / m_THM->Parameters().T * densities[i] * V);
1756 normweight *= exp(aij * totals[j] / m_THM->Parameters().T / m_THM->Volume() * totals[i] - aij * densities[j] / m_THM->Parameters().T * densities[i] * V);
1757 }
1758
1759 }
1760 if (!fl)
1761 return -1.;
1762
1763 m_LastWeight = weight;
1764 m_LastLogWeight = logweight;
1765 m_LastNormWeight = normweight;
1766
1767 ret = normweight;
1768 }
1769
1770 return ret;
1771 }
1772
1773 double EventGeneratorBase::ComputeWeightNew(const std::vector<int>& totals) const
1774 {
1775 // Compute the normlaized weight factor due to EV/vdW interactions
1776 // If V - bN < 0, returns -1.
1777 std::vector<double>* densitiesid = NULL;
1778 std::vector<double> tmpdens;
1779 const std::vector<double>& densities = m_THM->Densities();
1780 if (m_THM->InteractionModel() != ThermalModelBase::Ideal) {
1781 tmpdens = m_DensitiesIdeal;
1782 densitiesid = &tmpdens;
1783 }
1784
1785 double ret = 1.;
1786
1787 if (m_THM->InteractionModel() == ThermalModelBase::DiagonalEV) {
1788 ThermalModelEVDiagonal* model = static_cast<ThermalModelEVDiagonal*>(m_THM);
1789 double V = m_THM->Volume();
1790 double VVN = m_THM->Volume();
1791
1792 for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i)
1793 VVN -= model->ExcludedVolume(i) * totals[i];
1794
1795 if (VVN < 0.)
1796 return -1.;
1797
1798 double weight = 1.;
1799 double logweight = 0.;
1800
1801 double normweight = 1.;
1802 double weightev = 1.;
1803 double VVNev = m_THM->Volume();
1804 for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i)
1805 VVNev -= model->ExcludedVolume(i) * densities[i] * m_THM->Volume();
1806
1807 for (size_t i = 0; i < m_THM->TPS()->Particles().size(); ++i) {
1808 weight *= pow(VVN / m_THM->Volume() * densitiesid->operator[](i) / densities[i], totals[i]);
1809 if (densitiesid->operator[](i) > 0. && densities[i] > 0.)
1810 logweight += totals[i] * log(VVN / m_THM->Volume() * densitiesid->operator[](i) / densities[i]);
1811
1812 weightev *= pow(VVNev / V * densitiesid->operator[](i) / densities[i], densities[i] * V);
1813
1814 if (densitiesid->operator[](i) > 0. && densities[i] > 0.)
1815 //normweight *= pow(VVN / V, totals[i]) / pow(VVNev / V, densities[i] * V) * pow(densitiesid->operator[](i) / densities[i], totals[i] - (densities[i] * V));
1816 normweight *= pow(VVN / VVNev, totals[i]) * pow(VVNev / V, totals[i] - densities[i] * V) * pow(densitiesid->operator[](i) / densities[i], totals[i] - (densities[i] * V));
1817 }
1818
1819 m_LastWeight = weight;
1820 m_LastLogWeight = logweight;
1821 m_LastNormWeight = normweight;
1822
1823 ret = normweight;
1824 }
1825
1826
1827 if (m_THM->InteractionModel() == ThermalModelBase::CrosstermsEV) {
1828 ThermalModelEVCrossterms* model = static_cast<ThermalModelEVCrossterms*>(m_THM);
1829 double V = m_THM->Volume();
1830
1831 double weight = 1.;
1832 double logweight = 0.;
1833 double normweight = 1.;
1834 double weightev = 1.;
1835 bool fl = true;
1836 int Nspecies = m_THM->TPS()->Particles().size();
1837
1838 int NEVcomp = model->EVComponentIndices().size();
1839 std::vector<int> Nscomp(NEVcomp, 0);
1840 std::vector<double> Nevscomp(NEVcomp, 0.);
1841 std::vector<double> bns(NEVcomp, 0.), bnevs(NEVcomp, 0.), dmuTs(NEVcomp, 0.);
1842 const std::vector< std::vector<double> >& virial = model->VirialMatrix();
1843
1844 for (size_t icomp = 0; icomp < NEVcomp; ++icomp) {
1845 const std::vector<int>& indis = model->EVComponentIndices()[icomp];
1846 int Nlocal = indis.size();
1847 for (size_t ilocal = 0; ilocal < Nlocal; ++ilocal) {
1848 int ip = indis[ilocal];
1849 Nscomp[icomp] += totals[ip];
1850 Nevscomp[icomp] += densities[ip] * V;
1851 }
1852
1853 if (indis.size()) {
1854 int i1 = indis[0];
1855
1856 for (size_t j = 0; j < Nspecies; ++j) {
1857 //bns[icomp] += model->VirialCoefficient(j, i1) * totals[j] / V;
1858 //bnevs[icomp] += model->VirialCoefficient(j, i1) * densities[j];
1859 bns[icomp] += virial[j][i1] * totals[j];// / V;
1860 bnevs[icomp] += virial[j][i1] * densities[j];
1861 }
1862 bns[icomp] /= V;
1863
1864 if (bns[icomp] > 1.)
1865 fl = false;
1866
1867 //dmuTs[icomp] = model->DeltaMu(i1) / model->Parameters().T;
1868 dmuTs[icomp] = log(densities[i1] / densitiesid->operator[](i1) / (1. - bnevs[icomp]));
1869 }
1870
1871 normweight *= pow((1. - bns[icomp]) / (1. - bnevs[icomp]), Nscomp[icomp]) * exp(-dmuTs[icomp] * (Nscomp[icomp] - Nevscomp[icomp]));
1872
1873 }
1874
1875 if (!fl)
1876 return -1.;
1877
1878 m_LastWeight = normweight;
1879 m_LastLogWeight = log(normweight);
1880 m_LastNormWeight = normweight;
1881
1882 ret = normweight;
1883 }
1884
1885 if (m_THM->InteractionModel() == ThermalModelBase::QvdW) {
1886 ThermalModelVDW* model = static_cast<ThermalModelVDW*>(m_THM);
1887 double V = m_THM->Volume();
1888
1889 double weight = 1.;
1890 double logweight = 0.;
1891 double normweight = 1.;
1892 double weightvdw = 1.;
1893 bool fl = true;
1894
1895 int Nspecies = m_THM->TPS()->Particles().size();
1896
1897 int NVDWcomp = model->VDWComponentIndices().size();
1898 std::vector<int> Nscomp(NVDWcomp, 0);
1899 std::vector<double> Nvdwscomp(NVDWcomp, 0.);
1900 std::vector<double> bns(NVDWcomp, 0.), bnvdws(NVDWcomp, 0.), dmuTs(NVDWcomp, 0.), aijs(NVDWcomp, 0.), aijvdws(NVDWcomp, 0.);
1901 const std::vector< std::vector<double> >& virial = model->VirialMatrix();
1902 const std::vector< std::vector<double> >& attr = model->AttractionMatrix();
1903
1904 for (size_t icomp = 0; icomp < NVDWcomp; ++icomp) {
1905 const std::vector<int>& indis = model->VDWComponentIndices()[icomp];
1906 int Nlocal = indis.size();
1907 for (size_t ilocal = 0; ilocal < Nlocal; ++ilocal) {
1908 int ip = indis[ilocal];
1909 Nscomp[icomp] += totals[ip];
1910 Nvdwscomp[icomp] += densities[ip] * V;
1911 }
1912
1913 if (indis.size()) {
1914 int i1 = indis[0];
1915
1916 for (size_t j = 0; j < Nspecies; ++j) {
1917 //bns[icomp] += model->VirialCoefficient(j, i1) * totals[j] / V;
1918 //bnevs[icomp] += model->VirialCoefficient(j, i1) * densities[j];
1919 bns[icomp] += virial[j][i1] * totals[j];// / V;
1920 bnvdws[icomp] += virial[j][i1] * densities[j];
1921
1922 aijs[icomp] += attr[i1][j] * totals[j];
1923 aijvdws[icomp] += attr[i1][j] * densities[j];
1924 }
1925 bns[icomp] /= V;
1926 aijs[icomp] /= V * model->Parameters().T;
1927 aijvdws[icomp] /= model->Parameters().T;
1928
1929
1930 if (bns[icomp] > 1.)
1931 fl = false;
1932
1933 //dmuTs[icomp] = model->DeltaMu(i1) / model->Parameters().T;
1934 dmuTs[icomp] = log(densities[i1] / densitiesid->operator[](i1) / (1. - bnvdws[icomp]));
1935 }
1936
1937 normweight *= pow((1. - bns[icomp]) / (1. - bnvdws[icomp]), Nscomp[icomp])
1938 * exp(-dmuTs[icomp] * (Nscomp[icomp] - Nvdwscomp[icomp]))
1939 * exp(aijs[icomp] * Nscomp[icomp] - aijvdws[icomp] * Nvdwscomp[icomp]);
1940 }
1941
1942 if (!fl)
1943 return -1.;
1944
1945 m_LastWeight = normweight;
1946 m_LastLogWeight = log(normweight);
1947 m_LastNormWeight = normweight;
1948
1949 ret = normweight;
1950 }
1951
1952 if (m_THM->InteractionModel() == ThermalModelBase::RealGas) {
1953 ThermalModelRealGas* model = static_cast<ThermalModelRealGas*>(m_THM);
1954 double V = m_THM->Volume();
1955
1956 double weight = 1.;
1957 double logweight = 0.;
1958 double normweight = 1.;
1959 double weightvdw = 1.;
1960 bool fl = true;
1961
1962 int Nspecies = m_THM->TPS()->Particles().size();
1963
1964 ExcludedVolumeModelMultiBase* evmod = model->ExcludedVolumeModel();
1965 int Nevcomp = evmod->ComponentsNumber();
1966 const std::vector<int>& evinds = evmod->ComponentIndices();
1967 const std::vector<int>& evindsfrom = evmod->ComponentIndicesFrom();
1968 std::vector<int> Nscomp(Nevcomp, 0);
1969 std::vector<double> Nvdwscomp(Nevcomp, 0.);
1970 //std::vector<double> dmuTs(Nspecies, 0.);
1971
1972 std::vector<double> ev_favs(Nevcomp, 0.);
1973 for (size_t icomp = 0; icomp < Nevcomp; ++icomp) {
1974 ev_favs[icomp] = evmod->f(evindsfrom[icomp]);
1975 }
1976 std::vector<double> densities_sampled(Nspecies);
1977 for (size_t j = 0; j < Nspecies; ++j) {
1978 densities_sampled[j] = totals[j] / V;
1979 }
1980 evmod->SetDensities(densities_sampled);
1981 std::vector<double> ev_fsampled(Nevcomp, 0.);
1982 for (size_t icomp = 0; icomp < Nevcomp; ++icomp) {
1983 ev_fsampled[icomp] = evmod->f(evindsfrom[icomp]);
1984 }
1985 evmod->SetDensities(m_THM->Densities());
1986
1987 MeanFieldModelMultiBase* mfmod = model->MeanFieldModel();
1988 double mf_vav = mfmod->v();
1989 mfmod->SetDensities(densities_sampled);
1990 double mf_vsampled = mfmod->v();
1991 mfmod->SetDensities(m_THM->Densities());
1992
1993 for (size_t j = 0; j < Nspecies; ++j) {
1994 int jcomp = evinds[j];
1995 normweight *= pow(ev_fsampled[jcomp] / ev_favs[jcomp], totals[j]);
1996
1997 double dmuT = 0.;
1998 if (densitiesid->operator[](j))
1999 dmuT = log(densities[j] / densitiesid->operator[](j) / ev_favs[jcomp]);
2000 normweight *= exp(-dmuT * (totals[j] - densities[j] * V));
2001 }
2002
2003 normweight *= exp(-V * (mf_vsampled - mf_vav) / m_THM->Parameters().T);
2004
2005 if (!fl)
2006 return -1.;
2007
2008 m_LastWeight = normweight;
2009 m_LastLogWeight = log(normweight);
2010 m_LastNormWeight = normweight;
2011
2012 ret = normweight;
2013 }
2014
2015 //printf("Weight: %lf\n", ret);
2016
2017 return ret;
2018 }
2019
2034
2035} // namespace thermalfist
ExcludedVolumeHelper.h
Contains some functions to deal with excluded volumes.
thermalfist::EventGeneratorBase::ComputeEVRadii
std::vector< std::vector< double > > ComputeEVRadii() const
Definition EventGeneratorBase.cpp:866
thermalfist::EventGeneratorBase::GCEMeanYields
virtual std::vector< double > GCEMeanYields() const
The grand-canonical mean yields.
Definition EventGeneratorBase.cpp:1579
thermalfist::EventGeneratorBase::SampleParticleByPdg
virtual SimpleParticle SampleParticleByPdg(long long pdgid) const
Samples the position and momentum of a particle species with given pdg code.
Definition EventGeneratorBase.cpp:1189
thermalfist::EventGeneratorBase::m_BWGens
std::vector< RandomGenerators::ThermalBreitWignerGenerator * > m_BWGens
Definition EventGeneratorBase.h:348
thermalfist::EventGeneratorBase::GenerateTotalsSCESubVolume
std::vector< int > GenerateTotalsSCESubVolume(double VolumeSC) const
Definition EventGeneratorBase.cpp:607
thermalfist::EventGeneratorBase::PerformDecays
static SimpleEvent PerformDecays(const SimpleEvent &evtin, const ThermalParticleSystem *TPS, const DecayerFlags &decayerFlags=DecayerFlags())
Performs decays of all unstable particles until only stable ones left.
Definition EventGeneratorBase.cpp:1426
thermalfist::EventGeneratorBase::SampleYields
virtual std::pair< std::vector< int >, double > SampleYields() const
Samples the primordial yields for each particle species.
Definition EventGeneratorBase.cpp:1160
thermalfist::EventGeneratorBase::GetEvent
virtual SimpleEvent GetEvent(bool PerformDecays=true) const
Generates a single event.
Definition EventGeneratorBase.cpp:1298
thermalfist::EventGeneratorBase::m_MomentumGens
std::vector< RandomGenerators::ParticleMomentumGenerator * > m_MomentumGens
Vector of momentum generators for each particle species.
Definition EventGeneratorBase.h:344
thermalfist::EventGeneratorBase::CheckEVOverlap
bool CheckEVOverlap(const std::vector< SimpleParticle > &evt, const SimpleParticle &cand, const std::vector< int > &ids=std::vector< int >(), const std::vector< std::vector< double > > &radii=std::vector< std::vector< double > >()) const
Definition EventGeneratorBase.cpp:884
thermalfist::EventGeneratorBase::SampleParticle
virtual SimpleParticle SampleParticle(int id) const
Samples the position and momentum of a particle species i.
Definition EventGeneratorBase.cpp:1168
thermalfist::EventGeneratorBase::ClearMomentumGenerators
void ClearMomentumGenerators()
Clears the momentum generators for all particles.
Definition EventGeneratorBase.cpp:72
thermalfist::EventGeneratorBase::m_Config
EventGeneratorConfiguration m_Config
Definition EventGeneratorBase.h:337
thermalfist::EventGeneratorBase::SetMomentumGenerators
virtual void SetMomentumGenerators()
Sets the momentum generators for all particles. Overloaded.
Definition EventGeneratorBase.h:133
thermalfist::EventGeneratorBase::GenerateTotalsCCESubVolume
std::vector< int > GenerateTotalsCCESubVolume(double VolumeSC) const
Definition EventGeneratorBase.cpp:801
thermalfist::EventGeneratorBase::CheckSetParameters
virtual void CheckSetParameters()
Sets the hypersurface parameters.
Definition EventGeneratorBase.h:268
thermalfist::EventGeneratorBase::ComputeWeight
double ComputeWeight(const std::vector< int > &totals) const
Definition EventGeneratorBase.cpp:1623
thermalfist::EventGeneratorBase::m_DensitiesIdeal
std::vector< double > m_DensitiesIdeal
Ideal gas densities used for sampling an interacting HRG.
Definition EventGeneratorBase.h:341
thermalfist::EventGeneratorBase::SampleMomenta
virtual SimpleEvent SampleMomenta(const std::vector< int > &yields) const
Samples the momenta of the particles and returns the sampled list of particles as an event.
Definition EventGeneratorBase.cpp:1198
thermalfist::EventGeneratorBase::GenerateTotalsGCE
std::vector< int > GenerateTotalsGCE() const
Definition EventGeneratorBase.cpp:463
thermalfist::EventGeneratorBase::GenerateTotalsSCE
std::vector< int > GenerateTotalsSCE() const
Definition EventGeneratorBase.cpp:483
thermalfist::EventGeneratorBase::GenerateTotalsCE
std::vector< int > GenerateTotalsCE() const
Definition EventGeneratorBase.cpp:924
thermalfist::EventGeneratorBase::SetConfiguration
void SetConfiguration(ThermalParticleSystem *TPS, const EventGeneratorConfiguration &config)
Sets the event generator configuration.
Definition EventGeneratorBase.cpp:88
thermalfist::EventGeneratorBase::ComputeWeightNew
double ComputeWeightNew(const std::vector< int > &totals) const
Definition EventGeneratorBase.cpp:1773
thermalfist::EventGeneratorBase::GenerateTotals
std::vector< int > GenerateTotals() const
Definition EventGeneratorBase.cpp:428
thermalfist::EventGeneratorBase::GenerateTotalsCCE
std::vector< int > GenerateTotalsCCE() const
Definition EventGeneratorBase.cpp:674
thermalfist::EventGeneratorBase::SampleMomentaWithShuffle
virtual SimpleEvent SampleMomentaWithShuffle(const std::vector< int > &yields) const
Samples the momenta of the particles and returns the sampled list of particles as an event.
Definition EventGeneratorBase.cpp:1235
thermalfist::EventGeneratorBase::SetVolume
void SetVolume(double V)
Set system volume.
Definition EventGeneratorBase.cpp:1588
thermalfist::EventGeneratorBase::RescaleCEMeans
void RescaleCEMeans(double Vmod)
Rescale the precalculated GCE means.
Definition EventGeneratorBase.cpp:1598
thermalfist::EventGeneratorBase::SetParameters
virtual void SetParameters()
Sets up the event generator ready for production.
Definition EventGeneratorBase.cpp:860
thermalfist::ExcludedVolumeModelBase
Base class implementing the ideal gas.
Definition ExcludedVolumeModels.h:26
thermalfist::ExcludedVolumeModelCS
Class implementing auxiliary functions for the Carnahan-Starling excluded volume model.
Definition ExcludedVolumeModels.h:66
thermalfist::ExcludedVolumeModelCrosstermsGeneralized
Implementation of a crossterms generalized excluded volume model.
Definition ExcludedVolumeModelsMulti.h:528
thermalfist::ExcludedVolumeModelMultiBase
Base class for multi-component excluded volume models.
Definition ExcludedVolumeModelsMulti.h:25
thermalfist::ExcludedVolumeModelMultiBase::f
virtual double f(int i) const
Calculates the suppression factor for species i.
Definition ExcludedVolumeModelsMulti.h:47
thermalfist::ExcludedVolumeModelMultiBase::ComponentsNumber
virtual const int ComponentsNumber() const
Gets the number of components.
Definition ExcludedVolumeModelsMulti.h:149
thermalfist::ExcludedVolumeModelMultiBase::ComponentIndices
virtual const std::vector< int > & ComponentIndices() const
Gets the component indices.
Definition ExcludedVolumeModelsMulti.h:135
thermalfist::ExcludedVolumeModelMultiBase::ComponentIndicesFrom
virtual const std::vector< int > & ComponentIndicesFrom() const
Gets the component indices from.
Definition ExcludedVolumeModelsMulti.h:142
thermalfist::ExcludedVolumeModelMultiBase::SetDensities
virtual void SetDensities(const std::vector< double > &n)
Sets the densities of particle species.
Definition ExcludedVolumeModelsMulti.h:128
thermalfist::ExcludedVolumeModelTVM
Class implementing auxiliary functions for the VDW excluded volume model truncated at n^3....
Definition ExcludedVolumeModels.h:135
thermalfist::ExcludedVolumeModelVDW
Class implementing auxiliary functions for the van der Waals excluded volume model.
Definition ExcludedVolumeModels.h:44
thermalfist::ExcludedVolumeModelVirial
Class implementing auxiliary functions for the VDW excluded volume model truncated at n^2.
Definition ExcludedVolumeModels.h:107
thermalfist::MeanFieldModelMultiBase
Base class for multi-component mean field models.
Definition MeanFieldModelsMulti.h:25
thermalfist::MeanFieldModelMultiBase::SetDensities
virtual void SetDensities(const std::vector< double > &n)
Sets the densities of particle species.
Definition MeanFieldModelsMulti.cpp:233
thermalfist::MeanFieldModelMultiBase::v
virtual double v() const
Calculates the mean field value.
Definition MeanFieldModelsMulti.h:44
thermalfist::MeanFieldModelMultiVDW
Implementation of the van der Waals mean field model for multiple components.
Definition MeanFieldModelsMulti.h:139
thermalfist::ThermalModelBase::QvdW
@ QvdW
Quantum van der Waals model.
Definition ThermalModelBase.h:50
thermalfist::ThermalModelBase::DiagonalEV
@ DiagonalEV
Diagonal excluded volume model.
Definition ThermalModelBase.h:48
thermalfist::ThermalModelBase::CrosstermsEV
@ CrosstermsEV
Crossterms excluded volume model.
Definition ThermalModelBase.h:49
thermalfist::ThermalModelBase::Parameters
const ThermalModelParameters & Parameters() const
Definition ThermalModelBase.h:121
thermalfist::ThermalModelEVCrossterms
Class implementing the crossterms excluded-volume model.
Definition ThermalModelEVCrossterms.h:36
thermalfist::ThermalModelEVCrossterms::EVComponentIndices
const std::vector< std::vector< int > > & EVComponentIndices() const
Definition ThermalModelEVCrossterms.h:71
thermalfist::ThermalModelEVDiagonal
Class implementing the diagonal excluded-volume model.
Definition ThermalModelEVDiagonal.h:46
thermalfist::ThermalModelEVDiagonal::ExcludedVolume
double ExcludedVolume(int i) const
The excluded volume parameter of particle species i.
Definition ThermalModelEVDiagonal.cpp:122
thermalfist::ThermalModelIdeal
Class implementing the Ideal HRG model.
Definition ThermalModelIdeal.h:20
thermalfist::ThermalModelRealGas
Class implementing the quantum real gas HRG model.
Definition ThermalModelRealGas.h:25
thermalfist::ThermalModelRealGas::ExcludedVolumeModel
ExcludedVolumeModelMultiBase * ExcludedVolumeModel() const
Get the excluded volume model used in the real gas HRG model.
Definition ThermalModelRealGas.h:60
thermalfist::ThermalModelRealGas::MeanFieldModel
MeanFieldModelMultiBase * MeanFieldModel() const
Get the mean field model used in the real gas HRG model.
Definition ThermalModelRealGas.h:67
thermalfist::ThermalModelVDW
Class implementing the quantum van der Waals HRG model.
Definition ThermalModelVDW.h:37
thermalfist::ThermalModelVDW::AttractionCoefficient
double AttractionCoefficient(int i, int j) const
QvdW mean field attraction coefficient .
Definition ThermalModelVDW.cpp:1155
thermalfist::ThermalModelVDW::VDWComponentIndices
const std::vector< std::vector< int > > & VDWComponentIndices() const
Definition ThermalModelVDW.h:207
thermalfist::ThermalModelVDW::VirialMatrix
const std::vector< std::vector< double > > & VirialMatrix() const
Definition ThermalModelVDW.h:209
thermalfist::ThermalModelVDW::AttractionMatrix
const std::vector< std::vector< double > > & AttractionMatrix() const
Definition ThermalModelVDW.h:210
thermalfist::ThermalModelVDW::VirialCoefficient
double VirialCoefficient(int i, int j) const
Excluded volume coefficient .
Definition ThermalModelVDW.cpp:1148
thermalfist::ThermalParticle
Class containing all information about a particle specie.
Definition ThermalParticle.h:64
thermalfist::ThermalParticle::PdgId
long long PdgId() const
Particle's Particle Data Group (PDG) ID number.
Definition ThermalParticle.h:348
thermalfist::ThermalParticle::Statistics
int Statistics() const
Particle's statistics.
Definition ThermalParticle.h:368
thermalfist::ThermalParticle::Mass
double Mass() const
Particle's mass [GeV].
Definition ThermalParticle.h:390
thermalfist::ThermalParticle::GetResonanceWidthIntegrationType
ResonanceWidthIntegration GetResonanceWidthIntegrationType() const
Resonance width integration scheme used to treat finite resonance widths.
Definition ThermalParticle.h:526
thermalfist::ThermalParticle::eBW
@ eBW
Energy-dependent Breit-Wigner scheme (eBW)
Definition ThermalParticle.h:87
thermalfist::ThermalParticle::ZeroWidth
@ ZeroWidth
Zero-width approximation.
Definition ThermalParticle.h:83
thermalfist::ThermalParticle::ZeroWidthEnforced
bool ZeroWidthEnforced() const
Whether zero-width approximation is enforced for this particle species.
Definition ThermalParticle.cpp:59
thermalfist::ThermalParticleSystem
Class containing the particle list.
Definition ThermalParticleSystem.h:30
thermalfist::ThermalParticleSystem::ResonanceWidthIntegrationType
ThermalParticle::ResonanceWidthIntegration ResonanceWidthIntegrationType() const
Definition ThermalParticleSystem.h:365
thermalfist::ThermalParticleSystem::PdgToId
int PdgToId(long long pdgid) const
Transforms PDG ID to a 0-based particle id number.
Definition ThermalParticleSystem.cpp:1348
thermalfist::ThermalParticleSystem::Particles
const std::vector< ThermalParticle > & Particles() const
Returns the vector of all particle species.
Definition ThermalParticleSystem.h:408
thermalfist::ThermalParticleSystem::ParticleByPDG
const ThermalParticle & ParticleByPDG(long long pdgid) const
ThermalParticle object corresponding to particle species with a provided PDG ID.
Definition ThermalParticleSystem.cpp:1330
thermalfist::CuteHRGHelper::rv
double rv(double v)
Computes the radius parameter from a given excluded volume parameter value.
Definition ExcludedVolumeHelper.h:41
thermalfist::ExtraParticles::ParticleByPdg
const ThermalParticle & ParticleByPdg(long long pdg)
Definition ExtraParticles.cpp:34
thermalfist::ParticleDecaysMC::ParticleDistanceSquared
double ParticleDistanceSquared(const SimpleParticle &part1, const SimpleParticle &part2)
Return the square of the distance between particles at equal time.
Definition ParticleDecaysMC.cpp:242
thermalfist::ParticleDecaysMC::ManyBodyDecay
std::vector< SimpleParticle > ManyBodyDecay(const SimpleParticle &Mother, std::vector< double > masses, std::vector< long long > pdgs)
Samples the decay products of a many-body decay.
Definition ParticleDecaysMC.cpp:167
thermalfist::RandomGenerators::RandomPoisson
int RandomPoisson(double mean)
Generates random integer distributed by Poisson with specified mean Uses randgenMT.
Definition RandomGenerators.cpp:27
thermalfist::RandomGenerators::randgenMT
MTRand randgenMT
The Mersenne Twister random number generator.
Definition RandomGenerators.cpp:19
thermalfist::RandomGenerators::rng_std
std::mt19937 rng_std
The Mersenne Twister random number generator from STD library.
thermalfist::xMath::GeVtoifm
constexpr double GeVtoifm()
A constant to transform GeV into fm .
Definition xMath.h:25
thermalfist
The main namespace where all classes and functions of the Thermal-FIST library reside.
Definition CosmicEoS.h:9
thermalfist::ThermalModelVDWFull
ThermalModelVDW ThermalModelVDWFull
For backward compatibility.
Definition ThermalModelVDW.h:346
thermalfist::LorentzBoost
std::vector< double > LorentzBoost(const std::vector< double > &fourvector, double vx, double vy, double vz)
Performs a Lorentz boost on a four-vector.
Definition EventGeneratorBase.cpp:32
thermalfist::EventGeneratorConfiguration
Structure containing the thermal event generator configuration.
Definition EventGeneratorBase.h:54
thermalfist::EventGeneratorConfiguration::fEnsemble
Ensemble fEnsemble
The statistical ensemble used.
Definition EventGeneratorBase.h:74
thermalfist::EventGeneratorConfiguration::B
int B
The total values of conserved charges in the CE.
Definition EventGeneratorBase.h:83
thermalfist::EventGeneratorConfiguration::fUseEVRejectionMultiplicity
bool fUseEVRejectionMultiplicity
Whether to use rejection sampling instead of importance sampling for the EV multiplicity sampling.
Definition EventGeneratorBase.h:107
thermalfist::EventGeneratorConfiguration::RealGasExcludedVolumePrescription
int RealGasExcludedVolumePrescription
The type of generalized excluded volume model prescription.
Definition EventGeneratorBase.h:98
thermalfist::EventGeneratorConfiguration::fUseEVUseSPRApproximation
bool fUseEVUseSPRApproximation
Whether to use the SPR (single-particle rejection) approximation for the EV effects in coordinate spa...
Definition EventGeneratorBase.h:113
thermalfist::EventGeneratorConfiguration::CFOParameters
ThermalModelParameters CFOParameters
The chemical freeze-out parameters.
Definition EventGeneratorBase.h:80
thermalfist::EventGeneratorConfiguration::CanonicalB
bool CanonicalB
Mixed-canonical configuration (full canonical by default)
Definition EventGeneratorBase.h:86
thermalfist::EventGeneratorConfiguration::fModelType
ModelType fModelType
The type of interaction model.
Definition EventGeneratorBase.h:77
thermalfist::EventGeneratorConfiguration::fUsePCE
bool fUsePCE
Whether partial chemical equilibrium (PCE) is used.
Definition EventGeneratorBase.h:101
thermalfist::EventGeneratorConfiguration::fUseEVRejectionCoordinates
bool fUseEVRejectionCoordinates
Whether to use rejection sampling in the coordinate space to model EV effects.
Definition EventGeneratorBase.h:110
thermalfist::EventGeneratorConfiguration::CrosstermsEV
@ CrosstermsEV
Crossterms excluded-volume.
Definition EventGeneratorBase.h:67
thermalfist::EventGeneratorConfiguration::fUseGCEConservedCharges
bool fUseGCEConservedCharges
Whether to calculate total conserved charge values from GCE.
Definition EventGeneratorBase.h:116
thermalfist::RandomGenerators::BesselDistributionGenerator::RandomBesselDevroye1
static int RandomBesselDevroye1(double a, int nu, MTRand &rangen)
Definition RandomGenerators.cpp:669
thermalfist::SimpleEvent
Structure holding information about a single event in the event generator.
Definition SimpleEvent.h:20
thermalfist::SimpleEvent::weight
double weight
Event weight factor.
Definition SimpleEvent.h:22
thermalfist::SimpleEvent::Particles
std::vector< SimpleParticle > Particles
Vector of all final particles in the event.
Definition SimpleEvent.h:28
thermalfist::SimpleEvent::logweight
double logweight
Log of the event weight factor.
Definition SimpleEvent.h:25
thermalfist::SimpleEvent::AllParticles
std::vector< SimpleParticle > AllParticles
Vector of all particles which ever appeared in the event (including those that decay and photons/lept...
Definition SimpleEvent.h:31
thermalfist::SimpleEvent::DecayMap
std::vector< int > DecayMap
Definition SimpleEvent.h:38
thermalfist::SimpleEvent::PhotonsLeptons
std::vector< SimpleParticle > PhotonsLeptons
Vector of all decay photons/leptons.
Definition SimpleEvent.h:34
thermalfist::SimpleEvent::DecayMapFinal
std::vector< int > DecayMapFinal
Vector for each Particles element pointing to the index of the primordial resonance from which this p...
Definition SimpleEvent.h:41
thermalfist::SimpleParticle
Structure holding information about a single particle in the event generator.
Definition SimpleParticle.h:15
thermalfist::ThermalModelParameters
Structure containing all thermal parameters of the model.
Definition ThermalModelParameters.h:18
xMath.h
Contains some extra mathematical functions used in the code.
Generated by doxygen 1.13.2