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

Last change on this file since 1201 was 1055, checked in by garnier, 17 years ago

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