73 G4cout <<
"CPA100 Elastic model is constructed " <<
G4endl
82 fpMolWaterDensity = 0;
95 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
96 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
115 if (verboseLevel > 3)
116 G4cout <<
"Calling G4DNACPA100ElasticModel::Initialise()" <<
G4endl;
121 G4Exception(
"*** WARNING: the G4DNACPA100ElasticModel is "
122 "not intented to be used with another particle than the electron",
130 G4cout <<
"G4DNACPA100ElasticModel: low energy limit increased from " <<
137 G4cout <<
"G4DNACPA100ElasticModel: high energy limit decreased from " <<
147 G4String fileElectron(
"dna/sigma_elastic_e_cpa100");
158 tableFile[electron] = fileElectron;
172 tableData[electron] = tableE;
176 char *path = getenv(
"G4LEDATA");
180 G4Exception(
"G4DNACPA100ElasticModel::Initialise",
"em0006",
185 std::ostringstream eFullFileName;
187 eFullFileName << path
188 <<
"/dna/sigmadiff_cumulated_elastic_e_cpa100.dat";
190 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
192 if (!eDiffCrossSection)
193 G4Exception(
"G4DNACPA100ElasticModel::Initialise",
"em0003",
195 "Missing data file:/dna/sigmadiff_cumulated_elastic_e_cpa100.dat");
202 eDiffCrossSectionData.clear();
206 eTdummyVec.push_back(0.);
208 while(!eDiffCrossSection.eof())
212 eDiffCrossSection>>tDummy>>eDummy;
216 if (tDummy != eTdummyVec.back())
218 eTdummyVec.push_back(tDummy);
219 eVecm[tDummy].push_back(0.);
222 eDiffCrossSection>>eDiffCrossSectionData[tDummy][eDummy];
224 if (eDummy != eVecm[tDummy].back()) eVecm[tDummy].push_back(eDummy);
230 if (verboseLevel > 2)
231 G4cout <<
"Loaded cross section files for CPA100 Elastic model" <<
G4endl;
237 G4cout <<
"CPA100 Elastic model is initialized " <<
G4endl
249 if (isInitialised) {
return; }
251 isInitialised =
true;
266 if (verboseLevel > 3)
268 "Calling CrossSectionPerVolume() of G4DNACPA100ElasticModel" <<
G4endl;
284 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
285 pos = tableData.find(particleName);
287 if (pos != tableData.end())
325 G4Exception(
"G4DNACPA100ElasticModel::ComputeCrossSectionPerVolume",
332 if (verboseLevel > 2)
334 G4cout <<
"__________________________________" <<
G4endl;
335 G4cout <<
"G4DNACPA100ElasticModel - XS INFO START" <<
G4endl;
336 G4cout <<
"Kinetic energy(eV)=" << ekin/eV <<
" particle : " << particleName <<
G4endl;
337 G4cout <<
"Cross section per water molecule (cm^2)=" << sigma/cm/cm <<
G4endl;
338 G4cout <<
"Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) <<
G4endl;
341 G4cout <<
"G4DNACPA100ElasticModel - XS INFO END" <<
G4endl;
345 return sigma*waterDensity;
357 if (verboseLevel > 3)
358 G4cout <<
"Calling SampleSecondaries() of G4DNACPA100ElasticModel" <<
G4endl;
363 G4double cosTheta = RandomizeCosTheta(electronEnergy0);
377 G4double CT1, ST1, CF1, SF1, CT2, ST2, CF2, SF2;
378 G4double sinTheta = std::sqrt (1-cosTheta*cosTheta);
390 ST1=std::sqrt(1.-CT1*CT1);
392 if (ST1!=0) CF1 = zVers.
x()/ST1;
else CF1 = std::cos(2. * pi *
G4UniformRand());
393 if (ST1!=0) SF1 = zVers.
y()/ST1;
else SF1 = std::sqrt(1.-CF1*CF1);
402 A3 = sinTheta*std::cos(phi);
403 A4 = A3*CT1 + ST1*cosTheta;
404 A5 = sinTheta * std::sin(phi);
405 A2 = A4 * SF1 + A5 * CF1;
406 A1 = A4 * CF1 - A5 * SF1;
408 CT2 = CT1*cosTheta - ST1*A3;
409 ST2 = std::sqrt(1.-CT2*CT2);
411 if (ST2==0) ST2=1E-6;
439 (electronEnergy0-1.214E-4*(1.-cosTheta)*electronEnergy0);
451G4double G4DNACPA100ElasticModel::Theta
468 if (k==eTdummyVec.back()) k=k*(1.-1e-12);
471 std::vector<G4double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
472 std::vector<G4double>::iterator t1 = t2-1;
474 std::vector<G4double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(),
476 std::vector<G4double>::iterator e11 = e12-1;
478 std::vector<G4double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(),
480 std::vector<G4double>::iterator e21 = e22-1;
489 xs11 = eDiffCrossSectionData[valueT1][valueE11];
490 xs12 = eDiffCrossSectionData[valueT1][valueE12];
491 xs21 = eDiffCrossSectionData[valueT2][valueE21];
492 xs22 = eDiffCrossSectionData[valueT2][valueE22];
497 if (xs11==0 && xs12==0 && xs21==0 && xs22==0)
return (0.);
499 theta = QuadInterpolator(
523 G4double value = std::exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
537 G4double value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
549 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
550 G4double b = std::log10(xs2) - a*std::log10(e2);
551 G4double sigma = a*std::log10(e) + b;
552 G4double value = (std::pow(10.,sigma));
580 G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12);
581 G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22);
582 G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
594 integrdiff = uniformRand;
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4GLOB_DLL std::ostream G4cout
virtual ~G4DNACPA100ElasticModel()
G4ParticleChangeForGamma * fParticleChangeForGamma
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4DNACPA100ElasticModel(const G4ParticleDefinition *p=0, const G4String &nam="DNACPA100ElasticModel")
virtual G4double FindValue(G4double e, G4int componentId=0) const
virtual G4bool LoadData(const G4String &argFileName)
const std::vector< G4double > * GetNumMolPerVolTableFor(const G4Material *) const
Retrieve a table of molecular densities (number of molecules per unit volume) in the G4 unit system f...
static G4DNAMolecularMaterial * Instance()
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
static G4Electron * ElectronDefinition()
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
const G4String & GetParticleName() const
void SetHighEnergyLimit(G4double)
G4ParticleChangeForGamma * GetParticleChangeForGamma()
G4double LowEnergyLimit() const
G4double HighEnergyLimit() const
void SetLowEnergyLimit(G4double)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)