Geant4 Cross Reference |
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 // $Id: G4PAIySection.cc,v 1.3 2007/10/01 18:38:10 vnivanch Exp $ >> 28 // GEANT4 tag $Name: geant4-09-01-patch-02 $ >> 29 // 27 // 30 // 28 // G4PAIySection.cc -- class implementation fi 31 // G4PAIySection.cc -- class implementation file 29 // 32 // 30 // GEANT 4 class implementation file 33 // GEANT 4 class implementation file 31 // 34 // 32 // For information related to this code, pleas 35 // For information related to this code, please, contact 33 // the Geant4 Collaboration. 36 // the Geant4 Collaboration. 34 // 37 // 35 // R&D: Vladimir.Grichine@cern.ch 38 // R&D: Vladimir.Grichine@cern.ch 36 // 39 // 37 // History: 40 // History: 38 // 41 // 39 // 01.10.07, V.Ivanchenko create using V.Grich 42 // 01.10.07, V.Ivanchenko create using V.Grichine G4PAIxSection class 40 // 26.07.09, V.Ivanchenko added protection for << 41 // low-density materials << 42 // 21.11.10 V. Grichine bug fixed in Initialis << 43 // material. Warning: the table is << 44 // 23.06.13 V.Grichine arrays->G4DataVectors << 45 // 43 // 46 44 47 #include "G4PAIySection.hh" 45 #include "G4PAIySection.hh" 48 46 49 #include "globals.hh" 47 #include "globals.hh" 50 #include "G4PhysicalConstants.hh" << 51 #include "G4SystemOfUnits.hh" << 52 #include "G4ios.hh" 48 #include "G4ios.hh" 53 #include "G4Poisson.hh" 49 #include "G4Poisson.hh" 54 #include "G4Material.hh" 50 #include "G4Material.hh" 55 #include "G4MaterialCutsCouple.hh" 51 #include "G4MaterialCutsCouple.hh" 56 #include "G4SandiaTable.hh" 52 #include "G4SandiaTable.hh" 57 #include "G4Exp.hh" << 58 #include "G4Log.hh" << 59 53 60 using namespace std; 54 using namespace std; 61 55 62 // Local class constants 56 // Local class constants 63 57 64 const G4double G4PAIySection::fDelta = 0.005; << 58 const G4double G4PAIySection::fDelta = 0.005 ; // energy shift from interval border 65 const G4double G4PAIySection::fError = 0.005; << 59 const G4double G4PAIySection::fError = 0.005 ; // error in lin-log approximation 66 60 67 const G4int G4PAIySection::fMaxSplineSize = 50 << 61 const G4int G4PAIySection::fMaxSplineSize = 500 ; // Max size of output spline 68 << 62 // arrays 69 63 70 ////////////////////////////////////////////// 64 ////////////////////////////////////////////////////////////////// 71 // 65 // 72 // Constructor 66 // Constructor 73 // 67 // 74 68 75 G4PAIySection::G4PAIySection() 69 G4PAIySection::G4PAIySection() 76 { << 70 {} 77 fSandia = nullptr; << 78 fDensity = fElectronDensity = fNormalization << 79 fIntervalNumber = fSplineNumber = 0; << 80 fVerbose = 0; << 81 << 82 betaBohr = fine_structure_const; << 83 G4double cofBetaBohr = 4.0; << 84 G4double betaBohr2 = fine_structure_const*fi << 85 betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr; << 86 << 87 fSplineEnergy = G4DataVector(fMaxSp << 88 fRePartDielectricConst = G4DataVector(fMaxSp << 89 fImPartDielectricConst = G4DataVector(fMaxSp << 90 fIntegralTerm = G4DataVector(fMaxSp << 91 fDifPAIySection = G4DataVector(fMaxSp << 92 fdNdxCerenkov = G4DataVector(fMaxSp << 93 fdNdxPlasmon = G4DataVector(fMaxSp << 94 fIntegralPAIySection = G4DataVector(fMaxSp << 95 fIntegralPAIdEdx = G4DataVector(fMaxSp << 96 fIntegralCerenkov = G4DataVector(fMaxSp << 97 fIntegralPlasmon = G4DataVector(fMaxSp << 98 << 99 for( G4int i = 0; i < 500; ++i ) << 100 { << 101 for( G4int j = 0; j < 112; ++j ) { fPAItab << 102 } << 103 } << 104 71 105 ////////////////////////////////////////////// 72 //////////////////////////////////////////////////////////////////////////// 106 // 73 // 107 // << 74 // Destructor 108 75 109 G4double G4PAIySection::GetLorentzFactor(G4int << 76 G4PAIySection::~G4PAIySection() 110 { << 77 {} 111 return fLorentzFactor[j]; << 112 } << 113 78 114 ////////////////////////////////////////////// 79 //////////////////////////////////////////////////////////////////////// 115 // 80 // 116 // Constructor with beta*gamma square value ca << 81 // Test Constructor with beta*gamma square value 117 82 118 void G4PAIySection::Initialize( const G4Materi 83 void G4PAIySection::Initialize( const G4Material* material, 119 G4double maxEn << 84 G4double maxEnergyTransfer, 120 G4double betaG << 85 G4double betaGammaSq) 121 G4SandiaTable* << 122 { 86 { 123 if(fVerbose > 0) << 87 G4int i, j, numberOfElements ; 124 { << 88 125 G4cout<<G4endl; << 89 fDensity = material->GetDensity(); 126 G4cout<<"G4PAIySection::Initialize(...,G4S << 90 fElectronDensity = material->GetElectronDensity() ; 127 G4cout<<G4endl; << 91 numberOfElements = material->GetNumberOfElements() ; 128 } << 129 G4int i, j; << 130 92 131 fSandia = sandia; << 93 fSandia = material->GetSandiaTable(); 132 fIntervalNumber = sandia->GetMaxInterval(); << 133 fDensity = material->GetDensity(); << 134 fElectronDensity = material->GetElectronDens << 135 94 136 // fIntervalNumber--; << 95 fIntervalNumber = fSandia->GetMaxInterval(); 137 96 138 if( fVerbose > 0 ) << 97 fIntervalNumber--; 139 { << 140 G4cout<<"fDensity = "<<fDensity<<"\t"<<fEl << 141 <<fIntervalNumber<< " (beta*gamma)^2 << 142 } << 143 fEnergyInterval = G4DataVector(fIntervalNumb << 144 fA1 = G4DataVector(fIntervalNumb << 145 fA2 = G4DataVector(fIntervalNumb << 146 fA3 = G4DataVector(fIntervalNumb << 147 fA4 = G4DataVector(fIntervalNumb << 148 98 149 for( i = 1; i <= fIntervalNumber; ++i ) << 99 for(i=1;i<=fIntervalNumber;i++) 150 { << 100 { 151 if ( sandia->GetSandiaMatTablePAI(i-1,0) < << 101 G4double e = fSandia->GetSandiaMatTablePAI(i,0); 152 { << 102 if(e >= maxEnergyTransfer || i > fIntervalNumber) 153 fIntervalNumber--; << 103 { 154 continue; << 104 fEnergyInterval[i] = maxEnergyTransfer ; 155 } << 105 fIntervalNumber = i ; 156 if( ( sandia->GetSandiaMatTablePAI(i-1,0) << 106 break; 157 || i >= fIntervalNumber ) << 107 } 158 { << 108 fEnergyInterval[i] = e; 159 fEnergyInterval[i] = maxEnergyTransfer; << 109 fA1[i] = fSandia->GetSandiaMatTablePAI(i,1); 160 fIntervalNumber = i; << 110 fA2[i] = fSandia->GetSandiaMatTablePAI(i,2); 161 break; << 111 fA3[i] = fSandia->GetSandiaMatTablePAI(i,3); 162 } << 112 fA4[i] = fSandia->GetSandiaMatTablePAI(i,4); 163 fEnergyInterval[i] = sandia->GetSandiaMatT << 113 164 fA1[i] = sandia->GetSandiaMatT << 114 } 165 fA2[i] = sandia->GetSandiaMatT << 115 if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer) 166 fA3[i] = sandia->GetSandiaMatT << 116 { 167 fA4[i] = sandia->GetSandiaMatT << 117 fIntervalNumber++; 168 << 118 fEnergyInterval[fIntervalNumber] = maxEnergyTransfer ; 169 if( fVerbose > 0 ) { << 119 fA1[fIntervalNumber] = fA1[fIntervalNumber-1] ; 170 G4cout<<i<<"\t"<<fEnergyInterval[i]/keV< << 120 fA2[fIntervalNumber] = fA2[fIntervalNumber-1] ; 171 <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4en << 121 fA3[fIntervalNumber] = fA3[fIntervalNumber-1] ; 172 } << 122 fA4[fIntervalNumber] = fA4[fIntervalNumber-1] ; 173 } << 123 } 174 if( fVerbose > 0 ) { << 124 175 G4cout<<"last i = "<<i<<"; "<<"fIntervalNu << 125 // Now checking, if two borders are too close together 176 <<fIntervalNumber<<G4endl; << 126 for(i=1;i<fIntervalNumber;i++) 177 } << 178 if( fEnergyInterval[fIntervalNumber] != maxE << 179 { << 180 fIntervalNumber++; << 181 fEnergyInterval[fIntervalNumber] = maxEn << 182 } << 183 if( fVerbose > 0 ) << 184 { << 185 for( i = 1; i <= fIntervalNumber; ++i ) << 186 { << 187 G4cout<<i<<"\t"<<fEnergyInterval[i]/keV< << 188 <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl; << 189 } << 190 } << 191 if( fVerbose > 0 ) { << 192 G4cout<<"Now checking, if two borders are << 193 } << 194 for( i = 1; i < fIntervalNumber; ++i ) << 195 { << 196 if( fEnergyInterval[i+1]-fEnergyInterval[i << 197 1.5*fDelta*(fEnergyInterval[i+1]+fEne << 198 else << 199 { << 200 for( j = i; j < fIntervalNumber; j++ ) << 201 { 127 { 202 fEnergyInterval[j] = fEnergyInte << 128 // G4cout<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t" 203 fA1[j] = fA1[j+1]; << 129 // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl ; 204 fA2[j] = fA2[j+1]; << 130 if(fEnergyInterval[i+1]-fEnergyInterval[i] > 205 fA3[j] = fA3[j+1]; << 131 1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i])) 206 fA4[j] = fA4[j+1]; << 132 { >> 133 continue ; >> 134 } >> 135 else >> 136 { >> 137 for(j=i;j<fIntervalNumber;j++) >> 138 { >> 139 fEnergyInterval[j] = fEnergyInterval[j+1] ; >> 140 fA1[j] = fA1[j+1] ; >> 141 fA2[j] = fA2[j+1] ; >> 142 fA3[j] = fA3[j+1] ; >> 143 fA4[j] = fA4[j+1] ; >> 144 } >> 145 fIntervalNumber-- ; >> 146 i-- ; >> 147 } 207 } 148 } 208 fIntervalNumber--; << 209 } << 210 } << 211 if( fVerbose > 0 ) << 212 { << 213 for( i = 1; i <= fIntervalNumber; ++i ) << 214 { << 215 G4cout<<i<<"\t"<<fEnergyInterval[i]/keV< << 216 <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl; << 217 } << 218 } << 219 // Preparation of fSplineEnergy array corres << 220 << 221 ComputeLowEnergyCof(material); << 222 << 223 G4double betaGammaSqRef = << 224 fLorentzFactor[fRefGammaNumber]*fLorentzFa << 225 149 226 NormShift(betaGammaSqRef); << 150 // Preparation of fSplineEnergy array corresponding to min ionisation, G~4 227 SplainPAI(betaGammaSqRef); << 228 151 229 // Preparation of integral PAI cross section << 152 G4double betaGammaSqRef = 230 << 153 fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1; 231 for( i = 1; i <= fSplineNumber; ++i ) << 232 { << 233 fDifPAIySection[i] = DifPAIySection(i,bet << 234 << 235 if( fVerbose > 0 ) G4cout<<i<<"; dNdxPAI << 236 } << 237 IntegralPAIySection(); << 238 } << 239 << 240 ////////////////////////////////////////////// << 241 // << 242 // Compute low energy cof. It reduces PAI xsc << 243 // << 244 << 245 void G4PAIySection::ComputeLowEnergyCof(const << 246 { << 247 G4int i, numberOfElements = (G4int)material- << 248 G4double sumZ = 0., sumCof = 0.; << 249 154 250 static const G4double p0 = 1.20923e+00; << 155 NormShift(betaGammaSqRef) ; 251 static const G4double p1 = 3.53256e-01; << 156 SplainPAI(betaGammaSqRef) ; 252 static const G4double p2 = -1.45052e-03; << 157 253 << 158 // Preparation of integral PAI cross section for input betaGammaSq 254 G4double* thisMaterialZ = new G4double[num << 255 G4double* thisMaterialCof = new G4double[num << 256 159 257 for( i = 0; i < numberOfElements; ++i ) << 160 for(i = 1 ; i <= fSplineNumber ; i++) 258 { << 161 { 259 thisMaterialZ[i] = material->GetElement(i) << 162 fDifPAIySection[i] = DifPAIySection(i,betaGammaSq); 260 sumZ += thisMaterialZ[i]; << 163 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq); 261 thisMaterialCof[i] = p0+p1*thisMaterialZ[i << 164 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq); 262 } << 165 } 263 for( i = 0; i < numberOfElements; ++i ) << 166 IntegralPAIySection() ; 264 { << 167 IntegralCerenkov() ; 265 sumCof += thisMaterialCof[i]*thisMaterialZ << 168 IntegralPlasmon() ; 266 } << 267 fLowEnergyCof = sumCof; << 268 delete [] thisMaterialZ; << 269 delete [] thisMaterialCof; << 270 // G4cout<<"fLowEnergyCof = "<<fLowEnergyCof << 271 } 169 } 272 170 273 ////////////////////////////////////////////// 171 ///////////////////////////////////////////////////////////////////////// 274 // 172 // 275 // General control function for class G4PAIySe 173 // General control function for class G4PAIySection 276 // 174 // 277 175 278 void G4PAIySection::InitPAI() 176 void G4PAIySection::InitPAI() 279 { 177 { 280 G4int i; << 178 G4int i ; 281 G4double betaGammaSq = fLorentzFactor[fRefG 179 G4double betaGammaSq = fLorentzFactor[fRefGammaNumber]* 282 fLorentzFactor[fRefG 180 fLorentzFactor[fRefGammaNumber] - 1; 283 181 284 // Preparation of integral PAI cross sectio 182 // Preparation of integral PAI cross section for reference gamma 285 183 286 NormShift(betaGammaSq); << 184 NormShift(betaGammaSq) ; 287 SplainPAI(betaGammaSq); << 185 SplainPAI(betaGammaSq) ; 288 186 289 IntegralPAIySection(); << 187 IntegralPAIySection() ; 290 IntegralCerenkov(); << 188 IntegralCerenkov() ; 291 IntegralPlasmon(); << 189 IntegralPlasmon() ; 292 190 293 for( i = 0; i<= fSplineNumber; ++i) << 191 for(i = 0 ; i<=fSplineNumber ; i++) 294 { 192 { 295 fPAItable[i][fRefGammaNumber] = fIntegral << 193 fPAItable[i][fRefGammaNumber] = fIntegralPAIySection[i] ; 296 << 194 if(i != 0) 297 if(i != 0) fPAItable[i][0] = fSplineEner << 195 { >> 196 fPAItable[i][0] = fSplineEnergy[i] ; >> 197 } 298 } 198 } 299 fPAItable[0][0] = fSplineNumber; << 199 fPAItable[0][0] = fSplineNumber ; 300 200 301 for( G4int j = 1; j < 112; ++j) // fo << 201 for(G4int j = 1 ; j < 112 ; j++) // for other gammas 302 { 202 { 303 if( j == fRefGammaNumber ) continue; << 203 if( j == fRefGammaNumber ) continue ; 304 204 305 betaGammaSq = fLorentzFactor[j]*fLorentz << 205 betaGammaSq = fLorentzFactor[j]*fLorentzFactor[j] - 1 ; 306 206 307 for(i = 1; i <= fSplineNumber; ++i) << 207 for(i = 1 ; i <= fSplineNumber ; i++) 308 { 208 { 309 fDifPAIySection[i] = DifPAIySection(i 209 fDifPAIySection[i] = DifPAIySection(i,betaGammaSq); 310 fdNdxCerenkov[i] = PAIdNdxCerenkov( 210 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq); 311 fdNdxPlasmon[i] = PAIdNdxPlasmon(i 211 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq); 312 } 212 } 313 IntegralPAIySection(); << 213 IntegralPAIySection() ; 314 IntegralCerenkov(); << 214 IntegralCerenkov() ; 315 IntegralPlasmon(); << 215 IntegralPlasmon() ; 316 216 317 for(i = 0; i <= fSplineNumber; ++i) << 217 for(i = 0 ; i <= fSplineNumber ; i++) 318 { 218 { 319 fPAItable[i][j] = fIntegralPAIySection << 219 fPAItable[i][j] = fIntegralPAIySection[i] ; 320 } 220 } 321 } 221 } >> 222 322 } 223 } 323 224 324 ////////////////////////////////////////////// 225 /////////////////////////////////////////////////////////////////////// 325 // 226 // 326 // Shifting from borders to intervals Creation 227 // Shifting from borders to intervals Creation of first energy points 327 // 228 // 328 229 329 void G4PAIySection::NormShift(G4double betaGam 230 void G4PAIySection::NormShift(G4double betaGammaSq) 330 { 231 { 331 G4int i, j; << 232 G4int i, j ; 332 233 333 for( i = 1; i <= fIntervalNumber-1; ++i) << 234 for( i = 1 ; i <= fIntervalNumber-1 ; i++ ) 334 { 235 { 335 for( j = 1; j <= 2; ++j) << 236 for( j = 1 ; j <= 2 ; j++ ) 336 { 237 { 337 fSplineNumber = (i-1)*2 + j; << 238 fSplineNumber = (i-1)*2 + j ; 338 239 339 if( j == 1 ) fSplineEnergy[fSplineNumber 240 if( j == 1 ) fSplineEnergy[fSplineNumber] = fEnergyInterval[i ]*(1+fDelta); 340 else fSplineEnergy[fSplineNumber 241 else fSplineEnergy[fSplineNumber] = fEnergyInterval[i+1]*(1-fDelta); 341 // G4cout<<"cn = "<<fSplineNumber<<"; 242 // G4cout<<"cn = "<<fSplineNumber<<"; "<<"energy = " 342 // <<fSplineEnergy[fSplineNumber]<<G4en 243 // <<fSplineEnergy[fSplineNumber]<<G4endl; 343 } 244 } 344 } 245 } 345 fIntegralTerm[1]=RutherfordIntegral(1,fEnerg 246 fIntegralTerm[1]=RutherfordIntegral(1,fEnergyInterval[1],fSplineEnergy[1]); 346 247 347 j = 1; << 248 j = 1 ; 348 249 349 for(i=2;i<=fSplineNumber;++i) << 250 for(i=2;i<=fSplineNumber;i++) 350 { 251 { 351 if(fSplineEnergy[i]<fEnergyInterval[j+1]) 252 if(fSplineEnergy[i]<fEnergyInterval[j+1]) 352 { 253 { 353 fIntegralTerm[i] = fIntegralTerm[i-1] 254 fIntegralTerm[i] = fIntegralTerm[i-1] + 354 RutherfordIntegral << 255 RutherfordIntegral(j,fSplineEnergy[i-1], 355 << 256 fSplineEnergy[i] ) ; 356 } 257 } 357 else 258 else 358 { 259 { 359 G4double x = RutherfordIntegral(j,fSpli 260 G4double x = RutherfordIntegral(j,fSplineEnergy[i-1], 360 fEn << 261 fEnergyInterval[j+1] ) ; 361 j++; 262 j++; 362 fIntegralTerm[i] = fIntegralTerm[i-1] 263 fIntegralTerm[i] = fIntegralTerm[i-1] + x + 363 RutherfordIntegral << 264 RutherfordIntegral(j,fEnergyInterval[j], 364 << 265 fSplineEnergy[i] ) ; 365 } 266 } 366 // G4cout<<i<<"\t"<<fSplineEnergy[i]<<"\t" 267 // G4cout<<i<<"\t"<<fSplineEnergy[i]<<"\t"<<fIntegralTerm[i]<<"\n"<<G4endl; 367 } 268 } 368 static const G4double nfactor = << 269 fNormalizationCof = 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2 ; 369 2*pi*pi*hbarc*hbarc*fine_structure_const/e << 270 fNormalizationCof *= fElectronDensity/fIntegralTerm[fSplineNumber] ; 370 fNormalizationCof = nfactor*fElectronDensity << 371 271 372 // G4cout<<"fNormalizationCof = "<<fNormaliz << 272 // G4cout<<"fNormalizationCof = "<<fNormalizationCof<<G4endl ; 373 273 374 // Calculation of PAI differrential cross-se << 274 // Calculation of PAI differrential cross-section (1/(keV*cm)) 375 // in the energy points near borders of ener << 275 // in the energy points near borders of energy intervals 376 276 377 for(G4int k=1; k<=fIntervalNumber-1; ++k) << 277 for(G4int k=1;k<=fIntervalNumber-1;k++) 378 { 278 { 379 for(j=1; j<=2; ++j) << 279 for(j=1;j<=2;j++) 380 { 280 { 381 i = (k-1)*2 + j; << 281 i = (k-1)*2 + j ; 382 fImPartDielectricConst[i] = fNormaliz 282 fImPartDielectricConst[i] = fNormalizationCof* 383 ImPartDie << 283 ImPartDielectricConst(k,fSplineEnergy[i]); 384 fRePartDielectricConst[i] = fNormaliz 284 fRePartDielectricConst[i] = fNormalizationCof* 385 RePartDie << 285 RePartDielectricConst(fSplineEnergy[i]); 386 fIntegralTerm[i] *= fNormalizationCof 286 fIntegralTerm[i] *= fNormalizationCof; 387 287 388 fDifPAIySection[i] = DifPAIySection(i 288 fDifPAIySection[i] = DifPAIySection(i,betaGammaSq); 389 fdNdxCerenkov[i] = PAIdNdxCerenkov( 289 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq); 390 fdNdxPlasmon[i] = PAIdNdxPlasmon(i 290 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq); 391 } 291 } 392 } 292 } 393 293 394 } // end of NormShift 294 } // end of NormShift 395 295 396 ////////////////////////////////////////////// 296 ///////////////////////////////////////////////////////////////////////// 397 // 297 // 398 // Creation of new energy points as geometrica 298 // Creation of new energy points as geometrical mean of existing 399 // one, calculation PAI_cs for them, while the 299 // one, calculation PAI_cs for them, while the error of logarithmic 400 // linear approximation would be smaller than 300 // linear approximation would be smaller than 'fError' 401 301 402 void G4PAIySection::SplainPAI(G4double betaGam 302 void G4PAIySection::SplainPAI(G4double betaGammaSq) 403 { 303 { 404 G4int k = 1; << 304 G4int k = 1 ; 405 G4int i = 1; << 305 G4int i = 1 ; 406 306 407 while ( (i < fSplineNumber) && (fSplineNumb 307 while ( (i < fSplineNumber) && (fSplineNumber < fMaxSplineSize-1) ) 408 { 308 { 409 if(fSplineEnergy[i+1] > fEnergyInterval[ 309 if(fSplineEnergy[i+1] > fEnergyInterval[k+1]) 410 { 310 { 411 k++; // Here next energy point is << 311 k++ ; // Here next energy point is in next energy interval 412 ++i; << 312 i++; 413 continue; 313 continue; 414 } 314 } 415 // Shifting of arrayes for inserting the << 315 // Shifting of arrayes for inserting the geometrical 416 // average of 'i' and 'i+1' energy point << 316 // average of 'i' and 'i+1' energy points to 'i+1' place 417 fSplineNumber++; 317 fSplineNumber++; 418 318 419 for(G4int j = fSplineNumber; j >= i+2; j << 319 for(G4int j = fSplineNumber; j >= i+2 ; j-- ) 420 { 320 { 421 fSplineEnergy[j] = fSplineEn 321 fSplineEnergy[j] = fSplineEnergy[j-1]; 422 fImPartDielectricConst[j] = fImPartDi 322 fImPartDielectricConst[j] = fImPartDielectricConst[j-1]; 423 fRePartDielectricConst[j] = fRePartDi << 323 fRePartDielectricConst[j] = fRePartDielectricConst[j-1]; 424 fIntegralTerm[j] = fIntegral << 324 fIntegralTerm[j] = fIntegralTerm[j-1]; 425 325 426 fDifPAIySection[j] = fDifPAIySection[ << 326 fDifPAIySection[j] = fDifPAIySection[j-1]; 427 fdNdxCerenkov[j] = fdNdxCerenkov[j- 327 fdNdxCerenkov[j] = fdNdxCerenkov[j-1]; 428 fdNdxPlasmon[j] = fdNdxPlasmon[j-1 328 fdNdxPlasmon[j] = fdNdxPlasmon[j-1]; 429 } 329 } 430 G4double x1 = fSplineEnergy[i]; 330 G4double x1 = fSplineEnergy[i]; 431 G4double x2 = fSplineEnergy[i+1]; 331 G4double x2 = fSplineEnergy[i+1]; 432 G4double yy1 = fDifPAIySection[i]; 332 G4double yy1 = fDifPAIySection[i]; 433 G4double y2 = fDifPAIySection[i+1]; 333 G4double y2 = fDifPAIySection[i+1]; 434 334 435 G4double en1 = sqrt(x1*x2); 335 G4double en1 = sqrt(x1*x2); 436 fSplineEnergy[i+1] = en1; 336 fSplineEnergy[i+1] = en1; 437 337 438 // Calculation of logarithmic linear app << 338 // Calculation of logarithmic linear approximation 439 // in this (enr) energy point, which num << 339 // in this (enr) energy point, which number is 'i+1' now 440 340 441 G4double a = log10(y2/yy1)/log10(x2/x1); 341 G4double a = log10(y2/yy1)/log10(x2/x1); 442 G4double b = log10(yy1) - a*log10(x1); 342 G4double b = log10(yy1) - a*log10(x1); 443 G4double y = a*log10(en1) + b; << 343 G4double y = a*log10(en1) + b ; 444 y = pow(10.,y); 344 y = pow(10.,y); 445 345 446 // Calculation of the PAI dif. cross-sec << 346 // Calculation of the PAI dif. cross-section at this point 447 347 448 fImPartDielectricConst[i+1] = fNormaliza 348 fImPartDielectricConst[i+1] = fNormalizationCof* 449 ImPartDiel << 349 ImPartDielectricConst(k,fSplineEnergy[i+1]); 450 fRePartDielectricConst[i+1] = fNormaliza 350 fRePartDielectricConst[i+1] = fNormalizationCof* 451 RePartDiel << 351 RePartDielectricConst(fSplineEnergy[i+1]); 452 fIntegralTerm[i+1] = fIntegralTerm[i] + 352 fIntegralTerm[i+1] = fIntegralTerm[i] + fNormalizationCof* 453 RutherfordIntegral( << 353 RutherfordIntegral(k,fSplineEnergy[i], 454 354 fSplineEnergy[i+1]); 455 355 456 fDifPAIySection[i+1] = DifPAIySection(i+ 356 fDifPAIySection[i+1] = DifPAIySection(i+1,betaGammaSq); 457 fdNdxCerenkov[i+1] = PAIdNdxCerenkov(i 357 fdNdxCerenkov[i+1] = PAIdNdxCerenkov(i+1,betaGammaSq); 458 fdNdxPlasmon[i+1] = PAIdNdxPlasmon(i+ 358 fdNdxPlasmon[i+1] = PAIdNdxPlasmon(i+1,betaGammaSq); 459 359 460 // Condition for next divisi << 360 // Condition for next division of this segment or to pass 461 // to higher energies << 361 // to higher energies 462 362 463 G4double x = 2*(fDifPAIySection[i+1] - y 363 G4double x = 2*(fDifPAIySection[i+1] - y)/(fDifPAIySection[i+1] + y); 464 364 465 G4double delta = 2.*(fSplineEnergy[i+1]- << 466 /(fSplineEnergy[i+1]+fSplineEnergy[i]) << 467 << 468 if( x < 0 ) 365 if( x < 0 ) 469 { 366 { 470 x = -x; << 367 x = -x ; 471 } 368 } 472 if( x > fError && fSplineNumber < fMaxSp << 369 if( x > fError && fSplineNumber < fMaxSplineSize-1 ) 473 { 370 { 474 continue; // next division << 371 continue; // next division 475 } 372 } 476 i += 2; // pass to next segment 373 i += 2; // pass to next segment 477 374 478 // Loop checking, 03-Aug-2015, Vladimir << 479 } // close 'while' 375 } // close 'while' 480 376 481 } // end of SplainPAI 377 } // end of SplainPAI 482 378 483 379 484 ////////////////////////////////////////////// 380 //////////////////////////////////////////////////////////////////// 485 // 381 // 486 // Integration over electrons that could be co 382 // Integration over electrons that could be considered 487 // quasi-free at energy transfer of interest 383 // quasi-free at energy transfer of interest 488 384 489 G4double G4PAIySection::RutherfordIntegral( G4 385 G4double G4PAIySection::RutherfordIntegral( G4int k, 490 G4 << 386 G4double x1, 491 << 387 G4double x2 ) 492 { 388 { 493 G4double c1, c2, c3; << 389 G4double c1, c2, c3 ; 494 // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<< << 390 // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl; 495 G4double x12 = x1*x2; << 391 c1 = (x2 - x1)/x1/x2 ; 496 c1 = (x2 - x1)/x12; << 392 c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2 ; 497 c2 = (x2 - x1)*(x2 + x1)/(x12*x12); << 393 c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ; 498 c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/(x12 << 499 // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "< 394 // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<<c2<<"; "<<"c3 = "<<c3<<G4endl; 500 395 501 return fA1[k]*log(x2/x1) + fA2[k]*c1 + fA3 << 396 return fA1[k]*log(x2/x1) + fA2[k]*c1 + fA3[k]*c2/2 + fA4[k]*c3/3 ; 502 397 503 } // end of RutherfordIntegral 398 } // end of RutherfordIntegral 504 399 505 400 506 ////////////////////////////////////////////// 401 ///////////////////////////////////////////////////////////////// 507 // 402 // 508 // Imaginary part of dielectric constant 403 // Imaginary part of dielectric constant 509 // (G4int k - interval number, G4double en1 - 404 // (G4int k - interval number, G4double en1 - energy point) 510 405 511 G4double G4PAIySection::ImPartDielectricConst( << 406 G4double G4PAIySection::ImPartDielectricConst( G4int k , >> 407 G4double energy1 ) 512 { 408 { 513 G4double energy2,energy3,energy4,result; 409 G4double energy2,energy3,energy4,result; 514 410 515 energy2 = energy1*energy1; 411 energy2 = energy1*energy1; 516 energy3 = energy2*energy1; 412 energy3 = energy2*energy1; 517 energy4 = energy3*energy1; 413 energy4 = energy3*energy1; 518 414 519 result = fA1[k]/energy1+fA2[k]/energy2+fA3[ << 415 result = fA1[k]/energy1+fA2[k]/energy2+fA3[k]/energy3+fA4[k]/energy4 ; 520 result *=hbarc/energy1; << 416 result *=hbarc/energy1 ; 521 417 522 return result; << 418 return result ; 523 419 524 } // end of ImPartDielectricConst 420 } // end of ImPartDielectricConst 525 421 526 422 527 ////////////////////////////////////////////// 423 ////////////////////////////////////////////////////////////////////////////// 528 // 424 // 529 // Real part of dielectric constant minus unit 425 // Real part of dielectric constant minus unit: epsilon_1 - 1 530 // (G4double enb - energy point) 426 // (G4double enb - energy point) 531 // 427 // 532 428 533 G4double G4PAIySection::RePartDielectricConst( 429 G4double G4PAIySection::RePartDielectricConst(G4double enb) 534 { 430 { 535 G4double x0, x02, x03, x04, x05, x1, x2, xx 431 G4double x0, x02, x03, x04, x05, x1, x2, xx1 ,xx2 , xx12, 536 c1, c2, c3, cof1, cof2, xln1, xln2 << 432 c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result ; 537 433 538 x0 = enb; << 434 x0 = enb ; 539 result = 0; << 435 result = 0 ; 540 436 541 for(G4int i=1;i<=fIntervalNumber-1;++i) << 437 for(G4int i=1;i<=fIntervalNumber-1;i++) 542 { 438 { 543 x1 = fEnergyInterval[i]; << 439 x1 = fEnergyInterval[i] ; 544 x2 = fEnergyInterval[i+1]; << 440 x2 = fEnergyInterval[i+1] ; 545 xx1 = x1 - x0; << 441 xx1 = x1 - x0 ; 546 xx2 = x2 - x0; << 442 xx2 = x2 - x0 ; 547 xx12 = xx2/xx1; << 443 xx12 = xx2/xx1 ; 548 444 549 if(xx12<0.) << 445 if(xx12<0) 550 { 446 { 551 xx12 = -xx12; << 447 xx12 = -xx12; 552 } 448 } 553 xln1 = log(x2/x1); << 449 xln1 = log(x2/x1) ; 554 xln2 = log(xx12); << 450 xln2 = log(xx12) ; 555 xln3 = log((x2 + x0)/(x1 + x0)); << 451 xln3 = log((x2 + x0)/(x1 + x0)) ; 556 x02 = x0*x0; << 452 x02 = x0*x0 ; 557 x03 = x02*x0; << 453 x03 = x02*x0 ; 558 x04 = x03*x0; << 454 x04 = x03*x0 ; 559 x05 = x04*x0; 455 x05 = x04*x0; 560 G4double x12 = x1*x2; << 456 c1 = (x2 - x1)/x1/x2 ; 561 c1 = (x2 - x1)/x12; << 457 c2 = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2 ; 562 c2 = (x2 - x1)*(x2 +x1)/(x12*x12); << 458 c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ; 563 c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/( << 459 564 << 460 result -= (fA1[i]/x02 + fA3[i]/x04)*xln1 ; 565 result -= (fA1[i]/x02 + fA3[i]/x04)*xln1 << 461 result -= (fA2[i]/x02 + fA4[i]/x04)*c1 ; 566 result -= (fA2[i]/x02 + fA4[i]/x04)*c1; << 462 result -= fA3[i]*c2/2/x02 ; 567 result -= fA3[i]*c2/2/x02; << 463 result -= fA4[i]*c3/3/x02 ; 568 result -= fA4[i]*c3/3/x02; << 569 464 570 cof1 = fA1[i]/x02 + fA3[i]/x04; << 465 cof1 = fA1[i]/x02 + fA3[i]/x04 ; 571 cof2 = fA2[i]/x03 + fA4[i]/x05; << 466 cof2 = fA2[i]/x03 + fA4[i]/x05 ; 572 467 573 result += 0.5*(cof1 +cof2)*xln2; << 468 result += 0.5*(cof1 +cof2)*xln2 ; 574 result += 0.5*(cof1 - cof2)*xln3; << 469 result += 0.5*(cof1 - cof2)*xln3 ; 575 } 470 } 576 result *= 2*hbarc/pi; << 471 result *= 2*hbarc/pi ; 577 472 578 return result; << 473 return result ; 579 474 580 } // end of RePartDielectricConst 475 } // end of RePartDielectricConst 581 476 582 ////////////////////////////////////////////// 477 ////////////////////////////////////////////////////////////////////// 583 // 478 // 584 // PAI differential cross-section in terms of 479 // PAI differential cross-section in terms of 585 // simplified Allison's equation 480 // simplified Allison's equation 586 // 481 // 587 482 588 G4double G4PAIySection::DifPAIySection( G4int 483 G4double G4PAIySection::DifPAIySection( G4int i , 589 G4doub 484 G4double betaGammaSq ) 590 { 485 { 591 G4double beta, be2,cof,x1,x2,x3,x4,x5,x6,x7 << 486 G4double be2,cof,x1,x2,x3,x4,x5,x6,x7,x8,result ; 592 be2 = betaGammaSq/(1 + betaGammaSq); << 487 //G4double beta, be4 ; 593 beta = std::sqrt(be2); << 488 G4double be4 ; 594 cof = 1; << 489 G4double betaBohr2 = fine_structure_const*fine_structure_const ; 595 x1 = log(2*electron_mass_c2/fSplineEnergy[i << 490 G4double betaBohr4 = betaBohr2*betaBohr2*4.0 ; >> 491 be2 = betaGammaSq/(1 + betaGammaSq) ; >> 492 be4 = be2*be2 ; >> 493 // beta = sqrt(be2) ; >> 494 cof = 1 ; >> 495 x1 = log(2*electron_mass_c2/fSplineEnergy[i]) ; 596 496 597 if( betaGammaSq < 0.01 ) x2 = log(be2); << 497 if( betaGammaSq < 0.01 ) x2 = log(be2) ; 598 else 498 else 599 { 499 { 600 x2 = -log( (1/betaGammaSq - fRePartDielec 500 x2 = -log( (1/betaGammaSq - fRePartDielectricConst[i])* 601 (1/betaGammaSq - fRePartDielec << 501 (1/betaGammaSq - fRePartDielectricConst[i]) + 602 fImPartDielectricConst[i]*fImP << 502 fImPartDielectricConst[i]*fImPartDielectricConst[i] )/2 ; 603 } 503 } 604 if( fImPartDielectricConst[i] == 0.0 ||beta 504 if( fImPartDielectricConst[i] == 0.0 ||betaGammaSq < 0.01 ) 605 { 505 { 606 x6=0; << 506 x6=0 ; 607 } 507 } 608 else 508 else 609 { 509 { 610 x3 = -fRePartDielectricConst[i] + 1/betaG << 510 x3 = -fRePartDielectricConst[i] + 1/betaGammaSq ; 611 x5 = -1 - fRePartDielectricConst[i] + 511 x5 = -1 - fRePartDielectricConst[i] + 612 be2*((1 +fRePartDielectricConst[i])* 512 be2*((1 +fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) + 613 fImPartDielectricConst[i]*fImPartDie << 513 fImPartDielectricConst[i]*fImPartDielectricConst[i]) ; 614 514 615 x7 = std::atan2(fImPartDielectricConst[i] << 515 x7 = atan2(fImPartDielectricConst[i],x3) ; 616 x6 = x5 * x7; << 516 x6 = x5 * x7 ; 617 } 517 } 618 x4 = ((x1 + x2)*fImPartDielectricConst[i] + << 518 // if(fImPartDielectricConst[i] == 0) x6 = 0 ; >> 519 >> 520 x4 = ((x1 + x2)*fImPartDielectricConst[i] + x6)/hbarc ; >> 521 // if( x4 < 0.0 ) x4 = 0.0 ; 619 x8 = (1 + fRePartDielectricConst[i])*(1 + f 522 x8 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) + 620 fImPartDielectricConst[i]*fImPartDiele << 523 fImPartDielectricConst[i]*fImPartDielectricConst[i] ; 621 << 622 result = (x4 + cof*fIntegralTerm[i]/fSpline << 623 result = std::max(result, 1.0e-8); << 624 result *= fine_structure_const/(be2*pi); << 625 // low energy correction << 626 << 627 G4double lowCof = fLowEnergyCof; // 6.0 ; / << 628 524 629 result *= (1 - std::exp(-beta/(betaBohr*low << 525 result = (x4 + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i]) ; 630 if(x8 > 0.) << 526 if(result < 1.0e-8) result = 1.0e-8 ; >> 527 result *= fine_structure_const/be2/pi ; >> 528 // result *= (1-exp(-beta/betaBohr))*(1-exp(-beta/betaBohr)) ; >> 529 // result *= (1-exp(-be2/betaBohr2)) ; >> 530 result *= (1-exp(-be4/betaBohr4)) ; >> 531 if(fDensity >= 0.1) 631 { 532 { 632 result /= x8; << 533 result /= x8 ; 633 } 534 } 634 return result; << 535 return result ; 635 536 636 } // end of DifPAIySection 537 } // end of DifPAIySection 637 538 638 ////////////////////////////////////////////// 539 ////////////////////////////////////////////////////////////////////////// 639 // 540 // 640 // Calculation od dN/dx of collisions with cre 541 // Calculation od dN/dx of collisions with creation of Cerenkov pseudo-photons 641 542 642 G4double G4PAIySection::PAIdNdxCerenkov( G4int << 543 G4double G4PAIySection::PAIdNdxCerenkov( G4int i , >> 544 G4double betaGammaSq ) 643 { 545 { 644 G4double logarithm, x3, x5, argument, modul << 546 G4double cof, logarithm, x3, x5, argument, modul2, dNdxC ; 645 G4double be2, be4; << 547 G4double be2, be4, betaBohr2,betaBohr4,cofBetaBohr ; 646 548 647 be2 = betaGammaSq/(1 + betaGammaSq); << 549 cof = 1.0 ; 648 be4 = be2*be2; << 550 cofBetaBohr = 4.0 ; >> 551 betaBohr2 = fine_structure_const*fine_structure_const ; >> 552 betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr ; 649 553 650 if( betaGammaSq < 0.01 ) logarithm = log(1. << 554 be2 = betaGammaSq/(1 + betaGammaSq) ; >> 555 be4 = be2*be2 ; >> 556 >> 557 if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq) ; // 0.0 ; 651 else 558 else 652 { 559 { 653 logarithm = -std::log( (1/betaGammaSq - f << 560 logarithm = -log( (1/betaGammaSq - fRePartDielectricConst[i])* 654 (1/betaGammaSq - fRePa << 561 (1/betaGammaSq - fRePartDielectricConst[i]) + 655 fImPartDielectricConst << 562 fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5 ; 656 logarithm += std::log(1+1.0/betaGammaSq); << 563 logarithm += log(1+1.0/betaGammaSq) ; 657 } 564 } 658 565 659 if( fImPartDielectricConst[i] == 0.0 || bet 566 if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 ) 660 { 567 { 661 argument = 0.0; << 568 argument = 0.0 ; 662 } 569 } 663 else 570 else 664 { 571 { 665 x3 = -fRePartDielectricConst[i] + 1.0/bet << 572 x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq ; 666 x5 = -1.0 - fRePartDielectricConst[i] + 573 x5 = -1.0 - fRePartDielectricConst[i] + 667 be2*((1.0 +fRePartDielectricConst[i] 574 be2*((1.0 +fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) + 668 fImPartDielectricConst[i]*fImPartDie << 575 fImPartDielectricConst[i]*fImPartDielectricConst[i]) ; 669 if( x3 == 0.0 ) argument = 0.5*pi; 576 if( x3 == 0.0 ) argument = 0.5*pi; 670 else argument = std::atan2(fIm << 577 else argument = atan2(fImPartDielectricConst[i],x3) ; 671 argument *= x5 ; << 578 argument *= x5 ; 672 } 579 } 673 dNdxC = ( logarithm*fImPartDielectricConst[ << 580 dNdxC = ( logarithm*fImPartDielectricConst[i] + argument )/hbarc ; 674 581 675 if(dNdxC < 1.0e-8) dNdxC = 1.0e-8; << 582 if(dNdxC < 1.0e-8) dNdxC = 1.0e-8 ; 676 583 677 dNdxC *= fine_structure_const/be2/pi; << 584 dNdxC *= fine_structure_const/be2/pi ; 678 585 679 dNdxC *= (1 - std::exp(-be4/betaBohr4)); << 586 dNdxC *= (1-exp(-be4/betaBohr4)) ; 680 587 681 modul2 = (1.0 + fRePartDielectricConst[i])* << 588 if(fDensity >= 0.1) 682 fImPartDielectricConst[i]* << 589 { 683 if(modul2 > 0.) << 590 modul2 = (1.0 + fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) + 684 { << 591 fImPartDielectricConst[i]*fImPartDielectricConst[i] ; 685 dNdxC /= modul2; << 592 dNdxC /= modul2 ; 686 } << 593 } 687 return dNdxC; << 594 return dNdxC ; 688 595 689 } // end of PAIdNdxCerenkov 596 } // end of PAIdNdxCerenkov 690 597 691 ////////////////////////////////////////////// 598 ////////////////////////////////////////////////////////////////////////// 692 // 599 // 693 // Calculation od dN/dx of collisions with cre 600 // Calculation od dN/dx of collisions with creation of longitudinal EM 694 // excitations (plasmons, delta-electrons) 601 // excitations (plasmons, delta-electrons) 695 602 696 G4double G4PAIySection::PAIdNdxPlasmon( G4int << 603 G4double G4PAIySection::PAIdNdxPlasmon( G4int i , >> 604 G4double betaGammaSq ) 697 { 605 { 698 G4double cof, resonance, modul2, dNdxP; << 606 G4double cof, resonance, modul2, dNdxP ; 699 G4double be2, be4; << 607 G4double be2, be4, betaBohr2, betaBohr4, cofBetaBohr ; 700 608 701 cof = 1; << 609 cof = 1 ; >> 610 cofBetaBohr = 4.0 ; >> 611 betaBohr2 = fine_structure_const*fine_structure_const ; >> 612 betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr ; 702 613 703 be2 = betaGammaSq/(1 + betaGammaSq); << 614 be2 = betaGammaSq/(1 + betaGammaSq) ; 704 be4 = be2*be2; << 615 be4 = be2*be2 ; 705 616 706 resonance = std::log(2*electron_mass_c2*be2 << 617 resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]) ; 707 resonance *= fImPartDielectricConst[i]/hbar << 618 resonance *= fImPartDielectricConst[i]/hbarc ; 708 619 709 dNdxP = ( resonance + cof*fIntegralTerm[i]/ << 710 620 711 dNdxP = std::max(dNdxP, 1.0e-8); << 621 dNdxP = ( resonance + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i] ) ; 712 622 713 dNdxP *= fine_structure_const/be2/pi; << 623 if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8 ; 714 dNdxP *= (1 - std::exp(-be4/betaBohr4)); << 715 624 716 modul2 = (1 + fRePartDielectricConst[i])*(1 << 625 dNdxP *= fine_structure_const/be2/pi ; 717 fImPartDielectricConst[i]*fImPartDielectr << 626 dNdxP *= (1-exp(-be4/betaBohr4)) ; 718 if(modul2 > 0.) << 627 719 { << 628 if( fDensity >= 0.1 ) 720 dNdxP /= modul2; << 629 { 721 } << 630 modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) + 722 return dNdxP; << 631 fImPartDielectricConst[i]*fImPartDielectricConst[i] ; >> 632 dNdxP /= modul2 ; >> 633 } >> 634 return dNdxP ; 723 635 724 } // end of PAIdNdxPlasmon 636 } // end of PAIdNdxPlasmon 725 637 726 ////////////////////////////////////////////// 638 //////////////////////////////////////////////////////////////////////// 727 // 639 // 728 // Calculation of the PAI integral cross-secti 640 // Calculation of the PAI integral cross-section 729 // fIntegralPAIySection[1] = specific primary 641 // fIntegralPAIySection[1] = specific primary ionisation, 1/cm 730 // and fIntegralPAIySection[0] = mean energy l 642 // and fIntegralPAIySection[0] = mean energy loss per cm in keV/cm 731 643 732 void G4PAIySection::IntegralPAIySection() 644 void G4PAIySection::IntegralPAIySection() 733 { 645 { 734 fIntegralPAIySection[fSplineNumber] = 0; << 646 fIntegralPAIySection[fSplineNumber] = 0 ; 735 fIntegralPAIdEdx[fSplineNumber] = 0; << 647 fIntegralPAIdEdx[fSplineNumber] = 0 ; 736 fIntegralPAIySection[0] = 0; << 648 fIntegralPAIySection[0] = 0 ; 737 G4int k = fIntervalNumber -1; << 649 G4int k = fIntervalNumber -1 ; 738 650 739 for(G4int i = fSplineNumber-1; i >= 1; i--) << 651 for(G4int i = fSplineNumber-1 ; i >= 1 ; i--) 740 { 652 { 741 if(fSplineEnergy[i] >= fEnergyInterval[k]) 653 if(fSplineEnergy[i] >= fEnergyInterval[k]) 742 { 654 { 743 fIntegralPAIySection[i] = fIntegralPAIyS << 655 fIntegralPAIySection[i] = fIntegralPAIySection[i+1] + SumOverInterval(i) ; 744 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i << 656 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] + SumOverIntervaldEdx(i) ; 745 } 657 } 746 else 658 else 747 { 659 { 748 fIntegralPAIySection[i] = fIntegralPAIyS 660 fIntegralPAIySection[i] = fIntegralPAIySection[i+1] + 749 SumOverBord << 661 SumOverBorder(i+1,fEnergyInterval[k]) ; 750 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i 662 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] + 751 SumOverBord << 663 SumOverBorderdEdx(i+1,fEnergyInterval[k]) ; 752 k--; << 664 k-- ; 753 } 665 } 754 } 666 } 755 } // end of IntegralPAIySection 667 } // end of IntegralPAIySection 756 668 757 ////////////////////////////////////////////// 669 //////////////////////////////////////////////////////////////////////// 758 // 670 // 759 // Calculation of the PAI Cerenkov integral cr 671 // Calculation of the PAI Cerenkov integral cross-section 760 // fIntegralCrenkov[1] = specific Crenkov ioni 672 // fIntegralCrenkov[1] = specific Crenkov ionisation, 1/cm 761 // and fIntegralCerenkov[0] = mean Cerenkov lo 673 // and fIntegralCerenkov[0] = mean Cerenkov loss per cm in keV/cm 762 674 763 void G4PAIySection::IntegralCerenkov() 675 void G4PAIySection::IntegralCerenkov() 764 { 676 { 765 G4int i, k; << 677 G4int i, k ; 766 fIntegralCerenkov[fSplineNumber] = 0; << 678 fIntegralCerenkov[fSplineNumber] = 0 ; 767 fIntegralCerenkov[0] = 0; << 679 fIntegralCerenkov[0] = 0 ; 768 k = fIntervalNumber -1; << 680 k = fIntervalNumber -1 ; 769 681 770 for( i = fSplineNumber-1; i >= 1; i-- ) << 682 for( i = fSplineNumber-1 ; i >= 1 ; i-- ) 771 { 683 { 772 if(fSplineEnergy[i] >= fEnergyInterval[k 684 if(fSplineEnergy[i] >= fEnergyInterval[k]) 773 { 685 { 774 fIntegralCerenkov[i] = fIntegralCerenk << 686 fIntegralCerenkov[i] = fIntegralCerenkov[i+1] + SumOverInterCerenkov(i) ; 775 // G4cout<<"int: i = "<<i<<"; sumC = " << 687 // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl; 776 } 688 } 777 else 689 else 778 { 690 { 779 fIntegralCerenkov[i] = fIntegralCerenk 691 fIntegralCerenkov[i] = fIntegralCerenkov[i+1] + 780 SumOverBord << 692 SumOverBordCerenkov(i+1,fEnergyInterval[k]) ; 781 k--; << 693 k-- ; 782 // G4cout<<"bord: i = "<<i<<"; sumC = << 694 // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl; 783 } 695 } 784 } 696 } 785 697 786 } // end of IntegralCerenkov 698 } // end of IntegralCerenkov 787 699 788 ////////////////////////////////////////////// 700 //////////////////////////////////////////////////////////////////////// 789 // 701 // 790 // Calculation of the PAI Plasmon integral cro 702 // Calculation of the PAI Plasmon integral cross-section 791 // fIntegralPlasmon[1] = splasmon primary ioni 703 // fIntegralPlasmon[1] = splasmon primary ionisation, 1/cm 792 // and fIntegralPlasmon[0] = mean plasmon loss 704 // and fIntegralPlasmon[0] = mean plasmon loss per cm in keV/cm 793 705 794 void G4PAIySection::IntegralPlasmon() 706 void G4PAIySection::IntegralPlasmon() 795 { 707 { 796 fIntegralPlasmon[fSplineNumber] = 0; << 708 fIntegralPlasmon[fSplineNumber] = 0 ; 797 fIntegralPlasmon[0] = 0; << 709 fIntegralPlasmon[0] = 0 ; 798 G4int k = fIntervalNumber -1; << 710 G4int k = fIntervalNumber -1 ; 799 for(G4int i=fSplineNumber-1;i>=1;i--) 711 for(G4int i=fSplineNumber-1;i>=1;i--) 800 { 712 { 801 if(fSplineEnergy[i] >= fEnergyInterval[k 713 if(fSplineEnergy[i] >= fEnergyInterval[k]) 802 { 714 { 803 fIntegralPlasmon[i] = fIntegralPlasmon << 715 fIntegralPlasmon[i] = fIntegralPlasmon[i+1] + SumOverInterPlasmon(i) ; 804 } 716 } 805 else 717 else 806 { 718 { 807 fIntegralPlasmon[i] = fIntegralPlasmon 719 fIntegralPlasmon[i] = fIntegralPlasmon[i+1] + 808 SumOverBord << 720 SumOverBordPlasmon(i+1,fEnergyInterval[k]) ; 809 k--; << 721 k-- ; 810 } 722 } 811 } 723 } >> 724 812 } // end of IntegralPlasmon 725 } // end of IntegralPlasmon 813 726 814 ////////////////////////////////////////////// 727 ////////////////////////////////////////////////////////////////////// 815 // 728 // 816 // Calculation the PAI integral cross-section 729 // Calculation the PAI integral cross-section inside 817 // of interval of continuous values of photo-i 730 // of interval of continuous values of photo-ionisation 818 // cross-section. Parameter 'i' is the number 731 // cross-section. Parameter 'i' is the number of interval. 819 732 820 G4double G4PAIySection::SumOverInterval( G4int 733 G4double G4PAIySection::SumOverInterval( G4int i ) 821 { 734 { 822 G4double x0,x1,y0,yy1,a,b,c,result; << 735 G4double x0,x1,y0,yy1,a,b,c,result ; 823 << 824 x0 = fSplineEnergy[i]; << 825 x1 = fSplineEnergy[i+1]; << 826 736 827 if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) << 737 x0 = fSplineEnergy[i] ; 828 << 738 x1 = fSplineEnergy[i+1] ; 829 y0 = fDifPAIySection[i]; << 739 y0 = fDifPAIySection[i] ; 830 yy1 = fDifPAIySection[i+1]; 740 yy1 = fDifPAIySection[i+1]; 831 //G4cout << "## x0= " << x0 << " x1= " << x << 832 c = x1/x0; 741 c = x1/x0; 833 //G4cout << "c= " << c << " y0= " << y0 << << 742 a = log10(yy1/y0)/log10(c) ; 834 a = log10(yy1/y0)/log10(c); << 743 // b = log10(y0) - a*log10(x0) ; 835 //G4cout << "a= " << a << G4endl; << 744 b = y0/pow(x0,a) ; 836 << 745 a += 1 ; 837 b = 0.0; << 838 if(a < 20.) b = y0/pow(x0,a); << 839 << 840 a += 1; << 841 if(a == 0) 746 if(a == 0) 842 { 747 { 843 result = b*log(x1/x0); << 748 result = b*log(x1/x0) ; 844 } 749 } 845 else 750 else 846 { 751 { 847 result = y0*(x1*pow(c,a-1) - x0)/a; << 752 result = y0*(x1*pow(c,a-1) - x0)/a ; 848 } 753 } 849 a++; 754 a++; 850 if(a == 0) 755 if(a == 0) 851 { 756 { 852 fIntegralPAIySection[0] += b*log(x1/x0); << 757 fIntegralPAIySection[0] += b*log(x1/x0) ; 853 } 758 } 854 else 759 else 855 { 760 { 856 fIntegralPAIySection[0] += y0*(x1*x1*pow << 761 fIntegralPAIySection[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a ; 857 } 762 } 858 return result; << 763 return result ; 859 764 860 } // end of SumOverInterval 765 } // end of SumOverInterval 861 766 862 ///////////////////////////////// 767 ///////////////////////////////// 863 768 864 G4double G4PAIySection::SumOverIntervaldEdx( G 769 G4double G4PAIySection::SumOverIntervaldEdx( G4int i ) 865 { 770 { 866 G4double x0,x1,y0,yy1,a,b,c,result; << 771 G4double x0,x1,y0,yy1,a,b,c,result ; 867 << 868 x0 = fSplineEnergy[i]; << 869 x1 = fSplineEnergy[i+1]; << 870 772 871 if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) << 773 x0 = fSplineEnergy[i] ; 872 << 774 x1 = fSplineEnergy[i+1] ; 873 y0 = fDifPAIySection[i]; << 775 y0 = fDifPAIySection[i] ; 874 yy1 = fDifPAIySection[i+1]; 776 yy1 = fDifPAIySection[i+1]; 875 c = x1/x0; 777 c = x1/x0; 876 a = log10(yy1/y0)/log10(c); << 778 a = log10(yy1/y0)/log10(c) ; 877 << 779 // b = log10(y0) - a*log10(x0) ; 878 b = 0.0; << 780 b = y0/pow(x0,a) ; 879 if(a < 20.) b = y0/pow(x0,a); << 781 a += 2 ; 880 << 881 a += 2; << 882 if(a == 0) 782 if(a == 0) 883 { 783 { 884 result = b*log(x1/x0); << 784 result = b*log(x1/x0) ; 885 } 785 } 886 else 786 else 887 { 787 { 888 result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a; << 788 result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a ; 889 } 789 } 890 return result; << 790 return result ; 891 791 892 } // end of SumOverInterval 792 } // end of SumOverInterval 893 793 894 ////////////////////////////////////////////// 794 ////////////////////////////////////////////////////////////////////// 895 // 795 // 896 // Calculation the PAI Cerenkov integral cross 796 // Calculation the PAI Cerenkov integral cross-section inside 897 // of interval of continuous values of photo-i 797 // of interval of continuous values of photo-ionisation Cerenkov 898 // cross-section. Parameter 'i' is the number 798 // cross-section. Parameter 'i' is the number of interval. 899 799 900 G4double G4PAIySection::SumOverInterCerenkov( 800 G4double G4PAIySection::SumOverInterCerenkov( G4int i ) 901 { 801 { 902 G4double x0,x1,y0,yy1,a,c,result; << 802 G4double x0,x1,y0,yy1,a,b,c,result ; 903 << 904 x0 = fSplineEnergy[i]; << 905 x1 = fSplineEnergy[i+1]; << 906 803 907 if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) << 804 x0 = fSplineEnergy[i] ; 908 << 805 x1 = fSplineEnergy[i+1] ; 909 y0 = fdNdxCerenkov[i]; << 806 y0 = fdNdxCerenkov[i] ; 910 yy1 = fdNdxCerenkov[i+1]; 807 yy1 = fdNdxCerenkov[i+1]; 911 // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<" 808 // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1 912 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4en 809 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl; 913 810 914 c = x1/x0; 811 c = x1/x0; 915 a = log10(yy1/y0)/log10(c); << 812 a = log10(yy1/y0)/log10(c) ; 916 G4double b = 0.0; << 813 b = y0/pow(x0,a) ; 917 if(a < 20.) b = y0/pow(x0,a); << 814 918 << 815 a += 1.0 ; 919 a += 1.0; << 816 if(a == 0) result = b*log(c) ; 920 if(a == 0) result = b*log(c); << 817 else result = y0*(x1*pow(c,a-1) - x0)/a ; 921 else result = y0*(x1*pow(c,a-1) - x0) << 818 a += 1.0 ; 922 a += 1.0; << 923 819 924 if( a == 0 ) fIntegralCerenkov[0] += b*log( << 820 if( a == 0 ) fIntegralCerenkov[0] += b*log(x1/x0) ; 925 else fIntegralCerenkov[0] += y0*(x1 << 821 else fIntegralCerenkov[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a ; 926 // G4cout<<"a = "<<a<<"; b = "<<b<<"; resu 822 // G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl; 927 return result; << 823 return result ; 928 824 929 } // end of SumOverInterCerenkov 825 } // end of SumOverInterCerenkov 930 826 931 ////////////////////////////////////////////// 827 ////////////////////////////////////////////////////////////////////// 932 // 828 // 933 // Calculation the PAI Plasmon integral cross- 829 // Calculation the PAI Plasmon integral cross-section inside 934 // of interval of continuous values of photo-i 830 // of interval of continuous values of photo-ionisation Plasmon 935 // cross-section. Parameter 'i' is the number 831 // cross-section. Parameter 'i' is the number of interval. 936 832 937 G4double G4PAIySection::SumOverInterPlasmon( G 833 G4double G4PAIySection::SumOverInterPlasmon( G4int i ) 938 { 834 { 939 G4double x0,x1,y0,yy1,a,c,result; << 835 G4double x0,x1,y0,yy1,a,b,c,result ; 940 836 941 x0 = fSplineEnergy[i]; << 837 x0 = fSplineEnergy[i] ; 942 x1 = fSplineEnergy[i+1]; << 838 x1 = fSplineEnergy[i+1] ; 943 << 839 y0 = fdNdxPlasmon[i] ; 944 if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) << 945 << 946 y0 = fdNdxPlasmon[i]; << 947 yy1 = fdNdxPlasmon[i+1]; 840 yy1 = fdNdxPlasmon[i+1]; 948 c = x1/x0; << 841 c =x1/x0; 949 a = log10(yy1/y0)/log10(c); << 842 a = log10(yy1/y0)/log10(c) ; 950 << 843 // b = log10(y0) - a*log10(x0) ; 951 G4double b = 0.0; << 844 b = y0/pow(x0,a) ; 952 if(a < 20.) b = y0/pow(x0,a); << 845 >> 846 a += 1.0 ; >> 847 if(a == 0) result = b*log(x1/x0) ; >> 848 else result = y0*(x1*pow(c,a-1) - x0)/a ; >> 849 a += 1.0 ; 953 850 954 a += 1.0; << 851 if( a == 0 ) fIntegralPlasmon[0] += b*log(x1/x0) ; 955 if(a == 0) result = b*log(x1/x0); << 852 else fIntegralPlasmon[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a ; 956 else result = y0*(x1*pow(c,a-1) - x0) << 957 a += 1.0; << 958 << 959 if( a == 0 ) fIntegralPlasmon[0] += b*log(x << 960 else fIntegralPlasmon[0] += y0*(x1* << 961 853 962 return result; << 854 return result ; 963 855 964 } // end of SumOverInterPlasmon 856 } // end of SumOverInterPlasmon 965 857 966 ////////////////////////////////////////////// 858 /////////////////////////////////////////////////////////////////////////////// 967 // 859 // 968 // Integration of PAI cross-section for the ca 860 // Integration of PAI cross-section for the case of 969 // passing across border between intervals 861 // passing across border between intervals 970 862 971 G4double G4PAIySection::SumOverBorder( G4int 863 G4double G4PAIySection::SumOverBorder( G4int i , 972 G4doubl 864 G4double en0 ) 973 { 865 { 974 G4double x0,x1,y0,yy1,a,d,e0,result; << 866 G4double x0,x1,y0,yy1,a,b,c,d,e0,result ; 975 867 976 e0 = en0; << 868 e0 = en0 ; 977 x0 = fSplineEnergy[i]; << 869 x0 = fSplineEnergy[i] ; 978 x1 = fSplineEnergy[i+1]; << 870 x1 = fSplineEnergy[i+1] ; 979 y0 = fDifPAIySection[i]; << 871 y0 = fDifPAIySection[i] ; 980 yy1 = fDifPAIySection[i+1]; << 872 yy1 = fDifPAIySection[i+1] ; 981 873 >> 874 c = x1/x0; 982 d = e0/x0; 875 d = e0/x0; 983 a = log10(yy1/y0)/log10(x1/x0); << 876 a = log10(yy1/y0)/log10(x1/x0) ; 984 << 877 // b0 = log10(y0) - a*log10(x0) ; 985 G4double b = 0.0; << 878 b = y0/pow(x0,a); // pow(10.,b) ; 986 if(a < 20.) b = y0/pow(x0,a); << 987 879 988 a += 1; << 880 a += 1 ; 989 if(a == 0) 881 if(a == 0) 990 { 882 { 991 result = b*log(x0/e0); << 883 result = b*log(x0/e0) ; 992 } 884 } 993 else 885 else 994 { 886 { 995 result = y0*(x0 - e0*pow(d,a-1))/a; << 887 result = y0*(x0 - e0*pow(d,a-1))/a ; 996 } 888 } 997 a++; << 889 a++ ; 998 if(a == 0) 890 if(a == 0) 999 { 891 { 1000 fIntegralPAIySection[0] += b*log(x0/e0) << 892 fIntegralPAIySection[0] += b*log(x0/e0) ; 1001 } 893 } 1002 else 894 else 1003 { 895 { 1004 fIntegralPAIySection[0] += y0*(x0*x0 - << 896 fIntegralPAIySection[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a ; 1005 } 897 } 1006 x0 = fSplineEnergy[i - 1]; << 898 x0 = fSplineEnergy[i - 1] ; 1007 x1 = fSplineEnergy[i - 2]; << 899 x1 = fSplineEnergy[i - 2] ; 1008 y0 = fDifPAIySection[i - 1]; << 900 y0 = fDifPAIySection[i - 1] ; 1009 yy1 = fDifPAIySection[i - 2]; << 901 yy1 = fDifPAIySection[i - 2] ; 1010 902 1011 //c = x1/x0; << 903 c = x1/x0; 1012 d = e0/x0; 904 d = e0/x0; 1013 a = log10(yy1/y0)/log10(x1/x0); << 905 a = log10(yy1/y0)/log10(x1/x0) ; 1014 << 906 // b0 = log10(y0) - a*log10(x0) ; 1015 b = 0.0; << 907 b = y0/pow(x0,a) ; 1016 if(a < 20.) b = y0/pow(x0,a); << 908 a += 1 ; 1017 << 1018 a += 1; << 1019 if(a == 0) 909 if(a == 0) 1020 { 910 { 1021 result += b*log(e0/x0); << 911 result += b*log(e0/x0) ; 1022 } 912 } 1023 else 913 else 1024 { 914 { 1025 result += y0*(e0*pow(d,a-1) - x0)/a; << 915 result += y0*(e0*pow(d,a-1) - x0)/a ; 1026 } 916 } 1027 a++; << 917 a++ ; 1028 if(a == 0) 918 if(a == 0) 1029 { 919 { 1030 fIntegralPAIySection[0] += b*log(e0/x0) << 920 fIntegralPAIySection[0] += b*log(e0/x0) ; 1031 } 921 } 1032 else 922 else 1033 { 923 { 1034 fIntegralPAIySection[0] += y0*(e0*e0*po << 924 fIntegralPAIySection[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a ; 1035 } 925 } 1036 return result; << 926 return result ; 1037 927 1038 } 928 } 1039 929 1040 ///////////////////////////////////////////// 930 /////////////////////////////////////////////////////////////////////// 1041 931 1042 G4double G4PAIySection::SumOverBorderdEdx( G4 932 G4double G4PAIySection::SumOverBorderdEdx( G4int i , 1043 G4doub 933 G4double en0 ) 1044 { 934 { 1045 G4double x0,x1,y0,yy1,a,/*c,*/d,e0,result; << 935 G4double x0,x1,y0,yy1,a,b,c,d,e0,result ; 1046 936 1047 e0 = en0; << 937 e0 = en0 ; 1048 x0 = fSplineEnergy[i]; << 938 x0 = fSplineEnergy[i] ; 1049 x1 = fSplineEnergy[i+1]; << 939 x1 = fSplineEnergy[i+1] ; 1050 y0 = fDifPAIySection[i]; << 940 y0 = fDifPAIySection[i] ; 1051 yy1 = fDifPAIySection[i+1]; << 941 yy1 = fDifPAIySection[i+1] ; 1052 942 >> 943 c = x1/x0; 1053 d = e0/x0; 944 d = e0/x0; 1054 a = log10(yy1/y0)/log10(x1/x0); << 945 a = log10(yy1/y0)/log10(x1/x0) ; 1055 << 946 // b0 = log10(y0) - a*log10(x0) ; 1056 G4double b = 0.0; << 947 b = y0/pow(x0,a); // pow(10.,b) ; 1057 if(a < 20.) b = y0/pow(x0,a); << 1058 948 1059 a += 2; << 949 a += 2 ; 1060 if(a == 0) 950 if(a == 0) 1061 { 951 { 1062 result = b*log(x0/e0); << 952 result = b*log(x0/e0) ; 1063 } 953 } 1064 else 954 else 1065 { 955 { 1066 result = y0*(x0*x0 - e0*e0*pow(d,a-2))/ << 956 result = y0*(x0*x0 - e0*e0*pow(d,a-2))/a ; 1067 } 957 } 1068 x0 = fSplineEnergy[i - 1]; << 958 x0 = fSplineEnergy[i - 1] ; 1069 x1 = fSplineEnergy[i - 2]; << 959 x1 = fSplineEnergy[i - 2] ; 1070 y0 = fDifPAIySection[i - 1]; << 960 y0 = fDifPAIySection[i - 1] ; 1071 yy1 = fDifPAIySection[i - 2]; << 961 yy1 = fDifPAIySection[i - 2] ; 1072 962 >> 963 c = x1/x0; 1073 d = e0/x0; 964 d = e0/x0; 1074 a = log10(yy1/y0)/log10(x1/x0); << 965 a = log10(yy1/y0)/log10(x1/x0) ; 1075 << 966 // b0 = log10(y0) - a*log10(x0) ; 1076 b = 0.0; << 967 b = y0/pow(x0,a) ; 1077 if(a < 20.) b = y0/pow(x0,a); << 968 a += 2 ; 1078 << 1079 a += 2; << 1080 if(a == 0) 969 if(a == 0) 1081 { 970 { 1082 result += b*log(e0/x0); << 971 result += b*log(e0/x0) ; 1083 } 972 } 1084 else 973 else 1085 { 974 { 1086 result += y0*(e0*e0*pow(d,a-2) - x0*x0) << 975 result += y0*(e0*e0*pow(d,a-2) - x0*x0)/a ; 1087 } 976 } 1088 return result; << 977 return result ; >> 978 1089 } 979 } 1090 980 1091 ///////////////////////////////////////////// 981 /////////////////////////////////////////////////////////////////////////////// 1092 // 982 // 1093 // Integration of Cerenkov cross-section for 983 // Integration of Cerenkov cross-section for the case of 1094 // passing across border between intervals 984 // passing across border between intervals 1095 985 1096 G4double G4PAIySection::SumOverBordCerenkov( 986 G4double G4PAIySection::SumOverBordCerenkov( G4int i , 1097 987 G4double en0 ) 1098 { 988 { 1099 G4double x0,x1,y0,yy1,a,e0,c,d,result; << 989 G4double x0,x1,y0,yy1,a,b,e0,c,d,result ; 1100 990 1101 e0 = en0; << 991 e0 = en0 ; 1102 x0 = fSplineEnergy[i]; << 992 x0 = fSplineEnergy[i] ; 1103 x1 = fSplineEnergy[i+1]; << 993 x1 = fSplineEnergy[i+1] ; 1104 y0 = fdNdxCerenkov[i]; << 994 y0 = fdNdxCerenkov[i] ; 1105 yy1 = fdNdxCerenkov[i+1]; << 995 yy1 = fdNdxCerenkov[i+1] ; 1106 996 1107 // G4cout<<G4endl; 997 // G4cout<<G4endl; 1108 //G4cout<<"SumBordC, i = "<<i<<"; en0 = "< << 998 // G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1 1109 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G 999 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl; 1110 c = x1/x0; << 1000 c = x1/x0 ; 1111 d = e0/x0; << 1001 d = e0/x0 ; 1112 a = log10(yy1/y0)/log10(c); << 1002 a = log10(yy1/y0)/log10(c) ; 1113 << 1003 // b0 = log10(y0) - a*log10(x0) ; 1114 G4double b = 0.0; << 1004 b = y0/pow(x0,a); // pow(10.,b0) ; 1115 if(a < 20.) b = y0/pow(x0,a); << 1116 1005 1117 a += 1.0; << 1006 a += 1.0 ; 1118 if( a == 0 ) result = b*log(x0/e0); << 1007 if( a == 0 ) result = b*log(x0/e0) ; 1119 else result = y0*(x0 - e0*pow(d,a- << 1008 else result = y0*(x0 - e0*pow(d,a-1))/a ; 1120 a += 1.0; << 1009 a += 1.0 ; 1121 1010 1122 if( a == 0 ) fIntegralCerenkov[0] += b*log << 1011 if( a == 0 ) fIntegralCerenkov[0] += b*log(x0/e0) ; 1123 else fIntegralCerenkov[0] += y0*(x << 1012 else fIntegralCerenkov[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a ; 1124 1013 1125 //G4cout<<"a = "<<a<<"; b = "<<b<<"; resul << 1014 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "<<b<<"; result = "<<result<<G4endl; 1126 1015 1127 x0 = fSplineEnergy[i - 1]; << 1016 x0 = fSplineEnergy[i - 1] ; 1128 x1 = fSplineEnergy[i - 2]; << 1017 x1 = fSplineEnergy[i - 2] ; 1129 y0 = fdNdxCerenkov[i - 1]; << 1018 y0 = fdNdxCerenkov[i - 1] ; 1130 yy1 = fdNdxCerenkov[i - 2]; << 1019 yy1 = fdNdxCerenkov[i - 2] ; 1131 1020 1132 //G4cout<<"x0 ="<<x0<<"; x1 = "<<x1 << 1021 // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1 1133 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4 1022 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl; 1134 1023 1135 c = x1/x0; << 1024 c = x1/x0 ; 1136 d = e0/x0; << 1025 d = e0/x0 ; 1137 a = log10(yy1/y0)/log10(x1/x0); << 1026 a = log10(yy1/y0)/log10(x1/x0) ; 1138 << 1027 // b0 = log10(y0) - a*log10(x0) ; 1139 // G4cout << "a= " << a << G4endl; << 1028 b = y0/pow(x0,a); // pow(10.,b0) ; 1140 if(a > 20.0) b = 0.0; << 1029 1141 else b = y0/pow(x0,a); << 1030 a += 1.0 ; 1142 << 1031 if( a == 0 ) result += b*log(e0/x0) ; 1143 //G4cout << "b= " << b << G4endl; << 1032 else result += y0*(e0*pow(d,a-1) - x0 )/a ; >> 1033 a += 1.0 ; 1144 1034 1145 a += 1.0; << 1035 if( a == 0 ) fIntegralCerenkov[0] += b*log(e0/x0) ; 1146 if( a == 0 ) result += b*log(e0/x0); << 1036 else fIntegralCerenkov[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a ; 1147 else result += y0*(e0*pow(d,a-1) - << 1148 a += 1.0; << 1149 //G4cout << "result= " << result << G4endl << 1150 1037 1151 if( a == 0 ) fIntegralCerenkov[0] += b*l << 1038 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = " 1152 else fIntegralCerenkov[0] += y0* << 1039 // <<b<<"; result = "<<result<<G4endl; 1153 1040 1154 //G4cout<<"a = "<<a<<"; b = "<<b<<"; resul << 1041 return result ; 1155 1042 1156 return result; << 1157 } 1043 } 1158 1044 1159 ///////////////////////////////////////////// 1045 /////////////////////////////////////////////////////////////////////////////// 1160 // 1046 // 1161 // Integration of Plasmon cross-section for t 1047 // Integration of Plasmon cross-section for the case of 1162 // passing across border between intervals 1048 // passing across border between intervals 1163 1049 1164 G4double G4PAIySection::SumOverBordPlasmon( G 1050 G4double G4PAIySection::SumOverBordPlasmon( G4int i , 1165 1051 G4double en0 ) 1166 { 1052 { 1167 G4double x0,x1,y0,yy1,a,c,d,e0,result; << 1053 G4double x0,x1,y0,yy1,a,b,c,d,e0,result ; 1168 << 1169 e0 = en0; << 1170 x0 = fSplineEnergy[i]; << 1171 x1 = fSplineEnergy[i+1]; << 1172 y0 = fdNdxPlasmon[i]; << 1173 yy1 = fdNdxPlasmon[i+1]; << 1174 1054 1175 c = x1/x0; << 1055 e0 = en0 ; 1176 d = e0/x0; << 1056 x0 = fSplineEnergy[i] ; 1177 a = log10(yy1/y0)/log10(c); << 1057 x1 = fSplineEnergy[i+1] ; 1178 << 1058 y0 = fdNdxPlasmon[i] ; 1179 G4double b = 0.0; << 1059 yy1 = fdNdxPlasmon[i+1] ; 1180 if(a < 20.) b = y0/pow(x0,a); << 1060 >> 1061 c = x1/x0 ; >> 1062 d = e0/x0 ; >> 1063 a = log10(yy1/y0)/log10(c) ; >> 1064 // b0 = log10(y0) - a*log10(x0) ; >> 1065 b = y0/pow(x0,a); //pow(10.,b) ; 1181 1066 1182 a += 1.0; << 1067 a += 1.0 ; 1183 if( a == 0 ) result = b*log(x0/e0); << 1068 if( a == 0 ) result = b*log(x0/e0) ; 1184 else result = y0*(x0 - e0*pow(d,a- << 1069 else result = y0*(x0 - e0*pow(d,a-1))/a ; 1185 a += 1.0; << 1070 a += 1.0 ; 1186 1071 1187 if( a == 0 ) fIntegralPlasmon[0] += b*log( << 1072 if( a == 0 ) fIntegralPlasmon[0] += b*log(x0/e0) ; 1188 else fIntegralPlasmon[0] += y0*(x0 << 1073 else fIntegralPlasmon[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a ; 1189 1074 1190 x0 = fSplineEnergy[i - 1]; << 1075 x0 = fSplineEnergy[i - 1] ; 1191 x1 = fSplineEnergy[i - 2]; << 1076 x1 = fSplineEnergy[i - 2] ; 1192 y0 = fdNdxPlasmon[i - 1]; << 1077 y0 = fdNdxPlasmon[i - 1] ; 1193 yy1 = fdNdxPlasmon[i - 2]; << 1078 yy1 = fdNdxPlasmon[i - 2] ; 1194 << 1079 1195 c = x1/x0; << 1080 c = x1/x0 ; 1196 d = e0/x0; << 1081 d = e0/x0 ; 1197 a = log10(yy1/y0)/log10(c); << 1082 a = log10(yy1/y0)/log10(c) ; 1198 << 1083 // b0 = log10(y0) - a*log10(x0) ; 1199 if(a < 20.) b = y0/pow(x0,a); << 1084 b = y0/pow(x0,a);// pow(10.,b0) ; >> 1085 >> 1086 a += 1.0 ; >> 1087 if( a == 0 ) result += b*log(e0/x0) ; >> 1088 else result += y0*(e0*pow(d,a-1) - x0)/a ; >> 1089 a += 1.0 ; 1200 1090 1201 a += 1.0; << 1091 if( a == 0 ) fIntegralPlasmon[0] += b*log(e0/x0) ; 1202 if( a == 0 ) result += b*log(e0/x0); << 1092 else fIntegralPlasmon[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a ; 1203 else result += y0*(e0*pow(d,a-1) - << 1204 a += 1.0; << 1205 << 1206 if( a == 0 ) fIntegralPlasmon[0] += b*lo << 1207 else fIntegralPlasmon[0] += y0*( << 1208 1093 1209 return result; << 1094 return result ; 1210 1095 1211 } 1096 } 1212 1097 1213 ///////////////////////////////////////////// 1098 ///////////////////////////////////////////////////////////////////////// 1214 // 1099 // 1215 // 1100 // 1216 1101 1217 G4double G4PAIySection::GetStepEnergyLoss( G4 1102 G4double G4PAIySection::GetStepEnergyLoss( G4double step ) 1218 { 1103 { 1219 G4int iTransfer ; << 1104 G4int iTransfer ; 1220 G4long numOfCollisions; << 1105 G4long numOfCollisions ; 1221 G4double loss = 0.0; << 1106 G4double loss = 0.0 ; 1222 G4double meanNumber, position; << 1107 G4double meanNumber, position ; 1223 1108 1224 // G4cout<<" G4PAIySection::GetStepEnergyLo << 1109 // G4cout<<" G4PAIySection::GetStepEnergyLoss "<<G4endl ; 1225 1110 1226 1111 1227 1112 1228 meanNumber = fIntegralPAIySection[1]*step; << 1113 meanNumber = fIntegralPAIySection[1]*step ; 1229 numOfCollisions = G4Poisson(meanNumber); << 1114 numOfCollisions = G4Poisson(meanNumber) ; 1230 1115 1231 // G4cout<<"numOfCollisions = "<<numOfCol << 1116 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl ; 1232 1117 1233 while(numOfCollisions) 1118 while(numOfCollisions) 1234 { 1119 { 1235 position = fIntegralPAIySection[1]*G4Unif << 1120 position = fIntegralPAIySection[1]*G4UniformRand() ; 1236 1121 1237 for( iTransfer=1; iTransfer<=fSplineNumbe << 1122 for( iTransfer=1 ; iTransfer<=fSplineNumber ; iTransfer++ ) 1238 { 1123 { 1239 if( position >= fIntegralPAIySection[ << 1124 if( position >= fIntegralPAIySection[iTransfer] ) break ; 1240 } 1125 } 1241 loss += fSplineEnergy[iTransfer] ; << 1126 loss += fSplineEnergy[iTransfer] ; 1242 numOfCollisions--; << 1127 numOfCollisions-- ; 1243 // Loop checking, 03-Aug-2015, Vladimir I << 1244 } 1128 } 1245 // G4cout<<"PAI energy loss = "<<loss/keV<< << 1129 // G4cout<<"PAI energy loss = "<<loss/keV<<" keV"<<G4endl ; 1246 1130 1247 return loss; << 1131 return loss ; 1248 } 1132 } 1249 1133 1250 ///////////////////////////////////////////// 1134 ///////////////////////////////////////////////////////////////////////// 1251 // 1135 // 1252 // 1136 // 1253 1137 1254 G4double G4PAIySection::GetStepCerenkovLoss( 1138 G4double G4PAIySection::GetStepCerenkovLoss( G4double step ) 1255 { 1139 { 1256 G4int iTransfer ; << 1140 G4int iTransfer ; 1257 G4long numOfCollisions; << 1141 G4long numOfCollisions ; 1258 G4double loss = 0.0; << 1142 G4double loss = 0.0 ; 1259 G4double meanNumber, position; << 1143 G4double meanNumber, position ; 1260 1144 1261 // G4cout<<" G4PAIySection::GetStepCreLosnk << 1145 // G4cout<<" G4PAIySection::GetStepCreLosnkovs "<<G4endl ; 1262 1146 1263 1147 1264 1148 1265 meanNumber = fIntegralCerenkov[1]*step; << 1149 meanNumber = fIntegralCerenkov[1]*step ; 1266 numOfCollisions = G4Poisson(meanNumber); << 1150 numOfCollisions = G4Poisson(meanNumber) ; 1267 1151 1268 // G4cout<<"numOfCollisions = "<<numOfCol << 1152 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl ; 1269 1153 1270 while(numOfCollisions) 1154 while(numOfCollisions) 1271 { 1155 { 1272 position = fIntegralCerenkov[1]*G4Uniform << 1156 position = fIntegralCerenkov[1]*G4UniformRand() ; 1273 1157 1274 for( iTransfer=1; iTransfer<=fSplineNumbe << 1158 for( iTransfer=1 ; iTransfer<=fSplineNumber ; iTransfer++ ) 1275 { 1159 { 1276 if( position >= fIntegralCerenkov[iTr << 1160 if( position >= fIntegralCerenkov[iTransfer] ) break ; 1277 } 1161 } 1278 loss += fSplineEnergy[iTransfer] ; << 1162 loss += fSplineEnergy[iTransfer] ; 1279 numOfCollisions--; << 1163 numOfCollisions-- ; 1280 // Loop checking, 03-Aug-2015, Vladimir I << 1281 } 1164 } 1282 // G4cout<<"PAI Cerenkov loss = "<<loss/keV << 1165 // G4cout<<"PAI Cerenkov loss = "<<loss/keV<<" keV"<<G4endl ; 1283 1166 1284 return loss; << 1167 return loss ; 1285 } 1168 } 1286 1169 1287 ///////////////////////////////////////////// 1170 ///////////////////////////////////////////////////////////////////////// 1288 // 1171 // 1289 // 1172 // 1290 1173 1291 G4double G4PAIySection::GetStepPlasmonLoss( G 1174 G4double G4PAIySection::GetStepPlasmonLoss( G4double step ) 1292 { 1175 { 1293 G4int iTransfer ; << 1176 G4int iTransfer ; 1294 G4long numOfCollisions; << 1177 G4long numOfCollisions ; 1295 G4double loss = 0.0; << 1178 G4double loss = 0.0 ; 1296 G4double meanNumber, position; << 1179 G4double meanNumber, position ; 1297 1180 1298 // G4cout<<" G4PAIySection::GetStepCreLosnk << 1181 // G4cout<<" G4PAIySection::GetStepCreLosnkovs "<<G4endl ; 1299 1182 1300 1183 1301 1184 1302 meanNumber = fIntegralPlasmon[1]*step; << 1185 meanNumber = fIntegralPlasmon[1]*step ; 1303 numOfCollisions = G4Poisson(meanNumber); << 1186 numOfCollisions = G4Poisson(meanNumber) ; 1304 1187 1305 // G4cout<<"numOfCollisions = "<<numOfCol << 1188 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl ; 1306 1189 1307 while(numOfCollisions) 1190 while(numOfCollisions) 1308 { 1191 { 1309 position = fIntegralPlasmon[1]*G4UniformR << 1192 position = fIntegralPlasmon[1]*G4UniformRand() ; 1310 1193 1311 for( iTransfer=1; iTransfer<=fSplineNumbe << 1194 for( iTransfer=1 ; iTransfer<=fSplineNumber ; iTransfer++ ) 1312 { 1195 { 1313 if( position >= fIntegralPlasmon[iTra << 1196 if( position >= fIntegralPlasmon[iTransfer] ) break ; 1314 } 1197 } 1315 loss += fSplineEnergy[iTransfer] ; << 1198 loss += fSplineEnergy[iTransfer] ; 1316 numOfCollisions--; << 1199 numOfCollisions-- ; 1317 // Loop checking, 03-Aug-2015, Vladimir I << 1318 } 1200 } 1319 // G4cout<<"PAI Plasmon loss = "<<loss/keV< << 1201 // G4cout<<"PAI Plasmon loss = "<<loss/keV<<" keV"<<G4endl ; 1320 1202 1321 return loss; << 1203 return loss ; 1322 } 1204 } 1323 1205 1324 ///////////////////////////////////////////// << 1325 // << 1326 1206 1327 void G4PAIySection::CallError(G4int i, const << 1328 { << 1329 G4String head = "G4PAIySection::" + methodN << 1330 G4ExceptionDescription ed; << 1331 ed << "Wrong index " << i << " fSplineNumbe << 1332 G4Exception(head,"pai001",FatalException,ed << 1333 } << 1334 1207 1335 ///////////////////////////////////////////// 1208 ///////////////////////////////////////////////////////////////////////////// 1336 // 1209 // 1337 // Init array of Lorentz factors 1210 // Init array of Lorentz factors 1338 // 1211 // 1339 1212 1340 G4int G4PAIySection::fNumberOfGammas = 111; << 1213 G4int G4PAIySection::fNumberOfGammas = 111 ; 1341 1214 1342 const G4double G4PAIySection::fLorentzFactor[ 1215 const G4double G4PAIySection::fLorentzFactor[112] = // fNumberOfGammas+1 1343 { 1216 { 1344 0.0, 1217 0.0, 1345 1.094989e+00, 1.107813e+00, 1.122369e+00, 1.1 1218 1.094989e+00, 1.107813e+00, 1.122369e+00, 1.138890e+00, 1.157642e+00, 1346 1.178925e+00, 1.203082e+00, 1.230500e+00, 1.2 1219 1.178925e+00, 1.203082e+00, 1.230500e+00, 1.261620e+00, 1.296942e+00, // 10 1347 1.337032e+00, 1.382535e+00, 1.434181e+00, 1.4 1220 1.337032e+00, 1.382535e+00, 1.434181e+00, 1.492800e+00, 1.559334e+00, 1348 1.634850e+00, 1.720562e+00, 1.817845e+00, 1.9 1221 1.634850e+00, 1.720562e+00, 1.817845e+00, 1.928263e+00, 2.053589e+00, // 20 1349 2.195835e+00, 2.357285e+00, 2.540533e+00, 2.7 1222 2.195835e+00, 2.357285e+00, 2.540533e+00, 2.748522e+00, 2.984591e+00, 1350 3.252533e+00, 3.556649e+00, 3.901824e+00, 4.2 1223 3.252533e+00, 3.556649e+00, 3.901824e+00, 4.293602e+00, 4.738274e+00, // 30 1351 5.242981e+00, 5.815829e+00, 6.466019e+00, 7.2 1224 5.242981e+00, 5.815829e+00, 6.466019e+00, 7.203990e+00, 8.041596e+00, 1352 8.992288e+00, 1.007133e+01, 1.129606e+01, 1.2 1225 8.992288e+00, 1.007133e+01, 1.129606e+01, 1.268614e+01, 1.426390e+01, // 40 1353 1.605467e+01, 1.808721e+01, 2.039417e+01, 2.3 1226 1.605467e+01, 1.808721e+01, 2.039417e+01, 2.301259e+01, 2.598453e+01, 1354 2.935771e+01, 3.318630e+01, 3.753180e+01, 4.2 1227 2.935771e+01, 3.318630e+01, 3.753180e+01, 4.246399e+01, 4.806208e+01, // 50 1355 5.441597e+01, 6.162770e+01, 6.981310e+01, 7.9 1228 5.441597e+01, 6.162770e+01, 6.981310e+01, 7.910361e+01, 8.964844e+01, 1356 1.016169e+02, 1.152013e+02, 1.306197e+02, 1.4 1229 1.016169e+02, 1.152013e+02, 1.306197e+02, 1.481198e+02, 1.679826e+02, // 60 1357 1.905270e+02, 2.161152e+02, 2.451581e+02, 2.7 1230 1.905270e+02, 2.161152e+02, 2.451581e+02, 2.781221e+02, 3.155365e+02, 1358 3.580024e+02, 4.062016e+02, 4.609081e+02, 5.2 1231 3.580024e+02, 4.062016e+02, 4.609081e+02, 5.230007e+02, 5.934765e+02, // 70 1359 6.734672e+02, 7.642575e+02, 8.673056e+02, 9.8 1232 6.734672e+02, 7.642575e+02, 8.673056e+02, 9.842662e+02, 1.117018e+03, 1360 1.267692e+03, 1.438709e+03, 1.632816e+03, 1.8 1233 1.267692e+03, 1.438709e+03, 1.632816e+03, 1.853128e+03, 2.103186e+03, // 80 1361 2.387004e+03, 2.709140e+03, 3.074768e+03, 3.4 1234 2.387004e+03, 2.709140e+03, 3.074768e+03, 3.489760e+03, 3.960780e+03, 1362 4.495394e+03, 5.102185e+03, 5.790900e+03, 6.5 1235 4.495394e+03, 5.102185e+03, 5.790900e+03, 6.572600e+03, 7.459837e+03, // 90 1363 8.466860e+03, 9.609843e+03, 1.090714e+04, 1.2 1236 8.466860e+03, 9.609843e+03, 1.090714e+04, 1.237959e+04, 1.405083e+04, 1364 1.594771e+04, 1.810069e+04, 2.054434e+04, 2.3 1237 1.594771e+04, 1.810069e+04, 2.054434e+04, 2.331792e+04, 2.646595e+04, // 100 1365 3.003901e+04, 3.409446e+04, 3.869745e+04, 4.3 1238 3.003901e+04, 3.409446e+04, 3.869745e+04, 4.392189e+04, 4.985168e+04, 1366 5.658206e+04, 6.422112e+04, 7.289153e+04, 8.2 1239 5.658206e+04, 6.422112e+04, 7.289153e+04, 8.273254e+04, 9.390219e+04, // 110 1367 1.065799e+05 1240 1.065799e+05 1368 }; << 1241 } ; 1369 1242 1370 ///////////////////////////////////////////// 1243 /////////////////////////////////////////////////////////////////////// 1371 // 1244 // 1372 // The number of gamma for creation of splin 1245 // The number of gamma for creation of spline (near ion-min , G ~ 4 ) 1373 // 1246 // 1374 1247 1375 const G4int G4PAIySection::fRefGammaNumber = << 1248 const >> 1249 G4int G4PAIySection::fRefGammaNumber = 29 ; 1376 1250 >> 1251 1377 // 1252 // 1378 // end of G4PAIySection implementation file 1253 // end of G4PAIySection implementation file 1379 // 1254 // 1380 ///////////////////////////////////////////// 1255 //////////////////////////////////////////////////////////////////////////// 1381 1256 1382 1257