212{
214 static const G4double expxl = -expxu;
215
219
220 static const G4int numMul = 1200;
221 static const G4int numSec = 60;
222
225
228
229 static G4bool first =
true;
230 static G4double protmul[numMul], protnorm[numSec];
231 static G4double neutmul[numMul], neutnorm[numSec];
232
233
234
235
236 G4int i, counter, nt, npos, nneg, nzero;
237
238 if (first) {
239
240 first = false;
241 for (i = 0; i < numMul; i++) protmul[i] = 0.0;
242 for (i = 0; i < numSec; i++) protnorm[i] = 0.0;
243 counter = -1;
244 for (npos = 0; npos < (numSec/3); npos++) {
245 for (nneg = std::max(0,npos-1); nneg <= (npos+1); nneg++) {
246 for( nzero=0; nzero<numSec/3; nzero++ )
247 {
248 if( ++counter < numMul )
249 {
250 nt = npos+nneg+nzero;
251 if( (nt>0) && (nt<=numSec) )
252 {
253 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 {
267 for( nneg=npos; nneg<=(npos+2); nneg++ )
268 {
269 for( nzero=0; nzero<numSec/3; nzero++ )
270 {
271 if( ++counter < numMul )
272 {
273 nt = npos+nneg+nzero;
274 if( (nt>0) && (nt<=numSec) )
275 {
276 neutmul[counter] =
pmltpc(npos,nneg,nzero,nt,neutb,c)
278 neutnorm[nt-1] += neutmul[counter];
279 }
280 }
281 }
282 }
283 }
284 for( i=0; i<numSec; i++ )
285 {
286 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
287 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
288 }
289 }
290
291
292
293
294 pv[0] = incidentParticle;
295 pv[1] = targetParticle;
296 vecLen = 2;
297
298 if( !inElastic )
299 {
300 if( targetCode == protonCode )
301 {
302 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
303 G4int iplab =
G4int( std::min( 9.0, incidentTotalMomentum*2.5) );
304 if(
G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
305 {
308 }
309 }
310 return;
311 }
313 return;
314
315
316
317 npos = 0, nneg = 0, nzero = 0;
320
321 G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98};
323 {
324
325
328 {
329 w0 = -
sqr(1.+neutb)/(2.*c*c);
330 wm = w0 = std::exp(w0);
331 w0 = w0/2.;
333 { npos = 0; nneg = 0; nzero = 1; }
334 else
335 { npos = 0; nneg = 1; nzero = 0; }
336 }
337 else
338 {
339 w0 = -
sqr(1.+protb)/(2.*c*c);
340 w0 = std::exp(w0);
341 wp = w0/2.;
342 wm = -
sqr(-1.+protb)/(2.*c*c);
343 wm = std::exp(wm)/2.;
344 wt = w0+wp+wm;
345 wp = w0+wp;
347 if( ran < w0/wt)
348 { npos = 0; nneg = 0; nzero = 1; }
349 else if( ran < wp/wt)
350 { npos = 1; nneg = 0; nzero = 0; }
351 else
352 { npos = 0; nneg = 1; nzero = 0; }
353 }
354 }
355 else
356 {
357
358
359 G4double aleab = std::log(availableEnergy);
360 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
361 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
362
363
364
366
367 for (nt=1; nt<=numSec; nt++) {
368 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
369 dum =
pi*nt/(2.0*
n*
n);
370 if (std::fabs(dum) < 1.0) {
371 if( test >= 1.0e-10 )anpn += dum*test;
372 } else {
373 anpn += dum*test;
374 }
375 }
376
379 if( targetCode == protonCode )
380 {
381 counter = -1;
382 for (npos=0; npos<numSec/3; npos++) {
383 for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
384 for (nzero=0; nzero<numSec/3; nzero++) {
385 if (++counter < numMul) {
386 nt = npos+nneg+nzero;
387 if ( (nt>0) && (nt<=numSec) ) {
388 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
389 dum = (
pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
390 if (std::fabs(dum) < 1.0) {
391 if( test >= 1.0e-10 )excs += dum*test;
392 } else {
393 excs += dum*test;
394 }
395 if (ran < excs) goto outOfLoop;
396 }
397 }
398 }
399 }
400 }
401
402
403 inElastic = false;
404 return;
405 }
406 else
407 {
408 counter = -1;
409 for (npos=0; npos<numSec/3; npos++) {
410 for (nneg=npos; nneg<=(npos+2); nneg++) {
411 for (nzero=0; nzero<numSec/3; nzero++) {
412 if (++counter < numMul) {
413 nt = npos+nneg+nzero;
414 if ( (nt>=1) && (nt<=numSec) ) {
415 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
416 dum = (
pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
417 if (std::fabs(dum) < 1.0) {
418 if( test >= 1.0e-10 )excs += dum*test;
419 } else {
420 excs += dum*test;
421 }
422 if (ran < excs) goto outOfLoop;
423 }
424 }
425 }
426 }
427 }
428
429 inElastic = false;
430 return;
431 }
432 }
433 outOfLoop:
434
436 {
437 if( npos == nneg)
438 {
439 }
440 else if (npos == (nneg-1))
441 {
443 {
445 }
446 else
447 {
449 }
450 }
451 else
452 {
455 }
456 }
457 else
458 {
459 if( npos == nneg)
460 {
462 {
465 }
466 else
467 {
468 }
469 }
470 else if ( npos == (1+nneg))
471 {
473 }
474 else
475 {
477 }
478 }
479
480
481 nt = npos + nneg + nzero;
482 while ( nt > 0)
483 {
486 {
487 if( npos > 0 )
489 npos--;
490 }
491 }
492 else if ( ran < (
G4double)(npos+nneg)/nt)
493 {
494 if( nneg > 0 )
495 {
497 nneg--;
498 }
499 }
500 else
501 {
502 if( nzero > 0 )
503 {
505 nzero--;
506 }
507 }
508 nt = npos + nneg + nzero;
509 }
511 {
512 G4cout <<
"Particles produced: " ;
515 for (i=2; i < vecLen; i++)
516 {
518 }
520 }
521 return;
522 }
G4double pmltpc(G4int np, G4int nm, G4int nz, G4int n, G4double b, G4double c)
G4double getTotalMomentum() const