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

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

update par rapport a CVS
File size: 32.7 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.10 2009/08/13 11:32:47 sincerti Exp $
27// GEANT4 tag $Name: emlowen-V09-02-64 $
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  // Test if water material
274
275  flagMaterialIsWater= false;
276  densityWater = 0;
277
278  const G4ProductionCutsTable* theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
279
280  if(theCoupleTable) 
281  {
282    G4int numOfCouples = theCoupleTable->GetTableSize();
283 
284    if(numOfCouples>0) 
285    {
286          for (G4int i=0; i<numOfCouples; i++) 
287          {
288            const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
289            const G4Material* material = couple->GetMaterial();
290
291            if (material->GetName() == "G4_WATER") 
292            {
293              G4double density = material->GetAtomicNumDensityVector()[1];
294              flagMaterialIsWater = true; 
295              densityWater = density; 
296             
297              if (verboseLevel > 3) 
298              G4cout << "****** Water material is found with density(cm^-3)=" << density/(cm*cm*cm) << G4endl;
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    // creating neutral water molechule...
566
567    G4DNAGenericMoleculeManager *instance;
568    instance = G4DNAGenericMoleculeManager::Instance();
569    G4ParticleDefinition* waterDef = NULL;
570    G4Molecule* water = instance->GetMolecule("H2O");
571    waterDef = (G4ParticleDefinition*)water;
572
573    direction.set(0.,0.,0.);
574
575    //G4DynamicParticle* dynamicWater = new G4DynamicParticle(waterDef, direction, bindingEnergy);
576    G4DynamicMolecule* dynamicWater = new G4DynamicMolecule(water, direction, bindingEnergy);
577    //dynamicWater->RemoveElectron(ionizationShell, 1);
578
579    G4DynamicMolecule* dynamicWater2 = new G4DynamicMolecule(water, direction, bindingEnergy);
580    G4DynamicMolecule* dynamicWater3 = new G4DynamicMolecule(water, direction, bindingEnergy);
581    // insertion inside secondaries
582
583    fvect->push_back(dynamicWater);
584    fvect->push_back(dynamicWater2);
585    fvect->push_back(dynamicWater3);
586      */
587  }
588
589  // SI - not useful since low energy of model is 0 eV
590 
591  if (k < lowLim) 
592  { 
593    fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
594    fParticleChangeForGamma->ProposeLocalEnergyDeposit(k);
595  }
596
597}
598
599//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
600
601G4double G4DNARuddIonisationModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition, 
602                                                                    G4double k, 
603                                                                    G4int shell)
604{
605  G4double maximumKineticEnergyTransfer = 0.;
606 
607  G4DNAGenericIonsManager *instance;
608  instance = G4DNAGenericIonsManager::Instance();
609
610  if (particleDefinition == G4Proton::ProtonDefinition() 
611      || particleDefinition == instance->GetIon("hydrogen"))
612  { 
613      maximumKineticEnergyTransfer= 4.* (electron_mass_c2 / proton_mass_c2) * k;
614  }
615
616  if (particleDefinition == instance->GetIon("helium") 
617      || particleDefinition == instance->GetIon("alpha+")
618      || particleDefinition == instance->GetIon("alpha++"))
619  { 
620      maximumKineticEnergyTransfer= 4.* (0.511 / 3728) * k;
621  }
622
623  G4double crossSectionMaximum = 0.;
624 
625  for(G4double value=waterStructure.IonisationEnergy(shell); value<=4.*waterStructure.IonisationEnergy(shell) ; value+=0.1*eV)
626  {
627      G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k, value, shell);
628      if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
629  }
630 
631  G4double secElecKinetic = 0.;
632 
633  do
634  {
635    secElecKinetic = G4UniformRand() * maximumKineticEnergyTransfer;
636  } while(G4UniformRand()*crossSectionMaximum > DifferentialCrossSection(particleDefinition, 
637                                                                         k,
638                                                                         secElecKinetic+waterStructure.IonisationEnergy(shell),
639                                                                         shell));
640
641  return(secElecKinetic);
642}
643
644//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
645
646
647void G4DNARuddIonisationModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition, 
648                                                                   G4double k, 
649                                                                   G4double secKinetic, 
650                                                                   G4double & cosTheta, 
651                                                                   G4double & phi )
652{
653  G4DNAGenericIonsManager *instance;
654  instance = G4DNAGenericIonsManager::Instance();
655
656  G4double maxSecKinetic = 0.;
657 
658  if (particleDefinition == G4Proton::ProtonDefinition() 
659      || particleDefinition == instance->GetIon("hydrogen")) 
660  { 
661      maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
662  }
663 
664  if (particleDefinition == instance->GetIon("helium") 
665      || particleDefinition == instance->GetIon("alpha+")
666      || particleDefinition == instance->GetIon("alpha++"))
667  { 
668      maxSecKinetic = 4.* (0.511 / 3728) * k;
669  }
670 
671  phi = twopi * G4UniformRand();
672  cosTheta = std::sqrt(secKinetic / maxSecKinetic);
673}
674
675//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
676
677
678G4double G4DNARuddIonisationModel::DifferentialCrossSection(G4ParticleDefinition* particleDefinition, 
679                                                              G4double k, 
680                                                              G4double energyTransfer, 
681                                                              G4int ionizationLevelIndex)
682{
683  // Shells ids are 0 1 2 3 4 (4 is k shell)
684  // !!Attention, "energyTransfer" here is the energy transfered to the electron which means
685  //             that the secondary kinetic energy is w = energyTransfer - bindingEnergy
686  //
687  //   ds            S                F1(nu) + w * F2(nu)
688  //  ---- = G(k) * ----     -------------------------------------------
689  //   dw            Bj       (1+w)^3 * [1 + exp{alpha * (w - wc) / nu}]
690  //
691  // w is the secondary electron kinetic Energy in eV
692  //
693  // All the other parameters can be found in Rudd's Papers
694  //
695  // M.Eugene Rudd, 1988, User-Friendly model for the energy distribution of
696  // electrons from protons or electron collisions. Nucl. Tracks Rad. Meas.Vol 16 N0 2/3 pp 219-218
697  //
698
699  const G4int j=ionizationLevelIndex;
700
701  G4double A1 ; 
702  G4double B1 ; 
703  G4double C1 ; 
704  G4double D1 ; 
705  G4double E1 ;
706  G4double A2 ; 
707  G4double B2 ; 
708  G4double C2 ; 
709  G4double D2 ; 
710  G4double alphaConst ;
711
712  if (j == 4) 
713  {
714      //Data For Liquid Water K SHELL from Dingfelder (Protons in Water)
715      A1 = 1.25; 
716      B1 = 0.5; 
717      C1 = 1.00; 
718      D1 = 1.00; 
719      E1 = 3.00; 
720      A2 = 1.10; 
721      B2 = 1.30;
722      C2 = 1.00; 
723      D2 = 0.00; 
724      alphaConst = 0.66;
725  }
726  else 
727  {
728      //Data For Liquid Water from Dingfelder (Protons in Water)
729      A1 = 1.02; 
730      B1 = 82.0; 
731      C1 = 0.45; 
732      D1 = -0.80; 
733      E1 = 0.38; 
734      A2 = 1.07; 
735      B2 = 14.6;
736      C2 = 0.60; 
737      D2 = 0.04; 
738      alphaConst = 0.64;
739  }
740 
741  const G4double n = 2.;
742  const G4double Gj[5] = {0.99, 1.11, 1.11, 0.52, 1.};
743
744  G4DNAGenericIonsManager* instance;
745  instance = G4DNAGenericIonsManager::Instance();
746
747  G4double wBig = (energyTransfer - waterStructure.IonisationEnergy(ionizationLevelIndex));
748  G4double w = wBig / waterStructure.IonisationEnergy(ionizationLevelIndex);
749  G4double Ry = 13.6*eV;
750
751  G4double tau = 0.;
752
753  if (particleDefinition == G4Proton::ProtonDefinition() 
754      || particleDefinition == instance->GetIon("hydrogen")) 
755  {
756      tau = (electron_mass_c2/proton_mass_c2) * k ;
757  }
758   
759  if ( particleDefinition == instance->GetIon("helium") 
760       || particleDefinition == instance->GetIon("alpha+")
761       || particleDefinition == instance->GetIon("alpha++")) 
762  {
763      tau = (0.511/3728.) * k ;
764  }
765 
766  G4double S = 4.*pi * Bohr_radius*Bohr_radius * n * std::pow((Ry/waterStructure.IonisationEnergy(ionizationLevelIndex)),2);
767  G4double v2 = tau / waterStructure.IonisationEnergy(ionizationLevelIndex);
768  G4double v = std::sqrt(v2);
769  G4double wc = 4.*v2 - 2.*v - (Ry/(4.*waterStructure.IonisationEnergy(ionizationLevelIndex)));
770
771  G4double L1 = (C1* std::pow(v,(D1))) / (1.+ E1*std::pow(v, (D1+4.)));
772  G4double L2 = C2*std::pow(v,(D2));
773  G4double H1 = (A1*std::log(1.+v2)) / (v2+(B1/v2));
774  G4double H2 = (A2/v2) + (B2/(v2*v2));
775
776  G4double F1 = L1+H1;
777  G4double F2 = (L2*H2)/(L2+H2);
778
779  G4double sigma = CorrectionFactor(particleDefinition, k/eV) 
780    * Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex)) 
781    * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
782
783  if (    particleDefinition == G4Proton::ProtonDefinition() 
784          || particleDefinition == instance->GetIon("hydrogen")
785          ) 
786  {
787      return(sigma);
788  }
789
790  if (particleDefinition == instance->GetIon("alpha++") ) 
791  {
792      slaterEffectiveCharge[0]=0.;
793      slaterEffectiveCharge[1]=0.;
794      slaterEffectiveCharge[2]=0.;
795      sCoefficient[0]=0.;
796      sCoefficient[1]=0.;
797      sCoefficient[2]=0.;
798  }
799
800  if (particleDefinition == instance->GetIon("alpha+") ) 
801  {
802      slaterEffectiveCharge[0]=2.0;
803      slaterEffectiveCharge[1]=1.15;
804      slaterEffectiveCharge[2]=1.15;
805      sCoefficient[0]=0.7;
806      sCoefficient[1]=0.15;
807      sCoefficient[2]=0.15;
808  }
809
810  if (particleDefinition == instance->GetIon("helium") ) 
811  {
812      slaterEffectiveCharge[0]=1.7;
813      slaterEffectiveCharge[1]=1.15;
814      slaterEffectiveCharge[2]=1.15;
815      sCoefficient[0]=0.5;
816      sCoefficient[1]=0.25;
817      sCoefficient[2]=0.25;
818  }
819 
820  if (    particleDefinition == instance->GetIon("helium") 
821          || particleDefinition == instance->GetIon("alpha+")
822          || particleDefinition == instance->GetIon("alpha++")
823          ) 
824  {
825      sigma = Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex)) * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
826   
827      G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
828 
829      zEff -= ( sCoefficient[0] * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.) +
830                sCoefficient[1] * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.) +
831                sCoefficient[2] * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.) );
832           
833      return zEff * zEff * sigma ;
834  } 
835 
836  return 0;
837}
838
839//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
840
841G4double G4DNARuddIonisationModel::S_1s(G4double t, 
842                                          G4double energyTransferred, 
843                                          G4double slaterEffectiveChg, 
844                                          G4double shellNumber)
845{
846  // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
847  // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
848 
849  G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
850  G4double value = 1. - std::exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
851 
852  return value;
853}
854
855//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
856
857G4double G4DNARuddIonisationModel::S_2s(G4double t,
858                                          G4double energyTransferred, 
859                                          G4double slaterEffectiveChg, 
860                                          G4double shellNumber)
861{
862  // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
863  // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
864
865  G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
866  G4double value =  1. - std::exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
867
868  return value;
869 
870}
871
872//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
873
874G4double G4DNARuddIonisationModel::S_2p(G4double t, 
875                                          G4double energyTransferred,
876                                          G4double slaterEffectiveChg, 
877                                          G4double shellNumber)
878{
879  // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
880  // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
881
882  G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
883  G4double value =  1. - std::exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r  + 1.);
884
885  return value;
886}
887
888//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
889
890G4double G4DNARuddIonisationModel::R(G4double t,
891                                       G4double energyTransferred,
892                                       G4double slaterEffectiveChg,
893                                       G4double shellNumber) 
894{
895  // tElectron = m_electron / m_alpha * t
896  // Dingfelder, in Chattanooga 2005 proceedings, p 4
897
898  G4double tElectron = 0.511/3728. * t;
899  G4double value = 2. * tElectron * slaterEffectiveChg / (energyTransferred * shellNumber);
900 
901  return value;
902}
903
904//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
905
906G4double G4DNARuddIonisationModel::CorrectionFactor(G4ParticleDefinition* particleDefinition, G4double k) 
907{
908  G4DNAGenericIonsManager *instance;
909  instance = G4DNAGenericIonsManager::Instance();
910
911  if (particleDefinition == G4Proton::Proton()) 
912  {
913      return(1.);
914  }
915  else 
916    if (particleDefinition == instance->GetIon("hydrogen")) 
917    { 
918        G4double value = (std::log(k/eV)-4.2)/0.5;
919        return((0.8/(1+std::exp(value))) + 0.9);
920    }
921    else 
922    {   
923        return(1.);
924    }
925}
926
927//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
928
929G4int G4DNARuddIonisationModel::RandomSelect(G4double k, const G4String& particle )
930{   
931 
932  // BEGIN PART 1/2 OF ELECTRON CORRECTION
933
934  // add ONE or TWO electron-water excitation for alpha+ and helium
935   
936  G4DNAGenericIonsManager *instance;
937  instance = G4DNAGenericIonsManager::Instance();
938  G4double kElectron(0);
939  G4double electronComponent(0);
940  G4DNACrossSectionDataSet * electronDataset = new G4DNACrossSectionDataSet (new G4LogLogInterpolation, eV, (1./3.343e22)*m*m);
941 
942  if ( particle == instance->GetIon("alpha+")->GetParticleName()
943       ||
944       particle == instance->GetIon("helium")->GetParticleName()
945       ) 
946  {     
947      electronDataset->LoadData("dna/sigma_ionisation_e_born");
948
949      kElectron = k * 0.511/3728;
950       
951      electronComponent = electronDataset->FindValue(kElectron);
952       
953  }     
954 
955  delete electronDataset;
956 
957  // END PART 1/2 OF ELECTRON CORRECTION
958 
959  G4int level = 0;
960
961  // Retrieve data table corresponding to the current particle type 
962
963  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
964  pos = tableData.find(particle);
965
966  if (pos != tableData.end())
967  {
968      G4DNACrossSectionDataSet* table = pos->second;
969
970      if (table != 0)
971      {
972          G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
973           
974          const size_t n(table->NumberOfComponents());
975          size_t i(n);
976          G4double value = 0.;
977           
978          while (i>0)
979          { 
980              i--;
981              valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
982
983              // BEGIN PART 2/2 OF ELECTRON CORRECTION
984
985              if (particle == instance->GetIon("alpha+")->GetParticleName()) 
986                {valuesBuffer[i]=table->GetComponent(i)->FindValue(k) + electronComponent; }
987     
988              if (particle == instance->GetIon("helium")->GetParticleName()) 
989                {valuesBuffer[i]=table->GetComponent(i)->FindValue(k) + 2*electronComponent; }
990     
991              // BEGIN PART 2/2 OF ELECTRON CORRECTION
992
993              value += valuesBuffer[i];
994          }
995           
996          value *= G4UniformRand();
997           
998          i = n;
999           
1000          while (i > 0)
1001          {
1002              i--;
1003               
1004              if (valuesBuffer[i] > value)
1005              {
1006                  delete[] valuesBuffer;
1007                  return i;
1008              }
1009              value -= valuesBuffer[i];
1010          }
1011           
1012          if (valuesBuffer) delete[] valuesBuffer;
1013           
1014      }
1015  }
1016  else
1017  {
1018    G4Exception("G4DNARuddIonisationModel::RandomSelect: attempting to calculate cross section for wrong particle");
1019  }
1020     
1021  return level;
1022}
1023
1024//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1025
1026G4double G4DNARuddIonisationModel::PartialCrossSection(const G4Track& track )
1027{
1028  G4double sigma = 0.;
1029
1030  const G4DynamicParticle* particle = track.GetDynamicParticle();
1031  G4double k = particle->GetKineticEnergy();
1032 
1033  G4double lowLim = 0;
1034  G4double highLim = 0;
1035
1036  const G4String& particleName = particle->GetDefinition()->GetParticleName();
1037
1038  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
1039  pos1 = lowEnergyLimit.find(particleName);
1040
1041  if (pos1 != lowEnergyLimit.end())
1042  {
1043    lowLim = pos1->second;
1044  }
1045
1046  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
1047  pos2 = highEnergyLimit.find(particleName);
1048
1049  if (pos2 != highEnergyLimit.end())
1050  {
1051    highLim = pos2->second;
1052  }
1053
1054  if (k >= lowLim && k <= highLim)
1055  {
1056      std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
1057      pos = tableData.find(particleName);
1058       
1059      if (pos != tableData.end())
1060      {
1061          G4DNACrossSectionDataSet* table = pos->second;
1062          if (table != 0)
1063          {
1064              sigma = table->FindValue(k);
1065          }
1066      }
1067      else
1068      {
1069          G4Exception("G4DNARuddIonisationModel::PartialCrossSection: attempting to calculate cross section for wrong particle");
1070      }
1071  }
1072
1073  return sigma;
1074}
1075
1076//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1077
1078G4double G4DNARuddIonisationModel::Sum(G4double /* energy */, const G4String& /* particle */)
1079{
1080  return 0;
1081}
1082
Note: See TracBrowser for help on using the repository browser.