| 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 // G4RDHadronIonisation 29 // G4RDHadronIonisation
30 // 30 //
>> 31 // $Id: G4hImpactIonisation.cc,v 1.4 2010/11/25 19:49:43 pia Exp $
>> 32 // GEANT4 tag $Name: geant4-09-04 $
31 // 33 //
32 // Author: Maria Grazia Pia (MariaGrazia.Pia@g 34 // Author: Maria Grazia Pia (MariaGrazia.Pia@ge.infn.it)
33 // 35 //
34 // 08 Sep 2008 - MGP - Created (initially base 36 // 08 Sep 2008 - MGP - Created (initially based on G4hLowEnergyIonisation)
35 // Added PIXE capabilities 37 // Added PIXE capabilities
36 // Partial clean-up of the 38 // Partial clean-up of the implementation (more needed)
37 // Calculation of Microsco 39 // Calculation of MicroscopicCrossSection delegated to specialised cla// Documentation available in:
38 // M.G. Pia et al., PIXE Simulation With Geant 40 // M.G. Pia et al., PIXE Simulation With Geant4,
39 // IEEE Trans. Nucl. Sci., vol. 56, no. 6, pp. 41 // IEEE Trans. Nucl. Sci., vol. 56, no. 6, pp. 3614-3649, Dec. 2009.
40 42
41 // 43 //
42 // ------------------------------------------- 44 // ------------------------------------------------------------
43 45
44 #include "G4hImpactIonisation.hh" 46 #include "G4hImpactIonisation.hh"
45 #include "globals.hh" 47 #include "globals.hh"
46 #include "G4ios.hh" 48 #include "G4ios.hh"
47 #include "Randomize.hh" 49 #include "Randomize.hh"
48 #include "G4PhysicalConstants.hh" <<
49 #include "G4SystemOfUnits.hh" <<
50 #include "G4Poisson.hh" 50 #include "G4Poisson.hh"
51 #include "G4UnitsTable.hh" 51 #include "G4UnitsTable.hh"
52 #include "G4EnergyLossTables.hh" 52 #include "G4EnergyLossTables.hh"
53 #include "G4Material.hh" 53 #include "G4Material.hh"
54 #include "G4DynamicParticle.hh" 54 #include "G4DynamicParticle.hh"
55 #include "G4ParticleDefinition.hh" 55 #include "G4ParticleDefinition.hh"
56 #include "G4AtomicDeexcitation.hh" 56 #include "G4AtomicDeexcitation.hh"
57 #include "G4AtomicTransitionManager.hh" 57 #include "G4AtomicTransitionManager.hh"
58 #include "G4PixeCrossSectionHandler.hh" 58 #include "G4PixeCrossSectionHandler.hh"
59 #include "G4IInterpolator.hh" 59 #include "G4IInterpolator.hh"
60 #include "G4LogLogInterpolator.hh" 60 #include "G4LogLogInterpolator.hh"
61 #include "G4Gamma.hh" 61 #include "G4Gamma.hh"
62 #include "G4Electron.hh" 62 #include "G4Electron.hh"
63 #include "G4Proton.hh" 63 #include "G4Proton.hh"
64 #include "G4ProcessManager.hh" 64 #include "G4ProcessManager.hh"
65 #include "G4ProductionCutsTable.hh" 65 #include "G4ProductionCutsTable.hh"
66 #include "G4PhysicsLogVector.hh" 66 #include "G4PhysicsLogVector.hh"
67 #include "G4PhysicsLinearVector.hh" 67 #include "G4PhysicsLinearVector.hh"
68 68
69 #include "G4VLowEnergyModel.hh" 69 #include "G4VLowEnergyModel.hh"
70 #include "G4hNuclearStoppingModel.hh" 70 #include "G4hNuclearStoppingModel.hh"
71 #include "G4hBetheBlochModel.hh" 71 #include "G4hBetheBlochModel.hh"
72 #include "G4hParametrisedLossModel.hh" 72 #include "G4hParametrisedLossModel.hh"
73 #include "G4QAOLowEnergyLoss.hh" 73 #include "G4QAOLowEnergyLoss.hh"
74 #include "G4hIonEffChargeSquare.hh" 74 #include "G4hIonEffChargeSquare.hh"
75 #include "G4IonChuFluctuationModel.hh" 75 #include "G4IonChuFluctuationModel.hh"
76 #include "G4IonYangFluctuationModel.hh" 76 #include "G4IonYangFluctuationModel.hh"
77 77
78 #include "G4MaterialCutsCouple.hh" 78 #include "G4MaterialCutsCouple.hh"
79 #include "G4Track.hh" 79 #include "G4Track.hh"
80 #include "G4Step.hh" 80 #include "G4Step.hh"
81 81
82 G4hImpactIonisation::G4hImpactIonisation(const 82 G4hImpactIonisation::G4hImpactIonisation(const G4String& processName)
83 : G4hRDEnergyLoss(processName), 83 : G4hRDEnergyLoss(processName),
84 betheBlochModel(0), 84 betheBlochModel(0),
85 protonModel(0), 85 protonModel(0),
86 antiprotonModel(0), 86 antiprotonModel(0),
87 theIonEffChargeModel(0), 87 theIonEffChargeModel(0),
88 theNuclearStoppingModel(0), 88 theNuclearStoppingModel(0),
89 theIonChuFluctuationModel(0), 89 theIonChuFluctuationModel(0),
90 theIonYangFluctuationModel(0), 90 theIonYangFluctuationModel(0),
91 protonTable("ICRU_R49p"), 91 protonTable("ICRU_R49p"),
92 antiprotonTable("ICRU_R49p"), 92 antiprotonTable("ICRU_R49p"),
93 theNuclearTable("ICRU_R49"), 93 theNuclearTable("ICRU_R49"),
94 nStopping(true), 94 nStopping(true),
95 theBarkas(true), 95 theBarkas(true),
96 theMeanFreePathTable(0), 96 theMeanFreePathTable(0),
97 paramStepLimit (0.005), 97 paramStepLimit (0.005),
98 pixeCrossSectionHandler(0) 98 pixeCrossSectionHandler(0)
99 { 99 {
100 InitializeMe(); 100 InitializeMe();
101 } 101 }
102 102
103 103
104 104
105 void G4hImpactIonisation::InitializeMe() 105 void G4hImpactIonisation::InitializeMe()
106 { 106 {
107 LowestKineticEnergy = 10.0*eV ; 107 LowestKineticEnergy = 10.0*eV ;
108 HighestKineticEnergy = 100.0*GeV ; 108 HighestKineticEnergy = 100.0*GeV ;
109 MinKineticEnergy = 10.0*eV ; 109 MinKineticEnergy = 10.0*eV ;
110 TotBin = 360 ; 110 TotBin = 360 ;
111 protonLowEnergy = 1.*keV ; 111 protonLowEnergy = 1.*keV ;
112 protonHighEnergy = 100.*MeV ; 112 protonHighEnergy = 100.*MeV ;
113 antiprotonLowEnergy = 25.*keV ; 113 antiprotonLowEnergy = 25.*keV ;
114 antiprotonHighEnergy = 2.*MeV ; 114 antiprotonHighEnergy = 2.*MeV ;
115 minGammaEnergy = 100 * eV; 115 minGammaEnergy = 100 * eV;
116 minElectronEnergy = 250.* eV; 116 minElectronEnergy = 250.* eV;
117 verboseLevel = 0; 117 verboseLevel = 0;
118 118
119 // Min and max energy of incident particle f 119 // Min and max energy of incident particle for the calculation of shell cross sections
120 // for PIXE generation 120 // for PIXE generation
121 eMinPixe = 1.* keV; 121 eMinPixe = 1.* keV;
122 eMaxPixe = 200. * MeV; 122 eMaxPixe = 200. * MeV;
123 123
124 G4String defaultPixeModel("ecpssr"); 124 G4String defaultPixeModel("ecpssr");
125 modelK = defaultPixeModel; 125 modelK = defaultPixeModel;
126 modelL = defaultPixeModel; 126 modelL = defaultPixeModel;
127 modelM = defaultPixeModel; 127 modelM = defaultPixeModel;
128 <<
129 // Calculated on the fly, but should be init <<
130 fdEdx = 0.0; <<
131 fRangeNow = 0.0; <<
132 charge = 0.0; <<
133 chargeSquare = 0.0; <<
134 initialMass = 0.0; <<
135 fBarkas = 0.0; <<
136 } 128 }
137 129
138 130
139 G4hImpactIonisation::~G4hImpactIonisation() 131 G4hImpactIonisation::~G4hImpactIonisation()
140 { 132 {
141 if (theMeanFreePathTable) 133 if (theMeanFreePathTable)
142 { 134 {
143 theMeanFreePathTable->clearAndDestroy(); 135 theMeanFreePathTable->clearAndDestroy();
144 delete theMeanFreePathTable; 136 delete theMeanFreePathTable;
145 } 137 }
146 138
147 if (betheBlochModel) delete betheBlochModel; 139 if (betheBlochModel) delete betheBlochModel;
148 if (protonModel) delete protonModel; 140 if (protonModel) delete protonModel;
149 if (antiprotonModel) delete antiprotonModel; 141 if (antiprotonModel) delete antiprotonModel;
150 if (theNuclearStoppingModel) delete theNucle 142 if (theNuclearStoppingModel) delete theNuclearStoppingModel;
151 if (theIonEffChargeModel) delete theIonEffCh 143 if (theIonEffChargeModel) delete theIonEffChargeModel;
152 if (theIonChuFluctuationModel) delete theIon 144 if (theIonChuFluctuationModel) delete theIonChuFluctuationModel;
153 if (theIonYangFluctuationModel) delete theIo 145 if (theIonYangFluctuationModel) delete theIonYangFluctuationModel;
154 146
155 delete pixeCrossSectionHandler; 147 delete pixeCrossSectionHandler;
156 148
157 // ---- MGP ---- The following is to be chec 149 // ---- MGP ---- The following is to be checked
158 // if (shellVacancy) delete shellVacancy; 150 // if (shellVacancy) delete shellVacancy;
159 151
160 cutForDelta.clear(); 152 cutForDelta.clear();
161 } 153 }
162 154
163 // ------------------------------------------- 155 // --------------------------------------------------------------------
164 void G4hImpactIonisation::SetElectronicStoppin 156 void G4hImpactIonisation::SetElectronicStoppingPowerModel(const G4ParticleDefinition* particle,
165 const G4String& dedxTable) 157 const G4String& dedxTable)
166 // This method defines the ionisation parametr 158 // This method defines the ionisation parametrisation method via its name
167 { 159 {
168 if (particle->GetPDGCharge() > 0 ) 160 if (particle->GetPDGCharge() > 0 )
169 { 161 {
170 // Positive charge 162 // Positive charge
171 SetProtonElectronicStoppingPowerModel(de 163 SetProtonElectronicStoppingPowerModel(dedxTable) ;
172 } 164 }
173 else 165 else
174 { 166 {
175 // Antiprotons 167 // Antiprotons
176 SetAntiProtonElectronicStoppingPowerMode 168 SetAntiProtonElectronicStoppingPowerModel(dedxTable) ;
177 } 169 }
178 } 170 }
179 171
180 172
181 // ------------------------------------------- 173 // --------------------------------------------------------------------
182 void G4hImpactIonisation::InitializeParametris 174 void G4hImpactIonisation::InitializeParametrisation()
183 175
184 { 176 {
185 // Define models for parametrisation of elec 177 // Define models for parametrisation of electronic energy losses
186 betheBlochModel = new G4hBetheBlochModel("Be 178 betheBlochModel = new G4hBetheBlochModel("Bethe-Bloch") ;
187 protonModel = new G4hParametrisedLossModel(p 179 protonModel = new G4hParametrisedLossModel(protonTable) ;
188 protonHighEnergy = std::min(protonHighEnergy 180 protonHighEnergy = std::min(protonHighEnergy,protonModel->HighEnergyLimit(0, 0));
189 antiprotonModel = new G4QAOLowEnergyLoss(ant 181 antiprotonModel = new G4QAOLowEnergyLoss(antiprotonTable) ;
190 theNuclearStoppingModel = new G4hNuclearStop 182 theNuclearStoppingModel = new G4hNuclearStoppingModel(theNuclearTable) ;
191 theIonEffChargeModel = new G4hIonEffChargeSq 183 theIonEffChargeModel = new G4hIonEffChargeSquare("Ziegler1988") ;
192 theIonChuFluctuationModel = new G4IonChuFluc 184 theIonChuFluctuationModel = new G4IonChuFluctuationModel("Chu") ;
193 theIonYangFluctuationModel = new G4IonYangFl 185 theIonYangFluctuationModel = new G4IonYangFluctuationModel("Yang") ;
194 } 186 }
195 187
196 188
197 // ------------------------------------------- 189 // --------------------------------------------------------------------
198 void G4hImpactIonisation::BuildPhysicsTable(co 190 void G4hImpactIonisation::BuildPhysicsTable(const G4ParticleDefinition& particleDef)
199 191
200 // just call BuildLossTable+BuildLambdaTable 192 // just call BuildLossTable+BuildLambdaTable
201 { 193 {
202 194
203 // Verbose print-out 195 // Verbose print-out
204 if(verboseLevel > 0) 196 if(verboseLevel > 0)
205 { 197 {
206 G4cout << "G4hImpactIonisation::BuildPhy 198 G4cout << "G4hImpactIonisation::BuildPhysicsTable for "
207 << particleDef.GetParticleName() 199 << particleDef.GetParticleName()
208 << " mass(MeV)= " << particleDef.GetPDG 200 << " mass(MeV)= " << particleDef.GetPDGMass()/MeV
209 << " charge= " << particleDef.GetPDGCha 201 << " charge= " << particleDef.GetPDGCharge()/eplus
210 << " type= " << particleDef.GetParticle 202 << " type= " << particleDef.GetParticleType()
211 << G4endl; 203 << G4endl;
212 204
213 if(verboseLevel > 1) 205 if(verboseLevel > 1)
214 { 206 {
215 G4ProcessVector* pv = particleDef.GetProce 207 G4ProcessVector* pv = particleDef.GetProcessManager()->GetProcessList();
216 208
217 G4cout << " 0: " << (*pv)[0]->GetProcessNa 209 G4cout << " 0: " << (*pv)[0]->GetProcessName() << " " << (*pv)[0]
218 << " 1: " << (*pv)[1]->GetProcessName() < 210 << " 1: " << (*pv)[1]->GetProcessName() << " " << (*pv)[1]
219 // << " 2: " << (*pv)[2]->GetProc 211 // << " 2: " << (*pv)[2]->GetProcessName() << " " << (*pv)[2]
220 << G4endl; 212 << G4endl;
221 G4cout << "ionModel= " << theIonEffChargeM 213 G4cout << "ionModel= " << theIonEffChargeModel
222 << " MFPtable= " << theMeanFreePathTable 214 << " MFPtable= " << theMeanFreePathTable
223 << " iniMass= " << initialMass 215 << " iniMass= " << initialMass
224 << G4endl; 216 << G4endl;
225 } 217 }
226 } 218 }
227 // End of verbose print-out 219 // End of verbose print-out
228 220
229 if (particleDef.GetParticleType() == "nucleu 221 if (particleDef.GetParticleType() == "nucleus" &&
230 particleDef.GetParticleName() != "Generi 222 particleDef.GetParticleName() != "GenericIon" &&
231 particleDef.GetParticleSubType() == "gen 223 particleDef.GetParticleSubType() == "generic")
232 { 224 {
233 225
234 G4EnergyLossTables::Register(&particleDe 226 G4EnergyLossTables::Register(&particleDef,
235 theDEDXpTable, 227 theDEDXpTable,
236 theRangepTable, 228 theRangepTable,
237 theInverseRangepTable, 229 theInverseRangepTable,
238 theLabTimepTable, 230 theLabTimepTable,
239 theProperTimepTable, 231 theProperTimepTable,
240 LowestKineticEnergy, HighestKinetic 232 LowestKineticEnergy, HighestKineticEnergy,
241 proton_mass_c2/particleDef.GetPDGMa 233 proton_mass_c2/particleDef.GetPDGMass(),
242 TotBin); 234 TotBin);
243 235
244 return; 236 return;
245 } 237 }
246 238
247 if( !CutsWhereModified() && theLossTable) re 239 if( !CutsWhereModified() && theLossTable) return;
248 240
249 InitializeParametrisation() ; 241 InitializeParametrisation() ;
250 G4Proton* proton = G4Proton::Proton(); 242 G4Proton* proton = G4Proton::Proton();
251 G4AntiProton* antiproton = G4AntiProton::Ant 243 G4AntiProton* antiproton = G4AntiProton::AntiProton();
252 244
253 charge = particleDef.GetPDGCharge() / eplus; 245 charge = particleDef.GetPDGCharge() / eplus;
254 chargeSquare = charge*charge ; 246 chargeSquare = charge*charge ;
255 247
256 const G4ProductionCutsTable* theCoupleTable= 248 const G4ProductionCutsTable* theCoupleTable=
257 G4ProductionCutsTable::GetProductionCutsTa 249 G4ProductionCutsTable::GetProductionCutsTable();
258 G4int numOfCouples = (G4int)theCoupleTable-> << 250 size_t numOfCouples = theCoupleTable->GetTableSize();
259 251
260 cutForDelta.clear(); 252 cutForDelta.clear();
261 cutForGamma.clear(); 253 cutForGamma.clear();
262 254
263 for (G4int j=0; j<numOfCouples; ++j) { << 255 for (size_t j=0; j<numOfCouples; j++) {
264 256
265 // get material parameters needed for the 257 // get material parameters needed for the energy loss calculation
266 const G4MaterialCutsCouple* couple = theCo 258 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
267 const G4Material* material= couple->GetMat 259 const G4Material* material= couple->GetMaterial();
268 260
269 // the cut cannot be below lowest limit 261 // the cut cannot be below lowest limit
270 G4double tCut = (*(theCoupleTable->GetEner 262 G4double tCut = (*(theCoupleTable->GetEnergyCutsVector(1)))[j];
271 if(tCut > HighestKineticEnergy) tCut = Hig 263 if(tCut > HighestKineticEnergy) tCut = HighestKineticEnergy;
272 264
273 G4double excEnergy = material->GetIonisati 265 G4double excEnergy = material->GetIonisation()->GetMeanExcitationEnergy();
274 266
275 tCut = std::max(tCut,excEnergy); 267 tCut = std::max(tCut,excEnergy);
276 cutForDelta.push_back(tCut); 268 cutForDelta.push_back(tCut);
277 269
278 // the cut cannot be below lowest limit 270 // the cut cannot be below lowest limit
279 tCut = (*(theCoupleTable->GetEnergyCutsVec 271 tCut = (*(theCoupleTable->GetEnergyCutsVector(0)))[j];
280 if(tCut > HighestKineticEnergy) tCut = Hig 272 if(tCut > HighestKineticEnergy) tCut = HighestKineticEnergy;
281 tCut = std::max(tCut,minGammaEnergy); 273 tCut = std::max(tCut,minGammaEnergy);
282 cutForGamma.push_back(tCut); 274 cutForGamma.push_back(tCut);
283 } 275 }
284 276
285 if(verboseLevel > 0) { 277 if(verboseLevel > 0) {
286 G4cout << "Cuts are defined " << G4endl; 278 G4cout << "Cuts are defined " << G4endl;
287 } 279 }
288 280
289 if(0.0 < charge) 281 if(0.0 < charge)
290 { 282 {
291 { 283 {
292 BuildLossTable(*proton) ; 284 BuildLossTable(*proton) ;
293 285
294 // The following vector has a fixed dim 286 // The following vector has a fixed dimension (see src/G4hImpactLoss.cc for more details)
295 // It happended in the past that caused 287 // It happended in the past that caused memory corruption errors. The problem is still pending, even if temporary solved
296 // G4cout << "[NOTE]: __LINE__=" << _ 288 // G4cout << "[NOTE]: __LINE__=" << __LINE__ << ", particleDef=" << particleDef.GetParticleName() << ", proton=" << proton << ", theLossTable=" << theLossTable << ", CounterOfpProcess=" << CounterOfpProcess << G4endl;
297 289
298 RecorderOfpProcess[CounterOfpProcess] 290 RecorderOfpProcess[CounterOfpProcess] = theLossTable ;
299 CounterOfpProcess++; 291 CounterOfpProcess++;
300 } 292 }
301 } else { 293 } else {
302 { 294 {
303 BuildLossTable(*antiproton) ; 295 BuildLossTable(*antiproton) ;
304 296
305 // The following vector has a fixed 297 // The following vector has a fixed dimension (see src/G4hImpactLoss.cc for more details)
306 // It happended in the past that ca 298 // It happended in the past that caused memory corruption errors. The problem is still pending, even if temporary solved
307 // G4cout << "[NOTE]: __LINE__=" 299 // G4cout << "[NOTE]: __LINE__=" << __LINE__ << ", particleDef=" << particleDef.GetParticleName() << ", antiproton=" << antiproton << ", theLossTable=" << theLossTable << ", CounterOfpbarProcess=" << CounterOfpbarProcess << G4endl;
308 300
309 RecorderOfpbarProcess[CounterOfpbarProce 301 RecorderOfpbarProcess[CounterOfpbarProcess] = theLossTable ;
310 CounterOfpbarProcess++; 302 CounterOfpbarProcess++;
311 } 303 }
312 } 304 }
313 305
314 if(verboseLevel > 0) { 306 if(verboseLevel > 0) {
315 G4cout << "G4hImpactIonisation::BuildPhysi 307 G4cout << "G4hImpactIonisation::BuildPhysicsTable: "
316 << "Loss table is built " 308 << "Loss table is built "
317 // << theLossTable 309 // << theLossTable
318 << G4endl; 310 << G4endl;
319 } 311 }
320 312
321 BuildLambdaTable(particleDef) ; 313 BuildLambdaTable(particleDef) ;
322 // BuildDataForFluorescence(particleDef); 314 // BuildDataForFluorescence(particleDef);
323 315
324 if(verboseLevel > 1) { 316 if(verboseLevel > 1) {
325 G4cout << (*theMeanFreePathTable) << G4end 317 G4cout << (*theMeanFreePathTable) << G4endl;
326 } 318 }
327 319
328 if(verboseLevel > 0) { 320 if(verboseLevel > 0) {
329 G4cout << "G4hImpactIonisation::BuildPhysi 321 G4cout << "G4hImpactIonisation::BuildPhysicsTable: "
330 << "DEDX table will be built " 322 << "DEDX table will be built "
331 // << theDEDXpTable << " " << theDED 323 // << theDEDXpTable << " " << theDEDXpbarTable
332 // << " " << theRangepTable << " " < 324 // << " " << theRangepTable << " " << theRangepbarTable
333 << G4endl; 325 << G4endl;
334 } 326 }
335 327
336 BuildDEDXTable(particleDef) ; 328 BuildDEDXTable(particleDef) ;
337 329
338 if(verboseLevel > 1) { 330 if(verboseLevel > 1) {
339 G4cout << (*theDEDXpTable) << G4endl; 331 G4cout << (*theDEDXpTable) << G4endl;
340 } 332 }
341 333
342 if((&particleDef == proton) || (&particleDe 334 if((&particleDef == proton) || (&particleDef == antiproton)) PrintInfoDefinition() ;
343 335
344 if(verboseLevel > 0) { 336 if(verboseLevel > 0) {
345 G4cout << "G4hImpactIonisation::BuildPhysi 337 G4cout << "G4hImpactIonisation::BuildPhysicsTable: end for "
346 << particleDef.GetParticleName() << 338 << particleDef.GetParticleName() << G4endl;
347 } 339 }
348 } 340 }
349 341
350 342
351 343
352 344
353 345
354 // ------------------------------------------- 346 // --------------------------------------------------------------------
355 void G4hImpactIonisation::BuildLossTable(const 347 void G4hImpactIonisation::BuildLossTable(const G4ParticleDefinition& particleDef)
356 { 348 {
>> 349
357 // Initialisation 350 // Initialisation
358 G4double lowEdgeEnergy , ionloss, ionlossBB, 351 G4double lowEdgeEnergy , ionloss, ionlossBB, paramB ;
359 //G4double lowEnergy, highEnergy; << 352 G4double lowEnergy, highEnergy;
360 G4double highEnergy; <<
361 G4Proton* proton = G4Proton::Proton(); 353 G4Proton* proton = G4Proton::Proton();
362 354
363 if (particleDef == *proton) 355 if (particleDef == *proton)
364 { 356 {
365 //lowEnergy = protonLowEnergy ; << 357 lowEnergy = protonLowEnergy ;
366 highEnergy = protonHighEnergy ; 358 highEnergy = protonHighEnergy ;
367 charge = 1. ; 359 charge = 1. ;
368 } 360 }
369 else 361 else
370 { 362 {
371 //lowEnergy = antiprotonLowEnergy ; << 363 lowEnergy = antiprotonLowEnergy ;
372 highEnergy = antiprotonHighEnergy ; 364 highEnergy = antiprotonHighEnergy ;
373 charge = -1. ; 365 charge = -1. ;
374 } 366 }
375 chargeSquare = 1. ; 367 chargeSquare = 1. ;
376 368
377 const G4ProductionCutsTable* theCoupleTable= 369 const G4ProductionCutsTable* theCoupleTable=
378 G4ProductionCutsTable::GetProductionCutsTa 370 G4ProductionCutsTable::GetProductionCutsTable();
379 G4int numOfCouples = (G4int)theCoupleTable-> << 371 size_t numOfCouples = theCoupleTable->GetTableSize();
380 372
381 if ( theLossTable) 373 if ( theLossTable)
382 { 374 {
383 theLossTable->clearAndDestroy(); 375 theLossTable->clearAndDestroy();
384 delete theLossTable; 376 delete theLossTable;
385 } 377 }
386 378
387 theLossTable = new G4PhysicsTable(numOfCoupl 379 theLossTable = new G4PhysicsTable(numOfCouples);
388 380
389 // loop for materials 381 // loop for materials
390 for (G4int j=0; j<numOfCouples; ++j) { << 382 for (size_t j=0; j<numOfCouples; j++) {
391 383
392 // create physics vector and fill it 384 // create physics vector and fill it
393 G4PhysicsLogVector* aVector = new G4Physic 385 G4PhysicsLogVector* aVector = new G4PhysicsLogVector(LowestKineticEnergy,
394 386 HighestKineticEnergy,
395 387 TotBin);
396 388
397 // get material parameters needed for the 389 // get material parameters needed for the energy loss calculation
398 const G4MaterialCutsCouple* couple = theCo 390 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
399 const G4Material* material= couple->GetMat 391 const G4Material* material= couple->GetMaterial();
400 392
401 if ( charge > 0.0 ) { 393 if ( charge > 0.0 ) {
402 ionloss = ProtonParametrisedDEDX(couple, 394 ionloss = ProtonParametrisedDEDX(couple,highEnergy) ;
403 } else { 395 } else {
404 ionloss = AntiProtonParametrisedDEDX(cou 396 ionloss = AntiProtonParametrisedDEDX(couple,highEnergy) ;
405 } 397 }
406 398
407 ionlossBB = betheBlochModel->TheValue(&par 399 ionlossBB = betheBlochModel->TheValue(&particleDef,material,highEnergy) ;
408 ionlossBB -= DeltaRaysEnergy(couple,highEn 400 ionlossBB -= DeltaRaysEnergy(couple,highEnergy,proton_mass_c2) ;
409 401
410 402
411 paramB = ionloss/ionlossBB - 1.0 ; 403 paramB = ionloss/ionlossBB - 1.0 ;
412 404
413 // now comes the loop for the kinetic ener 405 // now comes the loop for the kinetic energy values
414 for (G4int i = 0 ; i < TotBin ; i++) { 406 for (G4int i = 0 ; i < TotBin ; i++) {
415 lowEdgeEnergy = aVector->GetLowEdgeEnerg 407 lowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ;
416 408
417 // low energy part for this material, pa 409 // low energy part for this material, parametrised energy loss formulae
418 if ( lowEdgeEnergy < highEnergy ) { 410 if ( lowEdgeEnergy < highEnergy ) {
419 411
420 if ( charge > 0.0 ) { 412 if ( charge > 0.0 ) {
421 ionloss = ProtonParametrisedDEDX(cou 413 ionloss = ProtonParametrisedDEDX(couple,lowEdgeEnergy) ;
422 } else { 414 } else {
423 ionloss = AntiProtonParametrisedDEDX 415 ionloss = AntiProtonParametrisedDEDX(couple,lowEdgeEnergy) ;
424 } 416 }
425 417
426 } else { 418 } else {
427 419
428 // high energy part for this material, 420 // high energy part for this material, Bethe-Bloch formula
429 ionloss = betheBlochModel->TheValue(pr 421 ionloss = betheBlochModel->TheValue(proton,material,
430 lowEdgeEnergy) ; 422 lowEdgeEnergy) ;
431 423
432 ionloss -= DeltaRaysEnergy(couple,lowE 424 ionloss -= DeltaRaysEnergy(couple,lowEdgeEnergy,proton_mass_c2) ;
433 425
434 ionloss *= (1.0 + paramB*highEnergy/lowEdgeE 426 ionloss *= (1.0 + paramB*highEnergy/lowEdgeEnergy) ;
435 } 427 }
436 428
437 // now put the loss into the vector 429 // now put the loss into the vector
438 if(verboseLevel > 1) { 430 if(verboseLevel > 1) {
439 G4cout << "E(MeV)= " << lowEdgeEnergy/ 431 G4cout << "E(MeV)= " << lowEdgeEnergy/MeV
440 << " dE/dx(MeV/mm)= " << ionlo 432 << " dE/dx(MeV/mm)= " << ionloss*mm/MeV
441 << " in " << material->GetName( 433 << " in " << material->GetName() << G4endl;
442 } 434 }
443 aVector->PutValue(i,ionloss) ; 435 aVector->PutValue(i,ionloss) ;
444 } 436 }
445 // Insert vector for this material into th 437 // Insert vector for this material into the table
446 theLossTable->insert(aVector) ; 438 theLossTable->insert(aVector) ;
447 } 439 }
448 } 440 }
449 441
450 442
451 443
452 // ------------------------------------------- 444 // --------------------------------------------------------------------
453 void G4hImpactIonisation::BuildLambdaTable(con 445 void G4hImpactIonisation::BuildLambdaTable(const G4ParticleDefinition& particleDef)
454 446
455 { 447 {
456 // Build mean free path tables for the delta 448 // Build mean free path tables for the delta ray production process
457 // tables are built for MATERIALS 449 // tables are built for MATERIALS
458 450
459 if(verboseLevel > 1) { 451 if(verboseLevel > 1) {
460 G4cout << "G4hImpactIonisation::BuildLambd 452 G4cout << "G4hImpactIonisation::BuildLambdaTable for "
461 << particleDef.GetParticleName() << 453 << particleDef.GetParticleName() << " is started" << G4endl;
462 } 454 }
463 455
464 456
465 G4double lowEdgeEnergy, value; 457 G4double lowEdgeEnergy, value;
466 charge = particleDef.GetPDGCharge()/eplus ; 458 charge = particleDef.GetPDGCharge()/eplus ;
467 chargeSquare = charge*charge ; 459 chargeSquare = charge*charge ;
468 initialMass = particleDef.GetPDGMass(); 460 initialMass = particleDef.GetPDGMass();
469 461
470 const G4ProductionCutsTable* theCoupleTable= 462 const G4ProductionCutsTable* theCoupleTable=
471 G4ProductionCutsTable::GetProductionCutsTa 463 G4ProductionCutsTable::GetProductionCutsTable();
472 G4int numOfCouples = (G4int)theCoupleTable-> << 464 size_t numOfCouples = theCoupleTable->GetTableSize();
473 465
474 466
475 if (theMeanFreePathTable) { 467 if (theMeanFreePathTable) {
476 theMeanFreePathTable->clearAndDestroy(); 468 theMeanFreePathTable->clearAndDestroy();
477 delete theMeanFreePathTable; 469 delete theMeanFreePathTable;
478 } 470 }
479 471
480 theMeanFreePathTable = new G4PhysicsTable(nu 472 theMeanFreePathTable = new G4PhysicsTable(numOfCouples);
481 473
482 // loop for materials 474 // loop for materials
483 475
484 for (G4int j=0 ; j < numOfCouples; ++j) { << 476 for (size_t j=0 ; j < numOfCouples; j++) {
485 477
486 //create physics vector then fill it .... 478 //create physics vector then fill it ....
487 G4PhysicsLogVector* aVector = new G4Physic 479 G4PhysicsLogVector* aVector = new G4PhysicsLogVector(LowestKineticEnergy,
488 480 HighestKineticEnergy,
489 481 TotBin);
490 482
491 // compute the (macroscopic) cross section 483 // compute the (macroscopic) cross section first
492 const G4MaterialCutsCouple* couple = theCo 484 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
493 const G4Material* material= couple->GetMat 485 const G4Material* material= couple->GetMaterial();
494 486
495 const G4ElementVector* theElementVector = 487 const G4ElementVector* theElementVector = material->GetElementVector() ;
496 const G4double* theAtomicNumDensityVector 488 const G4double* theAtomicNumDensityVector = material->GetAtomicNumDensityVector();
497 const G4int numberOfElements = (G4int)mate << 489 const G4int numberOfElements = material->GetNumberOfElements() ;
498 490
499 // get the electron kinetic energy cut for 491 // get the electron kinetic energy cut for the actual material,
500 // it will be used in ComputeMicroscopicC 492 // it will be used in ComputeMicroscopicCrossSection
501 // ( it is the SAME for ALL the ELEMENTS i 493 // ( it is the SAME for ALL the ELEMENTS in THIS MATERIAL )
502 // ------------------------------------- 494 // ------------------------------------------------------
503 495
504 G4double deltaCut = cutForDelta[j]; 496 G4double deltaCut = cutForDelta[j];
505 497
506 for ( G4int i = 0 ; i < TotBin ; i++ ) { 498 for ( G4int i = 0 ; i < TotBin ; i++ ) {
507 lowEdgeEnergy = aVector->GetLowEdgeEnerg 499 lowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ;
508 G4double sigma = 0.0 ; 500 G4double sigma = 0.0 ;
509 G4int Z; 501 G4int Z;
510 502
511 for (G4int iel=0; iel<numberOfElements; 503 for (G4int iel=0; iel<numberOfElements; iel++ )
512 { 504 {
513 Z = (G4int) (*theElementVector)[iel]->GetZ 505 Z = (G4int) (*theElementVector)[iel]->GetZ();
514 // ---- MGP --- Corrected duplicated cross 506 // ---- MGP --- Corrected duplicated cross section calculation here from
515 // G4hLowEnergyIonisation original 507 // G4hLowEnergyIonisation original
516 G4double microCross = MicroscopicCrossSect 508 G4double microCross = MicroscopicCrossSection( particleDef,
517 lowEdgeEnergy, 509 lowEdgeEnergy,
518 Z, 510 Z,
519 deltaCut ) ; 511 deltaCut ) ;
520 //totalCrossSectionMap [Z] = microCross; 512 //totalCrossSectionMap [Z] = microCross;
521 sigma += theAtomicNumDensityVector[iel] * 513 sigma += theAtomicNumDensityVector[iel] * microCross ;
522 } 514 }
523 515
524 // mean free path = 1./macroscopic cross 516 // mean free path = 1./macroscopic cross section
525 517
526 value = sigma<=0 ? DBL_MAX : 1./sigma ; 518 value = sigma<=0 ? DBL_MAX : 1./sigma ;
527 519
528 aVector->PutValue(i, value) ; 520 aVector->PutValue(i, value) ;
529 } 521 }
530 522
531 theMeanFreePathTable->insert(aVector); 523 theMeanFreePathTable->insert(aVector);
532 } 524 }
533 525
534 } 526 }
535 527
536 528
537 // ------------------------------------------- 529 // --------------------------------------------------------------------
538 G4double G4hImpactIonisation::MicroscopicCross 530 G4double G4hImpactIonisation::MicroscopicCrossSection(const G4ParticleDefinition& particleDef,
539 G4double kineticEnergy, 531 G4double kineticEnergy,
540 G4double atomicNumber, 532 G4double atomicNumber,
541 G4double deltaCutInEnergy) c 533 G4double deltaCutInEnergy) const
542 { 534 {
543 //****************************************** 535 //******************************************************************
544 // cross section formula is OK for spin=0, 536 // cross section formula is OK for spin=0, 1/2, 1 only !
545 // *************************************** 537 // *****************************************************************
546 538
547 // Calculates the microscopic cross sectio 539 // Calculates the microscopic cross section in GEANT4 internal units
548 // Formula documented in Geant4 Phys. Ref. 540 // Formula documented in Geant4 Phys. Ref. Manual
549 // ( it is called for elements, AtomicNumb 541 // ( it is called for elements, AtomicNumber = z )
550 542
551 G4double totalCrossSection = 0.; 543 G4double totalCrossSection = 0.;
552 544
553 // Particle mass and energy 545 // Particle mass and energy
554 G4double particleMass = initialMass; 546 G4double particleMass = initialMass;
555 G4double energy = kineticEnergy + particleMa 547 G4double energy = kineticEnergy + particleMass;
556 548
557 // Some kinematics 549 // Some kinematics
558 G4double gamma = energy / particleMass; 550 G4double gamma = energy / particleMass;
559 G4double beta2 = 1. - 1. / (gamma * gamma); 551 G4double beta2 = 1. - 1. / (gamma * gamma);
560 G4double var = electron_mass_c2 / particleMa 552 G4double var = electron_mass_c2 / particleMass;
561 G4double tMax = 2. * electron_mass_c2 * (g 553 G4double tMax = 2. * electron_mass_c2 * (gamma*gamma - 1.) / (1. + 2.* gamma*var + var*var);
562 554
563 // Calculate the total cross section 555 // Calculate the total cross section
564 556
565 if ( tMax > deltaCutInEnergy ) 557 if ( tMax > deltaCutInEnergy )
566 { 558 {
567 var = deltaCutInEnergy / tMax; 559 var = deltaCutInEnergy / tMax;
568 totalCrossSection = (1. - var * (1. - be 560 totalCrossSection = (1. - var * (1. - beta2 * std::log(var))) / deltaCutInEnergy ;
569 561
570 G4double spin = particleDef.GetPDGSpin() 562 G4double spin = particleDef.GetPDGSpin() ;
571 563
572 // +term for spin=1/2 particle 564 // +term for spin=1/2 particle
573 if (spin == 0.5) 565 if (spin == 0.5)
574 { 566 {
575 totalCrossSection += 0.5 * (tMax - deltaC 567 totalCrossSection += 0.5 * (tMax - deltaCutInEnergy) / (energy*energy);
576 } 568 }
577 // +term for spin=1 particle 569 // +term for spin=1 particle
578 else if (spin > 0.9 ) 570 else if (spin > 0.9 )
579 { 571 {
580 totalCrossSection += -std::log(var) / 572 totalCrossSection += -std::log(var) /
581 (3. * deltaCutInEnergy) + (tMax - deltaC 573 (3. * deltaCutInEnergy) + (tMax - deltaCutInEnergy) * ( (5. + 1. /var)*0.25 / (energy*energy) -
582 beta2 / (tMax * deltaCutIn 574 beta2 / (tMax * deltaCutInEnergy) ) / 3. ;
583 } 575 }
584 totalCrossSection *= twopi_mc2_rcl2 * at 576 totalCrossSection *= twopi_mc2_rcl2 * atomicNumber / beta2 ;
585 } 577 }
586 578
587 //std::cout << "Microscopic = " << totalCros 579 //std::cout << "Microscopic = " << totalCrossSection/barn
588 // << ", e = " << kineticEnergy/MeV <<std 580 // << ", e = " << kineticEnergy/MeV <<std:: endl;
589 581
590 return totalCrossSection ; 582 return totalCrossSection ;
591 } 583 }
592 584
593 585
594 586
595 // ------------------------------------------- 587 // --------------------------------------------------------------------
596 G4double G4hImpactIonisation::GetMeanFreePath( 588 G4double G4hImpactIonisation::GetMeanFreePath(const G4Track& track,
597 G4double, // previousStepSize 589 G4double, // previousStepSize
598 enum G4ForceCondition* conditi 590 enum G4ForceCondition* condition)
599 { 591 {
600 const G4DynamicParticle* dynamicParticle = t 592 const G4DynamicParticle* dynamicParticle = track.GetDynamicParticle();
601 const G4MaterialCutsCouple* couple = track.G 593 const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
602 const G4Material* material = couple->GetMate 594 const G4Material* material = couple->GetMaterial();
603 595
604 G4double meanFreePath = DBL_MAX; 596 G4double meanFreePath = DBL_MAX;
605 // ---- MGP ---- What is the meaning of the 597 // ---- MGP ---- What is the meaning of the local variable isOutOfRange?
606 G4bool isOutRange = false; 598 G4bool isOutRange = false;
607 599
608 *condition = NotForced ; 600 *condition = NotForced ;
609 601
610 G4double kineticEnergy = (dynamicParticle->G 602 G4double kineticEnergy = (dynamicParticle->GetKineticEnergy())*initialMass/(dynamicParticle->GetMass());
611 charge = dynamicParticle->GetCharge()/eplus; 603 charge = dynamicParticle->GetCharge()/eplus;
612 chargeSquare = theIonEffChargeModel->TheValu 604 chargeSquare = theIonEffChargeModel->TheValue(dynamicParticle, material);
613 605
614 if (kineticEnergy < LowestKineticEnergy) 606 if (kineticEnergy < LowestKineticEnergy)
615 { 607 {
616 meanFreePath = DBL_MAX; 608 meanFreePath = DBL_MAX;
617 } 609 }
618 else 610 else
619 { 611 {
620 if (kineticEnergy > HighestKineticEnergy 612 if (kineticEnergy > HighestKineticEnergy) kineticEnergy = HighestKineticEnergy;
621 meanFreePath = (((*theMeanFreePathTable) 613 meanFreePath = (((*theMeanFreePathTable)(couple->GetIndex()))->
622 GetValue(kineticEnergy,isOutRange))/ 614 GetValue(kineticEnergy,isOutRange))/chargeSquare;
623 } 615 }
624 616
625 return meanFreePath ; 617 return meanFreePath ;
626 } 618 }
627 619
628 620
629 // ------------------------------------------- 621 // --------------------------------------------------------------------
630 G4double G4hImpactIonisation::GetConstraints(c 622 G4double G4hImpactIonisation::GetConstraints(const G4DynamicParticle* particle,
631 const G4MaterialCutsCouple* cou 623 const G4MaterialCutsCouple* couple)
632 { 624 {
633 // returns the Step limit 625 // returns the Step limit
634 // dEdx is calculated as well as the range 626 // dEdx is calculated as well as the range
635 // based on Effective Charge Approach 627 // based on Effective Charge Approach
636 628
637 const G4Material* material = couple->GetMate 629 const G4Material* material = couple->GetMaterial();
638 G4Proton* proton = G4Proton::Proton(); 630 G4Proton* proton = G4Proton::Proton();
639 G4AntiProton* antiproton = G4AntiProton::Ant 631 G4AntiProton* antiproton = G4AntiProton::AntiProton();
640 632
641 G4double stepLimit = 0.; 633 G4double stepLimit = 0.;
642 G4double dx, highEnergy; 634 G4double dx, highEnergy;
643 635
644 G4double massRatio = proton_mass_c2/(particl 636 G4double massRatio = proton_mass_c2/(particle->GetMass()) ;
645 G4double kineticEnergy = particle->GetKineti 637 G4double kineticEnergy = particle->GetKineticEnergy() ;
646 638
647 // Scale the kinetic energy 639 // Scale the kinetic energy
648 640
649 G4double tScaled = kineticEnergy*massRatio ; 641 G4double tScaled = kineticEnergy*massRatio ;
650 fBarkas = 0.; 642 fBarkas = 0.;
651 643
652 if (charge > 0.) 644 if (charge > 0.)
653 { 645 {
654 highEnergy = protonHighEnergy ; 646 highEnergy = protonHighEnergy ;
655 fRangeNow = G4EnergyLossTables::GetRange 647 fRangeNow = G4EnergyLossTables::GetRange(proton, tScaled, couple);
656 dx = G4EnergyLossTables::GetRange(proton 648 dx = G4EnergyLossTables::GetRange(proton, highEnergy, couple);
657 fdEdx = G4EnergyLossTables::GetDEDX(prot 649 fdEdx = G4EnergyLossTables::GetDEDX(proton, tScaled, couple)
658 * chargeSquare ; 650 * chargeSquare ;
659 651
660 // Correction for positive ions 652 // Correction for positive ions
661 if (theBarkas && tScaled > highEnergy) 653 if (theBarkas && tScaled > highEnergy)
662 { 654 {
663 fBarkas = BarkasTerm(material,tScaled)*std 655 fBarkas = BarkasTerm(material,tScaled)*std::sqrt(chargeSquare)*chargeSquare
664 + BlochTerm(material,tScaled,chargeSquar 656 + BlochTerm(material,tScaled,chargeSquare);
665 } 657 }
666 // Antiprotons and negative hadrons 658 // Antiprotons and negative hadrons
667 } 659 }
668 else 660 else
669 { 661 {
670 highEnergy = antiprotonHighEnergy ; 662 highEnergy = antiprotonHighEnergy ;
671 fRangeNow = G4EnergyLossTables::GetRange 663 fRangeNow = G4EnergyLossTables::GetRange(antiproton, tScaled, couple);
672 dx = G4EnergyLossTables::GetRange(antipr 664 dx = G4EnergyLossTables::GetRange(antiproton, highEnergy, couple);
673 fdEdx = G4EnergyLossTables::GetDEDX(anti 665 fdEdx = G4EnergyLossTables::GetDEDX(antiproton, tScaled, couple) * chargeSquare ;
674 666
675 if (theBarkas && tScaled > highEnergy) 667 if (theBarkas && tScaled > highEnergy)
676 { 668 {
677 fBarkas = -BarkasTerm(material,tScaled)*st 669 fBarkas = -BarkasTerm(material,tScaled)*std::sqrt(chargeSquare)*chargeSquare
678 + BlochTerm(material,tScaled,chargeSquar 670 + BlochTerm(material,tScaled,chargeSquare);
679 } 671 }
680 } 672 }
681 /* 673 /*
682 const G4Material* mat = couple->GetMateria 674 const G4Material* mat = couple->GetMaterial();
683 G4double fac = gram/(MeV*cm2*mat->GetDensi 675 G4double fac = gram/(MeV*cm2*mat->GetDensity());
684 G4cout << particle->GetDefinition()->GetPa 676 G4cout << particle->GetDefinition()->GetParticleName()
685 << " in " << mat->GetName() 677 << " in " << mat->GetName()
686 << " E(MeV)= " << kineticEnergy/MeV 678 << " E(MeV)= " << kineticEnergy/MeV
687 << " dedx(MeV*cm^2/g)= " << fdEdx*fac 679 << " dedx(MeV*cm^2/g)= " << fdEdx*fac
688 << " barcas(MeV*cm^2/gram)= " << fBarkas*f 680 << " barcas(MeV*cm^2/gram)= " << fBarkas*fac
689 << " Q^2= " << chargeSquare 681 << " Q^2= " << chargeSquare
690 << G4endl; 682 << G4endl;
691 */ 683 */
692 // scaling back 684 // scaling back
693 fRangeNow /= (chargeSquare*massRatio) ; 685 fRangeNow /= (chargeSquare*massRatio) ;
694 dx /= (chargeSquare*massRatio) ; 686 dx /= (chargeSquare*massRatio) ;
695 687
696 stepLimit = fRangeNow ; 688 stepLimit = fRangeNow ;
697 G4double r = std::min(finalRange, couple->Ge 689 G4double r = std::min(finalRange, couple->GetProductionCuts()
698 ->GetProductionCut(idxG4ElectronCut)); 690 ->GetProductionCut(idxG4ElectronCut));
699 691
700 if (fRangeNow > r) 692 if (fRangeNow > r)
701 { 693 {
702 stepLimit = dRoverRange*fRangeNow + r*(1 694 stepLimit = dRoverRange*fRangeNow + r*(1.0 - dRoverRange)*(2.0 - r/fRangeNow);
703 if (stepLimit > fRangeNow) stepLimit = f 695 if (stepLimit > fRangeNow) stepLimit = fRangeNow;
704 } 696 }
705 // compute the (random) Step limit in stan 697 // compute the (random) Step limit in standard energy range
706 if(tScaled > highEnergy ) 698 if(tScaled > highEnergy )
707 { 699 {
708 // add Barkas correction directly to ded 700 // add Barkas correction directly to dedx
709 fdEdx += fBarkas; 701 fdEdx += fBarkas;
710 702
711 if(stepLimit > fRangeNow - dx*0.9) stepL 703 if(stepLimit > fRangeNow - dx*0.9) stepLimit = fRangeNow - dx*0.9 ;
712 704
713 // Step limit in low energy range 705 // Step limit in low energy range
714 } 706 }
715 else 707 else
716 { 708 {
717 G4double x = dx*paramStepLimit; 709 G4double x = dx*paramStepLimit;
718 if (stepLimit > x) stepLimit = x; 710 if (stepLimit > x) stepLimit = x;
719 } 711 }
720 return stepLimit; 712 return stepLimit;
721 } 713 }
722 714
723 715
724 // ------------------------------------------- 716 // --------------------------------------------------------------------
725 G4VParticleChange* G4hImpactIonisation::AlongS 717 G4VParticleChange* G4hImpactIonisation::AlongStepDoIt(const G4Track& track,
726 const G4Step& step) 718 const G4Step& step)
727 { 719 {
728 // compute the energy loss after a step 720 // compute the energy loss after a step
729 G4Proton* proton = G4Proton::Proton(); 721 G4Proton* proton = G4Proton::Proton();
730 G4AntiProton* antiproton = G4AntiProton::Ant 722 G4AntiProton* antiproton = G4AntiProton::AntiProton();
731 G4double finalT = 0.; 723 G4double finalT = 0.;
732 724
733 aParticleChange.Initialize(track) ; 725 aParticleChange.Initialize(track) ;
734 726
735 const G4MaterialCutsCouple* couple = track.G 727 const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
736 const G4Material* material = couple->GetMate 728 const G4Material* material = couple->GetMaterial();
737 729
738 // get the actual (true) Step length from st 730 // get the actual (true) Step length from step
739 const G4double stepLength = step.GetStepLeng 731 const G4double stepLength = step.GetStepLength() ;
740 732
741 const G4DynamicParticle* particle = track.Ge 733 const G4DynamicParticle* particle = track.GetDynamicParticle() ;
742 734
743 G4double kineticEnergy = particle->GetKineti 735 G4double kineticEnergy = particle->GetKineticEnergy() ;
744 G4double massRatio = proton_mass_c2/(particl 736 G4double massRatio = proton_mass_c2/(particle->GetMass()) ;
745 G4double tScaled = kineticEnergy * massRatio 737 G4double tScaled = kineticEnergy * massRatio ;
746 G4double eLoss = 0.0 ; 738 G4double eLoss = 0.0 ;
747 G4double nLoss = 0.0 ; 739 G4double nLoss = 0.0 ;
748 740
749 741
750 // very small particle energy 742 // very small particle energy
751 if (kineticEnergy < MinKineticEnergy) 743 if (kineticEnergy < MinKineticEnergy)
752 { 744 {
753 eLoss = kineticEnergy ; 745 eLoss = kineticEnergy ;
754 // particle energy outside tabulated ene 746 // particle energy outside tabulated energy range
755 } 747 }
756 748
757 else if( kineticEnergy > HighestKineticEnerg 749 else if( kineticEnergy > HighestKineticEnergy)
758 { 750 {
759 eLoss = stepLength * fdEdx ; 751 eLoss = stepLength * fdEdx ;
760 // big step 752 // big step
761 } 753 }
762 else if (stepLength >= fRangeNow ) 754 else if (stepLength >= fRangeNow )
763 { 755 {
764 eLoss = kineticEnergy ; 756 eLoss = kineticEnergy ;
765 757
766 // tabulated range 758 // tabulated range
767 } 759 }
768 else 760 else
769 { 761 {
770 // step longer than linear step limit 762 // step longer than linear step limit
771 if(stepLength > linLossLimit * fRangeNow 763 if(stepLength > linLossLimit * fRangeNow)
772 { 764 {
773 G4double rScaled = fRangeNow * massRatio * 765 G4double rScaled = fRangeNow * massRatio * chargeSquare ;
774 G4double sScaled = stepLength * massRatio 766 G4double sScaled = stepLength * massRatio * chargeSquare ;
775 767
776 if(charge > 0.0) 768 if(charge > 0.0)
777 { 769 {
778 eLoss = G4EnergyLossTables::GetPrecise 770 eLoss = G4EnergyLossTables::GetPreciseEnergyFromRange(proton,rScaled, couple) -
779 G4EnergyLossTables::GetPreciseEnergyFromRa 771 G4EnergyLossTables::GetPreciseEnergyFromRange(proton,rScaled-sScaled,couple) ;
780 772
781 } 773 }
782 else 774 else
783 { 775 {
784 // Antiproton 776 // Antiproton
785 eLoss = G4EnergyLossTables::GetPrecise 777 eLoss = G4EnergyLossTables::GetPreciseEnergyFromRange(antiproton,rScaled,couple) -
786 G4EnergyLossTables::GetPreciseEnergyFromRa 778 G4EnergyLossTables::GetPreciseEnergyFromRange(antiproton,rScaled-sScaled,couple) ;
787 } 779 }
788 eLoss /= massRatio ; 780 eLoss /= massRatio ;
789 781
790 // Barkas correction at big step 782 // Barkas correction at big step
791 eLoss += fBarkas * stepLength; 783 eLoss += fBarkas * stepLength;
792 784
793 // step shorter than linear step limit 785 // step shorter than linear step limit
794 } 786 }
795 else 787 else
796 { 788 {
797 eLoss = stepLength *fdEdx ; 789 eLoss = stepLength *fdEdx ;
798 } 790 }
799 if (nStopping && tScaled < protonHighEne 791 if (nStopping && tScaled < protonHighEnergy)
800 { 792 {
801 nLoss = (theNuclearStoppingModel->TheValue 793 nLoss = (theNuclearStoppingModel->TheValue(particle, material)) * stepLength;
802 } 794 }
803 } 795 }
804 796
805 if (eLoss < 0.0) eLoss = 0.0; 797 if (eLoss < 0.0) eLoss = 0.0;
806 798
807 finalT = kineticEnergy - eLoss - nLoss; 799 finalT = kineticEnergy - eLoss - nLoss;
808 800
809 if ( EnlossFlucFlag && 0.0 < eLoss && finalT 801 if ( EnlossFlucFlag && 0.0 < eLoss && finalT > MinKineticEnergy)
810 { 802 {
811 803
812 // now the electron loss with fluctuati 804 // now the electron loss with fluctuation
813 eLoss = ElectronicLossFluctuation(partic 805 eLoss = ElectronicLossFluctuation(particle, couple, eLoss, stepLength) ;
814 if (eLoss < 0.0) eLoss = 0.0; 806 if (eLoss < 0.0) eLoss = 0.0;
815 finalT = kineticEnergy - eLoss - nLoss; 807 finalT = kineticEnergy - eLoss - nLoss;
816 } 808 }
817 809
818 // stop particle if the kinetic energy <= M 810 // stop particle if the kinetic energy <= MinKineticEnergy
819 if (finalT*massRatio <= MinKineticEnergy ) 811 if (finalT*massRatio <= MinKineticEnergy )
820 { 812 {
821 813
822 finalT = 0.0; 814 finalT = 0.0;
823 if (!particle->GetDefinition()->GetProce 815 if (!particle->GetDefinition()->GetProcessManager()->GetAtRestProcessVector()->size())
824 aParticleChange.ProposeTrackStatus(fSt 816 aParticleChange.ProposeTrackStatus(fStopAndKill);
825 else 817 else
826 aParticleChange.ProposeTrackStatus(fSt 818 aParticleChange.ProposeTrackStatus(fStopButAlive);
827 } 819 }
828 820
829 aParticleChange.ProposeEnergy( finalT ); 821 aParticleChange.ProposeEnergy( finalT );
830 eLoss = kineticEnergy-finalT; 822 eLoss = kineticEnergy-finalT;
831 823
832 aParticleChange.ProposeLocalEnergyDeposit(eL 824 aParticleChange.ProposeLocalEnergyDeposit(eLoss);
833 return &aParticleChange ; 825 return &aParticleChange ;
834 } 826 }
835 827
836 828
837 829
838 // ------------------------------------------- 830 // --------------------------------------------------------------------
839 G4double G4hImpactIonisation::ProtonParametris 831 G4double G4hImpactIonisation::ProtonParametrisedDEDX(const G4MaterialCutsCouple* couple,
840 G4double kineticEnergy) const 832 G4double kineticEnergy) const
841 { 833 {
842 const G4Material* material = couple->GetMate 834 const G4Material* material = couple->GetMaterial();
843 G4Proton* proton = G4Proton::Proton(); 835 G4Proton* proton = G4Proton::Proton();
844 G4double eLoss = 0.; 836 G4double eLoss = 0.;
845 837
846 // Free Electron Gas Model 838 // Free Electron Gas Model
847 if(kineticEnergy < protonLowEnergy) { 839 if(kineticEnergy < protonLowEnergy) {
848 eLoss = (protonModel->TheValue(proton, mat 840 eLoss = (protonModel->TheValue(proton, material, protonLowEnergy))
849 * std::sqrt(kineticEnergy/protonLowEnerg 841 * std::sqrt(kineticEnergy/protonLowEnergy) ;
850 842
851 // Parametrisation 843 // Parametrisation
852 } else { 844 } else {
853 eLoss = protonModel->TheValue(proton, mate 845 eLoss = protonModel->TheValue(proton, material, kineticEnergy) ;
854 } 846 }
855 847
856 // Delta rays energy 848 // Delta rays energy
857 eLoss -= DeltaRaysEnergy(couple,kineticEnerg 849 eLoss -= DeltaRaysEnergy(couple,kineticEnergy,proton_mass_c2) ;
858 850
859 if(verboseLevel > 2) { 851 if(verboseLevel > 2) {
860 G4cout << "p E(MeV)= " << kineticEnergy/Me 852 G4cout << "p E(MeV)= " << kineticEnergy/MeV
861 << " dE/dx(MeV/mm)= " << eLoss*mm/M 853 << " dE/dx(MeV/mm)= " << eLoss*mm/MeV
862 << " for " << material->GetName() 854 << " for " << material->GetName()
863 << " model: " << protonModel << G4e 855 << " model: " << protonModel << G4endl;
864 } 856 }
865 857
866 if(eLoss < 0.0) eLoss = 0.0 ; 858 if(eLoss < 0.0) eLoss = 0.0 ;
867 859
868 return eLoss ; 860 return eLoss ;
869 } 861 }
870 862
871 863
872 864
873 // ------------------------------------------- 865 // --------------------------------------------------------------------
874 G4double G4hImpactIonisation::AntiProtonParame 866 G4double G4hImpactIonisation::AntiProtonParametrisedDEDX(const G4MaterialCutsCouple* couple,
875 G4double kineticEnergy) const 867 G4double kineticEnergy) const
876 { 868 {
877 const G4Material* material = couple->GetMate 869 const G4Material* material = couple->GetMaterial();
878 G4AntiProton* antiproton = G4AntiProton::Ant 870 G4AntiProton* antiproton = G4AntiProton::AntiProton();
879 G4double eLoss = 0.0 ; 871 G4double eLoss = 0.0 ;
880 872
881 // Antiproton model is used 873 // Antiproton model is used
882 if(antiprotonModel->IsInCharge(antiproton,ma 874 if(antiprotonModel->IsInCharge(antiproton,material)) {
883 if(kineticEnergy < antiprotonLowEnergy) { 875 if(kineticEnergy < antiprotonLowEnergy) {
884 eLoss = antiprotonModel->TheValue(antipr 876 eLoss = antiprotonModel->TheValue(antiproton,material,antiprotonLowEnergy)
885 * std::sqrt(kineticEnergy/antiprotonLowEnerg 877 * std::sqrt(kineticEnergy/antiprotonLowEnergy) ;
886 878
887 // Parametrisation 879 // Parametrisation
888 } else { 880 } else {
889 eLoss = antiprotonModel->TheValue(antipr 881 eLoss = antiprotonModel->TheValue(antiproton,material,
890 kineticEnergy); 882 kineticEnergy);
891 } 883 }
892 884
893 // The proton model is used + Barkas corre 885 // The proton model is used + Barkas correction
894 } else { 886 } else {
895 if(kineticEnergy < protonLowEnergy) { 887 if(kineticEnergy < protonLowEnergy) {
896 eLoss = protonModel->TheValue(G4Proton:: 888 eLoss = protonModel->TheValue(G4Proton::Proton(),material,protonLowEnergy)
897 * std::sqrt(kineticEnergy/protonLowEnergy) ; 889 * std::sqrt(kineticEnergy/protonLowEnergy) ;
898 890
899 // Parametrisation 891 // Parametrisation
900 } else { 892 } else {
901 eLoss = protonModel->TheValue(G4Proton:: 893 eLoss = protonModel->TheValue(G4Proton::Proton(),material,
902 kineticEnergy); 894 kineticEnergy);
903 } 895 }
904 //if(theBarkas) eLoss -= 2.0*BarkasTerm(ma 896 //if(theBarkas) eLoss -= 2.0*BarkasTerm(material, kineticEnergy);
905 } 897 }
906 898
907 // Delta rays energy 899 // Delta rays energy
908 eLoss -= DeltaRaysEnergy(couple,kineticEnerg 900 eLoss -= DeltaRaysEnergy(couple,kineticEnergy,proton_mass_c2) ;
909 901
910 if(verboseLevel > 2) { 902 if(verboseLevel > 2) {
911 G4cout << "pbar E(MeV)= " << kineticEnergy 903 G4cout << "pbar E(MeV)= " << kineticEnergy/MeV
912 << " dE/dx(MeV/mm)= " << eLoss*mm/M 904 << " dE/dx(MeV/mm)= " << eLoss*mm/MeV
913 << " for " << material->GetName() 905 << " for " << material->GetName()
914 << " model: " << protonModel << G4e 906 << " model: " << protonModel << G4endl;
915 } 907 }
916 908
917 if(eLoss < 0.0) eLoss = 0.0 ; 909 if(eLoss < 0.0) eLoss = 0.0 ;
918 910
919 return eLoss ; 911 return eLoss ;
920 } 912 }
921 913
922 914
923 // ------------------------------------------- 915 // --------------------------------------------------------------------
924 G4double G4hImpactIonisation::DeltaRaysEnergy( 916 G4double G4hImpactIonisation::DeltaRaysEnergy(const G4MaterialCutsCouple* couple,
925 G4double kineticEnergy, 917 G4double kineticEnergy,
926 G4double particleMass) const 918 G4double particleMass) const
927 { 919 {
928 G4double dLoss = 0.; 920 G4double dLoss = 0.;
929 921
930 G4double deltaCutNow = cutForDelta[(couple-> 922 G4double deltaCutNow = cutForDelta[(couple->GetIndex())] ;
931 const G4Material* material = couple->GetMate 923 const G4Material* material = couple->GetMaterial();
932 G4double electronDensity = material->GetElec 924 G4double electronDensity = material->GetElectronDensity();
933 G4double excitationEnergy = material->GetIon 925 G4double excitationEnergy = material->GetIonisation()->GetMeanExcitationEnergy();
934 926
935 G4double tau = kineticEnergy / particleMass 927 G4double tau = kineticEnergy / particleMass ;
936 G4double rateMass = electron_mass_c2/particl 928 G4double rateMass = electron_mass_c2/particleMass ;
937 929
938 // some local variables 930 // some local variables
939 931
940 G4double gamma = tau + 1.0 ; 932 G4double gamma = tau + 1.0 ;
941 G4double bg2 = tau*(tau+2.0) ; 933 G4double bg2 = tau*(tau+2.0) ;
942 G4double beta2 = bg2/(gamma*gamma) ; 934 G4double beta2 = bg2/(gamma*gamma) ;
943 G4double tMax = 2.*electron_mass_c2*bg2/(1.0 935 G4double tMax = 2.*electron_mass_c2*bg2/(1.0+2.0*gamma*rateMass+rateMass*rateMass) ;
944 936
945 // Validity range for delta electron cross s 937 // Validity range for delta electron cross section
946 G4double deltaCut = std::max(deltaCutNow, ex 938 G4double deltaCut = std::max(deltaCutNow, excitationEnergy);
947 939
948 if ( deltaCut < tMax) 940 if ( deltaCut < tMax)
949 { 941 {
950 G4double x = deltaCut / tMax ; 942 G4double x = deltaCut / tMax ;
951 dLoss = ( beta2 * (x-1.) - std::log(x) ) 943 dLoss = ( beta2 * (x-1.) - std::log(x) ) * twopi_mc2_rcl2 * electronDensity / beta2 ;
952 } 944 }
953 return dLoss ; 945 return dLoss ;
954 } 946 }
955 947
956 948
957 // ------------------------------------------- 949 // -------------------------------------------------------------------------
958 950
959 G4VParticleChange* G4hImpactIonisation::PostSt 951 G4VParticleChange* G4hImpactIonisation::PostStepDoIt(const G4Track& track,
960 const G4Step& step) 952 const G4Step& step)
961 { 953 {
962 // Units are expressed in GEANT4 internal un 954 // Units are expressed in GEANT4 internal units.
963 955
964 // std::cout << "----- Calling PostStepDoIt 956 // std::cout << "----- Calling PostStepDoIt ----- " << std::endl;
965 957
966 aParticleChange.Initialize(track) ; 958 aParticleChange.Initialize(track) ;
967 const G4MaterialCutsCouple* couple = track.G 959 const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
968 const G4DynamicParticle* aParticle = track.G 960 const G4DynamicParticle* aParticle = track.GetDynamicParticle() ;
969 961
970 // Some kinematics 962 // Some kinematics
971 963
972 G4ParticleDefinition* definition = track.Get 964 G4ParticleDefinition* definition = track.GetDefinition();
973 G4double mass = definition->GetPDGMass(); 965 G4double mass = definition->GetPDGMass();
974 G4double kineticEnergy = aParticle->GetKinet 966 G4double kineticEnergy = aParticle->GetKineticEnergy();
975 G4double totalEnergy = kineticEnergy + mass 967 G4double totalEnergy = kineticEnergy + mass ;
976 G4double pSquare = kineticEnergy *( totalEne 968 G4double pSquare = kineticEnergy *( totalEnergy + mass) ;
977 G4double eSquare = totalEnergy * totalEnergy 969 G4double eSquare = totalEnergy * totalEnergy;
978 G4double betaSquare = pSquare / eSquare; 970 G4double betaSquare = pSquare / eSquare;
979 G4ThreeVector particleDirection = aParticle- 971 G4ThreeVector particleDirection = aParticle->GetMomentumDirection() ;
980 972
981 G4double gamma = kineticEnergy / mass + 1.; 973 G4double gamma = kineticEnergy / mass + 1.;
982 G4double r = electron_mass_c2 / mass; 974 G4double r = electron_mass_c2 / mass;
983 G4double tMax = 2. * electron_mass_c2 *(gamm 975 G4double tMax = 2. * electron_mass_c2 *(gamma*gamma - 1.) / (1. + 2.*gamma*r + r*r);
984 976
985 // Validity range for delta electron cross s 977 // Validity range for delta electron cross section
986 G4double deltaCut = cutForDelta[couple->GetI 978 G4double deltaCut = cutForDelta[couple->GetIndex()];
987 979
988 // This should not be a case 980 // This should not be a case
989 if (deltaCut >= tMax) 981 if (deltaCut >= tMax)
990 return G4VContinuousDiscreteProcess::PostS 982 return G4VContinuousDiscreteProcess::PostStepDoIt(track,step);
991 983
992 G4double xc = deltaCut / tMax; 984 G4double xc = deltaCut / tMax;
993 G4double rate = tMax / totalEnergy; 985 G4double rate = tMax / totalEnergy;
994 rate = rate*rate ; 986 rate = rate*rate ;
995 G4double spin = aParticle->GetDefinition()-> 987 G4double spin = aParticle->GetDefinition()->GetPDGSpin() ;
996 988
997 // Sampling follows ... 989 // Sampling follows ...
998 G4double x = 0.; 990 G4double x = 0.;
999 G4double gRej = 0.; 991 G4double gRej = 0.;
1000 992
1001 do { 993 do {
1002 x = xc / (1. - (1. - xc) * G4UniformRand( 994 x = xc / (1. - (1. - xc) * G4UniformRand());
1003 995
1004 if (0.0 == spin) 996 if (0.0 == spin)
1005 { 997 {
1006 gRej = 1.0 - betaSquare * x ; 998 gRej = 1.0 - betaSquare * x ;
1007 } 999 }
1008 else if (0.5 == spin) 1000 else if (0.5 == spin)
1009 { 1001 {
1010 gRej = (1. - betaSquare * x + 0.5 * x*x * r 1002 gRej = (1. - betaSquare * x + 0.5 * x*x * rate) / (1. + 0.5 * rate) ;
1011 } 1003 }
1012 else 1004 else
1013 { 1005 {
1014 gRej = (1. - betaSquare * x ) * (1. + x/(3. 1006 gRej = (1. - betaSquare * x ) * (1. + x/(3.*xc)) +
1015 x*x * rate * (1. + 0.5*x/xc) / 3.0 / 1007 x*x * rate * (1. + 0.5*x/xc) / 3.0 /
1016 (1. + 1./(3.*xc) + rate *(1.+ 0.5/xc) / 3 1008 (1. + 1./(3.*xc) + rate *(1.+ 0.5/xc) / 3.);
1017 } 1009 }
1018 1010
1019 } while( G4UniformRand() > gRej ); 1011 } while( G4UniformRand() > gRej );
1020 1012
1021 G4double deltaKineticEnergy = x * tMax; 1013 G4double deltaKineticEnergy = x * tMax;
1022 G4double deltaTotalMomentum = std::sqrt(del 1014 G4double deltaTotalMomentum = std::sqrt(deltaKineticEnergy *
1023 (deltaKineticEnergy + 2. * electr 1015 (deltaKineticEnergy + 2. * electron_mass_c2 ));
1024 G4double totalMomentum = std::sqrt(pSquare) 1016 G4double totalMomentum = std::sqrt(pSquare) ;
1025 G4double cosTheta = deltaKineticEnergy * (t 1017 G4double cosTheta = deltaKineticEnergy * (totalEnergy + electron_mass_c2) / (deltaTotalMomentum*totalMomentum);
1026 1018
1027 // protection against cosTheta > 1 or < -1 1019 // protection against cosTheta > 1 or < -1
1028 if ( cosTheta < -1. ) cosTheta = -1.; 1020 if ( cosTheta < -1. ) cosTheta = -1.;
1029 if ( cosTheta > 1. ) cosTheta = 1.; 1021 if ( cosTheta > 1. ) cosTheta = 1.;
1030 1022
1031 // direction of the delta electron 1023 // direction of the delta electron
1032 G4double phi = twopi * G4UniformRand(); 1024 G4double phi = twopi * G4UniformRand();
1033 G4double sinTheta = std::sqrt(1. - cosTheta 1025 G4double sinTheta = std::sqrt(1. - cosTheta*cosTheta);
1034 G4double dirX = sinTheta * std::cos(phi); 1026 G4double dirX = sinTheta * std::cos(phi);
1035 G4double dirY = sinTheta * std::sin(phi); 1027 G4double dirY = sinTheta * std::sin(phi);
1036 G4double dirZ = cosTheta; 1028 G4double dirZ = cosTheta;
1037 1029
1038 G4ThreeVector deltaDirection(dirX,dirY,dirZ 1030 G4ThreeVector deltaDirection(dirX,dirY,dirZ);
1039 deltaDirection.rotateUz(particleDirection); 1031 deltaDirection.rotateUz(particleDirection);
1040 1032
1041 // create G4DynamicParticle object for delt 1033 // create G4DynamicParticle object for delta ray
1042 G4DynamicParticle* deltaRay = new G4Dynamic 1034 G4DynamicParticle* deltaRay = new G4DynamicParticle;
1043 deltaRay->SetKineticEnergy( deltaKineticEne 1035 deltaRay->SetKineticEnergy( deltaKineticEnergy );
1044 deltaRay->SetMomentumDirection(deltaDirecti 1036 deltaRay->SetMomentumDirection(deltaDirection.x(),
1045 deltaDirection.y(), 1037 deltaDirection.y(),
1046 deltaDirection.z()); 1038 deltaDirection.z());
1047 deltaRay->SetDefinition(G4Electron::Electro 1039 deltaRay->SetDefinition(G4Electron::Electron());
1048 1040
1049 // fill aParticleChange 1041 // fill aParticleChange
1050 G4double finalKineticEnergy = kineticEnergy 1042 G4double finalKineticEnergy = kineticEnergy - deltaKineticEnergy;
1051 std::size_t totalNumber = 1; << 1043 size_t totalNumber = 1;
1052 1044
1053 // Atomic relaxation 1045 // Atomic relaxation
1054 1046
1055 // ---- MGP ---- Temporary limitation: curr 1047 // ---- MGP ---- Temporary limitation: currently PIXE only for incident protons
1056 1048
1057 std::size_t nSecondaries = 0; << 1049 size_t nSecondaries = 0;
1058 std::vector<G4DynamicParticle*>* secondaryV 1050 std::vector<G4DynamicParticle*>* secondaryVector = 0;
1059 1051
1060 if (definition == G4Proton::ProtonDefinitio 1052 if (definition == G4Proton::ProtonDefinition())
1061 { 1053 {
1062 const G4Material* material = couple->Ge 1054 const G4Material* material = couple->GetMaterial();
1063 1055
1064 // Lazy initialization of pixeCrossSect 1056 // Lazy initialization of pixeCrossSectionHandler
1065 if (pixeCrossSectionHandler == 0) 1057 if (pixeCrossSectionHandler == 0)
1066 { 1058 {
1067 // Instantiate pixeCrossSectionHandler wi 1059 // Instantiate pixeCrossSectionHandler with selected shell cross section models
1068 // Ownership of interpolation is transfer 1060 // Ownership of interpolation is transferred to pixeCrossSectionHandler
1069 G4IInterpolator* interpolation = new G4Lo 1061 G4IInterpolator* interpolation = new G4LogLogInterpolator();
1070 pixeCrossSectionHandler = 1062 pixeCrossSectionHandler =
1071 new G4PixeCrossSectionHandler(interpola 1063 new G4PixeCrossSectionHandler(interpolation,modelK,modelL,modelM,eMinPixe,eMaxPixe);
1072 G4String fileName("proton"); 1064 G4String fileName("proton");
1073 pixeCrossSectionHandler->LoadShellData(fi 1065 pixeCrossSectionHandler->LoadShellData(fileName);
1074 // pixeCrossSectionHandler->PrintData( 1066 // pixeCrossSectionHandler->PrintData();
1075 } 1067 }
1076 1068
1077 // Select an atom in the current materi 1069 // Select an atom in the current material based on the total shell cross sections
1078 G4int Z = pixeCrossSectionHandler->Sele 1070 G4int Z = pixeCrossSectionHandler->SelectRandomAtom(material,kineticEnergy);
1079 // std::cout << "G4hImpactIonisati 1071 // std::cout << "G4hImpactIonisation::PostStepDoIt - Z = " << Z << std::endl;
1080 1072
1081 // G4double microscopicCross = Mic 1073 // G4double microscopicCross = MicroscopicCrossSection(*definition,
1082 // kineticEnergy 1074 // kineticEnergy,
1083 // Z, deltaCut); 1075 // Z, deltaCut);
1084 // G4double crossFromShells = pixeCrossS 1076 // G4double crossFromShells = pixeCrossSectionHandler->FindValue(Z,kineticEnergy);
1085 1077
1086 //std::cout << "G4hImpactIonisation: Z= 1078 //std::cout << "G4hImpactIonisation: Z= "
1087 // << Z 1079 // << Z
1088 // << ", energy = " 1080 // << ", energy = "
1089 // << kineticEnergy/MeV 1081 // << kineticEnergy/MeV
1090 // <<" MeV, microscopic = " 1082 // <<" MeV, microscopic = "
1091 // << microscopicCross/barn 1083 // << microscopicCross/barn
1092 // << " barn, from shells = " 1084 // << " barn, from shells = "
1093 // << crossFromShells/barn 1085 // << crossFromShells/barn
1094 // << " barn" 1086 // << " barn"
1095 // << std::endl; 1087 // << std::endl;
1096 1088
1097 // Select a shell in the target atom ba 1089 // Select a shell in the target atom based on the individual shell cross sections
1098 G4int shellIndex = pixeCrossSectionHand 1090 G4int shellIndex = pixeCrossSectionHandler->SelectRandomShell(Z,kineticEnergy);
1099 1091
1100 G4AtomicTransitionManager* transitionMa 1092 G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
1101 const G4AtomicShell* atomicShell = tran 1093 const G4AtomicShell* atomicShell = transitionManager->Shell(Z,shellIndex);
1102 G4double bindingEnergy = atomicShell->B 1094 G4double bindingEnergy = atomicShell->BindingEnergy();
1103 1095
1104 // if (verboseLevel > 1) 1096 // if (verboseLevel > 1)
1105 // { 1097 // {
1106 // G4cout << "G4hImpactIonisation::Po 1098 // G4cout << "G4hImpactIonisation::PostStepDoIt - Z = "
1107 // << Z 1099 // << Z
1108 // << ", shell = " 1100 // << ", shell = "
1109 // << shellIndex 1101 // << shellIndex
1110 // << ", bindingE (keV) = " 1102 // << ", bindingE (keV) = "
1111 // << bindingEnergy/keV 1103 // << bindingEnergy/keV
1112 // << G4endl; 1104 // << G4endl;
1113 // } 1105 // }
1114 1106
1115 // Generate PIXE if binding energy larg 1107 // Generate PIXE if binding energy larger than cut for photons or electrons
1116 1108
1117 G4ParticleDefinition* type = 0; 1109 G4ParticleDefinition* type = 0;
1118 1110
1119 if (finalKineticEnergy >= bindingEnergy 1111 if (finalKineticEnergy >= bindingEnergy)
1120 // && (bindingEnergy >= minGammaEnergy | 1112 // && (bindingEnergy >= minGammaEnergy || bindingEnergy >= minElectronEnergy) )
1121 { 1113 {
1122 // Vacancy in subshell shellIndex; shellI 1114 // Vacancy in subshell shellIndex; shellId is the subshell identifier in EADL jargon
1123 G4int shellId = atomicShell->ShellId(); 1115 G4int shellId = atomicShell->ShellId();
1124 // Atomic relaxation: generate secondarie 1116 // Atomic relaxation: generate secondaries
1125 secondaryVector = atomicDeexcitation.Gene 1117 secondaryVector = atomicDeexcitation.GenerateParticles(Z, shellId);
1126 1118
1127 // ---- Debug ---- 1119 // ---- Debug ----
1128 //std::cout << "ShellId = " 1120 //std::cout << "ShellId = "
1129 // <<shellId << " ---- Atomic relaxati 1121 // <<shellId << " ---- Atomic relaxation secondaries: ---- "
1130 // << secondaryVector->size() 1122 // << secondaryVector->size()
1131 // << std::endl; 1123 // << std::endl;
1132 1124
1133 // ---- End debug --- 1125 // ---- End debug ---
1134 1126
1135 if (secondaryVector != 0) 1127 if (secondaryVector != 0)
1136 { 1128 {
1137 nSecondaries = secondaryVector->size( 1129 nSecondaries = secondaryVector->size();
1138 for (std::size_t i = 0; i<nSecondarie << 1130 for (size_t i = 0; i<nSecondaries; i++)
1139 { 1131 {
1140 G4DynamicParticle* aSecondary = (*secon 1132 G4DynamicParticle* aSecondary = (*secondaryVector)[i];
1141 if (aSecondary) 1133 if (aSecondary)
1142 { 1134 {
1143 G4double e = aSecondary->GetKinetic 1135 G4double e = aSecondary->GetKineticEnergy();
1144 type = aSecondary->GetDefinition(); 1136 type = aSecondary->GetDefinition();
1145 1137
1146 // ---- Debug ---- 1138 // ---- Debug ----
1147 //if (type == G4Gamma::GammaDefinit 1139 //if (type == G4Gamma::GammaDefinition())
1148 // { 1140 // {
1149 // std::cout << "Z = " << Z 1141 // std::cout << "Z = " << Z
1150 // << ", shell: " << shellId 1142 // << ", shell: " << shellId
1151 // << ", PIXE photon energy 1143 // << ", PIXE photon energy (keV) = " << e/keV
1152 // << std::endl; 1144 // << std::endl;
1153 // } 1145 // }
1154 // ---- End debug --- 1146 // ---- End debug ---
1155 1147
1156 if (e < finalKineticEnergy && 1148 if (e < finalKineticEnergy &&
1157 ((type == G4Gamma::Gamma() && e > min 1149 ((type == G4Gamma::Gamma() && e > minGammaEnergy ) ||
1158 (type == G4Electron::Electron() && e 1150 (type == G4Electron::Electron() && e > minElectronEnergy )))
1159 { 1151 {
1160 // Subtract the energy of the emitted 1152 // Subtract the energy of the emitted secondary from the primary
1161 finalKineticEnergy -= e; 1153 finalKineticEnergy -= e;
1162 totalNumber++; 1154 totalNumber++;
1163 // ---- Debug ---- 1155 // ---- Debug ----
1164 //if (type == G4Gamma::GammaDefinitio 1156 //if (type == G4Gamma::GammaDefinition())
1165 // { 1157 // {
1166 // std::cout << "Z = " << Z 1158 // std::cout << "Z = " << Z
1167 // << ", shell: " << shellId 1159 // << ", shell: " << shellId
1168 // << ", PIXE photon energy (k 1160 // << ", PIXE photon energy (keV) = " << e/keV
1169 // << std::endl; 1161 // << std::endl;
1170 // } 1162 // }
1171 // ---- End debug --- 1163 // ---- End debug ---
1172 } 1164 }
1173 else 1165 else
1174 { 1166 {
1175 // The atomic relaxation product has 1167 // The atomic relaxation product has energy below the cut
1176 // ---- Debug ---- 1168 // ---- Debug ----
1177 // if (type == G4Gamma::GammaDefiniti 1169 // if (type == G4Gamma::GammaDefinition())
1178 // { 1170 // {
1179 // std::cout << "Z = " << Z 1171 // std::cout << "Z = " << Z
1180 // 1172 //
1181 // << ", PIXE photon energy = 1173 // << ", PIXE photon energy = " << e/keV
1182 // << " keV below threshold 1174 // << " keV below threshold " << minGammaEnergy/keV << " keV"
1183 // << std::endl; 1175 // << std::endl;
1184 // } 1176 // }
1185 // ---- End debug --- 1177 // ---- End debug ---
1186 1178
1187 delete aSecondary; 1179 delete aSecondary;
1188 (*secondaryVector)[i] = 0; 1180 (*secondaryVector)[i] = 0;
1189 } 1181 }
1190 } 1182 }
1191 } 1183 }
1192 } 1184 }
1193 } 1185 }
1194 } 1186 }
1195 1187
1196 1188
1197 // Save delta-electrons 1189 // Save delta-electrons
1198 1190
1199 G4double eDeposit = 0.; 1191 G4double eDeposit = 0.;
1200 1192
1201 if (finalKineticEnergy > MinKineticEnergy) 1193 if (finalKineticEnergy > MinKineticEnergy)
1202 { 1194 {
1203 G4double finalPx = totalMomentum*partic 1195 G4double finalPx = totalMomentum*particleDirection.x() - deltaTotalMomentum*deltaDirection.x();
1204 G4double finalPy = totalMomentum*partic 1196 G4double finalPy = totalMomentum*particleDirection.y() - deltaTotalMomentum*deltaDirection.y();
1205 G4double finalPz = totalMomentum*partic 1197 G4double finalPz = totalMomentum*particleDirection.z() - deltaTotalMomentum*deltaDirection.z();
1206 G4double finalMomentum = std::sqrt(fina 1198 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz) ;
1207 finalPx /= finalMomentum; 1199 finalPx /= finalMomentum;
1208 finalPy /= finalMomentum; 1200 finalPy /= finalMomentum;
1209 finalPz /= finalMomentum; 1201 finalPz /= finalMomentum;
1210 1202
1211 aParticleChange.ProposeMomentumDirectio 1203 aParticleChange.ProposeMomentumDirection( finalPx,finalPy,finalPz );
1212 } 1204 }
1213 else 1205 else
1214 { 1206 {
1215 eDeposit = finalKineticEnergy; 1207 eDeposit = finalKineticEnergy;
1216 finalKineticEnergy = 0.; 1208 finalKineticEnergy = 0.;
1217 aParticleChange.ProposeMomentumDirectio 1209 aParticleChange.ProposeMomentumDirection(particleDirection.x(),
1218 particleDirection.y(), 1210 particleDirection.y(),
1219 particleDirection.z()); 1211 particleDirection.z());
1220 if(!aParticle->GetDefinition()->GetProc 1212 if(!aParticle->GetDefinition()->GetProcessManager()->
1221 GetAtRestProcessVector()->size()) 1213 GetAtRestProcessVector()->size())
1222 aParticleChange.ProposeTrackStatus(fS 1214 aParticleChange.ProposeTrackStatus(fStopAndKill);
1223 else 1215 else
1224 aParticleChange.ProposeTrackStatus(fS 1216 aParticleChange.ProposeTrackStatus(fStopButAlive);
1225 } 1217 }
1226 1218
1227 aParticleChange.ProposeEnergy(finalKineticE 1219 aParticleChange.ProposeEnergy(finalKineticEnergy);
1228 aParticleChange.ProposeLocalEnergyDeposit ( 1220 aParticleChange.ProposeLocalEnergyDeposit (eDeposit);
1229 aParticleChange.SetNumberOfSecondaries((G4i << 1221 aParticleChange.SetNumberOfSecondaries(totalNumber);
1230 aParticleChange.AddSecondary(deltaRay); 1222 aParticleChange.AddSecondary(deltaRay);
1231 1223
1232 // ---- Debug ---- 1224 // ---- Debug ----
1233 // std::cout << "RDHadronIonisation - fina 1225 // std::cout << "RDHadronIonisation - finalKineticEnergy (MeV) = "
1234 // << finalKineticEnergy/MeV 1226 // << finalKineticEnergy/MeV
1235 // << ", delta KineticEnergy (keV) = " 1227 // << ", delta KineticEnergy (keV) = "
1236 // << deltaKineticEnergy/keV 1228 // << deltaKineticEnergy/keV
1237 // << ", energy deposit (MeV) = " 1229 // << ", energy deposit (MeV) = "
1238 // << eDeposit/MeV 1230 // << eDeposit/MeV
1239 // << std::endl; 1231 // << std::endl;
1240 // ---- End debug --- 1232 // ---- End debug ---
1241 1233
1242 // Save Fluorescence and Auger 1234 // Save Fluorescence and Auger
1243 1235
1244 if (secondaryVector != 0) 1236 if (secondaryVector != 0)
1245 { 1237 {
1246 for (std::size_t l = 0; l < nSecondarie << 1238 for (size_t l = 0; l < nSecondaries; l++)
1247 { 1239 {
1248 G4DynamicParticle* secondary = (*secondar 1240 G4DynamicParticle* secondary = (*secondaryVector)[l];
1249 if (secondary) aParticleChange.AddSeconda 1241 if (secondary) aParticleChange.AddSecondary(secondary);
1250 1242
1251 // ---- Debug ---- 1243 // ---- Debug ----
1252 //if (secondary != 0) 1244 //if (secondary != 0)
1253 // { 1245 // {
1254 // if (secondary->GetDefinition() == G4 1246 // if (secondary->GetDefinition() == G4Gamma::GammaDefinition())
1255 // { 1247 // {
1256 // G4double eX = secondary->GetKinetic 1248 // G4double eX = secondary->GetKineticEnergy();
1257 // std::cout << " PIXE photon of energ 1249 // std::cout << " PIXE photon of energy " << eX/keV
1258 // << " keV added to ParticleChang 1250 // << " keV added to ParticleChange; total number of secondaries is " << totalNumber
1259 // << std::endl; 1251 // << std::endl;
1260 // } 1252 // }
1261 //} 1253 //}
1262 // ---- End debug --- 1254 // ---- End debug ---
1263 1255
1264 } 1256 }
1265 delete secondaryVector; 1257 delete secondaryVector;
1266 } 1258 }
1267 1259
1268 return G4VContinuousDiscreteProcess::PostSt 1260 return G4VContinuousDiscreteProcess::PostStepDoIt(track,step);
1269 } 1261 }
1270 1262
1271 // ------------------------------------------ 1263 // -------------------------------------------------------------------------
1272 1264
1273 G4double G4hImpactIonisation::ComputeDEDX(con 1265 G4double G4hImpactIonisation::ComputeDEDX(const G4ParticleDefinition* aParticle,
1274 const G4MaterialCutsCouple* coupl 1266 const G4MaterialCutsCouple* couple,
1275 G4double kineticEnergy) 1267 G4double kineticEnergy)
1276 { 1268 {
1277 const G4Material* material = couple->GetMat 1269 const G4Material* material = couple->GetMaterial();
1278 G4Proton* proton = G4Proton::Proton(); 1270 G4Proton* proton = G4Proton::Proton();
1279 G4AntiProton* antiproton = G4AntiProton::An 1271 G4AntiProton* antiproton = G4AntiProton::AntiProton();
1280 G4double dedx = 0.; 1272 G4double dedx = 0.;
1281 1273
1282 G4double tScaled = kineticEnergy * proton_m 1274 G4double tScaled = kineticEnergy * proton_mass_c2 / (aParticle->GetPDGMass()) ;
1283 charge = aParticle->GetPDGCharge() ; 1275 charge = aParticle->GetPDGCharge() ;
1284 1276
1285 if( charge > 0.) 1277 if( charge > 0.)
1286 { 1278 {
1287 if (tScaled > protonHighEnergy) 1279 if (tScaled > protonHighEnergy)
1288 { 1280 {
1289 dedx = G4EnergyLossTables::GetDEDX(proton 1281 dedx = G4EnergyLossTables::GetDEDX(proton,tScaled,couple) ;
1290 } 1282 }
1291 else 1283 else
1292 { 1284 {
1293 dedx = ProtonParametrisedDEDX(couple,tSca 1285 dedx = ProtonParametrisedDEDX(couple,tScaled) ;
1294 } 1286 }
1295 } 1287 }
1296 else 1288 else
1297 { 1289 {
1298 if (tScaled > antiprotonHighEnergy) 1290 if (tScaled > antiprotonHighEnergy)
1299 { 1291 {
1300 dedx = G4EnergyLossTables::GetDEDX(antipr 1292 dedx = G4EnergyLossTables::GetDEDX(antiproton,tScaled,couple);
1301 } 1293 }
1302 else 1294 else
1303 { 1295 {
1304 dedx = AntiProtonParametrisedDEDX(couple, 1296 dedx = AntiProtonParametrisedDEDX(couple,tScaled) ;
1305 } 1297 }
1306 } 1298 }
1307 dedx *= theIonEffChargeModel->TheValue(aPar 1299 dedx *= theIonEffChargeModel->TheValue(aParticle, material, kineticEnergy) ;
1308 1300
1309 return dedx ; 1301 return dedx ;
1310 } 1302 }
1311 1303
1312 1304
1313 G4double G4hImpactIonisation::BarkasTerm(cons 1305 G4double G4hImpactIonisation::BarkasTerm(const G4Material* material,
1314 G4double kineticEnergy) const 1306 G4double kineticEnergy) const
1315 //Function to compute the Barkas term for pro 1307 //Function to compute the Barkas term for protons:
1316 // 1308 //
1317 //Ref. Z_1^3 effect in the stopping power of 1309 //Ref. Z_1^3 effect in the stopping power of matter for charged particles
1318 // J.C Ashley and R.H.Ritchie 1310 // J.C Ashley and R.H.Ritchie
1319 // Physical review B Vol.5 No.7 1 April 1 1311 // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
1320 // 1312 //
1321 { 1313 {
1322 static G4ThreadLocal G4double FTable[47][2] << 1314 static double FTable[47][2] = {
1323 { 0.02, 21.5}, 1315 { 0.02, 21.5},
1324 { 0.03, 20.0}, 1316 { 0.03, 20.0},
1325 { 0.04, 18.0}, 1317 { 0.04, 18.0},
1326 { 0.05, 15.6}, 1318 { 0.05, 15.6},
1327 { 0.06, 15.0}, 1319 { 0.06, 15.0},
1328 { 0.07, 14.0}, 1320 { 0.07, 14.0},
1329 { 0.08, 13.5}, 1321 { 0.08, 13.5},
1330 { 0.09, 13.}, 1322 { 0.09, 13.},
1331 { 0.1, 12.2}, 1323 { 0.1, 12.2},
1332 { 0.2, 9.25}, 1324 { 0.2, 9.25},
1333 { 0.3, 7.0}, 1325 { 0.3, 7.0},
1334 { 0.4, 6.0}, 1326 { 0.4, 6.0},
1335 { 0.5, 4.5}, 1327 { 0.5, 4.5},
1336 { 0.6, 3.5}, 1328 { 0.6, 3.5},
1337 { 0.7, 3.0}, 1329 { 0.7, 3.0},
1338 { 0.8, 2.5}, 1330 { 0.8, 2.5},
1339 { 0.9, 2.0}, 1331 { 0.9, 2.0},
1340 { 1.0, 1.7}, 1332 { 1.0, 1.7},
1341 { 1.2, 1.2}, 1333 { 1.2, 1.2},
1342 { 1.3, 1.0}, 1334 { 1.3, 1.0},
1343 { 1.4, 0.86}, 1335 { 1.4, 0.86},
1344 { 1.5, 0.7}, 1336 { 1.5, 0.7},
1345 { 1.6, 0.61}, 1337 { 1.6, 0.61},
1346 { 1.7, 0.52}, 1338 { 1.7, 0.52},
1347 { 1.8, 0.5}, 1339 { 1.8, 0.5},
1348 { 1.9, 0.43}, 1340 { 1.9, 0.43},
1349 { 2.0, 0.42}, 1341 { 2.0, 0.42},
1350 { 2.1, 0.3}, 1342 { 2.1, 0.3},
1351 { 2.4, 0.2}, 1343 { 2.4, 0.2},
1352 { 3.0, 0.13}, 1344 { 3.0, 0.13},
1353 { 3.08, 0.1}, 1345 { 3.08, 0.1},
1354 { 3.1, 0.09}, 1346 { 3.1, 0.09},
1355 { 3.3, 0.08}, 1347 { 3.3, 0.08},
1356 { 3.5, 0.07}, 1348 { 3.5, 0.07},
1357 { 3.8, 0.06}, 1349 { 3.8, 0.06},
1358 { 4.0, 0.051}, 1350 { 4.0, 0.051},
1359 { 4.1, 0.04}, 1351 { 4.1, 0.04},
1360 { 4.8, 0.03}, 1352 { 4.8, 0.03},
1361 { 5.0, 0.024}, 1353 { 5.0, 0.024},
1362 { 5.1, 0.02}, 1354 { 5.1, 0.02},
1363 { 6.0, 0.013}, 1355 { 6.0, 0.013},
1364 { 6.5, 0.01}, 1356 { 6.5, 0.01},
1365 { 7.0, 0.009}, 1357 { 7.0, 0.009},
1366 { 7.1, 0.008}, 1358 { 7.1, 0.008},
1367 { 8.0, 0.006}, 1359 { 8.0, 0.006},
1368 { 9.0, 0.0032}, 1360 { 9.0, 0.0032},
1369 { 10.0, 0.0025} }; 1361 { 10.0, 0.0025} };
1370 1362
1371 // Information on particle and material 1363 // Information on particle and material
1372 G4double kinE = kineticEnergy ; 1364 G4double kinE = kineticEnergy ;
1373 if(0.5*MeV > kinE) kinE = 0.5*MeV ; 1365 if(0.5*MeV > kinE) kinE = 0.5*MeV ;
1374 G4double gamma = 1.0 + kinE / proton_mass_c 1366 G4double gamma = 1.0 + kinE / proton_mass_c2 ;
1375 G4double beta2 = 1.0 - 1.0/(gamma*gamma) ; 1367 G4double beta2 = 1.0 - 1.0/(gamma*gamma) ;
1376 if(0.0 >= beta2) return 0.0; 1368 if(0.0 >= beta2) return 0.0;
1377 1369
1378 G4double BTerm = 0.0; << 1370 G4double BarkasTerm = 0.0;
1379 //G4double AMaterial = 0.0; << 1371 G4double AMaterial = 0.0;
1380 G4double ZMaterial = 0.0; 1372 G4double ZMaterial = 0.0;
1381 const G4ElementVector* theElementVector = m 1373 const G4ElementVector* theElementVector = material->GetElementVector();
1382 G4int numberOfElements = (G4int)material->G << 1374 G4int numberOfElements = material->GetNumberOfElements();
1383 1375
1384 for (G4int i = 0; i<numberOfElements; i++) 1376 for (G4int i = 0; i<numberOfElements; i++) {
1385 1377
1386 //AMaterial = (*theElementVector)[i]->Get << 1378 AMaterial = (*theElementVector)[i]->GetA()*mole/g;
1387 ZMaterial = (*theElementVector)[i]->GetZ( 1379 ZMaterial = (*theElementVector)[i]->GetZ();
1388 1380
1389 G4double X = 137.0 * 137.0 * beta2 / ZMat 1381 G4double X = 137.0 * 137.0 * beta2 / ZMaterial;
1390 1382
1391 // Variables to compute L_1 1383 // Variables to compute L_1
1392 G4double Eta0Chi = 0.8; 1384 G4double Eta0Chi = 0.8;
1393 G4double EtaChi = Eta0Chi * ( 1.0 + 6.02* 1385 G4double EtaChi = Eta0Chi * ( 1.0 + 6.02*std::pow( ZMaterial,-1.19 ) );
1394 G4double W = ( EtaChi * std::pow( ZMateri 1386 G4double W = ( EtaChi * std::pow( ZMaterial,1.0/6.0 ) ) / std::sqrt(X);
1395 G4double FunctionOfW = FTable[46][1]*FTab 1387 G4double FunctionOfW = FTable[46][1]*FTable[46][0]/W ;
1396 1388
1397 for(G4int j=0; j<47; j++) { 1389 for(G4int j=0; j<47; j++) {
1398 1390
1399 if( W < FTable[j][0] ) { 1391 if( W < FTable[j][0] ) {
1400 1392
1401 if(0 == j) { 1393 if(0 == j) {
1402 FunctionOfW = FTable[0][1] ; 1394 FunctionOfW = FTable[0][1] ;
1403 1395
1404 } else { 1396 } else {
1405 FunctionOfW = (FTable[j][1] - FTabl 1397 FunctionOfW = (FTable[j][1] - FTable[j-1][1]) * (W - FTable[j-1][0])
1406 / (FTable[j][0] - FTable[j-1][0]) 1398 / (FTable[j][0] - FTable[j-1][0])
1407 + FTable[j-1][1] ; 1399 + FTable[j-1][1] ;
1408 } 1400 }
1409 1401
1410 break; 1402 break;
1411 } 1403 }
1412 1404
1413 } 1405 }
1414 1406
1415 BTerm += FunctionOfW /( std::sqrt(ZMateri << 1407 BarkasTerm += FunctionOfW /( std::sqrt(ZMaterial * X) * X);
1416 } 1408 }
1417 1409
1418 BTerm *= twopi_mc2_rcl2 * (material->GetEle << 1410 BarkasTerm *= twopi_mc2_rcl2 * (material->GetElectronDensity()) / beta2 ;
1419 1411
1420 return BTerm; << 1412 return BarkasTerm;
1421 } 1413 }
1422 1414
1423 1415
1424 1416
1425 G4double G4hImpactIonisation::BlochTerm(const 1417 G4double G4hImpactIonisation::BlochTerm(const G4Material* material,
1426 G4double kineticEnergy, 1418 G4double kineticEnergy,
1427 G4double cSquare) const 1419 G4double cSquare) const
1428 //Function to compute the Bloch term for prot 1420 //Function to compute the Bloch term for protons:
1429 // 1421 //
1430 //Ref. Z_1^3 effect in the stopping power of 1422 //Ref. Z_1^3 effect in the stopping power of matter for charged particles
1431 // J.C Ashley and R.H.Ritchie 1423 // J.C Ashley and R.H.Ritchie
1432 // Physical review B Vol.5 No.7 1 April 1 1424 // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
1433 // 1425 //
1434 { 1426 {
1435 G4double eLoss = 0.0 ; 1427 G4double eLoss = 0.0 ;
1436 G4double gamma = 1.0 + kineticEnergy / prot 1428 G4double gamma = 1.0 + kineticEnergy / proton_mass_c2 ;
1437 G4double beta2 = 1.0 - 1.0/(gamma*gamma) ; 1429 G4double beta2 = 1.0 - 1.0/(gamma*gamma) ;
1438 G4double y = cSquare / (137.0*137.0*beta2) 1430 G4double y = cSquare / (137.0*137.0*beta2) ;
1439 1431
1440 if(y < 0.05) { 1432 if(y < 0.05) {
1441 eLoss = 1.202 ; 1433 eLoss = 1.202 ;
1442 1434
1443 } else { 1435 } else {
1444 eLoss = 1.0 / (1.0 + y) ; 1436 eLoss = 1.0 / (1.0 + y) ;
1445 G4double de = eLoss ; 1437 G4double de = eLoss ;
1446 1438
1447 for(G4int i=2; de>eLoss*0.01; i++) { 1439 for(G4int i=2; de>eLoss*0.01; i++) {
1448 de = 1.0/( i * (i*i + y)) ; 1440 de = 1.0/( i * (i*i + y)) ;
1449 eLoss += de ; 1441 eLoss += de ;
1450 } 1442 }
1451 } 1443 }
1452 eLoss *= -1.0 * y * cSquare * twopi_mc2_rcl 1444 eLoss *= -1.0 * y * cSquare * twopi_mc2_rcl2 *
1453 (material->GetElectronDensity()) / beta2 1445 (material->GetElectronDensity()) / beta2 ;
1454 1446
1455 return eLoss; 1447 return eLoss;
1456 } 1448 }
1457 1449
1458 1450
1459 1451
1460 G4double G4hImpactIonisation::ElectronicLossF 1452 G4double G4hImpactIonisation::ElectronicLossFluctuation(
1461 const G4DynamicParticle* partic 1453 const G4DynamicParticle* particle,
1462 const G4MaterialCutsCouple* cou 1454 const G4MaterialCutsCouple* couple,
1463 G4double meanLoss, 1455 G4double meanLoss,
1464 G4double step) const 1456 G4double step) const
1465 // calculate actual loss from the mean loss 1457 // calculate actual loss from the mean loss
1466 // The model used to get the fluctuation is 1458 // The model used to get the fluctuation is essentially the same
1467 // as in Glandz in Geant3. 1459 // as in Glandz in Geant3.
1468 { 1460 {
1469 // data members to speed up the fluctuation 1461 // data members to speed up the fluctuation calculation
1470 // G4int imat ; 1462 // G4int imat ;
1471 // G4double f1Fluct,f2Fluct,e1Fluct,e2Fluc 1463 // G4double f1Fluct,f2Fluct,e1Fluct,e2Fluct,rateFluct,ipotFluct;
1472 // G4double e1LogFluct,e2LogFluct,ipotLogF 1464 // G4double e1LogFluct,e2LogFluct,ipotLogFluct;
1473 1465
1474 static const G4double minLoss = 1.*eV ; 1466 static const G4double minLoss = 1.*eV ;
1475 static const G4double kappa = 10. ; 1467 static const G4double kappa = 10. ;
1476 static const G4double theBohrBeta2 = 50.0 * 1468 static const G4double theBohrBeta2 = 50.0 * keV/proton_mass_c2 ;
1477 1469
1478 const G4Material* material = couple->GetMat 1470 const G4Material* material = couple->GetMaterial();
1479 G4int imaterial = couple->GetIndex() ; 1471 G4int imaterial = couple->GetIndex() ;
1480 G4double ipotFluct = material->GetIonisat 1472 G4double ipotFluct = material->GetIonisation()->GetMeanExcitationEnergy() ;
1481 G4double electronDensity = material->GetEle 1473 G4double electronDensity = material->GetElectronDensity() ;
1482 G4double zeff = electronDensity/(material-> 1474 G4double zeff = electronDensity/(material->GetTotNbOfAtomsPerVolume()) ;
1483 1475
1484 // get particle data 1476 // get particle data
1485 G4double tkin = particle->GetKineticEnerg 1477 G4double tkin = particle->GetKineticEnergy();
1486 G4double particleMass = particle->GetMass() 1478 G4double particleMass = particle->GetMass() ;
1487 G4double deltaCutInKineticEnergyNow = cutFo 1479 G4double deltaCutInKineticEnergyNow = cutForDelta[imaterial];
1488 1480
1489 // shortcut for very very small loss 1481 // shortcut for very very small loss
1490 if(meanLoss < minLoss) return meanLoss ; 1482 if(meanLoss < minLoss) return meanLoss ;
1491 1483
1492 // Validity range for delta electron cross 1484 // Validity range for delta electron cross section
1493 G4double threshold = std::max(deltaCutInKin 1485 G4double threshold = std::max(deltaCutInKineticEnergyNow,ipotFluct);
1494 G4double loss, siga; 1486 G4double loss, siga;
1495 1487
1496 G4double rmass = electron_mass_c2/particleM 1488 G4double rmass = electron_mass_c2/particleMass;
1497 G4double tau = tkin/particleMass; 1489 G4double tau = tkin/particleMass;
1498 G4double tau1 = tau+1.0; 1490 G4double tau1 = tau+1.0;
1499 G4double tau2 = tau*(tau+2.); 1491 G4double tau2 = tau*(tau+2.);
1500 G4double tMax = 2.*electron_mass_c2*tau2 1492 G4double tMax = 2.*electron_mass_c2*tau2/(1.+2.*tau1*rmass+rmass*rmass);
1501 1493
1502 1494
1503 if(tMax > threshold) tMax = threshold; 1495 if(tMax > threshold) tMax = threshold;
1504 G4double beta2 = tau2/(tau1*tau1); 1496 G4double beta2 = tau2/(tau1*tau1);
1505 1497
1506 // Gaussian fluctuation 1498 // Gaussian fluctuation
1507 if(meanLoss > kappa*tMax || tMax < kappa*ip 1499 if(meanLoss > kappa*tMax || tMax < kappa*ipotFluct )
1508 { 1500 {
1509 siga = tMax * (1.0-0.5*beta2) * step * 1501 siga = tMax * (1.0-0.5*beta2) * step * twopi_mc2_rcl2
1510 * electronDensity / beta2 ; 1502 * electronDensity / beta2 ;
1511 1503
1512 // High velocity or negatively charged 1504 // High velocity or negatively charged particle
1513 if( beta2 > 3.0*theBohrBeta2*zeff || ch 1505 if( beta2 > 3.0*theBohrBeta2*zeff || charge < 0.0) {
1514 siga = std::sqrt( siga * chargeSquare ) ; 1506 siga = std::sqrt( siga * chargeSquare ) ;
1515 1507
1516 // Low velocity - additional ion charge flu 1508 // Low velocity - additional ion charge fluctuations according to
1517 // Q.Yang et al., NIM B61(1991)149-155. 1509 // Q.Yang et al., NIM B61(1991)149-155.
1518 } else { 1510 } else {
1519 G4double chu = theIonChuFluctuationModel->T 1511 G4double chu = theIonChuFluctuationModel->TheValue(particle, material);
1520 G4double yang = theIonYangFluctuationModel- 1512 G4double yang = theIonYangFluctuationModel->TheValue(particle, material);
1521 siga = std::sqrt( siga * (chargeSquare * ch 1513 siga = std::sqrt( siga * (chargeSquare * chu + yang)) ;
1522 } 1514 }
1523 1515
1524 do { 1516 do {
1525 loss = G4RandGauss::shoot(meanLoss,si 1517 loss = G4RandGauss::shoot(meanLoss,siga);
1526 } while (loss < 0. || loss > 2.0*meanLo 1518 } while (loss < 0. || loss > 2.0*meanLoss);
1527 return loss; 1519 return loss;
1528 } 1520 }
1529 1521
1530 // Non Gaussian fluctuation 1522 // Non Gaussian fluctuation
1531 static const G4double probLim = 0.01 ; 1523 static const G4double probLim = 0.01 ;
1532 static const G4double sumaLim = -std::log(p 1524 static const G4double sumaLim = -std::log(probLim) ;
1533 static const G4double alim = 10.; 1525 static const G4double alim = 10.;
1534 1526
1535 G4double suma,w1,w2,C,e0,lossc,w; 1527 G4double suma,w1,w2,C,e0,lossc,w;
1536 G4double a1,a2,a3; 1528 G4double a1,a2,a3;
1537 G4int p1,p2,p3; 1529 G4int p1,p2,p3;
1538 G4int nb; 1530 G4int nb;
1539 G4double corrfac, na,alfa,rfac,namean,sa,al 1531 G4double corrfac, na,alfa,rfac,namean,sa,alfa1,ea,sea;
1540 G4double dp3; 1532 G4double dp3;
1541 1533
1542 G4double f1Fluct = material->GetIonisat 1534 G4double f1Fluct = material->GetIonisation()->GetF1fluct();
1543 G4double f2Fluct = material->GetIonisat 1535 G4double f2Fluct = material->GetIonisation()->GetF2fluct();
1544 G4double e1Fluct = material->GetIonisat 1536 G4double e1Fluct = material->GetIonisation()->GetEnergy1fluct();
1545 G4double e2Fluct = material->GetIonisat 1537 G4double e2Fluct = material->GetIonisation()->GetEnergy2fluct();
1546 G4double e1LogFluct = material->GetIonisat 1538 G4double e1LogFluct = material->GetIonisation()->GetLogEnergy1fluct();
1547 G4double e2LogFluct = material->GetIonisat 1539 G4double e2LogFluct = material->GetIonisation()->GetLogEnergy2fluct();
1548 G4double rateFluct = material->GetIonisat 1540 G4double rateFluct = material->GetIonisation()->GetRateionexcfluct();
1549 G4double ipotLogFluct= material->GetIonisat 1541 G4double ipotLogFluct= material->GetIonisation()->GetLogMeanExcEnergy();
1550 1542
1551 w1 = tMax/ipotFluct; 1543 w1 = tMax/ipotFluct;
1552 w2 = std::log(2.*electron_mass_c2*tau2); 1544 w2 = std::log(2.*electron_mass_c2*tau2);
1553 1545
1554 C = meanLoss*(1.-rateFluct)/(w2-ipotLogFluc 1546 C = meanLoss*(1.-rateFluct)/(w2-ipotLogFluct-beta2);
1555 1547
1556 a1 = C*f1Fluct*(w2-e1LogFluct-beta2)/e1Fluc 1548 a1 = C*f1Fluct*(w2-e1LogFluct-beta2)/e1Fluct;
1557 a2 = C*f2Fluct*(w2-e2LogFluct-beta2)/e2Fluc 1549 a2 = C*f2Fluct*(w2-e2LogFluct-beta2)/e2Fluct;
1558 a3 = rateFluct*meanLoss*(tMax-ipotFluct)/(i 1550 a3 = rateFluct*meanLoss*(tMax-ipotFluct)/(ipotFluct*tMax*std::log(w1));
1559 if(a1 < 0.0) a1 = 0.0; 1551 if(a1 < 0.0) a1 = 0.0;
1560 if(a2 < 0.0) a2 = 0.0; 1552 if(a2 < 0.0) a2 = 0.0;
1561 if(a3 < 0.0) a3 = 0.0; 1553 if(a3 < 0.0) a3 = 0.0;
1562 1554
1563 suma = a1+a2+a3; 1555 suma = a1+a2+a3;
1564 1556
1565 loss = 0.; 1557 loss = 0.;
1566 1558
1567 1559
1568 if(suma < sumaLim) // very smal 1560 if(suma < sumaLim) // very small Step
1569 { 1561 {
1570 e0 = material->GetIonisation()->GetEner 1562 e0 = material->GetIonisation()->GetEnergy0fluct();
1571 1563
1572 if(tMax == ipotFluct) 1564 if(tMax == ipotFluct)
1573 { 1565 {
1574 a3 = meanLoss/e0; 1566 a3 = meanLoss/e0;
1575 1567
1576 if(a3>alim) 1568 if(a3>alim)
1577 { 1569 {
1578 siga=std::sqrt(a3) ; 1570 siga=std::sqrt(a3) ;
1579 p3 = std::max(0,G4int(G4RandGauss::sh 1571 p3 = std::max(0,G4int(G4RandGauss::shoot(a3,siga)+0.5));
1580 } 1572 }
1581 else 1573 else
1582 p3 = (G4int)G4Poisson(a3); << 1574 p3 = G4Poisson(a3);
1583 1575
1584 loss = p3*e0 ; 1576 loss = p3*e0 ;
1585 1577
1586 if(p3 > 0) 1578 if(p3 > 0)
1587 loss += (1.-2.*G4UniformRand())*e0 ; 1579 loss += (1.-2.*G4UniformRand())*e0 ;
1588 1580
1589 } 1581 }
1590 else 1582 else
1591 { 1583 {
1592 tMax = tMax-ipotFluct+e0 ; 1584 tMax = tMax-ipotFluct+e0 ;
1593 a3 = meanLoss*(tMax-e0)/(tMax*e0*std::log 1585 a3 = meanLoss*(tMax-e0)/(tMax*e0*std::log(tMax/e0));
1594 1586
1595 if(a3>alim) 1587 if(a3>alim)
1596 { 1588 {
1597 siga=std::sqrt(a3) ; 1589 siga=std::sqrt(a3) ;
1598 p3 = std::max(0,int(G4RandGauss::shoo 1590 p3 = std::max(0,int(G4RandGauss::shoot(a3,siga)+0.5));
1599 } 1591 }
1600 else 1592 else
1601 p3 = (G4int)G4Poisson(a3); << 1593 p3 = G4Poisson(a3);
1602 1594
1603 if(p3 > 0) 1595 if(p3 > 0)
1604 { 1596 {
1605 w = (tMax-e0)/tMax ; 1597 w = (tMax-e0)/tMax ;
1606 if(p3 > nmaxCont2) 1598 if(p3 > nmaxCont2)
1607 { 1599 {
1608 dp3 = G4float(p3) ; 1600 dp3 = G4float(p3) ;
1609 corrfac = dp3/G4float(nmaxCont2) ; 1601 corrfac = dp3/G4float(nmaxCont2) ;
1610 p3 = G4int(nmaxCont2) ; << 1602 p3 = nmaxCont2 ;
1611 } 1603 }
1612 else 1604 else
1613 corrfac = 1. ; 1605 corrfac = 1. ;
1614 1606
1615 for(G4int i=0; i<p3; i++) loss += 1./ 1607 for(G4int i=0; i<p3; i++) loss += 1./(1.-w*G4UniformRand()) ;
1616 loss *= e0*corrfac ; 1608 loss *= e0*corrfac ;
1617 } 1609 }
1618 } 1610 }
1619 } 1611 }
1620 1612
1621 else // not so 1613 else // not so small Step
1622 { 1614 {
1623 // excitation type 1 1615 // excitation type 1
1624 if(a1>alim) 1616 if(a1>alim)
1625 { 1617 {
1626 siga=std::sqrt(a1) ; 1618 siga=std::sqrt(a1) ;
1627 p1 = std::max(0,G4int(G4RandGauss::shoot( 1619 p1 = std::max(0,G4int(G4RandGauss::shoot(a1,siga)+0.5));
1628 } 1620 }
1629 else 1621 else
1630 p1 = (G4int)G4Poisson(a1); << 1622 p1 = G4Poisson(a1);
1631 1623
1632 // excitation type 2 1624 // excitation type 2
1633 if(a2>alim) 1625 if(a2>alim)
1634 { 1626 {
1635 siga=std::sqrt(a2) ; 1627 siga=std::sqrt(a2) ;
1636 p2 = std::max(0,G4int(G4RandGauss::shoot( 1628 p2 = std::max(0,G4int(G4RandGauss::shoot(a2,siga)+0.5));
1637 } 1629 }
1638 else 1630 else
1639 p2 = (G4int)G4Poisson(a2); << 1631 p2 = G4Poisson(a2);
1640 1632
1641 loss = p1*e1Fluct+p2*e2Fluct; 1633 loss = p1*e1Fluct+p2*e2Fluct;
1642 1634
1643 // smearing to avoid unphysical peaks 1635 // smearing to avoid unphysical peaks
1644 if(p2 > 0) 1636 if(p2 > 0)
1645 loss += (1.-2.*G4UniformRand())*e2Flu 1637 loss += (1.-2.*G4UniformRand())*e2Fluct;
1646 else if (loss>0.) 1638 else if (loss>0.)
1647 loss += (1.-2.*G4UniformRand())*e1Flu 1639 loss += (1.-2.*G4UniformRand())*e1Fluct;
1648 1640
1649 // ionisation ......................... 1641 // ionisation .......................................
1650 if(a3 > 0.) 1642 if(a3 > 0.)
1651 { 1643 {
1652 if(a3>alim) 1644 if(a3>alim)
1653 { 1645 {
1654 siga=std::sqrt(a3) ; 1646 siga=std::sqrt(a3) ;
1655 p3 = std::max(0,G4int(G4RandGauss::sh 1647 p3 = std::max(0,G4int(G4RandGauss::shoot(a3,siga)+0.5));
1656 } 1648 }
1657 else 1649 else
1658 p3 = (G4int)G4Poisson(a3); << 1650 p3 = G4Poisson(a3);
1659 1651
1660 lossc = 0.; 1652 lossc = 0.;
1661 if(p3 > 0) 1653 if(p3 > 0)
1662 { 1654 {
1663 na = 0.; 1655 na = 0.;
1664 alfa = 1.; 1656 alfa = 1.;
1665 if (p3 > nmaxCont2) 1657 if (p3 > nmaxCont2)
1666 { 1658 {
1667 dp3 = G4float(p3); 1659 dp3 = G4float(p3);
1668 rfac = dp3/(G4float(nmaxCont2)+dp 1660 rfac = dp3/(G4float(nmaxCont2)+dp3);
1669 namean = G4float(p3)*rfac; 1661 namean = G4float(p3)*rfac;
1670 sa = G4float(nmaxCont1)*rfac; 1662 sa = G4float(nmaxCont1)*rfac;
1671 na = G4RandGauss::shoot(namean, 1663 na = G4RandGauss::shoot(namean,sa);
1672 if (na > 0.) 1664 if (na > 0.)
1673 { 1665 {
1674 alfa = w1*G4float(nmaxCont2+p3)/ 1666 alfa = w1*G4float(nmaxCont2+p3)/
1675 (w1*G4float(nmaxCont2)+G4float(p3)); 1667 (w1*G4float(nmaxCont2)+G4float(p3));
1676 alfa1 = alfa*std::log(alfa)/(alfa- 1668 alfa1 = alfa*std::log(alfa)/(alfa-1.);
1677 ea = na*ipotFluct*alfa1; 1669 ea = na*ipotFluct*alfa1;
1678 sea = ipotFluct*std::sqrt(na*(al 1670 sea = ipotFluct*std::sqrt(na*(alfa-alfa1*alfa1));
1679 lossc += G4RandGauss::shoot(ea,sea) 1671 lossc += G4RandGauss::shoot(ea,sea);
1680 } 1672 }
1681 } 1673 }
1682 1674
1683 nb = G4int(G4float(p3)-na); 1675 nb = G4int(G4float(p3)-na);
1684 if (nb > 0) 1676 if (nb > 0)
1685 { 1677 {
1686 w2 = alfa*ipotFluct; 1678 w2 = alfa*ipotFluct;
1687 w = (tMax-w2)/tMax; 1679 w = (tMax-w2)/tMax;
1688 for (G4int k=0; k<nb; k++) lossc += w2/ 1680 for (G4int k=0; k<nb; k++) lossc += w2/(1.-w*G4UniformRand());
1689 } 1681 }
1690 } 1682 }
1691 loss += lossc; 1683 loss += lossc;
1692 } 1684 }
1693 } 1685 }
1694 1686
1695 return loss ; 1687 return loss ;
1696 } 1688 }
1697 1689
1698 1690
1699 1691
1700 void G4hImpactIonisation::SetCutForSecondaryP 1692 void G4hImpactIonisation::SetCutForSecondaryPhotons(G4double cut)
1701 { 1693 {
1702 minGammaEnergy = cut; 1694 minGammaEnergy = cut;
1703 } 1695 }
1704 1696
1705 1697
1706 1698
1707 void G4hImpactIonisation::SetCutForAugerElect 1699 void G4hImpactIonisation::SetCutForAugerElectrons(G4double cut)
1708 { 1700 {
1709 minElectronEnergy = cut; 1701 minElectronEnergy = cut;
1710 } 1702 }
1711 1703
1712 1704
1713 1705
1714 void G4hImpactIonisation::ActivateAugerElectr 1706 void G4hImpactIonisation::ActivateAugerElectronProduction(G4bool val)
1715 { 1707 {
1716 atomicDeexcitation.ActivateAugerElectronPro 1708 atomicDeexcitation.ActivateAugerElectronProduction(val);
1717 } 1709 }
1718 1710
1719 1711
1720 1712
1721 void G4hImpactIonisation::PrintInfoDefinition 1713 void G4hImpactIonisation::PrintInfoDefinition() const
1722 { 1714 {
1723 G4String comments = " Knock-on electron cr 1715 G4String comments = " Knock-on electron cross sections . ";
1724 comments += "\n Good description abo 1716 comments += "\n Good description above the mean excitation energy.\n";
1725 comments += " Delta ray energy sampl 1717 comments += " Delta ray energy sampled from differential Xsection.";
1726 1718
1727 G4cout << G4endl << GetProcessName() << ": 1719 G4cout << G4endl << GetProcessName() << ": " << comments
1728 << "\n PhysicsTables from " < 1720 << "\n PhysicsTables from " << LowestKineticEnergy / eV << " eV "
1729 << " to " << HighestKineticEnergy / 1721 << " to " << HighestKineticEnergy / TeV << " TeV "
1730 << " in " << TotBin << " bins." 1722 << " in " << TotBin << " bins."
1731 << "\n Electronic stopping power mo 1723 << "\n Electronic stopping power model is "
1732 << protonTable 1724 << protonTable
1733 << "\n from " << protonLowEne 1725 << "\n from " << protonLowEnergy / keV << " keV "
1734 << " to " << protonHighEnergy / MeV 1726 << " to " << protonHighEnergy / MeV << " MeV " << "." << G4endl ;
1735 G4cout << "\n Parametrisation model 1727 G4cout << "\n Parametrisation model for antiprotons is "
1736 << antiprotonTable 1728 << antiprotonTable
1737 << "\n from " << antiprotonLo 1729 << "\n from " << antiprotonLowEnergy / keV << " keV "
1738 << " to " << antiprotonHighEnergy / 1730 << " to " << antiprotonHighEnergy / MeV << " MeV " << "." << G4endl ;
1739 if(theBarkas){ 1731 if(theBarkas){
1740 G4cout << " Parametrization of the 1732 G4cout << " Parametrization of the Barkas effect is switched on."
1741 << G4endl ; 1733 << G4endl ;
1742 } 1734 }
1743 if(nStopping) { 1735 if(nStopping) {
1744 G4cout << " Nuclear stopping power 1736 G4cout << " Nuclear stopping power model is " << theNuclearTable
1745 << G4endl ; 1737 << G4endl ;
1746 } 1738 }
1747 1739
1748 G4bool printHead = true; 1740 G4bool printHead = true;
1749 1741
1750 const G4ProductionCutsTable* theCoupleTable 1742 const G4ProductionCutsTable* theCoupleTable=
1751 G4ProductionCutsTable::GetProductionCutsT 1743 G4ProductionCutsTable::GetProductionCutsTable();
1752 G4int numOfCouples = (G4int)theCoupleTable- << 1744 size_t numOfCouples = theCoupleTable->GetTableSize();
1753 1745
1754 // loop for materials 1746 // loop for materials
1755 1747
1756 for (G4int j=0 ; j < numOfCouples; ++j) { << 1748 for (size_t j=0 ; j < numOfCouples; j++) {
1757 1749
1758 const G4MaterialCutsCouple* couple = theC 1750 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
1759 const G4Material* material= couple->GetMa 1751 const G4Material* material= couple->GetMaterial();
1760 G4double deltaCutNow = cutForDelta[(coupl 1752 G4double deltaCutNow = cutForDelta[(couple->GetIndex())] ;
1761 G4double excitationEnergy = material->Get 1753 G4double excitationEnergy = material->GetIonisation()->GetMeanExcitationEnergy();
1762 1754
1763 if(excitationEnergy > deltaCutNow) { 1755 if(excitationEnergy > deltaCutNow) {
1764 if(printHead) { 1756 if(printHead) {
1765 printHead = false ; 1757 printHead = false ;
1766 1758
1767 G4cout << " material 1759 G4cout << " material min.delta energy(keV) " << G4endl;
1768 G4cout << G4endl; 1760 G4cout << G4endl;
1769 } 1761 }
1770 1762
1771 G4cout << std::setw(20) << material->Ge 1763 G4cout << std::setw(20) << material->GetName()
1772 << std::setw(15) << excitationEnergy/k 1764 << std::setw(15) << excitationEnergy/keV << G4endl;
1773 } 1765 }
1774 } 1766 }
1775 } 1767 }
1776 1768