Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4VEmissionProbability.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 //
27 // $Id: G4VEmissionProbability.cc 66241 2012-12-13 18:34:42Z gunter $
28 //
29 // Hadronic Process: Nuclear De-excitations
30 // by V. Lara (Oct 1998)
31 //
32 // Modifications:
33 // 28.10.2010 V.Ivanchenko defined members in constructor and cleaned up
34 
36 #include "G4NuclearLevelData.hh"
37 #include "G4DeexPrecoParameters.hh"
38 #include "Randomize.hh"
39 
41  :OPTxs(3),fVerbose(0),theZ(Z),theA(A),
42  LevelDensity(0.1),elimit(CLHEP::MeV),accuracy(0.02)
43 {
46  length = nfilled = 0;
48 }
49 
51 {}
52 
54 {
56  OPTxs = param->GetDeexModelType();
57  LevelDensity = param->GetLevelDensity();
58 }
59 
61 {
62  if(nbin > 0) {
63  fProb.clear();
64  fEner.clear();
65  fEner.resize(nbin+1, 0.0);
66  fProb.resize(nbin+1, 0.0);
67  length = nbin;
68  }
69  if(de > 0.0) { elimit = de; }
70  if(accuracy > 0.0) { accuracy = eps; }
71 }
72 
74 {
75  return 0.0;
76 }
77 
78 G4double
80 {
81  totProbability = 0.0;
82  if(elow >= ehigh) { return totProbability; }
83 
84  emin = elow;
85  emax = ehigh;
86  eCoulomb = cb;
87 
88  size_t nbin = (size_t)((emax - emin)/elimit) + 1;
89  G4double edelta = elimit;
90  if(nbin < 3) {
91  nbin = 3;
92  edelta = (emax - emin)/(G4double)nbin;
93  }
94  if(nbin > length) {
95  fEner.resize(nbin + 1);
96  fProb.resize(nbin + 1);
97  length = nbin;
98  }
99 
100  G4double x(emin), del, y;
101  eprobmax = x + edelta*0.1;
103  probmax = problast;
104  fEner[0] = emin;
105  fProb[0] = probmax;
106  size_t i(0);
107  if(fVerbose > 1) {
108  G4cout << "### G4VEmissionProbability::IntegrateProbability: "
109  << " Emax= " << emax << " QB= " << cb << " nbin= " << nbin
110  << G4endl;
111  G4cout << " 0. E= " << emin << " prob= " << problast << G4endl;
112  }
113  static const G4double edeltamin = 0.2*CLHEP::MeV;
114  static const G4double edeltamax = 2*CLHEP::MeV;
115  G4bool peak = true;
116  do {
117  ++i;
118  x += edelta;
119  if(x > emax) {
120  x = emax;
121  edelta = emax - fEner[i-1];
122  }
124  if(fVerbose > 1) {
125  G4cout << " " << i << ". E= " << x << " prob= " << y
126  << " Edel= " << edelta << G4endl;
127  }
128  fEner[i] = x;
129  fProb[i] = y;
130  del = (y + problast)*edelta*0.5;
131  totProbability += del;
132 
133  // smart step definition
134  if(del < accuracy*totProbability) { break; }
135  else if(del != totProbability && del > 0.8*totProbability
136  && edelta > edeltamin)
137  { edelta *= 0.5; }
138  else if(del < 0.1*totProbability && edelta < edeltamax)
139  { edelta *= 2.0; }
140  if(y > probmax) {
141  probmax = y;
142  eprobmax = x;
143  } else if(peak && y < probmax) {
144  edelta *= 2.0;
145  peak = false;
146  }
147  problast = y;
148  // Loop checking, 10-Mar-2017, Vladimir Ivanchenko
149  } while(i < nbin && x < emax);
150 
151  nfilled = i;
152  if(fVerbose > 1) { G4cout << " Probability= " << totProbability << G4endl; }
153  return totProbability;
154 }
155 
157 {
158  static const G4double fact = 1.1;
159  probmax *= fact;
160  if(0.0 == fProb[nfilled] && nfilled > 2) { --nfilled; }
161  for(size_t i=0; i<=nfilled; ++i) {
162  fProb[i] *= fact;
163  }
164  G4double ekin(0.0), s1(1.0), s2(0.0), ksi(0.0), psi(1.0), z(0.0);
165 
166  if(fVerbose > 1) {
167  G4cout << "### G4VEmissionProbability::SampleEnergy: "
168  << " Emin= " << emin << " Emax= " << emax
169  << " Nf= " << nfilled << G4endl;
170  }
171  G4double x0 = emax;
172  G4double x1 = fEner[nfilled-1];
174  G4double y0 = probmax;
175  G4double y1 = fProb[nfilled-1];
177  G4bool islog(false);
178  // a condition if a special treatment of falling down
179  // distribution is needed
180  //
181  G4double ee = 0.5*(emax - emin);
182  if(nfilled > 5 && y2 > 0.0 && y2 < 0.1*y0 && y1 > 0.0 && x1 - emin < ee) {
183  ksi = G4Log(y1/y2)/G4Log(x2/x1);
184  x0 = x2*G4Exp(-G4Log(y0/y2)/ksi);
185  // general condition to have two sampling area
186  if(x0 < x1 && x0 > eprobmax && x0 - emin < ee) {
187  // condition when the first area does not exist
188  if(x0 <= emin) {
189  s1 = 0.0;
190  s2 = 1.0;
191  y0 *= G4Exp(G4Log(x0/emin)*ksi);
192  x0 = emin;
193  } else {
194  s1 = (x0 - emin)*y0;
195  }
196  // parameters of the second area
197  if(std::abs(1.0 - ksi) < 0.1) {
198  z = G4Log(emax/x0);
199  if(s1 > 0.0) { s2 = x0*y0*z; }
200  islog = true;
201  } else {
202  psi = 1.0/(1.0 - ksi);
203  z = G4Exp(G4Log(emax/x0)*(1. - ksi)) - 1.;
204  if(s1 > 0.0) { s2 = std::max(y0*x0*z*psi, 0.0); }
205  }
206  }
207  }
208  G4double sum = s1 + s2;
209  if(fVerbose > 1) {
210  G4cout << " Epeak= " << eprobmax << " e0= " << x0 << " e1= " << x1 << " e2= " << x2
211  << " psi= " << psi << G4endl;
212  G4cout << " s1= " << s1 << " s2= " << s2
213  << " y0= " << y0 << " y1= " << y1 << " y2= " << y2 << " z= " << z
214  << G4endl;
215  for(size_t i=0; i<=nfilled; ++i) {
216  G4cout << i << ". E= " << fEner[i] << " P= " << fProb[i] << G4endl;
217  }
218  }
219 
220  CLHEP::HepRandomEngine* rndm = G4Random::getTheEngine();
221  static const G4int nmax = 100;
222  G4double gr, g;
223  for(size_t i=0; i<nmax; ++i) {
224  G4double q = sum*rndm->flat();
225  if(s2 <= 0.0 || q <= s1) {
226  gr = y0;
227  ekin = emin + (x0 - emin)*q/s1;
228  } else if(islog) {
229  ekin = x0*G4Exp((q - s1)*z/s2);
230  gr = y0*x0/ekin;
231  } else {
232  ekin = x0*G4Exp(G4Log(1.0 + (q - s1)*z/s2)*psi);
233  gr = y0*G4Exp(G4Log(x0/ekin)*ksi);
234  }
235  g = ComputeProbability(ekin, eCoulomb);
236  if(fVerbose > 1) {
237  G4cout << " " << i
238  << ". prob= " << g << " probmax= " << gr
239  << " Ekin= " << ekin << G4endl;
240  }
241  if(g > gr && fVerbose > 0) {
242  G4cout << "### G4VEmissionProbability::SampleEnergy Warning i= " << i
243  << " prob/probmax= " << g/gr << " rndm= " << q
244  << "\n prob= " << g << " probmax= " << gr
245  << " z= " << z << " ksi= " << ksi
246  << "\n Ekin= " << ekin << " Emin= " << emin
247  << " Emax= " << emax << " E0= " << x0 << G4endl;
248  G4cout << " s1= " << s1 << " s2= " << s2
249  << " y0= " << y0 << " y1= " << y1 << " y2= " << y2 << G4endl;
250  for(size_t j=0; j<=nfilled; ++j) {
251  G4cout << j << ". E= " << fEner[j] << " P= " << fProb[j] << G4endl;
252  }
253  }
254  if(gr*rndm->flat() <= g) { break; }
255  }
256  return ekin;
257 }
258 
Float_t x
Definition: compare.C:6
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
static constexpr double MeV
Definition: G4SIunits.hh:214
Float_t y1[n_points_granero]
Definition: compare.C:5
G4DeexPrecoParameters * GetParameters()
Float_t x1[n_points_granero]
Definition: compare.C:5
#define G4endl
Definition: G4ios.hh:61
Float_t y
Definition: compare.C:6
static G4NuclearLevelData * GetInstance()
Double_t z
G4PairingCorrection * fPairCorr
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double IntegrateProbability(G4double elow, G4double ehigh, G4double CB)
Float_t y2[n_points_geant4]
Definition: compare.C:26
const G4int nmax
Float_t Z
static constexpr double g
Definition: G4SIunits.hh:183
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
static G4Pow * GetInstance()
Definition: G4Pow.cc:57
void ResetIntegrator(size_t nbin, G4double de, G4double eps)
virtual double flat()=0
static constexpr double MeV
double A(double temperature)
virtual G4double ComputeProbability(G4double anEnergy, G4double CB)
G4VEmissionProbability(G4int Z, G4int A)
G4double GetLevelDensity() const
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
static const G4double eps
G4PairingCorrection * GetPairingCorrection()
Float_t x2[n_points_geant4]
Definition: compare.C:26
std::vector< G4double > fEner
std::vector< G4double > fProb