216{
218 static const G4double expxl = -expxu;
219
223
224 static const G4int numMul = 1200;
225 static const G4int numSec = 60;
226
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
245 first = false;
246 for (i=0; i<numMul; i++) protmul[i] = 0.0;
247 for (i=0; i<numSec; i++) protnorm[i] = 0.0;
248 counter = -1;
249 for (npos=0; npos<(numSec/3); npos++) {
250 for (nneg=
Imax(0,npos-2); nneg<=npos; nneg++) {
251 for (nzero=0; nzero<numSec/3; nzero++) {
252 if (++counter < numMul) {
253 nt = npos+nneg+nzero;
254 if ( (nt>0) && (nt<=numSec) ) {
255 protmul[counter] =
pmltpc(npos,nneg,nzero,nt,protb,c)
257 protnorm[nt-1] += protmul[counter];
258 }
259 }
260 }
261 }
262 }
263
264 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
265 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
266 counter = -1;
267 for (npos=0; npos<numSec/3; npos++) {
268 for (nneg=
Imax(0,npos-1); nneg<=(npos+1); nneg++) {
269 for (nzero=0; nzero<numSec/3; nzero++) {
270 if (++counter < numMul) {
271 nt = npos+nneg+nzero;
272 if ( (nt>0) && (nt<=numSec) ) {
273 neutmul[counter] =
pmltpc(npos,nneg,nzero,nt,neutb,c)
275 neutnorm[nt-1] += neutmul[counter];
276 }
277 }
278 }
279 }
280 }
281 for (i=0; i<numSec; i++) {
282 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
283 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
284 }
285 }
286
287
288
289
290 pv[0] = incidentParticle;
291 pv[1] = targetParticle;
292 vecLen = 2;
293
294 if (!inElastic) {
296 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
298 if (
G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) ) {
299
302 }
303 }
304 return;
305 } else if (availableEnergy <= pionMass) return;
306
307
308 npos = 0, nneg = 0, nzero = 0;
311
312 G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98};
314
315
317 if (targetCode == protonCode) {
318 w0 = -
sqr(1.+protb)/(2.*c*c);
319 wp = w0 = std::exp(w0);
321 npos = 0;
322 nneg = 0;
323 nzero = 1;
324 } else {
325 npos = 1;
326 nneg = 0;
327 nzero = 0; }
328 } else {
329 w0 = -
sqr(1.+neutb)/(2.*c*c);
330 w0 = std::exp(w0);
331 wp = w0/2.;
332 wm = -
sqr(-1.+neutb)/(2.*c*c);
333 wm = std::exp(wm)/2.;
334 wt = w0+wp+wm;
335 wp = w0+wp;
337 if( ran < w0/wt)
338 { npos = 0; nneg = 0; nzero = 1; }
339 else if( ran < wp/wt)
340 { npos = 1; nneg = 0; nzero = 0; }
341 else
342 { npos = 0; nneg = 1; nzero = 0; }
343 }
344 } else {
345
346
347 G4double aleab = std::log(availableEnergy);
348 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
349 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
350
351
352
354
355 for (nt=1; nt<=numSec; nt++) {
356 test = std::exp(
Amin( expxu,
Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
357 dum =
pi*nt/(2.0*
n*
n);
358 if (std::fabs(dum) < 1.0) {
359 if( test >= 1.0e-10 )anpn += dum*test;
360 } else {
361 anpn += dum*test;
362 }
363 }
364
367 if( targetCode == protonCode )
368 {
369 counter = -1;
370 for (npos=0; npos<numSec/3; npos++) {
371 for (nneg=
Imax(0,npos-2); nneg<=npos; nneg++) {
372 for (nzero=0; nzero<numSec/3; nzero++) {
373 if (++counter < numMul) {
374 nt = npos+nneg+nzero;
375 if ( (nt>0) && (nt<=numSec) ) {
376 test = std::exp(
Amin( expxu,
Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
377 dum = (
pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
378 if (std::fabs(dum) < 1.0) {
379 if( test >= 1.0e-10 )excs += dum*test;
380 } else {
381 excs += dum*test;
382 }
383 if (ran < excs) goto outOfLoop;
384 }
385 }
386 }
387 }
388 }
389
390
391 inElastic = false;
392 return;
393 }
394 else
395 {
396 counter = -1;
397 for (npos=0; npos<numSec/3; npos++) {
398 for (nneg=
Imax(0,npos-1); nneg<=(npos+1); nneg++) {
399 for (nzero=0; nzero<numSec/3; nzero++) {
400 if (++counter < numMul) {
401 nt = npos+nneg+nzero;
402 if ( (nt>=1) && (nt<=numSec) ) {
403 test = std::exp(
Amin( expxu,
Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
404 dum = (
pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
405 if (std::fabs(dum) < 1.0) {
406 if( test >= 1.0e-10 )excs += dum*test;
407 } else {
408 excs += dum*test;
409 }
410 if (ran < excs) goto outOfLoop;
411 }
412 }
413 }
414 }
415 }
416
417 inElastic = false;
418 return;
419 }
420 }
421 outOfLoop:
422
423 if( targetCode == protonCode)
424 {
425 if( npos == nneg)
426 {
427 }
428 else if (npos == (1+nneg))
429 {
431 {
433 }
434 else
435 {
437 }
438 }
439 else
440 {
443 }
444 }
445 else
446 {
447 if( npos == nneg)
448 {
450 {
453 }
454 else
455 {
456 }
457 }
458 else if ( npos == (1+nneg))
459 {
461 }
462 else
463 {
465 }
466 }
467
468
469 nt = npos + nneg + nzero;
470 while ( nt > 0)
471 {
474 {
475 if( npos > 0 )
477 npos--;
478 }
479 }
480 else if ( ran < (
G4double)(npos+nneg)/nt)
481 {
482 if( nneg > 0 )
483 {
485 nneg--;
486 }
487 }
488 else
489 {
490 if( nzero > 0 )
491 {
493 nzero--;
494 }
495 }
496 nt = npos + nneg + nzero;
497 }
499 {
500 G4cout <<
"Particles produced: " ;
503 for (i=2; i < vecLen; i++)
504 {
506 }
508 }
509 return;
510 }
G4double pmltpc(G4int np, G4int nm, G4int nz, G4int n, G4double b, G4double c)
G4double Amin(G4double a, G4double b)
G4double Amax(G4double a, G4double b)
G4int Imax(G4int a, G4int b)
G4double getTotalMomentum() const