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