217{
219 static const G4double expxl = -expxu;
220
224
225 static const G4int numMul = 1200;
226 static const G4int numSec = 60;
227
231
232 static G4bool first =
true;
233 static G4double protmul[numMul], protnorm[numSec];
234 static G4double neutmul[numMul], neutnorm[numSec];
235
236
237
238
239 G4int i, counter, nt, npos, nneg, nzero;
240
241 if (first) {
242 first = false;
243 for( i=0; i<numMul; i++ )protmul[i] = 0.0;
244 for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
245 counter = -1;
246 for (npos = 0; npos < (numSec/3); npos++) {
247 for (nneg = std::max(0,npos-2); nneg <= (npos+1); nneg++) {
248 for (nzero = 0; nzero < numSec/3; nzero++) {
249 if (++counter < numMul) {
250 nt = npos+nneg+nzero;
251 if ((nt>0) && (nt<=numSec) ) {
252 protmul[counter] =
pmltpc(npos,nneg,nzero,nt,protb,c);
253 protnorm[nt-1] += protmul[counter];
254 }
255 }
256 }
257 }
258 }
259
260 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
261 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
262 counter = -1;
263 for (npos = 0; npos < numSec/3; npos++) {
264 for( nneg=std::max(0,npos-1); nneg<=(npos+2); nneg++ )
265 {
266 for( nzero=0; nzero<numSec/3; nzero++ )
267 {
268 if( ++counter < numMul )
269 {
270 nt = npos+nneg+nzero;
271 if( (nt>0) && (nt<=numSec) )
272 {
273 neutmul[counter] =
pmltpc(npos,nneg,nzero,nt,neutb,c);
274 neutnorm[nt-1] += neutmul[counter];
275 }
276 }
277 }
278 }
279 }
280 for( i=0; i<numSec; i++ )
281 {
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 pv[0] = incidentParticle;
288 pv[1] = targetParticle;
289 vecLen = 2;
290
291 if (!inElastic) {
292 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
293 G4int iplab =
G4int( std::min( 9.0, incidentTotalMomentum*2.5 ) );
294 if (
G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) ) {
296 if (targetCode == protonCode) {
297 if( ran < 0.2)
298 {
301 }
302 else if(ran < 0.4)
303 {
305 }
306 else if(ran < 0.6)
307 {
310 }
311 else if(ran < 0.8)
312 {
315 }
316 else
317 {
320 }
321 } else {
322 if(ran < 0.2)
323 {
325 }
326 else if(ran < 0.4)
327 {
330 }
331 else if(ran < 0.6)
332 {
335 }
336 else if(ran < 0.8)
337 {
340 }
341 else
342 {
345 }
346 }
347 }
348 return;
349 }
351 return;
352
353
354 npos = 0; nneg = 0; nzero = 0;
355
356
357
358 G4double aleab = std::log(availableEnergy);
359 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
360 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
361
362
363
365
366 for (nt=1; nt<=numSec; nt++) {
367 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
368 dum =
pi*nt/(2.0*
n*
n);
369 if (std::fabs(dum) < 1.0) {
370 if( test >= 1.0e-10 )anpn += dum*test;
371 } else {
372 anpn += dum*test;
373 }
374 }
375
378 if (targetCode == protonCode) {
379 counter = -1;
380 for (npos = 0; npos < numSec/3; npos++) {
381 for (nneg=std::max(0,npos-2); nneg<=(npos+1); nneg++) {
382 for (nzero=0; nzero<numSec/3; nzero++) {
383 if (++counter < numMul) {
384 nt = npos+nneg+nzero;
385 if ( (nt>0) && (nt<=numSec) ) {
386 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
387 dum = (
pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
388 if (std::fabs(dum) < 1.0) {
389 if( test >= 1.0e-10 )excs += dum*test;
390 } else {
391 excs += dum*test;
392 }
393 if (ran < excs) goto outOfLoop;
394 }
395 }
396 }
397 }
398 }
399
400
401 inElastic = false;
402 return;
403
404 } else {
405 counter = -1;
406 for (npos=0; npos<numSec/3; npos++) {
407 for (nneg=std::max(0,npos-1); nneg<=(npos+2); nneg++) {
408 for (nzero=0; nzero<numSec/3; nzero++) {
409 if (++counter < numMul) {
410 nt = npos+nneg+nzero;
411 if ( (nt>=1) && (nt<=numSec) ) {
412 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
413 dum = (
pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
414 if (std::fabs(dum) < 1.0) {
415 if( test >= 1.0e-10 )excs += dum*test;
416 } else {
417 excs += dum*test;
418 }
419 if (ran < excs) goto outOfLoop;
420 }
421 }
422 }
423 }
424 }
425
426 inElastic = false;
427 return;
428 }
429
430 outOfLoop:
431
433 if (targetCode == protonCode) {
434 if (npos == nneg) {
435 if (ran < 0.25)
436 {
437 }
438 else if(ran < 0.5)
439 {
441 }
442 else
443 {
446 }
447 } else if (npos == (nneg+1)) {
449 {
451 }
452 else if(ran < 0.5)
453 {
456 }
457 else
458 {
460 }
461 } else if (npos == (nneg-1)) {
463 } else {
466 }
467
468 } else {
469 if (npos == nneg) {
470 if(ran < 0.5)
471 {
472 }
473 else
474 {
477 }
478 } else if (npos == (nneg-1)) {
479 if( ran < 0.25)
480 {
482 }
483 else if(ran < 0.5)
484 {
487 }
488 else
489 {
491 }
492 } else if (npos == (1+nneg)) {
494 } else {
497 }
498 }
499
500 nt = npos + nneg + nzero;
501 while (nt > 0) {
504 if (npos > 0) {
506 npos--;
507 }
508 }
else if (rnd < (
G4double)(npos+nneg)/nt) {
509 if (nneg > 0) {
511 nneg--;
512 }
513 } else {
514 if (nzero > 0) {
516 nzero--;
517 }
518 }
519 nt = npos + nneg + nzero;
520 }
521
523 G4cout <<
"Particles produced: " ;
526 for (i = 2; i < vecLen; i++)
G4cout << pv[i].getName() <<
" " ;
528 }
529 return;
530}
G4double pmltpc(G4int np, G4int nm, G4int nz, G4int n, G4double b, G4double c)
G4double getTotalMomentum() const