| 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 // GEANT4 Class class file 29 // GEANT4 Class class file
30 // 30 //
31 // 31 //
32 // File name: G4VAtomDeexcitation 32 // File name: G4VAtomDeexcitation
33 // 33 //
34 // Author: Alfonso Mantero & Vladimir I 34 // Author: Alfonso Mantero & Vladimir Ivanchenko
35 // 35 //
36 // Creation date: 21.04.2010 36 // Creation date: 21.04.2010
37 // 37 //
38 // Modifications: 38 // Modifications:
39 // 39 //
40 // Class Description: 40 // Class Description:
41 // 41 //
42 // Abstract interface to energy loss models 42 // Abstract interface to energy loss models
43 43
44 // ------------------------------------------- 44 // -------------------------------------------------------------------
45 // 45 //
46 46
47 #include "G4VAtomDeexcitation.hh" 47 #include "G4VAtomDeexcitation.hh"
48 #include "G4SystemOfUnits.hh" 48 #include "G4SystemOfUnits.hh"
49 #include "G4EmParameters.hh" <<
50 #include "G4ParticleDefinition.hh" 49 #include "G4ParticleDefinition.hh"
51 #include "G4DynamicParticle.hh" 50 #include "G4DynamicParticle.hh"
52 #include "G4Step.hh" 51 #include "G4Step.hh"
53 #include "G4Region.hh" 52 #include "G4Region.hh"
54 #include "G4RegionStore.hh" 53 #include "G4RegionStore.hh"
55 #include "G4MaterialCutsCouple.hh" 54 #include "G4MaterialCutsCouple.hh"
56 #include "G4MaterialCutsCouple.hh" 55 #include "G4MaterialCutsCouple.hh"
57 #include "G4Material.hh" 56 #include "G4Material.hh"
58 #include "G4Element.hh" 57 #include "G4Element.hh"
59 #include "G4ElementVector.hh" 58 #include "G4ElementVector.hh"
60 #include "Randomize.hh" 59 #include "Randomize.hh"
61 #include "G4VParticleChange.hh" 60 #include "G4VParticleChange.hh"
62 #include "G4EmSecondaryParticleType.hh" << 61 #include "G4PhysicsModelCatalog.hh"
63 #include "G4Gamma.hh" 62 #include "G4Gamma.hh"
64 #include "G4Log.hh" <<
65 63
66 #ifdef G4MULTITHREADED 64 #ifdef G4MULTITHREADED
67 G4Mutex G4VAtomDeexcitation::atomDeexcitatio 65 G4Mutex G4VAtomDeexcitation::atomDeexcitationMutex = G4MUTEX_INITIALIZER;
68 #endif 66 #endif
69 67
70 //....oooOO0OOooo........oooOO0OOooo........oo 68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71 69
>> 70 G4int G4VAtomDeexcitation::pixeIDg = -1;
>> 71 G4int G4VAtomDeexcitation::pixeIDe = -1;
>> 72
72 G4VAtomDeexcitation::G4VAtomDeexcitation(const 73 G4VAtomDeexcitation::G4VAtomDeexcitation(const G4String& modname)
73 : name(modname) << 74 : verbose(1), name(modname), isActive(false), flagAuger(false),
>> 75 flagAugerCascade(false), flagPIXE(false), ignoreCuts(false),
>> 76 isActiveLocked(false), isAugerLocked(false),
>> 77 isAugerCascadeLocked(false), isPIXELocked(false)
74 { 78 {
>> 79 theParameters = G4EmParameters::Instance();
75 vdyn.reserve(5); 80 vdyn.reserve(5);
76 theCoupleTable = nullptr; 81 theCoupleTable = nullptr;
>> 82 G4String gg = "gammaPIXE";
>> 83 G4String ee = "e-PIXE";
>> 84 if(pixeIDg < 0) {
>> 85 #ifdef G4MULTITHREADED
>> 86 G4MUTEXLOCK(&atomDeexcitationMutex);
>> 87 if(pixeIDg < 0) {
>> 88 #endif
>> 89 pixeIDg = G4PhysicsModelCatalog::Register(gg);
>> 90 pixeIDe = G4PhysicsModelCatalog::Register(ee);
>> 91 #ifdef G4MULTITHREADED
>> 92 }
>> 93 G4MUTEXUNLOCK(&atomDeexcitationMutex);
>> 94 #endif
>> 95 }
77 gamma = G4Gamma::Gamma(); 96 gamma = G4Gamma::Gamma();
78 } 97 }
79 98
80 //....oooOO0OOooo........oooOO0OOooo........oo 99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
81 100
82 G4VAtomDeexcitation::~G4VAtomDeexcitation() = << 101 G4VAtomDeexcitation::~G4VAtomDeexcitation()
>> 102 {}
83 103
84 //....oooOO0OOooo........oooOO0OOooo........oo 104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
85 105
86 void G4VAtomDeexcitation::InitialiseAtomicDeex 106 void G4VAtomDeexcitation::InitialiseAtomicDeexcitation()
87 { 107 {
88 G4EmParameters* theParameters = G4EmParamete <<
89 theParameters->DefineRegParamForDeex(this); 108 theParameters->DefineRegParamForDeex(this);
90 109
91 // Define list of couples 110 // Define list of couples
92 theCoupleTable = G4ProductionCutsTable::GetP 111 theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
93 nCouples = (G4int)theCoupleTable->GetTableSi << 112 G4int numOfCouples = theCoupleTable->GetTableSize();
94 113
95 // needed for unit tests 114 // needed for unit tests
96 std::size_t nn = std::max(nCouples, 1); << 115 size_t nn = std::max(numOfCouples, 1);
97 if(activeDeexcitationMedia.size() != nn) { 116 if(activeDeexcitationMedia.size() != nn) {
98 activeDeexcitationMedia.resize(nn, false); 117 activeDeexcitationMedia.resize(nn, false);
99 activeAugerMedia.resize(nn, false); 118 activeAugerMedia.resize(nn, false);
100 activePIXEMedia.resize(nn, false); 119 activePIXEMedia.resize(nn, false);
101 } 120 }
102 if(activeZ.size() != 93) { activeZ.resize(93 121 if(activeZ.size() != 93) { activeZ.resize(93, false); }
103 122
104 // initialisation of flags and options 123 // initialisation of flags and options
105 // normally there is no locksed flags 124 // normally there is no locksed flags
106 if(!isActiveLocked) { isActive = theParamet << 125 if(!isActiveLocked) { isActive = theParameters->Fluo(); }
107 if(!isAugerLocked) { flagAuger = theParamet << 126 if(!isAugerLocked) { flagAuger = theParameters->Auger(); }
108 if(!isPIXELocked) { flagPIXE = theParamet << 127 if(!isAugerCascadeLocked) { flagAugerCascade = theParameters->AugerCascade(); }
>> 128 if(!isPIXELocked) { flagPIXE = theParameters->Pixe(); }
109 ignoreCuts = theParameters->DeexcitationIgno 129 ignoreCuts = theParameters->DeexcitationIgnoreCut();
110 130
111 // Define list of regions 131 // Define list of regions
112 std::size_t nRegions = deRegions.size(); << 132 size_t nRegions = deRegions.size();
113 // check if deexcitation is active for the g 133 // check if deexcitation is active for the given run
114 if(!isActive && 0 == nRegions) { return; } 134 if(!isActive && 0 == nRegions) { return; }
115 135
116 // if no active regions add a world 136 // if no active regions add a world
117 if(0 == nRegions) { 137 if(0 == nRegions) {
118 SetDeexcitationActiveRegion("World",isActi 138 SetDeexcitationActiveRegion("World",isActive,flagAuger,flagPIXE);
119 nRegions = deRegions.size(); 139 nRegions = deRegions.size();
120 } 140 }
121 141
122 if(0 < verbose) { 142 if(0 < verbose) {
123 G4cout << G4endl; 143 G4cout << G4endl;
124 G4cout << "### === Deexcitation model " < 144 G4cout << "### === Deexcitation model " << name
125 << " is activated for " << nRegions 145 << " is activated for " << nRegions;
126 if(1 == nRegions) { G4cout << " region:" < 146 if(1 == nRegions) { G4cout << " region:" << G4endl; }
127 else { G4cout << " regions:" 147 else { G4cout << " regions:" << G4endl;}
128 } 148 }
129 149
130 // Identify active media 150 // Identify active media
131 const G4RegionStore* regionStore = G4RegionS << 151 G4RegionStore* regionStore = G4RegionStore::GetInstance();
132 for(std::size_t j=0; j<nRegions; ++j) { << 152 for(size_t j=0; j<nRegions; ++j) {
133 const G4Region* reg = regionStore->GetRegi 153 const G4Region* reg = regionStore->GetRegion(activeRegions[j], false);
134 if(nullptr != reg && 0 < nCouples) { << 154 if(reg && 0 < numOfCouples) {
135 const G4ProductionCuts* rpcuts = reg->Ge 155 const G4ProductionCuts* rpcuts = reg->GetProductionCuts();
136 if(0 < verbose) { 156 if(0 < verbose) {
137 G4cout << " " << activeRegion 157 G4cout << " " << activeRegions[j]
138 << " " << deRegions[j] << " << 158 << " " << deRegions[j] << " " << AugerRegions[j]
139 << " " << PIXERegions[j] << G4 << 159 << " " << PIXERegions[j] << G4endl;
140 } 160 }
141 for(G4int i=0; i<nCouples; ++i) { << 161 for(G4int i=0; i<numOfCouples; ++i) {
142 const G4MaterialCutsCouple* couple = 162 const G4MaterialCutsCouple* couple =
143 theCoupleTable->GetMaterialCutsCoupl 163 theCoupleTable->GetMaterialCutsCouple(i);
144 if (couple->GetProductionCuts() == rpc 164 if (couple->GetProductionCuts() == rpcuts) {
145 activeDeexcitationMedia[i] = deRegio 165 activeDeexcitationMedia[i] = deRegions[j];
146 activeAugerMedia[i] = AugerRegions[j 166 activeAugerMedia[i] = AugerRegions[j];
147 activePIXEMedia[i] = PIXERegions[j]; 167 activePIXEMedia[i] = PIXERegions[j];
148 } 168 }
149 } 169 }
150 } 170 }
151 } 171 }
152 std::size_t nelm = G4Element::GetNumberOfEle << 172 G4int nelm = G4Element::GetNumberOfElements();
153 //G4cout << nelm << G4endl; 173 //G4cout << nelm << G4endl;
154 for(std::size_t k=0; k<nelm; ++k) { << 174 for(G4int k=0; k<nelm; ++k) {
155 G4int Z = (*(G4Element::GetElementTable()) 175 G4int Z = (*(G4Element::GetElementTable()))[k]->GetZasInt();
156 if(Z > 5 && Z < 93) { 176 if(Z > 5 && Z < 93) {
157 activeZ[Z] = true; 177 activeZ[Z] = true;
158 //G4cout << "!!! Active de-excitation Z= 178 //G4cout << "!!! Active de-excitation Z= " << Z << G4endl;
159 } 179 }
160 } 180 }
161 181
162 // Initialise derived class 182 // Initialise derived class
163 InitialiseForNewRun(); 183 InitialiseForNewRun();
164 184
165 if(0 < verbose && flagAuger) { 185 if(0 < verbose && flagAuger) {
166 G4cout << "### === Auger flag: " << flagA << 186 G4cout << "### === Auger cascade flag: " << flagAugerCascade
167 << G4endl; << 187 << G4endl;
168 } 188 }
169 if(0 < verbose) { 189 if(0 < verbose) {
170 G4cout << "### === Ignore cuts flag: " 190 G4cout << "### === Ignore cuts flag: " << ignoreCuts
171 << G4endl; << 191 << G4endl;
172 } 192 }
173 if(0 < verbose && flagPIXE) { 193 if(0 < verbose && flagPIXE) {
174 G4cout << "### === PIXE model for hadrons 194 G4cout << "### === PIXE model for hadrons: "
175 << theParameters->PIXECrossSectionM 195 << theParameters->PIXECrossSectionModel()
176 << G4endl; 196 << G4endl;
177 G4cout << "### === PIXE model for e+-: 197 G4cout << "### === PIXE model for e+-: "
178 << theParameters->PIXEElectronCross 198 << theParameters->PIXEElectronCrossSectionModel()
179 << G4endl; 199 << G4endl;
180 } 200 }
181 } 201 }
182 202
183 //....oooOO0OOooo........oooOO0OOooo........oo 203 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
184 204
185 void 205 void
186 G4VAtomDeexcitation::SetDeexcitationActiveRegi 206 G4VAtomDeexcitation::SetDeexcitationActiveRegion(const G4String& rname,
187 207 G4bool valDeexcitation,
188 208 G4bool valAuger,
189 209 G4bool valPIXE)
190 { 210 {
191 // no PIXE in parallel world 211 // no PIXE in parallel world
192 if(rname == "DefaultRegionForParallelWorld") 212 if(rname == "DefaultRegionForParallelWorld") { return; }
193 213
194 G4String ss = rname; 214 G4String ss = rname;
195 /* 215 /*
196 G4cout << "### G4VAtomDeexcitation::SetDeexc 216 G4cout << "### G4VAtomDeexcitation::SetDeexcitationActiveRegion " << ss
197 << " " << valDeexcitation << " " << 217 << " " << valDeexcitation << " " << valAuger
198 << " " << valPIXE << G4endl; 218 << " " << valPIXE << G4endl;
199 */ 219 */
200 if(ss == "world" || ss == "World" || ss == " 220 if(ss == "world" || ss == "World" || ss == "WORLD") {
201 ss = "DefaultRegionForTheWorld"; 221 ss = "DefaultRegionForTheWorld";
202 } 222 }
203 std::size_t n = deRegions.size(); << 223 size_t n = deRegions.size();
204 for(std::size_t i=0; i<n; ++i) { << 224 for(size_t i=0; i<n; ++i) {
205 225
206 // Region already exist 226 // Region already exist
207 if(ss == activeRegions[i]) { 227 if(ss == activeRegions[i]) {
208 deRegions[i] = valDeexcitation; 228 deRegions[i] = valDeexcitation;
209 AugerRegions[i] = valAuger; 229 AugerRegions[i] = valAuger;
210 PIXERegions[i] = valPIXE; 230 PIXERegions[i] = valPIXE;
211 return; 231 return;
212 } 232 }
213 } 233 }
214 // New region 234 // New region
215 activeRegions.push_back(ss); 235 activeRegions.push_back(ss);
216 deRegions.push_back(valDeexcitation); 236 deRegions.push_back(valDeexcitation);
217 AugerRegions.push_back(valAuger); 237 AugerRegions.push_back(valAuger);
218 PIXERegions.push_back(valPIXE); 238 PIXERegions.push_back(valPIXE);
219 239
220 // if de-excitation defined for the world vo 240 // if de-excitation defined for the world volume
221 // it should be active for all G4Regions 241 // it should be active for all G4Regions
222 if(ss == "DefaultRegionForTheWorld") { 242 if(ss == "DefaultRegionForTheWorld") {
223 G4RegionStore* regions = G4RegionStore::Ge 243 G4RegionStore* regions = G4RegionStore::GetInstance();
224 std::size_t nn = regions->size(); << 244 G4int nn = regions->size();
225 for(std::size_t i=0; i<nn; ++i) { << 245 for(G4int i=0; i<nn; ++i) {
226 if(ss == (*regions)[i]->GetName()) { con 246 if(ss == (*regions)[i]->GetName()) { continue; }
227 SetDeexcitationActiveRegion((*regions)[i 247 SetDeexcitationActiveRegion((*regions)[i]->GetName(), valDeexcitation,
228 valAuger, va 248 valAuger, valPIXE);
229 249
230 } 250 }
231 } 251 }
232 } 252 }
233 253
234 void G4VAtomDeexcitation::GenerateParticles(st << 254 void G4VAtomDeexcitation::GenerateParticles(std::vector<G4DynamicParticle*>* v,
235 co << 255 const G4AtomicShell* as,
236 G4 << 256 G4int Z, G4int idx)
237 { 257 {
238 G4double gCut = DBL_MAX; 258 G4double gCut = DBL_MAX;
239 if(ignoreCuts) { 259 if(ignoreCuts) {
240 gCut = 0.0; 260 gCut = 0.0;
241 } else if (nullptr != theCoupleTable) { << 261 } else if (theCoupleTable) {
242 gCut = (*(theCoupleTable->GetEnergyCutsVec 262 gCut = (*(theCoupleTable->GetEnergyCutsVector(0)))[idx];
243 } 263 }
244 if(gCut < as->BindingEnergy()) { 264 if(gCut < as->BindingEnergy()) {
245 G4double eCut = DBL_MAX; 265 G4double eCut = DBL_MAX;
246 if(CheckAugerActiveRegion(idx)) { 266 if(CheckAugerActiveRegion(idx)) {
247 if(ignoreCuts) { 267 if(ignoreCuts) {
248 eCut = 0.0; 268 eCut = 0.0;
249 } else if (nullptr != theCoupleTable) { << 269 } else if (theCoupleTable) {
250 eCut = (*(theCoupleTable->GetEnergyCut 270 eCut = (*(theCoupleTable->GetEnergyCutsVector(1)))[idx];
251 } 271 }
252 } 272 }
253 GenerateParticles(v, as, Z, gCut, eCut); 273 GenerateParticles(v, as, Z, gCut, eCut);
254 } 274 }
255 } 275 }
256 276
257 //....oooOO0OOooo........oooOO0OOooo........oo 277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
258 278
259 void 279 void
260 G4VAtomDeexcitation::AlongStepDeexcitation(std 280 G4VAtomDeexcitation::AlongStepDeexcitation(std::vector<G4Track*>& tracks,
261 con 281 const G4Step& step,
262 G4d 282 G4double& eLossMax,
263 G4i 283 G4int coupleIndex)
264 { 284 {
265 G4double truelength = step.GetStepLength(); 285 G4double truelength = step.GetStepLength();
266 if(!flagPIXE && !activePIXEMedia[coupleIndex 286 if(!flagPIXE && !activePIXEMedia[coupleIndex]) { return; }
267 if(eLossMax <= 0.0 || truelength <= 0.0) 287 if(eLossMax <= 0.0 || truelength <= 0.0) { return; }
268 288
269 // step parameters 289 // step parameters
270 const G4StepPoint* preStep = step.GetPreStep 290 const G4StepPoint* preStep = step.GetPreStepPoint();
271 const G4ThreeVector prePos = preStep->GetPos << 291 G4ThreeVector prePos = preStep->GetPosition();
272 const G4ThreeVector delta = step.GetPostStep << 292 G4ThreeVector delta = step.GetPostStepPoint()->GetPosition() - prePos;
273 const G4double preTime = preStep->GetGlobalT << 293 G4double preTime = preStep->GetGlobalTime();
274 const G4double dt = step.GetPostStepPoint()- << 294 G4double dt = step.GetPostStepPoint()->GetGlobalTime() - preTime;
275 295
276 // particle parameters 296 // particle parameters
277 const G4Track* track = step.GetTrack(); 297 const G4Track* track = step.GetTrack();
278 const G4ParticleDefinition* part = track->Ge 298 const G4ParticleDefinition* part = track->GetDefinition();
279 G4double ekin = preStep->GetKineticEnergy(); 299 G4double ekin = preStep->GetKineticEnergy();
280 300
281 // media parameters 301 // media parameters
282 G4double gCut = (*theCoupleTable->GetEnergyC 302 G4double gCut = (*theCoupleTable->GetEnergyCutsVector(0))[coupleIndex];
283 if(ignoreCuts) { gCut = 0.0; } 303 if(ignoreCuts) { gCut = 0.0; }
284 G4double eCut = DBL_MAX; 304 G4double eCut = DBL_MAX;
285 if(CheckAugerActiveRegion(coupleIndex)) { 305 if(CheckAugerActiveRegion(coupleIndex)) {
286 eCut = (*theCoupleTable->GetEnergyCutsVect 306 eCut = (*theCoupleTable->GetEnergyCutsVector(1))[coupleIndex];
287 if(ignoreCuts) { eCut = 0.0; } 307 if(ignoreCuts) { eCut = 0.0; }
288 } 308 }
289 309
290 //G4cout<<"!Sample PIXE gCut(MeV)= "<<gCut<< 310 //G4cout<<"!Sample PIXE gCut(MeV)= "<<gCut<<" eCut(MeV)= "<<eCut
291 // <<" Ekin(MeV)= " << ekin/MeV << G4 311 // <<" Ekin(MeV)= " << ekin/MeV << G4endl;
292 312
293 const G4Material* material = preStep->GetMat 313 const G4Material* material = preStep->GetMaterial();
294 const G4ElementVector* theElementVector = ma 314 const G4ElementVector* theElementVector = material->GetElementVector();
295 const G4double* theAtomNumDensityVector = 315 const G4double* theAtomNumDensityVector =
296 material->GetVecNbOfAtomsPerVolume(); 316 material->GetVecNbOfAtomsPerVolume();
297 const std::size_t nelm = material->GetNumber << 317 G4int nelm = material->GetNumberOfElements();
298 318
299 // loop over deexcitations 319 // loop over deexcitations
300 for(std::size_t i=0; i<nelm; ++i) { << 320 for(G4int i=0; i<nelm; ++i) {
301 G4int Z = (*theElementVector)[i]->GetZasIn 321 G4int Z = (*theElementVector)[i]->GetZasInt();
302 if(activeZ[Z] && Z < 93) { 322 if(activeZ[Z] && Z < 93) {
303 G4int nshells = 323 G4int nshells =
304 std::min(9,(*theElementVector)[i]->Get 324 std::min(9,(*theElementVector)[i]->GetNbOfAtomicShells());
305 G4double rho = truelength*theAtomNumDens 325 G4double rho = truelength*theAtomNumDensityVector[i];
306 //G4cout<<" Z "<< Z <<" is active x(m 326 //G4cout<<" Z "<< Z <<" is active x(mm)= " << truelength/mm << G4endl;
307 for(G4int ii=0; ii<nshells; ++ii) { 327 for(G4int ii=0; ii<nshells; ++ii) {
308 auto as = (G4AtomicShellEnumerator)(ii << 328 G4AtomicShellEnumerator as = G4AtomicShellEnumerator(ii);
309 const G4AtomicShell* shell = GetAtomic 329 const G4AtomicShell* shell = GetAtomicShell(Z, as);
310 const G4double bindingEnergy = shell-> << 330 G4double bindingEnergy = shell->BindingEnergy();
311 331
312 if(gCut > bindingEnergy) { break; } 332 if(gCut > bindingEnergy) { break; }
313 333
314 if(eLossMax > bindingEnergy) { 334 if(eLossMax > bindingEnergy) {
315 G4double sig = rho* 335 G4double sig = rho*
316 GetShellIonisationCrossSectionPerA 336 GetShellIonisationCrossSectionPerAtom(part, Z, as, ekin, material);
317 337
318 // mfp is mean free path in units of 338 // mfp is mean free path in units of step size
319 if(sig > 0.0) { 339 if(sig > 0.0) {
320 G4double mfp = 1.0/sig; 340 G4double mfp = 1.0/sig;
321 G4double stot = 0.0; 341 G4double stot = 0.0;
322 //G4cout << " Shell " << ii << " m 342 //G4cout << " Shell " << ii << " mfp(mm)= " << mfp/mm << G4endl;
323 // sample ionisation points 343 // sample ionisation points
324 do { 344 do {
325 stot -= mfp*G4Log(G4UniformRand( << 345 stot -= mfp*std::log(G4UniformRand());
326 if( stot > 1.0 || eLossMax < bin 346 if( stot > 1.0 || eLossMax < bindingEnergy) { break; }
327 // sample deexcitation 347 // sample deexcitation
328 vdyn.clear(); 348 vdyn.clear();
329 GenerateParticles(&vdyn, shell, 349 GenerateParticles(&vdyn, shell, Z, gCut, eCut);
330 std::size_t nsec = vdyn.size(); << 350 G4int nsec = vdyn.size();
331 if(nsec > 0) { 351 if(nsec > 0) {
332 G4ThreeVector r = prePos + st 352 G4ThreeVector r = prePos + stot*delta;
333 G4double time = preTime + st 353 G4double time = preTime + stot*dt;
334 for(std::size_t j=0; j<nsec; + << 354 for(G4int j=0; j<nsec; ++j) {
335 G4DynamicParticle* dp = vdyn 355 G4DynamicParticle* dp = vdyn[j];
336 G4double e = dp->GetKineticE 356 G4double e = dp->GetKineticEnergy();
337 357
338 // save new secondary if the 358 // save new secondary if there is enough energy
339 if(eLossMax >= e) { 359 if(eLossMax >= e) {
340 eLossMax -= e; 360 eLossMax -= e;
341 G4Track* t = new G4Track(d 361 G4Track* t = new G4Track(dp, time, r);
342 362
343 // defined secondary type 363 // defined secondary type
344 if(dp->GetDefinition() == 364 if(dp->GetDefinition() == gamma) {
345 t->SetCreatorModelID(_Ga << 365 t->SetCreatorModelIndex(pixeIDg);
346 } else { 366 } else {
347 t->SetCreatorModelID(_eP << 367 t->SetCreatorModelIndex(pixeIDe);
348 } 368 }
>> 369
349 tracks.push_back(t); 370 tracks.push_back(t);
350 } else { 371 } else {
351 delete dp; 372 delete dp;
352 } 373 }
353 } 374 }
354 } 375 }
355 // Loop checking, 03-Aug-2015, V 376 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
356 } while (stot < 1.0); 377 } while (stot < 1.0);
357 } 378 }
358 } 379 }
359 } 380 }
360 } 381 }
361 } 382 }
362 return; 383 return;
363 } 384 }
364 385
365 //....oooOO0OOooo........oooOO0OOooo........oo 386 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
366 387