Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros 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 96060 2016-03-11 12:58:04Z gcosmo $
27 //
28 
30 #include "G4PhysicalConstants.hh"
31 #include "G4SystemOfUnits.hh"
33 #include "G4DNAChemistryManager.hh"
34 
35 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
36 
37 using namespace std;
38 
39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
40 
42  const G4String& nam) :
43  G4VEmModel(nam), isInitialised(false)
44 {
45  // nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
46  fpMolWaterDensity = 0;
47  numberOfPartialCrossSections[0] = 0;
48  numberOfPartialCrossSections[1] = 0;
49  numberOfPartialCrossSections[2] = 0;
50 
51  verboseLevel = 0;
52  // Verbosity scale:
53  // 0 = nothing
54  // 1 = warning for energy non-conservation
55  // 2 = details of energy budget
56  // 3 = calculation of cross sections, file openings, sampling of atoms
57  // 4 = entering in methods
58 
59  if (verboseLevel > 0)
60  {
61  G4cout << "Dingfelder charge decrease model is constructed " << G4endl;
62  }
64 
65  // Selection of stationary mode
66 
67  statCode = false;
68 }
69 
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71 
73 {
74 }
75 
76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77 
79  const G4DataVector& /*cuts*/)
80 {
81 
82  if (verboseLevel > 3)
83  {
84  G4cout << "Calling G4DNADingfelderChargeDecreaseModel::Initialise()"
85  << G4endl;
86  }
87 
88  // Energy limits
89 
90  G4DNAGenericIonsManager *instance;
93  G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
94  G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
95 
97  G4String alphaPlusPlus;
98  G4String alphaPlus;
99 
100  // LIMITS
101 
102  proton = protonDef->GetParticleName();
103  lowEnergyLimit[proton] = 100. * eV;
104  highEnergyLimit[proton] = 100. * MeV;
105 
106  alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
107  lowEnergyLimit[alphaPlusPlus] = 1. * keV;
108  highEnergyLimit[alphaPlusPlus] = 400. * MeV;
109 
110  alphaPlus = alphaPlusDef->GetParticleName();
111  lowEnergyLimit[alphaPlus] = 1. * keV;
112  highEnergyLimit[alphaPlus] = 400. * MeV;
113 
114  //
115 
116  if (particle==protonDef)
117  {
118  SetLowEnergyLimit(lowEnergyLimit[proton]);
119  SetHighEnergyLimit(highEnergyLimit[proton]);
120  }
121 
122  if (particle==alphaPlusPlusDef)
123  {
124  SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
125  SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
126  }
127 
128  if (particle==alphaPlusDef)
129  {
130  SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
131  SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
132  }
133 
134  // Final state
135 
136  //PROTON
137  f0[0][0]=1.;
138  a0[0][0]=-0.180;
139  a1[0][0]=-3.600;
140  b0[0][0]=-18.22;
141  b1[0][0]=-1.997;
142  c0[0][0]=0.215;
143  d0[0][0]=3.550;
144  x0[0][0]=3.450;
145  x1[0][0]=5.251;
146 
147  numberOfPartialCrossSections[0] = 1;
148 
149  //ALPHA++
150  f0[0][1]=1.; a0[0][1]=0.95;
151  a1[0][1]=-2.75;
152  b0[0][1]=-23.00;
153  c0[0][1]=0.215;
154  d0[0][1]=2.95;
155  x0[0][1]=3.50;
156 
157  f0[1][1]=1.;
158  a0[1][1]=0.95;
159  a1[1][1]=-2.75;
160  b0[1][1]=-23.73;
161  c0[1][1]=0.250;
162  d0[1][1]=3.55;
163  x0[1][1]=3.72;
164 
165  x1[0][1]=-1.;
166  b1[0][1]=-1.;
167 
168  x1[1][1]=-1.;
169  b1[1][1]=-1.;
170 
171  numberOfPartialCrossSections[1] = 2;
172 
173  // ALPHA+
174  f0[0][2]=1.;
175  a0[0][2]=0.65;
176  a1[0][2]=-2.75;
177  b0[0][2]=-21.81;
178  c0[0][2]=0.232;
179  d0[0][2]=2.95;
180  x0[0][2]=3.53;
181 
182  x1[0][2]=-1.;
183  b1[0][2]=-1.;
184 
185  numberOfPartialCrossSections[2] = 1;
186 
187  //
188 
189  if( verboseLevel>0 )
190  {
191  G4cout << "Dingfelder charge decrease model is initialized " << G4endl
192  << "Energy range: "
193  << LowEnergyLimit() / keV << " keV - "
194  << HighEnergyLimit() / MeV << " MeV for "
195  << particle->GetParticleName()
196  << G4endl;
197  }
198 
199  // Initialize water density pointer
201 
202  if (isInitialised)
203  { return;}
205  isInitialised = true;
206 
207 }
208 
209 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
210 
212  const G4ParticleDefinition* particleDefinition,
213  G4double k,
214  G4double,
215  G4double)
216 {
217  if (verboseLevel > 3)
218  {
219  G4cout
220  << "Calling CrossSectionPerVolume() of G4DNADingfelderChargeDecreaseModel"
221  << G4endl;
222  }
223 
224  // Calculate total cross section for model
225 
226  G4DNAGenericIonsManager *instance;
228 
229  if (
230  particleDefinition != G4Proton::ProtonDefinition()
231  &&
232  particleDefinition != instance->GetIon("alpha++")
233  &&
234  particleDefinition != instance->GetIon("alpha+")
235  )
236 
237  return 0;
238 
239  G4double lowLim = 0;
240  G4double highLim = 0;
241  G4double crossSection = 0.;
242 
243  G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
244 
245  if(waterDensity!= 0.0)
246  // if (material == nistwater || material->GetBaseMaterial() == nistwater)
247  {
248  const G4String& particleName = particleDefinition->GetParticleName();
249 
250  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
251  pos1 = lowEnergyLimit.find(particleName);
252 
253  if (pos1 != lowEnergyLimit.end())
254  {
255  lowLim = pos1->second;
256  }
257 
258  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
259  pos2 = highEnergyLimit.find(particleName);
260 
261  if (pos2 != highEnergyLimit.end())
262  {
263  highLim = pos2->second;
264  }
265 
266  if (k >= lowLim && k < highLim)
267  {
268  crossSection = Sum(k,particleDefinition);
269  }
270 
271  if (verboseLevel > 2)
272  {
273  G4cout << "_______________________________________" << G4endl;
274  G4cout << "G4DNADingfelderChargeDecreaeModel" << G4endl;
275  G4cout << "Kinetic energy(eV)=" << k/eV << "particle :" << particleName << G4endl;
276  G4cout << "Cross section per water molecule (cm^2)=" << crossSection/cm/cm << G4endl;
277  G4cout << "Cross section per water molecule (cm^-1)=" << crossSection*
278  waterDensity/(1./cm) << G4endl;
279 // material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
280  }
281  }
282 
283  return crossSection*waterDensity;
284 // return crossSection*material->GetAtomicNumDensityVector()[1];
285 
286 }
287 
288 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
289 
291  G4DynamicParticle*>* fvect,
292  const G4MaterialCutsCouple* /*couple*/,
293  const G4DynamicParticle* aDynamicParticle,
294  G4double,
295  G4double)
296 {
297  if (verboseLevel > 3)
298  {
299  G4cout
300  << "Calling SampleSecondaries() of G4DNADingfelderChargeDecreaseModel"
301  << G4endl;
302  }
303 
304  G4double inK = aDynamicParticle->GetKineticEnergy();
305 
306  G4ParticleDefinition* definition = aDynamicParticle->GetDefinition();
307 
308  G4double particleMass = definition->GetPDGMass();
309 
310  G4int finalStateIndex = RandomSelect(inK,definition);
311 
312  G4int n = NumberOfFinalStates(definition, finalStateIndex);
313  G4double waterBindingEnergy = WaterBindingEnergyConstant(definition, finalStateIndex);
314  G4double outgoingParticleBindingEnergy = OutgoingParticleBindingEnergyConstant(definition, finalStateIndex);
315 
316  G4double outK = 0.;
317 
318  if (definition==G4Proton::Proton())
319  {
320  if (!statCode) outK = inK - n*(inK*electron_mass_c2/proton_mass_c2)
321  - waterBindingEnergy + outgoingParticleBindingEnergy;
322  else outK = inK;
323  }
324 
325  else
326  {
327  if (!statCode) outK = inK - n*(inK*electron_mass_c2/particleMass)
328  - waterBindingEnergy + outgoingParticleBindingEnergy;
329  else outK = inK;
330  }
331 
332  if (outK<0)
333  {
334  G4Exception("G4DNADingfelderChargeDecreaseModel::SampleSecondaries","em0004",
335  FatalException,"Final kinetic energy is negative.");
336  }
337 
339 
340  if (!statCode) fParticleChangeForGamma->ProposeLocalEnergyDeposit(waterBindingEnergy);
341 
342  else
343 
344  {
345  if (definition==G4Proton::Proton())
347  + waterBindingEnergy - outgoingParticleBindingEnergy);
348  else
350  + waterBindingEnergy - outgoingParticleBindingEnergy);
351  }
352 
353  G4DynamicParticle* dp = new G4DynamicParticle (OutgoingParticleDefinition(definition, finalStateIndex),
354  aDynamicParticle->GetMomentumDirection(),
355  outK);
356  fvect->push_back(dp);
357 
358  const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
360  1,
361  theIncomingTrack);
362 
363 }
364 
365 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
366 
367 G4int G4DNADingfelderChargeDecreaseModel::NumberOfFinalStates(G4ParticleDefinition* particleDefinition,
368  G4int finalStateIndex)
369 
370 {
371  if (particleDefinition == G4Proton::Proton())
372  return 1;
373 
374  G4DNAGenericIonsManager*instance;
376 
377  if (particleDefinition == instance->GetIon("alpha++"))
378  {
379  if (finalStateIndex == 0)
380  return 1;
381  return 2;
382  }
383 
384  if (particleDefinition == instance->GetIon("alpha+"))
385  return 1;
386 
387  return 0;
388 }
389 
390 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
391 
392 G4ParticleDefinition* G4DNADingfelderChargeDecreaseModel::OutgoingParticleDefinition(G4ParticleDefinition* particleDefinition,
393  G4int finalStateIndex)
394 {
396 
397  if (particleDefinition == G4Proton::Proton())
398  return instance->GetIon("hydrogen");
399 
400  if (particleDefinition == instance->GetIon("alpha++"))
401  {
402  if (finalStateIndex == 0)
403  return instance->GetIon("alpha+");
404  return instance->GetIon("helium");
405  }
406 
407  if (particleDefinition == instance->GetIon("alpha+"))
408  return instance->GetIon("helium");
409 
410  return 0;
411 }
412 
413 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
414 
415 G4double G4DNADingfelderChargeDecreaseModel::WaterBindingEnergyConstant(G4ParticleDefinition* particleDefinition,
416  G4int finalStateIndex)
417 {
418  // Ionization energy of first water shell
419  // Rad. Phys. Chem. 59 p.267 by Dingf. et al.
420  // W + 10.79 eV -> W+ + e-
421 
423 
424  if (particleDefinition == G4Proton::Proton())
425  return 10.79 * eV;
426 
427  if (particleDefinition == instance->GetIon("alpha++"))
428  {
429  // Binding energy for W+ -> W++ + e- 10.79 eV
430  // Binding energy for W -> W+ + e- 10.79 eV
431  //
432  // Ionization energy of first water shell
433  // Rad. Phys. Chem. 59 p.267 by Dingf. et al.
434 
435  if (finalStateIndex == 0)
436  return 10.79 * eV;
437 
438  return 10.79 * 2 * eV;
439  }
440 
441  if (particleDefinition == instance->GetIon("alpha+"))
442  {
443  // Binding energy for W+ -> W++ + e- 10.79 eV
444  // Binding energy for W -> W+ + e- 10.79 eV
445  //
446  // Ionization energy of first water shell
447  // Rad. Phys. Chem. 59 p.267 by Dingf. et al.
448 
449  return 10.79 * eV;
450  }
451 
452  return 0;
453 }
454 
455 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
456 
457 G4double G4DNADingfelderChargeDecreaseModel::OutgoingParticleBindingEnergyConstant(G4ParticleDefinition* particleDefinition,
458  G4int finalStateIndex)
459 {
461 
462  if (particleDefinition == G4Proton::Proton())
463  return 13.6 * eV;
464 
465  if (particleDefinition == instance->GetIon("alpha++"))
466  {
467  // Binding energy for He+ -> He++ + e- 54.509 eV
468  // Binding energy for He -> He+ + e- 24.587 eV
469 
470  if (finalStateIndex == 0)
471  return 54.509 * eV;
472 
473  return (54.509 + 24.587) * eV;
474  }
475 
476  if (particleDefinition == instance->GetIon("alpha+"))
477  {
478  // Binding energy for He+ -> He++ + e- 54.509 eV
479  // Binding energy for He -> He+ + e- 24.587 eV
480 
481  return 24.587 * eV;
482  }
483 
484  return 0;
485 }
486 
487 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
488 
489 G4double G4DNADingfelderChargeDecreaseModel::PartialCrossSection(G4double k,
490  G4int index,
491  const G4ParticleDefinition* particleDefinition)
492 {
493  G4int particleTypeIndex = 0;
494  G4DNAGenericIonsManager* instance;
496 
497  if (particleDefinition == G4Proton::ProtonDefinition())
498  particleTypeIndex = 0;
499  if (particleDefinition == instance->GetIon("alpha++"))
500  particleTypeIndex = 1;
501  if (particleDefinition == instance->GetIon("alpha+"))
502  particleTypeIndex = 2;
503 
504  //
505  // sigma(T) = f0 10 ^ y(log10(T/eV))
506  //
507  // / a0 x + b0 x < x0
508  // |
509  // y(x) = < a0 x + b0 - c0 (x - x0)^d0 x0 <= x < x1
510  // |
511  // \ a1 x + b1 x >= x1
512  //
513  //
514  // f0, a0, a1, b0, b1, c0, d0, x0, x1 are parameters that change for protons and helium (0, +, ++)
515  //
516  // 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)
517  //
518  // From Rad. Phys. and Chem. 59 (2000) 255-275, M. Dingfelder et al.
519  // Inelastic-collision cross sections of liquid water for interactions of energetic proton
520  //
521 
522  if (x1[index][particleTypeIndex] < x0[index][particleTypeIndex])
523  {
524  //
525  // 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)
526  //
527  // x1 = x0 + ((a0 - a1)/(c0 * d0)) ^ (1 / (d0 - 1))
528  //
529  // b1 = (a0 - a1) * x1 + b0 - c0 * (x1 - x0) ^ d0
530  //
531 
532  x1[index][particleTypeIndex] = x0[index][particleTypeIndex]
533  + std::pow((a0[index][particleTypeIndex] - a1[index][particleTypeIndex])
534  / (c0[index][particleTypeIndex]
535  * d0[index][particleTypeIndex]),
536  1. / (d0[index][particleTypeIndex] - 1.));
537  b1[index][particleTypeIndex] = (a0[index][particleTypeIndex]
538  - a1[index][particleTypeIndex]) * x1[index][particleTypeIndex]
539  + b0[index][particleTypeIndex]
540  - c0[index][particleTypeIndex]
541  * std::pow(x1[index][particleTypeIndex]
542  - x0[index][particleTypeIndex],
543  d0[index][particleTypeIndex]);
544  }
545 
546  G4double x(std::log10(k / eV));
547  G4double y;
548 
549  if (x < x0[index][particleTypeIndex])
550  y = a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex];
551  else if (x < x1[index][particleTypeIndex])
552  y = a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex]
553  - c0[index][particleTypeIndex]
554  * std::pow(x - x0[index][particleTypeIndex],
555  d0[index][particleTypeIndex]);
556  else
557  y = a1[index][particleTypeIndex] * x + b1[index][particleTypeIndex];
558 
559  return f0[index][particleTypeIndex] * std::pow(10., y) * m * m;
560 
561 }
562 
563 G4int G4DNADingfelderChargeDecreaseModel::RandomSelect(G4double k,
564  const G4ParticleDefinition* particleDefinition)
565 {
566  G4int particleTypeIndex = 0;
567  G4DNAGenericIonsManager *instance;
569 
570  if (particleDefinition == G4Proton::ProtonDefinition())
571  particleTypeIndex = 0;
572  if (particleDefinition == instance->GetIon("alpha++"))
573  particleTypeIndex = 1;
574  if (particleDefinition == instance->GetIon("alpha+"))
575  particleTypeIndex = 2;
576 
577  const G4int n = numberOfPartialCrossSections[particleTypeIndex];
578  G4double* values(new G4double[n]);
579  G4double value(0);
580  G4int i = n;
581 
582  while (i > 0)
583  {
584  i--;
585  values[i] = PartialCrossSection(k, i, particleDefinition);
586  value += values[i];
587  }
588 
589  value *= G4UniformRand();
590 
591  i = n;
592  while (i > 0)
593  {
594  i--;
595 
596  if (values[i] > value)
597  break;
598 
599  value -= values[i];
600  }
601 
602  delete[] values;
603 
604  return i;
605 }
606 
607 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
608 
609 G4double G4DNADingfelderChargeDecreaseModel::Sum(G4double k,
610  const G4ParticleDefinition* particleDefinition)
611 {
612  G4int particleTypeIndex = 0;
613  G4DNAGenericIonsManager* instance;
615 
616  if (particleDefinition == G4Proton::ProtonDefinition())
617  particleTypeIndex = 0;
618  if (particleDefinition == instance->GetIon("alpha++"))
619  particleTypeIndex = 1;
620  if (particleDefinition == instance->GetIon("alpha+"))
621  particleTypeIndex = 2;
622 
623  G4double totalCrossSection = 0.;
624 
625  for (G4int i = 0; i < numberOfPartialCrossSections[particleTypeIndex]; i++)
626  {
627  totalCrossSection += PartialCrossSection(k, i, particleDefinition);
628  }
629  return totalCrossSection;
630 }
631 
G4DNADingfelderChargeDecreaseModel(const G4ParticleDefinition *p=0, const G4String &nam="DNADingfelderChargeDecreaseModel")
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:640
static constexpr double proton_mass_c2
G4double GetKineticEnergy() const
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:633
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:602
size_t GetIndex() const
Definition: G4Material.hh:262
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
G4ParticleDefinition * GetDefinition() const
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
static constexpr double electron_mass_c2
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:724
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
const std::vector< double > * GetNumMolPerVolTableFor(const G4Material *) const
static constexpr double m
Definition: G4SIunits.hh:129
const XML_Char int const XML_Char * value
Definition: expat.h:331
const G4ThreeVector & GetMomentumDirection() const
static constexpr double cm
Definition: G4SIunits.hh:119
static G4Proton * Proton()
Definition: G4Proton.cc:93
static constexpr double eV
Definition: G4SIunits.hh:215
static G4DNAGenericIonsManager * Instance(void)
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 G4DNAChemistryManager * Instance()
static G4DNAMolecularMaterial * Instance()
G4double GetPDGMass() const
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
const G4Track * GetCurrentTrack() const
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:731
static constexpr double keV
Definition: G4SIunits.hh:216
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:132
G4ParticleDefinition * GetIon(const G4String &name)