| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // Authors: S. Meylan and C. Villagrasa (IRSN, 26 // Authors: S. Meylan and C. Villagrasa (IRSN, France)
27 // Models come from 27 // Models come from
28 // M. Bug et al, Rad. Phys and Chem. 130, 459- 28 // M. Bug et al, Rad. Phys and Chem. 130, 459-479 (2017)
29 // 29 //
30 30
31 #include "G4DNAPTBExcitationModel.hh" 31 #include "G4DNAPTBExcitationModel.hh"
32 << 32 #include "G4SystemOfUnits.hh"
33 #include "G4DNAChemistryManager.hh" 33 #include "G4DNAChemistryManager.hh"
34 #include "G4DNAMaterialManager.hh" <<
35 #include "G4DNAMolecularMaterial.hh" 34 #include "G4DNAMolecularMaterial.hh"
36 #include "G4SystemOfUnits.hh" <<
37 35
38 G4DNAPTBExcitationModel::G4DNAPTBExcitationMod << 36 G4DNAPTBExcitationModel::G4DNAPTBExcitationModel(const G4String& applyToMaterial, const G4ParticleDefinition*,
39 << 37 const G4String& nam)
40 : G4VDNAModel(nam, applyToMaterial) << 38 : G4VDNAModel(nam, applyToMaterial)
41 { 39 {
42 fpTHF = G4Material::GetMaterial("THF", false << 40 verboseLevel= 0;
43 fpPY = G4Material::GetMaterial("PY", false); << 41 // Verbosity scale:
44 fpPU = G4Material::GetMaterial("PU", false); << 42 // 0 = nothing
45 fpTMP = G4Material::GetMaterial("TMP", false << 43 // 1 = warning for energy non-conservation
46 fpG4_WATER = G4Material::GetMaterial("G4_WAT << 44 // 2 = details of energy budget
47 fpBackbone_THF = G4Material::GetMaterial("ba << 45 // 3 = calculation of cross sections, file openings, sampling of atoms
48 fpCytosine_PY = G4Material::GetMaterial("cyt << 46 // 4 = entering in methods
49 fpThymine_PY = G4Material::GetMaterial("thym << 47
50 fpAdenine_PU = G4Material::GetMaterial("aden << 48 // initialisation of mean energy loss for each material
51 fpBackbone_TMP = G4Material::GetMaterial("ba << 49 tableMeanEnergyPTB["THF"] = 8.01*eV;
52 fpGuanine_PU = G4Material::GetMaterial("guan << 50 tableMeanEnergyPTB["PY"] = 7.61*eV;
53 fpN2 = G4Material::GetMaterial("N2", false); << 51 tableMeanEnergyPTB["PU"] = 7.61*eV;
54 // initialisation of mean energy loss for ea << 52 tableMeanEnergyPTB["TMP"] = 8.01*eV;
55 << 53
56 if (fpTHF != nullptr) { << 54 if( verboseLevel>0 )
57 fTableMeanEnergyPTB[fpTHF->GetIndex()] = 8 << 55 {
58 } << 56 G4cout << "PTB excitation model is constructed " << G4endl;
59 << 57 }
60 if (fpPY != nullptr) { << 58 }
61 fTableMeanEnergyPTB[fpPY->GetIndex()] = 7. << 59
62 } << 60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63 << 61
64 if (fpPU != nullptr) { << 62 G4DNAPTBExcitationModel::~G4DNAPTBExcitationModel()
65 fTableMeanEnergyPTB[fpPU->GetIndex()] = 7. << 63 {
66 } << 64
67 if (fpTMP != nullptr) { <<
68 fTableMeanEnergyPTB[fpTMP->GetIndex()] = 8 <<
69 } <<
70 <<
71 if (verboseLevel > 0) { <<
72 G4cout << "PTB excitation model is constru <<
73 } <<
74 } 65 }
75 66
76 //....oooOO0OOooo........oooOO0OOooo........oo 67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77 68
78 void G4DNAPTBExcitationModel::Initialise(const 69 void G4DNAPTBExcitationModel::Initialise(const G4ParticleDefinition* particle,
79 const << 70 const G4DataVector& /*cuts*/, G4ParticleChangeForGamma*)
80 { 71 {
81 if (isInitialised) { << 72 if (verboseLevel > 3)
82 return; << 73 G4cout << "Calling G4DNAPTBExcitationModel::Initialise()" << G4endl;
83 } << 74
84 if (verboseLevel > 3) << 75 G4double scaleFactor = 1e-16*cm*cm;
85 { << 76 G4double scaleFactorBorn = (1.e-22 / 3.343) * m*m;
86 G4cout << "Calling G4DNAPTBExcitationModel << 77
87 } << 78 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
88 << 79
89 if (particle != G4Electron::ElectronDefiniti << 80 //*******************************************************
90 std::ostringstream oss; << 81 // Cross section data
91 oss << " Model is not applied for this par << 82 //*******************************************************
92 G4Exception("G4DNAPTBExcitationModel::Init << 83
93 } << 84 if(particle == electronDef)
94 << 85 {
95 G4double scaleFactor = 1e-16 * cm * cm; << 86 G4String particleName = particle->GetParticleName();
96 G4double scaleFactorBorn = (1.e-22 / 3.343) << 87
97 << 88 AddCrossSectionData("THF",
98 //****************************************** << 89 particleName,
99 // Cross section data << 90 "dna/sigma_excitation_e-_PTB_THF",
100 //****************************************** << 91 scaleFactor);
101 std::size_t index; << 92 SetLowELimit("THF", particleName, 9.*eV);
102 if (fpTHF != nullptr) { << 93 SetHighELimit("THF", particleName, 1.*keV);
103 index = fpTHF->GetIndex(); << 94
104 AddCrossSectionData(index, particle, "dna/ << 95 AddCrossSectionData("PY",
105 SetLowELimit(index, particle, 9. * eV); << 96 particleName,
106 SetHighELimit(index, particle, 1. * keV); << 97 "dna/sigma_excitation_e-_PTB_PY",
107 } << 98 scaleFactor);
108 if (fpPY != nullptr) { << 99 SetLowELimit("PY", particleName, 9.*eV);
109 index = fpPY->GetIndex(); << 100 SetHighELimit("PY", particleName, 1.*keV);
110 AddCrossSectionData(index, particle, "dna/ << 101
111 SetLowELimit(index, particle, 9. * eV); << 102 AddCrossSectionData("PU",
112 SetHighELimit(index, particle, 1. * keV); << 103 particleName,
113 } << 104 "dna/sigma_excitation_e-_PTB_PU",
114 << 105 scaleFactor);
115 if (fpPU != nullptr) { << 106 SetLowELimit("PU", particleName, 9.*eV);
116 index = fpPU->GetIndex(); << 107 SetHighELimit("PU", particleName, 1.*keV);
117 AddCrossSectionData(index, particle, "dna/ << 108
118 SetLowELimit(index, particle, 9. * eV); << 109 AddCrossSectionData("TMP",
119 SetHighELimit(index, particle, 1. * keV); << 110 particleName,
120 } << 111 "dna/sigma_excitation_e-_PTB_TMP",
121 << 112 scaleFactor);
122 if (fpTMP != nullptr) { << 113 SetLowELimit("TMP", particleName, 9.*eV);
123 index = fpTMP->GetIndex(); << 114 SetHighELimit("TMP", particleName, 1.*keV);
124 AddCrossSectionData(index, particle, "dna/ << 115
125 SetLowELimit(index, particle, 9. * eV); << 116 AddCrossSectionData("G4_WATER",
126 SetHighELimit(index, particle, 1. * keV); << 117 particleName,
127 } << 118 "dna/sigma_excitation_e_born",
128 if (fpG4_WATER != nullptr) { << 119 scaleFactorBorn);
129 index = fpG4_WATER->GetIndex(); << 120 SetLowELimit("G4_WATER", particleName, 9.*eV);
130 AddCrossSectionData(index, particle, "dna/ << 121 SetHighELimit("G4_WATER", particleName, 1.*keV);
131 SetLowELimit(index, particle, 9. * eV); << 122
132 SetHighELimit(index, particle, 1. * keV); << 123 // DNA materials
133 } << 124 //
134 // DNA materials << 125 AddCrossSectionData("backbone_THF",
135 // << 126 particleName,
136 if (fpBackbone_THF != nullptr) { << 127 "dna/sigma_excitation_e-_PTB_THF",
137 index = fpBackbone_THF->GetIndex(); << 128 scaleFactor*33./30);
138 AddCrossSectionData(index, particle, "dna/ << 129 SetLowELimit("backbone_THF", particleName, 9.*eV);
139 SetLowELimit(index, particle, 9. * eV); << 130 SetHighELimit("backbone_THF", particleName, 1.*keV);
140 SetHighELimit(index, particle, 1. * keV); << 131
141 } << 132 AddCrossSectionData("cytosine_PY",
142 if (fpCytosine_PY != nullptr) { << 133 particleName,
143 index = fpCytosine_PY->GetIndex(); << 134 "dna/sigma_excitation_e-_PTB_PY",
144 AddCrossSectionData(index, particle, "dna/ << 135 scaleFactor*42./30);
145 SetLowELimit(index, particle, 9. * eV); << 136 SetLowELimit("cytosine_PY", particleName, 9.*eV);
146 SetHighELimit(index, particle, 1. * keV); << 137 SetHighELimit("cytosine_PY", particleName, 1.*keV);
147 } << 138
148 if (fpThymine_PY != nullptr) { << 139 AddCrossSectionData("thymine_PY",
149 index = fpThymine_PY->GetIndex(); << 140 particleName,
150 AddCrossSectionData(index, particle, "dna/ << 141 "dna/sigma_excitation_e-_PTB_PY",
151 SetLowELimit(index, particle, 9. * eV); << 142 scaleFactor*48./30);
152 SetHighELimit(index, particle, 1. * keV); << 143 SetLowELimit("thymine_PY", particleName, 9.*eV);
153 } << 144 SetHighELimit("thymine_PY", particleName, 1.*keV);
154 if (fpAdenine_PU != nullptr) { << 145
155 index = fpAdenine_PU->GetIndex(); << 146 AddCrossSectionData("adenine_PU",
156 AddCrossSectionData(index, particle, "dna/ << 147 particleName,
157 SetLowELimit(index, particle, 9. * eV); << 148 "dna/sigma_excitation_e-_PTB_PU",
158 SetHighELimit(index, particle, 1. * keV); << 149 scaleFactor*50./44);
159 } << 150 SetLowELimit("adenine_PU", particleName, 9.*eV);
160 if (fpGuanine_PU != nullptr) { << 151 SetHighELimit("adenine_PU", particleName, 1.*keV);
161 index = fpGuanine_PU->GetIndex(); << 152
162 AddCrossSectionData(index, particle, "dna/ << 153 AddCrossSectionData("guanine_PU",
163 SetLowELimit(index, particle, 9. * eV); << 154 particleName,
164 SetHighELimit(index, particle, 1. * keV); << 155 "dna/sigma_excitation_e-_PTB_PU",
165 } << 156 scaleFactor*56./44);
166 if (fpBackbone_TMP != nullptr) { << 157 SetLowELimit("guanine_PU", particleName, 9.*eV);
167 index = fpBackbone_TMP->GetIndex(); << 158 SetHighELimit("guanine_PU", particleName, 1.*keV);
168 AddCrossSectionData(index, particle, "dna/ << 159
169 SetLowELimit(index, particle, 9. * eV); << 160 AddCrossSectionData("backbone_TMP",
170 SetHighELimit(index, particle, 1. * keV); << 161 particleName,
171 } << 162 "dna/sigma_excitation_e-_PTB_TMP",
172 // MPietrzak, adding paths for N2 << 163 scaleFactor*33./50);
173 if (fpN2 != nullptr) { << 164 SetLowELimit("backbone_TMP", particleName, 9.*eV);
174 index = fpN2->GetIndex(); << 165 SetHighELimit("backbone_TMP", particleName, 1.*keV);
175 AddCrossSectionData(index, particle, "dna/ << 166
176 SetLowELimit(index, particle, 13. * eV); << 167 // MPietrzak, adding paths for N2
177 SetHighELimit(index, particle, 1.02 * MeV) << 168 AddCrossSectionData("N2",
178 } << 169 particleName,
179 if (!G4DNAMaterialManager::Instance()->IsLoc << 170 "dna/sigma_excitation_e-_PTB_N2",
180 // Load data << 171 scaleFactor);
181 LoadCrossSectionData(particle); << 172 SetLowELimit("N2", particleName, 13.*eV);
182 G4DNAMaterialManager::Instance()->SetMaste << 173 SetHighELimit("N2", particleName, 1.02*MeV);
183 fpModelData = this; << 174 // MPietrzak
184 } <<
185 else { <<
186 auto dataModel = dynamic_cast<G4DNAPTBExci <<
187 G4DNAMaterialManager::Instance()->GetMod <<
188 if (dataModel == nullptr) { <<
189 G4cout << "G4DNAPTBExcitationModel::Init <<
190 G4Exception("G4DNAPTBExcitationModel::In <<
191 "not good modelData"); <<
192 } <<
193 else { <<
194 fpModelData = dataModel; <<
195 } 175 }
196 } <<
197 176
198 fParticleChangeForGamma = GetParticleChangeF << 177 //*******************************************************
199 isInitialised = true; << 178 // Load data
>> 179 //*******************************************************
>> 180
>> 181 LoadCrossSectionData(particle->GetParticleName() );
>> 182
>> 183 //*******************************************************
>> 184 // Verbose
>> 185 //*******************************************************
>> 186
>> 187 if( verboseLevel>0 )
>> 188 {
>> 189 G4cout << "PTB excitation model is initialized " << G4endl;
>> 190 }
200 } 191 }
201 192
202 //....oooOO0OOooo........oooOO0OOooo........oo 193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
203 194
204 G4double G4DNAPTBExcitationModel::CrossSection << 195 G4double G4DNAPTBExcitationModel::CrossSectionPerVolume(const G4Material* /*material*/,
205 << 196 const G4String& materialName,
206 << 197 const G4ParticleDefinition* particleDefinition,
>> 198 G4double ekin,
>> 199 G4double /*emin*/,
207 200 G4double /*emax*/)
208 { 201 {
209 // Get the name of the current particle << 202 if (verboseLevel > 3)
210 G4String particleName = p->GetParticleName() << 203 G4cout << "Calling CrossSectionPerVolume() of G4DNAPTBExcitationModel" << G4endl;
211 const std::size_t& MatID = material->GetInde << 204
212 // initialise variables << 205 // Get the name of the current particle
213 G4double lowLim; << 206 G4String particleName = particleDefinition->GetParticleName();
214 G4double highLim; << 207
215 G4double sigma = 0; << 208 // initialise variables
216 << 209 G4double lowLim = 0;
217 // Get the low energy limit for the current << 210 G4double highLim = 0;
218 lowLim = fpModelData->GetLowELimit(MatID, p) << 211 G4double sigma=0;
219 << 212
220 // Get the high energy limit for the current << 213 // Get the low energy limit for the current particle
221 highLim = fpModelData->GetHighELimit(MatID, << 214 lowLim = GetLowELimit(materialName, particleName);
222 << 215
223 // Check that we are in the correct energy r << 216 // Get the high energy limit for the current particle
224 if (ekin >= lowLim && ekin < highLim) { << 217 highLim = GetHighELimit(materialName, particleName);
225 // Get the map with all the data tables << 218
226 auto Data = fpModelData->GetData(); << 219 // Check that we are in the correct energy range
227 << 220 if (ekin >= lowLim && ekin < highLim)
228 if ((*Data)[MatID][p] == nullptr) { << 221 {
229 G4Exception("G4DNAPTBExcitationModel::Cr << 222 // Get the map with all the data tables
230 "No model is registered"); << 223 TableMapData* tableData = GetTableData();
231 } << 224
232 // Retrieve the cross section value << 225 // Retrieve the cross section value
233 sigma = (*Data)[MatID][p]->FindValue(ekin) << 226 sigma = (*tableData)[materialName][particleName]->FindValue(ekin);
>> 227
>> 228 if (verboseLevel > 2)
>> 229 {
>> 230 G4cout << "__________________________________" << G4endl;
>> 231 G4cout << "°°° G4DNAPTBExcitationModel - XS INFO START" << G4endl;
>> 232 G4cout << "°°° Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
>> 233 G4cout << "°°° Cross section per "<< materialName <<" molecule (cm^2)=" << sigma/cm/cm << G4endl;
>> 234 G4cout << "°°° G4DNAPTBExcitationModel - XS INFO END" << G4endl;
>> 235 }
234 236
235 if (verboseLevel > 2) { <<
236 G4cout << "_____________________________ <<
237 G4cout << "°°° G4DNAPTBExcitationMode <<
238 G4cout << "°°° Kinetic energy(eV)=" < <<
239 G4cout << "°°° Cross section per " << <<
240 << G4endl; <<
241 G4cout << "°°° G4DNAPTBExcitationMode <<
242 } 237 }
243 } <<
244 238
245 // Return the cross section value << 239 // Return the cross section value
246 auto MolDensity = << 240 return sigma;
247 (*G4DNAMolecularMaterial::Instance()->GetN <<
248 return sigma * MolDensity; <<
249 } 241 }
250 242
251 //....oooOO0OOooo........oooOO0OOooo........oo 243 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
252 244
253 void G4DNAPTBExcitationModel::SampleSecondarie 245 void G4DNAPTBExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
254 << 246 const G4MaterialCutsCouple* /*couple*/,
>> 247 const G4String& materialName,
255 248 const G4DynamicParticle* aDynamicParticle,
256 << 249 G4ParticleChangeForGamma* particleChangeForGamma,
>> 250 G4double /*tmin*/,
>> 251 G4double /*tmax*/)
257 { 252 {
258 const std::size_t& materialID = (std::size_t << 253 if (verboseLevel > 3)
>> 254 G4cout << "Calling SampleSecondaries() of G4DNAPTBExcitationModel" << G4endl;
259 255
260 // Get the incident particle kinetic energy << 256 // Get the incident particle kinetic energy
261 G4double k = aDynamicParticle->GetKineticEne << 257 G4double k = aDynamicParticle->GetKineticEnergy();
262 // Get the particle name << 258 //Get the particle name
263 const auto& particle = aDynamicParticle->Get << 259 const G4String& particleName = aDynamicParticle->GetDefinition()->GetParticleName();
264 // Get the energy limits << 260 // Get the energy limits
265 G4double lowLim = fpModelData->GetLowELimit( << 261 G4double lowLim = GetLowELimit(materialName, particleName);
266 G4double highLim = fpModelData->GetHighELimi << 262 G4double highLim = GetHighELimit(materialName, particleName);
267 << 263
268 // Check if we are in the correct energy ran << 264 // Check if we are in the correct energy range
269 if (k >= lowLim && k < highLim) { << 265 if (k >= lowLim && k < highLim)
270 if (fpN2 != nullptr && materialID == fpN2- << 266 {
271 // Retrieve the excitation energy for th << 267 if(materialName=="N2")
272 G4int level = fpModelData->RandomSelectS << 268 {
273 G4double excitationEnergy = ptbExcitatio << 269
274 << 270 // Retrieve the excitation energy for the current material
275 // Calculate the new energy of the parti << 271 G4int level = RandomSelectShell(k,particleName,materialName);
276 G4double newEnergy = k - excitationEnerg << 272 G4double excitationEnergy = ptbExcitationStructure.ExcitationEnergy(level, materialName);
277 << 273
278 // Check that the new energy is above ze << 274 // Calculate the new energy of the particle
279 // Otherwise, do nothing. << 275 G4double newEnergy = k - excitationEnergy;
280 if (newEnergy > 0) { << 276
281 fParticleChangeForGamma->ProposeMoment << 277 // Check that the new energy is above zero before applying it the particle.
282 fParticleChangeForGamma->SetProposedKi << 278 // Otherwise, do nothing.
283 fParticleChangeForGamma->ProposeLocalE << 279 if (newEnergy > 0)
284 G4double ioniThres = ptbIonisationStru << 280 {
285 // if excitation energy greater than i << 281 particleChangeForGamma->ProposeMomentumDirection(aDynamicParticle->GetMomentumDirection());
286 if ((excitationEnergy > ioniThres) && << 282 particleChangeForGamma->SetProposedKineticEnergy(newEnergy);
287 fParticleChangeForGamma->ProposeLoca << 283 particleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
288 // energy of ejected electron << 284 G4double ioniThres = ptbIonisationStructure.IonisationEnergy(0,materialName);
289 G4double secondaryKinetic = excitati << 285 // if excitation energy greater than ionisation threshold, then autoionisaiton
290 // random direction << 286 if((excitationEnergy>ioniThres)&&(G4UniformRand()<0.5))
291 G4double cosTheta = 2 * G4UniformRan << 287 {
292 G4double sinTheta = std::sqrt(1. - c << 288 particleChangeForGamma->ProposeLocalEnergyDeposit(ioniThres);
293 G4double ux = sinTheta * std::cos(ph << 289 // energy of ejected electron
294 G4ThreeVector deltaDirection(ux, uy, << 290 G4double secondaryKinetic = excitationEnergy - ioniThres;
295 // Create the new particle with its << 291 // random direction
296 auto dp = new G4DynamicParticle(G4El << 292 G4double cosTheta = 2*G4UniformRand() - 1., phi = CLHEP::twopi*G4UniformRand();
297 fvect->push_back(dp); << 293 G4double sinTheta = std::sqrt(1. - cosTheta*cosTheta);
>> 294 G4double ux = sinTheta*std::cos(phi),
>> 295 uy = sinTheta*std::sin(phi),
>> 296 uz = cosTheta;
>> 297 G4ThreeVector deltaDirection(ux,uy,uz);
>> 298 // Create the new particle with its characteristics
>> 299 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
>> 300 fvect->push_back(dp);
>> 301 }
>> 302 } else {
>> 303 G4ExceptionDescription description;
>> 304 description<<"Kinetic energy <= 0 at "<<materialName<<" material !!!";
>> 305 G4Exception("G4DNAPTBExcitationModel::SampleSecondaries","",FatalException,description);
>> 306 }
>> 307 } else if(materialName!="G4_WATER"){
>> 308 // Retrieve the excitation energy for the current material
>> 309 G4double excitationEnergy = tableMeanEnergyPTB[materialName];
>> 310 // Calculate the new energy of the particle
>> 311 G4double newEnergy = k - excitationEnergy;
>> 312 // Check that the new energy is above zero before applying it the particle.
>> 313 // Otherwise, do nothing.
>> 314 if (newEnergy > 0){
>> 315 particleChangeForGamma->ProposeMomentumDirection(aDynamicParticle->GetMomentumDirection());
>> 316 particleChangeForGamma->SetProposedKineticEnergy(newEnergy);
>> 317 particleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
>> 318 } else {
>> 319 G4ExceptionDescription description;
>> 320 description<<"Kinetic energy <= 0 at "<<materialName<<" material !!!";
>> 321 G4Exception("G4DNAPTBExcitationModel::SampleSecondaries","",FatalException,description);
>> 322 }
>> 323 } else {
>> 324 G4int level = RandomSelectShell(k,particleName, materialName);
>> 325 G4double excitationEnergy = waterStructure.ExcitationEnergy(level);
>> 326 G4double newEnergy = k - excitationEnergy;
>> 327
>> 328 if (newEnergy > 0){
>> 329 particleChangeForGamma->ProposeMomentumDirection(aDynamicParticle->GetMomentumDirection());
>> 330 particleChangeForGamma->SetProposedKineticEnergy(newEnergy);
>> 331 particleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
>> 332 const G4Track * theIncomingTrack = particleChangeForGamma->GetCurrentTrack();
>> 333 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eExcitedMolecule,
>> 334 level,
>> 335 theIncomingTrack);
>> 336 } else {
>> 337 G4ExceptionDescription description;
>> 338 description<<"Kinetic energy <= 0 at "<<materialName<<" material !!!";
>> 339 G4Exception("G4DNAPTBExcitationModel::SampleSecondaries","",FatalException,description);
>> 340 }
298 } 341 }
299 } <<
300 else { <<
301 G4ExceptionDescription description; <<
302 description << "Kinetic energy <= 0 at <<
303 G4Exception("G4DNAPTBExcitationModel:: <<
304 } <<
305 } <<
306 else if (fpG4_WATER == nullptr || material <<
307 // Retrieve the excitation energy for th <<
308 G4double excitationEnergy = fTableMeanEn <<
309 // Calculate the new energy of the parti <<
310 G4double newEnergy = k - excitationEnerg <<
311 // Check that the new energy is above ze <<
312 // Otherwise, do nothing. <<
313 if (newEnergy > 0) { <<
314 fParticleChangeForGamma->ProposeMoment <<
315 fParticleChangeForGamma->SetProposedKi <<
316 fParticleChangeForGamma->ProposeLocalE <<
317 } <<
318 else { <<
319 G4ExceptionDescription description; <<
320 description << "Kinetic energy <= 0 at <<
321 G4Exception("G4DNAPTBExcitationModel:: <<
322 } <<
323 } <<
324 else { <<
325 G4int level = RandomSelectShell(k, parti <<
326 G4double excitationEnergy = waterStructu <<
327 G4double newEnergy = k - excitationEnerg <<
328 <<
329 if (newEnergy > 0) { <<
330 fParticleChangeForGamma->ProposeMoment <<
331 fParticleChangeForGamma->SetProposedKi <<
332 fParticleChangeForGamma->ProposeLocalE <<
333 const G4Track* theIncomingTrack = fPar <<
334 G4DNAChemistryManager::Instance()->Cre <<
335 <<
336 } <<
337 else { <<
338 G4ExceptionDescription description; <<
339 description << "Kinetic energy <= 0 at <<
340 G4Exception("G4DNAPTBExcitationModel:: <<
341 } <<
342 } 342 }
343 } << 343
344 } 344 }
345 345