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