| 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 "G4DNABornExcitationModel2.hh" 28 #include "G4DNABornExcitationModel2.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 "G4PhysicsTable.hh" 32 #include "G4PhysicsTable.hh"
33 #include "G4PhysicsVector.hh" 33 #include "G4PhysicsVector.hh"
34 #include "G4UnitsTable.hh" 34 #include "G4UnitsTable.hh"
35 #include <map> 35 #include <map>
36 36
37 //....oooOO0OOooo........oooOO0OOooo........oo 37 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
38 38
39 using namespace std; 39 using namespace std;
40 40
41 //....oooOO0OOooo........oooOO0OOooo........oo 41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42 42
43 G4DNABornExcitationModel2::G4DNABornExcitation 43 G4DNABornExcitationModel2::G4DNABornExcitationModel2(const G4ParticleDefinition*,
44 44 const G4String& nam) :
45 G4VEmModel(nam) << 45 G4VEmModel(nam), isInitialised(false), fTableData(0)
46 { 46 {
47 fpMolWaterDensity = nullptr; << 47 fpMolWaterDensity = 0;
48 fHighEnergy = 0; 48 fHighEnergy = 0;
49 fLowEnergy = 0; 49 fLowEnergy = 0;
50 fParticleDefinition = nullptr; << 50 fParticleDefinition = 0;
51 51
52 verboseLevel = 0; 52 verboseLevel = 0;
53 // Verbosity scale: 53 // Verbosity scale:
54 // 0 = nothing 54 // 0 = nothing
55 // 1 = warning for energy non-conservation 55 // 1 = warning for energy non-conservation
56 // 2 = details of energy budget 56 // 2 = details of energy budget
57 // 3 = calculation of cross sections, file o 57 // 3 = calculation of cross sections, file openings, sampling of atoms
58 // 4 = entering in methods 58 // 4 = entering in methods
59 59
60 if (verboseLevel > 0) 60 if (verboseLevel > 0)
61 { 61 {
62 G4cout << "Born excitation model is constr 62 G4cout << "Born excitation model is constructed " << G4endl;
63 } 63 }
64 fParticleChangeForGamma = nullptr; << 64 fParticleChangeForGamma = 0;
65 fLastBinCallForFinalXS = 0; 65 fLastBinCallForFinalXS = 0;
66 fTotalXS = nullptr; << 66 fTotalXS = 0;
67 fTableData = nullptr; << 67 fTableData = 0;
68 68
69 // Selection of stationary mode 69 // Selection of stationary mode
70 70
71 statCode = false; 71 statCode = false;
72 } 72 }
73 73
74 //....oooOO0OOooo........oooOO0OOooo........oo 74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75 75
76 G4DNABornExcitationModel2::~G4DNABornExcitatio 76 G4DNABornExcitationModel2::~G4DNABornExcitationModel2()
77 { 77 {
78 // Cross section 78 // Cross section
79 << 79 if (fTableData)
80 delete fTableData; 80 delete fTableData;
81 } 81 }
82 82
83 //....oooOO0OOooo........oooOO0OOooo........oo 83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
84 84
85 void G4DNABornExcitationModel2::Initialise(con 85 void G4DNABornExcitationModel2::Initialise(const G4ParticleDefinition* particle,
86 con 86 const G4DataVector& /*cuts*/)
87 { 87 {
88 88
89 if (verboseLevel > 3) 89 if (verboseLevel > 3)
90 { 90 {
91 G4cout << "Calling G4DNABornExcitationMode 91 G4cout << "Calling G4DNABornExcitationModel2::Initialise()" << G4endl;
92 } 92 }
93 93
94 if(fParticleDefinition != nullptr && fPartic << 94 if(fParticleDefinition != 0 && fParticleDefinition != particle)
95 { 95 {
96 G4Exception("G4DNABornExcitationModel2::In 96 G4Exception("G4DNABornExcitationModel2::Initialise","em0001",
97 FatalException,"Model already initiali 97 FatalException,"Model already initialized for another particle type.");
98 } 98 }
99 99
100 fParticleDefinition = particle; 100 fParticleDefinition = particle;
101 101
102 std::ostringstream fullFileName; 102 std::ostringstream fullFileName;
103 const char* path = G4FindDataDir("G4LEDATA") << 103 char *path = getenv("G4LEDATA");
104 104
105 if(G4String(path).empty()) << 105 if(G4String(path) == "")
106 { 106 {
107 G4Exception("G4DNABornExcitationModel2::In 107 G4Exception("G4DNABornExcitationModel2::Initialise","G4LEDATA-CHECK",
108 FatalException, "G4LEDATA not defined 108 FatalException, "G4LEDATA not defined in environment variables");
109 } 109 }
110 110
111 fullFileName << path; 111 fullFileName << path;
112 112
113 if(particle->GetParticleName() == "e-") 113 if(particle->GetParticleName() == "e-")
114 { 114 {
115 fullFileName << "/dna/bornExcitation-e.dat 115 fullFileName << "/dna/bornExcitation-e.dat";
116 fLowEnergy = 9*eV; 116 fLowEnergy = 9*eV;
117 fHighEnergy = 1*MeV; 117 fHighEnergy = 1*MeV;
118 } 118 }
119 else if(particle->GetParticleName() == "prot 119 else if(particle->GetParticleName() == "proton")
120 { 120 {
121 fullFileName << "/dna/bornExcitation-p.dat 121 fullFileName << "/dna/bornExcitation-p.dat";
122 fLowEnergy = 500. * keV; 122 fLowEnergy = 500. * keV;
123 fHighEnergy = 100. * MeV; 123 fHighEnergy = 100. * MeV;
124 } 124 }
125 125
126 SetLowEnergyLimit(fLowEnergy); 126 SetLowEnergyLimit(fLowEnergy);
127 SetHighEnergyLimit(fHighEnergy); 127 SetHighEnergyLimit(fHighEnergy);
128 128
129 //G4double scaleFactor = (1.e-22 / 3.343) * 129 //G4double scaleFactor = (1.e-22 / 3.343) * m*m;
130 130
131 fTableData = new G4PhysicsTable(); 131 fTableData = new G4PhysicsTable();
132 fTableData->RetrievePhysicsTable(fullFileNam 132 fTableData->RetrievePhysicsTable(fullFileName.str().c_str(), true);
133 /* << 133 for(size_t level = 0; level<fTableData->size(); ++level)
134 for(std::size_t level = 0; level<fTableData- <<
135 { 134 {
136 //(*fTableData)(level)->ScaleVector(1,scal 135 //(*fTableData)(level)->ScaleVector(1,scaleFactor);
>> 136 (*fTableData)(level)->SetSpline(true);
137 } 137 }
138 */ << 138
139 std::size_t finalBin_i = 2000; << 139 size_t finalBin_i = 2000;
140 G4double E_min = fLowEnergy; 140 G4double E_min = fLowEnergy;
141 G4double E_max = fHighEnergy; 141 G4double E_max = fHighEnergy;
142 fTotalXS = new G4PhysicsLogVector(E_min, E_m << 142 fTotalXS = new G4PhysicsLogVector(E_min, E_max, finalBin_i);
>> 143 fTotalXS->SetSpline(true);
143 G4double energy; 144 G4double energy;
144 G4double finalXS; 145 G4double finalXS;
145 146
146 for(std::size_t energy_i = 0; energy_i < fin << 147 for(size_t energy_i = 0; energy_i < finalBin_i; ++energy_i)
147 { 148 {
148 energy = fTotalXS->Energy(energy_i); 149 energy = fTotalXS->Energy(energy_i);
149 finalXS = 0; 150 finalXS = 0;
150 151
151 for(std::size_t level = 0; level<fTableDat << 152 for(size_t level = 0; level<fTableData->size(); ++level)
152 { 153 {
153 finalXS += (*fTableData)(level)->Value(e 154 finalXS += (*fTableData)(level)->Value(energy);
154 } 155 }
155 fTotalXS->PutValue(energy_i, finalXS); 156 fTotalXS->PutValue(energy_i, finalXS);
156 //G4cout << "energy = " << energy << " " < 157 //G4cout << "energy = " << energy << " " << fTotalXS->Value(energy)
157 // << " " << energy_i << " " << fina 158 // << " " << energy_i << " " << finalXS << G4endl;
158 } 159 }
159 160
160 // for(energy = LowEnergyLimit() ; energy 161 // for(energy = LowEnergyLimit() ; energy < HighEnergyLimit() ; energy += 1*pow(10,log10(energy)))
161 // { 162 // {
162 // G4cout << "energy = " << energy << " 163 // G4cout << "energy = " << energy << " " << fTotalXS->Value(energy) << G4endl;
163 // } 164 // }
164 165
165 if( verboseLevel>0 ) 166 if( verboseLevel>0 )
166 { 167 {
167 G4cout << "Born excitation model is initia 168 G4cout << "Born excitation model is initialized " << G4endl
168 << "Energy range: " 169 << "Energy range: "
169 << LowEnergyLimit() / eV << " eV - " 170 << LowEnergyLimit() / eV << " eV - "
170 << HighEnergyLimit() / keV << " keV for " 171 << HighEnergyLimit() / keV << " keV for "
171 << particle->GetParticleName() 172 << particle->GetParticleName()
172 << G4endl; 173 << G4endl;
173 } 174 }
174 175
175 // Initialize water density pointer 176 // Initialize water density pointer
176 fpMolWaterDensity = G4DNAMolecularMaterial:: 177 fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
177 178
178 if (isInitialised) 179 if (isInitialised)
179 { return;} 180 { return;}
180 fParticleChangeForGamma = GetParticleChangeF 181 fParticleChangeForGamma = GetParticleChangeForGamma();
181 isInitialised = true; 182 isInitialised = true;
182 } 183 }
183 184
184 //....oooOO0OOooo........oooOO0OOooo........oo 185 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
185 186
186 G4double G4DNABornExcitationModel2::CrossSecti 187 G4double G4DNABornExcitationModel2::CrossSectionPerVolume(const G4Material* material,
187 188 const G4ParticleDefinition* particleDefinition,
188 189 G4double ekin,
189 190 G4double,
190 191 G4double)
191 { 192 {
192 if (verboseLevel > 3) 193 if (verboseLevel > 3)
193 { 194 {
194 G4cout << "Calling CrossSectionPerVolume() 195 G4cout << "Calling CrossSectionPerVolume() of G4DNABornExcitationModel2"
195 << G4endl; 196 << G4endl;
196 } 197 }
197 198
198 if(particleDefinition != fParticleDefinition 199 if(particleDefinition != fParticleDefinition) return 0;
199 200
200 // Calculate total cross section for model 201 // Calculate total cross section for model
201 202
202 G4double sigma=0; 203 G4double sigma=0;
203 204
204 G4double waterDensity = (*fpMolWaterDensity) 205 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
205 206
206 if (ekin >= fLowEnergy && ekin <= fHighEnerg 207 if (ekin >= fLowEnergy && ekin <= fHighEnergy)
207 { 208 {
208 sigma = fTotalXS->Value(ekin, fLastBinCall 209 sigma = fTotalXS->Value(ekin, fLastBinCallForFinalXS);
209 210
210 // for(std::size_t i = 0; i < 5; ++i) << 211 // for(size_t i = 0; i < 5; ++i)
211 // sigma += (*fTableData)[i]->Value(ekin); 212 // sigma += (*fTableData)[i]->Value(ekin);
212 213
213 if(sigma == 0) 214 if(sigma == 0)
214 { 215 {
215 G4cerr << "PROBLEM SIGMA = 0 at " << G4B 216 G4cerr << "PROBLEM SIGMA = 0 at " << G4BestUnit(ekin, "Energy")<< G4endl;
216 } 217 }
217 } 218 }
218 219
219 if (verboseLevel > 2) 220 if (verboseLevel > 2)
220 { 221 {
221 G4cout << "_______________________________ 222 G4cout << "__________________________________" << G4endl;
222 G4cout << "G4DNABornExcitationModel2 - XS 223 G4cout << "G4DNABornExcitationModel2 - XS INFO START" << G4endl;
223 G4cout << "Kinetic energy(eV)=" << ekin/eV 224 G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
224 G4cout << "Cross section per water molecul 225 G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
225 G4cout << "Cross section per water molecul 226 G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
226 G4cout << "G4DNABornExcitationModel2 - XS 227 G4cout << "G4DNABornExcitationModel2 - XS INFO END" << G4endl;
227 } 228 }
228 229
229 return sigma*waterDensity; 230 return sigma*waterDensity;
230 } 231 }
231 232
232 //....oooOO0OOooo........oooOO0OOooo........oo 233 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
233 234
234 void G4DNABornExcitationModel2::SampleSecondar 235 void G4DNABornExcitationModel2::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
235 236 const G4MaterialCutsCouple* /*couple*/,
236 237 const G4DynamicParticle* aDynamicParticle,
237 238 G4double,
238 239 G4double)
239 { 240 {
240 241
241 if (verboseLevel > 3) 242 if (verboseLevel > 3)
242 { 243 {
243 G4cout << "Calling SampleSecondaries() of 244 G4cout << "Calling SampleSecondaries() of G4DNABornExcitationModel2"
244 << G4endl; 245 << G4endl;
245 } 246 }
246 247
247 G4double k = aDynamicParticle->GetKineticEne 248 G4double k = aDynamicParticle->GetKineticEnergy();
248 249
249 G4int level = RandomSelect(k); 250 G4int level = RandomSelect(k);
250 G4double excitationEnergy = waterStructure.E 251 G4double excitationEnergy = waterStructure.ExcitationEnergy(level);
251 G4double newEnergy = k - excitationEnergy; 252 G4double newEnergy = k - excitationEnergy;
252 253
253 if (newEnergy > 0) 254 if (newEnergy > 0)
254 { 255 {
255 fParticleChangeForGamma->ProposeMomentumDi 256 fParticleChangeForGamma->ProposeMomentumDirection(aDynamicParticle->GetMomentumDirection());
256 257
257 if (!statCode) fParticleChangeForGamma->Se 258 if (!statCode) fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
258 else fParticleChangeForGamma->SetProposedK 259 else fParticleChangeForGamma->SetProposedKineticEnergy(k);
259 260
260 fParticleChangeForGamma->ProposeLocalEnerg 261 fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
261 } 262 }
262 263
263 const G4Track * theIncomingTrack = fParticle 264 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
264 G4DNAChemistryManager::Instance()->CreateWat 265 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eExcitedMolecule,
265 level, 266 level,
266 theIncomingTrack); 267 theIncomingTrack);
267 } 268 }
268 269
269 //....oooOO0OOooo........oooOO0OOooo........oo 270 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
270 271
271 G4double G4DNABornExcitationModel2::GetPartial 272 G4double G4DNABornExcitationModel2::GetPartialCrossSection(const G4Material*,
272 273 G4int level,
273 274 const G4ParticleDefinition* particle,
274 275 G4double kineticEnergy)
275 { 276 {
276 if (fParticleDefinition != particle) 277 if (fParticleDefinition != particle)
277 { 278 {
278 G4Exception("G4DNABornExcitationModel2::Ge 279 G4Exception("G4DNABornExcitationModel2::GetPartialCrossSection",
279 "bornParticleType", 280 "bornParticleType",
280 FatalException, 281 FatalException,
281 "Model initialized for another 282 "Model initialized for another particle type.");
282 } 283 }
283 284
284 return (*fTableData)(level)->Value(kineticEn 285 return (*fTableData)(level)->Value(kineticEnergy);
285 } 286 }
286 287
287 //....oooOO0OOooo........oooOO0OOooo........oo 288 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
288 289
289 G4int G4DNABornExcitationModel2::RandomSelect( 290 G4int G4DNABornExcitationModel2::RandomSelect(G4double k)
290 { 291 {
291 const std::size_t n(fTableData->size()); << 292 const size_t n(fTableData->size());
292 std::size_t i(n); << 293 size_t i(n);
293 294
294 G4double value = fTotalXS->Value(k, fLastBin 295 G4double value = fTotalXS->Value(k, fLastBinCallForFinalXS);
295 296
296 value *= G4UniformRand(); 297 value *= G4UniformRand();
297 i = n; 298 i = n;
298 299
299 G4double partialXS; 300 G4double partialXS;
300 301
301 while (i > 0) 302 while (i > 0)
302 { 303 {
303 i--; 304 i--;
304 305
305 partialXS = (*fTableData)(i)->Value(k); 306 partialXS = (*fTableData)(i)->Value(k);
306 if (partialXS > value) 307 if (partialXS > value)
307 { 308 {
308 return (G4int)i; << 309 return i;
309 } 310 }
310 value -= partialXS; 311 value -= partialXS;
311 } 312 }
312 313
313 return 0; 314 return 0;
314 } 315 }
315 316
316 317