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