Ignore:
Timestamp:
Apr 6, 2009, 12:21:12 PM (15 years ago)
Author:
garnier
Message:

update processes

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/electromagnetic/lowenergy/src/G4FinalStateIonisationRudd.cc

    r819 r961  
    2424// ********************************************************************
    2525//
    26 //
    27 // $Id: G4FinalStateIonisationRudd.cc,v 1.5 2007/11/26 17:27:09 pia Exp $
    28 // GEANT4 tag $Name:  $
    29 //
    30 // Contact Author: Sebastien Incerti (incerti@cenbg.in2p3.fr)
    31 //                 Maria Grazia Pia  (Maria.Grazia.Pia@cern.ch)
    32 //
    33 ///
    34 // Reference: TNS Geant4-DNA paper
    35 // Reference for implementation model: NIM. 155, pp. 145-156, 1978
    36 //
    37 // History:
    38 // -----------
    39 // Date         Name              Modification
    40 // 28 Apr 2007  M.G. Pia          Created in compliance with design described in TNS paper
    41 //    Nov 2007  S. Incerti        Implementation
    42 // 26 Nov 2007  MGP               Cleaned up std::
    43 //
    44 // -------------------------------------------------------------------
    45 
    46 // Class description:
    47 // Reference: TNS Geant4-DNA paper
    48 // S. Chauvie et al., Geant4 physics processes for microdosimetry simulation:
    49 // design foundation and implementation of the first set of models,
    50 // IEEE Trans. Nucl. Sci., vol. 54, no. 6, Dec. 2007.
    51 // Further documentation available from http://www.ge.infn.it/geant4/dna
    52 
    53 // -------------------------------------------------------------------
    54 
     26// $Id: G4FinalStateIonisationRudd.cc,v 1.8 2008/08/20 14:51:48 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    5528
    5629#include "G4FinalStateIonisationRudd.hh"
    57 #include "G4Track.hh"
    58 #include "G4Step.hh"
    59 #include "G4DynamicParticle.hh"
    60 #include "Randomize.hh"
    61 
    62 #include "G4ParticleTypes.hh"
    63 #include "G4ParticleDefinition.hh"
    64 #include "G4Electron.hh"
    65 #include "G4Proton.hh"
    66 #include "G4SystemOfUnits.hh"
    67 #include "G4ParticleMomentum.hh"
    68 #include "G4DNAGenericIonsManager.hh"
    69 
     30
     31//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    7032
    7133G4FinalStateIonisationRudd::G4FinalStateIonisationRudd()
    7234{
    73   name = "IonisationBorn";
    74   // Default energy limits (defined for protection against anomalous behaviour only)
    7535  lowEnergyLimitDefault = 100 * eV;
    7636  highEnergyLimitDefault = 100 * MeV;
    … …  
    11272}
    11373
     74//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    11475
    11576G4FinalStateIonisationRudd::~G4FinalStateIonisationRudd()
    116 { }
    117 
    118 
     77{}
     78
     79//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    11980
    12081const G4FinalStateProduct& G4FinalStateIonisationRudd::GenerateFinalState(const G4Track& track, const G4Step& /* step */)
    12182{
    122   // Clear previous secondaries, energy deposit and particle kill status
    12383  product.Clear();
    12484
    … …  
    13292  const G4String& particleName = particle->GetDefinition()->GetParticleName();
    13393
    134   // Retrieve energy limits for the current particle type
    135 
    13694  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
    13795  pos1 = lowEnergyLimit.find(particleName);
    13896
    139   // Lower limit
    14097  if (pos1 != lowEnergyLimit.end())
    141     {
    142       lowLim = pos1->second;
    143     }
    144 
    145   // Upper limit
     98  {
     99    lowLim = pos1->second;
     100  }
     101
    146102  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
    147103  pos2 = highEnergyLimit.find(particleName);
    148104
    149105  if (pos2 != highEnergyLimit.end())
    150     {
    151       highLim = pos2->second;
    152     }
    153 
    154   // Verify that the current track is within the energy limits of validity of the cross section model
     106  {
     107    highLim = pos2->second;
     108  }
    155109
    156110  if (k >= lowLim && k <= highLim)
    157     {
    158       // Kinetic energy of primary particle
    159 
     111  {
    160112      G4ParticleDefinition* definition = particle->GetDefinition();
    161113      G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
    … …  
    186138      G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
    187139
    188       // Primary Particle Direction
    189140      G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
    190141      G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
    … …  
    200151      G4DynamicParticle* aElectron = new G4DynamicParticle(G4Electron::Electron(),deltaDirection,secondaryKinetic);
    201152      product.AddSecondary(aElectron);
    202     }
    203 
    204   if (k < lowLim) {product.KillPrimaryParticle();product.AddEnergyDeposit(k);}
     153  }
     154
     155  if (k < lowLim)
     156  { 
     157    product.KillPrimaryParticle();
     158  }
    205159 
    206160  return product;
    207161}
    208162
    209 
     163//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    210164
    211165G4double G4FinalStateIonisationRudd::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
    … …  
    220174  if (particleDefinition == G4Proton::ProtonDefinition()
    221175      || particleDefinition == instance->GetIon("hydrogen"))
    222          
    223     {
     176  {
    224177      maximumKineticEnergyTransfer= 4.* (electron_mass_c2 / proton_mass_c2) * k;
    225     }
     178  }
    226179
    227180  if (particleDefinition == instance->GetIon("helium")
    228181      || particleDefinition == instance->GetIon("alpha+")
    229182      || particleDefinition == instance->GetIon("alpha++"))
    230     {
     183  {
    231184      maximumKineticEnergyTransfer= 4.* (0.511 / 3728) * k;
    232     }
     185  }
    233186
    234187  G4double crossSectionMaximum = 0.;
     188 
    235189  for(G4double value=waterStructure.IonisationEnergy(shell); value<=4.*waterStructure.IonisationEnergy(shell) ; value+=0.1*eV)
    236     {
     190  {
    237191      G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k, value, shell);
    238192      if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
    239     }
     193  }
     194 
    240195  G4double secElecKinetic = 0.;
    241   do{
     196 
     197  do
     198  {
    242199    secElecKinetic = G4UniformRand() * maximumKineticEnergyTransfer;
    243200  } while(G4UniformRand()*crossSectionMaximum > DifferentialCrossSection(particleDefinition,
    … …  
    249206}
    250207
     208//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     209
    251210
    252211void G4FinalStateIonisationRudd::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
    253212                                                                   G4double k,
    254213                                                                   G4double secKinetic,
    255                                                                    G4double cosTheta,
    256                                                                    G4double phi )
     214                                                                   G4double & cosTheta,
     215                                                                   G4double & phi )
    257216{
    258217  G4DNAGenericIonsManager *instance;
    … …  
    263222  if (particleDefinition == G4Proton::ProtonDefinition()
    264223      || particleDefinition == instance->GetIon("hydrogen"))
    265     {
     224  {
    266225      maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
    267     }
     226  }
    268227 
    269228  if (particleDefinition == instance->GetIon("helium")
    270229      || particleDefinition == instance->GetIon("alpha+")
    271230      || particleDefinition == instance->GetIon("alpha++"))
    272     {
     231  {
    273232      maxSecKinetic = 4.* (0.511 / 3728) * k;
    274     }
     233  }
    275234 
    276235  phi = twopi * G4UniformRand();
    277236  cosTheta = std::sqrt(secKinetic / maxSecKinetic);
    278237}
     238
     239//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    279240
    280241
    … …  
    314275
    315276  if (j == 4)
    316     {
     277  {
    317278      //Data For Liquid Water K SHELL from Dingfelder (Protons in Water)
    318279      A1 = 1.25;
    … …  
    326287      D2 = 0.00;
    327288      alphaConst = 0.66;
    328     }
     289  }
    329290  else
    330     {
     291  {
    331292      //Data For Liquid Water from Dingfelder (Protons in Water)
    332293      A1 = 1.02;
    … …  
    340301      D2 = 0.04;
    341302      alphaConst = 0.64;
    342     }
     303  }
    343304 
    344305  const G4double n = 2.;
    345306  const G4double Gj[5] = {0.99, 1.11, 1.11, 0.52, 1.};
    346 
    347   //const G4double I[5]={12.61*eV, 14.73*eV, 18.55*eV, 32.2*eV, 539.7*eV}; // for water Vapor
    348   //const G4double energyConstant[]={10.79*eV, 13.39*eV, 16.05*eV, 32.30*eV, 539.*eV};
    349307
    350308  G4DNAGenericIonsManager* instance;
    … …  
    359317  if (particleDefinition == G4Proton::ProtonDefinition()
    360318      || particleDefinition == instance->GetIon("hydrogen"))
    361     {
     319  {
    362320      tau = (electron_mass_c2/proton_mass_c2) * k ;
    363     }
     321  }
    364322   
    365323  if ( particleDefinition == instance->GetIon("helium")
    366324       || particleDefinition == instance->GetIon("alpha+")
    367325       || particleDefinition == instance->GetIon("alpha++"))
    368     {
     326  {
    369327      tau = (0.511/3728.) * k ;
    370     }
     328  }
    371329 
    372330  G4double S = 4.*pi * Bohr_radius*Bohr_radius * n * std::pow((Ry/waterStructure.IonisationEnergy(ionizationLevelIndex)),2);
    … …  
    390348          || particleDefinition == instance->GetIon("hydrogen")
    391349          )
    392     {
     350  {
    393351      return(sigma);
    394     }
    395 
    396   // ------------
    397  
     352  }
     353
    398354  if (particleDefinition == instance->GetIon("alpha++") )
    399     {
     355  {
    400356      slaterEffectiveCharge[0]=0.;
    401357      slaterEffectiveCharge[1]=0.;
    … …  
    404360      sCoefficient[1]=0.;
    405361      sCoefficient[2]=0.;
    406     }
     362  }
    407363
    408364  if (particleDefinition == instance->GetIon("alpha+") )
    409     {
     365  {
    410366      slaterEffectiveCharge[0]=2.0;
    411367      slaterEffectiveCharge[1]=1.15;
    … …  
    414370      sCoefficient[1]=0.15;
    415371      sCoefficient[2]=0.15;
    416     }
     372  }
    417373
    418374  if (particleDefinition == instance->GetIon("helium") )
    419     {
     375  {
    420376      slaterEffectiveCharge[0]=1.7;
    421377      slaterEffectiveCharge[1]=1.15;
    … …  
    424380      sCoefficient[1]=0.25;
    425381      sCoefficient[2]=0.25;
    426     }
     382  }
    427383 
    428384  if (    particleDefinition == instance->GetIon("helium")
    … …  
    430386          || particleDefinition == instance->GetIon("alpha++")
    431387          )
    432     {
     388  {
    433389      sigma = Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex)) * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
    434390   
    … …  
    440396           
    441397      return zEff * zEff * sigma ;
    442     } 
     398  } 
    443399 
    444400  return 0;
    445401}
     402
     403//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    446404
    447405G4double G4FinalStateIonisationRudd::S_1s(G4double t,
    … …  
    459417}
    460418
    461 
     419//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    462420
    463421G4double G4FinalStateIonisationRudd::S_2s(G4double t,
    … …  
    476434}
    477435
    478 
     436//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    479437
    480438G4double G4FinalStateIonisationRudd::S_2p(G4double t,
    … …  
    492450}
    493451
    494 
     452//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    495453
    496454G4double G4FinalStateIonisationRudd::R(G4double t,
    … …  
    500458{
    501459  // tElectron = m_electron / m_alpha * t
    502   // Hardcoded in Riccardo's implementation; to be corrected
    503460  // Dingfelder, in Chattanooga 2005 proceedings, p 4
    504461
    … …  
    509466}
    510467
    511 
     468//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    512469
    513470G4double G4FinalStateIonisationRudd::CorrectionFactor(G4ParticleDefinition* particleDefinition, G4double k)
    … …  
    517474
    518475  if (particleDefinition == G4Proton::Proton())
    519     {
     476  {
    520477      return(1.);
    521     }
     478  }
    522479  else
    523480    if (particleDefinition == instance->GetIon("hydrogen"))
    524       {
     481    {
    525482        G4double value = (std::log(k/eV)-4.2)/0.5;
    526483        return((0.8/(1+std::exp(value))) + 0.9);
    527       }
     484    }
    528485    else
    529       {   
     486    {   
    530487        return(1.);
    531       }
    532 }
     488    }
     489}
Note: See TracChangeset for help on using the changeset viewer.