206{
208 static const G4double expxl = -expxu;
209
213
214 static const G4int numMul = 1200;
215 static const G4int numSec = 60;
216
219
222
223 static G4bool first =
true;
224 static G4double protmul[numMul], protnorm[numSec];
225 static G4double neutmul[numMul], neutnorm[numSec];
226
227
228
229
230 G4int i, counter, nt, npos, nneg, nzero;
231
232 if( first )
233 {
234 first = false;
235 for( i=0; i<numMul; i++ )protmul[i] = 0.0;
236 for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
237 counter = -1;
238 for( npos=0; npos<(numSec/3); npos++ )
239 {
240 for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ )
241 {
242 for( nzero=0; nzero<numSec/3; nzero++ )
243 {
244 if( ++counter < numMul )
245 {
246 nt = npos+nneg+nzero;
247 if( (nt>0) && (nt<=numSec) )
248 {
249 protmul[counter] =
250 pmltpc(npos,nneg,nzero,nt,protb,c) ;
251 protnorm[nt-1] += protmul[counter];
252 }
253 }
254 }
255 }
256 }
257 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
258 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
259 counter = -1;
260 for( npos=0; npos<numSec/3; npos++ )
261 {
262 for( nneg=npos; nneg<=(npos+2); nneg++ )
263 {
264 for( nzero=0; nzero<numSec/3; nzero++ )
265 {
266 if( ++counter < numMul )
267 {
268 nt = npos+nneg+nzero;
269 if( (nt>0) && (nt<=numSec) )
270 {
271 neutmul[counter] =
272 pmltpc(npos,nneg,nzero,nt,neutb,c);
273 neutnorm[nt-1] += neutmul[counter];
274 }
275 }
276 }
277 }
278 }
279 for( i=0; i<numSec; i++ )
280 {
281 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
282 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
283 }
284 }
285
286
287
288
289 pv[0] = incidentParticle;
290 pv[1] = targetParticle;
291 vecLen = 2;
292
293 if( !inElastic )
294 {
295 if( targetCode == protonCode )
296 {
297 G4double cech[] = {0.33,0.27,0.29,0.31,0.27,0.18,0.13,0.10,0.09,0.07};
298 G4int iplab =
G4int( std::min( 9.0, incidentTotalMomentum*5. ) );
299 if(
G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
300 {
303 }
304 }
305 return;
306 }
308 return;
309
310
311
312 npos = 0, nneg = 0, nzero = 0;
315
316 G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98};
318 {
319
320
323 {
324 w0 = -
sqr(1.+protb)/(2.*c*c);
325 w0 = std::exp(w0);
326 wm = -
sqr(-1.+protb)/(2.*c*c);
327 wm = std::exp(wm);
328 w0 = w0/2.;
329 wm = wm*1.5;
330 if(
G4UniformRand() < w0/(w0+wm) ) { npos = 0; nneg = 0; nzero = 1; }
331 else
332 { npos = 0; nneg = 1; nzero = 0; }
333 }
334 else
335 {
336 w0 = -
sqr(1.+neutb)/(2.*c*c);
337 wp = w0 = std::exp(w0);
338 wm = -
sqr(-1.+neutb)/(2.*c*c);
339 wm = std::exp(wm);
340 wt = w0+wp+wm;
341 wp = w0+wp;
343 if( ran < w0/wt)
344 { npos = 0; nneg = 0; nzero = 1; }
345 else if( ran < wp/wt)
346 { npos = 1; nneg = 0; nzero = 0; }
347 else
348 { npos = 0; nneg = 1; nzero = 0; }
349 }
350 }
351 else
352 {
353
354
355 G4double aleab = std::log(availableEnergy);
356 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
357 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
358
359
360
362
363 for (nt=1; nt<=numSec; nt++) {
364 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
365 dum =
pi*nt/(2.0*
n*
n);
366 if (std::fabs(dum) < 1.0) {
367 if( test >= 1.0e-10 )anpn += dum*test;
368 } else {
369 anpn += dum*test;
370 }
371 }
372
375 if( targetCode == protonCode )
376 {
377 counter = -1;
378 for( npos=0; npos<numSec/3; npos++ )
379 {
380 for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ )
381 {
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 {
406 counter = -1;
407 for( npos=0; npos<numSec/3; npos++ )
408 {
409 for( nneg=npos; nneg<=(npos+2); nneg++ )
410 {
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 == (nneg+1))
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; }
G4double pmltpc(G4int np, G4int nm, G4int nz, G4int n, G4double b, G4double c)
G4double getTotalMomentum() const