Geant4  10.01.p03
G4PreCompoundTransitions.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: G4PreCompoundTransitions.cc 90591 2015-06-04 13:45:29Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4PreCompoundTransitions
34 //
35 // Author: V.Lara
36 //
37 // Modified:
38 // 16.02.2008 J.M.Quesada fixed bugs
39 // 06.09.2008 J.M.Quesada added external choices for:
40 // - "never go back" hipothesis (useNGB=true)
41 // - CEM transition probabilities (useCEMtr=true)
42 // 30.10.2009 J.M.Quesada: CEM transition probabilities have been renormalized
43 // (IAEA benchmark)
44 // 20.08.2010 V.Ivanchenko move constructor and destructor to the source and
45 // optimise the code
46 // 30.08.2011 M.Kelsey - Skip CalculateProbability if no excitons
47 
49 #include "G4PhysicalConstants.hh"
50 #include "G4SystemOfUnits.hh"
51 #include "Randomize.hh"
53 #include "G4Fragment.hh"
54 #include "G4Proton.hh"
55 #include "G4Exp.hh"
56 #include "G4Log.hh"
57 
58 static const G4double sixdpi2 = 6.0/CLHEP::pi2;
59 
61 {
64  FermiEnergy = param.GetFermiEnergy();
65  r0 = param.GetTransitionsr0();
66  aLDP = param.GetLevelDensity();
67 }
68 
70 {}
71 
72 // Calculates transition probabilities with
73 // DeltaN = +2 (Trans1) -2 (Trans2) and 0 (Trans3)
75 CalculateProbability(const G4Fragment & aFragment)
76 {
77  // Number of holes
78  G4int H = aFragment.GetNumberOfHoles();
79  // Number of Particles
80  G4int P = aFragment.GetNumberOfParticles();
81  // Number of Excitons
82  G4int N = P+H;
83  // Nucleus
84  G4int A = aFragment.GetA_asInt();
85  G4int Z = aFragment.GetZ_asInt();
86  G4double U = aFragment.GetExcitationEnergy();
87  TransitionProb2 = 0.0;
88  TransitionProb3 = 0.0;
89 
90  //G4cout << aFragment << G4endl;
91 
92  if(U < 10*eV || 0==N) { return 0.0; }
93 
94  //J. M. Quesada (Feb. 08) new physics
95  // OPT=1 Transitions are calculated according to Gudima's paper
96  // (original in G4PreCompound from VL)
97  // OPT=2 Transitions are calculated according to Gupta's formulae
98  //
99  if (useCEMtr) {
100  // Relative Energy (T_{rel})
101  G4double RelativeEnergy = 1.6*FermiEnergy + U/G4double(N);
102 
103  // Sample kind of nucleon-projectile
104  G4bool ChargedNucleon(false);
105  if(G4int(P*G4UniformRand()) <= aFragment.GetNumberOfCharged()) {
106  ChargedNucleon = true;
107  }
108 
109  // Relative Velocity:
110  // <V_{rel}>^2
111  G4double RelativeVelocitySqr;
112  if (ChargedNucleon) {
113  RelativeVelocitySqr = 2*RelativeEnergy/CLHEP::proton_mass_c2;
114  } else {
115  RelativeVelocitySqr = 2*RelativeEnergy/CLHEP::neutron_mass_c2;
116  }
117  // <V_{rel}>
118  G4double RelativeVelocity = std::sqrt(RelativeVelocitySqr);
119 
120  // Proton-Proton Cross Section
121  G4double ppXSection =
122  (10.63/RelativeVelocitySqr - 29.92/RelativeVelocity + 42.9)
124  // Proton-Neutron Cross Section
125  G4double npXSection =
126  (34.10/RelativeVelocitySqr - 82.20/RelativeVelocity + 82.2)
128 
129  // Averaged Cross Section: \sigma(V_{rel})
130  G4double AveragedXSection;
131  if (ChargedNucleon)
132  {
133  //JMQ: small bug fixed
134  AveragedXSection = ((Z-1)*ppXSection + (A-Z)*npXSection)/G4double(A-1);
135  }
136  else
137  {
138  AveragedXSection = ((A-Z-1)*ppXSection + Z*npXSection)/G4double(A-1);
139  }
140 
141  // Fermi relative energy ratio
142  G4double FermiRelRatio = FermiEnergy/RelativeEnergy;
143 
144  // This factor is introduced to take into account the Pauli principle
145  G4double PauliFactor = 1.0 - 1.4*FermiRelRatio;
146  if (FermiRelRatio > 0.5) {
147  G4double x = 2.0 - 1.0/FermiRelRatio;
148  PauliFactor += 0.4*FermiRelRatio*x*x*std::sqrt(x);
149  }
150  // Interaction volume
151  G4double xx = 2*r0 + CLHEP::hbarc/(CLHEP::proton_mass_c2*RelativeVelocity);
152  G4double Vint = CLHEP::pi*xx*xx*xx/0.75;
153 
154  // Transition probability for \Delta n = +2
155 
156  TransitionProb1 = std::max(0.0, AveragedXSection*PauliFactor
157  *std::sqrt(2.0*RelativeEnergy/CLHEP::proton_mass_c2)/Vint);
158 
159  //JMQ 281009 phenomenological factor in order to increase
160  // equilibrium contribution
161  // G4double factor=5.0;
162  // TransitionProb1 *= factor;
163 
164  // GE = g*E where E is Excitation Energy
165  G4double GE = sixdpi2*aLDP*A*U;
166  G4double Fph = G4double(P*P+H*H+P-3*H)*0.25;
167 
168  if(!useNGB) {
169 
170  // F(p+1,h+1)
171  G4double Fph1 = Fph + N*0.5;
172 
173  static const G4double plimit = 100;
174 
175  //JMQ/AH bug fixed: if (U-Fph < 0.0)
176  if (GE-Fph1 > 0.0) {
177  G4double x0 = GE-Fph;
178  G4double x1 = (N+1)*G4Log(x0/(GE-Fph1));
179  if(x1 < plimit) {
180  x1 = G4Exp(x1)*TransitionProb1/x0;
181 
182  // Transition probability for \Delta n = -2 (at F(p,h) = 0)
183  TransitionProb2 = std::max(0.0, (P*H*(N+1)*(N-2))*x1/x0);
184 
185  // Transition probability for \Delta n = 0 (at F(p,h) = 0)
186  TransitionProb3 = std::max(0.0,((N+1)*(P*(P-1) + 4*P*H + H*(H-1)))*x1
187  /G4double(N));
188  }
189  }
190  }
191 
192  } else {
193  //JMQ: Transition probabilities from Gupta's work
194  // GE = g*E where E is Excitation Energy
195  TransitionProb1 = std::max(0.0, U*(4.2e+12 - 3.6e+10*U/G4double(N+1)))
196  /(16*CLHEP::c_light);
197 
198  if (!useNGB && N > 1) {
199  G4double GE = sixdpi2*aLDP*A*U;
200  TransitionProb2 = ((N-1)*(N-2)*P*H)*TransitionProb1/(GE*GE);
201  }
202  }
203  // G4cout<<"U = "<<U<<G4endl;
204  // G4cout<<"N="<<N<<" P="<<P<<" H="<<H<<G4endl;
205  // G4cout<<"l+ ="<<TransitionProb1<<" l- ="<< TransitionProb2
206  // <<" l0 ="<< TransitionProb3<<G4endl;
208 }
209 
211 {
212  G4double ChosenTransition =
214  G4int deltaN = 0;
215  G4int Npart = result.GetNumberOfParticles();
216  G4int Ncharged = result.GetNumberOfCharged();
217  G4int Nholes = result.GetNumberOfHoles();
218  if (ChosenTransition <= TransitionProb1)
219  {
220  // Number of excitons is increased on \Delta n = +2
221  deltaN = 2;
222  }
223  else if (ChosenTransition <= TransitionProb1+TransitionProb2)
224  {
225  // Number of excitons is increased on \Delta n = -2
226  deltaN = -2;
227  }
228 
229  // AH/JMQ: Randomly decrease the number of charges if deltaN is -2 and
230  // in proportion to the number charges w.r.t. number of particles,
231  // PROVIDED that there are charged particles
232  deltaN /= 2;
233 
234  //G4cout << "deltaN= " << deltaN << G4endl;
235 
236  // JMQ the following lines have to be before SetNumberOfCharged,
237  // otherwise the check on number of charged vs. number of particles fails
238  result.SetNumberOfParticles(Npart+deltaN);
239  result.SetNumberOfHoles(Nholes+deltaN);
240 
241  if(deltaN < 0) {
242  if( (Ncharged == Npart) ||
243  (Ncharged >= 1 && G4int(Npart*G4UniformRand()) <= Ncharged))
244  {
245  result.SetNumberOfCharged(Ncharged+deltaN); // deltaN is negative!
246  }
247 
248  } else if ( deltaN > 0 ) {
249  // With weight Z/A, number of charged particles is increased with +1
250  G4int A = result.GetA_asInt() - Npart;
251  G4int Z = result.GetZ_asInt() - Ncharged;
252  if((Z == A) || (Z > 0 && G4int(A*G4UniformRand()) <= Z))
253  {
254  result.SetNumberOfCharged(Ncharged+deltaN);
255  }
256  }
257 
258  // Number of charged can not be greater that number of particles
259  if ( Npart < Ncharged )
260  {
261  result.SetNumberOfCharged(Npart);
262  }
263  //G4cout << "### After transition" << G4endl;
264  //G4cout << result << G4endl;
265 }
266 
const G4double pi
void SetNumberOfHoles(G4int valueTot, G4int valueP=0)
Definition: G4Fragment.hh:360
G4int GetNumberOfParticles() const
Definition: G4Fragment.hh:330
Definition: Evaluator.cc:66
static const G4double sixdpi2
int G4int
Definition: G4Types.hh:78
G4int GetNumberOfHoles() const
Definition: G4Fragment.hh:350
#define G4UniformRand()
Definition: Randomize.hh:93
virtual G4double CalculateProbability(const G4Fragment &aFragment)
G4int GetA_asInt() const
Definition: G4Fragment.hh:243
bool G4bool
Definition: G4Types.hh:79
static G4Proton * Proton()
Definition: G4Proton.cc:93
const G4ParticleDefinition * proton
virtual void PerformTransition(G4Fragment &aFragment)
void SetNumberOfParticles(G4int value)
Definition: G4Fragment.hh:369
static const G4double A[nN]
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static const double eV
Definition: G4SIunits.hh:194
T max(const T t1, const T t2)
brief Return the largest of the two arguments
static const double millibarn
Definition: G4SIunits.hh:96
G4int GetZ_asInt() const
Definition: G4Fragment.hh:248
void SetNumberOfCharged(G4int value)
Definition: G4Fragment.hh:374
double G4double
Definition: G4Types.hh:76
G4int GetNumberOfCharged() const
Definition: G4Fragment.hh:335
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:260