Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4DNAELSEPAElasticModel.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// Created on 2016/01/18
27//
28// Authors: D. Sakata, W.G. Shin, S. Incerti
29//
30// Based on a recent release of the ELSEPA code
31// developed and provided kindly by F. Salvat et al.
32// See
33// Computer Physics Communications, 165(2), 157-190. (2005)
34// http://dx.doi.org/10.1016/j.cpc.2004.09.006
35//
36
37#include "G4DNAELSEPAElasticModel.hh"
38#include "G4PhysicalConstants.hh"
39#include "G4SystemOfUnits.hh"
40#include "G4DNAMolecularMaterial.hh"
41
42//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
43
44using namespace std;
45
46//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47
48G4DNAELSEPAElasticModel::G4DNAELSEPAElasticModel(const G4ParticleDefinition*,
49const G4String& nam) :
50G4VEmModel(nam)
51{
52 verboseLevel = 0;
53
54 G4ProductionCutsTable* theCoupleTable =
55 G4ProductionCutsTable::GetProductionCutsTable();
56 auto numOfCouples = (G4int)theCoupleTable->GetTableSize();
57
58 for(G4int i=0; i<numOfCouples; ++i)
59 {
60 const G4MaterialCutsCouple* couple =
61 theCoupleTable->GetMaterialCutsCouple(i);
62 const G4Material* material = couple->GetMaterial();
63 auto nelm = (G4int)material->GetNumberOfElements();
64 const G4ElementVector* theElementVector = material->GetElementVector();
65
66 if(nelm==1)
67 {// Protection: only for single element
68 G4int Z = 79;
69 Z = G4lrint((*theElementVector)[0]->GetZ());
70 // Protection: only for GOLD
71 if (Z==79){
72 fkillBelowEnergy_Au = 10. * eV; // Kills e- tracking
73 flowEnergyLimit = 0 * eV; // Must stay at zero for killing
74 fhighEnergyLimit = 1 * GeV; // Default
75 SetLowEnergyLimit (flowEnergyLimit);
76 SetHighEnergyLimit(fhighEnergyLimit);
77 }else{
78 //continue;
79 }
80 }else{// Protection: H2O only is available
81 if(material->GetName()=="G4_WATER"){
82 flowEnergyLimit = 10. * eV;
83 fhighEnergyLimit = 1 * MeV;
84 SetLowEnergyLimit (flowEnergyLimit);
85 SetHighEnergyLimit(fhighEnergyLimit);
86 }else{
87 //continue;
88 }
89 }
90
91 if (verboseLevel > 0)
92 {
93 G4cout << "ELSEPA Elastic model is constructed for "
94 << material->GetName() << G4endl
95 << "Energy range: "
96 << flowEnergyLimit / eV << " eV - "
97 << fhighEnergyLimit / MeV << " MeV"
98 << G4endl;
99 }
100 }
101
102
103 fParticleChangeForGamma = nullptr;
104 fpMolDensity = nullptr;
105
106 fpData_Au=nullptr;
107 fpData_H2O=nullptr;
108}
109
110//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111
112G4DNAELSEPAElasticModel::~G4DNAELSEPAElasticModel()
113{
114 //std::map<G4int,G4DNACrossSectionDataSet*,
115 // std::less<G4String>>::iterator posZ;
116 //for (posZ = tableZData.begin(); posZ != tableZData.end(); ++posZ)
117 //{
118 // G4DNACrossSectionDataSet* table = posZ->second;
119 // delete table;
120 //}
121 //for (posZ = tableZData_Au.begin(); posZ != tableZData_Au.end(); ++posZ)
122 //{
123 // G4DNACrossSectionDataSet* table = posZ->second;
124 // delete table;
125 //}
126 //for (posZ = tableZData_H2O.begin(); posZ != tableZData_H2O.end(); ++posZ)
127 //{
128 // G4DNACrossSectionDataSet* table = posZ->second;
129 // delete table;
130 //}
131
132 delete fpData_Au;
133 delete fpData_H2O;
134
135 //eEdummyVecZ.clear();
136 //eCumZ.clear();
137 //fAngleDataZ.clear();
138
139 eEdummyVec_Au.clear();
140 eEdummyVec_H2O.clear();
141 eCum_Au.clear();
142 eCum_H2O.clear();
143 fAngleData_Au.clear();
144 fAngleData_H2O.clear();
145}
146
147//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
148
149void G4DNAELSEPAElasticModel::Initialise(const G4ParticleDefinition* particle,
150const G4DataVector& )
151{
152 if (verboseLevel > 3)
153 G4cout << "Calling G4DNAELSEPAElasticModel::Initialise()" << G4endl;
154
155 if (isInitialised) {return;}
156
157 if(particle->GetParticleName() != "e-")
158 {
159 G4Exception("G4DNAELSEPAElasticModel::Initialise","em0001",
160 FatalException,"Model not applicable to particle type.");
161 return;
162 }
163
164 G4ProductionCutsTable* theCoupleTable =
165 G4ProductionCutsTable::GetProductionCutsTable();
166 auto numOfCouples = (G4int)theCoupleTable->GetTableSize();
167
168 // UNIT OF TCS
169 G4double scaleFactor = 1.*cm*cm;
170
171 //tableZData.clear();
172 //tableZData_Au.clear();
173 //tableZData_H2O.clear();
174
175 fpData_Au=nullptr;
176 fpData_H2O=nullptr;
177
178 for(G4int i=0; i<numOfCouples; ++i)
179 {
180 const G4MaterialCutsCouple* couple =
181 theCoupleTable->GetMaterialCutsCouple(i);
182 const G4Material* material = couple->GetMaterial();
183 const G4ElementVector* theElementVector = material->GetElementVector();
184
185 auto nelm = (G4int)material->GetNumberOfElements();
186 if (nelm==1){// Protection: only for single element
187 G4int Z = G4lrint((*theElementVector)[0]->GetZ());
188 if (Z!=79)// Protection: only for GOLD
189 {
190 continue;
191 }
192
193 if (Z>0)
194 {
195 G4String fileZElectron("dna/sigma_elastic_e_elsepa_Z");
196 std::ostringstream oss;
197 oss.str("");
198 oss.clear(stringstream::goodbit);
199 oss << Z;
200 fileZElectron += oss.str()+"_muffintin";
201
202 //G4DNACrossSectionDataSet* tableZE =
203 // new G4DNACrossSectionDataSet
204 // (new G4LogLogInterpolation, eV,scaleFactor );
205 //tableZE->LoadData(fileZElectron);
206 ////tableZData_Au[0] = tableZE;
207 //tableZData[Z] = tableZE;
208
209 fpData_Au = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
210 eV,
211 scaleFactor );
212 fpData_Au->LoadData(fileZElectron);
213
214 std::ostringstream eFullFileNameZ;
215 const char *path = G4FindDataDir("G4LEDATA");
216 if (path == nullptr)
217 {
218 G4Exception("G4DNAELSEPAElasticModel::Initialise","em0002",
219 FatalException,"G4LEDATA environment variable not set.");
220 return;
221 }
222
223 eFullFileNameZ.str("");
224 eFullFileNameZ.clear(stringstream::goodbit);
225
226 eFullFileNameZ
227 << path
228 << "/dna/sigmadiff_cumulated_elastic_e_elsepa_Z"
229 << Z << "_muffintin.dat";
230
231 std::ifstream eDiffCrossSectionZ(eFullFileNameZ.str().c_str());
232
233 if (!eDiffCrossSectionZ)
234 {
235 G4Exception("G4DNAELSEPAElasticModel::Initialise","em0003",
236 FatalException,"Missing data file for cumulated DCS");
237 return;
238 }
239
240 //eEdummyVecZ.clear();
241 //eCumZ.clear();
242 //fAngleDataZ.clear();
243
244 eEdummyVec_Au.clear();
245 eCum_Au.clear();
246 fAngleData_Au.clear();
247
248 //eEdummyVecZ[Z].push_back(0.);
249 eEdummyVec_Au.push_back(0.);
250 do
251 {
252 G4double eDummy;
253 G4double cumDummy;
254 eDiffCrossSectionZ>>eDummy>>cumDummy;
255 //if (eDummy != eEdummyVecZ[Z].back())
256 if (eDummy != eEdummyVec_Au.back())
257 {
258
259 //eEdummyVecZ[Z].push_back(eDummy);
260 eEdummyVec_Au.push_back(eDummy);
261 //eCumZ[Z][eDummy].push_back(0.);
262 eCum_Au[eDummy].push_back(0.);
263 }
264 //eDiffCrossSectionZ>>fAngleDataZ[Z][eDummy][cumDummy];
265 eDiffCrossSectionZ>>fAngleData_Au[eDummy][cumDummy];
266 //if (cumDummy != eCumZ[Z][eDummy].back())
267 if (cumDummy != eCum_Au[eDummy].back())
268 {
269 //eCumZ[Z][eDummy].push_back(cumDummy);
270 eCum_Au[eDummy].push_back(cumDummy);
271 }
272 }while(!eDiffCrossSectionZ.eof());
273 }
274
275 }else{// Protection: H2O only is available
276 if(material->GetName()=="G4_WATER"){
277 if (LowEnergyLimit() < 10*eV)
278 {
279 G4cout<<"G4DNAELSEPAElasticModel: low energy limit increased from "
280 << LowEnergyLimit()/eV << " eV to " << 10 << " eV"
281 << G4endl;
282 SetLowEnergyLimit(10.*eV);
283 }
284
285 if (HighEnergyLimit() > 1.*MeV)
286 {
287 G4cout<<"G4DNAELSEPAElasticModel: high energy limit decreased from "
288 << HighEnergyLimit()/MeV << " MeV to " << 1. << " MeV"
289 << G4endl;
290 SetHighEnergyLimit(1.*MeV);
291 }
292
293 G4String fileZElectron("dna/sigma_elastic_e_elsepa_muffin");
294
295 //G4DNACrossSectionDataSet* tableZE =
296 // new G4DNACrossSectionDataSet(
297 // new G4LogLogInterpolation, eV,scaleFactor );
298 //tableZE->LoadData(fileZElectron);
299 ////tableZData_H2O[0] = tableZE;
300 //tableZData[0] = tableZE;
301
302 fpData_H2O = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
303 eV,
304 scaleFactor );
305 fpData_H2O->LoadData(fileZElectron);
306
307 std::ostringstream eFullFileNameZ;
308
309 const char *path = G4FindDataDir("G4LEDATA");
310 if (path == nullptr)
311 {
312 G4Exception("G4DNAELSEPAElasticModel::Initialise","em0004",
313 FatalException,"G4LEDATA environment variable not set.");
314 return;
315 }
316
317 eFullFileNameZ.str("");
318 eFullFileNameZ.clear(stringstream::goodbit);
319
320 eFullFileNameZ
321 << path
322 << "/dna/sigmadiff_cumulated_elastic_e_elsepa_muffin.dat";
323
324 std::ifstream eDiffCrossSectionZ(eFullFileNameZ.str().c_str());
325
326 if (!eDiffCrossSectionZ)
327 G4Exception("G4DNAELSEPAElasticModel::Initialise","em0005",
328 FatalException,
329 "Missing data file for cumulated DCS");
330
331 //eEdummyVecZ.clear();
332 //eCumZ.clear();
333 //fAngleDataZ.clear();
334
335 eEdummyVec_H2O.clear();
336 eCum_H2O.clear();
337 fAngleData_H2O.clear();
338
339 //eEdummyVecZ[0].push_back(0.);
340 eEdummyVec_H2O.push_back(0.);
341
342 do
343 {
344 G4double eDummy;
345 G4double cumDummy;
346 eDiffCrossSectionZ>>eDummy>>cumDummy;
347 //if (eDummy != eEdummyVecZ[0].back())
348 if (eDummy != eEdummyVec_H2O.back())
349 {
350 //eEdummyVecZ[0].push_back(eDummy);
351 eEdummyVec_H2O.push_back(eDummy);
352 //eCumZ[0][eDummy].push_back(0.);
353 eCum_H2O[eDummy].push_back(0.);
354 }
355 //eDiffCrossSectionZ>>fAngleDataZ[0][eDummy][cumDummy];
356 eDiffCrossSectionZ>>fAngleData_H2O[eDummy][cumDummy];
357 //if (cumDummy != eCumZ[0][eDummy].back()){
358 if (cumDummy != eCum_H2O[eDummy].back()){
359 //eCumZ[0][eDummy].push_back(cumDummy);
360 eCum_H2O[eDummy].push_back(cumDummy);
361 }
362 }while(!eDiffCrossSectionZ.eof());
363 }
364 }
365 if (verboseLevel > 2)
366 G4cout << "Loaded cross section files of ELSEPA Elastic model for"
367 << material->GetName() << G4endl;
368
369 if( verboseLevel>0 )
370 {
371 G4cout << "ELSEPA elastic model is initialized " << G4endl
372 << "Energy range: "
373 << LowEnergyLimit() / eV << " eV - "
374 << HighEnergyLimit()/ MeV << " MeV"
375 << G4endl;
376 }
377 } // Loop on couples
378
379
380 fParticleChangeForGamma = GetParticleChangeForGamma();
381 fpMolDensity = nullptr;
382
383 isInitialised = true;
384}
385
386//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
387
388G4double G4DNAELSEPAElasticModel::CrossSectionPerVolume
389(const G4Material* material,
390 const G4ParticleDefinition* particle,
391 G4double ekin,
392 G4double,
393 G4double)
394{
395
396 if (verboseLevel > 3)
397 {
398 G4cout <<
399 "Calling CrossSectionPerVolume() of G4DNAELSEPAElasticModel"
400 << G4endl;
401 }
402
403 G4double atomicNDensity=0.0;
404 G4double sigma=0;
405
406 const G4ElementVector* theElementVector = material->GetElementVector();
407 std::size_t nelm = material->GetNumberOfElements();
408 if (nelm==1) // Protection: only for single element
409 {
410 // Protection: only for GOLD
411 if (material->GetZ()!=79) return 0.0;
412
413 G4int Z = G4lrint((*theElementVector)[0]->GetZ());
414
415 const G4String& particleName = particle->GetParticleName();
416 atomicNDensity = material->GetAtomicNumDensityVector()[0];
417 if(atomicNDensity!= 0.0)
418 {
419 if (ekin < fhighEnergyLimit)
420 {
421 if (ekin < fkillBelowEnergy_Au) return DBL_MAX;
422
423 //std::map< G4int,G4DNACrossSectionDataSet*,
424 // std::less<G4String> >::iterator pos;
425 ////pos = tableZData_Au.find(0);
426 //pos = tableZData.find(Z);
427 //
428 ////if (pos != tableZData_Au.end())
429 //if (pos != tableZData.end())
430 //{
431 // G4DNACrossSectionDataSet* table = pos->second;
432 // if (table != 0)
433 // {
434 // // XS takes its 10 eV value below 10 eV for GOLD
435 // if (ekin < 10*eV) sigma = table->FindValue(10*eV);
436 // else sigma = table->FindValue(ekin);
437 // }
438 //}
439 //else
440 //{
441 // G4Exception("G4DNAELSEPAElasticModel::ComputeCrossSectionPerVolume",
442 // "em0006",FatalException,"Model not applicable to particle type.");
443 //}
444
445 if (ekin < 10*eV) sigma = fpData_Au->FindValue(10*eV);
446 else sigma = fpData_Au->FindValue(ekin);
447 }
448 }
449 if (verboseLevel > 2)
450 {
451 G4cout << "__________________________________" << G4endl;
452 G4cout << "=== G4DNAELSEPAElasticModel - XS INFO START" << G4endl;
453 G4cout << "=== Material is made of one element with Z =" << Z << G4endl;
454 G4cout << "=== Kinetic energy(eV)=" << ekin/eV << " particle : "
455 << particleName << G4endl;
456 G4cout << "=== Cross section per atom for Z="<<Z<<" is (cm^2)"
457 << sigma/cm/cm << G4endl;
458 G4cout << "=== Cross section per atom for Z="<<Z<<" is (cm^-1)="
459 << sigma*atomicNDensity/(1./cm) << G4endl;
460 G4cout << "=== G4DNAELSEPAElasticModel - XS INFO END" << G4endl;
461 }
462 }
463 else
464 {
465 fpMolDensity =
466 G4DNAMolecularMaterial::Instance()->
467 GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
468 atomicNDensity = (*fpMolDensity)[material->GetIndex()];
469 if(atomicNDensity!= 0.0)
470 {
471 if (ekin < HighEnergyLimit() && ekin >= LowEnergyLimit())
472 {
473 //std::map< G4int,G4DNACrossSectionDataSet*,
474 //std::less<G4String> >::iterator pos;
475 ////pos = tableZData_H2O.find(0); // the data is stored as Z=0
476 //pos = tableZData.find(0); // the data is stored as Z=0
477 ////SI : XS must not be zero
478 //// otherwise sampling of secondaries method ignored
479 ////if (pos != tableZData_H2O.end())
480 //if (pos != tableZData.end())
481 //{
482 // G4DNACrossSectionDataSet* table = pos->second;
483 // if (table != 0)
484 // {
485 // sigma = table->FindValue(ekin);
486 // }
487 //}
488
489 sigma = fpData_H2O->FindValue(ekin);
490 }
491 }
492 if (verboseLevel > 2)
493 {
494 G4cout << "__________________________________" << G4endl;
495 G4cout << "=== G4DNAELSEPAElasticModel - XS INFO START" << G4endl;
496 G4cout << "=== Kinetic energy(eV)=" << ekin/eV
497 << " particle : " << particle->GetParticleName() << G4endl;
498 G4cout << "=== Cross section per water molecule (cm^2)="
499 << sigma/cm/cm << G4endl;
500 G4cout << "=== Cross section per water molecule (cm^-1)="
501 << sigma*atomicNDensity/(1./cm) << G4endl;
502 G4cout << "=== G4DNAELSEPAElasticModel - XS INFO END" << G4endl;
503 }
504 }
505
506 return sigma*atomicNDensity;
507}
508
509//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
510
511void G4DNAELSEPAElasticModel::SampleSecondaries(
512 std::vector<G4DynamicParticle*>*,
513 const G4MaterialCutsCouple* couple,
514 const G4DynamicParticle* aDynamicElectron,
515 G4double,
516 G4double)
517{
518
519 if (verboseLevel > 3){
520 G4cout <<
521 "Calling SampleSecondaries() of G4DNAELSEPAElasticModel"
522 << G4endl;
523 }
524
525 G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
526
527 const G4Material* material = couple->GetMaterial();
528 const G4ElementVector* theElementVector = material->GetElementVector();
529 std::size_t nelm = material->GetNumberOfElements();
530 if (nelm==1) // Protection: only for single element
531 {
532 G4int Z = G4lrint((*theElementVector)[0]->GetZ());
533 if (Z!=79) return;
534 if (electronEnergy0 < fkillBelowEnergy_Au)
535 {
536 fParticleChangeForGamma->SetProposedKineticEnergy(0.);
537 fParticleChangeForGamma->ProposeMomentumDirection(G4ThreeVector(0,0,0));
538 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
539 fParticleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0);
540 return;
541 }
542
543 if(electronEnergy0>= fkillBelowEnergy_Au && electronEnergy0 < fhighEnergyLimit)
544 {
545 G4double cosTheta = 0;
546 if (electronEnergy0>=10*eV)
547 {
548 cosTheta = RandomizeCosTheta(Z,electronEnergy0);
549 }
550 else
551 {
552 cosTheta = RandomizeCosTheta(Z,10*eV);
553 }
554
555 G4double phi = 2. * CLHEP::pi * G4UniformRand();
556
557 G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
558 G4ThreeVector xVers = zVers.orthogonal();
559 G4ThreeVector yVers = zVers.cross(xVers);
560
561 G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
562 G4double yDir = xDir;
563 xDir *= std::cos(phi);
564 yDir *= std::sin(phi);
565
566 G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
567 fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit());
568 fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
569
570 }
571 }
572 else
573 {
574 if(material->GetName()=="G4_WATER")
575 {
576 //The data for water is stored as Z=0
577 G4double cosTheta = RandomizeCosTheta(0,electronEnergy0);
578
579 G4double phi = 2. * pi * G4UniformRand();
580
581 G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
582 G4ThreeVector xVers = zVers.orthogonal();
583 G4ThreeVector yVers = zVers.cross(xVers);
584
585 G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
586 G4double yDir = xDir;
587 xDir *= std::cos(phi);
588 yDir *= std::sin(phi);
589
590 G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
591 fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit());
592 fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
593 }
594 }
595}
596
597//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
598
599G4double G4DNAELSEPAElasticModel::Theta(G4int Z,
600 G4ParticleDefinition * particleDefinition,
601 G4double k,
602 G4double integrDiff)
603{
604
605 G4double theta = 0.;
606 G4double valueE1 = 0.;
607 G4double valueE2 = 0.;
608 G4double valuecum21 = 0.;
609 G4double valuecum22 = 0.;
610 G4double valuecum12 = 0.;
611 G4double valuecum11 = 0.;
612 G4double a11 = 0.;
613 G4double a12 = 0.;
614 G4double a21 = 0.;
615 G4double a22 = 0.;
616
617 if (particleDefinition == G4Electron::ElectronDefinition())
618 {
619 //std::vector<G4double>::iterator e2
620 // = std::upper_bound(eEdummyVecZ[Z].begin(),
621 // eEdummyVecZ[Z].end(), k);
622 std::vector<G4double>::iterator e2;
623 if(Z==0){
624 e2 = std::upper_bound(eEdummyVec_H2O.begin(),
625 eEdummyVec_H2O.end(), k);
626 }else if (Z==79){
627 e2 = std::upper_bound(eEdummyVec_Au.begin(),
628 eEdummyVec_Au.end(), k);
629 }
630
631 auto e1 = e2 - 1;
632
633 //std::vector<G4double>::iterator cum12
634 // = std::upper_bound(eCumZ[Z][(*e1)].begin(),
635 // eCumZ[Z][(*e1)].end(),integrDiff);
636 std::vector<G4double>::iterator cum12;
637 if(Z==0){
638 cum12 = std::upper_bound(eCum_H2O[(*e1)].begin(),
639 eCum_H2O[(*e1)].end(),integrDiff);
640 }else if (Z==79){
641 cum12 = std::upper_bound(eCum_Au[(*e1)].begin(),
642 eCum_Au[(*e1)].end(),integrDiff);
643 }
644
645 auto cum11 = cum12 - 1;
646
647 //std::vector<G4double>::iterator cum22
648 // = std::upper_bound(eCumZ[Z][(*e2)].begin(),
649 // eCumZ[Z][(*e2)].end(),integrDiff);
650 std::vector<G4double>::iterator cum22;
651 if(Z==0){
652 cum22 = std::upper_bound(eCum_H2O[(*e2)].begin(),
653 eCum_H2O[(*e2)].end(),integrDiff);
654 }else if(Z==79){
655 cum22 = std::upper_bound(eCum_Au[(*e2)].begin(),
656 eCum_Au[(*e2)].end(),integrDiff);
657 }
658
659 auto cum21 = cum22 - 1;
660
661 valueE1 = *e1;
662 valueE2 = *e2;
663 valuecum11 = *cum11;
664 valuecum12 = *cum12;
665 valuecum21 = *cum21;
666 valuecum22 = *cum22;
667
668
669 //a11 = fAngleDataZ[Z][valueE1][valuecum11];
670 //a12 = fAngleDataZ[Z][valueE1][valuecum12];
671 //a21 = fAngleDataZ[Z][valueE2][valuecum21];
672 //a22 = fAngleDataZ[Z][valueE2][valuecum22];
673 if(Z==0){
674 a11 = fAngleData_H2O[valueE1][valuecum11];
675 a12 = fAngleData_H2O[valueE1][valuecum12];
676 a21 = fAngleData_H2O[valueE2][valuecum21];
677 a22 = fAngleData_H2O[valueE2][valuecum22];
678 }else if (Z==79){
679 a11 = fAngleData_Au[valueE1][valuecum11];
680 a12 = fAngleData_Au[valueE1][valuecum12];
681 a21 = fAngleData_Au[valueE2][valuecum21];
682 a22 = fAngleData_Au[valueE2][valuecum22];
683 }
684
685 }
686
687 if (a11 == 0 && a12 == 0 && a21 == 0 && a22 == 0) return (0.);
688
689 theta = QuadInterpolator(valuecum11, valuecum12, valuecum21, valuecum22,
690 a11, a12,a21, a22, valueE1, valueE2, k, integrDiff);
691 return theta;
692}
693
694//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
695//
696G4double G4DNAELSEPAElasticModel::LogLinInterpolate(G4double e1,
697G4double e2,
698G4double e,
699G4double xs1,
700G4double xs2)
701{
702 G4double value=0.;
703 if(e1!=0){
704 G4double a = std::log10(e) - std::log10(e1);
705 G4double b = std::log10(e2) - std::log10(e);
706 value = xs1 + a/(a+b)*(xs2-xs1);
707 }
708 else{
709 G4double d1 = xs1;
710 G4double d2 = xs2;
711 value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
712 }
713
714 return value;
715}
716
717//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
718
719G4double G4DNAELSEPAElasticModel::LinLogInterpolate(G4double e1,
720G4double e2,
721G4double e,
722G4double xs1,
723G4double xs2)
724{
725 G4double d1 = std::log10(xs1);
726 G4double d2 = std::log10(xs2);
727 G4double value = std::pow(10,(d1 + (d2 - d1) * (e - e1) / (e2 - e1)));
728 return value;
729}
730
731//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
732
733G4double G4DNAELSEPAElasticModel::LinLinInterpolate(G4double e1,
734G4double e2,
735G4double e,
736G4double xs1,
737G4double xs2)
738{
739 G4double d1 = xs1;
740 G4double d2 = xs2;
741 G4double value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
742 return value;
743}
744
745//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
746
747G4double G4DNAELSEPAElasticModel::LogLogInterpolate(G4double e1,
748G4double e2,
749G4double e,
750G4double xs1,
751G4double xs2)
752{
753 G4double a = (std::log10(xs2) - std::log10(xs1))
754 / (std::log10(e2) - std::log10(e1));
755 G4double b = std::log10(xs2) - a * std::log10(e2);
756 G4double sigma = a * std::log10(e) + b;
757 G4double value = (std::pow(10., sigma));
758 return value;
759}
760
761//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
762
763G4double G4DNAELSEPAElasticModel::QuadInterpolator(
764G4double cum11,
765G4double cum12,
766G4double cum21,
767G4double cum22,
768G4double a11,
769G4double a12,
770G4double a21,
771G4double a22,
772G4double e1,
773G4double e2,
774G4double t,
775G4double cum)
776{
777 G4double value=0;
778 G4double interpolatedvalue1=0;
779 G4double interpolatedvalue2=0;
780
781 if(cum11!=0){
782 interpolatedvalue1 = LinLogInterpolate(cum11, cum12, cum, a11, a12);
783 }
784 else{
785 interpolatedvalue1 = LinLinInterpolate(cum11, cum12, cum, a11, a12);
786 }
787 if(cum21!=0){
788 interpolatedvalue2 = LinLogInterpolate(cum21, cum22, cum, a21, a22);
789 }
790 else{
791 interpolatedvalue2 = LinLinInterpolate(cum21, cum22, cum, a21, a22);
792 }
793
794 value = LogLinInterpolate(e1,e2,t,interpolatedvalue1,interpolatedvalue2);
795
796 return value;
797}
798
799//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
800
801G4double G4DNAELSEPAElasticModel::RandomizeCosTheta(G4int Z, G4double k)
802{
803
804 G4double integrdiff = 0.;
805 G4double uniformRand = G4UniformRand();
806 integrdiff = uniformRand;
807
808 G4double theta = 0.;
809 G4double cosTheta = 0.;
810 theta = Theta(Z, G4Electron::ElectronDefinition(), k / eV, integrdiff);
811
812 cosTheta = std::cos(theta * CLHEP::pi / 180.);
813
814 return cosTheta;
815}
816
817//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
818
819void G4DNAELSEPAElasticModel::SetKillBelowThreshold(G4double threshold)
820{
821 fkillBelowEnergy_Au = threshold;
822
823 if (threshold < 10 * eV)
824 {
825 G4cout<< "*** WARNING : the G4DNAELSEPAElasticModel model is not "
826 "defined below 10 eV !" << G4endl;
827 }
828}
G4ElementVector
std::vector< const G4Element * > G4ElementVector
Definition G4ElementVector.hh:33
G4FindDataDir
const char * G4FindDataDir(const char *)
Definition G4FindDataDir.cc:94
FatalException
@ FatalException
Definition G4ExceptionSeverity.hh:67
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition G4Exception.cc:59
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition G4ThreeVector.hh:36
fStopAndKill
@ fStopAndKill
Definition G4TrackStatus.hh:46
G4double
double G4double
Definition G4Types.hh:83
G4int
int G4int
Definition G4Types.hh:85
G4endl
#define G4endl
Definition G4ios.hh:67
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition Randomize.hh:52
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
CLHEP::Hep3Vector::orthogonal
Hep3Vector orthogonal() const
CLHEP::Hep3Vector::cross
Hep3Vector cross(const Hep3Vector &) const
G4DNACrossSectionDataSet::LoadData
G4bool LoadData(const G4String &argFileName) override
Definition G4DNACrossSectionDataSet.cc:83
G4DNACrossSectionDataSet::FindValue
G4double FindValue(G4double e, G4int componentId=0) const override
Definition G4DNACrossSectionDataSet.cc:483
G4DNAELSEPAElasticModel::G4DNAELSEPAElasticModel
G4DNAELSEPAElasticModel(const G4ParticleDefinition *particle=nullptr, const G4String &nam="DNAELSEPAElasticModel")
Definition G4DNAELSEPAElasticModel.cc:48
G4DNAELSEPAElasticModel::CrossSectionPerVolume
G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *particle, G4double ekin, G4double emin, G4double emax) override
Definition G4DNAELSEPAElasticModel.cc:389
G4DNAELSEPAElasticModel::Initialise
void Initialise(const G4ParticleDefinition *particle, const G4DataVector &) override
Definition G4DNAELSEPAElasticModel.cc:149
G4DNAELSEPAElasticModel::SampleSecondaries
void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
Definition G4DNAELSEPAElasticModel.cc:511
G4DNAELSEPAElasticModel::SetKillBelowThreshold
void SetKillBelowThreshold(G4double threshold)
Definition G4DNAELSEPAElasticModel.cc:819
G4DNAELSEPAElasticModel::fParticleChangeForGamma
G4ParticleChangeForGamma * fParticleChangeForGamma
Definition G4DNAELSEPAElasticModel.hh:86
G4DNAMolecularMaterial::Instance
static G4DNAMolecularMaterial * Instance()
Definition G4DNAMolecularMaterial.cc:80
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4Electron::ElectronDefinition
static G4Electron * ElectronDefinition()
Definition G4Electron.cc:86
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition G4MaterialCutsCouple.hh:166
G4Material::GetElementVector
const G4ElementVector * GetElementVector() const
Definition G4Material.hh:183
G4Material::GetZ
G4double GetZ() const
Definition G4Material.cc:700
G4Material::GetAtomicNumDensityVector
const G4double * GetAtomicNumDensityVector() const
Definition G4Material.hh:202
G4Material::GetIndex
std::size_t GetIndex() const
Definition G4Material.hh:239
G4Material::GetNumberOfElements
std::size_t GetNumberOfElements() const
Definition G4Material.hh:180
G4Material::GetName
const G4String & GetName() const
Definition G4Material.hh:171
G4Material::GetMaterial
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition G4Material.cc:653
G4ParticleChangeForGamma::SetProposedKineticEnergy
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Definition G4ParticleChangeForGamma.hh:106
G4ParticleChangeForGamma::ProposeMomentumDirection
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)
Definition G4ParticleChangeForGamma.hh:120
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition G4ParticleDefinition.hh:96
G4ProductionCutsTable::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
Definition G4ProductionCutsTable.hh:256
G4ProductionCutsTable::GetTableSize
std::size_t GetTableSize() const
Definition G4ProductionCutsTable.hh:249
G4ProductionCutsTable::GetProductionCutsTable
static G4ProductionCutsTable * GetProductionCutsTable()
Definition G4ProductionCutsTable.cc:58
G4VEmModel::SetHighEnergyLimit
void SetHighEnergyLimit(G4double)
Definition G4VEmModel.hh:738
G4VEmModel::GetParticleChangeForGamma
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition G4VEmModel.cc:120
G4VEmModel::LowEnergyLimit
G4double LowEnergyLimit() const
Definition G4VEmModel.hh:640
G4VEmModel::HighEnergyLimit
G4double HighEnergyLimit() const
Definition G4VEmModel.hh:633
G4VEmModel::SetLowEnergyLimit
void SetLowEnergyLimit(G4double)
Definition G4VEmModel.hh:745
G4VParticleChange::ProposeTrackStatus
void ProposeTrackStatus(G4TrackStatus status)
G4VParticleChange::ProposeLocalEnergyDeposit
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4lrint
int G4lrint(double ad)
Definition templates.hh:134
DBL_MAX
#define DBL_MAX
Definition templates.hh:62