ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNAChampionElasticModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 //
27 
28 #include "G4DNAChampionElasticModel.hh"
29 #include "G4PhysicalConstants.hh"
30 #include "G4SystemOfUnits.hh"
31 #include "G4DNAMolecularMaterial.hh"
32 #include "G4Exp.hh"
33 
34 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
35 
36 using namespace std;
37 
38 #define CHAMPION_VERBOSE
39 
40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
41 
42 G4DNAChampionElasticModel::
43 G4DNAChampionElasticModel(const G4ParticleDefinition*, const G4String& nam) :
44  G4VEmModel(nam), isInitialised(false)
45 {
46  SetLowEnergyLimit(7.4 * eV);
47  SetHighEnergyLimit(1. * MeV);
48 
49  verboseLevel = 0;
50  // Verbosity scale:
51  // 0 = nothing
52  // 1 = warning for energy non-conservation
53  // 2 = details of energy budget
54  // 3 = calculation of cross sections, file openings, sampling of atoms
55  // 4 = entering in methods
56 
57 #ifdef CHAMPION_VERBOSE
58  if (verboseLevel > 0)
59  {
60  G4cout << "Champion Elastic model is constructed "
61  << G4endl
62  << "Energy range: "
63  << LowEnergyLimit() / eV << " eV - "
64  << HighEnergyLimit() / MeV << " MeV"
65  << G4endl;
66  }
67 #endif
68 
69  fParticleChangeForGamma = 0;
70  fpMolWaterDensity = 0;
71  fpData = 0;
72 }
73 
74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75 
76 G4DNAChampionElasticModel::~G4DNAChampionElasticModel()
77 {
78  // For total cross section
79  if(fpData) delete fpData;
80 
81  // For final state
82  eVecm.clear();
83 }
84 
85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
86 
87 void G4DNAChampionElasticModel::Initialise(const G4ParticleDefinition* particle,
88  const G4DataVector& /*cuts*/)
89 {
90 #ifdef CHAMPION_VERBOSE
91  if (verboseLevel > 3)
92  {
93  G4cout << "Calling G4DNAChampionElasticModel::Initialise()" << G4endl;
94  }
95 #endif
96 
97  if(particle->GetParticleName() != "e-")
98  {
99  G4Exception("G4DNAChampionElasticModel::Initialise",
100  "em0002",
101  FatalException,
102  "Model not applicable to particle type.");
103  }
104 
105  // Energy limits
106 
107  if (LowEnergyLimit() < 7.4*eV)
108  {
109  G4cout << "G4DNAChampionElasticModel: low energy limit increased from "
110  << LowEnergyLimit()/eV << " eV to " << 7.4 << " eV"
111  << G4endl;
112  SetLowEnergyLimit(7.4*eV);
113  }
114 
115  if (HighEnergyLimit() > 1.*MeV)
116  {
117  G4cout << "G4DNAChampionElasticModel: high energy limit decreased from "
118  << HighEnergyLimit()/MeV << " MeV to " << 1. << " MeV"
119  << G4endl;
120  SetHighEnergyLimit(1.*MeV);
121  }
122 
123  if (isInitialised) { return; }
124 
125  // *** ELECTRON
126  // For total cross section
127  // Reading of data files
128 
129  G4double scaleFactor = 1e-16*cm*cm;
130 
131  G4String fileElectron("dna/sigma_elastic_e_champion");
132 
133  fpData = new G4DNACrossSectionDataSet(new G4LogLogInterpolation(),
134  eV,
135  scaleFactor );
136  fpData->LoadData(fileElectron);
137 
138  // For final state
139 
140  char *path = getenv("G4LEDATA");
141 
142  if (!path)
143  {
144  G4Exception("G4ChampionElasticModel::Initialise",
145  "em0006",
146  FatalException,
147  "G4LEDATA environment variable not set.");
148  return;
149  }
150 
151  std::ostringstream eFullFileName;
152  eFullFileName << path << "/dna/sigmadiff_cumulated_elastic_e_champion_hp.dat";
153  std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
154 
155  if (!eDiffCrossSection)
156  {
157  G4ExceptionDescription errMsg;
158  errMsg << "Missing data file:/dna/sigmadiff_cumulated_elastic_e_champion_hp.dat; "
159  << "please use G4EMLOW6.36 and above.";
160 
161  G4Exception("G4DNAChampionElasticModel::Initialise",
162  "em0003",
163  FatalException,
164  errMsg);
165  }
166 
167  // March 25th, 2014 - Vaclav Stepan, Sebastien Incerti
168  // Added clear for MT
169 
170  eTdummyVec.clear();
171  eVecm.clear();
172  eDiffCrossSectionData.clear();
173 
174  //
175 
176  eTdummyVec.push_back(0.);
177 
178  while(!eDiffCrossSection.eof())
179  {
180  G4double tDummy;
181  G4double eDummy;
182  eDiffCrossSection >> tDummy >> eDummy;
183 
184  // SI : mandatory eVecm initialization
185 
186  if (tDummy != eTdummyVec.back())
187  {
188  eTdummyVec.push_back(tDummy);
189  eVecm[tDummy].push_back(0.);
190  }
191 
192  eDiffCrossSection >> eDiffCrossSectionData[tDummy][eDummy];
193 
194  if (eDummy != eVecm[tDummy].back()) eVecm[tDummy].push_back(eDummy);
195  }
196 
197  // End final state
198 
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
216  G4DNAMolecularMaterial::Instance()->Initialize();
217 
218  fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->
219  GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
220 
221  fParticleChangeForGamma = GetParticleChangeForGamma();
222  isInitialised = true;
223 
224 }
225 
226 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
227 
228 G4double
229 G4DNAChampionElasticModel::
230 CrossSectionPerVolume(const G4Material* material,
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 (ekin <= HighEnergyLimit() && ekin >= LowEnergyLimit())
254  {
255  //SI : XS must not be zero otherwise sampling of secondaries method ignored
256  //
257  sigma = fpData->FindValue(ekin);
258  }
259 
260 #ifdef CHAMPION_VERBOSE
261  if (verboseLevel > 2)
262  {
263  G4cout << "__________________________________" << G4endl;
264  G4cout << "=== G4DNAChampionElasticModel - XS INFO START" << G4endl;
265  G4cout << "=== Kinetic energy(eV)=" << ekin/eV << " particle : " << p->GetParticleName() << G4endl;
266  G4cout << "=== Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
267  G4cout << "=== Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
268  G4cout << "=== G4DNAChampionElasticModel - XS INFO END" << G4endl;
269  }
270 #endif
271 
272  return sigma*waterDensity;
273 }
274 
275 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
276 
277 void G4DNAChampionElasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
278  const G4MaterialCutsCouple* /*couple*/,
279  const G4DynamicParticle* aDynamicElectron,
280  G4double,
281  G4double)
282 {
283 
284 #ifdef CHAMPION_VERBOSE
285  if (verboseLevel > 3)
286  {
287  G4cout << "Calling SampleSecondaries() of G4DNAChampionElasticModel" << G4endl;
288  }
289 #endif
290 
291  G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
292 
293  G4double cosTheta = RandomizeCosTheta(electronEnergy0);
294 
295  G4double phi = 2. * pi * G4UniformRand();
296 
297  G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
298  G4ThreeVector xVers = zVers.orthogonal();
299  G4ThreeVector yVers = zVers.cross(xVers);
300 
301  G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
302  G4double yDir = xDir;
303  xDir *= std::cos(phi);
304  yDir *= std::sin(phi);
305 
306  G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
307 
308  fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit());
309 
310  fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
311 
312 }
313 
314 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
315 
316 G4double G4DNAChampionElasticModel::Theta(G4double k,
317  G4double integrDiff)
318 {
319  G4double theta = 0.;
320  G4double valueT1 = 0;
321  G4double valueT2 = 0;
322  G4double valueE21 = 0;
323  G4double valueE22 = 0;
324  G4double valueE12 = 0;
325  G4double valueE11 = 0;
326  G4double xs11 = 0;
327  G4double xs12 = 0;
328  G4double xs21 = 0;
329  G4double xs22 = 0;
330 
331  std::vector<G4double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),
332  eTdummyVec.end(), k);
333  std::vector<G4double>::iterator t1 = t2 - 1;
334 
335  std::vector<G4double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),
336  eVecm[(*t1)].end(),
337  integrDiff);
338  std::vector<G4double>::iterator e11 = e12 - 1;
339 
340  std::vector<G4double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),
341  eVecm[(*t2)].end(),
342  integrDiff);
343  std::vector<G4double>::iterator e21 = e22 - 1;
344 
345  valueT1 = *t1;
346  valueT2 = *t2;
347  valueE21 = *e21;
348  valueE22 = *e22;
349  valueE12 = *e12;
350  valueE11 = *e11;
351 
352  xs11 = eDiffCrossSectionData[valueT1][valueE11];
353  xs12 = eDiffCrossSectionData[valueT1][valueE12];
354  xs21 = eDiffCrossSectionData[valueT2][valueE21];
355  xs22 = eDiffCrossSectionData[valueT2][valueE22];
356 
357  if (xs11 == 0 && xs12 == 0 && xs21 == 0 && xs22 == 0) return (0.);
358 
359  theta = QuadInterpolator(valueE11, valueE12, valueE21, valueE22, xs11, xs12,
360  xs21, xs22, valueT1, valueT2, k, integrDiff);
361 
362  return theta;
363 }
364 
365 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
366 
367 G4double G4DNAChampionElasticModel::LinLogInterpolate(G4double e1,
368  G4double e2,
369  G4double e,
370  G4double xs1,
371  G4double xs2)
372 {
373  G4double d1 = std::log(xs1);
374  G4double d2 = std::log(xs2);
375  G4double value = G4Exp(d1 + (d2 - d1) * (e - e1) / (e2 - e1));
376  return value;
377 }
378 
379 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
380 
381 G4double G4DNAChampionElasticModel::LinLinInterpolate(G4double e1,
382  G4double e2,
383  G4double e,
384  G4double xs1,
385  G4double xs2)
386 {
387  G4double d1 = xs1;
388  G4double d2 = xs2;
389  G4double value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
390  return value;
391 }
392 
393 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
394 
395 G4double G4DNAChampionElasticModel::LogLogInterpolate(G4double e1,
396  G4double e2,
397  G4double e,
398  G4double xs1,
399  G4double xs2)
400 {
401  G4double a = (std::log10(xs2) - std::log10(xs1))
402  / (std::log10(e2) - std::log10(e1));
403  G4double b = std::log10(xs2) - a * std::log10(e2);
404  G4double sigma = a * std::log10(e) + b;
405  G4double value = (std::pow(10., sigma));
406  return value;
407 }
408 
409 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
410 
411 G4double G4DNAChampionElasticModel::QuadInterpolator(G4double e11,
412  G4double e12,
413  G4double e21,
414  G4double e22,
415  G4double xs11,
416  G4double xs12,
417  G4double xs21,
418  G4double xs22,
419  G4double t1,
420  G4double t2,
421  G4double t,
422  G4double e)
423 {
424  // Log-Log
425  /*
426  G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
427  G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
428  G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
429 
430 
431  // Lin-Log
432  G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
433  G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
434  G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
435  */
436 
437  // Lin-Lin
438  G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12);
439  G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22);
440  G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1,
441  interpolatedvalue2);
442 
443  return value;
444 }
445 
446 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
447 
448 G4double G4DNAChampionElasticModel::RandomizeCosTheta(G4double k)
449 {
450 
451  G4double integrdiff = 0;
452  G4double uniformRand = G4UniformRand();
453  integrdiff = uniformRand;
454 
455  G4double theta = 0.;
456  G4double cosTheta = 0.;
457  theta = Theta(k / eV, integrdiff);
458 
459  cosTheta = std::cos(theta * pi / 180);
460 
461  return cosTheta;
462 }
463 
464 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
465 
466 void G4DNAChampionElasticModel::SetKillBelowThreshold(G4double)
467 {
468  G4ExceptionDescription errMsg;
469  errMsg << "The method G4DNAChampionElasticModel::SetKillBelowThreshold is deprecated";
470 
471  G4Exception("G4DNAChampionElasticModel::SetKillBelowThreshold",
472  "deprecated",
473  JustWarning,
474  errMsg);
475 }