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