| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // Reference:
27 // A.D. Dominguez-Munoz, M.I. Gallardo, M.C
28 // Z. Francis, S. Incerti, M.A. Cortes-Gira
29 // Radiat. Phys. Chem. 199 (2022) 110363.
30 //
31 // Class authors:
32 // A.D. Dominguez-Munoz
33 // M.A. Cortes-Giraldo (miancortes -at- us.
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........oo
47
48 G4DNARPWBAExcitationModel::G4DNARPWBAExcitatio
49 const G4ParticleDefinition*, const G4String&
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 o
57 // 4 = entering in methods
58 if(verboseLevel > 0)
59 {
60 G4cout << "RPWBA excitation model is const
61 }
62 }
63
64 //....oooOO0OOooo........oooOO0OOooo........oo
65
66 G4DNARPWBAExcitationModel::~G4DNARPWBAExcitati
67
68 //....oooOO0OOooo........oooOO0OOooo........oo
69
70 void G4DNARPWBAExcitationModel::Initialise(con
71 con
72 {
73 if(isInitialised)
74 {
75 return;
76 }
77 if(verboseLevel > 3)
78 {
79 G4cout << "Calling G4DNARPWBAExcitationMod
80 }
81
82 if(fParticleDefinition != nullptr && fPartic
83 {
84 G4Exception("G4DNARPWBAExcitationModel::In
85 FatalException,
86 "Model already initialized for
87 }
88
89 fTableFile = "dna/sigma_excitation_p_RPWBA"
90 fLowEnergy = 100. * MeV;
91 fHighEnergy = 300. * MeV;
92
93 if(LowEnergyLimit() < fLowEnergy || HighEner
94 {
95 G4ExceptionDescription ed;
96 ed << "Model is applicable from "<<fLowEne
97 G4Exception("G4DNARPWBAExcitationModel::In
98 FatalException, ed);
99 }
100
101 G4double scaleFactor = 1 * cm * cm;
102 fTableData = make_unique<G4DNACrossSectionDa
103
104 fTableData->LoadData(fTableFile);
105
106 if(verboseLevel > 0)
107 {
108 G4cout << "RPWBA excitation model is initi
109 << "Energy range: " << LowEnergyLim
110 << HighEnergyLimit() / keV << " keV
111 << particle->GetParticleName() << G
112 }
113
114 // Initialize water density pointer
115 if(G4Material::GetMaterial("G4_WATER") != nu
116 fpMolWaterDensity =
117 G4DNAMolecularMaterial::Instance()->GetN
118 G4Material::GetMaterial("G4_WATER"));
119 }else{
120 G4ExceptionDescription exceptionDescriptio
121 exceptionDescription << "G4_WATER does not
122 G4Exception("G4DNARPWBAIonisationModel::In
123 FatalException, exceptionDescr
124 }
125 fParticleChangeForGamma = GetParticleChangeF
126 isInitialised = true;
127 }
128
129 //....oooOO0OOooo........oooOO0OOooo........oo
130
131 G4double G4DNARPWBAExcitationModel::CrossSecti
132 const G4Material* material, const G4Particle
133 G4double ekin, G4double, G4double)
134 {
135 if(verboseLevel > 3)
136 {
137 G4cout << "Calling CrossSectionPerVolume()
138 << G4endl;
139 }
140
141 if(fTableData == nullptr)
142 {
143 G4ExceptionDescription exceptionDescriptio
144 exceptionDescription << "No cross section
145 G4Exception("G4DNARPWBAIonisationModel::Cr
146 FatalException, exceptionDescr
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)
157
158 if(ekin >= fLowEnergy && ekin <= fHighEnergy
159 {
160 sigma = fTableData->FindValue(ekin);
161 }
162
163 if(verboseLevel > 2)
164 {
165 G4cout << "_______________________________
166 G4cout << "G4DNARPWBAExcitationModel - XS
167 G4cout << "Kinetic energy(eV)=" << ekin /
168 << " particle : " << particleDefini
169 G4cout << "Cross section per water molecul
170 << G4endl;
171 G4cout << "Cross section per water molecul
172 << sigma * waterDensity / (1. / cm)
173 G4cout << "G4DNARPWBAExcitationModel - XS
174 }
175
176 return sigma * waterDensity;
177 }
178
179 //....oooOO0OOooo........oooOO0OOooo........oo
180
181 void G4DNARPWBAExcitationModel::SampleSecondar
182 std::vector<G4DynamicParticle*>* /*fvect*/,
183 const G4MaterialCutsCouple* /*couple*/,
184 const G4DynamicParticle* aDynamicParticle, G
185 {
186 if(verboseLevel > 3)
187 {
188 G4cout << "Calling SampleSecondaries() of
189 << G4endl;
190 }
191
192 G4double k = aDynamicParticle->GetKineticEne
193
194 G4int level = RandomSelect(k);
195 G4double excitationEnergy = waterStructure.E
196 G4double newEnergy = k - excitationEn
197
198 if(newEnergy > 0)
199 {
200 fParticleChangeForGamma->ProposeMomentumDi
201 aDynamicParticle->GetMomentumDirection()
202
203 if(!statCode){
204 fParticleChangeForGamma->SetProposedKine
205 }
206 else{
207 fParticleChangeForGamma->SetProposedKine
208 }
209 fParticleChangeForGamma->ProposeLocalEnerg
210 }
211
212 const G4Track* theIncomingTrack = fParticleC
213 G4DNAChemistryManager::Instance()->CreateWat
214 eExcitedMolecule, level, theIncomingTrack)
215 }
216
217 //....oooOO0OOooo........oooOO0OOooo........oo
218
219 G4double G4DNARPWBAExcitationModel::GetPartial
220 const G4Material*, G4int level, const G4Part
221 G4double kineticEnergy)
222 {
223 if(fParticleDefinition != particle)
224 {
225 G4Exception("G4DNARPWBAExcitationModel::Ge
226 "RPWBAParticleType", FatalExce
227 "Model initialized for another
228 }
229
230 return fTableData->GetComponent(level)->Find
231 }
232
233 //....oooOO0OOooo........oooOO0OOooo........oo
234
235 G4int G4DNARPWBAExcitationModel::RandomSelect(
236 {
237 G4int level = 0;
238
239 auto valuesBuffer = new G4double[fTableData
240 const auto n = (G4int)fTableData->NumberOfC
241 G4int i(n);
242 G4double value = 0.;
243
244 while(i > 0)
245 {
246 --i;
247 valuesBuffer[i] = fTableData->GetComponent
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