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 9.4.p2)


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