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