ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PAIPhotModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 //
27 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class
30 // File name: G4PAIPhotModel.cc
31 //
32 // Author: Vladimir.Grichine@cern.ch on base of G4PAIModel MT interface
33 //
34 // Creation date: 07.10.2013
35 //
36 // Modifications:
37 //
38 //
39 
40 #include "G4PAIPhotModel.hh"
41 
42 #include "G4SystemOfUnits.hh"
43 #include "G4PhysicalConstants.hh"
44 #include "G4Region.hh"
45 #include "G4PhysicsLogVector.hh"
46 #include "G4PhysicsFreeVector.hh"
47 #include "G4PhysicsTable.hh"
48 #include "G4ProductionCutsTable.hh"
49 #include "G4MaterialCutsCouple.hh"
50 #include "G4MaterialTable.hh"
51 #include "G4SandiaTable.hh"
52 #include "G4OrderedTable.hh"
53 #include "G4RegionStore.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 
69 
70 using namespace std;
71 
74  fVerbose(0),
75  fModelData(nullptr),
76  fParticle(nullptr)
77 {
80 
81  fParticleChange = nullptr;
82 
83  if(p) { SetParticle(p); }
84  else { SetParticle(fElectron); }
85 
86  // default generator
88  fLowestTcut = 12.5*CLHEP::eV;
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);
110 
111  if( IsMaster() )
112  {
113 
115 
116  delete fModelData;
117  fMaterialCutsCoupleVector.clear();
118 
119  G4double tmin = LowEnergyLimit()*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");
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 G4MaterialCutsCouple*)
198 {
199  return fLowestTcut;
200 }
201 
203 
205  const G4ParticleDefinition* p,
206  G4double kineticEnergy,
207  G4double cutEnergy)
208 {
209  G4int coupleIndex = FindCoupleIndex(CurrentCouple());
210  if(0 > coupleIndex) { return 0.0; }
211 
212  G4double cut = std::min(MaxSecondaryEnergy(p, kineticEnergy), cutEnergy);
213 
214  G4double scaledTkin = kineticEnergy*fRatio;
215 
216  return fChargeSquare*fModelData->DEDXPerVolume(coupleIndex, scaledTkin, cut);
217 }
218 
220 
222  const G4ParticleDefinition* p,
223  G4double kineticEnergy,
224  G4double cutEnergy,
225  G4double maxEnergy )
226 {
227  G4int coupleIndex = FindCoupleIndex(CurrentCouple());
228 
229  if(0 > coupleIndex) return 0.0;
230 
231  G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy);
232 
233  if(tmax <= cutEnergy) return 0.0;
234 
235  G4double scaledTkin = kineticEnergy*fRatio;
237  scaledTkin,
238  cutEnergy, tmax);
239  return xsc;
240 }
241 
243 //
244 // It is analog of PostStepDoIt in terms of secondary electron.
245 //
246 
247 void G4PAIPhotModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
248  const G4MaterialCutsCouple* matCC,
249  const G4DynamicParticle* dp,
250  G4double tmin,
251  G4double maxEnergy)
252 {
253  G4int coupleIndex = FindCoupleIndex(matCC);
254  if(0 > coupleIndex) { return; }
255 
256  SetParticle(dp->GetDefinition());
257 
258  G4double kineticEnergy = dp->GetKineticEnergy();
259 
260  G4double tmax = MaxSecondaryEnergy(fParticle, kineticEnergy);
261 
262  if( maxEnergy < tmax) tmax = maxEnergy;
263  if( tmin >= tmax) return;
264 
265  G4ThreeVector direction = dp->GetMomentumDirection();
266  G4double scaledTkin = kineticEnergy*fRatio;
267  G4double totalEnergy = kineticEnergy + fMass;
268  G4double totalMomentum = sqrt(kineticEnergy*(totalEnergy + fMass));
269  G4double plRatio = fModelData->GetPlasmonRatio(coupleIndex, scaledTkin);
270 
271  if( G4UniformRand() <= plRatio )
272  {
273  G4double deltaTkin = fModelData->SamplePostStepPlasmonTransfer(coupleIndex, scaledTkin);
274 
275  // G4cout<<"G4PAIPhotModel::SampleSecondaries; dp "<<dp->GetParticleDefinition()->GetParticleName()
276  // <<"; Tkin = "<<kineticEnergy/keV<<" keV; transfer = "<<deltaTkin/keV<<" keV "<<G4endl;
277 
278  if( deltaTkin <= 0. && fVerbose > 0)
279  {
280  G4cout<<"G4PAIPhotModel::SampleSecondary e- deltaTkin = "<<deltaTkin<<G4endl;
281  }
282  if( deltaTkin <= 0.) return;
283 
284  if( deltaTkin > tmax) deltaTkin = tmax;
285 
286  const G4Element* anElement = SelectTargetAtom(matCC,fParticle,kineticEnergy,
287  dp->GetLogKineticEnergy());
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
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,
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 //