ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4MicroElecInelasticModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4MicroElecInelasticModel.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // 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 = std::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