Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PAIPhotModel.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: 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 {
79  fElectron = G4Electron::Electron();
80  fPositron = G4Positron::Positron();
81 
82  fParticleChange = nullptr;
83 
84  if(p) { SetParticle(p); }
85  else { SetParticle(fElectron); }
86 
87  // default generator
89 }
90 
92 
94 {
95  //G4cout << "G4PAIPhotModel::~G4PAIPhotModel() " << this << G4endl;
96  if(IsMaster()) { delete fModelData; fModelData = nullptr; }
97 }
98 
100 
102  const G4DataVector& cuts)
103 {
104  if(fVerbose > 0)
105  {
106  G4cout<<"G4PAIPhotModel::Initialise for "<<p->GetParticleName()<<G4endl;
107  }
108  SetParticle(p);
109  fParticleChange = GetParticleChangeForLoss();
110 
111  if( IsMaster() )
112  {
113 
115 
116  delete fModelData;
117  fMaterialCutsCoupleVector.clear();
118 
119  G4double tmin = LowEnergyLimit()*fRatio;
120  G4double tmax = HighEnergyLimit()*fRatio;
121  fModelData = new G4PAIPhotData(tmin, tmax, fVerbose);
122 
123  // Prepare initialization
124  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
125  size_t numOfMat = G4Material::GetNumberOfMaterials();
126  size_t numRegions = fPAIRegionVector.size();
127 
128  // protect for unit tests
129  if(0 == numRegions) {
130  G4Exception("G4PAIModel::Initialise()","em0106",JustWarning,
131  "no G4Regions are registered for the PAI model - World is used");
132  fPAIRegionVector.push_back(G4RegionStore::GetInstance()
133  ->GetRegion("DefaultRegionForTheWorld", false));
134  numRegions = 1;
135  }
136 
137  for( size_t iReg = 0; iReg < numRegions; ++iReg )
138  {
139  const G4Region* curReg = fPAIRegionVector[iReg];
140  G4Region* reg = const_cast<G4Region*>(curReg);
141 
142  for(size_t jMat = 0; jMat < numOfMat; ++jMat)
143  {
144  G4Material* mat = (*theMaterialTable)[jMat];
145  const G4MaterialCutsCouple* cutCouple = reg->FindCouple(mat);
146  //G4cout << "Couple <" << fCutCouple << G4endl;
147  if(cutCouple)
148  {
149  if(fVerbose>0)
150  {
151  G4cout << "Reg <" <<curReg->GetName() << "> mat <"
152  << mat->GetName() << "> fCouple= "
153  << cutCouple << ", idx= " << cutCouple->GetIndex()
154  <<" " << p->GetParticleName()
155  <<", cuts.size() = " << cuts.size() << G4endl;
156  }
157  // check if this couple is not already initialized
158 
159  size_t n = fMaterialCutsCoupleVector.size();
160 
161  G4bool isnew = true;
162  if( 0 < n )
163  {
164  for(size_t i=0; i<fMaterialCutsCoupleVector.size(); ++i)
165  {
166  if(cutCouple == fMaterialCutsCoupleVector[i]) {
167  isnew = false;
168  break;
169  }
170  }
171  }
172  // initialise data banks
173  if(isnew) {
174  fMaterialCutsCoupleVector.push_back(cutCouple);
175  G4double deltaCutInKinEnergy = cuts[cutCouple->GetIndex()];
176  fModelData->Initialise(cutCouple, deltaCutInKinEnergy, this);
177  }
178  }
179  }
180  }
181  }
182 }
183 
185 
187  G4VEmModel* masterModel)
188 {
189  fModelData = static_cast<G4PAIPhotModel*>(masterModel)->GetPAIPhotData();
190  fMaterialCutsCoupleVector = static_cast<G4PAIPhotModel*>(masterModel)->GetVectorOfCouples();
191  SetElementSelectors( masterModel->GetElementSelectors() );
192 }
193 
195 
197  const G4ParticleDefinition* p,
198  G4double kineticEnergy,
199  G4double cutEnergy)
200 {
201  G4int coupleIndex = FindCoupleIndex(CurrentCouple());
202  if(0 > coupleIndex) { return 0.0; }
203 
204  G4double cut = std::min(MaxSecondaryEnergy(p, kineticEnergy), cutEnergy);
205 
206  G4double scaledTkin = kineticEnergy*fRatio;
207 
208  return fChargeSquare*fModelData->DEDXPerVolume(coupleIndex, scaledTkin, cut);
209 }
210 
212 
214  const G4ParticleDefinition* p,
215  G4double kineticEnergy,
216  G4double cutEnergy,
217  G4double maxEnergy )
218 {
219  G4int coupleIndex = FindCoupleIndex(CurrentCouple());
220 
221  if(0 > coupleIndex) return 0.0;
222 
223  G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy);
224 
225  if(tmax <= cutEnergy) return 0.0;
226 
227  G4double scaledTkin = kineticEnergy*fRatio;
228  G4double xsc = fChargeSquare*fModelData->CrossSectionPerVolume(coupleIndex,
229  scaledTkin,
230  cutEnergy, tmax);
231  return xsc;
232 }
233 
235 //
236 // It is analog of PostStepDoIt in terms of secondary electron.
237 //
238 
239 void G4PAIPhotModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
240  const G4MaterialCutsCouple* matCC,
241  const G4DynamicParticle* dp,
242  G4double tmin,
243  G4double maxEnergy)
244 {
245  G4int coupleIndex = FindCoupleIndex(matCC);
246  if(0 > coupleIndex) { return; }
247 
248  SetParticle(dp->GetDefinition());
249 
250  G4double kineticEnergy = dp->GetKineticEnergy();
251 
252  G4double tmax = MaxSecondaryEnergy(fParticle, kineticEnergy);
253 
254  if( maxEnergy < tmax) tmax = maxEnergy;
255  if( tmin >= tmax) return;
256 
257  G4ThreeVector direction = dp->GetMomentumDirection();
258  G4double scaledTkin = kineticEnergy*fRatio;
259  G4double totalEnergy = kineticEnergy + fMass;
260  G4double totalMomentum = sqrt(kineticEnergy*(totalEnergy + fMass));
261  G4double plRatio = fModelData->GetPlasmonRatio(coupleIndex, scaledTkin);
262 
263  if( G4UniformRand() <= plRatio )
264  {
265  G4double deltaTkin = fModelData->SamplePostStepPlasmonTransfer(coupleIndex, scaledTkin);
266 
267  // G4cout<<"G4PAIPhotModel::SampleSecondaries; dp "<<dp->GetParticleDefinition()->GetParticleName()
268  // <<"; Tkin = "<<kineticEnergy/keV<<" keV; transfer = "<<deltaTkin/keV<<" keV "<<G4endl;
269 
270  if( deltaTkin <= 0. && fVerbose > 0)
271  {
272  G4cout<<"G4PAIPhotModel::SampleSecondary e- deltaTkin = "<<deltaTkin<<G4endl;
273  }
274  if( deltaTkin <= 0.) return;
275 
276  if( deltaTkin > tmax) deltaTkin = tmax;
277 
278  const G4Element* anElement = SelectRandomAtom(matCC,fParticle,kineticEnergy);
279  G4int Z = G4lrint(anElement->GetZ());
280 
281  G4DynamicParticle* deltaRay = new G4DynamicParticle(fElectron,
282  GetAngularDistribution()->SampleDirection(dp, deltaTkin,
283  Z, matCC->GetMaterial()),
284  deltaTkin);
285 
286  // primary change
287 
288  kineticEnergy -= deltaTkin;
289 
290  if( kineticEnergy <= 0. ) // kill primary as local? energy deposition
291  {
292  fParticleChange->SetProposedKineticEnergy(0.0);
293  fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy+deltaTkin);
294  // fParticleChange->ProposeTrackStatus(fStopAndKill);
295  return;
296  }
297  else
298  {
299  G4ThreeVector dir = totalMomentum*direction - deltaRay->GetMomentum();
300  direction = dir.unit();
301  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
302  fParticleChange->SetProposedMomentumDirection(direction);
303  vdp->push_back(deltaRay);
304  }
305  }
306  else // secondary X-ray CR photon
307  {
308  G4double deltaTkin = fModelData->SamplePostStepPhotonTransfer(coupleIndex, scaledTkin);
309 
310  // G4cout<<"PAIPhotonModel PhotonPostStepTransfer = "<<deltaTkin/keV<<" keV"<<G4endl;
311 
312  if( deltaTkin <= 0. )
313  {
314  G4cout<<"G4PAIPhotonModel::SampleSecondary gamma deltaTkin = "<<deltaTkin<<G4endl;
315  }
316  if( deltaTkin <= 0.) return;
317 
318  if( deltaTkin >= kineticEnergy ) // stop primary
319  {
320  deltaTkin = kineticEnergy;
321  kineticEnergy = 0.0;
322  }
323  G4double costheta = 0.; // G4UniformRand(); // VG: ??? for start only
324  G4double sintheta = sqrt((1.+costheta)*(1.-costheta));
325 
326  // direction of the 'Cherenkov' photon
327  G4double phi = twopi*G4UniformRand();
328  G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta;
329 
330  G4ThreeVector deltaDirection(dirx,diry,dirz);
331  deltaDirection.rotateUz(direction);
332 
333  if( kineticEnergy > 0.) // primary change
334  {
335  kineticEnergy -= deltaTkin;
336  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
337  }
338  else // stop primary, but pass X-ray CR
339  {
340  // fParticleChange->ProposeLocalEnergyDeposit(deltaTkin);
341  fParticleChange->SetProposedKineticEnergy(0.0);
342  }
343  // create G4DynamicParticle object for photon ray
344 
345  G4DynamicParticle* photonRay = new G4DynamicParticle;
346  photonRay->SetDefinition( G4Gamma::Gamma() );
347  photonRay->SetKineticEnergy( deltaTkin );
348  photonRay->SetMomentumDirection(deltaDirection);
349 
350  vdp->push_back(photonRay);
351  }
352  return;
353 }
354 
356 
358  const G4DynamicParticle* aParticle,
359  G4double, G4double step,
360  G4double eloss)
361 {
362  // return 0.;
363  G4int coupleIndex = FindCoupleIndex(matCC);
364  if(0 > coupleIndex) { return eloss; }
365 
366  SetParticle(aParticle->GetDefinition());
367 
368 
369  // G4cout << "G4PAIPhotModel::SampleFluctuations step(mm)= "<< step/mm
370  // << " Eloss(keV)= " << eloss/keV << " in "
371  // << matCC->GetMaterial()->GetName() << G4endl;
372 
373 
374  G4double Tkin = aParticle->GetKineticEnergy();
375  G4double scaledTkin = Tkin*fRatio;
376 
377  G4double loss = fModelData->SampleAlongStepPhotonTransfer(coupleIndex, Tkin,
378  scaledTkin,
379  step*fChargeSquare);
380  loss += fModelData->SampleAlongStepPlasmonTransfer(coupleIndex, Tkin,
381  scaledTkin,
382  step*fChargeSquare);
383 
384 
385  // G4cout<<" PAIPhotModel::SampleFluctuations loss = "<<loss/keV<<" keV, on step = "
386  // <<step/mm<<" mm"<<G4endl;
387  return loss;
388 
389 }
390 
392 //
393 // Returns the statistical estimation of the energy loss distribution variance
394 //
395 
396 
398  const G4DynamicParticle* aParticle,
399  G4double tmax,
400  G4double step )
401 {
402  G4double particleMass = aParticle->GetMass();
403  G4double electronDensity = material->GetElectronDensity();
404  G4double kineticEnergy = aParticle->GetKineticEnergy();
405  G4double q = aParticle->GetCharge()/eplus;
406  G4double etot = kineticEnergy + particleMass;
407  G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*particleMass)/(etot*etot);
408  G4double siga = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * step
409  * electronDensity * q * q;
410 
411  return siga;
412  /*
413  G4double loss, sumLoss=0., sumLoss2=0., sigma2, meanLoss=0.;
414  for(G4int i = 0; i < fMeanNumber; i++)
415  {
416  loss = SampleFluctuations(material,aParticle,tmax,step,meanLoss);
417  sumLoss += loss;
418  sumLoss2 += loss*loss;
419  }
420  meanLoss = sumLoss/fMeanNumber;
421  sigma2 = meanLoss*meanLoss + (sumLoss2-2*sumLoss*meanLoss)/fMeanNumber;
422  return sigma2;
423  */
424 }
425 
427 
429  G4double kinEnergy)
430 {
431  SetParticle(p);
432  G4double tmax = kinEnergy;
433  if(p == fElectron) { tmax *= 0.5; }
434  else if(p != fPositron) {
435  G4double ratio= electron_mass_c2/fMass;
436  G4double gamma= kinEnergy/fMass + 1.0;
437  tmax = 2.0*electron_mass_c2*(gamma*gamma - 1.) /
438  (1. + 2.0*gamma*ratio + ratio*ratio);
439  }
440  return tmax;
441 }
442 
444 
446 {
447  fPAIRegionVector.push_back(r);
448 }
449 
450 //
451 //
virtual G4double SampleFluctuations(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double, G4double, G4double) final
G4double SamplePostStepPhotonTransfer(G4int coupleIndex, G4double scaledTkin) const
virtual ~G4PAIPhotModel()
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:643
const G4String & GetName() const
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:118
G4double GetKineticEnergy() const
virtual void DefineForRegion(const G4Region *r) final
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4VEmModel.cc:146
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel) final
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:636
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) final
G4PAIPhotData * GetPAIPhotData()
const char * p
Definition: xmltok.h:285
const G4String & GetName() const
Definition: G4Material.hh:178
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) final
G4double GetPlasmonRatio(G4int coupleIndex, G4double scaledTkin) const
G4double GetZ() const
Definition: G4Element.hh:131
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:587
G4double SampleAlongStepPhotonTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
std::vector< G4Material * > G4MaterialTable
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:619
G4bool IsMaster() const
Definition: G4VEmModel.hh:725
void SetMomentumDirection(const G4ThreeVector &aDirection)
G4ParticleDefinition * GetDefinition() const
const std::vector< const G4MaterialCutsCouple * > & GetVectorOfCouples()
virtual G4double Dispersion(const G4Material *, const G4DynamicParticle *, G4double, G4double) final
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
void Initialise(const G4MaterialCutsCouple *, G4double cut, G4PAIPhotModel *)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
static const G4double reg
static constexpr double twopi
Definition: G4SIunits.hh:76
static G4RegionStore * GetInstance()
string material
Definition: eplot.py:19
virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kinEnergy) final
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy) final
const G4MaterialCutsCouple * CurrentCouple() const
Definition: G4VEmModel.hh:454
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy) final
G4double GetElectronDensity() const
Definition: G4Material.hh:217
G4double GetMass() const
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
G4PAIPhotModel(const G4ParticleDefinition *p=nullptr, const G4String &nam="PAI")
G4double DEDXPerVolume(G4int coupleIndex, G4double scaledTkin, G4double cut) const
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
G4double GetCharge() const
static constexpr double eplus
Definition: G4SIunits.hh:199
void SetProposedKineticEnergy(G4double proposedKinEnergy)
float electron_mass_c2
Definition: hepunit.py:274
G4MaterialCutsCouple * FindCouple(G4Material *mat)
const G4int n
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
std::vector< G4EmElementSelector * > * GetElementSelectors()
Definition: G4VEmModel.hh:809
static size_t GetNumberOfMaterials()
Definition: G4Material.cc:594
void SetProposedMomentumDirection(const G4ThreeVector &dir)
void SetKineticEnergy(G4double aEnergy)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4Positron * Positron()
Definition: G4Positron.cc:94
void SetElementSelectors(std::vector< G4EmElementSelector * > *)
Definition: G4VEmModel.hh:817
int G4lrint(double ad)
Definition: templates.hh:163
Hep3Vector unit() const
G4double SamplePostStepPlasmonTransfer(G4int coupleIndex, G4double scaledTkin) const
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:626
G4double CrossSectionPerVolume(G4int coupleIndex, G4double scaledTkin, G4double tcut, G4double tmax) const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double SampleAlongStepPlasmonTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:544
const G4Material * GetMaterial() const