| 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 //
27 //
28 // ===========================================
29 // GEANT4 class source file
30 //
31 // Class: G4IonDEDXHandler
32 //
33 // Author: Anton Lechner (Anton.
34 //
35 // First implementation: 11. 03. 2009
36 //
37 // Modifications: 12. 11 .2009 - Function Buil
38 // methods of st
39 // to interface
40 // Function Upda
41 // ScalingFactor
42 // interface cha
43 // Algorithm (AL
44 //
45 // Class description:
46 // Ion dE/dx table handler.
47 //
48 // Comments:
49 //
50 // ===========================================
51
52 #include <iomanip>
53
54 #include "G4IonDEDXHandler.hh"
55 #include "G4SystemOfUnits.hh"
56 #include "G4VIonDEDXTable.hh"
57 #include "G4VIonDEDXScalingAlgorithm.hh"
58 #include "G4ParticleDefinition.hh"
59 #include "G4Material.hh"
60 #include "G4PhysicsFreeVector.hh"
61 #include "G4Exp.hh"
62
63 // ###########################################
64
65 G4IonDEDXHandler::G4IonDEDXHandler(
66 G4VIonDEDXTable* ionTable,
67 G4VIonDEDXScalingAlgorithm* ionAlgo
68 const G4String& name,
69 G4int maxCacheSize,
70 G4bool splines) :
71 table(ionTable),
72 algorithm(ionAlgorithm),
73 tableName(name),
74 useSplines(splines),
75 maxCacheEntries(maxCacheSize) {
76
77 if(table == nullptr) {
78 G4cerr << "G4IonDEDXHandler::G4IonDEDXHan
79 << " Pointer to G4VIonDEDXTable ob
80 << G4endl;
81 }
82
83 if(algorithm == nullptr) {
84 G4cerr << "G4IonDEDXHandler::G4IonDEDXHan
85 << " Pointer to G4VIonDEDXScalingA
86 << G4endl;
87 }
88
89 if(maxCacheEntries <= 0) {
90 G4cerr << "G4IonDEDXHandler::G4IonDEDXHan
91 << " Cache size <=0. Resetting to
92 << G4endl;
93 maxCacheEntries = 5;
94 }
95 }
96
97 // ###########################################
98
99 G4IonDEDXHandler::~G4IonDEDXHandler() {
100
101 ClearCache();
102
103 // All stopping power vectors built accordin
104 // are deleted. All other stopping power vec
105 // deleted by their creator class (sub-class
106 // DEDXTableBraggRule::iterator iter = stopp
107 // DEDXTableBraggRule::iterator iter_end = s
108
109 // for(;iter != iter_end; iter++) delete it
110 stoppingPowerTableBragg.clear();
111
112 stoppingPowerTable.clear();
113
114 if(table != nullptr)
115 delete table;
116 if(algorithm != nullptr)
117 delete algorithm;
118 }
119
120 // ###########################################
121
122 G4bool G4IonDEDXHandler::IsApplicable(
123 const G4ParticleDefinition* p
124 const G4Material* material) {
125
126 G4bool isApplicable = true;
127
128 if(table == nullptr || algorithm == nullptr)
129 isApplicable = false;
130 }
131 else {
132
133 G4int atomicNumberIon = particle -> GetAt
134 G4int atomicNumberBase =
135 algorithm -> AtomicNumberBaseI
136
137 G4IonKey key = std::make_pair(atomicNumbe
138
139 DEDXTable::iterator iter = stoppingPowerT
140 if(iter == stoppingPowerTable.end()) isAp
141 }
142
143 return isApplicable;
144 }
145
146 // ###########################################
147
148 G4double G4IonDEDXHandler::GetDEDX(
149 const G4ParticleDefinition* p
150 const G4Material* material,
151 G4double kineticEnergy) {
152
153 G4double dedx = 0.0;
154
155 G4CacheValue value = GetCacheValue(particle,
156
157 if(kineticEnergy <= 0.0) dedx = 0.0;
158 else if(value.dedxVector != 0) {
159
160 G4bool b;
161 G4double factor = value.density;
162
163 factor *= algorithm -> ScalingFactorDEDX(
164 m
165 k
166 G4double scaledKineticEnergy = kineticEne
167
168 if(scaledKineticEnergy < value.lowerEnerg
169
170 factor *= std::sqrt(scaledKineticEnerg
171 scaledKineticEnergy = value.lowerEnerg
172 }
173
174 dedx = factor * value.dedxVector -> GetVa
175
176 if(dedx < 0.0) dedx = 0.0;
177 }
178 else dedx = 0.0;
179
180 #ifdef PRINT_DEBUG
181 G4cout << "G4IonDEDXHandler::GetDEDX() E
182 << kineticEnergy / MeV << " MeV *
183 << value.energyScaling << " = "
184 << kineticEnergy * value.energySca
185 << " MeV, dE/dx = " << dedx / MeV
186 << ", material = " << material ->
187 << G4endl;
188 #endif
189
190 return dedx;
191 }
192
193 // ###########################################
194
195 G4bool G4IonDEDXHandler::BuildDEDXTable(
196 const G4ParticleDefinition* p
197 const G4Material* material) {
198
199 G4int atomicNumberIon = particle -> GetAtomi
200
201 G4bool isApplicable = BuildDEDXTable(atomicN
202
203 return isApplicable;
204 }
205
206
207 // ###########################################
208
209 G4bool G4IonDEDXHandler::BuildDEDXTable(
210 G4int atomicNumberIon,
211 const G4Material* material) {
212
213 G4bool isApplicable = true;
214
215 if(table == 0 || algorithm == 0) {
216 isApplicable = false;
217 return isApplicable;
218 }
219
220 G4int atomicNumberBase =
221 algorithm -> AtomicNumberBaseI
222
223 // Checking if vector is already built, and
224 // the case
225 G4IonKey key = std::make_pair(atomicNumberBa
226
227 auto iter = stoppingPowerTable.find(key);
228 if(iter != stoppingPowerTable.end()) return
229
230 // Checking if table contains stopping power
231 // or chemical formula
232 const G4String& chemFormula = material -> Ge
233 const G4String& materialName = material -> G
234
235 isApplicable = table -> BuildPhysicsVector(a
236
237 if(isApplicable) {
238 stoppingPowerTable[key] =
239 table -> GetPhysicsVector(atomic
240 return isApplicable;
241 }
242
243 isApplicable = table -> BuildPhysicsVector(a
244 if(isApplicable) {
245 stoppingPowerTable[key] =
246 table -> GetPhysicsVector(atomic
247 return isApplicable;
248 }
249
250 // Building the stopping power vector based
251 const G4ElementVector* elementVector = mater
252
253 std::vector<G4PhysicsVector*> dEdxTable;
254
255 size_t nmbElements = material -> GetNumberOf
256
257 for(size_t i = 0; i < nmbElements; i++) {
258
259 G4int atomicNumberMat = G4int((*elementV
260
261 isApplicable = table -> BuildPhysicsVect
262
263 if(isApplicable) {
264
265 G4PhysicsVector* dEdx =
266 table -> GetPhysicsVector(at
267 dEdxTable.push_back(dEdx);
268 }
269 else {
270
271 dEdxTable.clear();
272 break;
273 }
274 }
275
276 if(isApplicable) {
277
278 if(dEdxTable.size() > 0) {
279
280 size_t nmbdEdxBins = dEdxTable[0] -> G
281 G4double lowerEdge = dEdxTable[0] -> G
282 G4double upperEdge = dEdxTable[0] -> G
283
284 G4PhysicsFreeVector* dEdxBragg =
285 new G4PhysicsFreeVector(nmbdEdxBins,
286 lowerEdge,
287 upperEdge,
288 useSplines);
289
290 const G4double* massFractionVector = m
291
292 G4bool b;
293 for(size_t j = 0; j < nmbdEdxBins; j++
294
295 G4double edge = dEdxTable[0] -> Ge
296
297 G4double value = 0.0;
298 for(size_t i = 0; i < nmbElements; i++
299
300 value += (dEdxTable[i] -> GetV
301
302 }
303
304 dEdxBragg -> PutValues(j, edge, va
305 }
306 if (useSplines)
307 dEdxBragg -> FillSecondDerivatives();
308
309 #ifdef PRINT_DEBUG
310 G4cout << "G4IonDEDXHandler::BuildPhys
311 << atomicNumberBase << " in "
312 << material -> GetName()
313 << G4endl;
314
315 G4cout << *dEdxBragg;
316 #endif
317
318 stoppingPowerTable[key] = dEdxBragg;
319 stoppingPowerTableBragg[key] = dEdxBragg;
320 }
321 }
322
323 ClearCache();
324
325 return isApplicable;
326 }
327
328 // ###########################################
329
330 G4CacheValue G4IonDEDXHandler::UpdateCacheValu
331 const G4ParticleDefinition* part
332 const G4Material* material) {
333
334 G4CacheValue value;
335
336 G4int atomicNumberIon = particle -> GetAtomi
337 G4int atomicNumberBase =
338 algorithm -> AtomicNumberBaseI
339
340 G4IonKey key = std::make_pair(atomicNumberBa
341
342 DEDXTable::iterator iter = stoppingPowerTabl
343
344 if(iter != stoppingPowerTable.end()) {
345 value.dedxVector = iter -> second;
346
347 G4double nmbNucleons = G4double(particle
348 value.energyScaling =
349 algorithm -> ScalingFactorEnergy(pa
350
351 size_t nmbdEdxBins = value.dedxVector ->
352 value.lowerEnergyEdge = value.dedxVector
353
354 value.upperEnergyEdge =
355 value.dedxVector -> Get
356 value.density = material -> GetDensity();
357 }
358 else {
359 value.dedxVector = 0;
360 value.energyScaling = 0.0;
361 value.lowerEnergyEdge = 0.0;
362 value.upperEnergyEdge = 0.0;
363 value.density = 0.0;
364 }
365
366 #ifdef PRINT_DEBUG
367 G4cout << "G4IonDEDXHandler::UpdateCacheValu
368 << particle -> GetParticleName() << "
369 << material -> GetName()
370 << G4endl;
371 #endif
372
373 return value;
374 }
375
376 // ###########################################
377
378 G4CacheValue G4IonDEDXHandler::GetCacheValue(
379 const G4ParticleDefinition* part
380 const G4Material* material) {
381
382 G4CacheKey key = std::make_pair(particle, ma
383
384 G4CacheEntry entry;
385 CacheEntryList::iterator* pointerIter =
386 (CacheEntryList::iterator*)
387
388 if(!pointerIter) {
389 entry.value = UpdateCacheValue(particle,
390
391 entry.key = key;
392 cacheEntries.push_front(entry);
393
394 CacheEntryList::iterator* pointerIter1 =
395 new CacheEn
396 *pointerIter1 = cacheEntries.begin();
397 cacheKeyPointers[key] = pointerIter1;
398
399 if(G4int(cacheEntries.size()) > maxCache
400
401 G4CacheEntry lastEntry = cacheEntries.bac
402
403 void* pointerIter2 = cacheKeyPointers
404 CacheEntryList::iterator* listPointer
405 (CacheEntryList::ite
406
407 cacheEntries.erase(*listPointerIter);
408
409 delete listPointerIter;
410 cacheKeyPointers.erase(lastEntry.key)
411 }
412 }
413 else {
414 entry = *(*pointerIter);
415 // Cache entries are currently not re-or
416 // Uncomment for activating re-ordering:
417 // cacheEntries.erase(*pointerIter)
418 // cacheEntries.push_front(entry);
419 // *pointerIter = cacheEntries.begi
420 }
421 return entry.value;
422 }
423
424 // ###########################################
425
426 void G4IonDEDXHandler::ClearCache() {
427
428 CacheIterPointerMap::iterator iter = cacheKe
429 CacheIterPointerMap::iterator iter_end = cac
430
431 for(;iter != iter_end; iter++) {
432 void* pointerIter = iter -> second;
433 CacheEntryList::iterator* listPointerIte
434 (CacheEntryList::ite
435
436 delete listPointerIter;
437 }
438
439 cacheEntries.clear();
440 cacheKeyPointers.clear();
441 }
442
443 // ###########################################
444
445 void G4IonDEDXHandler::PrintDEDXTable(
446 const G4ParticleDefinition*
447 const G4Material* material,
448 G4double lowerBoundary,
449 G4double upperBoundary,
450 G4int nmbBins,
451 G4bool logScaleEnergy) {
452
453 G4double atomicMassNumber = particle -> GetA
454 G4double materialDensity = material -> GetDe
455
456 G4cout << "# dE/dx table for " << particle -
457 << " in material " << material -> Get
458 << " of density " << materialDensity
459 << " g/cm3"
460 << G4endl
461 << "# Projectile mass number A1 = " <
462 << G4endl
463 << "# Energy range (per nucleon) of t
464 << GetLowerEnergyEdge(particle, mater
465 << " - "
466 << GetUpperEnergyEdge(particle, mater
467 << " MeV"
468 << G4endl
469 << "# -------------------------------
470 << G4endl;
471 G4cout << "#"
472 << std::setw(13) << std::right << "E"
473 << std::setw(14) << "E/A1"
474 << std::setw(14) << "dE/dx"
475 << std::setw(14) << "1/rho*dE/dx"
476 << G4endl;
477 G4cout << "#"
478 << std::setw(13) << std::right << "(M
479 << std::setw(14) << "(MeV)"
480 << std::setw(14) << "(MeV/cm)"
481 << std::setw(14) << "(MeV*cm2/mg)"
482 << G4endl
483 << "# -------------------------------
484 << G4endl;
485
486 //G4CacheValue value = GetCacheValue(particl
487
488 G4double energyLowerBoundary = lowerBoundary
489 G4double energyUpperBoundary = upperBoundary
490
491 if(logScaleEnergy) {
492
493 energyLowerBoundary = std::log(energyLowe
494 energyUpperBoundary = std::log(energyUppe
495 }
496
497 G4double deltaEnergy = (energyUpperBoundary
498
499
500 G4cout.precision(6);
501 for(int i = 0; i < nmbBins + 1; i++) {
502
503 G4double energy = energyLowerBoundary +
504 if(logScaleEnergy) energy = G4Exp(energy
505
506 G4double loss = GetDEDX(particle, materi
507
508 G4cout << std::setw(14) << std::right <<
509 << std::setw(14) << energy / atom
510 << std::setw(14) << loss / MeV *
511 << std::setw(14) << loss / materi
512 << G4endl;
513 }
514 }
515
516 // ###########################################
517
518 G4double G4IonDEDXHandler::GetLowerEnergyEdge(
519 const G4ParticleDefinition* p
520 const G4Material* material) {
521
522 G4double edge = 0.0;
523
524 G4CacheValue value = GetCacheValue(particle,
525
526 if(value.energyScaling > 0)
527 edge = value.lowerEnergyEdge / value
528
529 return edge;
530 }
531
532 // ###########################################
533
534 G4double G4IonDEDXHandler::GetUpperEnergyEdge(
535 const G4ParticleDefinition* p
536 const G4Material* material) {
537
538 G4double edge = 0.0;
539
540 G4CacheValue value = GetCacheValue(particle,
541
542 if(value.energyScaling > 0)
543 edge = value.upperEnergyEdge / value
544
545 return edge;
546 }
547
548 // ###########################################
549
550 G4String G4IonDEDXHandler::GetName() {
551
552 return tableName;
553 }
554
555 // ###########################################
556