ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNABornIonisationModel1.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DNABornIonisationModel1.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 
28 #include "G4DNABornIonisationModel1.hh"
29 #include "G4PhysicalConstants.hh"
30 #include "G4SystemOfUnits.hh"
31 #include "G4UAtomicDeexcitation.hh"
32 #include "G4LossTableManager.hh"
33 #include "G4DNAChemistryManager.hh"
34 #include "G4DNAMolecularMaterial.hh"
35 #include "G4DNABornAngle.hh"
36 #include "G4DeltaAngle.hh"
37 #include "G4Exp.hh"
38 
39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
40 
41 using namespace std;
42 
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44 
45 G4DNABornIonisationModel1::G4DNABornIonisationModel1(const G4ParticleDefinition*,
46  const G4String& nam) :
47 G4VEmModel(nam), isInitialised(false)
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  // Define default angular generator
69  SetAngularDistribution(new G4DNABornAngle());
70 
71  // Selection of computation method
72 
73  fasterCode = false;
74 
75  // Selection of stationary mode
76 
77  statCode = false;
78 
79  // Selection of SP scaling
80 
81  spScaling = true;
82 }
83 
84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
85 
86 G4DNABornIonisationModel1::~G4DNABornIonisationModel1()
87 {
88  // Cross section
89 
90  std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
91  for (pos = tableData.begin(); pos != tableData.end(); ++pos)
92  {
93  G4DNACrossSectionDataSet* table = pos->second;
94  delete table;
95  }
96 
97  // Final state
98 
99  eVecm.clear();
100  pVecm.clear();
101 }
102 
103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
104 
105 void G4DNABornIonisationModel1::Initialise(const G4ParticleDefinition* particle,
106  const G4DataVector& /*cuts*/)
107 {
108 
109  if (verboseLevel > 3)
110  {
111  G4cout << "Calling G4DNABornIonisationModel1::Initialise()" << G4endl;
112  }
113 
114  // Energy limits
115 
116  G4String fileElectron("dna/sigma_ionisation_e_born");
117  G4String fileProton("dna/sigma_ionisation_p_born");
118 
119  G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
120  G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
121 
122  G4String electron;
123  G4String proton;
124 
125  G4double scaleFactor = (1.e-22 / 3.343) * m*m;
126 
127  char *path = getenv("G4LEDATA");
128 
129  // *** ELECTRON
130 
131  electron = electronDef->GetParticleName();
132 
133  tableFile[electron] = fileElectron;
134 
135  lowEnergyLimit[electron] = 11. * eV;
136  highEnergyLimit[electron] = 1. * MeV;
137 
138  // Cross section
139 
140  G4DNACrossSectionDataSet* tableE = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
141  tableE->LoadData(fileElectron);
142 
143  tableData[electron] = tableE;
144 
145  // Final state
146 
147  std::ostringstream eFullFileName;
148 
149  if (fasterCode) eFullFileName << path << "/dna/sigmadiff_cumulated_ionisation_e_born_hp.dat";
150  if (!fasterCode) eFullFileName << path << "/dna/sigmadiff_ionisation_e_born.dat";
151 
152  std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
153 
154  if (!eDiffCrossSection)
155  {
156  if (fasterCode) G4Exception("G4DNABornIonisationModel1::Initialise","em0003",
157  FatalException,"Missing data file:/dna/sigmadiff_cumulated_ionisation_e_born_hp.dat");
158 
159  if (!fasterCode) G4Exception("G4DNABornIonisationModel1::Initialise","em0003",
160  FatalException,"Missing data file:/dna/sigmadiff_ionisation_e_born.dat");
161  }
162 
163  // Clear the arrays for re-initialization case (MT mode)
164  // March 25th, 2014 - Vaclav Stepan, Sebastien Incerti
165 
166  eTdummyVec.clear();
167  pTdummyVec.clear();
168 
169  eVecm.clear();
170  pVecm.clear();
171 
172  for (G4int j=0; j<5; j++)
173  {
174  eProbaShellMap[j].clear();
175  pProbaShellMap[j].clear();
176 
177  eDiffCrossSectionData[j].clear();
178  pDiffCrossSectionData[j].clear();
179 
180  eNrjTransfData[j].clear();
181  pNrjTransfData[j].clear();
182  }
183 
184  //
185 
186  eTdummyVec.push_back(0.);
187  while(!eDiffCrossSection.eof())
188  {
189  G4double tDummy;
190  G4double eDummy;
191  eDiffCrossSection>>tDummy>>eDummy;
192  if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
193 
194  G4double tmp;
195  for (G4int j=0; j<5; j++)
196  {
197  eDiffCrossSection>> tmp;
198 
199  eDiffCrossSectionData[j][tDummy][eDummy] = tmp;
200 
201  if (fasterCode)
202  {
203  eNrjTransfData[j][tDummy][eDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
204  eProbaShellMap[j][tDummy].push_back(eDiffCrossSectionData[j][tDummy][eDummy]);
205  }
206 
207  // SI - only if eof is not reached
208  if (!eDiffCrossSection.eof() && !fasterCode) eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
209 
210  if (!fasterCode) eVecm[tDummy].push_back(eDummy);
211 
212  }
213  }
214 
215  // *** PROTON
216 
217  proton = protonDef->GetParticleName();
218 
219  tableFile[proton] = fileProton;
220 
221  lowEnergyLimit[proton] = 500. * keV;
222  highEnergyLimit[proton] = 100. * MeV;
223 
224  // Cross section
225 
226  G4DNACrossSectionDataSet* tableP = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
227  tableP->LoadData(fileProton);
228 
229  tableData[proton] = tableP;
230 
231  // Final state
232 
233  std::ostringstream pFullFileName;
234 
235  if (fasterCode) pFullFileName << path << "/dna/sigmadiff_cumulated_ionisation_p_born_hp.dat";
236 
237  if (!fasterCode) pFullFileName << path << "/dna/sigmadiff_ionisation_p_born.dat";
238 
239  std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
240 
241  if (!pDiffCrossSection)
242  {
243  if (fasterCode) G4Exception("G4DNABornIonisationModel1::Initialise","em0003",
244  FatalException,"Missing data file:/dna/sigmadiff_cumulated_ionisation_p_born_hp.dat");
245 
246  if (!fasterCode) G4Exception("G4DNABornIonisationModel1::Initialise","em0003",
247  FatalException,"Missing data file:/dna/sigmadiff_ionisation_p_born.dat");
248  }
249 
250  pTdummyVec.push_back(0.);
251  while(!pDiffCrossSection.eof())
252  {
253  G4double tDummy;
254  G4double eDummy;
255  pDiffCrossSection>>tDummy>>eDummy;
256  if (tDummy != pTdummyVec.back()) pTdummyVec.push_back(tDummy);
257  for (G4int j=0; j<5; j++)
258  {
259  pDiffCrossSection>>pDiffCrossSectionData[j][tDummy][eDummy];
260 
261  if (fasterCode)
262  {
263  pNrjTransfData[j][tDummy][pDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
264  pProbaShellMap[j][tDummy].push_back(pDiffCrossSectionData[j][tDummy][eDummy]);
265  }
266 
267  // SI - only if eof is not reached !
268  if (!pDiffCrossSection.eof() && !fasterCode) pDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
269 
270  if (!fasterCode) pVecm[tDummy].push_back(eDummy);
271  }
272  }
273 
274  //
275 
276  if (particle==electronDef)
277  {
278  SetLowEnergyLimit(lowEnergyLimit[electron]);
279  SetHighEnergyLimit(highEnergyLimit[electron]);
280  }
281 
282  if (particle==protonDef)
283  {
284  SetLowEnergyLimit(lowEnergyLimit[proton]);
285  SetHighEnergyLimit(highEnergyLimit[proton]);
286  }
287 
288  if( verboseLevel>0 )
289  {
290  G4cout << "Born ionisation model is initialized " << G4endl
291  << "Energy range: "
292  << LowEnergyLimit() / eV << " eV - "
293  << HighEnergyLimit() / keV << " keV for "
294  << particle->GetParticleName()
295  << G4endl;
296  }
297 
298  // Initialize water density pointer
299 
300  fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->
301  GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
302 
303  // AD
304 
305  fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
306 
307  //
308 
309  if (isInitialised)
310  { return;}
311  fParticleChangeForGamma = GetParticleChangeForGamma();
312  isInitialised = true;
313 }
314 
315 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
316 
317 G4double G4DNABornIonisationModel1::CrossSectionPerVolume(const G4Material* material,
318  const G4ParticleDefinition* particleDefinition,
319  G4double ekin,
320  G4double,
321  G4double)
322 {
323  if (verboseLevel > 3)
324  {
325  G4cout << "Calling CrossSectionPerVolume() of G4DNABornIonisationModel1"
326  << G4endl;
327  }
328 
329  if (
330  particleDefinition != G4Proton::ProtonDefinition()
331  &&
332  particleDefinition != G4Electron::ElectronDefinition()
333  )
334 
335  return 0;
336 
337  // Calculate total cross section for model
338 
339  G4double lowLim = 0;
340  G4double highLim = 0;
341  G4double sigma=0;
342 
343  G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
344 
345  const G4String& particleName = particleDefinition->GetParticleName();
346 
347  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
348  pos1 = lowEnergyLimit.find(particleName);
349  if (pos1 != lowEnergyLimit.end())
350  {
351  lowLim = pos1->second;
352  }
353 
354  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
355  pos2 = highEnergyLimit.find(particleName);
356  if (pos2 != highEnergyLimit.end())
357  {
358  highLim = pos2->second;
359  }
360 
361  if (ekin >= lowLim && ekin <= highLim)
362  {
363  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
364  pos = tableData.find(particleName);
365 
366  if (pos != tableData.end())
367  {
368  G4DNACrossSectionDataSet* table = pos->second;
369  if (table != 0)
370  {
371  sigma = table->FindValue(ekin);
372 
373  // ICRU49 electronic SP scaling - ZF, SI
374 
375  if (particleDefinition == G4Proton::ProtonDefinition() && ekin < 70*MeV && spScaling)
376  {
377  G4double A = 1.39241700556072800000E-009 ;
378  G4double B = -8.52610412942622630000E-002 ;
379  sigma = sigma * G4Exp(A*(ekin/eV)+B);
380  }
381  //
382 
383  }
384  }
385  else
386  {
387  G4Exception("G4DNABornIonisationModel1::CrossSectionPerVolume","em0002",
388  FatalException,"Model not applicable to particle type.");
389  }
390  }
391 
392  if (verboseLevel > 2)
393  {
394  G4cout << "__________________________________" << G4endl;
395  G4cout << "G4DNABornIonisationModel1 - XS INFO START" << G4endl;
396  G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
397  G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
398  G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
399  G4cout << "G4DNABornIonisationModel1 - XS INFO END" << G4endl;
400  }
401 
402  return sigma*waterDensity;
403 }
404 
405 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
406 
407 void G4DNABornIonisationModel1::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
408  const G4MaterialCutsCouple* couple,
409  const G4DynamicParticle* particle,
410  G4double,
411  G4double)
412 {
413 
414  if (verboseLevel > 3)
415  {
416  G4cout << "Calling SampleSecondaries() of G4DNABornIonisationModel1"
417  << G4endl;
418  }
419 
420  G4double lowLim = 0;
421  G4double highLim = 0;
422 
423  G4double k = particle->GetKineticEnergy();
424 
425  const G4String& particleName = particle->GetDefinition()->GetParticleName();
426 
427  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
428  pos1 = lowEnergyLimit.find(particleName);
429 
430  if (pos1 != lowEnergyLimit.end())
431  {
432  lowLim = pos1->second;
433  }
434 
435  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
436  pos2 = highEnergyLimit.find(particleName);
437 
438  if (pos2 != highEnergyLimit.end())
439  {
440  highLim = pos2->second;
441  }
442 
443  if (k >= lowLim && k <= highLim)
444  {
445  G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
446  G4double particleMass = particle->GetDefinition()->GetPDGMass();
447  G4double totalEnergy = k + particleMass;
448  G4double pSquare = k * (totalEnergy + particleMass);
449  G4double totalMomentum = std::sqrt(pSquare);
450 
451  G4int ionizationShell = 0;
452 
453  if (!fasterCode) ionizationShell = RandomSelect(k,particleName);
454 
455  // SI: The following protection is necessary to avoid infinite loops :
456  // sigmadiff_ionisation_e_born.dat has non zero partial xs at 18 eV for shell 3 (ionizationShell ==2)
457  // sigmadiff_cumulated_ionisation_e_born.dat has zero cumulated partial xs at 18 eV for shell 3 (ionizationShell ==2)
458  // this is due to the fact that the max allowed transfered energy is (18+10.79)/2=17.025 eV and only transfered energies
459  // strictly above this value have non zero partial xs in sigmadiff_ionisation_e_born.dat (starting at trans = 17.12 eV)
460 
461  if (fasterCode)
462  do
463  {
464  ionizationShell = RandomSelect(k,particleName);
465  } while (k<19*eV && ionizationShell==2 && particle->GetDefinition()==G4Electron::ElectronDefinition());
466 
467  G4double bindingEnergy = 0;
468  bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
469 
470  // SI: additional protection if tcs interpolation method is modified
471  if (k<bindingEnergy) return;
472  //
473 
474  G4double secondaryKinetic=-1000*eV;
475 
476  if (fasterCode == false)
477  {
478  secondaryKinetic = RandomizeEjectedElectronEnergy(particle->GetDefinition(),k,ionizationShell);
479  }
480  else
481  {
482  secondaryKinetic = RandomizeEjectedElectronEnergyFromCumulatedDcs(particle->GetDefinition(),k,ionizationShell);
483  }
484  //
485 
486  G4int Z = 8;
487 
488  G4ThreeVector deltaDirection =
489  GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
490  Z, ionizationShell,
491  couple->GetMaterial());
492 
493  if (secondaryKinetic>0)
494  {
495  G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic);
496  fvect->push_back(dp);
497  }
498 
499  if (particle->GetDefinition() == G4Electron::ElectronDefinition())
500  {
501  G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
502 
503  G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
504  G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
505  G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
506  G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
507  finalPx /= finalMomentum;
508  finalPy /= finalMomentum;
509  finalPz /= finalMomentum;
510 
511  G4ThreeVector direction;
512  direction.set(finalPx,finalPy,finalPz);
513 
514  fParticleChangeForGamma->ProposeMomentumDirection(direction.unit());
515  }
516 
517  else fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
518 
519  // AM: sample deexcitation
520  // here we assume that H_{2}O electronic levels are the same as Oxygen.
521  // this can be considered true with a rough 10% error in energy on K-shell,
522 
523  size_t secNumberInit = 0;// need to know at a certain point the energy of secondaries
524  size_t secNumberFinal = 0;// So I'll make the diference and then sum the energies
525 
526  G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
527 
528  // SI: only atomic deexcitation from K shell is considered
529  if(fAtomDeexcitation && ionizationShell == 4)
530  {
531  const G4AtomicShell* shell =
532  fAtomDeexcitation->GetAtomicShell(Z, G4AtomicShellEnumerator(0));
533  secNumberInit = fvect->size();
534  fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
535  secNumberFinal = fvect->size();
536 
537  //TEST
538  //G4cout << "ionizationShell=" << ionizationShell<< G4endl;
539  //G4cout << "bindingEnergy=" << bindingEnergy/eV<< G4endl;
540 
541  if(secNumberFinal > secNumberInit)
542  {
543  for (size_t i=secNumberInit; i<secNumberFinal; ++i)
544  {
545  //Check if there is enough residual energy
546  if (bindingEnergy >= ((*fvect)[i])->GetKineticEnergy())
547  {
548  //Ok, this is a valid secondary: keep it
549  bindingEnergy -= ((*fvect)[i])->GetKineticEnergy();
550  //G4cout << "--deex nrj=" << ((*fvect)[i])->GetKineticEnergy()/eV
551  //<< G4endl;
552  }
553  else
554  {
555  //Invalid secondary: not enough energy to create it!
556  //Keep its energy in the local deposit
557  delete (*fvect)[i];
558  (*fvect)[i]=0;
559  }
560  }
561  }
562 
563  //TEST
564  //G4cout << "k=" << k/eV<< G4endl;
565  //G4cout << "secondaryKinetic=" << secondaryKinetic/eV<< G4endl;
566  //G4cout << "scatteredEnergy=" << scatteredEnergy/eV<< G4endl;
567  //G4cout << "deposited energy=" << bindingEnergy/eV<< G4endl;
568  //
569 
570  }
571 
572  //This should never happen
573  if(bindingEnergy < 0.0)
574  G4Exception("G4DNABornIonisatioModel1::SampleSecondaries()",
575  "em2050",FatalException,"Negative local energy deposit");
576 
577  //bindingEnergy has been decreased
578  //by the amount of energy taken away by deexc. products
579  if (!statCode)
580  {
581  fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
582  fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy);
583  }
584  else
585  {
586  fParticleChangeForGamma->SetProposedKineticEnergy(k);
587  fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy);
588  }
589 
590  // TEST //////////////////////////
591  // if (secondaryKinetic<0) abort();
592  // if (scatteredEnergy<0) abort();
593  // if (k-scatteredEnergy-secondaryKinetic-deexSecEnergy<0) abort();
594  // if (k-scatteredEnergy<0) abort();
596 
597  const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
598  G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
599  ionizationShell,
600  theIncomingTrack);
601  }
602 }
603 
604 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
605 
606 G4double G4DNABornIonisationModel1::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
607  G4double k,
608  G4int shell)
609 {
610  // G4cout << "*** SLOW computation for " << " " << particleDefinition->GetParticleName() << G4endl;
611 
612  if (particleDefinition == G4Electron::ElectronDefinition())
613  {
614  G4double maximumEnergyTransfer = 0.;
615  if ((k + waterStructure.IonisationEnergy(shell)) / 2. > k)
616  maximumEnergyTransfer = k;
617  else
618  maximumEnergyTransfer = (k + waterStructure.IonisationEnergy(shell)) / 2.;
619 
620  // SI : original method
621  /*
622  G4double crossSectionMaximum = 0.;
623  for(G4double value=waterStructure.IonisationEnergy(shell); value<=maximumEnergyTransfer; value+=0.1*eV)
624  {
625  G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
626  if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
627  }
628  */
629 
630  // SI : alternative method
631  G4double crossSectionMaximum = 0.;
632 
633  G4double minEnergy = waterStructure.IonisationEnergy(shell);
634  G4double maxEnergy = maximumEnergyTransfer;
635  G4int nEnergySteps = 50;
636 
637  G4double value(minEnergy);
638  G4double stpEnergy(std::pow(maxEnergy / value,
639  1. / static_cast<G4double>(nEnergySteps - 1)));
640  G4int step(nEnergySteps);
641  while (step > 0)
642  {
643  step--;
644  G4double differentialCrossSection =
645  DifferentialCrossSection(particleDefinition,
646  k / eV,
647  value / eV,
648  shell);
649  if (differentialCrossSection >= crossSectionMaximum)
650  crossSectionMaximum = differentialCrossSection;
651  value *= stpEnergy;
652  }
653  //
654 
655  G4double secondaryElectronKineticEnergy = 0.;
656  do
657  {
658  secondaryElectronKineticEnergy = G4UniformRand()* (maximumEnergyTransfer-waterStructure.IonisationEnergy(shell));
659  } while(G4UniformRand()*crossSectionMaximum >
660  DifferentialCrossSection(particleDefinition, k/eV,
661  (secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
662 
663  return secondaryElectronKineticEnergy;
664 
665  }
666 
667  else if (particleDefinition == G4Proton::ProtonDefinition())
668  {
669  G4double maximumKineticEnergyTransfer = 4.
670  * (electron_mass_c2 / proton_mass_c2) * k;
671 
672  G4double crossSectionMaximum = 0.;
673  for (G4double value = waterStructure.IonisationEnergy(shell);
674  value <= 4. * waterStructure.IonisationEnergy(shell); value += 0.1 * eV)
675  {
676  G4double differentialCrossSection =
677  DifferentialCrossSection(particleDefinition,
678  k / eV,
679  value / eV,
680  shell);
681  if (differentialCrossSection >= crossSectionMaximum)
682  crossSectionMaximum = differentialCrossSection;
683  }
684 
685  G4double secondaryElectronKineticEnergy = 0.;
686  do
687  {
688  secondaryElectronKineticEnergy = G4UniformRand()* maximumKineticEnergyTransfer;
689  } while(G4UniformRand()*crossSectionMaximum >=
690  DifferentialCrossSection(particleDefinition, k/eV,
691  (secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
692 
693  return secondaryElectronKineticEnergy;
694  }
695 
696  return 0;
697 }
698 
699 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
700 
701 // The following section is not used anymore but is kept for memory
702 // GetAngularDistribution()->SampleDirectionForShell is used instead
703 
704 /*
705  void G4DNABornIonisationModel1::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
706  G4double k,
707  G4double secKinetic,
708  G4double & cosTheta,
709  G4double & phi )
710  {
711  if (particleDefinition == G4Electron::ElectronDefinition())
712  {
713  phi = twopi * G4UniformRand();
714  if (secKinetic < 50.*eV) cosTheta = (2.*G4UniformRand())-1.;
715  else if (secKinetic <= 200.*eV)
716  {
717  if (G4UniformRand() <= 0.1) cosTheta = (2.*G4UniformRand())-1.;
718  else cosTheta = G4UniformRand()*(std::sqrt(2.)/2);
719  }
720  else
721  {
722  G4double sin2O = (1.-secKinetic/k) / (1.+secKinetic/(2.*electron_mass_c2));
723  cosTheta = std::sqrt(1.-sin2O);
724  }
725  }
726 
727  else if (particleDefinition == G4Proton::ProtonDefinition())
728  {
729  G4double maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
730  phi = twopi * G4UniformRand();
731 
732  // cosTheta = std::sqrt(secKinetic / maxSecKinetic);
733 
734  // Restriction below 100 eV from Emfietzoglou (2000)
735 
736  if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
737  else cosTheta = (2.*G4UniformRand())-1.;
738 
739  }
740  }
741  */
742 
743 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
744 G4double G4DNABornIonisationModel1::DifferentialCrossSection(G4ParticleDefinition * particleDefinition,
745  G4double k,
746  G4double energyTransfer,
747  G4int ionizationLevelIndex)
748 {
749  G4double sigma = 0.;
750 
751  if (energyTransfer >= waterStructure.IonisationEnergy(ionizationLevelIndex)/eV)
752  {
753  G4double valueT1 = 0;
754  G4double valueT2 = 0;
755  G4double valueE21 = 0;
756  G4double valueE22 = 0;
757  G4double valueE12 = 0;
758  G4double valueE11 = 0;
759 
760  G4double xs11 = 0;
761  G4double xs12 = 0;
762  G4double xs21 = 0;
763  G4double xs22 = 0;
764 
765  if (particleDefinition == G4Electron::ElectronDefinition())
766  {
767 
768  // Protection against out of boundary access - electron case : 1 MeV
769  if (k==eTdummyVec.back()) k=k*(1.-1e-12);
770  //
771 
772  // k should be in eV and energy transfer eV also
773 
774  std::vector<G4double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),
775  eTdummyVec.end(),
776  k);
777 
778  std::vector<G4double>::iterator t1 = t2 - 1;
779 
780  // SI : the following condition avoids situations where energyTransfer >last vector element
781  if (energyTransfer <= eVecm[(*t1)].back()
782  && energyTransfer <= eVecm[(*t2)].back())
783  {
784  std::vector<G4double>::iterator e12 =
785  std::upper_bound(eVecm[(*t1)].begin(),
786  eVecm[(*t1)].end(),
787  energyTransfer);
788  std::vector<G4double>::iterator e11 = e12 - 1;
789 
790  std::vector<G4double>::iterator e22 =
791  std::upper_bound(eVecm[(*t2)].begin(),
792  eVecm[(*t2)].end(),
793  energyTransfer);
794  std::vector<G4double>::iterator e21 = e22 - 1;
795 
796  valueT1 = *t1;
797  valueT2 = *t2;
798  valueE21 = *e21;
799  valueE22 = *e22;
800  valueE12 = *e12;
801  valueE11 = *e11;
802 
803  xs11 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
804  xs12 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
805  xs21 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
806  xs22 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
807 
808  }
809 
810  }
811 
812  if (particleDefinition == G4Proton::ProtonDefinition())
813  {
814  // Protection against out of boundary access - proton case : 100 MeV
815  if (k==pTdummyVec.back()) k=k*(1.-1e-12);
816  //
817 
818  // k should be in eV and energy transfer eV also
819  std::vector<G4double>::iterator t2 = std::upper_bound(pTdummyVec.begin(),
820  pTdummyVec.end(),
821  k);
822  std::vector<G4double>::iterator t1 = t2 - 1;
823 
824  std::vector<G4double>::iterator e12 = std::upper_bound(pVecm[(*t1)].begin(),
825  pVecm[(*t1)].end(),
826  energyTransfer);
827  std::vector<G4double>::iterator e11 = e12 - 1;
828 
829  std::vector<G4double>::iterator e22 = std::upper_bound(pVecm[(*t2)].begin(),
830  pVecm[(*t2)].end(),
831  energyTransfer);
832  std::vector<G4double>::iterator e21 = e22 - 1;
833 
834  valueT1 = *t1;
835  valueT2 = *t2;
836  valueE21 = *e21;
837  valueE22 = *e22;
838  valueE12 = *e12;
839  valueE11 = *e11;
840 
841  xs11 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
842  xs12 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
843  xs21 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
844  xs22 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
845 
846  }
847 
848  G4double xsProduct = xs11 * xs12 * xs21 * xs22;
849  if (xsProduct != 0.)
850  {
851  sigma = QuadInterpolator(valueE11,
852  valueE12,
853  valueE21,
854  valueE22,
855  xs11,
856  xs12,
857  xs21,
858  xs22,
859  valueT1,
860  valueT2,
861  k,
862  energyTransfer);
863  }
864  }
865 
866  return sigma;
867 }
868 
869 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
870 
871 G4double G4DNABornIonisationModel1::Interpolate(G4double e1,
872  G4double e2,
873  G4double e,
874  G4double xs1,
875  G4double xs2)
876 {
877  G4double value = 0.;
878 
879  // Log-log interpolation by default
880 
881  if (e1 != 0 && e2 != 0 && (std::log10(e2) - std::log10(e1)) != 0
882  && !fasterCode)
883  {
884  G4double a = (std::log10(xs2) - std::log10(xs1))
885  / (std::log10(e2) - std::log10(e1));
886  G4double b = std::log10(xs2) - a * std::log10(e2);
887  G4double sigma = a * std::log10(e) + b;
888  value = (std::pow(10., sigma));
889  }
890 
891  // Switch to lin-lin interpolation
892  /*
893  if ((e2-e1)!=0)
894  {
895  G4double d1 = xs1;
896  G4double d2 = xs2;
897  value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
898  }
899  */
900 
901  // Switch to log-lin interpolation for faster code
902  if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 && fasterCode)
903  {
904  G4double d1 = std::log10(xs1);
905  G4double d2 = std::log10(xs2);
906  value = std::pow(10., (d1 + (d2 - d1) * (e - e1) / (e2 - e1)));
907  }
908 
909  // Switch to lin-lin interpolation for faster code
910  // in case one of xs1 or xs2 (=cum proba) value is zero
911 
912  if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0) && fasterCode)
913  {
914  G4double d1 = xs1;
915  G4double d2 = xs2;
916  value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
917  }
918 
919  /*
920  G4cout
921  << e1 << " "
922  << e2 << " "
923  << e << " "
924  << xs1 << " "
925  << xs2 << " "
926  << value
927  << G4endl;
928  */
929 
930  return value;
931 }
932 
933 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
934 
935 G4double G4DNABornIonisationModel1::QuadInterpolator(G4double e11,
936  G4double e12,
937  G4double e21,
938  G4double e22,
939  G4double xs11,
940  G4double xs12,
941  G4double xs21,
942  G4double xs22,
943  G4double t1,
944  G4double t2,
945  G4double t,
946  G4double e)
947 {
948  G4double interpolatedvalue1 = Interpolate(e11, e12, e, xs11, xs12);
949  G4double interpolatedvalue2 = Interpolate(e21, e22, e, xs21, xs22);
950  G4double value = Interpolate(t1,
951  t2,
952  t,
953  interpolatedvalue1,
954  interpolatedvalue2);
955 
956  return value;
957 }
958 
959 G4double G4DNABornIonisationModel1::GetPartialCrossSection(const G4Material* /*material*/,
960  G4int level,
961  const G4ParticleDefinition* particle,
962  G4double kineticEnergy)
963 {
964  std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
965  pos = tableData.find(particle->GetParticleName());
966 
967  if (pos != tableData.end())
968  {
969  G4DNACrossSectionDataSet* table = pos->second;
970  return table->GetComponent(level)->FindValue(kineticEnergy);
971  }
972 
973  return 0;
974 }
975 
976 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
977 
978 G4int G4DNABornIonisationModel1::RandomSelect(G4double k,
979  const G4String& particle)
980 {
981  G4int level = 0;
982 
983  std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
984  pos = tableData.find(particle);
985 
986  if (pos != tableData.end())
987  {
988  G4DNACrossSectionDataSet* table = pos->second;
989 
990  if (table != 0)
991  {
992  G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
993  const size_t n(table->NumberOfComponents());
994  size_t i(n);
995  G4double value = 0.;
996 
997  while (i > 0)
998  {
999  i--;
1000  valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
1001  value += valuesBuffer[i];
1002  }
1003 
1004  value *= G4UniformRand();
1005 
1006  i = n;
1007 
1008  while (i > 0)
1009  {
1010  i--;
1011 
1012  if (valuesBuffer[i] > value)
1013  {
1014  delete[] valuesBuffer;
1015  return i;
1016  }
1017  value -= valuesBuffer[i];
1018  }
1019 
1020  if (valuesBuffer)
1021  delete[] valuesBuffer;
1022 
1023  }
1024  } else
1025  {
1026  G4Exception("G4DNABornIonisationModel1::RandomSelect",
1027  "em0002",
1028  FatalException,
1029  "Model not applicable to particle type.");
1030  }
1031 
1032  return level;
1033 }
1034 
1035 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1036 
1037 G4double G4DNABornIonisationModel1::RandomizeEjectedElectronEnergyFromCumulatedDcs(G4ParticleDefinition* particleDefinition,
1038  G4double k,
1039  G4int shell)
1040 {
1041  //G4cout << "*** FAST computation for " << " " << particleDefinition->GetParticleName() << G4endl;
1042 
1043  G4double secondaryElectronKineticEnergy = 0.;
1044 
1045  G4double random = G4UniformRand();
1046 
1047  secondaryElectronKineticEnergy = TransferedEnergy(particleDefinition,
1048  k / eV,
1049  shell,
1050  random) * eV
1051  - waterStructure.IonisationEnergy(shell);
1052 
1053  //G4cout << RandomTransferedEnergy(particleDefinition, k/eV, shell) << G4endl;
1054 
1055  if (secondaryElectronKineticEnergy < 0.)
1056  return 0.;
1057  //
1058 
1059  return secondaryElectronKineticEnergy;
1060 }
1061 
1062 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1063 
1064 G4double G4DNABornIonisationModel1::TransferedEnergy(G4ParticleDefinition* particleDefinition,
1065  G4double k,
1066  G4int ionizationLevelIndex,
1067  G4double random)
1068 {
1069  G4double nrj = 0.;
1070 
1071  G4double valueK1 = 0;
1072  G4double valueK2 = 0;
1073  G4double valuePROB21 = 0;
1074  G4double valuePROB22 = 0;
1075  G4double valuePROB12 = 0;
1076  G4double valuePROB11 = 0;
1077 
1078  G4double nrjTransf11 = 0;
1079  G4double nrjTransf12 = 0;
1080  G4double nrjTransf21 = 0;
1081  G4double nrjTransf22 = 0;
1082 
1083  if (particleDefinition == G4Electron::ElectronDefinition())
1084  {
1085  // Protection against out of boundary access - electron case : 1 MeV
1086  if (k==eTdummyVec.back()) k=k*(1.-1e-12);
1087  //
1088 
1089  // k should be in eV
1090  std::vector<G4double>::iterator k2 = std::upper_bound(eTdummyVec.begin(),
1091  eTdummyVec.end(),
1092  k);
1093  std::vector<G4double>::iterator k1 = k2 - 1;
1094 
1095  /*
1096  G4cout << "----> k=" << k
1097  << " " << *k1
1098  << " " << *k2
1099  << " " << random
1100  << " " << ionizationLevelIndex
1101  << " " << eProbaShellMap[ionizationLevelIndex][(*k1)].back()
1102  << " " << eProbaShellMap[ionizationLevelIndex][(*k2)].back()
1103  << G4endl;
1104  */
1105 
1106  // SI : the following condition avoids situations where random >last vector element
1107  if (random <= eProbaShellMap[ionizationLevelIndex][(*k1)].back()
1108  && random <= eProbaShellMap[ionizationLevelIndex][(*k2)].back())
1109  {
1110  std::vector<G4double>::iterator prob12 =
1111  std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
1112  eProbaShellMap[ionizationLevelIndex][(*k1)].end(),
1113  random);
1114 
1115  std::vector<G4double>::iterator prob11 = prob12 - 1;
1116 
1117  std::vector<G4double>::iterator prob22 =
1118  std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1119  eProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1120  random);
1121 
1122  std::vector<G4double>::iterator prob21 = prob22 - 1;
1123 
1124  valueK1 = *k1;
1125  valueK2 = *k2;
1126  valuePROB21 = *prob21;
1127  valuePROB22 = *prob22;
1128  valuePROB12 = *prob12;
1129  valuePROB11 = *prob11;
1130 
1131  /*
1132  G4cout << " " << random << " " << valuePROB11 << " "
1133  << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1134  */
1135 
1136  nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1137  nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1138  nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1139  nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1140 
1141  /*
1142  G4cout << " " << ionizationLevelIndex << " "
1143  << random << " " <<valueK1 << " " << valueK2 << G4endl;
1144 
1145  G4cout << " " << random << " " << nrjTransf11 << " "
1146  << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1147  */
1148 
1149  }
1150 
1151  // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
1152  if (random > eProbaShellMap[ionizationLevelIndex][(*k1)].back())
1153  {
1154  std::vector<G4double>::iterator prob22 =
1155  std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1156  eProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1157  random);
1158 
1159  std::vector<G4double>::iterator prob21 = prob22 - 1;
1160 
1161  valueK1 = *k1;
1162  valueK2 = *k2;
1163  valuePROB21 = *prob21;
1164  valuePROB22 = *prob22;
1165 
1166  //G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1167 
1168  nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1169  nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1170 
1171  G4double interpolatedvalue2 = Interpolate(valuePROB21,
1172  valuePROB22,
1173  random,
1174  nrjTransf21,
1175  nrjTransf22);
1176 
1177  // zeros are explicitely set
1178 
1179  G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
1180 
1181  /*
1182  G4cout << " " << ionizationLevelIndex << " "
1183  << random << " " <<valueK1 << " " << valueK2 << G4endl;
1184 
1185  G4cout << " " << random << " " << nrjTransf11 << " "
1186  << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1187 
1188  G4cout << "ici" << " " << value << G4endl;
1189  */
1190 
1191  return value;
1192  }
1193  }
1194  //
1195  else if (particleDefinition == G4Proton::ProtonDefinition())
1196  {
1197  // Protection against out of boundary access - proton case : 100 MeV
1198  if (k==pTdummyVec.back()) k=k*(1.-1e-12);
1199  //
1200 
1201  // k should be in eV
1202 
1203  std::vector<G4double>::iterator k2 = std::upper_bound(pTdummyVec.begin(),
1204  pTdummyVec.end(),
1205  k);
1206 
1207  std::vector<G4double>::iterator k1 = k2 - 1;
1208 
1209  /*
1210  G4cout << "----> k=" << k
1211  << " " << *k1
1212  << " " << *k2
1213  << " " << random
1214  << " " << ionizationLevelIndex
1215  << " " << pProbaShellMap[ionizationLevelIndex][(*k1)].back()
1216  << " " << pProbaShellMap[ionizationLevelIndex][(*k2)].back()
1217  << G4endl;
1218  */
1219 
1220  // SI : the following condition avoids situations where random > last vector element,
1221  // for eg. when the last element is zero
1222  if (random <= pProbaShellMap[ionizationLevelIndex][(*k1)].back()
1223  && random <= pProbaShellMap[ionizationLevelIndex][(*k2)].back())
1224  {
1225  std::vector<G4double>::iterator prob12 =
1226  std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
1227  pProbaShellMap[ionizationLevelIndex][(*k1)].end(),
1228  random);
1229 
1230  std::vector<G4double>::iterator prob11 = prob12 - 1;
1231 
1232  std::vector<G4double>::iterator prob22 =
1233  std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1234  pProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1235  random);
1236 
1237  std::vector<G4double>::iterator prob21 = prob22 - 1;
1238 
1239  valueK1 = *k1;
1240  valueK2 = *k2;
1241  valuePROB21 = *prob21;
1242  valuePROB22 = *prob22;
1243  valuePROB12 = *prob12;
1244  valuePROB11 = *prob11;
1245 
1246  /*
1247  G4cout << " " << random << " " << valuePROB11 << " "
1248  << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1249  */
1250 
1251  nrjTransf11 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1252  nrjTransf12 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1253  nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1254  nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1255 
1256  /*
1257  G4cout << " " << ionizationLevelIndex << " "
1258  << random << " " <<valueK1 << " " << valueK2 << G4endl;
1259 
1260  G4cout << " " << random << " " << nrjTransf11 << " "
1261  << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1262  */
1263  }
1264 
1265  // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
1266 
1267  if (random > pProbaShellMap[ionizationLevelIndex][(*k1)].back())
1268  {
1269  std::vector<G4double>::iterator prob22 =
1270  std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1271  pProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1272  random);
1273 
1274  std::vector<G4double>::iterator prob21 = prob22 - 1;
1275 
1276  valueK1 = *k1;
1277  valueK2 = *k2;
1278  valuePROB21 = *prob21;
1279  valuePROB22 = *prob22;
1280 
1281  //G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1282 
1283  nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1284  nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1285 
1286  G4double interpolatedvalue2 = Interpolate(valuePROB21,
1287  valuePROB22,
1288  random,
1289  nrjTransf21,
1290  nrjTransf22);
1291 
1292  // zeros are explicitely set
1293 
1294  G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
1295 
1296  /*
1297  G4cout << " " << ionizationLevelIndex << " "
1298  << random << " " <<valueK1 << " " << valueK2 << G4endl;
1299 
1300  G4cout << " " << random << " " << nrjTransf11 << " "
1301  << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1302 
1303  G4cout << "ici" << " " << value << G4endl;
1304  */
1305 
1306  return value;
1307  }
1308  }
1309  // End electron and proton cases
1310 
1311  G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21
1312  * nrjTransf22;
1313  //G4cout << "nrjTransfProduct=" << nrjTransfProduct << G4endl;
1314 
1315  if (nrjTransfProduct != 0.)
1316  {
1317  nrj = QuadInterpolator(valuePROB11,
1318  valuePROB12,
1319  valuePROB21,
1320  valuePROB22,
1321  nrjTransf11,
1322  nrjTransf12,
1323  nrjTransf21,
1324  nrjTransf22,
1325  valueK1,
1326  valueK2,
1327  k,
1328  random);
1329  }
1330  //G4cout << nrj << endl;
1331 
1332  return nrj;
1333 }