| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
>> 27 // $Id: G4CompetitiveFission.cc,v 1.5 2006/06/29 20:13:33 gunter Exp $
>> 28 // GEANT4 tag $Name: geant4-09-01-patch-01 $
>> 29 //
27 // Hadronic Process: Nuclear De-excitations 30 // Hadronic Process: Nuclear De-excitations
28 // by V. Lara (Oct 1998) 31 // by V. Lara (Oct 1998)
29 // <<
30 // J. M. Quesada (March 2009). Bugs fixed: <<
31 // - Full relativistic calculation (L <<
32 // - Fission pairing energy is includ <<
33 // Now Energy and momentum are conserved in fi <<
34 32
35 #include "G4CompetitiveFission.hh" 33 #include "G4CompetitiveFission.hh"
36 #include "G4PairingCorrection.hh" 34 #include "G4PairingCorrection.hh"
37 #include "G4ParticleMomentum.hh" << 35
38 #include "G4NuclearLevelData.hh" << 36 G4CompetitiveFission::G4CompetitiveFission() : G4VEvaporationChannel("fission")
39 #include "G4VFissionBarrier.hh" << 37 {
40 #include "G4FissionBarrier.hh" << 38 theFissionBarrierPtr = new G4FissionBarrier;
41 #include "G4FissionProbability.hh" << 39 MyOwnFissionBarrier = true;
42 #include "G4VLevelDensityParameter.hh" << 40
43 #include "G4FissionLevelDensityParameter.hh" << 41 theFissionProbabilityPtr = new G4FissionProbability;
44 #include "G4Pow.hh" << 42 MyOwnFissionProbability = true;
45 #include "Randomize.hh" << 43
46 #include "G4RandomDirection.hh" << 44 theLevelDensityPtr = new G4FissionLevelDensityParameter;
47 #include "G4PhysicalConstants.hh" << 45 MyOwnLevelDensity = true;
48 #include "G4PhysicsModelCatalog.hh" << 46
49 << 47 MaximalKineticEnergy = -1000.0*MeV;
50 G4CompetitiveFission::G4CompetitiveFission() : << 48 FissionBarrier = 0.0;
51 { << 49 FissionProbability = 0.0;
52 theFissionBarrierPtr = new G4FissionBarrier; << 50 LevelDensityParameter = 0.0;
53 theFissionProbabilityPtr = new G4FissionProb << 51 }
54 theLevelDensityPtr = new G4FissionLevelDensi << 52
55 pairingCorrection = G4NuclearLevelData::GetI << 53 G4CompetitiveFission::G4CompetitiveFission(const G4CompetitiveFission &) : G4VEvaporationChannel()
56 theSecID = G4PhysicsModelCatalog::GetModelID << 54 {
57 } 55 }
58 56
59 G4CompetitiveFission::~G4CompetitiveFission() 57 G4CompetitiveFission::~G4CompetitiveFission()
60 { 58 {
61 if (myOwnFissionBarrier) delete theFissionBa << 59 if (MyOwnFissionBarrier) delete theFissionBarrierPtr;
62 if (myOwnFissionProbability) delete theFissi << 60
63 if (myOwnLevelDensity) delete theLevelDensit << 61 if (MyOwnFissionProbability) delete theFissionProbabilityPtr;
64 } << 62
65 << 63 if (MyOwnLevelDensity) delete theLevelDensityPtr;
66 void G4CompetitiveFission::Initialise() << 64 }
67 { << 65
68 if (!isInitialised) { << 66 const G4CompetitiveFission & G4CompetitiveFission::operator=(const G4CompetitiveFission &)
69 isInitialised = true; << 67 {
70 G4VEvaporationChannel::Initialise(); << 68 throw G4HadronicException(__FILE__, __LINE__, "G4CompetitiveFission::operator= meant to not be accessable");
71 if (OPTxs == 1) { fFactor = 0.5; } << 69 return *this;
72 } << 70 }
73 } << 71
74 << 72 G4bool G4CompetitiveFission::operator==(const G4CompetitiveFission &right) const
75 G4double G4CompetitiveFission::GetEmissionProb << 73 {
76 { << 74 return (this == (G4CompetitiveFission *) &right);
77 if (!isInitialised) { Initialise(); } <<
78 G4int Z = fragment->GetZ_asInt(); <<
79 G4int A = fragment->GetA_asInt(); <<
80 fissionProbability = 0.0; <<
81 // Saddle point excitation energy ---> A = 6 <<
82 if (A >= 65 && Z > 16) { <<
83 G4double exEnergy = fragment->GetExcitatio <<
84 pairingCorrection->GetFissionPairingCorr <<
85 <<
86 if (exEnergy > 0.0) { <<
87 fissionBarrier = theFissionBarrierPtr->F <<
88 maxKineticEnergy = exEnergy - fissionBar <<
89 fissionProbability = <<
90 theFissionProbabilityPtr->EmissionProbabilit <<
91 maxKineticEnergy); <<
92 } <<
93 } <<
94 return fissionProbability*fFactor; <<
95 } 75 }
96 76
97 G4Fragment* G4CompetitiveFission::EmittedFragm << 77 G4bool G4CompetitiveFission::operator!=(const G4CompetitiveFission &right) const
98 { 78 {
99 G4Fragment * Fragment1 = nullptr; << 79 return (this != (G4CompetitiveFission *) &right);
100 // Nucleus data << 80 }
101 // Atomic number of nucleus << 81
102 G4int A = theNucleus->GetA_asInt(); <<
103 // Charge of nucleus <<
104 G4int Z = theNucleus->GetZ_asInt(); <<
105 // Excitation energy (in MeV) <<
106 G4double U = theNucleus->GetExcitationEnergy <<
107 G4double pcorr = pairingCorrection->GetFissi <<
108 if (U <= pcorr) { return Fragment1; } <<
109 82
110 // Atomic Mass of Nucleus (in MeV) <<
111 G4double M = theNucleus->GetGroundStateMass( <<
112 83
113 // Nucleus Momentum <<
114 G4LorentzVector theNucleusMomentum = theNucl <<
115 84
116 // Calculate fission parameters << 85 void G4CompetitiveFission::Initialize(const G4Fragment & fragment)
117 theParam.DefineParameters(A, Z, U-pcorr, fis << 86 {
>> 87 G4int anA = static_cast<G4int>(fragment.GetA());
>> 88 G4int aZ = static_cast<G4int>(fragment.GetZ());
>> 89 G4double ExEnergy = fragment.GetExcitationEnergy() -
>> 90 G4PairingCorrection::GetInstance()->GetFissionPairingCorrection(anA,aZ);
118 91
119 // First fragment << 92
120 G4int A1 = 0; << 93 // Saddle point excitation energy ---> A = 65
121 G4int Z1 = 0; << 94 // Fission is excluded for A < 65
122 G4double M1 = 0.0; << 95 if (anA >= 65 && ExEnergy > 0.0) {
123 << 96 FissionBarrier = theFissionBarrierPtr->FissionBarrier(anA,aZ,ExEnergy);
124 // Second fragment << 97 MaximalKineticEnergy = ExEnergy - FissionBarrier;
125 G4int A2 = 0; << 98 LevelDensityParameter = theLevelDensityPtr->LevelDensityParameter(anA,aZ,ExEnergy);
126 G4int Z2 = 0; << 99 FissionProbability = theFissionProbabilityPtr->EmissionProbability(fragment,MaximalKineticEnergy);
127 G4double M2 = 0.0; <<
128 <<
129 G4double FragmentsExcitationEnergy = 0.0; <<
130 G4double FragmentsKineticEnergy = 0.0; <<
131 <<
132 G4int Trials = 0; <<
133 do { <<
134 <<
135 // First fragment <<
136 A1 = FissionAtomicNumber(A); <<
137 Z1 = FissionCharge(A, Z, A1); <<
138 M1 = G4NucleiProperties::GetNuclearMass(A1 <<
139 <<
140 // Second Fragment <<
141 A2 = A - A1; <<
142 Z2 = Z - Z1; <<
143 if (A2 < 1 || Z2 < 0 || Z2 > A2) { <<
144 FragmentsExcitationEnergy = -1.0; <<
145 continue; <<
146 } 100 }
147 M2 = G4NucleiProperties::GetNuclearMass(A2 << 101 else {
148 // Maximal Kinetic Energy (available energ << 102 MaximalKineticEnergy = -1000.0*MeV;
149 G4double Tmax = M + U - M1 - M2 - pcorr; << 103 LevelDensityParameter = 0.0;
150 << 104 FissionProbability = 0.0;
151 // Check that fragment masses are less or <<
152 if (Tmax < 0.0) { <<
153 FragmentsExcitationEnergy = -1.0; <<
154 continue; <<
155 } 105 }
156 106
157 FragmentsKineticEnergy = FissionKineticEne << 107 return;
158 A1, Z1, << 108 }
159 A2, Z2, <<
160 U , Tmax); <<
161 <<
162 // Excitation Energy <<
163 // FragmentsExcitationEnergy = Tmax - Frag <<
164 // JMQ 04/03/09 BUG FIXED: in order to ful <<
165 // fragments carry the fission pairing ene <<
166 // excitation energy <<
167 109
168 FragmentsExcitationEnergy = <<
169 Tmax - FragmentsKineticEnergy + pcorr; <<
170 110
171 // Loop checking, 05-Aug-2015, Vladimir Iv <<
172 } while (FragmentsExcitationEnergy < 0.0 && <<
173 <<
174 if (FragmentsExcitationEnergy <= 0.0) { <<
175 throw G4HadronicException(__FILE__, __LINE <<
176 "G4CompetitiveFission::BreakItUp: Excita <<
177 } <<
178 111
179 // Fragment 1 << 112 G4FragmentVector * G4CompetitiveFission::BreakUp(const G4Fragment & theNucleus)
180 M1 += FragmentsExcitationEnergy * A1/static_ << 113 {
181 // Fragment 2 << 114 // Nucleus data
182 M2 += FragmentsExcitationEnergy * A2/static_ << 115 // Atomic number of nucleus
183 // primary << 116 G4int A = static_cast<G4int>(theNucleus.GetA());
184 M += U; << 117 // Charge of nucleus
185 << 118 G4int Z = static_cast<G4int>(theNucleus.GetZ());
186 G4double etot1 = ((M - M2)*(M + M2) + M1*M1) << 119 // Excitation energy (in MeV)
187 G4ParticleMomentum Momentum1 = << 120 G4double U = theNucleus.GetExcitationEnergy() -
188 std::sqrt((etot1 - M1)*(etot1+M1))*G4Rando << 121 G4PairingCorrection::GetInstance()->GetFissionPairingCorrection(A,Z);
189 G4LorentzVector FourMomentum1(Momentum1, eto << 122 // Check that U > 0
190 FourMomentum1.boost(theNucleusMomentum.boost << 123 if (U <= 0.0) {
>> 124 G4FragmentVector * theResult = new G4FragmentVector;
>> 125 theResult->push_back(new G4Fragment(theNucleus));
>> 126 return theResult;
>> 127 }
>> 128
>> 129 // Atomic Mass of Nucleus (in MeV)
>> 130 G4double M = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z,A)/MeV;
>> 131 // Nucleus Momentum
>> 132 G4LorentzVector theNucleusMomentum = theNucleus.GetMomentum();
>> 133
>> 134 // Calculate fission parameters
>> 135 G4FissionParameters theParameters(A,Z,U,FissionBarrier);
>> 136
>> 137 // First fragment
>> 138 G4int A1 = 0;
>> 139 G4int Z1 = 0;
>> 140 G4double M1 = 0.0;
>> 141
>> 142 // Second fragment
>> 143 G4int A2 = 0;
>> 144 G4int Z2 = 0;
>> 145 G4double M2 = 0.0;
>> 146
>> 147 G4double FragmentsExcitationEnergy = 0.0;
>> 148 G4double FragmentsKineticEnergy = 0.0;
>> 149
>> 150 G4int Trials = 0;
>> 151 do {
>> 152
>> 153 // First fragment
>> 154 A1 = FissionAtomicNumber(A,theParameters);
>> 155 Z1 = FissionCharge(A,Z,A1);
>> 156 M1 = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z1,A1);
>> 157
>> 158
>> 159 // Second Fragment
>> 160 A2 = A - A1;
>> 161 Z2 = Z - Z1;
>> 162 if (A2 < 1 || Z2 < 0)
>> 163 throw G4HadronicException(__FILE__, __LINE__, "G4CompetitiveFission::BreakUp: Can't define second fragment! ");
>> 164 M2 = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z2,A2)/MeV;
>> 165
>> 166 // Check that fragment masses are less or equal than total energy
>> 167 // if (M1 + M2 > theNucleusMomentum.mag()/MeV)
>> 168 if (M1 + M2 > theNucleusMomentum.e()/MeV)
>> 169 throw G4HadronicException(__FILE__, __LINE__, "G4CompetitiveFission::BreakUp: Fragments Mass > Total Energy");
>> 170
>> 171 // Maximal Kinetic Energy (available energy for fragments)
>> 172 // G4double Tmax = theNucleusMomentum.mag()/MeV - M1 - M2;
>> 173 G4double Tmax = M + U - M1 - M2;
>> 174
>> 175 FragmentsKineticEnergy = FissionKineticEnergy( A , Z,
>> 176 A1, Z1,
>> 177 A2, Z2,
>> 178 U , Tmax,
>> 179 theParameters);
>> 180
>> 181 // Excitation Energy
>> 182 FragmentsExcitationEnergy = Tmax - FragmentsKineticEnergy;
191 183
192 // Create Fragments << 184 } while (FragmentsExcitationEnergy < 0.0 && Trials++ < 100);
193 Fragment1 = new G4Fragment( A1, Z1, FourMome <<
194 if (Fragment1 != nullptr) { Fragment1->SetCr <<
195 theNucleusMomentum -= FourMomentum1; <<
196 theNucleus->SetZandA_asInt(Z2, A2); <<
197 theNucleus->SetMomentum(theNucleusMomentum); <<
198 theNucleus->SetCreatorModelID(theSecID); <<
199 return Fragment1; <<
200 } <<
201 <<
202 G4int <<
203 G4CompetitiveFission::FissionAtomicNumber(G4in <<
204 // Calculates the atomic number of a fission <<
205 { <<
206 <<
207 // For Simplicity reading code <<
208 G4int A1 = theParam.GetA1(); <<
209 G4int A2 = theParam.GetA2(); <<
210 G4double As = theParam.GetAs(); <<
211 G4double Sigma2 = theParam.GetSigma2(); <<
212 G4double SigmaS = theParam.GetSigmaS(); <<
213 G4double w = theParam.GetW(); <<
214 185
215 G4double C2A = A2 + 3.72*Sigma2; <<
216 G4double C2S = As + 3.72*SigmaS; <<
217 186
218 G4double C2 = 0.0; <<
219 if (w > 1000.0 ) { C2 = C2S; } <<
220 else if (w < 0.001) { C2 = C2A; } <<
221 else { C2 = std::max(C2A,C2S <<
222 <<
223 G4double C1 = A-C2; <<
224 if (C1 < 30.0) { <<
225 C2 = A-30.0; <<
226 C1 = 30.0; <<
227 } <<
228 <<
229 G4double Am1 = (As + A1)*0.5; <<
230 G4double Am2 = (A1 + A2)*0.5; <<
231 <<
232 // Get Mass distributions as sum of symmetri <<
233 G4double Mass1 = MassDistribution(As,A); <<
234 G4double Mass2 = MassDistribution(Am1,A); <<
235 G4double Mass3 = MassDistribution(G4double(A <<
236 G4double Mass4 = MassDistribution(Am2,A); <<
237 G4double Mass5 = MassDistribution(G4double(A <<
238 // get maximal value among Mass1,...,Mass5 <<
239 G4double MassMax = Mass1; <<
240 if (Mass2 > MassMax) { MassMax = Mass2; } <<
241 if (Mass3 > MassMax) { MassMax = Mass3; } <<
242 if (Mass4 > MassMax) { MassMax = Mass4; } <<
243 if (Mass5 > MassMax) { MassMax = Mass5; } <<
244 <<
245 // Sample a fragment mass number, which lies <<
246 G4double xm; <<
247 G4double Pm; <<
248 do { <<
249 xm = C1+G4UniformRand()*(C2-C1); <<
250 Pm = MassDistribution(xm,A); <<
251 // Loop checking, 05-Aug-2015, Vladimir Iv <<
252 } while (MassMax*G4UniformRand() > Pm); <<
253 G4int ires = G4lrint(xm); <<
254 <<
255 return ires; <<
256 } <<
257 <<
258 G4double <<
259 G4CompetitiveFission::MassDistribution(G4doubl <<
260 // This method gives mass distribution F(x) <<
261 // which consist of symmetric and asymmetric <<
262 { <<
263 G4double y0 = (x-theParam.GetAs())/theParam. <<
264 G4double Xsym = LocalExp(y0); <<
265 <<
266 G4double y1 = (x - theParam.GetA1())/thePara <<
267 G4double y2 = (x - theParam.GetA2())/thePara <<
268 G4double z1 = (x - A + theParam.GetA1())/the <<
269 G4double z2 = (x - A + theParam.GetA2())/the <<
270 G4double Xasym = LocalExp(y1) + LocalExp(y2) <<
271 + 0.5*(LocalExp(z1) + LocalExp(z2)); <<
272 <<
273 G4double res; <<
274 G4double w = theParam.GetW(); <<
275 if (w > 1000) { res = Xsym; } <<
276 else if (w < 0.001) { res = Xasym; } <<
277 else { res = w*Xsym+Xasym; } <<
278 return res; <<
279 } <<
280 <<
281 G4int G4CompetitiveFission::FissionCharge(G4in <<
282 // Calculates the charge of a fission produc <<
283 { <<
284 static const G4double sigma = 0.6; <<
285 G4double DeltaZ = 0.0; <<
286 if (Af >= 134.0) { DeltaZ = -0.45; <<
287 else if (Af <= (A-134.0)) { DeltaZ = 0.45; } <<
288 else { DeltaZ = -0.45*( <<
289 187
290 G4double Zmean = (Af/A)*Z + DeltaZ; << 188 if (FragmentsExcitationEnergy <= 0.0)
291 << 189 throw G4HadronicException(__FILE__, __LINE__, "G4CompetitiveFission::BreakItUp: Excitation energy for fragments < 0.0!");
292 G4double theZ; << 190
293 do { << 191
294 theZ = G4RandGauss::shoot(Zmean,sigma); << 192 // while (FragmentsExcitationEnergy < 0 && Trials < 100);
295 // Loop checking, 05-Aug-2015, Vladimir Iv <<
296 } while (theZ < 1.0 || theZ > (Z-1.0) || th <<
297 193
298 return G4lrint(theZ); << 194 // Fragment 1
>> 195 G4double U1 = FragmentsExcitationEnergy * (static_cast<G4double>(A1)/static_cast<G4double>(A));
>> 196 // Fragment 2
>> 197 G4double U2 = FragmentsExcitationEnergy * (static_cast<G4double>(A2)/static_cast<G4double>(A));
>> 198
>> 199
>> 200 G4double Pmax = std::sqrt( 2 * ( ( (M1+U1)*(M2+U2) ) /
>> 201 ( (M1+U1)+(M2+U2) ) ) * FragmentsKineticEnergy);
>> 202
>> 203 G4ParticleMomentum momentum1 = IsotropicVector( Pmax );
>> 204 G4ParticleMomentum momentum2( -momentum1 );
>> 205
>> 206 // Perform a Galileo boost for fragments
>> 207 momentum1 += (theNucleusMomentum.boostVector() * (M1+U1));
>> 208 momentum2 += (theNucleusMomentum.boostVector() * (M2+U2));
>> 209
>> 210
>> 211 // Create 4-momentum for first fragment
>> 212 // Warning!! Energy conservation is broken
>> 213 G4LorentzVector FourMomentum1( momentum1 , std::sqrt(momentum1.mag2() + (M1+U1)*(M1+U1)));
>> 214
>> 215 // Create 4-momentum for second fragment
>> 216 // Warning!! Energy conservation is broken
>> 217 G4LorentzVector FourMomentum2( momentum2 , std::sqrt(momentum2.mag2() + (M2+U2)*(M2+U2)));
>> 218
>> 219 // Create Fragments
>> 220 G4Fragment * Fragment1 = new G4Fragment( A1, Z1, FourMomentum1);
>> 221 if (!Fragment1) throw G4HadronicException(__FILE__, __LINE__, "G4CompetitiveFission::BreakItUp: Can't create Fragment1! ");
>> 222 G4Fragment * Fragment2 = new G4Fragment( A2, Z2, FourMomentum2);
>> 223 if (!Fragment2) throw G4HadronicException(__FILE__, __LINE__, "G4CompetitiveFission::BreakItUp: Can't create Fragment2! ");
>> 224
>> 225 #ifdef PRECOMPOUND_TEST
>> 226 Fragment1->SetCreatorModel(G4String("G4CompetitiveFission"));
>> 227 Fragment2->SetCreatorModel(G4String("G4CompetitiveFission"));
>> 228 #endif
>> 229 // Create Fragment Vector
>> 230 G4FragmentVector * theResult = new G4FragmentVector;
>> 231
>> 232 theResult->push_back(Fragment1);
>> 233 theResult->push_back(Fragment2);
>> 234
>> 235 #ifdef debug
>> 236 CheckConservation(theNucleus,theResult);
>> 237 #endif
>> 238
>> 239 return theResult;
299 } 240 }
300 241
301 G4double << 242
302 G4CompetitiveFission::FissionKineticEnergy(G4i << 243
303 G4int Af1, G4int /*Zf1*/, << 244 G4int G4CompetitiveFission::FissionAtomicNumber(const G4int A, const G4FissionParameters & theParam)
304 G4int Af2, G4int /*Zf2*/, << 245 // Calculates the atomic number of a fission product
305 G4double /*U*/, G4double Tmax) << 246 {
306 // Gives the kinetic energy of fission produ << 247
307 { << 248 // For Simplicity reading code
308 // Find maximal value of A for fragments << 249 const G4double A1 = theParam.GetA1();
309 G4int AfMax = std::max(Af1,Af2); << 250 const G4double A2 = theParam.GetA2();
310 << 251 const G4double As = theParam.GetAs();
311 // Weights for symmetric and asymmetric comp << 252 // const G4double Sigma1 = theParam.GetSigma1();
312 G4double Pas = 0.0; << 253 const G4double Sigma2 = theParam.GetSigma2();
313 if (theParam.GetW() <= 1000) { << 254 const G4double SigmaS = theParam.GetSigmaS();
314 G4double x1 = (AfMax-theParam.GetA1())/the << 255 const G4double w = theParam.GetW();
315 G4double x2 = (AfMax-theParam.GetA2())/the << 256
316 Pas = 0.5*LocalExp(x1) + LocalExp(x2); << 257
317 } << 258 // G4double FasymAsym = 2.0*std::exp(-((A2-As)*(A2-As))/(2.0*Sigma2*Sigma2)) +
318 << 259 // std::exp(-((A1-As)*(A1-As))/(2.0*Sigma1*Sigma1));
319 G4double Ps = 0.0; << 260
320 if (theParam.GetW() >= 0.001) { << 261 // G4double FsymA1A2 = std::exp(-((As-(A1+A2))*(As-(A1+A2)))/(2.0*SigmaS*SigmaS));
321 G4double xs = (AfMax-theParam.GetAs())/the << 262
322 Ps = theParam.GetW()*LocalExp(xs); << 263
323 } << 264 G4double C2A = A2 + 3.72*Sigma2;
324 G4double Psy = (Pas + Ps > 0.0) ? Ps/(Pas+Ps << 265 G4double C2S = As + 3.72*SigmaS;
325 <<
326 // Fission fractions Xsy and Xas formed in s <<
327 G4double PPas = theParam.GetSigma1() + 2.0 * <<
328 G4double PPsy = theParam.GetW() * theParam.G <<
329 G4double Xas = (PPas + PPsy > 0.0) ? PPas/(P <<
330 G4double Xsy = 1.0 - Xas; <<
331 <<
332 // Average kinetic energy for symmetric and <<
333 G4double Eaverage = (0.1071*(Z*Z)/G4Pow::Get <<
334 <<
335 // Compute maximal average kinetic energy of <<
336 G4double TaverageAfMax; <<
337 G4double ESigma = 10*CLHEP::MeV; <<
338 // Select randomly fission mode (symmetric o <<
339 if (G4UniformRand() > Psy) { // Asymmetric M <<
340 G4double A11 = theParam.GetA1()-0.7979*the <<
341 G4double A12 = theParam.GetA1()+0.7979*the <<
342 G4double A21 = theParam.GetA2()-0.7979*the <<
343 G4double A22 = theParam.GetA2()+0.7979*the <<
344 // scale factor <<
345 G4double ScaleFactor = 0.5*theParam.GetSig <<
346 (AsymmetricRatio(A,A11)+AsymmetricRatio( <<
347 theParam.GetSigma2()*(AsymmetricRatio(A, <<
348 // Compute average kinetic energy for frag <<
349 TaverageAfMax = (Eaverage + 12.5 * Xsy) * <<
350 AsymmetricRatio(A,G4double(AfMax)); <<
351 <<
352 } else { // Symmetric Mode <<
353 G4double As0 = theParam.GetAs() + 0.7979*t <<
354 // Compute average kinetic energy for frag <<
355 TaverageAfMax = (Eaverage - 12.5*CLHEP::Me <<
356 *SymmetricRatio(A, G4double(AfMax))/Symm <<
357 ESigma = 8.0*CLHEP::MeV; <<
358 } <<
359 <<
360 // Select randomly, in accordance with Gauss <<
361 // fragment kinetic energy <<
362 G4double KineticEnergy; <<
363 G4int i = 0; <<
364 do { <<
365 KineticEnergy = G4RandGauss::shoot(Taverag <<
366 if (++i > 100) return Eaverage; <<
367 // Loop checking, 05-Aug-2015, Vladimir Iv <<
368 } while (KineticEnergy < Eaverage-3.72*ESigm <<
369 KineticEnergy > Eaverage+3.72*ESigma || <<
370 KineticEnergy > Tmax); <<
371 266
372 return KineticEnergy; << 267 G4double C2 = 0.0;
>> 268 if (w > 1000.0 ) C2 = C2S;
>> 269 else if (w < 0.001) C2 = C2A;
>> 270 else C2 = std::max(C2A,C2S);
>> 271
>> 272 G4double C1 = A-C2;
>> 273 if (C1 < 30.0) {
>> 274 C2 = A-30.0;
>> 275 C1 = 30.0;
>> 276 }
>> 277
>> 278 G4double Am1 = (As + A1)/2.0;
>> 279 G4double Am2 = (A1 + A2)/2.0;
>> 280
>> 281 // Get Mass distributions as sum of symmetric and asymmetric Gasussians
>> 282 G4double Mass1 = MassDistribution(As,A,theParam);
>> 283 G4double Mass2 = MassDistribution(Am1,A,theParam);
>> 284 G4double Mass3 = MassDistribution(A1,A,theParam);
>> 285 G4double Mass4 = MassDistribution(Am2,A,theParam);
>> 286 G4double Mass5 = MassDistribution(A2,A,theParam);
>> 287 // get maximal value among Mass1,...,Mass5
>> 288 G4double MassMax = Mass1;
>> 289 if (Mass2 > MassMax) MassMax = Mass2;
>> 290 if (Mass3 > MassMax) MassMax = Mass3;
>> 291 if (Mass4 > MassMax) MassMax = Mass4;
>> 292 if (Mass5 > MassMax) MassMax = Mass5;
>> 293
>> 294 // Sample a fragment mass number, which lies between C1 and C2
>> 295 G4double m;
>> 296 G4double Pm;
>> 297 do {
>> 298 m = C1+G4UniformRand()*(C2-C1);
>> 299 Pm = MassDistribution(m,A,theParam);
>> 300 } while (G4UniformRand() > Pm/MassMax);
>> 301
>> 302 return static_cast<G4int>(m+0.5);
373 } 303 }
374 304
375 void G4CompetitiveFission::SetFissionBarrier(G << 305
>> 306
>> 307
>> 308
>> 309
>> 310 G4double G4CompetitiveFission::MassDistribution(const G4double x, const G4double A,
>> 311 const G4FissionParameters & theParam)
>> 312 // This method gives mass distribution F(x) = F_{asym}(x)+w*F_{sym}(x)
>> 313 // which consist of symmetric and asymmetric sum of gaussians components.
376 { 314 {
377 if (myOwnFissionBarrier) delete theFissionBa << 315 G4double Xsym = std::exp(-0.5*(x-theParam.GetAs())*(x-theParam.GetAs())/
378 theFissionBarrierPtr = aBarrier; << 316 (theParam.GetSigmaS()*theParam.GetSigmaS()));
379 myOwnFissionBarrier = false; << 317
>> 318 G4double Xasym = std::exp(-0.5*(x-theParam.GetA2())*(x-theParam.GetA2())/
>> 319 (theParam.GetSigma2()*theParam.GetSigma2())) +
>> 320 std::exp(-0.5*(x-(A-theParam.GetA2()))*(x-(A-theParam.GetA2()))/
>> 321 (theParam.GetSigma2()*theParam.GetSigma2())) +
>> 322 0.5*std::exp(-0.5*(x-theParam.GetA1())*(x-theParam.GetA1())/
>> 323 (theParam.GetSigma1()*theParam.GetSigma1())) +
>> 324 0.5*std::exp(-0.5*(x-(A-theParam.GetA1()))*(x-(A-theParam.GetA1()))/
>> 325 (theParam.GetSigma1()*theParam.GetSigma1()));
>> 326
>> 327 if (theParam.GetW() > 1000) return Xsym;
>> 328 else if (theParam.GetW() < 0.001) return Xasym;
>> 329 else return theParam.GetW()*Xsym+Xasym;
380 } 330 }
381 331
382 void << 332
383 G4CompetitiveFission::SetEmissionStrategy(G4VE << 333
>> 334
>> 335 G4int G4CompetitiveFission::FissionCharge(const G4double A,
>> 336 const G4double Z,
>> 337 const G4double Af)
>> 338 // Calculates the charge of a fission product for a given atomic number Af
384 { 339 {
385 if (myOwnFissionProbability) delete theFissi << 340 const G4double sigma = 0.6;
386 theFissionProbabilityPtr = aFissionProb; << 341 G4double DeltaZ = 0.0;
387 myOwnFissionProbability = false; << 342 if (Af >= 134.0) DeltaZ = -0.45; // 134 <= Af
>> 343 else if (A <= (A-134.0)) DeltaZ = 0.45; // Af <= (A-134)
>> 344 else DeltaZ = -0.45*(Af-(A/2.0))/(134.0-(A/2.0)); // (A-134) < Af < 134
>> 345
>> 346 G4double Zmean = (Af/A)*Z + DeltaZ;
>> 347
>> 348 G4double theZ;
>> 349 do {
>> 350 theZ = G4RandGauss::shoot(Zmean,sigma);
>> 351 } while (theZ < 1.0 || theZ > (Z-1.0) || theZ > Af);
>> 352 // return static_cast<G4int>(theZ+0.5);
>> 353 return static_cast<G4int>(theZ+0.5);
388 } 354 }
389 355
390 void << 356
391 G4CompetitiveFission::SetLevelDensityParameter << 357
392 { << 358
393 if (myOwnLevelDensity) delete theLevelDensit << 359 G4double G4CompetitiveFission::FissionKineticEnergy(const G4double A, const G4double Z,
394 theLevelDensityPtr = aLevelDensity; << 360 const G4double Af1, const G4double /*Zf1*/,
395 myOwnLevelDensity = false; << 361 const G4double Af2, const G4double /*Zf2*/,
>> 362 const G4double /*U*/, const G4double Tmax,
>> 363 const G4FissionParameters & theParam)
>> 364 // Gives the kinetic energy of fission products
>> 365 {
>> 366 // Find maximal value of A for fragments
>> 367 G4double AfMax = std::max(Af1,Af2);
>> 368 if (AfMax < (A/2.0)) AfMax = A - AfMax;
>> 369
>> 370 // Weights for symmetric and asymmetric components
>> 371 G4double Pas;
>> 372 if (theParam.GetW() > 1000) Pas = 0.0;
>> 373 else {
>> 374 G4double P1 = 0.5*std::exp(-0.5*(AfMax-theParam.GetA1())*(AfMax-theParam.GetA1())/
>> 375 (theParam.GetSigma1()*theParam.GetSigma1()));
>> 376
>> 377 G4double P2 = std::exp(-0.5*(AfMax-theParam.GetA2())*(AfMax-theParam.GetA2())/
>> 378 (theParam.GetSigma2()*theParam.GetSigma2()));
>> 379
>> 380 Pas = P1+P2;
>> 381 }
>> 382
>> 383
>> 384 G4double Ps;
>> 385 if (theParam.GetW() < 0.001) Ps = 0.0;
>> 386 else
>> 387 Ps = theParam.GetW()*std::exp(-0.5*(AfMax-theParam.GetAs())*(AfMax-theParam.GetAs())/
>> 388 (theParam.GetSigmaS()*theParam.GetSigmaS()));
>> 389
>> 390
>> 391
>> 392 G4double Psy = Ps/(Pas+Ps);
>> 393
>> 394
>> 395 // Fission fractions Xsy and Xas formed in symmetric and asymmetric modes
>> 396 G4double PPas = theParam.GetSigma1() + 2.0 * theParam.GetSigma2();
>> 397 G4double PPsy = theParam.GetW() * theParam.GetSigmaS();
>> 398 G4double Xas = PPas / (PPas+PPsy);
>> 399 G4double Xsy = PPsy / (PPas+PPsy);
>> 400
>> 401
>> 402 // Average kinetic energy for symmetric and asymmetric components
>> 403 G4double Eaverage = 0.1071*MeV*(Z*Z)/std::pow(A,1.0/3.0) + 22.2*MeV;
>> 404
>> 405
>> 406 // Compute maximal average kinetic energy of fragments and Energy Dispersion (sqrt)
>> 407 G4double TaverageAfMax;
>> 408 G4double ESigma;
>> 409 // Select randomly fission mode (symmetric or asymmetric)
>> 410 if (G4UniformRand() > Psy) { // Asymmetric Mode
>> 411 G4double A11 = theParam.GetA1()-0.7979*theParam.GetSigma1();
>> 412 G4double A12 = theParam.GetA1()+0.7979*theParam.GetSigma1();
>> 413 G4double A21 = theParam.GetA2()-0.7979*theParam.GetSigma2();
>> 414 G4double A22 = theParam.GetA2()+0.7979*theParam.GetSigma2();
>> 415 // scale factor
>> 416 G4double ScaleFactor = 0.5*theParam.GetSigma1()*(AsymmetricRatio(A,A11)+AsymmetricRatio(A,A12))+
>> 417 theParam.GetSigma2()*(AsymmetricRatio(A,A21)+AsymmetricRatio(A,A22));
>> 418 // Compute average kinetic energy for fragment with AfMax
>> 419 TaverageAfMax = (Eaverage + 12.5 * Xsy) * (PPas/ScaleFactor) * AsymmetricRatio(A,AfMax);
>> 420 ESigma = 10.0*MeV; // MeV
>> 421
>> 422 } else { // Symmetric Mode
>> 423 G4double As0 = theParam.GetAs() + 0.7979*theParam.GetSigmaS();
>> 424 // scale factor
>> 425 G4double ScaleFactor = theParam.GetW()*theParam.GetSigmaS()*SymmetricRatio(A,As0);
>> 426 // Compute average kinetic energy for fragment with AfMax
>> 427 TaverageAfMax = (Eaverage - 12.5*MeV*Xas) * (PPsy/ScaleFactor) * SymmetricRatio(A,AfMax);
>> 428 ESigma = 8.0*MeV;
>> 429 }
>> 430
>> 431
>> 432 // Select randomly, in accordance with Gaussian distribution, fragment kinetic energy
>> 433 G4double KineticEnergy;
>> 434 G4int i = 0;
>> 435 do {
>> 436 KineticEnergy = G4RandGauss::shoot(TaverageAfMax,ESigma);
>> 437 if (i++ > 100) return Eaverage;
>> 438 } while (KineticEnergy < Eaverage-3.72*ESigma ||
>> 439 KineticEnergy > Eaverage+3.72*ESigma ||
>> 440 KineticEnergy > Tmax);
>> 441
>> 442 return KineticEnergy;
396 } 443 }
>> 444
>> 445
>> 446
>> 447
>> 448
>> 449 G4double G4CompetitiveFission::AsymmetricRatio(const G4double A,const G4double A11)
>> 450 {
>> 451 const G4double B1 = 23.5;
>> 452 const G4double A00 = 134.0;
>> 453 return Ratio(A,A11,B1,A00);
>> 454 }
>> 455
>> 456 G4double G4CompetitiveFission::SymmetricRatio(const G4double A,const G4double A11)
>> 457 {
>> 458 const G4double B1 = 5.32;
>> 459 const G4double A00 = A/2.0;
>> 460 return Ratio(A,A11,B1,A00);
>> 461 }
>> 462
>> 463 G4double G4CompetitiveFission::Ratio(const G4double A,const G4double A11,
>> 464 const G4double B1,const G4double A00)
>> 465 {
>> 466 if (A == 0) throw G4HadronicException(__FILE__, __LINE__, "G4CompetitiveFission::Ratio: A == 0!");
>> 467 if (A11 >= A/2.0 && A11 <= (A00+10.0)) return 1.0-B1*((A11-A00)/A)*((A11-A00)/A);
>> 468 else return 1.0-B1*(10.0/A)*(10.0/A)-2.0*(10.0/A)*B1*((A11-A00-10.0)/A);
>> 469 }
>> 470
>> 471
>> 472
>> 473
>> 474
>> 475 G4ThreeVector G4CompetitiveFission::IsotropicVector(const G4double Magnitude)
>> 476 // Samples a isotropic random vectorwith a magnitud given by Magnitude.
>> 477 // By default Magnitude = 1.0
>> 478 {
>> 479 G4double CosTheta = 1.0 - 2.0*G4UniformRand();
>> 480 G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
>> 481 G4double Phi = twopi*G4UniformRand();
>> 482 G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
>> 483 Magnitude*std::sin(Phi)*SinTheta,
>> 484 Magnitude*CosTheta);
>> 485 return Vector;
>> 486 }
>> 487
>> 488
>> 489 #ifdef debug
>> 490 void G4CompetitiveFission::CheckConservation(const G4Fragment & theInitialState,
>> 491 G4FragmentVector * Result) const
>> 492 {
>> 493 G4double ProductsEnergy =0;
>> 494 G4ThreeVector ProductsMomentum;
>> 495 G4int ProductsA = 0;
>> 496 G4int ProductsZ = 0;
>> 497 G4FragmentVector::iterator h;
>> 498 for (h = Result->begin(); h != Result->end(); h++) {
>> 499 G4LorentzVector tmp = (*h)->GetMomentum();
>> 500 ProductsEnergy += tmp.e();
>> 501 ProductsMomentum += tmp.vect();
>> 502 ProductsA += static_cast<G4int>((*h)->GetA());
>> 503 ProductsZ += static_cast<G4int>((*h)->GetZ());
>> 504 }
>> 505
>> 506 if (ProductsA != theInitialState.GetA()) {
>> 507 G4cout << "!!!!!!!!!! Baryonic Number Conservation Violation !!!!!!!!!!" << G4endl;
>> 508 G4cout << "G4CompetitiveFission.cc: Barionic Number Conservation test for fission fragments"
>> 509 << G4endl;
>> 510 G4cout << "Initial A = " << theInitialState.GetA()
>> 511 << " Fragments A = " << ProductsA << " Diference --> "
>> 512 << theInitialState.GetA() - ProductsA << G4endl;
>> 513 }
>> 514 if (ProductsZ != theInitialState.GetZ()) {
>> 515 G4cout << "!!!!!!!!!! Charge Conservation Violation !!!!!!!!!!" << G4endl;
>> 516 G4cout << "G4CompetitiveFission.cc: Charge Conservation test for fission fragments"
>> 517 << G4endl;
>> 518 G4cout << "Initial Z = " << theInitialState.GetZ()
>> 519 << " Fragments Z = " << ProductsZ << " Diference --> "
>> 520 << theInitialState.GetZ() - ProductsZ << G4endl;
>> 521 }
>> 522 if (std::abs(ProductsEnergy-theInitialState.GetMomentum().e()) > 1.0*keV) {
>> 523 G4cout << "!!!!!!!!!! Energy Conservation Violation !!!!!!!!!!" << G4endl;
>> 524 G4cout << "G4CompetitiveFission.cc: Energy Conservation test for fission fragments"
>> 525 << G4endl;
>> 526 G4cout << "Initial E = " << theInitialState.GetMomentum().e()/MeV << " MeV"
>> 527 << " Fragments E = " << ProductsEnergy/MeV << " MeV Diference --> "
>> 528 << (theInitialState.GetMomentum().e() - ProductsEnergy)/MeV << " MeV" << G4endl;
>> 529 }
>> 530 if (std::abs(ProductsMomentum.x()-theInitialState.GetMomentum().x()) > 1.0*keV ||
>> 531 std::abs(ProductsMomentum.y()-theInitialState.GetMomentum().y()) > 1.0*keV ||
>> 532 std::abs(ProductsMomentum.z()-theInitialState.GetMomentum().z()) > 1.0*keV) {
>> 533 G4cout << "!!!!!!!!!! Momentum Conservation Violation !!!!!!!!!!" << G4endl;
>> 534 G4cout << "G4CompetitiveFission.cc: Momentum Conservation test for fission fragments"
>> 535 << G4endl;
>> 536 G4cout << "Initial P = " << theInitialState.GetMomentum().vect() << " MeV"
>> 537 << " Fragments P = " << ProductsMomentum << " MeV Diference --> "
>> 538 << theInitialState.GetMomentum().vect() - ProductsMomentum << " MeV" << G4endl;
>> 539 }
>> 540 return;
>> 541 }
>> 542 #endif
>> 543
>> 544
>> 545
397 546
398 547