Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/utils/src/G4VAtomDeexcitation.cc

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

Diff markup

Differences between /processes/electromagnetic/utils/src/G4VAtomDeexcitation.cc (Version 11.3.0) and /processes/electromagnetic/utils/src/G4VAtomDeexcitation.cc (Version 11.0.p4)


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