Geant4  10.02.p01
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 91834 2015-08-07 07:24:22Z 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;
64 }
65 
67 {
68  if (MyOwnLevelDensity) { delete theLevelDensityPtr; }
69 }
70 
72 {
73  G4int anA = fragment->GetA_asInt();
74  G4int aZ = fragment->GetZ_asInt();
75  ResidualA = anA - A;
76  ResidualZ = aZ - Z;
77  /*
78  G4cout << "G4GEMChannel: Z= " << Z << " A= " << A
79  << " FragmentZ= " << aZ << " FragmentA= " << anA
80  << " Zres= " << ResidualZ << " Ares= " << ResidualA
81  << G4endl;
82  */
83  // We only take into account channels which are physically allowed
84  EmissionProbability = 0.0;
85 
86  // Only channels which are physically allowed are taken into account
87  if (ResidualA >= ResidualZ && ResidualZ > 0 && ResidualA >= A) {
88 
89  //Effective excitation energy
90  G4double ExEnergy = fragment->GetExcitationEnergy()
92  if(ExEnergy > 0.0) {
94  G4double FragmentMass = fragment->GetGroundStateMass();
95  G4double Etot = FragmentMass + ExEnergy;
96  // Coulomb Barrier calculation
99  /*
100  G4cout << "Eexc(MeV)= " << ExEnergy/MeV
101  << " CoulBarrier(MeV)= " << CoulombBarrier/MeV << G4endl;
102  */
104 
105  // Maximal Kinetic Energy
107  + EvaporatedMass*EvaporatedMass)/(2.0*Etot)
109 
110  //G4cout << "CBarrier(MeV)= " << CoulombBarrier/MeV << G4endl;
111 
112  if (MaximalKineticEnergy > 0.0) {
113  // Total emission probability for this channel
116  }
117  }
118  }
119  }
120  //G4cout << "Prob= " << EmissionProbability << G4endl;
121  return EmissionProbability;
122 }
123 
125 {
126  G4Fragment* evFragment = 0;
127  G4double evEnergy = SampleKineticEnergy(*theNucleus) + EvaporatedMass;
128 
130  (std::sqrt((evEnergy - EvaporatedMass)*(evEnergy + EvaporatedMass))));
131 
132  G4LorentzVector EvaporatedMomentum(momentum, evEnergy);
133  G4LorentzVector ResidualMomentum = theNucleus->GetMomentum();
134  EvaporatedMomentum.boost(ResidualMomentum.boostVector());
135 
136  evFragment = new G4Fragment(A, Z, EvaporatedMomentum);
137  ResidualMomentum -= EvaporatedMomentum;
138  theNucleus->SetZandA_asInt(ResidualZ, ResidualA);
139  theNucleus->SetMomentum(ResidualMomentum);
140 
141  return evFragment;
142 }
143 
145 {
146  G4FragmentVector * theResult = new G4FragmentVector();
147  G4Fragment* frag0 = new G4Fragment(theNucleus);
148  G4Fragment* frag1 = EmittedFragment(frag0);
149  if(frag1) { theResult->push_back(frag1); }
150  theResult->push_back(frag0);
151  return theResult;
152 }
153 
155 // Samples fragment kinetic energy.
156 {
157  G4double U = fragment.GetExcitationEnergy();
158 
161 
162  // ***RESIDUAL***
163  //JMQ (September 2009) the following quantities refer to the RESIDUAL:
165  G4double Ux = (2.5 + 150.0/ResidualA)*MeV;
166  G4double Ex = Ux + delta0;
167  G4double InitialLevelDensity;
168  // ***end RESIDUAL ***
169 
170  // ***PARENT***
171  //JMQ (September 2009) the following quantities refer to the PARENT:
172 
174  fragment.GetZ_asInt());
176  fragment.GetZ_asInt(),
177  U-deltaCN);
178  G4double UxCN = (2.5 + 150.0/G4double(fragment.GetA_asInt()))*MeV;
179  G4double ExCN = UxCN + deltaCN;
180  G4double TCN = 1.0/(std::sqrt(aCN/UxCN) - 1.5/UxCN);
181  // ***end PARENT***
182 
183  //JMQ quantities calculated for CN in InitialLevelDensity
184  if ( U < ExCN )
185  {
186  G4double E0CN = ExCN - TCN*(G4Log(TCN/MeV) - G4Log(aCN*MeV)/4.0
187  - 1.25*G4Log(UxCN/MeV) + 2.0*std::sqrt(aCN*UxCN));
188  InitialLevelDensity = (pi/12.0)*G4Exp((U-E0CN)/TCN)/TCN;
189  }
190  else
191  {
192  G4double x = U-deltaCN;
193  G4double x1 = std::sqrt(aCN*x);
194  InitialLevelDensity = (pi/12.0)*G4Exp(2*x1)/(x*std::sqrt(x1));
195  }
196 
198  //JMQ BIG BUG fixed: hbarc instead of hbar_Planck !!!!
199  // it was fixed in total probability (for this channel) but remained still here!!
200  // G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbar_Planck*hbar_Planck);
201  G4double gg = (2.0*Spin+1.0)*EvaporatedMass/(pi2* hbarc*hbarc);
202  //
203  //JMQ fix on Rb and geometrical cross sections according to Furihata's paper
204  // (JAERI-Data/Code 2001-105, p6)
205  G4double Rb = 0.0;
206  if (A > 4)
207  {
208  G4double Ad = fG4pow->Z13(ResidualA);
209  G4double Aj = fG4pow->Z13(A);
210  Rb = (1.12*(Aj + Ad) - 0.86*((Aj+Ad)/(Aj*Ad))+2.85)*fermi;
211  }
212  else if (A>1)
213  {
214  G4double Ad = fG4pow->Z13(ResidualA);
215  G4double Aj = fG4pow->Z13(A);
216  Rb=1.5*(Aj+Ad)*fermi;
217  }
218  else
219  {
220  G4double Ad = fG4pow->Z13(ResidualA);
221  Rb = 1.5*Ad*fermi;
222  }
223  G4double GeometricalXS = pi*Rb*Rb;
224 
225  G4double ConstantFactor = gg*GeometricalXS*Alpha*pi/(InitialLevelDensity*12);
226  //JMQ : this is the assimptotic maximal kinetic energy of the ejectile
227  // (obtained by energy-momentum conservation).
228  // In general smaller than U-theSeparationEnergy
230  G4double KineticEnergy;
231  G4double Probability;
232 
233  for(G4int i=0; i<100; ++i) {
234  KineticEnergy = CoulombBarrier + G4UniformRand()*(MaximalKineticEnergy);
235  G4double edelta = theEnergy-KineticEnergy-delta0;
236  Probability = ConstantFactor*(KineticEnergy + Beta);
237  G4double a =
239  G4double T = 1.0/(std::sqrt(a/Ux) - 1.5/Ux);
240  //JMQ fix in units
241 
242  if (theEnergy - KineticEnergy < Ex) {
243  G4double E0 = Ex - T*(G4Log(T) - G4Log(a)*0.25
244  - 1.25*G4Log(Ux) + 2.0*std::sqrt(a*Ux));
245  Probability *= G4Exp((theEnergy-KineticEnergy-E0)/T)/T;
246  } else {
247  G4double e2 = edelta*edelta;
248  Probability *=
249  G4Exp(2*std::sqrt(a*edelta) - 0.25*G4Log(a*edelta*e2*e2));
250  }
251  if(EmissionProbability*G4UniformRand() <= Probability) { break; }
252  }
253 
254  return KineticEnergy;
255 }
256 
258  // Samples a isotropic random vectorwith a magnitude given by Magnitude.
259  // By default Magnitude = 1.0
260 {
261  G4double CosTheta = 1.0 - 2.0*G4UniformRand();
262  G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
263  G4double Phi = twopi*G4UniformRand();
264  G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
265  Magnitude*std::sin(Phi)*SinTheta,
266  Magnitude*CosTheta);
267  return Vector;
268 }
269 
270 void G4GEMChannel::Dump() const
271 {
273 }
274 
275 
276 
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4ThreeVector IsotropicVector(G4double Magnitude=1.0)
G4Pow * fG4pow
static const double MeV
Definition: G4SIunits.hh:211
static G4double GetNuclearMass(const G4double A, const G4double Z)
virtual ~G4GEMChannel()
Definition: G4GEMChannel.cc:66
CLHEP::Hep3Vector G4ThreeVector
G4double CalcAlphaParam(const G4Fragment &) const
static const G4double e2
G4double CoulombBarrier
void Dump() const
G4double a
Definition: TRTMaterials.hh:39
virtual G4Fragment * EmittedFragment(G4Fragment *theNucleus)
G4bool MyOwnLevelDensity
static const double pi2
Definition: G4SIunits.hh:77
int G4int
Definition: G4Types.hh:78
G4double SampleKineticEnergy(const G4Fragment &fragment)
G4VCoulombBarrier * theCoulombBarrierPtr
G4double ResidualMass
Definition: G4GEMChannel.hh:99
#define G4UniformRand()
Definition: Randomize.hh:97
G4double Z13(G4int Z) const
Definition: G4Pow.hh:127
double A(double temperature)
G4int GetA_asInt() const
Definition: G4Fragment.hh:251
G4double GetSpin(void) const
G4PairingCorrection * pairingCorrection
G4double EvaporatedMass
Definition: G4GEMChannel.hh:98
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:284
void SetMomentum(const G4LorentzVector &value)
Definition: G4Fragment.hh:289
static const double twopi
Definition: G4SIunits.hh:75
G4double GetPairingCorrection(G4int A, G4int Z) const
std::vector< G4Fragment * > G4FragmentVector
Definition: G4Fragment.hh:63
static const double GeV
Definition: G4SIunits.hh:214
G4double GetGroundStateMass() const
Definition: G4Fragment.hh:273
G4double EmissionProbability
G4double CalcBetaParam(const G4Fragment &) const
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:71
static const double pi
Definition: G4SIunits.hh:74
virtual void Dump() const
void SetZandA_asInt(G4int Znew, G4int Anew)
Definition: G4Fragment.hh:261
const G4double x[NPOINTSGL]
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:256
virtual G4double GetCoulombBarrier(G4int ARes, G4int ZRes, G4double U) const =0
G4GEMProbability * theEvaporationProbabilityPtr
double G4double
Definition: G4Types.hh:76
G4double MaximalKineticEnergy
static const double fermi
Definition: G4SIunits.hh:102
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:268
CLHEP::HepLorentzVector G4LorentzVector