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