00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043 #include "G4MuElecInelasticModel.hh"
00044
00045 #include "globals.hh"
00046 #include "G4PhysicalConstants.hh"
00047 #include "G4SystemOfUnits.hh"
00048 #include "G4ios.hh"
00049 #include "G4UnitsTable.hh"
00050 #include "G4UAtomicDeexcitation.hh"
00051 #include "G4LossTableManager.hh"
00052 #include "G4ionEffectiveCharge.hh"
00053
00054
00055
00056 using namespace std;
00057
00058
00059
00060 G4MuElecInelasticModel::G4MuElecInelasticModel(const G4ParticleDefinition*,
00061 const G4String& nam)
00062 :G4VEmModel(nam),fAtomDeexcitation(0),isInitialised(false)
00063 {
00064 nistSi = G4NistManager::Instance()->FindOrBuildMaterial("G4_Si");
00065
00066 verboseLevel= 0;
00067
00068
00069
00070
00071
00072
00073
00074 if( verboseLevel>0 )
00075 {
00076 G4cout << "MuElec inelastic model is constructed " << G4endl;
00077 }
00078
00079
00080 SetDeexcitationFlag(true);
00081 fParticleChangeForGamma = 0;
00082 }
00083
00084
00085
00086 G4MuElecInelasticModel::~G4MuElecInelasticModel()
00087 {
00088
00089
00090 std::map< G4String,G4MuElecCrossSectionDataSet*,std::less<G4String> >::iterator pos;
00091 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
00092 {
00093 G4MuElecCrossSectionDataSet* table = pos->second;
00094 delete table;
00095 }
00096
00097
00098
00099 eVecm.clear();
00100 pVecm.clear();
00101
00102 }
00103
00104
00105
00106 void G4MuElecInelasticModel::Initialise(const G4ParticleDefinition* particle,
00107 const G4DataVector& )
00108 {
00109
00110 if (verboseLevel > 3)
00111 G4cout << "Calling G4MuElecInelasticModel::Initialise()" << G4endl;
00112
00113
00114
00115 G4String fileElectron("muelec/sigma_inelastic_e_Si");
00116 G4String fileProton("muelec/sigma_inelastic_p_Si");
00117
00118 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
00119 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
00120
00121 G4String electron;
00122 G4String proton;
00123
00124 G4double scaleFactor = 1e-18 * cm *cm;
00125
00126 char *path = getenv("G4LEDATA");
00127
00128
00129 electron = electronDef->GetParticleName();
00130
00131 tableFile[electron] = fileElectron;
00132
00133 lowEnergyLimit[electron] = 16.7 * eV;
00134 highEnergyLimit[electron] = 100.0 * MeV;
00135
00136
00137
00138 G4MuElecCrossSectionDataSet* tableE = new G4MuElecCrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
00139 tableE->LoadData(fileElectron);
00140
00141 tableData[electron] = tableE;
00142
00143
00144
00145 std::ostringstream eFullFileName;
00146 eFullFileName << path << "/muelec/sigmadiff_inelastic_e_Si.dat";
00147 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
00148
00149 if (!eDiffCrossSection)
00150 {
00151 G4Exception("G4MuElecInelasticModel::Initialise","em0003",FatalException,"Missing data file:/muelec/sigmadiff_inelastic_e_Si.dat");
00152 }
00153
00154 eTdummyVec.push_back(0.);
00155 while(!eDiffCrossSection.eof())
00156 {
00157 double tDummy;
00158 double eDummy;
00159 eDiffCrossSection>>tDummy>>eDummy;
00160 if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
00161 for (int j=0; j<6; j++)
00162 {
00163 eDiffCrossSection>>eDiffCrossSectionData[j][tDummy][eDummy];
00164
00165
00166 if (!eDiffCrossSection.eof()) eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
00167
00168 eVecm[tDummy].push_back(eDummy);
00169
00170 }
00171 }
00172
00173
00174
00175
00176 proton = protonDef->GetParticleName();
00177
00178 tableFile[proton] = fileProton;
00179
00180 lowEnergyLimit[proton] = 50. * keV;
00181 highEnergyLimit[proton] = 1. * GeV;
00182
00183
00184
00185 G4MuElecCrossSectionDataSet* tableP = new G4MuElecCrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
00186 tableP->LoadData(fileProton);
00187
00188 tableData[proton] = tableP;
00189
00190
00191
00192 std::ostringstream pFullFileName;
00193 pFullFileName << path << "/muelec/sigmadiff_inelastic_p_Si.dat";
00194 std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
00195
00196 if (!pDiffCrossSection)
00197 {
00198 G4Exception("G4MuElecInelasticModel::Initialise","em0003",FatalException,"Missing data file:/muelec/sigmadiff_inelastic_p_Si.dat");
00199 }
00200
00201 pTdummyVec.push_back(0.);
00202 while(!pDiffCrossSection.eof())
00203 {
00204 double tDummy;
00205 double eDummy;
00206 pDiffCrossSection>>tDummy>>eDummy;
00207 if (tDummy != pTdummyVec.back()) pTdummyVec.push_back(tDummy);
00208 for (int j=0; j<6; j++)
00209 {
00210 pDiffCrossSection>>pDiffCrossSectionData[j][tDummy][eDummy];
00211
00212
00213 if (!pDiffCrossSection.eof()) pDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
00214
00215 pVecm[tDummy].push_back(eDummy);
00216 }
00217 }
00218
00219
00220 if (particle==electronDef)
00221 {
00222 SetLowEnergyLimit(lowEnergyLimit[electron]);
00223 SetHighEnergyLimit(highEnergyLimit[electron]);
00224 }
00225
00226 if (particle==protonDef)
00227 {
00228 SetLowEnergyLimit(lowEnergyLimit[proton]);
00229 SetHighEnergyLimit(highEnergyLimit[proton]);
00230 }
00231
00232 if( verboseLevel>0 )
00233 {
00234 G4cout << "MuElec Inelastic model is initialized " << G4endl
00235 << "Energy range: "
00236 << LowEnergyLimit() / eV << " eV - "
00237 << HighEnergyLimit() / keV << " keV for "
00238 << particle->GetParticleName()
00239 << " with mass (amu) " << particle->GetPDGMass()/proton_mass_c2
00240 << " and charge " << particle->GetPDGCharge()
00241 << G4endl << G4endl ;
00242 }
00243
00244
00245
00246 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
00247
00248 if (isInitialised) { return; }
00249 fParticleChangeForGamma = GetParticleChangeForGamma();
00250 isInitialised = true;
00251
00252 }
00253
00254
00255
00256 G4double G4MuElecInelasticModel::CrossSectionPerVolume(const G4Material* material,
00257 const G4ParticleDefinition* particleDefinition,
00258 G4double ekin,
00259 G4double,
00260 G4double)
00261 {
00262 if (verboseLevel > 3)
00263 G4cout << "Calling CrossSectionPerVolume() of G4MuElecInelasticModel" << G4endl;
00264
00265 G4double density = material->GetTotNbOfAtomsPerVolume();
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279 G4double lowLim = 0;
00280 G4double highLim = 0;
00281 G4double sigma=0;
00282
00283 const G4String& particleName = particleDefinition->GetParticleName();
00284 G4String nameLocal = particleName ;
00285
00286 G4double Zeff2 = 1.0;
00287 G4double Mion_c2 = particleDefinition->GetPDGMass();
00288
00289 if (Mion_c2 > proton_mass_c2)
00290 {
00291 G4ionEffectiveCharge EffCharge ;
00292 G4double Zeff = EffCharge.EffectiveCharge(particleDefinition, material,ekin);
00293 Zeff2 = Zeff*Zeff;
00294
00295 if (verboseLevel > 3)
00296 G4cout << "Before scaling : " << G4endl
00297 << "Particle : " << nameLocal << ", mass : " << Mion_c2/proton_mass_c2 << "*mp, charge " << Zeff
00298 << ", Ekin (eV) = " << ekin/eV << G4endl ;
00299
00300 ekin *= proton_mass_c2/Mion_c2 ;
00301 nameLocal = "proton" ;
00302
00303 if (verboseLevel > 3)
00304 G4cout << "After scaling : " << G4endl
00305 << "Particle : " << nameLocal << ", Ekin (eV) = " << ekin/eV << G4endl ;
00306 }
00307
00308 if (material == nistSi || material->GetBaseMaterial() == nistSi)
00309 {
00310
00311 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
00312 pos1 = lowEnergyLimit.find(nameLocal);
00313 if (pos1 != lowEnergyLimit.end())
00314 {
00315 lowLim = pos1->second;
00316 }
00317
00318 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
00319 pos2 = highEnergyLimit.find(nameLocal);
00320 if (pos2 != highEnergyLimit.end())
00321 {
00322 highLim = pos2->second;
00323 }
00324
00325 if (ekin >= lowLim && ekin < highLim)
00326 {
00327 std::map< G4String,G4MuElecCrossSectionDataSet*,std::less<G4String> >::iterator pos;
00328 pos = tableData.find(nameLocal);
00329
00330 if (pos != tableData.end())
00331 {
00332 G4MuElecCrossSectionDataSet* table = pos->second;
00333 if (table != 0)
00334 {
00335 sigma = table->FindValue(ekin);
00336 }
00337 }
00338 else
00339 {
00340 G4Exception("G4MuElecInelasticModel::CrossSectionPerVolume","em0002",FatalException,"Model not applicable to particle type.");
00341 }
00342 }
00343 else
00344 {
00345 if (nameLocal!="e-")
00346 {
00347
00348
00349 }
00350 }
00351
00352 if (verboseLevel > 3)
00353 {
00354 G4cout << "---> Kinetic energy (eV)=" << ekin/eV << G4endl;
00355 G4cout << " - Cross section per Si atom (cm^2)=" << sigma*Zeff2/cm2 << G4endl;
00356 G4cout << " - Cross section per Si atom (cm^-1)=" << sigma*density*Zeff2/(1./cm) << G4endl;
00357 }
00358
00359 }
00360 return sigma*density*Zeff2;
00361
00362
00363 }
00364
00365
00366
00367 void G4MuElecInelasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
00368 const G4MaterialCutsCouple* ,
00369 const G4DynamicParticle* particle,
00370 G4double,
00371 G4double)
00372 {
00373
00374 if (verboseLevel > 3)
00375 G4cout << "Calling SampleSecondaries() of G4MuElecInelasticModel" << G4endl;
00376
00377 G4double lowLim = 0;
00378 G4double highLim = 0;
00379
00380 G4double ekin = particle->GetKineticEnergy();
00381 G4double k = ekin ;
00382
00383 G4ParticleDefinition* PartDef = particle->GetDefinition();
00384 const G4String& particleName = PartDef->GetParticleName();
00385 G4String nameLocal2 = particleName ;
00386 G4double particleMass = particle->GetDefinition()->GetPDGMass();
00387
00388 if (particleMass > proton_mass_c2)
00389 {
00390 k *= proton_mass_c2/particleMass ;
00391 PartDef = G4Proton::ProtonDefinition();
00392 nameLocal2 = "proton" ;
00393 }
00394
00395 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
00396 pos1 = lowEnergyLimit.find(nameLocal2);
00397
00398 if (pos1 != lowEnergyLimit.end())
00399 {
00400 lowLim = pos1->second;
00401 }
00402
00403 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
00404 pos2 = highEnergyLimit.find(nameLocal2);
00405
00406 if (pos2 != highEnergyLimit.end())
00407 {
00408 highLim = pos2->second;
00409 }
00410
00411 if (k >= lowLim && k < highLim)
00412 {
00413 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
00414 G4double totalEnergy = ekin + particleMass;
00415 G4double pSquare = ekin * (totalEnergy + particleMass);
00416 G4double totalMomentum = std::sqrt(pSquare);
00417
00418 G4int Shell = RandomSelect(k,nameLocal2);
00419 G4double bindingEnergy = SiStructure.Energy(Shell);
00420 if (verboseLevel > 3)
00421 {
00422 G4cout << "---> Kinetic energy (eV)=" << k/eV << G4endl ;
00423 G4cout << "Shell: " << Shell << ", energy: " << bindingEnergy/eV << G4endl;
00424 }
00425
00426
00427
00428 G4int secNumberInit = 0;
00429 G4int secNumberFinal = 0;
00430
00431 if(fAtomDeexcitation && Shell > 3) {
00432 G4int Z = 14;
00433 G4AtomicShellEnumerator as = fKShell;
00434
00435 if (Shell == 5)
00436 {
00437 as = G4AtomicShellEnumerator(1);
00438 }
00439 else if (Shell == 4)
00440 {
00441 as = G4AtomicShellEnumerator(3);
00442 }
00443
00444 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
00445 secNumberInit = fvect->size();
00446 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
00447 secNumberFinal = fvect->size();
00448 }
00449
00450 G4double secondaryKinetic = RandomizeEjectedElectronEnergy(PartDef,k,Shell);
00451
00452 if (verboseLevel > 3)
00453 {
00454 G4cout << "Ionisation process" << G4endl;
00455 G4cout << "Shell: " << Shell << " Kin. energy (eV)=" << k/eV
00456 << " Sec. energy (eV)=" << secondaryKinetic/eV << G4endl;
00457 }
00458
00459 G4double cosTheta = 0.;
00460 G4double phi = 0.;
00461 RandomizeEjectedElectronDirection(PartDef, k, secondaryKinetic, cosTheta, phi);
00462
00463 G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
00464 G4double dirX = sinTheta*std::cos(phi);
00465 G4double dirY = sinTheta*std::sin(phi);
00466 G4double dirZ = cosTheta;
00467 G4ThreeVector deltaDirection(dirX,dirY,dirZ);
00468 deltaDirection.rotateUz(primaryDirection);
00469
00470
00471
00472 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
00473
00474 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
00475 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
00476 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
00477 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
00478 finalPx /= finalMomentum;
00479 finalPy /= finalMomentum;
00480 finalPz /= finalMomentum;
00481
00482 G4ThreeVector direction;
00483 direction.set(finalPx,finalPy,finalPz);
00484
00485 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
00486
00487
00488
00489
00490 G4double deexSecEnergy = 0;
00491 for (G4int j=secNumberInit; j < secNumberFinal; j++) {
00492 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();}
00493
00494 fParticleChangeForGamma->SetProposedKineticEnergy(ekin-bindingEnergy-secondaryKinetic);
00495 fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy-deexSecEnergy);
00496
00497 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
00498 fvect->push_back(dp);
00499
00500 }
00501
00502 }
00503
00504
00505
00506 G4double G4MuElecInelasticModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
00507 G4double k, G4int shell)
00508 {
00509 if (particleDefinition == G4Electron::ElectronDefinition())
00510 {
00511 G4double maximumEnergyTransfer=0.;
00512 if ((k+SiStructure.Energy(shell))/2. > k) maximumEnergyTransfer=k;
00513 else maximumEnergyTransfer = (k+SiStructure.Energy(shell))/2.;
00514
00515 G4double crossSectionMaximum = 0.;
00516
00517 G4double minEnergy = SiStructure.Energy(shell);
00518 G4double maxEnergy = maximumEnergyTransfer;
00519 G4int nEnergySteps = 100;
00520
00521 G4double value(minEnergy);
00522 G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
00523 G4int step(nEnergySteps);
00524 while (step>0)
00525 {
00526 step--;
00527 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
00528 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
00529 value*=stpEnergy;
00530 }
00531
00532
00533 G4double secondaryElectronKineticEnergy=0.;
00534 do
00535 {
00536 secondaryElectronKineticEnergy = G4UniformRand() * (maximumEnergyTransfer-SiStructure.Energy(shell));
00537 } while(G4UniformRand()*crossSectionMaximum >
00538 DifferentialCrossSection(particleDefinition, k/eV,(secondaryElectronKineticEnergy+SiStructure.Energy(shell))/eV,shell));
00539
00540 return secondaryElectronKineticEnergy;
00541
00542 }
00543
00544 if (particleDefinition == G4Proton::ProtonDefinition())
00545 {
00546 G4double maximumEnergyTransfer = 4.* (electron_mass_c2 / proton_mass_c2) * k;
00547 G4double crossSectionMaximum = 0.;
00548
00549 G4double minEnergy = SiStructure.Energy(shell);
00550 G4double maxEnergy = maximumEnergyTransfer;
00551 G4int nEnergySteps = 100;
00552
00553 G4double value(minEnergy);
00554 G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
00555 G4int step(nEnergySteps);
00556 while (step>0)
00557 {
00558 step--;
00559 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
00560 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
00561 value*=stpEnergy;
00562 }
00563 G4double secondaryElectronKineticEnergy = 0.;
00564 do
00565 {
00566 secondaryElectronKineticEnergy = G4UniformRand() * (maximumEnergyTransfer-SiStructure.Energy(shell));
00567
00568 } while(G4UniformRand()*crossSectionMaximum >=
00569 DifferentialCrossSection(particleDefinition, k/eV,(secondaryElectronKineticEnergy+SiStructure.Energy(shell))/eV,shell));
00570 return secondaryElectronKineticEnergy;
00571 }
00572
00573 return 0;
00574 }
00575
00576
00577
00578 void G4MuElecInelasticModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
00579 G4double k,
00580 G4double secKinetic,
00581 G4double & cosTheta,
00582 G4double & phi )
00583 {
00584 if (particleDefinition == G4Electron::ElectronDefinition())
00585 {
00586 phi = twopi * G4UniformRand();
00587 G4double sin2O = (1.-secKinetic/k) / (1.+secKinetic/(2.*electron_mass_c2));
00588 cosTheta = std::sqrt(1.-sin2O);
00589 }
00590
00591 if (particleDefinition == G4Proton::ProtonDefinition())
00592 {
00593 G4double maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
00594 phi = twopi * G4UniformRand();
00595 cosTheta = std::sqrt(secKinetic / maxSecKinetic);
00596 }
00597
00598 else
00599 {
00600 G4double maxSecKinetic = 4.* (electron_mass_c2 / particleDefinition->GetPDGMass()) * k;
00601 phi = twopi * G4UniformRand();
00602 cosTheta = std::sqrt(secKinetic / maxSecKinetic);
00603 }
00604 }
00605
00606
00607
00608 double G4MuElecInelasticModel::DifferentialCrossSection(G4ParticleDefinition * particleDefinition,
00609 G4double k,
00610 G4double energyTransfer,
00611 G4int LevelIndex)
00612 {
00613 G4double sigma = 0.;
00614
00615 if (energyTransfer >= SiStructure.Energy(LevelIndex))
00616 {
00617 G4double valueT1 = 0;
00618 G4double valueT2 = 0;
00619 G4double valueE21 = 0;
00620 G4double valueE22 = 0;
00621 G4double valueE12 = 0;
00622 G4double valueE11 = 0;
00623
00624 G4double xs11 = 0;
00625 G4double xs12 = 0;
00626 G4double xs21 = 0;
00627 G4double xs22 = 0;
00628
00629 if (particleDefinition == G4Electron::ElectronDefinition())
00630 {
00631
00632
00633 std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
00634 std::vector<double>::iterator t1 = t2-1;
00635
00636 if (energyTransfer <= eVecm[(*t1)].back() && energyTransfer <= eVecm[(*t2)].back() )
00637 {
00638 std::vector<double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), energyTransfer);
00639 std::vector<double>::iterator e11 = e12-1;
00640
00641 std::vector<double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), energyTransfer);
00642 std::vector<double>::iterator e21 = e22-1;
00643
00644 valueT1 =*t1;
00645 valueT2 =*t2;
00646 valueE21 =*e21;
00647 valueE22 =*e22;
00648 valueE12 =*e12;
00649 valueE11 =*e11;
00650
00651 xs11 = eDiffCrossSectionData[LevelIndex][valueT1][valueE11];
00652 xs12 = eDiffCrossSectionData[LevelIndex][valueT1][valueE12];
00653 xs21 = eDiffCrossSectionData[LevelIndex][valueT2][valueE21];
00654 xs22 = eDiffCrossSectionData[LevelIndex][valueT2][valueE22];
00655 }
00656
00657 }
00658
00659 if (particleDefinition == G4Proton::ProtonDefinition())
00660 {
00661
00662 std::vector<double>::iterator t2 = std::upper_bound(pTdummyVec.begin(),pTdummyVec.end(), k);
00663 std::vector<double>::iterator t1 = t2-1;
00664 if (energyTransfer <= pVecm[(*t1)].back() && energyTransfer <= pVecm[(*t2)].back() )
00665 {
00666 std::vector<double>::iterator e12 = std::upper_bound(pVecm[(*t1)].begin(),pVecm[(*t1)].end(), energyTransfer);
00667 std::vector<double>::iterator e11 = e12-1;
00668
00669 std::vector<double>::iterator e22 = std::upper_bound(pVecm[(*t2)].begin(),pVecm[(*t2)].end(), energyTransfer);
00670 std::vector<double>::iterator e21 = e22-1;
00671
00672 valueT1 =*t1;
00673 valueT2 =*t2;
00674 valueE21 =*e21;
00675 valueE22 =*e22;
00676 valueE12 =*e12;
00677 valueE11 =*e11;
00678
00679 xs11 = pDiffCrossSectionData[LevelIndex][valueT1][valueE11];
00680 xs12 = pDiffCrossSectionData[LevelIndex][valueT1][valueE12];
00681 xs21 = pDiffCrossSectionData[LevelIndex][valueT2][valueE21];
00682 xs22 = pDiffCrossSectionData[LevelIndex][valueT2][valueE22];
00683 }
00684 }
00685
00686 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
00687 if (xsProduct != 0.)
00688 {
00689 sigma = QuadInterpolator( valueE11, valueE12,
00690 valueE21, valueE22,
00691 xs11, xs12,
00692 xs21, xs22,
00693 valueT1, valueT2,
00694 k, energyTransfer);
00695 }
00696
00697 }
00698
00699 return sigma;
00700 }
00701
00702
00703
00704 G4double G4MuElecInelasticModel::LogLogInterpolate(G4double e1,
00705 G4double e2,
00706 G4double e,
00707 G4double xs1,
00708 G4double xs2)
00709 {
00710 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
00711 G4double b = std::log10(xs2) - a*std::log10(e2);
00712 G4double sigma = a*std::log10(e) + b;
00713 G4double value = (std::pow(10.,sigma));
00714 return value;
00715 }
00716
00717
00718
00719 G4double G4MuElecInelasticModel::QuadInterpolator(G4double e11, G4double e12,
00720 G4double e21, G4double e22,
00721 G4double xs11, G4double xs12,
00722 G4double xs21, G4double xs22,
00723 G4double t1, G4double t2,
00724 G4double t, G4double e)
00725 {
00726 G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
00727 G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
00728 G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
00729 return value;
00730 }
00731
00732
00733
00734 G4int G4MuElecInelasticModel::RandomSelect(G4double k, const G4String& particle )
00735 {
00736 G4int level = 0;
00737
00738 std::map< G4String,G4MuElecCrossSectionDataSet*,std::less<G4String> >::iterator pos;
00739 pos = tableData.find(particle);
00740
00741 if (pos != tableData.end())
00742 {
00743 G4MuElecCrossSectionDataSet* table = pos->second;
00744
00745 if (table != 0)
00746 {
00747 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
00748 const size_t n(table->NumberOfComponents());
00749 size_t i(n);
00750 G4double value = 0.;
00751
00752 while (i>0)
00753 {
00754 i--;
00755 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
00756 value += valuesBuffer[i];
00757 }
00758
00759 value *= G4UniformRand();
00760
00761 i = n;
00762
00763 while (i > 0)
00764 {
00765 i--;
00766
00767 if (valuesBuffer[i] > value)
00768 {
00769 delete[] valuesBuffer;
00770 return i;
00771 }
00772 value -= valuesBuffer[i];
00773 }
00774
00775 if (valuesBuffer) delete[] valuesBuffer;
00776
00777 }
00778 }
00779 else
00780 {
00781 G4Exception("G4MuElecInelasticModel::RandomSelect","em0002",FatalException,"Model not applicable to particle type.");
00782 }
00783
00784 return level;
00785 }
00786
00787
00788
00789