ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNAPTBElasticModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DNAPTBElasticModel.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 // Authors: S. Meylan and C. Villagrasa (IRSN, France)
27 // Models come from
28 // M. Bug et al, Rad. Phys and Chem. 130, 459-479 (2017)
29 //
30 
31 #include "G4DNAPTBElasticModel.hh"
32 #include "G4DNAChampionElasticModel.hh"
33 #include "G4PhysicalConstants.hh"
34 #include "G4SystemOfUnits.hh"
35 #include "G4DNAMolecularMaterial.hh"
36 #include "G4Proton.hh"
37 
38 G4DNAPTBElasticModel::G4DNAPTBElasticModel(const G4String& applyToMaterial, const G4ParticleDefinition*,
39  const G4String& nam)
40  : G4VDNAModel(nam, applyToMaterial)
41 {
42  fKillBelowEnergy = 10*eV; // will be override by the limits defined for each material
43 
44  verboseLevel= 0;
45  // Verbosity scale:
46  // 0 = nothing
47  // 1 = warning for energy non-conservation
48  // 2 = details of energy budget
49  // 3 = calculation of cross sections, file openings, sampling of atoms
50  // 4 = entering in methods
51 
52  if( verboseLevel>0 )
53  {
54  G4cout << "PTB Elastic model is constructed " << G4endl;
55  }
56 }
57 
58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
59 
60 G4DNAPTBElasticModel::~G4DNAPTBElasticModel()
61 {
62 
63 }
64 
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66 
67 void G4DNAPTBElasticModel::Initialise(const G4ParticleDefinition* particle,
68  const G4DataVector& /*cuts*/, G4ParticleChangeForGamma*)
69 {
70  if (verboseLevel > 3)
71  G4cout << "Calling G4DNAPTBElasticModel::Initialise()" << G4endl;
72 
73  G4double scaleFactor = 1e-16*cm*cm;
74 
75  G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
76 
77  //*******************************************************
78  // Cross section data
79  //*******************************************************
80 
81  if(particle == electronDef)
82  {
83  G4String particleName = particle->GetParticleName();
84 
85  AddCrossSectionData("THF",
86  particleName,
87  "dna/sigma_elastic_e-_PTB_THF",
88  "dna/sigmadiff_cumulated_elastic_e-_PTB_THF",
89  scaleFactor);
90  SetLowELimit("THF", particleName, 10*eV);
91  SetHighELimit("THF", particleName, 1*keV);
92 
93  AddCrossSectionData("PY",
94  particleName,
95  "dna/sigma_elastic_e-_PTB_PY",
96  "dna/sigmadiff_cumulated_elastic_e-_PTB_PY",
97  scaleFactor);
98  SetLowELimit("PY", particleName, 10*eV);
99  SetHighELimit("PY", particleName, 1*keV);
100 
101  AddCrossSectionData("PU",
102  particleName,
103  "dna/sigma_elastic_e-_PTB_PU",
104  "dna/sigmadiff_cumulated_elastic_e-_PTB_PU",
105  scaleFactor);
106  SetLowELimit("PU", particleName, 10*eV);
107  SetHighELimit("PU", particleName, 1*keV);
108 
109  AddCrossSectionData("TMP",
110  particleName,
111  "dna/sigma_elastic_e-_PTB_TMP",
112  "dna/sigmadiff_cumulated_elastic_e-_PTB_TMP",
113  scaleFactor);
114  SetLowELimit("TMP", particleName, 10*eV);
115  SetHighELimit("TMP", particleName, 1*keV);
116 
117  AddCrossSectionData("G4_WATER",
118  particleName,
119  "dna/sigma_elastic_e_champion",
120  "dna/sigmadiff_cumulated_elastic_e_champion",
121  scaleFactor);
122  SetLowELimit("G4_WATER", particleName, 10*eV);
123  SetHighELimit("G4_WATER", particleName, 1*keV);
124 
125  // DNA materials
126  //
127  AddCrossSectionData("backbone_THF",
128  particleName,
129  "dna/sigma_elastic_e-_PTB_THF",
130  "dna/sigmadiff_cumulated_elastic_e-_PTB_THF",
131  scaleFactor*33./30);
132  SetLowELimit("backbone_THF", particleName, 10*eV);
133  SetHighELimit("backbone_THF", particleName, 1*keV);
134 
135  AddCrossSectionData("cytosine_PY",
136  particleName,
137  "dna/sigma_elastic_e-_PTB_PY",
138  "dna/sigmadiff_cumulated_elastic_e-_PTB_PY",
139  scaleFactor*42./30);
140  SetLowELimit("cytosine_PY", particleName, 10*eV);
141  SetHighELimit("cytosine_PY", particleName, 1*keV);
142 
143  AddCrossSectionData("thymine_PY",
144  particleName,
145  "dna/sigma_elastic_e-_PTB_PY",
146  "dna/sigmadiff_cumulated_elastic_e-_PTB_PY",
147  scaleFactor*48./30);
148  SetLowELimit("thymine_PY", particleName, 10*eV);
149  SetHighELimit("thymine_PY", particleName, 1*keV);
150 
151  AddCrossSectionData("adenine_PU",
152  particleName,
153  "dna/sigma_elastic_e-_PTB_PU",
154  "dna/sigmadiff_cumulated_elastic_e-_PTB_PU",
155  scaleFactor*50./44);
156  SetLowELimit("adenine_PU", particleName, 10*eV);
157  SetHighELimit("adenine_PU", particleName, 1*keV);
158 
159  AddCrossSectionData("guanine_PU",
160  particleName,
161  "dna/sigma_elastic_e-_PTB_PU",
162  "dna/sigmadiff_cumulated_elastic_e-_PTB_PU",
163  scaleFactor*56./44);
164  SetLowELimit("guanine_PU", particleName, 10*eV);
165  SetHighELimit("guanine_PU", particleName, 1*keV);
166 
167  AddCrossSectionData("backbone_TMP",
168  particleName,
169  "dna/sigma_elastic_e-_PTB_TMP",
170  "dna/sigmadiff_cumulated_elastic_e-_PTB_TMP",
171  scaleFactor*33./50);
172  SetLowELimit("backbone_TMP", particleName, 10*eV);
173  SetHighELimit("backbone_TMP", particleName, 1*keV);
174  }
175 
176  //*******************************************************
177  // Load the data
178  //*******************************************************
179 
180  LoadCrossSectionData(particle->GetParticleName() );
181 
182  //*******************************************************
183  // Verbose output
184  //*******************************************************
185 
186  if (verboseLevel > 2)
187  G4cout << "Loaded cross section files for PTB Elastic model" << G4endl;
188 
189  if( verboseLevel>0 )
190  {
191  G4cout << "PTB Elastic model is initialized " << G4endl;
192  }
193 }
194 
195 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
196 
197 void G4DNAPTBElasticModel::ReadDiffCSFile(const G4String& materialName,
198  const G4String& particleName,
199  const G4String& file,
200  const G4double)
201 {
202  // Method to read and save the information contained within the differential cross section files.
203  // This method is not yet standard.
204 
205  // get the path of the G4LEDATA data folder
206  char *path = std::getenv("G4LEDATA");
207  // if it is not found then quit and print error message
208  if(!path)
209  {
210  G4Exception("G4DNAPTBElasticModel::ReadAllDiffCSFiles","em0006",
211  FatalException,"G4LEDATA environment variable not set.");
212  return;
213  }
214 
215  // build the fullFileName path of the data file
216  std::ostringstream fullFileName;
217  fullFileName << path <<"/"<< file<<".dat";
218 
219  // open the data file
220  std::ifstream diffCrossSection (fullFileName.str().c_str());
221  // error if file is not there
222  std::stringstream endPath;
223  if (!diffCrossSection)
224  {
225  endPath << "Missing data file: "<<file;
226  G4Exception("G4DNAPTBElasticModel::Initialise","em0003",
227  FatalException, endPath.str().c_str());
228  }
229 
230  tValuesVec[materialName][particleName].push_back(0.);
231 
232  G4String line;
233 
234  // read the file line by line until we reach the end of file point
235  while(std::getline(diffCrossSection, line))
236  {
237  // check if the line is comment or empty
238  //
239  std::istringstream testIss(line);
240  G4String test;
241  testIss >> test;
242  // check first caracter to determine if following information is data or comments
243  if(test=="#")
244  {
245  // skip the line by beginning a new while loop.
246  continue;
247  }
248  // check if line is empty
249  else if(line.empty())
250  {
251  // skip the line by beginning a new while loop.
252  continue;
253  }
254  //
255  // end of the check
256 
257  // transform the line into a iss
258  std::istringstream iss(line);
259 
260  // Variables to be filled by the input file
261  double tDummy;
262  double eDummy;
263 
264  // fill the variables with the content of the line
265  iss>>tDummy>>eDummy;
266 
267  // SI : mandatory Vecm initialization
268 
269  // Fill two vectors contained in maps of types:
270  // [materialName][particleName]=vector
271  // [materialName][particleName][T]=vector
272  // to list all the incident energies (tValues) and all the output energies (eValues) within the file
273  //
274  // Check if we already have the current T value in the vector.
275  // If not then add it
276  if (tDummy != tValuesVec[materialName][particleName].back())
277  {
278  // Add the current T value
279  tValuesVec[materialName][particleName].push_back(tDummy);
280 
281  // Make it correspond to a default zero E value
282  eValuesVect[materialName][particleName][tDummy].push_back(0.);
283  }
284 
285  // Put the differential cross section value of the input file within the diffCrossSectionData map
286  iss>>diffCrossSectionData[materialName][particleName][tDummy][eDummy];
287 
288  // If the current E value (eDummy) is different from the one already registered in the eVector then add it to the vector
289  if (eDummy != eValuesVect[materialName][particleName][tDummy].back()) eValuesVect[materialName][particleName][tDummy].push_back(eDummy);
290  }
291 }
292 
293 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
294 
295 G4double G4DNAPTBElasticModel::CrossSectionPerVolume(const G4Material* /*material*/,
296  const G4String& materialName,
297  const G4ParticleDefinition* p,
298  G4double ekin,
299  G4double /*emin*/,
300  G4double /*emax*/)
301 {
302  if (verboseLevel > 3)
303  G4cout << "Calling CrossSectionPerVolume() of G4DNAPTBElasticModel" << G4endl;
304 
305  // Get the name of the current particle
306  const G4String& particleName = p->GetParticleName();
307 
308  // set killBelowEnergy value for current material
309  fKillBelowEnergy = GetLowELimit(materialName, particleName);
310 
311  // initialise the return value (cross section) to zero
312  G4double sigma(0);
313 
314  // check if we are below the high energy limit
315  if (ekin < GetHighELimit(materialName, particleName) )
316  {
317  // This is used to kill the particle if its kinetic energy is below fKillBelowEnergy.
318  // If the energy is lower then we return a maximum cross section and thus the SampleSecondaries method will be called for sure.
319  // SampleSecondaries will remove the particle from the simulation.
320  //
321  //SI : XS must not be zero otherwise sampling of secondaries method ignored
322  if (ekin < fKillBelowEnergy) return DBL_MAX;
323 
324  // Get the tables with the cross section data
325  TableMapData* tableData = GetTableData();
326 
327  // Retrieve the cross section value
328  sigma = (*tableData)[materialName][particleName]->FindValue(ekin);
329  }
330 
331  if (verboseLevel > 2)
332  {
333  G4cout << "__________________________________" << G4endl;
334  G4cout << "°°° G4DNAPTBElasticModel - XS INFO START" << G4endl;
335  G4cout << "°°° Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
336  G4cout << "°°° Cross section per molecule (cm^2)=" << sigma/cm/cm << G4endl;
337  G4cout << "°°° G4DNAPTBElasticModel - XS INFO END" << G4endl;
338  }
339 
340  // Return the cross section
341  return sigma;
342 }
343 
344 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
345 
346 void G4DNAPTBElasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
347  const G4MaterialCutsCouple* /*couple*/,
348  const G4String& materialName,
349  const G4DynamicParticle* aDynamicElectron,
350  G4ParticleChangeForGamma* particleChangeForGamma,
351  G4double /*tmin*/,
352  G4double /*tmax*/)
353 {
354  if (verboseLevel > 3)
355  G4cout << "Calling SampleSecondaries() of G4DNAPTBElasticModel" << G4endl;
356 
357  G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
358 
359  const G4String& particleName = aDynamicElectron->GetParticleDefinition()->GetParticleName();
360 
361  // set killBelowEnergy value for material
362  fKillBelowEnergy = GetLowELimit(materialName, particleName);
363 
364  // If the particle (electron here) energy is below the kill limit then we remove it from the simulation
365  if (electronEnergy0 < fKillBelowEnergy)
366  {
367  particleChangeForGamma->SetProposedKineticEnergy(0.);
368  particleChangeForGamma->ProposeTrackStatus(fStopAndKill);
369  particleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0);
370  }
371  // If we are above the kill limite and below the high limit then we proceed
372  else if (electronEnergy0>= fKillBelowEnergy && electronEnergy0 < GetHighELimit(materialName, particleName) )
373  {
374  // Random sampling of the cosTheta
375  G4double cosTheta = RandomizeCosTheta(electronEnergy0, materialName);
376 
377  // Random sampling of phi
378  G4double phi = 2. * pi * G4UniformRand();
379 
380  G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
381  G4ThreeVector xVers = zVers.orthogonal();
382  G4ThreeVector yVers = zVers.cross(xVers);
383 
384  G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
385  G4double yDir = xDir;
386  xDir *= std::cos(phi);
387  yDir *= std::sin(phi);
388 
389  // Particle direction after ModelInterface
390  G4ThreeVector zPrikeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
391 
392  // Give the new direction
393  particleChangeForGamma->ProposeMomentumDirection(zPrikeVers.unit()) ;
394 
395  // Update the energy which does not change here
396  particleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
397  }
398 }
399 
400 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
401 
402 G4double G4DNAPTBElasticModel::Theta
403 (G4ParticleDefinition * particleDefinition, G4double k, G4double integrDiff, const G4String& materialName)
404 {
405  G4double theta = 0.;
406  G4double valueT1 = 0;
407  G4double valueT2 = 0;
408  G4double valueE21 = 0;
409  G4double valueE22 = 0;
410  G4double valueE12 = 0;
411  G4double valueE11 = 0;
412  G4double xs11 = 0;
413  G4double xs12 = 0;
414  G4double xs21 = 0;
415  G4double xs22 = 0;
416  G4String particleName = particleDefinition->GetParticleName();
417 
418  if (particleDefinition == G4Electron::ElectronDefinition())
419  {
420  std::vector<double>::iterator t2 = std::upper_bound(tValuesVec[materialName][particleName].begin(),tValuesVec[materialName][particleName].end(), k);
421  std::vector<double>::iterator t1 = t2-1;
422 
423  std::vector<double>::iterator e12 = std::upper_bound(eValuesVect[materialName][particleName][(*t1)].begin(),eValuesVect[materialName][particleName][(*t1)].end(), integrDiff);
424  std::vector<double>::iterator e11 = e12-1;
425 
426  std::vector<double>::iterator e22 = std::upper_bound(eValuesVect[materialName][particleName][(*t2)].begin(),eValuesVect[materialName][particleName][(*t2)].end(), integrDiff);
427  std::vector<double>::iterator e21 = e22-1;
428 
429  valueT1 =*t1;
430  valueT2 =*t2;
431  valueE21 =*e21;
432  valueE22 =*e22;
433  valueE12 =*e12;
434  valueE11 =*e11;
435 
436  xs11 = diffCrossSectionData[materialName][particleName][valueT1][valueE11];
437  xs12 = diffCrossSectionData[materialName][particleName][valueT1][valueE12];
438  xs21 = diffCrossSectionData[materialName][particleName][valueT2][valueE21];
439  xs22 = diffCrossSectionData[materialName][particleName][valueT2][valueE22];
440  }
441 
442  if (xs11==0 && xs12==0 && xs21==0 && xs22==0) return (0.);
443 
444  theta = QuadInterpolator ( valueE11, valueE12,
445  valueE21, valueE22,
446  xs11, xs12,
447  xs21, xs22,
448  valueT1, valueT2,
449  k, integrDiff );
450 
451  return theta;
452 }
453 
454 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
455 
456 G4double G4DNAPTBElasticModel::LinLogInterpolate(G4double e1,
457  G4double e2,
458  G4double e,
459  G4double xs1,
460  G4double xs2)
461 {
462  G4double d1 = std::log(xs1);
463  G4double d2 = std::log(xs2);
464  G4double value = std::exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
465  return value;
466 }
467 
468 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
469 
470 G4double G4DNAPTBElasticModel::LinLinInterpolate(G4double e1,
471  G4double e2,
472  G4double e,
473  G4double xs1,
474  G4double xs2)
475 {
476  G4double d1 = xs1;
477  G4double d2 = xs2;
478  G4double value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
479  return value;
480 }
481 
482 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
483 
484 G4double G4DNAPTBElasticModel::LogLogInterpolate(G4double e1,
485  G4double e2,
486  G4double e,
487  G4double xs1,
488  G4double xs2)
489 {
490  G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
491  G4double b = std::log10(xs2) - a*std::log10(e2);
492  G4double sigma = a*std::log10(e) + b;
493  G4double value = (std::pow(10.,sigma));
494  return value;
495 }
496 
497 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
498 
499 G4double G4DNAPTBElasticModel::QuadInterpolator(G4double e11, G4double e12,
500  G4double e21, G4double e22,
501  G4double xs11, G4double xs12,
502  G4double xs21, G4double xs22,
503  G4double t1, G4double t2,
504  G4double t, G4double e)
505 {
506  // Log-Log
507  /*
508  G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
509  G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
510  G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
511 
512 
513  // Lin-Log
514  G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
515  G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
516  G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
517 */
518 
519  // Lin-Lin
520  G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12);
521  G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22);
522  G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
523 
524  return value;
525 }
526 
527 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
528 
529 G4double G4DNAPTBElasticModel::RandomizeCosTheta(G4double k, const G4String& materialName)
530 {
531  G4double integrdiff=0;
532  G4double uniformRand=G4UniformRand();
533  integrdiff = uniformRand;
534 
535  G4double theta=0.;
536  G4double cosTheta=0.;
537  theta = Theta(G4Electron::ElectronDefinition(),k/eV,integrdiff, materialName);
538 
539  cosTheta= std::cos(theta*pi/180);
540 
541  return cosTheta;
542 }
543 
544 
545 
546