Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNADingfelderChargeDecreaseModel.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4DNADingfelderChargeDecreaseModel.cc 70171 2013-05-24 13:34:18Z gcosmo $
27 //
28 
30 #include "G4PhysicalConstants.hh"
31 #include "G4SystemOfUnits.hh"
33 
34 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
35 
36 using namespace std;
37 
38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
39 
41  const G4String& nam)
42  :G4VEmModel(nam),isInitialised(false)
43 {
44  // nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
45  fpMolWaterDensity = 0;
46  numberOfPartialCrossSections[0] = 0;
47  numberOfPartialCrossSections[1] = 0;
48  numberOfPartialCrossSections[2] = 0;
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 << "Dingfelder charge decrease model is constructed " << G4endl;
61  }
63 }
64 
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66 
68 {}
69 
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71 
73  const G4DataVector& /*cuts*/)
74 {
75 
76  if (verboseLevel > 3)
77  G4cout << "Calling G4DNADingfelderChargeDecreaseModel::Initialise()" << G4endl;
78 
79  // Energy limits
80 
81  G4DNAGenericIonsManager *instance;
84  G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
85  G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
86 
88  G4String alphaPlusPlus;
89  G4String alphaPlus;
90 
91  // LIMITS
92 
93  proton = protonDef->GetParticleName();
94  lowEnergyLimit[proton] = 100. * eV;
95  highEnergyLimit[proton] = 100. * MeV;
96 
97  alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
98  lowEnergyLimit[alphaPlusPlus] = 1. * keV;
99  highEnergyLimit[alphaPlusPlus] = 400. * MeV;
100 
101  alphaPlus = alphaPlusDef->GetParticleName();
102  lowEnergyLimit[alphaPlus] = 1. * keV;
103  highEnergyLimit[alphaPlus] = 400. * MeV;
104 
105  //
106 
107  if (particle==protonDef)
108  {
109  SetLowEnergyLimit(lowEnergyLimit[proton]);
110  SetHighEnergyLimit(highEnergyLimit[proton]);
111  }
112 
113  if (particle==alphaPlusPlusDef)
114  {
115  SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
116  SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
117  }
118 
119  if (particle==alphaPlusDef)
120  {
121  SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
122  SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
123  }
124 
125  // Final state
126 
127  //PROTON
128  f0[0][0]=1.;
129  a0[0][0]=-0.180;
130  a1[0][0]=-3.600;
131  b0[0][0]=-18.22;
132  b1[0][0]=-1.997;
133  c0[0][0]=0.215;
134  d0[0][0]=3.550;
135  x0[0][0]=3.450;
136  x1[0][0]=5.251;
137 
138  numberOfPartialCrossSections[0] = 1;
139 
140  //ALPHA++
141  f0[0][1]=1.; a0[0][1]=0.95;
142  a1[0][1]=-2.75;
143  b0[0][1]=-23.00;
144  c0[0][1]=0.215;
145  d0[0][1]=2.95;
146  x0[0][1]=3.50;
147 
148  f0[1][1]=1.;
149  a0[1][1]=0.95;
150  a1[1][1]=-2.75;
151  b0[1][1]=-23.73;
152  c0[1][1]=0.250;
153  d0[1][1]=3.55;
154  x0[1][1]=3.72;
155 
156  x1[0][1]=-1.;
157  b1[0][1]=-1.;
158 
159  x1[1][1]=-1.;
160  b1[1][1]=-1.;
161 
162  numberOfPartialCrossSections[1] = 2;
163 
164  // ALPHA+
165  f0[0][2]=1.;
166  a0[0][2]=0.65;
167  a1[0][2]=-2.75;
168  b0[0][2]=-21.81;
169  c0[0][2]=0.232;
170  d0[0][2]=2.95;
171  x0[0][2]=3.53;
172 
173  x1[0][2]=-1.;
174  b1[0][2]=-1.;
175 
176  numberOfPartialCrossSections[2] = 1;
177 
178  //
179 
180  if( verboseLevel>0 )
181  {
182  G4cout << "Dingfelder charge decrease model is initialized " << G4endl
183  << "Energy range: "
184  << LowEnergyLimit() / keV << " keV - "
185  << HighEnergyLimit() / MeV << " MeV for "
186  << particle->GetParticleName()
187  << G4endl;
188  }
189 
190  // Initialize water density pointer
192 
193  if (isInitialised) { return; }
195  isInitialised = true;
196 
197 }
198 
199 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
200 
202  const G4ParticleDefinition* particleDefinition,
203  G4double k,
204  G4double,
205  G4double)
206 {
207  if (verboseLevel > 3)
208  G4cout << "Calling CrossSectionPerVolume() of G4DNADingfelderChargeDecreaseModel" << G4endl;
209 
210  // Calculate total cross section for model
211 
212  G4DNAGenericIonsManager *instance;
214 
215  if (
216  particleDefinition != G4Proton::ProtonDefinition()
217  &&
218  particleDefinition != instance->GetIon("alpha++")
219  &&
220  particleDefinition != instance->GetIon("alpha+")
221  )
222 
223  return 0;
224 
225  G4double lowLim = 0;
226  G4double highLim = 0;
227  G4double crossSection = 0.;
228 
229  G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
230 
231  if(waterDensity!= 0.0)
232  // if (material == nistwater || material->GetBaseMaterial() == nistwater)
233  {
234  const G4String& particleName = particleDefinition->GetParticleName();
235 
236  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
237  pos1 = lowEnergyLimit.find(particleName);
238 
239  if (pos1 != lowEnergyLimit.end())
240  {
241  lowLim = pos1->second;
242  }
243 
244  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
245  pos2 = highEnergyLimit.find(particleName);
246 
247  if (pos2 != highEnergyLimit.end())
248  {
249  highLim = pos2->second;
250  }
251 
252  if (k >= lowLim && k < highLim)
253  {
254  crossSection = Sum(k,particleDefinition);
255  }
256 
257  if (verboseLevel > 2)
258  {
259  G4cout << "_______________________________________" << G4endl;
260  G4cout << "°°° G4DNADingfelderChargeDecreaeModel" << G4endl;
261  G4cout << "---> Kinetic energy(eV)=" << k/eV << "particle :" << particleName << G4endl;
262  G4cout << " - Cross section per water molecule (cm^2)=" << crossSection/cm/cm << G4endl;
263  G4cout << " - Cross section per water molecule (cm^-1)=" << crossSection*
264  waterDensity/(1./cm) << G4endl;
265 // material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
266  }
267  }
268 
269  return crossSection*waterDensity;
270 // return crossSection*material->GetAtomicNumDensityVector()[1];
271 
272 }
273 
274 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
275 
276 void G4DNADingfelderChargeDecreaseModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
277  const G4MaterialCutsCouple* /*couple*/,
278  const G4DynamicParticle* aDynamicParticle,
279  G4double,
280  G4double)
281 {
282  if (verboseLevel > 3)
283  G4cout << "Calling SampleSecondaries() of G4DNADingfelderChargeDecreaseModel" << G4endl;
284 
285  G4double inK = aDynamicParticle->GetKineticEnergy();
286 
287  G4ParticleDefinition* definition = aDynamicParticle->GetDefinition();
288 
289  G4double particleMass = definition->GetPDGMass();
290 
291  G4int finalStateIndex = RandomSelect(inK,definition);
292 
293  G4int n = NumberOfFinalStates(definition, finalStateIndex);
294  G4double waterBindingEnergy = WaterBindingEnergyConstant(definition, finalStateIndex);
295  G4double outgoingParticleBindingEnergy = OutgoingParticleBindingEnergyConstant(definition, finalStateIndex);
296 
297  G4double outK = 0.;
298  if (definition==G4Proton::Proton())
299  outK = inK - n*(inK*electron_mass_c2/proton_mass_c2) - waterBindingEnergy + outgoingParticleBindingEnergy;
300  else
301  outK = inK - n*(inK*electron_mass_c2/particleMass) - waterBindingEnergy + outgoingParticleBindingEnergy;
302 
303  if (outK<0)
304  {
305  G4Exception("G4DNADingfelderChargeDecreaseModel::SampleSecondaries","em0004",
306  FatalException,"Final kinetic energy is negative.");
307  }
308 
311 
312  G4DynamicParticle* dp = new G4DynamicParticle (OutgoingParticleDefinition(definition, finalStateIndex),
313  aDynamicParticle->GetMomentumDirection(),
314  outK) ;
315  fvect->push_back(dp);
316 
317 }
318 
319 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
320 
321 G4int G4DNADingfelderChargeDecreaseModel::NumberOfFinalStates(G4ParticleDefinition* particleDefinition,
322  G4int finalStateIndex )
323 
324 {
325  if (particleDefinition == G4Proton::Proton()) return 1;
326 
327  G4DNAGenericIonsManager*instance;
329 
330  if (particleDefinition == instance->GetIon("alpha++") )
331  {
332  if (finalStateIndex==0) return 1;
333  return 2;
334  }
335 
336  if (particleDefinition == instance->GetIon("alpha+") ) return 1;
337 
338  return 0;
339 }
340 
341 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
342 
343 G4ParticleDefinition* G4DNADingfelderChargeDecreaseModel::OutgoingParticleDefinition (G4ParticleDefinition* particleDefinition,
344  G4int finalStateIndex)
345 {
347 
348  if (particleDefinition == G4Proton::Proton()) return instance->GetIon("hydrogen");
349 
350  if (particleDefinition == instance->GetIon("alpha++") )
351  {
352  if (finalStateIndex == 0) return instance->GetIon("alpha+");
353  return instance->GetIon("helium");
354  }
355 
356  if (particleDefinition == instance->GetIon("alpha+") ) return instance->GetIon("helium");
357 
358  return 0;
359 }
360 
361 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
362 
363 G4double G4DNADingfelderChargeDecreaseModel::WaterBindingEnergyConstant(G4ParticleDefinition* particleDefinition,
364  G4int finalStateIndex )
365 {
366  // Ionization energy of first water shell
367  // Rad. Phys. Chem. 59 p.267 by Dingf. et al.
368  // W + 10.79 eV -> W+ + e-
369 
371 
372  if(particleDefinition == G4Proton::Proton()) return 10.79*eV;
373 
374  if (particleDefinition == instance->GetIon("alpha++") )
375  {
376  // Binding energy for W+ -> W++ + e- 10.79 eV
377  // Binding energy for W -> W+ + e- 10.79 eV
378  //
379  // Ionization energy of first water shell
380  // Rad. Phys. Chem. 59 p.267 by Dingf. et al.
381 
382  if (finalStateIndex == 0) return 10.79*eV;
383 
384  return 10.79*2*eV;
385  }
386 
387  if (particleDefinition == instance->GetIon("alpha+") )
388  {
389  // Binding energy for W+ -> W++ + e- 10.79 eV
390  // Binding energy for W -> W+ + e- 10.79 eV
391  //
392  // Ionization energy of first water shell
393  // Rad. Phys. Chem. 59 p.267 by Dingf. et al.
394 
395  return 10.79*eV;
396  }
397 
398  return 0;
399 }
400 
401 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
402 
403 G4double G4DNADingfelderChargeDecreaseModel::OutgoingParticleBindingEnergyConstant(G4ParticleDefinition* particleDefinition,
404  G4int finalStateIndex)
405 {
407 
408  if(particleDefinition == G4Proton::Proton()) return 13.6*eV;
409 
410  if (particleDefinition == instance->GetIon("alpha++") )
411  {
412  // Binding energy for He+ -> He++ + e- 54.509 eV
413  // Binding energy for He -> He+ + e- 24.587 eV
414 
415  if (finalStateIndex==0) return 54.509*eV;
416 
417  return (54.509 + 24.587)*eV;
418  }
419 
420  if (particleDefinition == instance->GetIon("alpha+") )
421  {
422  // Binding energy for He+ -> He++ + e- 54.509 eV
423  // Binding energy for He -> He+ + e- 24.587 eV
424 
425  return 24.587*eV;
426  }
427 
428  return 0;
429 }
430 
431 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
432 
433 G4double G4DNADingfelderChargeDecreaseModel::PartialCrossSection(G4double k, G4int index,
434  const G4ParticleDefinition* particleDefinition)
435 {
436  G4int particleTypeIndex = 0;
437  G4DNAGenericIonsManager* instance;
439 
440  if (particleDefinition == G4Proton::ProtonDefinition()) particleTypeIndex=0;
441  if (particleDefinition == instance->GetIon("alpha++")) particleTypeIndex=1;
442  if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=2;
443 
444  //
445  // sigma(T) = f0 10 ^ y(log10(T/eV))
446  //
447  // / a0 x + b0 x < x0
448  // |
449  // y(x) = < a0 x + b0 - c0 (x - x0)^d0 x0 <= x < x1
450  // |
451  // \ a1 x + b1 x >= x1
452  //
453  //
454  // f0, a0, a1, b0, b1, c0, d0, x0, x1 are parameters that change for protons and helium (0, +, ++)
455  //
456  // f0 has been added to the code in order to manage partial (shell-dependent) cross sections (if no shell dependence is present. f0=1. Sum of f0 over the considered shells should give 1)
457  //
458  // From Rad. Phys. and Chem. 59 (2000) 255-275, M. Dingfelder et al.
459  // Inelastic-collision cross sections of liquid water for interactions of energetic proton
460  //
461 
462  if (x1[index][particleTypeIndex]<x0[index][particleTypeIndex])
463  {
464  //
465  // if x1 < x0 means that x1 and b1 will be calculated with the following formula (this piece of code is run on all alphas and not on protons)
466  //
467  // x1 = x0 + ((a0 - a1)/(c0 * d0)) ^ (1 / (d0 - 1))
468  //
469  // b1 = (a0 - a1) * x1 + b0 - c0 * (x1 - x0) ^ d0
470  //
471 
472  x1[index][particleTypeIndex]=x0[index][particleTypeIndex] + std::pow((a0[index][particleTypeIndex] - a1[index][particleTypeIndex])
473  / (c0[index][particleTypeIndex] * d0[index][particleTypeIndex]), 1. / (d0[index][particleTypeIndex] - 1.));
474  b1[index][particleTypeIndex]=(a0[index][particleTypeIndex] - a1[index][particleTypeIndex]) * x1[index][particleTypeIndex]
475  + b0[index][particleTypeIndex] - c0[index][particleTypeIndex] * std::pow(x1[index][particleTypeIndex]
476  - x0[index][particleTypeIndex], d0[index][particleTypeIndex]);
477  }
478 
479  G4double x(std::log10(k/eV));
480  G4double y;
481 
482  if (x<x0[index][particleTypeIndex])
483  y=a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex];
484  else if (x<x1[index][particleTypeIndex])
485  y=a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex] - c0[index][particleTypeIndex]
486  * std::pow(x - x0[index][particleTypeIndex], d0[index][particleTypeIndex]);
487  else
488  y=a1[index][particleTypeIndex] * x + b1[index][particleTypeIndex];
489 
490  return f0[index][particleTypeIndex] * std::pow(10., y)*m*m;
491 
492 }
493 
494 G4int G4DNADingfelderChargeDecreaseModel::RandomSelect(G4double k,
495  const G4ParticleDefinition* particleDefinition)
496 {
497  G4int particleTypeIndex = 0;
498  G4DNAGenericIonsManager *instance;
500 
501  if (particleDefinition == G4Proton::ProtonDefinition()) particleTypeIndex = 0;
502  if (particleDefinition == instance->GetIon("alpha++")) particleTypeIndex = 1;
503  if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex = 2;
504 
505  const G4int n = numberOfPartialCrossSections[particleTypeIndex];
506  G4double* values(new G4double[n]);
507  G4double value(0);
508  G4int i = n;
509 
510  while (i>0)
511  {
512  i--;
513  values[i]=PartialCrossSection(k, i, particleDefinition);
514  value+=values[i];
515  }
516 
517  value*=G4UniformRand();
518 
519  i=n;
520  while (i>0)
521  {
522  i--;
523 
524  if (values[i]>value)
525  break;
526 
527  value-=values[i];
528  }
529 
530  delete[] values;
531 
532  return i;
533 }
534 
535 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
536 
537 G4double G4DNADingfelderChargeDecreaseModel::Sum(G4double k, const G4ParticleDefinition* particleDefinition)
538 {
539  G4int particleTypeIndex = 0;
540  G4DNAGenericIonsManager* instance;
542 
543  if (particleDefinition == G4Proton::ProtonDefinition()) particleTypeIndex=0;
544  if (particleDefinition == instance->GetIon("alpha++")) particleTypeIndex=1;
545  if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=2;
546 
547  G4double totalCrossSection = 0.;
548 
549  for (G4int i=0; i<numberOfPartialCrossSections[particleTypeIndex]; i++)
550  {
551  totalCrossSection += PartialCrossSection(k,i,particleDefinition);
552  }
553  return totalCrossSection;
554 }
555 
G4DNADingfelderChargeDecreaseModel(const G4ParticleDefinition *p=0, const G4String &nam="DNADingfelderChargeDecreaseModel")
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:599
G4double GetKineticEnergy() const
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:592
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:578
size_t GetIndex() const
Definition: G4Material.hh:260
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
G4ParticleDefinition * GetDefinition() const
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:683
string material
Definition: eplot.py:19
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
const std::vector< double > * GetNumMolPerVolTableFor(const G4Material *) const
const G4ThreeVector & GetMomentumDirection() const
static G4Proton * Proton()
Definition: G4Proton.cc:93
float proton_mass_c2
Definition: hepunit.py:275
float electron_mass_c2
Definition: hepunit.py:274
static G4DNAGenericIonsManager * Instance(void)
const G4int n
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4DNAMolecularMaterial * Instance()
G4double GetPDGMass() const
const XML_Char int const XML_Char * value
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:690
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:121
G4ParticleDefinition * GetIon(const G4String &name)