Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/dna/models/src/G4DNAScreenedRutherfordElasticModel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/dna/models/src/G4DNAScreenedRutherfordElasticModel.cc (Version 11.3.0) and /processes/electromagnetic/dna/models/src/G4DNAScreenedRutherfordElasticModel.cc (Version 10.6.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 //                                                 26 //
 27                                                    27 
 28 #include "G4DNAScreenedRutherfordElasticModel.     28 #include "G4DNAScreenedRutherfordElasticModel.hh"
 29 #include "G4PhysicalConstants.hh"                  29 #include "G4PhysicalConstants.hh"
 30 #include "G4SystemOfUnits.hh"                      30 #include "G4SystemOfUnits.hh"
 31 #include "G4DNAMolecularMaterial.hh"               31 #include "G4DNAMolecularMaterial.hh"
 32 #include "G4Exp.hh"                                32 #include "G4Exp.hh"
 33 #include "G4Log.hh"                                33 #include "G4Log.hh"
 34                                                    34 
 35 //....oooOO0OOooo........oooOO0OOooo........oo     35 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 36                                                    36 
 37 using namespace std;                               37 using namespace std;
 38                                                    38 
 39 //....oooOO0OOooo........oooOO0OOooo........oo     39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 40                                                    40 
 41 G4DNAScreenedRutherfordElasticModel::              41 G4DNAScreenedRutherfordElasticModel::
 42 G4DNAScreenedRutherfordElasticModel(const G4Pa     42 G4DNAScreenedRutherfordElasticModel(const G4ParticleDefinition*,
 43                                     const G4St     43                                     const G4String& nam) :
 44 G4VEmModel(nam)                                <<  44 G4VEmModel(nam), isInitialised(false)
 45 {                                                  45 {
 46   fpWaterDensity = nullptr;                    <<  46   fpWaterDensity = 0;
 47                                                    47 
 48   lowEnergyLimit = 0 * eV;                         48   lowEnergyLimit = 0 * eV;
 49   intermediateEnergyLimit = 200 * eV; // Switc     49   intermediateEnergyLimit = 200 * eV; // Switch between two final state models
 50   highEnergyLimit = 1. * MeV;                      50   highEnergyLimit = 1. * MeV;
 51                                                    51   
 52   SetLowEnergyLimit(lowEnergyLimit);               52   SetLowEnergyLimit(lowEnergyLimit);
 53   SetHighEnergyLimit(highEnergyLimit);             53   SetHighEnergyLimit(highEnergyLimit);
 54                                                    54 
 55   verboseLevel = 0;                                55   verboseLevel = 0;
 56   // Verbosity scale:                              56   // Verbosity scale:
 57   // 0 = nothing                                   57   // 0 = nothing
 58   // 1 = warning for energy non-conservation       58   // 1 = warning for energy non-conservation
 59   // 2 = details of energy budget                  59   // 2 = details of energy budget
 60   // 3 = calculation of cross sections, file o     60   // 3 = calculation of cross sections, file openings, sampling of atoms
 61   // 4 = entering in methods                       61   // 4 = entering in methods
 62                                                    62 
 63 #ifdef SR_VERBOSE                                  63 #ifdef SR_VERBOSE
 64   if (verboseLevel > 0)                            64   if (verboseLevel > 0)
 65   {                                                65   {
 66     G4cout << "Screened Rutherford Elastic mod     66     G4cout << "Screened Rutherford Elastic model is constructed "
 67            << G4endl                               67            << G4endl
 68            << "Energy range: "                     68            << "Energy range: "
 69            << lowEnergyLimit / eV << " eV - "      69            << lowEnergyLimit / eV << " eV - "
 70            << highEnergyLimit / MeV << " MeV"      70            << highEnergyLimit / MeV << " MeV"
 71            << G4endl;                              71            << G4endl;
 72   }                                                72   }
 73 #endif                                             73 #endif
 74   fParticleChangeForGamma = nullptr;           <<  74   fParticleChangeForGamma = 0;
 75                                                    75 
 76   // Selection of computation method               76   // Selection of computation method
 77   // We do not recommend "true" usage with the     77   // We do not recommend "true" usage with the current cumul. proba. settings
 78   fasterCode = false;                              78   fasterCode = false; 
 79 }                                                  79 }
 80                                                    80 
 81 //....oooOO0OOooo........oooOO0OOooo........oo     81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 82                                                    82 
 83 G4DNAScreenedRutherfordElasticModel::~G4DNAScr     83 G4DNAScreenedRutherfordElasticModel::~G4DNAScreenedRutherfordElasticModel()
 84 = default;                                     <<  84 {
                                                   >>  85 }
 85                                                    86 
 86 //....oooOO0OOooo........oooOO0OOooo........oo     87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 87                                                    88 
 88 void G4DNAScreenedRutherfordElasticModel::         89 void G4DNAScreenedRutherfordElasticModel::
 89 Initialise(const G4ParticleDefinition* particl     90 Initialise(const G4ParticleDefinition* particle,
 90            const G4DataVector& /*cuts*/)           91            const G4DataVector& /*cuts*/)
 91 {                                                  92 {
 92 #ifdef SR_VERBOSE                                  93 #ifdef SR_VERBOSE
 93   if (verboseLevel > 3)                            94   if (verboseLevel > 3)
 94   {                                                95   {
 95     G4cout << "Calling G4DNAScreenedRutherford     96     G4cout << "Calling G4DNAScreenedRutherfordElasticModel::Initialise()"
 96            << G4endl;                              97            << G4endl;
 97   }                                                98   }
 98 #endif                                             99 #endif
 99                                                   100   
100   if(particle->GetParticleName() != "e-")         101   if(particle->GetParticleName() != "e-")
101   {                                               102   {
102     G4Exception ("*** WARNING: the G4DNAScreen    103     G4Exception ("*** WARNING: the G4DNAScreenedRutherfordElasticModel is not "
103                  "intented to be used with ano    104                  "intented to be used with another particle than the electron",
104                  "",FatalException,"") ;          105                  "",FatalException,"") ;
105   }                                               106   }
106                                                   107   
107   // Energy limits                                108   // Energy limits
108   if (LowEnergyLimit() < 9*eV)                    109   if (LowEnergyLimit() < 9*eV)
109   {                                               110   {
110     G4Exception("*** WARNING: the G4DNAScreene    111     G4Exception("*** WARNING: the G4DNAScreenedRutherfordElasticModel class is "
111                 "not validated below 9 eV",       112                 "not validated below 9 eV",
112                 "",JustWarning,"") ;              113                 "",JustWarning,"") ;
113   }                                               114   }
114                                                   115 
115   if (HighEnergyLimit() > 1*MeV)                  116   if (HighEnergyLimit() > 1*MeV)
116   {                                               117   {
117     G4Exception("*** WARNING: the G4DNAScreene    118     G4Exception("*** WARNING: the G4DNAScreenedRutherfordElasticModel class is "
118                 "not validated above 1 MeV",      119                 "not validated above 1 MeV",
119                 "",JustWarning,"") ;              120                 "",JustWarning,"") ;
120   }                                               121   }
121                                                   122 
122   //                                              123   //
123 #ifdef SR_VERBOSE                                 124 #ifdef SR_VERBOSE
124   if( verboseLevel>0 )                            125   if( verboseLevel>0 )
125   {                                               126   {
126     G4cout << "Screened Rutherford elastic mod    127     G4cout << "Screened Rutherford elastic model is initialized " << G4endl
127            << "Energy range: "                    128            << "Energy range: "
128            << LowEnergyLimit() / eV << " eV -     129            << LowEnergyLimit() / eV << " eV - "
129            << HighEnergyLimit() / MeV << " MeV    130            << HighEnergyLimit() / MeV << " MeV"
130            << G4endl;                             131            << G4endl;
131   }                                               132   }
132 #endif                                            133 #endif
133                                                   134 
134   if (isInitialised) { return; } // return her    135   if (isInitialised) { return; } // return here, prevent reinit consts + pointer
135                                                   136   
136   // Initialize water density pointer             137   // Initialize water density pointer
137   fpWaterDensity = G4DNAMolecularMaterial::Ins    138   fpWaterDensity = G4DNAMolecularMaterial::Instance()->
138                   GetNumMolPerVolTableFor(G4Ma    139                   GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
139                                                   140   
140   fParticleChangeForGamma = GetParticleChangeF    141   fParticleChangeForGamma = GetParticleChangeForGamma();
141   isInitialised = true;                           142   isInitialised = true;
142                                                   143   
143   // Constants for final state by Brenner & Za    144   // Constants for final state by Brenner & Zaider
144   // note: if called after if(isInitialised) n    145   // note: if called after if(isInitialised) no need for clear and resetting
145   // the values at every call                     146   // the values at every call
146                                                   147 
147   betaCoeff=                                      148   betaCoeff=
148   {                                               149   {
149     7.51525,                                      150     7.51525,
150     -0.41912,                                     151     -0.41912,
151     7.2017E-3,                                    152     7.2017E-3,
152     -4.646E-5,                                    153     -4.646E-5,
153     1.02897E-7};                                  154     1.02897E-7};
154                                                   155 
155   deltaCoeff=                                     156   deltaCoeff=
156   {                                               157   {
157     2.9612,                                       158     2.9612,
158     -0.26376,                                     159     -0.26376,
159     4.307E-3,                                     160     4.307E-3,
160     -2.6895E-5,                                   161     -2.6895E-5,
161     5.83505E-8};                                  162     5.83505E-8};
162                                                   163 
163   gamma035_10Coeff =                              164   gamma035_10Coeff =
164   {                                               165   {
165     -1.7013,                                      166     -1.7013,
166     -1.48284,                                     167     -1.48284,
167     0.6331,                                       168     0.6331,
168     -0.10911,                                     169     -0.10911,
169     8.358E-3,                                     170     8.358E-3,
170     -2.388E-4};                                   171     -2.388E-4};
171                                                   172 
172   gamma10_100Coeff =                              173   gamma10_100Coeff =
173   {                                               174   {
174     -3.32517,                                     175     -3.32517,
175     0.10996,                                      176     0.10996,
176     -4.5255E-3,                                   177     -4.5255E-3,
177     5.8372E-5,                                    178     5.8372E-5,
178     -2.4659E-7};                                  179     -2.4659E-7};
179                                                   180 
180   gamma100_200Coeff =                             181   gamma100_200Coeff =
181   {                                               182   {
182     2.4775E-2,                                    183     2.4775E-2,
183     -2.96264E-5,                                  184     -2.96264E-5,
184     -1.20655E-7};                                 185     -1.20655E-7};
185 }                                                 186 }
186                                                   187 
187 //....oooOO0OOooo........oooOO0OOooo........oo    188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
188                                                   189 
189 G4double G4DNAScreenedRutherfordElasticModel::    190 G4double G4DNAScreenedRutherfordElasticModel::
190 CrossSectionPerVolume(const G4Material* materi    191 CrossSectionPerVolume(const G4Material* material,
191 #ifdef SR_VERBOSE                                 192 #ifdef SR_VERBOSE
192                       const G4ParticleDefiniti    193                       const G4ParticleDefinition* particleDefinition,
193 #else                                             194 #else
194                       const G4ParticleDefiniti    195                       const G4ParticleDefinition*,
195 #endif                                            196 #endif
196                       G4double ekin,              197                       G4double ekin,
197                       G4double,                   198                       G4double,
198                       G4double)                   199                       G4double)
199 {                                                 200 {
200 #ifdef SR_VERBOSE                                 201 #ifdef SR_VERBOSE
201   if (verboseLevel > 3)                           202   if (verboseLevel > 3)
202   {                                               203   {
203     G4cout << "Calling CrossSectionPerVolume()    204     G4cout << "Calling CrossSectionPerVolume() of "
204             "G4DNAScreenedRutherfordElasticMod    205             "G4DNAScreenedRutherfordElasticModel"
205            << G4endl;                             206            << G4endl;
206   }                                               207   }
207 #endif                                            208 #endif
208                                                   209   
209   // Calculate total cross section for model      210   // Calculate total cross section for model
210                                                   211 
211   G4double sigma=0.;                              212   G4double sigma=0.;
212   G4double waterDensity = (*fpWaterDensity)[ma    213   G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
213                                                   214 
214   if(ekin <= HighEnergyLimit() && ekin >= LowE    215   if(ekin <= HighEnergyLimit() && ekin >= LowEnergyLimit())
215   {                                               216   {
216     G4double z = 10.;                             217     G4double z = 10.;
217     G4double n = ScreeningFactor(ekin,z);         218     G4double n = ScreeningFactor(ekin,z);
218     G4double crossSection = RutherfordCrossSec    219     G4double crossSection = RutherfordCrossSection(ekin, z);
219     sigma = pi * crossSection / (n * (n + 1.))    220     sigma = pi * crossSection / (n * (n + 1.));
220   }                                               221   }
221                                                   222 
222 #ifdef SR_VERBOSE                                 223 #ifdef SR_VERBOSE 
223   if (verboseLevel > 2)                           224   if (verboseLevel > 2)
224   {                                               225   {
225     G4cout << "_______________________________    226     G4cout << "__________________________________" << G4endl;
226     G4cout << "=== G4DNAScreenedRutherfordElas    227     G4cout << "=== G4DNAScreenedRutherfordElasticModel - XS INFO START"
227            << G4endl;                             228            << G4endl;
228     G4cout << "=== Kinetic energy(eV)=" << eki    229     G4cout << "=== Kinetic energy(eV)=" << ekin/eV
229            << " particle : " << particleDefini    230            << " particle : " << particleDefinition->GetParticleName()
230            << G4endl;                             231            << G4endl;
231     G4cout << "=== Cross section per water mol    232     G4cout << "=== Cross section per water molecule (cm^2)=" << sigma/cm/cm
232            << G4endl;                             233            << G4endl;
233     G4cout << "=== Cross section per water mol    234     G4cout << "=== Cross section per water molecule (cm^-1)="
234            << sigma*waterDensity/(1./cm) << G4    235            << sigma*waterDensity/(1./cm) << G4endl;
235     G4cout << "=== G4DNAScreenedRutherfordElas    236     G4cout << "=== G4DNAScreenedRutherfordElasticModel - XS INFO END"
236            << G4endl;                             237            << G4endl;
237   }                                               238   }
238 #endif                                            239 #endif
239                                                   240 
240   return sigma*waterDensity;                      241   return sigma*waterDensity;
241 }                                                 242 }
242                                                   243 
243 //....oooOO0OOooo........oooOO0OOooo........oo    244 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
244                                                   245 
245 G4double G4DNAScreenedRutherfordElasticModel::    246 G4double G4DNAScreenedRutherfordElasticModel::RutherfordCrossSection(G4double k,
246                                                   247                                                                      G4double z)
247 {                                                 248 {
248   //                                              249   //
249   //                               e^4            250   //                               e^4         /      K + m_e c^2      \^2
250   // sigma_Ruth(K) = Z (Z+1) -----------------    251   // sigma_Ruth(K) = Z (Z+1) -------------------- | --------------------- |
251   //                          (4 pi epsilon_0)    252   //                          (4 pi epsilon_0)^2  \  K * (K + 2 m_e c^2)  /
252   //                                              253   //
253   // Where K is the electron non-relativistic     254   // Where K is the electron non-relativistic kinetic energy
254   //                                              255   //
255   // NIM 155, pp. 145-156, 1978                   256   // NIM 155, pp. 145-156, 1978
256                                                   257 
257   G4double length = (e_squared * (k + electron    258   G4double length = (e_squared * (k + electron_mass_c2))
258       / (4 * pi * epsilon0 * k * (k + 2 * elec    259       / (4 * pi * epsilon0 * k * (k + 2 * electron_mass_c2));
259   G4double cross = z * (z + 1) * length * leng    260   G4double cross = z * (z + 1) * length * length;
260                                                   261 
261   return cross;                                   262   return cross;
262 }                                                 263 }
263                                                   264 
264 //....oooOO0OOooo........oooOO0OOooo........oo    265 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
265                                                   266 
266 G4double G4DNAScreenedRutherfordElasticModel::    267 G4double G4DNAScreenedRutherfordElasticModel::ScreeningFactor(G4double k,
267                                                   268                                                               G4double z)
268 {                                                 269 {
269   //                                              270   //
270   //         alpha_1 + beta_1 ln(K/eV)   const    271   //         alpha_1 + beta_1 ln(K/eV)   constK Z^(2/3)
271   // n(T) = -------------------------- -------    272   // n(T) = -------------------------- -----------------
272   //              K/(m_e c^2)            2 + K    273   //              K/(m_e c^2)            2 + K/(m_e c^2)
273   //                                              274   //
274   // Where K is the electron non-relativistic     275   // Where K is the electron non-relativistic kinetic energy
275   //                                              276   //
276   // n(T) > 0 for T < ~ 400 MeV                   277   // n(T) > 0 for T < ~ 400 MeV
277   //                                              278   //
278   // NIM 155, pp. 145-156, 1978                   279   // NIM 155, pp. 145-156, 1978
279   // Formulae (2) and (5)                         280   // Formulae (2) and (5)
280                                                   281 
281   const G4double alpha_1(1.64);                   282   const G4double alpha_1(1.64);
282   const G4double beta_1(-0.0825);                 283   const G4double beta_1(-0.0825);
283   const G4double constK(1.7E-5);                  284   const G4double constK(1.7E-5);
284                                                   285 
285   G4double numerator = (alpha_1 + beta_1 * G4L    286   G4double numerator = (alpha_1 + beta_1 * G4Log(k / eV)) * constK
286                        * std::pow(z, 2. / 3.);    287                        * std::pow(z, 2. / 3.);
287                                                   288 
288   k /= electron_mass_c2;                          289   k /= electron_mass_c2;
289                                                   290 
290   G4double denominator = k * (2 + k);             291   G4double denominator = k * (2 + k);
291                                                   292 
292   G4double value = 0.;                            293   G4double value = 0.;
293   if (denominator > 0.) value = numerator / de    294   if (denominator > 0.) value = numerator / denominator;
294                                                   295 
295   return value;                                   296   return value;
296 }                                                 297 }
297                                                   298 
298 //....oooOO0OOooo........oooOO0OOooo........oo    299 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
299                                                   300 
300 void G4DNAScreenedRutherfordElasticModel::        301 void G4DNAScreenedRutherfordElasticModel::
301 SampleSecondaries(std::vector<G4DynamicParticl    302 SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
302                   const G4MaterialCutsCouple*     303                   const G4MaterialCutsCouple* /*couple*/,
303                   const G4DynamicParticle* aDy    304                   const G4DynamicParticle* aDynamicElectron,
304                   G4double,                       305                   G4double,
305                   G4double)                       306                   G4double)
306 {                                                 307 {
307 #ifdef SR_VERBOSE                                 308 #ifdef SR_VERBOSE
308   if (verboseLevel > 3)                           309   if (verboseLevel > 3)
309   {                                               310   {
310     G4cout << "Calling SampleSecondaries() of     311     G4cout << "Calling SampleSecondaries() of "
311               "G4DNAScreenedRutherfordElasticM    312               "G4DNAScreenedRutherfordElasticModel"
312            << G4endl;                             313            << G4endl;
313   }                                               314   }
314 #endif                                            315 #endif
315                                                   316   
316   G4double electronEnergy0 = aDynamicElectron-    317   G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
317   G4double cosTheta = 0.;                         318   G4double cosTheta = 0.;
318                                                   319 
319   if (electronEnergy0<intermediateEnergyLimit)    320   if (electronEnergy0<intermediateEnergyLimit)
320   {                                               321   {
321 #ifdef SR_VERBOSE                                 322 #ifdef SR_VERBOSE
322     if (verboseLevel > 3)                         323     if (verboseLevel > 3)
323       {G4cout << "---> Using Brenner & Zaider     324       {G4cout << "---> Using Brenner & Zaider model" << G4endl;}
324 #endif                                            325 #endif
325     cosTheta = BrennerZaiderRandomizeCosTheta(    326     cosTheta = BrennerZaiderRandomizeCosTheta(electronEnergy0);
326   }                                               327   }
327                                                   328 
328   if (electronEnergy0>=intermediateEnergyLimit    329   if (electronEnergy0>=intermediateEnergyLimit)
329   {                                               330   {
330 #ifdef SR_VERBOSE                                 331 #ifdef SR_VERBOSE
331     if (verboseLevel > 3)                         332     if (verboseLevel > 3)
332       {G4cout << "---> Using Screened Rutherfo    333       {G4cout << "---> Using Screened Rutherford model" << G4endl;}
333 #endif                                            334 #endif
334     G4double z = 10.;                             335     G4double z = 10.;
335     cosTheta = ScreenedRutherfordRandomizeCosT    336     cosTheta = ScreenedRutherfordRandomizeCosTheta(electronEnergy0,z);
336   }                                               337   }
337                                                   338 
338   G4double phi = 2. * pi * G4UniformRand();       339   G4double phi = 2. * pi * G4UniformRand();
339                                                   340 
340   G4ThreeVector zVers = aDynamicElectron->GetM    341   G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
341   G4ThreeVector xVers = zVers.orthogonal();       342   G4ThreeVector xVers = zVers.orthogonal();
342   G4ThreeVector yVers = zVers.cross(xVers);       343   G4ThreeVector yVers = zVers.cross(xVers);
343                                                   344 
344   G4double xDir = std::sqrt(1. - cosTheta*cosT    345   G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
345   G4double yDir = xDir;                           346   G4double yDir = xDir;
346   xDir *= std::cos(phi);                          347   xDir *= std::cos(phi);
347   yDir *= std::sin(phi);                          348   yDir *= std::sin(phi);
348                                                   349 
349   G4ThreeVector zPrimeVers((xDir*xVers + yDir*    350   G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
350                                                   351 
351   fParticleChangeForGamma->ProposeMomentumDire    352   fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit());
352                                                   353 
353   fParticleChangeForGamma->SetProposedKineticE    354   fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
354 }                                                 355 }
355                                                   356 
356 //....oooOO0OOooo........oooOO0OOooo........oo    357 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
357                                                   358 
358 G4double G4DNAScreenedRutherfordElasticModel::    359 G4double G4DNAScreenedRutherfordElasticModel::
359 BrennerZaiderRandomizeCosTheta(G4double k)        360 BrennerZaiderRandomizeCosTheta(G4double k)
360 {                                                 361 {
361   //  d sigma_el                         1        362   //  d sigma_el                         1                                 beta(K)
362   // ------------ (K) ~ ----------------------    363   // ------------ (K) ~ --------------------------------- + ---------------------------------
363   //   d Omega           (1 + 2 gamma(K) - cos    364   //   d Omega           (1 + 2 gamma(K) - cos(theta))^2     (1 + 2 delta(K) + cos(theta))^2
364   //                                              365   //
365   // Maximum is < 1/(4 gamma(K)^2) + beta(K)/(    366   // Maximum is < 1/(4 gamma(K)^2) + beta(K)/((2+2delta(K))^2)
366   //                                              367   //
367   // Phys. Med. Biol. 29 N.4 (1983) 443-447       368   // Phys. Med. Biol. 29 N.4 (1983) 443-447
368                                                   369 
369   // gamma(K), beta(K) and delta(K) are polyno    370   // gamma(K), beta(K) and delta(K) are polynomials with coefficients for energy measured in eV
370                                                   371 
371   k /= eV;                                        372   k /= eV;
372                                                   373 
373   G4double beta = G4Exp(CalculatePolynomial(k,    374   G4double beta = G4Exp(CalculatePolynomial(k, betaCoeff));
374   G4double delta = G4Exp(CalculatePolynomial(k    375   G4double delta = G4Exp(CalculatePolynomial(k, deltaCoeff));
375   G4double gamma;                                 376   G4double gamma;
376                                                   377 
377   if (k > 100.)                                   378   if (k > 100.)
378   {                                               379   {
379     gamma = CalculatePolynomial(k, gamma100_20    380     gamma = CalculatePolynomial(k, gamma100_200Coeff);
380     // Only in this case it is not the exponen    381     // Only in this case it is not the exponent of the polynomial
381   }                                               382   }
382   else                                            383   else
383   {                                               384   {
384     if (k > 10)                                   385     if (k > 10)
385     {                                             386     {
386       gamma = G4Exp(CalculatePolynomial(k, gam    387       gamma = G4Exp(CalculatePolynomial(k, gamma10_100Coeff));
387     }                                             388     }
388     else                                          389     else
389     {                                             390     {
390       gamma = G4Exp(CalculatePolynomial(k, gam    391       gamma = G4Exp(CalculatePolynomial(k, gamma035_10Coeff));
391     }                                             392     }
392   }                                               393   }
393                                                   394 
394   // ***** Original method                        395   // ***** Original method
395                                                   396   
396   if (!fasterCode)                                397   if (!fasterCode)
397   {                                               398   {
398                                                   399 
399    G4double oneOverMax = 1.                       400    G4double oneOverMax = 1.
400       / (1. / (4. * gamma * gamma) + beta         401       / (1. / (4. * gamma * gamma) + beta
401           / ((2. + 2. * delta) * (2. + 2. * de    402           / ((2. + 2. * delta) * (2. + 2. * delta)));
402                                                   403 
403    G4double cosTheta = 0.;                        404    G4double cosTheta = 0.;
404    G4double leftDenominator = 0.;                 405    G4double leftDenominator = 0.;
405    G4double rightDenominator = 0.;                406    G4double rightDenominator = 0.;
406    G4double fCosTheta = 0.;                       407    G4double fCosTheta = 0.;
407                                                   408 
408    do                                             409    do
409    {                                              410    {
410      cosTheta = 2. * G4UniformRand()- 1.;         411      cosTheta = 2. * G4UniformRand()- 1.;
411                                                   412 
412      leftDenominator = (1. + 2.*gamma - cosThe    413      leftDenominator = (1. + 2.*gamma - cosTheta);
413      rightDenominator = (1. + 2.*delta + cosTh    414      rightDenominator = (1. + 2.*delta + cosTheta);
414      if ( (leftDenominator * rightDenominator)    415      if ( (leftDenominator * rightDenominator) != 0. )
415      {                                            416      {
416        fCosTheta = oneOverMax * (1./(leftDenom    417        fCosTheta = oneOverMax * (1./(leftDenominator*leftDenominator)
417                                  + beta/(right    418                                  + beta/(rightDenominator*rightDenominator));
418      }                                            419      }
419    }                                              420    }
420    while (fCosTheta < G4UniformRand());           421    while (fCosTheta < G4UniformRand());
421                                                   422 
422    return cosTheta;                               423    return cosTheta;
423   }                                               424   }
424                                                   425 
425   // ***** Alternative method using cumulative    426   // ***** Alternative method using cumulative probability
426                                                   427   
427   if (fasterCode)                                 428   if (fasterCode)
428   {                                               429   {
429                                                   430    
430    //                                             431    // 
431    // modified by Shogo OKADA @ KEK, JP, 2016.    432    // modified by Shogo OKADA @ KEK, JP, 2016.2.27(Sat.)
432    //                                             433    // 
433    // An integral of differential cross-sectio    434    // An integral of differential cross-section formula shown above this member function
434    // (integral variable: cos(theta), integral    435    // (integral variable: cos(theta), integral interval: [-1, x]) is as follows:
435    //                                             436    // 
436    //          1.0 + x                beta * (    437    //          1.0 + x                beta * (1 + x)
437    // I = --------------------- + ------------    438    // I = --------------------- + ----------------------   (1)
438    //      (a - x) * (a + 1.0)      (b + x) *     439    //      (a - x) * (a + 1.0)      (b + x) * (b - 1.0)
439    //                                             440    //
440    // where a = 1.0 + 2.0 * gamma(K), b = 1.0     441    // where a = 1.0 + 2.0 * gamma(K), b = 1.0 + 2.0 * delta(K)
441    //                                             442    // 
442    // Then, a cumulative probability (cp) is a    443    // Then, a cumulative probability (cp) is as follows:
443    //                                             444    // 
444    //  cp       1.0 + x                beta *     445    //  cp       1.0 + x                beta * (1 + x)      
445    // ---- = --------------------- + ---------    446    // ---- = --------------------- + ----------------------  (2)
446    //  S      (a - x) * (a + 1.0)      (b + x)    447    //  S      (a - x) * (a + 1.0)      (b + x) * (b - 1.0) 
447    //                                             448    //
448    // where 1/S is the integral of differnetic    449    // where 1/S is the integral of differnetical cross-section (1) on interval [-1, 1]
449    //                                             450    // 
450    //   1           2.0                     2.    451    //   1           2.0                     2.0 * beta
451    //  --- = ----------------------- + -------    452    //  --- = ----------------------- + -----------------------   (3)
452    //   S     (a - 1.0) * (a + 1.0)     (b + 1    453    //   S     (a - 1.0) * (a + 1.0)     (b + 1.0) * (b - 1.0)
453    //                                             454    //
454    // x is calculated from the quadratic equat    455    // x is calculated from the quadratic equation derived from (2) and (3):
455    //                                             456    //
456    // A * x^2 + B * x + C = 0                     457    // A * x^2 + B * x + C = 0
457    //                                             458    //
458    // where A, B, anc C are coefficients of th    459    // where A, B, anc C are coefficients of the equation:
459    //  A = S * {(b - 1.0) - beta * (a + 1.0)}     460    //  A = S * {(b - 1.0) - beta * (a + 1.0)} + cp * (a + 1.0) * (b - 1.0),
460    //  B = S * {(b - 1.0) * (b + 1.0) + beta *    461    //  B = S * {(b - 1.0) * (b + 1.0) + beta * (a - 1.0) * (a + 1.0)} - cp * (a + 1.0) * (b - 1.0) * (a - b)
461    //  C = S * {b * (b - 1.0) + beta * a * (a     462    //  C = S * {b * (b - 1.0) + beta * a * (a + 1.0)} - cp * (a + 1.0) * (b - 1.0) * ab
462    //                                             463    //
463                                                   464    
464    // sampling cumulative probability             465    // sampling cumulative probability
465    G4double cp = G4UniformRand();                 466    G4double cp = G4UniformRand();
466                                                   467   
467    G4double a = 1.0 + 2.0 * gamma;                468    G4double a = 1.0 + 2.0 * gamma;
468    G4double b = 1.0 + 2.0 * delta;                469    G4double b = 1.0 + 2.0 * delta;
469    G4double a1 = a - 1.0;                         470    G4double a1 = a - 1.0;
470    G4double a2 = a + 1.0;                         471    G4double a2 = a + 1.0;
471    G4double b1 = b - 1.0;                         472    G4double b1 = b - 1.0;
472    G4double b2 = b + 1.0;                         473    G4double b2 = b + 1.0;
473    G4double c1 = a - b;                           474    G4double c1 = a - b;
474    G4double c2 = a * b;                           475    G4double c2 = a * b;
475                                                   476 
476    G4double S = 2.0 / (a1 * a2) + 2.0 * beta /    477    G4double S = 2.0 / (a1 * a2) + 2.0 * beta / (b1 * b2); S = 1.0 / S;
477                                                   478  
478    // coefficients for the quadratic equation     479    // coefficients for the quadratic equation
479    G4double A = S * (b1 - beta * a2) + cp * a2    480    G4double A = S * (b1 - beta * a2) + cp * a2 * b1;
480    G4double B = S * (b1 * b2 + beta * a1 * a2)    481    G4double B = S * (b1 * b2 + beta * a1 * a2) - cp * a2 * b1 * c1;
481    G4double C = S * (b * b1 + beta * a * a2) -    482    G4double C = S * (b * b1 + beta * a * a2) - cp * a2 * b1 * c2;
482                                                   483   
483    // calculate cos(theta)                        484    // calculate cos(theta)
484    return (-1.0 * B + std::sqrt(B * B - 4.0 *     485    return (-1.0 * B + std::sqrt(B * B - 4.0 * A * C)) / (2.0 * A);
485                                                   486    
486    /*                                             487    /*
487    G4double cosTheta = -1;                        488    G4double cosTheta = -1;
488    G4double cumul = 0;                            489    G4double cumul = 0;
489    G4double value = 0;                            490    G4double value = 0;
490    G4double leftDenominator = 0.;                 491    G4double leftDenominator = 0.;
491    G4double rightDenominator = 0.;                492    G4double rightDenominator = 0.;
492                                                   493 
493    // Number of integration steps in the -1,1     494    // Number of integration steps in the -1,1 range
494    G4int iMax=200;                                495    G4int iMax=200;
495                                                   496 
496    G4double random = G4UniformRand();             497    G4double random = G4UniformRand();
497                                                   498 
498    // Cumulate differential cross section         499    // Cumulate differential cross section
499    for (G4int i=0; i<iMax; i++)                   500    for (G4int i=0; i<iMax; i++)
500    {                                              501    {
501     cosTheta = -1 + i*2./(iMax-1);                502     cosTheta = -1 + i*2./(iMax-1);
502     leftDenominator = (1. + 2.*gamma - cosThet    503     leftDenominator = (1. + 2.*gamma - cosTheta);
503     rightDenominator = (1. + 2.*delta + cosThe    504     rightDenominator = (1. + 2.*delta + cosTheta);
504     if ( (leftDenominator * rightDenominator)     505     if ( (leftDenominator * rightDenominator) != 0. )
505     {                                             506     {
506       cumul = cumul + (1./(leftDenominator*lef    507       cumul = cumul + (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator));
507     }                                             508     }
508    }                                              509    }
509                                                   510 
510    // Select cosTheta                             511    // Select cosTheta
511    for (G4int i=0; i<iMax; i++)                   512    for (G4int i=0; i<iMax; i++)
512    {                                              513    {
513      cosTheta = -1 + i*2./(iMax-1);               514      cosTheta = -1 + i*2./(iMax-1);
514      leftDenominator = (1. + 2.*gamma - cosThe    515      leftDenominator = (1. + 2.*gamma - cosTheta);
515      rightDenominator = (1. + 2.*delta + cosTh    516      rightDenominator = (1. + 2.*delta + cosTheta);
516      if (cumul !=0 && (leftDenominator * right    517      if (cumul !=0 && (leftDenominator * rightDenominator) != 0.)
517       value = value + (1./(leftDenominator*lef    518       value = value + (1./(leftDenominator*leftDenominator) + beta/(rightDenominator*rightDenominator)) / cumul;
518      if (random < value) break;                   519      if (random < value) break;
519    }                                              520    }
520                                                   521 
521    return cosTheta;                               522    return cosTheta;
522    */                                             523    */
523   }                                               524   }
524                                                   525  
525   return 0.;                                      526   return 0.;
526 }                                                 527 }
527                                                   528 
528 //....oooOO0OOooo........oooOO0OOooo........oo    529 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
529                                                   530 
530 G4double G4DNAScreenedRutherfordElasticModel::    531 G4double G4DNAScreenedRutherfordElasticModel::
531 CalculatePolynomial(G4double k,                   532 CalculatePolynomial(G4double k,
532                     std::vector<G4double>& vec    533                     std::vector<G4double>& vec)
533 {                                                 534 {
534   // Sum_{i=0}^{size-1} vector_i k^i              535   // Sum_{i=0}^{size-1} vector_i k^i
535   //                                              536   //
536   // Phys. Med. Biol. 29 N.4 (1983) 443-447       537   // Phys. Med. Biol. 29 N.4 (1983) 443-447
537                                                   538 
538   G4double result = 0.;                           539   G4double result = 0.;
539   size_t size = vec.size();                       540   size_t size = vec.size();
540                                                   541 
541   while (size > 0)                                542   while (size > 0)
542   {                                               543   {
543     size--;                                       544     size--;
544                                                   545 
545     result *= k;                                  546     result *= k;
546     result += vec[size];                          547     result += vec[size];
547   }                                               548   }
548                                                   549 
549   return result;                                  550   return result;
550 }                                                 551 }
551                                                   552 
552 //....oooOO0OOooo........oooOO0OOooo........oo    553 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
553                                                   554 
554 G4double G4DNAScreenedRutherfordElasticModel::    555 G4double G4DNAScreenedRutherfordElasticModel::
555 ScreenedRutherfordRandomizeCosTheta(G4double k    556 ScreenedRutherfordRandomizeCosTheta(G4double k,
556                                     G4double z    557                                     G4double z)
557 {                                                 558 {
558                                                   559 
559   //  d sigma_el                sigma_Ruth(K)     560   //  d sigma_el                sigma_Ruth(K)
560   // ------------ (K) ~ ----------------------    561   // ------------ (K) ~ -----------------------------
561   //   d Omega           (1 + 2 n(K) - cos(the    562   //   d Omega           (1 + 2 n(K) - cos(theta))^2
562   //                                              563   //
563   // We extract cos(theta) distributed as (1 +    564   // We extract cos(theta) distributed as (1 + 2 n(K) - cos(theta))^-2
564   //                                              565   //
565   // Maximum is for theta=0: 1/(4 n(K)^2) (Whe    566   // Maximum is for theta=0: 1/(4 n(K)^2) (When n(K) is positive, that is always satisfied within the validity of the process)
566   //                                              567   //
567   // Phys. Med. Biol. 45 (2000) 3171-3194         568   // Phys. Med. Biol. 45 (2000) 3171-3194
568                                                   569 
569   // ***** Original method                        570   // ***** Original method
570                                                   571 
571   if (!fasterCode)                                572   if (!fasterCode)
572   {                                               573   {
573                                                   574 
574    G4double n = ScreeningFactor(k, z);            575    G4double n = ScreeningFactor(k, z);
575                                                   576 
576    G4double oneOverMax = (4. * n * n);            577    G4double oneOverMax = (4. * n * n);
577                                                   578 
578    G4double cosTheta = 0.;                        579    G4double cosTheta = 0.;
579    G4double fCosTheta;                            580    G4double fCosTheta;
580                                                   581 
581    do                                             582    do
582    {                                              583    {
583      cosTheta = 2. * G4UniformRand()- 1.;         584      cosTheta = 2. * G4UniformRand()- 1.;
584      fCosTheta = (1 + 2.*n - cosTheta);           585      fCosTheta = (1 + 2.*n - cosTheta);
585      if (fCosTheta !=0.) fCosTheta = oneOverMa    586      if (fCosTheta !=0.) fCosTheta = oneOverMax / (fCosTheta*fCosTheta);
586    }                                              587    }
587    while (fCosTheta < G4UniformRand());           588    while (fCosTheta < G4UniformRand());
588                                                   589 
589    return cosTheta;                               590    return cosTheta;
590   }                                               591   }
591                                                   592 
592   // ***** Alternative method using cumulative    593   // ***** Alternative method using cumulative probability
                                                   >> 594   if (fasterCode)
                                                   >> 595   {
593                                                   596   
594                                                << 
595    //                                             597    //
596    // modified by Shogo OKADA @ KEK, JP, 2016.    598    // modified by Shogo OKADA @ KEK, JP, 2016.2.27(Sat.)
597    //                                             599    //
598    // The cumulative probability (cp) is calcu    600    // The cumulative probability (cp) is calculated by integrating 
599    // the differential cross-section fomula wi    601    // the differential cross-section fomula with cos(theta):
600    //                                             602    // 
601    //         n(K) * (1.0 + cos(theta))           603    //         n(K) * (1.0 + cos(theta))
602    //  cp = ---------------------------------     604    //  cp = ---------------------------------
603    //         1.0 + 2.0 * n(K) - cos(theta)       605    //         1.0 + 2.0 * n(K) - cos(theta)
604    //                                             606    //
605    // Then, cos(theta) is as follows:             607    // Then, cos(theta) is as follows:
606    //                                             608    //
607    //               cp * (1.0 + 2.0 * n(K)) -     609    //               cp * (1.0 + 2.0 * n(K)) - n(K)
608    // cos(theta) = ---------------------------    610    // cos(theta) = --------------------------------
609    //                       n(k) + cp             611    //                       n(k) + cp
610    //                                             612    // 
611    // where, K is kinetic energy, n(K) is scre    613    // where, K is kinetic energy, n(K) is screeing factor, and cp is cumulative probability 
612    //                                             614    // 
613                                                   615   
614    G4double n = ScreeningFactor(k, z);            616    G4double n = ScreeningFactor(k, z);
615    G4double cp = G4UniformRand();                 617    G4double cp = G4UniformRand();
616    G4double numerator = cp * (1.0 + 2.0 * n) -    618    G4double numerator = cp * (1.0 + 2.0 * n) - n;
617    G4double denominator = n + cp;                 619    G4double denominator = n + cp;
618    return numerator / denominator;                620    return numerator / denominator;
619                                                   621       
620    /*                                             622    /*
621    G4double cosTheta = -1;                        623    G4double cosTheta = -1;
622    G4double cumul = 0;                            624    G4double cumul = 0;
623    G4double value = 0;                            625    G4double value = 0;
624    G4double n = ScreeningFactor(k, z);            626    G4double n = ScreeningFactor(k, z);
625    G4double fCosTheta;                            627    G4double fCosTheta;
626                                                   628 
627    // Number of integration steps in the -1,1     629    // Number of integration steps in the -1,1 range
628    G4int iMax=200;                                630    G4int iMax=200;
629                                                   631 
630    G4double random = G4UniformRand();             632    G4double random = G4UniformRand();
631                                                   633 
632    // Cumulate differential cross section         634    // Cumulate differential cross section
633    for (G4int i=0; i<iMax; i++)                   635    for (G4int i=0; i<iMax; i++)
634    {                                              636    {
635      cosTheta = -1 + i*2./(iMax-1);               637      cosTheta = -1 + i*2./(iMax-1);
636      fCosTheta = (1 + 2.*n - cosTheta);           638      fCosTheta = (1 + 2.*n - cosTheta);
637      if (fCosTheta !=0.) cumul = cumul + 1./(f    639      if (fCosTheta !=0.) cumul = cumul + 1./(fCosTheta*fCosTheta);
638    }                                              640    }
639                                                   641 
640    // Select cosTheta                             642    // Select cosTheta
641    for (G4int i=0; i<iMax; i++)                   643    for (G4int i=0; i<iMax; i++)
642    {                                              644    {
643      cosTheta = -1 + i*2./(iMax-1);               645      cosTheta = -1 + i*2./(iMax-1);
644      fCosTheta = (1 + 2.*n - cosTheta);           646      fCosTheta = (1 + 2.*n - cosTheta);
645      if (cumul !=0.) value = value + (1./(fCos    647      if (cumul !=0.) value = value + (1./(fCosTheta*fCosTheta)) / cumul;
646      if (random < value) break;                   648      if (random < value) break;
647    }                                              649    }
648    return cosTheta;                               650    return cosTheta;
649    */                                             651    */
650                                                << 652   }
651                                                   653 
652   //return 0.;                                 << 654   return 0.;
653 }                                                 655 }
654                                                   656 
655                                                   657 
656                                                   658