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