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