Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNABornIonisationModel.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4DNABornIonisationModel.cc 70823 2013-06-06 08:26:50Z gcosmo $
27 //
28 
30 #include "G4PhysicalConstants.hh"
31 #include "G4SystemOfUnits.hh"
32 #include "G4UAtomicDeexcitation.hh"
33 #include "G4LossTableManager.hh"
34 #include "G4DNAChemistryManager.hh"
36 
37 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
38 
39 using namespace std;
40 
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42 
44  const G4String& nam)
45  :G4VEmModel(nam),isInitialised(false)
46 {
47  verboseLevel= 0;
48  // Verbosity scale:
49  // 0 = nothing
50  // 1 = warning for energy non-conservation
51  // 2 = details of energy budget
52  // 3 = calculation of cross sections, file openings, sampling of atoms
53  // 4 = entering in methods
54 
55  if( verboseLevel>0 )
56  {
57  G4cout << "Born ionisation model is constructed " << G4endl;
58  }
59 
60  //Mark this model as "applicable" for atomic deexcitation
61  SetDeexcitationFlag(true);
62  fAtomDeexcitation = 0;
64  fpMolWaterDensity = 0;
65 
66  // Selection of computation method
67 
68  fasterCode = false;
69 }
70 
71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
72 
74 {
75  // Cross section
76 
77  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
78  for (pos = tableData.begin(); pos != tableData.end(); ++pos)
79  {
80  G4DNACrossSectionDataSet* table = pos->second;
81  delete table;
82  }
83 
84  // Final state
85 
86  eVecm.clear();
87  pVecm.clear();
88 
89 }
90 
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92 
94  const G4DataVector& /*cuts*/)
95 {
96 
97  if (verboseLevel > 3)
98  G4cout << "Calling G4DNABornIonisationModel::Initialise()" << G4endl;
99 
100  // Energy limits
101 
102  G4String fileElectron("dna/sigma_ionisation_e_born");
103  G4String fileProton("dna/sigma_ionisation_p_born");
104 
107 
110 
111  G4double scaleFactor = (1.e-22 / 3.343) * m*m;
112 
113  char *path = getenv("G4LEDATA");
114 
115  // *** ELECTRON
116 
117  electron = electronDef->GetParticleName();
118 
119  tableFile[electron] = fileElectron;
120 
121  lowEnergyLimit[electron] = 11. * eV;
122  highEnergyLimit[electron] = 1. * MeV;
123 
124  // Cross section
125 
127  tableE->LoadData(fileElectron);
128 
129  tableData[electron] = tableE;
130 
131  // Final state
132 
133  std::ostringstream eFullFileName;
134 
135  if (fasterCode) eFullFileName << path << "/dna/sigmadiff_cumulated_ionisation_e_born.dat";
136  if (!fasterCode) eFullFileName << path << "/dna/sigmadiff_ionisation_e_born.dat";
137 
138  std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
139 
140  if (!eDiffCrossSection)
141  {
142  if (fasterCode) G4Exception("G4DNABornIonisationModel::Initialise","em0003",
143  FatalException,"Missing data file:/dna/sigmadiff_cumulated_ionisation_e_born.dat");
144 
145  if (!fasterCode) G4Exception("G4DNABornIonisationModel::Initialise","em0003",
146  FatalException,"Missing data file:/dna/sigmadiff_ionisation_e_born.dat");
147  }
148 
149  eTdummyVec.push_back(0.);
150  while(!eDiffCrossSection.eof())
151  {
152  double tDummy;
153  double eDummy;
154  eDiffCrossSection>>tDummy>>eDummy;
155  if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
156  for (int j=0; j<5; j++)
157  {
158  eDiffCrossSection>>eDiffCrossSectionData[j][tDummy][eDummy];
159 
160  if (fasterCode)
161  {
162  eNrjTransfData[j][tDummy][eDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
163  eProbaShellMap[j][tDummy].push_back(eDiffCrossSectionData[j][tDummy][eDummy]);
164  }
165 
166  // SI - only if eof is not reached
167  if (!eDiffCrossSection.eof() && !fasterCode) eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
168 
169  if (!fasterCode) eVecm[tDummy].push_back(eDummy);
170 
171  }
172  }
173 
174  // *** PROTON
175 
176  proton = protonDef->GetParticleName();
177 
178  tableFile[proton] = fileProton;
179 
180  lowEnergyLimit[proton] = 500. * keV;
181  highEnergyLimit[proton] = 100. * MeV;
182 
183  // Cross section
184 
186  tableP->LoadData(fileProton);
187 
188  tableData[proton] = tableP;
189 
190  // Final state
191 
192  std::ostringstream pFullFileName;
193 
194  if (fasterCode) pFullFileName << path << "/dna/sigmadiff_cumulated_ionisation_p_born.dat";
195 
196  if (!fasterCode) pFullFileName << path << "/dna/sigmadiff_ionisation_p_born.dat";
197 
198  std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
199 
200  if (!pDiffCrossSection)
201  {
202  if (fasterCode) G4Exception("G4DNABornIonisationModel::Initialise","em0003",
203  FatalException,"Missing data file:/dna/sigmadiff_cumulated_ionisation_p_born.dat");
204 
205  if (!fasterCode) G4Exception("G4DNABornIonisationModel::Initialise","em0003",
206  FatalException,"Missing data file:/dna/sigmadiff_ionisation_p_born.dat");
207  }
208 
209  pTdummyVec.push_back(0.);
210  while(!pDiffCrossSection.eof())
211  {
212  double tDummy;
213  double eDummy;
214  pDiffCrossSection>>tDummy>>eDummy;
215  if (tDummy != pTdummyVec.back()) pTdummyVec.push_back(tDummy);
216  for (int j=0; j<5; j++)
217  {
218  pDiffCrossSection>>pDiffCrossSectionData[j][tDummy][eDummy];
219 
220  if (fasterCode)
221  {
222  pNrjTransfData[j][tDummy][pDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
223  pProbaShellMap[j][tDummy].push_back(pDiffCrossSectionData[j][tDummy][eDummy]);
224  }
225 
226  // SI - only if eof is not reached !
227  if (!pDiffCrossSection.eof() && !fasterCode) pDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
228 
229  if (!fasterCode) pVecm[tDummy].push_back(eDummy);
230  }
231  }
232 
233  //
234 
235  if (particle==electronDef)
236  {
237  SetLowEnergyLimit(lowEnergyLimit[electron]);
238  SetHighEnergyLimit(highEnergyLimit[electron]);
239  }
240 
241  if (particle==protonDef)
242  {
243  SetLowEnergyLimit(lowEnergyLimit[proton]);
244  SetHighEnergyLimit(highEnergyLimit[proton]);
245  }
246 
247  if( verboseLevel>0 )
248  {
249  G4cout << "Born ionisation model is initialized " << G4endl
250  << "Energy range: "
251  << LowEnergyLimit() / eV << " eV - "
252  << HighEnergyLimit() / keV << " keV for "
253  << particle->GetParticleName()
254  << G4endl;
255  }
256 
257  // Initialize water density pointer
259 
260  //
261  fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
262 
263  if (isInitialised) { return; }
265  isInitialised = true;
266 }
267 
268 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
269 
271  const G4ParticleDefinition* particleDefinition,
272  G4double ekin,
273  G4double,
274  G4double)
275 {
276  if (verboseLevel > 3)
277  G4cout << "Calling CrossSectionPerVolume() of G4DNABornIonisationModel" << G4endl;
278 
279  if (
280  particleDefinition != G4Proton::ProtonDefinition()
281  &&
282  particleDefinition != G4Electron::ElectronDefinition()
283  )
284 
285  return 0;
286 
287  // Calculate total cross section for model
288 
289  G4double lowLim = 0;
290  G4double highLim = 0;
291  G4double sigma=0;
292 
293  G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
294 
295  if(waterDensity!= 0.0)
296  // if (material == nistwater || material->GetBaseMaterial() == nistwater)
297  {
298  const G4String& particleName = particleDefinition->GetParticleName();
299 
300  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
301  pos1 = lowEnergyLimit.find(particleName);
302  if (pos1 != lowEnergyLimit.end())
303  {
304  lowLim = pos1->second;
305  }
306 
307  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
308  pos2 = highEnergyLimit.find(particleName);
309  if (pos2 != highEnergyLimit.end())
310  {
311  highLim = pos2->second;
312  }
313 
314  if (ekin >= lowLim && ekin < highLim)
315  {
316  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
317  pos = tableData.find(particleName);
318 
319  if (pos != tableData.end())
320  {
321  G4DNACrossSectionDataSet* table = pos->second;
322  if (table != 0)
323  {
324  sigma = table->FindValue(ekin);
325  }
326  }
327  else
328  {
329  G4Exception("G4DNABornIonisationModel::CrossSectionPerVolume","em0002",
330  FatalException,"Model not applicable to particle type.");
331  }
332  }
333 
334  if (verboseLevel > 2)
335  {
336  G4cout << "__________________________________" << G4endl;
337  G4cout << "G4DNABornIonisationModel - XS INFO START" << G4endl;
338  G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
339  G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
340  G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
341  G4cout << "G4DNABornIonisationModel - XS INFO END" << G4endl;
342  }
343 
344  } // if (waterMaterial)
345 
346  return sigma*waterDensity;
347 }
348 
349 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
350 
351 void G4DNABornIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
352  const G4MaterialCutsCouple* ,//must be set!
353  const G4DynamicParticle* particle,
354  G4double,
355  G4double)
356 {
357 
358  if (verboseLevel > 3)
359  G4cout << "Calling SampleSecondaries() of G4DNABornIonisationModel" << G4endl;
360 
361  G4double lowLim = 0;
362  G4double highLim = 0;
363 
364  G4double k = particle->GetKineticEnergy();
365 
366  const G4String& particleName = particle->GetDefinition()->GetParticleName();
367 
368  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
369  pos1 = lowEnergyLimit.find(particleName);
370 
371  if (pos1 != lowEnergyLimit.end())
372  {
373  lowLim = pos1->second;
374  }
375 
376  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
377  pos2 = highEnergyLimit.find(particleName);
378 
379  if (pos2 != highEnergyLimit.end())
380  {
381  highLim = pos2->second;
382  }
383 
384  if (k >= lowLim && k < highLim)
385  {
386  G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
387  G4double particleMass = particle->GetDefinition()->GetPDGMass();
388  G4double totalEnergy = k + particleMass;
389  G4double pSquare = k * (totalEnergy + particleMass);
390  G4double totalMomentum = std::sqrt(pSquare);
391 
392  G4int ionizationShell = -1;
393 
394  if (!fasterCode) ionizationShell = RandomSelect(k,particleName);
395 
396  // SI: The following protection is necessary to avoid infinite loops :
397  // sigmadiff_ionisation_e_born.dat has non zero partial xs at 18 eV for shell 3 (ionizationShell ==2)
398  // sigmadiff_cumulated_ionisation_e_born.dat has zero cumulated partial xs at 18 eV for shell 3 (ionizationShell ==2)
399  // this is due to the fact that the max allowed transfered energy is (18+10.79)/2=17.025 eV and only transfered energies
400  // strictly above this value have non zero partial xs in sigmadiff_ionisation_e_born.dat (starting at trans = 17.12 eV)
401 
402  if (fasterCode)
403  do
404  {
405  ionizationShell = RandomSelect(k,particleName);
406  } while (k<19*eV && ionizationShell==2 && particle->GetDefinition()==G4Electron::ElectronDefinition());
407 
408 
409  // AM: sample deexcitation
410  // here we assume that H_{2}O electronic levels are the same of Oxigen.
411  // this can be considered true with a rough 10% error in energy on K-shell,
412 
413  G4int secNumberInit = 0; // need to know at a certain point the enrgy of secondaries
414  G4int secNumberFinal = 0; // So I'll make the diference and then sum the energies
415 
417  bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
418 
419  if(fAtomDeexcitation) {
420  G4int Z = 8;
422 
423  if (ionizationShell <5 && ionizationShell >1)
424  {
425  as = G4AtomicShellEnumerator(4-ionizationShell);
426  }
427  else if (ionizationShell <2)
428  {
429  as = G4AtomicShellEnumerator(3);
430  }
431 
432  // FOR DEBUG ONLY
433  // if (ionizationShell == 4) {
434  //
435  // G4cout << "Z: " << Z << " as: " << as
436  // << " ionizationShell: " << ionizationShell << " bindingEnergy: "<< bindingEnergy/eV << G4endl;
437  // G4cout << "Press <Enter> key to continue..." << G4endl;
438  // G4cin.ignore();
439  // }
440 
441  const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
442  secNumberInit = fvect->size();
443  fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
444  secNumberFinal = fvect->size();
445  }
446 
447  G4double secondaryKinetic=-1000*eV;
448 
449  if (!fasterCode) secondaryKinetic = RandomizeEjectedElectronEnergy(particle->GetDefinition(),k,ionizationShell);
450 
451  if (fasterCode)
452  do
453  {
454  secondaryKinetic = RandomizeEjectedElectronEnergyFromCumulatedDcs(particle->GetDefinition(),k,ionizationShell);
455 
456  }
457  while (secondaryKinetic<0) ;
458 
459  G4double cosTheta = 0.;
460  G4double phi = 0.;
461  RandomizeEjectedElectronDirection(particle->GetDefinition(), k,secondaryKinetic, cosTheta, phi);
462 
463  G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
464  G4double dirX = sinTheta*std::cos(phi);
465  G4double dirY = sinTheta*std::sin(phi);
466  G4double dirZ = cosTheta;
467  G4ThreeVector deltaDirection(dirX,dirY,dirZ);
468  deltaDirection.rotateUz(primaryDirection);
469 
470  if (particle->GetDefinition() == G4Electron::ElectronDefinition())
471  {
472  G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
473 
474  G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
475  G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
476  G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
477  G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
478  finalPx /= finalMomentum;
479  finalPy /= finalMomentum;
480  finalPz /= finalMomentum;
481 
482  G4ThreeVector direction;
483  direction.set(finalPx,finalPy,finalPz);
484 
486  }
487 
488  else fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection) ;
489 
490  // note thta secondaryKinetic is the nergy of the delta ray, not of all secondaries.
491  G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
492  G4double deexSecEnergy = 0;
493  for (G4int j=secNumberInit; j < secNumberFinal; j++) {
494 
495  deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
496 
497  }
498 
500 
501  fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy-secondaryKinetic-deexSecEnergy);
502 
503  G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
504  fvect->push_back(dp);
505 
506  const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
508  ionizationShell,
509  theIncomingTrack);
510  }
511 
512 }
513 
514 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
515 
516 G4double G4DNABornIonisationModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
517  G4double k, G4int shell)
518 {
519  // G4cout << "*** SLOW computation for " << " " << particleDefinition->GetParticleName() << G4endl;
520 
521  if (particleDefinition == G4Electron::ElectronDefinition())
522  {
523  G4double maximumEnergyTransfer=0.;
524  if ((k+waterStructure.IonisationEnergy(shell))/2. > k) maximumEnergyTransfer=k;
525  else maximumEnergyTransfer = (k+waterStructure.IonisationEnergy(shell))/2.;
526 
527  // SI : original method
528  /*
529  G4double crossSectionMaximum = 0.;
530  for(G4double value=waterStructure.IonisationEnergy(shell); value<=maximumEnergyTransfer; value+=0.1*eV)
531  {
532  G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
533  if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
534  }
535  */
536 
537  // SI : alternative method
538 
539  G4double crossSectionMaximum = 0.;
540 
541  G4double minEnergy = waterStructure.IonisationEnergy(shell);
542  G4double maxEnergy = maximumEnergyTransfer;
543  G4int nEnergySteps = 50;
544 
545  G4double value(minEnergy);
546  G4double stpEnergy(std::pow(maxEnergy/value, 1./static_cast<G4double>(nEnergySteps-1)));
547  G4int step(nEnergySteps);
548  while (step>0)
549  {
550  step--;
551  G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
552  if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
553  value*=stpEnergy;
554  }
555  //
556 
557  G4double secondaryElectronKineticEnergy=0.;
558  do
559  {
560  secondaryElectronKineticEnergy = G4UniformRand() * (maximumEnergyTransfer-waterStructure.IonisationEnergy(shell));
561  } while(G4UniformRand()*crossSectionMaximum >
562  DifferentialCrossSection(particleDefinition, k/eV,
563  (secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
564 
565  return secondaryElectronKineticEnergy;
566 
567  }
568 
569  else if (particleDefinition == G4Proton::ProtonDefinition())
570  {
571  G4double maximumKineticEnergyTransfer = 4.* (electron_mass_c2 / proton_mass_c2) * k;
572 
573  G4double crossSectionMaximum = 0.;
574  for (G4double value = waterStructure.IonisationEnergy(shell);
575  value<=4.*waterStructure.IonisationEnergy(shell) ;
576  value+=0.1*eV)
577  {
578  G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
579  if (differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
580  }
581 
582  G4double secondaryElectronKineticEnergy = 0.;
583  do
584  {
585  secondaryElectronKineticEnergy = G4UniformRand() * maximumKineticEnergyTransfer;
586  } while(G4UniformRand()*crossSectionMaximum >=
587  DifferentialCrossSection(particleDefinition, k/eV,
588  (secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
589 
590  return secondaryElectronKineticEnergy;
591  }
592 
593  return 0;
594 }
595 
596 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
597 
598 void G4DNABornIonisationModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
599  G4double k,
600  G4double secKinetic,
601  G4double & cosTheta,
602  G4double & phi )
603 {
604  if (particleDefinition == G4Electron::ElectronDefinition())
605  {
606  phi = twopi * G4UniformRand();
607  if (secKinetic < 50.*eV) cosTheta = (2.*G4UniformRand())-1.;
608  else if (secKinetic <= 200.*eV)
609  {
610  if (G4UniformRand() <= 0.1) cosTheta = (2.*G4UniformRand())-1.;
611  else cosTheta = G4UniformRand()*(std::sqrt(2.)/2);
612  }
613  else
614  {
615  G4double sin2O = (1.-secKinetic/k) / (1.+secKinetic/(2.*electron_mass_c2));
616  cosTheta = std::sqrt(1.-sin2O);
617  }
618  }
619 
620  else if (particleDefinition == G4Proton::ProtonDefinition())
621  {
622  G4double maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
623  phi = twopi * G4UniformRand();
624 
625  // cosTheta = std::sqrt(secKinetic / maxSecKinetic);
626 
627  // Restriction below 100 eV from Emfietzoglou (2000)
628 
629  if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
630  else cosTheta = (2.*G4UniformRand())-1.;
631 
632  }
633 }
634 
635 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
636 
638  G4double k,
639  G4double energyTransfer,
640  G4int ionizationLevelIndex)
641 {
642  G4double sigma = 0.;
643 
644  if (energyTransfer >= waterStructure.IonisationEnergy(ionizationLevelIndex))
645  {
646  G4double valueT1 = 0;
647  G4double valueT2 = 0;
648  G4double valueE21 = 0;
649  G4double valueE22 = 0;
650  G4double valueE12 = 0;
651  G4double valueE11 = 0;
652 
653  G4double xs11 = 0;
654  G4double xs12 = 0;
655  G4double xs21 = 0;
656  G4double xs22 = 0;
657 
658  if (particleDefinition == G4Electron::ElectronDefinition())
659  {
660  // k should be in eV and energy transfer eV also
661 
662  std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
663 
664  std::vector<double>::iterator t1 = t2-1;
665 
666  // SI : the following condition avoids situations where energyTransfer >last vector element
667  if (energyTransfer <= eVecm[(*t1)].back() && energyTransfer <= eVecm[(*t2)].back() )
668  {
669  std::vector<double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), energyTransfer);
670  std::vector<double>::iterator e11 = e12-1;
671 
672  std::vector<double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), energyTransfer);
673  std::vector<double>::iterator e21 = e22-1;
674 
675  valueT1 =*t1;
676  valueT2 =*t2;
677  valueE21 =*e21;
678  valueE22 =*e22;
679  valueE12 =*e12;
680  valueE11 =*e11;
681 
682  xs11 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
683  xs12 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
684  xs21 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
685  xs22 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
686 
687  }
688 
689  }
690 
691  if (particleDefinition == G4Proton::ProtonDefinition())
692  {
693  // k should be in eV and energy transfer eV also
694  std::vector<double>::iterator t2 = std::upper_bound(pTdummyVec.begin(),pTdummyVec.end(), k);
695  std::vector<double>::iterator t1 = t2-1;
696 
697  std::vector<double>::iterator e12 = std::upper_bound(pVecm[(*t1)].begin(),pVecm[(*t1)].end(), energyTransfer);
698  std::vector<double>::iterator e11 = e12-1;
699 
700  std::vector<double>::iterator e22 = std::upper_bound(pVecm[(*t2)].begin(),pVecm[(*t2)].end(), energyTransfer);
701  std::vector<double>::iterator e21 = e22-1;
702 
703  valueT1 =*t1;
704  valueT2 =*t2;
705  valueE21 =*e21;
706  valueE22 =*e22;
707  valueE12 =*e12;
708  valueE11 =*e11;
709 
710  xs11 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
711  xs12 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
712  xs21 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
713  xs22 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
714 
715  }
716 
717  G4double xsProduct = xs11 * xs12 * xs21 * xs22;
718  if (xsProduct != 0.)
719  {
720  sigma = QuadInterpolator( valueE11, valueE12,
721  valueE21, valueE22,
722  xs11, xs12,
723  xs21, xs22,
724  valueT1, valueT2,
725  k, energyTransfer);
726  }
727 
728  }
729 
730  return sigma;
731 }
732 
733 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
734 
735 G4double G4DNABornIonisationModel::Interpolate( G4double e1,
736  G4double e2,
737  G4double e,
738  G4double xs1,
739  G4double xs2)
740 {
741 
742  G4double value = 0.;
743 
744  // Log-log interpolation by default
745 
746  if (e1!=0 && e2!=0 && (std::log10(e2)-std::log10(e1)) !=0 && !fasterCode)
747  {
748  G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
749  G4double b = std::log10(xs2) - a*std::log10(e2);
750  G4double sigma = a*std::log10(e) + b;
751  value = (std::pow(10.,sigma));
752  }
753 
754  // Switch to lin-lin interpolation
755  /*
756  if ((e2-e1)!=0)
757  {
758  G4double d1 = xs1;
759  G4double d2 = xs2;
760  value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
761  }
762  */
763 
764  // Switch to log-lin interpolation for faster code
765 
766  if ((e2-e1)!=0 && xs1 !=0 && xs2 !=0 && fasterCode )
767  {
768  G4double d1 = std::log10(xs1);
769  G4double d2 = std::log10(xs2);
770  value = std::pow(10.,(d1 + (d2 - d1)*(e - e1)/ (e2 - e1)) );
771  }
772 
773  // Switch to lin-lin interpolation for faster code
774  // in case one of xs1 or xs2 (=cum proba) value is zero
775 
776  if ((e2-e1)!=0 && (xs1 ==0 || xs2 ==0) && fasterCode )
777  {
778  G4double d1 = xs1;
779  G4double d2 = xs2;
780  value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
781  }
782 
783  /*
784  G4cout
785  << e1 << " "
786  << e2 << " "
787  << e << " "
788  << xs1 << " "
789  << xs2 << " "
790  << value
791  << G4endl;
792  */
793 
794  return value;
795 }
796 
797 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
798 
799 G4double G4DNABornIonisationModel::QuadInterpolator(G4double e11, G4double e12,
800  G4double e21, G4double e22,
801  G4double xs11, G4double xs12,
802  G4double xs21, G4double xs22,
803  G4double t1, G4double t2,
804  G4double t, G4double e)
805 {
806  G4double interpolatedvalue1 = Interpolate(e11, e12, e, xs11, xs12);
807  G4double interpolatedvalue2 = Interpolate(e21, e22, e, xs21, xs22);
808  G4double value = Interpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
809 
810  return value;
811 }
812 
813 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
814 
815 G4int G4DNABornIonisationModel::RandomSelect(G4double k, const G4String& particle )
816 {
817  G4int level = 0;
818 
819  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
820  pos = tableData.find(particle);
821 
822  if (pos != tableData.end())
823  {
824  G4DNACrossSectionDataSet* table = pos->second;
825 
826  if (table != 0)
827  {
828  G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
829  const size_t n(table->NumberOfComponents());
830  size_t i(n);
831  G4double value = 0.;
832 
833  while (i>0)
834  {
835  i--;
836  valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
837  value += valuesBuffer[i];
838  }
839 
840  value *= G4UniformRand();
841 
842  i = n;
843 
844  while (i > 0)
845  {
846  i--;
847 
848  if (valuesBuffer[i] > value)
849  {
850  delete[] valuesBuffer;
851  return i;
852  }
853  value -= valuesBuffer[i];
854  }
855 
856  if (valuesBuffer) delete[] valuesBuffer;
857 
858  }
859  }
860  else
861  {
862  G4Exception("G4DNABornIonisationModel::RandomSelect","em0002",
863  FatalException,"Model not applicable to particle type.");
864  }
865 
866  return level;
867 }
868 
869 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
870 
871 G4double G4DNABornIonisationModel::RandomizeEjectedElectronEnergyFromCumulatedDcs
872 (G4ParticleDefinition* particleDefinition, G4double k, G4int shell)
873 {
874  //G4cout << "*** FAST computation for " << " " << particleDefinition->GetParticleName() << G4endl;
875 
876  G4double secondaryElectronKineticEnergy = 0.;
877 
878  secondaryElectronKineticEnergy=
879  RandomTransferedEnergy(particleDefinition, k/eV, shell)*eV-waterStructure.IonisationEnergy(shell);
880 
881  return secondaryElectronKineticEnergy;
882 }
883 
884 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
885 
886 G4double G4DNABornIonisationModel::RandomTransferedEnergy
887 (G4ParticleDefinition* particleDefinition,G4double k, G4int ionizationLevelIndex)
888 {
889 
890  G4double random = G4UniformRand();
891 
892  G4double nrj = 0.;
893 
894  G4double valueK1 = 0;
895  G4double valueK2 = 0;
896  G4double valuePROB21 = 0;
897  G4double valuePROB22 = 0;
898  G4double valuePROB12 = 0;
899  G4double valuePROB11 = 0;
900 
901  G4double nrjTransf11 = 0;
902  G4double nrjTransf12 = 0;
903  G4double nrjTransf21 = 0;
904  G4double nrjTransf22 = 0;
905 
906  if (particleDefinition == G4Electron::ElectronDefinition())
907  {
908  // k should be in eV
909 
910  std::vector<double>::iterator k2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
911 
912  std::vector<double>::iterator k1 = k2-1;
913 
914  /*
915  G4cout << "----> k=" << k
916  << " " << *k1
917  << " " << *k2
918  << " " << random
919  << " " << ionizationLevelIndex
920  << " " << eProbaShellMap[ionizationLevelIndex][(*k1)].back()
921  << " " << eProbaShellMap[ionizationLevelIndex][(*k2)].back()
922  << G4endl;
923  */
924 
925  // SI : the following condition avoids situations where random >last vector element
926 
927  if ( random <= eProbaShellMap[ionizationLevelIndex][(*k1)].back()
928  && random <= eProbaShellMap[ionizationLevelIndex][(*k2)].back() )
929 
930  {
931  std::vector<double>::iterator prob12 = std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
932  eProbaShellMap[ionizationLevelIndex][(*k1)].end(), random);
933 
934  std::vector<double>::iterator prob11 = prob12-1;
935 
936 
937  std::vector<double>::iterator prob22 = std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
938  eProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
939 
940  std::vector<double>::iterator prob21 = prob22-1;
941 
942  valueK1 =*k1;
943  valueK2 =*k2;
944  valuePROB21 =*prob21;
945  valuePROB22 =*prob22;
946  valuePROB12 =*prob12;
947  valuePROB11 =*prob11;
948 
949 
950  /*
951  G4cout << " " << random << " " << valuePROB11 << " "
952  << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
953  */
954 
955  nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
956  nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
957  nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
958  nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
959 
960  /*
961  G4cout << " " << ionizationLevelIndex << " "
962  << random << " " <<valueK1 << " " << valueK2 << G4endl;
963 
964  G4cout << " " << random << " " << nrjTransf11 << " "
965  << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
966  */
967 
968  }
969 
970 
971  // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
972 
973  if ( random > eProbaShellMap[ionizationLevelIndex][(*k1)].back() )
974 
975  {
976  std::vector<double>::iterator prob22 = std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
977  eProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
978 
979  std::vector<double>::iterator prob21 = prob22-1;
980 
981  valueK1 =*k1;
982  valueK2 =*k2;
983  valuePROB21 =*prob21;
984  valuePROB22 =*prob22;
985 
986  //G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
987 
988  nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
989  nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
990 
991  G4double interpolatedvalue2 = Interpolate(valuePROB21, valuePROB22, random, nrjTransf21, nrjTransf22);
992 
993  // zeros are explicitely set
994 
995  G4double value = Interpolate(0., valueK2, k, 0., interpolatedvalue2);
996 
997  /*
998  G4cout << " " << ionizationLevelIndex << " "
999  << random << " " <<valueK1 << " " << valueK2 << G4endl;
1000 
1001  G4cout << " " << random << " " << nrjTransf11 << " "
1002  << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1003 
1004  G4cout << "ici" << " " << value << G4endl;
1005  */
1006 
1007  return value;
1008  }
1009 
1010  }
1011 
1012 
1013  //
1014 
1015  if (particleDefinition == G4Proton::ProtonDefinition())
1016  {
1017  // k should be in eV
1018 
1019  std::vector<double>::iterator k2 = std::upper_bound(pTdummyVec.begin(),pTdummyVec.end(), k);
1020 
1021  std::vector<double>::iterator k1 = k2-1;
1022 
1023  /*
1024  G4cout << "----> k=" << k
1025  << " " << *k1
1026  << " " << *k2
1027  << " " << random
1028  << " " << ionizationLevelIndex
1029  << " " << pProbaShellMap[ionizationLevelIndex][(*k1)].back()
1030  << " " << pProbaShellMap[ionizationLevelIndex][(*k2)].back()
1031  << G4endl;
1032  */
1033 
1034  // SI : the following condition avoids situations where random > last vector element,
1035  // for eg. when the last element is zero
1036 
1037  if ( random <= pProbaShellMap[ionizationLevelIndex][(*k1)].back()
1038  && random <= pProbaShellMap[ionizationLevelIndex][(*k2)].back() )
1039  {
1040  std::vector<double>::iterator prob12 = std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
1041  pProbaShellMap[ionizationLevelIndex][(*k1)].end(), random);
1042 
1043  std::vector<double>::iterator prob11 = prob12-1;
1044 
1045 
1046  std::vector<double>::iterator prob22 = std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1047  pProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
1048 
1049  std::vector<double>::iterator prob21 = prob22-1;
1050 
1051  valueK1 =*k1;
1052  valueK2 =*k2;
1053  valuePROB21 =*prob21;
1054  valuePROB22 =*prob22;
1055  valuePROB12 =*prob12;
1056  valuePROB11 =*prob11;
1057 
1058  /*
1059  G4cout << " " << random << " " << valuePROB11 << " "
1060  << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1061  */
1062 
1063  nrjTransf11 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1064  nrjTransf12 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1065  nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1066  nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1067 
1068  /*
1069  G4cout << " " << ionizationLevelIndex << " "
1070  << random << " " <<valueK1 << " " << valueK2 << G4endl;
1071 
1072  G4cout << " " << random << " " << nrjTransf11 << " "
1073  << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1074  */
1075  }
1076 
1077  // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
1078 
1079  if ( random > pProbaShellMap[ionizationLevelIndex][(*k1)].back() )
1080 
1081  {
1082  std::vector<double>::iterator prob22 = std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1083  pProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
1084 
1085  std::vector<double>::iterator prob21 = prob22-1;
1086 
1087  valueK1 =*k1;
1088  valueK2 =*k2;
1089  valuePROB21 =*prob21;
1090  valuePROB22 =*prob22;
1091 
1092  //G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1093 
1094  nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1095  nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1096 
1097  G4double interpolatedvalue2 = Interpolate(valuePROB21, valuePROB22, random, nrjTransf21, nrjTransf22);
1098 
1099  // zeros are explicitely set
1100 
1101  G4double value = Interpolate(0., valueK2, k, 0., interpolatedvalue2);
1102 
1103  /*
1104  G4cout << " " << ionizationLevelIndex << " "
1105  << random << " " <<valueK1 << " " << valueK2 << G4endl;
1106 
1107  G4cout << " " << random << " " << nrjTransf11 << " "
1108  << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1109 
1110  G4cout << "ici" << " " << value << G4endl;
1111  */
1112 
1113  return value;
1114  }
1115 
1116  }
1117 
1118  // End electron and proton cases
1119 
1120  G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21 * nrjTransf22;
1121 
1122  //G4cout << "nrjTransfProduct=" << nrjTransfProduct << G4endl;
1123 
1124  if (nrjTransfProduct != 0.)
1125  {
1126  nrj = QuadInterpolator( valuePROB11, valuePROB12,
1127  valuePROB21, valuePROB22,
1128  nrjTransf11, nrjTransf12,
1129  nrjTransf21, nrjTransf22,
1130  valueK1, valueK2,
1131  k, random);
1132  }
1133 
1134  //G4cout << nrj << endl;
1135 
1136  return nrj ;
1137 }
void set(double x, double y, double z)
static G4Electron * ElectronDefinition()
Definition: G4Electron.cc:89
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:599
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
static G4LossTableManager * Instance()
G4double GetKineticEnergy() const
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:592
double x() const
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:578
size_t GetIndex() const
Definition: G4Material.hh:260
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
virtual G4bool LoadData(const G4String &argFileName)
G4ParticleDefinition * GetDefinition() const
int G4int
Definition: G4Types.hh:78
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
const G4String & GetParticleName() const
double z() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:683
string material
Definition: eplot.py:19
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
const std::vector< double > * GetNumMolPerVolTableFor(const G4Material *) const
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
float proton_mass_c2
Definition: hepunit.py:275
G4DNABornIonisationModel(const G4ParticleDefinition *p=0, const G4String &nam="DNABornIonisationModel")
float electron_mass_c2
Definition: hepunit.py:274
const G4int n
virtual G4double FindValue(G4double e, G4int componentId=0) const
virtual size_t NumberOfComponents(void) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4DNAChemistryManager * Instance()
static G4DNAMolecularMaterial * Instance()
G4double GetPDGMass() const
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
Hep3Vector unit() const
tuple t1
Definition: plottest35.py:33
double y() const
G4ParticleChangeForGamma * fParticleChangeForGamma
double DifferentialCrossSection(G4ParticleDefinition *aParticleDefinition, G4double k, G4double energyTransfer, G4int shell)
const G4Track * GetCurrentTrack() const
const XML_Char int const XML_Char * value
static G4Electron * Electron()
Definition: G4Electron.cc:94
void SetProposedKineticEnergy(G4double proposedKinEnergy)
#define G4endl
Definition: G4ios.hh:61
G4VAtomDeexcitation * AtomDeexcitation()
double G4double
Definition: G4Types.hh:76
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:690
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:739
G4double bindingEnergy(G4int A, G4int Z)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &=*(new G4DataVector()))
G4AtomicShellEnumerator
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:121