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

Last change on this file since 1250 was 1228, checked in by garnier, 15 years ago

update geant4.9.3 tag
File size: 17.2 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.6 2009/08/13 11:32:47 sincerti Exp $
27// GEANT4 tag $Name: geant4-09-03 $
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] = 1. * keV;
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  // Test if water material
203
204  flagMaterialIsWater= false;
205  densityWater = 0;
206
207  const G4ProductionCutsTable* theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
208
209  if(theCoupleTable) 
210  {
211    G4int numOfCouples = theCoupleTable->GetTableSize();
212 
213    if(numOfCouples>0) 
214    {
215          for (G4int i=0; i<numOfCouples; i++) 
216          {
217            const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
218            const G4Material* material = couple->GetMaterial();
219
220            if (material->GetName() == "G4_WATER") 
221            {
222              G4double density = material->GetAtomicNumDensityVector()[1];
223              flagMaterialIsWater = true; 
224              densityWater = density; 
225             
226              if (verboseLevel > 3) 
227              G4cout << "****** Water material is found with density(cm^-3)=" << density/(cm*cm*cm) << G4endl;
228            }
229 
230          }
231
232    } // if(numOfCouples>0)
233
234  } // if (theCoupleTable)
235
236}
237
238//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
239
240G4double G4DNADingfelderChargeIncreaseModel::CrossSectionPerVolume(const G4Material*,
241                                           const G4ParticleDefinition* particleDefinition,
242                                           G4double k,
243                                           G4double,
244                                           G4double)
245{
246  if (verboseLevel > 3)
247    G4cout << "Calling CrossSectionPerVolume() of G4DNADingfelderChargeIncreaseModel" << G4endl;
248
249 // Calculate total cross section for model
250
251  G4DNAGenericIonsManager *instance;
252  instance = G4DNAGenericIonsManager::Instance();
253
254  if (
255      particleDefinition != instance->GetIon("hydrogen")
256      &&
257      particleDefinition != instance->GetIon("alpha+")
258      &&
259      particleDefinition != instance->GetIon("helium")
260      )
261                   
262    return 0;
263
264  G4double lowLim = 0;
265  G4double highLim = 0;
266  G4double totalCrossSection = 0.;
267 
268  if (flagMaterialIsWater)
269  {
270    const G4String& particleName = particleDefinition->GetParticleName();
271
272    std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
273    pos1 = lowEnergyLimit.find(particleName);
274
275    if (pos1 != lowEnergyLimit.end())
276    {
277      lowLim = pos1->second;
278    }
279
280    std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
281    pos2 = highEnergyLimit.find(particleName);
282
283    if (pos2 != highEnergyLimit.end())
284    {
285      highLim = pos2->second;
286    }
287
288    if (k >= lowLim && k < highLim)
289    {
290      //HYDROGEN
291      if (particleDefinition == instance->GetIon("hydrogen"))
292      {
293          const  G4double aa = 2.835;
294          const  G4double bb = 0.310;
295          const  G4double cc = 2.100;
296          const  G4double dd = 0.760;
297          const  G4double fac = 1.0e-18;
298          const  G4double rr = 13.606 * eV; 
299         
300          G4double t = k / (proton_mass_c2/electron_mass_c2); 
301          G4double x = t / rr; 
302          G4double temp = 4.0 * pi * Bohr_radius/nm * Bohr_radius/nm * fac;
303          G4double sigmal =  temp * cc * (std::pow(x,dd));
304          G4double sigmah = temp * (aa * std::log(1.0 + x) + bb) / x;
305          totalCrossSection = 1.0/(1.0/sigmal + 1.0/sigmah) *m*m;
306      }
307      else
308      {
309          totalCrossSection = Sum(k,particleDefinition);
310      }
311    }
312 
313    if (verboseLevel > 3)
314    {
315      G4cout << "---> Kinetic energy(eV)=" << k/eV << G4endl;
316      G4cout << " - Cross section per water molecule (cm^2)=" << totalCrossSection/cm/cm << G4endl;
317      G4cout << " - Cross section per water molecule (cm^-1)=" << totalCrossSection*densityWater/(1./cm) << G4endl;
318    } 
319
320  } // if (flagMaterialIsWater)
321
322  return totalCrossSection*densityWater;                   
323
324}
325
326//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
327
328void G4DNADingfelderChargeIncreaseModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
329                                              const G4MaterialCutsCouple* /*couple*/,
330                                              const G4DynamicParticle* aDynamicParticle,
331                                              G4double,
332                                              G4double)
333{
334  if (verboseLevel > 3)
335    G4cout << "Calling SampleSecondaries() of G4DNADingfelderChargeIncreaseModel" << G4endl;
336
337  fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
338  fParticleChangeForGamma->ProposeLocalEnergyDeposit(0.);
339
340  G4ParticleDefinition* definition = aDynamicParticle->GetDefinition();
341 
342  G4double inK = aDynamicParticle->GetKineticEnergy();
343
344  G4int finalStateIndex = RandomSelect(inK,definition);
345
346  G4int n = NumberOfFinalStates(definition,finalStateIndex);
347 
348  G4double outK = inK - IncomingParticleBindingEnergyConstant(definition,finalStateIndex);
349
350  G4DNAGenericIonsManager* instance;
351  instance = G4DNAGenericIonsManager::Instance();
352
353  G4double electronK;
354  if (definition == instance->GetIon("hydrogen")) electronK = inK*electron_mass_c2/proton_mass_c2;
355  else electronK = inK*electron_mass_c2/(3728*MeV);
356 
357  if (outK<0)
358  {
359    G4String message;
360    message="G4DNADingfelderChargeIncreaseModel::SampleSecondaries: final kinetic energy is below 0! Process ";
361  }
362 
363  G4DynamicParticle* dp = new G4DynamicParticle
364 
365  (OutgoingParticleDefinition(definition,finalStateIndex), 
366   aDynamicParticle->GetMomentumDirection(), 
367   outK) ;
368
369  fvect->push_back(dp);
370
371  n = n - 1;
372 
373  while (n>0)
374  {
375    n--;
376    fvect->push_back(new G4DynamicParticle
377    (G4Electron::Electron(), aDynamicParticle->GetMomentumDirection(), electronK) );
378  }
379       
380}
381
382//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
383
384G4int G4DNADingfelderChargeIncreaseModel::NumberOfFinalStates(G4ParticleDefinition* particleDefinition, 
385                                                      G4int finalStateIndex )
386 
387{
388  G4DNAGenericIonsManager* instance;
389  instance = G4DNAGenericIonsManager::Instance();
390  if (particleDefinition == instance->GetIon("hydrogen")) return 2;
391  if (particleDefinition == instance->GetIon("alpha+")) return 2;
392 
393  if (particleDefinition == instance->GetIon("helium")) 
394  {    if (finalStateIndex==0) return 2;
395       return 3;
396  }
397  return 0;
398}
399
400//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
401
402G4ParticleDefinition* G4DNADingfelderChargeIncreaseModel::OutgoingParticleDefinition (G4ParticleDefinition* particleDefinition, 
403                                                                              G4int finalStateIndex) 
404{
405  G4DNAGenericIonsManager * instance(G4DNAGenericIonsManager::Instance());
406  if (particleDefinition == instance->GetIon("hydrogen")) return G4Proton::Proton(); 
407  if (particleDefinition == instance->GetIon("alpha+")) return instance->GetIon("alpha++");
408
409  if (particleDefinition == instance->GetIon("helium"))
410  {
411    if (finalStateIndex==0) return instance->GetIon("alpha+"); 
412    return instance->GetIon("alpha++");
413  }
414  return 0;
415}
416
417//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
418
419G4double G4DNADingfelderChargeIncreaseModel::IncomingParticleBindingEnergyConstant(G4ParticleDefinition* particleDefinition, 
420                                                                           G4int finalStateIndex )
421{
422  G4DNAGenericIonsManager * instance(G4DNAGenericIonsManager::Instance());
423  if(particleDefinition == instance->GetIon("hydrogen")) return 13.6*eV;
424 
425  if(particleDefinition == instance->GetIon("alpha+"))
426  {
427      // Binding energy for    He+ -> He++ + e-    54.509 eV
428      // Binding energy for    He  -> He+  + e-    24.587 eV
429      return 54.509*eV;
430  }
431   
432  if(particleDefinition == instance->GetIon("helium"))
433  {
434      // Binding energy for    He+ -> He++ + e-    54.509 eV
435      // Binding energy for    He  -> He+  + e-    24.587 eV
436
437      if (finalStateIndex==0) return 24.587*eV;
438      return (54.509 + 24.587)*eV;
439  } 
440
441  return 0;
442}
443
444//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
445
446//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
447
448G4double G4DNADingfelderChargeIncreaseModel::PartialCrossSection(G4double k, G4int index, 
449                                                           const G4ParticleDefinition* particleDefinition)
450{
451
452  G4int particleTypeIndex = 0;
453  G4DNAGenericIonsManager *instance;
454  instance = G4DNAGenericIonsManager::Instance();
455
456  if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=0;
457  if (particleDefinition == instance->GetIon("helium")) particleTypeIndex=1;
458
459  //
460  // sigma(T) = f0 10 ^ y(log10(T/eV))
461  //
462  //         /  a0 x + b0                    x < x0
463  //         |
464  // y(x) = <   a0 x + b0 - c0 (x - x0)^d0   x0 <= x < x1
465  //         |
466  //         \  a1 x + b1                    x >= x1
467  //
468  //
469  // f0, a0, a1, b0, b1, c0, d0, x0, x1 are parameters that change for protons and helium (0, +, ++)
470  //
471  // 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)
472  //
473  // From Rad. Phys. and Chem. 59 (2000) 255-275, M. Dingfelder et al.
474  // Inelastic-collision cross sections of liquid water for interactions of energetic proton
475  //
476 
477  if (x1[index][particleTypeIndex]<x0[index][particleTypeIndex])
478  {
479      //
480      // 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)
481      //
482      // x1 = x0 + ((a0 - a1)/(c0 * d0)) ^ (1 / (d0 - 1))
483      //
484      // b1 = (a0 - a1) * x1 + b0 - c0 * (x1 - x0) ^ d0
485      //
486 
487      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.));
488      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]);
489  }
490
491  G4double x(std::log10(k/eV));
492  G4double y;
493 
494  if (x<x0[index][particleTypeIndex])
495    y=a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex];
496  else if (x<x1[index][particleTypeIndex])
497    y=a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex] - c0[index][particleTypeIndex] * std::pow(x - x0[index][particleTypeIndex], d0[index][particleTypeIndex]);
498  else
499    y=a1[index][particleTypeIndex] * x + b1[index][particleTypeIndex];
500
501  return f0[index][particleTypeIndex] * std::pow(10., y)*m*m;
502 
503}
504
505//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
506
507G4int G4DNADingfelderChargeIncreaseModel::RandomSelect(G4double k, 
508                                                        const G4ParticleDefinition* particleDefinition)
509{
510
511  G4int particleTypeIndex = 0;
512  G4DNAGenericIonsManager *instance;
513  instance = G4DNAGenericIonsManager::Instance();
514
515  if (particleDefinition == instance->GetIon("hydrogen")) return 0;
516  if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=0;
517  if (particleDefinition == instance->GetIon("helium")) particleTypeIndex=1;
518
519  const G4int n = numberOfPartialCrossSections[particleTypeIndex];
520  G4double* values(new G4double[n]);
521  G4double value = 0;
522  G4int i = n;
523 
524  while (i>0)
525  {
526      i--;
527      values[i]=PartialCrossSection(k, i, particleDefinition);
528      value+=values[i];
529  }
530 
531  value*=G4UniformRand();
532 
533  i=n;
534  while (i>0)
535  {
536      i--;
537   
538      if (values[i]>value)
539        break;
540 
541      value-=values[i];
542  }
543 
544  delete[] values;
545 
546  return i;
547}
548
549//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
550
551G4double G4DNADingfelderChargeIncreaseModel::Sum(G4double k, const G4ParticleDefinition* particleDefinition)
552{
553  G4int particleTypeIndex = 0;
554  G4DNAGenericIonsManager *instance;
555  instance = G4DNAGenericIonsManager::Instance();
556
557  if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=0;
558  if (particleDefinition == instance->GetIon("helium")) particleTypeIndex=1;
559
560  G4double totalCrossSection = 0.;
561
562  for (G4int i=0; i<numberOfPartialCrossSections[particleTypeIndex]; i++)
563  {
564    totalCrossSection += PartialCrossSection(k,i,particleDefinition);
565  }
566  return totalCrossSection;
567}
Note: See TracBrowser for help on using the repository browser.