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