Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4GEMChannelVI.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: G4GEMChannelVI.cc 98577 2016-07-25 13:05:12Z vnivanch $
27 //
28 // GEM de-excitation model
29 // by V. Ivanchenko (July 2016)
30 //
31 
32 #include "G4GEMChannelVI.hh"
33 #include "G4VCoulombBarrier.hh"
34 #include "G4CoulombBarrier.hh"
35 #include "G4PairingCorrection.hh"
36 #include "G4NuclearLevelData.hh"
37 #include "G4LevelManager.hh"
38 #include "G4RandomDirection.hh"
39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "Randomize.hh"
42 #include "G4Pow.hh"
43 #include "G4Log.hh"
44 
45 // 10-Points Gauss-Legendre abcisas and weights
46 const G4double G4GEMChannelVI::ws[] = {
47  0.0666713443086881,
48  0.149451349150581,
49  0.219086362515982,
50  0.269266719309996,
51  0.295524224714753,
52  0.295524224714753,
53  0.269266719309996,
54  0.219086362515982,
55  0.149451349150581,
56  0.0666713443086881
57  };
58 const G4double G4GEMChannelVI::xs[] = {
59  -0.973906528517172,
60  -0.865063366688985,
61  -0.679409568299024,
62  -0.433395394129247,
63  -0.148874338981631,
64  0.148874338981631,
65  0.433395394129247,
66  0.679409568299024,
67  0.865063366688985,
68  0.973906528517172
69 };
70 
72  : A(theA), Z(theZ), levelDensity(0.1)
73 {
74  fG4pow = G4Pow::GetInstance();
75  Z13 = fG4pow->Z13(Z);
76  A13 = fG4pow->Z13(A);
77 
78  cBarrier = new G4CoulombBarrier(A, Z);
79  pairingCorrection = G4PairingCorrection::GetInstance();
80 
82  levelManager = nData->GetLevelManager(Z, A);
83  maxLevelE = levelManager->MaxLevelEnergy();
84 
85  massGround = G4NucleiProperties::GetNuclearMass(A, Z);
86 
87  resA = resZ = fragZ = fragA = nWarn = 0;
88  massGround = maxLevelE = Z13 = A13 = massFrag = eCBarrier
89  = resMassGround = resZ13 = resA13 = delta0
90  = delta1 = maxExc = maxProb = alphaP = betaP = maxKinEnergy = 0.0;
92 }
93 
95 {
96  delete cBarrier;
97 }
98 
100 {
101  G4DeexPrecoParameters* param =
103  levelDensity = param->GetLevelDensity();
104 }
105 
107 {
108  fragZ = fragment->GetZ_asInt();
109  fragA = fragment->GetA_asInt();
110  resZ = fragZ - Z;
111  resA = fragA - A;
112  G4double prob = 0.0;
113  if(resA < A || resA < resZ || resZ < 0 || (resA == A && resZ < Z)) {
114  return prob;
115  }
116 
117  resMassGround = G4NucleiProperties::GetNuclearMass(resA, resZ);
118  G4double exc = fragment->GetExcitationEnergy();
119  massFrag = fragment->GetGroundStateMass() + exc;
120  delta0 = pairingCorrection->GetPairingCorrection(fragA, fragZ);
121  eCBarrier = cBarrier->GetCoulombBarrier(resA, resZ, exc);
122 
123  maxExc = massFrag - massGround - resMassGround - eCBarrier - delta0;
124  if(maxExc < 0.0) { return prob; }
125 
126  resZ13 = fG4pow->Z13(resZ);
127  resA13 = fG4pow->Z13(resA);
128  delta1 = pairingCorrection->GetPairingCorrection(resA, resZ);
129 
130  G4double C = 0.0;
131  if(resA >= 50) {
132  C = -0.10/G4double(A);
133  } else if(resZ > 20) {
134  C = (0.123482-0.00534691*Z-0.0000610624*(Z*Z)+5.93719*1e-7*(Z*Z*Z)+
135  1.95687*1e-8*(Z*Z*Z*Z))/G4double(A);
136  }
137  if(0 == Z) {
138  alphaP = 0.76+1.93/resA13;
139  betaP = (1.66/(resA13*resA13)-0.05)*CLHEP::MeV/alphaP;
140  } else {
141  alphaP = 1.0 + C;
142  betaP = - eCBarrier;
143  }
144 
145  maxProb = 0.0;
146  G4double e0 = maxExc*0.5;
147 
148  // e is an excitation of emitted fragment
149  for (G4int i=0; i<NPOINTSGEM; ++i) {
150  prob += ws[i]*IntegratedProbability(e0*(xs[i] + 1.0));
151  }
152  prob *= coeff*e0*e0;
153 
154  /*
155  G4cout << "G4GEMChannelVI: Z= " << Z << " A= " << A
156  << " FragmentZ= " << aZ << " FragmentA= " << anA
157  << " Zres= " << ResidualZ << " Ares= " << ResidualA
158  << G4endl;
159  */
160 
161  //G4cout << "Prob= " << prob << G4endl;
162  return prob;
163 }
164 
165 G4double G4GEMChannelVI::IntegratedProbability(G4double exc)
166 {
167  G4double e0 = (maxExc - exc)*0.5;
168 
169  G4double y;
170  G4double sum = 0.0;
171 
172  for (G4int i=0; i<NPOINTSGEM; ++i) {
173  y = ProbabilityDistributionFunction(exc, e0*(xs[i] + 1.0));
174  maxProb = std::max(maxProb, y);
175  sum += ws[i]*y;
176  }
177  return sum;
178 }
179 
181 {
182  G4double exc;
183  G4double resExc;
184 
185  static G4double factor = 1.2;
186  maxProb *= factor;
187 
188  CLHEP::HepRandomEngine* rndm = G4Random::getTheEngine();
189  for(G4int i=0; i<100; ++i) {
190  do {
191  exc = maxExc*rndm->flat();
192  resExc = maxExc*rndm->flat();
193  } while (exc + resExc > maxExc);
194 
195  G4double prob = ProbabilityDistributionFunction(exc, resExc);
196  if(prob > maxProb && nWarn < 10) {
197  ++nWarn;
198  G4cout << "### G4GEMChannelVI::EmittedFragment WARNING: majoranta "
199  << maxProb << " is exceeded " << prob << "\n"
200  << " fragZ= " << fragZ << " fragA= " << fragA
201  << " Z= " << Z << " A= " << A
202  << " resZ= " << resZ << " resA= " << resA << "\n"
203  << " exc(MeV)= " << exc << " resExc(MeV)= " << resExc
204  << " maxExc(MeV)= " << maxExc << G4endl;
205  }
206  if(maxProb*rndm->flat() <= prob) { break; }
207  }
208  if(exc <= maxLevelE) {
209  exc = FindLevel(levelManager, exc, maxExc - resExc);
210  }
211  if(resA >= nData->GetMinA(resZ) && resA <= nData->GetMaxA(resZ)
212  && resExc < nData->GetMaxLevelEnergy(Z, A)) {
213  const G4LevelManager* lman = nData->GetLevelManager(Z, A);
214  if(lman) { resExc = FindLevel(lman, resExc, maxExc - exc); }
215  }
216 
217  G4LorentzVector lv0 = theNucleus->GetMomentum();
218  G4double mass1 = massGround + exc;
219  G4double mass2 = resMassGround + resExc;
220 
221  G4double e1 = 0.5*((massFrag - mass2)*(massFrag + mass2)
222  + mass1*mass1)/massFrag;
223 
224  G4double p1(0.0);
225  if(e1 > mass1) {
226  p1 = std::sqrt((e1 - mass1)*(e1 + mass1));
227  } else {
228  e1 = mass1;
229  }
231  G4LorentzVector lv1(p1*v, e1);
232 
233  G4ThreeVector boostVector = lv0.boostVector();
234  lv1.boost(boostVector);
235 
236  G4Fragment* frag = new G4Fragment(A, Z, lv1);
237 
238  G4double e2 = massFrag - e1;
239  if(e2 < mass2) {
240  e2 = mass2;
241  p1 = 0.0;
242  }
243  lv0.set(-v*p1, e2);
244  lv0.boost(boostVector);
245 
246  theNucleus->SetZandA_asInt(resZ, resA);
247  theNucleus->SetMomentum(lv0);
248 
249  return frag;
250 }
251 
252 G4double
253 G4GEMChannelVI::ProbabilityDistributionFunction(G4double exc, G4double resExc)
254 {
255  G4double Ux = (2.5 + 150.0/G4double(resA))*CLHEP::MeV;
256  G4double Ex = Ux + delta1;
257  G4double T = 1.0/(std::sqrt(levelDensity/Ux) - 1.5/Ux);
258  G4double E0 = Ex - T*(G4Log(T) - G4Log(levelDensity)*0.25
259  - 1.25*G4Log(Ux) + 2.0*std::sqrt(levelDensity*Ux));
260 
261  G4double UxCN = (2.5 + 150.0/G4double(A))*CLHEP::MeV;
262  G4double ExCN = UxCN + delta0;
263  G4double TCN = 1.0/(std::sqrt(levelDensity/UxCN) - 1.5/UxCN);
264 
265  G4double mass1 = massGround + exc;
266  G4double mass2 = resMassGround + resExc;
267 
268  maxKinEnergy = 0.5*((massFrag - mass2)*(massFrag + mass2)
269  + mass1*mass1)/massFrag - mass1;
270  maxKinEnergy = std::max(maxKinEnergy, 0.0);
271 
272  G4double Width = 0.0;
273  G4double t = maxKinEnergy/T;
274  if ( maxKinEnergy < Ex ) {
275  Width = (I1(t,t)*T + (betaP+eCBarrier)*I0(t))/G4Exp(E0/T);
276 
277  } else {
278 
279  G4double tx = Ex/T;
280  G4double s0 = 2.0*std::sqrt(levelDensity*(maxKinEnergy-delta0));
281  G4double sx = 2.0*std::sqrt(levelDensity*(Ex-delta0));
282 
283  // VI: protection against FPE exception
284  if(s0 > 350.) { s0 = 350.; }
285 
286  G4double expE0T = G4Exp(E0/T);
287  G4double exps0 = G4Exp(s0);
288  static const G4double sqrt2 = std::sqrt(2.0);
289 
290  Width = I1(t,tx)*T/expE0T + I3(s0,sx)*exps0/(sqrt2*levelDensity);
291 
292  if (0 == Z) {
293  Width += (betaP+eCBarrier)*(I0(tx)/expE0T + 2.0*sqrt2*I2(s0,sx)*exps0);
294  }
295  }
296  Width *= alphaP*massFrag;
297 
298 
299  //JMQ 190709 fix on Rb and geometrical cross sections according to
300  // Furihata's paper (JAERI-Data/Code 2001-105, p6)
301  G4double Rb = 0.0;
302  if (A > 4) {
303  Rb = 1.12*(resA13 + A13) - 0.86*((resA13 + A13)/(resA13*A13))+2.85;
304  } else if (A > 1) {
305  Rb=1.5*(resA13 + A13);
306  } else {
307  Rb = 1.5*resA13;
308  }
309 
310  G4double ild;
311  if (exc < ExCN ) {
312  G4double E0CN = ExCN - TCN*(G4Log(TCN) - 0.25*G4Log(levelDensity)
313  - 1.25*G4Log(UxCN)
314  + 2.0*std::sqrt(levelDensity*UxCN));
315  ild = G4Exp((exc-E0CN)/TCN)/TCN;
316  } else {
317  G4double x = exc - delta0;
318  G4double x1 = std::sqrt(levelDensity*x);
319  ild = G4Exp(2*x1)/(x*std::sqrt(x1));
320  }
321 
322  Width *= (Rb*Rb/ild);
323  return Width;
324 }
325 
326 G4double G4GEMChannelVI::FindLevel(const G4LevelManager* man,
327  G4double exc, G4double exclim)
328 {
329  size_t idx = man->NearestLowEdgeLevelIndex(exc);
330  size_t idxm = man->NumberOfTransitions();
331  G4double e1 = (G4double)man->LevelEnergy(idx);
332  if(idx + 1 < idxm) {
333  G4double e2 = (G4double)man->LevelEnergy(idx+1);
334  if(e2 <= exclim) {
335  G4int s1 = std::abs(man->SpinParity(idx))+1;
336  G4int s2 = std::abs(man->SpinParity(idx+1))+1;
337  G4double pr = (G4double)s1/(G4double)(s1 + s2);
338  pr = (exc - e1 <= e2 - exc) ? 1.0 - (1.0 - pr)*2*(exc - e1)/(e2 - e1) :
339  2*pr*(e2 - exc)/(e2 - e1);
340  exc = (G4UniformRand() < pr) ? e1 : e2;
341  } else {
342  exc = e1;
343  }
344  } else { exc = e1; }
345  return exc;
346 }
347 
348 G4double G4GEMChannelVI::I3(G4double s0, G4double sx)
349 {
350  G4double s2 = s0*s0;
351  G4double sx2 = sx*sx;
352  G4double S = 1.0/std::sqrt(s0);
353  G4double S2 = S*S;
354  G4double Sx = 1.0/std::sqrt(sx);
355  G4double Sx2 = Sx*Sx;
356 
357  G4double p1 = S *(2.0 + S2 *( 4.0 + S2 *( 13.5 + S2 *( 60.0 + S2 * 325.125 ))));
358  G4double p2 = Sx*Sx2 *((s2-sx2)
359  + Sx2 *((1.5*s2+0.5*sx2)
360  + Sx2 *((3.75*s2+0.25*sx2)
361  + Sx2 *((12.875*s2+0.625*sx2)
362  + Sx2 *((59.0625*s2+0.9375*sx2)
363  + Sx2 *(324.8*s2+3.28*sx2))))));
364  p2 *= G4Exp(sx-s0);
365  return p1-p2;
366 }
367 
369 {
370 }
371 
372 
373 
size_t NearestLowEdgeLevelIndex(G4double energy) const
G4int GetMinA(G4int Z) const
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
Hep3Vector boostVector() const
static G4double GetNuclearMass(const G4double A, const G4double Z)
double S(double temp)
G4int SpinParity(size_t i) const
G4GEMChannelVI(G4int theA, G4int theZ)
const G4LevelManager * GetLevelManager(G4int Z, G4int A)
static constexpr double hbarc
G4ThreeVector G4RandomDirection()
tuple x
Definition: test.py:50
virtual double flat()=0
double C(double temp)
virtual G4Fragment * EmittedFragment(G4Fragment *theNucleus) final
virtual void Initialise() final
G4double GetLevelDensity() const
int G4int
Definition: G4Types.hh:78
G4float LevelEnergy(size_t i) const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
G4double Z13(G4int Z) const
Definition: G4Pow.hh:127
double A(double temperature)
G4int GetA_asInt() const
Definition: G4Fragment.hh:266
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:307
HepLorentzVector & boost(double, double, double)
void SetMomentum(const G4LorentzVector &value)
Definition: G4Fragment.hh:312
static constexpr double MeV
G4double GetPairingCorrection(G4int A, G4int Z) const
G4DeexPrecoParameters * GetParameters()
size_t NumberOfTransitions() const
G4double GetGroundStateMass() const
Definition: G4Fragment.hh:288
G4int GetMaxA(G4int Z) const
tuple v
Definition: test.py:18
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static G4PairingCorrection * GetInstance()
void set(double x, double y, double z, double t)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
void SetZandA_asInt(G4int Znew, G4int Anew)
Definition: G4Fragment.hh:276
virtual ~G4GEMChannelVI()
static constexpr double fermi
Definition: SystemOfUnits.h:83
G4int GetZ_asInt() const
Definition: G4Fragment.hh:271
G4double GetMaxLevelEnergy(G4int Z, G4int A) const
virtual void Dump() const
virtual G4double GetCoulombBarrier(G4int ARes, G4int ZRes, G4double U) const =0
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4float MaxLevelEnergy() const
static G4NuclearLevelData * GetInstance()
static constexpr double pi
Definition: SystemOfUnits.h:54
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:283
const G4int NPOINTSGEM
virtual G4double GetEmissionProbability(G4Fragment *theNucleus) final