| Geant4 Cross Reference |
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 // Hadronic Process: Nuclear De-excitations
28 // by V. Lara (Oct 1998)
29 //
30 // J. M. Quesada (September 2009) bugs fixed in probability distribution for kinetic
31 // energy sampling:
32 // -hbarc instead of hbar_Planck (BIG BUG)
33 // -quantities for initial nucleus and residual are calculated separately
34 // V.Ivanchenko (September 2009) Added proper protection for unphysical final state
35 // J. M. Quesada (October 2009) fixed bug in CoulombBarrier calculation
36
37 #include "G4GEMChannel.hh"
38 #include "G4VCoulombBarrier.hh"
39 #include "G4GEMCoulombBarrier.hh"
40 #include "G4PairingCorrection.hh"
41 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4Pow.hh"
44 #include "G4Log.hh"
45 #include "G4Exp.hh"
46 #include "G4RandomDirection.hh"
47 #include "G4PhysicsModelCatalog.hh"
48
49 G4GEMChannel::G4GEMChannel(G4int theA, G4int theZ, const G4String & aName,
50 G4GEMProbability * aEmissionStrategy) :
51 G4VEvaporationChannel(aName),
52 A(theA),
53 Z(theZ),
54 EmissionProbability(0.0),
55 MaximalKineticEnergy(-CLHEP::GeV),
56 theEvaporationProbabilityPtr(aEmissionStrategy),
57 secID(-1)
58 {
59 theCoulombBarrierPtr = new G4GEMCoulombBarrier(theA, theZ);
60 theEvaporationProbabilityPtr->SetCoulomBarrier(theCoulombBarrierPtr);
61 theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
62 MyOwnLevelDensity = true;
63 EvaporatedMass = G4NucleiProperties::GetNuclearMass(A, Z);
64 ResidualMass = CoulombBarrier = 0.0;
65 fG4pow = G4Pow::GetInstance();
66 ResidualZ = ResidualA = 0;
67 fNucData = G4NuclearLevelData::GetInstance();
68 secID = G4PhysicsModelCatalog::GetModelID("model_G4GEMChannel");
69 }
70
71 G4GEMChannel::~G4GEMChannel()
72 {
73 if (MyOwnLevelDensity) { delete theLevelDensityPtr; }
74 delete theCoulombBarrierPtr;
75 }
76
77 G4double G4GEMChannel::GetEmissionProbability(G4Fragment* fragment)
78 {
79 G4int anA = fragment->GetA_asInt();
80 G4int aZ = fragment->GetZ_asInt();
81 ResidualA = anA - A;
82 ResidualZ = aZ - Z;
83 /*
84 G4cout << "G4GEMChannel: Z= " << Z << " A= " << A
85 << " FragmentZ= " << aZ << " FragmentA= " << anA
86 << " Zres= " << ResidualZ << " Ares= " << ResidualA
87 << G4endl;
88 */
89 // We only take into account channels which are physically allowed
90 EmissionProbability = 0.0;
91
92 // Only channels which are physically allowed are taken into account
93 if (ResidualA >= ResidualZ && ResidualZ >= 0 && ResidualA >= A) {
94
95 //Effective excitation energy
96 G4double ExEnergy = fragment->GetExcitationEnergy()
97 - fNucData->GetPairingCorrection(aZ, anA);
98 if(ExEnergy > 0.0) {
99 ResidualMass = G4NucleiProperties::GetNuclearMass(ResidualA, ResidualZ);
100 G4double FragmentMass = fragment->GetGroundStateMass();
101 G4double Etot = FragmentMass + ExEnergy;
102 // Coulomb Barrier calculation
103 CoulombBarrier =
104 theCoulombBarrierPtr->GetCoulombBarrier(ResidualA,ResidualZ,ExEnergy);
105 /*
106 G4cout << "Eexc(MeV)= " << ExEnergy/MeV
107 << " CoulBarrier(MeV)= " << CoulombBarrier/MeV << G4endl;
108 */
109 if(Etot > ResidualMass + EvaporatedMass + CoulombBarrier) {
110
111 // Maximal Kinetic Energy
112 MaximalKineticEnergy = ((Etot-ResidualMass)*(Etot+ResidualMass)
113 + EvaporatedMass*EvaporatedMass)/(2.0*Etot)
114 - EvaporatedMass - CoulombBarrier;
115
116 //G4cout << "CBarrier(MeV)= " << CoulombBarrier/MeV << G4endl;
117
118 if (MaximalKineticEnergy > 0.0) {
119 // Total emission probability for this channel
120 EmissionProbability = theEvaporationProbabilityPtr->
121 EmissionProbability(*fragment, MaximalKineticEnergy);
122 }
123 }
124 }
125 }
126 //G4cout << "Prob= " << EmissionProbability << G4endl;
127 return EmissionProbability;
128 }
129
130 G4Fragment* G4GEMChannel::EmittedFragment(G4Fragment* theNucleus)
131 {
132 G4Fragment* evFragment = 0;
133 G4double evEnergy = SampleKineticEnergy(*theNucleus) + EvaporatedMass;
134
135 G4ThreeVector momentum = G4RandomDirection()*
136 std::sqrt((evEnergy - EvaporatedMass)*(evEnergy + EvaporatedMass));
137
138 G4LorentzVector EvaporatedMomentum(momentum, evEnergy);
139 G4LorentzVector ResidualMomentum = theNucleus->GetMomentum();
140 EvaporatedMomentum.boost(ResidualMomentum.boostVector());
141
142 evFragment = new G4Fragment(A, Z, EvaporatedMomentum);
143 if ( evFragment != nullptr ) { evFragment->SetCreatorModelID(secID); }
144 ResidualMomentum -= EvaporatedMomentum;
145 theNucleus->SetZandA_asInt(ResidualZ, ResidualA);
146 theNucleus->SetMomentum(ResidualMomentum);
147 theNucleus->SetCreatorModelID(secID);
148
149 return evFragment;
150 }
151
152 G4double G4GEMChannel::SampleKineticEnergy(const G4Fragment & fragment)
153 // Samples fragment kinetic energy.
154 {
155 G4double U = fragment.GetExcitationEnergy();
156
157 G4double Alpha = theEvaporationProbabilityPtr->CalcAlphaParam(fragment);
158 G4double Beta = theEvaporationProbabilityPtr->CalcBetaParam(fragment);
159
160 // ***RESIDUAL***
161 //JMQ (September 2009) the following quantities refer to the RESIDUAL:
162 G4double delta0 = fNucData->GetPairingCorrection(ResidualZ,ResidualA);
163 G4double Ux = (2.5 + 150.0/ResidualA)*MeV;
164 G4double Ex = Ux + delta0;
165 G4double InitialLevelDensity;
166 // ***end RESIDUAL ***
167
168 // ***PARENT***
169 //JMQ (September 2009) the following quantities refer to the PARENT:
170
171 G4double deltaCN = fNucData->GetPairingCorrection(fragment.GetZ_asInt(),
172 fragment.GetA_asInt());
173 G4double aCN = theLevelDensityPtr->LevelDensityParameter(fragment.GetA_asInt(),
174 fragment.GetZ_asInt(),
175 U-deltaCN);
176 G4double UxCN = (2.5 + 150.0/G4double(fragment.GetA_asInt()))*MeV;
177 G4double ExCN = UxCN + deltaCN;
178 G4double TCN = 1.0/(std::sqrt(aCN/UxCN) - 1.5/UxCN);
179 // ***end PARENT***
180
181 //JMQ quantities calculated for CN in InitialLevelDensity
182 if ( U < ExCN )
183 {
184 G4double E0CN = ExCN - TCN*(G4Log(TCN/MeV) - G4Log(aCN*MeV)/4.0
185 - 1.25*G4Log(UxCN/MeV) + 2.0*std::sqrt(aCN*UxCN));
186 InitialLevelDensity = (pi/12.0)*G4Exp((U-E0CN)/TCN)/TCN;
187 }
188 else
189 {
190 G4double x = U-deltaCN;
191 G4double x1 = std::sqrt(aCN*x);
192 InitialLevelDensity = (pi/12.0)*G4Exp(2*x1)/(x*std::sqrt(x1));
193 }
194
195 G4double Spin = theEvaporationProbabilityPtr->GetSpin();
196 //JMQ BIG BUG fixed: hbarc instead of hbar_Planck !!!!
197 // it was also fixed in total probability
198 G4double gg = (2.0*Spin+1.0)*EvaporatedMass/(pi2* hbarc*hbarc);
199 //
200 //JMQ fix on Rb and geometrical cross sections according to Furihata's paper
201 // (JAERI-Data/Code 2001-105, p6)
202 G4double Rb = 0.0;
203 if (A > 4)
204 {
205 G4double Ad = fG4pow->Z13(ResidualA);
206 G4double Aj = fG4pow->Z13(A);
207 Rb = (1.12*(Aj + Ad) - 0.86*((Aj+Ad)/(Aj*Ad))+2.85)*fermi;
208 }
209 else if (A>1)
210 {
211 G4double Ad = fG4pow->Z13(ResidualA);
212 G4double Aj = fG4pow->Z13(A);
213 Rb=1.5*(Aj+Ad)*fermi;
214 }
215 else
216 {
217 G4double Ad = fG4pow->Z13(ResidualA);
218 Rb = 1.5*Ad*fermi;
219 }
220 G4double GeometricalXS = pi*Rb*Rb;
221
222 G4double ConstantFactor = gg*GeometricalXS*Alpha*pi/(InitialLevelDensity*12);
223 //JMQ : this is the assimptotic maximal kinetic energy of the ejectile
224 // (obtained by energy-momentum conservation).
225 // In general smaller than U-theSeparationEnergy
226 G4double theEnergy = MaximalKineticEnergy + CoulombBarrier;
227 G4double KineticEnergy;
228 G4double Probability;
229
230 for(G4int i=0; i<100; ++i) {
231 KineticEnergy = CoulombBarrier + G4UniformRand()*(MaximalKineticEnergy);
232 G4double edelta = theEnergy-KineticEnergy-delta0;
233 Probability = ConstantFactor*(KineticEnergy + Beta);
234 G4double a =
235 theLevelDensityPtr->LevelDensityParameter(ResidualA,ResidualZ,edelta);
236 G4double T = 1.0/(std::sqrt(a/Ux) - 1.5/Ux);
237 //JMQ fix in units
238
239 if (theEnergy - KineticEnergy < Ex) {
240 G4double E0 = Ex - T*(G4Log(T) - G4Log(a)*0.25
241 - 1.25*G4Log(Ux) + 2.0*std::sqrt(a*Ux));
242 Probability *= G4Exp((theEnergy-KineticEnergy-E0)/T)/T;
243 } else {
244 G4double e2 = edelta*edelta;
245 Probability *=
246 G4Exp(2*std::sqrt(a*edelta) - 0.25*G4Log(a*edelta*e2*e2));
247 }
248 if(EmissionProbability*G4UniformRand() <= Probability) { break; }
249 }
250
251 return KineticEnergy;
252 }
253
254 void G4GEMChannel::Dump() const
255 {
256 theEvaporationProbabilityPtr->Dump();
257 }
258
259
260
261