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