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