Geant4  10.02.p01
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 
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 
63  const G4String& nam)
64 :G4VEmModel(nam),fAtomDeexcitation(0),isInitialised(false)
65 {
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);
84 
85  // default generator
87 }
88 
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
90 
92 {
93  // Cross section
94 
95  std::map< G4String,G4MicroElecCrossSectionDataSet*,std::less<G4String> >::iterator pos;
96  for (pos = tableData.begin(); pos != tableData.end(); ++pos)
97  {
98  G4MicroElecCrossSectionDataSet* table = pos->second;
99  delete table;
100  }
101 
102  // Final state
103 
104  eVecm.clear();
105  pVecm.clear();
106 
107 }
108 
109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
110 
112  const G4DataVector& /*cuts*/)
113 {
114 
115  if (verboseLevel > 3)
116  G4cout << "Calling G4MicroElecInelasticModel::Initialise()" << G4endl;
117 
118  // Energy limits
119 
120  G4String fileElectron("microelec/sigma_inelastic_e_Si");
121  G4String fileProton("microelec/sigma_inelastic_p_Si");
122 
125 
128 
129  G4double scaleFactor = 1e-18 * cm *cm;
130 
131  char *path = getenv("G4LEDATA");
132 
133  // *** ELECTRON
134  electron = electronDef->GetParticleName();
135 
136  tableFile[electron] = fileElectron;
137 
138  lowEnergyLimit[electron] = 16.7 * eV;
139  highEnergyLimit[electron] = 100.0 * MeV;
140 
141  // Cross section
142 
144  tableE->LoadData(fileElectron);
145 
146  tableData[electron] = tableE;
147 
148  // Final state
149 
150  std::ostringstream eFullFileName;
151  eFullFileName << path << "/microelec/sigmadiff_inelastic_e_Si.dat";
152  std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
153 
154  if (!eDiffCrossSection)
155  {
156  G4Exception("G4MicroElecInelasticModel::Initialise","em0003",FatalException,"Missing data file:/microelec/sigmadiff_inelastic_e_Si.dat");
157  }
158 
159  //
160 
161  // Clear the arrays for re-initialization case (MT mode)
162  // Octobre 22nd, 2014 - Melanie Raine
163 
164  eTdummyVec.clear();
165  pTdummyVec.clear();
166 
167  eVecm.clear();
168  pVecm.clear();
169 
170  eDiffCrossSectionData->clear();
171  pDiffCrossSectionData->clear();
172 
173  //
174 
175 
176  eTdummyVec.push_back(0.);
177  while(!eDiffCrossSection.eof())
178  {
179  double tDummy;
180  double eDummy;
181  eDiffCrossSection>>tDummy>>eDummy;
182  if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
183  for (int j=0; j<6; j++)
184  {
185  eDiffCrossSection>>eDiffCrossSectionData[j][tDummy][eDummy];
186 
187  // SI - only if eof is not reached !
188  if (!eDiffCrossSection.eof()) eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
189 
190  eVecm[tDummy].push_back(eDummy);
191 
192  }
193  }
194  //
195 
196  // *** PROTON
197 
198  proton = protonDef->GetParticleName();
199 
200  tableFile[proton] = fileProton;
201 
202  lowEnergyLimit[proton] = 50. * keV;
203  highEnergyLimit[proton] = 10. * GeV;
204 
205  // Cross section
206 
208  tableP->LoadData(fileProton);
209 
210  tableData[proton] = tableP;
211 
212  // Final state
213 
214  std::ostringstream pFullFileName;
215  pFullFileName << path << "/microelec/sigmadiff_inelastic_p_Si.dat";
216  std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
217 
218  if (!pDiffCrossSection)
219  {
220  G4Exception("G4MicroElecInelasticModel::Initialise","em0003",FatalException,"Missing data file:/microelec/sigmadiff_inelastic_p_Si.dat");
221  }
222 
223  pTdummyVec.push_back(0.);
224  while(!pDiffCrossSection.eof())
225  {
226  double tDummy;
227  double eDummy;
228  pDiffCrossSection>>tDummy>>eDummy;
229  if (tDummy != pTdummyVec.back()) pTdummyVec.push_back(tDummy);
230  for (int j=0; j<6; j++)
231  {
232  pDiffCrossSection>>pDiffCrossSectionData[j][tDummy][eDummy];
233 
234  // SI - only if eof is not reached !
235  if (!pDiffCrossSection.eof()) pDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
236 
237  pVecm[tDummy].push_back(eDummy);
238  }
239  }
240 
241 
242  if (particle==electronDef)
243  {
246  }
247 
248  if (particle==protonDef)
249  {
252  }
253 
254  if( verboseLevel>0 )
255  {
256  G4cout << "MicroElec Inelastic model is initialized " << G4endl
257  << "Energy range: "
258  << LowEnergyLimit() / keV << " keV - "
259  << HighEnergyLimit() / MeV << " MeV for "
260  << particle->GetParticleName()
261  << " with mass (amu) " << particle->GetPDGMass()/proton_mass_c2
262  << " and charge " << particle->GetPDGCharge()
263  << G4endl << G4endl ;
264  }
265 
266  //
267 
269 
270  if (isInitialised) { return; }
272  isInitialised = true;
273 
274 }
275 
276 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
277 
279  const G4ParticleDefinition* particleDefinition,
280  G4double ekin,
281  G4double,
282  G4double)
283 {
284  if (verboseLevel > 3)
285  G4cout << "Calling CrossSectionPerVolume() of G4MicroElecInelasticModel" << G4endl;
286 
288 
289  /* if (
290  particleDefinition != G4Proton::ProtonDefinition()
291  &&
292  particleDefinition != G4Electron::ElectronDefinition()
293  &&
294  particleDefinition != G4GenericIon::GenericIonDefinition()
295  )
296 
297  return 0;*/
298 
299  // Calculate total cross section for model
300 
301  G4double lowLim = 0;
302  G4double highLim = 0;
303  G4double sigma=0;
304 
305  const G4String& particleName = particleDefinition->GetParticleName();
306  G4String nameLocal = particleName ;
307 
308  G4double Zeff2 = 1.0;
309  G4double Mion_c2 = particleDefinition->GetPDGMass();
310 
311  if (Mion_c2 > proton_mass_c2)
312  {
313  G4ionEffectiveCharge EffCharge ;
314  G4double Zeff = EffCharge.EffectiveCharge(particleDefinition, material,ekin);
315  Zeff2 = Zeff*Zeff;
316 
317  if (verboseLevel > 3)
318  G4cout << "Before scaling : " << G4endl
319  << "Particle : " << nameLocal << ", mass : " << Mion_c2/proton_mass_c2 << "*mp, charge " << Zeff
320  << ", Ekin (eV) = " << ekin/eV << G4endl ;
321 
322  ekin *= proton_mass_c2/Mion_c2 ;
323  nameLocal = "proton" ;
324 
325  if (verboseLevel > 3)
326  G4cout << "After scaling : " << G4endl
327  << "Particle : " << nameLocal << ", Ekin (eV) = " << ekin/eV << G4endl ;
328  }
329 
330  if (material == nistSi || material->GetBaseMaterial() == nistSi)
331  {
332 
333  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
334  pos1 = lowEnergyLimit.find(nameLocal);
335  if (pos1 != lowEnergyLimit.end())
336  {
337  lowLim = pos1->second;
338  }
339 
340  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
341  pos2 = highEnergyLimit.find(nameLocal);
342  if (pos2 != highEnergyLimit.end())
343  {
344  highLim = pos2->second;
345  }
346 
347  if (ekin >= lowLim && ekin < highLim)
348  {
349  std::map< G4String,G4MicroElecCrossSectionDataSet*,std::less<G4String> >::iterator pos;
350  pos = tableData.find(nameLocal);
351 
352  if (pos != tableData.end())
353  {
354  G4MicroElecCrossSectionDataSet* table = pos->second;
355  if (table != 0)
356  {
357  sigma = table->FindValue(ekin);
358  }
359  }
360  else
361  {
362  G4Exception("G4MicroElecInelasticModel::CrossSectionPerVolume","em0002",FatalException,"Model not applicable to particle type.");
363  }
364  }
365  else
366  {
367  if (nameLocal!="e-")
368  {
369  // G4cout << "Particle : " << nameLocal << ", Ekin (eV) = " << ekin/eV << G4endl;
370  // G4cout << "### Warning: particle energy out of bounds! ###" << G4endl;
371  }
372  }
373 
374  if (verboseLevel > 3)
375  {
376  G4cout << "---> Kinetic energy (eV)=" << ekin/eV << G4endl;
377  G4cout << " - Cross section per Si atom (cm^2)=" << sigma*Zeff2/cm2 << G4endl;
378  G4cout << " - Cross section per Si atom (cm^-1)=" << sigma*density*Zeff2/(1./cm) << G4endl;
379  }
380 
381  } // if (SiMaterial)
382  return sigma*density*Zeff2;
383 
384 
385 }
386 
387 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
388 
389 void G4MicroElecInelasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
390  const G4MaterialCutsCouple* couple,
391  const G4DynamicParticle* particle,
392  G4double,
393  G4double)
394 {
395 
396  if (verboseLevel > 3)
397  G4cout << "Calling SampleSecondaries() of G4MicroElecInelasticModel" << G4endl;
398 
399  G4double lowLim = 0;
400  G4double highLim = 0;
401 
402  G4double ekin = particle->GetKineticEnergy();
403  G4double k = ekin ;
404 
405  G4ParticleDefinition* PartDef = particle->GetDefinition();
406  const G4String& particleName = PartDef->GetParticleName();
407  G4String nameLocal2 = particleName ;
408  G4double particleMass = particle->GetDefinition()->GetPDGMass();
409 
410  if (particleMass > proton_mass_c2)
411  {
412  k *= proton_mass_c2/particleMass ;
413  PartDef = G4Proton::ProtonDefinition();
414  nameLocal2 = "proton" ;
415  }
416 
417  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
418  pos1 = lowEnergyLimit.find(nameLocal2);
419 
420  if (pos1 != lowEnergyLimit.end())
421  {
422  lowLim = pos1->second;
423  }
424 
425  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
426  pos2 = highEnergyLimit.find(nameLocal2);
427 
428  if (pos2 != highEnergyLimit.end())
429  {
430  highLim = pos2->second;
431  }
432 
433  if (k >= lowLim && k < highLim)
434  {
435  G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
436  G4double totalEnergy = ekin + particleMass;
437  G4double pSquare = ekin * (totalEnergy + particleMass);
438  G4double totalMomentum = std::sqrt(pSquare);
439 
440  G4int Shell = RandomSelect(k,nameLocal2);
442 
443  if (verboseLevel > 3)
444  {
445  G4cout << "---> Kinetic energy (eV)=" << k/eV << G4endl ;
446  G4cout << "Shell: " << Shell << ", energy: " << bindingEnergy/eV << G4endl;
447  }
448 
449  // sample deexcitation
450 
451  G4int secNumberInit = 0; // need to know at a certain point the energy of secondaries
452  G4int secNumberFinal = 0; // So I'll make the difference and then sum the energies
453 
454  G4int Z = 14;
455 
456  if(fAtomDeexcitation && Shell > 2) {
457 
459 
460  if (Shell == 4)
461  {
462  as = G4AtomicShellEnumerator(1);
463  }
464  else if (Shell == 3)
465  {
466  as = G4AtomicShellEnumerator(3);
467  }
468 
469  const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
470  secNumberInit = fvect->size();
471  fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
472  secNumberFinal = fvect->size();
473  }
474 
475  G4double secondaryKinetic = RandomizeEjectedElectronEnergy(PartDef,k,Shell);
476 
477  if (verboseLevel > 3)
478  {
479  G4cout << "Ionisation process" << G4endl;
480  G4cout << "Shell: " << Shell << " Kin. energy (eV)=" << k/eV
481  << " Sec. energy (eV)=" << secondaryKinetic/eV << G4endl;
482  }
483 
484  G4ThreeVector deltaDirection =
485  GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
486  Z, Shell,
487  couple->GetMaterial());
488 
489  //if (particle->GetDefinition() == G4Electron::ElectronDefinition())
490  //{
491  G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
492 
493  G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
494  G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
495  G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
496  G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
497  finalPx /= finalMomentum;
498  finalPy /= finalMomentum;
499  finalPz /= finalMomentum;
500 
501  G4ThreeVector direction;
502  direction.set(finalPx,finalPy,finalPz);
503 
505  //}
506  //else fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection) ;
507 
508  // note that secondaryKinetic is the energy of the delta ray, not of all secondaries.
509  G4double deexSecEnergy = 0;
510  for (G4int j=secNumberInit; j < secNumberFinal; j++) {
511  deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();}
512 
513  fParticleChangeForGamma->SetProposedKineticEnergy(ekin-bindingEnergy-secondaryKinetic);
514  fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy-deexSecEnergy);
515 
516  G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
517  fvect->push_back(dp);
518 
519  }
520 
521 }
522 
523 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
524 
526  G4double k, G4int shell)
527 {
528  if (particleDefinition == G4Electron::ElectronDefinition())
529  {
530  G4double maximumEnergyTransfer=0.;
531  if ((k+SiStructure.Energy(shell))/2. > k) maximumEnergyTransfer=k;
532  else maximumEnergyTransfer = (k+SiStructure.Energy(shell))/2.;
533 
534  G4double crossSectionMaximum = 0.;
535 
536  G4double minEnergy = SiStructure.Energy(shell);
537  G4double maxEnergy = maximumEnergyTransfer;
538  G4int nEnergySteps = 100;
539 
540  G4double value(minEnergy);
541  G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
542  G4int step(nEnergySteps);
543  while (step>0)
544  {
545  step--;
546  G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
547  if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
548  value*=stpEnergy;
549  }
550 
551 
552  G4double secondaryElectronKineticEnergy=0.;
553  do
554  {
555  secondaryElectronKineticEnergy = G4UniformRand() * (maximumEnergyTransfer-SiStructure.Energy(shell));
556  } while(G4UniformRand()*crossSectionMaximum >
557  DifferentialCrossSection(particleDefinition, k/eV,(secondaryElectronKineticEnergy+SiStructure.Energy(shell))/eV,shell));
558 
559  return secondaryElectronKineticEnergy;
560 
561  }
562 
563  if (particleDefinition == G4Proton::ProtonDefinition())
564  {
565  G4double maximumEnergyTransfer = 4.* (electron_mass_c2 / proton_mass_c2) * k;
566  G4double crossSectionMaximum = 0.;
567 
568  G4double minEnergy = SiStructure.Energy(shell);
569  G4double maxEnergy = maximumEnergyTransfer;
570  G4int nEnergySteps = 100;
571 
572  G4double value(minEnergy);
573  G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
574  G4int step(nEnergySteps);
575  while (step>0)
576  {
577  step--;
578  G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
579  if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
580  value*=stpEnergy;
581  }
582  G4double secondaryElectronKineticEnergy = 0.;
583  do
584  {
585  secondaryElectronKineticEnergy = G4UniformRand() * (maximumEnergyTransfer-SiStructure.Energy(shell));
586 
587  } while(G4UniformRand()*crossSectionMaximum >
588  DifferentialCrossSection(particleDefinition, k/eV,(secondaryElectronKineticEnergy+SiStructure.Energy(shell))/eV,shell));
589  return secondaryElectronKineticEnergy;
590  }
591 
592  return 0;
593 }
594 
595 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
596 
597 // The following section is not used anymore but is kept for memory
598 // GetAngularDistribution()->SampleDirectionForShell is used instead
599 
600 /*void G4MicroElecInelasticModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
601  G4double k,
602  G4double secKinetic,
603  G4double & cosTheta,
604  G4double & phi )
605  {
606  if (particleDefinition == G4Electron::ElectronDefinition())
607  {
608  phi = twopi * G4UniformRand();
609  G4double sin2O = (1.-secKinetic/k) / (1.+secKinetic/(2.*electron_mass_c2));
610  cosTheta = std::sqrt(1.-sin2O);
611  }
612 
613  if (particleDefinition == G4Proton::ProtonDefinition())
614  {
615  G4double maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
616  phi = twopi * G4UniformRand();
617  cosTheta = std::sqrt(secKinetic / maxSecKinetic);
618  }
619 
620  else
621  {
622  G4double maxSecKinetic = 4.* (electron_mass_c2 / particleDefinition->GetPDGMass()) * k;
623  phi = twopi * G4UniformRand();
624  cosTheta = std::sqrt(secKinetic / maxSecKinetic);
625  }
626  }
627  */
628 
629 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
630 
632  G4double k,
633  G4double energyTransfer,
634  G4int LevelIndex)
635 {
636  G4double sigma = 0.;
637 
638  if (energyTransfer >= SiStructure.Energy(LevelIndex))
639  {
640  G4double valueT1 = 0;
641  G4double valueT2 = 0;
642  G4double valueE21 = 0;
643  G4double valueE22 = 0;
644  G4double valueE12 = 0;
645  G4double valueE11 = 0;
646 
647  G4double xs11 = 0;
648  G4double xs12 = 0;
649  G4double xs21 = 0;
650  G4double xs22 = 0;
651 
652  if (particleDefinition == G4Electron::ElectronDefinition())
653  {
654  // k should be in eV and energy transfer eV also
655 
656  std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
657  std::vector<double>::iterator t1 = t2-1;
658  // SI : the following condition avoids situations where energyTransfer >last vector element
659  if (energyTransfer <= eVecm[(*t1)].back() && energyTransfer <= eVecm[(*t2)].back() )
660  {
661  std::vector<double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), energyTransfer);
662  std::vector<double>::iterator e11 = e12-1;
663 
664  std::vector<double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), energyTransfer);
665  std::vector<double>::iterator e21 = e22-1;
666 
667  valueT1 =*t1;
668  valueT2 =*t2;
669  valueE21 =*e21;
670  valueE22 =*e22;
671  valueE12 =*e12;
672  valueE11 =*e11;
673 
674  xs11 = eDiffCrossSectionData[LevelIndex][valueT1][valueE11];
675  xs12 = eDiffCrossSectionData[LevelIndex][valueT1][valueE12];
676  xs21 = eDiffCrossSectionData[LevelIndex][valueT2][valueE21];
677  xs22 = eDiffCrossSectionData[LevelIndex][valueT2][valueE22];
678  }
679 
680  }
681 
682  if (particleDefinition == G4Proton::ProtonDefinition())
683  {
684  // k should be in eV and energy transfer eV also
685  std::vector<double>::iterator t2 = std::upper_bound(pTdummyVec.begin(),pTdummyVec.end(), k);
686  std::vector<double>::iterator t1 = t2-1;
687  if (energyTransfer <= pVecm[(*t1)].back() && energyTransfer <= pVecm[(*t2)].back() )
688  {
689  std::vector<double>::iterator e12 = std::upper_bound(pVecm[(*t1)].begin(),pVecm[(*t1)].end(), energyTransfer);
690  std::vector<double>::iterator e11 = e12-1;
691 
692  std::vector<double>::iterator e22 = std::upper_bound(pVecm[(*t2)].begin(),pVecm[(*t2)].end(), energyTransfer);
693  std::vector<double>::iterator e21 = e22-1;
694 
695  valueT1 =*t1;
696  valueT2 =*t2;
697  valueE21 =*e21;
698  valueE22 =*e22;
699  valueE12 =*e12;
700  valueE11 =*e11;
701 
702  xs11 = pDiffCrossSectionData[LevelIndex][valueT1][valueE11];
703  xs12 = pDiffCrossSectionData[LevelIndex][valueT1][valueE12];
704  xs21 = pDiffCrossSectionData[LevelIndex][valueT2][valueE21];
705  xs22 = pDiffCrossSectionData[LevelIndex][valueT2][valueE22];
706  }
707  }
708 
709  G4double xsProduct = xs11 * xs12 * xs21 * xs22;
710  if (xsProduct != 0.)
711  {
712  sigma = QuadInterpolator( valueE11, valueE12,
713  valueE21, valueE22,
714  xs11, xs12,
715  xs21, xs22,
716  valueT1, valueT2,
717  k, energyTransfer);
718  }
719 
720  }
721 
722  return sigma;
723 }
724 
725 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
726 
728  G4double e2,
729  G4double e,
730  G4double xs1,
731  G4double xs2)
732 {
733  G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
734  G4double b = std::log10(xs2) - a*std::log10(e2);
735  G4double sigma = a*std::log10(e) + b;
736  G4double value = (std::pow(10.,sigma));
737  return value;
738 }
739 
740 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
741 
743  G4double e21, G4double e22,
744  G4double xs11, G4double xs12,
745  G4double xs21, G4double xs22,
746  G4double t1, G4double t2,
747  G4double t, G4double e)
748 {
749  G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
750  G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
751  G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
752  return value;
753 }
754 
755 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
756 
758 {
759  G4int level = 0;
760 
761  std::map< G4String,G4MicroElecCrossSectionDataSet*,std::less<G4String> >::iterator pos;
762  pos = tableData.find(particle);
763 
764  if (pos != tableData.end())
765  {
766  G4MicroElecCrossSectionDataSet* table = pos->second;
767 
768  if (table != 0)
769  {
770  G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
771  const size_t n(table->NumberOfComponents());
772  size_t i(n);
773  G4double value = 0.;
774 
775  while (i>0)
776  {
777  i--;
778  valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
779  value += valuesBuffer[i];
780  }
781 
782  value *= G4UniformRand();
783 
784  i = n;
785 
786  while (i > 0)
787  {
788  i--;
789 
790  if (valuesBuffer[i] > value)
791  {
792  delete[] valuesBuffer;
793  return i;
794  }
795  value -= valuesBuffer[i];
796  }
797 
798  if (valuesBuffer) delete[] valuesBuffer;
799 
800  }
801  }
802  else
803  {
804  G4Exception("G4MicroElecInelasticModel::RandomSelect","em0002",FatalException,"Model not applicable to particle type.");
805  }
806 
807  return level;
808 }
809 
810 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
811 
812 
813 
static const double cm
Definition: G4SIunits.hh:118
static G4Electron * ElectronDefinition()
Definition: G4Electron.cc:89
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:641
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
G4Material * FindOrBuildMaterial(const G4String &name, G4bool isotopes=true, G4bool warning=false)
static const double MeV
Definition: G4SIunits.hh:211
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
static G4LossTableManager * Instance()
static const double cm2
Definition: G4SIunits.hh:119
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:634
G4MicroElecInelasticModel(const G4ParticleDefinition *p=0, const G4String &nam="MicroElecInelasticModel")
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:617
static const G4double e2
G4ParticleDefinition * GetDefinition() const
G4double a
Definition: TRTMaterials.hh:39
double DifferentialCrossSection(G4ParticleDefinition *aParticleDefinition, G4double k, G4double energyTransfer, G4int shell)
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
const G4String & GetParticleName() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4double Energy(G4int level)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:725
virtual G4ThreeVector & SampleDirectionForShell(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, G4int shellID, const G4Material *)
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
G4double density
Definition: TRTMaterials.hh:39
G4double RandomizeEjectedElectronEnergy(G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell)
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
virtual size_t NumberOfComponents(void) const
const G4ThreeVector & GetMomentumDirection() const
G4int RandomSelect(G4double energy, const G4String &particle)
static const double GeV
Definition: G4SIunits.hh:214
G4double QuadInterpolator(G4double e11, G4double e12, G4double e21, G4double e22, G4double x11, G4double x12, G4double x21, G4double x22, G4double t1, G4double t2, G4double t, G4double e)
const G4int n
G4ParticleChangeForGamma * fParticleChangeForGamma
virtual G4double FindValue(G4double e, G4int componentId=0) const
virtual G4bool LoadData(const G4String &argFileName)
static const G4double e1
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double EffectiveCharge(const G4ParticleDefinition *p, const G4Material *material, G4double kineticEnergy)
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:209
std::map< G4String, G4double, std::less< G4String > > highEnergyLimit
static const double eV
Definition: G4SIunits.hh:212
G4double GetPDGMass() const
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:624
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
const G4Material * GetBaseMaterial() const
Definition: G4Material.hh:233
static G4Electron * Electron()
Definition: G4Electron.cc:94
void SetProposedKineticEnergy(G4double proposedKinEnergy)
#define G4endl
Definition: G4ios.hh:61
G4VAtomDeexcitation * AtomDeexcitation()
static const double keV
Definition: G4SIunits.hh:213
double G4double
Definition: G4Types.hh:76
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:732
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:781
G4double GetPDGCharge() const
G4ThreeVector G4ParticleMomentum
G4double bindingEnergy(G4int A, G4int Z)
G4VAtomDeexcitation * fAtomDeexcitation
std::map< G4String, G4double, std::less< G4String > > lowEnergyLimit
static const G4double pos
G4AtomicShellEnumerator
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:134
const G4Material * GetMaterial() const
G4double LogLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)