Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4DNAChampionElasticModel.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: G4DNAChampionElasticModel.cc 105719 2017-08-16 12:36:37Z gcosmo $
27 //
28 
30 #include "G4PhysicalConstants.hh"
31 #include "G4SystemOfUnits.hh"
33 #include "G4Exp.hh"
34 
35 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
36 
37 using namespace std;
38 
39 #define CHAMPION_VERBOSE
40 
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42 
45  G4VEmModel(nam), isInitialised(false)
46 {
47  SetLowEnergyLimit(7.4 * eV);
48  SetHighEnergyLimit(1. * MeV);
49 
50  verboseLevel = 0;
51  // Verbosity scale:
52  // 0 = nothing
53  // 1 = warning for energy non-conservation
54  // 2 = details of energy budget
55  // 3 = calculation of cross sections, file openings, sampling of atoms
56  // 4 = entering in methods
57 
58 #ifdef CHAMPION_VERBOSE
59  if (verboseLevel > 0)
60  {
61  G4cout << "Champion Elastic model is constructed "
62  << G4endl
63  << "Energy range: "
64  << LowEnergyLimit() / eV << " eV - "
65  << HighEnergyLimit() / MeV << " MeV"
66  << G4endl;
67  }
68 #endif
69 
72  fpData = 0;
73 }
74 
75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76 
78 {
79  // For total cross section
80  if(fpData) delete fpData;
81 
82  // For final state
83  eVecm.clear();
84 }
85 
86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
87 
89  const G4DataVector& /*cuts*/)
90 {
91 #ifdef CHAMPION_VERBOSE
92  if (verboseLevel > 3)
93  {
94  G4cout << "Calling G4DNAChampionElasticModel::Initialise()" << G4endl;
95  }
96 #endif
97 
98  if(particle->GetParticleName() != "e-")
99  {
100  G4Exception("G4DNAChampionElasticModel::Initialise",
101  "em0002",
103  "Model not applicable to particle type.");
104  }
105 
106  // Energy limits
107 
108  if (LowEnergyLimit() < 7.4*eV)
109  {
110  G4cout << "G4DNAChampionElasticModel: low energy limit increased from "
111  << LowEnergyLimit()/eV << " eV to " << 7.4 << " eV"
112  << G4endl;
113  SetLowEnergyLimit(7.4*eV);
114  }
115 
116  if (HighEnergyLimit() > 1.*MeV)
117  {
118  G4cout << "G4DNAChampionElasticModel: high energy limit decreased from "
119  << HighEnergyLimit()/MeV << " MeV to " << 1. << " MeV"
120  << G4endl;
122  }
123 
124  if (isInitialised) { return; }
125 
126  // *** ELECTRON
127  // For total cross section
128  // Reading of data files
129 
130  G4double scaleFactor = 1e-16*cm*cm;
131 
132  G4String fileElectron("dna/sigma_elastic_e_champion");
133 
135  eV,
136  scaleFactor );
137  fpData->LoadData(fileElectron);
138 
139  // For final state
140 
141  char *path = getenv("G4LEDATA");
142 
143  if (!path)
144  {
145  G4Exception("G4ChampionElasticModel::Initialise",
146  "em0006",
148  "G4LEDATA environment variable not set.");
149  return;
150  }
151 
152  std::ostringstream eFullFileName;
153  eFullFileName << path << "/dna/sigmadiff_cumulated_elastic_e_champion_hp.dat";
154  std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
155 
156  if (!eDiffCrossSection)
157  {
158  G4ExceptionDescription errMsg;
159  errMsg << "Missing data file:/dna/sigmadiff_cumulated_elastic_e_champion_hp.dat; "
160  << "please use G4EMLOW6.36 and above.";
161 
162  G4Exception("G4DNAChampionElasticModel::Initialise",
163  "em0003",
165  errMsg);
166  }
167 
168  // March 25th, 2014 - Vaclav Stepan, Sebastien Incerti
169  // Added clear for MT
170 
171  eTdummyVec.clear();
172  eVecm.clear();
173  eDiffCrossSectionData.clear();
174 
175  //
176 
177  eTdummyVec.push_back(0.);
178 
179  while(!eDiffCrossSection.eof())
180  {
181  G4double tDummy;
182  G4double eDummy;
183  eDiffCrossSection >> tDummy >> eDummy;
184 
185  // SI : mandatory eVecm initialization
186 
187  if (tDummy != eTdummyVec.back())
188  {
189  eTdummyVec.push_back(tDummy);
190  eVecm[tDummy].push_back(0.);
191  }
192 
193  eDiffCrossSection >> eDiffCrossSectionData[tDummy][eDummy];
194 
195  if (eDummy != eVecm[tDummy].back()) eVecm[tDummy].push_back(eDummy);
196  }
197 
198  // End final state
199 #ifdef CHAMPION_VERBOSE
200  if (verboseLevel>0)
201  {
202  if (verboseLevel > 2)
203  {
204  G4cout << "Loaded cross section files for Champion Elastic model" << G4endl;
205  }
206 
207  G4cout << "Champion Elastic model is initialized " << G4endl
208  << "Energy range: "
209  << LowEnergyLimit() / eV << " eV - "
210  << HighEnergyLimit() / MeV << " MeV"
211  << G4endl;
212  }
213 #endif
214 
215  // Initialize water density pointer
217 
219  GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
220 
222  isInitialised = true;
223 
224 }
225 
226 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
227 
228 G4double
231 #ifdef CHAMPION_VERBOSE
232  const G4ParticleDefinition* p,
233 #else
234  const G4ParticleDefinition*,
235 #endif
236  G4double ekin,
237  G4double,
238  G4double)
239 {
240 #ifdef CHAMPION_VERBOSE
241  if (verboseLevel > 3)
242  {
243  G4cout << "Calling CrossSectionPerVolume() of G4DNAChampionElasticModel"
244  << G4endl;
245  }
246 #endif
247 
248  // Calculate total cross section for model
249 
250  G4double sigma = 0.;
251  G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
252 
253  if(waterDensity!= 0.0)
254  {
255  if (ekin < HighEnergyLimit() && ekin >= LowEnergyLimit())
256  {
257  //SI : XS must not be zero otherwise sampling of secondaries method ignored
258  //
259  sigma = fpData->FindValue(ekin);
260  }
261 
262 #ifdef CHAMPION_VERBOSE
263  if (verboseLevel > 2)
264  {
265  G4cout << "__________________________________" << G4endl;
266  G4cout << "=== G4DNAChampionElasticModel - XS INFO START" << G4endl;
267  G4cout << "=== Kinetic energy(eV)=" << ekin/eV << " particle : " << p->GetParticleName() << G4endl;
268  G4cout << "=== Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
269  G4cout << "=== Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
270  G4cout << "=== G4DNAChampionElasticModel - XS INFO END" << G4endl;
271  }
272 #endif
273  }
274 
275  return sigma*waterDensity;
276 }
277 
278 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
279 
280 void G4DNAChampionElasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
281  const G4MaterialCutsCouple* /*couple*/,
282  const G4DynamicParticle* aDynamicElectron,
283  G4double,
284  G4double)
285 {
286 
287 #ifdef CHAMPION_VERBOSE
288  if (verboseLevel > 3)
289  {
290  G4cout << "Calling SampleSecondaries() of G4DNAChampionElasticModel" << G4endl;
291  }
292 #endif
293 
294  G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
295 
296 // if (electronEnergy0 < HighEnergyLimit()) // necessaire ?
297  {
298  G4double cosTheta = RandomizeCosTheta(electronEnergy0);
299 
300  G4double phi = 2. * pi * G4UniformRand();
301 
302  G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
303  G4ThreeVector xVers = zVers.orthogonal();
304  G4ThreeVector yVers = zVers.cross(xVers);
305 
306  G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
307  G4double yDir = xDir;
308  xDir *= std::cos(phi);
309  yDir *= std::sin(phi);
310 
311  G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
312 
314 
316  // necessaire ?
317  }
318 }
319 
320 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
321 
322 G4double G4DNAChampionElasticModel::Theta(//G4ParticleDefinition * particleDefinition,
323  G4double k,
324  G4double integrDiff)
325 {
326  G4double theta = 0.;
327  G4double valueT1 = 0;
328  G4double valueT2 = 0;
329  G4double valueE21 = 0;
330  G4double valueE22 = 0;
331  G4double valueE12 = 0;
332  G4double valueE11 = 0;
333  G4double xs11 = 0;
334  G4double xs12 = 0;
335  G4double xs21 = 0;
336  G4double xs22 = 0;
337 
338 // if (particleDefinition == G4Electron::ElectronDefinition()) // necessaire ?
339  {
340  std::vector<G4double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),
341  eTdummyVec.end(), k);
342  std::vector<G4double>::iterator t1 = t2 - 1;
343 
344  std::vector<G4double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),
345  eVecm[(*t1)].end(),
346  integrDiff);
347  std::vector<G4double>::iterator e11 = e12 - 1;
348 
349  std::vector<G4double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),
350  eVecm[(*t2)].end(),
351  integrDiff);
352  std::vector<G4double>::iterator e21 = e22 - 1;
353 
354  valueT1 = *t1;
355  valueT2 = *t2;
356  valueE21 = *e21;
357  valueE22 = *e22;
358  valueE12 = *e12;
359  valueE11 = *e11;
360 
361  xs11 = eDiffCrossSectionData[valueT1][valueE11];
362  xs12 = eDiffCrossSectionData[valueT1][valueE12];
363  xs21 = eDiffCrossSectionData[valueT2][valueE21];
364  xs22 = eDiffCrossSectionData[valueT2][valueE22];
365  }
366 
367  if (xs11 == 0 && xs12 == 0 && xs21 == 0 && xs22 == 0) return (0.);
368 
369  theta = QuadInterpolator(valueE11, valueE12, valueE21, valueE22, xs11, xs12,
370  xs21, xs22, valueT1, valueT2, k, integrDiff);
371 
372  return theta;
373 }
374 
375 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
376 
378  G4double e2,
379  G4double e,
380  G4double xs1,
381  G4double xs2)
382 {
383  G4double d1 = std::log(xs1);
384  G4double d2 = std::log(xs2);
385  G4double value = G4Exp(d1 + (d2 - d1) * (e - e1) / (e2 - e1));
386  return value;
387 }
388 
389 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
390 
392  G4double e2,
393  G4double e,
394  G4double xs1,
395  G4double xs2)
396 {
397  G4double d1 = xs1;
398  G4double d2 = xs2;
399  G4double value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
400  return value;
401 }
402 
403 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
404 
406  G4double e2,
407  G4double e,
408  G4double xs1,
409  G4double xs2)
410 {
411  G4double a = (std::log10(xs2) - std::log10(xs1))
412  / (std::log10(e2) - std::log10(e1));
413  G4double b = std::log10(xs2) - a * std::log10(e2);
414  G4double sigma = a * std::log10(e) + b;
415  G4double value = (std::pow(10., sigma));
416  return value;
417 }
418 
419 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
420 
422  G4double e12,
423  G4double e21,
424  G4double e22,
425  G4double xs11,
426  G4double xs12,
427  G4double xs21,
428  G4double xs22,
429  G4double t1,
430  G4double t2,
431  G4double t,
432  G4double e)
433 {
434  // Log-Log
435  /*
436  G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
437  G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
438  G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
439 
440 
441  // Lin-Log
442  G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
443  G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
444  G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
445  */
446 
447  // Lin-Lin
448  G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12);
449  G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22);
450  G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1,
451  interpolatedvalue2);
452 
453  return value;
454 }
455 
456 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
457 
459 {
460 
461  G4double integrdiff = 0;
462  G4double uniformRand = G4UniformRand();
463  integrdiff = uniformRand;
464 
465  G4double theta = 0.;
466  G4double cosTheta = 0.;
467  theta = Theta(//G4Electron::ElectronDefinition(),
468  k / eV, integrdiff);
469 
470  cosTheta = std::cos(theta * pi / 180);
471 
472  return cosTheta;
473 }
474 
475 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
476 
478 {
479  G4ExceptionDescription errMsg;
480  errMsg << "The method G4DNAChampionElasticModel::SetKillBelowThreshold is deprecated";
481 
482  G4Exception("G4DNAChampionElasticModel::SetKillBelowThreshold",
483  "deprecated",
484  JustWarning,
485  errMsg);
486 }
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
virtual G4bool LoadData(const G4String &argFileName)
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:45
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4ParticleChangeForGamma * fParticleChangeForGamma
TTree * t1
Definition: plottest35.C:26
static constexpr double MeV
Definition: G4SIunits.hh:214
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:707
size_t GetIndex() const
Definition: G4Material.hh:262
#define G4endl
Definition: G4ios.hh:61
const G4ThreeVector & GetMomentumDirection() const
const char * p
Definition: xmltok.h:285
const G4String & GetParticleName() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:608
void SetKillBelowThreshold(G4double threshold)
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:616
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4DNAChampionElasticModel(const G4ParticleDefinition *p=0, const G4String &nam="DNAChampionElasticModel")
double G4double
Definition: G4Types.hh:76
static const G4double d2
const XML_Char int const XML_Char * value
Definition: expat.h:331
static G4DNAMolecularMaterial * Instance()
static const G4double d1
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
const std::vector< G4double > * fpMolWaterDensity
#define G4UniformRand()
Definition: Randomize.hh:53
static constexpr double eV
Definition: G4SIunits.hh:215
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
G4double LinLinInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
Hep3Vector unit() const
Hep3Vector cross(const Hep3Vector &) const
TTree * t2
Definition: plottest35.C:36
G4double LinLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
G4double QuadInterpolator(G4double e11, G4double e12, G4double e21, G4double e22, G4double x11, G4double x12, G4double x21, G4double x22, G4double t1, G4double t2, G4double t, G4double e)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
G4double LogLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:714
G4DNACrossSectionDataSet * fpData
G4double GetKineticEnergy() const
virtual G4double FindValue(G4double e, G4int componentId=0) const
static constexpr double cm
Definition: G4SIunits.hh:119
G4GLOB_DLL std::ostream G4cout
std::vector< G4double > eTdummyVec
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
Hep3Vector orthogonal() const
static constexpr double pi
Definition: G4SIunits.hh:75
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:131
G4double Theta(G4double k, G4double integrDiff)
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:609
#define CHAMPION_VERBOSE