ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNADingfelderChargeIncreaseModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DNADingfelderChargeIncreaseModel.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 //
27 
28 #include "G4DNADingfelderChargeIncreaseModel.hh"
29 #include "G4PhysicalConstants.hh"
30 #include "G4SystemOfUnits.hh"
31 #include "G4DNAMolecularMaterial.hh"
32 
33 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
34 
35 using namespace std;
36 
37 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
38 
39 G4DNADingfelderChargeIncreaseModel::G4DNADingfelderChargeIncreaseModel(const G4ParticleDefinition*,
40  const G4String& nam) :
41 G4VEmModel(nam), isInitialised(false)
42 {
43  fpMolWaterDensity = 0;
44 
45  numberOfPartialCrossSections[0] = 0;
46  numberOfPartialCrossSections[1] = 0;
47 
48  verboseLevel = 0;
49  // Verbosity scale:
50  // 0 = nothing
51  // 1 = warning for energy non-conservation
52  // 2 = details of energy budget
53  // 3 = calculation of cross sections, file openings, sampling of atoms
54  // 4 = entering in methods
55 
56  if (verboseLevel > 0)
57  {
58  G4cout << "Dingfelder charge increase model is constructed " << G4endl;
59  }
60  fParticleChangeForGamma = 0;
61 
62  // Selection of stationary mode
63 
64  statCode = false;
65 }
66 
67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
68 
69 G4DNADingfelderChargeIncreaseModel::~G4DNADingfelderChargeIncreaseModel()
70 {}
71 
72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73 
74 void G4DNADingfelderChargeIncreaseModel::Initialise(const G4ParticleDefinition* particle,
75  const G4DataVector& /*cuts*/)
76 {
77 
78  if (verboseLevel > 3)
79  {
80  G4cout << "Calling G4DNADingfelderChargeIncreaseModel::Initialise()"
81  << G4endl;
82  }
83 
84  // Energy limits
85 
86  G4DNAGenericIonsManager *instance;
87  instance = G4DNAGenericIonsManager::Instance();
88  G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
89  G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
90  G4ParticleDefinition* heliumDef = instance->GetIon("helium");
91 
92  G4String hydrogen;
93  G4String alphaPlus;
94  G4String helium;
95 
96  // Limits
97 
98  hydrogen = hydrogenDef->GetParticleName();
99  lowEnergyLimit[hydrogen] = 100. * eV;
100  highEnergyLimit[hydrogen] = 100. * MeV;
101 
102  alphaPlus = alphaPlusDef->GetParticleName();
103  lowEnergyLimit[alphaPlus] = 1. * keV;
104  highEnergyLimit[alphaPlus] = 400. * MeV;
105 
106  helium = heliumDef->GetParticleName();
107  lowEnergyLimit[helium] = 1. * keV;
108  highEnergyLimit[helium] = 400. * MeV;
109 
110  //
111 
112  if (particle==hydrogenDef)
113  {
114  SetLowEnergyLimit(lowEnergyLimit[hydrogen]);
115  SetHighEnergyLimit(highEnergyLimit[hydrogen]);
116  }
117 
118  if (particle==alphaPlusDef)
119  {
120  SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
121  SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
122  }
123 
124  if (particle==heliumDef)
125  {
126  SetLowEnergyLimit(lowEnergyLimit[helium]);
127  SetHighEnergyLimit(highEnergyLimit[helium]);
128  }
129 
130  // Final state
131 
132  //ALPHA+
133 
134  f0[0][0]=1.;
135  a0[0][0]=2.25;
136  a1[0][0]=-0.75;
137  b0[0][0]=-32.10;
138  c0[0][0]=0.600;
139  d0[0][0]=2.40;
140  x0[0][0]=4.60;
141 
142  x1[0][0]=-1.;
143  b1[0][0]=-1.;
144 
145  numberOfPartialCrossSections[0]=1;
146 
147  //HELIUM
148 
149  f0[0][1]=1.;
150  a0[0][1]=2.25;
151  a1[0][1]=-0.75;
152  b0[0][1]=-30.93;
153  c0[0][1]=0.590;
154  d0[0][1]=2.35;
155  x0[0][1]=4.29;
156 
157  f0[1][1]=1.;
158  a0[1][1]=2.25;
159  a1[1][1]=-0.75;
160  b0[1][1]=-32.61;
161  c0[1][1]=0.435;
162  d0[1][1]=2.70;
163  x0[1][1]=4.45;
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  //
174 
175  if( verboseLevel>0 )
176  {
177  G4cout << "Dingfelder charge increase model is initialized " << G4endl
178  << "Energy range: "
179  << LowEnergyLimit() / keV << " keV - "
180  << HighEnergyLimit() / MeV << " MeV for "
181  << particle->GetParticleName()
182  << G4endl;
183  }
184 
185  // Initialize water density pointer
186  fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
187 
188  if (isInitialised) return;
189 
190  fParticleChangeForGamma = GetParticleChangeForGamma();
191  isInitialised = true;
192 }
193 
194 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
195 
196 G4double G4DNADingfelderChargeIncreaseModel::CrossSectionPerVolume(const G4Material* material,
197  const G4ParticleDefinition* particleDefinition,
198  G4double k,
199  G4double,
200  G4double)
201 {
202  if (verboseLevel > 3)
203  {
204  G4cout
205  << "Calling CrossSectionPerVolume() of G4DNADingfelderChargeIncreaseModel"
206  << G4endl;
207  }
208 
209  // Calculate total cross section for model
210 
211  G4DNAGenericIonsManager *instance;
212  instance = G4DNAGenericIonsManager::Instance();
213 
214  if (
215  particleDefinition != instance->GetIon("hydrogen")
216  &&
217  particleDefinition != instance->GetIon("alpha+")
218  &&
219  particleDefinition != instance->GetIon("helium")
220  )
221 
222  return 0;
223 
224  G4double lowLim = 0;
225  G4double highLim = 0;
226  G4double totalCrossSection = 0.;
227 
228  G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
229 
230  const G4String& particleName = particleDefinition->GetParticleName();
231 
232  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
233  pos1 = lowEnergyLimit.find(particleName);
234 
235  if (pos1 != lowEnergyLimit.end())
236  {
237  lowLim = pos1->second;
238  }
239 
240  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
241  pos2 = highEnergyLimit.find(particleName);
242 
243  if (pos2 != highEnergyLimit.end())
244  {
245  highLim = pos2->second;
246  }
247 
248  if (k >= lowLim && k <= highLim)
249  {
250  //HYDROGEN
251  if (particleDefinition == instance->GetIon("hydrogen"))
252  {
253  const G4double aa = 2.835;
254  const G4double bb = 0.310;
255  const G4double cc = 2.100;
256  const G4double dd = 0.760;
257  const G4double fac = 1.0e-18;
258  const G4double rr = 13.606 * eV;
259 
260  G4double t = k / (proton_mass_c2/electron_mass_c2);
261  G4double x = t / rr;
262  G4double temp = 4.0 * pi * Bohr_radius/nm * Bohr_radius/nm * fac;
263  G4double sigmal = temp * cc * (std::pow(x,dd));
264  G4double sigmah = temp * (aa * std::log(1.0 + x) + bb) / x;
265  totalCrossSection = 1.0/(1.0/sigmal + 1.0/sigmah) *m*m;
266  }
267  else
268  {
269  totalCrossSection = Sum(k,particleDefinition);
270  }
271  }
272 
273  if (verboseLevel > 2)
274  {
275  G4cout << "__________________________________" << G4endl;
276  G4cout << "G4DNADingfelderChargeIncreaseModel - XS INFO START" << G4endl;
277  G4cout << "Kinetic energy(eV)=" << k/eV << " particle : " << particleName << G4endl;
278  G4cout << "Cross section per water molecule (cm^2)=" << totalCrossSection/cm/cm << G4endl;
279  G4cout << "Cross section per water molecule (cm^-1)=" << totalCrossSection*waterDensity/(1./cm) << G4endl;
280  // G4cout << " - Cross section per water molecule (cm^-1)="
281  // << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
282  G4cout << "G4DNADingfelderChargeIncreaseModel - XS INFO END" << G4endl;
283  }
284 
285  return totalCrossSection*waterDensity;
286  // return totalCrossSection*material->GetAtomicNumDensityVector()[1];
287 
288 }
289 
290 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
291 
292 void G4DNADingfelderChargeIncreaseModel::SampleSecondaries(std::vector<
293  G4DynamicParticle*>* fvect,
294  const G4MaterialCutsCouple* /*couple*/,
295  const G4DynamicParticle* aDynamicParticle,
296  G4double,
297  G4double)
298 {
299  if (verboseLevel > 3)
300  {
301  G4cout
302  << "Calling SampleSecondaries() of G4DNADingfelderChargeIncreaseModel"
303  << G4endl;
304  }
305 
306  if (!statCode) fParticleChangeForGamma->ProposeLocalEnergyDeposit(0.);
307 
308  G4ParticleDefinition* definition = aDynamicParticle->GetDefinition();
309 
310  G4double particleMass = definition->GetPDGMass();
311 
312  G4double inK = aDynamicParticle->GetKineticEnergy();
313 
314  G4int finalStateIndex = RandomSelect(inK,definition);
315 
316  G4int n = NumberOfFinalStates(definition,finalStateIndex);
317 
318  G4double outK = 0.;
319 
320  if (!statCode) outK = inK - IncomingParticleBindingEnergyConstant(definition,finalStateIndex);
321 
322  else outK = inK;
323 
324  if (statCode)
325  fParticleChangeForGamma->
326  ProposeLocalEnergyDeposit(IncomingParticleBindingEnergyConstant(definition,finalStateIndex));
327 
328  fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
329 
330  G4DNAGenericIonsManager* instance;
331  instance = G4DNAGenericIonsManager::Instance();
332 
333  G4double electronK;
334  if (definition == instance->GetIon("hydrogen")) electronK = inK*electron_mass_c2/proton_mass_c2;
335  else electronK = inK*electron_mass_c2/(particleMass);
336 
337  if (outK<0)
338  {
339  G4Exception("G4DNADingfelderChargeIncreaseModel::SampleSecondaries","em0004",
340  FatalException,"Final kinetic energy is negative.");
341  }
342 
343  G4DynamicParticle* dp = new G4DynamicParticle(OutgoingParticleDefinition(definition,finalStateIndex),
344  aDynamicParticle->GetMomentumDirection(),
345  outK);
346 
347  fvect->push_back(dp);
348 
349  n = n - 1;
350 
351  while (n>0)
352  {
353  n--;
354  fvect->push_back(new G4DynamicParticle
355  (G4Electron::Electron(), aDynamicParticle->GetMomentumDirection(), electronK) );
356  }
357 }
358 
359 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
360 
361 G4int G4DNADingfelderChargeIncreaseModel::NumberOfFinalStates(G4ParticleDefinition* particleDefinition,
362  G4int finalStateIndex)
363 
364 {
365  G4DNAGenericIonsManager* instance;
366  instance = G4DNAGenericIonsManager::Instance();
367 
368  if (particleDefinition == instance->GetIon("hydrogen"))
369  return 2;
370 
371  if (particleDefinition == instance->GetIon("alpha+"))
372  return 2;
373 
374  if (particleDefinition == instance->GetIon("helium"))
375  {
376  if (finalStateIndex == 0)
377  return 2;
378  return 3;
379  }
380 
381  return 0;
382 }
383 
384 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
385 
386 G4ParticleDefinition* G4DNADingfelderChargeIncreaseModel::OutgoingParticleDefinition(G4ParticleDefinition* particleDefinition,
387  G4int finalStateIndex)
388 {
389  G4DNAGenericIonsManager * instance(G4DNAGenericIonsManager::Instance());
390 
391  if (particleDefinition == instance->GetIon("hydrogen"))
392  return G4Proton::Proton();
393 
394  if (particleDefinition == instance->GetIon("alpha+"))
395  return instance->GetIon("alpha++");
396 
397  if (particleDefinition == instance->GetIon("helium"))
398  {
399  if (finalStateIndex == 0)
400  return instance->GetIon("alpha+");
401  return instance->GetIon("alpha++");
402  }
403 
404  return 0;
405 }
406 
407 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
408 
409 G4double G4DNADingfelderChargeIncreaseModel::IncomingParticleBindingEnergyConstant(G4ParticleDefinition* particleDefinition,
410  G4int finalStateIndex)
411 {
412  G4DNAGenericIonsManager * instance(G4DNAGenericIonsManager::Instance());
413 
414  if (particleDefinition == instance->GetIon("hydrogen"))
415  return 13.6 * eV;
416 
417  if (particleDefinition == instance->GetIon("alpha+"))
418  {
419  // Binding energy for He+ -> He++ + e- 54.509 eV
420  // Binding energy for He -> He+ + e- 24.587 eV
421  return 54.509 * eV;
422  }
423 
424  if (particleDefinition == instance->GetIon("helium"))
425  {
426  // Binding energy for He+ -> He++ + e- 54.509 eV
427  // Binding energy for He -> He+ + e- 24.587 eV
428 
429  if (finalStateIndex == 0)
430  return 24.587 * eV;
431  return (54.509 + 24.587) * eV;
432  }
433 
434  return 0;
435 }
436 
437 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
438 
439 G4double G4DNADingfelderChargeIncreaseModel::PartialCrossSection(G4double k,
440  G4int index,
441  const G4ParticleDefinition* particleDefinition)
442 {
443  G4int particleTypeIndex = 0;
444  G4DNAGenericIonsManager *instance;
445  instance = G4DNAGenericIonsManager::Instance();
446 
447  if (particleDefinition == instance->GetIon("alpha+"))
448  particleTypeIndex = 0;
449 
450  if (particleDefinition == instance->GetIon("helium"))
451  particleTypeIndex = 1;
452 
453  //
454  // sigma(T) = f0 10 ^ y(log10(T/eV))
455  //
456  // / a0 x + b0 x < x0
457  // |
458  // y(x) = < a0 x + b0 - c0 (x - x0)^d0 x0 <= x < x1
459  // |
460  // \ a1 x + b1 x >= x1
461  //
462  //
463  // f0, a0, a1, b0, b1, c0, d0, x0, x1 are parameters that change for protons and helium (0, +, ++)
464  //
465  // f0 has been added to the code in order to manage partial (shell-dependent) cross sections
466  // (if no shell dependence is present. f0=1. Sum of f0 over the considered shells should give 1)
467  //
468  // From Rad. Phys. and Chem. 59 (2000) 255-275, M. Dingfelder et al.
469  // Inelastic-collision cross sections of liquid water for interactions of energetic proton
470  //
471 
472  if (x1[index][particleTypeIndex] < x0[index][particleTypeIndex])
473  {
474  //
475  // if x1 < x0 means that x1 and b1 will be calculated with the following formula
476  // (this piece of code is run on all alphas and not on protons)
477  //
478  // x1 = x0 + ((a0 - a1)/(c0 * d0)) ^ (1 / (d0 - 1))
479  //
480  // b1 = (a0 - a1) * x1 + b0 - c0 * (x1 - x0) ^ d0
481  //
482 
483  x1[index][particleTypeIndex] = x0[index][particleTypeIndex]
484  + std::pow((a0[index][particleTypeIndex] - a1[index][particleTypeIndex])
485  / (c0[index][particleTypeIndex]
486  * d0[index][particleTypeIndex]),
487  1. / (d0[index][particleTypeIndex] - 1.));
488  b1[index][particleTypeIndex] = (a0[index][particleTypeIndex]
489  - a1[index][particleTypeIndex]) * x1[index][particleTypeIndex]
490  + b0[index][particleTypeIndex]
491  - c0[index][particleTypeIndex]
492  * std::pow(x1[index][particleTypeIndex]
493  - x0[index][particleTypeIndex],
494  d0[index][particleTypeIndex]);
495  }
496 
497  G4double x(std::log10(k / eV));
498  G4double y;
499 
500  if (x < x0[index][particleTypeIndex])
501  y = a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex];
502  else if (x < x1[index][particleTypeIndex])
503  y = a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex]
504  - c0[index][particleTypeIndex]
505  * std::pow(x - x0[index][particleTypeIndex],
506  d0[index][particleTypeIndex]);
507  else
508  y = a1[index][particleTypeIndex] * x + b1[index][particleTypeIndex];
509 
510  return f0[index][particleTypeIndex] * std::pow(10., y) * m * m;
511 
512 }
513 
514 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
515 
516 G4int G4DNADingfelderChargeIncreaseModel::RandomSelect(G4double k,
517  const G4ParticleDefinition* particleDefinition)
518 {
519  G4int particleTypeIndex = 0;
520  G4DNAGenericIonsManager *instance;
521  instance = G4DNAGenericIonsManager::Instance();
522 
523  if (particleDefinition == instance->GetIon("hydrogen"))
524  return 0;
525 
526  if (particleDefinition == instance->GetIon("alpha+"))
527  particleTypeIndex = 0;
528 
529  if (particleDefinition == instance->GetIon("helium"))
530  particleTypeIndex = 1;
531 
532  const G4int n = numberOfPartialCrossSections[particleTypeIndex];
533  G4double* values(new G4double[n]);
534  G4double value = 0;
535  G4int i = n;
536 
537  while (i > 0)
538  {
539  i--;
540  values[i] = PartialCrossSection(k, i, particleDefinition);
541  value += values[i];
542  }
543 
544  value *= G4UniformRand();
545 
546  i = n;
547  while (i > 0)
548  {
549  i--;
550 
551  if (values[i] > value)
552  break;
553 
554  value -= values[i];
555  }
556 
557  delete[] values;
558 
559  return i;
560 }
561 
562 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
563 
564 G4double G4DNADingfelderChargeIncreaseModel::Sum(G4double k,
565  const G4ParticleDefinition* particleDefinition)
566 {
567  G4int particleTypeIndex = 0;
568  G4DNAGenericIonsManager *instance;
569  instance = G4DNAGenericIonsManager::Instance();
570 
571  if (particleDefinition == instance->GetIon("alpha+"))
572  particleTypeIndex = 0;
573 
574  if (particleDefinition == instance->GetIon("helium"))
575  particleTypeIndex = 1;
576 
577  G4double totalCrossSection = 0.;
578 
579  for (G4int i = 0; i < numberOfPartialCrossSections[particleTypeIndex]; i++)
580  {
581  totalCrossSection += PartialCrossSection(k, i, particleDefinition);
582  }
583 
584  return totalCrossSection;
585 }