Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4GEMChannel.cc
이 파일의 문서화 페이지로 가기
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 107060 2017-11-01 15:00:04Z 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 "G4NuclearLevelData.hh"
42 #include "G4PairingCorrection.hh"
43 #include "G4PhysicalConstants.hh"
44 #include "G4SystemOfUnits.hh"
45 #include "G4Pow.hh"
46 #include "G4Log.hh"
47 #include "G4Exp.hh"
48 #include "G4RandomDirection.hh"
49 
50 G4GEMChannel::G4GEMChannel(G4int theA, G4int theZ, const G4String & aName,
51  G4GEMProbability * aEmissionStrategy) :
52  G4VEvaporationChannel(aName),
53  A(theA),
54  Z(theZ),
55  theEvaporationProbabilityPtr(aEmissionStrategy),
56  EmissionProbability(0.0),
57  MaximalKineticEnergy(-CLHEP::GeV)
58 {
59  theCoulombBarrierPtr = new G4GEMCoulombBarrier(theA, theZ);
62  MyOwnLevelDensity = true;
66  ResidualZ = ResidualA = 0;
69 }
70 
72 {
73  if (MyOwnLevelDensity) { delete theLevelDensityPtr; }
74  delete theCoulombBarrierPtr;
75 }
76 
78 {
79  G4int anA = fragment->GetA_asInt();
80  G4int aZ = fragment->GetZ_asInt();
81  ResidualA = anA - A;
82  ResidualZ = aZ - Z;
83  /*
84  G4cout << "G4GEMChannel: Z= " << Z << " A= " << A
85  << " FragmentZ= " << aZ << " FragmentA= " << anA
86  << " Zres= " << ResidualZ << " Ares= " << ResidualA
87  << G4endl;
88  */
89  // We only take into account channels which are physically allowed
90  EmissionProbability = 0.0;
91 
92  // Only channels which are physically allowed are taken into account
93  if (ResidualA >= ResidualZ && ResidualZ > 0 && ResidualA >= A) {
94 
95  //Effective excitation energy
96  G4double ExEnergy = fragment->GetExcitationEnergy()
98  if(ExEnergy > 0.0) {
100  G4double FragmentMass = fragment->GetGroundStateMass();
101  G4double Etot = FragmentMass + ExEnergy;
102  // Coulomb Barrier calculation
103  CoulombBarrier =
105  /*
106  G4cout << "Eexc(MeV)= " << ExEnergy/MeV
107  << " CoulBarrier(MeV)= " << CoulombBarrier/MeV << G4endl;
108  */
110 
111  // Maximal Kinetic Energy
113  + EvaporatedMass*EvaporatedMass)/(2.0*Etot)
115 
116  //G4cout << "CBarrier(MeV)= " << CoulombBarrier/MeV << G4endl;
117 
118  if (MaximalKineticEnergy > 0.0) {
119  // Total emission probability for this channel
122  }
123  }
124  }
125  }
126  //G4cout << "Prob= " << EmissionProbability << G4endl;
127  return EmissionProbability;
128 }
129 
131 {
132  G4Fragment* evFragment = 0;
133  G4double evEnergy = SampleKineticEnergy(*theNucleus) + EvaporatedMass;
134 
135  G4ThreeVector momentum = G4RandomDirection()*
136  std::sqrt((evEnergy - EvaporatedMass)*(evEnergy + EvaporatedMass));
137 
138  G4LorentzVector EvaporatedMomentum(momentum, evEnergy);
139  G4LorentzVector ResidualMomentum = theNucleus->GetMomentum();
140  EvaporatedMomentum.boost(ResidualMomentum.boostVector());
141 
142  evFragment = new G4Fragment(A, Z, EvaporatedMomentum);
143  ResidualMomentum -= EvaporatedMomentum;
144  theNucleus->SetZandA_asInt(ResidualZ, ResidualA);
145  theNucleus->SetMomentum(ResidualMomentum);
146 
147  return evFragment;
148 }
149 
151 // Samples fragment kinetic energy.
152 {
153  G4double U = fragment.GetExcitationEnergy();
154 
157 
158  // ***RESIDUAL***
159  //JMQ (September 2009) the following quantities refer to the RESIDUAL:
161  G4double Ux = (2.5 + 150.0/ResidualA)*MeV;
162  G4double Ex = Ux + delta0;
163  G4double InitialLevelDensity;
164  // ***end RESIDUAL ***
165 
166  // ***PARENT***
167  //JMQ (September 2009) the following quantities refer to the PARENT:
168 
170  fragment.GetZ_asInt());
172  fragment.GetZ_asInt(),
173  U-deltaCN);
174  G4double UxCN = (2.5 + 150.0/G4double(fragment.GetA_asInt()))*MeV;
175  G4double ExCN = UxCN + deltaCN;
176  G4double TCN = 1.0/(std::sqrt(aCN/UxCN) - 1.5/UxCN);
177  // ***end PARENT***
178 
179  //JMQ quantities calculated for CN in InitialLevelDensity
180  if ( U < ExCN )
181  {
182  G4double E0CN = ExCN - TCN*(G4Log(TCN/MeV) - G4Log(aCN*MeV)/4.0
183  - 1.25*G4Log(UxCN/MeV) + 2.0*std::sqrt(aCN*UxCN));
184  InitialLevelDensity = (pi/12.0)*G4Exp((U-E0CN)/TCN)/TCN;
185  }
186  else
187  {
188  G4double x = U-deltaCN;
189  G4double x1 = std::sqrt(aCN*x);
190  InitialLevelDensity = (pi/12.0)*G4Exp(2*x1)/(x*std::sqrt(x1));
191  }
192 
194  //JMQ BIG BUG fixed: hbarc instead of hbar_Planck !!!!
195  // it was fixed in total probability (for this channel) but remained still here!!
196  // G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbar_Planck*hbar_Planck);
197  G4double gg = (2.0*Spin+1.0)*EvaporatedMass/(pi2* hbarc*hbarc);
198  //
199  //JMQ fix on Rb and geometrical cross sections according to Furihata's paper
200  // (JAERI-Data/Code 2001-105, p6)
201  G4double Rb = 0.0;
202  if (A > 4)
203  {
204  G4double Ad = fG4pow->Z13(ResidualA);
205  G4double Aj = fG4pow->Z13(A);
206  Rb = (1.12*(Aj + Ad) - 0.86*((Aj+Ad)/(Aj*Ad))+2.85)*fermi;
207  }
208  else if (A>1)
209  {
210  G4double Ad = fG4pow->Z13(ResidualA);
211  G4double Aj = fG4pow->Z13(A);
212  Rb=1.5*(Aj+Ad)*fermi;
213  }
214  else
215  {
216  G4double Ad = fG4pow->Z13(ResidualA);
217  Rb = 1.5*Ad*fermi;
218  }
219  G4double GeometricalXS = pi*Rb*Rb;
220 
221  G4double ConstantFactor = gg*GeometricalXS*Alpha*pi/(InitialLevelDensity*12);
222  //JMQ : this is the assimptotic maximal kinetic energy of the ejectile
223  // (obtained by energy-momentum conservation).
224  // In general smaller than U-theSeparationEnergy
226  G4double KineticEnergy;
227  G4double Probability;
228 
229  for(G4int i=0; i<100; ++i) {
230  KineticEnergy = CoulombBarrier + G4UniformRand()*(MaximalKineticEnergy);
231  G4double edelta = theEnergy-KineticEnergy-delta0;
232  Probability = ConstantFactor*(KineticEnergy + Beta);
233  G4double a =
235  G4double T = 1.0/(std::sqrt(a/Ux) - 1.5/Ux);
236  //JMQ fix in units
237 
238  if (theEnergy - KineticEnergy < Ex) {
239  G4double E0 = Ex - T*(G4Log(T) - G4Log(a)*0.25
240  - 1.25*G4Log(Ux) + 2.0*std::sqrt(a*Ux));
241  Probability *= G4Exp((theEnergy-KineticEnergy-E0)/T)/T;
242  } else {
243  G4double e2 = edelta*edelta;
244  Probability *=
245  G4Exp(2*std::sqrt(a*edelta) - 0.25*G4Log(a*edelta*e2*e2));
246  }
247  if(EmissionProbability*G4UniformRand() <= Probability) { break; }
248  }
249 
250  return KineticEnergy;
251 }
252 
253 void G4GEMChannel::Dump() const
254 {
256 }
257 
258 
259 
Float_t x
Definition: compare.C:6
G4double EvaporatedMass
Definition: G4GEMChannel.hh:90
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:300
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4double Z13(G4int Z) const
Definition: G4Pow.hh:132
G4double MaximalKineticEnergy
static constexpr double MeV
Definition: G4SIunits.hh:214
virtual G4double LevelDensityParameter(G4int A, G4int Z, G4double U) const =0
Float_t x1[n_points_granero]
Definition: compare.C:5
void SetZandA_asInt(G4int Znew, G4int Anew)
Definition: G4Fragment.hh:269
G4double ResidualMass
Definition: G4GEMChannel.hh:91
static G4NuclearLevelData * GetInstance()
G4PairingCorrection * pairingCorrection
static constexpr double hbarc
G4int GetA_asInt() const
Definition: G4Fragment.hh:259
G4double CoulombBarrier
G4GEMChannel(G4int theA, G4int theZ, const G4String &aName, G4GEMProbability *aEmissionStrategy)
Definition: G4GEMChannel.cc:50
G4ThreeVector G4RandomDirection()
G4double GetGroundStateMass() const
Definition: G4Fragment.hh:281
void Dump() const
G4double CalcBetaParam(const G4Fragment &) const
G4double SampleKineticEnergy(const G4Fragment &fragment)
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4GEMProbability * theEvaporationProbabilityPtr
Definition: G4GEMChannel.hh:96
Float_t Z
double G4double
Definition: G4Types.hh:76
static constexpr double fermi
Definition: G4SIunits.hh:103
G4double CalcAlphaParam(const G4Fragment &) const
static G4Pow * GetInstance()
Definition: G4Pow.cc:57
#define G4UniformRand()
Definition: Randomize.hh:53
virtual void Dump() const
virtual ~G4GEMChannel()
Definition: G4GEMChannel.cc:71
double A(double temperature)
void SetMomentum(const G4LorentzVector &value)
Definition: G4Fragment.hh:305
static G4double GetNuclearMass(const G4double A, const G4double Z)
static constexpr double pi2
Definition: G4SIunits.hh:78
G4VLevelDensityParameter * theLevelDensityPtr
G4int GetZ_asInt() const
Definition: G4Fragment.hh:264
G4double GetSpin(void) const
virtual G4double GetCoulombBarrier(G4int ARes, G4int ZRes, G4double U) const =0
int G4int
Definition: G4Types.hh:78
G4VCoulombBarrier * theCoulombBarrierPtr
G4double GetPairingCorrection(G4int A, G4int Z) const
void SetCoulomBarrier(const G4VCoulombBarrier *aCoulombBarrierStrategy)
Hep3Vector boostVector() const
virtual G4Fragment * EmittedFragment(G4Fragment *theNucleus)
static constexpr double pi
Definition: G4SIunits.hh:75
G4bool MyOwnLevelDensity
Definition: G4GEMChannel.hh:99
G4PairingCorrection * GetPairingCorrection()
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:276
virtual G4double GetEmissionProbability(G4Fragment *theNucleus)
Definition: G4GEMChannel.cc:77
static constexpr double GeV
Definition: G4SIunits.hh:217
G4double EmissionProbability
HepLorentzVector & boost(double, double, double)
G4Pow * fG4pow
Definition: G4GEMChannel.hh:93