| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35
36 #include "globals.hh"
37
38 #include "G4INCLPhaseSpaceRauboldLynch.hh"
39 #include "G4INCLRandom.hh"
40 #include "G4INCLKinematicsUtils.hh"
41 #include <algorithm>
42 #include <numeric>
43 #include <functional>
44
45 namespace G4INCL {
46
47 const G4double PhaseSpaceRauboldLynch::wMaxMasslessX[PhaseSpaceRauboldLynch::wMaxNE] = {
48 0.01,
49 0.017433288222,
50 0.0303919538231,
51 0.0529831690628,
52 0.0923670857187,
53 0.161026202756,
54 0.280721620394,
55 0.489390091848,
56 0.853167852417,
57 1.48735210729,
58 2.5929437974,
59 4.52035365636,
60 7.88046281567,
61 13.7382379588,
62 23.9502661999,
63 41.7531893656,
64 72.7895384398,
65 126.896100317,
66 221.221629107,
67 385.662042116,
68 672.33575365,
69 1172.10229753,
70 2043.35971786,
71 3562.24789026,
72 6210.16941892,
73 10826.3673387,
74 18873.9182214,
75 32903.4456231,
76 57361.5251045,
77 100000.0
78 };
79
80 const G4double PhaseSpaceRauboldLynch::wMaxMasslessY[PhaseSpaceRauboldLynch::wMaxNE] = {
81 -4.69180212056,
82 -4.13600582212,
83 -3.58020940928,
84 -3.02441303018,
85 -2.46861663471,
86 -1.91282025644,
87 -1.35702379603,
88 -0.801227342882,
89 -0.245430866403,
90 0.310365269464,
91 0.866161720461,
92 1.42195839972,
93 1.97775459763,
94 2.5335510251,
95 3.08934734235,
96 3.64514378357,
97 4.20094013304,
98 4.75673663205,
99 5.3125329676,
100 5.86832950014,
101 6.42412601468,
102 6.97992197797,
103 7.53571856765,
104 8.09151503031,
105 8.64731144398,
106 9.20310774148,
107 9.7589041009,
108 10.314700625,
109 10.8704971207,
110 11.4262935304
111 };
112
113 const G4double PhaseSpaceRauboldLynch::wMaxCorrectionX[PhaseSpaceRauboldLynch::wMaxNE] = {
114 8.77396538453e-07,
115 1.52959067398e-06,
116 2.66657950812e-06,
117 4.6487249132e-06,
118 8.10425612766e-06,
119 1.41283832898e-05,
120 2.46304178003e-05,
121 4.2938917254e-05,
122 7.4856652043e-05,
123 0.00013049975904,
124 0.000227503991225,
125 0.000396614265067,
126 0.000691429079588,
127 0.00120538824295,
128 0.00210138806588,
129 0.00366341038188,
130 0.00638652890627,
131 0.0111338199161,
132 0.0194099091609,
133 0.0338378540766,
134 0.0589905062931,
135 0.102839849857,
136 0.179283674326,
137 0.312550396803,
138 0.544878115136,
139 0.949901722703,
140 1.65599105145,
141 2.88693692929,
142 5.03288035671,
143 8.77396538453
144 };
145 const G4double PhaseSpaceRauboldLynch::wMaxCorrectionY[PhaseSpaceRauboldLynch::wMaxNE] = {
146 7.28682764024,
147 7.0089298602,
148 6.7310326046,
149 6.45313436085,
150 6.17523713298,
151 5.89734080335,
152 5.61944580818,
153 5.34155326611,
154 5.06366530628,
155 4.78578514804,
156 4.50791742491,
157 4.23007259554,
158 3.97566018117,
159 3.72116551935,
160 3.44355099732,
161 3.1746329047,
162 2.89761210229,
163 2.62143335535,
164 2.34649707964,
165 2.07366126445,
166 1.8043078447,
167 1.54056992953,
168 1.27913996411,
169 0.981358535322,
170 0.73694629682,
171 0.587421056631,
172 0.417178268911,
173 0.280881750176,
174 0.187480311508,
175 0.116321254846
176 };
177
178 const G4double PhaseSpaceRauboldLynch::wMaxInterpolationMargin = std::log(1.5);
179
180 PhaseSpaceRauboldLynch::PhaseSpaceRauboldLynch() :
181 nParticles(0),
182 sqrtS(0.),
183 availableEnergy(0.),
184 maxGeneratedWeight(0.)
185 {
186 std::vector<G4double> wMaxMasslessXV(wMaxMasslessX, wMaxMasslessX + wMaxNE);
187 std::vector<G4double> wMaxMasslessYV(wMaxMasslessY, wMaxMasslessY + wMaxNE);
188 wMaxMassless = new InterpolationTable(wMaxMasslessXV, wMaxMasslessYV);
189 std::vector<G4double> wMaxCorrectionXV(wMaxCorrectionX, wMaxCorrectionX + wMaxNE);
190 std::vector<G4double> wMaxCorrectionYV(wMaxCorrectionY, wMaxCorrectionY + wMaxNE);
191 wMaxCorrection = new InterpolationTable(wMaxCorrectionXV, wMaxCorrectionYV);
192
193 // initialize the precalculated coefficients
194 prelog[0] = 0.;
195 for(size_t i=1; i<wMaxNP; ++i) {
196 prelog[i] = -std::log(G4double(i));
197 }
198 }
199
200 PhaseSpaceRauboldLynch::~PhaseSpaceRauboldLynch() {
201 delete wMaxMassless;
202 delete wMaxCorrection;
203 }
204
205 void PhaseSpaceRauboldLynch::generate(const G4double sqs, ParticleList &particles) {
206 maxGeneratedWeight = 0.;
207
208 sqrtS = sqs;
209
210 // initialize the structures containing the particle masses
211 initialize(particles);
212
213 // initialize the maximum weight
214 const G4double weightMax = computeMaximumWeightParam();
215
216 const G4int maxIter = 500;
217 G4int iter = 0;
218 G4double weight, r;
219 do {
220 weight = computeWeight();
221 maxGeneratedWeight = std::max(weight, maxGeneratedWeight);
222 // assert(weight<=weightMax);
223 r = Random::shoot();
224 } while(++iter<maxIter && r*weightMax>weight); /* Loop checking, 10.07.2015, D.Mancusi */
225
226 #ifndef INCLXX_IN_GEANT4_MODE
227 if(iter>=maxIter) {
228 INCL_WARN("Number of tries exceeded in PhaseSpaceRauboldLynch::generate()\n"
229 << "nParticles=" << nParticles << ", sqrtS=" << sqrtS << '\n');
230 }
231 #endif
232
233 generateEvent(particles);
234 }
235
236 void PhaseSpaceRauboldLynch::initialize(ParticleList &particles) {
237 nParticles = particles.size();
238 // assert(nParticles>2);
239
240 // masses and sum of masses
241 masses.resize(nParticles);
242 sumMasses.resize(nParticles);
243 std::transform(particles.begin(), particles.end(), masses.begin(), std::mem_fn(&Particle::getMass));
244 std::partial_sum(masses.begin(), masses.end(), sumMasses.begin());
245
246 // sanity check
247 availableEnergy = sqrtS-sumMasses[nParticles-1];
248 // assert(availableEnergy>-1.e-5);
249 if(availableEnergy<0.)
250 availableEnergy = 0.;
251
252 rnd.resize(nParticles);
253 invariantMasses.resize(nParticles);
254 momentaCM.resize(nParticles-1);
255 }
256
257 G4double PhaseSpaceRauboldLynch::computeMaximumWeightNaive() {
258 // compute the maximum weight
259 G4double eMMax = sqrtS + masses[0];
260 G4double eMMin = 0.;
261 G4double wMax = 1.;
262 for(size_t i=1; i<nParticles; i++) {
263 eMMin += masses[i-1];
264 eMMax += masses[i];
265 wMax *= KinematicsUtils::momentumInCM(eMMax, eMMin, masses[i]);
266 }
267 return wMax;
268 }
269
270 G4double PhaseSpaceRauboldLynch::computeMaximumWeightParam() {
271 #ifndef INCLXX_IN_GEANT4_MODE
272 if(nParticles>=wMaxNP) {
273 INCL_WARN("The requested number of particles (" << nParticles << ") requires extrapolation the tables in PhaseSpaceRauboldLynch. YMMV." << '\n');
274 }
275
276 if(availableEnergy>wMaxMasslessX[wMaxNE-1] || availableEnergy<wMaxMasslessX[0] ) {
277 INCL_WARN("The requested available energy (" << availableEnergy << " MeV) requires extrapolation the tables in PhaseSpaceRauboldLynch. YMMV." << '\n');
278 }
279 #endif
280
281 // compute the maximum weight
282 const G4double logMassless = ((*wMaxMassless)(availableEnergy)+prelog[nParticles])*(nParticles-1);
283 const G4double reducedSqrtS = availableEnergy/sumMasses[nParticles-1];
284 const G4double correction = (*wMaxCorrection)(reducedSqrtS);
285 const G4double wMax = std::exp(correction*G4double(nParticles-1) + logMassless + wMaxInterpolationMargin);
286 if(wMax>0.)
287 return wMax;
288 else
289 return computeMaximumWeightNaive();
290 }
291
292 G4double PhaseSpaceRauboldLynch::computeWeight() {
293 // Generate nParticles-2 sorted random numbers, bracketed by 0 and 1
294 rnd[0] = 0.;
295 for(size_t i=1; i<nParticles-1; ++i)
296 rnd[i] = Random::shoot();
297 rnd[nParticles-1] = 1.;
298 std::sort(rnd.begin()+1, rnd.begin()+nParticles-1);
299
300 // invariant masses
301 for(size_t i=0; i<nParticles; ++i)
302 invariantMasses[i] = rnd[i]*availableEnergy + sumMasses[i];
303
304 // compute the CM momenta and the weight for the current event
305 G4double weight=KinematicsUtils::momentumInCM(invariantMasses[1], invariantMasses[0], masses[1]);
306 momentaCM[0] = weight;
307 for(size_t i=1; i<nParticles-1; ++i) {
308 G4double momentumCM;
309 // assert(invariantMasses[i+1]-invariantMasses[i]-masses[i+1]> - 1.E-8);
310 if(invariantMasses[i+1]-invariantMasses[i]-masses[i+1] < 0.) momentumCM = 0.;
311 else momentumCM = KinematicsUtils::momentumInCM(invariantMasses[i+1], invariantMasses[i], masses[i+1]);
312 momentaCM[i] = momentumCM;
313 weight *= momentumCM;
314 }
315
316 return weight;
317 }
318
319 void PhaseSpaceRauboldLynch::generateEvent(ParticleList &particles) {
320 Particle *p = particles[0];
321 ThreeVector momentum = Random::normVector(momentaCM[0]);
322 p->setMomentum(momentum);
323 p->adjustEnergyFromMomentum();
324
325 ThreeVector boostV;
326
327 for(size_t i=1; i<nParticles; ++i) {
328 p = particles[i];
329 p->setMomentum(-momentum);
330 p->adjustEnergyFromMomentum();
331
332 if(i==nParticles-1)
333 break;
334
335 momentum = Random::normVector(momentaCM[i]);
336
337 const G4double iM = invariantMasses[i];
338 const G4double recoilE = std::sqrt(momentum.mag2() + iM*iM);
339 boostV = -momentum/recoilE;
340 for(size_t j=0; j<=i; ++j)
341 particles[j]->boost(boostV);
342 }
343 }
344
345 G4double PhaseSpaceRauboldLynch::getMaxGeneratedWeight() const {
346 return maxGeneratedWeight;
347 }
348 }
349