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