ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNAPTBIonisationModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DNAPTBIonisationModel.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 // Authors: S. Meylan and C. Villagrasa (IRSN, France)
27 // Models come from
28 // M. Bug et al, Rad. Phys and Chem. 130, 459-479 (2017)
29 //
30 
32 #include "G4PhysicalConstants.hh"
33 #include "G4SystemOfUnits.hh"
34 #include "G4UAtomicDeexcitation.hh"
35 #include "G4LossTableManager.hh"
36 #include "G4DNAChemistryManager.hh"
37 
39  const G4ParticleDefinition*,
40  const G4String& nam, const G4bool isAuger)
41  : G4VDNAModel(nam, applyToMaterial)
42 {
43  verboseLevel= 0;
44  // Verbosity scale:
45  // 0 = nothing
46  // 1 = warning for energy non-conservation
47  // 2 = details of energy budget
48  // 3 = calculation of cross sections, file openings, sampling of atoms
49  // 4 = entering in methods
50 
51  if( verboseLevel>0 )
52  {
53  G4cout << "PTB ionisation model is constructed " << G4endl;
54  }
55 
56  if(isAuger)
57  {
58  // create the PTB Auger model
59  fDNAPTBAugerModel = new G4DNAPTBAugerModel("e-_G4DNAPTBAugerModel");
60  }
61  else
62  {
63  // no PTB Auger model
65  }
66 }
67 
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
69 
71 {
72  // To delete the DNAPTBAugerModel created at initialisation of the ionisation class
74 }
75 
76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77 
79  const G4DataVector& /*cuts*/, G4ParticleChangeForGamma*)
80 {
81 
82  if (verboseLevel > 3)
83  G4cout << "Calling G4DNAPTBIonisationModel::Initialise()" << G4endl;
84 
85  G4double scaleFactor = 1e-16 * cm*cm;
86  G4double scaleFactorBorn = (1.e-22 / 3.343) * m*m;
87 
90 
91  //*******************************************************
92  // Cross section data
93  //*******************************************************
94 
95  if(particle == electronDef)
96  {
97  G4String particleName = particle->GetParticleName();
98 
99  // Raw materials
100  //
101  AddCrossSectionData("THF",
102  particleName,
103  "dna/sigma_ionisation_e-_PTB_THF",
104  "dna/sigmadiff_cumulated_ionisation_e-_PTB_THF",
105  scaleFactor);
106  SetLowELimit("THF", particleName, 12.*eV);
107  SetHighELimit("THF", particleName, 1.*keV);
108 
109  AddCrossSectionData("PY",
110  particleName,
111  "dna/sigma_ionisation_e-_PTB_PY",
112  "dna/sigmadiff_cumulated_ionisation_e-_PTB_PY",
113  scaleFactor);
114  SetLowELimit("PY", particleName, 12.*eV);
115  SetHighELimit("PY", particleName, 1.*keV);
116 
117  AddCrossSectionData("PU",
118  particleName,
119  "dna/sigma_ionisation_e-_PTB_PU",
120  "dna/sigmadiff_cumulated_ionisation_e-_PTB_PU",
121  scaleFactor);
122  SetLowELimit("PU", particleName, 12.*eV);
123  SetHighELimit("PU", particleName, 1.*keV);
124 
125  AddCrossSectionData("TMP",
126  particleName,
127  "dna/sigma_ionisation_e-_PTB_TMP",
128  "dna/sigmadiff_cumulated_ionisation_e-_PTB_TMP",
129  scaleFactor);
130  SetLowELimit("TMP", particleName, 12.*eV);
131  SetHighELimit("TMP", particleName, 1.*keV);
132 
133  AddCrossSectionData("G4_WATER",
134  particleName,
135  "dna/sigma_ionisation_e_born",
136  "dna/sigmadiff_ionisation_e_born",
137  scaleFactorBorn);
138  SetLowELimit("G4_WATER", particleName, 12.*eV);
139  SetHighELimit("G4_WATER", particleName, 1.*keV);
140 
141  // DNA materials
142  //
143  AddCrossSectionData("backbone_THF",
144  particleName,
145  "dna/sigma_ionisation_e-_PTB_THF",
146  "dna/sigmadiff_cumulated_ionisation_e-_PTB_THF",
147  scaleFactor*33./30);
148  SetLowELimit("backbone_THF", particleName, 12.*eV);
149  SetHighELimit("backbone_THF", particleName, 1.*keV);
150 
151  AddCrossSectionData("cytosine_PY",
152  particleName,
153  "dna/sigma_ionisation_e-_PTB_PY",
154  "dna/sigmadiff_cumulated_ionisation_e-_PTB_PY",
155  scaleFactor*42./30);
156  SetLowELimit("cytosine_PY", particleName, 12.*eV);
157  SetHighELimit("cytosine_PY", particleName, 1.*keV);
158 
159  AddCrossSectionData("thymine_PY",
160  particleName,
161  "dna/sigma_ionisation_e-_PTB_PY",
162  "dna/sigmadiff_cumulated_ionisation_e-_PTB_PY",
163  scaleFactor*48./30);
164  SetLowELimit("thymine_PY", particleName, 12.*eV);
165  SetHighELimit("thymine_PY", particleName, 1.*keV);
166 
167  AddCrossSectionData("adenine_PU",
168  particleName,
169  "dna/sigma_ionisation_e-_PTB_PU",
170  "dna/sigmadiff_cumulated_ionisation_e-_PTB_PU",
171  scaleFactor*50./44);
172  SetLowELimit("adenine_PU", particleName, 12.*eV);
173  SetHighELimit("adenine_PU", particleName, 1.*keV);
174 
175  AddCrossSectionData("guanine_PU",
176  particleName,
177  "dna/sigma_ionisation_e-_PTB_PU",
178  "dna/sigmadiff_cumulated_ionisation_e-_PTB_PU",
179  scaleFactor*56./44);
180  SetLowELimit("guanine_PU", particleName, 12.*eV);
181  SetHighELimit("guanine_PU", particleName, 1.*keV);
182 
183  AddCrossSectionData("backbone_TMP",
184  particleName,
185  "dna/sigma_ionisation_e-_PTB_TMP",
186  "dna/sigmadiff_cumulated_ionisation_e-_PTB_TMP",
187  scaleFactor*33./50);
188  SetLowELimit("backbone_TMP", particleName, 12.*eV);
189  SetHighELimit("backbone_TMP", particleName, 1.*keV);
190  }
191 
192  else if (particle == protonDef)
193  {
194  G4String particleName = particle->GetParticleName();
195 
196  // Raw materials
197  //
198  AddCrossSectionData("THF",
199  particleName,
200  "dna/sigma_ionisation_p_HKS_THF",
201  "dna/sigmadiff_cumulated_ionisation_p_PTB_THF",
202  scaleFactor);
203  SetLowELimit("THF", particleName, 70.*keV);
204  SetHighELimit("THF", particleName, 10.*MeV);
205 
206 
207  AddCrossSectionData("PY",
208  particleName,
209  "dna/sigma_ionisation_p_HKS_PY",
210  "dna/sigmadiff_cumulated_ionisation_p_PTB_PY",
211  scaleFactor);
212  SetLowELimit("PY", particleName, 70.*keV);
213  SetHighELimit("PY", particleName, 10.*MeV);
214 
215  /*
216  AddCrossSectionData("PU",
217  particleName,
218  "dna/sigma_ionisation_e-_PTB_PU",
219  "dna/sigmadiff_cumulated_ionisation_e-_PTB_PU",
220  scaleFactor);
221  SetLowELimit("PU", particleName2, 70.*keV);
222  SetHighELimit("PU", particleName2, 10.*keV);
223 */
224 
225  AddCrossSectionData("TMP",
226  particleName,
227  "dna/sigma_ionisation_p_HKS_TMP",
228  "dna/sigmadiff_cumulated_ionisation_p_PTB_TMP",
229  scaleFactor);
230  SetLowELimit("TMP", particleName, 70.*keV);
231  SetHighELimit("TMP", particleName, 10.*MeV);
232  }
233 
234  // *******************************************************
235  // deal with composite materials
236  // *******************************************************
237 
239 
240  // *******************************************************
241  // Verbose
242  // *******************************************************
243 
244  // initialise DNAPTBAugerModel
246 }
247 
248 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
249 
251  const G4String& materialName,
252  const G4ParticleDefinition* p,
253  G4double ekin,
254  G4double /*emin*/,
255  G4double /*emax*/)
256 {
257  if(verboseLevel > 3)
258  G4cout << "Calling CrossSectionPerVolume() of G4DNAPTBIonisationModel" << G4endl;
259 
260  // initialise the cross section value (output value)
261  G4double sigma(0);
262 
263  // Get the current particle name
264  const G4String& particleName = p->GetParticleName();
265 
266  // Set the low and high energy limits
267  G4double lowLim = GetLowELimit(materialName, particleName);
268  G4double highLim = GetHighELimit(materialName, particleName);
269 
270  // Check that we are in the correct energy range
271  if (ekin >= lowLim && ekin < highLim)
272  {
273  // Get the map with all the model data tables
274  TableMapData* tableData = GetTableData();
275 
276  // Retrieve the cross section value for the current material, particle and energy values
277  sigma = (*tableData)[materialName][particleName]->FindValue(ekin);
278 
279  if (verboseLevel > 2)
280  {
281  G4cout << "__________________________________" << G4endl;
282  G4cout << "°°° G4DNAPTBIonisationModel - XS INFO START" << G4endl;
283  G4cout << "°°° Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
284  G4cout << "°°° Cross section per "<< materialName <<" molecule (cm^2)=" << sigma/cm/cm << G4endl;
285  G4cout << "°°° G4DNAPTBIonisationModel - XS INFO END" << G4endl;
286  }
287  }
288 
289  // Return the cross section value
290  return sigma;
291 }
292 
293 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
294 
295 void G4DNAPTBIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
296  const G4MaterialCutsCouple* /*couple*/,
297  const G4String& materialName,
298  const G4DynamicParticle* aDynamicParticle,
299  G4ParticleChangeForGamma* particleChangeForGamma,
300  G4double /*tmin*/,
301  G4double /*tmax*/)
302 {
303  if (verboseLevel > 3)
304  G4cout << "Calling SampleSecondaries() of G4DNAPTBIonisationModel" << G4endl;
305 
306  // Get the current particle energy
307  G4double k = aDynamicParticle->GetKineticEnergy();
308 
309  // Get the current particle name
310  const G4String& particleName = aDynamicParticle->GetDefinition()->GetParticleName();
311 
312  // Get the energy limits
313  G4double lowLim = GetLowELimit(materialName, particleName);
314  G4double highLim = GetHighELimit(materialName, particleName);
315 
316  // Check if we are in the correct energy range
317  if (k >= lowLim && k < highLim)
318  {
319  G4ParticleMomentum primaryDirection = aDynamicParticle->GetMomentumDirection();
320  G4double particleMass = aDynamicParticle->GetDefinition()->GetPDGMass();
321  G4double totalEnergy = k + particleMass;
322  G4double pSquare = k * (totalEnergy + particleMass);
323  G4double totalMomentum = std::sqrt(pSquare);
324 
325  // Get the ionisation shell from a random sampling
326  G4int ionizationShell = RandomSelectShell(k, particleName, materialName);
327 
328  // Get the binding energy from the ptbStructure class
329  G4double bindingEnergy = ptbStructure.IonisationEnergy(ionizationShell, materialName);
330 
331  // Initialize the secondary kinetic energy to a negative value.
332  G4double secondaryKinetic (-1000*eV);
333 
334  if(materialName!="G4_WATER")
335  {
336  // Get the energy of the secondary particle
337  secondaryKinetic = RandomizeEjectedElectronEnergyFromCumulated(aDynamicParticle->GetDefinition(),k/eV,ionizationShell, materialName);
338  }
339  else
340  {
341  secondaryKinetic = RandomizeEjectedElectronEnergy(aDynamicParticle->GetDefinition(),k,ionizationShell, materialName);
342  }
343 
344  if(secondaryKinetic<=0)
345  {
346  G4cout<<"Fatal error *************************************** "<<secondaryKinetic/eV<<G4endl;
347  G4cout<<"secondaryKinetic: "<<secondaryKinetic/eV<<G4endl;
348  G4cout<<"k: "<<k/eV<<G4endl;
349  G4cout<<"shell: "<<ionizationShell<<G4endl;
350  G4cout<<"material:"<<materialName<<G4endl;
351  exit(EXIT_FAILURE);
352  }
353 
354  G4double cosTheta = 0.;
355  G4double phi = 0.;
356  RandomizeEjectedElectronDirection(aDynamicParticle->GetDefinition(), k, secondaryKinetic, cosTheta, phi);
357 
358  G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
359  G4double dirX = sinTheta*std::cos(phi);
360  G4double dirY = sinTheta*std::sin(phi);
361  G4double dirZ = cosTheta;
362  G4ThreeVector deltaDirection(dirX,dirY,dirZ);
363  deltaDirection.rotateUz(primaryDirection);
364 
365  // The model is written only for electron and thus we want the change the direction of the incident electron
366  // after each ionization. However, if other particle are going to be introduced within this model the following should be added:
367  //
368  // Check if the particle is an electron
369  if(aDynamicParticle->GetDefinition() == G4Electron::ElectronDefinition() )
370  {
371  // If yes do the following code until next commented "else" statement
372 
373  G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
374  G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
375  G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
376  G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
377  G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
378  finalPx /= finalMomentum;
379  finalPy /= finalMomentum;
380  finalPz /= finalMomentum;
381 
382  G4ThreeVector direction(finalPx,finalPy,finalPz);
383  if(direction.unit().getX()>1||direction.unit().getY()>1||direction.unit().getZ()>1)
384  {
385  G4cout<<"Fatal error ****************************"<<G4endl;
386  G4cout<<"direction problem "<<direction.unit()<<G4endl;
387  exit(EXIT_FAILURE);
388  }
389 
390  // Give the new direction to the particle
391  particleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
392  }
393  // If the particle is not an electron
394  else particleChangeForGamma->ProposeMomentumDirection(primaryDirection) ;
395 
396  // note that secondaryKinetic is the energy of the delta ray, not of all secondaries.
397  G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
398 
399  if(scatteredEnergy<=0)
400  {
401  G4cout<<"Fatal error ****************************"<<G4endl;
402  G4cout<<"k: "<<k/eV<<G4endl;
403  G4cout<<"secondaryKinetic: "<<secondaryKinetic/eV<<G4endl;
404  G4cout<<"shell: "<<ionizationShell<<G4endl;
405  G4cout<<"bindingEnergy: "<<bindingEnergy/eV<<G4endl;
406  G4cout<<"scatteredEnergy: "<<scatteredEnergy/eV<<G4endl;
407  G4cout<<"material: "<<materialName<<G4endl;
408  exit(EXIT_FAILURE);
409  }
410 
411  // Set the new energy of the particle
412  particleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
413 
414  // Set the energy deposited by the ionization
415  particleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy-secondaryKinetic);
416 
417  // Create the new particle with its characteristics
418  G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
419  fvect->push_back(dp);
420 
421  // Check if the auger model is activated (ie instanciated)
423  {
424  // run the PTB Auger model
425  if(materialName!="G4_WATER") fDNAPTBAugerModel->ComputeAugerEffect(fvect, materialName, bindingEnergy);
426  }
427  }
428 }
429 
430 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
431 
433  const G4String& particleName,
434  const G4String& file,
435  const G4double scaleFactor)
436 {
437  // To read and save the informations contained within the differential cross section files
438 
439  // get the path of the G4LEDATA data folder
440  char *path = std::getenv("G4LEDATA");
441  // if it is not found then quit and print error message
442  if(!path)
443  {
444  G4Exception("G4DNAPTBIonisationModel::ReadAllDiffCSFiles","em0006",
445  FatalException,"G4LEDATA environment variable not set.");
446  return;
447  }
448 
449  // build the fullFileName path of the data file
450  std::ostringstream fullFileName;
451  fullFileName << path <<"/"<< file<<".dat";
452 
453  // open the data file
454  std::ifstream diffCrossSection (fullFileName.str().c_str());
455  // error if file is not there
456  std::stringstream endPath;
457  if (!diffCrossSection)
458  {
459  endPath << "Missing data file: "<<file;
460  G4Exception("G4DNAPTBIonisationModel::Initialise","em0003",
461  FatalException, endPath.str().c_str());
462  }
463 
464  // load data from the file
465  fTMapWithVec[materialName][particleName].push_back(0.);
466 
467  G4String line;
468 
469  // read the file until we reach the end of file point
470  // fill fTMapWithVec, diffCrossSectionData, fEnergyTransferData, fProbaShellMap and fEMapWithVector
471  while(std::getline(diffCrossSection, line))
472  {
473  // check if the line is comment or empty
474  //
475  std::istringstream testIss(line);
476  G4String test;
477  testIss >> test;
478  // check first caracter to determine if following information is data or comments
479  if(test=="#")
480  {
481  // skip the line by beginning a new while loop.
482  continue;
483  }
484  // check if line is empty
485  else if(line.empty())
486  {
487  // skip the line by beginning a new while loop.
488  continue;
489  }
490  //
491  // end of the check
492 
493  // transform the line into a iss
494  std::istringstream iss(line);
495 
496  // Initialise the variables to be filled
497  double T;
498  double E;
499 
500  // Filled T and E with the first two numbers of each file line
501  iss>>T>>E;
502 
503  // Fill the fTMapWithVec container with all the different T values contained within the file.
504  // Duplicate must be avoided and this is the purpose of the if statement
505  if (T != fTMapWithVec[materialName][particleName].back()) fTMapWithVec[materialName][particleName].push_back(T);
506 
507  // iterate on each shell of the corresponding material
508  for (int shell=0, eshell=ptbStructure.NumberOfLevels(materialName); shell<eshell; ++shell)
509  {
510  // map[material][particle][shell][T][E]=diffCrossSectionValue
511  // Fill the map with the informations of the input file
512  iss>>diffCrossSectionData[materialName][particleName][shell][T][E];
513 
514  if(materialName!="G4_WATER")
515  {
516  // map[material][particle][shell][T][CS]=E
517  // Fill the map
518  fEnergySecondaryData[materialName][particleName][shell][T][diffCrossSectionData[materialName][particleName][shell][T][E] ]=E;
519 
520  // map[material][particle][shell][T]=CS_vector
521  // Fill the vector within the map
522  fProbaShellMap[materialName][particleName][shell][T].push_back(diffCrossSectionData[materialName][particleName][shell][T][E]);
523  }
524  else
525  {
526  diffCrossSectionData[materialName][particleName][shell][T][E]*=scaleFactor;
527 
528  fEMapWithVector[materialName][particleName][T].push_back(E);
529  }
530  }
531  }
532 }
533 
534 
535 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
536 
538  G4double k, G4int shell, const G4String& materialName)
539 {
540  if (particleDefinition == G4Electron::ElectronDefinition())
541  {
542  //G4double Tcut=25.0E-6;
543  G4double maximumEnergyTransfer=0.;
544  if ((k+ptbStructure.IonisationEnergy(shell, materialName))/2. > k) maximumEnergyTransfer=k;
545  else maximumEnergyTransfer = (k+ptbStructure.IonisationEnergy(shell,materialName))/2.;
546 
547  // SI : original method
548  /*
549  G4double crossSectionMaximum = 0.;
550  for(G4double value=waterStructure.IonisationEnergy(shell); value<=maximumEnergyTransfer; value+=0.1*eV)
551  {
552  G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
553  if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
554  }
555  */
556 
557 
558  // SI : alternative method
559 
560  //if (k > Tcut)
561  //{
562  G4double crossSectionMaximum = 0.;
563 
564  G4double minEnergy = ptbStructure.IonisationEnergy(shell, materialName);
565  G4double maxEnergy = maximumEnergyTransfer;
566  G4int nEnergySteps = 50;
567  G4double value(minEnergy);
568  G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
569  G4int step(nEnergySteps);
570  while (step>0)
571  {
572  step--;
573  G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell, materialName);
574  if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
575  value *= stpEnergy;
576 
577  }
578  //
579 
580 
581  G4double secondaryElectronKineticEnergy=0.;
582 
583  do
584  {
585  secondaryElectronKineticEnergy = G4UniformRand() * (maximumEnergyTransfer-ptbStructure.IonisationEnergy(shell, materialName));
586 
587  } while(G4UniformRand()*crossSectionMaximum >
588  DifferentialCrossSection(particleDefinition, k/eV,(secondaryElectronKineticEnergy+ptbStructure.IonisationEnergy(shell, materialName))/eV,shell, materialName));
589 
590  return secondaryElectronKineticEnergy;
591 
592  // }
593 
594  // else if (k < Tcut)
595  // {
596 
597  // G4double bindingEnergy = ptbStructure.IonisationEnergy(shell, materialName);
598  // G4double maxEnergy = ((k-bindingEnergy)/2.);
599 
600  // G4double secondaryElectronKineticEnergy = G4UniformRand()*maxEnergy;
601  // return secondaryElectronKineticEnergy;
602  // }
603  }
604 
605 
606  else if (particleDefinition == G4Proton::ProtonDefinition())
607  {
608  G4double maximumKineticEnergyTransfer = 4.* (electron_mass_c2 / proton_mass_c2) * k;
609 
610  G4double crossSectionMaximum = 0.;
611  for (G4double value = ptbStructure.IonisationEnergy(shell, materialName);
612  value<=4.*ptbStructure.IonisationEnergy(shell, materialName) ;
613  value+=0.1*eV)
614  {
615  G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell, materialName);
616  if (differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
617  }
618 
619  G4double secondaryElectronKineticEnergy = 0.;
620  do
621  {
622  secondaryElectronKineticEnergy = G4UniformRand() * maximumKineticEnergyTransfer;
623  } while(G4UniformRand()*crossSectionMaximum >=
624  DifferentialCrossSection(particleDefinition, k/eV,(secondaryElectronKineticEnergy+ptbStructure.IonisationEnergy(shell, materialName))/eV,shell, materialName));
625 
626  return secondaryElectronKineticEnergy;
627  }
628 
629  return 0;
630 }
631 
632 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
633 
635  G4double k,
636  G4double secKinetic,
637  G4double & cosTheta,
638  G4double & phi)
639 {
640  if (particleDefinition == G4Electron::ElectronDefinition())
641  {
642  phi = twopi * G4UniformRand();
643  if (secKinetic < 50.*eV) cosTheta = (2.*G4UniformRand())-1.;
644  else if (secKinetic <= 200.*eV)
645  {
646  if (G4UniformRand() <= 0.1) cosTheta = (2.*G4UniformRand())-1.;
647  else cosTheta = G4UniformRand()*(std::sqrt(2.)/2);
648  }
649  else
650  {
651  G4double sin2O = (1.-secKinetic/k) / (1.+secKinetic/(2.*electron_mass_c2));
652  cosTheta = std::sqrt(1.-sin2O);
653  }
654  }
655 
656  else if (particleDefinition == G4Proton::ProtonDefinition())
657  {
658  G4double maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
659  phi = twopi * G4UniformRand();
660 
661  // cosTheta = std::sqrt(secKinetic / maxSecKinetic);
662 
663  // Restriction below 100 eV from Emfietzoglou (2000)
664 
665  if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
666  else cosTheta = (2.*G4UniformRand())-1.;
667 
668  }
669 }
670 
671 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
672 
674  G4double k,
675  G4double energyTransfer,
676  G4int ionizationLevelIndex,
677  const G4String& materialName)
678 {
679  G4double sigma = 0.;
680  const G4String& particleName = particleDefinition->GetParticleName();
681 
682  G4double shellEnergy (ptbStructure.IonisationEnergy(ionizationLevelIndex, materialName));
683  G4double kSE (energyTransfer-shellEnergy);
684 
685  if (energyTransfer >= shellEnergy)
686  {
687  G4double valueT1 = 0;
688  G4double valueT2 = 0;
689  G4double valueE21 = 0;
690  G4double valueE22 = 0;
691  G4double valueE12 = 0;
692  G4double valueE11 = 0;
693 
694  G4double xs11 = 0;
695  G4double xs12 = 0;
696  G4double xs21 = 0;
697  G4double xs22 = 0;
698 
699  if (particleDefinition == G4Electron::ElectronDefinition())
700  {
701  // k should be in eV and energy transfer eV also
702  std::vector<double>::iterator t2 = std::upper_bound(fTMapWithVec[materialName][particleName].begin(),fTMapWithVec[materialName][particleName].end(), k);
703  std::vector<double>::iterator t1 = t2-1;
704 
705  // SI : the following condition avoids situations where energyTransfer >last vector element
706  if (kSE <= fEMapWithVector[materialName][particleName][(*t1)].back() && kSE <= fEMapWithVector[materialName][particleName][(*t2)].back() )
707  {
708  std::vector<double>::iterator e12 = std::upper_bound(fEMapWithVector[materialName][particleName][(*t1)].begin(),fEMapWithVector[materialName][particleName][(*t1)].end(), kSE);
709  std::vector<double>::iterator e11 = e12-1;
710 
711  std::vector<double>::iterator e22 = std::upper_bound(fEMapWithVector[materialName][particleName][(*t2)].begin(),fEMapWithVector[materialName][particleName][(*t2)].end(), kSE);
712  std::vector<double>::iterator e21 = e22-1;
713 
714  valueT1 =*t1;
715  valueT2 =*t2;
716  valueE21 =*e21;
717  valueE22 =*e22;
718  valueE12 =*e12;
719  valueE11 =*e11;
720 
721  xs11 = diffCrossSectionData[materialName][particleName][ionizationLevelIndex][valueT1][valueE11];
722  xs12 = diffCrossSectionData[materialName][particleName][ionizationLevelIndex][valueT1][valueE12];
723  xs21 = diffCrossSectionData[materialName][particleName][ionizationLevelIndex][valueT2][valueE21];
724  xs22 = diffCrossSectionData[materialName][particleName][ionizationLevelIndex][valueT2][valueE22];
725  }
726  }
727 
728  if (particleDefinition == G4Proton::ProtonDefinition())
729  {
730  // k should be in eV and energy transfer eV also
731  std::vector<double>::iterator t2 = std::upper_bound(fTMapWithVec[materialName][particleName].begin(),fTMapWithVec[materialName][particleName].end(), k);
732  std::vector<double>::iterator t1 = t2-1;
733 
734  std::vector<double>::iterator e12 = std::upper_bound(fEMapWithVector[materialName][particleName][(*t1)].begin(),fEMapWithVector[materialName][particleName][(*t1)].end(), kSE);
735  std::vector<double>::iterator e11 = e12-1;
736 
737  std::vector<double>::iterator e22 = std::upper_bound(fEMapWithVector[materialName][particleName][(*t2)].begin(),fEMapWithVector[materialName][particleName][(*t2)].end(), kSE);
738  std::vector<double>::iterator e21 = e22-1;
739 
740  valueT1 =*t1;
741  valueT2 =*t2;
742  valueE21 =*e21;
743  valueE22 =*e22;
744  valueE12 =*e12;
745  valueE11 =*e11;
746 
747  xs11 = diffCrossSectionData[materialName][particleName][ionizationLevelIndex][valueT1][valueE11];
748  xs12 = diffCrossSectionData[materialName][particleName][ionizationLevelIndex][valueT1][valueE12];
749  xs21 = diffCrossSectionData[materialName][particleName][ionizationLevelIndex][valueT2][valueE21];
750  xs22 = diffCrossSectionData[materialName][particleName][ionizationLevelIndex][valueT2][valueE22];
751  }
752 
753  G4double xsProduct = xs11 * xs12 * xs21 * xs22;
754 
755  if (xsProduct != 0.)
756  {
757  sigma = QuadInterpolator(valueE11, valueE12,
758  valueE21, valueE22,
759  xs11, xs12,
760  xs21, xs22,
761  valueT1, valueT2,
762  k, kSE);
763  }
764  }
765 
766 
767  return sigma;
768 }
769 
770 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
771 
773 {
774  // k should be in eV
775 
776  // Schematic explanation.
777  // We will do an interpolation to get a final E value (ejected electron energy).
778  // 1/ We choose a random number between 0 and 1 (ie we select a cumulated cross section).
779  // 2/ We look for T_lower and T_upper.
780  // 3/ We look for the cumulated corresponding cross sections and their associated E values.
781  //
782  // T_low | CS_low_1 -> E_low_1
783  // | CS_low_2 -> E_low_2
784  // T_up | CS_up_1 -> E_up_1
785  // | CS_up_2 -> E_up_2
786  //
787  // 4/ We interpolate to get our E value.
788  //
789  // T_low | CS_low_1 -> E_low_1 -----
790  // | |----> E_low --
791  // | CS_low_2 -> E_low_2 ----- |
792  // | ---> E_final
793  // T_up | CS_up_1 -> E_up_1 ------- |
794  // | |----> E_up ---
795  // | CS_up_2 -> E_up_2 -------
796 
797  // Initialize some values
798  //
799  G4double ejectedElectronEnergy = 0.;
800  G4double valueK1 = 0;
801  G4double valueK2 = 0;
802  G4double valueCumulCS21 = 0;
803  G4double valueCumulCS22 = 0;
804  G4double valueCumulCS12 = 0;
805  G4double valueCumulCS11 = 0;
806  G4double secElecE11 = 0;
807  G4double secElecE12 = 0;
808  G4double secElecE21 = 0;
809  G4double secElecE22 = 0;
810  G4String particleName = particleDefinition->GetParticleName();
811 
812  // ***************************************************************************
813  // Get a random number between 0 and 1 to compare with the cumulated CS
814  // ***************************************************************************
815  //
816  // It will allow us to choose an ejected electron energy with respect to the CS.
817  G4double random = G4UniformRand();
818 
819  // **********************************************
820  // Take the input from the data tables
821  // **********************************************
822 
823  // Cumulated tables are like this: T E cumulatedCS1 cumulatedCS2 cumulatedCS3
824  // We have two sets of loaded data: fTMapWithVec which contains data about T (incident particle energy)
825  // and fProbaShellMap which contains cumulated cross section data.
826  // Since we already have a specific T energy value which could not be explicitly in the table, we must interpolate all the values.
827 
828  // First, we select the upper and lower T data values surrounding our T value (ie "k").
829  std::vector<double>::iterator k2 = std::upper_bound(fTMapWithVec[materialName][particleName].begin(),fTMapWithVec[materialName][particleName].end(), k);
830  std::vector<double>::iterator k1 = k2-1;
831 
832  // Check if we have found a k2 value (0 if we did not found it).
833  // A missing k2 value can be caused by a energy to high for the data table,
834  // Ex : table done for 12*eV -> 1000*eV and k=2000*eV
835  // then k2 = 0 and k1 = max of the table.
836  // To detect this, we check that k1 is not superior to k2.
837  if(*k1 > *k2)
838  {
839  // Error
840  G4cerr<<"**************** Fatal error ******************"<<G4endl;
841  G4cerr<<"G4DNAPTBIonisationModel::RandomizeEjectedElectronEnergyFromCumulated"<<G4endl;
842  G4cerr<<"You have *k1 > *k2 with k1 "<<*k1<<" and k2 "<<*k2<<G4endl;
843  G4cerr<<"This may be because the energy of the incident particle is to high for the data table."<<G4endl;
844  G4cerr<<"Particle energy (eV): "<<k<<G4endl;
845  exit(EXIT_FAILURE);
846  }
847 
848 
849  // We have a random number and we select the cumulated cross section data values surrounding our random number.
850  // But we need to do that for each T value (ie two T values) previously selected.
851  //
852  // First one.
853  std::vector<double>::iterator cumulCS12 = std::upper_bound(fProbaShellMap[materialName][particleName][ionizationLevelIndex][(*k1)].begin(),
854  fProbaShellMap[materialName][particleName][ionizationLevelIndex][(*k1)].end(), random);
855  std::vector<double>::iterator cumulCS11 = cumulCS12-1;
856  // Second one.
857  std::vector<double>::iterator cumulCS22 = std::upper_bound(fProbaShellMap[materialName][particleName][ionizationLevelIndex][(*k2)].begin(),
858  fProbaShellMap[materialName][particleName][ionizationLevelIndex][(*k2)].end(), random);
859  std::vector<double>::iterator cumulCS21 = cumulCS22-1;
860 
861  // Now that we have the "values" through pointers, we access them.
862  valueK1 = *k1;
863  valueK2 = *k2;
864  valueCumulCS11 = *cumulCS11;
865  valueCumulCS12 = *cumulCS12;
866  valueCumulCS21 = *cumulCS21;
867  valueCumulCS22 = *cumulCS22;
868 
869  // *************************************************************
870  // Do the interpolation to get the ejected electron energy
871  // *************************************************************
872 
873  // Here we will get four E values corresponding to our four cumulated cross section values previously selected.
874  // But we need to take into account a specific case: we have selected a shell by using the ionisation cross section table
875  // and, since we get two T values, we could have differential cross sections (or cumulated) equal to 0 for the lower T
876  // and not for the upper T. When looking for the cumulated cross section values which surround the selected random number (for the lower T),
877  // the upper_bound method will only found 0 values. Thus, the upper_bound method will return the last E value present in the table for the
878  // selected T. The last E value being the highest, we will later perform an interpolation between a high E value (for the lower T) and
879  // a small E value (for the upper T). This is inconsistent because if the cross section are equal to zero for the lower T then it
880  // means it is not possible to ionize and, thus, to have a secondary electron. But, in our situation, it is possible to ionize for the upper T
881  // AND for an interpolate T value between Tupper Tlower. That's why the final E value should be interpolate between 0 and the E value (upper T).
882  //
883  if(cumulCS12==fProbaShellMap[materialName][particleName][ionizationLevelIndex][(*k1)].end())
884  {
885  // Here we are in the special case and we force Elower1 and Elower2 to be equal at 0 for the interpolation.
886  secElecE11 = 0;
887  secElecE12 = 0;
888  secElecE21 = fEnergySecondaryData[materialName][particleName][ionizationLevelIndex][valueK2][valueCumulCS21];
889  secElecE22 = fEnergySecondaryData[materialName][particleName][ionizationLevelIndex][valueK2][valueCumulCS22];
890 
891  valueCumulCS11 = 0;
892  valueCumulCS12 = 0;
893  }
894  else
895  {
896  // No special case, interpolation will happen as usual.
897  secElecE11 = fEnergySecondaryData[materialName][particleName][ionizationLevelIndex][valueK1][valueCumulCS11];
898  secElecE12 = fEnergySecondaryData[materialName][particleName][ionizationLevelIndex][valueK1][valueCumulCS12];
899  secElecE21 = fEnergySecondaryData[materialName][particleName][ionizationLevelIndex][valueK2][valueCumulCS21];
900  secElecE22 = fEnergySecondaryData[materialName][particleName][ionizationLevelIndex][valueK2][valueCumulCS22];
901  }
902 
903  ejectedElectronEnergy = QuadInterpolator(valueCumulCS11, valueCumulCS12,
904  valueCumulCS21, valueCumulCS22,
905  secElecE11, secElecE12,
906  secElecE21, secElecE22,
907  valueK1, valueK2,
908  k, random);
909 
910  // **********************************************
911  // Some tests for debugging
912  // **********************************************
913 
914  G4double bindingEnergy (ptbStructure.IonisationEnergy(ionizationLevelIndex, materialName)/eV);
915  if(k-ejectedElectronEnergy-bindingEnergy<=0 || ejectedElectronEnergy<=0)
916  {
917  G4cout<<"k "<<k<<G4endl;
918  G4cout<<"material "<<materialName<<G4endl;
919  G4cout<<"secondaryKin "<<ejectedElectronEnergy<<G4endl;
920  G4cout<<"shell "<<ionizationLevelIndex<<G4endl;
921  G4cout<<"bindingEnergy "<<bindingEnergy<<G4endl;
922  G4cout<<"scatteredEnergy "<<k-ejectedElectronEnergy-bindingEnergy<<G4endl;
923  G4cout<<"rand "<<random<<G4endl;
924  G4cout<<"surrounding k values: valueK1 valueK2\n"<<valueK1<<" "<<valueK2<<G4endl;
925  G4cout<<"surrounding E values: secElecE11 secElecE12 secElecE21 secElecE22\n"
926  <<secElecE11<<" "<<secElecE12<<" "<<secElecE21<<" "<<secElecE22<<" "<<G4endl;
927  G4cout<<"surrounding cumulCS values: valueCumulCS11 valueCumulCS12 valueCumulCS21 valueCumulCS22\n"
928  <<valueCumulCS11<<" "<<valueCumulCS12<<" "<<valueCumulCS21<<" "<<valueCumulCS22<<" "<<G4endl;
929  G4cerr<<"*****************************"<<G4endl;
930  G4cerr<<"Fatal error, EXIT."<<G4endl;
931  exit(EXIT_FAILURE);
932  }
933 
934  return ejectedElectronEnergy*eV;
935 }
936 
937 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
938 
940  G4double e2,
941  G4double e,
942  G4double xs1,
943  G4double xs2)
944 {
945  G4double value (0);
946 
947  // Switch to log-lin interpolation for faster code
948 
949  if ((e2-e1)!=0 && xs1 !=0 && xs2 !=0)
950  {
951  G4double d1 = std::log10(xs1);
952  G4double d2 = std::log10(xs2);
953  value = std::pow(10.,(d1 + (d2 - d1)*(e - e1)/ (e2 - e1)) );
954  }
955 
956  // Switch to lin-lin interpolation for faster code
957  // in case one of xs1 or xs2 (=cum proba) value is zero
958 
959  if ((e2-e1)!=0 && (xs1 ==0 || xs2 ==0))
960  {
961  G4double d1 = xs1;
962  G4double d2 = xs2;
963  value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
964  }
965 
966  return value;
967 }
968 
969 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
970 
972  G4double e21, G4double e22,
973  G4double xs11, G4double xs12,
974  G4double xs21, G4double xs22,
976  G4double t, G4double e)
977 {
978  G4double interpolatedvalue1 (-1);
979  if(xs11!=xs12) interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
980  else interpolatedvalue1 = xs11;
981 
982  G4double interpolatedvalue2 (-1);
983  if(xs21!=xs22) interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
984  else interpolatedvalue2 = xs21;
985 
986  G4double value (-1);
987  if(interpolatedvalue1!=interpolatedvalue2) value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
988  else value = interpolatedvalue1;
989 
990  return value;
991 
992  // G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
993  // G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
994  // G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
995  // return value;
996 }