205{
207 static const G4double expxl = -expxu;
208
212
213 static const G4int numMul = 1200;
214 static const G4int numSec = 60;
215
217
220
221 static G4bool first =
true;
222 static G4double protmul[numMul], protnorm[numSec];
223 static G4double neutmul[numMul], neutnorm[numSec];
224
225
226
227
228 G4int i, counter, nt, npos, nneg, nzero;
229
230 if (first) {
231
232 first = false;
233 for( i=0; i<numMul; i++ )protmul[i] = 0.0;
234 for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
235 counter = -1;
236 for( npos=0; npos<(numSec/3); npos++ )
237 {
238 for( nneg=std::max(0,npos-1); nneg<=npos; nneg++ )
239 {
240 for( nzero=0; nzero<numSec/3; nzero++ )
241 {
242 if( ++counter < numMul )
243 {
244 nt = npos+nneg+nzero;
245 if( (nt>0) && (nt<=numSec) )
246 {
247 protmul[counter] =
pmltpc(npos,nneg,nzero,nt,protb,c);
248 protnorm[nt-1] += protmul[counter];
249 }
250 }
251 }
252 }
253 }
254 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
255 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
256 counter = -1;
257 for( npos=0; npos<numSec/3; npos++ )
258 {
259 for( nneg=npos; nneg<=(npos+1); nneg++ )
260 {
261 for( nzero=0; nzero<numSec/3; nzero++ )
262 {
263 if( ++counter < numMul )
264 {
265 nt = npos+nneg+nzero;
266 if( (nt>0) && (nt<=numSec) )
267 {
268 neutmul[counter] =
pmltpc(npos,nneg,nzero,nt,neutb,c);
269 neutnorm[nt-1] += neutmul[counter];
270 }
271 }
272 }
273 }
274 }
275 for( i=0; i<numSec; i++ )
276 {
277 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
278 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
279 }
280 }
281
282
283
284
285 pv[0] = incidentParticle;
286 pv[1] = targetParticle;
287 vecLen = 2;
288
289 if( !inElastic )
290 {
291 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
292 G4int iplab =
G4int( std::min( 9.0, incidentTotalMomentum*2.5 ) );
293 if(
G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
294 {
296 if( targetCode == protonCode)
297 {
298 if (ran < 0.2)
299 {
302 }
303 else if (ran < 0.4)
304 {
307 }
308 else if (ran < 0.6)
309 {
312 }
313 else if (ran < 0.8)
314 {
317 }
318 else
319 {
322 }
323 }
324 else
325 {
326 if (ran < 0.2)
327 {
330 }
331 else if (ran < 0.4)
332 {
335 }
336 else if (ran < 0.6)
337 {
340 }
341 else if (ran < 0.8)
342 {
345 }
346 else
347 {
350 }
351 }
352 }
353 return;
354 }
356 return;
357
358
359 npos = 0; nneg = 0; nzero = 0;
360
361
362 G4double aleab = std::log(availableEnergy);
363 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
364 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
365
366
367
369
370 for (nt=1; nt<=numSec; nt++) {
371 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
372 dum =
pi*nt/(2.0*
n*
n);
373 if (std::fabs(dum) < 1.0) {
374 if( test >= 1.0e-10 )anpn += dum*test;
375 } else {
376 anpn += dum*test;
377 }
378 }
379
382 if (targetCode == protonCode) {
383 counter = -1;
384 for (npos = 0; npos < numSec/3; npos++) {
385 for( nneg=std::max(0,npos-1); nneg<=npos; nneg++ )
386 {
387 for( nzero=0; nzero<numSec/3; nzero++ )
388 {
389 if( ++counter < numMul )
390 {
391 nt = npos+nneg+nzero;
392 if ( (nt>0) && (nt<=numSec) ) {
393 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
394 dum = (
pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
395 if (std::fabs(dum) < 1.0) {
396 if( test >= 1.0e-10 )excs += dum*test;
397 } else {
398 excs += dum*test;
399 }
400 if (ran < excs) goto outOfLoop;
401 }
402 }
403 }
404 }
405 }
406
407
408 inElastic = false;
409 return;
410 }
411 else
412 {
413 counter = -1;
414 for( npos=0; npos<numSec/3; npos++ )
415 {
416 for( nneg=npos; nneg<=(npos+1); nneg++ )
417 {
418 for( nzero=0; nzero<numSec/3; nzero++ )
419 {
420 if( ++counter < numMul )
421 {
422 nt = npos+nneg+nzero;
423 if ( (nt>=1) && (nt<=numSec) ) {
424 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
425 dum = (
pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
426 if (std::fabs(dum) < 1.0) {
427 if( test >= 1.0e-10 )excs += dum*test;
428 } else {
429 excs += dum*test;
430 }
431 if (ran < excs) goto outOfLoop;
432 }
433 }
434 }
435 }
436 }
437
438 inElastic = false;
439 return;
440 }
441
442 outOfLoop:
443
444
445
446
448 if (targetCode == protonCode) {
449 if( npos == nneg)
450 {
451 }
452 else
453 {
455 }
456 } else {
457 if (npos == nneg)
458 {
459 }
460 else
461 {
463 }
464 }
465
466 nt = npos + nneg + nzero;
467 while (nt > 0) {
470 if (npos > 0) {
472 npos--;
473 }
474 }
else if (rnd < (
G4double)(npos+nneg)/nt) {
475 if (nneg > 0) {
477 nneg--;
478 }
479 } else {
480 if (nzero > 0) {
482 nzero--;
483 }
484 }
485 nt = npos + nneg + nzero;
486 }
487
489 G4cout <<
"Particles produced: " ;
492 for (i=2; i < vecLen; i++)
G4cout << pv[i].getName() <<
" " ;
494 }
495 return;
496}
G4double pmltpc(G4int np, G4int nm, G4int nz, G4int n, G4double b, G4double c)
G4double getTotalMomentum() const