source: trunk/source/processes/hadronic/stopping/src/G4AntiNeutronAnnihilationAtRest.cc@ 1036

Last change on this file since 1036 was 1007, checked in by garnier, 17 years ago

update to geant4.9.2
File size: 19.9 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// G4AntiNeutronAnnihilationAtRest physics process
27// Larry Felawka (TRIUMF), April 1998
28//---------------------------------------------------------------------
29
30#include "G4AntiNeutronAnnihilationAtRest.hh"
31#include "G4DynamicParticle.hh"
32#include "G4ParticleTypes.hh"
33#include "Randomize.hh"
34#include <string.h>
35#include <cmath>
36#include <stdio.h>
37
38#define MAX_SECONDARIES 100
39
40// constructor
41
42G4AntiNeutronAnnihilationAtRest::G4AntiNeutronAnnihilationAtRest(const G4String& processName,
43 G4ProcessType aType) :
44 G4VRestProcess (processName, aType), // initialization
45 massPionMinus(G4PionMinus::PionMinus()->GetPDGMass()/GeV),
46 massPionZero(G4PionZero::PionZero()->GetPDGMass()/GeV),
47 massPionPlus(G4PionPlus::PionPlus()->GetPDGMass()/GeV),
48 massGamma(G4Gamma::Gamma()->GetPDGMass()/GeV),
49 massAntiNeutron(G4AntiNeutron::AntiNeutron()->GetPDGMass()/GeV),
50 massNeutron(G4Neutron::Neutron()->GetPDGMass()/GeV),
51 pdefGamma(G4Gamma::Gamma()),
52 pdefPionPlus(G4PionPlus::PionPlus()),
53 pdefPionZero(G4PionZero::PionZero()),
54 pdefPionMinus(G4PionMinus::PionMinus()),
55 pdefProton(G4Proton::Proton()),
56 pdefNeutron(G4Neutron::Neutron()),
57 pdefAntiNeutron(G4AntiNeutron::AntiNeutron()),
58 pdefDeuteron(G4Deuteron::Deuteron()),
59 pdefTriton(G4Triton::Triton()),
60 pdefAlpha(G4Alpha::Alpha())
61{
62 if (verboseLevel>0) {
63 G4cout << GetProcessName() << " is created "<< G4endl;
64 }
65 SetProcessSubType(fHadronAtRest);
66 pv = new G4GHEKinematicsVector [MAX_SECONDARIES+1];
67 eve = new G4GHEKinematicsVector [MAX_SECONDARIES];
68 gkin = new G4GHEKinematicsVector [MAX_SECONDARIES];
69
70}
71
72// destructor
73
74G4AntiNeutronAnnihilationAtRest::~G4AntiNeutronAnnihilationAtRest()
75{
76 delete [] pv;
77 delete [] eve;
78 delete [] gkin;
79}
80
81
82// methods.............................................................................
83
84G4bool G4AntiNeutronAnnihilationAtRest::IsApplicable(
85 const G4ParticleDefinition& particle
86 )
87{
88 return ( &particle == pdefAntiNeutron );
89
90}
91
92// Warning - this method may be optimized away if made "inline"
93G4int G4AntiNeutronAnnihilationAtRest::GetNumberOfSecondaries()
94{
95 return ( ngkine );
96
97}
98
99// Warning - this method may be optimized away if made "inline"
100G4GHEKinematicsVector* G4AntiNeutronAnnihilationAtRest::GetSecondaryKinematics()
101{
102 return ( &gkin[0] );
103
104}
105
106G4double G4AntiNeutronAnnihilationAtRest::AtRestGetPhysicalInteractionLength(
107 const G4Track& track,
108 G4ForceCondition* condition
109 )
110{
111 // beggining of tracking
112 ResetNumberOfInteractionLengthLeft();
113
114 // condition is set to "Not Forced"
115 *condition = NotForced;
116
117 // get mean life time
118 currentInteractionLength = GetMeanLifeTime(track, condition);
119
120 if ((currentInteractionLength <0.0) || (verboseLevel>2)){
121 G4cout << "G4AntiNeutronAnnihilationAtRestProcess::AtRestGetPhysicalInteractionLength ";
122 G4cout << "[ " << GetProcessName() << "]" <<G4endl;
123 track.GetDynamicParticle()->DumpInfo();
124 G4cout << " in Material " << track.GetMaterial()->GetName() <<G4endl;
125 G4cout << "MeanLifeTime = " << currentInteractionLength/ns << "[ns]" <<G4endl;
126 }
127
128 return theNumberOfInteractionLengthLeft * currentInteractionLength;
129
130}
131
132G4VParticleChange* G4AntiNeutronAnnihilationAtRest::AtRestDoIt(
133 const G4Track& track,
134 const G4Step&
135 )
136//
137// Handles AntiNeutrons at rest; an AntiNeutron can either create secondaries
138// or do nothing (in which case it should be sent back to decay-handling
139// section
140//
141{
142
143// Initialize ParticleChange
144// all members of G4VParticleChange are set to equal to
145// corresponding member in G4Track
146
147 aParticleChange.Initialize(track);
148
149// Store some global quantities that depend on current material and particle
150
151 globalTime = track.GetGlobalTime()/s;
152 G4Material * aMaterial = track.GetMaterial();
153 const G4int numberOfElements = aMaterial->GetNumberOfElements();
154 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
155
156 const G4double* theAtomicNumberDensity = aMaterial->GetAtomicNumDensityVector();
157 G4double normalization = 0;
158 for ( G4int i1=0; i1 < numberOfElements; i1++ )
159 {
160 normalization += theAtomicNumberDensity[i1] ; // change when nucleon specific
161 // probabilities are included.
162 }
163 G4double runningSum= 0.;
164 G4double random = G4UniformRand()*normalization;
165 for ( G4int i2=0; i2 < numberOfElements; i2++ )
166 {
167 runningSum += theAtomicNumberDensity[i2]; // change when nucleon specific
168 // probabilities are included.
169 if (random<=runningSum)
170 {
171 targetCharge = G4double( ((*theElementVector)[i2])->GetZ());
172 targetAtomicMass = (*theElementVector)[i2]->GetN();
173 }
174 }
175 if (random>runningSum)
176 {
177 targetCharge = G4double( ((*theElementVector)[numberOfElements-1])->GetZ());
178 targetAtomicMass = (*theElementVector)[numberOfElements-1]->GetN();
179 }
180
181 if (verboseLevel>1) {
182 G4cout << "G4AntiNeutronAnnihilationAtRest::AtRestDoIt is invoked " <<G4endl;
183 }
184
185 G4ParticleMomentum momentum;
186 G4float localtime;
187
188 G4ThreeVector position = track.GetPosition();
189
190 GenerateSecondaries(); // Generate secondaries
191
192 aParticleChange.SetNumberOfSecondaries( ngkine );
193
194 for ( G4int isec = 0; isec < ngkine; isec++ ) {
195 G4DynamicParticle* aNewParticle = new G4DynamicParticle;
196 aNewParticle->SetDefinition( gkin[isec].GetParticleDef() );
197 aNewParticle->SetMomentum( gkin[isec].GetMomentum() * GeV );
198
199 localtime = globalTime + gkin[isec].GetTOF();
200
201 G4Track* aNewTrack = new G4Track( aNewParticle, localtime*s, position );
202 aNewTrack->SetTouchableHandle(track.GetTouchableHandle());
203 aParticleChange.AddSecondary( aNewTrack );
204
205 }
206
207 aParticleChange.ProposeLocalEnergyDeposit( 0.0*GeV );
208
209 aParticleChange.ProposeTrackStatus(fStopAndKill); // Kill the incident AntiNeutron
210
211// clear InteractionLengthLeft
212
213 ResetNumberOfInteractionLengthLeft();
214
215 return &aParticleChange;
216
217}
218
219
220void G4AntiNeutronAnnihilationAtRest::GenerateSecondaries()
221{
222 static G4int index;
223 static G4int l;
224 static G4int nopt;
225 static G4int i;
226 static G4ParticleDefinition* jnd;
227
228 for (i = 1; i <= MAX_SECONDARIES; ++i) {
229 pv[i].SetZero();
230 }
231
232
233 ngkine = 0; // number of generated secondary particles
234 ntot = 0;
235 result.SetZero();
236 result.SetMass( massAntiNeutron );
237 result.SetKineticEnergyAndUpdate( 0. );
238 result.SetTOF( 0. );
239 result.SetParticleDef( pdefAntiNeutron );
240
241 // *** SELECT PROCESS FOR CURRENT PARTICLE ***
242
243 AntiNeutronAnnihilation(&nopt);
244
245 // *** CHECK WHETHER THERE ARE NEW PARTICLES GENERATED ***
246 if (ntot != 0 || result.GetParticleDef() != pdefAntiNeutron) {
247 // *** CURRENT PARTICLE IS NOT THE SAME AS IN THE BEGINNING OR/AND ***
248 // *** ONE OR MORE SECONDARIES HAVE BEEN GENERATED ***
249
250 // --- INITIAL PARTICLE TYPE HAS BEEN CHANGED ==> PUT NEW TYPE ON ---
251 // --- THE GEANT TEMPORARY STACK ---
252
253 // --- PUT PARTICLE ON THE STACK ---
254 gkin[0] = result;
255 gkin[0].SetTOF( result.GetTOF() * 5e-11 );
256 ngkine = 1;
257
258 // --- ALL QUANTITIES ARE TAKEN FROM THE GHEISHA STACK WHERE THE ---
259 // --- CONVENTION IS THE FOLLOWING ---
260
261 // --- ONE OR MORE SECONDARIES HAVE BEEN GENERATED ---
262 for (l = 1; l <= ntot; ++l) {
263 index = l - 1;
264 jnd = eve[index].GetParticleDef();
265
266 // --- ADD PARTICLE TO THE STACK IF STACK NOT YET FULL ---
267 if (ngkine < MAX_SECONDARIES) {
268 gkin[ngkine] = eve[index];
269 gkin[ngkine].SetTOF( eve[index].GetTOF() * 5e-11 );
270 ++ngkine;
271 }
272 }
273 }
274 else {
275 // --- NO SECONDARIES GENERATED AND PARTICLE IS STILL THE SAME ---
276 // --- ==> COPY EVERYTHING BACK IN THE CURRENT GEANT STACK ---
277 ngkine = 0;
278 ntot = 0;
279 globalTime += result.GetTOF() * G4float(5e-11);
280 }
281
282 // --- LIMIT THE VALUE OF NGKINE IN CASE OF OVERFLOW ---
283 ngkine = G4int(std::min(ngkine,G4int(MAX_SECONDARIES)));
284
285} // GenerateSecondaries
286
287
288void G4AntiNeutronAnnihilationAtRest::Poisso(G4float xav, G4int *iran)
289{
290 static G4int i;
291 static G4float r, p1, p2, p3;
292 static G4int mm;
293 static G4float rr, ran, rrr, ran1;
294
295 // *** GENERATION OF POISSON DISTRIBUTION ***
296 // *** NVE 16-MAR-1988 CERN GENEVA ***
297 // ORIGIN : H.FESEFELDT (27-OCT-1983)
298
299 // --- USE NORMAL DISTRIBUTION FOR <X> > 9.9 ---
300 if (xav > G4float(9.9)) {
301 // ** NORMAL DISTRIBUTION WITH SIGMA**2 = <X>
302 Normal(&ran1);
303 ran1 = xav + ran1 * std::sqrt(xav);
304 *iran = G4int(ran1);
305 if (*iran < 0) {
306 *iran = 0;
307 }
308 }
309 else {
310 mm = G4int(xav * G4float(5.));
311 *iran = 0;
312 if (mm > 0) {
313 r = std::exp(-G4double(xav));
314 ran1 = G4UniformRand();
315 if (ran1 > r) {
316 rr = r;
317 for (i = 1; i <= mm; ++i) {
318 ++(*iran);
319 if (i <= 5) {
320 rrr = std::pow(xav, G4float(i)) / NFac(i);
321 }
322 // ** STIRLING' S FORMULA FOR LARGE NUMBERS
323 if (i > 5) {
324 rrr = std::exp(i * std::log(xav) -
325 (i + G4float(.5)) * std::log(i * G4float(1.)) +
326 i - G4float(.9189385));
327 }
328 rr += r * rrr;
329 if (ran1 <= rr) {
330 break;
331 }
332 }
333 }
334 }
335 else {
336 // ** FOR VERY SMALL XAV TRY IRAN=1,2,3
337 p1 = xav * std::exp(-G4double(xav));
338 p2 = xav * p1 / G4float(2.);
339 p3 = xav * p2 / G4float(3.);
340 ran = G4UniformRand();
341 if (ran >= p3) {
342 if (ran >= p2) {
343 if (ran >= p1) {
344 *iran = 0;
345 }
346 else {
347 *iran = 1;
348 }
349 }
350 else {
351 *iran = 2;
352 }
353 }
354 else {
355 *iran = 3;
356 }
357 }
358 }
359
360} // Poisso
361
362
363G4int G4AntiNeutronAnnihilationAtRest::NFac(G4int n)
364{
365 G4int ret_val;
366
367 static G4int i, m;
368
369 // *** NVE 16-MAR-1988 CERN GENEVA ***
370 // ORIGIN : H.FESEFELDT (27-OCT-1983)
371
372 ret_val = 1;
373 m = n;
374 if (m > 1) {
375 if (m > 10) {
376 m = 10;
377 }
378 for (i = 2; i <= m; ++i) {
379 ret_val *= i;
380 }
381 }
382 return ret_val;
383
384} // NFac
385
386
387void G4AntiNeutronAnnihilationAtRest::Normal(G4float *ran)
388{
389 static G4int i;
390
391 // *** NVE 14-APR-1988 CERN GENEVA ***
392 // ORIGIN : H.FESEFELDT (27-OCT-1983)
393
394 *ran = G4float(-6.);
395 for (i = 1; i <= 12; ++i) {
396 *ran += G4UniformRand();
397 }
398
399} // Normal
400
401
402void G4AntiNeutronAnnihilationAtRest::AntiNeutronAnnihilation(G4int *nopt)
403{
404 static G4float brr[3] = { G4float(.125),G4float(.25),G4float(.5) };
405
406 G4float r__1;
407
408 static G4int i, ii, kk;
409 static G4int nt;
410 static G4float cfa, eka;
411 static G4int ika, nbl;
412 static G4float ran, pcm;
413 static G4int isw;
414 static G4float tex;
415 static G4ParticleDefinition* ipa1;
416 static G4float ran1, ran2, ekin, tkin;
417 static G4float targ;
418 static G4ParticleDefinition* inve;
419 static G4float ekin1, ekin2, black;
420 static G4float pnrat, rmnve1, rmnve2;
421 static G4float ek, en;
422
423 // *** ANTI NEUTRON ANNIHILATION AT REST ***
424 // *** NVE 04-MAR-1988 CERN GENEVA ***
425 // ORIGIN : H.FESEFELDT (09-JULY-1987)
426
427 // NOPT=0 NO ANNIHILATION
428 // NOPT=1 ANNIH.IN PI+ PI-
429 // NOPT=2 ANNIH.IN PI0 PI0
430 // NOPT=3 ANNIH.IN PI+ PI0
431 // NOPT=4 ANNIH.IN GAMMA GAMMA
432
433 pv[1].SetZero();
434 pv[1].SetMass( massAntiNeutron );
435 pv[1].SetKineticEnergyAndUpdate( 0. );
436 pv[1].SetTOF( result.GetTOF() );
437 pv[1].SetParticleDef( result.GetParticleDef() );
438 isw = 1;
439 ran = G4UniformRand();
440 if (ran > brr[0]) {
441 isw = 2;
442 }
443 if (ran > brr[1]) {
444 isw = 3;
445 }
446 if (ran > brr[2]) {
447 isw = 4;
448 }
449 *nopt = isw;
450 // **
451 // ** EVAPORATION
452 // **
453 rmnve1 = massPionPlus;
454 rmnve2 = massPionMinus;
455 if (isw == 2) {
456 rmnve1 = massPionZero;
457 }
458 if (isw == 2) {
459 rmnve2 = massPionZero;
460 }
461 if (isw == 3) {
462 rmnve2 = massPionZero;
463 }
464 if (isw == 4) {
465 rmnve1 = massGamma;
466 }
467 if (isw == 4) {
468 rmnve2 = massGamma;
469 }
470 ek = massNeutron + massAntiNeutron - rmnve1 - rmnve2;
471 tkin = ExNu(ek);
472 ek -= tkin;
473 if (ek < G4float(1e-4)) {
474 ek = G4float(1e-4);
475 }
476 ek /= G4float(2.);
477 en = ek + (rmnve1 + rmnve2) / G4float(2.);
478 r__1 = en * en - rmnve1 * rmnve2;
479 pcm = r__1 > 0 ? std::sqrt(r__1) : 0;
480 pv[2].SetZero();
481 pv[2].SetMass( rmnve1 );
482 pv[3].SetZero();
483 pv[3].SetMass( rmnve2 );
484 if (isw > 3) {
485 pv[2].SetMass( 0. );
486 pv[3].SetMass( 0. );
487 }
488 pv[2].SetEnergyAndUpdate( std::sqrt(pv[2].GetMass()*pv[2].GetMass()+pcm*pcm) );
489 pv[2].SetTOF( result.GetTOF() );
490 pv[3].SetEnergy( std::sqrt(pv[3].GetMass()*pv[3].GetMass()+pcm*pcm) );
491 pv[3].SetMomentumAndUpdate( -pv[2].GetMomentum().x(), -pv[2].GetMomentum().y(), -pv[2].GetMomentum().z() );
492 pv[3].SetTOF( result.GetTOF() );
493 switch ((int)isw) {
494 case 1:
495 pv[2].SetParticleDef( pdefPionPlus );
496 pv[3].SetParticleDef( pdefPionMinus );
497 break;
498 case 2:
499 pv[2].SetParticleDef( pdefPionZero );
500 pv[3].SetParticleDef( pdefPionZero );
501 break;
502 case 3:
503 pv[2].SetParticleDef( pdefPionPlus );
504 pv[3].SetParticleDef( pdefPionZero );
505 break;
506 case 4:
507 pv[2].SetParticleDef( pdefGamma );
508 pv[3].SetParticleDef( pdefGamma );
509 break;
510 default:
511 break;
512 }
513 nt = 3;
514 if (targetAtomicMass >= G4float(1.5)) {
515 cfa = (targetAtomicMass - G4float(1.)) / G4float(120.) *
516 G4float(.025) * std::exp(-G4double(targetAtomicMass - G4float(1.)) /
517 G4float(120.));
518 targ = G4float(1.);
519 tex = evapEnergy1;
520 if (tex >= G4float(.001)) {
521 black = (targ * G4float(1.25) +
522 G4float(1.5)) * evapEnergy1 / (evapEnergy1 + evapEnergy3);
523 Poisso(black, &nbl);
524 if (G4float(G4int(targ) + nbl) > targetAtomicMass) {
525 nbl = G4int(targetAtomicMass - targ);
526 }
527 if (nt + nbl > (MAX_SECONDARIES - 2)) {
528 nbl = (MAX_SECONDARIES - 2) - nt;
529 }
530 if (nbl > 0) {
531 ekin = tex / nbl;
532 ekin2 = G4float(0.);
533 for (i = 1; i <= nbl; ++i) {
534 if (nt == (MAX_SECONDARIES - 2)) {
535 continue;
536 }
537 if (ekin2 > tex) {
538 break;
539 }
540 ran1 = G4UniformRand();
541 Normal(&ran2);
542 ekin1 = -G4double(ekin) * std::log(ran1) -
543 cfa * (ran2 * G4float(.5) + G4float(1.));
544 if (ekin1 < G4float(0.)) {
545 ekin1 = std::log(ran1) * G4float(-.01);
546 }
547 ekin1 *= G4float(1.);
548 ekin2 += ekin1;
549 if (ekin2 > tex) {
550 ekin1 = tex - (ekin2 - ekin1);
551 }
552 if (ekin1 < G4float(0.)) {
553 ekin1 = G4float(.001);
554 }
555 ipa1 = pdefNeutron;
556 pnrat = G4float(1.) - targetCharge / targetAtomicMass;
557 if (G4UniformRand() > pnrat) {
558 ipa1 = pdefProton;
559 }
560 ++nt;
561 pv[nt].SetZero();
562 pv[nt].SetMass( ipa1->GetPDGMass()/GeV );
563 pv[nt].SetKineticEnergyAndUpdate( ekin1 );
564 pv[nt].SetTOF( result.GetTOF() );
565 pv[nt].SetParticleDef( ipa1 );
566 }
567 if (targetAtomicMass >= G4float(230.) && ek <= G4float(2.)) {
568 ii = nt + 1;
569 kk = 0;
570 eka = ek;
571 if (eka > G4float(1.)) {
572 eka *= eka;
573 }
574 if (eka < G4float(.1)) {
575 eka = G4float(.1);
576 }
577 ika = G4int(G4float(3.6) / eka);
578 for (i = 1; i <= nt; ++i) {
579 --ii;
580 if (pv[ii].GetParticleDef() != pdefProton) {
581 continue;
582 }
583 ipa1 = pdefNeutron;
584 pv[ii].SetMass( ipa1->GetPDGMass()/GeV );
585 pv[ii].SetParticleDef( ipa1 );
586 ++kk;
587 if (kk > ika) {
588 break;
589 }
590 }
591 }
592 }
593 }
594 // **
595 // ** THEN ALSO DEUTERONS, TRITONS AND ALPHAS
596 // **
597 tex = evapEnergy3;
598 if (tex >= G4float(.001)) {
599 black = (targ * G4float(1.25) + G4float(1.5)) * evapEnergy3 /
600 (evapEnergy1 + evapEnergy3);
601 Poisso(black, &nbl);
602 if (nt + nbl > (MAX_SECONDARIES - 2)) {
603 nbl = (MAX_SECONDARIES - 2) - nt;
604 }
605 if (nbl > 0) {
606 ekin = tex / nbl;
607 ekin2 = G4float(0.);
608 for (i = 1; i <= nbl; ++i) {
609 if (nt == (MAX_SECONDARIES - 2)) {
610 continue;
611 }
612 if (ekin2 > tex) {
613 break;
614 }
615 ran1 = G4UniformRand();
616 Normal(&ran2);
617 ekin1 = -G4double(ekin) * std::log(ran1) -
618 cfa * (ran2 * G4float(.5) + G4float(1.));
619 if (ekin1 < G4float(0.)) {
620 ekin1 = std::log(ran1) * G4float(-.01);
621 }
622 ekin1 *= G4float(1.);
623 ekin2 += ekin1;
624 if (ekin2 > tex) {
625 ekin1 = tex - (ekin2 - ekin1);
626 }
627 if (ekin1 < G4float(0.)) {
628 ekin1 = G4float(.001);
629 }
630 ran = G4UniformRand();
631 inve = pdefDeuteron;
632 if (ran > G4float(.6)) {
633 inve = pdefTriton;
634 }
635 if (ran > G4float(.9)) {
636 inve = pdefAlpha;
637 }
638 ++nt;
639 pv[nt].SetZero();
640 pv[nt].SetMass( inve->GetPDGMass()/GeV );
641 pv[nt].SetKineticEnergyAndUpdate( ekin1 );
642 pv[nt].SetTOF( result.GetTOF() );
643 pv[nt].SetParticleDef( inve );
644 }
645 }
646 }
647 }
648 result = pv[2];
649 if (nt == 2) {
650 return;
651 }
652 for (i = 3; i <= nt; ++i) {
653 if (ntot >= MAX_SECONDARIES) {
654 return;
655 }
656 eve[ntot++] = pv[i];
657 }
658
659} // AntiNeutronAnnihilation
660
661
662G4double G4AntiNeutronAnnihilationAtRest::ExNu(G4float ek1)
663{
664 G4float ret_val, r__1;
665
666 static G4float cfa, gfa, ran1, ran2, ekin1, atno3;
667 static G4int magic;
668 static G4float fpdiv;
669
670 // *** NUCLEAR EVAPORATION AS FUNCTION OF ATOMIC NUMBER ATNO ***
671 // *** AND KINETIC ENERGY EKIN OF PRIMARY PARTICLE ***
672 // *** NVE 04-MAR-1988 CERN GENEVA ***
673 // ORIGIN : H.FESEFELDT (10-DEC-1986)
674
675 ret_val = G4float(0.);
676 if (targetAtomicMass >= G4float(1.5)) {
677 magic = 0;
678 if (G4int(targetCharge + G4float(.1)) == 82) {
679 magic = 1;
680 }
681 ekin1 = ek1;
682 if (ekin1 < G4float(.1)) {
683 ekin1 = G4float(.1);
684 }
685 if (ekin1 > G4float(4.)) {
686 ekin1 = G4float(4.);
687 }
688 // ** 0.35 VALUE AT 1 GEV
689 // ** 0.05 VALUE AT 0.1 GEV
690 cfa = G4float(.13043478260869565);
691 cfa = cfa * std::log(ekin1) + G4float(.35);
692 if (cfa < G4float(.15)) {
693 cfa = G4float(.15);
694 }
695 ret_val = cfa * G4float(7.716) * std::exp(-G4double(cfa));
696 atno3 = targetAtomicMass;
697 if (atno3 > G4float(120.)) {
698 atno3 = G4float(120.);
699 }
700 cfa = (atno3 - G4float(1.)) /
701 G4float(120.) * std::exp(-G4double(atno3 - G4float(1.)) / G4float(120.));
702 ret_val *= cfa;
703 r__1 = ekin1;
704 fpdiv = G4float(1.) - r__1 * r__1 * G4float(.25);
705 if (fpdiv < G4float(.5)) {
706 fpdiv = G4float(.5);
707 }
708 gfa = (targetAtomicMass - G4float(1.)) /
709 G4float(70.) * G4float(2.) *
710 std::exp(-G4double(targetAtomicMass - G4float(1.)) / G4float(70.));
711 evapEnergy1 = ret_val * fpdiv;
712 evapEnergy3 = ret_val - evapEnergy1;
713 Normal(&ran1);
714 Normal(&ran2);
715 if (magic == 1) {
716 ran1 = G4float(0.);
717 ran2 = G4float(0.);
718 }
719 evapEnergy1 *= ran1 * gfa + G4float(1.);
720 if (evapEnergy1 < G4float(0.)) {
721 evapEnergy1 = G4float(0.);
722 }
723 evapEnergy3 *= ran2 * gfa + G4float(1.);
724 if (evapEnergy3 < G4float(0.)) {
725 evapEnergy3 = G4float(0.);
726 }
727 while ((ret_val = evapEnergy1 + evapEnergy3) >= ek1) {
728 evapEnergy1 *= G4float(1.) - G4UniformRand() * G4float(.5);
729 evapEnergy3 *= G4float(1.) - G4UniformRand() * G4float(.5);
730 }
731 }
732 return ret_val;
733
734} // ExNu
Note: See TracBrowser for help on using the repository browser.