source: trunk/source/processes/electromagnetic/lowenergy/src/G4DNABornIonisationModel.cc @ 968

Last change on this file since 968 was 968, checked in by garnier, 15 years ago

fichier ajoutes
File size: 21.9 KB
Line 
1// ********************************************************************
2// * License and Disclaimer                                           *
3// *                                                                  *
4// * The  Geant4 software  is  copyright of the Copyright Holders  of *
5// * the Geant4 Collaboration.  It is provided  under  the terms  and *
6// * conditions of the Geant4 Software License,  included in the file *
7// * LICENSE and available at  http://cern.ch/geant4/license .  These *
8// * include a list of copyright holders.                             *
9// *                                                                  *
10// * Neither the authors of this software system, nor their employing *
11// * institutes,nor the agencies providing financial support for this *
12// * work  make  any representation or  warranty, express or implied, *
13// * regarding  this  software system or assume any liability for its *
14// * use.  Please see the license in the file  LICENSE  and URL above *
15// * for the full disclaimer and the limitation of liability.         *
16// *                                                                  *
17// * This  code  implementation is the result of  the  scientific and *
18// * technical work of the GEANT4 collaboration.                      *
19// * By using,  copying,  modifying or  distributing the software (or *
20// * any work based  on the software)  you  agree  to acknowledge its *
21// * use  in  resulting  scientific  publications,  and indicate your *
22// * acceptance of all terms of the Geant4 Software license.          *
23// ********************************************************************
24//
25// $Id: G4DNABornIonisationModel.cc,v 1.4 2009/02/16 11:00:11 sincerti Exp $
26// GEANT4 tag $Name: geant4-09-02-ref-02 $
27//
28
29#include "G4DNABornIonisationModel.hh"
30
31//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
32
33using namespace std;
34
35//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
36
37G4DNABornIonisationModel::G4DNABornIonisationModel(const G4ParticleDefinition*,
38                                             const G4String& nam)
39:G4VEmModel(nam),isInitialised(false)
40{
41  verboseLevel= 0;
42  // Verbosity scale:
43  // 0 = nothing
44  // 1 = warning for energy non-conservation
45  // 2 = details of energy budget
46  // 3 = calculation of cross sections, file openings, sampling of atoms
47  // 4 = entering in methods
48 
49  G4cout << "Born ionisation model is constructed " << G4endl;
50}
51
52//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
53
54G4DNABornIonisationModel::~G4DNABornIonisationModel()
55{ 
56  // Cross section
57 
58  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
59  for (pos = tableData.begin(); pos != tableData.end(); ++pos)
60  {
61    G4DNACrossSectionDataSet* table = pos->second;
62    delete table;
63  }
64 
65  // Final state
66 
67  eVecm.clear();
68  pVecm.clear();
69
70}
71
72//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73
74void G4DNABornIonisationModel::Initialise(const G4ParticleDefinition* particle,
75                                       const G4DataVector& /*cuts*/)
76{
77
78  if (verboseLevel > 3)
79    G4cout << "Calling G4DNABornIonisationModel::Initialise()" << G4endl;
80
81  // Energy limits
82
83  G4String fileElectron("dna/sigma_ionisation_e_born");
84  G4String fileProton("dna/sigma_ionisation_p_born");
85
86  G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
87  G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
88
89  G4String electron;
90  G4String proton;
91 
92  G4double scaleFactor = (1.e-22 / 3.343) * m*m;
93
94  char *path = getenv("G4LEDATA");
95
96  if (electronDef != 0)
97  {
98    electron = electronDef->GetParticleName();
99
100    tableFile[electron] = fileElectron;
101
102    lowEnergyLimit[electron] = 12.61 * eV; 
103    highEnergyLimit[electron] = 30. * keV;
104
105    // Cross section
106   
107    G4DNACrossSectionDataSet* tableE = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
108    tableE->LoadData(fileElectron);
109     
110    tableData[electron] = tableE;
111   
112    // Final state
113   
114    std::ostringstream eFullFileName;
115    eFullFileName << path << "/dna/sigmadiff_ionisation_e_born.dat";
116    std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
117
118    if (!eDiffCrossSection)
119    { 
120      G4Exception("G4DNABornIonisationModel::ERROR OPENING electron DATA FILE");
121    }
122     
123    eTdummyVec.push_back(0.);
124    while(!eDiffCrossSection.eof())
125    {
126      double tDummy;
127      double eDummy;
128      eDiffCrossSection>>tDummy>>eDummy;
129      if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
130      for (int j=0; j<5; j++)
131      {
132        eDiffCrossSection>>eDiffCrossSectionData[j][tDummy][eDummy];
133
134        // SI - only if eof is not reached !
135        if (!eDiffCrossSection.eof()) eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
136
137        eVecm[tDummy].push_back(eDummy);
138
139      }
140    }
141
142    //
143  }
144  else
145  {
146    G4Exception("G4DNABornIonisationModel::Initialise(): electron is not defined");
147  }
148
149  if (protonDef != 0)
150  {
151    proton = protonDef->GetParticleName();
152
153    tableFile[proton] = fileProton;
154
155    lowEnergyLimit[proton] = 500. * keV;
156    highEnergyLimit[proton] = 10. * MeV;
157
158    // Cross section
159   
160    G4DNACrossSectionDataSet* tableP = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
161    tableP->LoadData(fileProton);
162     
163    tableData[proton] = tableP;
164   
165    // Final state
166
167    std::ostringstream pFullFileName;
168    pFullFileName << path << "/dna/sigmadiff_ionisation_p_born.dat";
169    std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
170   
171    if (!pDiffCrossSection)
172    { 
173      G4Exception("G4DNABornIonisationModel::ERROR OPENING proton DATA FILE");
174    }
175     
176    pTdummyVec.push_back(0.);
177    while(!pDiffCrossSection.eof())
178    {
179      double tDummy;
180      double eDummy;
181      pDiffCrossSection>>tDummy>>eDummy;
182      if (tDummy != pTdummyVec.back()) pTdummyVec.push_back(tDummy);
183      for (int j=0; j<5; j++)
184      {
185        pDiffCrossSection>>pDiffCrossSectionData[j][tDummy][eDummy];
186
187        // SI - only if eof is not reached !
188        if (!pDiffCrossSection.eof()) pDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
189
190        pVecm[tDummy].push_back(eDummy);
191      }
192    }
193 
194  }
195  else
196  {
197    G4Exception("G4DNABornIonisationModel::Initialise(): proton is not defined");
198  }
199
200  if (particle==electronDef) 
201  {
202    SetLowEnergyLimit(lowEnergyLimit[electron]);
203    SetHighEnergyLimit(highEnergyLimit[electron]);
204  }
205
206  if (particle==protonDef) 
207  {
208    SetLowEnergyLimit(lowEnergyLimit[proton]);
209    SetHighEnergyLimit(highEnergyLimit[proton]);
210  }
211
212  G4cout << "Born ionisation model is initialized " << G4endl
213         << "Energy range: "
214         << LowEnergyLimit() / eV << " eV - "
215         << HighEnergyLimit() / keV << " keV for "
216         << particle->GetParticleName()
217         << G4endl;
218
219  //
220 
221  if(!isInitialised) 
222  {
223    isInitialised = true;
224 
225    if(pParticleChange)
226      fParticleChangeForGamma = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
227    else
228      fParticleChangeForGamma = new G4ParticleChangeForGamma();
229  }   
230
231  // InitialiseElementSelectors(particle,cuts);
232
233  // Test if water material
234
235  flagMaterialIsWater= false;
236  densityWater = 0;
237
238  const G4ProductionCutsTable* theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
239
240  if(theCoupleTable) 
241  {
242    G4int numOfCouples = theCoupleTable->GetTableSize();
243 
244    if(numOfCouples>0) 
245    {
246          for (G4int i=0; i<numOfCouples; i++) 
247          {
248            const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
249            const G4Material* material = couple->GetMaterial();
250
251            size_t j = material->GetNumberOfElements();
252            while (j>0)
253            {
254               j--;
255               const G4Element* element(material->GetElement(j));
256               if (element->GetZ() == 8.)
257               {
258                  G4double density = material->GetAtomicNumDensityVector()[j];
259                  if (density > 0.) 
260                  { 
261                    flagMaterialIsWater = true; 
262                    densityWater = density; 
263                   
264                    if (verboseLevel > 3) 
265                    G4cout << "Water material is found with density(cm^-3)=" << density/(cm*cm*cm) << G4endl;
266                  }
267               }
268            }
269 
270          }
271   } // if(numOfCouples>0)
272
273  } // if (theCoupleTable)
274
275}
276
277//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
278
279G4double G4DNABornIonisationModel::CrossSectionPerVolume(const G4Material*,
280                                           const G4ParticleDefinition* particleDefinition,
281                                           G4double ekin,
282                                           G4double,
283                                           G4double)
284{
285  if (verboseLevel > 3)
286    G4cout << "Calling CrossSectionPerVolume() of G4DNABornIonisationModel" << G4endl;
287
288  if (
289      particleDefinition != G4Proton::ProtonDefinition()
290      &&
291      particleDefinition != G4Electron::ElectronDefinition()
292     )
293           
294    return 0;
295 
296  // Calculate total cross section for model
297
298  G4double lowLim = 0;
299  G4double highLim = 0;
300  G4double sigma=0;
301
302  if (flagMaterialIsWater)
303  {
304    const G4String& particleName = particleDefinition->GetParticleName();
305 
306    std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
307    pos1 = lowEnergyLimit.find(particleName);
308    if (pos1 != lowEnergyLimit.end())
309    {
310      lowLim = pos1->second;
311    }
312 
313    std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
314    pos2 = highEnergyLimit.find(particleName);
315    if (pos2 != highEnergyLimit.end())
316    {
317      highLim = pos2->second;
318    }
319
320    if (ekin >= lowLim && ekin < highLim)
321    {
322      std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
323      pos = tableData.find(particleName);
324       
325      if (pos != tableData.end())
326      {
327        G4DNACrossSectionDataSet* table = pos->second;
328        if (table != 0)
329        {
330          sigma = table->FindValue(ekin);
331        }
332      }
333      else
334      {
335        G4Exception("G4DNABornIonisationModel::CrossSectionPerVolume: attempting to calculate cross section for wrong particle");
336      }
337    }
338
339    if (verboseLevel > 3)
340    {
341      G4cout << "---> Kinetic energy(eV)=" << ekin/eV << G4endl;
342      G4cout << " - Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
343      G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*densityWater/(1./cm) << G4endl;
344    } 
345 
346  } // if (waterMaterial)
347 
348 return sigma*densityWater;               
349
350}
351
352//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
353
354void G4DNABornIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
355                                              const G4MaterialCutsCouple* /*couple*/,
356                                              const G4DynamicParticle* particle,
357                                              G4double,
358                                              G4double)
359{
360
361  if (verboseLevel > 3)
362    G4cout << "Calling SampleSecondaries() of G4DNABornIonisationModel" << G4endl;
363
364  G4double lowLim = 0;
365  G4double highLim = 0;
366
367  G4double k = particle->GetKineticEnergy();
368
369  const G4String& particleName = particle->GetDefinition()->GetParticleName();
370
371  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
372  pos1 = lowEnergyLimit.find(particleName);
373
374  if (pos1 != lowEnergyLimit.end())
375  {
376    lowLim = pos1->second;
377  }
378
379  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
380  pos2 = highEnergyLimit.find(particleName);
381
382  if (pos2 != highEnergyLimit.end())
383  {
384    highLim = pos2->second;
385  }
386
387  if (k >= lowLim && k < highLim)
388  {
389    G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
390    G4double particleMass = particle->GetDefinition()->GetPDGMass();
391    G4double totalEnergy = k + particleMass;
392    G4double pSquare = k * (totalEnergy + particleMass);
393    G4double totalMomentum = std::sqrt(pSquare);
394
395    G4int ionizationShell = RandomSelect(k,particleName);
396 
397    G4double secondaryKinetic = RandomizeEjectedElectronEnergy(particle->GetDefinition(),k,ionizationShell);
398 
399    G4double bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
400
401    G4double cosTheta = 0.;
402    G4double phi = 0.; 
403    RandomizeEjectedElectronDirection(particle->GetDefinition(), k,secondaryKinetic, cosTheta, phi);
404
405    G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
406    G4double dirX = sinTheta*std::cos(phi);
407    G4double dirY = sinTheta*std::sin(phi);
408    G4double dirZ = cosTheta;
409    G4ThreeVector deltaDirection(dirX,dirY,dirZ);
410    deltaDirection.rotateUz(primaryDirection);
411
412    G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
413
414    G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
415    G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
416    G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
417    G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
418    finalPx /= finalMomentum;
419    finalPy /= finalMomentum;
420    finalPz /= finalMomentum;
421   
422    G4ThreeVector direction;
423    direction.set(finalPx,finalPy,finalPz);
424   
425    fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
426    fParticleChangeForGamma->SetProposedKineticEnergy(k-bindingEnergy-secondaryKinetic);
427    fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy);
428
429    G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
430    fvect->push_back(dp);
431  }
432
433}
434
435//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
436
437G4double G4DNABornIonisationModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition, 
438G4double k, G4int shell)
439{
440  if (particleDefinition == G4Electron::ElectronDefinition()) 
441  {
442    G4double maximumEnergyTransfer=0.;
443    if ((k+waterStructure.IonisationEnergy(shell))/2. > k) maximumEnergyTransfer=k;
444    else maximumEnergyTransfer = (k+waterStructure.IonisationEnergy(shell))/2.;
445   
446    G4double crossSectionMaximum = 0.;
447    for(G4double value=waterStructure.IonisationEnergy(shell); value<=maximumEnergyTransfer; value+=0.1*eV)
448    {
449      G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
450      if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
451    }
452 
453    G4double secondaryElectronKineticEnergy=0.;
454    do 
455    {
456      secondaryElectronKineticEnergy = G4UniformRand() * (maximumEnergyTransfer-waterStructure.IonisationEnergy(shell));
457    } while(G4UniformRand()*crossSectionMaximum >
458      DifferentialCrossSection(particleDefinition, k/eV,(secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
459
460    return secondaryElectronKineticEnergy;
461 
462  }
463 
464  if (particleDefinition == G4Proton::ProtonDefinition()) 
465  {
466    G4double maximumKineticEnergyTransfer = 4.* (electron_mass_c2 / proton_mass_c2) * k - (waterStructure.IonisationEnergy(shell));
467
468    G4double crossSectionMaximum = 0.;
469    for (G4double value = waterStructure.IonisationEnergy(shell); 
470         value<=4.*waterStructure.IonisationEnergy(shell) ; 
471         value+=0.1*eV)
472    {
473      G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
474      if (differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
475    }
476
477    G4double secondaryElectronKineticEnergy = 0.;
478    do
479    {
480      secondaryElectronKineticEnergy = G4UniformRand() * maximumKineticEnergyTransfer;
481    } while(G4UniformRand()*crossSectionMaximum >= 
482              DifferentialCrossSection(particleDefinition, k/eV,(secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
483
484    return secondaryElectronKineticEnergy;
485  }
486
487  return 0;
488}
489
490//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
491
492void G4DNABornIonisationModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition, 
493                                                                   G4double k, 
494                                                                   G4double secKinetic, 
495                                                                   G4double & cosTheta, 
496                                                                   G4double & phi )
497{
498  if (particleDefinition == G4Electron::ElectronDefinition()) 
499  {
500    phi = twopi * G4UniformRand();
501    if (secKinetic < 50.*eV) cosTheta = (2.*G4UniformRand())-1.;
502    else if (secKinetic <= 200.*eV)     
503    {
504      if (G4UniformRand() <= 0.1) cosTheta = (2.*G4UniformRand())-1.;
505      else cosTheta = G4UniformRand()*(std::sqrt(2.)/2);
506    }
507    else       
508    {
509      G4double sin2O = (1.-secKinetic/k) / (1.+secKinetic/(2.*electron_mass_c2));
510      cosTheta = std::sqrt(1.-sin2O);
511    }
512  }
513 
514  if (particleDefinition == G4Proton::ProtonDefinition()) 
515  {
516    G4double maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
517    phi = twopi * G4UniformRand();
518    cosTheta = std::sqrt(secKinetic / maxSecKinetic);
519  }                     
520}
521
522//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
523
524double G4DNABornIonisationModel::DifferentialCrossSection(G4ParticleDefinition * particleDefinition, 
525                                                            G4double k, 
526                                                            G4double energyTransfer, 
527                                                            G4int ionizationLevelIndex)
528{
529  G4double sigma = 0.;
530
531  if (energyTransfer >= waterStructure.IonisationEnergy(ionizationLevelIndex))
532  {
533    G4double valueT1 = 0;
534    G4double valueT2 = 0;
535    G4double valueE21 = 0;
536    G4double valueE22 = 0;
537    G4double valueE12 = 0;
538    G4double valueE11 = 0;
539
540    G4double xs11 = 0;   
541    G4double xs12 = 0; 
542    G4double xs21 = 0; 
543    G4double xs22 = 0; 
544
545    if (particleDefinition == G4Electron::ElectronDefinition()) 
546    {
547      // k should be in eV and energy transfer eV also
548
549      std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
550
551      std::vector<double>::iterator t1 = t2-1;
552
553      // SI : the following condition avoids situations where energyTransfer >last vector element
554      if (energyTransfer <= eVecm[(*t1)].back())
555      {
556        std::vector<double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), energyTransfer);
557        std::vector<double>::iterator e11 = e12-1;
558
559        std::vector<double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), energyTransfer);
560        std::vector<double>::iterator e21 = e22-1;
561
562        valueT1  =*t1;
563        valueT2  =*t2;
564        valueE21 =*e21;
565        valueE22 =*e22;
566        valueE12 =*e12;
567        valueE11 =*e11;
568
569        xs11 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
570        xs12 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
571        xs21 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
572        xs22 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
573      }
574
575    }
576 
577   if (particleDefinition == G4Proton::ProtonDefinition()) 
578   {
579      // k should be in eV and energy transfer eV also
580      std::vector<double>::iterator t2 = std::upper_bound(pTdummyVec.begin(),pTdummyVec.end(), k);
581      std::vector<double>::iterator t1 = t2-1;
582     
583        std::vector<double>::iterator e12 = std::upper_bound(pVecm[(*t1)].begin(),pVecm[(*t1)].end(), energyTransfer);
584        std::vector<double>::iterator e11 = e12-1;
585
586        std::vector<double>::iterator e22 = std::upper_bound(pVecm[(*t2)].begin(),pVecm[(*t2)].end(), energyTransfer);
587        std::vector<double>::iterator e21 = e22-1;
588 
589        valueT1  =*t1;
590        valueT2  =*t2;
591        valueE21 =*e21;
592        valueE22 =*e22;
593        valueE12 =*e12;
594        valueE11 =*e11;
595
596        xs11 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
597        xs12 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
598        xs21 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
599        xs22 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
600
601   }
602 
603   G4double xsProduct = xs11 * xs12 * xs21 * xs22;
604   if (xsProduct != 0.)
605   {
606     sigma = QuadInterpolator(     valueE11, valueE12, 
607                                   valueE21, valueE22, 
608                                   xs11, xs12, 
609                                   xs21, xs22, 
610                                   valueT1, valueT2, 
611                                   k, energyTransfer);
612   }
613 
614 }
615 
616  return sigma;
617}
618
619//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
620
621G4double G4DNABornIonisationModel::LogLogInterpolate(G4double e1, 
622                                                       G4double e2, 
623                                                       G4double e, 
624                                                       G4double xs1, 
625                                                       G4double xs2)
626{
627  G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
628  G4double b = std::log10(xs2) - a*std::log10(e2);
629  G4double sigma = a*std::log10(e) + b;
630  G4double value = (std::pow(10.,sigma));
631  return value;
632}
633
634//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
635
636G4double G4DNABornIonisationModel::QuadInterpolator(G4double e11, G4double e12, 
637                                                      G4double e21, G4double e22, 
638                                                      G4double xs11, G4double xs12, 
639                                                      G4double xs21, G4double xs22, 
640                                                      G4double t1, G4double t2, 
641                                                      G4double t, G4double e)
642{
643  G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
644  G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
645  G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
646  return value;
647}
648
649//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
650
651G4int G4DNABornIonisationModel::RandomSelect(G4double k, const G4String& particle )
652{   
653  G4int level = 0;
654
655  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
656  pos = tableData.find(particle);
657
658  if (pos != tableData.end())
659  {
660    G4DNACrossSectionDataSet* table = pos->second;
661
662    if (table != 0)
663    {
664      G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
665      const size_t n(table->NumberOfComponents());
666      size_t i(n);
667      G4double value = 0.;
668           
669      while (i>0)
670      { 
671        i--;
672        valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
673        value += valuesBuffer[i];
674      }
675           
676      value *= G4UniformRand();
677   
678      i = n;
679           
680      while (i > 0)
681      {
682        i--;
683
684        if (valuesBuffer[i] > value)
685        {
686          delete[] valuesBuffer;
687          return i;
688        }
689        value -= valuesBuffer[i];
690      }
691           
692      if (valuesBuffer) delete[] valuesBuffer;
693   
694    }
695  }
696  else
697  {
698    G4Exception("G4DNABornIonisationModel::RandomSelect attempting to calculate cross section for wrong particle");
699  }
700     
701  return level;
702}
703
Note: See TracBrowser for help on using the repository browser.