Geant4  10.01.p02
G4VAtomDeexcitation.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4VAtomDeexcitation.cc 86805 2014-11-18 13:25:55Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class class file
31 //
32 //
33 // File name: G4VAtomDeexcitation
34 //
35 // Author: Alfonso Mantero & Vladimir Ivanchenko
36 //
37 // Creation date: 21.04.2010
38 //
39 // Modifications:
40 //
41 // Class Description:
42 //
43 // Abstract interface to energy loss models
44 
45 // -------------------------------------------------------------------
46 //
47 
48 #include "G4VAtomDeexcitation.hh"
49 #include "G4SystemOfUnits.hh"
50 #include "G4ParticleDefinition.hh"
51 #include "G4DynamicParticle.hh"
52 #include "G4Step.hh"
53 #include "G4Region.hh"
54 #include "G4RegionStore.hh"
55 #include "G4MaterialCutsCouple.hh"
56 #include "G4MaterialCutsCouple.hh"
57 #include "G4Material.hh"
58 #include "G4Element.hh"
59 #include "G4ElementVector.hh"
60 #include "Randomize.hh"
61 #include "G4VParticleChange.hh"
62 #include "G4PhysicsModelCatalog.hh"
63 #include "G4Gamma.hh"
64 
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66 
69 
71  const G4String& pname)
72  : lowestKinEnergy(keV), verbose(1), name(modname), namePIXE(pname),
73  nameElectronPIXE(""), isActive(false), flagAuger(false),
74  flagPIXE(false), ignoreCuts(false)
75 {
76  vdyn.reserve(5);
77  theCoupleTable = 0;
78  G4String gg = "gammaPIXE";
79  G4String ee = "e-PIXE";
84 }
85 
86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
87 
89 {}
90 
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92 
94 {
95  // Define list of couples
97  G4int numOfCouples = theCoupleTable->GetTableSize();
98 
99  // needed for unit tests
100  G4int nn = std::max(numOfCouples, 1);
101  activeDeexcitationMedia.resize(nn, false);
102  activeAugerMedia.resize(nn, false);
103  activePIXEMedia.resize(nn, false);
104  activeZ.resize(93, false);
105 
106  // check if deexcitation is active for the given run
107  if( !isActive ) { return; }
108 
109  // Define list of regions
110  size_t nRegions = deRegions.size();
111 
112  if(0 == nRegions) {
114  nRegions = deRegions.size();
115  }
116 
117  if(0 < verbose) {
118  G4cout << G4endl;
119  G4cout << "### === Deexcitation model " << name
120  << " is activated for " << nRegions;
121  if(1 == nRegions) { G4cout << " region:" << G4endl; }
122  else { G4cout << " regions:" << G4endl;}
123  }
124 
125  // Identify active media
126  G4RegionStore* regionStore = G4RegionStore::GetInstance();
127  for(size_t j=0; j<nRegions; ++j) {
128  const G4Region* reg = regionStore->GetRegion(activeRegions[j], false);
129  if(reg && 0 < numOfCouples) {
130  const G4ProductionCuts* rpcuts = reg->GetProductionCuts();
131  if(0 < verbose) {
132  G4cout << " " << activeRegions[j] << G4endl;
133  }
134  for(G4int i=0; i<numOfCouples; ++i) {
135  const G4MaterialCutsCouple* couple =
137  if (couple->GetProductionCuts() == rpcuts) {
141  }
142  }
143  }
144  }
146  //G4cout << nelm << G4endl;
147  for(G4int k=0; k<nelm; ++k) {
148  G4int Z = G4lrint((*(G4Element::GetElementTable()))[k]->GetZ());
149  if(Z > 5 && Z < 93) {
150  activeZ[Z] = true;
151  //G4cout << "!!! Active de-excitation Z= " << Z << G4endl;
152  }
153  }
154 
155  // Initialise derived class
157 
158  if(0 < verbose && flagPIXE) {
159  G4cout << "### === PIXE model for hadrons: " << namePIXE
160  << " " << IsPIXEActive()
161  << G4endl;
162  G4cout << "### === PIXE model for e+-: " << nameElectronPIXE
163  << " " << IsPIXEActive()
164  << G4endl;
165  }
166 }
167 
168 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
169 
170 void
172  G4bool valDeexcitation,
173  G4bool valAuger,
174  G4bool valPIXE)
175 {
176  // no PIXE in parallel world
177  if(rname == "DefaultRegionForParallelWorld") { return; }
178 
179  G4String ss = rname;
180  //G4cout << "### G4VAtomDeexcitation::SetDeexcitationActiveRegion " << ss
181  // << " " << valDeexcitation << " " << valAuger
182  // << " " << valPIXE << G4endl;
183  if(ss == "world" || ss == "World" || ss == "WORLD") {
184  ss = "DefaultRegionForTheWorld";
185  }
186  size_t n = deRegions.size();
187  if(n > 0) {
188  for(size_t i=0; i<n; ++i) {
189 
190  // Region already exist
191  if(ss == activeRegions[i]) {
192  deRegions[i] = valDeexcitation;
193  AugerRegions[i] = valAuger;
194  PIXERegions[i] = valPIXE;
195  return;
196  }
197  }
198  }
199  // New region
200  activeRegions.push_back(ss);
201  deRegions.push_back(valDeexcitation);
202  AugerRegions.push_back(valAuger);
203  PIXERegions.push_back(valPIXE);
204 
205  // if de-excitation defined fo rthe world volume
206  // it should be active everywhere
207  if(ss == "DefaultRegionForTheWorld") {
209  G4int nn = regions->size();
210  for(G4int i=0; i<nn; ++i) {
211  SetDeexcitationActiveRegion((*regions)[i]->GetName(), valDeexcitation,
212  valAuger, valPIXE);
213  }
214  }
215 }
216 
217 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
218 
219 void
220 G4VAtomDeexcitation::AlongStepDeexcitation(std::vector<G4Track*>& tracks,
221  const G4Step& step,
222  G4double& eLossMax,
223  G4int coupleIndex)
224 {
225  G4double truelength = step.GetStepLength();
226  if(!flagPIXE && !activePIXEMedia[coupleIndex]) { return; }
227  if(eLossMax <= 0.0 || truelength <= 0.0) { return; }
228 
229  // step parameters
230  const G4StepPoint* preStep = step.GetPreStepPoint();
231  G4ThreeVector prePos = preStep->GetPosition();
232  G4ThreeVector delta = step.GetPostStepPoint()->GetPosition() - prePos;
233  G4double preTime = preStep->GetGlobalTime();
234  G4double dt = step.GetPostStepPoint()->GetGlobalTime() - preTime;
235 
236  // particle parameters
237  const G4Track* track = step.GetTrack();
238  const G4ParticleDefinition* part = track->GetDefinition();
239  G4double ekin = preStep->GetKineticEnergy();
240 
241  // media parameters
242  G4double gCut = (*theCoupleTable->GetEnergyCutsVector(0))[coupleIndex];
243  if(ignoreCuts) { gCut = 0.0; }
244  G4double eCut = DBL_MAX;
245  if(CheckAugerActiveRegion(coupleIndex)) {
246  eCut = (*theCoupleTable->GetEnergyCutsVector(1))[coupleIndex];
247  if(ignoreCuts) { eCut = 0.0; }
248  }
249 
250  //G4cout<<"!Sample PIXE gCut(MeV)= "<<gCut<<" eCut(MeV)= "<<eCut
251  // <<" Ekin(MeV)= " << ekin/MeV << G4endl;
252 
253  const G4Material* material = preStep->GetMaterial();
254  const G4ElementVector* theElementVector = material->GetElementVector();
255  const G4double* theAtomNumDensityVector =
256  material->GetVecNbOfAtomsPerVolume();
257  G4int nelm = material->GetNumberOfElements();
258 
259  // loop over deexcitations
260  for(G4int i=0; i<nelm; ++i) {
261  G4int Z = G4lrint((*theElementVector)[i]->GetZ());
262  if(activeZ[Z] && Z < 93) {
263  G4int nshells =
264  std::min(9,(*theElementVector)[i]->GetNbOfAtomicShells());
265  G4double rho = truelength*theAtomNumDensityVector[i];
266  //G4cout<<" Z "<< Z <<" is active x(mm)= " << truelength/mm << G4endl;
267  for(G4int ii=0; ii<nshells; ++ii) {
269  const G4AtomicShell* shell = GetAtomicShell(Z, as);
271 
272  if(gCut > bindingEnergy) { break; }
273 
274  if(eLossMax > bindingEnergy) {
275  G4double sig = rho*
276  GetShellIonisationCrossSectionPerAtom(part, Z, as, ekin, material);
277 
278  // mfp is mean free path in units of step size
279  if(sig > 0.0) {
280  G4double mfp = 1.0/sig;
281  G4double stot = 0.0;
282  //G4cout << " Shell " << ii << " mfp(mm)= " << mfp/mm << G4endl;
283  // sample ionisation points
284  do {
285  stot -= mfp*std::log(G4UniformRand());
286  if( stot > 1.0 || eLossMax < bindingEnergy) { break; }
287  // sample deexcitation
288  vdyn.clear();
289  GenerateParticles(&vdyn, shell, Z, gCut, eCut);
290  G4int nsec = vdyn.size();
291  if(nsec > 0) {
292  G4ThreeVector r = prePos + stot*delta;
293  G4double time = preTime + stot*dt;
294  for(G4int j=0; j<nsec; ++j) {
295  G4DynamicParticle* dp = vdyn[j];
296  G4double e = dp->GetKineticEnergy();
297 
298  // save new secondary if there is enough energy
299  if(eLossMax >= e) {
300  eLossMax -= e;
301  G4Track* t = new G4Track(dp, time, r);
302 
303  // defined secondary type
304  if(dp->GetDefinition() == gamma) {
306  } else {
308  }
309 
310  tracks.push_back(t);
311  } else {
312  delete dp;
313  }
314  }
315  }
316  } while (stot < 1.0);
317  }
318  }
319  }
320  }
321  }
322  return;
323 }
324 
325 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
G4ParticleDefinition * GetDefinition() const
G4ProductionCuts * GetProductionCuts() const
G4Region * GetRegion(const G4String &name, G4bool verbose=true) const
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
void AlongStepDeexcitation(std::vector< G4Track * > &tracks, const G4Step &step, G4double &eLoss, G4int coupleIndex)
std::vector< G4String > activeRegions
std::vector< G4Element * > G4ElementVector
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double GetStepLength() const
G4bool IsPIXEActive() const
std::vector< G4DynamicParticle * > vdyn
G4String name
Definition: TRTMaterials.hh:40
G4Material * GetMaterial() const
virtual G4double GetShellIonisationCrossSectionPerAtom(const G4ParticleDefinition *, G4int Z, G4AtomicShellEnumerator shell, G4double kinE, const G4Material *mat=0)=0
void SetDeexcitationActiveRegion(const G4String &rname, G4bool valDeexcitation, G4bool valAuger, G4bool valPIXE)
G4ParticleDefinition * GetDefinition() const
G4double BindingEnergy() const
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
std::vector< G4bool > activeAugerMedia
std::vector< G4bool > AugerRegions
void SetCreatorModelIndex(G4int idx)
static const G4double reg
static G4RegionStore * GetInstance()
G4StepPoint * GetPreStepPoint() const
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:206
#define G4UniformRand()
Definition: Randomize.hh:93
G4GLOB_DLL std::ostream G4cout
static size_t GetNumberOfElements()
Definition: G4Element.cc:410
const G4ThreeVector & GetPosition() const
const G4ParticleDefinition * gamma
bool G4bool
Definition: G4Types.hh:79
G4VAtomDeexcitation(const G4String &modname="Deexcitation", const G4String &pixename="")
Definition: G4Step.hh:76
G4EmParameters * theParameters
const G4int n
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4ProductionCutsTable * theCoupleTable
static G4ProductionCutsTable * GetProductionCutsTable()
G4bool CheckAugerActiveRegion(G4int coupleIndex)
std::vector< G4bool > activeDeexcitationMedia
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
int G4lrint(double ad)
Definition: templates.hh:163
T max(const T t1, const T t2)
brief Return the largest of the two arguments
static G4int Register(const G4String &)
const G4String & GetName() const
G4StepPoint * GetPostStepPoint() const
static G4EmParameters * Instance()
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
virtual void InitialiseForNewRun()=0
std::vector< G4bool > deRegions
#define G4endl
Definition: G4ios.hh:61
G4double GetGlobalTime() const
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
std::vector< G4bool > activeZ
static const double keV
Definition: G4SIunits.hh:195
G4double GetKineticEnergy() const
double G4double
Definition: G4Types.hh:76
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
G4Track * GetTrack() const
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:403
G4double bindingEnergy(G4int A, G4int Z)
G4ProductionCuts * GetProductionCuts() const
#define DBL_MAX
Definition: templates.hh:83
std::vector< G4bool > PIXERegions
G4AtomicShellEnumerator
std::vector< G4bool > activePIXEMedia