Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4VAtomDeexcitation.cc
이 파일의 문서화 페이지로 가기
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 108386 2018-02-09 15:38:32Z 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 #ifdef G4MULTITHREADED
66  G4Mutex G4VAtomDeexcitation::atomDeexcitationMutex = G4MUTEX_INITIALIZER;
67 #endif
68 
69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
70 
73 
75  : verbose(1), name(modname), isActive(false), flagAuger(false),
76  flagAugerCascade(false), flagPIXE(false), ignoreCuts(false),
77  isActiveLocked(false), isAugerLocked(false),
78  isAugerCascadeLocked(false), isPIXELocked(false)
79 {
81  vdyn.reserve(5);
82  theCoupleTable = nullptr;
83  G4String gg = "gammaPIXE";
84  G4String ee = "e-PIXE";
85  if(pixeIDg < 0) {
86 #ifdef G4MULTITHREADED
87  G4MUTEXLOCK(&atomDeexcitationMutex);
88  if(pixeIDg < 0) {
89 #endif
92 #ifdef G4MULTITHREADED
93  }
94  G4MUTEXUNLOCK(&atomDeexcitationMutex);
95 #endif
96  }
98 }
99 
100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101 
103 {}
104 
105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
106 
108 {
110 
111  // Define list of couples
113  G4int numOfCouples = theCoupleTable->GetTableSize();
114 
115  // needed for unit tests
116  G4int nn = std::max(numOfCouples, 1);
117  activeDeexcitationMedia.resize(nn, false);
118  activeAugerMedia.resize(nn, false);
119  activePIXEMedia.resize(nn, false);
120  activeZ.resize(93, false);
121 
122  // initialisation of flags and options
123  // normally there is no locksed flags
127  if(!isPIXELocked) { flagPIXE = theParameters->Pixe(); }
129 
130  // Define list of regions
131  size_t nRegions = deRegions.size();
132  // check if deexcitation is active for the given run
133  if(!isActive && 0 == nRegions) { return; }
134 
135  // if no active regions add a world
136  if(0 == nRegions) {
138  nRegions = deRegions.size();
139  }
140 
141  if(0 < verbose) {
142  G4cout << G4endl;
143  G4cout << "### === Deexcitation model " << name
144  << " is activated for " << nRegions;
145  if(1 == nRegions) { G4cout << " region:" << G4endl; }
146  else { G4cout << " regions:" << G4endl;}
147  }
148 
149  // Identify active media
150  G4RegionStore* regionStore = G4RegionStore::GetInstance();
151  for(size_t j=0; j<nRegions; ++j) {
152  const G4Region* reg = regionStore->GetRegion(activeRegions[j], false);
153  if(reg && 0 < numOfCouples) {
154  const G4ProductionCuts* rpcuts = reg->GetProductionCuts();
155  if(0 < verbose) {
156  G4cout << " " << activeRegions[j]
157  << " " << deRegions[j] << " " << AugerRegions[j]
158  << " " << PIXERegions[j] << G4endl;
159  }
160  for(G4int i=0; i<numOfCouples; ++i) {
161  const G4MaterialCutsCouple* couple =
163  if (couple->GetProductionCuts() == rpcuts) {
167  }
168  }
169  }
170  }
172  //G4cout << nelm << G4endl;
173  for(G4int k=0; k<nelm; ++k) {
174  G4int Z = (*(G4Element::GetElementTable()))[k]->GetZasInt();
175  if(Z > 5 && Z < 93) {
176  activeZ[Z] = true;
177  //G4cout << "!!! Active de-excitation Z= " << Z << G4endl;
178  }
179  }
180 
181  // Initialise derived class
183 
184  if(0 < verbose && flagAuger) {
185  G4cout << "### === Auger cascade flag: " << flagAugerCascade
186  << G4endl;
187  }
188  if(0 < verbose) {
189  G4cout << "### === Ignore cuts flag: " << ignoreCuts
190  << G4endl;
191  }
192  if(0 < verbose && flagPIXE) {
193  G4cout << "### === PIXE model for hadrons: "
195  << G4endl;
196  G4cout << "### === PIXE model for e+-: "
198  << G4endl;
199  }
200 }
201 
202 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
203 
204 void
206  G4bool valDeexcitation,
207  G4bool valAuger,
208  G4bool valPIXE)
209 {
210  // no PIXE in parallel world
211  if(rname == "DefaultRegionForParallelWorld") { return; }
212 
213  G4String ss = rname;
214  /*
215  G4cout << "### G4VAtomDeexcitation::SetDeexcitationActiveRegion " << ss
216  << " " << valDeexcitation << " " << valAuger
217  << " " << valPIXE << G4endl;
218  */
219  if(ss == "world" || ss == "World" || ss == "WORLD") {
220  ss = "DefaultRegionForTheWorld";
221  }
222  size_t n = deRegions.size();
223  for(size_t i=0; i<n; ++i) {
224 
225  // Region already exist
226  if(ss == activeRegions[i]) {
227  deRegions[i] = valDeexcitation;
228  AugerRegions[i] = valAuger;
229  PIXERegions[i] = valPIXE;
230  return;
231  }
232  }
233  // New region
234  activeRegions.push_back(ss);
235  deRegions.push_back(valDeexcitation);
236  AugerRegions.push_back(valAuger);
237  PIXERegions.push_back(valPIXE);
238 
239  // if de-excitation defined for the world volume
240  // it should be active for all G4Regions
241  if(ss == "DefaultRegionForTheWorld") {
243  G4int nn = regions->size();
244  for(G4int i=0; i<nn; ++i) {
245  if(ss == (*regions)[i]->GetName()) { continue; }
246  SetDeexcitationActiveRegion((*regions)[i]->GetName(), valDeexcitation,
247  valAuger, valPIXE);
248 
249  }
250  }
251 }
252 
253 void G4VAtomDeexcitation::GenerateParticles(std::vector<G4DynamicParticle*>* v,
254  const G4AtomicShell* as,
255  G4int Z, G4int idx)
256 {
257  G4double gCut = DBL_MAX;
258  if(ignoreCuts) {
259  gCut = 0.0;
260  } else if (theCoupleTable) {
261  gCut = (*(theCoupleTable->GetEnergyCutsVector(0)))[idx];
262  }
263  if(gCut < as->BindingEnergy()) {
264  G4double eCut = DBL_MAX;
265  if(CheckAugerActiveRegion(idx)) {
266  if(ignoreCuts) {
267  eCut = 0.0;
268  } else if (theCoupleTable) {
269  eCut = (*(theCoupleTable->GetEnergyCutsVector(1)))[idx];
270  }
271  }
272  GenerateParticles(v, as, Z, gCut, eCut);
273  }
274 }
275 
276 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
277 
278 void
279 G4VAtomDeexcitation::AlongStepDeexcitation(std::vector<G4Track*>& tracks,
280  const G4Step& step,
281  G4double& eLossMax,
282  G4int coupleIndex)
283 {
284  G4double truelength = step.GetStepLength();
285  if(!flagPIXE && !activePIXEMedia[coupleIndex]) { return; }
286  if(eLossMax <= 0.0 || truelength <= 0.0) { return; }
287 
288  // step parameters
289  const G4StepPoint* preStep = step.GetPreStepPoint();
290  G4ThreeVector prePos = preStep->GetPosition();
291  G4ThreeVector delta = step.GetPostStepPoint()->GetPosition() - prePos;
292  G4double preTime = preStep->GetGlobalTime();
293  G4double dt = step.GetPostStepPoint()->GetGlobalTime() - preTime;
294 
295  // particle parameters
296  const G4Track* track = step.GetTrack();
297  const G4ParticleDefinition* part = track->GetDefinition();
298  G4double ekin = preStep->GetKineticEnergy();
299 
300  // media parameters
301  G4double gCut = (*theCoupleTable->GetEnergyCutsVector(0))[coupleIndex];
302  if(ignoreCuts) { gCut = 0.0; }
303  G4double eCut = DBL_MAX;
304  if(CheckAugerActiveRegion(coupleIndex)) {
305  eCut = (*theCoupleTable->GetEnergyCutsVector(1))[coupleIndex];
306  if(ignoreCuts) { eCut = 0.0; }
307  }
308 
309  //G4cout<<"!Sample PIXE gCut(MeV)= "<<gCut<<" eCut(MeV)= "<<eCut
310  // <<" Ekin(MeV)= " << ekin/MeV << G4endl;
311 
312  const G4Material* material = preStep->GetMaterial();
313  const G4ElementVector* theElementVector = material->GetElementVector();
314  const G4double* theAtomNumDensityVector =
315  material->GetVecNbOfAtomsPerVolume();
316  G4int nelm = material->GetNumberOfElements();
317 
318  // loop over deexcitations
319  for(G4int i=0; i<nelm; ++i) {
320  G4int Z = (*theElementVector)[i]->GetZasInt();
321  if(activeZ[Z] && Z < 93) {
322  G4int nshells =
323  std::min(9,(*theElementVector)[i]->GetNbOfAtomicShells());
324  G4double rho = truelength*theAtomNumDensityVector[i];
325  //G4cout<<" Z "<< Z <<" is active x(mm)= " << truelength/mm << G4endl;
326  for(G4int ii=0; ii<nshells; ++ii) {
328  const G4AtomicShell* shell = GetAtomicShell(Z, as);
330 
331  if(gCut > bindingEnergy) { break; }
332 
333  if(eLossMax > bindingEnergy) {
334  G4double sig = rho*
335  GetShellIonisationCrossSectionPerAtom(part, Z, as, ekin, material);
336 
337  // mfp is mean free path in units of step size
338  if(sig > 0.0) {
339  G4double mfp = 1.0/sig;
340  G4double stot = 0.0;
341  //G4cout << " Shell " << ii << " mfp(mm)= " << mfp/mm << G4endl;
342  // sample ionisation points
343  do {
344  stot -= mfp*std::log(G4UniformRand());
345  if( stot > 1.0 || eLossMax < bindingEnergy) { break; }
346  // sample deexcitation
347  vdyn.clear();
348  GenerateParticles(&vdyn, shell, Z, gCut, eCut);
349  G4int nsec = vdyn.size();
350  if(nsec > 0) {
351  G4ThreeVector r = prePos + stot*delta;
352  G4double time = preTime + stot*dt;
353  for(G4int j=0; j<nsec; ++j) {
354  G4DynamicParticle* dp = vdyn[j];
355  G4double e = dp->GetKineticEnergy();
356 
357  // save new secondary if there is enough energy
358  if(eLossMax >= e) {
359  eLossMax -= e;
360  G4Track* t = new G4Track(dp, time, r);
361 
362  // defined secondary type
363  if(dp->GetDefinition() == gamma) {
365  } else {
367  }
368 
369  tracks.push_back(t);
370  } else {
371  delete dp;
372  }
373  }
374  }
375  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
376  } while (stot < 1.0);
377  }
378  }
379  }
380  }
381  }
382  return;
383 }
384 
385 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
T max(const T t1, const T t2)
brief Return the largest of the two arguments
const XML_Char * name
Definition: expat.h:151
std::vector< G4bool > AugerRegions
G4EmParameters * theParameters
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
std::vector< G4bool > deRegions
G4StepPoint * GetPreStepPoint() const
G4ProductionCutsTable * theCoupleTable
void AlongStepDeexcitation(std::vector< G4Track * > &tracks, const G4Step &step, G4double &eLoss, G4int coupleIndex)
const G4String & GetName() const
std::vector< G4bool > activeZ
void SetCreatorModelIndex(G4int idx)
void SetDeexcitationActiveRegion(const G4String &rname, G4bool valDeexcitation, G4bool valAuger, G4bool valPIXE)
#define G4endl
Definition: G4ios.hh:61
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:191
std::vector< G4String > activeRegions
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:207
G4AtomicShellEnumerator
G4ParticleDefinition * GetDefinition() const
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
G4Region * GetRegion(const G4String &name, G4bool verbose=true) const
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
Float_t Z
std::vector< G4bool > PIXERegions
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
TString part[npart]
Definition: Style.C:32
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:88
G4bool CheckAugerActiveRegion(G4int coupleIndex)
G4ParticleDefinition * GetDefinition() const
G4Material * GetMaterial() const
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
G4double GetStepLength() const
const G4ThreeVector & GetPosition() const
std::vector< G4DynamicParticle * > vdyn
std::vector< G4bool > activePIXEMedia
#define G4UniformRand()
Definition: Randomize.hh:53
std::vector< G4bool > activeAugerMedia
G4Track * GetTrack() const
G4double GetGlobalTime() const
const G4ParticleDefinition * gamma
Definition: G4Step.hh:76
G4StepPoint * GetPostStepPoint() const
G4double GetKineticEnergy() const
G4bool Pixe() const
#define G4MUTEXUNLOCK(mutex)
Definition: G4Threading.hh:234
G4ProductionCuts * GetProductionCuts() const
G4bool AugerCascade() const
const G4String & PIXECrossSectionModel()
static G4ProductionCutsTable * GetProductionCutsTable()
G4VAtomDeexcitation(const G4String &modname="Deexcitation")
static G4RegionStore * GetInstance()
std::vector< G4Element * > G4ElementVector
int G4int
Definition: G4Types.hh:78
G4bool DeexcitationIgnoreCut() const
const G4String & PIXEElectronCrossSectionModel()
#define G4MUTEXLOCK(mutex)
Definition: G4Threading.hh:233
G4double GetKineticEnergy() const
G4GLOB_DLL std::ostream G4cout
void DefineRegParamForDeex(G4VAtomDeexcitation *) const
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
G4ProductionCuts * GetProductionCuts() const
G4bool Auger() const
Char_t n[5]
static size_t GetNumberOfElements()
Definition: G4Element.cc:405
static const G4double reg
std::vector< G4bool > activeDeexcitationMedia
G4double bindingEnergy(G4int A, G4int Z)
#define DBL_MAX
Definition: templates.hh:83
virtual G4double GetShellIonisationCrossSectionPerAtom(const G4ParticleDefinition *, G4int Z, G4AtomicShellEnumerator shell, G4double kinE, const G4Material *mat=nullptr)=0
virtual void InitialiseForNewRun()=0
static G4EmParameters * Instance()
size_t GetNumberOfElements() const
Definition: G4Material.hh:187
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:398
G4double BindingEnergy() const
std::mutex G4Mutex
Definition: G4Threading.hh:84
static G4int Register(const G4String &)
G4bool Fluo() const