215{
217 static const G4double expxl = -expxu;
218
222
223 static const G4int numMul = 1200;
224 static const G4int numSec = 60;
225
229
230 static G4bool first =
true;
231 static G4double protmul[numMul], protnorm[numSec];
232 static G4double neutmul[numMul], neutnorm[numSec];
233
234
235
236
237 G4int i, counter, nt, npos, nneg, nero;
238
239 if (first) {
240
241 first = false;
242 for( i=0; i<numMul; i++ )protmul[i] = 0.0;
243 for( i=0; i<numSec; i++ )protnorm[i] = 0.0;
244 counter = -1;
245 for( npos=0; npos<(numSec/3); npos++ )
246 {
247 for( nneg=
Imax(0,npos-1); nneg<=(npos+1); nneg++ )
248 {
249 for( nero=0; nero<numSec/3; nero++ )
250 {
251 if( ++counter < numMul )
252 {
253 nt = npos+nneg+nero;
254 if( (nt>0) && (nt<=numSec) )
255 {
256 protmul[counter] =
257 pmltpc(npos,nneg,nero,nt,protb,c) ;
258 protnorm[nt-1] += protmul[counter];
259 }
260 }
261 }
262 }
263 }
264 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
265 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
266 counter = -1;
267 for( npos=0; npos<numSec/3; npos++ )
268 {
269 for( nneg=npos; nneg<=(npos+2); nneg++ )
270 {
271 for( nero=0; nero<numSec/3; nero++ )
272 {
273 if( ++counter < numMul )
274 {
275 nt = npos+nneg+nero;
276 if( (nt>0) && (nt<=numSec) )
277 {
278 neutmul[counter] =
279 pmltpc(npos,nneg,nero,nt,neutb,c);
280 neutnorm[nt-1] += neutmul[counter];
281 }
282 }
283 }
284 }
285 }
286 for( i=0; i<numSec; i++ )
287 {
288 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
289 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
290 }
291 }
292
293
294
295
296 pv[0] = incidentParticle;
297 pv[1] = targetParticle;
298 vecLen = 2;
299
301 return;
302
303
304
305
306 npos = 0, nneg = 0, nero = 0;
309
310 G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98};
312 {
313
314
315 G4double cech[] = {1., 0.95, 0.79, 0.32, 0.19, 0.16, 0.14, 0.12, 0.10, 0.08};
318 {
319 if( targetCode == protonCode )
320 {
323 return;
324 }
325 else
326 {
327 return;
328 }
329 }
330
332 if( targetCode == protonCode )
333 {
334 w0 = -
sqr(1.+protb)/(2.*c*c);
335 wp = w0 = std::exp(w0);
336 wm = -
sqr(-1.+protb)/(2.*c*c);
337 wm = std::exp(wm);
338 wp *= 10.;
339 wt = w0+wp+wm;
340 wp = w0+wp;
342 if( ran < w0/wt )
343 { npos = 0; nneg = 0; nero = 1; }
344 else if ( ran < wp/wt )
345 { npos = 1; nneg = 0; nero = 0; }
346 else
347 { npos = 0; nneg = 1; nero = 0; }
348 }
349 else
350 {
351 w0 = -
sqr(1.+neutb)/(2.*c*c);
352 wm = -
sqr(-1.+neutb)/(2.*c*c);
353 w0 = std::exp(w0);
354 wm = std::exp(wm);
356 { npos = 0; nneg = 0; nero = 1; }
357 else
358 { npos = 0; nneg = 1; nero = 0; }
359 }
360 }
361 else
362 {
363
364
365 G4double aleab = std::log(availableEnergy);
366 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
367 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
368
369
370
372
373 for (nt=1; nt<=numSec; nt++) {
374 test = std::exp(
Amin( expxu,
Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
375 dum =
pi*nt/(2.0*
n*
n);
376 if (std::fabs(dum) < 1.0) {
377 if( test >= 1.0e-10 )anpn += dum*test;
378 } else {
379 anpn += dum*test;
380 }
381 }
382
385 if( targetCode == protonCode )
386 {
387 counter = -1;
388 for (npos=0; npos<numSec/3; npos++) {
389 for (nneg=
Imax(0,npos-1); nneg<=(npos+1); nneg++) {
390 for (nero=0; nero<numSec/3; nero++) {
391 if (++counter < numMul) {
392 nt = npos+nneg+nero;
393 if ( (nt>0) && (nt<=numSec) ) {
394 test = std::exp(
Amin( expxu,
Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
395 dum = (
pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
396 if (std::fabs(dum) < 1.0) {
397 if( test >= 1.0e-10 )excs += dum*test;
398 } else {
399 excs += dum*test;
400 }
401 if (ran < excs) goto outOfLoop;
402 }
403 }
404 }
405 }
406 }
407
408
409 inElastic = false;
410 return;
411 }
412 else
413 {
414 counter = -1;
415 for (npos=0; npos<numSec/3; npos++) {
416 for (nneg=npos; nneg<=(npos+2); nneg++) {
417 for (nero=0; nero<numSec/3; nero++) {
418 if (++counter < numMul) {
419 nt = npos+nneg+nero;
420 if ( (nt>=1) && (nt<=numSec) ) {
421 test = std::exp(
Amin( expxu,
Amax( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
422 dum = (
pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
423 if (std::fabs(dum) < 1.0) {
424 if( test >= 1.0e-10 )excs += dum*test;
425 } else {
426 excs += dum*test;
427 }
428 if (ran < excs) goto outOfLoop;
429 }
430 }
431 }
432 }
433 }
434
435 inElastic = false;
436 return;
437 }
438 }
439 outOfLoop:
440
441 if( targetCode == protonCode)
442 {
443 if( npos == (1+nneg))
444 {
446 }
447 else if (npos == nneg)
448 {
450 {
451 }
452 else
453 {
456 }
457 }
458 else
459 {
461 }
462 }
463 else
464 {
465 if( npos == (nneg-1))
466 {
468 {
470 }
471 else
472 {
474 }
475 }
476 else if ( npos == nneg)
477 {
478 }
479 else
480 {
483 }
484 }
485
486
487 nt = npos + nneg + nero;
488 while ( nt > 0)
489 {
492 {
493 if( npos > 0 )
495 npos--;
496 }
497 }
498 else if ( ran < (
G4double)(npos+nneg)/nt)
499 {
500 if( nneg > 0 )
501 {
503 nneg--;
504 }
505 }
506 else
507 {
508 if( nero > 0 )
509 {
511 nero--;
512 }
513 }
514 nt = npos + nneg + nero;
515 }
517 {
518 G4cout <<
"Particles produced: " ;
521 for (i=2; i < vecLen; i++)
522 {
524 }
526 }
527 return;
528 }
G4double pmltpc(G4int np, G4int nm, G4int nz, G4int n, G4double b, G4double c)
G4double Amin(G4double a, G4double b)
G4int Imin(G4int a, G4int b)
G4double Amax(G4double a, G4double b)
G4int Imax(G4int a, G4int b)
G4double getTotalMomentum() const