Geant4  10.00.p02
G4EvaporationChannel.cc
Go to the documentation of this file.
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 // $Id: G4EvaporationChannel.cc 74869 2013-10-23 09:26:17Z gcosmo $
27 //
28 //J.M. Quesada (August2008). Based on:
29 //
30 // Hadronic Process: Nuclear De-excitations
31 // by V. Lara (Oct 1998)
32 //
33 // Modified:
34 // 03-09-2008 J.M. Quesada for external choice of inverse cross section option
35 // 06-09-2008 J.M. Quesada Also external choices have been added for superimposed
36 // Coulomb barrier (if useSICB is set true, by default is false)
37 // 17-11-2010 V.Ivanchenko in constructor replace G4VEmissionProbability by
38 // G4EvaporationProbability and do not new and delete probability
39 // object at each call; use G4Pow
40 
41 #include "G4EvaporationChannel.hh"
42 #include "G4PairingCorrection.hh"
43 #include "G4NucleiProperties.hh"
44 #include "G4Pow.hh"
45 #include "G4Log.hh"
46 #include "G4Exp.hh"
48 #include "G4PhysicalConstants.hh"
49 #include "G4SystemOfUnits.hh"
50 #include "Randomize.hh"
51 #include "G4Alpha.hh"
52 
54  const G4String & aName,
55  G4EvaporationProbability* aEmissionStrategy,
56  G4VCoulombBarrier* aCoulombBarrier):
57  G4VEvaporationChannel(aName),
58  theA(anA),
59  theZ(aZ),
60  theEvaporationProbabilityPtr(aEmissionStrategy),
61  theCoulombBarrierPtr(aCoulombBarrier),
62  EmissionProbability(0.0),
63  MaximalKineticEnergy(-1000.0)
64 {
65  ResidualA = 0;
66  ResidualZ = 0;
70 }
71 
74  theA(0),
75  theZ(0),
76  theEvaporationProbabilityPtr(0),
77  theCoulombBarrierPtr(0),
78  EmissionProbability(0.0),
79  MaximalKineticEnergy(-1000.0)
80 {
81  ResidualA = 0;
82  ResidualZ = 0;
85 }
86 
88 {
89  delete theLevelDensityPtr;
90 }
91 
92 //void G4EvaporationChannel::Initialize(const G4Fragment &)
93 //{}
94 
96 {
97  //for inverse cross section choice
99  // for superimposed Coulomb Barrier for inverse cross sections
101 
102  G4int FragmentA = fragment->GetA_asInt();
103  G4int FragmentZ = fragment->GetZ_asInt();
104  ResidualA = FragmentA - theA;
105  ResidualZ = FragmentZ - theZ;
106  //G4cout << "G4EvaporationChannel::Initialize Z= " << theZ << " A= " << theA
107  // << " FragZ= " << FragmentZ << " FragA= " << FragmentA << G4endl;
108 
109  //Effective excitation energy
110  G4double ExEnergy = fragment->GetExcitationEnergy() -
112 
113  // Only channels which are physically allowed are taken into account
114  if (ResidualA <= 0 || ResidualZ <= 0 || ResidualA < ResidualZ ||
115  (ResidualA == ResidualZ && ResidualA > 1) || ExEnergy <= 0.0) {
117  MaximalKineticEnergy = -1000.0*MeV;
118  EmissionProbability = 0.0;
119  } else {
121  G4double FragmentMass = fragment->GetGroundStateMass();
122  CoulombBarrier =
124  // Maximal Kinetic Energy
125  MaximalKineticEnergy = CalcMaximalKineticEnergy(FragmentMass + ExEnergy);
126  //MaximalKineticEnergy = ExEnergy + fragment.GetGroundStateMass()
127  // - EvaporatedMass - ResidualMass;
128 
129  // Emission probability
130  // Protection for the case Tmax<V. If not set in this way we could end up in an
131  // infinite loop in the method GetKineticEnergy if OPTxs!=0 && useSICB=true.
132  // Of course for OPTxs=0 we have the Coulomb barrier
133 
134  G4double limit;
135  if (OPTxs==0 || (OPTxs!=0 && useSICB))
136  limit= CoulombBarrier;
137  else limit=0.;
138  limit = std::max(limit, FragmentMass - ResidualMass - EvaporatedMass);
139 
140  // The threshold for charged particle emission must be set to 0 if Coulomb
141  //cutoff is included in the cross sections
142  if (MaximalKineticEnergy <= limit) EmissionProbability = 0.0;
143  else {
144  // Total emission probability for this channel
147  }
148  }
149  //G4cout << "G4EvaporationChannel:: probability= " << EmissionProbability << G4endl;
150 
151  return EmissionProbability;
152 }
153 
155 {
156  /*
157  G4double Ecm = GetKineticEnergy(theNucleus) + ResidualMass + EvaporatedMass;
158 
159  G4double EvaporatedEnergy =
160  ((Ecm-ResidualMass)*(Ecm+ResidualMass) + EvaporatedMass*EvaporatedMass)/(2*Ecm);
161  */
162  G4double EvaporatedEnergy = GetKineticEnergy(theNucleus) + EvaporatedMass;
163 
165  (std::sqrt((EvaporatedEnergy - EvaporatedMass)*
166  (EvaporatedEnergy + EvaporatedMass))));
167 
168  G4LorentzVector EvaporatedMomentum(momentum,EvaporatedEnergy);
169  G4LorentzVector ResidualMomentum = theNucleus.GetMomentum();
170  EvaporatedMomentum.boost(ResidualMomentum.boostVector());
171 
172  G4Fragment * EvaporatedFragment = new G4Fragment(theA,theZ,EvaporatedMomentum);
173  ResidualMomentum -= EvaporatedMomentum;
174 
175  G4Fragment * ResidualFragment = new G4Fragment(ResidualA, ResidualZ, ResidualMomentum);
176 
177  G4FragmentVector * theResult = new G4FragmentVector;
178 
179 #ifdef debug
180  G4double Efinal = ResidualMomentum.e() + EvaporatedMomentum.e();
181  G4ThreeVector Pfinal = ResidualMomentum.vect() + EvaporatedMomentum.vect();
182  if (std::abs(Efinal-theNucleus.GetMomentum().e()) > 1.0*keV) {
183  G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4Evaporation Chanel: ENERGY @@@@@@@@@@@@@@@@"
184  << G4endl;
185  G4cout << "Initial : " << theNucleus.GetMomentum().e()/MeV << " MeV Final :"
186  <<Efinal/MeV << " MeV Delta: "
187  << (Efinal-theNucleus.GetMomentum().e())/MeV
188  << " MeV" << G4endl;
189  }
190  if (std::abs(Pfinal.x()-theNucleus.GetMomentum().x()) > 1.0*keV ||
191  std::abs(Pfinal.y()-theNucleus.GetMomentum().y()) > 1.0*keV ||
192  std::abs(Pfinal.z()-theNucleus.GetMomentum().z()) > 1.0*keV ) {
193  G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4Evaporation Chanel: MOMENTUM @@@@@@@@@@@@@@@@"
194  << G4endl;
195  G4cout << "Initial : " << theNucleus.GetMomentum().vect() << " MeV Final :"
196  <<Pfinal/MeV << " MeV Delta: " << Pfinal-theNucleus.GetMomentum().vect()
197  << " MeV" << G4endl;
198 
199  }
200 #endif
201  theResult->push_back(EvaporatedFragment);
202  theResult->push_back(ResidualFragment);
203  return theResult;
204 }
205 
207 // Calculates the maximal kinetic energy that can be carried by fragment.
209 {
210  // This is the "true" assimptotic kinetic energy (from energy conservation)
211  G4double Tmax =
212  ((NucleusTotalE-ResidualMass)*(NucleusTotalE+ResidualMass)
213  + EvaporatedMass*EvaporatedMass)/(2.0*NucleusTotalE) - EvaporatedMass;
214 
215  //JMQ (13-09-08) bug fixed: in the original version the Tmax is calculated
216  //at the Coulomb barrier
217  //IMPORTANT: meaning of Tmax differs in OPTxs=0 and OPTxs!=0
218  //When OPTxs!=0 Tmax is the TRUE (assimptotic) maximal kinetic energy
219 
220  if(OPTxs==0)
221  Tmax=Tmax- CoulombBarrier;
222 
223  return Tmax;
224 }
225 
227 //JMQ: New method for MC sampling of kinetic energy. Substitutes old CalcKineticEnergy
229 {
230 
231  if (OPTxs==0) {
232  // It uses Dostrovsky's approximation for the inverse reaction cross
233  // in the probability for fragment emission
234  // MaximalKineticEnergy energy in the original version (V.Lara) was calculated at
235  //the Coulomb barrier.
236 
237 
238  if (MaximalKineticEnergy < 0.0) {
239  throw G4HadronicException(__FILE__, __LINE__,
240  "G4EvaporationChannel::CalcKineticEnergy: maximal kinetic at the Coulomb barrier is less than 0");
241  }
242  G4double Rb = 4.0*theLevelDensityPtr->
243  LevelDensityParameter(ResidualA+theA,ResidualZ+theZ,MaximalKineticEnergy)*
245  G4double RbSqrt = std::sqrt(Rb);
246  G4double PEX1 = 0.0;
247  if (RbSqrt < 160.0) PEX1 = G4Exp(-RbSqrt);
248  G4double Rk = 0.0;
249  G4double FRk = 0.0;
250  do {
251  G4double RandNumber = G4UniformRand();
252  Rk = 1.0 + (1./RbSqrt)*G4Log(RandNumber + (1.0-RandNumber)*PEX1);
253  G4double Q1 = 1.0;
254  G4double Q2 = 1.0;
255  if (theZ == 0) { // for emitted neutron
256  G4double Beta = (2.12/G4Pow::GetInstance()->Z23(ResidualA) - 0.05)*MeV/
257  (0.76 + 2.2/G4Pow::GetInstance()->Z13(ResidualA));
258  Q1 = 1.0 + Beta/(MaximalKineticEnergy);
259  Q2 = Q1*std::sqrt(Q1);
260  }
261 
262  FRk = (3.0*std::sqrt(3.0)/2.0)/Q2 * Rk * (Q1 - Rk*Rk);
263 
264  } while (FRk < G4UniformRand());
265 
266  G4double result = MaximalKineticEnergy * (1.0-Rk*Rk) + CoulombBarrier;
267  return result;
268  } else if (OPTxs==1 || OPTxs==2 || OPTxs==3 || OPTxs==4) {
269 
270  // Coulomb barrier is just included in the cross sections
271  G4double V = 0;
272  if(useSICB) { V= CoulombBarrier; }
273 
275 
277  G4double T(0.0);
278  G4double NormalizedProbability(1.0);
279 
280  // VI: This is very ineffective - create new objects at each call to the method
281  /*
282  // A pointer is created in order to access the distribution function.
283  G4EvaporationProbability * G4EPtemp = 0;
284 
285  if (theA==1 && theZ==0) G4EPtemp=new G4NeutronEvaporationProbability();
286  else if (theA==1 && theZ==1) G4EPtemp=new G4ProtonEvaporationProbability();
287  else if (theA==2 && theZ==1 ) G4EPtemp=new G4DeuteronEvaporationProbability();
288  else if (theA==3 && theZ==1 ) G4EPtemp=new G4TritonEvaporationProbability();
289  else if (theA==3 && theZ==2 ) G4EPtemp=new G4He3EvaporationProbability();
290  else if (theA==4 && theZ==2) G4EPtemp=new G4AlphaEvaporationProbability();
291  else {
292  std::ostringstream errOs;
293  errOs << "ejected particle out of range in G4EvaporationChannel" << G4endl;
294  throw G4HadronicException(__FILE__, __LINE__, errOs.str());
295  }
296 
297  //for cross section selection and superimposed Coulom Barrier for xs
298  G4EPtemp->SetOPTxs(OPTxs);
299  G4EPtemp->UseSICB(useSICB);
300  */
301 
302  // use local pointer and not create a new one
303  do
304  {
305  T=V+G4UniformRand()*(Tmax-V);
306  NormalizedProbability =
308  GetEmissionProbability(const_cast<G4Fragment*>(&aFragment));
309 
310  }
311  while (G4UniformRand() > NormalizedProbability);
312  // delete G4EPtemp;
313  return T;
314  } else{
315  std::ostringstream errOs;
316  errOs << "Bad option for energy sampling in evaporation" <<G4endl;
317  throw G4HadronicException(__FILE__, __LINE__, errOs.str());
318  }
319 }
320 
322  // Samples a isotropic random vectorwith a magnitud given by Magnitude.
323  // By default Magnitude = 1.0
324 {
325  G4double CosTheta = 1.0 - 2.0*G4UniformRand();
326  G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
327  G4double Phi = twopi*G4UniformRand();
328  G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
329  Magnitude*std::sin(Phi)*SinTheta,
330  Magnitude*CosTheta);
331  return Vector;
332 }
333 
334 
335 
336 
337 
338 
339 
340 
341 
342 
static G4Pow * GetInstance()
Definition: G4Pow.cc:53
static const double MeV
Definition: G4SIunits.hh:193
static G4double GetNuclearMass(const G4double A, const G4double Z)
CLHEP::Hep3Vector G4ThreeVector
G4VCoulombBarrier * theCoulombBarrierPtr
G4EvaporationProbability * theEvaporationProbabilityPtr
virtual G4double GetEmissionProbability(G4Fragment *fragment)
int G4int
Definition: G4Types.hh:78
G4double ProbabilityDistributionFunction(const G4Fragment &aFragment, G4double K)
G4VLevelDensityParameter * theLevelDensityPtr
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
G4double Z13(G4int Z) const
Definition: G4Pow.hh:129
G4int GetA_asInt() const
Definition: G4Fragment.hh:238
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:271
G4double GetPairingCorrection(G4int A, G4int Z) const
std::vector< G4Fragment * > G4FragmentVector
Definition: G4Fragment.hh:65
G4double GetGroundStateMass() const
Definition: G4Fragment.hh:260
G4double CalcMaximalKineticEnergy(G4double U)
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static G4PairingCorrection * GetInstance()
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4FragmentVector * BreakUp(const G4Fragment &theNucleus)
G4int GetZ_asInt() const
Definition: G4Fragment.hh:243
G4double Z23(G4int Z) const
Definition: G4Pow.hh:153
virtual G4double GetCoulombBarrier(G4int ARes, G4int ZRes, G4double U) const =0
#define G4endl
Definition: G4ios.hh:61
G4ThreeVector IsotropicVector(G4double Magnitude=1.0)
static const double keV
Definition: G4SIunits.hh:195
double G4double
Definition: G4Types.hh:76
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:255
G4double GetKineticEnergy(const G4Fragment &aFragment)
CLHEP::HepLorentzVector G4LorentzVector