Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4MicroElecInelasticModel.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // G4MicroElecInelasticModel.cc, 2011/08/29 A.Valentin, M. Raine
28 //
29 // Based on the following publications
30 //
31 // - Inelastic cross-sections of low energy electrons in silicon
32 // for the simulation of heavy ion tracks with theGeant4-DNA toolkit,
33 // NSS Conf. Record 2010, pp. 80-85.
34 // - Geant4 physics processes for microdosimetry simulation:
35 // very low energy electromagnetic models for electrons in Si,
36 // NIM B, vol. 288, pp. 66 - 73, 2012.
37 // - Geant4 physics processes for microdosimetry simulation:
38 // very low energy electromagnetic models for protons and
39 // heavy ions in Si, NIM B, vol. 287, pp. 124 - 129, 2012.
40 //
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
42 
43 #include "G4MicroElecInelasticModel.hh"
44 
45 #include "globals.hh"
46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "G4ios.hh"
49 #include "G4UnitsTable.hh"
50 #include "G4UAtomicDeexcitation.hh"
51 #include "G4LossTableManager.hh"
52 #include "G4ionEffectiveCharge.hh"
53 
54 #include "G4DeltaAngle.hh"
55 
56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
57 
58 using namespace std;
59 
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
61 
62 G4MicroElecInelasticModel::G4MicroElecInelasticModel(const G4ParticleDefinition*,
63  const G4String& nam)
64 :G4VEmModel(nam),isInitialised(false)
65 {
66  nistSi = G4NistManager::Instance()->FindOrBuildMaterial("G4_Si");
67 
68  verboseLevel= 0;
69  // Verbosity scale:
70  // 0 = nothing
71  // 1 = warning for energy non-conservation
72  // 2 = details of energy budget
73  // 3 = calculation of cross sections, file openings, sampling of atoms
74  // 4 = entering in methods
75 
76  if( verboseLevel>0 )
77  {
78  G4cout << "MicroElec inelastic model is constructed " << G4endl;
79  }
80 
81  //Mark this model as "applicable" for atomic deexcitation
82  SetDeexcitationFlag(true);
83  fAtomDeexcitation = 0;
84  fParticleChangeForGamma = 0;
85 
86  // default generator
87  SetAngularDistribution(new G4DeltaAngle());
88 
89  // Selection of computation method
90  fasterCode = true; //false;
91 }
92 
93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
94 
95 G4MicroElecInelasticModel::~G4MicroElecInelasticModel()
96 {
97  // Cross section
98 
99  std::map< G4String,G4MicroElecCrossSectionDataSet*,std::less<G4String> >::iterator pos;
100  for (pos = tableData.begin(); pos != tableData.end(); ++pos)
101  {
102  G4MicroElecCrossSectionDataSet* table = pos->second;
103  delete table;
104  }
105 
106  // Final state
107 
108  eVecm.clear();
109  pVecm.clear();
110 
111 }
112 
113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
114 
115 void G4MicroElecInelasticModel::Initialise(const G4ParticleDefinition* particle,
116  const G4DataVector& /*cuts*/)
117 {
118 
119  if (verboseLevel > 3)
120  G4cout << "Calling G4MicroElecInelasticModel::Initialise()" << G4endl;
121 
122  // Energy limits
123 
124  G4String fileElectron("microelec/sigma_inelastic_e_Si");
125  G4String fileProton("microelec/sigma_inelastic_p_Si");
126 
127  G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
128  G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
129 
130  G4String electron;
131  G4String proton;
132 
133  G4double scaleFactor = 1e-18 * cm *cm;
134 
135  char *path = getenv("G4LEDATA");
136 
137  // *** ELECTRON
138  electron = electronDef->GetParticleName();
139 
140  tableFile[electron] = fileElectron;
141 
142  lowEnergyLimit[electron] = 16.7 * eV;
143  highEnergyLimit[electron] = 100.0 * MeV;
144 
145  // Cross section
146 
147  G4MicroElecCrossSectionDataSet* tableE = new G4MicroElecCrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
148  tableE->LoadData(fileElectron);
149 
150  tableData[electron] = tableE;
151 
152  // Final state
153 
154  std::ostringstream eFullFileName;
155 
156  if (fasterCode) eFullFileName << path << "/microelec/sigmadiff_cumulated_inelastic_e_Si.dat";
157  else eFullFileName << path << "/microelec/sigmadiff_inelastic_e_Si.dat";
158 
159  std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
160 
161  if (!eDiffCrossSection)
162  {
163  if (fasterCode) G4Exception("G4MicroElecInelasticModel::Initialise","em0003",
164  FatalException,"Missing data file:/microelec/sigmadiff_cumulated_inelastic_e_Si.dat");
165 
166  else G4Exception("G4MicroElecInelasticModel::Initialise","em0003",
167  FatalException,"Missing data file:/microelec/sigmadiff_inelastic_e_Si.dat");
168 
169  }
170 
171  //
172 
173  // Clear the arrays for re-initialization case (MT mode)
174  // Octobre 22nd, 2014 - Melanie Raine
175 
176  eTdummyVec.clear();
177  pTdummyVec.clear();
178 
179  eVecm.clear();
180  pVecm.clear();
181 
182  for (int j=0; j<6; j++)
183  {
184  eProbaShellMap[j].clear();
185  pProbaShellMap[j].clear();
186 
187  eDiffCrossSectionData[j].clear();
188  pDiffCrossSectionData[j].clear();
189 
190  eNrjTransfData[j].clear();
191  pNrjTransfData[j].clear();
192  }
193 
194  //
195 
196 
197  eTdummyVec.push_back(0.);
198  while(!eDiffCrossSection.eof())
199  {
200  double tDummy;
201  double eDummy;
202  eDiffCrossSection>>tDummy>>eDummy;
203  if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
204 
205  double tmp;
206  for (int j=0; j<6; j++)
207  {
208  eDiffCrossSection>> tmp;
209 
210  eDiffCrossSectionData[j][tDummy][eDummy] = tmp;
211 
212  if (fasterCode)
213  {
214  eNrjTransfData[j][tDummy][eDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
215  eProbaShellMap[j][tDummy].push_back(eDiffCrossSectionData[j][tDummy][eDummy]);
216  }
217  else
218  {
219  // SI - only if eof is not reached !
220  if (!eDiffCrossSection.eof()) eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
221  eVecm[tDummy].push_back(eDummy);
222  }
223 
224  }
225  }
226  //
227 
228  // *** PROTON
229 
230  proton = protonDef->GetParticleName();
231 
232  tableFile[proton] = fileProton;
233 
234  lowEnergyLimit[proton] = 50. * keV;
235  highEnergyLimit[proton] = 10. * GeV;
236 
237  // Cross section
238 
239  G4MicroElecCrossSectionDataSet* tableP = new G4MicroElecCrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
240  tableP->LoadData(fileProton);
241 
242  tableData[proton] = tableP;
243 
244  // Final state
245 
246  std::ostringstream pFullFileName;
247 
248  if (fasterCode) pFullFileName << path << "/microelec/sigmadiff_cumulated_inelastic_p_Si.dat";
249  else pFullFileName << path << "/microelec/sigmadiff_inelastic_p_Si.dat";
250 
251  std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
252 
253  if (!pDiffCrossSection)
254  {
255  if (fasterCode) G4Exception("G4MicroElecInelasticModel::Initialise","em0003",
256  FatalException,"Missing data file:/microelec/sigmadiff_cumulated_inelastic_p_Si.dat");
257 
258  else G4Exception("G4MicroElecInelasticModel::Initialise","em0003",
259  FatalException,"Missing data file:/microelec/sigmadiff_inelastic_p_Si.dat");
260  }
261 
262  pTdummyVec.push_back(0.);
263  while(!pDiffCrossSection.eof())
264  {
265  double tDummy;
266  double eDummy;
267  pDiffCrossSection>>tDummy>>eDummy;
268  if (tDummy != pTdummyVec.back()) pTdummyVec.push_back(tDummy);
269  for (int j=0; j<6; j++)
270  {
271  pDiffCrossSection>>pDiffCrossSectionData[j][tDummy][eDummy];
272 
273  if (fasterCode)
274  {
275  pNrjTransfData[j][tDummy][pDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
276  pProbaShellMap[j][tDummy].push_back(pDiffCrossSectionData[j][tDummy][eDummy]);
277  }
278  else
279  {
280  // SI - only if eof is not reached !
281  if (!pDiffCrossSection.eof()) pDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
282  pVecm[tDummy].push_back(eDummy);
283  }
284  }
285  }
286 
287 
288  if (particle==electronDef)
289  {
290  SetLowEnergyLimit(lowEnergyLimit[electron]);
291  SetHighEnergyLimit(highEnergyLimit[electron]);
292  }
293 
294  if (particle==protonDef)
295  {
296  SetLowEnergyLimit(lowEnergyLimit[proton]);
297  SetHighEnergyLimit(highEnergyLimit[proton]);
298  }
299 
300  if( verboseLevel>0 )
301  {
302  G4cout << "MicroElec Inelastic model is initialized " << G4endl
303  << "Energy range: "
304  << LowEnergyLimit() / keV << " keV - "
305  << HighEnergyLimit() / MeV << " MeV for "
306  << particle->GetParticleName()
307  << " with mass (amu) " << particle->GetPDGMass()/proton_mass_c2
308  << " and charge " << particle->GetPDGCharge()
309  << G4endl << G4endl ;
310  }
311 
312  //
313 
314  fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
315 
316  if (isInitialised) { return; }
317  fParticleChangeForGamma = GetParticleChangeForGamma();
318  isInitialised = true;
319 
320 }
321 
322 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
323 
324 G4double G4MicroElecInelasticModel::CrossSectionPerVolume(const G4Material* material,
325  const G4ParticleDefinition* particleDefinition,
326  G4double ekin,
327  G4double,
328  G4double)
329 {
330  if (verboseLevel > 3)
331  G4cout << "Calling CrossSectionPerVolume() of G4MicroElecInelasticModel" << G4endl;
332 
333  G4double density = material->GetTotNbOfAtomsPerVolume();
334 
335  /* if (
336  particleDefinition != G4Proton::ProtonDefinition()
337  &&
338  particleDefinition != G4Electron::ElectronDefinition()
339  &&
340  particleDefinition != G4GenericIon::GenericIonDefinition()
341  )
342 
343  return 0;*/
344 
345  // Calculate total cross section for model
346 
347  G4double lowLim = 0;
348  G4double highLim = 0;
349  G4double sigma=0;
350 
351  const G4String& particleName = particleDefinition->GetParticleName();
352  G4String nameLocal = particleName ;
353 
354  G4double Zeff2 = 1.0;
355  G4double Mion_c2 = particleDefinition->GetPDGMass();
356 
357  if (Mion_c2 > proton_mass_c2)
358  {
359  G4ionEffectiveCharge EffCharge ;
360  G4double Zeff = EffCharge.EffectiveCharge(particleDefinition, material,ekin);
361  Zeff2 = Zeff*Zeff;
362 
363  if (verboseLevel > 3)
364  G4cout << "Before scaling : " << G4endl
365  << "Particle : " << nameLocal << ", mass : " << Mion_c2/proton_mass_c2 << "*mp, charge " << Zeff
366  << ", Ekin (eV) = " << ekin/eV << G4endl ;
367 
368  ekin *= proton_mass_c2/Mion_c2 ;
369  nameLocal = "proton" ;
370 
371  if (verboseLevel > 3)
372  G4cout << "After scaling : " << G4endl
373  << "Particle : " << nameLocal << ", Ekin (eV) = " << ekin/eV << G4endl ;
374  }
375 
376  if (material == nistSi || material->GetBaseMaterial() == nistSi)
377  {
378 
379  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
380  pos1 = lowEnergyLimit.find(nameLocal);
381  if (pos1 != lowEnergyLimit.end())
382  {
383  lowLim = pos1->second;
384  }
385 
386  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
387  pos2 = highEnergyLimit.find(nameLocal);
388  if (pos2 != highEnergyLimit.end())
389  {
390  highLim = pos2->second;
391  }
392 
393  if (ekin >= lowLim && ekin < highLim)
394  {
395  std::map< G4String,G4MicroElecCrossSectionDataSet*,std::less<G4String> >::iterator pos;
396  pos = tableData.find(nameLocal);
397 
398  if (pos != tableData.end())
399  {
400  G4MicroElecCrossSectionDataSet* table = pos->second;
401  if (table != 0)
402  {
403  sigma = table->FindValue(ekin);
404  }
405  }
406  else
407  {
408  G4Exception("G4MicroElecInelasticModel::CrossSectionPerVolume","em0002",FatalException,"Model not applicable to particle type.");
409  }
410  }
411  else
412  {
413  if (nameLocal!="e-")
414  {
415  // G4cout << "Particle : " << nameLocal << ", Ekin (eV) = " << ekin/eV << G4endl;
416  // G4cout << "### Warning: particle energy out of bounds! ###" << G4endl;
417  }
418  }
419 
420  if (verboseLevel > 3)
421  {
422  G4cout << "---> Kinetic energy (eV)=" << ekin/eV << G4endl;
423  G4cout << " - Cross section per Si atom (cm^2)=" << sigma*Zeff2/cm2 << G4endl;
424  G4cout << " - Cross section per Si atom (cm^-1)=" << sigma*density*Zeff2/(1./cm) << G4endl;
425  }
426 
427  } // if (SiMaterial)
428  return sigma*density*Zeff2;
429 
430 
431 }
432 
433 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
434 
435 void G4MicroElecInelasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
436  const G4MaterialCutsCouple* couple,
437  const G4DynamicParticle* particle,
438  G4double,
439  G4double)
440 {
441 
442  if (verboseLevel > 3)
443  G4cout << "Calling SampleSecondaries() of G4MicroElecInelasticModel" << G4endl;
444 
445  G4double lowLim = 0;
446  G4double highLim = 0;
447 
448  G4double ekin = particle->GetKineticEnergy();
449  G4double k = ekin ;
450 
451  G4ParticleDefinition* PartDef = particle->GetDefinition();
452  const G4String& particleName = PartDef->GetParticleName();
453  G4String nameLocal2 = particleName ;
454  G4double particleMass = particle->GetDefinition()->GetPDGMass();
455 
456  if (particleMass > proton_mass_c2)
457  {
458  k *= proton_mass_c2/particleMass ;
459  PartDef = G4Proton::ProtonDefinition();
460  nameLocal2 = "proton" ;
461  }
462 
463  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
464  pos1 = lowEnergyLimit.find(nameLocal2);
465 
466  if (pos1 != lowEnergyLimit.end())
467  {
468  lowLim = pos1->second;
469  }
470 
471  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
472  pos2 = highEnergyLimit.find(nameLocal2);
473 
474  if (pos2 != highEnergyLimit.end())
475  {
476  highLim = pos2->second;
477  }
478 
479  if (k >= lowLim && k < highLim)
480  {
481  G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
482  G4double totalEnergy = ekin + particleMass;
483  G4double pSquare = ekin * (totalEnergy + particleMass);
484  G4double totalMomentum = std::sqrt(pSquare);
485 
486  G4int Shell = 0;
487 
488  /* if (!fasterCode)*/ Shell = RandomSelect(k,nameLocal2);
489 
490  // SI: The following protection is necessary to avoid infinite loops :
491  // sigmadiff_ionisation_e_born.dat has non zero partial xs at 18 eV for shell 3 (ionizationShell ==2)
492  // sigmadiff_cumulated_ionisation_e_born.dat has zero cumulated partial xs at 18 eV for shell 3 (ionizationShell ==2)
493  // this is due to the fact that the max allowed transfered energy is (18+10.79)/2=17.025 eV and only transfered energies
494  // strictly above this value have non zero partial xs in sigmadiff_ionisation_e_born.dat (starting at trans = 17.12 eV)
495 
496  /*if (fasterCode)
497  do
498  {
499  Shell = RandomSelect(k,nameLocal2);
500  }while (k<19*eV && ionizationShell==2 && particle->GetDefinition()==G4Electron::ElectronDefinition());*/
501 
502  G4double bindingEnergy = SiStructure.Energy(Shell);
503 
504  if (verboseLevel > 3)
505  {
506  G4cout << "---> Kinetic energy (eV)=" << k/eV << G4endl ;
507  G4cout << "Shell: " << Shell << ", energy: " << bindingEnergy/eV << G4endl;
508  }
509 
510  // sample deexcitation
511 
512  G4int secNumberInit = 0; // need to know at a certain point the energy of secondaries
513  G4int secNumberFinal = 0; // So I'll make the difference and then sum the energies
514 
515  //SI: additional protection if tcs interpolation method is modified
516  if (k<bindingEnergy) return;
517 
518  G4int Z = 14;
519 
520  if(fAtomDeexcitation && Shell > 2) {
521 
522  G4AtomicShellEnumerator as = fKShell;
523 
524  if (Shell == 4)
525  {
526  as = G4AtomicShellEnumerator(1);
527  }
528  else if (Shell == 3)
529  {
530  as = G4AtomicShellEnumerator(3);
531  }
532 
533  const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
534  secNumberInit = fvect->size();
535  fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
536  secNumberFinal = fvect->size();
537  }
538 
539  G4double secondaryKinetic=-1000*eV;
540 
541  if (!fasterCode)
542  {
543  secondaryKinetic = RandomizeEjectedElectronEnergy(PartDef,k,Shell);
544  }
545  else
546  {
547  secondaryKinetic = RandomizeEjectedElectronEnergyFromCumulatedDcs(PartDef,k,Shell);
548  }
549 
550 
551  if (verboseLevel > 3)
552  {
553  G4cout << "Ionisation process" << G4endl;
554  G4cout << "Shell: " << Shell << " Kin. energy (eV)=" << k/eV
555  << " Sec. energy (eV)=" << secondaryKinetic/eV << G4endl;
556  }
557 
558  G4ThreeVector deltaDirection =
559  GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
560  Z, Shell,
561  couple->GetMaterial());
562 
563  if (particle->GetDefinition() == G4Electron::ElectronDefinition())
564  {
565  G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
566 
567  G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
568  G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
569  G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
570  G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
571  finalPx /= finalMomentum;
572  finalPy /= finalMomentum;
573  finalPz /= finalMomentum;
574 
575  G4ThreeVector direction;
576  direction.set(finalPx,finalPy,finalPz);
577 
578  fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
579  }
580  else fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection) ;
581 
582  // note that secondaryKinetic is the energy of the delta ray, not of all secondaries.
583  G4double deexSecEnergy = 0;
584  for (G4int j=secNumberInit; j < secNumberFinal; j++) {
585  deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();}
586 
587  fParticleChangeForGamma->SetProposedKineticEnergy(ekin-bindingEnergy-secondaryKinetic);
588  fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy-deexSecEnergy);
589 
590  if (secondaryKinetic>0)
591  {
592  G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
593  fvect->push_back(dp);
594  }
595 
596  }
597 }
598 
599 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
600 
601 G4double G4MicroElecInelasticModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
602  G4double k, G4int shell)
603 {
604  if (particleDefinition == G4Electron::ElectronDefinition())
605  {
606  G4double maximumEnergyTransfer=0.;
607  if ((k+SiStructure.Energy(shell))/2. > k) maximumEnergyTransfer=k;
608  else maximumEnergyTransfer = (k+SiStructure.Energy(shell))/2.;
609 
610  G4double crossSectionMaximum = 0.;
611 
612  G4double minEnergy = SiStructure.Energy(shell);
613  G4double maxEnergy = maximumEnergyTransfer;
614  G4int nEnergySteps = 100;
615 
616  G4double value(minEnergy);
617  G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
618  G4int step(nEnergySteps);
619  while (step>0)
620  {
621  step--;
622  G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
623  if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
624  value*=stpEnergy;
625  }
626 
627 
628  G4double secondaryElectronKineticEnergy=0.;
629  do
630  {
631  secondaryElectronKineticEnergy = G4UniformRand() * (maximumEnergyTransfer-SiStructure.Energy(shell));
632  } while(G4UniformRand()*crossSectionMaximum >
633  DifferentialCrossSection(particleDefinition, k/eV,(secondaryElectronKineticEnergy+SiStructure.Energy(shell))/eV,shell));
634 
635  return secondaryElectronKineticEnergy;
636 
637  }
638 
639  if (particleDefinition == G4Proton::ProtonDefinition())
640  {
641  G4double maximumEnergyTransfer = 4.* (electron_mass_c2 / proton_mass_c2) * k;
642  G4double crossSectionMaximum = 0.;
643 
644  G4double minEnergy = SiStructure.Energy(shell);
645  G4double maxEnergy = maximumEnergyTransfer;
646  G4int nEnergySteps = 100;
647 
648  G4double value(minEnergy);
649  G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
650  G4int step(nEnergySteps);
651  while (step>0)
652  {
653  step--;
654  G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
655  if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
656  value*=stpEnergy;
657  }
658 
659  G4double secondaryElectronKineticEnergy = 0.;
660  do
661  {
662  secondaryElectronKineticEnergy = G4UniformRand() * (maximumEnergyTransfer-SiStructure.Energy(shell));
663 
664  } while(G4UniformRand()*crossSectionMaximum >
665  DifferentialCrossSection(particleDefinition, k/eV,(secondaryElectronKineticEnergy+SiStructure.Energy(shell))/eV,shell));
666  return secondaryElectronKineticEnergy;
667  }
668 
669  return 0;
670 }
671 
672 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
673 
674 // The following section is not used anymore but is kept for memory
675 // GetAngularDistribution()->SampleDirectionForShell is used instead
676 
677 /*void G4MicroElecInelasticModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
678  G4double k,
679  G4double secKinetic,
680  G4double & cosTheta,
681  G4double & phi )
682  {
683  if (particleDefinition == G4Electron::ElectronDefinition())
684  {
685  phi = twopi * G4UniformRand();
686  G4double sin2O = (1.-secKinetic/k) / (1.+secKinetic/(2.*electron_mass_c2));
687  cosTheta = std::sqrt(1.-sin2O);
688  }
689 
690  if (particleDefinition == G4Proton::ProtonDefinition())
691  {
692  G4double maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
693  phi = twopi * G4UniformRand();
694  cosTheta = std::sqrt(secKinetic / maxSecKinetic);
695  }
696 
697  else
698  {
699  G4double maxSecKinetic = 4.* (electron_mass_c2 / particleDefinition->GetPDGMass()) * k;
700  phi = twopi * G4UniformRand();
701  cosTheta = std::sqrt(secKinetic / maxSecKinetic);
702  }
703  }
704  */
705 
706 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
707 
708 double G4MicroElecInelasticModel::DifferentialCrossSection(G4ParticleDefinition * particleDefinition,
709  G4double k,
710  G4double energyTransfer,
711  G4int LevelIndex)
712 {
713  G4double sigma = 0.;
714 
715  if (energyTransfer >= SiStructure.Energy(LevelIndex))
716  {
717  G4double valueT1 = 0;
718  G4double valueT2 = 0;
719  G4double valueE21 = 0;
720  G4double valueE22 = 0;
721  G4double valueE12 = 0;
722  G4double valueE11 = 0;
723 
724  G4double xs11 = 0;
725  G4double xs12 = 0;
726  G4double xs21 = 0;
727  G4double xs22 = 0;
728 
729  if (particleDefinition == G4Electron::ElectronDefinition())
730  {
731  // k should be in eV and energy transfer eV also
732 
733  std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
734  std::vector<double>::iterator t1 = t2-1;
735  // SI : the following condition avoids situations where energyTransfer >last vector element
736  if (energyTransfer <= eVecm[(*t1)].back() && energyTransfer <= eVecm[(*t2)].back() )
737  {
738  std::vector<double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), energyTransfer);
739  std::vector<double>::iterator e11 = e12-1;
740 
741  std::vector<double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), energyTransfer);
742  std::vector<double>::iterator e21 = e22-1;
743 
744  valueT1 =*t1;
745  valueT2 =*t2;
746  valueE21 =*e21;
747  valueE22 =*e22;
748  valueE12 =*e12;
749  valueE11 =*e11;
750 
751  xs11 = eDiffCrossSectionData[LevelIndex][valueT1][valueE11];
752  xs12 = eDiffCrossSectionData[LevelIndex][valueT1][valueE12];
753  xs21 = eDiffCrossSectionData[LevelIndex][valueT2][valueE21];
754  xs22 = eDiffCrossSectionData[LevelIndex][valueT2][valueE22];
755  }
756 
757  }
758 
759  if (particleDefinition == G4Proton::ProtonDefinition())
760  {
761  // k should be in eV and energy transfer eV also
762  std::vector<double>::iterator t2 = std::upper_bound(pTdummyVec.begin(),pTdummyVec.end(), k);
763  std::vector<double>::iterator t1 = t2-1;
764  if (energyTransfer <= pVecm[(*t1)].back() && energyTransfer <= pVecm[(*t2)].back() )
765  {
766  std::vector<double>::iterator e12 = std::upper_bound(pVecm[(*t1)].begin(),pVecm[(*t1)].end(), energyTransfer);
767  std::vector<double>::iterator e11 = e12-1;
768 
769  std::vector<double>::iterator e22 = std::upper_bound(pVecm[(*t2)].begin(),pVecm[(*t2)].end(), energyTransfer);
770  std::vector<double>::iterator e21 = e22-1;
771 
772  valueT1 =*t1;
773  valueT2 =*t2;
774  valueE21 =*e21;
775  valueE22 =*e22;
776  valueE12 =*e12;
777  valueE11 =*e11;
778 
779  xs11 = pDiffCrossSectionData[LevelIndex][valueT1][valueE11];
780  xs12 = pDiffCrossSectionData[LevelIndex][valueT1][valueE12];
781  xs21 = pDiffCrossSectionData[LevelIndex][valueT2][valueE21];
782  xs22 = pDiffCrossSectionData[LevelIndex][valueT2][valueE22];
783  }
784  }
785 
786  // G4double xsProduct = xs11 * xs12 * xs21 * xs22;
787  // if (xsProduct != 0.)
788  // {
789  sigma = QuadInterpolator( valueE11, valueE12,
790  valueE21, valueE22,
791  xs11, xs12,
792  xs21, xs22,
793  valueT1, valueT2,
794  k, energyTransfer);
795  // }
796 
797  }
798 
799  return sigma;
800 }
801 
802 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
803 
804 G4double G4MicroElecInelasticModel::Interpolate(G4double e1,
805  G4double e2,
806  G4double e,
807  G4double xs1,
808  G4double xs2)
809 {
810  G4double value = 0.;
811 
812  // Log-log interpolation by default
813  if (e1 != 0 && e2 != 0 && (std::log10(e2) - std::log10(e1)) != 0
814  && !fasterCode)
815  {
816  G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
817  G4double b = std::log10(xs2) - a*std::log10(e2);
818  G4double sigma = a*std::log10(e) + b;
819  value = (std::pow(10.,sigma));
820 
821  }
822 
823  // Switch to log-lin interpolation for faster code
824  if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 && fasterCode)
825  {
826  G4double d1 = std::log10(xs1);
827  G4double d2 = std::log10(xs2);
828  value = std::pow(10., (d1 + (d2 - d1) * (e - e1) / (e2 - e1)));
829  }
830 
831  // Switch to lin-lin interpolation for faster code
832  // in case one of xs1 or xs2 (=cum proba) value is zero
833 
834  if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0)) // && fasterCode)
835  {
836  G4double d1 = xs1;
837  G4double d2 = xs2;
838  value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
839  }
840 
841 
842  return value;
843 }
844 
845 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
846 
847 G4double G4MicroElecInelasticModel::QuadInterpolator(G4double e11, G4double e12,
848  G4double e21, G4double e22,
849  G4double xs11, G4double xs12,
850  G4double xs21, G4double xs22,
851  G4double t1, G4double t2,
852  G4double t, G4double e)
853 {
854  G4double interpolatedvalue1 = Interpolate(e11, e12, e, xs11, xs12);
855  G4double interpolatedvalue2 = Interpolate(e21, e22, e, xs21, xs22);
856  G4double value = Interpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
857  return value;
858 }
859 
860 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
861 
862 G4int G4MicroElecInelasticModel::RandomSelect(G4double k, const G4String& particle )
863 {
864  G4int level = 0;
865 
866  std::map< G4String,G4MicroElecCrossSectionDataSet*,std::less<G4String> >::iterator pos;
867  pos = tableData.find(particle);
868 
869  if (pos != tableData.end())
870  {
871  G4MicroElecCrossSectionDataSet* table = pos->second;
872 
873  if (table != 0)
874  {
875  G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
876  const size_t n(table->NumberOfComponents());
877  size_t i(n);
878  G4double value = 0.;
879 
880  while (i>0)
881  {
882  i--;
883  valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
884  value += valuesBuffer[i];
885  }
886 
887  value *= G4UniformRand();
888 
889  i = n;
890 
891  while (i > 0)
892  {
893  i--;
894 
895  if (valuesBuffer[i] > value)
896  {
897  delete[] valuesBuffer;
898  return i;
899  }
900  value -= valuesBuffer[i];
901  }
902 
903  if (valuesBuffer) delete[] valuesBuffer;
904 
905  }
906  }
907  else
908  {
909  G4Exception("G4MicroElecInelasticModel::RandomSelect","em0002",FatalException,"Model not applicable to particle type.");
910  }
911 
912  return level;
913 }
914 
915 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
916 
917 G4double G4MicroElecInelasticModel::RandomizeEjectedElectronEnergyFromCumulatedDcs(G4ParticleDefinition* particleDefinition,
918  G4double k,
919  G4int shell)
920 {
921 
922  G4double secondaryElectronKineticEnergy = 0.;
923 
924  G4double random = G4UniformRand();
925 
926  secondaryElectronKineticEnergy = TransferedEnergy(particleDefinition,
927  k / eV,
928  shell,
929  random) * eV
930  - SiStructure.Energy(shell);
931 
932  if (secondaryElectronKineticEnergy < 0.)
933  return 0.;
934  //
935 
936  return secondaryElectronKineticEnergy;
937 }
938 
939 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
940 
941 G4double G4MicroElecInelasticModel::TransferedEnergy(G4ParticleDefinition* particleDefinition,
942  G4double k,
943  G4int ionizationLevelIndex,
944  G4double random)
945 {
946  G4double nrj = 0.;
947 
948  G4double valueK1 = 0;
949  G4double valueK2 = 0;
950  G4double valuePROB21 = 0;
951  G4double valuePROB22 = 0;
952  G4double valuePROB12 = 0;
953  G4double valuePROB11 = 0;
954 
955  G4double nrjTransf11 = 0;
956  G4double nrjTransf12 = 0;
957  G4double nrjTransf21 = 0;
958  G4double nrjTransf22 = 0;
959 
960  G4double maximumEnergyTransfer1 = 0;
961  G4double maximumEnergyTransfer2 = 0;
962  G4double maximumEnergyTransferP = 4.* (electron_mass_c2 / proton_mass_c2) * k;
963  G4double bindingEnergy = SiStructure.Energy(ionizationLevelIndex)*1e6;
964 
965  if (particleDefinition == G4Electron::ElectronDefinition())
966  {
967  // k should be in eV
968  std::vector<double>::iterator k2 = std::upper_bound(eTdummyVec.begin(),
969  eTdummyVec.end(),
970  k);
971  std::vector<double>::iterator k1 = k2 - 1;
972 
973  /*
974  G4cout << "----> k=" << k
975  << " " << *k1
976  << " " << *k2
977  << " " << random
978  << " " << ionizationLevelIndex
979  << " " << eProbaShellMap[ionizationLevelIndex][(*k1)].back()
980  << " " << eProbaShellMap[ionizationLevelIndex][(*k2)].back()
981  << G4endl;
982  */
983 
984  // SI : the following condition avoids situations where random >last vector element
985  if (random <= eProbaShellMap[ionizationLevelIndex][(*k1)].back()
986  && random <= eProbaShellMap[ionizationLevelIndex][(*k2)].back())
987  {
988  std::vector<double>::iterator prob12 =
989  std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
990  eProbaShellMap[ionizationLevelIndex][(*k1)].end(),
991  random);
992 
993  std::vector<double>::iterator prob11 = prob12 - 1;
994 
995  std::vector<double>::iterator prob22 =
996  std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
997  eProbaShellMap[ionizationLevelIndex][(*k2)].end(),
998  random);
999 
1000  std::vector<double>::iterator prob21 = prob22 - 1;
1001 
1002  valueK1 = *k1;
1003  valueK2 = *k2;
1004  valuePROB21 = *prob21;
1005  valuePROB22 = *prob22;
1006  valuePROB12 = *prob12;
1007  valuePROB11 = *prob11;
1008 
1009  /*
1010  G4cout << " " << random << " " << valuePROB11 << " "
1011  << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1012  */
1013 
1014  // The following condition avoid getting transfered energy < binding energy and forces cumxs = 1 for maximum energy transfer.
1015  if(valuePROB11 == 0) nrjTransf11 = bindingEnergy;
1016  else nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1017  if(valuePROB12 == 1)
1018  {
1019  if ((valueK1+bindingEnergy)/2. > valueK1) maximumEnergyTransfer1=valueK1;
1020  else maximumEnergyTransfer1 = (valueK1+bindingEnergy)/2.;
1021 
1022  nrjTransf12 = maximumEnergyTransfer1;
1023  }
1024  else nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1025 
1026  if(valuePROB21 == 0) nrjTransf21 = bindingEnergy;
1027  else nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1028  if(valuePROB22 == 1)
1029  {
1030  if ((valueK2+bindingEnergy)/2. > valueK2) maximumEnergyTransfer2=valueK2;
1031  else maximumEnergyTransfer2 = (valueK2+bindingEnergy)/2.;
1032 
1033  nrjTransf22 = maximumEnergyTransfer2;
1034  }
1035  else nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1036 
1037 
1038  /*nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1039  nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1040  nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1041  nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];*/
1042 
1043  /*
1044  G4cout << " " << ionizationLevelIndex << " "
1045  << random << " " <<valueK1 << " " << valueK2 << G4endl;
1046 
1047  G4cout << " " << random << " " << nrjTransf11 << " "
1048  << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1049  */
1050 
1051  }
1052  // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
1053  if (random > eProbaShellMap[ionizationLevelIndex][(*k1)].back())
1054  {
1055  std::vector<double>::iterator prob22 =
1056  std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1057  eProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1058  random);
1059 
1060  std::vector<double>::iterator prob21 = prob22 - 1;
1061 
1062  valueK1 = *k1;
1063  valueK2 = *k2;
1064  valuePROB21 = *prob21;
1065  valuePROB22 = *prob22;
1066 
1067  //G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1068 
1069  nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1070  nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1071 
1072  G4double interpolatedvalue2 = Interpolate(valuePROB21,
1073  valuePROB22,
1074  random,
1075  nrjTransf21,
1076  nrjTransf22);
1077 
1078  // zeros are explicitely set
1079 
1080  G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
1081 
1082  /*
1083  G4cout << " " << ionizationLevelIndex << " "
1084  << random << " " <<valueK1 << " " << valueK2 << G4endl;
1085 
1086  G4cout << " " << random << " " << nrjTransf11 << " "
1087  << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1088 
1089  G4cout << "ici" << " " << value << G4endl;
1090  */
1091 
1092  return value;
1093  }
1094  }
1095  //
1096  else if (particleDefinition == G4Proton::ProtonDefinition())
1097  {
1098  // k should be in eV
1099 
1100  std::vector<double>::iterator k2 = std::upper_bound(pTdummyVec.begin(),
1101  pTdummyVec.end(),
1102  k);
1103 
1104  std::vector<double>::iterator k1 = k2 - 1;
1105 
1106  /*
1107  G4cout << "----> k=" << k
1108  << " " << *k1
1109  << " " << *k2
1110  << " " << random
1111  << " " << ionizationLevelIndex
1112  << " " << pProbaShellMap[ionizationLevelIndex][(*k1)].back()
1113  << " " << pProbaShellMap[ionizationLevelIndex][(*k2)].back()
1114  << G4endl;
1115  */
1116 
1117  // SI : the following condition avoids situations where random > last vector element,
1118  // for eg. when the last element is zero
1119  if (random <= pProbaShellMap[ionizationLevelIndex][(*k1)].back()
1120  && random <= pProbaShellMap[ionizationLevelIndex][(*k2)].back())
1121  {
1122  std::vector<double>::iterator prob12 =
1123  std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
1124  pProbaShellMap[ionizationLevelIndex][(*k1)].end(),
1125  random);
1126 
1127  std::vector<double>::iterator prob11 = prob12 - 1;
1128 
1129  std::vector<double>::iterator prob22 =
1130  std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1131  pProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1132  random);
1133 
1134  std::vector<double>::iterator prob21 = prob22 - 1;
1135 
1136  valueK1 = *k1;
1137  valueK2 = *k2;
1138  valuePROB21 = *prob21;
1139  valuePROB22 = *prob22;
1140  valuePROB12 = *prob12;
1141  valuePROB11 = *prob11;
1142 
1143  /*
1144  G4cout << " " << random << " " << valuePROB11 << " "
1145  << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1146  */
1147 
1148  // The following condition avoid getting transfered energy < binding energy and forces cumxs = 1 for maximum energy transfer.
1149  if(valuePROB11 == 0) nrjTransf11 = bindingEnergy;
1150  else nrjTransf11 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1151  if(valuePROB12 == 1) nrjTransf12 = maximumEnergyTransferP;
1152  else nrjTransf12 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1153  if(valuePROB21 == 0) nrjTransf21 = bindingEnergy;
1154  else nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1155  if(valuePROB22 == 1) nrjTransf22 = maximumEnergyTransferP;
1156  else nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1157 
1158 
1159  /* nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1160  nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1161  nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1162  nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];*/
1163 
1164  /*
1165  G4cout << " " << ionizationLevelIndex << " "
1166  << random << " " <<valueK1 << " " << valueK2 << G4endl;
1167 
1168  G4cout << " " << random << " " << nrjTransf11 << " "
1169  << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1170  */
1171  }
1172 
1173  // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
1174 
1175  if (random > pProbaShellMap[ionizationLevelIndex][(*k1)].back())
1176  {
1177  std::vector<double>::iterator prob22 =
1178  std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1179  pProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1180  random);
1181 
1182  std::vector<double>::iterator prob21 = prob22 - 1;
1183 
1184  valueK1 = *k1;
1185  valueK2 = *k2;
1186  valuePROB21 = *prob21;
1187  valuePROB22 = *prob22;
1188 
1189  //G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1190 
1191  nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1192  nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1193 
1194  G4double interpolatedvalue2 = Interpolate(valuePROB21,
1195  valuePROB22,
1196  random,
1197  nrjTransf21,
1198  nrjTransf22);
1199 
1200  // zeros are explicitely set
1201 
1202  G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
1203 
1204  /*
1205  G4cout << " " << ionizationLevelIndex << " "
1206  << random << " " <<valueK1 << " " << valueK2 << G4endl;
1207 
1208  G4cout << " " << random << " " << nrjTransf11 << " "
1209  << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1210 
1211  G4cout << "ici" << " " << value << G4endl;
1212  */
1213 
1214  return value;
1215  }
1216  }
1217  // End electron and proton cases
1218 
1219  G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21
1220  * nrjTransf22;
1221  //G4cout << "nrjTransfProduct=" << nrjTransfProduct << G4endl;
1222 
1223  if (nrjTransfProduct != 0.)
1224  {
1225  nrj = QuadInterpolator(valuePROB11,
1226  valuePROB12,
1227  valuePROB21,
1228  valuePROB22,
1229  nrjTransf11,
1230  nrjTransf12,
1231  nrjTransf21,
1232  nrjTransf22,
1233  valueK1,
1234  valueK2,
1235  k,
1236  random);
1237  }
1238  //G4cout << nrj << endl;
1239 
1240  return nrj;
1241 }
1242 
1243 
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
G4Electron::ElectronDefinition
static G4Electron * ElectronDefinition()
Definition: G4Electron.cc:89
G4VEmModel::LowEnergyLimit
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:640
G4MicroElecInelasticModel::CrossSectionPerVolume
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
Definition: G4MicroElecInelasticModel.cc:324
G4NistManager::FindOrBuildMaterial
G4Material * FindOrBuildMaterial(const G4String &name, G4bool isotopes=true, G4bool warning=false)
Definition: G4NistManager.hh:474
G4MicroElecCrossSectionDataSet::GetComponent
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
Definition: G4MicroElecCrossSectionDataSet.hh:51
G4VEMDataSet::FindValue
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
G4LossTableManager::Instance
static G4LossTableManager * Instance()
Definition: G4LossTableManager.cc:112
CLHEP::proton_mass_c2
static constexpr double proton_mass_c2
Definition: PhysicalConstants.h:80
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
cm2
static constexpr double cm2
Definition: G4SIunits.hh:120
G4ionEffectiveCharge.hh
G4VEmModel::HighEnergyLimit
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:633
CLHEP::Hep3Vector::x
double x() const
G4MicroElecInelasticModel::G4MicroElecInelasticModel
G4MicroElecInelasticModel(const G4ParticleDefinition *p=0, const G4String &nam="MicroElecInelasticModel")
Definition: G4MicroElecInelasticModel.cc:62
d2
static const G4double d2
Definition: G4ElectroNuclearCrossSection.cc:89
G4MicroElecInelasticModel::Initialise
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4MicroElecInelasticModel.cc:115
G4Proton::ProtonDefinition
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
G4UnitsTable.hh
G4VEmModel::GetAngularDistribution
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:616
G4LossTableManager.hh
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4MicroElecInelasticModel::DifferentialCrossSection
double DifferentialCrossSection(G4ParticleDefinition *aParticleDefinition, G4double k, G4double energyTransfer, G4int shell)
Definition: G4MicroElecInelasticModel.cc:708
G4int
int G4int
Definition: G4Types.hh:78
G4NistManager::Instance
static G4NistManager * Instance()
Definition: G4NistManager.cc:68
G4ParticleChangeForGamma::ProposeMomentumDirection
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
Definition: G4ParticleChangeForGamma.hh:147
fKShell
Definition: G4AtomicShellEnumerator.hh:53
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:120
CLHEP::Hep3Vector::z
double z() const
G4VParticleChange::ProposeLocalEnergyDeposit
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4MicroElecSiStructure::Energy
G4double Energy(G4int level)
Definition: G4MicroElecSiStructure.cc:63
CLHEP::electron_mass_c2
static constexpr double electron_mass_c2
Definition: PhysicalConstants.h:79
G4VEmModel::SetHighEnergyLimit
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:724
G4VEmAngularDistribution::SampleDirectionForShell
virtual G4ThreeVector & SampleDirectionForShell(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, G4int shellID, const G4Material *)
Definition: G4VEmAngularDistribution.cc:60
G4VAtomDeexcitation::GetAtomicShell
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4cout
G4GLOB_DLL std::ostream G4cout
value
const XML_Char int const XML_Char * value
Definition: expat.h:331
G4MicroElecCrossSectionDataSet::NumberOfComponents
virtual size_t NumberOfComponents(void) const
Definition: G4MicroElecCrossSectionDataSet.hh:57
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
cm
static constexpr double cm
Definition: G4SIunits.hh:119
eV
static constexpr double eV
Definition: G4SIunits.hh:215
G4MicroElecInelasticModel::fParticleChangeForGamma
G4ParticleChangeForGamma * fParticleChangeForGamma
Definition: G4MicroElecInelasticModel.hh:97
G4MicroElecCrossSectionDataSet::FindValue
virtual G4double FindValue(G4double e, G4int componentId=0) const
Definition: G4MicroElecCrossSectionDataSet.cc:474
G4MicroElecCrossSectionDataSet::LoadData
virtual G4bool LoadData(const G4String &argFileName)
Definition: G4MicroElecCrossSectionDataSet.cc:58
d1
static const G4double d1
Definition: G4ElectroNuclearCrossSection.cc:88
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4PhysicalConstants.hh
G4UAtomicDeexcitation.hh
G4MicroElecInelasticModel.hh
G4ionEffectiveCharge::EffectiveCharge
G4double EffectiveCharge(const G4ParticleDefinition *p, const G4Material *material, G4double kineticEnergy)
Definition: G4ionEffectiveCharge.cc:89
G4Material::GetTotNbOfAtomsPerVolume
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:209
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:122
G4DeltaAngle.hh
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
G4VEmModel::SetAngularDistribution
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:623
G4ios.hh
CLHEP::Hep3Vector::y
double y() const
FatalException
Definition: G4ExceptionSeverity.hh:60
G4MicroElecInelasticModel::SampleSecondaries
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
Definition: G4MicroElecInelasticModel.cc:435
G4Material::GetBaseMaterial
const G4Material * GetBaseMaterial() const
Definition: G4Material.hh:233
GeV
static constexpr double GeV
Definition: G4SIunits.hh:217
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4ParticleChangeForGamma::SetProposedKineticEnergy
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Definition: G4ParticleChangeForGamma.hh:129
G4endl
#define G4endl
Definition: G4ios.hh:61
MeV
static constexpr double MeV
Definition: G4SIunits.hh:214
G4LossTableManager::AtomDeexcitation
G4VAtomDeexcitation * AtomDeexcitation()
Definition: G4LossTableManager.hh:444
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
G4VEmModel::SetLowEnergyLimit
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:731
G4VAtomDeexcitation::GenerateParticles
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
Definition: G4VAtomDeexcitation.hh:274
G4VEmModel::SetDeexcitationFlag
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:780
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:124
G4InuclSpecialFunctions::bindingEnergy
G4double bindingEnergy(G4int A, G4int Z)
Definition: bindingEnergy.cc:37
keV
static constexpr double keV
Definition: G4SIunits.hh:216
pos
static const G4double pos
Definition: G4ElectroNuclearCrossSection.cc:66
G4AtomicShellEnumerator
G4AtomicShellEnumerator
Definition: G4AtomicShellEnumerator.hh:50
G4MicroElecInelasticModel::TransferedEnergy
G4double TransferedEnergy(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell, G4double random)
Definition: G4MicroElecInelasticModel.cc:941
G4VEmModel::GetParticleChangeForGamma
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:132
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:150