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