| 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 // ------------------------------------------- 28 // -----------------------------------------------------------
29 // GEANT 4 class implementation file 29 // GEANT 4 class implementation file
30 // 30 //
31 // History: based on object model of 31 // History: based on object model of
32 // 2nd December 1995, G.Cosmo 32 // 2nd December 1995, G.Cosmo
33 // ---------- G4hEnergyLoss physics proce 33 // ---------- G4hEnergyLoss physics process -----------
34 // by Laszlo Urban, 30 May 1997 34 // by Laszlo Urban, 30 May 1997
35 // 35 //
36 // ******************************************* 36 // **************************************************************
37 // It is the first implementation of the NEW U 37 // It is the first implementation of the NEW UNIFIED ENERGY LOSS PROCESS.
38 // It calculates the energy loss of charged ha 38 // It calculates the energy loss of charged hadrons.
39 // ******************************************* 39 // **************************************************************
40 // 40 //
41 // 7/10/98 bug fixes + some cleanup , L.Urb 41 // 7/10/98 bug fixes + some cleanup , L.Urban
42 // 22/10/98 cleanup , L.Urban 42 // 22/10/98 cleanup , L.Urban
43 // 07/12/98 works for ions as well+ bug corr 43 // 07/12/98 works for ions as well+ bug corrected, L.Urban
44 // 02/02/99 several bugs fixed, L.Urban 44 // 02/02/99 several bugs fixed, L.Urban
45 // 31/03/00 rename to lowenergy as G4hRDEner 45 // 31/03/00 rename to lowenergy as G4hRDEnergyLoss.cc V.Ivanchenko
46 // 05/11/00 new method to calculate particle 46 // 05/11/00 new method to calculate particle ranges
47 // 10/05/01 V.Ivanchenko Clean up againist L 47 // 10/05/01 V.Ivanchenko Clean up againist Linux compilation with -Wall
48 // 23/11/01 V.Ivanchenko Move static member- 48 // 23/11/01 V.Ivanchenko Move static member-functions from header to source
49 // 28/10/02 V.Ivanchenko Optimal binning for 49 // 28/10/02 V.Ivanchenko Optimal binning for dE/dx
50 // 21/01/03 V.Ivanchenko Cut per region 50 // 21/01/03 V.Ivanchenko Cut per region
51 // 23/01/03 V.Ivanchenko Fix in table build 51 // 23/01/03 V.Ivanchenko Fix in table build
52 52
53 // 31 Jul 2008 MGP Short term supply of en 53 // 31 Jul 2008 MGP Short term supply of energy loss of hadrons through clone of
54 // former G4hLowEnergyLoss 54 // former G4hLowEnergyLoss (with some initial cleaning)
55 // To be replaced by rewor 55 // To be replaced by reworked class to deal with condensed/discrete
56 // issues properly 56 // issues properly
57 57
58 // 25 Nov 2010 MGP Added some protections 58 // 25 Nov 2010 MGP Added some protections for FPE (mostly division by 0)
59 // The whole energy loss d 59 // The whole energy loss domain would profit from a design iteration
60 60
61 // 20 Jun 2011 MGP Corrected some compilat 61 // 20 Jun 2011 MGP Corrected some compilation warnings. The whole class will be heavily refactored anyway.
62 62
63 // ------------------------------------------- 63 // --------------------------------------------------------------
64 64
65 #include "G4hRDEnergyLoss.hh" 65 #include "G4hRDEnergyLoss.hh"
66 #include "G4PhysicalConstants.hh" 66 #include "G4PhysicalConstants.hh"
67 #include "G4SystemOfUnits.hh" 67 #include "G4SystemOfUnits.hh"
68 #include "G4EnergyLossTables.hh" 68 #include "G4EnergyLossTables.hh"
69 #include "G4Poisson.hh" 69 #include "G4Poisson.hh"
70 #include "G4ProductionCutsTable.hh" 70 #include "G4ProductionCutsTable.hh"
71 71
72 72
73 // Initialisation of static members ********** 73 // Initialisation of static members ******************************************
74 // contributing processes : ion.loss ->NumberO 74 // contributing processes : ion.loss ->NumberOfProcesses is initialized
75 // to 1 . YOU DO NOT HAVE TO CHANGE this var 75 // to 1 . YOU DO NOT HAVE TO CHANGE this variable for a 'normal' run.
76 // You have to change NumberOfProcesses 76 // You have to change NumberOfProcesses
77 // if you invent a new process contributing to 77 // if you invent a new process contributing to the cont. energy loss,
78 // NumberOfProcesses should be 2 in this cas 78 // NumberOfProcesses should be 2 in this case,
79 // or for debugging purposes. 79 // or for debugging purposes.
80 // The NumberOfProcesses data member can be c 80 // The NumberOfProcesses data member can be changed using the (public static)
81 // functions Get/Set/Plus/MinusNumberOfProces 81 // functions Get/Set/Plus/MinusNumberOfProcesses (see G4hRDEnergyLoss.hh)
82 82
83 83
84 // The following vectors have a fixed dimensio 84 // The following vectors have a fixed dimension this means that if they are
85 // filled in with more than 100 elements will 85 // filled in with more than 100 elements will corrupt the memory.
86 G4ThreadLocal G4int G4hRDEnergyLoss::NumberOfP 86 G4ThreadLocal G4int G4hRDEnergyLoss::NumberOfProcesses = 1 ;
87 87
88 G4ThreadLocal G4int G4hRDEnergyLoss 88 G4ThreadLocal G4int G4hRDEnergyLoss::CounterOfProcess = 0 ;
89 G4ThreadLocal G4PhysicsTable** G4hRDEnergyLoss 89 G4ThreadLocal G4PhysicsTable** G4hRDEnergyLoss::RecorderOfProcess = 0 ;
90 90
91 G4ThreadLocal G4int G4hRDEnergyLoss 91 G4ThreadLocal G4int G4hRDEnergyLoss::CounterOfpProcess = 0 ;
92 G4ThreadLocal G4PhysicsTable** G4hRDEnergyLoss 92 G4ThreadLocal G4PhysicsTable** G4hRDEnergyLoss::RecorderOfpProcess = 0 ;
93 93
94 G4ThreadLocal G4int G4hRDEnergyLoss 94 G4ThreadLocal G4int G4hRDEnergyLoss::CounterOfpbarProcess = 0 ;
95 G4ThreadLocal G4PhysicsTable** G4hRDEnergyLoss 95 G4ThreadLocal G4PhysicsTable** G4hRDEnergyLoss::RecorderOfpbarProcess = 0 ;
96 96
97 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss: 97 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss::theDEDXpTable = 0 ;
98 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss: 98 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss::theDEDXpbarTable = 0 ;
99 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss: 99 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss::theRangepTable = 0 ;
100 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss: 100 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss::theRangepbarTable = 0 ;
101 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss: 101 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss::theInverseRangepTable = 0 ;
102 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss: 102 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss::theInverseRangepbarTable = 0 ;
103 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss: 103 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss::theLabTimepTable = 0 ;
104 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss: 104 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss::theLabTimepbarTable = 0 ;
105 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss: 105 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss::theProperTimepTable = 0 ;
106 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss: 106 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss::theProperTimepbarTable = 0 ;
107 107
108 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss: 108 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss::thepRangeCoeffATable = 0 ;
109 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss: 109 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss::thepRangeCoeffBTable = 0 ;
110 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss: 110 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss::thepRangeCoeffCTable = 0 ;
111 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss: 111 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss::thepbarRangeCoeffATable = 0 ;
112 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss: 112 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss::thepbarRangeCoeffBTable = 0 ;
113 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss: 113 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss::thepbarRangeCoeffCTable = 0 ;
114 114
115 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss: 115 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss::theDEDXTable = 0 ;
116 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss: 116 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss::theRangeTable = 0 ;
117 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss: 117 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss::theInverseRangeTable = 0 ;
118 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss: 118 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss::theLabTimeTable = 0 ;
119 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss: 119 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss::theProperTimeTable = 0 ;
120 120
121 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss: 121 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss::theRangeCoeffATable = 0 ;
122 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss: 122 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss::theRangeCoeffBTable = 0 ;
123 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss: 123 G4ThreadLocal G4PhysicsTable* G4hRDEnergyLoss::theRangeCoeffCTable = 0 ;
124 124
125 //const G4Proton* G4hRDEnergyLoss::theProton=G 125 //const G4Proton* G4hRDEnergyLoss::theProton=G4Proton::Proton() ;
126 //const G4AntiProton* G4hRDEnergyLoss::theAnti 126 //const G4AntiProton* G4hRDEnergyLoss::theAntiProton=G4AntiProton::AntiProton() ;
127 127
128 G4ThreadLocal G4double G4hRDEnergyLoss::Partic 128 G4ThreadLocal G4double G4hRDEnergyLoss::ParticleMass;
129 G4ThreadLocal G4double G4hRDEnergyLoss::ptable 129 G4ThreadLocal G4double G4hRDEnergyLoss::ptableElectronCutInRange = 0.0 ;
130 G4ThreadLocal G4double G4hRDEnergyLoss::pbarta 130 G4ThreadLocal G4double G4hRDEnergyLoss::pbartableElectronCutInRange = 0.0 ;
131 131
132 G4ThreadLocal G4double G4hRDEnergyLoss::Mass, 132 G4ThreadLocal G4double G4hRDEnergyLoss::Mass,
133 G4hRDEnergyLoss::taulow 133 G4hRDEnergyLoss::taulow,
134 G4hRDEnergyLoss::tauhig 134 G4hRDEnergyLoss::tauhigh,
135 G4hRDEnergyLoss::ltaulo 135 G4hRDEnergyLoss::ltaulow,
136 G4hRDEnergyLoss::ltauhi 136 G4hRDEnergyLoss::ltauhigh;
137 137
138 G4ThreadLocal G4double G4hRDEnergyLoss::dRover 138 G4ThreadLocal G4double G4hRDEnergyLoss::dRoverRange = 0.20 ;
139 G4ThreadLocal G4double G4hRDEnergyLoss::finalR 139 G4ThreadLocal G4double G4hRDEnergyLoss::finalRange = 0.2; // 200.*micrometer
140 140
141 G4ThreadLocal G4double G4hRDEnergyLoss::c1lim 141 G4ThreadLocal G4double G4hRDEnergyLoss::c1lim = 0.20 ;
142 G4ThreadLocal G4double G4hRDEnergyLoss::c2lim 142 G4ThreadLocal G4double G4hRDEnergyLoss::c2lim = 0.32 ;
143 // 2.*(1.-(0.20))*(200.*micrometer) 143 // 2.*(1.-(0.20))*(200.*micrometer)
144 G4ThreadLocal G4double G4hRDEnergyLoss::c3lim 144 G4ThreadLocal G4double G4hRDEnergyLoss::c3lim = -0.032 ;
145 // -(1.-(0.20))*(200.*micrometer)*(200.*micro 145 // -(1.-(0.20))*(200.*micrometer)*(200.*micrometer)
146 146
147 G4ThreadLocal G4double G4hRDEnergyLoss::Charge 147 G4ThreadLocal G4double G4hRDEnergyLoss::Charge ;
148 148
149 G4ThreadLocal G4bool G4hRDEnergyLoss::rndmSt 149 G4ThreadLocal G4bool G4hRDEnergyLoss::rndmStepFlag = false ;
150 G4ThreadLocal G4bool G4hRDEnergyLoss::Enloss 150 G4ThreadLocal G4bool G4hRDEnergyLoss::EnlossFlucFlag = true ;
151 151
152 G4ThreadLocal G4double G4hRDEnergyLoss::Lowest 152 G4ThreadLocal G4double G4hRDEnergyLoss::LowestKineticEnergy = 1e-05 ; // 10.*eV;
153 G4ThreadLocal G4double G4hRDEnergyLoss::Highes 153 G4ThreadLocal G4double G4hRDEnergyLoss::HighestKineticEnergy= 1.e5; // 100.*GeV;
154 G4ThreadLocal G4int G4hRDEnergyLoss::TotBin 154 G4ThreadLocal G4int G4hRDEnergyLoss::TotBin = 360;
155 G4ThreadLocal G4double G4hRDEnergyLoss::RTable 155 G4ThreadLocal G4double G4hRDEnergyLoss::RTable, G4hRDEnergyLoss::LOGRTable;
156 156
157 //-------------------------------------------- 157 //--------------------------------------------------------------------------------
158 158
159 // constructor and destructor 159 // constructor and destructor
160 160
161 G4hRDEnergyLoss::G4hRDEnergyLoss(const G4Strin 161 G4hRDEnergyLoss::G4hRDEnergyLoss(const G4String& processName)
162 : G4VContinuousDiscreteProcess (processName) 162 : G4VContinuousDiscreteProcess (processName),
163 MaxExcitationNumber (1.e6), 163 MaxExcitationNumber (1.e6),
164 probLimFluct (0.01), 164 probLimFluct (0.01),
165 nmaxDirectFluct (100), 165 nmaxDirectFluct (100),
166 nmaxCont1(4), 166 nmaxCont1(4),
167 nmaxCont2(16), 167 nmaxCont2(16),
168 theLossTable(0), 168 theLossTable(0),
169 linLossLimit(0.05), 169 linLossLimit(0.05),
170 MinKineticEnergy(0.0) 170 MinKineticEnergy(0.0)
171 { 171 {
172 if (!RecorderOfpbarProcess) RecorderOfpbarPr 172 if (!RecorderOfpbarProcess) RecorderOfpbarProcess = new G4PhysicsTable*[100];
173 if (!RecorderOfpProcess) RecorderOfpProcess 173 if (!RecorderOfpProcess) RecorderOfpProcess = new G4PhysicsTable*[100];
174 if (!RecorderOfProcess) RecorderOfProcess = 174 if (!RecorderOfProcess) RecorderOfProcess = new G4PhysicsTable*[100];
175 } 175 }
176 176
177 //-------------------------------------------- 177 //--------------------------------------------------------------------------------
178 178
179 G4hRDEnergyLoss::~G4hRDEnergyLoss() 179 G4hRDEnergyLoss::~G4hRDEnergyLoss()
180 { 180 {
181 if(theLossTable) 181 if(theLossTable)
182 { 182 {
183 theLossTable->clearAndDestroy(); 183 theLossTable->clearAndDestroy();
184 delete theLossTable; 184 delete theLossTable;
185 } 185 }
186 } 186 }
187 187
188 //-------------------------------------------- 188 //--------------------------------------------------------------------------------
189 189
190 G4int G4hRDEnergyLoss::GetNumberOfProcesses() 190 G4int G4hRDEnergyLoss::GetNumberOfProcesses()
191 { 191 {
192 return NumberOfProcesses; 192 return NumberOfProcesses;
193 } 193 }
194 194
195 //-------------------------------------------- 195 //--------------------------------------------------------------------------------
196 196
197 void G4hRDEnergyLoss::SetNumberOfProcesses(G4i 197 void G4hRDEnergyLoss::SetNumberOfProcesses(G4int number)
198 { 198 {
199 NumberOfProcesses=number; 199 NumberOfProcesses=number;
200 } 200 }
201 201
202 //-------------------------------------------- 202 //--------------------------------------------------------------------------------
203 203
204 void G4hRDEnergyLoss::PlusNumberOfProcesses() 204 void G4hRDEnergyLoss::PlusNumberOfProcesses()
205 { 205 {
206 NumberOfProcesses++; 206 NumberOfProcesses++;
207 } 207 }
208 208
209 //-------------------------------------------- 209 //--------------------------------------------------------------------------------
210 210
211 void G4hRDEnergyLoss::MinusNumberOfProcesses() 211 void G4hRDEnergyLoss::MinusNumberOfProcesses()
212 { 212 {
213 NumberOfProcesses--; 213 NumberOfProcesses--;
214 } 214 }
215 215
216 //-------------------------------------------- 216 //--------------------------------------------------------------------------------
217 217
218 void G4hRDEnergyLoss::SetdRoverRange(G4double 218 void G4hRDEnergyLoss::SetdRoverRange(G4double value)
219 { 219 {
220 dRoverRange = value; 220 dRoverRange = value;
221 } 221 }
222 222
223 //-------------------------------------------- 223 //--------------------------------------------------------------------------------
224 224
225 void G4hRDEnergyLoss::SetRndmStep (G4bool va 225 void G4hRDEnergyLoss::SetRndmStep (G4bool value)
226 { 226 {
227 rndmStepFlag = value; 227 rndmStepFlag = value;
228 } 228 }
229 229
230 //-------------------------------------------- 230 //--------------------------------------------------------------------------------
231 231
232 void G4hRDEnergyLoss::SetEnlossFluc (G4bool va 232 void G4hRDEnergyLoss::SetEnlossFluc (G4bool value)
233 { 233 {
234 EnlossFlucFlag = value; 234 EnlossFlucFlag = value;
235 } 235 }
236 236
237 //-------------------------------------------- 237 //--------------------------------------------------------------------------------
238 238
239 void G4hRDEnergyLoss::SetStepFunction (G4doubl 239 void G4hRDEnergyLoss::SetStepFunction (G4double c1, G4double c2)
240 { 240 {
241 dRoverRange = c1; 241 dRoverRange = c1;
242 finalRange = c2; 242 finalRange = c2;
243 c1lim=dRoverRange; 243 c1lim=dRoverRange;
244 c2lim=2.*(1-dRoverRange)*finalRange; 244 c2lim=2.*(1-dRoverRange)*finalRange;
245 c3lim=-(1.-dRoverRange)*finalRange*finalRang 245 c3lim=-(1.-dRoverRange)*finalRange*finalRange;
246 } 246 }
247 247
248 //-------------------------------------------- 248 //--------------------------------------------------------------------------------
249 249
250 void G4hRDEnergyLoss::BuildDEDXTable(const G4P 250 void G4hRDEnergyLoss::BuildDEDXTable(const G4ParticleDefinition& aParticleType)
251 { 251 {
252 // calculate data members TotBin,LOGRTable, 252 // calculate data members TotBin,LOGRTable,RTable first
253 253
254 if (!RecorderOfpbarProcess) RecorderOfpbarPr 254 if (!RecorderOfpbarProcess) RecorderOfpbarProcess = new G4PhysicsTable*[100];
255 if (!RecorderOfpProcess) RecorderOfpProcess 255 if (!RecorderOfpProcess) RecorderOfpProcess = new G4PhysicsTable*[100];
256 if (!RecorderOfProcess) RecorderOfProcess = 256 if (!RecorderOfProcess) RecorderOfProcess = new G4PhysicsTable*[100];
257 257
258 const G4ProductionCutsTable* theCoupleTable= 258 const G4ProductionCutsTable* theCoupleTable=
259 G4ProductionCutsTable::GetProductionCutsTa 259 G4ProductionCutsTable::GetProductionCutsTable();
260 std::size_t numOfCouples = theCoupleTable->G 260 std::size_t numOfCouples = theCoupleTable->GetTableSize();
261 261
262 // create/fill proton or antiproton tables d 262 // create/fill proton or antiproton tables depending on the charge
263 Charge = aParticleType.GetPDGCharge()/eplus; 263 Charge = aParticleType.GetPDGCharge()/eplus;
264 ParticleMass = aParticleType.GetPDGMass() ; 264 ParticleMass = aParticleType.GetPDGMass() ;
265 265
266 if (Charge>0.) {theDEDXTable= theDEDXpTable; 266 if (Charge>0.) {theDEDXTable= theDEDXpTable;}
267 else {theDEDXTable= theDEDXpbarTab 267 else {theDEDXTable= theDEDXpbarTable;}
268 268
269 if( ((Charge>0.) && (theDEDXTable==0)) || 269 if( ((Charge>0.) && (theDEDXTable==0)) ||
270 ((Charge<0.) && (theDEDXTable==0)) 270 ((Charge<0.) && (theDEDXTable==0))
271 ) 271 )
272 { 272 {
273 273
274 // Build energy loss table as a sum of t 274 // Build energy loss table as a sum of the energy loss due to the
275 // different processes. 275 // different processes.
276 if( Charge >0.) 276 if( Charge >0.)
277 { 277 {
278 RecorderOfProcess=RecorderOfpProcess; 278 RecorderOfProcess=RecorderOfpProcess;
279 CounterOfProcess=CounterOfpProcess; 279 CounterOfProcess=CounterOfpProcess;
280 280
281 if(CounterOfProcess == NumberOfProcesses) 281 if(CounterOfProcess == NumberOfProcesses)
282 { 282 {
283 if(theDEDXpTable) 283 if(theDEDXpTable)
284 { theDEDXpTable->clearAndDestroy(); 284 { theDEDXpTable->clearAndDestroy();
285 delete theDEDXpTable; } 285 delete theDEDXpTable; }
286 theDEDXpTable = new G4PhysicsTable(num 286 theDEDXpTable = new G4PhysicsTable(numOfCouples);
287 theDEDXTable = theDEDXpTable; 287 theDEDXTable = theDEDXpTable;
288 } 288 }
289 } 289 }
290 else 290 else
291 { 291 {
292 RecorderOfProcess=RecorderOfpbarProcess; 292 RecorderOfProcess=RecorderOfpbarProcess;
293 CounterOfProcess=CounterOfpbarProcess; 293 CounterOfProcess=CounterOfpbarProcess;
294 294
295 if(CounterOfProcess == NumberOfProcesses) 295 if(CounterOfProcess == NumberOfProcesses)
296 { 296 {
297 if(theDEDXpbarTable) 297 if(theDEDXpbarTable)
298 { theDEDXpbarTable->clearAndDestroy(); 298 { theDEDXpbarTable->clearAndDestroy();
299 delete theDEDXpbarTable; } 299 delete theDEDXpbarTable; }
300 theDEDXpbarTable = new G4PhysicsTable( 300 theDEDXpbarTable = new G4PhysicsTable(numOfCouples);
301 theDEDXTable = theDEDXpbarTable; 301 theDEDXTable = theDEDXpbarTable;
302 } 302 }
303 } 303 }
304 304
305 if(CounterOfProcess == NumberOfProcesses 305 if(CounterOfProcess == NumberOfProcesses)
306 { 306 {
307 // loop for materials 307 // loop for materials
308 G4double LowEdgeEnergy , Value ; 308 G4double LowEdgeEnergy , Value ;
309 G4bool isOutRange ; 309 G4bool isOutRange ;
310 G4PhysicsTable* pointer ; 310 G4PhysicsTable* pointer ;
311 311
312 for (std::size_t J=0; J<numOfCouples; ++J) 312 for (std::size_t J=0; J<numOfCouples; ++J)
313 { 313 {
314 // create physics vector and fill it 314 // create physics vector and fill it
315 G4PhysicsLogVector* aVector = 315 G4PhysicsLogVector* aVector =
316 new G4PhysicsLogVector(LowestK 316 new G4PhysicsLogVector(LowestKineticEnergy,
317 Highest 317 HighestKineticEnergy, TotBin);
318 318
319 // loop for the kinetic energy 319 // loop for the kinetic energy
320 for (G4int i=0; i<TotBin; i++) 320 for (G4int i=0; i<TotBin; i++)
321 { 321 {
322 LowEdgeEnergy = aVector->GetLowEdgeEnerg 322 LowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ;
323 Value = 0. ; 323 Value = 0. ;
324 324
325 // loop for the contributing processes 325 // loop for the contributing processes
326 for (G4int process=0; process < NumberOf 326 for (G4int process=0; process < NumberOfProcesses; process++)
327 { 327 {
328 pointer= RecorderOfProcess[process]; 328 pointer= RecorderOfProcess[process];
329 Value += (*pointer)[J]-> 329 Value += (*pointer)[J]->
330 GetValue(LowEdgeEnergy,isOutRange) ; 330 GetValue(LowEdgeEnergy,isOutRange) ;
331 } 331 }
332 332
333 aVector->PutValue(i,Value) ; 333 aVector->PutValue(i,Value) ;
334 } 334 }
335 335
336 theDEDXTable->insert(aVector) ; 336 theDEDXTable->insert(aVector) ;
337 } 337 }
338 338
339 // reset counter to zero ................ 339 // reset counter to zero ..................
340 if( Charge >0.) 340 if( Charge >0.)
341 CounterOfpProcess=0 ; 341 CounterOfpProcess=0 ;
342 else 342 else
343 CounterOfpbarProcess=0 ; 343 CounterOfpbarProcess=0 ;
344 344
345 // Build range table 345 // Build range table
346 BuildRangeTable( aParticleType); 346 BuildRangeTable( aParticleType);
347 347
348 // Build lab/proper time tables 348 // Build lab/proper time tables
349 BuildTimeTables( aParticleType) ; 349 BuildTimeTables( aParticleType) ;
350 350
351 // Build coeff tables for the energy loss 351 // Build coeff tables for the energy loss calculation
352 BuildRangeCoeffATable( aParticleType); 352 BuildRangeCoeffATable( aParticleType);
353 BuildRangeCoeffBTable( aParticleType); 353 BuildRangeCoeffBTable( aParticleType);
354 BuildRangeCoeffCTable( aParticleType); 354 BuildRangeCoeffCTable( aParticleType);
355 355
356 // invert the range table 356 // invert the range table
357 357
358 BuildInverseRangeTable(aParticleType); 358 BuildInverseRangeTable(aParticleType);
359 } 359 }
360 } 360 }
361 // make the energy loss and the range table 361 // make the energy loss and the range table available
362 362
363 G4EnergyLossTables::Register(&aParticleType, 363 G4EnergyLossTables::Register(&aParticleType,
364 (Charge>0)? 364 (Charge>0)?
365 theDEDXpTable: theDEDXpbarTable, 365 theDEDXpTable: theDEDXpbarTable,
366 (Charge>0)? 366 (Charge>0)?
367 theRangepTable: theRangepbarTable 367 theRangepTable: theRangepbarTable,
368 (Charge>0)? 368 (Charge>0)?
369 theInverseRangepTable: theInverse 369 theInverseRangepTable: theInverseRangepbarTable,
370 (Charge>0)? 370 (Charge>0)?
371 theLabTimepTable: theLabTimepbarT 371 theLabTimepTable: theLabTimepbarTable,
372 (Charge>0)? 372 (Charge>0)?
373 theProperTimepTable: theProperTim 373 theProperTimepTable: theProperTimepbarTable,
374 LowestKineticEnergy, HighestKinet 374 LowestKineticEnergy, HighestKineticEnergy,
375 proton_mass_c2/aParticleType.GetP 375 proton_mass_c2/aParticleType.GetPDGMass(),
376 TotBin); 376 TotBin);
377 377
378 } 378 }
379 379
380 //-------------------------------------------- 380 //--------------------------------------------------------------------------------
381 381
382 void G4hRDEnergyLoss::BuildRangeTable(const G4 382 void G4hRDEnergyLoss::BuildRangeTable(const G4ParticleDefinition& aParticleType)
383 { 383 {
384 // Build range table from the energy loss ta 384 // Build range table from the energy loss table
385 385
386 Mass = aParticleType.GetPDGMass(); 386 Mass = aParticleType.GetPDGMass();
387 387
388 const G4ProductionCutsTable* theCoupleTable= 388 const G4ProductionCutsTable* theCoupleTable=
389 G4ProductionCutsTable::GetProductionCutsTa 389 G4ProductionCutsTable::GetProductionCutsTable();
390 G4int numOfCouples = (G4int)theCoupleTable-> 390 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
391 391
392 if( Charge >0.) 392 if( Charge >0.)
393 { 393 {
394 if(theRangepTable) 394 if(theRangepTable)
395 { theRangepTable->clearAndDestroy(); 395 { theRangepTable->clearAndDestroy();
396 delete theRangepTable; } 396 delete theRangepTable; }
397 theRangepTable = new G4PhysicsTable(numO 397 theRangepTable = new G4PhysicsTable(numOfCouples);
398 theRangeTable = theRangepTable ; 398 theRangeTable = theRangepTable ;
399 } 399 }
400 else 400 else
401 { 401 {
402 if(theRangepbarTable) 402 if(theRangepbarTable)
403 { theRangepbarTable->clearAndDestroy(); 403 { theRangepbarTable->clearAndDestroy();
404 delete theRangepbarTable; } 404 delete theRangepbarTable; }
405 theRangepbarTable = new G4PhysicsTable(n 405 theRangepbarTable = new G4PhysicsTable(numOfCouples);
406 theRangeTable = theRangepbarTable ; 406 theRangeTable = theRangepbarTable ;
407 } 407 }
408 408
409 // loop for materials 409 // loop for materials
410 410
411 for (G4int J=0; J<numOfCouples; ++J) 411 for (G4int J=0; J<numOfCouples; ++J)
412 { 412 {
413 G4PhysicsLogVector* aVector; 413 G4PhysicsLogVector* aVector;
414 aVector = new G4PhysicsLogVector(LowestK 414 aVector = new G4PhysicsLogVector(LowestKineticEnergy,
415 HighestKineticEnergy,TotBin); 415 HighestKineticEnergy,TotBin);
416 416
417 BuildRangeVector(J, aVector); 417 BuildRangeVector(J, aVector);
418 theRangeTable->insert(aVector); 418 theRangeTable->insert(aVector);
419 } 419 }
420 } 420 }
421 421
422 //-------------------------------------------- 422 //--------------------------------------------------------------------------------
423 423
424 void G4hRDEnergyLoss::BuildTimeTables(const G4 424 void G4hRDEnergyLoss::BuildTimeTables(const G4ParticleDefinition& aParticleType)
425 { 425 {
426 const G4ProductionCutsTable* theCoupleTable= 426 const G4ProductionCutsTable* theCoupleTable=
427 G4ProductionCutsTable::GetProductionCutsTa 427 G4ProductionCutsTable::GetProductionCutsTable();
428 G4int numOfCouples = (G4int)theCoupleTable-> 428 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
429 429
430 if(&aParticleType == G4Proton::Proton()) 430 if(&aParticleType == G4Proton::Proton())
431 { 431 {
432 if(theLabTimepTable) 432 if(theLabTimepTable)
433 { theLabTimepTable->clearAndDestroy(); 433 { theLabTimepTable->clearAndDestroy();
434 delete theLabTimepTable; } 434 delete theLabTimepTable; }
435 theLabTimepTable = new G4PhysicsTable(nu 435 theLabTimepTable = new G4PhysicsTable(numOfCouples);
436 theLabTimeTable = theLabTimepTable ; 436 theLabTimeTable = theLabTimepTable ;
437 437
438 if(theProperTimepTable) 438 if(theProperTimepTable)
439 { theProperTimepTable->clearAndDestroy(); 439 { theProperTimepTable->clearAndDestroy();
440 delete theProperTimepTable; } 440 delete theProperTimepTable; }
441 theProperTimepTable = new G4PhysicsTable 441 theProperTimepTable = new G4PhysicsTable(numOfCouples);
442 theProperTimeTable = theProperTimepTable 442 theProperTimeTable = theProperTimepTable ;
443 } 443 }
444 444
445 if(&aParticleType == G4AntiProton::AntiProto 445 if(&aParticleType == G4AntiProton::AntiProton())
446 { 446 {
447 if(theLabTimepbarTable) 447 if(theLabTimepbarTable)
448 { theLabTimepbarTable->clearAndDestroy(); 448 { theLabTimepbarTable->clearAndDestroy();
449 delete theLabTimepbarTable; } 449 delete theLabTimepbarTable; }
450 theLabTimepbarTable = new G4PhysicsTable 450 theLabTimepbarTable = new G4PhysicsTable(numOfCouples);
451 theLabTimeTable = theLabTimepbarTable ; 451 theLabTimeTable = theLabTimepbarTable ;
452 452
453 if(theProperTimepbarTable) 453 if(theProperTimepbarTable)
454 { theProperTimepbarTable->clearAndDestroy(); 454 { theProperTimepbarTable->clearAndDestroy();
455 delete theProperTimepbarTable; } 455 delete theProperTimepbarTable; }
456 theProperTimepbarTable = new G4PhysicsTa 456 theProperTimepbarTable = new G4PhysicsTable(numOfCouples);
457 theProperTimeTable = theProperTimepbarTa 457 theProperTimeTable = theProperTimepbarTable ;
458 } 458 }
459 459
460 for (G4int J=0; J<numOfCouples; ++J) 460 for (G4int J=0; J<numOfCouples; ++J)
461 { 461 {
462 G4PhysicsLogVector* aVector; 462 G4PhysicsLogVector* aVector;
463 G4PhysicsLogVector* bVector; 463 G4PhysicsLogVector* bVector;
464 464
465 aVector = new G4PhysicsLogVector(LowestK 465 aVector = new G4PhysicsLogVector(LowestKineticEnergy,
466 HighestKineticEnergy,TotBin); 466 HighestKineticEnergy,TotBin);
467 467
468 BuildLabTimeVector(J, aVector); 468 BuildLabTimeVector(J, aVector);
469 theLabTimeTable->insert(aVector); 469 theLabTimeTable->insert(aVector);
470 470
471 bVector = new G4PhysicsLogVector(LowestK 471 bVector = new G4PhysicsLogVector(LowestKineticEnergy,
472 HighestKineticEnergy,TotBin); 472 HighestKineticEnergy,TotBin);
473 473
474 BuildProperTimeVector(J, bVector); 474 BuildProperTimeVector(J, bVector);
475 theProperTimeTable->insert(bVector); 475 theProperTimeTable->insert(bVector);
476 } 476 }
477 } 477 }
478 478
479 //-------------------------------------------- 479 //--------------------------------------------------------------------------------
480 480
481 void G4hRDEnergyLoss::BuildRangeVector(G4int m 481 void G4hRDEnergyLoss::BuildRangeVector(G4int materialIndex,
482 G4PhysicsLogVector* rangeVector 482 G4PhysicsLogVector* rangeVector)
483 { 483 {
484 // create range vector for a material 484 // create range vector for a material
485 485
486 G4bool isOut; 486 G4bool isOut;
487 G4PhysicsVector* physicsVector= (*theDEDXTab 487 G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex];
488 G4double energy1 = rangeVector->GetLowEdgeEn 488 G4double energy1 = rangeVector->GetLowEdgeEnergy(0);
489 G4double dedx = physicsVector->GetValue(e 489 G4double dedx = physicsVector->GetValue(energy1,isOut);
490 G4double range = 0.5*energy1/dedx; 490 G4double range = 0.5*energy1/dedx;
491 rangeVector->PutValue(0,range); 491 rangeVector->PutValue(0,range);
492 G4int n = 100; 492 G4int n = 100;
493 G4double del = 1.0/(G4double)n ; 493 G4double del = 1.0/(G4double)n ;
494 494
495 for (G4int j=1; j<TotBin; j++) { 495 for (G4int j=1; j<TotBin; j++) {
496 496
497 G4double energy2 = rangeVector->GetLowEdge 497 G4double energy2 = rangeVector->GetLowEdgeEnergy(j);
498 G4double de = (energy2 - energy1) * del ; 498 G4double de = (energy2 - energy1) * del ;
499 G4double dedx1 = dedx ; 499 G4double dedx1 = dedx ;
500 500
501 for (G4int i=1; i<n; i++) { 501 for (G4int i=1; i<n; i++) {
502 G4double energy = energy1 + i*de ; 502 G4double energy = energy1 + i*de ;
503 G4double dedx2 = physicsVector->GetValu 503 G4double dedx2 = physicsVector->GetValue(energy,isOut);
504 range += 0.5*de*(1.0/dedx1 + 1.0/dedx2) 504 range += 0.5*de*(1.0/dedx1 + 1.0/dedx2);
505 dedx1 = dedx2; 505 dedx1 = dedx2;
506 } 506 }
507 rangeVector->PutValue(j,range); 507 rangeVector->PutValue(j,range);
508 dedx = dedx1 ; 508 dedx = dedx1 ;
509 energy1 = energy2 ; 509 energy1 = energy2 ;
510 } 510 }
511 } 511 }
512 512
513 //-------------------------------------------- 513 //--------------------------------------------------------------------------------
514 514
515 void G4hRDEnergyLoss::BuildLabTimeVector(G4int 515 void G4hRDEnergyLoss::BuildLabTimeVector(G4int materialIndex,
516 G4PhysicsLogVector* timeVector) 516 G4PhysicsLogVector* timeVector)
517 { 517 {
518 // create lab time vector for a material 518 // create lab time vector for a material
519 519
520 G4int nbin=100; 520 G4int nbin=100;
521 G4bool isOut; 521 G4bool isOut;
522 G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-p 522 G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-parlowen ;
523 //G4double losslim,clim,taulim,timelim,ltaul 523 //G4double losslim,clim,taulim,timelim,ltaulim,ltaumax,
524 G4double losslim,clim,taulim,timelim, 524 G4double losslim,clim,taulim,timelim,
525 LowEdgeEnergy,tau,Value ; 525 LowEdgeEnergy,tau,Value ;
526 526
527 G4PhysicsVector* physicsVector= (*theDEDXTab 527 G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex];
528 528
529 // low energy part first... 529 // low energy part first...
530 losslim = physicsVector->GetValue(tlim,isOut 530 losslim = physicsVector->GetValue(tlim,isOut);
531 taulim=tlim/ParticleMass ; 531 taulim=tlim/ParticleMass ;
532 clim=std::sqrt(ParticleMass*tlim/2.)/(c_ligh 532 clim=std::sqrt(ParticleMass*tlim/2.)/(c_light*losslim*ppar) ;
533 //ltaulim = std::log(taulim); 533 //ltaulim = std::log(taulim);
534 //ltaumax = std::log(HighestKineticEnergy/Pa 534 //ltaumax = std::log(HighestKineticEnergy/ParticleMass) ;
535 535
536 G4int i=-1; 536 G4int i=-1;
537 G4double oldValue = 0. ; 537 G4double oldValue = 0. ;
538 G4double tauold ; 538 G4double tauold ;
539 do 539 do
540 { 540 {
541 i += 1 ; 541 i += 1 ;
542 LowEdgeEnergy = timeVector->GetLowEdgeEn 542 LowEdgeEnergy = timeVector->GetLowEdgeEnergy(i);
543 tau = LowEdgeEnergy/ParticleMass ; 543 tau = LowEdgeEnergy/ParticleMass ;
544 if ( tau <= taulim ) 544 if ( tau <= taulim )
545 { 545 {
546 Value = clim*std::exp(ppar*std::log(tau/ta 546 Value = clim*std::exp(ppar*std::log(tau/taulim)) ;
547 } 547 }
548 else 548 else
549 { 549 {
550 timelim=clim ; 550 timelim=clim ;
551 ltaulow = std::log(taulim); 551 ltaulow = std::log(taulim);
552 ltauhigh = std::log(tau); 552 ltauhigh = std::log(tau);
553 Value = timelim+LabTimeIntLog(physicsVecto 553 Value = timelim+LabTimeIntLog(physicsVector,nbin);
554 } 554 }
555 timeVector->PutValue(i,Value); 555 timeVector->PutValue(i,Value);
556 oldValue = Value ; 556 oldValue = Value ;
557 tauold = tau ; 557 tauold = tau ;
558 } while (tau<=taulim) ; 558 } while (tau<=taulim) ;
559 559
560 i += 1 ; 560 i += 1 ;
561 for (G4int j=i; j<TotBin; j++) 561 for (G4int j=i; j<TotBin; j++)
562 { 562 {
563 LowEdgeEnergy = timeVector->GetLowEdgeEn 563 LowEdgeEnergy = timeVector->GetLowEdgeEnergy(j);
564 tau = LowEdgeEnergy/ParticleMass ; 564 tau = LowEdgeEnergy/ParticleMass ;
565 ltaulow = std::log(tauold); 565 ltaulow = std::log(tauold);
566 ltauhigh = std::log(tau); 566 ltauhigh = std::log(tau);
567 Value = oldValue+LabTimeIntLog(physicsVe 567 Value = oldValue+LabTimeIntLog(physicsVector,nbin);
568 timeVector->PutValue(j,Value); 568 timeVector->PutValue(j,Value);
569 oldValue = Value ; 569 oldValue = Value ;
570 tauold = tau ; 570 tauold = tau ;
571 } 571 }
572 } 572 }
573 573
574 //-------------------------------------------- 574 //--------------------------------------------------------------------------------
575 575
576 void G4hRDEnergyLoss::BuildProperTimeVector(G4 576 void G4hRDEnergyLoss::BuildProperTimeVector(G4int materialIndex,
577 G4PhysicsLogVector* timeVector) 577 G4PhysicsLogVector* timeVector)
578 { 578 {
579 // create proper time vector for a material 579 // create proper time vector for a material
580 580
581 G4int nbin=100; 581 G4int nbin=100;
582 G4bool isOut; 582 G4bool isOut;
583 G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-p 583 G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-parlowen ;
584 //G4double losslim,clim,taulim,timelim,ltaul 584 //G4double losslim,clim,taulim,timelim,ltaulim,ltaumax,
585 G4double losslim,clim,taulim,timelim, 585 G4double losslim,clim,taulim,timelim,
586 LowEdgeEnergy,tau,Value ; 586 LowEdgeEnergy,tau,Value ;
587 587
588 G4PhysicsVector* physicsVector= (*theDEDXTab 588 G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex];
589 589
590 // low energy part first... 590 // low energy part first...
591 losslim = physicsVector->GetValue(tlim,isOut 591 losslim = physicsVector->GetValue(tlim,isOut);
592 taulim=tlim/ParticleMass ; 592 taulim=tlim/ParticleMass ;
593 clim=std::sqrt(ParticleMass*tlim/2.)/(c_ligh 593 clim=std::sqrt(ParticleMass*tlim/2.)/(c_light*losslim*ppar) ;
594 //ltaulim = std::log(taulim); 594 //ltaulim = std::log(taulim);
595 //ltaumax = std::log(HighestKineticEnergy/Pa 595 //ltaumax = std::log(HighestKineticEnergy/ParticleMass) ;
596 596
597 G4int i=-1; 597 G4int i=-1;
598 G4double oldValue = 0. ; 598 G4double oldValue = 0. ;
599 G4double tauold ; 599 G4double tauold ;
600 do 600 do
601 { 601 {
602 i += 1 ; 602 i += 1 ;
603 LowEdgeEnergy = timeVector->GetLowEdgeEn 603 LowEdgeEnergy = timeVector->GetLowEdgeEnergy(i);
604 tau = LowEdgeEnergy/ParticleMass ; 604 tau = LowEdgeEnergy/ParticleMass ;
605 if ( tau <= taulim ) 605 if ( tau <= taulim )
606 { 606 {
607 Value = clim*std::exp(ppar*std::log(tau/ta 607 Value = clim*std::exp(ppar*std::log(tau/taulim)) ;
608 } 608 }
609 else 609 else
610 { 610 {
611 timelim=clim ; 611 timelim=clim ;
612 ltaulow = std::log(taulim); 612 ltaulow = std::log(taulim);
613 ltauhigh = std::log(tau); 613 ltauhigh = std::log(tau);
614 Value = timelim+ProperTimeIntLog(physicsVe 614 Value = timelim+ProperTimeIntLog(physicsVector,nbin);
615 } 615 }
616 timeVector->PutValue(i,Value); 616 timeVector->PutValue(i,Value);
617 oldValue = Value ; 617 oldValue = Value ;
618 tauold = tau ; 618 tauold = tau ;
619 } while (tau<=taulim) ; 619 } while (tau<=taulim) ;
620 620
621 i += 1 ; 621 i += 1 ;
622 for (G4int j=i; j<TotBin; j++) 622 for (G4int j=i; j<TotBin; j++)
623 { 623 {
624 LowEdgeEnergy = timeVector->GetLowEdgeEn 624 LowEdgeEnergy = timeVector->GetLowEdgeEnergy(j);
625 tau = LowEdgeEnergy/ParticleMass ; 625 tau = LowEdgeEnergy/ParticleMass ;
626 ltaulow = std::log(tauold); 626 ltaulow = std::log(tauold);
627 ltauhigh = std::log(tau); 627 ltauhigh = std::log(tau);
628 Value = oldValue+ProperTimeIntLog(physic 628 Value = oldValue+ProperTimeIntLog(physicsVector,nbin);
629 timeVector->PutValue(j,Value); 629 timeVector->PutValue(j,Value);
630 oldValue = Value ; 630 oldValue = Value ;
631 tauold = tau ; 631 tauold = tau ;
632 } 632 }
633 } 633 }
634 634
635 //-------------------------------------------- 635 //--------------------------------------------------------------------------------
636 636
637 G4double G4hRDEnergyLoss::RangeIntLin(G4Physic 637 G4double G4hRDEnergyLoss::RangeIntLin(G4PhysicsVector* physicsVector,
638 G4int nbin) 638 G4int nbin)
639 { 639 {
640 // num. integration, linear binning 640 // num. integration, linear binning
641 641
642 G4double dtau,Value,taui,ti,lossi,ci; 642 G4double dtau,Value,taui,ti,lossi,ci;
643 G4bool isOut; 643 G4bool isOut;
644 dtau = (tauhigh-taulow)/nbin; 644 dtau = (tauhigh-taulow)/nbin;
645 Value = 0.; 645 Value = 0.;
646 646
647 for (G4int i=0; i<=nbin; i++) 647 for (G4int i=0; i<=nbin; i++)
648 { 648 {
649 taui = taulow + dtau*i ; 649 taui = taulow + dtau*i ;
650 ti = Mass*taui; 650 ti = Mass*taui;
651 lossi = physicsVector->GetValue(ti,isOut 651 lossi = physicsVector->GetValue(ti,isOut);
652 if(i==0) 652 if(i==0)
653 ci=0.5; 653 ci=0.5;
654 else 654 else
655 { 655 {
656 if(i<nbin) 656 if(i<nbin)
657 ci=1.; 657 ci=1.;
658 else 658 else
659 ci=0.5; 659 ci=0.5;
660 } 660 }
661 Value += ci/lossi; 661 Value += ci/lossi;
662 } 662 }
663 Value *= Mass*dtau; 663 Value *= Mass*dtau;
664 return Value; 664 return Value;
665 } 665 }
666 666
667 //-------------------------------------------- 667 //--------------------------------------------------------------------------------
668 668
669 G4double G4hRDEnergyLoss::RangeIntLog(G4Physic 669 G4double G4hRDEnergyLoss::RangeIntLog(G4PhysicsVector* physicsVector,
670 G4int nbin) 670 G4int nbin)
671 { 671 {
672 // num. integration, logarithmic binning 672 // num. integration, logarithmic binning
673 673
674 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci 674 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
675 G4bool isOut; 675 G4bool isOut;
676 ltt = ltauhigh-ltaulow; 676 ltt = ltauhigh-ltaulow;
677 dltau = ltt/nbin; 677 dltau = ltt/nbin;
678 Value = 0.; 678 Value = 0.;
679 679
680 for (G4int i=0; i<=nbin; i++) 680 for (G4int i=0; i<=nbin; i++)
681 { 681 {
682 ui = ltaulow+dltau*i; 682 ui = ltaulow+dltau*i;
683 taui = std::exp(ui); 683 taui = std::exp(ui);
684 ti = Mass*taui; 684 ti = Mass*taui;
685 lossi = physicsVector->GetValue(ti,isOut 685 lossi = physicsVector->GetValue(ti,isOut);
686 if(i==0) 686 if(i==0)
687 ci=0.5; 687 ci=0.5;
688 else 688 else
689 { 689 {
690 if(i<nbin) 690 if(i<nbin)
691 ci=1.; 691 ci=1.;
692 else 692 else
693 ci=0.5; 693 ci=0.5;
694 } 694 }
695 Value += ci*taui/lossi; 695 Value += ci*taui/lossi;
696 } 696 }
697 Value *= Mass*dltau; 697 Value *= Mass*dltau;
698 return Value; 698 return Value;
699 } 699 }
700 700
701 //-------------------------------------------- 701 //--------------------------------------------------------------------------------
702 702
703 G4double G4hRDEnergyLoss::LabTimeIntLog(G4Phys 703 G4double G4hRDEnergyLoss::LabTimeIntLog(G4PhysicsVector* physicsVector,
704 G4int nbin) 704 G4int nbin)
705 { 705 {
706 // num. integration, logarithmic binning 706 // num. integration, logarithmic binning
707 707
708 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci 708 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
709 G4bool isOut; 709 G4bool isOut;
710 ltt = ltauhigh-ltaulow; 710 ltt = ltauhigh-ltaulow;
711 dltau = ltt/nbin; 711 dltau = ltt/nbin;
712 Value = 0.; 712 Value = 0.;
713 713
714 for (G4int i=0; i<=nbin; i++) 714 for (G4int i=0; i<=nbin; i++)
715 { 715 {
716 ui = ltaulow+dltau*i; 716 ui = ltaulow+dltau*i;
717 taui = std::exp(ui); 717 taui = std::exp(ui);
718 ti = ParticleMass*taui; 718 ti = ParticleMass*taui;
719 lossi = physicsVector->GetValue(ti,isOut 719 lossi = physicsVector->GetValue(ti,isOut);
720 if(i==0) 720 if(i==0)
721 ci=0.5; 721 ci=0.5;
722 else 722 else
723 { 723 {
724 if(i<nbin) 724 if(i<nbin)
725 ci=1.; 725 ci=1.;
726 else 726 else
727 ci=0.5; 727 ci=0.5;
728 } 728 }
729 Value += ci*taui*(ti+ParticleMass)/(std: 729 Value += ci*taui*(ti+ParticleMass)/(std::sqrt(ti*(ti+2.*ParticleMass))*lossi);
730 } 730 }
731 Value *= ParticleMass*dltau/c_light; 731 Value *= ParticleMass*dltau/c_light;
732 return Value; 732 return Value;
733 } 733 }
734 734
735 //-------------------------------------------- 735 //--------------------------------------------------------------------------------
736 736
737 G4double G4hRDEnergyLoss::ProperTimeIntLog(G4P 737 G4double G4hRDEnergyLoss::ProperTimeIntLog(G4PhysicsVector* physicsVector,
738 G4int nbin) 738 G4int nbin)
739 { 739 {
740 // num. integration, logarithmic binning 740 // num. integration, logarithmic binning
741 741
742 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci 742 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
743 G4bool isOut; 743 G4bool isOut;
744 ltt = ltauhigh-ltaulow; 744 ltt = ltauhigh-ltaulow;
745 dltau = ltt/nbin; 745 dltau = ltt/nbin;
746 Value = 0.; 746 Value = 0.;
747 747
748 for (G4int i=0; i<=nbin; i++) 748 for (G4int i=0; i<=nbin; i++)
749 { 749 {
750 ui = ltaulow+dltau*i; 750 ui = ltaulow+dltau*i;
751 taui = std::exp(ui); 751 taui = std::exp(ui);
752 ti = ParticleMass*taui; 752 ti = ParticleMass*taui;
753 lossi = physicsVector->GetValue(ti,isOut 753 lossi = physicsVector->GetValue(ti,isOut);
754 if(i==0) 754 if(i==0)
755 ci=0.5; 755 ci=0.5;
756 else 756 else
757 { 757 {
758 if(i<nbin) 758 if(i<nbin)
759 ci=1.; 759 ci=1.;
760 else 760 else
761 ci=0.5; 761 ci=0.5;
762 } 762 }
763 Value += ci*taui*ParticleMass/(std::sqrt 763 Value += ci*taui*ParticleMass/(std::sqrt(ti*(ti+2.*ParticleMass))*lossi);
764 } 764 }
765 Value *= ParticleMass*dltau/c_light; 765 Value *= ParticleMass*dltau/c_light;
766 return Value; 766 return Value;
767 } 767 }
768 768
769 //-------------------------------------------- 769 //--------------------------------------------------------------------------------
770 770
771 void G4hRDEnergyLoss::BuildRangeCoeffATable( c 771 void G4hRDEnergyLoss::BuildRangeCoeffATable( const G4ParticleDefinition& )
772 { 772 {
773 // Build tables of coefficients for the ener 773 // Build tables of coefficients for the energy loss calculation
774 // create table for coefficients "A" 774 // create table for coefficients "A"
775 775
776 G4int numOfCouples = (G4int)G4ProductionCuts 776 G4int numOfCouples = (G4int)G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
777 777
778 if(Charge>0.) 778 if(Charge>0.)
779 { 779 {
780 if(thepRangeCoeffATable) 780 if(thepRangeCoeffATable)
781 { thepRangeCoeffATable->clearAndDestroy(); 781 { thepRangeCoeffATable->clearAndDestroy();
782 delete thepRangeCoeffATable; } 782 delete thepRangeCoeffATable; }
783 thepRangeCoeffATable = new G4PhysicsTabl 783 thepRangeCoeffATable = new G4PhysicsTable(numOfCouples);
784 theRangeCoeffATable = thepRangeCoeffATab 784 theRangeCoeffATable = thepRangeCoeffATable ;
785 theRangeTable = theRangepTable ; 785 theRangeTable = theRangepTable ;
786 } 786 }
787 else 787 else
788 { 788 {
789 if(thepbarRangeCoeffATable) 789 if(thepbarRangeCoeffATable)
790 { thepbarRangeCoeffATable->clearAndDestroy() 790 { thepbarRangeCoeffATable->clearAndDestroy();
791 delete thepbarRangeCoeffATable; } 791 delete thepbarRangeCoeffATable; }
792 thepbarRangeCoeffATable = new G4PhysicsT 792 thepbarRangeCoeffATable = new G4PhysicsTable(numOfCouples);
793 theRangeCoeffATable = thepbarRangeCoeffA 793 theRangeCoeffATable = thepbarRangeCoeffATable ;
794 theRangeTable = theRangepbarTable ; 794 theRangeTable = theRangepbarTable ;
795 } 795 }
796 796
797 G4double R2 = RTable*RTable ; 797 G4double R2 = RTable*RTable ;
798 G4double R1 = RTable+1.; 798 G4double R1 = RTable+1.;
799 G4double w = R1*(RTable-1.)*(RTable-1.); 799 G4double w = R1*(RTable-1.)*(RTable-1.);
800 G4double w1 = RTable/w , w2 = -RTable*R1/w , 800 G4double w1 = RTable/w , w2 = -RTable*R1/w , w3 = R2/w ;
801 G4double Ti , Tim , Tip , Ri , Rim , Rip , V 801 G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
802 G4bool isOut; 802 G4bool isOut;
803 803
804 // loop for materials 804 // loop for materials
805 for (G4int J=0; J<numOfCouples; J++) 805 for (G4int J=0; J<numOfCouples; J++)
806 { 806 {
807 G4int binmax=TotBin ; 807 G4int binmax=TotBin ;
808 G4PhysicsLinearVector* aVector = 808 G4PhysicsLinearVector* aVector =
809 new G4PhysicsLinearVector(0.,binmax, TotBin) 809 new G4PhysicsLinearVector(0.,binmax, TotBin);
810 Ti = LowestKineticEnergy ; 810 Ti = LowestKineticEnergy ;
811 if (Ti < DBL_MIN) Ti = 1.e-8; 811 if (Ti < DBL_MIN) Ti = 1.e-8;
812 G4PhysicsVector* rangeVector= (*theRange 812 G4PhysicsVector* rangeVector= (*theRangeTable)[J];
813 813
814 for ( G4int i=0; i<TotBin; i++) 814 for ( G4int i=0; i<TotBin; i++)
815 { 815 {
816 Ri = rangeVector->GetValue(Ti,isOut) ; 816 Ri = rangeVector->GetValue(Ti,isOut) ;
817 if (Ti < DBL_MIN) Ti = 1.e-8; 817 if (Ti < DBL_MIN) Ti = 1.e-8;
818 if ( i==0 ) 818 if ( i==0 )
819 Rim = 0. ; 819 Rim = 0. ;
820 else 820 else
821 { 821 {
822 // ---- MGP ---- Modified to avoid a f 822 // ---- MGP ---- Modified to avoid a floating point exception
823 // The correction is just a temporary 823 // The correction is just a temporary patch, the whole class should be redesigned
824 // Original: Tim = Ti/RTable results 824 // Original: Tim = Ti/RTable results in 0./0.
825 if (RTable != 0.) 825 if (RTable != 0.)
826 { 826 {
827 Tim = Ti/RTable ; 827 Tim = Ti/RTable ;
828 } 828 }
829 else 829 else
830 { 830 {
831 Tim = 0.; 831 Tim = 0.;
832 } 832 }
833 Rim = rangeVector->GetValue(Tim,isOut) 833 Rim = rangeVector->GetValue(Tim,isOut);
834 } 834 }
835 if ( i==(TotBin-1)) 835 if ( i==(TotBin-1))
836 Rip = Ri ; 836 Rip = Ri ;
837 else 837 else
838 { 838 {
839 Tip = Ti*RTable ; 839 Tip = Ti*RTable ;
840 Rip = rangeVector->GetValue(Tip,isOut) 840 Rip = rangeVector->GetValue(Tip,isOut);
841 } 841 }
842 Value = (w1*Rip + w2*Ri + w3*Rim)/(Ti*Ti) 842 Value = (w1*Rip + w2*Ri + w3*Rim)/(Ti*Ti) ;
843 843
844 aVector->PutValue(i,Value); 844 aVector->PutValue(i,Value);
845 Ti = RTable*Ti ; 845 Ti = RTable*Ti ;
846 } 846 }
847 847
848 theRangeCoeffATable->insert(aVector); 848 theRangeCoeffATable->insert(aVector);
849 } 849 }
850 } 850 }
851 851
852 //-------------------------------------------- 852 //--------------------------------------------------------------------------------
853 853
854 void G4hRDEnergyLoss::BuildRangeCoeffBTable( c 854 void G4hRDEnergyLoss::BuildRangeCoeffBTable( const G4ParticleDefinition& )
855 { 855 {
856 // Build tables of coefficients for the ener 856 // Build tables of coefficients for the energy loss calculation
857 // create table for coefficients "B" 857 // create table for coefficients "B"
858 858
859 G4int numOfCouples = 859 G4int numOfCouples =
860 (G4int)G4ProductionCutsTable::GetProdu 860 (G4int)G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
861 861
862 if(Charge>0.) 862 if(Charge>0.)
863 { 863 {
864 if(thepRangeCoeffBTable) 864 if(thepRangeCoeffBTable)
865 { thepRangeCoeffBTable->clearAndDestroy(); 865 { thepRangeCoeffBTable->clearAndDestroy();
866 delete thepRangeCoeffBTable; } 866 delete thepRangeCoeffBTable; }
867 thepRangeCoeffBTable = new G4PhysicsTabl 867 thepRangeCoeffBTable = new G4PhysicsTable(numOfCouples);
868 theRangeCoeffBTable = thepRangeCoeffBTab 868 theRangeCoeffBTable = thepRangeCoeffBTable ;
869 theRangeTable = theRangepTable ; 869 theRangeTable = theRangepTable ;
870 } 870 }
871 else 871 else
872 { 872 {
873 if(thepbarRangeCoeffBTable) 873 if(thepbarRangeCoeffBTable)
874 { thepbarRangeCoeffBTable->clearAndDestroy() 874 { thepbarRangeCoeffBTable->clearAndDestroy();
875 delete thepbarRangeCoeffBTable; } 875 delete thepbarRangeCoeffBTable; }
876 thepbarRangeCoeffBTable = new G4PhysicsT 876 thepbarRangeCoeffBTable = new G4PhysicsTable(numOfCouples);
877 theRangeCoeffBTable = thepbarRangeCoeffB 877 theRangeCoeffBTable = thepbarRangeCoeffBTable ;
878 theRangeTable = theRangepbarTable ; 878 theRangeTable = theRangepbarTable ;
879 } 879 }
880 880
881 G4double R2 = RTable*RTable ; 881 G4double R2 = RTable*RTable ;
882 G4double R1 = RTable+1.; 882 G4double R1 = RTable+1.;
883 G4double w = R1*(RTable-1.)*(RTable-1.); 883 G4double w = R1*(RTable-1.)*(RTable-1.);
884 if (w < DBL_MIN) w = DBL_MIN; 884 if (w < DBL_MIN) w = DBL_MIN;
885 G4double w1 = -R1/w , w2 = R1*(R2+1.)/w , w3 885 G4double w1 = -R1/w , w2 = R1*(R2+1.)/w , w3 = -R2*R1/w ;
886 G4double Ti , Tim , Tip , Ri , Rim , Rip , V 886 G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
887 G4bool isOut; 887 G4bool isOut;
888 888
889 // loop for materials 889 // loop for materials
890 for (G4int J=0; J<numOfCouples; J++) 890 for (G4int J=0; J<numOfCouples; J++)
891 { 891 {
892 G4int binmax=TotBin ; 892 G4int binmax=TotBin ;
893 G4PhysicsLinearVector* aVector = 893 G4PhysicsLinearVector* aVector =
894 new G4PhysicsLinearVector(0.,binmax, TotBin) 894 new G4PhysicsLinearVector(0.,binmax, TotBin);
895 Ti = LowestKineticEnergy ; 895 Ti = LowestKineticEnergy ;
896 if (Ti < DBL_MIN) Ti = 1.e-8; 896 if (Ti < DBL_MIN) Ti = 1.e-8;
897 G4PhysicsVector* rangeVector= (*theRange 897 G4PhysicsVector* rangeVector= (*theRangeTable)[J];
898 898
899 for ( G4int i=0; i<TotBin; i++) 899 for ( G4int i=0; i<TotBin; i++)
900 { 900 {
901 Ri = rangeVector->GetValue(Ti,isOut) ; 901 Ri = rangeVector->GetValue(Ti,isOut) ;
902 if (Ti < DBL_MIN) Ti = 1.e-8; 902 if (Ti < DBL_MIN) Ti = 1.e-8;
903 if ( i==0 ) 903 if ( i==0 )
904 Rim = 0. ; 904 Rim = 0. ;
905 else 905 else
906 { 906 {
907 if (RTable < DBL_MIN) RTable = DBL_MIN 907 if (RTable < DBL_MIN) RTable = DBL_MIN;
908 Tim = Ti/RTable ; 908 Tim = Ti/RTable ;
909 Rim = rangeVector->GetValue(Tim,isOut) 909 Rim = rangeVector->GetValue(Tim,isOut);
910 } 910 }
911 if ( i==(TotBin-1)) 911 if ( i==(TotBin-1))
912 Rip = Ri ; 912 Rip = Ri ;
913 else 913 else
914 { 914 {
915 Tip = Ti*RTable ; 915 Tip = Ti*RTable ;
916 Rip = rangeVector->GetValue(Tip,isOut) 916 Rip = rangeVector->GetValue(Tip,isOut);
917 } 917 }
918 if (Ti < DBL_MIN) Ti = DBL_MIN; 918 if (Ti < DBL_MIN) Ti = DBL_MIN;
919 Value = (w1*Rip + w2*Ri + w3*Rim)/Ti; 919 Value = (w1*Rip + w2*Ri + w3*Rim)/Ti;
920 920
921 aVector->PutValue(i,Value); 921 aVector->PutValue(i,Value);
922 Ti = RTable*Ti ; 922 Ti = RTable*Ti ;
923 } 923 }
924 theRangeCoeffBTable->insert(aVector); 924 theRangeCoeffBTable->insert(aVector);
925 } 925 }
926 } 926 }
927 927
928 //-------------------------------------------- 928 //--------------------------------------------------------------------------------
929 929
930 void G4hRDEnergyLoss::BuildRangeCoeffCTable( c 930 void G4hRDEnergyLoss::BuildRangeCoeffCTable( const G4ParticleDefinition& )
931 { 931 {
932 // Build tables of coefficients for the ener 932 // Build tables of coefficients for the energy loss calculation
933 // create table for coefficients "C" 933 // create table for coefficients "C"
934 934
935 G4int numOfCouples = 935 G4int numOfCouples =
936 (G4int)G4ProductionCutsTable::GetProdu 936 (G4int)G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
937 937
938 if(Charge>0.) 938 if(Charge>0.)
939 { 939 {
940 if(thepRangeCoeffCTable) 940 if(thepRangeCoeffCTable)
941 { thepRangeCoeffCTable->clearAndDestroy(); 941 { thepRangeCoeffCTable->clearAndDestroy();
942 delete thepRangeCoeffCTable; } 942 delete thepRangeCoeffCTable; }
943 thepRangeCoeffCTable = new G4PhysicsTabl 943 thepRangeCoeffCTable = new G4PhysicsTable(numOfCouples);
944 theRangeCoeffCTable = thepRangeCoeffCTab 944 theRangeCoeffCTable = thepRangeCoeffCTable ;
945 theRangeTable = theRangepTable ; 945 theRangeTable = theRangepTable ;
946 } 946 }
947 else 947 else
948 { 948 {
949 if(thepbarRangeCoeffCTable) 949 if(thepbarRangeCoeffCTable)
950 { thepbarRangeCoeffCTable->clearAndDestroy() 950 { thepbarRangeCoeffCTable->clearAndDestroy();
951 delete thepbarRangeCoeffCTable; } 951 delete thepbarRangeCoeffCTable; }
952 thepbarRangeCoeffCTable = new G4PhysicsT 952 thepbarRangeCoeffCTable = new G4PhysicsTable(numOfCouples);
953 theRangeCoeffCTable = thepbarRangeCoeffC 953 theRangeCoeffCTable = thepbarRangeCoeffCTable ;
954 theRangeTable = theRangepbarTable ; 954 theRangeTable = theRangepbarTable ;
955 } 955 }
956 956
957 G4double R2 = RTable*RTable ; 957 G4double R2 = RTable*RTable ;
958 G4double R1 = RTable+1.; 958 G4double R1 = RTable+1.;
959 G4double w = R1*(RTable-1.)*(RTable-1.); 959 G4double w = R1*(RTable-1.)*(RTable-1.);
960 G4double w1 = 1./w , w2 = -RTable*R1/w , w3 960 G4double w1 = 1./w , w2 = -RTable*R1/w , w3 = RTable*R2/w ;
961 G4double Ti , Tim , Tip , Ri , Rim , Rip , V 961 G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
962 G4bool isOut; 962 G4bool isOut;
963 963
964 // loop for materials 964 // loop for materials
965 for (G4int J=0; J<numOfCouples; J++) 965 for (G4int J=0; J<numOfCouples; J++)
966 { 966 {
967 G4int binmax=TotBin ; 967 G4int binmax=TotBin ;
968 G4PhysicsLinearVector* aVector = 968 G4PhysicsLinearVector* aVector =
969 new G4PhysicsLinearVector(0.,binmax, TotBin) 969 new G4PhysicsLinearVector(0.,binmax, TotBin);
970 Ti = LowestKineticEnergy ; 970 Ti = LowestKineticEnergy ;
971 G4PhysicsVector* rangeVector= (*theRange 971 G4PhysicsVector* rangeVector= (*theRangeTable)[J];
972 972
973 for ( G4int i=0; i<TotBin; i++) 973 for ( G4int i=0; i<TotBin; i++)
974 { 974 {
975 Ri = rangeVector->GetValue(Ti,isOut) ; 975 Ri = rangeVector->GetValue(Ti,isOut) ;
976 if ( i==0 ) 976 if ( i==0 )
977 Rim = 0. ; 977 Rim = 0. ;
978 else 978 else
979 { 979 {
980 Tim = Ti/RTable ; 980 Tim = Ti/RTable ;
981 Rim = rangeVector->GetValue(Tim,isOut) 981 Rim = rangeVector->GetValue(Tim,isOut);
982 } 982 }
983 if ( i==(TotBin-1)) 983 if ( i==(TotBin-1))
984 Rip = Ri ; 984 Rip = Ri ;
985 else 985 else
986 { 986 {
987 Tip = Ti*RTable ; 987 Tip = Ti*RTable ;
988 Rip = rangeVector->GetValue(Tip,isOut) 988 Rip = rangeVector->GetValue(Tip,isOut);
989 } 989 }
990 Value = w1*Rip + w2*Ri + w3*Rim ; 990 Value = w1*Rip + w2*Ri + w3*Rim ;
991 991
992 aVector->PutValue(i,Value); 992 aVector->PutValue(i,Value);
993 Ti = RTable*Ti ; 993 Ti = RTable*Ti ;
994 } 994 }
995 theRangeCoeffCTable->insert(aVector); 995 theRangeCoeffCTable->insert(aVector);
996 } 996 }
997 } 997 }
998 998
999 //-------------------------------------------- 999 //--------------------------------------------------------------------------------
1000 1000
1001 void G4hRDEnergyLoss:: 1001 void G4hRDEnergyLoss::
1002 BuildInverseRangeTable(const G4ParticleDefini 1002 BuildInverseRangeTable(const G4ParticleDefinition& aParticleType)
1003 { 1003 {
1004 // Build inverse table of the range table 1004 // Build inverse table of the range table
1005 1005
1006 G4bool b; 1006 G4bool b;
1007 1007
1008 const G4ProductionCutsTable* theCoupleTable 1008 const G4ProductionCutsTable* theCoupleTable=
1009 G4ProductionCutsTable::GetProductionCutsT 1009 G4ProductionCutsTable::GetProductionCutsTable();
1010 std::size_t numOfCouples = theCoupleTable-> 1010 std::size_t numOfCouples = theCoupleTable->GetTableSize();
1011 1011
1012 if(&aParticleType == G4Proton::Proton()) 1012 if(&aParticleType == G4Proton::Proton())
1013 { 1013 {
1014 if(theInverseRangepTable) 1014 if(theInverseRangepTable)
1015 { theInverseRangepTable->clearAndDestroy(); 1015 { theInverseRangepTable->clearAndDestroy();
1016 delete theInverseRangepTable; } 1016 delete theInverseRangepTable; }
1017 theInverseRangepTable = new G4PhysicsTa 1017 theInverseRangepTable = new G4PhysicsTable(numOfCouples);
1018 theInverseRangeTable = theInverseRangep 1018 theInverseRangeTable = theInverseRangepTable ;
1019 theRangeTable = theRangepTable ; 1019 theRangeTable = theRangepTable ;
1020 theDEDXTable = theDEDXpTable ; 1020 theDEDXTable = theDEDXpTable ;
1021 theRangeCoeffATable = thepRangeCoeffATa 1021 theRangeCoeffATable = thepRangeCoeffATable ;
1022 theRangeCoeffBTable = thepRangeCoeffBTa 1022 theRangeCoeffBTable = thepRangeCoeffBTable ;
1023 theRangeCoeffCTable = thepRangeCoeffCTa 1023 theRangeCoeffCTable = thepRangeCoeffCTable ;
1024 } 1024 }
1025 1025
1026 if(&aParticleType == G4AntiProton::AntiProt 1026 if(&aParticleType == G4AntiProton::AntiProton())
1027 { 1027 {
1028 if(theInverseRangepbarTable) 1028 if(theInverseRangepbarTable)
1029 { theInverseRangepbarTable->clearAndDestroy 1029 { theInverseRangepbarTable->clearAndDestroy();
1030 delete theInverseRangepbarTable; } 1030 delete theInverseRangepbarTable; }
1031 theInverseRangepbarTable = new G4Physic 1031 theInverseRangepbarTable = new G4PhysicsTable(numOfCouples);
1032 theInverseRangeTable = theInverseRangep 1032 theInverseRangeTable = theInverseRangepbarTable ;
1033 theRangeTable = theRangepbarTable ; 1033 theRangeTable = theRangepbarTable ;
1034 theDEDXTable = theDEDXpbarTable ; 1034 theDEDXTable = theDEDXpbarTable ;
1035 theRangeCoeffATable = thepbarRangeCoeff 1035 theRangeCoeffATable = thepbarRangeCoeffATable ;
1036 theRangeCoeffBTable = thepbarRangeCoeff 1036 theRangeCoeffBTable = thepbarRangeCoeffBTable ;
1037 theRangeCoeffCTable = thepbarRangeCoeff 1037 theRangeCoeffCTable = thepbarRangeCoeffCTable ;
1038 } 1038 }
1039 1039
1040 // loop for materials 1040 // loop for materials
1041 for (std::size_t i=0; i<numOfCouples; ++i) 1041 for (std::size_t i=0; i<numOfCouples; ++i)
1042 { 1042 {
1043 1043
1044 G4PhysicsVector* pv = (*theRangeTable)[ 1044 G4PhysicsVector* pv = (*theRangeTable)[i];
1045 std::size_t nbins = pv->GetVectorLeng 1045 std::size_t nbins = pv->GetVectorLength();
1046 G4double elow = pv->GetLowEdgeEnergy(0 1046 G4double elow = pv->GetLowEdgeEnergy(0);
1047 G4double ehigh = pv->GetLowEdgeEnergy(n 1047 G4double ehigh = pv->GetLowEdgeEnergy(nbins-1);
1048 G4double rlow = pv->GetValue(elow, b); 1048 G4double rlow = pv->GetValue(elow, b);
1049 G4double rhigh = pv->GetValue(ehigh, b) 1049 G4double rhigh = pv->GetValue(ehigh, b);
1050 1050
1051 if (rlow <DBL_MIN) rlow = 1.e-8; 1051 if (rlow <DBL_MIN) rlow = 1.e-8;
1052 if (rhigh > 1.e16) rhigh = 1.e16; 1052 if (rhigh > 1.e16) rhigh = 1.e16;
1053 if (rhigh < 1.e-8) rhigh =1.e-8; 1053 if (rhigh < 1.e-8) rhigh =1.e-8;
1054 G4double tmpTrick = rhigh/rlow; 1054 G4double tmpTrick = rhigh/rlow;
1055 1055
1056 //std::cout << nbins << ", elow " << el 1056 //std::cout << nbins << ", elow " << elow << ", ehigh " << ehigh
1057 // << ", rlow " << rlow << ", rhigh 1057 // << ", rlow " << rlow << ", rhigh " << rhigh
1058 // << ", trick " << tmpT 1058 // << ", trick " << tmpTrick << std::endl;
1059 1059
1060 if (tmpTrick <= 0. || tmpTrick < DBL_MI 1060 if (tmpTrick <= 0. || tmpTrick < DBL_MIN) tmpTrick = 1.e-8;
1061 if (tmpTrick > 1.e16) tmpTrick = 1.e16; 1061 if (tmpTrick > 1.e16) tmpTrick = 1.e16;
1062 1062
1063 rhigh *= std::exp(std::log(tmpTrick)/(( 1063 rhigh *= std::exp(std::log(tmpTrick)/((G4double)(nbins-1)));
1064 1064
1065 G4PhysicsLogVector* v = new G4PhysicsLo 1065 G4PhysicsLogVector* v = new G4PhysicsLogVector(rlow, rhigh, nbins);
1066 1066
1067 v->PutValue(0,elow); 1067 v->PutValue(0,elow);
1068 G4double energy1 = elow; 1068 G4double energy1 = elow;
1069 G4double range1 = rlow; 1069 G4double range1 = rlow;
1070 G4double energy2 = elow; 1070 G4double energy2 = elow;
1071 G4double range2 = rlow; 1071 G4double range2 = rlow;
1072 std::size_t ilow = 0; 1072 std::size_t ilow = 0;
1073 std::size_t ihigh; 1073 std::size_t ihigh;
1074 1074
1075 for (std::size_t j=1; j<nbins; ++j) { 1075 for (std::size_t j=1; j<nbins; ++j) {
1076 1076
1077 G4double range = v->GetLowEdgeEnergy(j); 1077 G4double range = v->GetLowEdgeEnergy(j);
1078 1078
1079 for (ihigh=ilow+1; ihigh<nbins; ihigh++) { 1079 for (ihigh=ilow+1; ihigh<nbins; ihigh++) {
1080 energy2 = pv->GetLowEdgeEnergy(ihigh); 1080 energy2 = pv->GetLowEdgeEnergy(ihigh);
1081 range2 = pv->GetValue(energy2, b); 1081 range2 = pv->GetValue(energy2, b);
1082 if(range2 >= range || ihigh == nbins-1) { 1082 if(range2 >= range || ihigh == nbins-1) {
1083 ilow = ihigh - 1; 1083 ilow = ihigh - 1;
1084 energy1 = pv->GetLowEdgeEnergy(ilow); 1084 energy1 = pv->GetLowEdgeEnergy(ilow);
1085 range1 = pv->GetValue(energy1, b); 1085 range1 = pv->GetValue(energy1, b);
1086 break; 1086 break;
1087 } 1087 }
1088 } 1088 }
1089 1089
1090 G4double e = std::log(energy1) + std::log(e 1090 G4double e = std::log(energy1) + std::log(energy2/energy1)*std::log(range/range1)/std::log(range2/range1);
1091 1091
1092 v->PutValue(j,std::exp(e)); 1092 v->PutValue(j,std::exp(e));
1093 } 1093 }
1094 theInverseRangeTable->insert(v); 1094 theInverseRangeTable->insert(v);
1095 1095
1096 } 1096 }
1097 } 1097 }
1098 1098
1099 //------------------------------------------- 1099 //--------------------------------------------------------------------------------
1100 1100
1101 void G4hRDEnergyLoss::InvertRangeVector(G4int 1101 void G4hRDEnergyLoss::InvertRangeVector(G4int materialIndex,
1102 G4PhysicsLogVector* aVector) 1102 G4PhysicsLogVector* aVector)
1103 { 1103 {
1104 // invert range vector for a material 1104 // invert range vector for a material
1105 1105
1106 G4double LowEdgeRange,A,B,C,discr,KineticEn 1106 G4double LowEdgeRange,A,B,C,discr,KineticEnergy ;
1107 G4double Tbin = LowestKineticEnergy/RTable 1107 G4double Tbin = LowestKineticEnergy/RTable ;
1108 G4double rangebin = 0.0 ; 1108 G4double rangebin = 0.0 ;
1109 G4int binnumber = -1 ; 1109 G4int binnumber = -1 ;
1110 G4bool isOut ; 1110 G4bool isOut ;
1111 1111
1112 1112
1113 //loop for range values 1113 //loop for range values
1114 for( G4int i=0; i<TotBin; i++) 1114 for( G4int i=0; i<TotBin; i++)
1115 { 1115 {
1116 LowEdgeRange = aVector->GetLowEdgeEnerg 1116 LowEdgeRange = aVector->GetLowEdgeEnergy(i) ; //i.e. GetLowEdgeValue(i)
1117 1117
1118 if( rangebin < LowEdgeRange ) 1118 if( rangebin < LowEdgeRange )
1119 { 1119 {
1120 do 1120 do
1121 { 1121 {
1122 binnumber += 1 ; 1122 binnumber += 1 ;
1123 Tbin *= RTable ; 1123 Tbin *= RTable ;
1124 rangebin = (*theRangeTable)(materialI 1124 rangebin = (*theRangeTable)(materialIndex)->GetValue(Tbin,isOut) ;
1125 } 1125 }
1126 while ((rangebin < LowEdgeRange) && (binn 1126 while ((rangebin < LowEdgeRange) && (binnumber < TotBin )) ;
1127 } 1127 }
1128 1128
1129 if(binnumber == 0) 1129 if(binnumber == 0)
1130 KineticEnergy = LowestKineticEnergy ; 1130 KineticEnergy = LowestKineticEnergy ;
1131 else if(binnumber == TotBin-1) 1131 else if(binnumber == TotBin-1)
1132 KineticEnergy = HighestKineticEnergy ; 1132 KineticEnergy = HighestKineticEnergy ;
1133 else 1133 else
1134 { 1134 {
1135 A = (*(*theRangeCoeffATable)(materialInde 1135 A = (*(*theRangeCoeffATable)(materialIndex))(binnumber-1) ;
1136 B = (*(*theRangeCoeffBTable)(materialInde 1136 B = (*(*theRangeCoeffBTable)(materialIndex))(binnumber-1) ;
1137 C = (*(*theRangeCoeffCTable)(materialInde 1137 C = (*(*theRangeCoeffCTable)(materialIndex))(binnumber-1) ;
1138 if(A==0.) 1138 if(A==0.)
1139 KineticEnergy = (LowEdgeRange -C )/B ; 1139 KineticEnergy = (LowEdgeRange -C )/B ;
1140 else 1140 else
1141 { 1141 {
1142 discr = B*B - 4.*A*(C-LowEdgeRange); 1142 discr = B*B - 4.*A*(C-LowEdgeRange);
1143 discr = discr>0. ? std::sqrt(discr) : 1143 discr = discr>0. ? std::sqrt(discr) : 0.;
1144 KineticEnergy = 0.5*(discr-B)/A ; 1144 KineticEnergy = 0.5*(discr-B)/A ;
1145 } 1145 }
1146 } 1146 }
1147 1147
1148 aVector->PutValue(i,KineticEnergy) ; 1148 aVector->PutValue(i,KineticEnergy) ;
1149 } 1149 }
1150 } 1150 }
1151 1151
1152 //------------------------------------------- 1152 //------------------------------------------------------------------------------
1153 1153
1154 G4bool G4hRDEnergyLoss::CutsWhereModified() 1154 G4bool G4hRDEnergyLoss::CutsWhereModified()
1155 { 1155 {
1156 G4bool wasModified = false; 1156 G4bool wasModified = false;
1157 const G4ProductionCutsTable* theCoupleTable 1157 const G4ProductionCutsTable* theCoupleTable=
1158 G4ProductionCutsTable::GetProductionCutsT 1158 G4ProductionCutsTable::GetProductionCutsTable();
1159 G4int numOfCouples = (G4int)theCoupleTable- 1159 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
1160 1160
1161 for (G4int j=0; j<numOfCouples; ++j){ 1161 for (G4int j=0; j<numOfCouples; ++j){
1162 if (theCoupleTable->GetMaterialCutsCouple 1162 if (theCoupleTable->GetMaterialCutsCouple(j)->IsRecalcNeeded()) {
1163 wasModified = true; 1163 wasModified = true;
1164 break; 1164 break;
1165 } 1165 }
1166 } 1166 }
1167 return wasModified; 1167 return wasModified;
1168 } 1168 }
1169 1169
1170 //------------------------------------------- 1170 //-----------------------------------------------------------------------
1171 //--- --- --- --- --- --- --- --- --- --- --- 1171 //--- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
1172 1172
1173 //G4bool G4hRDEnergyLoss::IsApplicable(const 1173 //G4bool G4hRDEnergyLoss::IsApplicable(const G4ParticleDefinition&
1174 // 1174 // particle)
1175 //{ 1175 //{
1176 // return((particle.GetPDGCharge()!= 0.) && 1176 // return((particle.GetPDGCharge()!= 0.) && (particle.GetLeptonNumber() == 0));
1177 //} 1177 //}
1178 1178
1179 //G4double G4hRDEnergyLoss::GetContinuousStep 1179 //G4double G4hRDEnergyLoss::GetContinuousStepLimit(
1180 // 1180 // const G4Track& track,
1181 // 1181 // G4double,
1182 // 1182 // G4double currentMinimumStep,
1183 // 1183 // G4double&)
1184 //{ 1184 //{
1185 // 1185 //
1186 // G4double Step = 1186 // G4double Step =
1187 // GetConstraints(track.GetDynamicParticle 1187 // GetConstraints(track.GetDynamicParticle(),track.GetMaterial()) ;
1188 // 1188 //
1189 // if((Step>0.0)&&(Step<currentMinimumStep)) 1189 // if((Step>0.0)&&(Step<currentMinimumStep))
1190 // currentMinimumStep = Step ; 1190 // currentMinimumStep = Step ;
1191 // 1191 //
1192 // return Step ; 1192 // return Step ;
1193 //} 1193 //}
1194 1194