204{
206 static const G4double expxl = -expxu;
207
211
212 static const G4int numMul = 1200;
213 static const G4int numSec = 60;
214
218
219 static G4bool first =
true;
220 static G4double protmul[numMul], protnorm[numSec];
221 static G4double neutmul[numMul], neutnorm[numSec];
222
223
224
225
226 G4int i, counter, nt, npos, nneg, nzero;
227
228 if (first) {
229
230 first = false;
231 for (i = 0; i < numMul; i++) protmul[i] = 0.0;
232 for (i = 0; i < numSec; i++) protnorm[i] = 0.0;
233 counter = -1;
234 for (npos = 0; npos < (numSec/3); npos++) {
235 for (nneg = std::max(0,npos-2); nneg <= npos; nneg++) {
236 for (nzero = 0; nzero < numSec/3; nzero++) {
237 if (++counter < numMul) {
238 nt = npos+nneg+nzero;
239 if ((nt>0) && (nt<=numSec) ) {
240 protmul[counter] =
pmltpc(npos,nneg,nzero,nt,protb,c);
241 protnorm[nt-1] += protmul[counter];
242 }
243 }
244 }
245 }
246 }
247
248 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
249 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
250 counter = -1;
251 for( npos=0; npos<numSec/3; npos++ )
252 {
253 for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ )
254 {
255 for( nzero=0; nzero<numSec/3; nzero++ )
256 {
257 if( ++counter < numMul )
258 {
259 nt = npos+nneg+nzero;
260 if( (nt>0) && (nt<=numSec) )
261 {
262 neutmul[counter] =
pmltpc(npos,nneg,nzero,nt,neutb,c);
263 neutnorm[nt-1] += neutmul[counter];
264 }
265 }
266 }
267 }
268 }
269 for (i = 0; i < numSec; i++) {
270 if (protnorm[i] > 0.0) protnorm[i] = 1.0/protnorm[i];
271 if (neutnorm[i] > 0.0) neutnorm[i] = 1.0/neutnorm[i];
272 }
273 }
274
275
276 pv[0] = incidentParticle;
277 pv[1] = targetParticle;
278 vecLen = 2;
279
280 if( !inElastic )
281 {
282 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
283 G4int iplab =
G4int( std::min( 9.0, incidentTotalMomentum*2.5 ) );
284 if(
G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
285 {
287 if( targetCode == protonCode)
288 {
289 if (ran < 0.2)
290 {
293 }
294 else if (ran < 0.4)
295 {
298 }
299 else if (ran < 0.6)
300 {
303 }
304 else if (ran < 0.8)
305 {
308 }
309 else
310 {
313 }
314 }
315 else
316 {
317 if (ran < 0.2)
318 {
321 }
322 else if (ran < 0.3)
323 {
326 }
327 else if (ran < 0.4)
328 {
331 }
332 else if (ran < 0.5)
333 {
336 }
337 else if (ran < 0.6)
338 {
341 }
342 else if (ran < 0.7)
343 {
346 }
347 else if (ran < 0.8)
348 {
351 }
352 else if (ran < 0.9)
353 {
356 }
357 else
358 {
361 }
362 }
363 }
364 return;
365 }
367 return;
368
369
370
371 npos = 0; nneg = 0; nzero = 0;
372
373
374
375 G4double aleab = std::log(availableEnergy);
376 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
377 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
378
379
380
382
383 for (nt=1; nt<=numSec; nt++) {
384 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
385 dum =
pi*nt/(2.0*
n*
n);
386 if (std::fabs(dum) < 1.0) {
387 if (test >= 1.0e-10) anpn += dum*test;
388 } else {
389 anpn += dum*test;
390 }
391 }
392
395 if( targetCode == protonCode )
396 {
397 counter = -1;
398 for (npos=0; npos<numSec/3; npos++) {
399 for (nneg=std::max(0,npos-2); nneg<=npos; nneg++) {
400 for (nzero=0; nzero<numSec/3; nzero++) {
401 if (++counter < numMul) {
402 nt = npos+nneg+nzero;
403 if ( (nt>0) && (nt<=numSec) ) {
404 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
405 dum = (
pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
406 if (std::fabs(dum) < 1.0) {
407 if (test >= 1.0e-10) excs += dum*test;
408 } else {
409 excs += dum*test;
410 }
411 if (ran < excs) goto outOfLoop;
412 }
413 }
414 }
415 }
416 }
417
418
419
420 inElastic = false;
421 return;
422 }
423 else
424 {
425 counter = -1;
426 for (npos=0; npos<numSec/3; npos++) {
427 for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
428 for (nzero=0; nzero<numSec/3; nzero++) {
429 if (++counter < numMul) {
430 nt = npos+nneg+nzero;
431 if ( (nt>=1) && (nt<=numSec) ) {
432 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
433 dum = (
pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
434 if (std::fabs(dum) < 1.0) {
435 if (test >= 1.0e-10) excs += dum*test;
436 } else {
437 excs += dum*test;
438 }
439 if (ran < excs) goto outOfLoop;
440 }
441 }
442 }
443 }
444 }
445
446
447 inElastic = false;
448 return;
449 }
450
451 outOfLoop:
452
453
454
455
456
458 if (targetCode == protonCode) {
459 if( npos == nneg)
460 {
461 }
462 else if (npos == (nneg+1))
463 {
464 if( ran < 0.50)
465 {
467 }
468 else
469 {
471 }
472 }
473 else
474 {
477 }
478 } else {
479 if (npos == nneg)
480 {
481 if (ran < 0.5)
482 {
483 }
484 else
485 {
488 }
489 }
490 else if (npos == (nneg-1))
491 {
493 }
494 else
495 {
497 }
498 }
499
500 nt = npos + nneg + nzero;
501 while (nt > 0) {
504 if (npos > 0) {
506 npos--;
507 }
508 }
else if (rnd < (
G4double)(npos+nneg)/nt) {
509 if (nneg > 0) {
511 nneg--;
512 }
513 } else {
514 if (nzero > 0) {
516 nzero--;
517 }
518 }
519 nt = npos + nneg + nzero;
520 }
521
523 G4cout <<
"Particles produced: " ;
526 for ( i = 2; i < vecLen; i++)
G4cout << pv[i].getCode() <<
" " ;
528 }
529
530 return;
531}
G4double pmltpc(G4int np, G4int nm, G4int nz, G4int n, G4double b, G4double c)
G4double getTotalMomentum() const