Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationProbability.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 // J.M. Quesada (August2008). Based on:
 27 //
 28 // Hadronic Process: Nuclear De-excitations
 29 // by V. Lara (Oct 1998)
 30 //
 31 // Modif (03 September 2008) by J. M. Quesada for external choice of inverse 
 32 // cross section option
 33 // JMQ (06 September 2008) Also external choices have been added for 
 34 // superimposed Coulomb barrier (if useSICB is set true, by default is false) 
 35 //
 36 // JMQ (14 february 2009) bug fixed in emission width: hbarc instead of 
 37 //                        hbar_Planck in the denominator
 38 //
 39 // V.Ivanchenko general clean-up since 2010
 40 //
 41 #include "G4EvaporationProbability.hh"
 42 #include "G4NuclearLevelData.hh"
 43 #include "G4VCoulombBarrier.hh"
 44 #include "G4PhysicalConstants.hh"
 45 #include "G4SystemOfUnits.hh"
 46 #include "G4PairingCorrection.hh"
 47 #include "G4NucleiProperties.hh"
 48 #include "G4KalbachCrossSection.hh"
 49 #include "G4ChatterjeeCrossSection.hh"
 50 #include "G4InterfaceToXS.hh"
 51 #include "G4IsotopeList.hh"
 52 #include "G4Neutron.hh"
 53 #include "G4Proton.hh"
 54 #include "G4Deuteron.hh"
 55 #include "G4Triton.hh"
 56 #include "G4He3.hh"
 57 #include "G4Alpha.hh"
 58 #include "Randomize.hh"
 59 #include "G4Exp.hh"
 60 #include "G4Log.hh"
 61 #include "G4Pow.hh"
 62 
 63 namespace
 64 {
 65   const G4double explim = 160.;
 66 }
 67 
 68 G4EvaporationProbability::G4EvaporationProbability(G4int anA, G4int aZ, 
 69                G4double aGamma) 
 70   : G4VEmissionProbability(aZ, anA), fGamma(aGamma)
 71 {
 72   resA13 = lastA = muu = freeU = a0 = a1 = delta0 = delta1 = 0.0;
 73   pcoeff = fGamma*pEvapMass*CLHEP::millibarn
 74     /((CLHEP::pi*CLHEP::hbarc)*(CLHEP::pi*CLHEP::hbarc)); 
 75 
 76   if (1 == theZ && 1 == theA) { index = 1; }
 77   else if (1 == theZ && 2 == theA) { index = 2; }
 78   else if (1 == theZ && 3 == theA) { index = 3; }
 79   else if (2 == theZ && 3 == theA) { index = 4; }
 80   else if (2 == theZ && 4 == theA) { index = 5; }
 81 
 82   if (OPTxs == 1) {
 83     const G4ParticleDefinition* part = nullptr;
 84     if (index == 1) { part = G4Proton::Proton(); }
 85     else if (index == 2) { part = G4Deuteron::Deuteron(); }
 86     else if (index == 3) { part = G4Triton::Triton(); }
 87     else if (index == 4) { part = G4He3::He3(); }
 88     else if (index == 5) { part = G4Alpha::Alpha(); }
 89     else { part = G4Neutron::Neutron(); }
 90     fXSection = new G4InterfaceToXS(part, index);
 91   }
 92   
 93   if (0 == aZ) {
 94     ResetIntegrator(30, 0.15*CLHEP::MeV, 0.02);
 95   } else {
 96     ResetIntegrator(30, 0.25*CLHEP::MeV, 0.03);
 97   }
 98 }
 99 
100 G4EvaporationProbability::~G4EvaporationProbability()
101 {
102   delete fXSection;
103 }
104 
105 G4double G4EvaporationProbability::CalcAlphaParam(const G4Fragment&)
106 {
107   return 1.0;
108 }
109  
110 G4double G4EvaporationProbability::CalcBetaParam(const G4Fragment&)
111 {
112   return 1.0;
113 }
114 
115 G4double G4EvaporationProbability::TotalProbability(
116   const G4Fragment& fragment, G4double minEnergy, G4double maxEnergy, 
117   G4double CB, G4double exEnergy)
118 {
119   G4int fragA = fragment.GetA_asInt();
120   G4int fragZ = fragment.GetZ_asInt(); 
121   freeU = exEnergy;
122   a0 = pNuclearLevelData->GetLevelDensity(fragZ, fragA, freeU);
123   delta0 = pNuclearLevelData->GetPairingCorrection(fragZ, fragA);
124   delta1 = pNuclearLevelData->GetPairingCorrection(resZ, resA);
125   resA13 = pG4pow->Z13(resA);
126   /*      
127   G4cout << "G4EvaporationProbability: Z= " << theZ << " A= " << theA 
128    << " resZ= " << resZ << " resA= " << resA 
129    << " fragZ= " << fragZ << " fragA= " << fragA 
130    << "\n   freeU= " << freeU  
131    << " a0= " << a0 << " OPT= " << OPTxs << " emin= " 
132          << minEnergy << " emax= " << maxEnergy 
133    << " CB= " << CB << G4endl;
134   */
135   if (OPTxs==0) {
136 
137     G4double SystemEntropy = 2.0*std::sqrt(a0*freeU);
138     const G4double RN2 = 2.25*CLHEP::fermi*CLHEP::fermi
139       /(CLHEP::twopi*CLHEP::hbar_Planck*hbar_Planck);
140 
141     G4double Alpha = CalcAlphaParam(fragment);
142     G4double Beta = CalcBetaParam(fragment);
143 
144     // to be checked where to use a0, where - a1  
145     a1 = pNuclearLevelData->GetLevelDensity(resZ,resA,freeU);
146     G4double GlobalFactor = fGamma*Alpha*pEvapMass*RN2*resA13*resA13/(a1*a1);
147     
148     G4double maxea = maxEnergy*a1;
149     G4double Term1 = Beta*a1 - 1.5 + maxea;
150     G4double Term2 = (2.0*Beta*a1-3.0)*std::sqrt(maxea) + 2*maxea;
151   
152     G4double ExpTerm1 = (SystemEntropy <= explim) ? G4Exp(-SystemEntropy) : 0.0;
153   
154     G4double ExpTerm2 = 2.*std::sqrt(maxea) - SystemEntropy;
155     ExpTerm2 = std::min(ExpTerm2, explim);
156     ExpTerm2 = G4Exp(ExpTerm2);
157   
158     pProbability = GlobalFactor*(Term1*ExpTerm1 + Term2*ExpTerm2);
159              
160   } else {
161     // if Coulomb barrier cutoff is superimposed for all cross sections 
162     // then the limit is the Coulomb Barrier
163     pProbability = IntegrateProbability(minEnergy, maxEnergy, CB);
164   }
165   /*
166   G4cout << "TotalProbability:  Emin=" << minEnergy << " Emax= " << maxEnergy 
167    << " CB= " << CB << " prob=" << pProbability << G4endl;
168   */
169   return pProbability;
170 }
171 
172 G4double G4EvaporationProbability::ComputeProbability(G4double K, G4double CB)
173 {
174   // abnormal case - should never happens
175   if(pMass < pEvapMass + pResMass + K) { return 0.0; }
176     
177   G4double pEvapM2 = pEvapMass*pEvapMass;
178   G4double mres = std::sqrt(pMass*pMass + pEvapM2 - 2.*pMass*(pEvapMass + K));
179 
180   G4double excRes = mres - pResMass;
181   if (excRes < 0.0) { return 0.0; }
182   G4double K1 = (pMass*(K + pEvapMass) - pEvapM2)/mres - pEvapMass;
183   K1 = std::max(K1, 0.0); 
184   G4double xs = CrossSection(K1, CB);
185   if (xs <= 0.0) { return 0.0; }
186 
187   a1 = pNuclearLevelData->GetLevelDensity(resZ, resA, excRes);
188   G4double E0 = std::max(freeU - delta0, 0.0);
189   G4double E1 = std::max(excRes - delta1, 0.0);
190   G4double prob = pcoeff*G4Exp(2.0*(std::sqrt(a1*E1) - std::sqrt(a0*E0)))*K1*xs;
191   return prob;
192 }
193 
194 G4double 
195 G4EvaporationProbability::CrossSection(G4double K, G4double CB)
196 {
197   // compute power once
198   if (OPTxs > 1 && 0 < index && resA != lastA) {
199     lastA = resA;
200     muu = G4KalbachCrossSection::ComputePowerParameter(resA, index);
201   }
202   if (OPTxs == 1) {
203     recentXS = fXSection->GetElementCrossSection(K, resZ)/CLHEP::millibarn;
204 
205   } else if (OPTxs == 2) { 
206     recentXS = G4ChatterjeeCrossSection::ComputeCrossSection(K, CB, resA13, muu,
207                                                              index, theZ, resA); 
208   } else if (OPTxs == 3) {
209     recentXS = G4KalbachCrossSection::ComputeCrossSection(K, CB, resA13, muu,
210                                                           index, theZ, theA, resA);
211   }
212   return recentXS;
213 }
214