source: trunk/source/processes/electromagnetic/lowenergy/src/G4DNAMillerGreenExcitationModel.cc @ 968

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

fichier ajoutes
File size: 18.3 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: G4DNAMillerGreenExcitationModel.cc,v 1.3 2009/02/16 11:00:11 sincerti Exp $
27// GEANT4 tag $Name: geant4-09-02-ref-02 $
28//
29
30#include "G4DNAMillerGreenExcitationModel.hh"
31
32//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
33
34using namespace std;
35
36//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
37
38G4DNAMillerGreenExcitationModel::G4DNAMillerGreenExcitationModel(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  G4cout << "Miller & Green excitation model is constructed " << G4endl;
52
53}
54
55//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
56
57G4DNAMillerGreenExcitationModel::~G4DNAMillerGreenExcitationModel()
58{}
59
60//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
61
62void G4DNAMillerGreenExcitationModel::Initialise(const G4ParticleDefinition* particle,
63                                       const G4DataVector& /*cuts*/)
64{
65
66  if (verboseLevel > 3)
67    G4cout << "Calling G4DNAMillerGreenExcitationModel::Initialise()" << G4endl;
68
69  // Energy limits
70 
71  G4DNAGenericIonsManager *instance;
72  instance = G4DNAGenericIonsManager::Instance();
73  G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
74  G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
75  G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
76  G4ParticleDefinition* heliumDef = instance->GetIon("helium");
77
78  G4String proton;
79  G4String alphaPlusPlus;
80  G4String alphaPlus;
81  G4String helium;
82
83  if (protonDef != 0)
84  {
85    proton = protonDef->GetParticleName();
86    lowEnergyLimit[proton] = 10. * eV;
87    highEnergyLimit[proton] = 500. * keV;
88   
89    kineticEnergyCorrection[0] = 1.;
90    slaterEffectiveCharge[0][0] = 0.;
91    slaterEffectiveCharge[1][0] = 0.;
92    slaterEffectiveCharge[2][0] = 0.;
93    sCoefficient[0][0] = 0.;
94    sCoefficient[1][0] = 0.;
95    sCoefficient[2][0] = 0.;
96  }
97  else
98  {
99    G4Exception("G4DNAMillerGreenExcitationModel::Initialise: proton is not defined");
100  }
101
102  if (alphaPlusPlusDef != 0)
103  {
104    alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
105    lowEnergyLimit[alphaPlusPlus] = 1. * keV;
106    highEnergyLimit[alphaPlusPlus] = 10. * MeV;
107
108    kineticEnergyCorrection[1] = 0.9382723/3.727417;
109    slaterEffectiveCharge[0][1]=0.;
110    slaterEffectiveCharge[1][1]=0.;
111    slaterEffectiveCharge[2][1]=0.;
112    sCoefficient[0][1]=0.;
113    sCoefficient[1][1]=0.;
114    sCoefficient[2][1]=0.;
115  }
116  else
117  {
118      G4Exception("G4DNAMillerGreenExcitationModel::Initialise: alphaPlusPlus is not defined");
119  }
120
121  if (alphaPlusDef != 0)
122  {
123    alphaPlus = alphaPlusDef->GetParticleName();
124    lowEnergyLimit[alphaPlus] = 1. * keV;
125    highEnergyLimit[alphaPlus] = 10. * MeV;
126
127    kineticEnergyCorrection[2] = 0.9382723/3.727417;
128    slaterEffectiveCharge[0][2]=2.0;
129    slaterEffectiveCharge[1][2]=1.15;
130    slaterEffectiveCharge[2][2]=1.15;
131    sCoefficient[0][2]=0.7;
132    sCoefficient[1][2]=0.15;
133    sCoefficient[2][2]=0.15;
134  }
135  else
136  {
137    G4Exception("G4DNAMillerGreenExcitationModel::Initialise: alphaPlus is not defined");
138  }
139
140  if (heliumDef != 0)
141  {
142    helium = heliumDef->GetParticleName();
143    lowEnergyLimit[helium] = 1. * keV;
144    highEnergyLimit[helium] = 10. * MeV;
145   
146    kineticEnergyCorrection[3] = 0.9382723/3.727417;
147    slaterEffectiveCharge[0][3]=1.7;
148    slaterEffectiveCharge[1][3]=1.15;
149    slaterEffectiveCharge[2][3]=1.15;
150    sCoefficient[0][3]=0.5;
151    sCoefficient[1][3]=0.25;
152    sCoefficient[2][3]=0.25;
153  }
154  else
155  {
156    G4Exception("G4DNAMillerGreenExcitationModel::Initialise: helium is not defined");
157  }
158
159  if (particle==protonDef) 
160  {
161    SetLowEnergyLimit(lowEnergyLimit[proton]);
162    SetHighEnergyLimit(highEnergyLimit[proton]);
163  }
164
165  if (particle==alphaPlusPlusDef) 
166  {
167    SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
168    SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
169  }
170
171  if (particle==alphaPlusDef) 
172  {
173    SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
174    SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
175  }
176
177  if (particle==heliumDef) 
178  {
179    SetLowEnergyLimit(lowEnergyLimit[helium]);
180    SetHighEnergyLimit(highEnergyLimit[helium]);
181  }
182
183  //
184 
185  nLevels = waterExcitation.NumberOfLevels();
186
187  //
188 
189  G4cout << "Miller & Green excitation model is initialized " << G4endl
190         << "Energy range: "
191         << LowEnergyLimit() / eV << " eV - "
192         << HighEnergyLimit() / keV << " keV for "
193         << particle->GetParticleName()
194         << G4endl;
195
196  if(!isInitialised) 
197  {
198    isInitialised = true;
199 
200    if(pParticleChange)
201      fParticleChangeForGamma = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
202    else
203      fParticleChangeForGamma = new G4ParticleChangeForGamma();
204  }   
205
206  // InitialiseElementSelectors(particle,cuts);
207 
208  // Test if water material
209
210  flagMaterialIsWater= false;
211  densityWater = 0;
212
213  const G4ProductionCutsTable* theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
214
215  if(theCoupleTable) 
216  {
217    G4int numOfCouples = theCoupleTable->GetTableSize();
218 
219    if(numOfCouples>0) 
220    {
221          for (G4int i=0; i<numOfCouples; i++) 
222          {
223            const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
224            const G4Material* material = couple->GetMaterial();
225
226            size_t j = material->GetNumberOfElements();
227            while (j>0)
228            {
229               j--;
230               const G4Element* element(material->GetElement(j));
231               if (element->GetZ() == 8.)
232               {
233                  G4double density = material->GetAtomicNumDensityVector()[j];
234                  if (density > 0.) 
235                  { 
236                    flagMaterialIsWater = true; 
237                    densityWater = density; 
238                   
239                    if (verboseLevel > 3) 
240                    G4cout << "Water material is found with density(cm^-3)=" << density/(cm*cm*cm) << G4endl;
241                  }
242               }
243            }
244 
245          }
246   } // if(numOfCouples>0)
247
248  } // if (theCoupleTable)
249
250}
251
252//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
253
254G4double G4DNAMillerGreenExcitationModel::CrossSectionPerVolume(const G4Material* material,
255                                           const G4ParticleDefinition* particleDefinition,
256                                           G4double k,
257                                           G4double,
258                                           G4double)
259{
260  if (verboseLevel > 3)
261    G4cout << "Calling CrossSectionPerVolume() of G4DNAMillerGreenExcitationModel" << G4endl;
262
263 // Calculate total cross section for model
264
265  G4DNAGenericIonsManager *instance;
266  instance = G4DNAGenericIonsManager::Instance();
267
268  if (
269      particleDefinition != G4Proton::ProtonDefinition()
270      &&
271      particleDefinition != instance->GetIon("alpha++")
272      &&
273      particleDefinition != instance->GetIon("alpha+")
274      &&
275      particleDefinition != instance->GetIon("helium")
276     )
277           
278    return 0;
279
280  G4double lowLim = 0;
281  G4double highLim = 0;
282  G4double crossSection = 0.;
283
284  if (flagMaterialIsWater)
285  {
286    const G4String& particleName = particleDefinition->GetParticleName();
287
288    std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
289    pos1 = lowEnergyLimit.find(particleName);
290
291    if (pos1 != lowEnergyLimit.end())
292    {
293      lowLim = pos1->second;
294    }
295
296    std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
297    pos2 = highEnergyLimit.find(particleName);
298
299    if (pos2 != highEnergyLimit.end())
300    {
301      highLim = pos2->second;
302    }
303
304    if (k >= lowLim && k < highLim)
305    {
306      crossSection = Sum(k,particleDefinition);
307
308      G4DNAGenericIonsManager *instance;
309      instance = G4DNAGenericIonsManager::Instance();
310
311      // add ONE or TWO electron-water excitation for alpha+ and helium
312 
313      if ( particleDefinition == instance->GetIon("alpha+") 
314           ||
315           particleDefinition == instance->GetIon("helium")
316         ) 
317      {
318          G4DNAEmfietzoglouExcitationModel * excitationXS = new G4DNAEmfietzoglouExcitationModel();
319
320          G4double sigmaExcitation=0;
321          G4double tmp =0.;
322         
323          if (k*0.511/3728 > 7.4*eV && k*0.511/3728 < 10*keV) sigmaExcitation = 
324            excitationXS->CrossSectionPerVolume(material,particleDefinition,k*0.511/3728,tmp,tmp)/densityWater;
325       
326          if ( particleDefinition == instance->GetIon("alpha+") ) 
327            crossSection = crossSection +  sigmaExcitation ;
328         
329          if ( particleDefinition == instance->GetIon("helium") ) 
330            crossSection = crossSection + 2*sigmaExcitation ;
331         
332          delete excitationXS;
333      }     
334
335    }
336
337    if (verboseLevel > 3)
338    {
339      G4cout << "---> Kinetic energy(eV)=" << k/eV << G4endl;
340      G4cout << " - Cross section per water molecule (cm^2)=" << crossSection/cm/cm << G4endl;
341      G4cout << " - Cross section per water molecule (cm^-1)=" << crossSection*densityWater/(1./cm) << G4endl;
342    } 
343
344  } // if (flagMaterialIsWater)
345
346 return crossSection*densityWater;                 
347
348}
349
350//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
351
352void G4DNAMillerGreenExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
353                                              const G4MaterialCutsCouple* /*couple*/,
354                                              const G4DynamicParticle* aDynamicParticle,
355                                              G4double,
356                                              G4double)
357{
358
359  if (verboseLevel > 3)
360    G4cout << "Calling SampleSecondaries() of G4DNAMillerGreenExcitationModel" << G4endl;
361
362  G4double particleEnergy0 = aDynamicParticle->GetKineticEnergy();
363 
364  G4int level = RandomSelect(particleEnergy0,aDynamicParticle->GetDefinition());
365
366  G4double excitationEnergy = waterExcitation.ExcitationEnergy(level);
367  G4double newEnergy = particleEnergy0 - excitationEnergy;
368 
369  if (newEnergy>0)
370  {
371      fParticleChangeForGamma->ProposeMomentumDirection(aDynamicParticle->GetMomentumDirection());
372      fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
373      fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
374  }
375
376}
377
378//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
379
380G4double G4DNAMillerGreenExcitationModel::PartialCrossSection(G4double k, G4int excitationLevel, 
381                                                                  const G4ParticleDefinition* particleDefinition)
382{
383  //                               ( ( z * aj ) ^ omegaj ) * ( t - ej ) ^ nu
384  // sigma(t) = zEff^2 * sigma0 * --------------------------------------------
385  //                               jj ^ ( omegaj + nu ) + t ^ ( omegaj + nu )
386  //
387  // where t is the kinetic energy corrected by Helium mass over proton mass for Helium ions
388  //
389  // zEff is:
390  //  1 for protons
391  //  2 for alpha++
392  //  and  2 - c1 S_1s - c2 S_2s - c3 S_2p for alpha+ and He
393  //
394  // Dingfelder et al., RPC 59, 255-275, 2000 from Miller and Green (1973)
395  // Formula (34) and Table 2
396 
397  const G4double sigma0(1.E+8 * barn);
398  const G4double nu(1.);
399  const G4double aj[]={876.*eV, 2084.* eV, 1373.*eV, 692.*eV, 900.*eV};
400  const G4double jj[]={19820.*eV, 23490.*eV, 27770.*eV, 30830.*eV, 33080.*eV};
401  const G4double omegaj[]={0.85, 0.88, 0.88, 0.78, 0.78};
402 
403  G4int particleTypeIndex = 0;
404  G4DNAGenericIonsManager* instance;
405  instance = G4DNAGenericIonsManager::Instance();
406
407  if (particleDefinition == G4Proton::ProtonDefinition()) particleTypeIndex=0;
408  if (particleDefinition == instance->GetIon("alpha++")) particleTypeIndex=1;
409  if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=2;
410  if (particleDefinition == instance->GetIon("helium")) particleTypeIndex=3;
411 
412  G4double tCorrected;
413  tCorrected = k * kineticEnergyCorrection[particleTypeIndex];
414
415  // SI - added protection
416  if (tCorrected < waterExcitation.ExcitationEnergy(excitationLevel)) return 0;
417  //
418 
419  G4int z = 10;
420
421  G4double numerator;
422  numerator = std::pow(z * aj[excitationLevel], omegaj[excitationLevel]) * 
423    std::pow(tCorrected - waterExcitation.ExcitationEnergy(excitationLevel), nu);
424
425  G4double power;
426  power = omegaj[excitationLevel] + nu;
427
428  G4double denominator;
429  denominator = std::pow(jj[excitationLevel], power) + std::pow(tCorrected, power);
430
431  G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
432
433  zEff -= ( sCoefficient[0][particleTypeIndex] * S_1s(k, waterExcitation.ExcitationEnergy(excitationLevel), slaterEffectiveCharge[0][particleTypeIndex], 1.) +
434            sCoefficient[1][particleTypeIndex] * S_2s(k, waterExcitation.ExcitationEnergy(excitationLevel), slaterEffectiveCharge[1][particleTypeIndex], 2.) +
435            sCoefficient[2][particleTypeIndex] * S_2p(k, waterExcitation.ExcitationEnergy(excitationLevel), slaterEffectiveCharge[2][particleTypeIndex], 2.) );
436
437  G4double cross = sigma0 * zEff * zEff * numerator / denominator;
438
439  return cross;
440}
441
442//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
443
444G4int G4DNAMillerGreenExcitationModel::RandomSelect(G4double k,const G4ParticleDefinition* particle)
445{
446  G4int i = nLevels;
447  G4double value = 0.;
448  std::deque<double> values;
449 
450  G4DNAGenericIonsManager *instance;
451  instance = G4DNAGenericIonsManager::Instance();
452
453  if ( particle == instance->GetIon("alpha++") ||
454       particle == G4Proton::ProtonDefinition()  )
455  { 
456     while (i > 0)
457     {
458        i--;
459        G4double partial = PartialCrossSection(k,i,particle);
460        values.push_front(partial);
461        value += partial;
462     }
463
464     value *= G4UniformRand();
465   
466     i = nLevels;
467
468     while (i > 0)
469     {
470        i--;
471        if (values[i] > value) return i;
472        value -= values[i];
473     }
474  } 
475
476  // add ONE or TWO electron-water excitation for alpha+ and helium
477   
478  if ( particle == instance->GetIon("alpha+") 
479       ||
480       particle == instance->GetIon("helium")
481     ) 
482  {
483    while (i>0)
484    {
485          i--;
486         
487          G4DNAEmfietzoglouExcitationModel * excitationXS = new G4DNAEmfietzoglouExcitationModel();
488         
489          G4double sigmaExcitation=0;
490
491          if (k*0.511/3728 > 7.4*eV && k*0.511/3728 < 10*keV) sigmaExcitation = excitationXS->PartialCrossSection(k*0.511/3728,i);
492 
493          G4double partial = PartialCrossSection(k,i,particle);
494          if (particle == instance->GetIon("alpha+")) partial = PartialCrossSection(k,i,particle) + sigmaExcitation;
495          if (particle == instance->GetIon("helium")) partial = PartialCrossSection(k,i,particle) + 2*sigmaExcitation;
496          values.push_front(partial);
497          value += partial;
498          delete excitationXS;
499    }
500 
501    value*=G4UniformRand();
502 
503    i=5;
504    while (i>0)
505    {
506          i--;
507   
508          if (values[i]>value) return i;
509 
510          value-=values[i];
511    }
512  }     
513
514  return 0;
515}
516
517//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
518
519G4double G4DNAMillerGreenExcitationModel::Sum(G4double k, const G4ParticleDefinition* particle)
520{
521  G4double totalCrossSection = 0.;
522
523  for (G4int i=0; i<nLevels; i++)
524  {
525    totalCrossSection += PartialCrossSection(k,i,particle);
526  }
527  return totalCrossSection;
528}
529
530//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
531
532G4double G4DNAMillerGreenExcitationModel::S_1s(G4double t,
533                                                          G4double energyTransferred,
534                                                          G4double slaterEffectiveCharge,
535                                                          G4double shellNumber)
536{
537  // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
538  // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
539 
540  G4double r = R(t, energyTransferred, slaterEffectiveCharge, shellNumber);
541  G4double value = 1. - std::exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
542 
543  return value;
544}
545
546
547//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
548
549G4double G4DNAMillerGreenExcitationModel::S_2s(G4double t,
550                                                          G4double energyTransferred,
551                                                          G4double slaterEffectiveCharge,
552                                                          G4double shellNumber)
553{
554  // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
555  // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
556
557  G4double r = R(t, energyTransferred, slaterEffectiveCharge, shellNumber);
558  G4double value =  1. - std::exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
559
560  return value;
561 
562}
563
564//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
565
566G4double G4DNAMillerGreenExcitationModel::S_2p(G4double t,
567                                                          G4double energyTransferred,
568                                                          G4double slaterEffectiveCharge,
569                                                          G4double shellNumber)
570{
571  // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
572  // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
573
574  G4double r = R(t, energyTransferred, slaterEffectiveCharge, shellNumber);
575  G4double value =  1. - std::exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r  + 1.);
576
577  return value;
578}
579
580//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
581
582G4double G4DNAMillerGreenExcitationModel::R(G4double t,
583                                                       G4double energyTransferred,
584                                                       G4double slaterEffectiveCharge,
585                                                       G4double shellNumber) 
586{
587  // tElectron = m_electron / m_alpha * t
588  // Dingfelder, in Chattanooga 2005 proceedings, p 4
589
590  G4double tElectron = 0.511/3728. * t;
591  G4double value = 2. * tElectron * slaterEffectiveCharge / (energyTransferred * shellNumber);
592 
593  return value;
594}
595
Note: See TracBrowser for help on using the repository browser.