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

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

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