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=
Imax(0,npos-2); nneg<=npos; 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=
Imax(0,npos-1); nneg<=(npos+1); 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 {
296 {
297 G4double cech[] = {0.33,0.27,0.29,0.31,0.27,0.18,0.13,0.10,0.09,0.07};
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
322 if( targetCode == protonCode )
323 {
324 w0 = -
sqr(1.+protb)/(2.*c*c);
325 wp = w0 = std::exp(w0);
326 wp *= 2.;
328 { npos = 0; nneg = 0; nzero = 1; }
329 else
330 { npos = 1; nneg = 0; nzero = 0; }
331 }
332 else
333 {
334 w0 = -
sqr(1.+neutb)/(2.*c*c);
335 wp = w0 = std::exp(w0);
336 wm = -
sqr(-1.+neutb)/(2.*c*c);
337 wm = std::exp(wm);
338 wt = w0+wp+wm;
339 wp = w0+wp;
341 if( ran < w0/wt)
342 { npos = 0; nneg = 0; nzero = 1; }
343 else if( ran < wp/wt)
344 { npos = 1; nneg = 0; nzero = 0; }
345 else
346 { npos = 0; nneg = 1; nzero = 0; }
347 }
348 }
349 else
350 {
351
352
353 G4double aleab = std::log(availableEnergy);
354 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
355 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
356
357
358
360
361 for (nt=1; nt<=numSec; nt++) {
362 test = std::exp(
Amin( expxu,
Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
363 dum =
pi*nt/(2.0*
n*
n);
364 if (std::fabs(dum) < 1.0) {
365 if( test >= 1.0e-10 )anpn += dum*test;
366 } else {
367 anpn += dum*test;
368 }
369 }
370
373 if( targetCode == protonCode )
374 {
375 counter = -1;
376 for( npos=0; npos<numSec/3; npos++ )
377 {
378 for( nneg=
Imax(0,npos-2); nneg<=npos; nneg++ )
379 {
380 for (nzero=0; nzero<numSec/3; nzero++) {
381 if (++counter < numMul) {
382 nt = npos+nneg+nzero;
383 if ( (nt>0) && (nt<=numSec) ) {
384 test = std::exp(
Amin( expxu,
Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
385 dum = (
pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
386 if (std::fabs(dum) < 1.0) {
387 if( test >= 1.0e-10 )excs += dum*test;
388 } else {
389 excs += dum*test;
390 }
391 if (ran < excs) goto outOfLoop;
392 }
393 }
394 }
395 }
396 }
397
398
399 inElastic = false;
400 return;
401 }
402 else
403 {
404 counter = -1;
405 for( npos=0; npos<numSec/3; npos++ )
406 {
407 for( nneg=
Imax(0,npos-1); nneg<=(npos+1); nneg++ )
408 {
409 for (nzero=0; nzero<numSec/3; nzero++) {
410 if (++counter < numMul) {
411 nt = npos+nneg+nzero;
412 if ( (nt>=1) && (nt<=numSec) ) {
413 test = std::exp(
Amin( expxu,
Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
414 dum = (
pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
415 if (std::fabs(dum) < 1.0) {
416 if( test >= 1.0e-10 )excs += dum*test;
417 } else {
418 excs += dum*test;
419 }
420 if (ran < excs) goto outOfLoop;
421 }
422 }
423 }
424 }
425 }
426
427 inElastic = false;
428 return;
429 }
430 }
431 outOfLoop:
432
433 if( targetCode == protonCode)
434 {
435 if( npos == nneg)
436 {
437 }
438 else if (npos == (1+nneg))
439 {
441 {
443 }
444 else
445 {
447 }
448 }
449 else
450 {
453 }
454 }
455 else
456 {
457 if( npos == nneg)
458 {
460 {
463 }
464 else
465 {
466 }
467 }
468 else if ( npos == (1+nneg))
469 {
471 }
472 else
473 {
475 }
476 }
477
478
479 nt = npos + nneg + nzero;
480 while ( nt > 0)
481 {
484 {
485 if( npos > 0 )
487 npos--;
488 }
489 }
490 else if ( ran < (
G4double)(npos+nneg)/nt)
491 {
492 if( nneg > 0 )
493 {
495 nneg--;
496 }
497 }
498 else
499 {
500 if( nzero > 0 )
501 {
503 nzero--;
504 }
505 }
506 nt = npos + nneg + nzero;
507 }
509 {
510 G4cout <<
"Particles produced: " ;
513 for (i=2; i < vecLen; i++)
514 {
516 }
518 }
519 return;
520 }
G4double pmltpc(G4int np, G4int nm, G4int nz, G4int n, G4double b, G4double c)
G4double Amin(G4double a, G4double b)
G4double Amax(G4double a, G4double b)
G4int Imax(G4int a, G4int b)
G4double getTotalMomentum() const