source: trunk/source/processes/electromagnetic/lowenergy/src/G4DNARuddIonisationExtendedModel.cc @ 1350

Last change on this file since 1350 was 1350, checked in by garnier, 13 years ago

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