ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNACPA100IonisationModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DNACPA100IonisationModel.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 // CPA100 ionisation model class for electrons
27 //
28 // Based on the work of M. Terrissol and M. C. Bordage
29 //
30 // Users are requested to cite the following papers:
31 // - M. Terrissol, A. Baudre, Radiat. Prot. Dosim. 31 (1990) 175-177
32 // - M.C. Bordage, J. Bordes, S. Edel, M. Terrissol, X. Franceries,
33 // M. Bardies, N. Lampe, S. Incerti, Phys. Med. 32 (2016) 1833-1840
34 //
35 // Authors of this class:
36 // M.C. Bordage, M. Terrissol, S. Edel, J. Bordes, S. Incerti
37 //
38 // 15.01.2014: creation
39 //
40 
41 #include "G4DNACPA100IonisationModel.hh"
42 #include "G4PhysicalConstants.hh"
43 #include "G4SystemOfUnits.hh"
44 #include "G4UAtomicDeexcitation.hh"
45 #include "G4LossTableManager.hh"
46 #include "G4DNAChemistryManager.hh"
47 #include "G4DNAMolecularMaterial.hh"
48 
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
50 
51 using namespace std;
52 
53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
54 
55 G4DNACPA100IonisationModel::G4DNACPA100IonisationModel(const G4ParticleDefinition*,
56  const G4String& nam)
57 :G4VEmModel(nam),isInitialised(false)
58 {
59  verboseLevel= 0;
60  // Verbosity scale:
61  // 0 = nothing
62  // 1 = warning for energy non-conservation
63  // 2 = details of energy budget
64  // 3 = calculation of cross sections, file openings, sampling of atoms
65  // 4 = entering in methods
66 
67  if( verboseLevel>0 )
68  {
69  G4cout << "CPA100 ionisation model is constructed " << G4endl;
70  }
71 
72  SetLowEnergyLimit(11*eV);
73  SetHighEnergyLimit(255955*eV);
74 
75  // Mark this model as "applicable" for atomic deexcitation
76  SetDeexcitationFlag(true);
77  fAtomDeexcitation = 0;
78  fParticleChangeForGamma = 0;
79  fpMolWaterDensity = 0;
80 
81  // Selection of computation method
82 
83  // useDcs = true if usage of dcs for sampling of secondaries
84  // useDcs = false if usage of composition sampling (DEFAULT)
85 
86  useDcs = true;
87 
88  // if useDcs is true, one has the following choice
89  // fasterCode = true for usage of cumulated dcs (DEFAULT)
90  // fasterCode = false for usage of non-cumulated dcs
91 
92  fasterCode = true;
93 
94  // Selection of stationary mode
95 
96  statCode = false;
97 }
98 
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100 
101 G4DNACPA100IonisationModel::~G4DNACPA100IonisationModel()
102 {
103  // Cross section
104 
105  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
106  for (pos = tableData.begin(); pos != tableData.end(); ++pos)
107  {
108  G4DNACrossSectionDataSet* table = pos->second;
109  delete table;
110  }
111 
112  // Final state
113 
114  eVecm.clear();
115 
116 }
117 
118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
119 
120 void G4DNACPA100IonisationModel::Initialise(const G4ParticleDefinition* particle,
121  const G4DataVector& /*cuts*/)
122 {
123 
124  if (verboseLevel > 3)
125  G4cout << "Calling G4DNACPA100IonisationModel::Initialise()" << G4endl;
126 
127  // Energy limits
128 
129  // The following file is proved by M. Terrissol et al. (sigion3)
130 
131  G4String fileElectron("dna/sigma_ionisation_e_cpa100_form_rel");
132 
133  G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
134 
135  G4String electron;
136 
137  G4double scaleFactor = 1.e-20 * m*m;
138 
139  char *path = getenv("G4LEDATA");
140 
141  // *** ELECTRON
142 
143  electron = electronDef->GetParticleName();
144 
145  tableFile[electron] = fileElectron;
146 
147  // Cross section
148 
149  G4DNACrossSectionDataSet* tableE =
150  new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
151 
152  //G4DNACrossSectionDataSet* tableE =
153  // new G4DNACrossSectionDataSet(new G4DNACPA100LogLogInterpolation, eV,scaleFactor );
154 
155  tableE->LoadData(fileElectron);
156 
157  tableData[electron] = tableE;
158 
159  // Final state
160 
161  // ******************************
162 
163  if (useDcs)
164  {
165 
166  std::ostringstream eFullFileName;
167 
168  if (fasterCode) eFullFileName << path << "/dna/sigmadiff_cumulated_ionisation_e_cpa100_rel.dat";
169 
170  if (!fasterCode) eFullFileName << path << "/dna/sigmadiff_ionisation_e_cpa100_rel.dat";
171 
172  std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
173 
174  if (!eDiffCrossSection)
175  {
176  if (fasterCode) G4Exception("G4DNACPA100IonisationModel::Initialise","em0003",
177  FatalException,"Missing data file:/dna/sigmadiff_cumulated_ionisation_e_cpa100_rel.dat");
178 
179  if (!fasterCode) G4Exception("G4DNACPA100IonisationModel::Initialise","em0003",
180  FatalException,"Missing data file:/dna/sigmadiff_ionisation_e_cpa100_rel.dat");
181  }
182 
183  // Clear the arrays for re-initialization case (MT mode)
184  // March 25th, 2014 - Vaclav Stepan, Sebastien Incerti
185 
186  eTdummyVec.clear();
187  eVecm.clear();
188  eProbaShellMap->clear();
189  eDiffCrossSectionData->clear();
190  eNrjTransfData->clear();
191 
192  //
193 
194  eTdummyVec.push_back(0.);
195  while(!eDiffCrossSection.eof())
196  {
197  G4double tDummy;
198  G4double eDummy;
199  eDiffCrossSection>>tDummy>>eDummy;
200  if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
201  for (G4int j=0; j<5; j++)
202  {
203  eDiffCrossSection>>eDiffCrossSectionData[j][tDummy][eDummy];
204 
205  if (fasterCode)
206  {
207  eNrjTransfData[j][tDummy][eDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
208  eProbaShellMap[j][tDummy].push_back(eDiffCrossSectionData[j][tDummy][eDummy]);
209  }
210 
211  // SI - only if eof is not reached
212  if (!eDiffCrossSection.eof() && !fasterCode)
213  eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
214 
215  if (!fasterCode) eVecm[tDummy].push_back(eDummy);
216 
217  }
218  }
219 
220  //
221 
222  } // end of if (useDcs)
223 
224  // ******************************
225 
226  //
227 
228  if( verboseLevel>0 )
229  {
230  G4cout << "CPA100 ionisation model is initialized " << G4endl
231  << "Energy range: "
232  << LowEnergyLimit() / eV << " eV - "
233  << HighEnergyLimit() / keV << " keV for "
234  << particle->GetParticleName()
235  << G4endl;
236  }
237 
238  // Initialize water density pointer
239  fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
240 
241  // AD
242  fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
243 
244  if (isInitialised) return;
245  fParticleChangeForGamma = GetParticleChangeForGamma();
246  isInitialised = true;
247 }
248 
249 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
250 
251 G4double G4DNACPA100IonisationModel::CrossSectionPerVolume( const G4Material* material,
252  const G4ParticleDefinition* particleDefinition,
253  G4double ekin,
254  G4double,
255  G4double)
256 {
257 
258  if (verboseLevel > 3)
259  G4cout << "Calling CrossSectionPerVolume() of G4DNACPA100IonisationModel" << G4endl;
260 
261  if (particleDefinition != G4Electron::ElectronDefinition()) return 0;
262 
263  // Calculate total cross section for model
264 
265  G4double sigma=0;
266 
267  G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
268 
269  const G4String& particleName = particleDefinition->GetParticleName();
270 
271  if (ekin >= LowEnergyLimit() && ekin <= HighEnergyLimit())
272  {
273  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
274  pos = tableData.find(particleName);
275 
276  if (pos != tableData.end())
277  {
278  G4DNACrossSectionDataSet* table = pos->second;
279  if (table != 0) sigma = table->FindValue(ekin);
280  }
281  else
282  {
283  G4Exception("G4DNACPA100IonisationModel::CrossSectionPerVolume","em0002",
284  FatalException,"Model not applicable to particle type.");
285  }
286  }
287 
288  if (verboseLevel > 2)
289  {
290  G4cout << "__________________________________" << G4endl;
291  G4cout << "G4DNACPA100IonisationModel - XS INFO START" << G4endl;
292  G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
293  G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
294  G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
295  G4cout << "G4DNACPA100IonisationModel - XS INFO END" << G4endl;
296  }
297 
298  return sigma*waterDensity;
299 }
300 
301 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
302 
303 void G4DNACPA100IonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
304  const G4MaterialCutsCouple* ,//must be set!
305  const G4DynamicParticle* particle,
306  G4double,
307  G4double)
308 {
309  if (verboseLevel > 3)
310  G4cout << "Calling SampleSecondaries() of G4DNACPA100IonisationModel" << G4endl;
311 
312  G4double k = particle->GetKineticEnergy();
313 
314  const G4String& particleName = particle->GetDefinition()->GetParticleName();
315 
316  if (k >= LowEnergyLimit() && k <= HighEnergyLimit())
317  {
318  G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
319  G4double particleMass = particle->GetDefinition()->GetPDGMass();
320  G4double totalEnergy = k + particleMass;
321  G4double pSquare = k * (totalEnergy + particleMass);
322  G4double totalMomentum = std::sqrt(pSquare);
323 
324  G4int ionizationShell = -1;
325 
326  ionizationShell = RandomSelect(k,particleName);
327 
328  //SI: PROTECTION FOR G4LOGLOGINTERPOLATION ON UPPER VALUE
329  if (k<waterStructure.IonisationEnergy(ionizationShell)) { return; }
330 
331  G4double bindingEnergy = 0;
332  bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
333 
334  G4double secondaryKinetic=-1000*eV;
335 
336  if (useDcs && !fasterCode)
337  secondaryKinetic = RandomizeEjectedElectronEnergy(particle->GetDefinition(),k,ionizationShell);
338 
339  if (useDcs && fasterCode)
340  secondaryKinetic = RandomizeEjectedElectronEnergyFromCumulatedDcs(particle->GetDefinition(),k,ionizationShell);
341 
342  if (!useDcs)
343  secondaryKinetic = RandomizeEjectedElectronEnergyFromCompositionSampling(particle->GetDefinition(),k,ionizationShell);
344 
345  // Quick test
346  /*
347  FILE* myFile;
348  myFile=fopen("nrj.txt","a");
349  fprintf(myFile,"%e\n", secondaryKinetic/eV );
350  fclose(myFile);
351  */
352 
353  G4double cosTheta = 0.;
354  G4double phi = 0.;
355  RandomizeEjectedElectronDirection(particle->GetDefinition(), k,secondaryKinetic, cosTheta, phi);
356 
357  G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
358  G4double dirX = sinTheta*std::cos(phi);
359  G4double dirY = sinTheta*std::sin(phi);
360  G4double dirZ = cosTheta;
361  G4ThreeVector deltaDirection(dirX,dirY,dirZ);
362  deltaDirection.rotateUz(primaryDirection);
363 
364  // SI - For atom. deexc. tagging - 23/05/2017
365  if (secondaryKinetic>0)
366  {
367  G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
368  fvect->push_back(dp);
369  }
370  //
371 
372  if (particle->GetDefinition() == G4Electron::ElectronDefinition())
373  {
374  G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
375 
376  G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
377  G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
378  G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
379  G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
380  finalPx /= finalMomentum;
381  finalPy /= finalMomentum;
382  finalPz /= finalMomentum;
383 
384  G4ThreeVector direction;
385  direction.set(finalPx,finalPy,finalPz);
386 
387  fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
388  }
389 
390  else fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection) ;
391 
392  // SI - For atom. deexc. tagging - 23/05/2017
393 
394  // AM: sample deexcitation
395  // here we assume that H_{2}O electronic levels are the same of Oxigen.
396  // this can be considered true with a rough 10% error in energy on K-shell,
397 
398  size_t secNumberInit = 0;// need to know at a certain point the energy of secondaries
399  size_t secNumberFinal = 0;// So I'll make the diference and then sum the energies
400 
401  G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
402 
403  // SI: only atomic deexcitation from K shell is considered
404  if(fAtomDeexcitation && ionizationShell == 4)
405  {
406  G4int Z = 8;
407  const G4AtomicShell* shell =
408  fAtomDeexcitation->GetAtomicShell(Z, G4AtomicShellEnumerator(0));
409  secNumberInit = fvect->size();
410  fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
411  secNumberFinal = fvect->size();
412 
413  if(secNumberFinal > secNumberInit)
414  {
415  for (size_t i=secNumberInit; i<secNumberFinal; ++i)
416  {
417  //Check if there is enough residual energy
418  if (bindingEnergy >= ((*fvect)[i])->GetKineticEnergy())
419  {
420  //Ok, this is a valid secondary: keep it
421  bindingEnergy -= ((*fvect)[i])->GetKineticEnergy();
422  }
423  else
424  {
425  //Invalid secondary: not enough energy to create it!
426  //Keep its energy in the local deposit
427  delete (*fvect)[i];
428  (*fvect)[i]=0;
429  }
430  }
431  }
432 
433  }
434 
435  //This should never happen
436  if(bindingEnergy < 0.0)
437  G4Exception("G4DNACPA100IonisatioModel1::SampleSecondaries()",
438  "em2050",FatalException,"Negative local energy deposit");
439 
440  //bindingEnergy has been decreased
441  //by the amount of energy taken away by deexc. products
442 
443  if (!statCode)
444  {
445  fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
446  fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy);
447  }
448  else
449  {
450  fParticleChangeForGamma->SetProposedKineticEnergy(k);
451  fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy);
452  }
453 
454  // TEST //////////////////////////
455  // if (secondaryKinetic<0) abort();
456  // if (scatteredEnergy<0) abort();
457  // if (k-scatteredEnergy-secondaryKinetic-deexSecEnergy<0) abort();
458  // if (k-scatteredEnergy<0) abort();
460 
461  const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
462  G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
463  ionizationShell,
464  theIncomingTrack);
465  }
466 
467 }
468 
469 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
470 
471 G4double G4DNACPA100IonisationModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
472  G4double k, G4int shell)
473 {
474  // G4cout << "*** SLOW computation for " << " " << particleDefinition->GetParticleName() << G4endl;
475 
476  if (particleDefinition == G4Electron::ElectronDefinition())
477  {
478  G4double maximumEnergyTransfer=0.;
479  if ((k+waterStructure.IonisationEnergy(shell))/2. > k) maximumEnergyTransfer=k;
480  else maximumEnergyTransfer = (k+waterStructure.IonisationEnergy(shell))/2.;
481 
482  // SI : original method
483  /*
484  G4double crossSectionMaximum = 0.;
485  for(G4double value=waterStructure.IonisationEnergy(shell); value<=maximumEnergyTransfer; value+=0.1*eV)
486  {
487  G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
488  if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
489  }
490  */
491 
492  // SI : alternative method
493 
494  G4double crossSectionMaximum = 0.;
495 
496  G4double minEnergy = waterStructure.IonisationEnergy(shell);
497  G4double maxEnergy = maximumEnergyTransfer;
498 
499  // nEnergySteps can be optimized - 100 by default
500  G4int nEnergySteps = 50;
501 
502  // *** METHOD 1
503  // FOR SLOW COMPUTATION ONLY
504  /*
505  G4double value(minEnergy);
506  G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
507  G4int step(nEnergySteps);
508  while (step>0)
509  {
510  step--;
511  G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
512  if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
513  value*=stpEnergy;
514  }
515  */
516 
517  // *** METHOD 2 : Faster method for CPA100 only since DCS is monotonously decreasing
518  // FOR SLOW COMPUTATION ONLY
519 
520  G4double value(minEnergy);
521  G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
522  G4int step(nEnergySteps);
523  G4double differentialCrossSection = 0.;
524  while (step>0)
525  {
526  step--;
527  differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
528  if(differentialCrossSection >0)
529  {
530  crossSectionMaximum=differentialCrossSection;
531  break;
532  }
533  value*=stpEnergy;
534  }
535 
536  //
537 
538  G4double secondaryElectronKineticEnergy=0.;
539  do
540  {
541  secondaryElectronKineticEnergy = G4UniformRand() * (maximumEnergyTransfer-waterStructure.IonisationEnergy(shell));
542  } while(G4UniformRand()*crossSectionMaximum >
543  DifferentialCrossSection(particleDefinition, k/eV,
544  (secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
545 
546  return secondaryElectronKineticEnergy;
547 
548  }
549 
550  return 0;
551 }
552 
553 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
554 
555 void G4DNACPA100IonisationModel::RandomizeEjectedElectronDirection(G4ParticleDefinition*,
556  G4double k,
557  G4double secKinetic,
558  G4double & cosTheta,
559  G4double & phi )
560 {
561 
562  phi = twopi * G4UniformRand();
563  G4double sin2O = (1.-secKinetic/k) / (1.+secKinetic/(2.*electron_mass_c2));
564  cosTheta = std::sqrt(1.-sin2O);
565 
566 }
567 
568 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
569 
570 G4double G4DNACPA100IonisationModel::DifferentialCrossSection(G4ParticleDefinition * particleDefinition,
571  G4double k,
572  G4double energyTransfer,
573  G4int ionizationLevelIndex)
574 {
575  G4double sigma = 0.;
576 
577  if (energyTransfer >= waterStructure.IonisationEnergy(ionizationLevelIndex)/eV)
578  {
579  G4double valueT1 = 0;
580  G4double valueT2 = 0;
581  G4double valueE21 = 0;
582  G4double valueE22 = 0;
583  G4double valueE12 = 0;
584  G4double valueE11 = 0;
585 
586  G4double xs11 = 0;
587  G4double xs12 = 0;
588  G4double xs21 = 0;
589  G4double xs22 = 0;
590 
591  if (particleDefinition == G4Electron::ElectronDefinition())
592  {
593  // Protection against out of boundary access
594  if (k==eTdummyVec.back()) k=k*(1.-1e-12);
595  //
596 
597  // k should be in eV and energy transfer eV also
598 
599  std::vector<G4double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
600 
601  std::vector<G4double>::iterator t1 = t2-1;
602 
603  // SI : the following condition avoids situations where energyTransfer >last vector element
604 
605  if (energyTransfer <= eVecm[(*t1)].back() && energyTransfer <= eVecm[(*t2)].back() )
606  {
607  std::vector<G4double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), energyTransfer);
608  std::vector<G4double>::iterator e11 = e12-1;
609 
610  std::vector<G4double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), energyTransfer);
611  std::vector<G4double>::iterator e21 = e22-1;
612 
613  valueT1 =*t1;
614  valueT2 =*t2;
615  valueE21 =*e21;
616  valueE22 =*e22;
617  valueE12 =*e12;
618  valueE11 =*e11;
619 
620  xs11 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
621  xs12 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
622  xs21 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
623  xs22 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
624 
625  }
626 
627  }
628 
629  G4double xsProduct = xs11 * xs12 * xs21 * xs22;
630  if (xsProduct != 0.)
631  {
632  sigma = QuadInterpolator( valueE11, valueE12,
633  valueE21, valueE22,
634  xs11, xs12,
635  xs21, xs22,
636  valueT1, valueT2,
637  k, energyTransfer);
638  }
639 
640  }
641  return sigma;
642 }
643 
644 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
645 
646 G4double G4DNACPA100IonisationModel::Interpolate( G4double e1,
647  G4double e2,
648  G4double e,
649  G4double xs1,
650  G4double xs2)
651 {
652 
653  G4double value = 0.;
654 
655  // Log-log interpolation by default
656 
657  if (e1!=0 && e2!=0 && (std::log10(e2)-std::log10(e1)) !=0 && !fasterCode && useDcs)
658  {
659  G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
660  G4double b = std::log10(xs2) - a*std::log10(e2);
661  G4double sigma = a*std::log10(e) + b;
662  value = (std::pow(10.,sigma));
663  }
664 
665  // Switch to lin-lin interpolation
666  /*
667  if ((e2-e1)!=0)
668  {
669  G4double d1 = xs1;
670  G4double d2 = xs2;
671  value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
672  }
673  */
674 
675  // Switch to log-lin interpolation for faster code
676 
677  if ((e2-e1)!=0 && xs1 !=0 && xs2 !=0 && fasterCode && useDcs )
678  {
679  G4double d1 = std::log10(xs1);
680  G4double d2 = std::log10(xs2);
681  value = std::pow(10.,(d1 + (d2 - d1)*(e - e1)/ (e2 - e1)) );
682  }
683 
684  // Switch to lin-lin interpolation for faster code
685  // in case one of xs1 or xs2 (=cum proba) value is zero
686 
687  if ((e2-e1)!=0 && (xs1 ==0 || xs2 ==0) && fasterCode && useDcs )
688  {
689  G4double d1 = xs1;
690  G4double d2 = xs2;
691  value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
692  }
693 
694  /*
695  G4cout
696  << e1 << " "
697  << e2 << " "
698  << e << " "
699  << xs1 << " "
700  << xs2 << " "
701  << value
702  << G4endl;
703  */
704 
705  return value;
706 }
707 
708 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
709 
710 G4double G4DNACPA100IonisationModel::QuadInterpolator(G4double e11, G4double e12,
711  G4double e21, G4double e22,
712  G4double xs11, G4double xs12,
713  G4double xs21, G4double xs22,
714  G4double t1, G4double t2,
715  G4double t, G4double e)
716 {
717  G4double interpolatedvalue1 = Interpolate(e11, e12, e, xs11, xs12);
718  G4double interpolatedvalue2 = Interpolate(e21, e22, e, xs21, xs22);
719  G4double value = Interpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
720 
721  return value;
722 }
723 
724 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
725 
726 G4int G4DNACPA100IonisationModel::RandomSelect(G4double k, const G4String& particle )
727 {
728  G4int level = 0;
729 
730  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
731  pos = tableData.find(particle);
732 
733  if (pos != tableData.end())
734  {
735  G4DNACrossSectionDataSet* table = pos->second;
736 
737  if (table != 0)
738  {
739  G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
740  const size_t n(table->NumberOfComponents());
741  size_t i(n);
742  G4double value = 0.;
743 
744  //Verification
745  /*
746  G4double tmp=200*keV;
747  G4cout << table->GetComponent(0)->FindValue(tmp)/(1e-20*m*m) << G4endl;
748  G4cout << table->GetComponent(1)->FindValue(tmp)/(1e-20*m*m) << G4endl;
749  G4cout << table->GetComponent(2)->FindValue(tmp)/(1e-20*m*m) << G4endl;
750  G4cout << table->GetComponent(3)->FindValue(tmp)/(1e-20*m*m) << G4endl;
751  G4cout << table->GetComponent(4)->FindValue(tmp)/(1e-20*m*m) << G4endl;
752  G4cout <<
753  table->GetComponent(0)->FindValue(tmp)/(1e-20*m*m) +
754  table->GetComponent(1)->FindValue(tmp)/(1e-20*m*m) +
755  table->GetComponent(2)->FindValue(tmp)/(1e-20*m*m) +
756  table->GetComponent(3)->FindValue(tmp)/(1e-20*m*m)
757  << G4endl;
758  abort();
759  */
760  //
761  //Dump
762  //
763  /*
764  G4double minEnergy = 10.985 * eV;
765  G4double maxEnergy = 255955. * eV;
766  G4int nEnergySteps = 1000;
767  G4double energy(minEnergy);
768  G4double stpEnergy(std::pow(maxEnergy/energy, 1./static_cast<G4double>(nEnergySteps-1)));
769  G4int step(nEnergySteps);
770  system ("rm -rf ionisation-cpa100.out");
771  FILE* myFile=fopen("ionisation-cpa100.out","a");
772  while (step>0)
773  {
774  step--;
775  fprintf (myFile,"%16.9le %16.9le %16.9le %16.9le %16.9le %16.9le %16.9le \n",
776  energy/eV,
777  table->GetComponent(0)->FindValue(energy)/(1e-20*m*m),
778  table->GetComponent(1)->FindValue(energy)/(1e-20*m*m),
779  table->GetComponent(2)->FindValue(energy)/(1e-20*m*m),
780  table->GetComponent(3)->FindValue(energy)/(1e-20*m*m),
781  table->GetComponent(4)->FindValue(energy)/(1e-20*m*m),
782  table->GetComponent(0)->FindValue(energy)/(1e-20*m*m)+
783  table->GetComponent(1)->FindValue(energy)/(1e-20*m*m)+
784  table->GetComponent(2)->FindValue(energy)/(1e-20*m*m)+
785  table->GetComponent(3)->FindValue(energy)/(1e-20*m*m)+
786  table->GetComponent(4)->FindValue(energy)/(1e-20*m*m)
787  );
788  energy*=stpEnergy;
789  }
790  fclose (myFile);
791  abort();
792  */
793  //
794  // end of dump
795  //
796  // Test of diff XS
797  // G4double nrj1 = .26827E+04; // in eV
798  // G4double nrj2 = .57991E+03; // in eV
799  // Shells run from 0 to 4
800  // G4cout << DifferentialCrossSection(G4Electron::ElectronDefinition(), nrj1, nrj2, 0)/(1e-20*m*m) << G4endl;
801  // G4cout << DifferentialCrossSection(G4Electron::ElectronDefinition(), nrj1, nrj2, 1)/(1e-20*m*m) << G4endl;
802  // G4cout << DifferentialCrossSection(G4Electron::ElectronDefinition(), nrj1, nrj2, 2)/(1e-20*m*m) << G4endl;
803  // G4cout << DifferentialCrossSection(G4Electron::ElectronDefinition(), nrj1, nrj2, 3)/(1e-20*m*m) << G4endl;
804  // G4cout << DifferentialCrossSection(G4Electron::ElectronDefinition(), nrj1, nrj2, 4)/(1e-20*m*m) << G4endl;
805  // abort();
806  //
807 
808  while (i>0)
809  {
810  i--;
811  valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
812  value += valuesBuffer[i];
813  }
814 
815  value *= G4UniformRand();
816 
817  i = n;
818 
819  while (i > 0)
820  {
821  i--;
822  if (valuesBuffer[i] > value)
823  {
824  delete[] valuesBuffer;
825 
826  return i;
827  }
828  value -= valuesBuffer[i];
829  }
830 
831  if (valuesBuffer) delete[] valuesBuffer;
832 
833  }
834  }
835  else
836  {
837  G4Exception("G4DNACPA100IonisationModel::RandomSelect","em0002",
838  FatalException,"Model not applicable to particle type.");
839  }
840 
841  return level;
842 }
843 
844 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
845 
846 G4double G4DNACPA100IonisationModel::RandomizeEjectedElectronEnergyFromCumulatedDcs
847 (G4ParticleDefinition* particleDefinition, G4double k, G4int shell)
848 {
849  //G4cout << "*** FAST computation for " << " " << particleDefinition->GetParticleName() << G4endl;
850 
851  G4double secondaryElectronKineticEnergy = 0.;
852 
853  secondaryElectronKineticEnergy=
854  RandomTransferedEnergy(particleDefinition, k/eV, shell)*eV-waterStructure.IonisationEnergy(shell);
855 
856  //G4cout << RandomTransferedEnergy(particleDefinition, k/eV, shell) << G4endl;
857  if (secondaryElectronKineticEnergy<0.) return 0.;
858  //
859 
860  return secondaryElectronKineticEnergy;
861 }
862 
863 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
864 
865 G4double G4DNACPA100IonisationModel::RandomTransferedEnergy
866 (G4ParticleDefinition* particleDefinition,G4double k, G4int ionizationLevelIndex)
867 {
868 
869  G4double random = G4UniformRand();
870 
871  G4double nrj = 0.;
872 
873  G4double valueK1 = 0;
874  G4double valueK2 = 0;
875  G4double valuePROB21 = 0;
876  G4double valuePROB22 = 0;
877  G4double valuePROB12 = 0;
878  G4double valuePROB11 = 0;
879 
880  G4double nrjTransf11 = 0;
881  G4double nrjTransf12 = 0;
882  G4double nrjTransf21 = 0;
883  G4double nrjTransf22 = 0;
884 
885  if (particleDefinition == G4Electron::ElectronDefinition())
886  {
887 
888  // Protection against out of boundary access
889  if (k==eTdummyVec.back()) k=k*(1.-1e-12);
890  //
891 
892  // k should be in eV
893 
894  std::vector<G4double>::iterator k2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
895 
896  std::vector<G4double>::iterator k1 = k2-1;
897 
898  /*
899  G4cout << "----> k=" << k
900  << " " << *k1
901  << " " << *k2
902  << " " << random
903  << " " << ionizationLevelIndex
904  << " " << eProbaShellMap[ionizationLevelIndex][(*k1)].back()
905  << " " << eProbaShellMap[ionizationLevelIndex][(*k2)].back()
906  << G4endl;
907  */
908 
909  // SI : the following condition avoids situations where random >last vector element
910 
911  if ( random <= eProbaShellMap[ionizationLevelIndex][(*k1)].back()
912  && random <= eProbaShellMap[ionizationLevelIndex][(*k2)].back() )
913 
914  {
915 
916  std::vector<G4double>::iterator prob12 =
917  std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
918  eProbaShellMap[ionizationLevelIndex][(*k1)].end(), random);
919 
920  std::vector<G4double>::iterator prob11 = prob12-1;
921 
922 
923  std::vector<G4double>::iterator prob22 =
924  std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
925  eProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
926 
927  std::vector<G4double>::iterator prob21 = prob22-1;
928 
929  valueK1 =*k1;
930  valueK2 =*k2;
931  valuePROB21 =*prob21;
932  valuePROB22 =*prob22;
933  valuePROB12 =*prob12;
934  valuePROB11 =*prob11;
935 
936 
937  /*
938  G4cout << " " << random << " " << valuePROB11 << " "
939  << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
940  */
941 
942  nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
943  nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
944  nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
945  nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
946 
947  /*
948  G4cout << " " << ionizationLevelIndex << " "
949  << random << " " <<valueK1 << " " << valueK2 << G4endl;
950 
951  G4cout << " " << random << " " << nrjTransf11 << " "
952  << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
953  */
954 
955  }
956 
957 
958  // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
959 
960  if ( random > eProbaShellMap[ionizationLevelIndex][(*k1)].back() )
961  {
962 
963  std::vector<G4double>::iterator prob22 =
964 
965  std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
966  eProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
967 
968  std::vector<G4double>::iterator prob21 = prob22-1;
969 
970  valueK1 =*k1;
971  valueK2 =*k2;
972  valuePROB21 =*prob21;
973  valuePROB22 =*prob22;
974 
975  //G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
976 
977  nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
978  nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
979 
980  G4double interpolatedvalue2 = Interpolate(valuePROB21, valuePROB22, random, nrjTransf21, nrjTransf22);
981 
982  // zero is explicitely set
983 
984  G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
985 
986  /*
987  G4cout << " " << ionizationLevelIndex << " "
988  << random << " " <<valueK1 << " " << valueK2 << G4endl;
989 
990  G4cout << " " << random << " " << nrjTransf11 << " "
991  << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
992 
993  G4cout << "ici" << " " << value << G4endl;
994  */
995 
996  return value;
997  }
998 
999  }
1000 
1001  // End electron case
1002 
1003  G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21 * nrjTransf22;
1004 
1005  //G4cout << "nrjTransfProduct=" << nrjTransfProduct << G4endl;
1006 
1007  if (nrjTransfProduct != 0.)
1008  {
1009  nrj = QuadInterpolator( valuePROB11, valuePROB12,
1010  valuePROB21, valuePROB22,
1011  nrjTransf11, nrjTransf12,
1012  nrjTransf21, nrjTransf22,
1013  valueK1, valueK2,
1014  k, random);
1015  }
1016 
1017  //G4cout << nrj << endl;
1018 
1019  return nrj ;
1020 }
1021 
1022 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1023 
1024 G4double G4DNACPA100IonisationModel::RandomizeEjectedElectronEnergyFromCompositionSampling
1025 (G4ParticleDefinition*, G4double tt, G4int shell)
1026 {
1027  //G4cout << "*** Rejection method for " << " " << particleDefinition->GetParticleName() << G4endl;
1028 
1029  // ***** METHOD 1 ***** (sequential)
1030  /*
1031 
1032  // ww is KINETIC ENERGY OF SECONDARY ELECTRON
1033  G4double un=1.;
1034  G4double deux=2.;
1035 
1036  G4double bb = waterStructure.IonisationEnergy(shell);
1037  G4double uu = waterStructure.UEnergy(shell);
1038 
1039  if (tt<=bb) return 0.;
1040 
1041  G4double t = tt/bb;
1042  G4double u = uu/bb;
1043  G4double tp1 = t + un;
1044  G4double tu1 = t + u + un;
1045  G4double tm1 = t - un;
1046  G4double tp12 = tp1 * tp1;
1047  G4double dlt = std::log(t);
1048 
1049  G4double a1 = t * tm1 / tu1 / tp12;
1050  G4double a2 = tm1 / tu1 / t / tp1 / deux;
1051  G4double a3 = dlt * (tp12 - deux * deux ) / tu1 / tp12;
1052  G4double ato = a1 + a2 + a3;
1053 
1054  // 15
1055 
1056  G4double r1 =G4UniformRand();
1057  G4double r2 =G4UniformRand();
1058  G4double r3 =G4UniformRand();
1059 
1060  while (r1<=a1/ato)
1061  {
1062  G4double fx1=r2*tm1/tp1;
1063  G4double wx1=un/(un-fx1)-un;
1064  G4double gx1=(t-wx1)/t;
1065  if(r3 <= gx1) return wx1*bb;
1066 
1067  r1 =G4UniformRand();
1068  r2 =G4UniformRand();
1069  r3 =G4UniformRand();
1070 
1071  }
1072 
1073  // 20
1074 
1075  while (r1<=(a1+a2)/ato)
1076  {
1077  G4double fx2=tp1+r2*tm1;
1078  G4double wx2=t-t*tp1/fx2;
1079  G4double gx2=deux*(un-(t-wx2)/tp1);
1080  if(r3 <= gx2) return wx2*bb;
1081 
1082  // REPEAT 15
1083  r1 =G4UniformRand();
1084  r2 =G4UniformRand();
1085  r3 =G4UniformRand();
1086 
1087  while (r1<=a1/ato)
1088  {
1089  G4double fx1=r2*tm1/tp1;
1090  G4double wx1=un/(un-fx1)-un;
1091  G4double gx1=(t-wx1)/t;
1092  if(r3 <= gx1) return wx1*bb;
1093  r1 =G4UniformRand();
1094  r2 =G4UniformRand();
1095  r3 =G4UniformRand();
1096  }
1097  // END 15
1098 
1099  }
1100 
1101  // 30
1102 
1103  G4double wx3=std::sqrt(un/(un-r2*(tp12-deux*deux)/tp12))-un;
1104  G4double gg3=(wx3+un)/(t-wx3);
1105  G4double gx3=(un+gg3*gg3*gg3)/deux;
1106 
1107  while (r3>gx3)
1108  {
1109 
1110  // 15
1111 
1112  r1 =G4UniformRand();
1113  r2 =G4UniformRand();
1114  r3 =G4UniformRand();
1115 
1116  while (r1<=a1/ato)
1117  {
1118  G4double fx1=r2*tm1/tp1;
1119  G4double wx1=un/(un-fx1)-un;
1120  G4double gx1=(t-wx1)/t;
1121  if(r3 <= gx1) return wx1*bb;
1122 
1123  r1 =G4UniformRand();
1124  r2 =G4UniformRand();
1125  r3 =G4UniformRand();
1126 
1127  }
1128 
1129  // 20
1130 
1131  while (r1<=(a1+a2)/ato)
1132  {
1133  G4double fx2=tp1+r2*tm1;
1134  G4double wx2=t-t*tp1/fx2;
1135  G4double gx2=deux*(un-(t-wx2)/tp1);
1136  if(r3 <= gx2)return wx2*bb;
1137 
1138  // REPEAT 15
1139  r1 =G4UniformRand();
1140  r2 =G4UniformRand();
1141  r3 =G4UniformRand();
1142 
1143  while (r1<=a1/ato)
1144  {
1145  G4double fx1=r2*tm1/tp1;
1146  G4double wx1=un/(un-fx1)-un;
1147  G4double gx1=(t-wx1)/t;
1148  if(r3 <= gx1) return wx1*bb;
1149 
1150  r1 =G4UniformRand();
1151  r2 =G4UniformRand();
1152  r3 =G4UniformRand();
1153  }
1154  //
1155 
1156  }
1157 
1158  wx3=std::sqrt(un/(un-r2*(tp12-deux*deux)/tp12))-un;
1159  gg3=(wx3+un)/(t-wx3);
1160  gx3=(un+gg3*gg3*gg3)/deux;
1161 
1162  }
1163 
1164  //
1165 
1166  return wx3*bb;
1167  */
1168 
1169  // ***** METHOD by M. C. Bordage ***** (optimized)
1170 
1171  G4double un=1.;
1172  G4double deux=2.;
1173 
1174  G4double bb = waterStructure.IonisationEnergy(shell);
1175  G4double uu = waterStructure.UEnergy(shell);
1176 
1177  if (tt<=bb) return 0.;
1178 
1179  G4double t = tt/bb;
1180  G4double u = uu/bb;
1181  G4double tp1 = t + un;
1182  G4double tu1 = t + u + un;
1183  G4double tm1 = t - un;
1184  G4double tp12 = tp1 * tp1;
1185  G4double dlt = std::log(t);
1186 
1187  G4double a1 = t * tm1 / tu1 / tp12;
1188  G4double a2 = tm1 / tu1 / t / tp1 / deux;
1189  G4double a3 = dlt * (tp12 - deux * deux ) / tu1 / tp12;
1190  G4double ato = a1 + a2 + a3;
1191 
1192  G4double A1 = a1/ato;
1193  G4double A2 = (a1+a2)/ato;
1194  G4int F = 0;
1195  G4double fx=0;
1196  G4double gx=0;
1197  G4double gg=0;
1198  G4double wx=0;
1199 
1200  G4double r1=0;
1201  G4double r2=0;
1202  G4double r3=0;
1203 
1204  //
1205 
1206  do
1207  {
1208  r1 =G4UniformRand();
1209  r2 =G4UniformRand();
1210  r3 =G4UniformRand();
1211 
1212  if (r1>A2)
1213  F=3;
1214  else if ((r1>A1) && (r1< A2))
1215  F=2;
1216  else
1217  F=1;
1218 
1219  switch (F)
1220  {
1221  case 1:
1222  {
1223  fx=r2*tm1/tp1;
1224  wx=un/(un-fx)-un;
1225  gx=(t-wx)/t;
1226  break;
1227  }
1228 
1229  case 2:
1230  {
1231  fx=tp1+r2*tm1;
1232  wx=t-t*tp1/fx;
1233  gx=deux*(un-(t-wx)/tp1);
1234  break;
1235  }
1236 
1237  case 3:
1238  {
1239  fx=un-r2*(tp12-deux*deux)/tp12;
1240  wx=sqrt(un/fx)-un;
1241  gg=(wx+un)/(t-wx);
1242  gx=(un+gg*gg*gg)/deux;
1243  break;
1244  }
1245  } // switch
1246 
1247  } while (r3>gx);
1248 
1249  return wx*bb;
1250 
1251 }