Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNASancheSolvatationModel.cc
Go to the documentation of this file.
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: G4DNASancheSolvatationModel.cc 64057 2012-10-30 15:04:49Z gcosmo $
27 //
28 // Author: Mathieu Karamitros (kara (AT) cenbg . in2p3 . fr)
29 //
30 // WARNING : This class is released as a prototype.
31 // It might strongly evolve or even disapear in the next releases.
32 //
33 // History:
34 // -----------
35 // 10 Oct 2011 M.Karamitros created
36 //
37 // -------------------------------------------------------------------
38 
40 #include "G4PhysicalConstants.hh"
41 #include "G4SystemOfUnits.hh"
44 #include "G4Electron.hh"
45 #include "G4NistManager.hh"
46 #include "G4DNAChemistryManager.hh"
48 
50  const G4String& nam):
51  G4VEmModel(nam),fIsInitialised(false)
52 {
53  fVerboseLevel = 0 ;
55  G4DNAWaterExcitationStructure exStructure ;
56  SetHighEnergyLimit(exStructure.ExcitationEnergy(0));
58  // fNistWater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
59  fpWaterDensity = 0;
60 }
61 
62 //______________________________________________________________________
64 {}
65 
66 //______________________________________________________________________
68  const G4DataVector&)
69 {
70 #ifdef G4VERBOSE
71  if (fVerboseLevel)
72  G4cout << "Calling G4SancheSolvatationModel::Initialise()" << G4endl;
73 #endif
74  if (particleDefinition != G4Electron::ElectronDefinition())
75  {
76  G4ExceptionDescription exceptionDescription ;
77  exceptionDescription << "Attempting to calculate cross section for wrong particle";
78  G4Exception("G4DNASancheSolvatationModel::CrossSectionPerVolume","G4DNASancheSolvatationModel001",
79  FatalErrorInArgument,exceptionDescription);
80  return;
81  }
82 
83  if(!fIsInitialised)
84  {
85  fIsInitialised = true;
87  }
88 
89  // Initialize water density pointer
91 }
92 
93 //______________________________________________________________________
95  const G4ParticleDefinition*,
96  G4double ekin,
97  G4double,
98  G4double)
99 {
100 #ifdef G4VERBOSE
101  if (fVerboseLevel > 1)
102  G4cout << "Calling CrossSectionPerVolume() of G4SancheSolvatationModel" << G4endl;
103 #endif
104 
105  if(ekin > HighEnergyLimit())
106  {
107  return 0.0;
108  }
109 
110  G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
111 
112  if(waterDensity!= 0.0)
113  // if (material == nistwater || material->GetBaseMaterial() == nistwater)
114  {
115  if (ekin <= HighEnergyLimit())
116  {
117  return DBL_MAX;
118  }
119  }
120  return 0. ;
121 }
122 
123 //______________________________________________________________________
125 {
126  G4double sigma = std::sqrt(1.57)/2*expectationValue;
127 
128  G4double XValueForfMax = std::sqrt(2.*sigma*sigma);
129  G4double fMaxValue = std::sqrt(2./3.14) * 1./(sigma*sigma*sigma) *
130  (XValueForfMax*XValueForfMax)*
131  std::exp(-1./2. * (XValueForfMax*XValueForfMax)/(sigma*sigma));
132 
133  G4double R;
134 
135  do
136  {
137  G4double aRandomfValue = fMaxValue*G4UniformRand();
138 
139  G4double sign;
140  if(G4UniformRand() > 0.5)
141  {
142  sign = +1.;
143  }
144  else
145  {
146  sign = -1;
147  }
148 
149  R = expectationValue + sign*3.*sigma* G4UniformRand();
150  G4double f = std::sqrt(2./3.14) * 1/std::pow(sigma, 3) * R*R * std::exp(-1./2. * R*R/(sigma*sigma));
151 
152  if(aRandomfValue < f)
153  {
154  break;
155  }
156  }
157  while(1);
158 
159  G4double costheta = (2.*G4UniformRand()-1.);
160  G4double theta = std::acos (costheta);
161  G4double phi = 2.*pi*G4UniformRand();
162 
163  G4double xDirection = R*std::cos(phi)* std::sin(theta);
164  G4double yDirection = R*std::sin(theta)*std::sin(phi);
165  G4double zDirection = R*costheta;
166  G4ThreeVector RandDirection = G4ThreeVector(xDirection, yDirection, zDirection);
167 
168  return RandDirection;
169 }
170 
171 //______________________________________________________________________
172 void G4DNASancheSolvatationModel::SampleSecondaries(std::vector<G4DynamicParticle*>*,
173  const G4MaterialCutsCouple*,
174  const G4DynamicParticle* particle,
175  G4double,
176  G4double)
177 {
178 #ifdef G4VERBOSE
179  if (fVerboseLevel)
180  G4cout << "Calling SampleSecondaries() of G4SancheSolvatationModel" << G4endl;
181 #endif
182  G4double k = particle->GetKineticEnergy();
183 
184  if (k <= HighEnergyLimit())
185  {
186  G4double r_mean =
187  (-0.003*std::pow(k/eV,6) + 0.0749*std::pow(k/eV,5) - 0.7197*std::pow(k/eV,4)
188  + 3.1384*std::pow(k/eV,3) - 5.6926*std::pow(k/eV,2) + 5.6237*k/eV - 0.7883)*nanometer;
189 
190  G4ThreeVector displacement = RadialDistributionOfProducts (r_mean);
191 
192  //______________________________________________________________
193  const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
194  G4ThreeVector finalPosition(theIncomingTrack->GetPosition()+displacement);
195 
196  G4DNAChemistryManager::Instance()->CreateSolvatedElectron(theIncomingTrack,&finalPosition);
197 
201  }
202 }