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 "G4DNABornExcitationModel.hh"
00030 #include "G4SystemOfUnits.hh"
00031 #include "G4DNAChemistryManager.hh"
00032 #include "G4DNAMolecularMaterial.hh"
00033
00034
00035
00036 using namespace std;
00037
00038
00039
00040 G4DNABornExcitationModel::G4DNABornExcitationModel(const G4ParticleDefinition*,
00041 const G4String& nam)
00042 :G4VEmModel(nam),isInitialised(false)
00043 {
00044
00045 fpMolWaterDensity = 0;
00046
00047 verboseLevel= 0;
00048
00049
00050
00051
00052
00053
00054
00055 if( verboseLevel>0 )
00056 {
00057 G4cout << "Born excitation model is constructed " << G4endl;
00058 }
00059 fParticleChangeForGamma = 0;
00060 }
00061
00062
00063
00064 G4DNABornExcitationModel::~G4DNABornExcitationModel()
00065 {
00066
00067
00068 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
00069 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
00070 {
00071 G4DNACrossSectionDataSet* table = pos->second;
00072 delete table;
00073 }
00074
00075 }
00076
00077
00078
00079 void G4DNABornExcitationModel::Initialise(const G4ParticleDefinition* particle,
00080 const G4DataVector& )
00081 {
00082
00083 if (verboseLevel > 3)
00084 G4cout << "Calling G4DNABornExcitationModel::Initialise()" << G4endl;
00085
00086 G4String fileElectron("dna/sigma_excitation_e_born");
00087 G4String fileProton("dna/sigma_excitation_p_born");
00088
00089 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
00090 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
00091
00092 G4String electron;
00093 G4String proton;
00094
00095 G4double scaleFactor = (1.e-22 / 3.343) * m*m;
00096
00097
00098
00099 electron = electronDef->GetParticleName();
00100
00101 tableFile[electron] = fileElectron;
00102
00103 lowEnergyLimit[electron] = 9. * eV;
00104 highEnergyLimit[electron] = 1. * MeV;
00105
00106
00107
00108 G4DNACrossSectionDataSet* tableE = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
00109 tableE->LoadData(fileElectron);
00110
00111 tableData[electron] = tableE;
00112
00113
00114
00115 proton = protonDef->GetParticleName();
00116
00117 tableFile[proton] = fileProton;
00118
00119 lowEnergyLimit[proton] = 500. * keV;
00120 highEnergyLimit[proton] = 100. * MeV;
00121
00122
00123
00124 G4DNACrossSectionDataSet* tableP = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
00125 tableP->LoadData(fileProton);
00126
00127 tableData[proton] = tableP;
00128
00129
00130
00131 if (particle==electronDef)
00132 {
00133 SetLowEnergyLimit(lowEnergyLimit[electron]);
00134 SetHighEnergyLimit(highEnergyLimit[electron]);
00135 }
00136
00137 if (particle==protonDef)
00138 {
00139 SetLowEnergyLimit(lowEnergyLimit[proton]);
00140 SetHighEnergyLimit(highEnergyLimit[proton]);
00141 }
00142
00143 if( verboseLevel>0 )
00144 {
00145 G4cout << "Born excitation model is initialized " << G4endl
00146 << "Energy range: "
00147 << LowEnergyLimit() / eV << " eV - "
00148 << HighEnergyLimit() / keV << " keV for "
00149 << particle->GetParticleName()
00150 << G4endl;
00151 }
00152
00153
00154 fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
00155
00156 if (isInitialised) { return; }
00157 fParticleChangeForGamma = GetParticleChangeForGamma();
00158 isInitialised = true;
00159 }
00160
00161
00162
00163 G4double G4DNABornExcitationModel::CrossSectionPerVolume(const G4Material* material,
00164 const G4ParticleDefinition* particleDefinition,
00165 G4double ekin,
00166 G4double,
00167 G4double)
00168 {
00169 if (verboseLevel > 3)
00170 G4cout << "Calling CrossSectionPerVolume() of G4DNABornExcitationModel" << G4endl;
00171
00172 if (
00173 particleDefinition != G4Proton::ProtonDefinition()
00174 &&
00175 particleDefinition != G4Electron::ElectronDefinition()
00176 )
00177
00178 return 0;
00179
00180
00181
00182 G4double lowLim = 0;
00183 G4double highLim = 0;
00184 G4double sigma=0;
00185
00186 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
00187
00188 if(waterDensity!= 0.0)
00189
00190 {
00191 const G4String& particleName = particleDefinition->GetParticleName();
00192
00193 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
00194 pos1 = lowEnergyLimit.find(particleName);
00195 if (pos1 != lowEnergyLimit.end())
00196 {
00197 lowLim = pos1->second;
00198 }
00199
00200 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
00201 pos2 = highEnergyLimit.find(particleName);
00202 if (pos2 != highEnergyLimit.end())
00203 {
00204 highLim = pos2->second;
00205 }
00206
00207 if (ekin >= lowLim && ekin < highLim)
00208 {
00209 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
00210 pos = tableData.find(particleName);
00211
00212 if (pos != tableData.end())
00213 {
00214 G4DNACrossSectionDataSet* table = pos->second;
00215 if (table != 0)
00216 {
00217 sigma = table->FindValue(ekin);
00218 }
00219 }
00220 else
00221 {
00222 G4Exception("G4DNABornExcitationModel::CrossSectionPerVolume","em0002",
00223 FatalException,"Model not applicable to particle type.");
00224 }
00225 }
00226
00227 if (verboseLevel > 2)
00228 {
00229 G4cout << "__________________________________" << G4endl;
00230 G4cout << "°°° G4DNABornExcitationModel - XS INFO START" << G4endl;
00231 G4cout << "°°° Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
00232 G4cout << "°°° Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
00233 G4cout << "°°° Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
00234
00235 G4cout << "°°° G4DNABornExcitationModel - XS INFO END" << G4endl;
00236 }
00237
00238 }
00239
00240 return sigma*waterDensity;
00241
00242 }
00243
00244
00245
00246 void G4DNABornExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* ,
00247 const G4MaterialCutsCouple* ,
00248 const G4DynamicParticle* aDynamicParticle,
00249 G4double,
00250 G4double)
00251 {
00252
00253 if (verboseLevel > 3)
00254 G4cout << "Calling SampleSecondaries() of G4DNABornExcitationModel" << G4endl;
00255
00256 G4double k = aDynamicParticle->GetKineticEnergy();
00257
00258 const G4String& particleName = aDynamicParticle->GetDefinition()->GetParticleName();
00259
00260 G4int level = RandomSelect(k,particleName);
00261 G4double excitationEnergy = waterStructure.ExcitationEnergy(level);
00262 G4double newEnergy = k - excitationEnergy;
00263
00264 if (newEnergy > 0)
00265 {
00266 fParticleChangeForGamma->ProposeMomentumDirection(aDynamicParticle->GetMomentumDirection());
00267 fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
00268 fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
00269 }
00270
00271 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
00272 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eExcitedMolecule,
00273 level,
00274 theIncomingTrack);
00275 }
00276
00277
00278
00279 G4int G4DNABornExcitationModel::RandomSelect(G4double k, const G4String& particle)
00280 {
00281 G4int level = 0;
00282
00283 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
00284 pos = tableData.find(particle);
00285
00286 if (pos != tableData.end())
00287 {
00288 G4DNACrossSectionDataSet* table = pos->second;
00289
00290 if (table != 0)
00291 {
00292 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
00293 const size_t n(table->NumberOfComponents());
00294 size_t i(n);
00295 G4double value = 0.;
00296
00297 while (i>0)
00298 {
00299 i--;
00300 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
00301 value += valuesBuffer[i];
00302 }
00303
00304 value *= G4UniformRand();
00305
00306 i = n;
00307
00308 while (i > 0)
00309 {
00310 i--;
00311
00312 if (valuesBuffer[i] > value)
00313 {
00314 delete[] valuesBuffer;
00315 return i;
00316 }
00317 value -= valuesBuffer[i];
00318 }
00319
00320 if (valuesBuffer) delete[] valuesBuffer;
00321
00322 }
00323 }
00324 else
00325 {
00326 G4Exception("G4DNABornExcitationModel::RandomSelect","em0002",
00327 FatalException,"Model not applicable to particle type.");
00328 }
00329 return level;
00330 }
00331
00332