source: trunk/source/processes/electromagnetic/lowenergy/src/G4DNADingfelderChargeIncreaseModel.cc @ 1347

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago

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