217{
219 static const G4double expxl = -expxu;
220
224
225 static const G4int numMul = 1200;
226 static const G4int numSec = 60;
227
230
233
234 static G4bool first =
true;
235 static G4double protmul[numMul], protnorm[numSec];
236 static G4double neutmul[numMul], neutnorm[numSec];
237
238
239
240
241 G4int i, counter, nt, npos, nneg, nzero;
242
243 if (first) {
244 first = false;
245 for( i=0; i<numMul; i++ )protmul[i] = 0.0;
246 for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
247 counter = -1;
248 for (npos = 0; npos < (numSec/3); npos++) {
249 for (nneg = std::max(0,npos-1); nneg <= (npos+1); nneg++) {
250 for (nzero = 0; nzero < numSec/3; nzero++) {
251 if (++counter < numMul) {
252 nt = npos+nneg+nzero;
253 if ((nt > 0) && (nt<=numSec) ) {
254 protmul[counter] =
pmltpc(npos,nneg,nzero,nt,protb,c);
255 protnorm[nt-1] += protmul[counter];
256 }
257 }
258 }
259 }
260 }
261
262 for (i = 0; i < numMul; i++) neutmul[i] = 0.0;
263 for (i = 0; i < numSec; i++) neutnorm[i] = 0.0;
264 counter = -1;
265 for (npos = 0; npos < numSec/3; npos++) {
266 for (nneg = npos; nneg <= (npos+2); nneg++) {
267 for (nzero = 0; nzero < numSec/3; nzero++) {
268 if (++counter < numMul) {
269 nt = npos+nneg+nzero;
270 if ((nt>0) && (nt<=numSec) ) {
271 neutmul[counter] =
pmltpc(npos,nneg,nzero,nt,neutb,c);
272 neutnorm[nt-1] += neutmul[counter];
273 }
274 }
275 }
276 }
277 }
278 for (i = 0; i < numSec; i++) {
279 if (protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
280 if (neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
281 }
282 }
283
284 pv[0] = incidentParticle;
285 pv[1] = targetParticle;
286 vecLen = 2;
287
288 if (!inElastic) {
289 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
290 G4int iplab =
G4int( std::min( 9.0, incidentTotalMomentum*2.5 ) );
291 if (
G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) ) {
296 } else {
297 if (ran < 0.2) {
300 } else if (ran < 0.4) {
303 } else if (ran < 0.6) {
306 } else if (ran < 0.8) {
309 } else {
312 }
313 }
314 }
315 return;
316
318
319
320 npos = 0, nneg = 0, nzero = 0;
321
322
323 G4double aleab = std::log(availableEnergy);
324 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
325 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
326
327
328
330
331 for (nt = 1; nt <= numSec; nt++) {
332 test = std::exp(std::min(expxu, std::max(expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
333 dum =
pi*nt/(2.0*
n*
n);
334 if (std::fabs(dum) < 1.0) {
335 if (test >= 1.0e-10) anpn += dum*test;
336 } else {
337 anpn += dum*test;
338 }
339 }
340
343 if (targetCode == protonCode) {
344 counter = -1;
345 for (npos=0; npos<numSec/3; npos++) {
346 for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
347 for (nzero=0; nzero<numSec/3; nzero++) {
348 if (++counter < numMul) {
349 nt = npos+nneg+nzero;
350 if ( (nt>0) && (nt<=numSec) ) {
351 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
352 dum = (
pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
353 if (std::fabs(dum) < 1.0) {
354 if( test >= 1.0e-10 )excs += dum*test;
355 } else {
356 excs += dum*test;
357 }
358 if (ran < excs) goto outOfLoop;
359 }
360 }
361 }
362 }
363 }
364
365
366 inElastic = false;
367 return;
368
369 } else {
370 counter = -1;
371 for (npos=0; npos<numSec/3; npos++) {
372 for (nneg=npos; nneg<=(npos+2); nneg++) {
373 for (nzero=0; nzero<numSec/3; nzero++) {
374 if (++counter < numMul) {
375 nt = npos+nneg+nzero;
376 if ( (nt>=1) && (nt<=numSec) ) {
377 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
378 dum = (
pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
379 if (std::fabs(dum) < 1.0) {
380 if( test >= 1.0e-10 )excs += dum*test;
381 } else {
382 excs += dum*test;
383 }
384 if (ran < excs) goto outOfLoop;
385 }
386 }
387 }
388 }
389 }
390
391
392 inElastic = false;
393 return;
394 }
395
396 outOfLoop:
397
400 if( npos == nneg)
401 {
402 }
403 else if (npos == (nneg-1))
404 {
405 if( ran < 0.25)
406 {
408 }
409 else if (ran < 0.5)
410 {
412 }
413 else
414 {
416 }
417 }
418 else
419 {
420 if(ran < 0.5)
421 {
424 }
425 else
426 {
429 }
430 }
431 } else {
432 if (npos == nneg)
433 {
434 if (ran < 0.5)
435 {
436 }
437 else if (ran < 0.75)
438 {
441 }
442 else
443 {
446 }
447 }
448 else if (npos == (nneg+1))
449 {
451 }
452 else
453 {
454 if (ran < 0.5)
455 {
457 }
458 else
459 {
461 }
462 }
463 }
464
465 nt = npos + nneg + nzero;
466 while (nt > 0) {
469 if (npos > 0) {
471 npos--;
472 }
473 }
else if (rnd < (
G4double)(npos+nneg)/nt) {
474 if (nneg > 0) {
476 nneg--;
477 }
478 } else {
479 if (nzero > 0) {
481 nzero--;
482 }
483 }
484 nt = npos + nneg + nzero;
485 }
486
488 G4cout <<
"Particles produced: " ;
491 for (i = 2; i < vecLen; i++)
G4cout << pv[i].getCode() <<
" " ;
493 }
494
495 return;
496}
G4double pmltpc(G4int np, G4int nm, G4int nz, G4int n, G4double b, G4double c)
G4double getTotalMomentum() const