ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4EmCalculator.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4EmCalculator.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 file
30 //
31 //
32 // File name: G4EmCalculator
33 //
34 // Author: Vladimir Ivanchenko
35 //
36 // Creation date: 28.06.2004
37 //
38 //
39 // Class Description: V.Ivanchenko & M.Novak
40 //
41 // -------------------------------------------------------------------
42 //
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45 
46 #include "G4EmCalculator.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "G4LossTableManager.hh"
49 #include "G4EmParameters.hh"
50 #include "G4NistManager.hh"
51 #include "G4VEmProcess.hh"
52 #include "G4VEnergyLossProcess.hh"
53 #include "G4VMultipleScattering.hh"
54 #include "G4Material.hh"
55 #include "G4MaterialCutsCouple.hh"
56 #include "G4ParticleDefinition.hh"
57 #include "G4ParticleTable.hh"
58 #include "G4IonTable.hh"
59 #include "G4PhysicsTable.hh"
60 #include "G4ProductionCutsTable.hh"
61 #include "G4ProcessManager.hh"
62 #include "G4ionEffectiveCharge.hh"
63 #include "G4RegionStore.hh"
64 #include "G4Element.hh"
65 #include "G4EmCorrections.hh"
66 #include "G4GenericIon.hh"
67 #include "G4ProcessVector.hh"
68 #include "G4Gamma.hh"
69 #include "G4Electron.hh"
70 #include "G4Positron.hh"
71 
72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73 
74 G4EmCalculator::G4EmCalculator()
75 {
76  manager = G4LossTableManager::Instance();
77  nist = G4NistManager::Instance();
78  theParameters = G4EmParameters::Instance();
79  corr = manager->EmCorrections();
80  nLocalMaterials = 0;
81  verbose = 0;
82  currentCoupleIndex = 0;
83  currentCouple = nullptr;
84  currentMaterial = cutMaterial = nullptr;
85  currentParticle = nullptr;
86  lambdaParticle = nullptr;
87  baseParticle = nullptr;
88  currentLambda = nullptr;
89  currentModel = nullptr;
90  currentProcess = nullptr;
91  curProcess = nullptr;
92  loweModel = nullptr;
93  chargeSquare = 1.0;
94  massRatio = 1.0;
95  mass = 0.0;
96  currentCut = 0.0;
97  cutenergy[0] = cutenergy[1] = cutenergy[2] = DBL_MAX;
98  currentParticleName= "";
99  currentMaterialName= "";
100  currentName = "";
101  lambdaName = "";
102  theGenericIon = G4GenericIon::GenericIon();
103  ionEffCharge = new G4ionEffectiveCharge();
104  ionTable = G4ParticleTable::GetParticleTable()->GetIonTable();
105  isIon = false;
106  isApplicable = false;
107 }
108 
109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
110 
111 G4EmCalculator::~G4EmCalculator()
112 {
113  delete ionEffCharge;
114  for (G4int i=0; i<nLocalMaterials; ++i) {
115  delete localCouples[i];
116  }
117 }
118 
119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
120 
121 G4double G4EmCalculator::GetDEDX(G4double kinEnergy,
122  const G4ParticleDefinition* p,
123  const G4Material* mat,
124  const G4Region* region)
125 {
126  G4double res = 0.0;
127  const G4MaterialCutsCouple* couple = FindCouple(mat, region);
128  if(couple && UpdateParticle(p, kinEnergy) ) {
129  res = manager->GetDEDX(p, kinEnergy, couple);
130 
131  if(isIon) {
132  if(FindEmModel(p, currentProcessName, kinEnergy)) {
133  G4double length = CLHEP::nm;
134  G4double eloss = res*length;
135  //G4cout << "### GetDEDX: E= " << kinEnergy << " dedx0= " << res
136  // << " de= " << eloss << G4endl;;
137  G4double niel = 0.0;
138  dynParticle.SetKineticEnergy(kinEnergy);
139  currentModel->GetChargeSquareRatio(p, mat, kinEnergy);
140  currentModel->CorrectionsAlongStep(couple,&dynParticle,eloss,niel,length);
141  res = eloss/length;
142  //G4cout << " de1= " << eloss << " res1= " << res
143  // << " " << p->GetParticleName() <<G4endl;;
144  }
145  }
146 
147  if(verbose>0) {
148  G4cout << "G4EmCalculator::GetDEDX: E(MeV)= " << kinEnergy/MeV
149  << " DEDX(MeV/mm)= " << res*mm/MeV
150  << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
151  << " " << p->GetParticleName()
152  << " in " << mat->GetName()
153  << " isIon= " << isIon
154  << G4endl;
155  }
156  }
157  return res;
158 }
159 
160 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
161 
162 G4double G4EmCalculator::GetRangeFromRestricteDEDX(G4double kinEnergy,
163  const G4ParticleDefinition* p,
164  const G4Material* mat,
165  const G4Region* region)
166 {
167  G4double res = 0.0;
168  const G4MaterialCutsCouple* couple = FindCouple(mat,region);
169  if(couple && UpdateParticle(p, kinEnergy)) {
170  res = manager->GetRangeFromRestricteDEDX(p, kinEnergy, couple);
171  if(verbose>1) {
172  G4cout << " G4EmCalculator::GetRangeFromRestrictedDEDX: E(MeV)= "
173  << kinEnergy/MeV
174  << " range(mm)= " << res/mm
175  << " " << p->GetParticleName()
176  << " in " << mat->GetName()
177  << G4endl;
178  }
179  }
180  return res;
181 }
182 
183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
184 
185 G4double G4EmCalculator::GetCSDARange(G4double kinEnergy,
186  const G4ParticleDefinition* p,
187  const G4Material* mat,
188  const G4Region* region)
189 {
190  G4double res = 0.0;
191  if(!theParameters->BuildCSDARange()) {
192  G4ExceptionDescription ed;
193  ed << "G4EmCalculator::GetCSDARange: CSDA table is not built; "
194  << " use UI command: /process/eLoss/CSDARange true";
195  G4Exception("G4EmCalculator::GetCSDARange", "em0077",
196  JustWarning, ed);
197  return res;
198  }
199 
200  const G4MaterialCutsCouple* couple = FindCouple(mat,region);
201  if(couple && UpdateParticle(p, kinEnergy)) {
202  res = manager->GetCSDARange(p, kinEnergy, couple);
203  if(verbose>1) {
204  G4cout << " G4EmCalculator::GetCSDARange: E(MeV)= " << kinEnergy/MeV
205  << " range(mm)= " << res/mm
206  << " " << p->GetParticleName()
207  << " in " << mat->GetName()
208  << G4endl;
209  }
210  }
211  return res;
212 }
213 
214 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
215 
216 G4double G4EmCalculator::GetRange(G4double kinEnergy,
217  const G4ParticleDefinition* p,
218  const G4Material* mat,
219  const G4Region* region)
220 {
221  G4double res = 0.0;
222  if(theParameters->BuildCSDARange()) {
223  res = GetCSDARange(kinEnergy, p, mat, region);
224  } else {
225  res = GetRangeFromRestricteDEDX(kinEnergy, p, mat, region);
226  }
227  return res;
228 }
229 
230 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
231 
232 G4double G4EmCalculator::GetKinEnergy(G4double range,
233  const G4ParticleDefinition* p,
234  const G4Material* mat,
235  const G4Region* region)
236 {
237  G4double res = 0.0;
238  const G4MaterialCutsCouple* couple = FindCouple(mat,region);
239  if(couple && UpdateParticle(p, 1.0*GeV)) {
240  res = manager->GetEnergy(p, range, couple);
241  if(verbose>0) {
242  G4cout << "G4EmCalculator::GetKinEnergy: Range(mm)= " << range/mm
243  << " KinE(MeV)= " << res/MeV
244  << " " << p->GetParticleName()
245  << " in " << mat->GetName()
246  << G4endl;
247  }
248  }
249  return res;
250 }
251 
252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
253 
254 G4double G4EmCalculator::GetCrossSectionPerVolume(G4double kinEnergy,
255  const G4ParticleDefinition* p,
256  const G4String& processName,
257  const G4Material* mat,
258  const G4Region* region)
259 {
260  G4double res = 0.0;
261  const G4MaterialCutsCouple* couple = FindCouple(mat,region);
262 
263  if(couple && UpdateParticle(p, kinEnergy)) {
264  if(FindEmModel(p, processName, kinEnergy)) {
265  G4int idx = couple->GetIndex();
266  G4int procType = -1;
267  FindLambdaTable(p, processName, kinEnergy, procType);
268 
269  G4VEmProcess* emproc = FindDiscreteProcess(p, processName);
270  if(emproc) {
271  res = emproc->CrossSectionPerVolume(kinEnergy, couple);
272  } else if(currentLambda) {
273  // special tables are built for Msc models (procType is set in FindLambdaTable
274  if(procType==2) {
275  G4VMscModel* mscM = static_cast<G4VMscModel*>(currentModel);
276  mscM->SetCurrentCouple(couple);
277  G4double tr1Mfp = mscM->GetTransportMeanFreePath(p, kinEnergy);
278  if (tr1Mfp<DBL_MAX) {
279  res = 1./tr1Mfp;
280  }
281  } else {
282  G4double e = kinEnergy*massRatio;
283  res = (((*currentLambda)[idx])->Value(e))*chargeSquare;
284  }
285  } else {
286  res = ComputeCrossSectionPerVolume(kinEnergy, p, processName, mat, kinEnergy);
287  }
288  if(verbose>0) {
289  G4cout << "G4EmCalculator::GetXSPerVolume: E(MeV)= " << kinEnergy/MeV
290  << " cross(cm-1)= " << res*cm
291  << " " << p->GetParticleName()
292  << " in " << mat->GetName();
293  if(verbose>1)
294  G4cout << " idx= " << idx << " Escaled((MeV)= "
295  << kinEnergy*massRatio
296  << " q2= " << chargeSquare;
297  G4cout << G4endl;
298  }
299  }
300  }
301  return res;
302 }
303 
304 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
305 
306 G4double G4EmCalculator::GetShellIonisationCrossSectionPerAtom(
307  const G4String& particle,
308  G4int Z,
309  G4AtomicShellEnumerator shell,
310  G4double kinEnergy)
311 {
312  G4double res = 0.0;
313  const G4ParticleDefinition* p = FindParticle(particle);
314  G4VAtomDeexcitation* ad = manager->AtomDeexcitation();
315  if(p && ad) {
316  res = ad->GetShellIonisationCrossSectionPerAtom(p, Z, shell, kinEnergy);
317  }
318  return res;
319 }
320 
321 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
322 
323 G4double G4EmCalculator::GetMeanFreePath(G4double kinEnergy,
324  const G4ParticleDefinition* p,
325  const G4String& processName,
326  const G4Material* mat,
327  const G4Region* region)
328 {
329  G4double res = DBL_MAX;
330  G4double x = GetCrossSectionPerVolume(kinEnergy,p, processName, mat,region);
331  if(x > 0.0) { res = 1.0/x; }
332  if(verbose>1) {
333  G4cout << "G4EmCalculator::GetMeanFreePath: E(MeV)= " << kinEnergy/MeV
334  << " MFP(mm)= " << res/mm
335  << " " << p->GetParticleName()
336  << " in " << mat->GetName()
337  << G4endl;
338  }
339  return res;
340 }
341 
342 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
343 
344 void G4EmCalculator::PrintDEDXTable(const G4ParticleDefinition* p)
345 {
346  const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p);
347  G4cout << "##### DEDX Table for " << p->GetParticleName() << G4endl;
348  if(elp) G4cout << *(elp->DEDXTable()) << G4endl;
349 }
350 
351 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
352 
353 void G4EmCalculator::PrintRangeTable(const G4ParticleDefinition* p)
354 {
355  const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p);
356  G4cout << "##### Range Table for " << p->GetParticleName() << G4endl;
357  if(elp) G4cout << *(elp->RangeTableForLoss()) << G4endl;
358 }
359 
360 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
361 
362 void G4EmCalculator::PrintInverseRangeTable(const G4ParticleDefinition* p)
363 {
364  const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p);
365  G4cout << "### G4EmCalculator: Inverse Range Table for "
366  << p->GetParticleName() << G4endl;
367  if(elp) G4cout << *(elp->InverseRangeTable()) << G4endl;
368 }
369 
370 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
371 
372 G4double G4EmCalculator::ComputeDEDX(G4double kinEnergy,
373  const G4ParticleDefinition* p,
374  const G4String& processName,
375  const G4Material* mat,
376  G4double cut)
377 {
378  SetupMaterial(mat);
379  G4double res = 0.0;
380  if(verbose > 1) {
381  G4cout << "### G4EmCalculator::ComputeDEDX: " << p->GetParticleName()
382  << " in " << currentMaterialName
383  << " e(MeV)= " << kinEnergy/MeV << " cut(MeV)= " << cut/MeV
384  << G4endl;
385  }
386  if(UpdateParticle(p, kinEnergy)) {
387  if(FindEmModel(p, processName, kinEnergy)) {
388 
389  // Special case of ICRU'73 model
390  const G4String& mname = currentModel->GetName();
391  if(mname == "ParamICRU73" || mname == "LinhardSorensen" || mname == "Atima") {
392  res = currentModel->ComputeDEDXPerVolume(mat, p, kinEnergy, cut);
393  if(verbose > 1) {
394  G4cout << " ICRU73 ion E(MeV)= " << kinEnergy << " ";
395  G4cout << currentModel->GetName() << ": DEDX(MeV/mm)= " << res*mm/MeV
396  << " DEDX(MeV*cm^2/g)= "
397  << res*gram/(MeV*cm2*mat->GetDensity())
398  << G4endl;
399  }
400  } else {
401 
402  G4double escaled = kinEnergy*massRatio;
403  if(baseParticle) {
404  res = currentModel->ComputeDEDXPerVolume(
405  mat, baseParticle, escaled, cut) * chargeSquare;
406  if(verbose > 1) {
407  G4cout << baseParticle->GetParticleName()
408  << " Escaled(MeV)= " << escaled;
409  }
410  } else {
411  res = currentModel->ComputeDEDXPerVolume(mat, p, kinEnergy, cut);
412  if(verbose > 1) { G4cout << " no basePart E(MeV)= " << kinEnergy << " "; }
413  }
414  if(verbose > 1) {
415  G4cout << currentModel->GetName() << ": DEDX(MeV/mm)= " << res*mm/MeV
416  << " DEDX(MeV*cm^2/g)= "
417  << res*gram/(MeV*cm2*mat->GetDensity())
418  << G4endl;
419  }
420 
421  // emulate smoothing procedure
422  G4double eth = currentModel->LowEnergyLimit();
423  // G4cout << "massRatio= " << massRatio << " eth= " << eth << G4endl;
424  if(loweModel) {
425  G4double res0 = 0.0;
426  G4double res1 = 0.0;
427  if(baseParticle) {
428  res1 = currentModel->ComputeDEDXPerVolume(mat, baseParticle, eth, cut)
429  * chargeSquare;
430  res0 = loweModel->ComputeDEDXPerVolume(mat, baseParticle, eth, cut)
431  * chargeSquare;
432  } else {
433  res1 = currentModel->ComputeDEDXPerVolume(mat, p, eth, cut);
434  res0 = loweModel->ComputeDEDXPerVolume(mat, p, eth, cut);
435  }
436  if(verbose > 1) {
437  G4cout << "At boundary energy(MeV)= " << eth/MeV
438  << " DEDX(MeV/mm)= " << res1*mm/MeV
439  << G4endl;
440  }
441 
442  //G4cout << "eth= " << eth << " escaled= " << escaled
443  // << " res0= " << res0 << " res1= "
444  // << res1 << " q2= " << chargeSquare << G4endl;
445 
446  if(res1 > 0.0 && escaled > 0.0) {
447  res *= (1.0 + (res0/res1 - 1.0)*eth/escaled);
448  }
449  }
450 
451  // low energy correction for ions
452  if(isIon) {
453  G4double length = CLHEP::nm;
454  const G4Region* r = 0;
455  const G4MaterialCutsCouple* couple = FindCouple(mat, r);
456  G4double eloss = res*length;
457  G4double niel = 0.0;
458  dynParticle.SetKineticEnergy(kinEnergy);
459  currentModel->GetChargeSquareRatio(p, mat, kinEnergy);
460  currentModel->CorrectionsAlongStep(couple,&dynParticle,eloss,niel,length);
461  res = eloss/length;
462 
463  if(verbose > 1) {
464  G4cout << "After Corrections: DEDX(MeV/mm)= " << res*mm/MeV
465  << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
466  << G4endl;
467  }
468  }
469  }
470  }
471  if(verbose > 0) {
472  G4cout << "Sum: E(MeV)= " << kinEnergy/MeV
473  << " DEDX(MeV/mm)= " << res*mm/MeV
474  << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
475  << " cut(MeV)= " << cut/MeV
476  << " " << p->GetParticleName()
477  << " in " << currentMaterialName
478  << " Zi^2= " << chargeSquare
479  << " isIon=" << isIon
480  << G4endl;
481  }
482  }
483  return res;
484 }
485 
486 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
487 
488 G4double G4EmCalculator::ComputeElectronicDEDX(G4double kinEnergy,
489  const G4ParticleDefinition* part,
490  const G4Material* mat,
491  G4double cut)
492 {
493  SetupMaterial(mat);
494  G4double dedx = 0.0;
495  if(UpdateParticle(part, kinEnergy)) {
496 
497  G4LossTableManager* lManager = G4LossTableManager::Instance();
498  const std::vector<G4VEnergyLossProcess*> vel =
499  lManager->GetEnergyLossProcessVector();
500  G4int n = vel.size();
501 
502  //G4cout << "ComputeElectronicDEDX for " << part->GetParticleName()
503  // << " n= " << n << G4endl;
504 
505  for(G4int i=0; i<n; ++i) {
506  if(vel[i]) {
507  G4VProcess* p = reinterpret_cast<G4VProcess*>(vel[i]);
508  if(ActiveForParticle(part, p)) {
509  //G4cout << "idx= " << i << " " << (vel[i])->GetProcessName()
510  // << " " << (vel[i])->Particle()->GetParticleName() << G4endl;
511  dedx += ComputeDEDX(kinEnergy,part,(vel[i])->GetProcessName(),mat,cut);
512  }
513  }
514  }
515  }
516  return dedx;
517 }
518 
519 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
520 
521 G4double G4EmCalculator::ComputeDEDXForCutInRange(G4double kinEnergy,
522  const G4ParticleDefinition* part,
523  const G4Material* mat,
524  G4double rangecut)
525 {
526  SetupMaterial(mat);
527  G4double dedx = 0.0;
528  if(UpdateParticle(part, kinEnergy)) {
529 
530  G4LossTableManager* lManager = G4LossTableManager::Instance();
531  const std::vector<G4VEnergyLossProcess*> vel =
532  lManager->GetEnergyLossProcessVector();
533  G4int n = vel.size();
534 
535  if(mat != cutMaterial) {
536  cutMaterial = mat;
537  cutenergy[0] = ComputeEnergyCutFromRangeCut(rangecut, G4Gamma::Gamma(), mat);
538  cutenergy[1] = ComputeEnergyCutFromRangeCut(rangecut, G4Electron::Electron(), mat);
539  cutenergy[2] = ComputeEnergyCutFromRangeCut(rangecut, G4Positron::Positron(), mat);
540  }
541 
542  //G4cout << "ComputeElectronicDEDX for " << part->GetParticleName()
543  // << " n= " << n << G4endl;
544 
545  for(G4int i=0; i<n; ++i) {
546  if(vel[i]) {
547  G4VProcess* p = reinterpret_cast<G4VProcess*>(vel[i]);
548  if(ActiveForParticle(part, p)) {
549  //G4cout << "idx= " << i << " " << (vel[i])->GetProcessName()
550  // << " " << (vel[i])->Particle()->GetParticleName() << G4endl;
551  const G4ParticleDefinition* sec = (vel[i])->SecondaryParticle();
552  G4int idx = 0;
553  if(sec == G4Electron::Electron()) { idx = 1; }
554  else if(sec == G4Positron::Positron()) { idx = 2; }
555 
556  dedx += ComputeDEDX(kinEnergy,part,(vel[i])->GetProcessName(),
557  mat,cutenergy[idx]);
558  }
559  }
560  }
561  }
562  return dedx;
563 }
564 
565 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
566 
567 G4double G4EmCalculator::ComputeTotalDEDX(G4double kinEnergy,
568  const G4ParticleDefinition* part,
569  const G4Material* mat,
570  G4double cut)
571 {
572  G4double dedx = ComputeElectronicDEDX(kinEnergy,part,mat,cut);
573  if(mass > 700.*MeV) { dedx += ComputeNuclearDEDX(kinEnergy,part,mat); }
574  return dedx;
575 }
576 
577 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
578 
579 G4double G4EmCalculator::ComputeNuclearDEDX(G4double kinEnergy,
580  const G4ParticleDefinition* p,
581  const G4Material* mat)
582 {
583  G4double res = 0.0;
584  G4VEmProcess* nucst = FindDiscreteProcess(p, "nuclearStopping");
585  if(nucst) {
586  G4VEmModel* mod = nucst->GetModelByIndex();
587  if(mod) {
588  mod->SetFluctuationFlag(false);
589  res = mod->ComputeDEDXPerVolume(mat, p, kinEnergy);
590  }
591  }
592 
593  if(verbose > 1) {
594  G4cout << p->GetParticleName() << " E(MeV)= " << kinEnergy/MeV
595  << " NuclearDEDX(MeV/mm)= " << res*mm/MeV
596  << " NuclearDEDX(MeV*cm^2/g)= "
597  << res*gram/(MeV*cm2*mat->GetDensity())
598  << G4endl;
599  }
600  return res;
601 }
602 
603 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
604 
605 G4double G4EmCalculator::ComputeCrossSectionPerVolume(
606  G4double kinEnergy,
607  const G4ParticleDefinition* p,
608  const G4String& processName,
609  const G4Material* mat,
610  G4double cut)
611 {
612  SetupMaterial(mat);
613  G4double res = 0.0;
614  if(UpdateParticle(p, kinEnergy)) {
615  if(FindEmModel(p, processName, kinEnergy)) {
616  G4double e = kinEnergy;
617  G4double aCut = std::max(cut, theParameters->LowestElectronEnergy());
618  if(baseParticle) {
619  e *= kinEnergy*massRatio;
620  res = currentModel->CrossSectionPerVolume(
621  mat, baseParticle, e, aCut, e) * chargeSquare;
622  } else {
623  res = currentModel->CrossSectionPerVolume(mat, p, e, aCut, e);
624  }
625  if(verbose>0) {
626  G4cout << "G4EmCalculator::ComputeXSPerVolume: E(MeV)= " << kinEnergy/MeV
627  << " cross(cm-1)= " << res*cm
628  << " cut(keV)= " << aCut/keV
629  << " " << p->GetParticleName()
630  << " in " << mat->GetName()
631  << G4endl;
632  }
633  }
634  }
635  return res;
636 }
637 
638 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
639 
640 G4double G4EmCalculator::ComputeCrossSectionPerAtom(
641  G4double kinEnergy,
642  const G4ParticleDefinition* p,
643  const G4String& processName,
644  G4double Z, G4double A,
645  G4double cut)
646 {
647  G4double res = 0.0;
648  if(UpdateParticle(p, kinEnergy)) {
649  G4int iz = G4lrint(Z);
650  CheckMaterial(iz);
651  if(FindEmModel(p, processName, kinEnergy)) {
652  G4double e = kinEnergy;
653  G4double aCut = std::max(cut, theParameters->LowestElectronEnergy());
654  if(baseParticle) {
655  e *= kinEnergy*massRatio;
656  currentModel->InitialiseForElement(baseParticle, iz);
657  res = currentModel->ComputeCrossSectionPerAtom(
658  baseParticle, e, Z, A, aCut) * chargeSquare;
659  } else {
660  currentModel->InitialiseForElement(p, iz);
661  res = currentModel->ComputeCrossSectionPerAtom(p, e, Z, A, aCut);
662  }
663  if(verbose>0) {
664  G4cout << "E(MeV)= " << kinEnergy/MeV
665  << " cross(barn)= " << res/barn
666  << " " << p->GetParticleName()
667  << " Z= " << Z << " A= " << A/(g/mole) << " g/mole"
668  << " cut(keV)= " << aCut/keV
669  << G4endl;
670  }
671  }
672  }
673  return res;
674 }
675 
676 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
677 
678 G4double G4EmCalculator::ComputeCrossSectionPerShell(G4double kinEnergy,
679  const G4ParticleDefinition* p,
680  const G4String& processName,
681  G4int Z, G4int shellIdx,
682  G4double cut)
683 {
684  G4double res = 0.0;
685  if(UpdateParticle(p, kinEnergy)) {
686  CheckMaterial(Z);
687  if(FindEmModel(p, processName, kinEnergy)) {
688  G4double e = kinEnergy;
689  G4double aCut = std::max(cut, theParameters->LowestElectronEnergy());
690  if(baseParticle) {
691  e *= kinEnergy*massRatio;
692  currentModel->InitialiseForElement(baseParticle, Z);
693  res = currentModel->ComputeCrossSectionPerShell(baseParticle, Z, shellIdx,
694  e, aCut) * chargeSquare;
695  } else {
696  currentModel->InitialiseForElement(p, Z);
697  res = currentModel->ComputeCrossSectionPerAtom(p, Z, shellIdx, e, aCut);
698  }
699  if(verbose>0) {
700  G4cout << "E(MeV)= " << kinEnergy/MeV
701  << " cross(barn)= " << res/barn
702  << " " << p->GetParticleName()
703  << " Z= " << Z << " shellIdx= " << shellIdx
704  << " cut(keV)= " << aCut/keV
705  << G4endl;
706  }
707  }
708  }
709  return res;
710 }
711 
712 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
713 
714 G4double
715 G4EmCalculator::ComputeGammaAttenuationLength(G4double kinEnergy,
716  const G4Material* mat)
717 {
718  G4double res = 0.0;
719  const G4ParticleDefinition* gamma = G4Gamma::Gamma();
720  res += ComputeCrossSectionPerVolume(kinEnergy, gamma, "conv", mat, 0.0);
721  res += ComputeCrossSectionPerVolume(kinEnergy, gamma, "compt", mat, 0.0);
722  res += ComputeCrossSectionPerVolume(kinEnergy, gamma, "phot", mat, 0.0);
723  res += ComputeCrossSectionPerVolume(kinEnergy, gamma, "Rayl", mat, 0.0);
724  if(res > 0.0) { res = 1.0/res; }
725  return res;
726 }
727 
728 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
729 
730 G4double G4EmCalculator::ComputeShellIonisationCrossSectionPerAtom(
731  const G4String& particle,
732  G4int Z,
733  G4AtomicShellEnumerator shell,
734  G4double kinEnergy,
735  const G4Material* mat)
736 {
737  G4double res = 0.0;
738  const G4ParticleDefinition* p = FindParticle(particle);
739  G4VAtomDeexcitation* ad = manager->AtomDeexcitation();
740  if(p && ad) {
741  res = ad->ComputeShellIonisationCrossSectionPerAtom(p, Z, shell,
742  kinEnergy, mat);
743  }
744  return res;
745 }
746 
747 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
748 
749 G4double G4EmCalculator::ComputeMeanFreePath(G4double kinEnergy,
750  const G4ParticleDefinition* p,
751  const G4String& processName,
752  const G4Material* mat,
753  G4double cut)
754 {
755  G4double mfp = DBL_MAX;
756  G4double x = ComputeCrossSectionPerVolume(kinEnergy, p, processName, mat, cut);
757  if(x > 0.0) { mfp = 1.0/x; }
758  if(verbose>1) {
759  G4cout << "E(MeV)= " << kinEnergy/MeV
760  << " MFP(mm)= " << mfp/mm
761  << " " << p->GetParticleName()
762  << " in " << mat->GetName()
763  << G4endl;
764  }
765  return mfp;
766 }
767 
768 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
769 
770 G4double G4EmCalculator::ComputeEnergyCutFromRangeCut(
771  G4double range,
772  const G4ParticleDefinition* part,
773  const G4Material* mat)
774 {
775  return G4ProductionCutsTable::GetProductionCutsTable()->
776  ConvertRangeToEnergy(part, mat, range);
777 }
778 
779 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
780 
781 G4bool G4EmCalculator::UpdateParticle(const G4ParticleDefinition* p,
782  G4double kinEnergy)
783 {
784  if(p != currentParticle) {
785 
786  // new particle
787  currentParticle = p;
788  dynParticle.SetDefinition(const_cast<G4ParticleDefinition*>(p));
789  dynParticle.SetKineticEnergy(kinEnergy);
790  baseParticle = 0;
791  currentParticleName = p->GetParticleName();
792  massRatio = 1.0;
793  mass = p->GetPDGMass();
794  chargeSquare = 1.0;
795  currentProcess = FindEnergyLossProcess(p);
796  currentProcessName = "";
797  isIon = false;
798 
799  // ionisation process exist
800  if(currentProcess) {
801  currentProcessName = currentProcess->GetProcessName();
802  baseParticle = currentProcess->BaseParticle();
803 
804  // base particle is used
805  if(baseParticle) {
806  massRatio = baseParticle->GetPDGMass()/p->GetPDGMass();
807  G4double q = p->GetPDGCharge()/baseParticle->GetPDGCharge();
808  chargeSquare = q*q;
809  }
810 
811  if(p->GetParticleType() == "nucleus"
812  && currentParticleName != "deuteron"
813  && currentParticleName != "triton"
814  && currentParticleName != "alpha+"
815  && currentParticleName != "helium"
816  && currentParticleName != "hydrogen"
817  ) {
818  isIon = true;
819  massRatio = theGenericIon->GetPDGMass()/p->GetPDGMass();
820  baseParticle = theGenericIon;
821  if(verbose>1) {
822  G4cout << "\n G4EmCalculator::UpdateParticle: isIon 1 "
823  << p->GetParticleName()
824  << " in " << currentMaterial->GetName()
825  << " e= " << kinEnergy << G4endl;
826  }
827  }
828  }
829  }
830 
831  // Effective charge for ions
832  if(isIon) {
833  chargeSquare =
834  corr->EffectiveChargeSquareRatio(p, currentMaterial, kinEnergy)
835  * corr->EffectiveChargeCorrection(p,currentMaterial,kinEnergy);
836  if(currentProcess) {
837  currentProcess->SetDynamicMassCharge(massRatio,chargeSquare);
838  if(verbose>1) {
839  G4cout <<"\n NewIon: massR= "<< massRatio << " q2= "
840  << chargeSquare << " " << currentProcess << G4endl;
841  }
842  }
843  }
844  return true;
845 }
846 
847 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
848 
849 const G4ParticleDefinition* G4EmCalculator::FindParticle(const G4String& name)
850 {
851  const G4ParticleDefinition* p = 0;
852  if(name != currentParticleName) {
853  p = G4ParticleTable::GetParticleTable()->FindParticle(name);
854  if(!p) {
855  G4cout << "### WARNING: G4EmCalculator::FindParticle fails to find "
856  << name << G4endl;
857  }
858  } else {
859  p = currentParticle;
860  }
861  return p;
862 }
863 
864 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
865 
866 const G4ParticleDefinition* G4EmCalculator::FindIon(G4int Z, G4int A)
867 {
868  const G4ParticleDefinition* p = ionTable->GetIon(Z,A,0);
869  return p;
870 }
871 
872 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
873 
874 const G4Material* G4EmCalculator::FindMaterial(const G4String& name)
875 {
876  if(name != currentMaterialName) {
877  SetupMaterial(G4Material::GetMaterial(name, false));
878  if(!currentMaterial) {
879  G4cout << "### WARNING: G4EmCalculator::FindMaterial fails to find "
880  << name << G4endl;
881  }
882  }
883  return currentMaterial;
884 }
885 
886 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
887 
888 const G4Region* G4EmCalculator::FindRegion(const G4String& reg)
889 {
890  const G4Region* r = 0;
891  if(reg != "" && reg != "world") {
892  r = G4RegionStore::GetInstance()->GetRegion(reg);
893  } else {
894  r = G4RegionStore::GetInstance()->GetRegion("DefaultRegionForTheWorld");
895  }
896  return r;
897 }
898 
899 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
900 
901 const G4MaterialCutsCouple* G4EmCalculator::FindCouple(
902  const G4Material* material,
903  const G4Region* region)
904 {
905  const G4MaterialCutsCouple* couple = nullptr;
906  SetupMaterial(material);
907  if(currentMaterial) {
908  // Access to materials
909  const G4ProductionCutsTable* theCoupleTable=
910  G4ProductionCutsTable::GetProductionCutsTable();
911  const G4Region* r = region;
912  if(r) {
913  couple = theCoupleTable->GetMaterialCutsCouple(material,
914  r->GetProductionCuts());
915  } else {
916  G4RegionStore* store = G4RegionStore::GetInstance();
917  size_t nr = store->size();
918  if(0 < nr) {
919  for(size_t i=0; i<nr; ++i) {
920  couple = theCoupleTable->GetMaterialCutsCouple(
921  material, ((*store)[i])->GetProductionCuts());
922  if(couple) { break; }
923  }
924  }
925  }
926  }
927  if(!couple) {
928  G4ExceptionDescription ed;
929  ed << "G4EmCalculator::FindCouple: fail for material <"
930  << currentMaterialName << ">";
931  if(region) { ed << " and region " << region->GetName(); }
932  G4Exception("G4EmCalculator::FindCouple", "em0078",
933  FatalException, ed);
934  }
935  return couple;
936 }
937 
938 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
939 
940 G4bool G4EmCalculator::UpdateCouple(const G4Material* material, G4double cut)
941 {
942  SetupMaterial(material);
943  if(!currentMaterial) { return false; }
944  for (G4int i=0; i<nLocalMaterials; ++i) {
945  if(material == localMaterials[i] && cut == localCuts[i]) {
946  currentCouple = localCouples[i];
947  currentCoupleIndex = currentCouple->GetIndex();
948  currentCut = cut;
949  return true;
950  }
951  }
952  const G4MaterialCutsCouple* cc = new G4MaterialCutsCouple(material);
953  localMaterials.push_back(material);
954  localCouples.push_back(cc);
955  localCuts.push_back(cut);
956  nLocalMaterials++;
957  currentCouple = cc;
958  currentCoupleIndex = currentCouple->GetIndex();
959  currentCut = cut;
960  return true;
961 }
962 
963 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
964 
965 void G4EmCalculator::FindLambdaTable(const G4ParticleDefinition* p,
966  const G4String& processName,
967  G4double kinEnergy, G4int& proctype)
968 {
969  // Search for the process
970  if (!currentLambda || p != lambdaParticle || processName != lambdaName) {
971  lambdaName = processName;
972  currentLambda = nullptr;
973  lambdaParticle = p;
974 
975  const G4ParticleDefinition* part = p;
976  if(isIon) { part = theGenericIon; }
977 
978  // Search for energy loss process
979  currentName = processName;
980  currentModel = nullptr;
981  loweModel = nullptr;
982 
983  G4VEnergyLossProcess* elproc = FindEnLossProcess(part, processName);
984  if(elproc) {
985  currentLambda = elproc->LambdaTable();
986  proctype = 0;
987  if(currentLambda) {
988  isApplicable = true;
989  if(verbose>1) {
990  G4cout << "G4VEnergyLossProcess is found out: " << currentName
991  << G4endl;
992  }
993  }
994  curProcess = elproc;
995  return;
996  }
997 
998  // Search for discrete process
999  G4VEmProcess* proc = FindDiscreteProcess(part, processName);
1000  if(proc) {
1001  currentLambda = proc->LambdaTable();
1002  proctype = 1;
1003  if(currentLambda) {
1004  isApplicable = true;
1005  if(verbose>1) {
1006  G4cout << "G4VEmProcess is found out: " << currentName << G4endl;
1007  }
1008  }
1009  curProcess = proc;
1010  return;
1011  }
1012 
1013  // Search for msc process
1014  G4VMultipleScattering* msc = FindMscProcess(part, processName);
1015  if(msc) {
1016  currentModel = msc->SelectModel(kinEnergy,0);
1017  proctype = 2;
1018  if(currentModel) {
1019  currentLambda = currentModel->GetCrossSectionTable();
1020  if(currentLambda) {
1021  isApplicable = true;
1022  if(verbose>1) {
1023  G4cout << "G4VMultipleScattering is found out: " << currentName
1024  << G4endl;
1025  }
1026  }
1027  }
1028  curProcess = msc;
1029  }
1030  }
1031 }
1032 
1033 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1034 
1035 G4bool G4EmCalculator::FindEmModel(const G4ParticleDefinition* p,
1036  const G4String& processName,
1037  G4double kinEnergy)
1038 {
1039  isApplicable = false;
1040  if(!p || !currentMaterial) {
1041  G4cout << "G4EmCalculator::FindEmModel WARNING: no particle"
1042  << " or materail defined; particle: " << p << G4endl;
1043  return isApplicable;
1044  }
1045  G4String partname = p->GetParticleName();
1046  const G4ParticleDefinition* part = p;
1047  G4double scaledEnergy = kinEnergy*massRatio;
1048  if(isIon) { part = theGenericIon; }
1049 
1050  if(verbose > 1) {
1051  G4cout << "## G4EmCalculator::FindEmModel for " << partname
1052  << " (type= " << p->GetParticleType()
1053  << ") and " << processName << " at E(MeV)= " << scaledEnergy
1054  << G4endl;
1055  if(p != part) { G4cout << " GenericIon is the base particle" << G4endl; }
1056  }
1057 
1058  // Search for energy loss process
1059  currentName = processName;
1060  currentModel = nullptr;
1061  loweModel = nullptr;
1062  size_t idx = 0;
1063 
1064  G4VEnergyLossProcess* elproc = FindEnLossProcess(part, processName);
1065  if(elproc) {
1066  currentModel = elproc->SelectModelForMaterial(scaledEnergy, idx);
1067  currentModel->InitialiseForMaterial(part, currentMaterial);
1068  currentModel->SetupForMaterial(part, currentMaterial, scaledEnergy);
1069  G4double eth = currentModel->LowEnergyLimit();
1070  if(eth > 0.0) {
1071  loweModel = elproc->SelectModelForMaterial(eth - CLHEP::eV, idx);
1072  if(loweModel == currentModel) { loweModel = nullptr; }
1073  else {
1074  loweModel->InitialiseForMaterial(part, currentMaterial);
1075  loweModel->SetupForMaterial(part, currentMaterial, eth - CLHEP::eV);
1076  }
1077  }
1078  }
1079 
1080  // Search for discrete process
1081  if(!currentModel) {
1082  G4VEmProcess* proc = FindDiscreteProcess(part, processName);
1083  if(proc) {
1084  currentModel = proc->SelectModelForMaterial(kinEnergy, idx);
1085  currentModel->InitialiseForMaterial(part, currentMaterial);
1086  currentModel->SetupForMaterial(part, currentMaterial, kinEnergy);
1087  G4double eth = currentModel->LowEnergyLimit();
1088  if(eth > 0.0) {
1089  loweModel = proc->SelectModelForMaterial(eth - CLHEP::eV, idx);
1090  if(loweModel == currentModel) { loweModel = nullptr; }
1091  else {
1092  loweModel->InitialiseForMaterial(part, currentMaterial);
1093  loweModel->SetupForMaterial(part, currentMaterial, eth - CLHEP::eV);
1094  }
1095  }
1096  }
1097  }
1098 
1099  // Search for msc process
1100  if(!currentModel) {
1101  G4VMultipleScattering* proc = FindMscProcess(part, processName);
1102  if(proc) {
1103  currentModel = proc->SelectModel(kinEnergy, idx);
1104  loweModel = nullptr;
1105  }
1106  }
1107  if(currentModel) {
1108  if(loweModel == currentModel) { loweModel = nullptr; }
1109  isApplicable = true;
1110  currentModel->InitialiseForMaterial(part, currentMaterial);
1111  if(loweModel) {
1112  loweModel->InitialiseForMaterial(part, currentMaterial);
1113  }
1114  if(verbose > 1) {
1115  G4cout << " Model <" << currentModel->GetName()
1116  << "> Emin(MeV)= " << currentModel->LowEnergyLimit()/MeV
1117  << " for " << part->GetParticleName();
1118  if(elproc) {
1119  G4cout << " and " << elproc->GetProcessName() << " " << elproc
1120  << G4endl;
1121  }
1122  if(loweModel) {
1123  G4cout << " LowEnergy model <" << loweModel->GetName() << ">";
1124  }
1125  G4cout << G4endl;
1126  }
1127  }
1128  return isApplicable;
1129 }
1130 
1131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1132 
1133 G4VEnergyLossProcess* G4EmCalculator::FindEnergyLossProcess(
1134  const G4ParticleDefinition* p)
1135 {
1136  G4VEnergyLossProcess* elp = nullptr;
1137  G4String partname = p->GetParticleName();
1138  const G4ParticleDefinition* part = p;
1139 
1140  if(p->GetParticleType() == "nucleus"
1141  && currentParticleName != "deuteron"
1142  && currentParticleName != "triton"
1143  && currentParticleName != "He3"
1144  && currentParticleName != "alpha"
1145  && currentParticleName != "alpha+"
1146  && currentParticleName != "helium"
1147  && currentParticleName != "hydrogen"
1148  ) { part = theGenericIon; }
1149 
1150  elp = manager->GetEnergyLossProcess(part);
1151  /*
1152  G4cout << "\n G4EmCalculator::FindEnergyLossProcess: for " << p->GetParticleName()
1153  << " found " << elp->GetProcessName() << " of "
1154  << elp->Particle()->GetParticleName() << " " << elp << G4endl;
1155  */
1156  return elp;
1157 }
1158 
1159 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1160 
1161 G4VEnergyLossProcess*
1162 G4EmCalculator::FindEnLossProcess(const G4ParticleDefinition* part,
1163  const G4String& processName)
1164 {
1165  G4VEnergyLossProcess* proc = 0;
1166  const std::vector<G4VEnergyLossProcess*> v =
1167  manager->GetEnergyLossProcessVector();
1168  G4int n = v.size();
1169  for(G4int i=0; i<n; ++i) {
1170  if((v[i])->GetProcessName() == processName) {
1171  G4VProcess* p = reinterpret_cast<G4VProcess*>(v[i]);
1172  if(ActiveForParticle(part, p)) {
1173  proc = v[i];
1174  break;
1175  }
1176  }
1177  }
1178  return proc;
1179 }
1180 
1181 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1182 
1183 G4VEmProcess*
1184 G4EmCalculator::FindDiscreteProcess(const G4ParticleDefinition* part,
1185  const G4String& processName)
1186 {
1187  G4VEmProcess* proc = nullptr;
1188  auto v = manager->GetEmProcessVector();
1189  G4int n = v.size();
1190  for(G4int i=0; i<n; ++i) {
1191  auto pName = v[i]->GetProcessName();
1192  if(pName == "GammaGeneralProc") {
1193  proc = v[i]->GetEmProcess(processName);
1194  break;
1195  } else if(pName == processName) {
1196  G4VProcess* p = reinterpret_cast<G4VProcess*>(v[i]);
1197  if(ActiveForParticle(part, p)) {
1198  proc = v[i];
1199  break;
1200  }
1201  }
1202  }
1203  return proc;
1204 }
1205 
1206 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1207 
1208 G4VMultipleScattering*
1209 G4EmCalculator::FindMscProcess(const G4ParticleDefinition* part,
1210  const G4String& processName)
1211 {
1212  G4VMultipleScattering* proc = 0;
1213  const std::vector<G4VMultipleScattering*> v =
1214  manager->GetMultipleScatteringVector();
1215  G4int n = v.size();
1216  for(G4int i=0; i<n; ++i) {
1217  if((v[i])->GetProcessName() == processName) {
1218  G4VProcess* p = reinterpret_cast<G4VProcess*>(v[i]);
1219  if(ActiveForParticle(part, p)) {
1220  proc = v[i];
1221  break;
1222  }
1223  }
1224  }
1225  return proc;
1226 }
1227 
1228 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1229 
1230 G4VProcess* G4EmCalculator::FindProcess(const G4ParticleDefinition* part,
1231  const G4String& processName)
1232 {
1233  G4VProcess* proc = 0;
1234  const G4ProcessManager* procman = part->GetProcessManager();
1235  G4ProcessVector* pv = procman->GetProcessList();
1236  G4int nproc = pv->size();
1237  for(G4int i=0; i<nproc; ++i) {
1238  if(processName == (*pv)[i]->GetProcessName()) {
1239  proc = (*pv)[i];
1240  break;
1241  }
1242  }
1243  return proc;
1244 }
1245 
1246 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1247 
1248 G4bool G4EmCalculator::ActiveForParticle(const G4ParticleDefinition* part,
1249  G4VProcess* proc)
1250 {
1251  G4ProcessManager* pm = part->GetProcessManager();
1252  G4ProcessVector* pv = pm->GetProcessList();
1253  G4int n = pv->size();
1254  G4bool res = false;
1255  for(G4int i=0; i<n; ++i) {
1256  if((*pv)[i] == proc) {
1257  if(pm->GetProcessActivation(i)) { res = true; }
1258  break;
1259  }
1260  }
1261  return res;
1262 }
1263 
1264 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1265 
1266 void G4EmCalculator::SetupMaterial(const G4Material* mat)
1267 {
1268  if(mat) {
1269  currentMaterial = mat;
1270  currentMaterialName = mat->GetName();
1271  } else {
1272  currentMaterial = 0;
1273  currentMaterialName = "";
1274  }
1275 }
1276 
1277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1278 
1279 void G4EmCalculator::SetupMaterial(const G4String& mname)
1280 {
1281  SetupMaterial(nist->FindOrBuildMaterial(mname));
1282 }
1283 
1284 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1285 
1286 void G4EmCalculator::CheckMaterial(G4int Z)
1287 {
1288  G4bool isFound = false;
1289  if(currentMaterial) {
1290  size_t nn = currentMaterial->GetNumberOfElements();
1291  for(size_t i=0; i<nn; ++i) {
1292  if(Z == currentMaterial->GetElement(i)->GetZasInt()) {
1293  isFound = true;
1294  break;
1295  }
1296  }
1297  }
1298  if(!isFound) {
1299  SetupMaterial(nist->FindOrBuildSimpleMaterial(Z));
1300  }
1301 }
1302 
1303 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1304 
1305 void G4EmCalculator::SetVerbose(G4int verb)
1306 {
1307  verbose = verb;
1308 }
1309 
1310 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1311