Geant4  10.00.p01
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 68028 2013-03-13 13:48:15Z 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"
52 #include "G4Pow.hh"
53 #include "G4HadronicException.hh"
55 #include "G4Proton.hh"
56 
58 {
61  FermiEnergy = param.GetFermiEnergy();
62  r0 = param.GetTransitionsr0();
63  aLDP = param.GetLevelDensity();
65 }
66 
68 {}
69 
70 // Calculates transition probabilities with
71 // DeltaN = +2 (Trans1) -2 (Trans2) and 0 (Trans3)
73 CalculateProbability(const G4Fragment & aFragment)
74 {
75  // Number of holes
76  G4int H = aFragment.GetNumberOfHoles();
77  // Number of Particles
78  G4int P = aFragment.GetNumberOfParticles();
79  // Number of Excitons
80  G4int N = P+H;
81  // Nucleus
82  G4int A = aFragment.GetA_asInt();
83  G4int Z = aFragment.GetZ_asInt();
84  G4double U = aFragment.GetExcitationEnergy();
85 
86  //G4cout << aFragment << G4endl;
87 
88  if(U < 10*eV || 0==N) { return 0.0; }
89 
90  //J. M. Quesada (Feb. 08) new physics
91  // OPT=1 Transitions are calculated according to Gudima's paper
92  // (original in G4PreCompound from VL)
93  // OPT=2 Transitions are calculated according to Gupta's formulae
94  //
95  if (useCEMtr){
96 
97  // Relative Energy (T_{rel})
98  G4double RelativeEnergy = 1.6*FermiEnergy + U/G4double(N);
99 
100  // Sample kind of nucleon-projectile
101  G4bool ChargedNucleon(false);
102  G4double chtest = 0.5;
103  if (P > 0) {
104  chtest = G4double(aFragment.GetNumberOfCharged())/G4double(P);
105  }
106  if (G4UniformRand() < chtest) { ChargedNucleon = true; }
107 
108  // Relative Velocity:
109  // <V_{rel}>^2
110  G4double RelativeVelocitySqr(0.0);
111  if (ChargedNucleon) {
112  RelativeVelocitySqr = 2.0*RelativeEnergy/CLHEP::proton_mass_c2;
113  } else {
114  RelativeVelocitySqr = 2.0*RelativeEnergy/CLHEP::neutron_mass_c2;
115  }
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 = (ppXSection+npXSection)/2.0;
131  G4double AveragedXSection(0.0);
132  if (ChargedNucleon)
133  {
134  //JMQ: small bug fixed
135  //AveragedXSection=((Z-1.0) * ppXSection + (A-Z-1.0) * npXSection)/(A-1.0);
136  AveragedXSection = ((Z-1)*ppXSection + (A-Z)*npXSection)/G4double(A-1);
137  }
138  else
139  {
140  AveragedXSection = ((A-Z-1)*ppXSection + Z*npXSection)/G4double(A-1);
141  //AveragedXSection = ((A-Z-1)*npXSection + Z*ppXSection)/G4double(A-1);
142  }
143 
144  // Fermi relative energy ratio
145  G4double FermiRelRatio = FermiEnergy/RelativeEnergy;
146 
147  // This factor is introduced to take into account the Pauli principle
148  G4double PauliFactor = 1.0 - 1.4*FermiRelRatio;
149  if (FermiRelRatio > 0.5) {
150  G4double x = 2.0 - 1.0/FermiRelRatio;
151  PauliFactor += 0.4*FermiRelRatio*x*x*std::sqrt(x);
152  //PauliFactor +=
153  //(2.0/5.0)*FermiRelRatio*std::pow(2.0 - (1.0/FermiRelRatio), 5.0/2.0);
154  }
155  // Interaction volume
156  // G4double Vint = (4.0/3.0)
157  //*pi*std::pow(2.0*r0 + hbarc/(proton_mass_c2*RelativeVelocity) , 3.0);
158  G4double xx = 2.0*r0 + hbarc/(CLHEP::proton_mass_c2*RelativeVelocity);
159  // G4double Vint = (4.0/3.0)*CLHEP::pi*xx*xx*xx;
160  G4double Vint = CLHEP::pi*xx*xx*xx/0.75;
161 
162  // Transition probability for \Delta n = +2
163 
164  TransitionProb1 = AveragedXSection*PauliFactor
165  *std::sqrt(2.0*RelativeEnergy/CLHEP::proton_mass_c2)/Vint;
166 
167  //JMQ 281009 phenomenological factor in order to increase
168  // equilibrium contribution
169  // G4double factor=5.0;
170  // TransitionProb1 *= factor;
171  //
172  if (TransitionProb1 < 0.0) { TransitionProb1 = 0.0; }
173 
174  // GE = g*E where E is Excitation Energy
175  G4double GE = (6.0/pi2)*aLDP*A*U;
176 
177  //G4double Fph = ((P*P+H*H+P-H)/4.0 - H/2.0);
178  G4double Fph = G4double(P*P+H*H+P-3*H)/4.0;
179 
180  G4bool NeverGoBack(false);
181  if(useNGB) { NeverGoBack=true; }
182 
183  //JMQ/AH bug fixed: if (U-Fph < 0.0) NeverGoBack = true;
184  if (GE-Fph < 0.0) { NeverGoBack = true; }
185 
186  // F(p+1,h+1)
187  G4double Fph1 = Fph + N/2.0;
188 
189  G4double ProbFactor = g4pow->powN((GE-Fph)/(GE-Fph1),N+1);
190 
191  if (NeverGoBack)
192  {
193  TransitionProb2 = 0.0;
194  TransitionProb3 = 0.0;
195  }
196  else
197  {
198  // Transition probability for \Delta n = -2 (at F(p,h) = 0)
199  TransitionProb2 =
200  TransitionProb1 * ProbFactor * (P*H*(N+1)*(N-2))/((GE-Fph)*(GE-Fph));
201  if (TransitionProb2 < 0.0) { TransitionProb2 = 0.0; }
202 
203  // Transition probability for \Delta n = 0 (at F(p,h) = 0)
204  TransitionProb3 = TransitionProb1*(N+1)* ProbFactor
205  * (P*(P-1) + 4.0*P*H + H*(H-1))/(N*(GE-Fph));
206  if (TransitionProb3 < 0.0) { TransitionProb3 = 0.0; }
207  }
208  } else {
209  //JMQ: Transition probabilities from Gupta's work
210  // GE = g*E where E is Excitation Energy
211  G4double GE = (6.0/pi2)*aLDP*A*U;
212 
213  G4double Kmfp=2.;
214 
215  //TransitionProb1=1./Kmfp*3./8.*1./c_light*1.0e-9*(1.4e+21*U-2./(N+1)*6.0e+18*U*U);
216  TransitionProb1 = 3.0e-9*(1.4e+21*U - 1.2e+19*U*U/G4double(N+1))
217  /(8*Kmfp*CLHEP::c_light);
218  if (TransitionProb1 < 0.0) { TransitionProb1 = 0.0; }
219 
220  TransitionProb2=0.;
221  TransitionProb3=0.;
222 
223  if (!useNGB && N > 1) {
224  // TransitionProb2=1./Kmfp*3./8.*1./c_light*1.0e-9*(N-1.)*(N-2.)*P*H/(GE*GE)*(1.4e+21*U - 2./(N-1)*6.0e+18*U*U);
225  TransitionProb2 =
226  3.0e-9*(N-2)*P*H*(1.4e+21*U*(N-1) - 1.2e+19*U*U)/(8*Kmfp*c_light*GE*GE);
227  if (TransitionProb2 < 0.0) TransitionProb2 = 0.0;
228  }
229  }
230  // G4cout<<"U = "<<U<<G4endl;
231  // G4cout<<"N="<<N<<" P="<<P<<" H="<<H<<G4endl;
232  // G4cout<<"l+ ="<<TransitionProb1<<" l- ="<< TransitionProb2
233  // <<" l0 ="<< TransitionProb3<<G4endl;
235 }
236 
238 {
239  G4double ChosenTransition =
241  G4int deltaN = 0;
242  // G4int Nexcitons = result.GetNumberOfExcitons();
243  G4int Npart = result.GetNumberOfParticles();
244  G4int Ncharged = result.GetNumberOfCharged();
245  G4int Nholes = result.GetNumberOfHoles();
246  if (ChosenTransition <= TransitionProb1)
247  {
248  // Number of excitons is increased on \Delta n = +2
249  deltaN = 2;
250  }
251  else if (ChosenTransition <= TransitionProb1+TransitionProb2)
252  {
253  // Number of excitons is increased on \Delta n = -2
254  deltaN = -2;
255  }
256 
257  // AH/JMQ: Randomly decrease the number of charges if deltaN is -2 and
258  // in proportion to the number charges w.r.t. number of particles,
259  // PROVIDED that there are charged particles
260  deltaN /= 2;
261 
262  //G4cout << "deltaN= " << deltaN << G4endl;
263 
264  // JMQ the following lines have to be before SetNumberOfCharged, otherwise the check on
265  // number of charged vs. number of particles fails
266  result.SetNumberOfParticles(Npart+deltaN);
267  result.SetNumberOfHoles(Nholes+deltaN);
268 
269  if(deltaN < 0) {
270  if( Ncharged >= 1 && G4int(Npart*G4UniformRand()) <= Ncharged)
271  {
272  result.SetNumberOfCharged(Ncharged+deltaN); // deltaN is negative!
273  }
274 
275  } else if ( deltaN > 0 ) {
276  // With weight Z/A, number of charged particles is increased with +1
277  G4int A = result.GetA_asInt();
278  G4int Z = result.GetZ_asInt();
279  if( G4int(std::max(1, A - Npart)*G4UniformRand()) <= Z)
280  {
281  result.SetNumberOfCharged(Ncharged+deltaN);
282  }
283  }
284 
285  // Number of charged can not be greater that number of particles
286  if ( Npart < Ncharged )
287  {
288  result.SetNumberOfCharged(Npart);
289  }
290  //G4cout << "### After transition" << G4endl;
291  //G4cout << result << G4endl;
292 }
293 
static G4Pow * GetInstance()
Definition: G4Pow.cc:53
G4double powN(G4double x, G4int n) const
Definition: G4Pow.cc:125
const G4double pi
void SetNumberOfHoles(G4int valueTot, G4int valueP=0)
Definition: G4Fragment.hh:355
G4int GetNumberOfParticles() const
Definition: G4Fragment.hh:325
int G4int
Definition: G4Types.hh:78
G4int GetNumberOfHoles() const
Definition: G4Fragment.hh:345
#define G4UniformRand()
Definition: Randomize.hh:87
virtual G4double CalculateProbability(const G4Fragment &aFragment)
G4int GetA_asInt() const
Definition: G4Fragment.hh:238
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:364
static const G4double A[nN]
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:243
void SetNumberOfCharged(G4int value)
Definition: G4Fragment.hh:369
double G4double
Definition: G4Types.hh:76
G4int GetNumberOfCharged() const
Definition: G4Fragment.hh:330
Definition: Evaluator.cc:66
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:255