Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/utils/src/G4EmCalculator.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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