00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include "G4DNAChampionElasticModel.hh"
00030 #include "G4PhysicalConstants.hh"
00031 #include "G4SystemOfUnits.hh"
00032 #include "G4DNAMolecularMaterial.hh"
00033
00034
00035
00036 using namespace std;
00037
00038
00039
00040 G4DNAChampionElasticModel::G4DNAChampionElasticModel(const G4ParticleDefinition*,
00041 const G4String& nam)
00042 :G4VEmModel(nam),isInitialised(false)
00043 {
00044
00045
00046 killBelowEnergy = 7.4*eV;
00047 lowEnergyLimit = 0 * eV;
00048 highEnergyLimit = 1. * MeV;
00049 SetLowEnergyLimit(lowEnergyLimit);
00050 SetHighEnergyLimit(highEnergyLimit);
00051
00052 verboseLevel= 0;
00053
00054
00055
00056
00057
00058
00059
00060 if( verboseLevel>0 )
00061 {
00062 G4cout << "Champion Elastic model is constructed " << G4endl
00063 << "Energy range: "
00064 << lowEnergyLimit / eV << " eV - "
00065 << highEnergyLimit / MeV << " MeV"
00066 << G4endl;
00067 }
00068 fParticleChangeForGamma = 0;
00069 fpMolWaterDensity = 0;
00070 }
00071
00072
00073
00074 G4DNAChampionElasticModel::~G4DNAChampionElasticModel()
00075 {
00076
00077
00078 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
00079 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
00080 {
00081 G4DNACrossSectionDataSet* table = pos->second;
00082 delete table;
00083 }
00084
00085
00086
00087 eVecm.clear();
00088
00089 }
00090
00091
00092
00093 void G4DNAChampionElasticModel::Initialise(const G4ParticleDefinition* ,
00094 const G4DataVector& )
00095 {
00096
00097 if (verboseLevel > 3)
00098 G4cout << "Calling G4DNAChampionElasticModel::Initialise()" << G4endl;
00099
00100
00101
00102 if (LowEnergyLimit() < lowEnergyLimit)
00103 {
00104 G4cout << "G4DNAChampionElasticModel: low energy limit increased from " <<
00105 LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
00106 SetLowEnergyLimit(lowEnergyLimit);
00107 }
00108
00109 if (HighEnergyLimit() > highEnergyLimit)
00110 {
00111 G4cout << "G4DNAChampionElasticModel: high energy limit decreased from " <<
00112 HighEnergyLimit()/MeV << " MeV to " << highEnergyLimit/MeV << " MeV" << G4endl;
00113 SetHighEnergyLimit(highEnergyLimit);
00114 }
00115
00116
00117
00118 G4double scaleFactor = 1e-16*cm*cm;
00119
00120 G4String fileElectron("dna/sigma_elastic_e_champion");
00121
00122 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
00123 G4String electron;
00124
00125
00126
00127
00128
00129 electron = electronDef->GetParticleName();
00130
00131 tableFile[electron] = fileElectron;
00132
00133 G4DNACrossSectionDataSet* tableE = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
00134 tableE->LoadData(fileElectron);
00135 tableData[electron] = tableE;
00136
00137
00138
00139 char *path = getenv("G4LEDATA");
00140
00141 if (!path)
00142 {
00143 G4Exception("G4ChampionElasticModel::Initialise","em0006",
00144 FatalException,"G4LEDATA environment variable not set.");
00145 return;
00146 }
00147
00148 std::ostringstream eFullFileName;
00149 eFullFileName << path << "/dna/sigmadiff_cumulated_elastic_e_champion.dat";
00150 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
00151
00152 if (!eDiffCrossSection)
00153 G4Exception("G4DNAChampionElasticModel::Initialise","em0003",
00154 FatalException,"Missing data file:/dna/sigmadiff_cumulated_elastic_e_champion.dat");
00155
00156 eTdummyVec.push_back(0.);
00157
00158 while(!eDiffCrossSection.eof())
00159 {
00160 double tDummy;
00161 double eDummy;
00162 eDiffCrossSection>>tDummy>>eDummy;
00163
00164
00165
00166 if (tDummy != eTdummyVec.back())
00167 {
00168 eTdummyVec.push_back(tDummy);
00169 eVecm[tDummy].push_back(0.);
00170 }
00171
00172 eDiffCrossSection>>eDiffCrossSectionData[tDummy][eDummy];
00173
00174 if (eDummy != eVecm[tDummy].back()) eVecm[tDummy].push_back(eDummy);
00175
00176 }
00177
00178
00179
00180 if (verboseLevel > 2)
00181 G4cout << "Loaded cross section files for Champion Elastic model" << G4endl;
00182
00183 if( verboseLevel>0 )
00184 {
00185 G4cout << "Champion Elastic model is initialized " << G4endl
00186 << "Energy range: "
00187 << LowEnergyLimit() / eV << " eV - "
00188 << HighEnergyLimit() / MeV << " MeV"
00189 << G4endl;
00190 }
00191
00192
00193 fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
00194
00195 if (isInitialised) { return; }
00196 fParticleChangeForGamma = GetParticleChangeForGamma();
00197 isInitialised = true;
00198
00199 }
00200
00201
00202
00203 G4double G4DNAChampionElasticModel::CrossSectionPerVolume(const G4Material* material,
00204 const G4ParticleDefinition* p,
00205 G4double ekin,
00206 G4double,
00207 G4double)
00208 {
00209 if (verboseLevel > 3)
00210 G4cout << "Calling CrossSectionPerVolume() of G4DNAChampionElasticModel" << G4endl;
00211
00212
00213
00214 G4double sigma=0;
00215
00216 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
00217
00218 if(waterDensity!= 0.0)
00219
00220 {
00221 const G4String& particleName = p->GetParticleName();
00222
00223 if (ekin < highEnergyLimit)
00224 {
00225
00226 if (ekin < killBelowEnergy) return DBL_MAX;
00227
00228
00229 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
00230 pos = tableData.find(particleName);
00231
00232 if (pos != tableData.end())
00233 {
00234 G4DNACrossSectionDataSet* table = pos->second;
00235 if (table != 0)
00236 {
00237 sigma = table->FindValue(ekin);
00238 }
00239 }
00240 else
00241 {
00242 G4Exception("G4DNAChampionElasticModel::ComputeCrossSectionPerVolume","em0002",
00243 FatalException,"Model not applicable to particle type.");
00244 }
00245 }
00246
00247 if (verboseLevel > 2)
00248 {
00249 G4cout << "__________________________________" << G4endl;
00250 G4cout << "°°° G4DNAChampionElasticModel - XS INFO START" << G4endl;
00251 G4cout << "°°° Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
00252 G4cout << "°°° Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
00253 G4cout << "°°° Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
00254
00255 G4cout << "°°° G4DNAChampionElasticModel - XS INFO END" << G4endl;
00256 }
00257
00258 }
00259
00260 return sigma*waterDensity;
00261
00262 }
00263
00264
00265
00266 void G4DNAChampionElasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* ,
00267 const G4MaterialCutsCouple* ,
00268 const G4DynamicParticle* aDynamicElectron,
00269 G4double,
00270 G4double)
00271 {
00272
00273 if (verboseLevel > 3)
00274 G4cout << "Calling SampleSecondaries() of G4DNAChampionElasticModel" << G4endl;
00275
00276 G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
00277
00278 if (electronEnergy0 < killBelowEnergy)
00279 {
00280 fParticleChangeForGamma->SetProposedKineticEnergy(0.);
00281 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
00282 fParticleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0);
00283 return ;
00284 }
00285
00286 if (electronEnergy0>= killBelowEnergy && electronEnergy0 < highEnergyLimit)
00287 {
00288
00289 G4double cosTheta = RandomizeCosTheta(electronEnergy0);
00290
00291 G4double phi = 2. * pi * G4UniformRand();
00292
00293 G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
00294 G4ThreeVector xVers = zVers.orthogonal();
00295 G4ThreeVector yVers = zVers.cross(xVers);
00296
00297 G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
00298 G4double yDir = xDir;
00299 xDir *= std::cos(phi);
00300 yDir *= std::sin(phi);
00301
00302 G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
00303
00304 fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit()) ;
00305
00306 fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
00307 }
00308
00309 }
00310
00311
00312
00313 G4double G4DNAChampionElasticModel::Theta
00314 (G4ParticleDefinition * particleDefinition, G4double k, G4double integrDiff)
00315 {
00316 G4double theta = 0.;
00317 G4double valueT1 = 0;
00318 G4double valueT2 = 0;
00319 G4double valueE21 = 0;
00320 G4double valueE22 = 0;
00321 G4double valueE12 = 0;
00322 G4double valueE11 = 0;
00323 G4double xs11 = 0;
00324 G4double xs12 = 0;
00325 G4double xs21 = 0;
00326 G4double xs22 = 0;
00327
00328 if (particleDefinition == G4Electron::ElectronDefinition())
00329 {
00330 std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
00331 std::vector<double>::iterator t1 = t2-1;
00332
00333 std::vector<double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), integrDiff);
00334 std::vector<double>::iterator e11 = e12-1;
00335
00336 std::vector<double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), integrDiff);
00337 std::vector<double>::iterator e21 = e22-1;
00338
00339 valueT1 =*t1;
00340 valueT2 =*t2;
00341 valueE21 =*e21;
00342 valueE22 =*e22;
00343 valueE12 =*e12;
00344 valueE11 =*e11;
00345
00346 xs11 = eDiffCrossSectionData[valueT1][valueE11];
00347 xs12 = eDiffCrossSectionData[valueT1][valueE12];
00348 xs21 = eDiffCrossSectionData[valueT2][valueE21];
00349 xs22 = eDiffCrossSectionData[valueT2][valueE22];
00350 }
00351
00352 if (xs11==0 && xs12==0 && xs21==0 && xs22==0) return (0.);
00353
00354 theta = QuadInterpolator ( valueE11, valueE12,
00355 valueE21, valueE22,
00356 xs11, xs12,
00357 xs21, xs22,
00358 valueT1, valueT2,
00359 k, integrDiff );
00360
00361 return theta;
00362 }
00363
00364
00365
00366 G4double G4DNAChampionElasticModel::LinLogInterpolate(G4double e1,
00367 G4double e2,
00368 G4double e,
00369 G4double xs1,
00370 G4double xs2)
00371 {
00372 G4double d1 = std::log(xs1);
00373 G4double d2 = std::log(xs2);
00374 G4double value = std::exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
00375 return value;
00376 }
00377
00378
00379
00380 G4double G4DNAChampionElasticModel::LinLinInterpolate(G4double e1,
00381 G4double e2,
00382 G4double e,
00383 G4double xs1,
00384 G4double xs2)
00385 {
00386 G4double d1 = xs1;
00387 G4double d2 = xs2;
00388 G4double value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
00389 return value;
00390 }
00391
00392
00393
00394 G4double G4DNAChampionElasticModel::LogLogInterpolate(G4double e1,
00395 G4double e2,
00396 G4double e,
00397 G4double xs1,
00398 G4double xs2)
00399 {
00400 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
00401 G4double b = std::log10(xs2) - a*std::log10(e2);
00402 G4double sigma = a*std::log10(e) + b;
00403 G4double value = (std::pow(10.,sigma));
00404 return value;
00405 }
00406
00407
00408
00409
00410 G4double G4DNAChampionElasticModel::QuadInterpolator(G4double e11, G4double e12,
00411 G4double e21, G4double e22,
00412 G4double xs11, G4double xs12,
00413 G4double xs21, G4double xs22,
00414 G4double t1, G4double t2,
00415 G4double t, G4double e)
00416 {
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431 G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12);
00432 G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22);
00433 G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
00434
00435 return value;
00436 }
00437
00438
00439
00440 G4double G4DNAChampionElasticModel::RandomizeCosTheta(G4double k)
00441 {
00442
00443 G4double integrdiff=0;
00444 G4double uniformRand=G4UniformRand();
00445 integrdiff = uniformRand;
00446
00447 G4double theta=0.;
00448 G4double cosTheta=0.;
00449 theta = Theta(G4Electron::ElectronDefinition(),k/eV,integrdiff);
00450
00451 cosTheta= std::cos(theta*pi/180);
00452
00453 return cosTheta;
00454 }