Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4PAIPhotModel.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: G4PAIPhotModel.cc 73607 2013-09-02 10:04:03Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class
31 // File name: G4PAIPhotModel.cc
32 //
33 // Author: Vladimir.Grichine@cern.ch on base of G4PAIModel MT interface
34 //
35 // Creation date: 07.10.2013
36 //
37 // Modifications:
38 //
39 //
40 
41 #include "G4PAIPhotModel.hh"
42 
43 #include "G4SystemOfUnits.hh"
44 #include "G4PhysicalConstants.hh"
45 #include "G4Region.hh"
46 #include "G4PhysicsLogVector.hh"
47 #include "G4PhysicsFreeVector.hh"
48 #include "G4PhysicsTable.hh"
49 #include "G4ProductionCutsTable.hh"
50 #include "G4MaterialCutsCouple.hh"
51 #include "G4MaterialTable.hh"
52 #include "G4SandiaTable.hh"
53 #include "G4OrderedTable.hh"
54 #include "G4RegionStore.hh"
55 
56 #include "Randomize.hh"
57 #include "G4Electron.hh"
58 #include "G4Positron.hh"
59 #include "G4Gamma.hh"
60 #include "G4Poisson.hh"
61 #include "G4Step.hh"
62 #include "G4Material.hh"
63 #include "G4DynamicParticle.hh"
64 #include "G4ParticleDefinition.hh"
66 #include "G4PAIPhotData.hh"
67 #include "G4DeltaAngle.hh"
68 
70 
71 using namespace std;
72 
75  fVerbose(0),
76  fModelData(nullptr),
77  fParticle(nullptr)
78 {
81 
82  fParticleChange = nullptr;
83 
84  if(p) { SetParticle(p); }
85  else { SetParticle(fElectron); }
86 
87  // default generator
89  fLowestTcut = 12.5*CLHEP::eV;
90 }
91 
93 
95 {
96  //G4cout << "G4PAIPhotModel::~G4PAIPhotModel() " << this << G4endl;
97  if(IsMaster()) { delete fModelData; fModelData = nullptr; }
98 }
99 
101 
103  const G4DataVector& cuts)
104 {
105  if(fVerbose > 0)
106  {
107  G4cout<<"G4PAIPhotModel::Initialise for "<<p->GetParticleName()<<G4endl;
108  }
109  SetParticle(p);
111 
112  if( IsMaster() )
113  {
114 
116 
117  delete fModelData;
118  fMaterialCutsCoupleVector.clear();
119 
120  G4double tmin = LowEnergyLimit()*fRatio;
122  fModelData = new G4PAIPhotData(tmin, tmax, fVerbose);
123 
124  // Prepare initialization
125  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
126  size_t numOfMat = G4Material::GetNumberOfMaterials();
127  size_t numRegions = fPAIRegionVector.size();
128 
129  // protect for unit tests
130  if(0 == numRegions) {
131  G4Exception("G4PAIModel::Initialise()","em0106",JustWarning,
132  "no G4Regions are registered for the PAI model - World is used");
134  ->GetRegion("DefaultRegionForTheWorld", false));
135  numRegions = 1;
136  }
137 
138  for( size_t iReg = 0; iReg < numRegions; ++iReg )
139  {
140  const G4Region* curReg = fPAIRegionVector[iReg];
141  G4Region* reg = const_cast<G4Region*>(curReg);
142 
143  for(size_t jMat = 0; jMat < numOfMat; ++jMat)
144  {
145  G4Material* mat = (*theMaterialTable)[jMat];
146  const G4MaterialCutsCouple* cutCouple = reg->FindCouple(mat);
147  //G4cout << "Couple <" << fCutCouple << G4endl;
148  if(cutCouple)
149  {
150  if(fVerbose>0)
151  {
152  G4cout << "Reg <" <<curReg->GetName() << "> mat <"
153  << mat->GetName() << "> fCouple= "
154  << cutCouple << ", idx= " << cutCouple->GetIndex()
155  <<" " << p->GetParticleName()
156  <<", cuts.size() = " << cuts.size() << G4endl;
157  }
158  // check if this couple is not already initialized
159 
160  size_t n = fMaterialCutsCoupleVector.size();
161 
162  G4bool isnew = true;
163  if( 0 < n )
164  {
165  for(size_t i=0; i<fMaterialCutsCoupleVector.size(); ++i)
166  {
167  if(cutCouple == fMaterialCutsCoupleVector[i]) {
168  isnew = false;
169  break;
170  }
171  }
172  }
173  // initialise data banks
174  if(isnew) {
175  fMaterialCutsCoupleVector.push_back(cutCouple);
176  G4double deltaCutInKinEnergy = cuts[cutCouple->GetIndex()];
177  fModelData->Initialise(cutCouple, deltaCutInKinEnergy, this);
178  }
179  }
180  }
181  }
182  }
183 }
184 
186 
188  G4VEmModel* masterModel)
189 {
190  fModelData = static_cast<G4PAIPhotModel*>(masterModel)->GetPAIPhotData();
191  fMaterialCutsCoupleVector = static_cast<G4PAIPhotModel*>(masterModel)->GetVectorOfCouples();
192  SetElementSelectors( masterModel->GetElementSelectors() );
193 }
194 
196 
198  const G4MaterialCutsCouple*)
199 {
200  return fLowestTcut;
201 }
202 
204 
206  const G4ParticleDefinition* p,
208  G4double cutEnergy)
209 {
210  G4int coupleIndex = FindCoupleIndex(CurrentCouple());
211  if(0 > coupleIndex) { return 0.0; }
212 
213  G4double cut = std::min(MaxSecondaryEnergy(p, kineticEnergy), cutEnergy);
214 
215  G4double scaledTkin = kineticEnergy*fRatio;
216 
217  return fChargeSquare*fModelData->DEDXPerVolume(coupleIndex, scaledTkin, cut);
218 }
219 
221 
223  const G4ParticleDefinition* p,
225  G4double cutEnergy,
226  G4double maxEnergy )
227 {
228  G4int coupleIndex = FindCoupleIndex(CurrentCouple());
229 
230  if(0 > coupleIndex) return 0.0;
231 
232  G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy);
233 
234  if(tmax <= cutEnergy) return 0.0;
235 
236  G4double scaledTkin = kineticEnergy*fRatio;
238  scaledTkin,
239  cutEnergy, tmax);
240  return xsc;
241 }
242 
244 //
245 // It is analog of PostStepDoIt in terms of secondary electron.
246 //
247 
248 void G4PAIPhotModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
249  const G4MaterialCutsCouple* matCC,
250  const G4DynamicParticle* dp,
251  G4double tmin,
252  G4double maxEnergy)
253 {
254  G4int coupleIndex = FindCoupleIndex(matCC);
255  if(0 > coupleIndex) { return; }
256 
257  SetParticle(dp->GetDefinition());
258 
260 
261  G4double tmax = MaxSecondaryEnergy(fParticle, kineticEnergy);
262 
263  if( maxEnergy < tmax) tmax = maxEnergy;
264  if( tmin >= tmax) return;
265 
266  G4ThreeVector direction = dp->GetMomentumDirection();
267  G4double scaledTkin = kineticEnergy*fRatio;
268  G4double totalEnergy = kineticEnergy + fMass;
269  G4double totalMomentum = sqrt(kineticEnergy*(totalEnergy + fMass));
270  G4double plRatio = fModelData->GetPlasmonRatio(coupleIndex, scaledTkin);
271 
272  if( G4UniformRand() <= plRatio )
273  {
274  G4double deltaTkin = fModelData->SamplePostStepPlasmonTransfer(coupleIndex, scaledTkin);
275 
276  // G4cout<<"G4PAIPhotModel::SampleSecondaries; dp "<<dp->GetParticleDefinition()->GetParticleName()
277  // <<"; Tkin = "<<kineticEnergy/keV<<" keV; transfer = "<<deltaTkin/keV<<" keV "<<G4endl;
278 
279  if( deltaTkin <= 0. && fVerbose > 0)
280  {
281  G4cout<<"G4PAIPhotModel::SampleSecondary e- deltaTkin = "<<deltaTkin<<G4endl;
282  }
283  if( deltaTkin <= 0.) return;
284 
285  if( deltaTkin > tmax) deltaTkin = tmax;
286 
287  const G4Element* anElement = SelectRandomAtom(matCC,fParticle,kineticEnergy);
288  G4int Z = G4lrint(anElement->GetZ());
289 
291  GetAngularDistribution()->SampleDirection(dp, deltaTkin,
292  Z, matCC->GetMaterial()),
293  deltaTkin);
294 
295  // primary change
296 
297  kineticEnergy -= deltaTkin;
298 
299  if( kineticEnergy <= 0. ) // kill primary as local? energy deposition
300  {
302  fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy+deltaTkin);
303  // fParticleChange->ProposeTrackStatus(fStopAndKill);
304  return;
305  }
306  else
307  {
308  G4ThreeVector dir = totalMomentum*direction - deltaRay->GetMomentum();
309  direction = dir.unit();
312  vdp->push_back(deltaRay);
313  }
314  }
315  else // secondary X-ray CR photon
316  {
317  G4double deltaTkin = fModelData->SamplePostStepPhotonTransfer(coupleIndex, scaledTkin);
318 
319  // G4cout<<"PAIPhotonModel PhotonPostStepTransfer = "<<deltaTkin/keV<<" keV"<<G4endl;
320 
321  if( deltaTkin <= 0. )
322  {
323  G4cout<<"G4PAIPhotonModel::SampleSecondary gamma deltaTkin = "<<deltaTkin<<G4endl;
324  }
325  if( deltaTkin <= 0.) return;
326 
327  if( deltaTkin >= kineticEnergy ) // stop primary
328  {
329  deltaTkin = kineticEnergy;
330  kineticEnergy = 0.0;
331  }
332  G4double costheta = 0.; // G4UniformRand(); // VG: ??? for start only
333  G4double sintheta = sqrt((1.+costheta)*(1.-costheta));
334 
335  // direction of the 'Cherenkov' photon
336  G4double phi = twopi*G4UniformRand();
337  G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta;
338 
339  G4ThreeVector deltaDirection(dirx,diry,dirz);
340  deltaDirection.rotateUz(direction);
341 
342  if( kineticEnergy > 0.) // primary change
343  {
344  kineticEnergy -= deltaTkin;
346  }
347  else // stop primary, but pass X-ray CR
348  {
349  // fParticleChange->ProposeLocalEnergyDeposit(deltaTkin);
351  }
352  // create G4DynamicParticle object for photon ray
353 
354  G4DynamicParticle* photonRay = new G4DynamicParticle;
355  photonRay->SetDefinition( G4Gamma::Gamma() );
356  photonRay->SetKineticEnergy( deltaTkin );
357  photonRay->SetMomentumDirection(deltaDirection);
358 
359  vdp->push_back(photonRay);
360  }
361  return;
362 }
363 
365 
367  const G4DynamicParticle* aParticle,
368  G4double, G4double step,
369  G4double eloss)
370 {
371  // return 0.;
372  G4int coupleIndex = FindCoupleIndex(matCC);
373  if(0 > coupleIndex) { return eloss; }
374 
375  SetParticle(aParticle->GetDefinition());
376 
377 
378  // G4cout << "G4PAIPhotModel::SampleFluctuations step(mm)= "<< step/mm
379  // << " Eloss(keV)= " << eloss/keV << " in "
380  // << matCC->GetMaterial()->GetName() << G4endl;
381 
382 
383  G4double Tkin = aParticle->GetKineticEnergy();
384  G4double scaledTkin = Tkin*fRatio;
385 
386  G4double loss = fModelData->SampleAlongStepPhotonTransfer(coupleIndex, Tkin,
387  scaledTkin,
388  step*fChargeSquare);
389  loss += fModelData->SampleAlongStepPlasmonTransfer(coupleIndex, Tkin,
390  scaledTkin,
391  step*fChargeSquare);
392 
393 
394  // G4cout<<" PAIPhotModel::SampleFluctuations loss = "<<loss/keV<<" keV, on step = "
395  // <<step/mm<<" mm"<<G4endl;
396  return loss;
397 
398 }
399 
401 //
402 // Returns the statistical estimation of the energy loss distribution variance
403 //
404 
405 
407  const G4DynamicParticle* aParticle,
408  G4double tmax,
409  G4double step )
410 {
411  G4double particleMass = aParticle->GetMass();
412  G4double electronDensity = material->GetElectronDensity();
413  G4double kineticEnergy = aParticle->GetKineticEnergy();
414  G4double q = aParticle->GetCharge()/eplus;
415  G4double etot = kineticEnergy + particleMass;
416  G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*particleMass)/(etot*etot);
417  G4double siga = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * step
418  * electronDensity * q * q;
419 
420  return siga;
421  /*
422  G4double loss, sumLoss=0., sumLoss2=0., sigma2, meanLoss=0.;
423  for(G4int i = 0; i < fMeanNumber; i++)
424  {
425  loss = SampleFluctuations(material,aParticle,tmax,step,meanLoss);
426  sumLoss += loss;
427  sumLoss2 += loss*loss;
428  }
429  meanLoss = sumLoss/fMeanNumber;
430  sigma2 = meanLoss*meanLoss + (sumLoss2-2*sumLoss*meanLoss)/fMeanNumber;
431  return sigma2;
432  */
433 }
434 
436 
438  G4double kinEnergy)
439 {
440  SetParticle(p);
441  G4double tmax = kinEnergy;
442  if(p == fElectron) { tmax *= 0.5; }
443  else if(p != fPositron) {
445  G4double gamma= kinEnergy/fMass + 1.0;
446  tmax = 2.0*electron_mass_c2*(gamma*gamma - 1.) /
447  (1. + 2.0*gamma*ratio + ratio*ratio);
448  }
449  return tmax;
450 }
451 
453 
455 {
456  fPAIRegionVector.push_back(r);
457 }
458 
459 //
460 //
const G4ParticleDefinition * fParticle
G4double SamplePostStepPlasmonTransfer(G4int coupleIndex, G4double scaledTkin) const
void SetElementSelectors(std::vector< G4EmElementSelector * > *)
Definition: G4VEmModel.hh:792
G4PAIPhotData * GetPAIPhotData()
static size_t GetNumberOfMaterials()
Definition: G4Material.cc:600
G4double fChargeSquare
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:593
void SetKineticEnergy(G4double aEnergy)
G4double Dispersion(const G4Material *, const G4DynamicParticle *, G4double, G4double) final
G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy) final
const std::vector< const G4MaterialCutsCouple * > & GetVectorOfCouples()
static constexpr double twopi_mc2_rcl2
const G4String & GetName() const
#define G4endl
Definition: G4ios.hh:61
const G4ThreeVector & GetMomentumDirection() const
const char * p
Definition: xmltok.h:285
const G4String & GetParticleName() const
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:552
void SetProposedMomentumDirection(const G4ThreeVector &dir)
std::vector< G4EmElementSelector * > * GetElementSelectors()
Definition: G4VEmModel.hh:784
G4int FindCoupleIndex(const G4MaterialCutsCouple *)
G4double GetPlasmonRatio(G4int coupleIndex, G4double scaledTkin) const
G4double fLowestTcut
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:616
void SetMomentumDirection(const G4ThreeVector &aDirection)
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy) final
const G4String & GetName() const
Definition: G4Material.hh:179
Float_t Z
G4double SampleFluctuations(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double, G4double, G4double) final
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:116
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
G4ParticleDefinition * GetDefinition() const
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4VEmModel.cc:146
G4ParticleChangeForLoss * fParticleChange
G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kinEnergy) final
static constexpr double electron_mass_c2
const G4ParticleDefinition * fElectron
void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) final
#define G4UniformRand()
Definition: Randomize.hh:53
static constexpr double twopi
Definition: G4SIunits.hh:76
G4double GetCharge() const
static G4Positron * Positron()
Definition: G4Positron.cc:94
static constexpr double eV
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:582
Float_t mat
static G4Electron * Electron()
Definition: G4Electron.cc:94
void Initialise(const G4MaterialCutsCouple *, G4double cut, G4PAIPhotModel *)
std::vector< const G4MaterialCutsCouple * > fMaterialCutsCoupleVector
const G4MaterialCutsCouple * CurrentCouple() const
Definition: G4VEmModel.hh:462
Hep3Vector unit() const
std::vector< G4Material * > G4MaterialTable
G4PAIPhotData * fModelData
G4double MinEnergyCut(const G4ParticleDefinition *, const G4MaterialCutsCouple *couple) final
static constexpr double eplus
Definition: G4SIunits.hh:199
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
G4double SamplePostStepPhotonTransfer(G4int coupleIndex, G4double scaledTkin) const
G4double SampleAlongStepPlasmonTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
G4PAIPhotModel(const G4ParticleDefinition *p=nullptr, const G4String &nam="PAI")
int G4lrint(double ad)
Definition: templates.hh:151
const G4ParticleDefinition * fPositron
static G4RegionStore * GetInstance()
int G4int
Definition: G4Types.hh:78
void DefineForRegion(const G4Region *r) final
void SetProposedKineticEnergy(G4double proposedKinEnergy)
TDirectory * dir
G4double GetKineticEnergy() const
G4double SampleAlongStepPhotonTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
G4double CrossSectionPerVolume(G4int coupleIndex, G4double scaledTkin, G4double tcut, G4double tmax) const
G4GLOB_DLL std::ostream G4cout
~G4PAIPhotModel() final
G4bool IsMaster() const
Definition: G4VEmModel.hh:700
void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel) final
G4MaterialCutsCouple * FindCouple(G4Material *mat)
void Initialise(const G4ParticleDefinition *, const G4DataVector &) final
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:575
G4double GetZ() const
Definition: G4Element.hh:131
Char_t n[5]
const G4Material * GetMaterial() const
G4double GetMass() const
static const G4double reg
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetParticle(const G4ParticleDefinition *p)
std::vector< const G4Region * > fPAIRegionVector
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:609
G4double GetElectronDensity() const
Definition: G4Material.hh:218
G4double DEDXPerVolume(G4int coupleIndex, G4double scaledTkin, G4double cut) const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments