G4DNASancheSolvatationModel.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id: G4DNASancheSolvatationModel.cc 64057 2012-10-30 15:04:49Z gcosmo $
00027 //
00028 // Author: Mathieu Karamitros (kara (AT) cenbg . in2p3 . fr) 
00029 //
00030 // WARNING : This class is released as a prototype.
00031 // It might strongly evolve or even disapear in the next releases.
00032 //
00033 // History:
00034 // -----------
00035 // 10 Oct 2011 M.Karamitros created
00036 //
00037 // -------------------------------------------------------------------
00038 
00039 #include "G4DNASancheSolvatationModel.hh"
00040 #include "G4PhysicalConstants.hh"
00041 #include "G4SystemOfUnits.hh"
00042 #include "G4DNAWaterExcitationStructure.hh"
00043 #include "G4ParticleChangeForGamma.hh"
00044 #include "G4Electron.hh"
00045 #include "G4NistManager.hh"
00046 #include "G4DNAChemistryManager.hh"
00047 #include "G4DNAMolecularMaterial.hh"
00048 
00049 G4DNASancheSolvatationModel::G4DNASancheSolvatationModel(const G4ParticleDefinition*,
00050                                                          const G4String& nam):
00051     G4VEmModel(nam),fIsInitialised(false)
00052 {
00053     fVerboseLevel = 0 ;
00054     SetLowEnergyLimit(0.);
00055     G4DNAWaterExcitationStructure exStructure ;
00056     SetHighEnergyLimit(exStructure.ExcitationEnergy(0));
00057     fParticleChangeForGamma = 0;
00058   //  fNistWater  = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
00059     fpWaterDensity = 0;
00060 }
00061 
00062 //______________________________________________________________________
00063 G4DNASancheSolvatationModel::~G4DNASancheSolvatationModel()
00064 {}
00065 
00066 //______________________________________________________________________
00067 void G4DNASancheSolvatationModel::Initialise(const G4ParticleDefinition* particleDefinition,
00068                                              const G4DataVector&)
00069 {
00070 #ifdef G4VERBOSE
00071     if (fVerboseLevel)
00072         G4cout << "Calling G4SancheSolvatationModel::Initialise()" << G4endl;
00073 #endif
00074     if (particleDefinition != G4Electron::ElectronDefinition())
00075     {
00076         G4ExceptionDescription exceptionDescription ;
00077         exceptionDescription << "Attempting to calculate cross section for wrong particle";
00078         G4Exception("G4DNASancheSolvatationModel::CrossSectionPerVolume","G4DNASancheSolvatationModel001",
00079                     FatalErrorInArgument,exceptionDescription);
00080         return;
00081     }
00082 
00083     if(!fIsInitialised)
00084     {
00085         fIsInitialised = true;
00086         fParticleChangeForGamma = GetParticleChangeForGamma();
00087     }
00088 
00089     // Initialize water density pointer
00090     fpWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
00091 }
00092 
00093 //______________________________________________________________________
00094 G4double G4DNASancheSolvatationModel::CrossSectionPerVolume(const G4Material* material,
00095                                                             const G4ParticleDefinition*,
00096                                                             G4double ekin,
00097                                                             G4double,
00098                                                             G4double)
00099 {
00100 #ifdef G4VERBOSE
00101     if (fVerboseLevel > 1)
00102         G4cout << "Calling CrossSectionPerVolume() of G4SancheSolvatationModel" << G4endl;
00103 #endif
00104 
00105     if(ekin > HighEnergyLimit())
00106     {
00107         return 0.0;
00108     }
00109 
00110     G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
00111 
00112     if(waterDensity!= 0.0)
00113    //  if (material == nistwater || material->GetBaseMaterial() == nistwater)
00114     {
00115         if (ekin <= HighEnergyLimit())
00116         {
00117             return DBL_MAX;
00118         }
00119     }
00120     return 0. ;
00121 }
00122 
00123 //______________________________________________________________________
00124 G4ThreeVector G4DNASancheSolvatationModel::RadialDistributionOfProducts(G4double expectationValue) const
00125 {
00126     G4double sigma = std::sqrt(1.57)/2*expectationValue;
00127 
00128     G4double XValueForfMax = std::sqrt(2.*sigma*sigma);
00129     G4double fMaxValue = std::sqrt(2./3.14) * 1./(sigma*sigma*sigma) *
00130             (XValueForfMax*XValueForfMax)*
00131             std::exp(-1./2. * (XValueForfMax*XValueForfMax)/(sigma*sigma));
00132 
00133     G4double R;
00134 
00135     do
00136     {
00137         G4double aRandomfValue = fMaxValue*G4UniformRand();
00138 
00139         G4double sign;
00140         if(G4UniformRand() > 0.5)
00141         {
00142             sign = +1.;
00143         }
00144         else
00145         {
00146             sign = -1;
00147         }
00148 
00149         R = expectationValue + sign*3.*sigma* G4UniformRand();
00150         G4double f = std::sqrt(2./3.14) * 1/std::pow(sigma, 3) * R*R * std::exp(-1./2. * R*R/(sigma*sigma));
00151 
00152         if(aRandomfValue < f)
00153         {
00154             break;
00155         }
00156     }
00157     while(1);
00158 
00159     G4double costheta = (2.*G4UniformRand()-1.);
00160     G4double theta = std::acos (costheta);
00161     G4double phi = 2.*pi*G4UniformRand();
00162 
00163     G4double xDirection = R*std::cos(phi)* std::sin(theta);
00164     G4double yDirection = R*std::sin(theta)*std::sin(phi);
00165     G4double zDirection = R*costheta;
00166     G4ThreeVector RandDirection = G4ThreeVector(xDirection, yDirection, zDirection);
00167 
00168     return RandDirection;
00169 }
00170 
00171 //______________________________________________________________________
00172 void G4DNASancheSolvatationModel::SampleSecondaries(std::vector<G4DynamicParticle*>*,
00173                                                     const G4MaterialCutsCouple*,
00174                                                     const G4DynamicParticle* particle,
00175                                                     G4double,
00176                                                     G4double)
00177 {
00178 #ifdef G4VERBOSE
00179     if (fVerboseLevel)
00180         G4cout << "Calling SampleSecondaries() of G4SancheSolvatationModel" << G4endl;
00181 #endif
00182     G4double k = particle->GetKineticEnergy();
00183 
00184     if (k <= HighEnergyLimit())
00185     {
00186         G4double r_mean =
00187                 (-0.003*std::pow(k/eV,6) + 0.0749*std::pow(k/eV,5) - 0.7197*std::pow(k/eV,4)
00188                  + 3.1384*std::pow(k/eV,3) - 5.6926*std::pow(k/eV,2) + 5.6237*k/eV - 0.7883)*nanometer;
00189 
00190         G4ThreeVector displacement = RadialDistributionOfProducts (r_mean);
00191 
00192         //______________________________________________________________
00193         const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
00194         G4ThreeVector finalPosition(theIncomingTrack->GetPosition()+displacement);
00195 
00196         G4DNAChemistryManager::Instance()->CreateSolvatedElectron(theIncomingTrack,&finalPosition);
00197 
00198         fParticleChangeForGamma->SetProposedKineticEnergy(25.e-3*eV);
00199         fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
00200         fParticleChangeForGamma->ProposeLocalEnergyDeposit(k);
00201     }
00202 }

Generated on Mon May 27 17:48:04 2013 for Geant4 by  doxygen 1.4.7