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