ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PreCompoundTransitions.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4PreCompoundTransitions.cc
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 //
27 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4PreCompoundTransitions
33 //
34 // Author: V.Lara
35 //
36 // Modified:
37 // 16.02.2008 J.M.Quesada fixed bugs
38 // 06.09.2008 J.M.Quesada added external choices for:
39 // - "never go back" hipothesis (useNGB=true)
40 // - CEM transition probabilities (useCEMtr=true)
41 // 30.10.2009 J.M.Quesada: CEM transition probabilities have been renormalized
42 // (IAEA benchmark)
43 // 20.08.2010 V.Ivanchenko move constructor and destructor to the source and
44 // optimise the code
45 // 30.08.2011 M.Kelsey - Skip CalculateProbability if no excitons
46 
48 #include "G4PhysicalConstants.hh"
49 #include "G4SystemOfUnits.hh"
50 #include "Randomize.hh"
51 #include "G4NuclearLevelData.hh"
52 #include "G4DeexPrecoParameters.hh"
53 #include "G4Fragment.hh"
54 #include "G4Proton.hh"
55 #include "G4Exp.hh"
56 #include "G4Log.hh"
57 
59 {
63  FermiEnergy = param->GetFermiEnergy();
64  r0 = param->GetTransitionsR0();
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  TransitionProb2 = 0.0;
86  TransitionProb3 = 0.0;
87  /*
88  G4cout << "G4PreCompoundTransitions::CalculateProbability H/P/N/Z/A= "
89  << H << " " << P << " " << N << " " << Z << " " << A <<G4endl;
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  static const G4double sixdpi2 = 6.0/CLHEP::pi2;
100  G4double GE = sixdpi2*U*fNuclData->GetLevelDensity(Z,A,U);
101  if (useCEMtr) {
102  // Relative Energy (T_{rel})
103  G4double RelativeEnergy = 1.6*FermiEnergy + U/G4double(N);
104 
105  // Sample kind of nucleon-projectile
106  G4bool ChargedNucleon(false);
107  if(G4int(P*G4UniformRand()) <= aFragment.GetNumberOfCharged()) {
108  ChargedNucleon = true;
109  }
110 
111  // Relative Velocity:
112  // <V_{rel}>^2
113  G4double RelativeVelocitySqr;
114  if (ChargedNucleon) {
115  RelativeVelocitySqr = 2*RelativeEnergy/CLHEP::proton_mass_c2;
116  } else {
117  RelativeVelocitySqr = 2*RelativeEnergy/CLHEP::neutron_mass_c2;
118  }
119  // <V_{rel}>
120  G4double RelativeVelocity = std::sqrt(RelativeVelocitySqr);
121 
122  // Proton-Proton Cross Section
123  G4double ppXSection =
124  (10.63/RelativeVelocitySqr - 29.92/RelativeVelocity + 42.9)
126  // Proton-Neutron Cross Section
127  G4double npXSection =
128  (34.10/RelativeVelocitySqr - 82.20/RelativeVelocity + 82.2)
130 
131  // Averaged Cross Section: \sigma(V_{rel})
132  G4double AveragedXSection;
133  if (ChargedNucleon)
134  {
135  //JMQ: small bug fixed
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  }
142 
143  // Fermi relative energy ratio
144  G4double FermiRelRatio = FermiEnergy/RelativeEnergy;
145 
146  // This factor is introduced to take into account the Pauli principle
147  G4double PauliFactor = 1.0 - 1.4*FermiRelRatio;
148  if (FermiRelRatio > 0.5) {
149  G4double x = 2.0 - 1.0/FermiRelRatio;
150  PauliFactor += 0.4*FermiRelRatio*x*x*std::sqrt(x);
151  }
152  // Interaction volume
153  G4double xx = 2*r0 + CLHEP::hbarc/(CLHEP::proton_mass_c2*RelativeVelocity);
154  G4double Vint = CLHEP::pi*xx*xx*xx/0.75;
155 
156  // Transition probability for \Delta n = +2
157 
158  TransitionProb1 = std::max(0.0, AveragedXSection*PauliFactor
159  *std::sqrt(2.0*RelativeEnergy/CLHEP::proton_mass_c2)/Vint);
160 
161  //JMQ 281009 phenomenological factor in order to increase
162  // equilibrium contribution
163  // G4double factor=5.0;
164  // TransitionProb1 *= factor;
165 
166  // GE = g*E where E is Excitation Energy
167  G4double Fph = G4double(P*P+H*H+P-3*H)*0.25;
168 
169  if(!useNGB) {
170 
171  // F(p+1,h+1)
172  G4double Fph1 = Fph + N*0.5;
173 
174  static const G4double plimit = 100;
175 
176  //JMQ/AH bug fixed: if (U-Fph < 0.0)
177  if (GE > Fph1) {
178  G4double x0 = GE-Fph;
179  G4double x1 = (N+1)*G4Log(x0/(GE-Fph1));
180  if(x1 < plimit) {
181  x1 = G4Exp(x1)*TransitionProb1/x0;
182 
183  // Transition probability for \Delta n = -2 (at F(p,h) = 0)
184  TransitionProb2 = std::max(0.0, (P*H*(N+1)*(N-2))*x1/x0);
185 
186  // Transition probability for \Delta n = 0 (at F(p,h) = 0)
187  TransitionProb3 = std::max(0.0,((N+1)*(P*(P-1) + 4*P*H + H*(H-1)))*x1
188  /G4double(N));
189  }
190  }
191  }
192 
193  } else {
194  //JMQ: Transition probabilities from Gupta's work
195  // GE = g*E where E is Excitation Energy
196  TransitionProb1 = std::max(0.0, U*(4.2e+12 - 3.6e+10*U/G4double(N+1)))
197  /(16*CLHEP::c_light);
198 
199  if (!useNGB && N > 1) {
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