Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PreCompoundTransitions Class Reference

#include <G4PreCompoundTransitions.hh>

Inheritance diagram for G4PreCompoundTransitions:
Collaboration diagram for G4PreCompoundTransitions:

Public Member Functions

 G4PreCompoundTransitions ()
 
virtual ~G4PreCompoundTransitions ()
 
virtual G4double CalculateProbability (const G4Fragment &aFragment)
 
virtual void PerformTransition (G4Fragment &aFragment)
 
- Public Member Functions inherited from G4VPreCompoundTransitions
 G4VPreCompoundTransitions ()
 
virtual ~G4VPreCompoundTransitions ()
 
G4double GetTransitionProb1 () const
 
G4double GetTransitionProb2 () const
 
G4double GetTransitionProb3 () const
 
void UseNGB (G4bool use)
 
void UseCEMtr (G4bool use)
 

Additional Inherited Members

- Protected Attributes inherited from G4VPreCompoundTransitions
G4bool useNGB
 
G4bool useCEMtr
 
G4double TransitionProb1
 
G4double TransitionProb2
 
G4double TransitionProb3
 

Detailed Description

Definition at line 50 of file G4PreCompoundTransitions.hh.

Constructor & Destructor Documentation

G4PreCompoundTransitions::G4PreCompoundTransitions ( )

Definition at line 59 of file G4PreCompoundTransitions.cc.

60 {
61  proton = G4Proton::Proton();
62  G4DeexPrecoParameters* param =
64  FermiEnergy = param->GetFermiEnergy();
65  r0 = param->GetTransitionsR0();
66  aLDP = param->GetLevelDensity();
67 }
G4double GetLevelDensity() const
G4double GetFermiEnergy() const
G4DeexPrecoParameters * GetParameters()
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4double GetTransitionsR0() const
static G4NuclearLevelData * GetInstance()

Here is the call graph for this function:

G4PreCompoundTransitions::~G4PreCompoundTransitions ( )
virtual

Definition at line 69 of file G4PreCompoundTransitions.cc.

70 {}

Member Function Documentation

G4double G4PreCompoundTransitions::CalculateProbability ( const G4Fragment aFragment)
virtual

Implements G4VPreCompoundTransitions.

Definition at line 75 of file G4PreCompoundTransitions.cc.

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 << "G4PreCompoundTransitions::CalculateProbability H/P/N/Z/A= "
91  << H << " " << P << " " << N << " " << Z << " " << A <<G4endl;
92  G4cout << aFragment << G4endl;
93  */
94  if(U < 10*eV || 0==N) { return 0.0; }
95 
96  //J. M. Quesada (Feb. 08) new physics
97  // OPT=1 Transitions are calculated according to Gudima's paper
98  // (original in G4PreCompound from VL)
99  // OPT=2 Transitions are calculated according to Gupta's formulae
100  //
101  static const G4double sixdpi2 = 6.0/CLHEP::pi2;
102  if (useCEMtr) {
103  // Relative Energy (T_{rel})
104  G4double RelativeEnergy = 1.6*FermiEnergy + U/G4double(N);
105 
106  // Sample kind of nucleon-projectile
107  G4bool ChargedNucleon(false);
108  if(G4int(P*G4UniformRand()) <= aFragment.GetNumberOfCharged()) {
109  ChargedNucleon = true;
110  }
111 
112  // Relative Velocity:
113  // <V_{rel}>^2
114  G4double RelativeVelocitySqr;
115  if (ChargedNucleon) {
116  RelativeVelocitySqr = 2*RelativeEnergy/CLHEP::proton_mass_c2;
117  } else {
118  RelativeVelocitySqr = 2*RelativeEnergy/CLHEP::neutron_mass_c2;
119  }
120  // <V_{rel}>
121  G4double RelativeVelocity = std::sqrt(RelativeVelocitySqr);
122 
123  // Proton-Proton Cross Section
124  G4double ppXSection =
125  (10.63/RelativeVelocitySqr - 29.92/RelativeVelocity + 42.9)
127  // Proton-Neutron Cross Section
128  G4double npXSection =
129  (34.10/RelativeVelocitySqr - 82.20/RelativeVelocity + 82.2)
131 
132  // Averaged Cross Section: \sigma(V_{rel})
133  G4double AveragedXSection;
134  if (ChargedNucleon)
135  {
136  //JMQ: small bug fixed
137  AveragedXSection = ((Z-1)*ppXSection + (A-Z)*npXSection)/G4double(A-1);
138  }
139  else
140  {
141  AveragedXSection = ((A-Z-1)*ppXSection + Z*npXSection)/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  }
153  // Interaction volume
154  G4double xx = 2*r0 + CLHEP::hbarc/(CLHEP::proton_mass_c2*RelativeVelocity);
155  G4double Vint = CLHEP::pi*xx*xx*xx/0.75;
156 
157  // Transition probability for \Delta n = +2
158 
159  TransitionProb1 = std::max(0.0, AveragedXSection*PauliFactor
160  *std::sqrt(2.0*RelativeEnergy/CLHEP::proton_mass_c2)/Vint);
161 
162  //JMQ 281009 phenomenological factor in order to increase
163  // equilibrium contribution
164  // G4double factor=5.0;
165  // TransitionProb1 *= factor;
166 
167  // GE = g*E where E is Excitation Energy
168  G4double GE = sixdpi2*aLDP*A*U;
169  G4double Fph = G4double(P*P+H*H+P-3*H)*0.25;
170 
171  if(!useNGB) {
172 
173  // F(p+1,h+1)
174  G4double Fph1 = Fph + N*0.5;
175 
176  static const G4double plimit = 100;
177 
178  //JMQ/AH bug fixed: if (U-Fph < 0.0)
179  if (GE-Fph1 > 0.0) {
180  G4double x0 = GE-Fph;
181  G4double x1 = (N+1)*G4Log(x0/(GE-Fph1));
182  if(x1 < plimit) {
183  x1 = G4Exp(x1)*TransitionProb1/x0;
184 
185  // Transition probability for \Delta n = -2 (at F(p,h) = 0)
186  TransitionProb2 = std::max(0.0, (P*H*(N+1)*(N-2))*x1/x0);
187 
188  // Transition probability for \Delta n = 0 (at F(p,h) = 0)
189  TransitionProb3 = std::max(0.0,((N+1)*(P*(P-1) + 4*P*H + H*(H-1)))*x1
190  /G4double(N));
191  }
192  }
193  }
194 
195  } else {
196  //JMQ: Transition probabilities from Gupta's work
197  // GE = g*E where E is Excitation Energy
198  TransitionProb1 = std::max(0.0, U*(4.2e+12 - 3.6e+10*U/G4double(N+1)))
199  /(16*CLHEP::c_light);
200 
201  if (!useNGB && N > 1) {
202  G4double GE = sixdpi2*aLDP*A*U;
203  TransitionProb2 = ((N-1)*(N-2)*P*H)*TransitionProb1/(GE*GE);
204  }
205  }
206  // G4cout<<"U = "<<U<<G4endl;
207  // G4cout<<"N="<<N<<" P="<<P<<" H="<<H<<G4endl;
208  // G4cout<<"l+ ="<<TransitionProb1<<" l- ="<< TransitionProb2
209  // <<" l0 ="<< TransitionProb3<<G4endl;
211 }
static constexpr double proton_mass_c2
static constexpr double pi2
Definition: SystemOfUnits.h:57
static constexpr double millibarn
Definition: SystemOfUnits.h:86
static constexpr double hbarc
G4int GetNumberOfParticles() const
Definition: G4Fragment.hh:345
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
static double P[]
G4int GetNumberOfHoles() const
Definition: G4Fragment.hh:365
#define G4UniformRand()
Definition: Randomize.hh:97
double A(double temperature)
G4int GetA_asInt() const
Definition: G4Fragment.hh:266
bool G4bool
Definition: G4Types.hh:79
Definition: Evaluator.cc:66
static constexpr double neutron_mass_c2
static constexpr double eV
Definition: G4SIunits.hh:215
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
T max(const T t1, const T t2)
brief Return the largest of the two arguments
static constexpr double c_light
G4int GetZ_asInt() const
Definition: G4Fragment.hh:271
**D E S C R I P T I O N
Definition: HEPEvtcom.cc:77
double G4double
Definition: G4Types.hh:76
G4int GetNumberOfCharged() const
Definition: G4Fragment.hh:350
static constexpr double pi
Definition: SystemOfUnits.h:54
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:283

Here is the call graph for this function:

void G4PreCompoundTransitions::PerformTransition ( G4Fragment aFragment)
virtual

Implements G4VPreCompoundTransitions.

Definition at line 213 of file G4PreCompoundTransitions.cc.

214 {
215  G4double ChosenTransition =
217  G4int deltaN = 0;
218  G4int Npart = result.GetNumberOfParticles();
219  G4int Ncharged = result.GetNumberOfCharged();
220  G4int Nholes = result.GetNumberOfHoles();
221  if (ChosenTransition <= TransitionProb1)
222  {
223  // Number of excitons is increased on \Delta n = +2
224  deltaN = 2;
225  }
226  else if (ChosenTransition <= TransitionProb1+TransitionProb2)
227  {
228  // Number of excitons is increased on \Delta n = -2
229  deltaN = -2;
230  }
231 
232  // AH/JMQ: Randomly decrease the number of charges if deltaN is -2 and
233  // in proportion to the number charges w.r.t. number of particles,
234  // PROVIDED that there are charged particles
235  deltaN /= 2;
236 
237  //G4cout << "deltaN= " << deltaN << G4endl;
238 
239  // JMQ the following lines have to be before SetNumberOfCharged,
240  // otherwise the check on number of charged vs. number of particles fails
241  result.SetNumberOfParticles(Npart+deltaN);
242  result.SetNumberOfHoles(Nholes+deltaN);
243 
244  if(deltaN < 0) {
245  if( (Ncharged == Npart) ||
246  (Ncharged >= 1 && G4int(Npart*G4UniformRand()) <= Ncharged))
247  {
248  result.SetNumberOfCharged(Ncharged+deltaN); // deltaN is negative!
249  }
250 
251  } else if ( deltaN > 0 ) {
252  // With weight Z/A, number of charged particles is increased with +1
253  G4int A = result.GetA_asInt() - Npart;
254  G4int Z = result.GetZ_asInt() - Ncharged;
255  if((Z == A) || (Z > 0 && G4int(A*G4UniformRand()) <= Z))
256  {
257  result.SetNumberOfCharged(Ncharged+deltaN);
258  }
259  }
260 
261  // Number of charged can not be greater that number of particles
262  if ( Npart < Ncharged )
263  {
264  result.SetNumberOfCharged(Npart);
265  }
266  //G4cout << "### After transition" << G4endl;
267  //G4cout << result << G4endl;
268 }
G4double G4ParticleHPJENDLHEData::G4double result
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:97
double A(double temperature)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:


The documentation for this class was generated from the following files: