source: trunk/source/processes/electromagnetic/lowenergy/src/G4DNARuddIonisationModel.cc @ 1315

Last change on this file since 1315 was 1315, checked in by garnier, 14 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 31.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: G4DNARuddIonisationModel.cc,v 1.17 2010/04/07 20:08:31 sincerti Exp $
27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
28//
29
30#include "G4DNARuddIonisationModel.hh"
31//#include "G4DynamicMolecule.hh"
32
33//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
34
35using namespace std;
36
37//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
38
39G4DNARuddIonisationModel::G4DNARuddIonisationModel(const G4ParticleDefinition*,
40                                             const G4String& nam)
41:G4VEmModel(nam),isInitialised(false)
42{
43  lowEnergyLimitForZ1 = 0 * eV; 
44  lowEnergyLimitForZ2 = 0 * eV; 
45  lowEnergyLimitOfModelForZ1 = 100 * eV; 
46  lowEnergyLimitOfModelForZ2 = 1 * keV; 
47  killBelowEnergyForZ1 = lowEnergyLimitOfModelForZ1; 
48  killBelowEnergyForZ2 = lowEnergyLimitOfModelForZ2; 
49
50  verboseLevel= 0;
51  // Verbosity scale:
52  // 0 = nothing
53  // 1 = warning for energy non-conservation
54  // 2 = details of energy budget
55  // 3 = calculation of cross sections, file openings, sampling of atoms
56  // 4 = entering in methods
57 
58  if( verboseLevel>0 ) 
59  { 
60    G4cout << "Rudd ionisation model is constructed " << G4endl;
61  }
62}
63
64//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
65
66G4DNARuddIonisationModel::~G4DNARuddIonisationModel()
67{ 
68  // Cross section
69 
70  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
71  for (pos = tableData.begin(); pos != tableData.end(); ++pos)
72  {
73    G4DNACrossSectionDataSet* table = pos->second;
74    delete table;
75  }
76 
77}
78
79//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
80
81void G4DNARuddIonisationModel::Initialise(const G4ParticleDefinition* particle,
82                                       const G4DataVector& /*cuts*/)
83{
84
85  if (verboseLevel > 3)
86    G4cout << "Calling G4DNARuddIonisationModel::Initialise()" << G4endl;
87
88  // Energy limits
89
90  G4String fileProton("dna/sigma_ionisation_p_rudd");
91  G4String fileHydrogen("dna/sigma_ionisation_h_rudd");
92  G4String fileAlphaPlusPlus("dna/sigma_ionisation_alphaplusplus_rudd");
93  G4String fileAlphaPlus("dna/sigma_ionisation_alphaplus_rudd");
94  G4String fileHelium("dna/sigma_ionisation_he_rudd");
95
96  G4DNAGenericIonsManager *instance;
97  instance = G4DNAGenericIonsManager::Instance();
98  G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
99  G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
100  G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
101  G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
102  G4ParticleDefinition* heliumDef = instance->GetIon("helium");
103
104  G4String proton;
105  G4String hydrogen;
106  G4String alphaPlusPlus;
107  G4String alphaPlus;
108  G4String helium;
109
110  G4double scaleFactor = 1 * m*m;
111
112  if (protonDef != 0)
113  {
114    proton = protonDef->GetParticleName();
115    tableFile[proton] = fileProton;
116
117    lowEnergyLimit[proton] = lowEnergyLimitForZ1;
118    highEnergyLimit[proton] = 500. * keV;
119
120    // Cross section
121
122    G4DNACrossSectionDataSet* tableProton = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, 
123                                                                           eV,
124                                                                           scaleFactor );
125    tableProton->LoadData(fileProton);
126    tableData[proton] = tableProton;
127  }
128  else
129  {
130    G4Exception("G4DNARuddIonisationModel::Initialise: proton is not defined");
131  }
132
133  if (hydrogenDef != 0)
134  {
135    hydrogen = hydrogenDef->GetParticleName();
136    tableFile[hydrogen] = fileHydrogen;
137
138    lowEnergyLimit[hydrogen] = lowEnergyLimitForZ1;
139    highEnergyLimit[hydrogen] = 100. * MeV;
140
141    // Cross section
142
143    G4DNACrossSectionDataSet* tableHydrogen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
144                                                                             eV,
145                                                                             scaleFactor );
146    tableHydrogen->LoadData(fileHydrogen);
147     
148    tableData[hydrogen] = tableHydrogen;
149  }
150  else
151  {
152    G4Exception("G4DNARuddIonisationModel::Initialise: hydrogen is not defined");
153  }
154
155  if (alphaPlusPlusDef != 0)
156  {
157    alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
158    tableFile[alphaPlusPlus] = fileAlphaPlusPlus;
159
160    lowEnergyLimit[alphaPlusPlus] = lowEnergyLimitForZ2;
161    highEnergyLimit[alphaPlusPlus] = 10. * MeV;
162
163    // Cross section
164
165    G4DNACrossSectionDataSet* tableAlphaPlusPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, 
166                                                                                  eV,
167                                                                                  scaleFactor );
168    tableAlphaPlusPlus->LoadData(fileAlphaPlusPlus);
169     
170    tableData[alphaPlusPlus] = tableAlphaPlusPlus;
171  }
172  else
173  {
174    G4Exception("G4DNARuddIonisationModel::Initialise: alphaPlusPlus is not defined");
175  }
176
177  if (alphaPlusDef != 0)
178  {
179    alphaPlus = alphaPlusDef->GetParticleName();
180    tableFile[alphaPlus] = fileAlphaPlus;
181
182    lowEnergyLimit[alphaPlus] = lowEnergyLimitForZ2;
183    highEnergyLimit[alphaPlus] = 10. * MeV;
184
185    // Cross section
186
187    G4DNACrossSectionDataSet* tableAlphaPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, 
188                                                                              eV,
189                                                                              scaleFactor );
190    tableAlphaPlus->LoadData(fileAlphaPlus);
191    tableData[alphaPlus] = tableAlphaPlus;
192  }
193  else
194  {
195    G4Exception("G4DNARuddIonisationModel::Initialise: alphaPlus is not defined");
196  }
197
198  if (heliumDef != 0)
199  {
200    helium = heliumDef->GetParticleName();
201    tableFile[helium] = fileHelium;
202
203    lowEnergyLimit[helium] = lowEnergyLimitForZ2;
204    highEnergyLimit[helium] = 10. * MeV;
205
206    // Cross section
207
208    G4DNACrossSectionDataSet* tableHelium = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, 
209                                                                           eV,
210                                                                           scaleFactor );
211    tableHelium->LoadData(fileHelium);
212    tableData[helium] = tableHelium;
213    }
214  else
215  {
216    G4Exception("G4DNARuddIonisationModel::Initialise: helium is not defined");
217  }
218
219  if (particle==protonDef) 
220  {
221    SetLowEnergyLimit(lowEnergyLimit[proton]);
222    SetHighEnergyLimit(highEnergyLimit[proton]);
223  }
224
225  if (particle==hydrogenDef) 
226  {
227    SetLowEnergyLimit(lowEnergyLimit[hydrogen]);
228    SetHighEnergyLimit(highEnergyLimit[hydrogen]);
229  }
230
231  if (particle==heliumDef) 
232  {
233    SetLowEnergyLimit(lowEnergyLimit[helium]);
234    SetHighEnergyLimit(highEnergyLimit[helium]);
235  }
236
237  if (particle==alphaPlusDef) 
238  {
239    SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
240    SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
241  }
242
243  if (particle==alphaPlusPlusDef) 
244  {
245    SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
246    SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
247  }
248
249  if( verboseLevel>0 ) 
250  { 
251    G4cout << "Rudd ionisation model is initialized " << G4endl
252           << "Energy range: "
253           << LowEnergyLimit() / eV << " eV - "
254           << HighEnergyLimit() / keV << " keV for "
255           << particle->GetParticleName()
256           << G4endl;
257  }
258
259  //
260 
261  if(!isInitialised) 
262  {
263    isInitialised = true;
264 
265    if(pParticleChange)
266      fParticleChangeForGamma = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
267    else
268      fParticleChangeForGamma = new G4ParticleChangeForGamma();
269  }   
270
271  // InitialiseElementSelectors(particle,cuts);
272
273}
274
275//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
276
277G4double G4DNARuddIonisationModel::CrossSectionPerVolume(const G4Material* material,
278                                           const G4ParticleDefinition* particleDefinition,
279                                           G4double k,
280                                           G4double,
281                                           G4double)
282{
283  if (verboseLevel > 3)
284    G4cout << "Calling CrossSectionPerVolume() of G4DNARuddIonisationModel" << G4endl;
285
286 // Calculate total cross section for model
287
288  G4DNAGenericIonsManager *instance;
289  instance = G4DNAGenericIonsManager::Instance();
290
291  if (
292      particleDefinition != G4Proton::ProtonDefinition()
293      &&
294      particleDefinition != instance->GetIon("hydrogen")
295      &&
296      particleDefinition != instance->GetIon("alpha++")
297      &&
298      particleDefinition != instance->GetIon("alpha+")
299      &&
300      particleDefinition != instance->GetIon("helium")
301     )
302           
303    return 0;
304
305  G4double lowLim = 0;
306 
307  if (     particleDefinition == G4Proton::ProtonDefinition()
308       ||  particleDefinition == instance->GetIon("hydrogen")
309     )
310       
311       lowLim = lowEnergyLimitOfModelForZ1;
312       
313  if (     particleDefinition == instance->GetIon("alpha++")
314       ||  particleDefinition == instance->GetIon("alpha+")
315       ||  particleDefinition == instance->GetIon("helium")
316     )
317       
318       lowLim = lowEnergyLimitOfModelForZ2;
319
320  G4double highLim = 0;
321  G4double sigma=0;
322
323  if (material->GetName() == "G4_WATER")
324  {
325    const G4String& particleName = particleDefinition->GetParticleName();
326   
327    // SI - the following is useless since lowLim is already defined
328    /*
329    std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
330    pos1 = lowEnergyLimit.find(particleName);
331
332    if (pos1 != lowEnergyLimit.end())
333    {
334      lowLim = pos1->second;
335    }
336    */
337
338    std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
339    pos2 = highEnergyLimit.find(particleName);
340
341    if (pos2 != highEnergyLimit.end())
342    {
343      highLim = pos2->second;
344    }
345
346    if (k < highLim)
347    {
348     
349      //SI : XS must not be zero otherwise sampling of secondaries method ignored
350
351      if (k < lowLim) k = lowLim;
352
353      //     
354     
355      std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
356      pos = tableData.find(particleName);
357       
358      if (pos != tableData.end())
359      {
360         G4DNACrossSectionDataSet* table = pos->second;
361         if (table != 0)
362         {
363              sigma = table->FindValue(k);
364
365              // BEGIN ELECTRON CORRECTION
366              // add ONE or TWO electron-water excitation for alpha+ and helium
367   
368              if ( particleDefinition == instance->GetIon("alpha+") 
369                   ||
370                   particleDefinition == instance->GetIon("helium")
371                   ) 
372              {
373     
374                  G4DNACrossSectionDataSet* electronDataset = new G4DNACrossSectionDataSet
375                    (new G4LogLogInterpolation, eV, (1./3.343e22)*m*m);
376       
377                  electronDataset->LoadData("dna/sigma_ionisation_e_born");
378
379                  G4double kElectron = k * 0.511/3728;
380
381                  if ( particleDefinition == instance->GetIon("alpha+") ) 
382                  {
383                      G4double tmp1 = table->FindValue(k) + electronDataset->FindValue(kElectron);
384                      delete electronDataset;
385                      if (verboseLevel > 3)
386                      {
387                        G4cout << "---> Kinetic energy(eV)=" << k/eV << G4endl;
388                        G4cout << " - Cross section per water molecule (cm^2)=" << tmp1/cm/cm << G4endl;
389                        G4cout << " - Cross section per water molecule (cm^-1)=" << 
390                        tmp1*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
391                      } 
392                      return tmp1*material->GetAtomicNumDensityVector()[1];
393                  }
394
395                  if ( particleDefinition == instance->GetIon("helium") ) 
396                  {
397                      G4double tmp2 = table->FindValue(k) +  2. * electronDataset->FindValue(kElectron);
398                      delete electronDataset;
399                      if (verboseLevel > 3)
400                      {
401                        G4cout << "---> Kinetic energy(eV)=" << k/eV << G4endl;
402                        G4cout << " - Cross section per water molecule (cm^2)=" << tmp2/cm/cm << G4endl;
403                        G4cout << " - Cross section per water molecule (cm^-1)=" << tmp2*
404                        material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
405                      } 
406                      return tmp2*material->GetAtomicNumDensityVector()[1];
407                  }
408              }     
409
410              // END ELECTRON CORRECTION
411         }
412      }
413      else
414      {
415        G4Exception("G4DNARuddIonisationModel::CrossSectionPerVolume: attempting to calculate cross section for wrong particle");
416      }
417 
418    } // if (k >= lowLim && k < highLim)
419   
420    if (verboseLevel > 3)
421    {
422      G4cout << "---> Kinetic energy(eV)=" << k/eV << G4endl;
423      G4cout << " - Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
424      G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*
425      material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
426    } 
427 
428  } // if (waterMaterial)
429 
430 return sigma*material->GetAtomicNumDensityVector()[1];           
431
432}
433
434//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
435
436void G4DNARuddIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
437                                              const G4MaterialCutsCouple* /*couple*/,
438                                              const G4DynamicParticle* particle,
439                                              G4double,
440                                              G4double)
441{
442  if (verboseLevel > 3)
443    G4cout << "Calling SampleSecondaries() of G4DNARuddIonisationModel" << G4endl;
444
445  G4double lowLim = 0;
446  G4double highLim = 0;
447
448  G4DNAGenericIonsManager *instance;
449  instance = G4DNAGenericIonsManager::Instance();
450
451  if (     particle->GetDefinition() == G4Proton::ProtonDefinition()
452       ||  particle->GetDefinition() == instance->GetIon("hydrogen")
453     )
454       
455       lowLim = killBelowEnergyForZ1;
456       
457  if (     particle->GetDefinition() == instance->GetIon("alpha++")
458       ||  particle->GetDefinition() == instance->GetIon("alpha+")
459       ||  particle->GetDefinition() == instance->GetIon("helium")
460     )
461       
462       lowLim = killBelowEnergyForZ2;
463
464  G4double k = particle->GetKineticEnergy();
465
466  const G4String& particleName = particle->GetDefinition()->GetParticleName();
467
468  // SI - the following is useless since lowLim is already defined
469  /*
470  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
471  pos1 = lowEnergyLimit.find(particleName);
472
473  if (pos1 != lowEnergyLimit.end())
474  {
475    lowLim = pos1->second;
476  }
477  */
478
479  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
480  pos2 = highEnergyLimit.find(particleName);
481
482  if (pos2 != highEnergyLimit.end())
483  {
484    highLim = pos2->second;
485  }
486
487  if (k >= lowLim && k < highLim)
488  {
489      G4ParticleDefinition* definition = particle->GetDefinition();
490      G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
491      G4double particleMass = definition->GetPDGMass();
492      G4double totalEnergy = k + particleMass;
493      G4double pSquare = k*(totalEnergy+particleMass);
494      G4double totalMomentum = std::sqrt(pSquare);
495
496      G4int ionizationShell = RandomSelect(k,particleName);
497 
498      G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell);
499 
500      G4double bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
501
502      G4double cosTheta = 0.;
503      G4double phi = 0.; 
504      RandomizeEjectedElectronDirection(definition, k,secondaryKinetic, cosTheta, phi);
505
506      G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
507      G4double dirX = sinTheta*std::cos(phi);
508      G4double dirY = sinTheta*std::sin(phi);
509      G4double dirZ = cosTheta;
510      G4ThreeVector deltaDirection(dirX,dirY,dirZ);
511      deltaDirection.rotateUz(primaryDirection);
512
513      G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
514
515      G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
516      G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
517      G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
518      G4double finalMomentum = std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
519      finalPx /= finalMomentum;
520      finalPy /= finalMomentum;
521      finalPz /= finalMomentum;
522
523      G4ThreeVector direction;
524      direction.set(finalPx,finalPy,finalPz);
525
526      fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
527      fParticleChangeForGamma->SetProposedKineticEnergy(k-bindingEnergy-secondaryKinetic);
528      fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy);
529
530      G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
531      fvect->push_back(dp);
532
533  }
534
535  // SI - not useful since low energy of model is 0 eV
536 
537  if (k < lowLim) 
538  { 
539    fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
540    fParticleChangeForGamma->ProposeLocalEnergyDeposit(k);
541  }
542
543}
544
545//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
546
547G4double G4DNARuddIonisationModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition, 
548                                                                    G4double k, 
549                                                                    G4int shell)
550{
551  G4double maximumKineticEnergyTransfer = 0.;
552 
553  G4DNAGenericIonsManager *instance;
554  instance = G4DNAGenericIonsManager::Instance();
555
556  if (particleDefinition == G4Proton::ProtonDefinition() 
557      || particleDefinition == instance->GetIon("hydrogen"))
558  { 
559      maximumKineticEnergyTransfer= 4.* (electron_mass_c2 / proton_mass_c2) * k;
560  }
561
562  if (particleDefinition == instance->GetIon("helium") 
563      || particleDefinition == instance->GetIon("alpha+")
564      || particleDefinition == instance->GetIon("alpha++"))
565  { 
566      maximumKineticEnergyTransfer= 4.* (0.511 / 3728) * k;
567  }
568
569  G4double crossSectionMaximum = 0.;
570 
571  for(G4double value=waterStructure.IonisationEnergy(shell); value<=5.*waterStructure.IonisationEnergy(shell) && k>=value ; value+=0.1*eV)
572  {
573      G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k, value, shell);
574      if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
575  }
576 
577  G4double secElecKinetic = 0.;
578 
579  do
580  {
581    secElecKinetic = G4UniformRand() * maximumKineticEnergyTransfer;
582  } while(G4UniformRand()*crossSectionMaximum > DifferentialCrossSection(particleDefinition, 
583                                                                         k,
584                                                                         secElecKinetic+waterStructure.IonisationEnergy(shell),
585                                                                         shell));
586
587  return(secElecKinetic);
588}
589
590//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
591
592
593void G4DNARuddIonisationModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition, 
594                                                                   G4double k, 
595                                                                   G4double secKinetic, 
596                                                                   G4double & cosTheta, 
597                                                                   G4double & phi )
598{
599  G4DNAGenericIonsManager *instance;
600  instance = G4DNAGenericIonsManager::Instance();
601
602  G4double maxSecKinetic = 0.;
603 
604  if (particleDefinition == G4Proton::ProtonDefinition() 
605      || particleDefinition == instance->GetIon("hydrogen")) 
606  { 
607      maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
608  }
609 
610  if (particleDefinition == instance->GetIon("helium") 
611      || particleDefinition == instance->GetIon("alpha+")
612      || particleDefinition == instance->GetIon("alpha++"))
613  { 
614      maxSecKinetic = 4.* (0.511 / 3728) * k;
615  }
616 
617  phi = twopi * G4UniformRand();
618  cosTheta = std::sqrt(secKinetic / maxSecKinetic);
619}
620
621//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
622
623
624G4double G4DNARuddIonisationModel::DifferentialCrossSection(G4ParticleDefinition* particleDefinition, 
625                                                              G4double k, 
626                                                              G4double energyTransfer, 
627                                                              G4int ionizationLevelIndex)
628{
629  // Shells ids are 0 1 2 3 4 (4 is k shell)
630  // !!Attention, "energyTransfer" here is the energy transfered to the electron which means
631  //             that the secondary kinetic energy is w = energyTransfer - bindingEnergy
632  //
633  //   ds            S                F1(nu) + w * F2(nu)
634  //  ---- = G(k) * ----     -------------------------------------------
635  //   dw            Bj       (1+w)^3 * [1 + exp{alpha * (w - wc) / nu}]
636  //
637  // w is the secondary electron kinetic Energy in eV
638  //
639  // All the other parameters can be found in Rudd's Papers
640  //
641  // M.Eugene Rudd, 1988, User-Friendly model for the energy distribution of
642  // electrons from protons or electron collisions. Nucl. Tracks Rad. Meas.Vol 16 N0 2/3 pp 219-218
643  //
644
645  const G4int j=ionizationLevelIndex;
646
647  G4double A1 ; 
648  G4double B1 ; 
649  G4double C1 ; 
650  G4double D1 ; 
651  G4double E1 ;
652  G4double A2 ; 
653  G4double B2 ; 
654  G4double C2 ; 
655  G4double D2 ; 
656  G4double alphaConst ;
657
658  const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV};
659
660  if (j == 4) 
661  {
662      //Data For Liquid Water K SHELL from Dingfelder (Protons in Water)
663      A1 = 1.25; 
664      B1 = 0.5; 
665      C1 = 1.00; 
666      D1 = 1.00; 
667      E1 = 3.00; 
668      A2 = 1.10; 
669      B2 = 1.30;
670      C2 = 1.00; 
671      D2 = 0.00; 
672      alphaConst = 0.66;
673  }
674  else 
675  {
676      //Data For Liquid Water from Dingfelder (Protons in Water)
677      A1 = 1.02; 
678      B1 = 82.0; 
679      C1 = 0.45; 
680      D1 = -0.80; 
681      E1 = 0.38; 
682      A2 = 1.07; 
683      B2 = 14.6;
684      C2 = 0.60; 
685      D2 = 0.04; 
686      alphaConst = 0.64;
687  }
688 
689  const G4double n = 2.;
690  const G4double Gj[5] = {0.99, 1.11, 1.11, 0.52, 1.};
691
692  G4DNAGenericIonsManager* instance;
693  instance = G4DNAGenericIonsManager::Instance();
694
695  G4double wBig = (energyTransfer - waterStructure.IonisationEnergy(ionizationLevelIndex));
696  if (wBig<0) return 0.;
697
698  G4double w = wBig / Bj[ionizationLevelIndex];
699  G4double Ry = 13.6*eV;
700
701  G4double tau = 0.;
702
703  if (particleDefinition == G4Proton::ProtonDefinition() 
704      || particleDefinition == instance->GetIon("hydrogen")) 
705  {
706      tau = (electron_mass_c2/proton_mass_c2) * k ;
707  }
708   
709  if ( particleDefinition == instance->GetIon("helium") 
710       || particleDefinition == instance->GetIon("alpha+")
711       || particleDefinition == instance->GetIon("alpha++")) 
712  {
713      tau = (0.511/3728.) * k ;
714  }
715 
716  G4double S = 4.*pi * Bohr_radius*Bohr_radius * n * std::pow((Ry/Bj[ionizationLevelIndex]),2);
717  G4double v2 = tau / Bj[ionizationLevelIndex];
718  G4double v = std::sqrt(v2);
719  G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj[ionizationLevelIndex]));
720
721  G4double L1 = (C1* std::pow(v,(D1))) / (1.+ E1*std::pow(v, (D1+4.)));
722  G4double L2 = C2*std::pow(v,(D2));
723  G4double H1 = (A1*std::log(1.+v2)) / (v2+(B1/v2));
724  G4double H2 = (A2/v2) + (B2/(v2*v2));
725
726  G4double F1 = L1+H1;
727  G4double F2 = (L2*H2)/(L2+H2);
728
729  G4double sigma = CorrectionFactor(particleDefinition, k) 
730    * Gj[j] * (S/Bj[ionizationLevelIndex]) 
731    * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
732
733  if ( (particleDefinition == instance->GetIon("hydrogen")) && (ionizationLevelIndex==4)) 
734
735    sigma = Gj[j] * (S/Bj[ionizationLevelIndex]) 
736    * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
737
738  if (    particleDefinition == G4Proton::ProtonDefinition() 
739          || particleDefinition == instance->GetIon("hydrogen")
740          ) 
741  {
742      return(sigma);
743  }
744
745  if (particleDefinition == instance->GetIon("alpha++") ) 
746  {
747      slaterEffectiveCharge[0]=0.;
748      slaterEffectiveCharge[1]=0.;
749      slaterEffectiveCharge[2]=0.;
750      sCoefficient[0]=0.;
751      sCoefficient[1]=0.;
752      sCoefficient[2]=0.;
753  }
754
755  if (particleDefinition == instance->GetIon("alpha+") ) 
756  {
757      slaterEffectiveCharge[0]=2.0;
758      slaterEffectiveCharge[1]=1.15;
759      slaterEffectiveCharge[2]=1.15;
760      sCoefficient[0]=0.7;
761      sCoefficient[1]=0.15;
762      sCoefficient[2]=0.15;
763  }
764
765  if (particleDefinition == instance->GetIon("helium") ) 
766  {
767      slaterEffectiveCharge[0]=1.7;
768      slaterEffectiveCharge[1]=1.15;
769      slaterEffectiveCharge[2]=1.15;
770      sCoefficient[0]=0.5;
771      sCoefficient[1]=0.25;
772      sCoefficient[2]=0.25;
773  }
774 
775  if (    particleDefinition == instance->GetIon("helium") 
776          || particleDefinition == instance->GetIon("alpha+")
777          || particleDefinition == instance->GetIon("alpha++")
778          ) 
779  {
780      sigma = Gj[j] * (S/Bj[ionizationLevelIndex]) * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
781   
782      G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
783 
784      zEff -= ( sCoefficient[0] * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.) +
785                sCoefficient[1] * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.) +
786                sCoefficient[2] * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.) );
787           
788      return zEff * zEff * sigma ;
789  } 
790 
791  return 0;
792}
793
794//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
795
796G4double G4DNARuddIonisationModel::S_1s(G4double t, 
797                                          G4double energyTransferred, 
798                                          G4double slaterEffectiveChg, 
799                                          G4double shellNumber)
800{
801  // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
802  // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
803 
804  G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
805  G4double value = 1. - std::exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
806 
807  return value;
808}
809
810//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
811
812G4double G4DNARuddIonisationModel::S_2s(G4double t,
813                                          G4double energyTransferred, 
814                                          G4double slaterEffectiveChg, 
815                                          G4double shellNumber)
816{
817  // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
818  // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
819
820  G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
821  G4double value =  1. - std::exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
822
823  return value;
824 
825}
826
827//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
828
829G4double G4DNARuddIonisationModel::S_2p(G4double t, 
830                                          G4double energyTransferred,
831                                          G4double slaterEffectiveChg, 
832                                          G4double shellNumber)
833{
834  // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
835  // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
836
837  G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
838  G4double value =  1. - std::exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r  + 1.);
839
840  return value;
841}
842
843//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
844
845G4double G4DNARuddIonisationModel::R(G4double t,
846                                       G4double energyTransferred,
847                                       G4double slaterEffectiveChg,
848                                       G4double shellNumber) 
849{
850  // tElectron = m_electron / m_alpha * t
851  // Dingfelder, in Chattanooga 2005 proceedings, p 4
852
853  G4double tElectron = 0.511/3728. * t;
854  G4double value = 2. * tElectron * slaterEffectiveChg / (energyTransferred * shellNumber);
855 
856  return value;
857}
858
859//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
860
861G4double G4DNARuddIonisationModel::CorrectionFactor(G4ParticleDefinition* particleDefinition, G4double k) 
862{
863  G4DNAGenericIonsManager *instance;
864  instance = G4DNAGenericIonsManager::Instance();
865
866  if (particleDefinition == G4Proton::Proton()) 
867  {
868      return(1.);
869  }
870  else 
871    if (particleDefinition == instance->GetIon("hydrogen")) 
872    { 
873        G4double value = (std::log10(k/eV)-4.2)/0.5;
874        return((0.8/(1+std::exp(value))) + 0.9);
875    }
876    else 
877    {   
878        return(1.);
879    }
880}
881
882//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
883
884G4int G4DNARuddIonisationModel::RandomSelect(G4double k, const G4String& particle )
885{   
886 
887  // BEGIN PART 1/2 OF ELECTRON CORRECTION
888
889  // add ONE or TWO electron-water ionisation for alpha+ and helium
890   
891  G4DNAGenericIonsManager *instance;
892  instance = G4DNAGenericIonsManager::Instance();
893  G4double kElectron(0);
894
895  G4DNACrossSectionDataSet * electronDataset = new G4DNACrossSectionDataSet (new G4LogLogInterpolation, eV, (1./3.343e22)*m*m);
896 
897  if ( particle == instance->GetIon("alpha+")->GetParticleName()
898       ||
899       particle == instance->GetIon("helium")->GetParticleName()
900       ) 
901  {     
902      electronDataset->LoadData("dna/sigma_ionisation_e_born");
903
904      kElectron = k * 0.511/3728;
905       
906  }     
907 
908  // END PART 1/2 OF ELECTRON CORRECTION
909 
910  G4int level = 0;
911
912  // Retrieve data table corresponding to the current particle type 
913
914  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
915  pos = tableData.find(particle);
916
917  if (pos != tableData.end())
918  {
919      G4DNACrossSectionDataSet* table = pos->second;
920
921      if (table != 0)
922      {
923          G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
924           
925          const size_t n(table->NumberOfComponents());
926          size_t i(n);
927          G4double value = 0.;
928           
929          while (i>0)
930          { 
931              i--;
932              valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
933
934              // BEGIN PART 2/2 OF ELECTRON CORRECTION
935              // Use only electron partial cross sections
936             
937              if (particle == instance->GetIon("alpha+")->GetParticleName()) 
938                {valuesBuffer[i]=table->GetComponent(i)->FindValue(k) + electronDataset->GetComponent(i)->FindValue(kElectron); }
939
940              if (particle == instance->GetIon("helium")->GetParticleName()) 
941                {valuesBuffer[i]=table->GetComponent(i)->FindValue(k) + 2*electronDataset->GetComponent(i)->FindValue(kElectron); }
942
943              // BEGIN PART 2/2 OF ELECTRON CORRECTION
944
945              value += valuesBuffer[i];
946          }
947           
948          value *= G4UniformRand();
949           
950          i = n;
951           
952          while (i > 0)
953          {
954              i--;
955             
956               
957              if (valuesBuffer[i] > value)
958              {
959                  delete[] valuesBuffer;
960                 
961                  if (electronDataset) delete electronDataset;
962                 
963                  return i;
964              }
965              value -= valuesBuffer[i];
966          }
967           
968          if (valuesBuffer) delete[] valuesBuffer;
969           
970      }
971  }
972  else
973  {
974    G4Exception("G4DNARuddIonisationModel::RandomSelect: attempting to calculate cross section for wrong particle");
975  }
976
977  delete electronDataset;
978     
979  return level;
980}
981
982//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
983
984G4double G4DNARuddIonisationModel::PartialCrossSection(const G4Track& track )
985{
986  G4double sigma = 0.;
987
988  const G4DynamicParticle* particle = track.GetDynamicParticle();
989  G4double k = particle->GetKineticEnergy();
990 
991  G4double lowLim = 0;
992  G4double highLim = 0;
993
994  const G4String& particleName = particle->GetDefinition()->GetParticleName();
995
996  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
997  pos1 = lowEnergyLimit.find(particleName);
998
999  if (pos1 != lowEnergyLimit.end())
1000  {
1001    lowLim = pos1->second;
1002  }
1003
1004  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
1005  pos2 = highEnergyLimit.find(particleName);
1006
1007  if (pos2 != highEnergyLimit.end())
1008  {
1009    highLim = pos2->second;
1010  }
1011
1012  if (k >= lowLim && k <= highLim)
1013  {
1014      std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
1015      pos = tableData.find(particleName);
1016       
1017      if (pos != tableData.end())
1018      {
1019          G4DNACrossSectionDataSet* table = pos->second;
1020          if (table != 0)
1021          {
1022              sigma = table->FindValue(k);
1023          }
1024      }
1025      else
1026      {
1027          G4Exception("G4DNARuddIonisationModel::PartialCrossSection: attempting to calculate cross section for wrong particle");
1028      }
1029  }
1030
1031  return sigma;
1032}
1033
1034//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1035
1036G4double G4DNARuddIonisationModel::Sum(G4double /* energy */, const G4String& /* particle */)
1037{
1038  return 0;
1039}
1040
Note: See TracBrowser for help on using the repository browser.