| 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 // Reference:
27 // A.D. Dominguez-Munoz, M.I. Gallardo, M.C. Bordage,
28 // Z. Francis, S. Incerti, M.A. Cortes-Giraldo,
29 // Radiat. Phys. Chem. 199 (2022) 110363.
30 //
31 // Class authors:
32 // A.D. Dominguez-Munoz
33 // M.A. Cortes-Giraldo (miancortes -at- us.es)
34 //
35 // Class creation: 2022-03-03
36 //
37 //
38
39 #include "G4DNARPWBAExcitationModel.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4DNAChemistryManager.hh"
42 #include "G4DNAMolecularMaterial.hh"
43 #include <map>
44 using namespace std;
45
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47
48 G4DNARPWBAExcitationModel::G4DNARPWBAExcitationModel(
49 const G4ParticleDefinition*, const G4String& nam)
50 : G4VEmModel(nam)
51 {
52 // Verbosity scale:
53 // 0 = nothing
54 // 1 = warning for energy non-conservation
55 // 2 = details of energy budget
56 // 3 = calculation of cross sections, file openings, sampling of atoms
57 // 4 = entering in methods
58 if(verboseLevel > 0)
59 {
60 G4cout << "RPWBA excitation model is constructed " << G4endl;
61 }
62 }
63
64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
65
66 G4DNARPWBAExcitationModel::~G4DNARPWBAExcitationModel() = default;
67
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
69
70 void G4DNARPWBAExcitationModel::Initialise(const G4ParticleDefinition* particle,
71 const G4DataVector& /*cuts*/)
72 {
73 if(isInitialised)
74 {
75 return;
76 }
77 if(verboseLevel > 3)
78 {
79 G4cout << "Calling G4DNARPWBAExcitationModel::Initialise()" << G4endl;
80 }
81
82 if(fParticleDefinition != nullptr && fParticleDefinition != particle)
83 {
84 G4Exception("G4DNARPWBAExcitationModel::Initialise", "em0001",
85 FatalException,
86 "Model already initialized for another particle type.");
87 }
88
89 fTableFile = "dna/sigma_excitation_p_RPWBA";
90 fLowEnergy = 100. * MeV;
91 fHighEnergy = 300. * MeV;
92
93 if(LowEnergyLimit() < fLowEnergy || HighEnergyLimit() > fHighEnergy)
94 {
95 G4ExceptionDescription ed;
96 ed << "Model is applicable from "<<fLowEnergy<<" to "<<fHighEnergy;
97 G4Exception("G4DNARPWBAExcitationModel::Initialise", "em0004",
98 FatalException, ed);
99 }
100
101 G4double scaleFactor = 1 * cm * cm;
102 fTableData = make_unique<G4DNACrossSectionDataSet>(new G4LogLogInterpolation,
103 eV, scaleFactor);
104 fTableData->LoadData(fTableFile);
105
106 if(verboseLevel > 0)
107 {
108 G4cout << "RPWBA excitation model is initialized " << G4endl
109 << "Energy range: " << LowEnergyLimit() / eV << " eV - "
110 << HighEnergyLimit() / keV << " keV for "
111 << particle->GetParticleName() << G4endl;
112 }
113
114 // Initialize water density pointer
115 if(G4Material::GetMaterial("G4_WATER") != nullptr){
116 fpMolWaterDensity =
117 G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(
118 G4Material::GetMaterial("G4_WATER"));
119 }else{
120 G4ExceptionDescription exceptionDescription;
121 exceptionDescription << "G4_WATER does not exist :";
122 G4Exception("G4DNARPWBAIonisationModel::Initialise", "em00020",
123 FatalException, exceptionDescription);
124 }
125 fParticleChangeForGamma = GetParticleChangeForGamma();
126 isInitialised = true;
127 }
128
129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
130
131 G4double G4DNARPWBAExcitationModel::CrossSectionPerVolume(
132 const G4Material* material, const G4ParticleDefinition* particleDefinition,
133 G4double ekin, G4double, G4double)
134 {
135 if(verboseLevel > 3)
136 {
137 G4cout << "Calling CrossSectionPerVolume() of G4DNARPWBAExcitationModel"
138 << G4endl;
139 }
140
141 if(fTableData == nullptr)
142 {
143 G4ExceptionDescription exceptionDescription;
144 exceptionDescription << "No cross section data ";
145 G4Exception("G4DNARPWBAIonisationModel::CrossSectionPerVolume", "em00120",
146 FatalException, exceptionDescription);
147 }
148
149 if(particleDefinition != fParticleDefinition)
150 return 0;
151
152 // Calculate total cross section for model
153
154 G4double sigma = 0;
155
156 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
157
158 if(ekin >= fLowEnergy && ekin <= fHighEnergy)
159 {
160 sigma = fTableData->FindValue(ekin);
161 }
162
163 if(verboseLevel > 2)
164 {
165 G4cout << "__________________________________" << G4endl;
166 G4cout << "G4DNARPWBAExcitationModel - XS INFO START" << G4endl;
167 G4cout << "Kinetic energy(eV)=" << ekin / eV
168 << " particle : " << particleDefinition->GetParticleName() << G4endl;
169 G4cout << "Cross section per water molecule (cm^2)=" << sigma / cm / cm
170 << G4endl;
171 G4cout << "Cross section per water molecule (cm^-1)="
172 << sigma * waterDensity / (1. / cm) << G4endl;
173 G4cout << "G4DNARPWBAExcitationModel - XS INFO END" << G4endl;
174 }
175
176 return sigma * waterDensity;
177 }
178
179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
180
181 void G4DNARPWBAExcitationModel::SampleSecondaries(
182 std::vector<G4DynamicParticle*>* /*fvect*/,
183 const G4MaterialCutsCouple* /*couple*/,
184 const G4DynamicParticle* aDynamicParticle, G4double, G4double)
185 {
186 if(verboseLevel > 3)
187 {
188 G4cout << "Calling SampleSecondaries() of G4DNARPWBAExcitationModel"
189 << G4endl;
190 }
191
192 G4double k = aDynamicParticle->GetKineticEnergy();
193
194 G4int level = RandomSelect(k);
195 G4double excitationEnergy = waterStructure.ExcitationEnergy(level);
196 G4double newEnergy = k - excitationEnergy;
197
198 if(newEnergy > 0)
199 {
200 fParticleChangeForGamma->ProposeMomentumDirection(
201 aDynamicParticle->GetMomentumDirection());
202
203 if(!statCode){
204 fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
205 }
206 else{
207 fParticleChangeForGamma->SetProposedKineticEnergy(k);
208 }
209 fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
210 }
211
212 const G4Track* theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
213 G4DNAChemistryManager::Instance()->CreateWaterMolecule(
214 eExcitedMolecule, level, theIncomingTrack);
215 }
216
217 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
218
219 G4double G4DNARPWBAExcitationModel::GetPartialCrossSection(
220 const G4Material*, G4int level, const G4ParticleDefinition* particle,
221 G4double kineticEnergy)
222 {
223 if(fParticleDefinition != particle)
224 {
225 G4Exception("G4DNARPWBAExcitationModel::GetPartialCrossSection",
226 "RPWBAParticleType", FatalException,
227 "Model initialized for another particle type.");
228 }
229
230 return fTableData->GetComponent(level)->FindValue(kineticEnergy);
231 }
232
233 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
234
235 G4int G4DNARPWBAExcitationModel::RandomSelect(G4double k)
236 {
237 G4int level = 0;
238
239 auto valuesBuffer = new G4double[fTableData->NumberOfComponents()];
240 const auto n = (G4int)fTableData->NumberOfComponents();
241 G4int i(n);
242 G4double value = 0.;
243
244 while(i > 0)
245 {
246 --i;
247 valuesBuffer[i] = fTableData->GetComponent(i)->FindValue(k);
248 value += valuesBuffer[i];
249 }
250
251 value *= G4UniformRand();
252 i = n;
253
254 while(i > 0)
255 {
256 --i;
257
258 if(valuesBuffer[i] > value)
259 {
260 delete[] valuesBuffer;
261 return i;
262 }
263 value -= valuesBuffer[i];
264 }
265 delete[] valuesBuffer;
266
267 return level;
268 }
269