Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4DNAPTBElasticModel.cc
Go to the documentation of this file.
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
34#include "G4SystemOfUnits.hh"
36#include "G4Proton.hh"
37
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
61{
62
63}
64
65//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66
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
76
77 //*******************************************************
78 // Cross section data
79 //*******************************************************
80
81 if(particle == electronDef)
82 {
83 G4String particleName = particle->GetParticleName();
84
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
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
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
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
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
197void 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
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
346void 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
402G4double 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
456G4double 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
470G4double 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
484G4double 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
499G4double 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
529G4double 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
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
@ fStopAndKill
double G4double
Definition: G4Types.hh:83
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:52
Hep3Vector unit() const
Hep3Vector orthogonal() const
Hep3Vector cross(const Hep3Vector &) const
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4String &materialName, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
CrossSectionPerVolume This method is mandatory for any model class. It finds and return the cross sec...
G4DNAPTBElasticModel(const G4String &applyToMaterial="all", const G4ParticleDefinition *p=0, const G4String &nam="DNAPTBElasticModel")
G4DNAPTBElasticModel Constructor.
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4String &materialName, const G4DynamicParticle *, G4ParticleChangeForGamma *particleChangeForGamma, G4double tmin, G4double tmax)
SampleSecondaries Method called after CrossSectionPerVolume if the process is the one which is select...
virtual ~G4DNAPTBElasticModel()
~G4DNAPTBElasticModel Destructor
virtual void Initialise(const G4ParticleDefinition *particle, const G4DataVector &, G4ParticleChangeForGamma *fpChangeForGamme=nullptr)
Initialise Mandatory method for every model class. The material/particle for which the model can be u...
const G4ThreeVector & GetMomentumDirection() const
const G4ParticleDefinition * GetParticleDefinition() const
G4double GetKineticEnergy() const
static G4Electron * ElectronDefinition()
Definition: G4Electron.cc:88
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
const G4String & GetParticleName() const
The G4VDNAModel class.
Definition: G4VDNAModel.hh:50
TableMapData * GetTableData()
GetTableData.
Definition: G4VDNAModel.hh:193
void SetHighELimit(const G4String &material, const G4String &particle, G4double lim)
SetHighEnergyLimit.
Definition: G4VDNAModel.hh:169
std::map< G4String, std::map< G4String, G4DNACrossSectionDataSet *, std::less< G4String > > > TableMapData
Definition: G4VDNAModel.hh:183
void AddCrossSectionData(G4String materialName, G4String particleName, G4String fileCS, G4String fileDiffCS, G4double scaleFactor)
AddCrossSectionData Method used during the initialization of the model class to add a new material....
Definition: G4VDNAModel.cc:58
void SetLowELimit(const G4String &material, const G4String &particle, G4double lim)
SetLowEnergyLimit.
Definition: G4VDNAModel.hh:177
G4double GetHighELimit(const G4String &material, const G4String &particle)
GetHighEnergyLimit.
Definition: G4VDNAModel.hh:153
void LoadCrossSectionData(const G4String &particleName)
LoadCrossSectionData Method to loop on all the registered materials in the model and load the corresp...
Definition: G4VDNAModel.cc:75
G4double GetLowELimit(const G4String &material, const G4String &particle)
GetLowEnergyLimit.
Definition: G4VDNAModel.hh:161
void ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
#define DBL_MAX
Definition: templates.hh:62