Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/de_excitation/fission/src/G4CompetitiveFission.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 (March 2009). Bugs fixed:
 31 //          - Full relativistic calculation (Lorentz boosts)
 32 //          - Fission pairing energy is included in fragment excitation energies
 33 // Now Energy and momentum are conserved in fission 
 34 
 35 #include "G4CompetitiveFission.hh"
 36 #include "G4PairingCorrection.hh"
 37 #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"
 45 #include "Randomize.hh"
 46 #include "G4RandomDirection.hh"
 47 #include "G4PhysicalConstants.hh"
 48 #include "G4PhysicsModelCatalog.hh"
 49 
 50 G4CompetitiveFission::G4CompetitiveFission() : G4VEvaporationChannel("fission"), theSecID(-1)
 51 {
 52   theFissionBarrierPtr = new G4FissionBarrier;
 53   theFissionProbabilityPtr = new G4FissionProbability;
 54   theLevelDensityPtr = new G4FissionLevelDensityParameter;
 55   pairingCorrection = G4NuclearLevelData::GetInstance()->GetPairingCorrection();
 56   theSecID = G4PhysicsModelCatalog::GetModelID("model_G4CompetitiveFission");
 57 }
 58 
 59 G4CompetitiveFission::~G4CompetitiveFission()
 60 {
 61   if (myOwnFissionBarrier) delete theFissionBarrierPtr;
 62   if (myOwnFissionProbability) delete theFissionProbabilityPtr;
 63   if (myOwnLevelDensity) delete theLevelDensityPtr;
 64 }
 65 
 66 void G4CompetitiveFission::Initialise()
 67 {
 68   if (!isInitialised) {
 69     isInitialised = true;
 70     G4VEvaporationChannel::Initialise();  
 71     if (OPTxs == 1) { fFactor = 0.5; }
 72   }
 73 }
 74 
 75 G4double G4CompetitiveFission::GetEmissionProbability(G4Fragment* fragment)
 76 {
 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 = 65
 82   if (A >= 65 && Z > 16) {
 83     G4double exEnergy = fragment->GetExcitationEnergy() - 
 84       pairingCorrection->GetFissionPairingCorrection(A, Z);
 85   
 86     if (exEnergy > 0.0) {
 87       fissionBarrier = theFissionBarrierPtr->FissionBarrier(A, Z, exEnergy);
 88       maxKineticEnergy = exEnergy - fissionBarrier;
 89       fissionProbability = 
 90   theFissionProbabilityPtr->EmissionProbability(*fragment,
 91                   maxKineticEnergy);
 92     }
 93   }
 94   return fissionProbability*fFactor;
 95 }
 96 
 97 G4Fragment* G4CompetitiveFission::EmittedFragment(G4Fragment* theNucleus)
 98 {
 99   G4Fragment * Fragment1 = nullptr; 
100   // Nucleus data
101   // Atomic number of nucleus
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->GetFissionPairingCorrection(A,Z);
108   if (U <= pcorr) { return Fragment1; }
109 
110   // Atomic Mass of Nucleus (in MeV)
111   G4double M = theNucleus->GetGroundStateMass();
112 
113   // Nucleus Momentum
114   G4LorentzVector theNucleusMomentum = theNucleus->GetMomentum();
115 
116   // Calculate fission parameters
117   theParam.DefineParameters(A, Z, U-pcorr, fissionBarrier);
118   
119   // First fragment
120   G4int A1 = 0;
121   G4int Z1 = 0;
122   G4double M1 = 0.0;
123 
124   // Second fragment
125   G4int A2 = 0;
126   G4int Z2 = 0;
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, Z1);
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     }
147     M2 = G4NucleiProperties::GetNuclearMass(A2, Z2);
148     // Maximal Kinetic Energy (available energy for fragments)
149     G4double Tmax = M + U - M1 - M2 - pcorr;
150 
151     // Check that fragment masses are less or equal than total energy
152     if (Tmax < 0.0) {
153       FragmentsExcitationEnergy = -1.0;
154       continue;
155     }
156 
157     FragmentsKineticEnergy = FissionKineticEnergy( A , Z,
158                A1, Z1,
159                A2, Z2,
160                U , Tmax);
161     
162     // Excitation Energy
163     // FragmentsExcitationEnergy = Tmax - FragmentsKineticEnergy;
164     // JMQ 04/03/09 BUG FIXED: in order to fulfill energy conservation the
165     // fragments carry the fission pairing energy in form of 
166     // excitation energy
167 
168     FragmentsExcitationEnergy = 
169       Tmax - FragmentsKineticEnergy + pcorr;
170 
171     // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
172   } while (FragmentsExcitationEnergy < 0.0 && ++Trials < 100);
173     
174   if (FragmentsExcitationEnergy <= 0.0) { 
175     throw G4HadronicException(__FILE__, __LINE__, 
176       "G4CompetitiveFission::BreakItUp: Excitation energy for fragments < 0.0!");
177   }
178 
179   // Fragment 1
180   M1 += FragmentsExcitationEnergy * A1/static_cast<G4double>(A);
181   // Fragment 2
182   M2 += FragmentsExcitationEnergy * A2/static_cast<G4double>(A);
183   // primary
184   M += U;
185 
186   G4double etot1 = ((M - M2)*(M + M2) + M1*M1)/(2*M);
187   G4ParticleMomentum Momentum1 = 
188     std::sqrt((etot1 - M1)*(etot1+M1))*G4RandomDirection();
189   G4LorentzVector FourMomentum1(Momentum1, etot1);
190   FourMomentum1.boost(theNucleusMomentum.boostVector());
191     
192   // Create Fragments
193   Fragment1 = new G4Fragment( A1, Z1, FourMomentum1);
194   if (Fragment1 != nullptr) { Fragment1->SetCreatorModelID(theSecID); }
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(G4int A)
204   // Calculates the atomic number of a fission product
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   
215   G4double C2A = A2 + 3.72*Sigma2;
216   G4double C2S = As + 3.72*SigmaS;
217   
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 symmetric and asymmetric Gasussians
233   G4double Mass1 = MassDistribution(As,A); 
234   G4double Mass2 = MassDistribution(Am1,A); 
235   G4double Mass3 = MassDistribution(G4double(A1),A); 
236   G4double Mass4 = MassDistribution(Am2,A); 
237   G4double Mass5 = MassDistribution(G4double(A2),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 between C1 and C2
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 Ivanchenko
252   } while (MassMax*G4UniformRand() > Pm);
253   G4int ires = G4lrint(xm);
254 
255   return ires;
256 }
257 
258 G4double 
259 G4CompetitiveFission::MassDistribution(G4double x, G4int A)
260   // This method gives mass distribution F(x) = F_{asym}(x)+w*F_{sym}(x)
261   // which consist of symmetric and asymmetric sum of gaussians components.
262 {
263   G4double y0 = (x-theParam.GetAs())/theParam.GetSigmaS();
264   G4double Xsym = LocalExp(y0);
265 
266   G4double y1 = (x - theParam.GetA1())/theParam.GetSigma1();
267   G4double y2 = (x - theParam.GetA2())/theParam.GetSigma2();
268   G4double z1 = (x - A + theParam.GetA1())/theParam.GetSigma1();
269   G4double z2 = (x - A + theParam.GetA2())/theParam.GetSigma2();
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(G4int A, G4int Z, G4double Af)
282   // Calculates the charge of a fission product for a given atomic number Af
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*(Af-A*0.5)/(134.0-A*0.5); }
289 
290   G4double Zmean = (Af/A)*Z + DeltaZ;
291  
292   G4double theZ;
293   do {
294     theZ = G4RandGauss::shoot(Zmean,sigma);
295     // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
296   } while (theZ  < 1.0 || theZ > (Z-1.0) || theZ > Af);
297   
298   return G4lrint(theZ);
299 }
300 
301 G4double 
302 G4CompetitiveFission::FissionKineticEnergy(G4int A, G4int Z,
303              G4int Af1, G4int /*Zf1*/,
304              G4int Af2, G4int /*Zf2*/,
305              G4double /*U*/, G4double Tmax)
306   // Gives the kinetic energy of fission products
307 {
308   // Find maximal value of A for fragments
309   G4int AfMax = std::max(Af1,Af2);
310 
311   // Weights for symmetric and asymmetric components
312   G4double Pas = 0.0;
313   if (theParam.GetW() <= 1000) { 
314     G4double x1 = (AfMax-theParam.GetA1())/theParam.GetSigma1();
315     G4double x2 = (AfMax-theParam.GetA2())/theParam.GetSigma2();
316     Pas = 0.5*LocalExp(x1) + LocalExp(x2);
317   }
318 
319   G4double Ps = 0.0;
320   if (theParam.GetW() >= 0.001) {
321     G4double xs = (AfMax-theParam.GetAs())/theParam.GetSigmaS();
322     Ps = theParam.GetW()*LocalExp(xs);
323   }
324   G4double Psy = (Pas + Ps > 0.0) ? Ps/(Pas+Ps) : 0.5;
325 
326   // Fission fractions Xsy and Xas formed in symmetric and asymmetric modes
327   G4double PPas = theParam.GetSigma1() + 2.0 * theParam.GetSigma2();
328   G4double PPsy = theParam.GetW() * theParam.GetSigmaS();
329   G4double Xas = (PPas + PPsy > 0.0) ? PPas/(PPas+PPsy) : 0.5;
330   G4double Xsy = 1.0 - Xas;
331 
332   // Average kinetic energy for symmetric and asymmetric components
333   G4double Eaverage = (0.1071*(Z*Z)/G4Pow::GetInstance()->Z13(A) + 22.2)*CLHEP::MeV;
334 
335   // Compute maximal average kinetic energy of fragments and Energy Dispersion 
336   G4double TaverageAfMax;
337   G4double ESigma = 10*CLHEP::MeV;
338   // Select randomly fission mode (symmetric or asymmetric)
339   if (G4UniformRand() > Psy) { // Asymmetric Mode
340     G4double A11 = theParam.GetA1()-0.7979*theParam.GetSigma1();
341     G4double A12 = theParam.GetA1()+0.7979*theParam.GetSigma1();
342     G4double A21 = theParam.GetA2()-0.7979*theParam.GetSigma2();
343     G4double A22 = theParam.GetA2()+0.7979*theParam.GetSigma2();
344     // scale factor
345     G4double ScaleFactor = 0.5*theParam.GetSigma1()*
346       (AsymmetricRatio(A,A11)+AsymmetricRatio(A,A12))+
347       theParam.GetSigma2()*(AsymmetricRatio(A,A21)+AsymmetricRatio(A,A22));
348     // Compute average kinetic energy for fragment with AfMax
349     TaverageAfMax = (Eaverage + 12.5 * Xsy) * (PPas/ScaleFactor) * 
350       AsymmetricRatio(A,G4double(AfMax));
351 
352   } else { // Symmetric Mode
353     G4double As0 = theParam.GetAs() + 0.7979*theParam.GetSigmaS();
354     // Compute average kinetic energy for fragment with AfMax
355     TaverageAfMax = (Eaverage - 12.5*CLHEP::MeV*Xas)
356       *SymmetricRatio(A, G4double(AfMax))/SymmetricRatio(A, As0);
357     ESigma = 8.0*CLHEP::MeV;
358   }
359 
360   // Select randomly, in accordance with Gaussian distribution, 
361   // fragment kinetic energy
362   G4double KineticEnergy;
363   G4int i = 0;
364   do {
365     KineticEnergy = G4RandGauss::shoot(TaverageAfMax, ESigma);
366     if (++i > 100) return Eaverage;
367     // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
368   } while (KineticEnergy < Eaverage-3.72*ESigma || 
369      KineticEnergy > Eaverage+3.72*ESigma ||
370      KineticEnergy > Tmax);
371   
372   return KineticEnergy;
373 }
374 
375 void G4CompetitiveFission::SetFissionBarrier(G4VFissionBarrier * aBarrier)
376 {
377   if (myOwnFissionBarrier) delete theFissionBarrierPtr;
378   theFissionBarrierPtr = aBarrier;
379   myOwnFissionBarrier = false;
380 }
381 
382 void 
383 G4CompetitiveFission::SetEmissionStrategy(G4VEmissionProbability * aFissionProb)
384 {
385   if (myOwnFissionProbability) delete theFissionProbabilityPtr;
386   theFissionProbabilityPtr = aFissionProb;
387   myOwnFissionProbability = false;
388 }
389 
390 void 
391 G4CompetitiveFission::SetLevelDensityParameter(G4VLevelDensityParameter* aLevelDensity)
392 { 
393   if (myOwnLevelDensity) delete theLevelDensityPtr;
394   theLevelDensityPtr = aLevelDensity;
395   myOwnLevelDensity = false;
396 }
397 
398