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
00030
00031 #include "G4DNASancheExcitationModel.hh"
00032 #include "G4SystemOfUnits.hh"
00033 #include "G4DNAMolecularMaterial.hh"
00034
00035
00036
00037 using namespace std;
00038
00039
00040
00041 G4DNASancheExcitationModel::G4DNASancheExcitationModel(const G4ParticleDefinition*,
00042 const G4String& nam)
00043 :G4VEmModel(nam),isInitialised(false)
00044 {
00045
00046 fpWaterDensity = 0;
00047
00048 lowEnergyLimit = 2 * eV;
00049 highEnergyLimit = 100 * eV;
00050 SetLowEnergyLimit(lowEnergyLimit);
00051 SetHighEnergyLimit(highEnergyLimit);
00052 nLevels = 9;
00053
00054 verboseLevel= 0;
00055
00056
00057
00058
00059
00060
00061
00062 if (verboseLevel > 0)
00063 {
00064 G4cout << "Sanche Excitation model is constructed " << G4endl
00065 << "Energy range: "
00066 << lowEnergyLimit / eV << " eV - "
00067 << highEnergyLimit / eV << " eV"
00068 << G4endl;
00069 }
00070 fParticleChangeForGamma = 0;
00071 fpWaterDensity = 0;
00072 }
00073
00074
00075
00076 G4DNASancheExcitationModel::~G4DNASancheExcitationModel()
00077 {}
00078
00079
00080
00081 void G4DNASancheExcitationModel::Initialise(const G4ParticleDefinition* ,
00082 const G4DataVector& )
00083 {
00084
00085
00086 if (verboseLevel > 3)
00087 G4cout << "Calling G4DNASancheExcitationModel::Initialise()" << G4endl;
00088
00089
00090
00091 if (LowEnergyLimit() < lowEnergyLimit)
00092 {
00093 G4cout << "G4DNASancheExcitationModel: low energy limit increased from " <<
00094 LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
00095 SetLowEnergyLimit(lowEnergyLimit);
00096 }
00097
00098 if (HighEnergyLimit() > highEnergyLimit)
00099 {
00100 G4cout << "G4DNASancheExcitationModel: high energy limit decreased from " <<
00101 HighEnergyLimit()/eV << " eV to " << highEnergyLimit/eV << " eV" << G4endl;
00102 SetHighEnergyLimit(highEnergyLimit);
00103 }
00104
00105
00106
00107 if (verboseLevel > 0)
00108 {
00109 G4cout << "Sanche Excitation model is initialized " << G4endl
00110 << "Energy range: "
00111 << LowEnergyLimit() / eV << " eV - "
00112 << HighEnergyLimit() / eV << " eV"
00113 << G4endl;
00114 }
00115
00116
00117 fpWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
00118
00119 if (isInitialised) { return; }
00120 fParticleChangeForGamma = GetParticleChangeForGamma();
00121 isInitialised = true;
00122
00123 char *path = getenv("G4LEDATA");
00124 std::ostringstream eFullFileName;
00125 eFullFileName << path << "/dna/sigma_excitationvib_e_sanche.dat";
00126 std::ifstream input(eFullFileName.str().c_str());
00127
00128 if (!input)
00129 {
00130 G4Exception("G4DNASancheExcitationModel::Initialise","em0003",
00131 FatalException,"Missing data file:/dna/sigma_excitationvib_e_sanche.dat");
00132 }
00133
00134 while(!input.eof())
00135 {
00136 double t;
00137 input>>t;
00138 tdummyVec.push_back(t);
00139 input>>map1[t][0]>>map1[t][1]>>map1[t][2]>>map1[t][3]>>map1[t][4]>>map1[t][5]>>map1[t][6]>>map1[t][7]>>map1[t][8];
00140
00141 }
00142 }
00143
00144
00145
00146 G4double G4DNASancheExcitationModel::CrossSectionPerVolume(const G4Material* material,
00147 const G4ParticleDefinition* particleDefinition,
00148 G4double ekin,
00149 G4double,
00150 G4double)
00151 {
00152 if (verboseLevel > 3)
00153 G4cout << "Calling CrossSectionPerVolume() of G4DNASancheExcitationModel" << G4endl;
00154
00155
00156
00157 G4double sigma=0;
00158
00159 G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
00160
00161 if(waterDensity!= 0.0)
00162
00163 {
00164
00165 if (particleDefinition == G4Electron::ElectronDefinition())
00166 {
00167 if (ekin >= lowEnergyLimit && ekin < highEnergyLimit)
00168 {
00169 sigma = Sum(ekin);
00170 }
00171 }
00172
00173 if (verboseLevel > 2)
00174 {
00175 G4cout << "__________________________________" << G4endl;
00176 G4cout << "°°° G4DNASancheExcitationModel - XS INFO START" << G4endl;
00177 G4cout << "°°° Kinetic energy(eV)=" << ekin/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
00178 G4cout << "°°° Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
00179 G4cout << "°°° Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
00180
00181 G4cout << "°°° G4DNASancheExcitationModel - XS INFO END" << G4endl;
00182 }
00183
00184 }
00185
00186
00187
00188 return sigma*2*waterDensity;
00189
00190
00191 }
00192
00193
00194
00195 void G4DNASancheExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>*,
00196 const G4MaterialCutsCouple*,
00197 const G4DynamicParticle* aDynamicElectron,
00198 G4double,
00199 G4double)
00200 {
00201
00202 if (verboseLevel > 3)
00203 G4cout << "Calling SampleSecondaries() of G4DNASancheExcitationModel" << G4endl;
00204
00205 G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
00206 G4int level = RandomSelect(electronEnergy0);
00207 G4double excitationEnergy = VibrationEnergy(level);
00208 G4double newEnergy = electronEnergy0 - excitationEnergy;
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228 if (electronEnergy0 < highEnergyLimit && newEnergy>0.)
00229 {
00230 fParticleChangeForGamma->ProposeMomentumDirection(aDynamicElectron->GetMomentumDirection());
00231 fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
00232 fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
00233 }
00234
00235
00236 }
00237
00238
00239
00240 G4double G4DNASancheExcitationModel::PartialCrossSection(G4double t, G4int level)
00241 {
00242 std::vector<double>::iterator t2 = std::upper_bound(tdummyVec.begin(),tdummyVec.end(), t/eV);
00243 std::vector<double>::iterator t1 = t2-1;
00244
00245 double sigma = LinInterpolate((*t1), (*t2), t/eV, map1[*t1][level], map1[*t2][level]);
00246 sigma*=1e-16*cm*cm;
00247 if(sigma==0.)sigma=1e-30;
00248 return (sigma);
00249 }
00250
00251
00252
00253 G4double G4DNASancheExcitationModel::VibrationEnergy(G4int level)
00254 {
00255 G4double energies[9] = {0.01, 0.024, 0.061, 0.092, 0.204, 0.417, 0.460, 0.500, 0.835};
00256 return(energies[level]*eV);
00257 }
00258
00259
00260
00261 G4int G4DNASancheExcitationModel::RandomSelect(G4double k)
00262 {
00263
00264
00265
00266 G4int i = nLevels;
00267 G4double value = 0.;
00268 std::deque<double> values;
00269
00270 while (i > 0)
00271 {
00272 i--;
00273 G4double partial = PartialCrossSection(k,i);
00274 values.push_front(partial);
00275 value += partial;
00276 }
00277
00278 value *= G4UniformRand();
00279
00280 i = nLevels;
00281
00282 while (i > 0)
00283 {
00284 i--;
00285 if (values[i] > value)
00286 {
00287
00288 return i;
00289 }
00290 value -= values[i];
00291 }
00292
00293
00294
00295 return 0;
00296 }
00297
00298
00299
00300 G4double G4DNASancheExcitationModel::Sum(G4double k)
00301 {
00302 G4double totalCrossSection = 0.;
00303
00304 for (G4int i=0; i<nLevels; i++)
00305 {
00306 totalCrossSection += PartialCrossSection(k,i);
00307 }
00308 return totalCrossSection;
00309 }
00310
00311
00312
00313 G4double G4DNASancheExcitationModel::LinInterpolate(G4double e1,
00314 G4double e2,
00315 G4double e,
00316 G4double xs1,
00317 G4double xs2)
00318 {
00319 G4double a = (xs2 - xs1) / (e2 - e1);
00320 G4double b = xs2 - a*e2;
00321 G4double value = a*e + b;
00322
00323
00324 return value;
00325 }
00326