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