source: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundTransitions.cc @ 1228

Last change on this file since 1228 was 1228, checked in by garnier, 14 years ago

update geant4.9.3 tag
File size: 10.5 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// $Id: G4PreCompoundTransitions.cc,v 1.22 2009/11/21 18:03:13 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-03 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name:     G4PreCompoundIon
35//
36// Author:         V.Lara
37//
38// Modified: 
39// 16.02.2008 J. M. Quesada fixed bugs
40// 06.09.2008 J. M. Quesada added external choices for:
41//                      - "never go back"  hipothesis (useNGB=true)
42//                      -  CEM transition probabilities (useCEMtr=true)
43
44// 30.10.09 J.M.Quesada: CEM transition probabilities have been renormalized
45//                       (IAEA benchmark)
46//
47#include "G4PreCompoundTransitions.hh"
48#include "G4HadronicException.hh"
49
50const G4PreCompoundTransitions & G4PreCompoundTransitions::
51operator=(const G4PreCompoundTransitions &)
52{
53  throw G4HadronicException(__FILE__, __LINE__, "G4PreCompoundTransitions::operator= meant to not be accessable");
54  return *this;
55}
56
57
58G4bool G4PreCompoundTransitions::operator==(const G4PreCompoundTransitions &) const
59{
60  return false;
61}
62
63G4bool G4PreCompoundTransitions::operator!=(const G4PreCompoundTransitions &) const
64{
65  return true;
66}
67
68
69G4double G4PreCompoundTransitions::
70CalculateProbability(const G4Fragment & aFragment)
71{
72  //G4cout<<"In G4PreCompoundTransitions.cc  useNGB="<<useNGB<<G4endl;
73  //G4cout<<"In G4PreCompoundTransitions.cc  useCEMtr="<<useCEMtr<<G4endl;
74
75  // Fermi energy
76  const G4double FermiEnergy = G4PreCompoundParameters::GetAddress()->GetFermiEnergy();
77 
78  // Nuclear radius
79  const G4double r0 = G4PreCompoundParameters::GetAddress()->GetTransitionsr0();
80   
81  // In order to calculate the level density parameter
82  // G4EvaporationLevelDensityParameter theLDP;
83
84  // Number of holes
85  G4double H = aFragment.GetNumberOfHoles();
86  // Number of Particles
87  G4double P = aFragment.GetNumberOfParticles();
88  // Number of Excitons
89  G4double N = P+H;
90  // Nucleus
91  G4double A = aFragment.GetA();
92  G4double Z = static_cast<G4double>(aFragment.GetZ());
93  G4double U = aFragment.GetExcitationEnergy();
94 
95  if(U<10*eV) return 0.0;
96 
97  //J. M. Quesada (Feb. 08) new physics
98  // OPT=1 Transitions are calculated according to Gudima's paper (original in G4PreCompound from VL)
99  // OPT=2 Transitions are calculated according to Gupta's formulae
100  //
101 
102 
103 
104  if (useCEMtr){
105
106   
107    // Relative Energy (T_{rel})
108    G4double RelativeEnergy = (8.0/5.0)*FermiEnergy + U/N;
109   
110    // Sample kind of nucleon-projectile
111    G4bool ChargedNucleon(false);
112    G4double chtest = 0.5;
113    if (P > 0) chtest = aFragment.GetNumberOfCharged()/P;
114    if (G4UniformRand() < chtest) ChargedNucleon = true;
115   
116    // Relative Velocity:
117    // <V_{rel}>^2
118    G4double RelativeVelocitySqr(0.0);
119    if (ChargedNucleon) RelativeVelocitySqr = 2.0*RelativeEnergy/proton_mass_c2;
120    else RelativeVelocitySqr = 2.0*RelativeEnergy/neutron_mass_c2;
121   
122    // <V_{rel}>
123    G4double RelativeVelocity = std::sqrt(RelativeVelocitySqr);
124   
125    // Proton-Proton Cross Section
126    G4double ppXSection = (10.63/RelativeVelocitySqr - 29.92/RelativeVelocity + 42.9)*millibarn;
127    // Proton-Neutron Cross Section
128    G4double npXSection = (34.10/RelativeVelocitySqr - 82.20/RelativeVelocity + 82.2)*millibarn;
129   
130    // Averaged Cross Section: \sigma(V_{rel})
131    //  G4double AveragedXSection = (ppXSection+npXSection)/2.0;
132    G4double AveragedXSection(0.0);
133    if (ChargedNucleon)
134      {
135        //JMQ: small bug fixed
136        //      AveragedXSection = ((Z-1.0) * ppXSection + (A-Z-1.0) * npXSection) / (A-1.0);
137        AveragedXSection = ((Z-1.0) * ppXSection + (A-Z) * npXSection) / (A-1.0);
138      }
139    else 
140      {
141        AveragedXSection = ((A-Z-1.0) * ppXSection + Z * npXSection) / (A-1.0);
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 - (7.0/5.0)*FermiRelRatio;
149    if (FermiRelRatio > 0.5) PauliFactor += (2.0/5.0)*FermiRelRatio*std::pow(2.0 - (1.0/FermiRelRatio), 5.0/2.0);
150   
151    // Interaction volume
152    //  G4double Vint = (4.0/3.0)*pi*std::pow(2.0*r0 + hbarc/(proton_mass_c2*RelativeVelocity) , 3.0);
153    G4double xx=2.0*r0 + hbarc/(proton_mass_c2*RelativeVelocity);
154    G4double Vint = (4.0/3.0)*pi*xx*xx*xx;
155   
156    // Transition probability for \Delta n = +2
157   
158    TransitionProb1 = AveragedXSection*PauliFactor*std::sqrt(2.0*RelativeEnergy/proton_mass_c2)/Vint;
159
160//JMQ 281009  phenomenological factor in order to increase equilibrium contribution
161//   G4double factor=5.0;
162//   TransitionProb1 *= factor;
163//
164    if (TransitionProb1 < 0.0) TransitionProb1 = 0.0; 
165   
166    G4double a = G4PreCompoundParameters::GetAddress()->GetLevelDensity();
167    // GE = g*E where E is Excitation Energy
168    G4double GE = (6.0/pi2)*a*A*U;
169   
170    G4double Fph = ((P*P+H*H+P-H)/4.0 - H/2.0);
171   
172    //G4bool NeverGoBack(false);
173    G4bool NeverGoBack;
174    if(useNGB)  NeverGoBack=true;
175    else NeverGoBack=false;
176   
177   
178    //JMQ/AH  bug fixed: if (U-Fph < 0.0) NeverGoBack = true;
179    if (GE-Fph < 0.0) NeverGoBack = true;
180   
181    // F(p+1,h+1)
182    G4double Fph1 = Fph + N/2.0;
183   
184    G4double ProbFactor = std::pow((GE-Fph)/(GE-Fph1),N+1.0);
185   
186   
187    if (NeverGoBack)
188      {
189      TransitionProb2 = 0.0;
190      TransitionProb3 = 0.0;
191      }
192    else 
193      {
194        // Transition probability for \Delta n = -2 (at F(p,h) = 0)
195        TransitionProb2 = TransitionProb1 * ProbFactor * (P*H*(N+1.0)*(N-2.0))/((GE-Fph)*(GE-Fph));
196        if (TransitionProb2 < 0.0) TransitionProb2 = 0.0; 
197       
198        // Transition probability for \Delta n = 0 (at F(p,h) = 0)
199        TransitionProb3 = TransitionProb1* ((N+1.0)/N) * ProbFactor  * (P*(P-1.0) + 4.0*P*H + H*(H-1.0))/(GE-Fph);
200        if (TransitionProb3 < 0.0) TransitionProb3 = 0.0; 
201      }
202   
203    //  G4cout<<"U = "<<U<<G4endl;
204    //  G4cout<<"N="<<N<<"  P="<<P<<"  H="<<H<<G4endl;
205    //  G4cout<<"l+ ="<<TransitionProb1<<"  l- ="<< TransitionProb2<<"  l0 ="<< TransitionProb3<<G4endl;
206    return TransitionProb1 + TransitionProb2 + TransitionProb3;
207  }
208 
209  else  {
210    //JMQ: Transition probabilities from Gupta's work
211   
212    G4double a = G4PreCompoundParameters::GetAddress()->GetLevelDensity();
213    // GE = g*E where E is Excitation Energy
214    G4double GE = (6.0/pi2)*a*A*U;
215 
216    G4double Kmfp=2.;
217       
218    TransitionProb1=1./Kmfp*3./8.*1./c_light*1.0e-9*(1.4e+21*U-2./(N+1)*6.0e+18*U*U);
219    if (TransitionProb1 < 0.0) TransitionProb1 = 0.0;
220   
221    if (useNGB){
222      TransitionProb2=0.;
223      TransitionProb3=0.;
224    }
225    else{       
226      if (N<=1) TransitionProb2=0. ;   
227      else  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);     
228      if (TransitionProb2 < 0.0) TransitionProb2 = 0.0; 
229      TransitionProb3=0.;
230    }
231   
232      //  G4cout<<"U = "<<U<<G4endl;
233    //  G4cout<<"N="<<N<<"  P="<<P<<"  H="<<H<<G4endl;
234    //  G4cout<<"l+ ="<<TransitionProb1<<"  l- ="<< TransitionProb2<<"  l0 ="<< TransitionProb3<<G4endl;
235    return TransitionProb1 + TransitionProb2 + TransitionProb3;
236  }
237}
238
239
240G4Fragment G4PreCompoundTransitions::PerformTransition(const G4Fragment & aFragment)
241{
242  G4Fragment result(aFragment);
243  G4double ChosenTransition = G4UniformRand()*(TransitionProb1 + TransitionProb2 + TransitionProb3);
244  G4int deltaN = 0;
245  G4int Nexcitons = result.GetNumberOfExcitons();
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 in proportion
258  // to the number charges w.r.t. number of particles,  PROVIDED that there are charged particles
259  if(deltaN < 0 && G4UniformRand() <= 
260     static_cast<G4double>(result.GetNumberOfCharged())/static_cast<G4double>(result.GetNumberOfParticles()) 
261     && (result.GetNumberOfCharged() >= 1)) {
262    result.SetNumberOfCharged(result.GetNumberOfCharged()+deltaN/2); // deltaN is negative!
263  }
264
265  // JMQ the following lines have to be before SetNumberOfCharged, otherwise the check on
266  // number of charged vs. number of particles fails
267  result.SetNumberOfParticles(result.GetNumberOfParticles()+deltaN/2);
268  result.SetNumberOfHoles(result.GetNumberOfHoles()+deltaN/2); 
269
270  // With weight Z/A, number of charged particles is increased with +1
271  if ( ( deltaN > 0 ) &&
272      (G4UniformRand() <= static_cast<G4double>(result.GetZ()-result.GetNumberOfCharged())/
273       std::max(static_cast<G4double>(result.GetA()-Nexcitons),1.)))
274    {
275      result.SetNumberOfCharged(result.GetNumberOfCharged()+deltaN/2);
276    }
277 
278  // Number of charged can not be greater that number of particles
279  if ( result.GetNumberOfParticles() < result.GetNumberOfCharged() ) 
280    {
281      result.SetNumberOfCharged(result.GetNumberOfParticles());
282    }
283 
284  return result;
285}
286
Note: See TracBrowser for help on using the repository browser.