ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PAIModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4PAIModel.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: G4PAIModel.cc
31 //
32 // Author: Vladimir.Grichine@cern.ch on base of V.Ivanchenko model interface
33 //
34 // Creation date: 05.10.2003
35 //
36 // Modifications:
37 //
38 // 17.08.04 V.Grichine, bug fixed for Tkin<=0 in SampleSecondary
39 // 16.08.04 V.Grichine, bug fixed in massRatio for DEDX, CrossSection,
40 // SampleSecondary
41 // 08.04.05 Major optimisation of internal interfaces (V.Ivantchenko)
42 // 26.07.09 Fixed logic to work with several materials (V.Ivantchenko)
43 // 21.11.10 V. Grichine verbose flag for protons and G4PAYySection to
44 // check sandia table
45 // 12.06.13 V. Grichine Bug fixed in SampleSecondaries for scaled Tkin
46 // (fMass -> proton_mass_c2)
47 // 19.08.13 V.Ivanchenko extract data handling to G4PAIModelData class
48 // added sharing of internal data between threads (MT migration)
49 //
50 
51 #include "G4PAIModel.hh"
52 
53 #include "G4SystemOfUnits.hh"
54 #include "G4PhysicalConstants.hh"
55 #include "G4Region.hh"
56 #include "G4MaterialCutsCouple.hh"
57 #include "G4MaterialTable.hh"
58 #include "G4RegionStore.hh"
59 
60 #include "Randomize.hh"
61 #include "G4Electron.hh"
62 #include "G4Positron.hh"
63 #include "G4Poisson.hh"
64 #include "G4Step.hh"
65 #include "G4Material.hh"
66 #include "G4DynamicParticle.hh"
67 #include "G4ParticleDefinition.hh"
69 #include "G4PAIModelData.hh"
70 #include "G4DeltaAngle.hh"
71 
73 
74 using namespace std;
75 
78  fVerbose(0),
79  fModelData(nullptr),
80  fParticle(nullptr)
81 {
84 
85  fParticleChange = nullptr;
86 
87  if(p) { SetParticle(p); }
88  else { SetParticle(fElectron); }
89 
90  // default generator
92  fLowestTcut = 12.5*CLHEP::eV;
93 }
94 
96 
98 {
99  if(IsMaster()) { delete fModelData; }
100 }
101 
103 
105  const G4DataVector& cuts)
106 {
107  if(fVerbose > 0) {
108  G4cout<<"G4PAIModel::Initialise for "<<p->GetParticleName()<<G4endl;
109  }
110  SetParticle(p);
112 
113  if(IsMaster()) {
114 
115  delete fModelData;
116  fMaterialCutsCoupleVector.clear();
117  if(fVerbose > 0) {
118  G4cout << "G4PAIModel instantiates data for " << p->GetParticleName()
119  << G4endl;
120  }
121  G4double tmin = LowEnergyLimit()*fRatio;
123  fModelData = new G4PAIModelData(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  // protect for unit tests
131  if(0 == numRegions) {
132  G4Exception("G4PAIModel::Initialise()","em0106",JustWarning,
133  "no G4Regions are registered for the PAI model - World is used");
135  ->GetRegion("DefaultRegionForTheWorld", false));
136  numRegions = 1;
137  }
138 
139  if(fVerbose > 0) {
140  G4cout << "G4PAIModel is defined for " << numRegions << " regions "
141  << G4endl;
142  G4cout << " total number of materials " << numOfMat << G4endl;
143  }
144  for(size_t iReg = 0; iReg<numRegions; ++iReg) {
145  const G4Region* curReg = fPAIRegionVector[iReg];
146  G4Region* reg = const_cast<G4Region*>(curReg);
147 
148  for(size_t jMat = 0; jMat<numOfMat; ++jMat) {
149  G4Material* mat = (*theMaterialTable)[jMat];
150  const G4MaterialCutsCouple* cutCouple = reg->FindCouple(mat);
151  size_t n = fMaterialCutsCoupleVector.size();
152  /*
153  G4cout << "Region: " << reg->GetName() << " " << reg
154  << " Couple " << cutCouple
155  << " PAI defined for " << n << " couples"
156  << " jMat= " << jMat << " " << mat->GetName()
157  << G4endl;
158  */
159  if(cutCouple) {
160  if(fVerbose > 0) {
161  G4cout << "Region <" << curReg->GetName() << "> mat <"
162  << mat->GetName() << "> CoupleIndex= "
163  << cutCouple->GetIndex()
164  << " " << p->GetParticleName()
165  << " cutsize= " << cuts.size() << G4endl;
166  }
167  // check if this couple is not already initialized
168  G4bool isnew = true;
169  if(0 < n) {
170  for(size_t i=0; i<n; ++i) {
171  if(cutCouple == fMaterialCutsCoupleVector[i]) {
172  isnew = false;
173  break;
174  }
175  }
176  }
177  // initialise data banks
178  //G4cout << " isNew: " << isnew << " " << cutCouple << G4endl;
179  if(isnew) {
180  fMaterialCutsCoupleVector.push_back(cutCouple);
181  fModelData->Initialise(cutCouple, this);
182  }
183  }
184  }
185  }
187  }
188 }
189 
191 
193  G4VEmModel* masterModel)
194 {
195  SetParticle(p);
196  fModelData = static_cast<G4PAIModel*>(masterModel)->GetPAIModelData();
198  static_cast<G4PAIModel*>(masterModel)->GetVectorOfCouples();
200 }
201 
203 
205  const G4MaterialCutsCouple*)
206 {
207  return fLowestTcut;
208 }
209 
211 
213  const G4ParticleDefinition* p,
214  G4double kineticEnergy,
215  G4double cutEnergy)
216 {
217  //G4cout << "===1=== " << CurrentCouple()
218  // << " idx= " << CurrentCouple()->GetIndex()
219  // << " " << fMaterialCutsCoupleVector[0]
220  // << G4endl;
221  G4int coupleIndex = FindCoupleIndex(CurrentCouple());
222  //G4cout << "===2=== " << coupleIndex << G4endl;
223  if(0 > coupleIndex) { return 0.0; }
224 
225  G4double cut = std::min(MaxSecondaryEnergy(p, kineticEnergy), cutEnergy);
226 
227  G4double scaledTkin = kineticEnergy*fRatio;
228 
229  return fChargeSquare*fModelData->DEDXPerVolume(coupleIndex, scaledTkin,
230  cut);
231 }
232 
234 
236  const G4ParticleDefinition* p,
237  G4double kineticEnergy,
238  G4double cutEnergy,
239  G4double maxEnergy )
240 {
241  //G4cout << "===3=== " << CurrentCouple()
242  // << " idx= " << CurrentCouple()->GetIndex()
243  // << " " << fMaterialCutsCoupleVector[0]
244  // << G4endl;
245  G4int coupleIndex = FindCoupleIndex(CurrentCouple());
246  //G4cout << "===4=== " << coupleIndex << G4endl;
247  if(0 > coupleIndex) { return 0.0; }
248 
249  G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy);
250  if(tmax <= cutEnergy) { return 0.0; }
251 
252  G4double scaledTkin = kineticEnergy*fRatio;
253 
254  return fChargeSquare*fModelData->CrossSectionPerVolume(coupleIndex,
255  scaledTkin,
256  cutEnergy,
257  tmax);
258 }
259 
261 //
262 // It is analog of PostStepDoIt in terms of secondary electron.
263 //
264 
265 void G4PAIModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
266  const G4MaterialCutsCouple* matCC,
267  const G4DynamicParticle* dp,
268  G4double tmin,
269  G4double maxEnergy)
270 {
271  G4int coupleIndex = FindCoupleIndex(matCC);
272 
273  //G4cout << "G4PAIModel::SampleSecondaries: coupleIndex= "<<coupleIndex<<G4endl;
274  if(0 > coupleIndex) { return; }
275 
276  SetParticle(dp->GetDefinition());
277  G4double kineticEnergy = dp->GetKineticEnergy();
278 
279  G4double tmax = MaxSecondaryEnergy(fParticle, kineticEnergy);
280  if(maxEnergy < tmax) { tmax = maxEnergy; }
281  if(tmin >= tmax) { return; }
282 
283  G4ThreeVector direction= dp->GetMomentumDirection();
284  G4double scaledTkin = kineticEnergy*fRatio;
285  G4double totalEnergy = kineticEnergy + fMass;
286  G4double totalMomentum = sqrt(kineticEnergy*(totalEnergy+fMass));
287 
288  G4double deltaTkin =
289  fModelData->SamplePostStepTransfer(coupleIndex, scaledTkin, tmin, tmax);
290 
291  //G4cout<<"G4PAIModel::SampleSecondaries; deltaKIn = "<<deltaTkin/keV
292  // <<" keV "<< " Escaled(MeV)= " << scaledTkin << G4endl;
293 
294  if( !(deltaTkin <= 0.) && !(deltaTkin > 0)) {
295  G4cout<<"G4PAIModel::SampleSecondaries; deltaKIn = "<<deltaTkin/keV
296  <<" keV "<< " Escaled(MeV)= " << scaledTkin << G4endl;
297  return;
298  }
299  if( deltaTkin <= 0.) { return; }
300 
301  if( deltaTkin > tmax) { deltaTkin = tmax; }
302 
303  const G4Element* anElement = SelectTargetAtom(matCC, fParticle, kineticEnergy,
304  dp->GetLogKineticEnergy());
305 
306  G4int Z = G4lrint(anElement->GetZ());
307 
309  GetAngularDistribution()->SampleDirection(dp, deltaTkin,
310  Z, matCC->GetMaterial()),
311  deltaTkin);
312 
313  // primary change
314  kineticEnergy -= deltaTkin;
315  G4ThreeVector dir = totalMomentum*direction - deltaRay->GetMomentum();
316  direction = dir.unit();
319 
320  vdp->push_back(deltaRay);
321 }
322 
324 
326  const G4DynamicParticle* aParticle,
327  G4double tmax, G4double step,
328  G4double eloss)
329 {
330  G4int coupleIndex = FindCoupleIndex(matCC);
331  if(0 > coupleIndex) { return eloss; }
332 
333  SetParticle(aParticle->GetDefinition());
334 
335  /*
336  G4cout << "G4PAIModel::SampleFluctuations step(mm)= "<< step/mm
337  << " Eloss(keV)= " << eloss/keV << " in "
338  << matCC->Getmaterial()->GetName() << G4endl;
339  */
340 
341  G4double Tkin = aParticle->GetKineticEnergy();
342  G4double scaledTkin = Tkin*fRatio;
343 
344  G4double loss = fModelData->SampleAlongStepTransfer(coupleIndex, Tkin,
345  scaledTkin, tmax,
346  step*fChargeSquare);
347 
348  // G4cout<<"PAIModel AlongStepLoss = "<<loss/keV<<" keV, on step = "
349  //<<step/mm<<" mm"<<G4endl;
350  return loss;
351 
352 }
353 
355 //
356 // Returns the statistical estimation of the energy loss distribution variance
357 //
358 
359 
361  const G4DynamicParticle* aParticle,
362  G4double tmax,
363  G4double step )
364 {
365  G4double particleMass = aParticle->GetMass();
366  G4double electronDensity = material->GetElectronDensity();
367  G4double kineticEnergy = aParticle->GetKineticEnergy();
368  G4double q = aParticle->GetCharge()/eplus;
369  G4double etot = kineticEnergy + particleMass;
370  G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*particleMass)/(etot*etot);
371  G4double siga = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * step
372  * electronDensity * q * q;
373 
374  return siga;
375 }
376 
378 
380  G4double kinEnergy)
381 {
382  SetParticle(p);
383  G4double tmax = kinEnergy;
384  if(p == fElectron) { tmax *= 0.5; }
385  else if(p != fPositron) {
387  G4double gamma= kinEnergy/fMass + 1.0;
388  tmax = 2.0*electron_mass_c2*(gamma*gamma - 1.) /
389  (1. + 2.0*gamma*ratio + ratio*ratio);
390  }
391  return tmax;
392 }
393 
395 
397 {
398  fPAIRegionVector.push_back(r);
399 }
400 
402