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/G4FinalStateElasticBrennerZaider.cc

    r819 r961  
    2424// ********************************************************************
    2525//
    26 //
    27 // $Id: G4FinalStateElasticBrennerZaider.cc,v 1.1 2007/10/12 23:11:41 pia Exp $
    28 // GEANT4 tag $Name:  $
    29 //
    30 // Contact Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
    31 //
    32 // Reference: TNS Geant4-DNA paper
    33 // Reference for implementation model: NIM. 155, pp. 145-156, 1978
    34 
    35 // History:
    36 // -----------
    37 // Date         Name              Modification
    38 // 28 Apr 2007  M.G. Pia          Created in compliance with design described in TNS paper
    39 //
    40 // -------------------------------------------------------------------
    41 
    42 // Class description:
    43 // Reference: TNS Geant4-DNA paper
    44 // S. Chauvie et al., Geant4 physics processes for microdosimetry simulation:
    45 // design foundation and implementation of the first set of models,
    46 // IEEE Trans. Nucl. Sci., vol. 54, no. 6, Dec. 2007.
    47 // Further documentation available from http://www.ge.infn.it/geant4/dna
    48 
    49 // -------------------------------------------------------------------
    50 
     26// $Id: G4FinalStateElasticBrennerZaider.cc,v 1.8 2008/12/05 11:58:16 sincerti Exp $
     27// GEANT4 tag $Name: geant4-09-02-ref-02 $
    5128
    5229#include "G4FinalStateElasticBrennerZaider.hh"
    53 #include "G4Track.hh"
    54 #include "G4Step.hh"
    55 #include "G4DynamicParticle.hh"
    56 #include "Randomize.hh"
    5730
    58 #include "G4ParticleTypes.hh"
    59 #include "G4ParticleDefinition.hh"
    60 #include "G4Electron.hh"
    61 #include "G4SystemOfUnits.hh"
    62 #include "G4ParticleMomentum.hh"
     31//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    6332
    6433G4FinalStateElasticBrennerZaider::G4FinalStateElasticBrennerZaider()
    6534{
    66   // These data members will be used in the next implementation iteration,
    67   // when the enriched PhysicsModel policy is implemented
    68   name = "FinalStateElasticBrennerZaider";
    69   lowEnergyLimit = 7.4 * eV;
    70   highEnergyLimit = 10 * MeV;
     35  lowEnergyLimit = 8.23 * eV; // SI : i/o of 7.4 * eV;
     36  highEnergyLimit = 200 * eV;
    7137
    7238  betaCoeff.push_back(7.51525);
    … …  
    10066}
    10167
     68//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    10269
    10370G4FinalStateElasticBrennerZaider::~G4FinalStateElasticBrennerZaider()
    104 {
    105   // empty
    106   // G4DynamicParticle objects produced are owned by client
    107 }
    108  
     71{}
     72
     73//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    10974
    11075const G4FinalStateProduct& G4FinalStateElasticBrennerZaider::GenerateFinalState(const G4Track& track, const G4Step& /* step */)
    11176{
    112   // Clear previous secondaries, energy deposit and particle kill status
    11377  product.Clear();
    11478
    115   // Kinetic energy of primary particle
    11679  G4double k = track.GetDynamicParticle()->GetKineticEnergy();
     80 
     81  if ( k>= lowEnergyLimit && k<=highEnergyLimit )
     82  {
     83    G4double cosTheta = RandomizeCosTheta(k);
     84 
     85    G4double phi = 2. * pi * G4UniformRand();
    11786
    118   // Assume material = water; H2O number of electrons
    119   // ---- MGP ---- To be generalized later
    120   // const G4int z = 10;
     87    G4ThreeVector zVers = track.GetDynamicParticle()->GetMomentumDirection();
     88    G4ThreeVector xVers = zVers.orthogonal();
     89    G4ThreeVector yVers = zVers.cross(xVers);
    12190
    122   G4double cosTheta = RandomizeCosTheta(k);
     91    G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
     92    G4double yDir = xDir;
     93    xDir *= std::cos(phi);
     94    yDir *= std::sin(phi);
     95
     96    G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
     97
     98    product.ModifyPrimaryParticle(zPrimeVers,k);
     99  }
     100
     101  if (k<lowEnergyLimit)
     102  {
     103    product.KillPrimaryParticle();
     104  }
    123105 
    124   G4double phi = 2. * pi * G4UniformRand();
    125 
    126   // G4cout << "cosTheta in GenerateFinalState = " << cosTheta << ", phi = " << phi << G4endl;
    127 
    128   G4ThreeVector zVers = track.GetDynamicParticle()->GetMomentumDirection();
    129   G4ThreeVector xVers = zVers.orthogonal();
    130   G4ThreeVector yVers = zVers.cross(xVers);
    131 
    132   G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
    133   G4double yDir = xDir;
    134   xDir *= std::cos(phi);
    135   yDir *= std::sin(phi);
    136 
    137   // G4cout << "xDir, yDir = " << xDir <<", " << yDir << G4endl;
    138 
    139   // G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers).unit());
    140   G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
    141 
    142   // G4cout << "zPrimeVers = (" << zPrimeVers.x() << ", "<< zPrimeVers.y() << ", "<< zPrimeVers.z() << ") " << G4endl;
    143 
    144   //  product.ModifyPrimaryParticle(zPrimeVers.x(),zPrimeVers.y(),zPrimeVers.z(),k);
    145   product.ModifyPrimaryParticle(zPrimeVers,k);
    146 
    147   //  this->aParticleChange.ProposeEnergy(k);
    148   //  this->aParticleChange.ProposeMomentumDirection(zPrimeVers);
    149   //  this->aParticleChange.SetNumberOfSecondaries(0);
    150 
    151106  return product;
    152107}
     108
     109//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    153110
    154111G4double G4FinalStateElasticBrennerZaider::RandomizeCosTheta(G4double k)
    … …  
    158115  //   d Omega           (1 + 2 gamma(K) - cos(theta))^2     (1 + 2 delta(K) + cos(theta))^2
    159116  //
    160   // Maximum is < 1/(4 gamma(K)^2) + beta(K)/(4 delta(K)^2)
     117  // Maximum is < 1/(4 gamma(K)^2) + beta(K)/((2+2delta(K))^2)
    161118  //
    162119  // Phys. Med. Biol. 29 N.4 (1983) 443-447
    163120 
    164121  // gamma(K), beta(K) and delta(K) are polynomials with coefficients for energy measured in eV
     122
    165123  k /= eV;
    166124 
    167125  G4double beta = std::exp(CalculatePolynomial(k,betaCoeff));
    168126  G4double delta = std::exp(CalculatePolynomial(k,deltaCoeff));
     127  G4double gamma;
    169128 
    170   G4double gamma;
    171129  if (k > 100.)
    172     {
    173       gamma = CalculatePolynomial(k, gamma100_200Coeff); // Only in this case it is not the exponent of the polynomial
    174     } 
     130  {
     131      gamma = CalculatePolynomial(k, gamma100_200Coeff);
     132      // Only in this case it is not the exponent of the polynomial
     133  } 
    175134  else
    176     {
     135  {
     136      if (k>10)
     137      {
     138          gamma = std::exp(CalculatePolynomial(k, gamma10_100Coeff));
     139      }
     140      else
     141      {
     142          gamma = std::exp(CalculatePolynomial(k, gamma035_10Coeff));
     143      }
     144  }
    177145
    178       if (k>10)
    179         {
    180           gamma = std::exp(CalculatePolynomial(k, gamma10_100Coeff));
    181         }
    182       else
    183         {
    184           gamma = std::exp(CalculatePolynomial(k, gamma035_10Coeff));
    185         }
    186     }
    187  
    188   // G4cout << "beta = " << beta << ", gamma = " << gamma << ", delta = " << delta << G4endl;
    189 
    190   G4double oneOverMax = 1. / (1./(4.*gamma*gamma) + beta/(4.*delta*delta));
     146  G4double oneOverMax = 1. / (1./(4.*gamma*gamma) + beta/( (2.+2.*delta)*(2.+2.*delta) ));
    191147 
    192148  G4double cosTheta = 0.;
    … …  
    196152 
    197153  do
    198     {
     154  {
    199155      cosTheta = 2. * G4UniformRand() - 1.;
    200       leftDenominator = (1 + 2.*gamma - cosTheta);
    201       rightDenominator = (1 + 2.*delta + cosTheta);
    202       fCosTheta = oneOverMax * (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator));
    203     }
     156      leftDenominator = (1. + 2.*gamma - cosTheta);
     157      rightDenominator = (1. + 2.*delta + cosTheta);
     158      if ( (leftDenominator * rightDenominator) != 0. )
     159      {
     160          fCosTheta = oneOverMax * (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator));
     161      }
     162  }
    204163  while (fCosTheta < G4UniformRand());
    205164
    206   //  G4cout << "cosTheta = " << cosTheta << G4endl;
    207  
    208165  return cosTheta;
    209166}
     167
     168//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    210169
    211170G4double G4FinalStateElasticBrennerZaider::CalculatePolynomial(G4double k, std::vector<G4double>& vec)
Note: See TracChangeset for help on using the changeset viewer.