Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4IonDEDXHandler.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/lowenergy/src/G4IonDEDXHandler.cc (Version 11.3.0) and /processes/electromagnetic/lowenergy/src/G4IonDEDXHandler.cc (Version 10.6.p3)


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