42 fKillBelowEnergy = 10*eV;
71 G4cout <<
"Calling G4DNAPTBElasticModel::Initialise()" <<
G4endl;
81 if(particle == electronDef)
87 "dna/sigma_elastic_e-_PTB_THF",
88 "dna/sigmadiff_cumulated_elastic_e-_PTB_THF",
95 "dna/sigma_elastic_e-_PTB_PY",
96 "dna/sigmadiff_cumulated_elastic_e-_PTB_PY",
103 "dna/sigma_elastic_e-_PTB_PU",
104 "dna/sigmadiff_cumulated_elastic_e-_PTB_PU",
111 "dna/sigma_elastic_e-_PTB_TMP",
112 "dna/sigmadiff_cumulated_elastic_e-_PTB_TMP",
119 "dna/sigma_elastic_e_champion",
120 "dna/sigmadiff_cumulated_elastic_e_champion",
129 "dna/sigma_elastic_e-_PTB_THF",
130 "dna/sigmadiff_cumulated_elastic_e-_PTB_THF",
137 "dna/sigma_elastic_e-_PTB_PY",
138 "dna/sigmadiff_cumulated_elastic_e-_PTB_PY",
145 "dna/sigma_elastic_e-_PTB_PY",
146 "dna/sigmadiff_cumulated_elastic_e-_PTB_PY",
153 "dna/sigma_elastic_e-_PTB_PU",
154 "dna/sigmadiff_cumulated_elastic_e-_PTB_PU",
161 "dna/sigma_elastic_e-_PTB_PU",
162 "dna/sigmadiff_cumulated_elastic_e-_PTB_PU",
169 "dna/sigma_elastic_e-_PTB_TMP",
170 "dna/sigmadiff_cumulated_elastic_e-_PTB_TMP",
186 if (verboseLevel > 2)
187 G4cout <<
"Loaded cross section files for PTB Elastic model" <<
G4endl;
191 G4cout <<
"PTB Elastic model is initialized " <<
G4endl;
197void G4DNAPTBElasticModel::ReadDiffCSFile(
const G4String& materialName,
206 char *path = std::getenv(
"G4LEDATA");
210 G4Exception(
"G4DNAPTBElasticModel::ReadAllDiffCSFiles",
"em0006",
216 std::ostringstream fullFileName;
217 fullFileName << path <<
"/"<< file<<
".dat";
220 std::ifstream diffCrossSection (fullFileName.str().c_str());
222 std::stringstream endPath;
223 if (!diffCrossSection)
225 endPath <<
"Missing data file: "<<file;
226 G4Exception(
"G4DNAPTBElasticModel::Initialise",
"em0003",
230 tValuesVec[materialName][particleName].push_back(0.);
235 while(std::getline(diffCrossSection, line))
239 std::istringstream testIss(line);
249 else if(line.empty())
258 std::istringstream iss(line);
276 if (tDummy != tValuesVec[materialName][particleName].back())
279 tValuesVec[materialName][particleName].push_back(tDummy);
282 eValuesVect[materialName][particleName][tDummy].push_back(0.);
286 iss>>diffCrossSectionData[materialName][particleName][tDummy][eDummy];
289 if (eDummy != eValuesVect[materialName][particleName][tDummy].back()) eValuesVect[materialName][particleName][tDummy].push_back(eDummy);
302 if (verboseLevel > 3)
303 G4cout <<
"Calling CrossSectionPerVolume() of G4DNAPTBElasticModel" <<
G4endl;
309 fKillBelowEnergy =
GetLowELimit(materialName, particleName);
322 if (ekin < fKillBelowEnergy)
return DBL_MAX;
328 sigma = (*tableData)[materialName][particleName]->FindValue(ekin);
331 if (verboseLevel > 2)
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;
354 if (verboseLevel > 3)
355 G4cout <<
"Calling SampleSecondaries() of G4DNAPTBElasticModel" <<
G4endl;
362 fKillBelowEnergy =
GetLowELimit(materialName, particleName);
365 if (electronEnergy0 < fKillBelowEnergy)
372 else if (electronEnergy0>= fKillBelowEnergy && electronEnergy0 <
GetHighELimit(materialName, particleName) )
375 G4double cosTheta = RandomizeCosTheta(electronEnergy0, materialName);
384 G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
386 xDir *= std::cos(phi);
387 yDir *= std::sin(phi);
390 G4ThreeVector zPrikeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
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;
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;
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;
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];
442 if (xs11==0 && xs12==0 && xs21==0 && xs22==0)
return (0.);
444 theta = QuadInterpolator ( valueE11, valueE12,
464 G4double value = std::exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
478 G4double value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
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));
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);
533 integrdiff = uniformRand;
539 cosTheta= std::cos(theta*pi/180);
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4GLOB_DLL std::ostream G4cout
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()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
const G4String & GetParticleName() const
TableMapData * GetTableData()
GetTableData.
void SetHighELimit(const G4String &material, const G4String &particle, G4double lim)
SetHighEnergyLimit.
std::map< G4String, std::map< G4String, G4DNACrossSectionDataSet *, std::less< G4String > > > TableMapData
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....
void SetLowELimit(const G4String &material, const G4String &particle, G4double lim)
SetLowEnergyLimit.
G4double GetHighELimit(const G4String &material, const G4String &particle)
GetHighEnergyLimit.
void LoadCrossSectionData(const G4String &particleName)
LoadCrossSectionData Method to loop on all the registered materials in the model and load the corresp...
G4double GetLowELimit(const G4String &material, const G4String &particle)
GetLowEnergyLimit.
void ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)