204{
206 static const G4double expxl = -expxu;
207
211
212 static const G4int numMul = 1200;
213 static const G4int numSec = 60;
214
217
220
221 static G4bool first =
true;
222 static G4double protmul[numMul], protnorm[numSec];
223 static G4double neutmul[numMul], neutnorm[numSec];
224
225
226
227
228 G4int i, counter, nt, npos, nneg, nzero;
229
230 if (first) {
231
232 first = false;
233 for (i = 0; i < numMul; i++) protmul[i] = 0.0;
234 for (i = 0; i < numSec; i++) protnorm[i] = 0.0;
235 counter = -1;
236 for( npos=0; npos<(numSec/3); npos++ )
237 {
238 for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ )
239 {
240 for( nzero=0; nzero<numSec/3; nzero++ )
241 {
242 if( ++counter < numMul )
243 {
244 nt = npos+nneg+nzero;
245 if( (nt>0) && (nt<=numSec) )
246 {
247 protmul[counter] =
pmltpc(npos,nneg,nzero,nt,protb,c);
248 protnorm[nt-1] += protmul[counter];
249 }
250 }
251 }
252 }
253 }
254 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
255 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
256 counter = -1;
257 for( npos=0; npos<numSec/3; npos++ )
258 {
259 for( nneg=npos; nneg<=(npos+2); nneg++ )
260 {
261 for( nzero=0; nzero<numSec/3; nzero++ )
262 {
263 if( ++counter < numMul )
264 {
265 nt = npos+nneg+nzero;
266 if( (nt>0) && (nt<=numSec) )
267 {
268 neutmul[counter] =
pmltpc(npos,nneg,nzero,nt,neutb,c);
269 neutnorm[nt-1] += neutmul[counter];
270 }
271 }
272 }
273 }
274 }
275 for( i=0; i<numSec; i++ )
276 {
277 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
278 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
279 }
280 }
281
282
283
284
285 pv[0] = incidentParticle;
286 pv[1] = targetParticle;
287 vecLen = 2;
288
289 if( !inElastic )
290 {
291 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
292 G4int iplab =
G4int( std::min( 9.0, incidentTotalMomentum*2.5 ) );
293 if(
G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
294 {
297 {
298 if (ran < 0.2)
299 {
302 }
303 else if (ran < 0.4)
304 {
307 }
308 else if (ran < 0.6)
309 {
312 }
313 else if (ran < 0.8)
314 {
317 }
318 else
319 {
322 }
323 }
324 else
325 {
326 if (ran < 0.2)
327 {
330 }
331 else if (ran < 0.3)
332 {
335 }
336 else if (ran < 0.4)
337 {
340 }
341 else if (ran < 0.5)
342 {
345 }
346 else if (ran < 0.7)
347 {
350 }
351 else if (ran < 0.8)
352 {
355 }
356 else
357 {
360 }
361 }
362 }
363 return;
364 }
366 return;
367
368
369
370 npos = 0; nneg = 0; nzero = 0;
371
372
373
374 G4double aleab = std::log(availableEnergy);
375 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
376 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
377
378
379
381
382 for (nt=1; nt<=numSec; nt++) {
383 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
384 dum =
pi*nt/(2.0*
n*
n);
385 if (std::fabs(dum) < 1.0) {
386 if (test >= 1.0e-10) anpn += dum*test;
387 } else {
388 anpn += dum*test;
389 }
390 }
391
394 if( targetCode == protonCode )
395 {
396 counter = -1;
397 for (npos=0; npos<numSec/3; npos++) {
398 for (nneg=std::max(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>0) && (nt<=numSec) ) {
403 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
404 dum = (
pi/anpn)*nt*protmul[counter]*protnorm[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
418 inElastic = false;
419 return;
420 }
421 else
422 {
423 counter = -1;
424 for (npos=0; npos<numSec/3; npos++) {
425 for (nneg=npos; nneg<=(npos+2); nneg++) {
426 for( nzero=0; nzero<numSec/3; nzero++) {
427 if (++counter < numMul) {
428 nt = npos+nneg+nzero;
429 if ( (nt>=1) && (nt<=numSec) ) {
430 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
431 dum = (
pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
432 if (std::fabs(dum) < 1.0) {
433 if( test >= 1.0e-10 )excs += dum*test;
434 } else {
435 excs += dum*test;
436 }
437 if (ran < excs) goto outOfLoop;
438 }
439 }
440 }
441 }
442 }
443
444 inElastic = false;
445 return;
446 }
447
448 outOfLoop:
449
450
451
452
453
456 if( npos == nneg)
457 {
458 }
459 else if (npos == (nneg-1))
460 {
461 if( ran < 0.50)
462 {
464 }
465 else
466 {
468 }
469 }
470 else
471 {
474 }
475 } else {
476 if (npos == nneg)
477 {
478 if (ran < 0.5)
479 {
480 }
481 else
482 {
485 }
486 }
487 else if (npos == (nneg+1))
488 {
490 }
491 else
492 {
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