Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4eeToHadronsModel.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: G4eeToHadronsModel.cc 109567 2018-05-02 07:04:10Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class header file
31 //
32 //
33 // File name: G4eeToHadronsModel
34 //
35 // Author: Vladimir Ivanchenko
36 //
37 // Creation date: 12.08.2003
38 //
39 // Modifications:
40 // 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
41 // 18-05-05 Use optimized interfaces (V.Ivantchenko)
42 //
43 //
44 // -------------------------------------------------------------------
45 //
46 
47 
48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
50 
51 #include "G4eeToHadronsModel.hh"
52 #include "Randomize.hh"
53 #include "G4PhysicalConstants.hh"
54 #include "G4SystemOfUnits.hh"
55 #include "G4Electron.hh"
56 #include "G4Gamma.hh"
57 #include "G4Positron.hh"
58 #include "G4PionPlus.hh"
59 #include "Randomize.hh"
60 #include "G4Vee2hadrons.hh"
61 #include "G4PhysicsVector.hh"
62 #include "G4PhysicsLogVector.hh"
63 #include "G4Log.hh"
64 #include "G4Exp.hh"
65 
66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
67 
68 using namespace std;
69 
71  const G4String& nam)
72  : G4VEmModel(nam),
73  model(mod),
74  crossPerElectron(0),
75  crossBornPerElectron(0),
76  isInitialised(false),
77  nbins(100),
78  verbose(ver)
79 {
86  epeak = emax;
87  //verbose = 1;
88 }
89 
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
91 
93 {
94  delete model;
95 }
96 
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
98 
100  const G4DataVector&)
101 {
102  if(isInitialised) { return; }
103  isInitialised = true;
104 
105  // CM system
106  emin = model->LowEnergy();
107  emax = model->HighEnergy();
108 
109  // peak energy
110  epeak = std::min(model->PeakEnergy(), emax);
111 
112  if(verbose>0) {
113  G4cout << "G4eeToHadronsModel::Initialise: " << G4endl;
114  G4cout << "CM System: emin(MeV)= " << emin/MeV
115  << " epeak(MeV)= " << epeak/MeV
116  << " emax(MeV)= " << emax/MeV
117  << G4endl;
118  }
119 
120  crossBornPerElectron = model->PhysicsVector();
121  crossPerElectron = model->PhysicsVector();
123  for(G4int i=0; i<nbins; ++i) {
125  G4double cs = model->ComputeCrossSection(e);
127  }
129 
130  if(verbose>1) {
131  G4cout << "G4eeToHadronsModel: Cross sections per electron"
132  << " nbins= " << nbins
133  << " emin(MeV)= " << emin/MeV
134  << " emax(MeV)= " << emax/MeV
135  << G4endl;
136  for(G4int i=0; i<nbins; ++i) {
140  G4cout << "E(MeV)= " << e/MeV
141  << " cross(nb)= " << s1/nanobarn
142  << " crossBorn(nb)= " << s2/nanobarn
143  << G4endl;
144  }
145  }
146 }
147 
148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
149 
151  const G4Material* mat,
152  const G4ParticleDefinition* p,
155 {
156  return mat->GetElectronDensity()*
157  ComputeCrossSectionPerElectron(p, kineticEnergy);
158 }
159 
160 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
161 
163  const G4ParticleDefinition* p,
167 {
168  return Z*ComputeCrossSectionPerElectron(p, kineticEnergy);
169 }
170 
171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
172 
174  const G4ParticleDefinition*,
177 {
178  return (crossPerElectron) ? crossPerElectron->Value(energy) : 0.0;
179 }
180 
181 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
182 
183 void G4eeToHadronsModel::SampleSecondaries(std::vector<G4DynamicParticle*>* newp,
184  const G4MaterialCutsCouple*,
185  const G4DynamicParticle* dParticle,
186  G4double,
187  G4double)
188 {
189  if(crossPerElectron) {
190  G4double t = dParticle->GetKineticEnergy() + 2*electron_mass_c2;
191  G4LorentzVector inlv = dParticle->Get4Momentum() +
192  G4LorentzVector(0.0,0.0,0.0,electron_mass_c2);
193  G4double e = inlv.m();
194  G4ThreeVector inBoost = inlv.boostVector();
195  //G4cout << "G4eeToHadronsModel::SampleSecondaries e= " << e
196  // << " " << inlv << " " << inBoost <<G4endl;
197  if(e > emin) {
199  G4LorentzVector gLv = gamma->Get4Momentum();
200  G4LorentzVector lv(0.0,0.0,0.0,e);
201  lv -= gLv;
202  G4double mass = lv.m();
203  //G4cout << "mass= " << mass << " " << lv << G4endl;
204  G4ThreeVector boost = lv.boostVector();
205  //G4cout << "mass= " << mass << " " << boost << G4endl;
206  const G4ThreeVector dir = gamma->GetMomentumDirection();
207  model->SampleSecondaries(newp, mass, dir);
208  G4int np = newp->size();
209  for(G4int j=0; j<np; ++j) {
210  G4DynamicParticle* dp = (*newp)[j];
211  G4LorentzVector v = dp->Get4Momentum();
212  v.boost(boost);
213  //G4cout << j << ". " << v << G4endl;
214  v.boost(inBoost);
215  //G4cout << " " << v << G4endl;
216  dp->Set4Momentum(v);
217  t -= v.e();
218  }
219  //G4cout << "Gamma " << gLv << G4endl;
220  gLv.boost(inBoost);
221  //G4cout << " " << gLv << G4endl;
222  gamma->Set4Momentum(gLv);
223  t -= gLv.e();
224  newp->push_back(gamma);
225  if(std::abs(t) > MeV) {
226  G4cout << "G4eeToHadronsModel::SampleSecondaries: Ebalance(MeV)= "
227  << t/MeV << " primary 4-momentum: " << inlv << G4endl;
228  }
229  }
230  }
231 }
232 
233 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
234 
236 {
237  for(G4int i=0; i<nbins; i++) {
239  G4double cs = 0.0;
240  if(i > 0) {
242  G4double bt = 2.0*fine_structure_const*(LL - 1.0)/pi;
243  G4double btm1= bt - 1.0;
244  G4double del = 1. + fine_structure_const*(1.5*LL + pi*pi/3. -2.)/pi;
246  G4double e1 = crossPerElectron->Energy(i-1);
247  G4double x1 = 1. - e1/e;
248  cs += s1*(del*G4Exp(G4Log(x1)*bt) - bt*(x1 - 0.25*x1*x1));
249  if(i > 1) {
250  G4double e2 = e1;
251  G4double x2 = x1;
253  G4double w2 = bt*(del*G4Exp(G4Log(x2)*btm1) - 1.0 + 0.5*x2);
254  G4double w1;
255 
256  for(G4int j=i-2; j>=0; --j) {
257  e1 = crossPerElectron->Energy(j);
258  x1 = 1. - e1/e;
259  s1 = crossBornPerElectron->Value(e1);
260  w1 = bt*(del*G4Exp(G4Log(x1)*btm1) - 1.0 + 0.5*x1);
261  cs += 0.5*(x1 - x2)*(w2*s2 + w1*s1);
262  e2 = e1;
263  x2 = x1;
264  s2 = s1;
265  w2 = w1;
266  }
267  }
268  }
269  crossPerElectron->PutValue(i, cs);
270  }
271 }
272 
273 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
274 
276 {
277  G4double x;
278  G4DynamicParticle* gamma = nullptr;
280  G4double bt = 2.0*fine_structure_const*(LL - 1.)/pi;
281  G4double btm1= bt - 1.0;
282  G4double del = 1. + fine_structure_const*(1.5*LL + pi*pi/3. -2.)/pi;
283 
285  G4double de = (emax - emin)/(G4double)nbins;
286  G4double xmax = 0.5*(1.0 - (emin*emin)/(e*e));
287  G4double xmin = std::min(de/e, xmax);
288  G4double ds = s0*(del*G4Exp(G4Log(xmin)*bt) - bt*(xmin - 0.25*xmin*xmin));
289  G4double e1 = e*(1. - xmin);
290 
291  //G4cout << "e1= " << e1 << G4endl;
292  if(e1 < emax && s0*G4UniformRand()<ds) {
293  x = xmin*G4Exp(G4Log(G4UniformRand())/bt);
294  } else {
295 
296  x = xmin;
298  G4double w1 = bt*(del*G4Exp(G4Log(x)*btm1) - 1.0 + 0.5*x);
299  G4double grej = s1*w1;
300  G4double f;
301  /*
302  G4cout << "e(GeV)= " << e/GeV << " epeak(GeV)= " << epeak/GeV
303  << " s1= " << s1 << " w1= " << w1
304  << " grej= " << grej << G4endl;
305  */
306  // Above emax cross section is const
307  if(e1 > emax) {
308  x = 0.5*(1. - (emax*emax)/(e*e));
310  G4double w2 = bt*(del*G4Exp(G4Log(x)*btm1) - 1.0 + 0.5*x);
311  grej = s2*w2;
312  //G4cout << "emax= " << emax << " s2= " << s2 << " w2= " << w2
313  // << " grej= " << grej << G4endl;
314  }
315 
316  if(e1 > epeak) {
317  x = 0.5*(1.0 - (epeak*epeak)/(e*e));
319  G4double w2 = bt*(del*G4Exp(G4Log(x)*btm1) - 1.0 + 0.5*x);
320  grej = std::max(grej,s2*w2);
321  //G4cout << "epeak= " << epeak << " s2= " << s2 << " w2= " << w2
322  // << " grej= " << grej << G4endl;
323  }
324  G4int ii = 0;
325  const G4int iimax = 1000;
326  do {
327  x = xmin + G4UniformRand()*(xmax - xmin);
328 
329  G4double s2 = crossBornPerElectron->Value(sqrt(1.0 - 2*x)*e);
330  G4double w2 = bt*(del*G4Exp(G4Log(x)*btm1) - 1.0 + 0.5*x);
331  /*
332  G4cout << "x= " << x << " xmin= " << xmin << " xmax= " << xmax
333  << " s2= " << s2 << " w2= " << w2 << G4endl;
334  */
335  f = s2*w2;
336  if(f > grej) {
337  G4cout << "G4DynamicParticle* G4eeToHadronsModel:WARNING "
338  << f << " > " << grej << " majorant is`small!"
339  << G4endl;
340  }
341  if(++ii >= iimax) { break; }
342  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
343  } while (f < grej*G4UniformRand());
344  }
345 
346  G4ThreeVector dir(0.0,0.0,1.0);
347  if(G4UniformRand() > 0.5) { dir.set(0.0,0.0,-1.0); }
348  //G4cout << "Egamma(MeV)= " << x*e << " " << dir << G4endl;
349  gamma = new G4DynamicParticle(theGamma,dir,x*e);
350  return gamma;
351 }
352 
353 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
354 
Float_t x
Definition: compare.C:6
G4double Energy(size_t index) const
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
static constexpr double MeV
Definition: G4SIunits.hh:214
static const G4int LL[nN]
G4eeToHadronsModel(G4Vee2hadrons *, G4int ver=0, const G4String &nam="eeToHadrons")
Float_t x1[n_points_granero]
Definition: compare.C:5
#define G4endl
Definition: G4ios.hh:61
const G4ThreeVector & GetMomentumDirection() const
const char * p
Definition: xmltok.h:285
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kineticEnergy, G4double Z, G4double A, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX) override
G4double Value(G4double theEnergy, size_t &lastidx) const
static constexpr double nanobarn
Definition: G4SIunits.hh:108
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:616
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
Float_t Z
G4DynamicParticle * GenerateCMPhoton(G4double)
double G4double
Definition: G4Types.hh:76
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin=0.0, G4double maxEnergy=DBL_MAX) override
G4ParticleDefinition * theGamma
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy) override
static constexpr double electron_mass_c2
#define G4UniformRand()
Definition: Randomize.hh:53
virtual G4double ComputeCrossSectionPerElectron(const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
double energy
Definition: plottest35.C:25
void Set4Momentum(const G4LorentzVector &momentum)
G4LorentzVector Get4Momentum() const
Float_t mat
int G4int
Definition: G4Types.hh:78
TDirectory * dir
G4double GetKineticEnergy() const
G4GLOB_DLL std::ostream G4cout
G4PhysicsVector * crossBornPerElectron
Hep3Vector boostVector() const
static constexpr double pi
Definition: G4SIunits.hh:75
CLHEP::HepLorentzVector G4LorentzVector
static constexpr double fine_structure_const
Float_t x2[n_points_geant4]
Definition: compare.C:26
G4PhysicsVector * crossPerElectron
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:609
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
const XML_Char XML_Content * model
Definition: expat.h:151
G4double GetElectronDensity() const
Definition: G4Material.hh:218
void PutValue(size_t index, G4double theValue)
HepLorentzVector & boost(double, double, double)
size_t GetVectorLength() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments