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