| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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 H
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::wMaxM
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::wMaxM
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::wMaxC
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::wMaxC
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::wMaxI
179
180 PhaseSpaceRauboldLynch::PhaseSpaceRauboldLyn
181 nParticles(0),
182 sqrtS(0.),
183 availableEnergy(0.),
184 maxGeneratedWeight(0.)
185 {
186 std::vector<G4double> wMaxMasslessXV(wMaxM
187 std::vector<G4double> wMaxMasslessYV(wMaxM
188 wMaxMassless = new InterpolationTable(wMax
189 std::vector<G4double> wMaxCorrectionXV(wMa
190 std::vector<G4double> wMaxCorrectionYV(wMa
191 wMaxCorrection = new InterpolationTable(wM
192
193 // initialize the precalculated coefficien
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::~PhaseSpaceRauboldLy
201 delete wMaxMassless;
202 delete wMaxCorrection;
203 }
204
205 void PhaseSpaceRauboldLynch::generate(const
206 maxGeneratedWeight = 0.;
207
208 sqrtS = sqs;
209
210 // initialize the structures containing th
211 initialize(particles);
212
213 // initialize the maximum weight
214 const G4double weightMax = computeMaximumW
215
216 const G4int maxIter = 500;
217 G4int iter = 0;
218 G4double weight, r;
219 do {
220 weight = computeWeight();
221 maxGeneratedWeight = std::max(weight, ma
222 // assert(weight<=weightMax);
223 r = Random::shoot();
224 } while(++iter<maxIter && r*weightMax>weig
225
226 #ifndef INCLXX_IN_GEANT4_MODE
227 if(iter>=maxIter) {
228 INCL_WARN("Number of tries exceeded in P
229 << "nParticles=" << nParticles
230 }
231 #endif
232
233 generateEvent(particles);
234 }
235
236 void PhaseSpaceRauboldLynch::initialize(Part
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(), particle
244 std::partial_sum(masses.begin(), masses.en
245
246 // sanity check
247 availableEnergy = sqrtS-sumMasses[nPartic
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::computeMaxi
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(eM
266 }
267 return wMax;
268 }
269
270 G4double PhaseSpaceRauboldLynch::computeMaxi
271 #ifndef INCLXX_IN_GEANT4_MODE
272 if(nParticles>=wMaxNP) {
273 INCL_WARN("The requested number of parti
274 }
275
276 if(availableEnergy>wMaxMasslessX[wMaxNE-1]
277 INCL_WARN("The requested available energ
278 }
279 #endif
280
281 // compute the maximum weight
282 const G4double logMassless = ((*wMaxMassle
283 const G4double reducedSqrtS = availableEne
284 const G4double correction = (*wMaxCorrecti
285 const G4double wMax = std::exp(correction*
286 if(wMax>0.)
287 return wMax;
288 else
289 return computeMaximumWeightNaive();
290 }
291
292 G4double PhaseSpaceRauboldLynch::computeWeig
293 // Generate nParticles-2 sorted random num
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()+nPart
299
300 // invariant masses
301 for(size_t i=0; i<nParticles; ++i)
302 invariantMasses[i] = rnd[i]*availableEne
303
304 // compute the CM momenta and the weight f
305 G4double weight=KinematicsUtils::momentumI
306 momentaCM[0] = weight;
307 for(size_t i=1; i<nParticles-1; ++i) {
308 G4double momentumCM;
309 // assert(invariantMasses[i+1]-invariantMasses
310 if(invariantMasses[i+1]-invariantMasses[
311 else momentumCM = KinematicsUtils::momen
312 momentaCM[i] = momentumCM;
313 weight *= momentumCM;
314 }
315
316 return weight;
317 }
318
319 void PhaseSpaceRauboldLynch::generateEvent(P
320 Particle *p = particles[0];
321 ThreeVector momentum = Random::normVector(
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[
336
337 const G4double iM = invariantMasses[i];
338 const G4double recoilE = std::sqrt(momen
339 boostV = -momentum/recoilE;
340 for(size_t j=0; j<=i; ++j)
341 particles[j]->boost(boostV);
342 }
343 }
344
345 G4double PhaseSpaceRauboldLynch::getMaxGener
346 return maxGeneratedWeight;
347 }
348 }
349