Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4PreCompoundEmission.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: G4PreCompoundEmission.cc 108685 2018-02-27 07:58:38Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4PreCompoundEmission
34 //
35 // Author: V.Lara
36 //
37 // Modified:
38 // 15.01.2010 J.M.Quesada added protection against unphysical values of parameter an
39 // 19.01.2010 V.Ivanchenko simplified computation of parameter an, sample cosTheta
40 // instead of theta; protect all calls to sqrt
41 // 20.08.2010 V.Ivanchenko added G4Pow and G4PreCompoundParameters pointers
42 // use int Z and A and cleanup
43 //
44 
45 #include "G4PreCompoundEmission.hh"
46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "G4Pow.hh"
49 #include "G4Exp.hh"
50 #include "G4Log.hh"
51 #include "Randomize.hh"
52 #include "G4RandomDirection.hh"
54 #include "G4HETCEmissionFactory.hh"
55 #include "G4HadronicException.hh"
56 #include "G4NuclearLevelData.hh"
57 #include "G4DeexPrecoParameters.hh"
58 
60 {
65  G4DeexPrecoParameters* param =
67  fLevelDensity = param->GetLevelDensity();
68  fFermiEnergy = param->GetFermiEnergy();
70 }
71 
73 {
74  delete theFragmentsFactory;
75  delete theFragmentsVector;
76 }
77 
79 {
82  if (theFragmentsVector) {
84  } else {
87  }
88 }
89 
91 {
94  if (theFragmentsVector) {
96  } else {
99  }
100 }
101 
104 {
105  // Choose a Fragment for emission
106  G4VPreCompoundFragment * thePreFragment =
108  if (thePreFragment == nullptr)
109  {
110  G4cout << "G4PreCompoundEmission::PerformEmission : "
111  << "I couldn't choose a fragment\n"
112  << "while trying to de-excite\n"
113  << aFragment << G4endl;
114  throw G4HadronicException(__FILE__, __LINE__, "");
115  }
116 
117  //G4cout << "Chosen fragment: " << G4endl;
118  //G4cout << *thePreFragment << G4endl;
119 
120  // Kinetic Energy of emitted fragment
121  G4double kinEnergy = thePreFragment->SampleKineticEnergy(aFragment);
122  kinEnergy = std::max(kinEnergy, 0.0);
123 
124  // Calculate the fragment momentum (three vector)
126  AngularDistribution(thePreFragment,aFragment,kinEnergy);
127  } else {
128  G4double pmag =
129  std::sqrt(kinEnergy*(kinEnergy + 2.0*thePreFragment->GetNuclearMass()));
131  }
132 
133  // Mass of emittef fragment
134  G4double EmittedMass = thePreFragment->GetNuclearMass();
135  // Now we can calculate the four momentum
136  // both options are valid and give the same result but 2nd one is faster
137  G4LorentzVector Emitted4Momentum(theFinalMomentum,EmittedMass + kinEnergy);
138 
139  // Perform Lorentz boost
140  G4LorentzVector Rest4Momentum = aFragment.GetMomentum();
141  Emitted4Momentum.boost(Rest4Momentum.boostVector());
142 
143  // Set emitted fragment momentum
144  thePreFragment->SetMomentum(Emitted4Momentum);
145 
146  // NOW THE RESIDUAL NUCLEUS
147  // ------------------------
148 
149  Rest4Momentum -= Emitted4Momentum;
150 
151  // Update nucleus parameters:
152  // --------------------------
153 
154  // Z and A
155  aFragment.SetZandA_asInt(thePreFragment->GetRestZ(),
156  thePreFragment->GetRestA());
157 
158  // Number of excitons
159  aFragment.SetNumberOfParticles(aFragment.GetNumberOfParticles()-
160  thePreFragment->GetA());
161  // Number of charges
162  aFragment.SetNumberOfCharged(aFragment.GetNumberOfCharged()-
163  thePreFragment->GetZ());
164 
165  // Update nucleus momentum
166  // A check on consistence of Z, A, and mass will be performed
167  aFragment.SetMomentum(Rest4Momentum);
168 
169  // Create a G4ReactionProduct
170  G4ReactionProduct * MyRP = thePreFragment->GetReactionProduct();
171 
172  // if(kinEnergy < MeV) {
173  //G4cout << "G4PreCompoundEmission::Fragment emitted" << G4endl;
174  //G4cout << *thePreFragment << G4endl;
175  // }
176  return MyRP;
177 }
178 
180  G4VPreCompoundFragment* thePreFragment,
181  const G4Fragment& aFragment,
182  G4double ekin)
183 {
184  G4int p = aFragment.GetNumberOfParticles();
185  G4int h = aFragment.GetNumberOfHoles();
186  G4double U = aFragment.GetExcitationEnergy();
187 
188  // Emission particle separation energy
189  G4double Bemission = thePreFragment->GetBindingEnergy();
190 
191  G4double gg = (6.0/pi2)*aFragment.GetA_asInt()*fLevelDensity;
192 
193  // Average exciton energy relative to bottom of nuclear well
194  G4double Eav = 2*p*(p+1)/((p+h)*gg);
195 
196  // Excitation energy relative to the Fermi Level
197  G4double Uf = std::max(U - (p - h)*fFermiEnergy , 0.0);
198  // G4double Uf = U - KineticEnergyOfEmittedFragment - Bemission;
199 
200  G4double w_num = rho(p+1, h, gg, Uf, fFermiEnergy);
201  G4double w_den = rho(p, h, gg, Uf, fFermiEnergy);
202  if (w_num > 0.0 && w_den > 0.0)
203  {
204  Eav *= (w_num/w_den);
205  Eav += - Uf/(p+h) + fFermiEnergy;
206  }
207  else
208  {
209  Eav = fFermiEnergy;
210  }
211 
212  // VI + JMQ 19/01/2010 update computation of the parameter an
213  //
214  G4double an = 0.0;
215  G4double Eeff = ekin + Bemission + fFermiEnergy;
216  if(ekin > DBL_MIN && Eeff > DBL_MIN) {
217 
218  G4double zeta = std::max(1.0,9.3/std::sqrt(ekin/MeV));
219 
220  // This should be the projectile energy. If I would know which is
221  // the projectile (proton, neutron) I could remove the binding energy.
222  // But, what happens if INC precedes precompound? This approximation
223  // seems to work well enough
224  G4double ProjEnergy = aFragment.GetExcitationEnergy();
225 
226  an = 3*std::sqrt((ProjEnergy+fFermiEnergy)*Eeff)/(zeta*Eav);
227 
228  G4int ne = aFragment.GetNumberOfExcitons() - 1;
229  if ( ne > 1 ) { an /= (G4double)ne; }
230 
231  // protection of exponent
232  if ( an > 10. ) { an = 10.; }
233  }
234 
235  // sample cosine of theta and not theta as in old versions
236  G4double random = G4UniformRand();
237  G4double cost;
238 
239  if(an < 0.1) { cost = 1. - 2*random; }
240  else {
241  G4double exp2an = G4Exp(-2*an);
242  cost = 1. + G4Log(1-random*(1-exp2an))/an;
243  if(cost > 1.) { cost = 1.; }
244  else if(cost < -1.) {cost = -1.; }
245  }
246 
248 
249  // Calculate the momentum magnitude of emitted fragment
250  G4double pmag =
251  std::sqrt(ekin*(ekin + 2.0*thePreFragment->GetNuclearMass()));
252 
253  G4double sint = std::sqrt((1.0-cost)*(1.0+cost));
254 
255  theFinalMomentum.set(pmag*std::cos(phi)*sint,pmag*std::sin(phi)*sint,
256  pmag*cost);
257 
258  // theta is the angle wrt the incident direction
259  G4ThreeVector theIncidentDirection = aFragment.GetMomentum().vect().unit();
260  theFinalMomentum.rotateUz(theIncidentDirection);
261 }
262 
264  G4double E, G4double Ef) const
265 {
266  // 25.02.2010 V.Ivanchenko added more protections
267  G4double Aph = (p*p + h*h + p - 3.0*h)/(4.0*gg);
268  // G4double alpha = (p*p + h*h)/(2.0*gg);
269 
270  if ( E - Aph < 0.0) { return 0.0; }
271 
272  G4double logConst = (p+h)*G4Log(gg)
273  - g4calc->logfactorial(p+h-1) - g4calc->logfactorial(p)
274  - g4calc->logfactorial(h);
275 
276  // initialise values using j=0
277 
278  G4double t1=1;
279  G4double t2=1;
280  G4double logt3 = (p+h-1) * G4Log(E-Aph) + logConst;
281  const G4double logmax = 200.;
282  if(logt3 > logmax) { logt3 = logmax; }
283  G4double tot = G4Exp( logt3 );
284 
285  // and now sum rest of terms
286  // 25.02.2010 V.Ivanchenko change while to for loop and cleanup
287  G4double Eeff = E - Aph;
288  for(G4int j=1; j<=h; ++j)
289  {
290  Eeff -= Ef;
291  if(Eeff < 0.0) { break; }
292  t1 *= -1.;
293  t2 *= (G4double)(h+1-j)/(G4double)j;
294  logt3 = (p+h-1) * G4Log( Eeff) + logConst;
295  if(logt3 > logmax) { logt3 = logmax; }
296  tot += t1*t2*G4Exp(logt3);
297  }
298 
299  return tot;
300 }
G4VPreCompoundEmissionFactory * theFragmentsFactory
void set(double x, double y, double z)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double GetFermiEnergy() const
G4ReactionProduct * GetReactionProduct() const
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:300
G4ReactionProduct * PerformEmission(G4Fragment &aFragment)
TTree * t1
Definition: plottest35.C:26
void SetMomentum(const G4LorentzVector &value)
static constexpr double MeV
Definition: G4SIunits.hh:214
G4DeexPrecoParameters * GetParameters()
#define G4endl
Definition: G4ios.hh:61
void SetZandA_asInt(G4int Znew, G4int Anew)
Definition: G4Fragment.hh:269
virtual G4double SampleKineticEnergy(const G4Fragment &aFragment)=0
const char * p
Definition: xmltok.h:285
static G4NuclearLevelData * GetInstance()
std::vector< G4VPreCompoundFragment * > * GetFragmentVector()
G4int GetNumberOfHoles() const
Definition: G4Fragment.hh:358
G4int GetA_asInt() const
Definition: G4Fragment.hh:259
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
void AngularDistribution(G4VPreCompoundFragment *theFragment, const G4Fragment &aFragment, G4double KineticEnergy)
G4int GetA() const
G4ThreeVector G4RandomDirection()
G4int GetNumberOfParticles() const
Definition: G4Fragment.hh:338
G4double G4Log(G4double x)
Definition: G4Log.hh:230
double G4double
Definition: G4Types.hh:76
static G4Pow * GetInstance()
Definition: G4Pow.cc:57
G4PreCompoundFragmentVector * theFragmentsVector
G4int GetZ() const
G4int GetNumberOfCharged() const
Definition: G4Fragment.hh:343
#define G4UniformRand()
Definition: Randomize.hh:53
G4int GetRestA() const
G4double GetBindingEnergy() const
void SetMomentum(const G4LorentzVector &value)
Definition: G4Fragment.hh:305
static constexpr double pi2
Definition: G4SIunits.hh:78
void SetNumberOfCharged(G4int value)
Definition: G4Fragment.hh:382
G4double logfactorial(G4int Z) const
Definition: G4Pow.hh:252
Hep3Vector unit() const
TTree * t2
Definition: plottest35.C:36
static constexpr double twopi
Definition: SystemOfUnits.h:55
G4double GetLevelDensity() const
G4double GetNuclearMass() const
int G4int
Definition: G4Types.hh:78
#define DBL_MIN
Definition: templates.hh:75
void SetNumberOfParticles(G4int value)
Definition: G4Fragment.hh:377
G4GLOB_DLL std::ostream G4cout
Hep3Vector boostVector() const
Hep3Vector vect() const
G4int GetRestZ() const
G4int GetNumberOfExcitons() const
Definition: G4Fragment.hh:333
G4double rho(G4int p, G4int h, G4double gg, G4double E, G4double Ef) const
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:276
HepLorentzVector & boost(double, double, double)
G4VPreCompoundFragment * ChooseFragment()