206{
208 static const G4double expxl = -expxu;
209
213
214 static const G4int numMul = 1200;
215 static const G4int numSec = 60;
216
220
223
224 static G4bool first =
true;
225 static G4double protmul[numMul], protnorm[numSec];
226 static G4double neutmul[numMul], neutnorm[numSec];
227
228
229
230
231 G4int i, counter, nt, npos, nneg, nzero;
232
233 if( first )
234 {
235 first = false;
236 for( i=0; i<numMul; i++ )protmul[i] = 0.0;
237 for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
238 counter = -1;
239 for( npos=0; npos<(numSec/3); npos++ )
240 {
241 for( nneg=
Imax(0,npos-2); nneg<=npos; nneg++ )
242 {
243 for( nzero=0; nzero<numSec/3; nzero++ )
244 {
245 if( ++counter < numMul )
246 {
247 nt = npos+nneg+nzero;
248 if( (nt>0) && (nt<=numSec) )
249 {
250 protmul[counter] =
251 pmltpc(npos,nneg,nzero,nt,protb,c) ;
252 protnorm[nt-1] += protmul[counter];
253 }
254 }
255 }
256 }
257 }
258 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
259 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
260 counter = -1;
261 for( npos=0; npos<numSec/3; npos++ )
262 {
263 for( nneg=
Imax(0,npos-1); nneg<=(npos+1); nneg++ )
264 {
265 for( nzero=0; nzero<numSec/3; nzero++ )
266 {
267 if( ++counter < numMul )
268 {
269 nt = npos+nneg+nzero;
270 if( (nt>0) && (nt<=numSec) )
271 {
272 neutmul[counter] =
273 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
288
289
290 pv[0] = incidentParticle;
291 pv[1] = targetParticle;
292 vecLen = 2;
293
294 if( !inElastic )
295 {
297 {
298 G4double cech[] = {0.33,0.27,0.29,0.31,0.27,0.18,0.13,0.10,0.09,0.07};
300 if(
G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
301 {
304 }
305 }
306 return;
307 }
308 else if (availableEnergy <= pionMass)
309 return;
310
311
312
313 npos = 0, nneg = 0, nzero = 0;
316
317 G4double supp[] = {0., 0.2, 0.45, 0.55, 0.65, 0.75, 0.85, 0.90, 0.94, 0.98};
319 {
320
321
323 if( targetCode == protonCode )
324 {
325 w0 = -
sqr(1.+protb)/(2.*c*c);
326 wp = w0 = std::exp(w0);
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 for (nneg=
Imax(0,npos-2); nneg<=npos; nneg++) {
378 for (nzero=0; nzero<numSec/3; nzero++) {
379 if (++counter < numMul) {
380 nt = npos+nneg+nzero;
381 if ( (nt>0) && (nt<=numSec) ) {
382 test = std::exp(
Amin( expxu,
Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
383 dum = (
pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
384 if (std::fabs(dum) < 1.0) {
385 if( test >= 1.0e-10 )excs += dum*test;
386 } else {
387 excs += dum*test;
388 }
389 if (ran < excs) goto outOfLoop;
390 }
391 }
392 }
393 }
394 }
395
396
397 inElastic = false;
398 return;
399 }
400 else
401 {
402 counter = -1;
403 for (npos=0; npos<numSec/3; npos++) {
404 for (nneg=
Imax(0,npos-1); nneg<=(npos+1); nneg++) {
405 for (nzero=0; nzero<numSec/3; nzero++) {
406 if (++counter < numMul) {
407 nt = npos+nneg+nzero;
408 if ( (nt>=1) && (nt<=numSec) ) {
409 test = std::exp(
Amin( expxu,
Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
410 dum = (
pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
411 if (std::fabs(dum) < 1.0) {
412 if( test >= 1.0e-10 )excs += dum*test;
413 } else {
414 excs += dum*test;
415 }
416 if (ran < excs) goto outOfLoop;
417 }
418 }
419 }
420 }
421 }
422
423 inElastic = false;
424 return;
425 }
426 }
427 outOfLoop:
428
429 if( targetCode == protonCode)
430 {
431 if( npos == nneg)
432 {
433 }
434 else if (npos == (1+nneg))
435 {
437 {
439 }
440 else
441 {
443 }
444 }
445 else
446 {
449 }
450 }
451 else
452 {
453 if( npos == nneg)
454 {
456 {
459 }
460 else
461 {
462 }
463 }
464 else if ( npos == (1+nneg))
465 {
467 }
468 else
469 {
471 }
472 }
473
474
475 nt = npos + nneg + nzero;
476 while ( nt > 0)
477 {
480 {
481 if( npos > 0 )
483 npos--;
484 }
485 }
486 else if ( ran < (
G4double)(npos+nneg)/nt)
487 {
488 if( nneg > 0 )
489 {
491 nneg--;
492 }
493 }
494 else
495 {
496 if( nzero > 0 )
497 {
499 nzero--;
500 }
501 }
502 nt = npos + nneg + nzero;
503 }
505 {
506 G4cout <<
"Particles produced: " ;
509 for (i=2; i < vecLen; i++)
510 {
512 }
514 }
515 return;
516 }
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