Geant4  10.00.p02
G4GEMChannel.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: G4GEMChannel.cc 74869 2013-10-23 09:26:17Z gcosmo $
27 //
28 // Hadronic Process: Nuclear De-excitations
29 // by V. Lara (Oct 1998)
30 //
31 // J. M. Quesada (September 2009) bugs fixed in probability distribution for kinetic
32 // energy sampling:
33 // -hbarc instead of hbar_Planck (BIG BUG)
34 // -quantities for initial nucleus and residual are calculated separately
35 // V.Ivanchenko (September 2009) Added proper protection for unphysical final state
36 // J. M. Quesada (October 2009) fixed bug in CoulombBarrier calculation
37 
38 #include "G4GEMChannel.hh"
39 #include "G4PairingCorrection.hh"
40 #include "G4PhysicalConstants.hh"
41 #include "G4SystemOfUnits.hh"
42 #include "G4Pow.hh"
43 #include "G4Log.hh"
44 #include "G4Exp.hh"
45 
46 G4GEMChannel::G4GEMChannel(G4int theA, G4int theZ, const G4String & aName,
47  G4GEMProbability * aEmissionStrategy,
48  G4VCoulombBarrier * aCoulombBarrier) :
49  G4VEvaporationChannel(aName),
50  A(theA),
51  Z(theZ),
52  theEvaporationProbabilityPtr(aEmissionStrategy),
53  theCoulombBarrierPtr(aCoulombBarrier),
54  EmissionProbability(0.0),
55  MaximalKineticEnergy(-CLHEP::GeV)
56 {
58  MyOwnLevelDensity = true;
62  ResidualZ = ResidualA = 0;
63 }
64 
66 {
67  if (MyOwnLevelDensity) { delete theLevelDensityPtr; }
68 }
69 
71 {
72  G4int anA = fragment->GetA_asInt();
73  G4int aZ = fragment->GetZ_asInt();
74  ResidualA = anA - A;
75  ResidualZ = aZ - Z;
76  //G4cout << "G4GEMChannel::Initialize: Z= " << aZ << " A= " << anA
77  // << " Zres= " << ResidualZ << " Ares= " << ResidualA << G4endl;
78 
79  // We only take into account channels which are physically allowed
80  if (ResidualA <= 0 || ResidualZ <= 0 || ResidualA < ResidualZ ||
81  (ResidualA == ResidualZ && ResidualA > 1))
82  {
83  CoulombBarrier = 0.0;
85  EmissionProbability = 0.0;
86  }
87  else
88  {
89  // Effective excitation energy
90  // JMQ 071009: pairing in ExEnergy should be the one of parent compound nucleus
91  // FIXED the bug causing reported crash by VI (negative Probabilities
92  // due to inconsistency in Coulomb barrier calculation (CoulombBarrier and -Beta
93  // param for protons must be the same)
94  // G4double ExEnergy = fragment.GetExcitationEnergy() -
95  // G4PairingCorrection::GetInstance()->GetPairingCorrection(ResidualA,ResidualZ);
96  G4double ExEnergy = fragment->GetExcitationEnergy() -
98 
99  //G4cout << "Eexc(MeV)= " << ExEnergy/MeV << G4endl;
100 
101  if( ExEnergy <= 0.0) {
102  CoulombBarrier = 0.0;
103  MaximalKineticEnergy = -1000.0*MeV;
104  EmissionProbability = 0.0;
105 
106  } else {
107 
109 
110  // Coulomb Barrier calculation
112  //G4cout << "CBarrier(MeV)= " << CoulombBarrier/MeV << G4endl;
113 
114  //Maximal kinetic energy (JMQ : at the Coulomb barrier)
116  CalcMaximalKineticEnergy(fragment->GetGroundStateMass()+ExEnergy);
117  //G4cout << "MaxE(MeV)= " << MaximalKineticEnergy/MeV << G4endl;
118 
119  // Emission probability
120  if (MaximalKineticEnergy <= 0.0)
121  {
122  EmissionProbability = 0.0;
123  }
124  else
125  {
126  // Total emission probability for this channel
130  }
131  }
132  }
133  //G4cout << "Prob= " << EmissionProbability << G4endl;
134  return EmissionProbability;
135 }
136 
138 {
139  G4double EvaporatedKineticEnergy = CalcKineticEnergy(theNucleus);
140  G4double EvaporatedEnergy = EvaporatedKineticEnergy + EvaporatedMass;
141 
142  G4ThreeVector momentum(IsotropicVector(std::sqrt(EvaporatedKineticEnergy*
143  (EvaporatedKineticEnergy+2.0*EvaporatedMass))));
144 
145  momentum.rotateUz(theNucleus.GetMomentum().vect().unit());
146 
147  G4LorentzVector EvaporatedMomentum(momentum,EvaporatedEnergy);
148  EvaporatedMomentum.boost(theNucleus.GetMomentum().boostVector());
149  G4Fragment * EvaporatedFragment = new G4Fragment(A,Z,EvaporatedMomentum);
150  // ** And now the residual nucleus **
151  G4double theExEnergy = theNucleus.GetExcitationEnergy();
152  G4double theMass = theNucleus.GetGroundStateMass();
153  G4double ResidualEnergy =
154  theMass + (theExEnergy - EvaporatedKineticEnergy) - EvaporatedMass;
155 
156  G4LorentzVector ResidualMomentum(-momentum,ResidualEnergy);
157  ResidualMomentum.boost(theNucleus.GetMomentum().boostVector());
158 
159  G4Fragment * ResidualFragment = new G4Fragment( ResidualA, ResidualZ, ResidualMomentum );
160 
161  G4FragmentVector * theResult = new G4FragmentVector;
162 
163  theResult->push_back(EvaporatedFragment);
164  theResult->push_back(ResidualFragment);
165  return theResult;
166 }
167 
169 // Calculate maximal kinetic energy that can be carried by fragment.
170 //JMQ this is not the assimptotic kinetic energy but the one at the Coulomb barrier
171 {
172  G4double T = (NucleusTotalE*NucleusTotalE + EvaporatedMass*EvaporatedMass - ResidualMass*ResidualMass)/
173  (2.0*NucleusTotalE) - EvaporatedMass - CoulombBarrier;
174 
175  return T;
176 }
177 
179 // Samples fragment kinetic energy.
180 {
181  G4double U = fragment.GetExcitationEnergy();
182 
185 
187 
188  // ***RESIDUAL***
189  //JMQ (September 2009) the following quantities refer to the RESIDUAL:
191  G4double Ux = (2.5 + 150.0/ResidualA)*MeV;
192  G4double Ex = Ux + delta0;
193  G4double InitialLevelDensity;
194  // ***end RESIDUAL ***
195 
196  // ***PARENT***
197  //JMQ (September 2009) the following quantities refer to the PARENT:
198 
199  G4double deltaCN =
201  fragment.GetZ_asInt());
203  fragment.GetZ_asInt(),U-deltaCN);
204  G4double UxCN = (2.5 + 150.0/G4double(fragment.GetA_asInt()))*MeV;
205  G4double ExCN = UxCN + deltaCN;
206  G4double TCN = 1.0/(std::sqrt(aCN/UxCN) - 1.5/UxCN);
207  // ***end PARENT***
208 
209  //JMQ quantities calculated for CN in InitialLevelDensity
210  if ( U < ExCN )
211  {
212  G4double E0CN = ExCN - TCN*(G4Log(TCN/MeV) - G4Log(aCN*MeV)/4.0
213  - 1.25*G4Log(UxCN/MeV) + 2.0*std::sqrt(aCN*UxCN));
214  InitialLevelDensity = (pi/12.0)*G4Exp((U-E0CN)/TCN)/TCN;
215  }
216  else
217  {
218  G4double x = U-deltaCN;
219  G4double x1 = std::sqrt(aCN*x);
220  InitialLevelDensity = (pi/12.0)*G4Exp(2*x1)/(x*std::sqrt(x1));
221  //InitialLevelDensity =
222  //(pi/12.0)*std::exp(2*std::sqrt(aCN*(U-deltaCN)))/std::pow(aCN*std::pow(U-deltaCN,5.0),1.0/4.0);
223  }
224 
226  //JMQ BIG BUG fixed: hbarc instead of hbar_Planck !!!!
227  // it was fixed in total probability (for this channel) but remained still here!!
228  // G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbar_Planck*hbar_Planck);
229  G4double gg = (2.0*Spin+1.0)*EvaporatedMass/(pi2* hbarc*hbarc);
230  //
231  //JMQ fix on Rb and geometrical cross sections according to Furihata's paper
232  // (JAERI-Data/Code 2001-105, p6)
233  G4double Rb = 0.0;
234  if (A > 4)
235  {
236  G4double Ad = fG4pow->Z13(ResidualA);
237  G4double Aj = fG4pow->Z13(A);
238  // RN = 1.12*(R1 + R2) - 0.86*((R1+R2)/(R1*R2));
239  Rb = 1.12*(Aj + Ad) - 0.86*((Aj+Ad)/(Aj*Ad))+2.85;
240  Rb *= fermi;
241  }
242  else if (A>1)
243  {
244  G4double Ad = fG4pow->Z13(ResidualA);
245  G4double Aj = fG4pow->Z13(A);
246  Rb=1.5*(Aj+Ad)*fermi;
247  }
248  else
249  {
250  G4double Ad = fG4pow->Z13(ResidualA);
251  Rb = 1.5*Ad*fermi;
252  }
253  // G4double GeometricalXS = pi*RN*RN*std::pow(G4double(fragment.GetA()-A),2./3.);
254  G4double GeometricalXS = pi*Rb*Rb;
255 
256  G4double ConstantFactor = gg*GeometricalXS*Alpha/InitialLevelDensity;
257  ConstantFactor *= pi/12.0;
258  //JMQ : this is the assimptotic maximal kinetic energy of the ejectile
259  // (obtained by energy-momentum conservation).
260  // In general smaller than U-theSeparationEnergy
262  G4double KineticEnergy;
263  G4double Probability;
264 
265  do
266  {
267  KineticEnergy = CoulombBarrier + G4UniformRand()*MaximalKineticEnergy;
268  Probability = ConstantFactor*(KineticEnergy + Beta);
269  G4double a =
270  theLevelDensityPtr->LevelDensityParameter(ResidualA,ResidualZ,theEnergy-KineticEnergy-delta0);
271  G4double T = 1.0/(std::sqrt(a/Ux) - 1.5/Ux);
272  //JMQ fix in units
273 
274  if ( theEnergy-KineticEnergy < Ex)
275  {
276  G4double E0 = Ex - T*(G4Log(T/MeV) - G4Log(a*MeV)/4.0
277  - 1.25*G4Log(Ux/MeV) + 2.0*std::sqrt(a*Ux));
278  Probability *= G4Exp((theEnergy-KineticEnergy-E0)/T)/T;
279  }
280  else
281  {
282  Probability *= G4Exp(2*std::sqrt(a*(theEnergy-KineticEnergy-delta0)))/
283  std::pow(a*fG4pow->powN(theEnergy-KineticEnergy-delta0,5), 0.25);
284  }
285  }
286  while (Normalization*G4UniformRand() > Probability);
287 
288  return KineticEnergy;
289 }
290 
291 
293  // Samples a isotropic random vectorwith a magnitud given by Magnitude.
294  // By default Magnitude = 1.0
295 {
296  G4double CosTheta = 1.0 - 2.0*G4UniformRand();
297  G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
298  G4double Phi = twopi*G4UniformRand();
299  G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
300  Magnitude*std::sin(Phi)*SinTheta,
301  Magnitude*CosTheta);
302  return Vector;
303 }
304 
305 
306 
static G4Pow * GetInstance()
Definition: G4Pow.cc:53
G4ThreeVector IsotropicVector(G4double Magnitude=1.0)
G4Pow * fG4pow
static const double MeV
Definition: G4SIunits.hh:193
static G4double GetNuclearMass(const G4double A, const G4double Z)
virtual ~G4GEMChannel()
Definition: G4GEMChannel.cc:65
G4double CalcKineticEnergy(const G4Fragment &fragment)
G4double powN(G4double x, G4int n) const
Definition: G4Pow.cc:125
CLHEP::Hep3Vector G4ThreeVector
G4double EmissionProbability(const G4Fragment &fragment, G4double anEnergy)
const G4double pi
G4double CalcAlphaParam(const G4Fragment &) const
G4double CoulombBarrier
G4double a
Definition: TRTMaterials.hh:39
G4bool MyOwnLevelDensity
int G4int
Definition: G4Types.hh:78
G4VCoulombBarrier * theCoulombBarrierPtr
G4double ResidualMass
#define G4UniformRand()
Definition: Randomize.hh:87
G4double Z13(G4int Z) const
Definition: G4Pow.hh:129
G4int GetA_asInt() const
Definition: G4Fragment.hh:238
G4double GetSpin(void) const
G4double GetNormalization(void) const
G4double EvaporatedMass
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:271
G4double GetPairingCorrection(G4int A, G4int Z) const
std::vector< G4Fragment * > G4FragmentVector
Definition: G4Fragment.hh:65
static const double GeV
Definition: G4SIunits.hh:196
G4double GetGroundStateMass() const
Definition: G4Fragment.hh:260
G4double EmissionProbability
G4double CalcBetaParam(const G4Fragment &) const
static const G4double A[nN]
G4GEMChannel(const G4int theA, const G4int theZ, const G4String &aName, G4GEMProbability *aEmissionStrategy, G4VCoulombBarrier *aCoulombBarrier)
Definition: G4GEMChannel.cc:46
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static G4PairingCorrection * GetInstance()
virtual G4double GetEmissionProbability(G4Fragment *theNucleus)
Definition: G4GEMChannel.cc:70
virtual G4FragmentVector * BreakUp(const G4Fragment &theNucleus)
G4VLevelDensityParameter * theLevelDensityPtr
virtual G4double LevelDensityParameter(G4int A, G4int Z, G4double U) const =0
G4int GetZ_asInt() const
Definition: G4Fragment.hh:243
virtual G4double GetCoulombBarrier(G4int ARes, G4int ZRes, G4double U) const =0
G4GEMProbability * theEvaporationProbabilityPtr
G4double CalcMaximalKineticEnergy(G4double U)
double G4double
Definition: G4Types.hh:76
G4double MaximalKineticEnergy
static const double fermi
Definition: G4SIunits.hh:93
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:255
CLHEP::HepLorentzVector G4LorentzVector