218{
220 static const G4double expxl = -expxu;
221
225
226 static const G4int numMul = 1200;
227 static const G4int numSec = 60;
228
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 {
251 for( nneg=std::max(0,npos-2); nneg<=npos; nneg++ )
252 {
253 for( nzero=0; nzero<numSec/3; nzero++ )
254 {
255 if( ++counter < numMul )
256 {
257 nt = npos+nneg+nzero;
258 if( (nt>0) && (nt<=numSec) )
259 {
260 protmul[counter] =
pmltpc(npos,nneg,nzero,nt,protb,c);
261 protnorm[nt-1] += protmul[counter];
262 }
263 }
264 }
265 }
266 }
267 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
268 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
269 counter = -1;
270 for( npos=0; npos<numSec/3; npos++ )
271 {
272 for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ )
273 {
274 for( nzero=0; nzero<numSec/3; nzero++ )
275 {
276 if( ++counter < numMul )
277 {
278 nt = npos+nneg+nzero;
279 if( (nt>0) && (nt<=numSec) )
280 {
281 neutmul[counter] =
pmltpc(npos,nneg,nzero,nt,neutb,c);
282 neutnorm[nt-1] += neutmul[counter];
283 }
284 }
285 }
286 }
287 }
288 for( i=0; i<numSec; i++ )
289 {
290 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
291 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
292 }
293 }
294
295
296
297
298 pv[0] = incidentParticle;
299 pv[1] = targetParticle;
300 vecLen = 2;
301
302 if( !inElastic )
303 {
304 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
305 G4int iplab =
G4int( std::min( 9.0, incidentTotalMomentum*2.5 ) );
306 if(
G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
307 {
309 if( targetCode == protonCode)
310 {
313 }
314 else
315 {
316 if(ran < 0.2)
317 {
320 }
321 else if(ran < 0.4)
322 {
325 }
326 else if(ran < 0.6)
327 {
330 }
331 else if(ran < 0.8)
332 {
335 }
336 else
337 {
340 }
341 }
342 }
343 return;
344 }
346 return;
347
348
349
350 npos = 0; nneg = 0; nzero = 0;
351
352
353
354 G4double aleab = std::log(availableEnergy);
355 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
356 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
357
358
359
361
362 for (nt=1; nt<=numSec; nt++) {
363 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
364 dum =
pi*nt/(2.0*
n*
n);
365 if (std::fabs(dum) < 1.0) {
366 if (test >= 1.0e-10) anpn += dum*test;
367 } else {
368 anpn += dum*test;
369 }
370 }
371
374 if( targetCode == protonCode )
375 {
376 counter = -1;
377 for (npos=0; npos<numSec/3; npos++) {
378 for (nneg=std::max(0,npos-2); nneg<=npos; nneg++) {
379 for (nzero=0; nzero<numSec/3; nzero++) {
380 if (++counter < numMul) {
381 nt = npos+nneg+nzero;
382 if ( (nt>0) && (nt<=numSec) ) {
383 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
384 dum = (
pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
385 if (std::fabs(dum) < 1.0) {
386 if (test >= 1.0e-10) excs += dum*test;
387 } else {
388 excs += dum*test;
389 }
390 if (ran < excs) goto outOfLoop;
391 }
392 }
393 }
394 }
395 }
396
397
398 inElastic = false;
399 return;
400 }
401 else
402 {
403 counter = -1;
404 for (npos=0; npos<numSec/3; npos++) {
405 for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
406 for (nzero=0; nzero<numSec/3; nzero++) {
407 if (++counter < numMul) {
408 nt = npos+nneg+nzero;
409 if ( (nt>=1) && (nt<=numSec) ) {
410 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
411 dum = (
pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
412 if (std::fabs(dum) < 1.0) {
413 if( test >= 1.0e-10 )excs += dum*test;
414 } else {
415 excs += dum*test;
416 }
417 if (ran < excs) goto outOfLoop;
418 }
419 }
420 }
421 }
422 }
423
424 inElastic = false;
425 return;
426 }
427
428 outOfLoop:
429
431 if (targetCode == protonCode) {
432 if( npos == nneg)
433 {
434 }
435 else if (npos == (nneg+1))
436 {
437 if( ran < 0.25)
438 {
440 }
441 else if(ran < 0.5)
442 {
444 }
445 else
446 {
448 }
449 }
450 else
451 {
452 if(ran < 0.5)
453 {
456 }
457 else
458 {
461 }
462 }
463 } else {
464 if (npos == nneg)
465 {
466 if (ran < 0.5)
467 {
468 }
469 else if (ran < 0.75)
470 {
473 }
474 else
475 {
478 }
479 }
480 else if (npos == (nneg-1))
481 {
483 }
484 else
485 {
486 if (ran < 0.5)
487 {
489 }
490 else
491 {
493 }
494 }
495 }
496
497 nt = npos + nneg + nzero;
498 while (nt > 0) {
501 if (npos > 0) {
503 npos--;
504 }
505 }
else if (rnd < (
G4double)(npos+nneg)/nt) {
506 if (nneg > 0) {
508 nneg--;
509 }
510 } else {
511 if (nzero > 0) {
513 nzero--;
514 }
515 }
516 nt = npos + nneg + nzero;
517 }
518
520 G4cout <<
"Particles produced: " ;
523 for (i = 2; i < vecLen; i++)
G4cout << pv[i].getName() <<
" " ;
525 }
526
527 return;
528}
G4double pmltpc(G4int np, G4int nm, G4int nz, G4int n, G4double b, G4double c)
G4double getTotalMomentum() const