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