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: G4PAIxSection.cc 83662 2014-09-08 10:00:20Z gcosmo $ >> 28 // GEANT4 tag $Name: geant4-09-03-ref-06 $ 27 // 29 // 28 // 30 // 29 // G4PAIxSection.cc -- class implementation fi 31 // G4PAIxSection.cc -- class implementation file 30 // 32 // 31 // GEANT 4 class implementation file 33 // GEANT 4 class implementation file 32 // 34 // 33 // For information related to this code, pleas 35 // For information related to this code, please, contact 34 // the Geant4 Collaboration. 36 // the Geant4 Collaboration. 35 // 37 // 36 // R&D: Vladimir.Grichine@cern.ch 38 // R&D: Vladimir.Grichine@cern.ch 37 // 39 // 38 // History: 40 // History: 39 // 41 // 40 // 13.05.03 V. Grichine, bug fixed for maxEner 42 // 13.05.03 V. Grichine, bug fixed for maxEnergyTransfer > max interval energy 41 // 28.05.01 V.Ivanchenko minor changes to prov 43 // 28.05.01 V.Ivanchenko minor changes to provide ANSI -wall compilation 42 // 17.05.01 V. Grichine, low energy extension 44 // 17.05.01 V. Grichine, low energy extension down to 10*keV of proton 43 // 20.11.98 adapted to a new Material/SandiaTa 45 // 20.11.98 adapted to a new Material/SandiaTable interface, mma 44 // 11.06.97 V. Grichine, 1st version 46 // 11.06.97 V. Grichine, 1st version 45 // 47 // 46 48 47 #include "G4PAIxSection.hh" 49 #include "G4PAIxSection.hh" 48 50 49 #include "globals.hh" 51 #include "globals.hh" 50 #include "G4PhysicalConstants.hh" 52 #include "G4PhysicalConstants.hh" 51 #include "G4SystemOfUnits.hh" 53 #include "G4SystemOfUnits.hh" 52 #include "G4ios.hh" 54 #include "G4ios.hh" 53 #include "G4Poisson.hh" 55 #include "G4Poisson.hh" 54 #include "G4Material.hh" 56 #include "G4Material.hh" 55 #include "G4MaterialCutsCouple.hh" 57 #include "G4MaterialCutsCouple.hh" 56 #include "G4SandiaTable.hh" 58 #include "G4SandiaTable.hh" 57 59 58 using namespace std; 60 using namespace std; 59 61 60 /* ******************************************* 62 /* ****************************************************************** 61 63 62 // Init array of Lorentz factors 64 // Init array of Lorentz factors 63 65 64 const G4double G4PAIxSection::fLorentzFactor[2 66 const G4double G4PAIxSection::fLorentzFactor[22] = 65 { 67 { 66 0.0 , 1.1 , 1.2 , 1.3 , 1 68 0.0 , 1.1 , 1.2 , 1.3 , 1.5 , 1.8 , 2.0 , 67 2.5 , 3.0 , 4.0 , 7.0 , 10 69 2.5 , 3.0 , 4.0 , 7.0 , 10.0 , 20.0 , 40.0 , 68 70.0 , 100.0 , 300.0 , 600.0 , 1000 70 70.0 , 100.0 , 300.0 , 600.0 , 1000.0 , 3000.0 , 69 10000.0 , 50000.0 71 10000.0 , 50000.0 70 }; 72 }; 71 73 72 const G4int G4PAIxSection:: 74 const G4int G4PAIxSection:: 73 fRefGammaNumber = 29; // The number of 75 fRefGammaNumber = 29; // The number of gamma for creation of 74 // spline (9) 76 // spline (9) 75 77 76 ********************************************** 78 ***************************************************************** */ 77 79 78 // Local class constants 80 // Local class constants 79 81 80 const G4double G4PAIxSection::fDelta = 0.005; 82 const G4double G4PAIxSection::fDelta = 0.005; // 0.005 energy shift from interval border 81 const G4double G4PAIxSection::fError = 0.005; 83 const G4double G4PAIxSection::fError = 0.005; // 0.005 error in lin-log approximation 82 84 83 const G4int G4PAIxSection::fMaxSplineSize = 10 85 const G4int G4PAIxSection::fMaxSplineSize = 1000; // Max size of output spline 84 86 // arrays 85 ////////////////////////////////////////////// 87 ////////////////////////////////////////////////////////////////// 86 // 88 // 87 // Constructor 89 // Constructor 88 // 90 // 89 91 90 G4PAIxSection::G4PAIxSection() 92 G4PAIxSection::G4PAIxSection() 91 { 93 { 92 fSandia = nullptr; << 94 fSandia = 0; 93 fMatSandiaMatrix = nullptr; << 95 fMatSandiaMatrix = 0; 94 fDensity = fElectronDensity = fNormalization 96 fDensity = fElectronDensity = fNormalizationCof = fLowEnergyCof = 0.0; 95 fIntervalNumber = fSplineNumber = 0; 97 fIntervalNumber = fSplineNumber = 0; 96 fVerbose = 0; 98 fVerbose = 0; 97 99 98 fSplineEnergy = G4DataVector(fMaxSp 100 fSplineEnergy = G4DataVector(fMaxSplineSize,0.0); 99 fRePartDielectricConst = G4DataVector(fMaxSp 101 fRePartDielectricConst = G4DataVector(fMaxSplineSize,0.0); 100 fImPartDielectricConst = G4DataVector(fMaxSp 102 fImPartDielectricConst = G4DataVector(fMaxSplineSize,0.0); 101 fIntegralTerm = G4DataVector(fMaxSp 103 fIntegralTerm = G4DataVector(fMaxSplineSize,0.0); 102 fDifPAIxSection = G4DataVector(fMaxSp 104 fDifPAIxSection = G4DataVector(fMaxSplineSize,0.0); 103 fdNdxCerenkov = G4DataVector(fMaxSp 105 fdNdxCerenkov = G4DataVector(fMaxSplineSize,0.0); 104 fdNdxPlasmon = G4DataVector(fMaxSp 106 fdNdxPlasmon = G4DataVector(fMaxSplineSize,0.0); 105 fdNdxMM = G4DataVector(fMaxSp 107 fdNdxMM = G4DataVector(fMaxSplineSize,0.0); 106 fdNdxResonance = G4DataVector(fMaxSp 108 fdNdxResonance = G4DataVector(fMaxSplineSize,0.0); 107 fIntegralPAIxSection = G4DataVector(fMaxSp 109 fIntegralPAIxSection = G4DataVector(fMaxSplineSize,0.0); 108 fIntegralPAIdEdx = G4DataVector(fMaxSp 110 fIntegralPAIdEdx = G4DataVector(fMaxSplineSize,0.0); 109 fIntegralCerenkov = G4DataVector(fMaxSp 111 fIntegralCerenkov = G4DataVector(fMaxSplineSize,0.0); 110 fIntegralPlasmon = G4DataVector(fMaxSp 112 fIntegralPlasmon = G4DataVector(fMaxSplineSize,0.0); 111 fIntegralMM = G4DataVector(fMaxSp 113 fIntegralMM = G4DataVector(fMaxSplineSize,0.0); 112 fIntegralResonance = G4DataVector(fMaxSp 114 fIntegralResonance = G4DataVector(fMaxSplineSize,0.0); 113 115 114 fMaterialIndex = 0; 116 fMaterialIndex = 0; 115 117 116 for( G4int i = 0; i < 500; ++i ) 118 for( G4int i = 0; i < 500; ++i ) 117 { 119 { 118 for( G4int j = 0; j < 112; ++j ) fPAItabl 120 for( G4int j = 0; j < 112; ++j ) fPAItable[i][j] = 0.0; 119 } 121 } 120 } 122 } 121 123 122 ////////////////////////////////////////////// 124 ////////////////////////////////////////////////////////////////// 123 // 125 // 124 // Constructor 126 // Constructor 125 // 127 // 126 128 127 G4PAIxSection::G4PAIxSection(G4MaterialCutsCou 129 G4PAIxSection::G4PAIxSection(G4MaterialCutsCouple* matCC) 128 { 130 { 129 fDensity = matCC->GetMaterial()->GetDe 131 fDensity = matCC->GetMaterial()->GetDensity(); 130 G4int matIndex = (G4int)matCC->GetMaterial() << 132 G4int matIndex = matCC->GetMaterial()->GetIndex(); 131 fMaterialIndex = matIndex; 133 fMaterialIndex = matIndex; 132 134 133 const G4MaterialTable* theMaterialTable = G4 135 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable(); 134 fSandia = (*theMaterialTable)[matIndex]->Get 136 fSandia = (*theMaterialTable)[matIndex]->GetSandiaTable(); 135 137 136 fVerbose = 0; 138 fVerbose = 0; 137 139 138 G4int i, j; 140 G4int i, j; 139 fMatSandiaMatrix = new G4OrderedTable(); 141 fMatSandiaMatrix = new G4OrderedTable(); 140 142 141 for (i = 0; i < fSandia->GetMaxInterval()-1; << 143 for (i = 0; i < fSandia->GetMaxInterval()-1; i++) 142 { 144 { 143 fMatSandiaMatrix->push_back(new G4DataVec 145 fMatSandiaMatrix->push_back(new G4DataVector(5,0.)); 144 } << 146 } 145 for (i = 0; i < fSandia->GetMaxInterval()-1; << 147 for (i = 0; i < fSandia->GetMaxInterval()-1; i++) 146 { 148 { 147 (*(*fMatSandiaMatrix)[i])[0] = fSandia->Ge 149 (*(*fMatSandiaMatrix)[i])[0] = fSandia->GetSandiaMatTable(i,0); 148 150 149 for(j = 1; j < 5; ++j) << 151 for(j = 1; j < 5; j++) 150 { 152 { 151 (*(*fMatSandiaMatrix)[i])[j] = fSandia-> 153 (*(*fMatSandiaMatrix)[i])[j] = fSandia->GetSandiaMatTable(i,j)*fDensity; 152 } 154 } 153 } 155 } 154 ComputeLowEnergyCof(); << 156 ComputeLowEnergyCof(); 155 // fEnergyInterval = fA1 = fA2 = fA3 = fA4 157 // fEnergyInterval = fA1 = fA2 = fA3 = fA4 = 0; 156 } 158 } 157 159 158 ////////////////////////////////////////////// 160 //////////////////////////////////////////////////////////////// 159 161 160 G4PAIxSection::G4PAIxSection(G4int materialInd 162 G4PAIxSection::G4PAIxSection(G4int materialIndex, 161 G4double maxEnerg << 163 G4double maxEnergyTransfer) 162 { 164 { 163 fSandia = nullptr; << 165 fSandia = 0; 164 fMatSandiaMatrix = nullptr; << 166 fMatSandiaMatrix = 0; 165 fVerbose = 0; 167 fVerbose = 0; 166 const G4MaterialTable* theMaterialTable = G4 168 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable(); 167 G4int i, j; 169 G4int i, j; 168 170 169 fMaterialIndex = materialIndex; 171 fMaterialIndex = materialIndex; 170 fDensity = (*theMaterialTable 172 fDensity = (*theMaterialTable)[materialIndex]->GetDensity(); 171 fElectronDensity = (*theMaterialTable 173 fElectronDensity = (*theMaterialTable)[materialIndex]-> 172 GetElectronDensit 174 GetElectronDensity(); 173 fIntervalNumber = (*theMaterialTable 175 fIntervalNumber = (*theMaterialTable)[materialIndex]-> 174 GetSandiaTable()- 176 GetSandiaTable()->GetMatNbOfIntervals(); 175 fIntervalNumber--; 177 fIntervalNumber--; 176 // G4cout<<fDensity<<"\t"<<fElectronDensity< 178 // G4cout<<fDensity<<"\t"<<fElectronDensity<<"\t"<<fIntervalNumber<<G4endl; 177 179 178 fEnergyInterval = G4DataVector(fIntervalNumb 180 fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0); 179 fA1 = G4DataVector(fIntervalNumb 181 fA1 = G4DataVector(fIntervalNumber+2,0.0); 180 fA2 = G4DataVector(fIntervalNumb 182 fA2 = G4DataVector(fIntervalNumber+2,0.0); 181 fA3 = G4DataVector(fIntervalNumb 183 fA3 = G4DataVector(fIntervalNumber+2,0.0); 182 fA4 = G4DataVector(fIntervalNumb 184 fA4 = G4DataVector(fIntervalNumber+2,0.0); 183 185 184 for(i = 1; i <= fIntervalNumber; i++ ) 186 for(i = 1; i <= fIntervalNumber; i++ ) 185 { 187 { 186 if(((*theMaterialTable)[materialIndex]-> 188 if(((*theMaterialTable)[materialIndex]-> 187 GetSandiaTable()->GetSandiaCofForMaterial( 189 GetSandiaTable()->GetSandiaCofForMaterial(i-1,0) >= maxEnergyTransfer) || 188 i > fIntervalNumber 190 i > fIntervalNumber ) 189 { << 191 { 190 fEnergyInterval[i] = maxEnergyTransf << 192 fEnergyInterval[i] = maxEnergyTransfer; 191 fIntervalNumber = i; << 193 fIntervalNumber = i; 192 break; << 194 break; 193 } << 195 } 194 fEnergyInterval[i] = (*theMaterialTab 196 fEnergyInterval[i] = (*theMaterialTable)[materialIndex]-> 195 GetSandiaTable() << 197 GetSandiaTable()->GetSandiaCofForMaterial(i-1,0); 196 fA1[i] = (*theMaterialTab 198 fA1[i] = (*theMaterialTable)[materialIndex]-> 197 GetSandiaTable() << 199 GetSandiaTable()->GetSandiaCofForMaterial(i-1,1); 198 fA2[i] = (*theMaterialTab 200 fA2[i] = (*theMaterialTable)[materialIndex]-> 199 GetSandiaTable() << 201 GetSandiaTable()->GetSandiaCofForMaterial(i-1,2); 200 fA3[i] = (*theMaterialTab 202 fA3[i] = (*theMaterialTable)[materialIndex]-> 201 GetSandiaTable() << 203 GetSandiaTable()->GetSandiaCofForMaterial(i-1,3); 202 fA4[i] = (*theMaterialTab 204 fA4[i] = (*theMaterialTable)[materialIndex]-> 203 GetSandiaTable() << 205 GetSandiaTable()->GetSandiaCofForMaterial(i-1,4); 204 // G4cout<<i<<"\t"<<fEnergyInterval[i << 206 // G4cout<<i<<"\t"<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t" 205 // <<fA << 207 // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl; 206 } 208 } 207 if(fEnergyInterval[fIntervalNumber] != maxEn 209 if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer) 208 { 210 { 209 fIntervalNumber++; 211 fIntervalNumber++; 210 fEnergyInterval[fIntervalNumber] = ma 212 fEnergyInterval[fIntervalNumber] = maxEnergyTransfer; 211 } 213 } 212 214 213 // Now checking, if two borders are too clos 215 // Now checking, if two borders are too close together 214 216 215 for(i=1;i<fIntervalNumber;i++) 217 for(i=1;i<fIntervalNumber;i++) 216 { 218 { 217 if(fEnergyInterval[i+1]-fEnergyInterva 219 if(fEnergyInterval[i+1]-fEnergyInterval[i] > 218 1.5*fDelta*(fEnergyInterval[i+1]+fE 220 1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i])) 219 { << 221 { 220 continue; 222 continue; 221 } << 223 } 222 else 224 else 223 { << 225 { 224 for(j=i;j<fIntervalNumber;j++) 226 for(j=i;j<fIntervalNumber;j++) 225 { << 227 { 226 fEnergyInterval[j] = fEnergyInterv 228 fEnergyInterval[j] = fEnergyInterval[j+1]; 227 fA1[j] = fA1[j+1]; 229 fA1[j] = fA1[j+1]; 228 fA2[j] = fA2[j+1]; 230 fA2[j] = fA2[j+1]; 229 fA3[j] = fA3[j+1]; 231 fA3[j] = fA3[j+1]; 230 fA4[j] = fA4[j+1]; 232 fA4[j] = fA4[j+1]; 231 } << 233 } 232 fIntervalNumber--; 234 fIntervalNumber--; 233 i--; 235 i--; 234 } << 236 } 235 } 237 } 236 238 237 239 238 /* ********************************* 240 /* ********************************* 239 241 240 fSplineEnergy = new G4double[fM 242 fSplineEnergy = new G4double[fMaxSplineSize]; 241 fRePartDielectricConst = new G4double[fM 243 fRePartDielectricConst = new G4double[fMaxSplineSize]; 242 fImPartDielectricConst = new G4double[fM 244 fImPartDielectricConst = new G4double[fMaxSplineSize]; 243 fIntegralTerm = new G4double[fM 245 fIntegralTerm = new G4double[fMaxSplineSize]; 244 fDifPAIxSection = new G4double[fM 246 fDifPAIxSection = new G4double[fMaxSplineSize]; 245 fIntegralPAIxSection = new G4double[fM 247 fIntegralPAIxSection = new G4double[fMaxSplineSize]; 246 248 247 for(i=0;i<fMaxSplineSize;i++) 249 for(i=0;i<fMaxSplineSize;i++) 248 { 250 { 249 fSplineEnergy[i] = 0.0; 251 fSplineEnergy[i] = 0.0; 250 fRePartDielectricConst[i] = 0.0; 252 fRePartDielectricConst[i] = 0.0; 251 fImPartDielectricConst[i] = 0.0; 253 fImPartDielectricConst[i] = 0.0; 252 fIntegralTerm[i] = 0.0; 254 fIntegralTerm[i] = 0.0; 253 fDifPAIxSection[i] = 0.0; 255 fDifPAIxSection[i] = 0.0; 254 fIntegralPAIxSection[i] = 0.0; 256 fIntegralPAIxSection[i] = 0.0; 255 } 257 } 256 **************************************** 258 ************************************************** */ 257 ComputeLowEnergyCof(); 259 ComputeLowEnergyCof(); 258 InitPAI(); // create arrays allocated a 260 InitPAI(); // create arrays allocated above 259 /* 261 /* 260 delete[] fEnergyInterval; 262 delete[] fEnergyInterval; 261 delete[] fA1; 263 delete[] fA1; 262 delete[] fA2; 264 delete[] fA2; 263 delete[] fA3; 265 delete[] fA3; 264 delete[] fA4; 266 delete[] fA4; 265 */ 267 */ 266 } 268 } 267 269 268 ////////////////////////////////////////////// 270 //////////////////////////////////////////////////////////////////////// 269 // 271 // 270 // Constructor called from G4PAIPhotonModel !! 272 // Constructor called from G4PAIPhotonModel !!! 271 273 272 G4PAIxSection::G4PAIxSection( G4int materialIn 274 G4PAIxSection::G4PAIxSection( G4int materialIndex, 273 G4double maxEner << 275 G4double maxEnergyTransfer, 274 G4double betaGam << 276 G4double betaGammaSq, 275 G4double** photo 277 G4double** photoAbsCof, 276 G4int intNumber 278 G4int intNumber ) 277 { 279 { 278 fSandia = nullptr; << 280 fSandia = 0; 279 fDensity = fElectronDensity = fNormalization 281 fDensity = fElectronDensity = fNormalizationCof = fLowEnergyCof = 0.0; 280 fIntervalNumber = fSplineNumber = 0; 282 fIntervalNumber = fSplineNumber = 0; 281 fVerbose = 0; 283 fVerbose = 0; 282 284 283 fSplineEnergy = G4DataVector(500,0. 285 fSplineEnergy = G4DataVector(500,0.0); 284 fRePartDielectricConst = G4DataVector(500,0. 286 fRePartDielectricConst = G4DataVector(500,0.0); 285 fImPartDielectricConst = G4DataVector(500,0. 287 fImPartDielectricConst = G4DataVector(500,0.0); 286 fIntegralTerm = G4DataVector(500,0. 288 fIntegralTerm = G4DataVector(500,0.0); 287 fDifPAIxSection = G4DataVector(500,0. 289 fDifPAIxSection = G4DataVector(500,0.0); 288 fdNdxCerenkov = G4DataVector(500,0. 290 fdNdxCerenkov = G4DataVector(500,0.0); 289 fdNdxPlasmon = G4DataVector(500,0. 291 fdNdxPlasmon = G4DataVector(500,0.0); 290 fdNdxMM = G4DataVector(500,0. 292 fdNdxMM = G4DataVector(500,0.0); 291 fdNdxResonance = G4DataVector(500,0. 293 fdNdxResonance = G4DataVector(500,0.0); 292 fIntegralPAIxSection = G4DataVector(500,0. 294 fIntegralPAIxSection = G4DataVector(500,0.0); 293 fIntegralPAIdEdx = G4DataVector(500,0. 295 fIntegralPAIdEdx = G4DataVector(500,0.0); 294 fIntegralCerenkov = G4DataVector(500,0. 296 fIntegralCerenkov = G4DataVector(500,0.0); 295 fIntegralPlasmon = G4DataVector(500,0. 297 fIntegralPlasmon = G4DataVector(500,0.0); 296 fIntegralMM = G4DataVector(500,0. 298 fIntegralMM = G4DataVector(500,0.0); 297 fIntegralResonance = G4DataVector(500,0. 299 fIntegralResonance = G4DataVector(500,0.0); 298 300 299 for( G4int i = 0; i < 500; ++i ) 301 for( G4int i = 0; i < 500; ++i ) 300 { 302 { 301 for( G4int j = 0; j < 112; ++j ) fPAItabl 303 for( G4int j = 0; j < 112; ++j ) fPAItable[i][j] = 0.0; 302 } 304 } 303 305 304 fSandia = nullptr; << 306 fSandia = 0; 305 fMatSandiaMatrix = nullptr; << 307 fMatSandiaMatrix = 0; 306 const G4MaterialTable* theMaterialTable = G4 308 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable(); 307 G4int i, j; 309 G4int i, j; 308 310 309 fMaterialIndex = materialIndex; 311 fMaterialIndex = materialIndex; 310 fDensity = (*theMaterialTable)[mater 312 fDensity = (*theMaterialTable)[materialIndex]->GetDensity(); 311 fElectronDensity = (*theMaterialTable)[mater 313 fElectronDensity = (*theMaterialTable)[materialIndex]->GetElectronDensity(); 312 314 313 fIntervalNumber = intNumber; 315 fIntervalNumber = intNumber; 314 fIntervalNumber--; 316 fIntervalNumber--; 315 // G4cout<<fDensity<<"\t"<<fElectronDensit 317 // G4cout<<fDensity<<"\t"<<fElectronDensity<<"\t"<<fIntervalNumber<<G4endl; 316 318 317 fEnergyInterval = G4DataVector(fIntervalNumb 319 fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0); 318 fA1 = G4DataVector(fIntervalNumb 320 fA1 = G4DataVector(fIntervalNumber+2,0.0); 319 fA2 = G4DataVector(fIntervalNumb 321 fA2 = G4DataVector(fIntervalNumber+2,0.0); 320 fA3 = G4DataVector(fIntervalNumb 322 fA3 = G4DataVector(fIntervalNumber+2,0.0); 321 fA4 = G4DataVector(fIntervalNumb 323 fA4 = G4DataVector(fIntervalNumber+2,0.0); 322 324 323 325 324 /* 326 /* 325 fEnergyInterval = new G4double[fInterval 327 fEnergyInterval = new G4double[fIntervalNumber+2]; 326 fA1 = new G4double[fInterval 328 fA1 = new G4double[fIntervalNumber+2]; 327 fA2 = new G4double[fInterval 329 fA2 = new G4double[fIntervalNumber+2]; 328 fA3 = new G4double[fInterval 330 fA3 = new G4double[fIntervalNumber+2]; 329 fA4 = new G4double[fInterval 331 fA4 = new G4double[fIntervalNumber+2]; 330 */ 332 */ 331 for( i = 1; i <= fIntervalNumber; i++ ) 333 for( i = 1; i <= fIntervalNumber; i++ ) 332 { 334 { 333 if( ( photoAbsCof[i-1][0] >= maxEnerg 335 if( ( photoAbsCof[i-1][0] >= maxEnergyTransfer ) || 334 i > fIntervalNumber ) 336 i > fIntervalNumber ) 335 { 337 { 336 fEnergyInterval[i] = maxEnergyTran 338 fEnergyInterval[i] = maxEnergyTransfer; 337 fIntervalNumber = i; << 339 fIntervalNumber = i; 338 break; << 340 break; 339 } 341 } 340 fEnergyInterval[i] = photoAbsCof[i-1] 342 fEnergyInterval[i] = photoAbsCof[i-1][0]; 341 fA1[i] = photoAbsCof[i-1] 343 fA1[i] = photoAbsCof[i-1][1]; 342 fA2[i] = photoAbsCof[i-1] 344 fA2[i] = photoAbsCof[i-1][2]; 343 fA3[i] = photoAbsCof[i-1] 345 fA3[i] = photoAbsCof[i-1][3]; 344 fA4[i] = photoAbsCof[i-1] 346 fA4[i] = photoAbsCof[i-1][4]; 345 // G4cout<<i<<"\t"<<fEnergyInterval[i << 347 // G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t" 346 // <<fA3[i]<<"\t"<<fA4[i]<<"\t"< << 348 // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl; 347 } 349 } 348 // G4cout<<"i last = "<<i<<"; "<<"fInter 350 // G4cout<<"i last = "<<i<<"; "<<"fIntervalNumber = "<<fIntervalNumber<<G4endl; 349 351 350 if(fEnergyInterval[fIntervalNumber] != maxEn 352 if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer) 351 { 353 { 352 fIntervalNumber++; 354 fIntervalNumber++; 353 fEnergyInterval[fIntervalNumber] = ma 355 fEnergyInterval[fIntervalNumber] = maxEnergyTransfer; 354 } 356 } 355 // G4cout<<"after check of max energy tr 357 // G4cout<<"after check of max energy transfer"<<G4endl; 356 358 357 for( i = 1; i <= fIntervalNumber; i++ ) 359 for( i = 1; i <= fIntervalNumber; i++ ) 358 { 360 { 359 // G4cout<<i<<"\t"<<fEnergyInterval[i] << 361 // G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t" 360 // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4e << 362 // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl; 361 } 363 } 362 // Now checking, if two borders are too 364 // Now checking, if two borders are too close together 363 365 364 for( i = 1; i < fIntervalNumber; i++ ) 366 for( i = 1; i < fIntervalNumber; i++ ) 365 { 367 { 366 if(fEnergyInterval[i+1]-fEnergyInterva 368 if(fEnergyInterval[i+1]-fEnergyInterval[i] > 367 1.5*fDelta*(fEnergyInterval[i+1]+fE 369 1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i])) 368 { << 370 { 369 continue; 371 continue; 370 } << 372 } 371 else 373 else 372 { << 374 { 373 for(j=i;j<fIntervalNumber;j++) 375 for(j=i;j<fIntervalNumber;j++) 374 { << 376 { 375 fEnergyInterval[j] = fEnergyInterv 377 fEnergyInterval[j] = fEnergyInterval[j+1]; 376 fA1[j] = fA1[j+1]; 378 fA1[j] = fA1[j+1]; 377 fA2[j] = fA2[j+1]; 379 fA2[j] = fA2[j+1]; 378 fA3[j] = fA3[j+1]; 380 fA3[j] = fA3[j+1]; 379 fA4[j] = fA4[j+1]; 381 fA4[j] = fA4[j+1]; 380 } << 382 } 381 fIntervalNumber--; 383 fIntervalNumber--; 382 i--; 384 i--; 383 } << 385 } 384 } 386 } 385 // G4cout<<"after check of close borders"<<G 387 // G4cout<<"after check of close borders"<<G4endl; 386 388 387 for( i = 1; i <= fIntervalNumber; i++ ) 389 for( i = 1; i <= fIntervalNumber; i++ ) 388 { 390 { 389 // G4cout<<i<<"\t"<<fEnergyInterval[i] << 391 // G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t" 390 // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4en << 392 // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl; 391 } 393 } 392 394 393 // Preparation of fSplineEnergy array corres 395 // Preparation of fSplineEnergy array corresponding to min ionisation, G~4 394 396 395 ComputeLowEnergyCof(); 397 ComputeLowEnergyCof(); 396 G4double betaGammaSqRef = 398 G4double betaGammaSqRef = 397 fLorentzFactor[fRefGammaNumber]*fLorentzFa 399 fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1; 398 400 399 NormShift(betaGammaSqRef); 401 NormShift(betaGammaSqRef); 400 SplainPAI(betaGammaSqRef); 402 SplainPAI(betaGammaSqRef); 401 403 402 // Preparation of integral PAI cross section 404 // Preparation of integral PAI cross section for input betaGammaSq 403 405 404 for(i = 1; i <= fSplineNumber; i++) 406 for(i = 1; i <= fSplineNumber; i++) 405 { 407 { 406 fdNdxCerenkov[i] = PAIdNdxCerenkov( 408 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq); 407 fdNdxMM[i] = PAIdNdxMM(i,betaGammaS 409 fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq); 408 fdNdxPlasmon[i] = PAIdNdxPlasmon(i 410 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq); 409 fdNdxResonance[i] = PAIdNdxResonance 411 fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq); 410 fDifPAIxSection[i] = DifPAIxSection(i 412 fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq); 411 413 412 // G4cout<<i<<"; dNdxC = "<<fdNdxCere << 414 // G4cout<<i<<"; dNdxC = "<<fdNdxCerenkov[i]<<"; dNdxP = "<<fdNdxPlasmon[i] 413 // <<"; dNdxPAI = "<<fDifPAIxSecti << 415 // <<"; dNdxPAI = "<<fDifPAIxSection[i]<<G4endl; 414 } 416 } 415 IntegralCerenkov(); 417 IntegralCerenkov(); 416 IntegralMM(); 418 IntegralMM(); 417 IntegralPlasmon(); 419 IntegralPlasmon(); 418 IntegralResonance(); 420 IntegralResonance(); 419 IntegralPAIxSection(); 421 IntegralPAIxSection(); 420 /* 422 /* 421 delete[] fEnergyInterval; 423 delete[] fEnergyInterval; 422 delete[] fA1; 424 delete[] fA1; 423 delete[] fA2; 425 delete[] fA2; 424 delete[] fA3; 426 delete[] fA3; 425 delete[] fA4; 427 delete[] fA4; 426 */ 428 */ 427 } 429 } 428 430 429 ////////////////////////////////////////////// 431 //////////////////////////////////////////////////////////////////////// 430 // 432 // 431 // Test Constructor with beta*gamma square val 433 // Test Constructor with beta*gamma square value 432 434 433 G4PAIxSection::G4PAIxSection( G4int materialIn 435 G4PAIxSection::G4PAIxSection( G4int materialIndex, 434 G4double maxEner << 436 G4double maxEnergyTransfer, 435 G4double betaGam << 437 G4double betaGammaSq ) 436 { 438 { 437 fSandia = nullptr; << 439 fSandia = 0; 438 fMatSandiaMatrix = nullptr; << 440 fMatSandiaMatrix = 0; 439 fVerbose = 0; 441 fVerbose = 0; 440 const G4MaterialTable* theMaterialTable = G4 442 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable(); 441 443 442 G4int i, j, numberOfElements; 444 G4int i, j, numberOfElements; 443 445 444 fMaterialIndex = materialIndex; 446 fMaterialIndex = materialIndex; 445 fDensity = (*theMaterialTable)[mater 447 fDensity = (*theMaterialTable)[materialIndex]->GetDensity(); 446 fElectronDensity = (*theMaterialTable)[mater 448 fElectronDensity = (*theMaterialTable)[materialIndex]->GetElectronDensity(); 447 numberOfElements = (G4int)(*theMaterialTable << 449 numberOfElements = (*theMaterialTable)[materialIndex]->GetNumberOfElements(); 448 450 449 G4int* thisMaterialZ = new G4int[numberOfEle 451 G4int* thisMaterialZ = new G4int[numberOfElements]; 450 452 451 for( i = 0; i < numberOfElements; ++i ) << 453 for( i = 0; i < numberOfElements; i++ ) 452 { 454 { 453 thisMaterialZ[i] = (G4int)(*theMateri 455 thisMaterialZ[i] = (G4int)(*theMaterialTable)[materialIndex]-> 454 GetEleme 456 GetElement(i)->GetZ(); 455 } 457 } 456 // fSandia = new G4SandiaTable(materialIndex 458 // fSandia = new G4SandiaTable(materialIndex); 457 fSandia = (*theMaterialTable)[materialIndex] 459 fSandia = (*theMaterialTable)[materialIndex]->GetSandiaTable(); 458 G4SandiaTable thisMaterialSandiaTable(ma 460 G4SandiaTable thisMaterialSandiaTable(materialIndex); 459 fIntervalNumber = thisMaterialSandiaTable.Sa 461 fIntervalNumber = thisMaterialSandiaTable.SandiaIntervals(thisMaterialZ, 460 << 462 numberOfElements); 461 fIntervalNumber = thisMaterialSandiaTable.Sa 463 fIntervalNumber = thisMaterialSandiaTable.SandiaMixing 462 ( thisMaterialZ , 464 ( thisMaterialZ , 463 (*theMaterialTable)[mate 465 (*theMaterialTable)[materialIndex]->GetFractionVector() , 464 numberOfElements, << 466 numberOfElements,fIntervalNumber); 465 467 466 fIntervalNumber--; 468 fIntervalNumber--; 467 469 468 fEnergyInterval = G4DataVector(fIntervalNumb 470 fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0); 469 fA1 = G4DataVector(fIntervalNumb 471 fA1 = G4DataVector(fIntervalNumber+2,0.0); 470 fA2 = G4DataVector(fIntervalNumb 472 fA2 = G4DataVector(fIntervalNumber+2,0.0); 471 fA3 = G4DataVector(fIntervalNumb 473 fA3 = G4DataVector(fIntervalNumber+2,0.0); 472 fA4 = G4DataVector(fIntervalNumb 474 fA4 = G4DataVector(fIntervalNumber+2,0.0); 473 475 474 /* 476 /* 475 fEnergyInterval = new G4double[fInterval 477 fEnergyInterval = new G4double[fIntervalNumber+2]; 476 fA1 = new G4double[fInterval 478 fA1 = new G4double[fIntervalNumber+2]; 477 fA2 = new G4double[fInterval 479 fA2 = new G4double[fIntervalNumber+2]; 478 fA3 = new G4double[fInterval 480 fA3 = new G4double[fIntervalNumber+2]; 479 fA4 = new G4double[fInterval 481 fA4 = new G4double[fIntervalNumber+2]; 480 */ 482 */ 481 for( i = 1; i <= fIntervalNumber; i++ ) 483 for( i = 1; i <= fIntervalNumber; i++ ) 482 { 484 { 483 if((thisMaterialSandiaTable.GetPhotoAbsorpCo 485 if((thisMaterialSandiaTable.GetPhotoAbsorpCof(i,0) >= maxEnergyTransfer) || 484 i > fIntervalNumber) 486 i > fIntervalNumber) 485 { 487 { 486 fEnergyInterval[i] = maxEnergyTran 488 fEnergyInterval[i] = maxEnergyTransfer; 487 fIntervalNumber = i; << 489 fIntervalNumber = i; 488 break; << 490 break; 489 } 491 } 490 fEnergyInterval[i] = thisMaterialSandiaTabl 492 fEnergyInterval[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,0); 491 fA1[i] = thisMaterialSandiaTabl 493 fA1[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,1)*fDensity; 492 fA2[i] = thisMaterialSandiaTabl 494 fA2[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,2)*fDensity; 493 fA3[i] = thisMaterialSandiaTabl 495 fA3[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,3)*fDensity; 494 fA4[i] = thisMaterialSandiaTabl 496 fA4[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,4)*fDensity; 495 497 496 } 498 } 497 if(fEnergyInterval[fIntervalNumber] != maxEn 499 if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer) 498 { 500 { 499 fIntervalNumber++; 501 fIntervalNumber++; 500 fEnergyInterval[fIntervalNumber] = ma 502 fEnergyInterval[fIntervalNumber] = maxEnergyTransfer; 501 fA1[fIntervalNumber] = fA1[fIntervalN 503 fA1[fIntervalNumber] = fA1[fIntervalNumber-1]; 502 fA2[fIntervalNumber] = fA2[fIntervalN 504 fA2[fIntervalNumber] = fA2[fIntervalNumber-1]; 503 fA3[fIntervalNumber] = fA3[fIntervalN 505 fA3[fIntervalNumber] = fA3[fIntervalNumber-1]; 504 fA4[fIntervalNumber] = fA4[fIntervalN 506 fA4[fIntervalNumber] = fA4[fIntervalNumber-1]; 505 } 507 } 506 for(i=1;i<=fIntervalNumber;i++) 508 for(i=1;i<=fIntervalNumber;i++) 507 { 509 { 508 // G4cout<<fEnergyInterval[i]<<"\t"<<f << 510 // G4cout<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t" 509 // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4e << 511 // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl; 510 } 512 } 511 // Now checking, if two borders are too clos 513 // Now checking, if two borders are too close together 512 514 513 for( i = 1; i < fIntervalNumber; i++ ) 515 for( i = 1; i < fIntervalNumber; i++ ) 514 { 516 { 515 if(fEnergyInterval[i+1]-fEnergyInterva 517 if(fEnergyInterval[i+1]-fEnergyInterval[i] > 516 1.5*fDelta*(fEnergyInterval[i+1]+fE 518 1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i])) 517 { << 519 { 518 continue; 520 continue; 519 } << 521 } 520 else 522 else 521 { << 523 { 522 for( j = i; j < fIntervalNumber; j++ 524 for( j = i; j < fIntervalNumber; j++ ) 523 { << 525 { 524 fEnergyInterval[j] = fEnergyInterv 526 fEnergyInterval[j] = fEnergyInterval[j+1]; 525 fA1[j] = fA1[j+1]; 527 fA1[j] = fA1[j+1]; 526 fA2[j] = fA2[j+1]; 528 fA2[j] = fA2[j+1]; 527 fA3[j] = fA3[j+1]; 529 fA3[j] = fA3[j+1]; 528 fA4[j] = fA4[j+1]; 530 fA4[j] = fA4[j+1]; 529 } << 531 } 530 fIntervalNumber--; 532 fIntervalNumber--; 531 i--; 533 i--; 532 } << 534 } 533 } 535 } 534 536 535 /* ********************************* 537 /* ********************************* 536 fSplineEnergy = new G4double[fM 538 fSplineEnergy = new G4double[fMaxSplineSize]; 537 fRePartDielectricConst = new G4double[fM 539 fRePartDielectricConst = new G4double[fMaxSplineSize]; 538 fImPartDielectricConst = new G4double[fM 540 fImPartDielectricConst = new G4double[fMaxSplineSize]; 539 fIntegralTerm = new G4double[fM 541 fIntegralTerm = new G4double[fMaxSplineSize]; 540 fDifPAIxSection = new G4double[fM 542 fDifPAIxSection = new G4double[fMaxSplineSize]; 541 fIntegralPAIxSection = new G4double[fM 543 fIntegralPAIxSection = new G4double[fMaxSplineSize]; 542 544 543 for(i=0;i<fMaxSplineSize;i++) 545 for(i=0;i<fMaxSplineSize;i++) 544 { 546 { 545 fSplineEnergy[i] = 0.0; 547 fSplineEnergy[i] = 0.0; 546 fRePartDielectricConst[i] = 0.0; 548 fRePartDielectricConst[i] = 0.0; 547 fImPartDielectricConst[i] = 0.0; 549 fImPartDielectricConst[i] = 0.0; 548 fIntegralTerm[i] = 0.0; 550 fIntegralTerm[i] = 0.0; 549 fDifPAIxSection[i] = 0.0; 551 fDifPAIxSection[i] = 0.0; 550 fIntegralPAIxSection[i] = 0.0; 552 fIntegralPAIxSection[i] = 0.0; 551 } 553 } 552 */ //////////////////////// 554 */ //////////////////////// 553 555 554 // Preparation of fSplineEnergy array co 556 // Preparation of fSplineEnergy array corresponding to min ionisation, G~4 555 557 556 ComputeLowEnergyCof(); 558 ComputeLowEnergyCof(); 557 G4double betaGammaSqRef = 559 G4double betaGammaSqRef = 558 fLorentzFactor[fRefGammaNumber]*fLorentzFa 560 fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1; 559 561 560 NormShift(betaGammaSqRef); 562 NormShift(betaGammaSqRef); 561 SplainPAI(betaGammaSqRef); 563 SplainPAI(betaGammaSqRef); 562 564 563 // Preparation of integral PAI cross section 565 // Preparation of integral PAI cross section for input betaGammaSq 564 566 565 for(i = 1; i <= fSplineNumber; i++) 567 for(i = 1; i <= fSplineNumber; i++) 566 { 568 { 567 fDifPAIxSection[i] = DifPAIxSection(i 569 fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq); 568 fdNdxCerenkov[i] = PAIdNdxCerenkov( 570 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq); 569 fdNdxMM[i] = PAIdNdxMM(i,betaGammaS 571 fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq); 570 fdNdxPlasmon[i] = PAIdNdxPlasmon(i 572 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq); 571 fdNdxResonance[i] = PAIdNdxResonance 573 fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq); 572 } 574 } 573 IntegralPAIxSection(); 575 IntegralPAIxSection(); 574 IntegralCerenkov(); 576 IntegralCerenkov(); 575 IntegralMM(); 577 IntegralMM(); 576 IntegralPlasmon(); 578 IntegralPlasmon(); 577 IntegralResonance(); 579 IntegralResonance(); 578 } 580 } 579 581 580 ////////////////////////////////////////////// 582 //////////////////////////////////////////////////////////////////////////// 581 // 583 // 582 // Destructor 584 // Destructor 583 585 584 G4PAIxSection::~G4PAIxSection() 586 G4PAIxSection::~G4PAIxSection() 585 { 587 { 586 /* ************************ 588 /* ************************ 587 delete[] fSplineEnergy ; 589 delete[] fSplineEnergy ; 588 delete[] fRePartDielectricConst; 590 delete[] fRePartDielectricConst; 589 delete[] fImPartDielectricConst; 591 delete[] fImPartDielectricConst; 590 delete[] fIntegralTerm ; 592 delete[] fIntegralTerm ; 591 delete[] fDifPAIxSection ; 593 delete[] fDifPAIxSection ; 592 delete[] fIntegralPAIxSection ; 594 delete[] fIntegralPAIxSection ; 593 */ //////////////////////// 595 */ //////////////////////// 594 delete fMatSandiaMatrix; 596 delete fMatSandiaMatrix; 595 } 597 } 596 598 597 G4double G4PAIxSection::GetLorentzFactor(G4int << 599 598 { << 600 599 return fLorentzFactor[j]; << 600 } << 601 601 602 ////////////////////////////////////////////// 602 //////////////////////////////////////////////////////////////////////// 603 // 603 // 604 // Constructor with beta*gamma square value ca 604 // Constructor with beta*gamma square value called from G4PAIPhotModel/Data 605 605 606 void G4PAIxSection::Initialize( const G4Materi 606 void G4PAIxSection::Initialize( const G4Material* material, 607 G4double maxEn << 607 G4double maxEnergyTransfer, 608 G4double betaG << 608 G4double betaGammaSq, 609 G4SandiaTable* << 609 G4SandiaTable* sandia) 610 { 610 { 611 if(fVerbose > 0) 611 if(fVerbose > 0) 612 { 612 { 613 G4cout<<G4endl; 613 G4cout<<G4endl; 614 G4cout<<"G4PAIxSection::Initialize(...,G4S 614 G4cout<<"G4PAIxSection::Initialize(...,G4SandiaTable* sandia)"<<G4endl; 615 G4cout<<G4endl; 615 G4cout<<G4endl; 616 } 616 } 617 G4int i, j; 617 G4int i, j; 618 618 619 fSandia = sandia; 619 fSandia = sandia; 620 fIntervalNumber = sandia->GetMaxInterval(); 620 fIntervalNumber = sandia->GetMaxInterval(); 621 fDensity = material->GetDensity(); 621 fDensity = material->GetDensity(); 622 fElectronDensity = material->GetElectronDens 622 fElectronDensity = material->GetElectronDensity(); 623 623 624 // fIntervalNumber--; 624 // fIntervalNumber--; 625 625 626 if( fVerbose > 0 ) 626 if( fVerbose > 0 ) 627 { 627 { 628 G4cout<<"fDensity = "<<fDensity<<"\t"<<fEl 628 G4cout<<"fDensity = "<<fDensity<<"\t"<<fElectronDensity<<"\t fIntervalNumber = "<<fIntervalNumber<<G4endl; 629 } 629 } 630 fEnergyInterval = G4DataVector(fIntervalNumb 630 fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0); 631 fA1 = G4DataVector(fIntervalNumb 631 fA1 = G4DataVector(fIntervalNumber+2,0.0); 632 fA2 = G4DataVector(fIntervalNumb 632 fA2 = G4DataVector(fIntervalNumber+2,0.0); 633 fA3 = G4DataVector(fIntervalNumb 633 fA3 = G4DataVector(fIntervalNumber+2,0.0); 634 fA4 = G4DataVector(fIntervalNumb 634 fA4 = G4DataVector(fIntervalNumber+2,0.0); 635 635 636 for( i = 1; i <= fIntervalNumber; i++ ) 636 for( i = 1; i <= fIntervalNumber; i++ ) 637 { 637 { 638 if ( sandia->GetSandiaMatTablePAI(i-1,0) < 638 if ( sandia->GetSandiaMatTablePAI(i-1,0) < 1.*eV && sandia->GetLowerI1() == false) 639 { 639 { 640 fIntervalNumber--; 640 fIntervalNumber--; 641 continue; 641 continue; 642 } 642 } 643 if( ( sandia->GetSandiaMatTablePAI(i-1,0) 643 if( ( sandia->GetSandiaMatTablePAI(i-1,0) >= maxEnergyTransfer ) || i >= fIntervalNumber ) 644 { 644 { 645 fEnergyInterval[i] = maxEnergyTransfer; 645 fEnergyInterval[i] = maxEnergyTransfer; 646 fIntervalNumber = i; 646 fIntervalNumber = i; 647 break; 647 break; 648 } 648 } 649 fEnergyInterval[i] = sandia->GetSandiaMatT 649 fEnergyInterval[i] = sandia->GetSandiaMatTablePAI(i-1,0); 650 fA1[i] = sandia->GetSandiaMatT 650 fA1[i] = sandia->GetSandiaMatTablePAI(i-1,1); 651 fA2[i] = sandia->GetSandiaMatT 651 fA2[i] = sandia->GetSandiaMatTablePAI(i-1,2); 652 fA3[i] = sandia->GetSandiaMatT 652 fA3[i] = sandia->GetSandiaMatTablePAI(i-1,3); 653 fA4[i] = sandia->GetSandiaMatT 653 fA4[i] = sandia->GetSandiaMatTablePAI(i-1,4); 654 654 655 if( fVerbose > 0 ) 655 if( fVerbose > 0 ) 656 { 656 { 657 G4cout<<i<<"\t"<<fEnergyInterval[i]/ke 657 G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t" 658 <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4e 658 <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl; 659 } 659 } 660 } 660 } 661 if( fVerbose > 0 ) G4cout<<"last i = "<<i<<" 661 if( fVerbose > 0 ) G4cout<<"last i = "<<i<<"; "<<"fIntervalNumber = "<<fIntervalNumber<<G4endl; 662 662 663 if( fEnergyInterval[fIntervalNumber] != maxE 663 if( fEnergyInterval[fIntervalNumber] != maxEnergyTransfer ) 664 { 664 { 665 fIntervalNumber++; 665 fIntervalNumber++; 666 fEnergyInterval[fIntervalNumber] = maxEn 666 fEnergyInterval[fIntervalNumber] = maxEnergyTransfer; 667 } 667 } 668 if( fVerbose > 0 ) 668 if( fVerbose > 0 ) 669 { 669 { 670 for( i = 1; i <= fIntervalNumber; i++ ) 670 for( i = 1; i <= fIntervalNumber; i++ ) 671 { 671 { 672 G4cout<<i<<"\t"<<fEnergyInterval[i]/keV< 672 G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t" 673 <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl; 673 <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl; 674 } 674 } 675 } 675 } 676 if( fVerbose > 0 ) G4cout<<"Now checking, 676 if( fVerbose > 0 ) G4cout<<"Now checking, if two borders are too close together"<<G4endl; 677 677 678 for( i = 1; i < fIntervalNumber; i++ ) 678 for( i = 1; i < fIntervalNumber; i++ ) 679 { 679 { 680 if( fEnergyInterval[i+1]-fEnergyInterval[i 680 if( fEnergyInterval[i+1]-fEnergyInterval[i] > 681 1.5*fDelta*(fEnergyInterval[i+1]+fEne << 681 1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]) && fEnergyInterval[i] > 0.) continue; 682 else 682 else 683 { 683 { 684 if( fVerbose > 0 ) G4cout<<i<<"\t"<<fEn 684 if( fVerbose > 0 ) G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fEnergyInterval[i+1]/keV; 685 685 686 for( j = i; j < fIntervalNumber; j++ ) 686 for( j = i; j < fIntervalNumber; j++ ) 687 { 687 { 688 fEnergyInterval[j] = fEnergyInte << 688 fEnergyInterval[j] = fEnergyInterval[j+1]; 689 fA1[j] = fA1[j+1]; << 689 fA1[j] = fA1[j+1]; 690 fA2[j] = fA2[j+1]; << 690 fA2[j] = fA2[j+1]; 691 fA3[j] = fA3[j+1]; << 691 fA3[j] = fA3[j+1]; 692 fA4[j] = fA4[j+1]; << 692 fA4[j] = fA4[j+1]; 693 } 693 } 694 fIntervalNumber--; 694 fIntervalNumber--; 695 i--; 695 i--; 696 } 696 } 697 } 697 } 698 if( fVerbose > 0 ) 698 if( fVerbose > 0 ) 699 { 699 { 700 for( i = 1; i <= fIntervalNumber; i++ ) 700 for( i = 1; i <= fIntervalNumber; i++ ) 701 { 701 { 702 G4cout<<i<<"\t"<<fEnergyInterval[i]/keV< 702 G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t" 703 <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl; 703 <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl; 704 } 704 } 705 } 705 } 706 // Preparation of fSplineEnergy array corres 706 // Preparation of fSplineEnergy array corresponding to min ionisation, G~4 707 707 708 ComputeLowEnergyCof(material); 708 ComputeLowEnergyCof(material); 709 709 710 G4double betaGammaSqRef = 710 G4double betaGammaSqRef = 711 fLorentzFactor[fRefGammaNumber]*fLorentzFa 711 fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1; 712 712 713 NormShift(betaGammaSqRef); 713 NormShift(betaGammaSqRef); 714 SplainPAI(betaGammaSqRef); 714 SplainPAI(betaGammaSqRef); 715 715 716 // Preparation of integral PAI cross section 716 // Preparation of integral PAI cross section for input betaGammaSq 717 717 718 for( i = 1; i <= fSplineNumber; i++ ) 718 for( i = 1; i <= fSplineNumber; i++ ) 719 { 719 { 720 fDifPAIxSection[i] = DifPAIxSection(i,bet 720 fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq); 721 721 722 722 723 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,be 723 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq); 724 fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq); 724 fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq); 725 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,bet 725 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq); 726 fdNdxResonance[i] = PAIdNdxResonance(i,b 726 fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq); 727 } 727 } 728 IntegralPAIxSection(); 728 IntegralPAIxSection(); 729 IntegralCerenkov(); 729 IntegralCerenkov(); 730 IntegralMM(); 730 IntegralMM(); 731 IntegralPlasmon(); 731 IntegralPlasmon(); 732 IntegralResonance(); 732 IntegralResonance(); 733 733 734 for( i = 1; i <= fSplineNumber; i++ ) 734 for( i = 1; i <= fSplineNumber; i++ ) 735 { 735 { 736 if(fVerbose>0) G4cout<<i<<"; w = "<<fSplin 736 if(fVerbose>0) G4cout<<i<<"; w = "<<fSplineEnergy[i]/keV<<" keV; dN/dx_>w = "<<fIntegralPAIxSection[i]<<" 1/mm"<<G4endl; 737 } 737 } 738 } 738 } 739 739 740 740 741 ////////////////////////////////////////////// 741 ///////////////////////////////////////////////////////////////////////// 742 // 742 // 743 // Compute low energy cof. It reduces PAI xsc 743 // Compute low energy cof. It reduces PAI xsc for Lorentz factors less than 4. 744 // 744 // 745 745 746 void G4PAIxSection::ComputeLowEnergyCof(const 746 void G4PAIxSection::ComputeLowEnergyCof(const G4Material* material) 747 { 747 { 748 G4int i, numberOfElements = (G4int)material- << 748 G4int i, numberOfElements = material->GetNumberOfElements(); 749 G4double sumZ = 0., sumCof = 0.; 749 G4double sumZ = 0., sumCof = 0.; 750 750 751 static const G4double p0 = 1.20923e+00; 751 static const G4double p0 = 1.20923e+00; 752 static const G4double p1 = 3.53256e-01; 752 static const G4double p1 = 3.53256e-01; 753 static const G4double p2 = -1.45052e-03; 753 static const G4double p2 = -1.45052e-03; 754 754 755 G4double* thisMaterialZ = new G4double[num 755 G4double* thisMaterialZ = new G4double[numberOfElements]; 756 G4double* thisMaterialCof = new G4double[num 756 G4double* thisMaterialCof = new G4double[numberOfElements]; 757 757 758 for( i = 0; i < numberOfElements; ++i ) << 758 for( i = 0; i < numberOfElements; i++ ) 759 { 759 { 760 thisMaterialZ[i] = material->GetElement(i) 760 thisMaterialZ[i] = material->GetElement(i)->GetZ(); 761 sumZ += thisMaterialZ[i]; 761 sumZ += thisMaterialZ[i]; 762 thisMaterialCof[i] = p0+p1*thisMaterialZ[i 762 thisMaterialCof[i] = p0+p1*thisMaterialZ[i]+p2*thisMaterialZ[i]*thisMaterialZ[i]; 763 } 763 } 764 for( i = 0; i < numberOfElements; ++i ) << 764 for( i = 0; i < numberOfElements; i++ ) 765 { 765 { 766 sumCof += thisMaterialCof[i]*thisMaterialZ 766 sumCof += thisMaterialCof[i]*thisMaterialZ[i]/sumZ; 767 } 767 } 768 fLowEnergyCof = sumCof; 768 fLowEnergyCof = sumCof; 769 delete [] thisMaterialZ; 769 delete [] thisMaterialZ; 770 delete [] thisMaterialCof; 770 delete [] thisMaterialCof; 771 // G4cout<<"fLowEnergyCof = "<<fLowEnergyCof 771 // G4cout<<"fLowEnergyCof = "<<fLowEnergyCof<<G4endl; 772 } 772 } 773 773 774 ////////////////////////////////////////////// 774 ///////////////////////////////////////////////////////////////////////// 775 // 775 // 776 // Compute low energy cof. It reduces PAI xsc 776 // Compute low energy cof. It reduces PAI xsc for Lorentz factors less than 4. 777 // 777 // 778 778 779 void G4PAIxSection::ComputeLowEnergyCof() 779 void G4PAIxSection::ComputeLowEnergyCof() 780 { 780 { 781 const G4MaterialTable* theMaterialTable = G4 781 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable(); 782 G4int i, numberOfElements = (G4int)(*theMate << 782 G4int i, numberOfElements = (*theMaterialTable)[fMaterialIndex]->GetNumberOfElements(); 783 G4double sumZ = 0., sumCof = 0.; 783 G4double sumZ = 0., sumCof = 0.; 784 784 785 const G4double p0 = 1.20923e+00; 785 const G4double p0 = 1.20923e+00; 786 const G4double p1 = 3.53256e-01; 786 const G4double p1 = 3.53256e-01; 787 const G4double p2 = -1.45052e-03; 787 const G4double p2 = -1.45052e-03; 788 788 789 G4double* thisMaterialZ = new G4double[num 789 G4double* thisMaterialZ = new G4double[numberOfElements]; 790 G4double* thisMaterialCof = new G4double[num 790 G4double* thisMaterialCof = new G4double[numberOfElements]; 791 791 792 for( i = 0; i < numberOfElements; ++i ) << 792 for( i = 0; i < numberOfElements; i++ ) 793 { 793 { 794 thisMaterialZ[i] = (*theMaterialTable)[fMa 794 thisMaterialZ[i] = (*theMaterialTable)[fMaterialIndex]->GetElement(i)->GetZ(); 795 sumZ += thisMaterialZ[i]; 795 sumZ += thisMaterialZ[i]; 796 thisMaterialCof[i] = p0+p1*thisMaterialZ[i 796 thisMaterialCof[i] = p0+p1*thisMaterialZ[i]+p2*thisMaterialZ[i]*thisMaterialZ[i]; 797 } 797 } 798 for( i = 0; i < numberOfElements; ++i ) << 798 for( i = 0; i < numberOfElements; i++ ) 799 { 799 { 800 sumCof += thisMaterialCof[i]*thisMaterialZ 800 sumCof += thisMaterialCof[i]*thisMaterialZ[i]/sumZ; 801 } 801 } 802 fLowEnergyCof = sumCof; 802 fLowEnergyCof = sumCof; 803 // G4cout<<"fLowEnergyCof = "<<fLowEnergyCof 803 // G4cout<<"fLowEnergyCof = "<<fLowEnergyCof<<G4endl; 804 delete [] thisMaterialZ; 804 delete [] thisMaterialZ; 805 delete [] thisMaterialCof; 805 delete [] thisMaterialCof; 806 } 806 } 807 807 808 ////////////////////////////////////////////// 808 ///////////////////////////////////////////////////////////////////////// 809 // 809 // 810 // General control function for class G4PAIxSe 810 // General control function for class G4PAIxSection 811 // 811 // 812 812 813 void G4PAIxSection::InitPAI() 813 void G4PAIxSection::InitPAI() 814 { 814 { 815 G4int i; 815 G4int i; 816 G4double betaGammaSq = fLorentzFactor[fRefG 816 G4double betaGammaSq = fLorentzFactor[fRefGammaNumber]* 817 fLorentzFactor[fRefG 817 fLorentzFactor[fRefGammaNumber] - 1; 818 818 819 // Preparation of integral PAI cross sectio 819 // Preparation of integral PAI cross section for reference gamma 820 820 821 NormShift(betaGammaSq); 821 NormShift(betaGammaSq); 822 SplainPAI(betaGammaSq); 822 SplainPAI(betaGammaSq); 823 823 824 IntegralPAIxSection(); 824 IntegralPAIxSection(); 825 IntegralCerenkov(); 825 IntegralCerenkov(); 826 IntegralMM(); 826 IntegralMM(); 827 IntegralPlasmon(); 827 IntegralPlasmon(); 828 IntegralResonance(); 828 IntegralResonance(); 829 829 830 for(i = 0; i<= fSplineNumber; i++) 830 for(i = 0; i<= fSplineNumber; i++) 831 { 831 { 832 fPAItable[i][fRefGammaNumber] = fIntegra 832 fPAItable[i][fRefGammaNumber] = fIntegralPAIxSection[i]; 833 if(i != 0) 833 if(i != 0) 834 { 834 { 835 fPAItable[i][0] = fSplineEnergy[i]; << 835 fPAItable[i][0] = fSplineEnergy[i]; 836 } 836 } 837 } 837 } 838 fPAItable[0][0] = fSplineNumber; 838 fPAItable[0][0] = fSplineNumber; 839 839 840 for(G4int j = 1; j < 112; j++) // for 840 for(G4int j = 1; j < 112; j++) // for other gammas 841 { 841 { 842 if( j == fRefGammaNumber ) continue; 842 if( j == fRefGammaNumber ) continue; 843 843 844 betaGammaSq = fLorentzFactor[j]*fLorentz 844 betaGammaSq = fLorentzFactor[j]*fLorentzFactor[j] - 1; 845 845 846 for(i = 1; i <= fSplineNumber; i++) 846 for(i = 1; i <= fSplineNumber; i++) 847 { 847 { 848 fDifPAIxSection[i] = DifPAIxSection(i 848 fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq); 849 fdNdxCerenkov[i] = PAIdNdxCerenkov( 849 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq); 850 fdNdxMM[i] = PAIdNdxMM(i,betaGammaS 850 fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq); 851 fdNdxPlasmon[i] = PAIdNdxPlasmon(i 851 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq); 852 fdNdxResonance[i] = PAIdNdxResonance 852 fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq); 853 } 853 } 854 IntegralPAIxSection(); 854 IntegralPAIxSection(); 855 IntegralCerenkov(); 855 IntegralCerenkov(); 856 IntegralMM(); 856 IntegralMM(); 857 IntegralPlasmon(); 857 IntegralPlasmon(); 858 IntegralResonance(); 858 IntegralResonance(); 859 859 860 for(i = 0; i <= fSplineNumber; i++) 860 for(i = 0; i <= fSplineNumber; i++) 861 { 861 { 862 fPAItable[i][j] = fIntegralPAIxSectio 862 fPAItable[i][j] = fIntegralPAIxSection[i]; 863 } 863 } 864 } 864 } 865 865 866 } 866 } 867 867 868 ////////////////////////////////////////////// 868 /////////////////////////////////////////////////////////////////////// 869 // 869 // 870 // Shifting from borders to intervals Creation 870 // Shifting from borders to intervals Creation of first energy points 871 // 871 // 872 872 873 void G4PAIxSection::NormShift(G4double betaGam 873 void G4PAIxSection::NormShift(G4double betaGammaSq) 874 { 874 { 875 G4int i, j; 875 G4int i, j; 876 876 877 if(fVerbose>0) G4cout<<" G4PAIxSection: 877 if(fVerbose>0) G4cout<<" G4PAIxSection::NormShift call "<<G4endl; 878 878 879 879 880 for( i = 1; i <= fIntervalNumber-1; i++ ) 880 for( i = 1; i <= fIntervalNumber-1; i++ ) 881 { 881 { 882 for( j = 1; j <= 2; j++ ) 882 for( j = 1; j <= 2; j++ ) 883 { 883 { 884 fSplineNumber = (i-1)*2 + j; 884 fSplineNumber = (i-1)*2 + j; 885 885 886 if( j == 1 ) fSplineEnergy[fSplineNumber 886 if( j == 1 ) fSplineEnergy[fSplineNumber] = fEnergyInterval[i ]*(1+fDelta); 887 else fSplineEnergy[fSplineNumber 887 else fSplineEnergy[fSplineNumber] = fEnergyInterval[i+1]*(1-fDelta); 888 if(fVerbose>0) G4cout<<"cn = "<<fSplineN 888 if(fVerbose>0) G4cout<<"cn = "<<fSplineNumber<<"; "<<"w = "<<fSplineEnergy[fSplineNumber]/keV<<" keV"<<G4endl; 889 } 889 } 890 } 890 } 891 fIntegralTerm[1]=RutherfordIntegral(1,fEnerg 891 fIntegralTerm[1]=RutherfordIntegral(1,fEnergyInterval[1],fSplineEnergy[1]); 892 892 893 j = 1; 893 j = 1; 894 894 895 for( i = 2; i <= fSplineNumber; i++ ) 895 for( i = 2; i <= fSplineNumber; i++ ) 896 { 896 { 897 if( fSplineEnergy[i]<fEnergyInterval[j+1] 897 if( fSplineEnergy[i]<fEnergyInterval[j+1] ) 898 { 898 { 899 fIntegralTerm[i] = fIntegralTerm[i-1] 899 fIntegralTerm[i] = fIntegralTerm[i-1] + 900 RutherfordIntegral << 900 RutherfordIntegral(j,fSplineEnergy[i-1], 901 901 fSplineEnergy[i] ); 902 } 902 } 903 else 903 else 904 { 904 { 905 G4double x = RutherfordIntegral(j,fSpli 905 G4double x = RutherfordIntegral(j,fSplineEnergy[i-1], 906 fEn 906 fEnergyInterval[j+1] ); 907 j++; 907 j++; 908 fIntegralTerm[i] = fIntegralTerm[i-1] 908 fIntegralTerm[i] = fIntegralTerm[i-1] + x + 909 RutherfordIntegral << 909 RutherfordIntegral(j,fEnergyInterval[j], 910 910 fSplineEnergy[i] ); 911 } 911 } 912 if(fVerbose>0) G4cout<<i<<" Shift: w = "< 912 if(fVerbose>0) G4cout<<i<<" Shift: w = "<<fSplineEnergy[i]/keV<<" keV \t"<<fIntegralTerm[i]<<"\n"<<G4endl; 913 } 913 } 914 fNormalizationCof = 2*pi*pi*hbarc*hbarc*fine 914 fNormalizationCof = 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2; 915 fNormalizationCof *= fElectronDensity/fInteg 915 fNormalizationCof *= fElectronDensity/fIntegralTerm[fSplineNumber]; 916 916 917 // G4cout<<"fNormalizationCof = "<<fNormaliz 917 // G4cout<<"fNormalizationCof = "<<fNormalizationCof<<G4endl; 918 918 919 // Calculation of PAI differrential << 919 // Calculation of PAI differrential cross-section (1/(keV*cm)) 920 // in the energy points near borders << 920 // in the energy points near borders of energy intervals 921 921 922 for(G4int k = 1; k <= fIntervalNumber-1; k+ 922 for(G4int k = 1; k <= fIntervalNumber-1; k++ ) 923 { 923 { 924 for( j = 1; j <= 2; j++ ) 924 for( j = 1; j <= 2; j++ ) 925 { 925 { 926 i = (k-1)*2 + j; 926 i = (k-1)*2 + j; 927 fImPartDielectricConst[i] = fNormaliz 927 fImPartDielectricConst[i] = fNormalizationCof* 928 ImPartDie << 928 ImPartDielectricConst(k,fSplineEnergy[i]); 929 fRePartDielectricConst[i] = fNormaliz 929 fRePartDielectricConst[i] = fNormalizationCof* 930 RePartDie << 930 RePartDielectricConst(fSplineEnergy[i]); 931 fIntegralTerm[i] *= fNormalizationCof 931 fIntegralTerm[i] *= fNormalizationCof; 932 932 933 fDifPAIxSection[i] = DifPAIxSection(i 933 fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq); 934 fdNdxCerenkov[i] = PAIdNdxCerenkov( 934 fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq); 935 fdNdxMM[i] = PAIdNdxMM(i,betaGammaS 935 fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq); 936 fdNdxPlasmon[i] = PAIdNdxPlasmon(i 936 fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq); 937 fdNdxResonance[i] = PAIdNdxResonan 937 fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq); 938 if(fVerbose>0) G4cout<<i<<" Shift: w = "< 938 if(fVerbose>0) G4cout<<i<<" Shift: w = "<<fSplineEnergy[i]/keV<<" keV, xsc = "<<fDifPAIxSection[i]<<"\n"<<G4endl; 939 } 939 } 940 } 940 } 941 941 942 } // end of NormShift 942 } // end of NormShift 943 943 944 ////////////////////////////////////////////// 944 ///////////////////////////////////////////////////////////////////////// 945 // 945 // 946 // Creation of new energy points as geometrica 946 // Creation of new energy points as geometrical mean of existing 947 // one, calculation PAI_cs for them, while the 947 // one, calculation PAI_cs for them, while the error of logarithmic 948 // linear approximation would be smaller than 948 // linear approximation would be smaller than 'fError' 949 949 950 void G4PAIxSection::SplainPAI(G4double betaGam 950 void G4PAIxSection::SplainPAI(G4double betaGammaSq) 951 { 951 { 952 G4int j, k = 1, i = 1; 952 G4int j, k = 1, i = 1; 953 953 954 if(fVerbose>0) G4cout<<" G 954 if(fVerbose>0) G4cout<<" G4PAIxSection::SplainPAI call "<<G4endl; 955 955 956 while ( (i < fSplineNumber) && (fSplineNumbe 956 while ( (i < fSplineNumber) && (fSplineNumber < fMaxSplineSize-1) ) 957 { 957 { 958 // if( std::abs(fSplineEnergy[i+1]-fEnerg 958 // if( std::abs(fSplineEnergy[i+1]-fEnergyInterval[k+1]) > (fSplineEnergy[i+1]+fEnergyInterval[k+1])*5.e-7 ) 959 if( fSplineEnergy[i+1] > fEnergyInterval[ 959 if( fSplineEnergy[i+1] > fEnergyInterval[k+1] ) 960 { 960 { 961 k++; // Here next energy point is 961 k++; // Here next energy point is in next energy interval 962 i++; << 962 i++; 963 if(fVerbose>0) G4cout<<" << 963 if(fVerbose>0) G4cout<<" in if: i = "<<i<<"; k = "<<k<<G4endl; 964 continue; 964 continue; 965 } 965 } 966 if(fVerbose>0) G4cout<<" out if: i 966 if(fVerbose>0) G4cout<<" out if: i = "<<i<<"; k = "<<k<<G4endl; 967 967 968 // Shifting of arrayes << 968 // Shifting of arrayes for inserting the geometrical 969 // average of 'i' and ' << 969 // average of 'i' and 'i+1' energy points to 'i+1' place 970 fSplineNumber++; 970 fSplineNumber++; 971 971 972 for( j = fSplineNumber; j >= i+2; j-- ) 972 for( j = fSplineNumber; j >= i+2; j-- ) 973 { 973 { 974 fSplineEnergy[j] = fSplineEn 974 fSplineEnergy[j] = fSplineEnergy[j-1]; 975 fImPartDielectricConst[j] = fImPartDi 975 fImPartDielectricConst[j] = fImPartDielectricConst[j-1]; 976 fRePartDielectricConst[j] = fRePartDi << 976 fRePartDielectricConst[j] = fRePartDielectricConst[j-1]; 977 fIntegralTerm[j] = fIntegral << 977 fIntegralTerm[j] = fIntegralTerm[j-1]; 978 978 979 fDifPAIxSection[j] = fDifPAIxSection[ << 979 fDifPAIxSection[j] = fDifPAIxSection[j-1]; 980 fdNdxCerenkov[j] = fdNdxCerenkov[j- 980 fdNdxCerenkov[j] = fdNdxCerenkov[j-1]; 981 fdNdxMM[j] = fdNdxMM[j-1]; 981 fdNdxMM[j] = fdNdxMM[j-1]; 982 fdNdxPlasmon[j] = fdNdxPlasmon[j-1 982 fdNdxPlasmon[j] = fdNdxPlasmon[j-1]; 983 fdNdxResonance[j] = fdNdxResonance[j 983 fdNdxResonance[j] = fdNdxResonance[j-1]; 984 } 984 } 985 G4double x1 = fSplineEnergy[i]; 985 G4double x1 = fSplineEnergy[i]; 986 G4double x2 = fSplineEnergy[i+1]; 986 G4double x2 = fSplineEnergy[i+1]; 987 G4double yy1 = fDifPAIxSection[i]; 987 G4double yy1 = fDifPAIxSection[i]; 988 G4double y2 = fDifPAIxSection[i+1]; 988 G4double y2 = fDifPAIxSection[i+1]; 989 989 990 if(fVerbose>0) G4cout<<"Spline: x1 = "<< 990 if(fVerbose>0) G4cout<<"Spline: x1 = "<<x1<<"; x2 = "<<x2<<", yy1 = "<<yy1<<"; y2 = "<<y2<<G4endl; 991 991 992 992 993 G4double en1 = sqrt(x1*x2); 993 G4double en1 = sqrt(x1*x2); 994 // G4double en1 = 0.5*(x1 + x2); 994 // G4double en1 = 0.5*(x1 + x2); 995 995 996 996 997 fSplineEnergy[i+1] = en1; 997 fSplineEnergy[i+1] = en1; 998 998 999 // Calculation of logarithmic << 999 // Calculation of logarithmic linear approximation 1000 // in this (enr) energy poin << 1000 // in this (enr) energy point, which number is 'i+1' now 1001 1001 1002 G4double a = log10(y2/yy1)/log10(x2/x1) 1002 G4double a = log10(y2/yy1)/log10(x2/x1); 1003 G4double b = log10(yy1) - a*log10(x1); 1003 G4double b = log10(yy1) - a*log10(x1); 1004 G4double y = a*log10(en1) + b; 1004 G4double y = a*log10(en1) + b; 1005 1005 1006 y = pow(10.,y); 1006 y = pow(10.,y); 1007 1007 1008 // Calculation of the PAI di << 1008 // Calculation of the PAI dif. cross-section at this point 1009 1009 1010 fImPartDielectricConst[i+1] = fNormaliz 1010 fImPartDielectricConst[i+1] = fNormalizationCof* 1011 ImPartDie << 1011 ImPartDielectricConst(k,fSplineEnergy[i+1]); 1012 fRePartDielectricConst[i+1] = fNormaliz 1012 fRePartDielectricConst[i+1] = fNormalizationCof* 1013 RePartDie << 1013 RePartDielectricConst(fSplineEnergy[i+1]); 1014 fIntegralTerm[i+1] = fIntegralTerm[i] + 1014 fIntegralTerm[i+1] = fIntegralTerm[i] + fNormalizationCof* 1015 RutherfordIntegral << 1015 RutherfordIntegral(k,fSplineEnergy[i], 1016 1016 fSplineEnergy[i+1]); 1017 1017 1018 fDifPAIxSection[i+1] = DifPAIxSection(i 1018 fDifPAIxSection[i+1] = DifPAIxSection(i+1,betaGammaSq); 1019 fdNdxCerenkov[i+1] = PAIdNdxCerenkov( 1019 fdNdxCerenkov[i+1] = PAIdNdxCerenkov(i+1,betaGammaSq); 1020 fdNdxMM[i+1] = PAIdNdxMM(i+1,be 1020 fdNdxMM[i+1] = PAIdNdxMM(i+1,betaGammaSq); 1021 fdNdxPlasmon[i+1] = PAIdNdxPlasmon(i 1021 fdNdxPlasmon[i+1] = PAIdNdxPlasmon(i+1,betaGammaSq); 1022 fdNdxResonance[i+1] = PAIdNdxResonance 1022 fdNdxResonance[i+1] = PAIdNdxResonance(i+1,betaGammaSq); 1023 1023 1024 // Condition for next divis << 1024 // Condition for next division of this segment or to pass 1025 1025 1026 if(fVerbose>0) G4cout<<"Spline, a = "<<a< 1026 if(fVerbose>0) G4cout<<"Spline, a = "<<a<<"; b = "<<b<<"; new xsc = "<<y<<"; compxsc = "<<fDifPAIxSection[i+1]<<G4endl; 1027 1027 1028 // to higher energies << 1028 // to higher energies 1029 1029 1030 G4double x = 2*(fDifPAIxSection[i+1] - 1030 G4double x = 2*(fDifPAIxSection[i+1] - y)/(fDifPAIxSection[i+1] + y); 1031 1031 1032 G4double delta = 2.*(fSplineEnergy[i+1] 1032 G4double delta = 2.*(fSplineEnergy[i+1]-fSplineEnergy[i])/(fSplineEnergy[i+1]+fSplineEnergy[i]); 1033 1033 1034 if( x < 0 ) 1034 if( x < 0 ) 1035 { 1035 { 1036 x = -x; << 1036 x = -x; 1037 } 1037 } 1038 if( x > fError && fSplineNumber < fMaxS 1038 if( x > fError && fSplineNumber < fMaxSplineSize-1 && delta > 2.*fDelta ) 1039 { 1039 { 1040 continue; // next division << 1040 continue; // next division 1041 } 1041 } 1042 i += 2; // pass to next segment 1042 i += 2; // pass to next segment 1043 1043 1044 // Loop checking, 03-Aug-2015, Vladimir << 1045 } // close 'while' 1044 } // close 'while' 1046 1045 1047 } // end of SplainPAI 1046 } // end of SplainPAI 1048 1047 1049 1048 1050 ///////////////////////////////////////////// 1049 //////////////////////////////////////////////////////////////////// 1051 // 1050 // 1052 // Integration over electrons that could be c 1051 // Integration over electrons that could be considered 1053 // quasi-free at energy transfer of interest 1052 // quasi-free at energy transfer of interest 1054 1053 1055 G4double G4PAIxSection::RutherfordIntegral( G 1054 G4double G4PAIxSection::RutherfordIntegral( G4int k, 1056 G << 1055 G4double x1, 1057 << 1056 G4double x2 ) 1058 { 1057 { 1059 G4double c1, c2, c3; 1058 G4double c1, c2, c3; 1060 // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "< 1059 // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl; 1061 c1 = (x2 - x1)/x1/x2; 1060 c1 = (x2 - x1)/x1/x2; 1062 c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2; 1061 c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2; 1063 c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/ 1062 c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2; 1064 // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = " 1063 // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<<c2<<"; "<<"c3 = "<<c3<<G4endl; 1065 1064 1066 return fA1[k]*log(x2/x1) + fA2[k]*c1 + fA 1065 return fA1[k]*log(x2/x1) + fA2[k]*c1 + fA3[k]*c2/2 + fA4[k]*c3/3; 1067 1066 1068 } // end of RutherfordIntegral 1067 } // end of RutherfordIntegral 1069 1068 1070 1069 1071 ///////////////////////////////////////////// 1070 ///////////////////////////////////////////////////////////////// 1072 // 1071 // 1073 // Imaginary part of dielectric constant 1072 // Imaginary part of dielectric constant 1074 // (G4int k - interval number, G4double en1 - 1073 // (G4int k - interval number, G4double en1 - energy point) 1075 1074 1076 G4double G4PAIxSection::ImPartDielectricConst 1075 G4double G4PAIxSection::ImPartDielectricConst( G4int k , 1077 << 1076 G4double energy1 ) 1078 { 1077 { 1079 G4double energy2,energy3,energy4,result; 1078 G4double energy2,energy3,energy4,result; 1080 1079 1081 energy2 = energy1*energy1; 1080 energy2 = energy1*energy1; 1082 energy3 = energy2*energy1; 1081 energy3 = energy2*energy1; 1083 energy4 = energy3*energy1; 1082 energy4 = energy3*energy1; 1084 1083 1085 result = fA1[k]/energy1+fA2[k]/energy2+fA3 1084 result = fA1[k]/energy1+fA2[k]/energy2+fA3[k]/energy3+fA4[k]/energy4; 1086 result *=hbarc/energy1; 1085 result *=hbarc/energy1; 1087 1086 1088 return result; 1087 return result; 1089 1088 1090 } // end of ImPartDielectricConst 1089 } // end of ImPartDielectricConst 1091 1090 1092 ///////////////////////////////////////////// 1091 ///////////////////////////////////////////////////////////////// 1093 // 1092 // 1094 // Returns lambda of photon with energy1 in c 1093 // Returns lambda of photon with energy1 in current material 1095 1094 1096 G4double G4PAIxSection::GetPhotonRange( G4dou 1095 G4double G4PAIxSection::GetPhotonRange( G4double energy1 ) 1097 { 1096 { 1098 G4int i; 1097 G4int i; 1099 G4double energy2, energy3, energy4, result, 1098 G4double energy2, energy3, energy4, result, lambda; 1100 1099 1101 energy2 = energy1*energy1; 1100 energy2 = energy1*energy1; 1102 energy3 = energy2*energy1; 1101 energy3 = energy2*energy1; 1103 energy4 = energy3*energy1; 1102 energy4 = energy3*energy1; 1104 1103 >> 1104 // G4double* SandiaCof = fSandia->GetSandiaCofForMaterialPAI(energy1); >> 1105 // result = SandiaCof[0]/energy1+SandiaCof[1]/energy2+SandiaCof[2]/energy3+SandiaCof[3]/energy4; >> 1106 // result *= fDensity; >> 1107 1105 for( i = 1; i <= fIntervalNumber; i++ ) 1108 for( i = 1; i <= fIntervalNumber; i++ ) 1106 { 1109 { 1107 if( energy1 < fEnergyInterval[i]) break; 1110 if( energy1 < fEnergyInterval[i]) break; 1108 } 1111 } 1109 i--; 1112 i--; 1110 if(i == 0) i = 1; 1113 if(i == 0) i = 1; 1111 1114 1112 result = fA1[i]/energy1+fA2[i]/energy2+fA3[ 1115 result = fA1[i]/energy1+fA2[i]/energy2+fA3[i]/energy3+fA4[i]/energy4; 1113 1116 1114 if( result > DBL_MIN ) lambda = 1./result; 1117 if( result > DBL_MIN ) lambda = 1./result; 1115 else lambda = DBL_MAX; 1118 else lambda = DBL_MAX; 1116 1119 1117 return lambda; 1120 return lambda; 1118 } 1121 } 1119 1122 1120 ///////////////////////////////////////////// 1123 ///////////////////////////////////////////////////////////////// 1121 // 1124 // 1122 // Return lambda of electron with energy1 in 1125 // Return lambda of electron with energy1 in current material 1123 // parametrisation from NIM A554(2005)474-493 1126 // parametrisation from NIM A554(2005)474-493 1124 1127 1125 G4double G4PAIxSection::GetElectronRange( G4d 1128 G4double G4PAIxSection::GetElectronRange( G4double energy ) 1126 { 1129 { 1127 G4double range; 1130 G4double range; 1128 /* 1131 /* 1129 const G4MaterialTable* theMaterialTable = G 1132 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable(); 1130 1133 1131 G4double Z = (*theMaterialTable)[fMaterialI 1134 G4double Z = (*theMaterialTable)[fMaterialIndex]->GetIonisation()->GetZeffective(); 1132 G4double A = (*theMaterialTable)[fMaterialI 1135 G4double A = (*theMaterialTable)[fMaterialIndex]->GetA(); 1133 1136 1134 energy /= keV; // energy in keV in parametr 1137 energy /= keV; // energy in keV in parametrised formula 1135 1138 1136 if (energy < 10.) 1139 if (energy < 10.) 1137 { 1140 { 1138 range = 3.872e-3*A/Z; 1141 range = 3.872e-3*A/Z; 1139 range *= pow(energy,1.492); 1142 range *= pow(energy,1.492); 1140 } 1143 } 1141 else 1144 else 1142 { 1145 { 1143 range = 6.97e-3*pow(energy,1.6); 1146 range = 6.97e-3*pow(energy,1.6); 1144 } 1147 } 1145 */ 1148 */ 1146 // Blum&Rolandi Particle Detection with Dri 1149 // Blum&Rolandi Particle Detection with Drift Chambers, p. 7 1147 1150 1148 G4double cofA = 5.37e-4*g/cm2/keV; 1151 G4double cofA = 5.37e-4*g/cm2/keV; 1149 G4double cofB = 0.9815; 1152 G4double cofB = 0.9815; 1150 G4double cofC = 3.123e-3/keV; 1153 G4double cofC = 3.123e-3/keV; 1151 // energy /= keV; 1154 // energy /= keV; 1152 1155 1153 range = cofA*energy*( 1 - cofB/(1 + cofC*en 1156 range = cofA*energy*( 1 - cofB/(1 + cofC*energy) ); 1154 1157 1155 // range *= g/cm2; 1158 // range *= g/cm2; 1156 range /= fDensity; 1159 range /= fDensity; 1157 1160 1158 return range; 1161 return range; 1159 } 1162 } 1160 1163 1161 ///////////////////////////////////////////// 1164 ////////////////////////////////////////////////////////////////////////////// 1162 // 1165 // 1163 // Real part of dielectric constant minus uni 1166 // Real part of dielectric constant minus unit: epsilon_1 - 1 1164 // (G4double enb - energy point) 1167 // (G4double enb - energy point) 1165 // 1168 // 1166 1169 1167 G4double G4PAIxSection::RePartDielectricConst 1170 G4double G4PAIxSection::RePartDielectricConst(G4double enb) 1168 { 1171 { 1169 G4double x0, x02, x03, x04, x05, x1, x2, x 1172 G4double x0, x02, x03, x04, x05, x1, x2, xx1 ,xx2 , xx12, 1170 c1, c2, c3, cof1, cof2, xln1, xln 1173 c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result; 1171 1174 1172 x0 = enb; 1175 x0 = enb; 1173 result = 0; 1176 result = 0; 1174 1177 1175 for(G4int i=1;i<=fIntervalNumber-1;i++) 1178 for(G4int i=1;i<=fIntervalNumber-1;i++) 1176 { 1179 { 1177 x1 = fEnergyInterval[i]; 1180 x1 = fEnergyInterval[i]; 1178 x2 = fEnergyInterval[i+1]; 1181 x2 = fEnergyInterval[i+1]; 1179 xx1 = x1 - x0; 1182 xx1 = x1 - x0; 1180 xx2 = x2 - x0; 1183 xx2 = x2 - x0; 1181 xx12 = xx2/xx1; 1184 xx12 = xx2/xx1; 1182 1185 1183 if(xx12<0) 1186 if(xx12<0) 1184 { 1187 { 1185 xx12 = -xx12; << 1188 xx12 = -xx12; 1186 } 1189 } 1187 xln1 = log(x2/x1); 1190 xln1 = log(x2/x1); 1188 xln2 = log(xx12); 1191 xln2 = log(xx12); 1189 xln3 = log((x2 + x0)/(x1 + x0)); 1192 xln3 = log((x2 + x0)/(x1 + x0)); 1190 x02 = x0*x0; 1193 x02 = x0*x0; 1191 x03 = x02*x0; 1194 x03 = x02*x0; 1192 x04 = x03*x0; 1195 x04 = x03*x0; 1193 x05 = x04*x0; 1196 x05 = x04*x0; 1194 c1 = (x2 - x1)/x1/x2; 1197 c1 = (x2 - x1)/x1/x2; 1195 c2 = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2; 1198 c2 = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2; 1196 c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/ 1199 c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2; 1197 1200 1198 result -= (fA1[i]/x02 + fA3[i]/x04)*xln 1201 result -= (fA1[i]/x02 + fA3[i]/x04)*xln1; 1199 result -= (fA2[i]/x02 + fA4[i]/x04)*c1; 1202 result -= (fA2[i]/x02 + fA4[i]/x04)*c1; 1200 result -= fA3[i]*c2/2/x02; 1203 result -= fA3[i]*c2/2/x02; 1201 result -= fA4[i]*c3/3/x02; 1204 result -= fA4[i]*c3/3/x02; 1202 1205 1203 cof1 = fA1[i]/x02 + fA3[i]/x04; 1206 cof1 = fA1[i]/x02 + fA3[i]/x04; 1204 cof2 = fA2[i]/x03 + fA4[i]/x05; 1207 cof2 = fA2[i]/x03 + fA4[i]/x05; 1205 1208 1206 result += 0.5*(cof1 +cof2)*xln2; 1209 result += 0.5*(cof1 +cof2)*xln2; 1207 result += 0.5*(cof1 - cof2)*xln3; 1210 result += 0.5*(cof1 - cof2)*xln3; 1208 } 1211 } 1209 result *= 2*hbarc/pi; 1212 result *= 2*hbarc/pi; 1210 1213 1211 return result; 1214 return result; 1212 1215 1213 } // end of RePartDielectricConst 1216 } // end of RePartDielectricConst 1214 1217 1215 ///////////////////////////////////////////// 1218 ////////////////////////////////////////////////////////////////////// 1216 // 1219 // 1217 // PAI differential cross-section in terms of 1220 // PAI differential cross-section in terms of 1218 // simplified Allison's equation 1221 // simplified Allison's equation 1219 // 1222 // 1220 1223 1221 G4double G4PAIxSection::DifPAIxSection( G4int << 1224 G4double G4PAIxSection::DifPAIxSection( G4int i , >> 1225 G4double betaGammaSq ) 1222 { 1226 { 1223 G4double cof,x1,x2,x3,x4,x5,x6,x7,x8,resul 1227 G4double cof,x1,x2,x3,x4,x5,x6,x7,x8,result; 1224 1228 1225 G4double betaBohr = fine_structure_const; 1229 G4double betaBohr = fine_structure_const; >> 1230 // G4double betaBohr2 = fine_structure_const*fine_structure_const; >> 1231 // G4double betaBohr3 = betaBohr*betaBohr2; // *4.0; >> 1232 1226 G4double be2 = betaGammaSq/(1 + betaGamma 1233 G4double be2 = betaGammaSq/(1 + betaGammaSq); 1227 G4double beta = std::sqrt(be2); << 1234 G4double beta = sqrt(be2); >> 1235 // G4double be3 = beta*be2; 1228 1236 1229 cof = 1.; 1237 cof = 1.; 1230 x1 = std::log(2*electron_mass_c2/fSplineE << 1238 x1 = log(2*electron_mass_c2/fSplineEnergy[i]); 1231 1239 1232 if( betaGammaSq < 0.01 ) x2 = std::log(be2 << 1240 if( betaGammaSq < 0.01 ) x2 = log(be2); 1233 else 1241 else 1234 { 1242 { 1235 x2 = -log( (1/betaGammaSq - fRePartDiele 1243 x2 = -log( (1/betaGammaSq - fRePartDielectricConst[i])* 1236 (1/betaGammaSq - fRePartDiele << 1244 (1/betaGammaSq - fRePartDielectricConst[i]) + 1237 fImPartDielectricConst[i]*fIm << 1245 fImPartDielectricConst[i]*fImPartDielectricConst[i] )/2; 1238 } 1246 } 1239 if( fImPartDielectricConst[i] == 0.0 ||bet 1247 if( fImPartDielectricConst[i] == 0.0 ||betaGammaSq < 0.01 ) 1240 { 1248 { 1241 x6 = 0.; 1249 x6 = 0.; 1242 } 1250 } 1243 else 1251 else 1244 { 1252 { 1245 x3 = -fRePartDielectricConst[i] + 1/beta 1253 x3 = -fRePartDielectricConst[i] + 1/betaGammaSq; 1246 x5 = -1 - fRePartDielectricConst[i] + 1254 x5 = -1 - fRePartDielectricConst[i] + 1247 be2*((1 +fRePartDielectricConst[i]) 1255 be2*((1 +fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) + 1248 fImPartDielectricConst[i]*fImPartDi << 1256 fImPartDielectricConst[i]*fImPartDielectricConst[i]); 1249 1257 1250 x7 = atan2(fImPartDielectricConst[i],x3) 1258 x7 = atan2(fImPartDielectricConst[i],x3); 1251 x6 = x5 * x7; 1259 x6 = x5 * x7; 1252 } 1260 } >> 1261 // if(fImPartDielectricConst[i] == 0) x6 = 0.; 1253 1262 1254 x4 = ((x1 + x2)*fImPartDielectricConst[i] 1263 x4 = ((x1 + x2)*fImPartDielectricConst[i] + x6)/hbarc; 1255 1264 >> 1265 // if( x4 < 0.0 ) x4 = 0.0; >> 1266 1256 x8 = (1 + fRePartDielectricConst[i])*(1 + 1267 x8 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) + 1257 fImPartDielectricConst[i]*fImPartDiel 1268 fImPartDielectricConst[i]*fImPartDielectricConst[i]; 1258 1269 1259 result = (x4 + cof*fIntegralTerm[i]/fSplin 1270 result = (x4 + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i]); 1260 1271 1261 if( result < 1.0e-8 ) result = 1.0e-8; 1272 if( result < 1.0e-8 ) result = 1.0e-8; 1262 1273 1263 result *= fine_structure_const/be2/pi; 1274 result *= fine_structure_const/be2/pi; 1264 1275 1265 // low energy correction 1276 // low energy correction 1266 1277 1267 G4double lowCof = fLowEnergyCof; // 6.0 ; 1278 G4double lowCof = fLowEnergyCof; // 6.0 ; // Ar ~ 4.; -> fLowCof as f(Z1,Z2)? 1268 1279 1269 result *= (1 - std::exp(-beta/betaBohr/low << 1280 result *= (1 - exp(-beta/betaBohr/lowCof)); 1270 if(x8 >= 0.0) << 1281 >> 1282 >> 1283 // result *= (1 - exp(-be2/betaBohr2/lowCof)); >> 1284 >> 1285 // result *= (1 - exp(-be3/betaBohr3/lowCof)); // ~ be for be<<betaBohr >> 1286 >> 1287 // result *= (1 - exp(-be4/betaBohr4/lowCof)); >> 1288 >> 1289 if(fDensity >= 0.1) 1271 { 1290 { 1272 result /= x8; << 1291 result /= x8; 1273 } 1292 } 1274 return result; 1293 return result; 1275 1294 1276 } // end of DifPAIxSection 1295 } // end of DifPAIxSection 1277 1296 1278 ///////////////////////////////////////////// 1297 ////////////////////////////////////////////////////////////////////////// 1279 // 1298 // 1280 // Calculation od dN/dx of collisions with cr 1299 // Calculation od dN/dx of collisions with creation of Cerenkov pseudo-photons 1281 1300 1282 G4double G4PAIxSection::PAIdNdxCerenkov( G4in 1301 G4double G4PAIxSection::PAIdNdxCerenkov( G4int i , 1283 G4do 1302 G4double betaGammaSq ) 1284 { 1303 { 1285 G4double logarithm, x3, x5, argument, modu 1304 G4double logarithm, x3, x5, argument, modul2, dNdxC; 1286 G4double be2, betaBohr2, cofBetaBohr; 1305 G4double be2, betaBohr2, cofBetaBohr; 1287 1306 1288 cofBetaBohr = 4.0; 1307 cofBetaBohr = 4.0; 1289 betaBohr2 = fine_structure_const*fine_stru << 1308 betaBohr2 = fine_structure_const*fine_structure_const; 1290 G4double betaBohr4 = betaBohr2*betaBohr2*c << 1309 G4double betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr; 1291 1310 1292 be2 = betaGammaSq/(1 + betaGammaSq); 1311 be2 = betaGammaSq/(1 + betaGammaSq); 1293 G4double be4 = be2*be2; 1312 G4double be4 = be2*be2; 1294 1313 1295 if( betaGammaSq < 0.01 ) logarithm = std:: << 1314 if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq); // 0.0; 1296 else 1315 else 1297 { 1316 { 1298 logarithm = -log( (1/betaGammaSq - fReP 1317 logarithm = -log( (1/betaGammaSq - fRePartDielectricConst[i])* 1299 (1/betaGammaSq - fReP << 1318 (1/betaGammaSq - fRePartDielectricConst[i]) + 1300 fImPartDielectricCons << 1319 fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5; 1301 logarithm += log(1+1.0/betaGammaSq); 1320 logarithm += log(1+1.0/betaGammaSq); 1302 } 1321 } 1303 1322 1304 if( fImPartDielectricConst[i] == 0.0 || be 1323 if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 ) 1305 { 1324 { 1306 argument = 0.0; 1325 argument = 0.0; 1307 } 1326 } 1308 else 1327 else 1309 { 1328 { 1310 x3 = -fRePartDielectricConst[i] + 1.0/be 1329 x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq; 1311 x5 = -1.0 - fRePartDielectricConst[i] + 1330 x5 = -1.0 - fRePartDielectricConst[i] + 1312 be2*((1.0 +fRePartDielectricConst[i 1331 be2*((1.0 +fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) + 1313 fImPartDielectricConst[i]*fImPartDi << 1332 fImPartDielectricConst[i]*fImPartDielectricConst[i]); 1314 if( x3 == 0.0 ) argument = 0.5*pi; 1333 if( x3 == 0.0 ) argument = 0.5*pi; 1315 else argument = std::atan2(fI << 1334 else argument = atan2(fImPartDielectricConst[i],x3); 1316 argument *= x5 ; 1335 argument *= x5 ; 1317 } 1336 } 1318 dNdxC = ( logarithm*fImPartDielectricConst 1337 dNdxC = ( logarithm*fImPartDielectricConst[i] + argument )/hbarc; 1319 1338 1320 if(dNdxC < 1.0e-8) dNdxC = 1.0e-8; 1339 if(dNdxC < 1.0e-8) dNdxC = 1.0e-8; 1321 1340 1322 dNdxC *= fine_structure_const/be2/pi; 1341 dNdxC *= fine_structure_const/be2/pi; 1323 1342 1324 dNdxC *= (1-std::exp(-be4/betaBohr4)); << 1343 dNdxC *= (1-exp(-be4/betaBohr4)); 1325 1344 1326 modul2 = (1.0 + fRePartDielectricConst[i]) << 1345 if(fDensity >= 0.1) 1327 fImPartDielectricConst[i]*fImPartDielect << 1328 if(modul2 >= 0.0) << 1329 { 1346 { 1330 dNdxC /= modul2; << 1347 modul2 = (1.0 + fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) + >> 1348 fImPartDielectricConst[i]*fImPartDielectricConst[i]; >> 1349 dNdxC /= modul2; 1331 } 1350 } 1332 return dNdxC; 1351 return dNdxC; 1333 1352 1334 } // end of PAIdNdxCerenkov 1353 } // end of PAIdNdxCerenkov 1335 1354 1336 ///////////////////////////////////////////// 1355 ////////////////////////////////////////////////////////////////////////// 1337 // 1356 // 1338 // Calculation od dN/dx of collisions of MM w 1357 // Calculation od dN/dx of collisions of MM with creation of Cerenkov pseudo-photons 1339 1358 1340 G4double G4PAIxSection::PAIdNdxMM( G4int i 1359 G4double G4PAIxSection::PAIdNdxMM( G4int i , 1341 G4do 1360 G4double betaGammaSq ) 1342 { 1361 { 1343 G4double logarithm, x3, x5, argument, dNdx 1362 G4double logarithm, x3, x5, argument, dNdxC; 1344 G4double be2, be4, betaBohr2,betaBohr4,cof 1363 G4double be2, be4, betaBohr2,betaBohr4,cofBetaBohr; 1345 1364 1346 cofBetaBohr = 4.0; 1365 cofBetaBohr = 4.0; 1347 betaBohr2 = fine_structure_const*fine_st 1366 betaBohr2 = fine_structure_const*fine_structure_const; 1348 betaBohr4 = betaBohr2*betaBohr2*cofBetaB 1367 betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr; 1349 1368 1350 be2 = betaGammaSq/(1 + betaGammaSq); 1369 be2 = betaGammaSq/(1 + betaGammaSq); 1351 be4 = be2*be2; 1370 be4 = be2*be2; 1352 1371 1353 if( betaGammaSq < 0.01 ) logarithm = log(1 1372 if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq); // 0.0; 1354 else 1373 else 1355 { 1374 { 1356 logarithm = -log( (1/betaGammaSq - fReP 1375 logarithm = -log( (1/betaGammaSq - fRePartDielectricConst[i])* 1357 (1/betaGammaSq - fReP << 1376 (1/betaGammaSq - fRePartDielectricConst[i]) + 1358 fImPartDielectricCons << 1377 fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5; 1359 logarithm += log(1+1.0/betaGammaSq); 1378 logarithm += log(1+1.0/betaGammaSq); 1360 } 1379 } 1361 1380 1362 if( fImPartDielectricConst[i] == 0.0 || be 1381 if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 ) 1363 { 1382 { 1364 argument = 0.0; 1383 argument = 0.0; 1365 } 1384 } 1366 else 1385 else 1367 { 1386 { 1368 x3 = -fRePartDielectricConst[i] + 1.0/be 1387 x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq; 1369 x5 = be2*( 1.0 + fRePartDielectricConst[ 1388 x5 = be2*( 1.0 + fRePartDielectricConst[i] ) - 1.0; 1370 if( x3 == 0.0 ) argument = 0.5*pi; 1389 if( x3 == 0.0 ) argument = 0.5*pi; 1371 else argument = atan2(fImPart 1390 else argument = atan2(fImPartDielectricConst[i],x3); 1372 argument *= x5 ; 1391 argument *= x5 ; 1373 } 1392 } 1374 dNdxC = ( logarithm*fImPartDielectricConst 1393 dNdxC = ( logarithm*fImPartDielectricConst[i]*be2 + argument )/hbarc; 1375 1394 1376 if(dNdxC < 1.0e-8) dNdxC = 1.0e-8; 1395 if(dNdxC < 1.0e-8) dNdxC = 1.0e-8; 1377 1396 1378 dNdxC *= fine_structure_const/be2/pi; 1397 dNdxC *= fine_structure_const/be2/pi; 1379 1398 1380 dNdxC *= (1-std::exp(-be4/betaBohr4)); << 1399 dNdxC *= (1-exp(-be4/betaBohr4)); 1381 return dNdxC; 1400 return dNdxC; 1382 1401 1383 } // end of PAIdNdxMM 1402 } // end of PAIdNdxMM 1384 1403 1385 ///////////////////////////////////////////// 1404 ////////////////////////////////////////////////////////////////////////// 1386 // 1405 // 1387 // Calculation od dN/dx of collisions with cr 1406 // Calculation od dN/dx of collisions with creation of longitudinal EM 1388 // excitations (plasmons, delta-electrons) 1407 // excitations (plasmons, delta-electrons) 1389 1408 1390 G4double G4PAIxSection::PAIdNdxPlasmon( G4int 1409 G4double G4PAIxSection::PAIdNdxPlasmon( G4int i , 1391 G4dou 1410 G4double betaGammaSq ) 1392 { 1411 { 1393 G4double resonance, modul2, dNdxP, cof = 1 1412 G4double resonance, modul2, dNdxP, cof = 1.; 1394 G4double be2, betaBohr; 1413 G4double be2, betaBohr; 1395 1414 1396 betaBohr = fine_structure_const; 1415 betaBohr = fine_structure_const; 1397 be2 = betaGammaSq/(1 + betaGammaSq); 1416 be2 = betaGammaSq/(1 + betaGammaSq); 1398 1417 1399 G4double beta = std::sqrt(be2); << 1418 G4double beta = sqrt(be2); 1400 1419 1401 resonance = std::log(2*electron_mass_c2*be << 1420 resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]); 1402 resonance *= fImPartDielectricConst[i]/hba 1421 resonance *= fImPartDielectricConst[i]/hbarc; 1403 1422 >> 1423 1404 dNdxP = ( resonance + cof*fIntegralTerm[i] 1424 dNdxP = ( resonance + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i] ); 1405 1425 1406 if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8; 1426 if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8; 1407 1427 1408 dNdxP *= fine_structure_const/be2/pi; 1428 dNdxP *= fine_structure_const/be2/pi; 1409 1429 1410 dNdxP *= (1 - std::exp(-beta/betaBohr/fLo << 1430 dNdxP *= (1 - exp(-beta/betaBohr/fLowEnergyCof)); 1411 1431 1412 modul2 = (1 + fRePartDielectricConst[i])*( << 1432 // dNdxP *= (1-exp(-be4/betaBohr4)); 1413 fImPartDielectricConst[i]*fImPartDielect << 1433 1414 if( modul2 >= 0.0 ) << 1434 if( fDensity >= 0.1 ) 1415 { 1435 { >> 1436 modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) + >> 1437 fImPartDielectricConst[i]*fImPartDielectricConst[i]; 1416 dNdxP /= modul2; 1438 dNdxP /= modul2; 1417 } 1439 } 1418 return dNdxP; 1440 return dNdxP; 1419 1441 1420 } // end of PAIdNdxPlasmon 1442 } // end of PAIdNdxPlasmon 1421 1443 1422 ///////////////////////////////////////////// 1444 ////////////////////////////////////////////////////////////////////////// 1423 // 1445 // 1424 // Calculation od dN/dx of collisions with cr 1446 // Calculation od dN/dx of collisions with creation of longitudinal EM 1425 // resonance excitations (plasmons, delta-ele 1447 // resonance excitations (plasmons, delta-electrons) 1426 1448 1427 G4double G4PAIxSection::PAIdNdxResonance( G4i 1449 G4double G4PAIxSection::PAIdNdxResonance( G4int i , 1428 G4dou 1450 G4double betaGammaSq ) 1429 { 1451 { 1430 G4double resonance, modul2, dNdxP; 1452 G4double resonance, modul2, dNdxP; 1431 G4double be2, be4, betaBohr2, betaBohr4, c 1453 G4double be2, be4, betaBohr2, betaBohr4, cofBetaBohr; 1432 1454 1433 cofBetaBohr = 4.0; 1455 cofBetaBohr = 4.0; 1434 betaBohr2 = fine_structure_const*fine_st 1456 betaBohr2 = fine_structure_const*fine_structure_const; 1435 betaBohr4 = betaBohr2*betaBohr2*cofBetaB 1457 betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr; 1436 1458 1437 be2 = betaGammaSq/(1 + betaGammaSq); 1459 be2 = betaGammaSq/(1 + betaGammaSq); 1438 be4 = be2*be2; 1460 be4 = be2*be2; 1439 1461 1440 resonance = log(2*electron_mass_c2*be2/fSp 1462 resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]); 1441 resonance *= fImPartDielectricConst[i]/hba 1463 resonance *= fImPartDielectricConst[i]/hbarc; 1442 1464 >> 1465 1443 dNdxP = resonance; 1466 dNdxP = resonance; 1444 1467 1445 if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8; 1468 if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8; 1446 1469 1447 dNdxP *= fine_structure_const/be2/pi; 1470 dNdxP *= fine_structure_const/be2/pi; 1448 dNdxP *= (1 - std::exp(-be4/betaBohr4)); << 1471 dNdxP *= (1-exp(-be4/betaBohr4)); 1449 1472 1450 modul2 = (1 + fRePartDielectricConst[i])*( << 1473 if( fDensity >= 0.1 ) 1451 fImPartDielectricConst[i]*fImPartDielect << 1452 if( modul2 >= 0.0 ) << 1453 { 1474 { >> 1475 modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) + >> 1476 fImPartDielectricConst[i]*fImPartDielectricConst[i]; 1454 dNdxP /= modul2; 1477 dNdxP /= modul2; 1455 } 1478 } 1456 return dNdxP; 1479 return dNdxP; 1457 1480 1458 } // end of PAIdNdxResonance 1481 } // end of PAIdNdxResonance 1459 1482 1460 ///////////////////////////////////////////// 1483 //////////////////////////////////////////////////////////////////////// 1461 // 1484 // 1462 // Calculation of the PAI integral cross-sect 1485 // Calculation of the PAI integral cross-section 1463 // fIntegralPAIxSection[1] = specific primary 1486 // fIntegralPAIxSection[1] = specific primary ionisation, 1/cm 1464 // and fIntegralPAIxSection[0] = mean energy 1487 // and fIntegralPAIxSection[0] = mean energy loss per cm in keV/cm 1465 1488 1466 void G4PAIxSection::IntegralPAIxSection() 1489 void G4PAIxSection::IntegralPAIxSection() 1467 { 1490 { 1468 fIntegralPAIxSection[fSplineNumber] = 0; 1491 fIntegralPAIxSection[fSplineNumber] = 0; 1469 fIntegralPAIdEdx[fSplineNumber] = 0; 1492 fIntegralPAIdEdx[fSplineNumber] = 0; 1470 fIntegralPAIxSection[0] = 0; 1493 fIntegralPAIxSection[0] = 0; 1471 G4int i, k = fIntervalNumber -1; 1494 G4int i, k = fIntervalNumber -1; 1472 1495 1473 for( i = fSplineNumber-1; i >= 1; i--) 1496 for( i = fSplineNumber-1; i >= 1; i--) 1474 { 1497 { 1475 if(fSplineEnergy[i] >= fEnergyInterval[k] 1498 if(fSplineEnergy[i] >= fEnergyInterval[k]) 1476 { 1499 { 1477 fIntegralPAIxSection[i] = fIntegralPAIx 1500 fIntegralPAIxSection[i] = fIntegralPAIxSection[i+1] + SumOverInterval(i); 1478 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[ 1501 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] + SumOverIntervaldEdx(i); 1479 } 1502 } 1480 else 1503 else 1481 { 1504 { 1482 fIntegralPAIxSection[i] = fIntegralPAIx 1505 fIntegralPAIxSection[i] = fIntegralPAIxSection[i+1] + 1483 SumOverBor << 1506 SumOverBorder(i+1,fEnergyInterval[k]); 1484 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[ 1507 fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] + 1485 SumOverBor << 1508 SumOverBorderdEdx(i+1,fEnergyInterval[k]); 1486 k--; 1509 k--; 1487 } 1510 } 1488 if(fVerbose>0) G4cout<<"i = "<<i<<"; k = 1511 if(fVerbose>0) G4cout<<"i = "<<i<<"; k = "<<k<<"; intPAIxsc[i] = "<<fIntegralPAIxSection[i]<<G4endl; 1489 } 1512 } 1490 } // end of IntegralPAIxSection 1513 } // end of IntegralPAIxSection 1491 1514 1492 ///////////////////////////////////////////// 1515 //////////////////////////////////////////////////////////////////////// 1493 // 1516 // 1494 // Calculation of the PAI Cerenkov integral c 1517 // Calculation of the PAI Cerenkov integral cross-section 1495 // fIntegralCrenkov[1] = specific Crenkov ion 1518 // fIntegralCrenkov[1] = specific Crenkov ionisation, 1/cm 1496 // and fIntegralCerenkov[0] = mean Cerenkov l 1519 // and fIntegralCerenkov[0] = mean Cerenkov loss per cm in keV/cm 1497 1520 1498 void G4PAIxSection::IntegralCerenkov() 1521 void G4PAIxSection::IntegralCerenkov() 1499 { 1522 { 1500 G4int i, k; 1523 G4int i, k; 1501 fIntegralCerenkov[fSplineNumber] = 0; 1524 fIntegralCerenkov[fSplineNumber] = 0; 1502 fIntegralCerenkov[0] = 0; 1525 fIntegralCerenkov[0] = 0; 1503 k = fIntervalNumber -1; 1526 k = fIntervalNumber -1; 1504 1527 1505 for( i = fSplineNumber-1; i >= 1; i-- ) 1528 for( i = fSplineNumber-1; i >= 1; i-- ) 1506 { 1529 { 1507 if(fSplineEnergy[i] >= fEnergyInterval[ 1530 if(fSplineEnergy[i] >= fEnergyInterval[k]) 1508 { 1531 { 1509 fIntegralCerenkov[i] = fIntegralCeren 1532 fIntegralCerenkov[i] = fIntegralCerenkov[i+1] + SumOverInterCerenkov(i); 1510 // G4cout<<"int: i = "<<i<<"; sumC = << 1533 // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl; 1511 } 1534 } 1512 else 1535 else 1513 { 1536 { 1514 fIntegralCerenkov[i] = fIntegralCeren 1537 fIntegralCerenkov[i] = fIntegralCerenkov[i+1] + 1515 SumOverBor << 1538 SumOverBordCerenkov(i+1,fEnergyInterval[k]); 1516 k--; << 1539 k--; 1517 // G4cout<<"bord: i = "<<i<<"; sumC = << 1540 // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl; 1518 } 1541 } 1519 } 1542 } 1520 1543 1521 } // end of IntegralCerenkov 1544 } // end of IntegralCerenkov 1522 1545 1523 ///////////////////////////////////////////// 1546 //////////////////////////////////////////////////////////////////////// 1524 // 1547 // 1525 // Calculation of the PAI MM-Cerenkov integra 1548 // Calculation of the PAI MM-Cerenkov integral cross-section 1526 // fIntegralMM[1] = specific MM-Cerenkov ioni 1549 // fIntegralMM[1] = specific MM-Cerenkov ionisation, 1/cm 1527 // and fIntegralMM[0] = mean MM-Cerenkov loss 1550 // and fIntegralMM[0] = mean MM-Cerenkov loss per cm in keV/cm 1528 1551 1529 void G4PAIxSection::IntegralMM() 1552 void G4PAIxSection::IntegralMM() 1530 { 1553 { 1531 G4int i, k; 1554 G4int i, k; 1532 fIntegralMM[fSplineNumber] = 0; 1555 fIntegralMM[fSplineNumber] = 0; 1533 fIntegralMM[0] = 0; 1556 fIntegralMM[0] = 0; 1534 k = fIntervalNumber -1; 1557 k = fIntervalNumber -1; 1535 1558 1536 for( i = fSplineNumber-1; i >= 1; i-- ) 1559 for( i = fSplineNumber-1; i >= 1; i-- ) 1537 { 1560 { 1538 if(fSplineEnergy[i] >= fEnergyInterval[ 1561 if(fSplineEnergy[i] >= fEnergyInterval[k]) 1539 { 1562 { 1540 fIntegralMM[i] = fIntegralMM[i+1] + S 1563 fIntegralMM[i] = fIntegralMM[i+1] + SumOverInterMM(i); 1541 // G4cout<<"int: i = "<<i<<"; sumC = << 1564 // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralMM[i]<<G4endl; 1542 } 1565 } 1543 else 1566 else 1544 { 1567 { 1545 fIntegralMM[i] = fIntegralMM[i+1] + 1568 fIntegralMM[i] = fIntegralMM[i+1] + 1546 SumOverBor << 1569 SumOverBordMM(i+1,fEnergyInterval[k]); 1547 k--; << 1570 k--; 1548 // G4cout<<"bord: i = "<<i<<"; sumC = << 1571 // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralMM[i]<<G4endl; 1549 } 1572 } 1550 } 1573 } 1551 1574 1552 } // end of IntegralMM 1575 } // end of IntegralMM 1553 1576 1554 ///////////////////////////////////////////// 1577 //////////////////////////////////////////////////////////////////////// 1555 // 1578 // 1556 // Calculation of the PAI Plasmon integral cr 1579 // Calculation of the PAI Plasmon integral cross-section 1557 // fIntegralPlasmon[1] = splasmon primary ion 1580 // fIntegralPlasmon[1] = splasmon primary ionisation, 1/cm 1558 // and fIntegralPlasmon[0] = mean plasmon los 1581 // and fIntegralPlasmon[0] = mean plasmon loss per cm in keV/cm 1559 1582 1560 void G4PAIxSection::IntegralPlasmon() 1583 void G4PAIxSection::IntegralPlasmon() 1561 { 1584 { 1562 fIntegralPlasmon[fSplineNumber] = 0; 1585 fIntegralPlasmon[fSplineNumber] = 0; 1563 fIntegralPlasmon[0] = 0; 1586 fIntegralPlasmon[0] = 0; 1564 G4int k = fIntervalNumber -1; 1587 G4int k = fIntervalNumber -1; 1565 for(G4int i=fSplineNumber-1;i>=1;i--) 1588 for(G4int i=fSplineNumber-1;i>=1;i--) 1566 { 1589 { 1567 if(fSplineEnergy[i] >= fEnergyInterval[ 1590 if(fSplineEnergy[i] >= fEnergyInterval[k]) 1568 { 1591 { 1569 fIntegralPlasmon[i] = fIntegralPlasmo 1592 fIntegralPlasmon[i] = fIntegralPlasmon[i+1] + SumOverInterPlasmon(i); 1570 } 1593 } 1571 else 1594 else 1572 { 1595 { 1573 fIntegralPlasmon[i] = fIntegralPlasmo 1596 fIntegralPlasmon[i] = fIntegralPlasmon[i+1] + 1574 SumOverBor << 1597 SumOverBordPlasmon(i+1,fEnergyInterval[k]); 1575 k--; << 1598 k--; 1576 } 1599 } 1577 } 1600 } 1578 1601 1579 } // end of IntegralPlasmon 1602 } // end of IntegralPlasmon 1580 1603 1581 ///////////////////////////////////////////// 1604 //////////////////////////////////////////////////////////////////////// 1582 // 1605 // 1583 // Calculation of the PAI resonance integral 1606 // Calculation of the PAI resonance integral cross-section 1584 // fIntegralResonance[1] = resonance primary 1607 // fIntegralResonance[1] = resonance primary ionisation, 1/cm 1585 // and fIntegralResonance[0] = mean resonance 1608 // and fIntegralResonance[0] = mean resonance loss per cm in keV/cm 1586 1609 1587 void G4PAIxSection::IntegralResonance() 1610 void G4PAIxSection::IntegralResonance() 1588 { 1611 { 1589 fIntegralResonance[fSplineNumber] = 0; 1612 fIntegralResonance[fSplineNumber] = 0; 1590 fIntegralResonance[0] = 0; 1613 fIntegralResonance[0] = 0; 1591 G4int k = fIntervalNumber -1; 1614 G4int k = fIntervalNumber -1; 1592 for(G4int i=fSplineNumber-1;i>=1;i--) 1615 for(G4int i=fSplineNumber-1;i>=1;i--) 1593 { 1616 { 1594 if(fSplineEnergy[i] >= fEnergyInterval[ 1617 if(fSplineEnergy[i] >= fEnergyInterval[k]) 1595 { 1618 { 1596 fIntegralResonance[i] = fIntegralReso 1619 fIntegralResonance[i] = fIntegralResonance[i+1] + SumOverInterResonance(i); 1597 } 1620 } 1598 else 1621 else 1599 { 1622 { 1600 fIntegralResonance[i] = fIntegralReso 1623 fIntegralResonance[i] = fIntegralResonance[i+1] + 1601 SumOverBor << 1624 SumOverBordResonance(i+1,fEnergyInterval[k]); 1602 k--; << 1625 k--; 1603 } 1626 } 1604 } 1627 } 1605 1628 1606 } // end of IntegralResonance 1629 } // end of IntegralResonance 1607 1630 1608 ///////////////////////////////////////////// 1631 ////////////////////////////////////////////////////////////////////// 1609 // 1632 // 1610 // Calculation the PAI integral cross-section 1633 // Calculation the PAI integral cross-section inside 1611 // of interval of continuous values of photo- 1634 // of interval of continuous values of photo-ionisation 1612 // cross-section. Parameter 'i' is the numbe 1635 // cross-section. Parameter 'i' is the number of interval. 1613 1636 1614 G4double G4PAIxSection::SumOverInterval( G4in 1637 G4double G4PAIxSection::SumOverInterval( G4int i ) 1615 { 1638 { 1616 G4double x0,x1,y0,yy1,a,b,c,result; 1639 G4double x0,x1,y0,yy1,a,b,c,result; 1617 1640 1618 x0 = fSplineEnergy[i]; 1641 x0 = fSplineEnergy[i]; 1619 x1 = fSplineEnergy[i+1]; 1642 x1 = fSplineEnergy[i+1]; 1620 if(fVerbose>0) G4cout<<"SumOverInterval i= 1643 if(fVerbose>0) G4cout<<"SumOverInterval i= " << i << " x0 = "<<x0<<"; x1 = "<<x1<<G4endl; 1621 1644 1622 if( x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/( 1645 if( x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.; 1623 1646 1624 y0 = fDifPAIxSection[i]; 1647 y0 = fDifPAIxSection[i]; 1625 yy1 = fDifPAIxSection[i+1]; 1648 yy1 = fDifPAIxSection[i+1]; 1626 1649 1627 if(fVerbose>0) G4cout<<"x0 = "<<x0<<"; x1 1650 if(fVerbose>0) G4cout<<"x0 = "<<x0<<"; x1 = "<<x1<<", y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl; 1628 1651 1629 c = x1/x0; 1652 c = x1/x0; 1630 a = log10(yy1/y0)/log10(c); 1653 a = log10(yy1/y0)/log10(c); 1631 1654 1632 if(fVerbose>0) G4cout<<"SumOverInterval, a 1655 if(fVerbose>0) G4cout<<"SumOverInterval, a = "<<a<<"; c = "<<c<<G4endl; 1633 1656 1634 b = 0.0; << 1657 // b = log10(y0) - a*log10(x0); 1635 if(a < 20.) b = y0/pow(x0,a); << 1658 b = y0/pow(x0,a); 1636 << 1637 a += 1.; 1659 a += 1.; 1638 if( std::abs(a) < 1.e-6 ) << 1660 if( std::fabs(a) < 1.e-6 ) 1639 { 1661 { 1640 result = b*log(x1/x0); 1662 result = b*log(x1/x0); 1641 } 1663 } 1642 else 1664 else 1643 { 1665 { 1644 result = y0*(x1*pow(c,a-1) - x0)/a; 1666 result = y0*(x1*pow(c,a-1) - x0)/a; 1645 } 1667 } 1646 a += 1.; 1668 a += 1.; 1647 if( std::abs(a) < 1.e-6 ) << 1669 if( std::fabs(a) < 1.e-6 ) 1648 { 1670 { 1649 fIntegralPAIxSection[0] += b*log(x1/x0); << 1671 fIntegralPAIxSection[0] += b*log(x1/x0); 1650 } 1672 } 1651 else 1673 else 1652 { 1674 { 1653 fIntegralPAIxSection[0] += y0*(x1*x1*po 1675 fIntegralPAIxSection[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a; 1654 } 1676 } 1655 if(fVerbose>0) G4cout<<"SumOverInterval, r 1677 if(fVerbose>0) G4cout<<"SumOverInterval, result = "<<result<<G4endl; 1656 return result; 1678 return result; 1657 1679 1658 } // end of SumOverInterval 1680 } // end of SumOverInterval 1659 1681 1660 ///////////////////////////////// 1682 ///////////////////////////////// 1661 1683 1662 G4double G4PAIxSection::SumOverIntervaldEdx( 1684 G4double G4PAIxSection::SumOverIntervaldEdx( G4int i ) 1663 { 1685 { 1664 G4double x0,x1,y0,yy1,a,b,c,result; 1686 G4double x0,x1,y0,yy1,a,b,c,result; 1665 1687 1666 x0 = fSplineEnergy[i]; 1688 x0 = fSplineEnergy[i]; 1667 x1 = fSplineEnergy[i+1]; 1689 x1 = fSplineEnergy[i+1]; 1668 1690 1669 if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x 1691 if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.; 1670 1692 1671 y0 = fDifPAIxSection[i]; 1693 y0 = fDifPAIxSection[i]; 1672 yy1 = fDifPAIxSection[i+1]; 1694 yy1 = fDifPAIxSection[i+1]; 1673 c = x1/x0; 1695 c = x1/x0; 1674 a = log10(yy1/y0)/log10(c); 1696 a = log10(yy1/y0)/log10(c); 1675 << 1697 // b = log10(y0) - a*log10(x0); 1676 b = 0.0; << 1698 b = y0/pow(x0,a); 1677 if(a < 20.) b = y0/pow(x0,a); << 1678 << 1679 a += 2; 1699 a += 2; 1680 if(a == 0) 1700 if(a == 0) 1681 { 1701 { 1682 result = b*log(x1/x0); 1702 result = b*log(x1/x0); 1683 } 1703 } 1684 else 1704 else 1685 { 1705 { 1686 result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a 1706 result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a; 1687 } 1707 } 1688 return result; 1708 return result; 1689 1709 1690 } // end of SumOverInterval 1710 } // end of SumOverInterval 1691 1711 1692 ///////////////////////////////////////////// 1712 ////////////////////////////////////////////////////////////////////// 1693 // 1713 // 1694 // Calculation the PAI Cerenkov integral cros 1714 // Calculation the PAI Cerenkov integral cross-section inside 1695 // of interval of continuous values of photo- 1715 // of interval of continuous values of photo-ionisation Cerenkov 1696 // cross-section. Parameter 'i' is the numbe 1716 // cross-section. Parameter 'i' is the number of interval. 1697 1717 1698 G4double G4PAIxSection::SumOverInterCerenkov( 1718 G4double G4PAIxSection::SumOverInterCerenkov( G4int i ) 1699 { 1719 { 1700 G4double x0,x1,y0,yy1,a,b,c,result; 1720 G4double x0,x1,y0,yy1,a,b,c,result; 1701 1721 1702 x0 = fSplineEnergy[i]; 1722 x0 = fSplineEnergy[i]; 1703 x1 = fSplineEnergy[i+1]; 1723 x1 = fSplineEnergy[i+1]; 1704 1724 1705 if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x 1725 if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.; 1706 1726 1707 y0 = fdNdxCerenkov[i]; 1727 y0 = fdNdxCerenkov[i]; 1708 yy1 = fdNdxCerenkov[i+1]; 1728 yy1 = fdNdxCerenkov[i+1]; 1709 // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<< 1729 // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1 1710 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4e 1730 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl; 1711 1731 1712 c = x1/x0; 1732 c = x1/x0; 1713 a = log10(yy1/y0)/log10(c); 1733 a = log10(yy1/y0)/log10(c); 1714 << 1734 b = y0/pow(x0,a); 1715 if(a > 20.0) b = 0.0; << 1716 else b = y0/pow(x0,a); << 1717 1735 1718 a += 1.0; 1736 a += 1.0; 1719 if(a == 0) result = b*log(c); 1737 if(a == 0) result = b*log(c); 1720 else result = y0*(x1*pow(c,a-1) - x0 1738 else result = y0*(x1*pow(c,a-1) - x0)/a; 1721 a += 1.0; 1739 a += 1.0; 1722 1740 1723 if( a == 0 ) fIntegralCerenkov[0] += b*log 1741 if( a == 0 ) fIntegralCerenkov[0] += b*log(x1/x0); 1724 else fIntegralCerenkov[0] += y0*(x 1742 else fIntegralCerenkov[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a; 1725 // G4cout<<"a = "<<a<<"; b = "<<b<<"; res 1743 // G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl; 1726 return result; 1744 return result; 1727 1745 1728 } // end of SumOverInterCerenkov 1746 } // end of SumOverInterCerenkov 1729 1747 1730 ///////////////////////////////////////////// 1748 ////////////////////////////////////////////////////////////////////// 1731 // 1749 // 1732 // Calculation the PAI MM-Cerenkov integral c 1750 // Calculation the PAI MM-Cerenkov integral cross-section inside 1733 // of interval of continuous values of photo- 1751 // of interval of continuous values of photo-ionisation Cerenkov 1734 // cross-section. Parameter 'i' is the numbe 1752 // cross-section. Parameter 'i' is the number of interval. 1735 1753 1736 G4double G4PAIxSection::SumOverInterMM( G4int 1754 G4double G4PAIxSection::SumOverInterMM( G4int i ) 1737 { 1755 { 1738 G4double x0,x1,y0,yy1,a,b,c,result; 1756 G4double x0,x1,y0,yy1,a,b,c,result; 1739 1757 1740 x0 = fSplineEnergy[i]; 1758 x0 = fSplineEnergy[i]; 1741 x1 = fSplineEnergy[i+1]; 1759 x1 = fSplineEnergy[i+1]; 1742 1760 1743 if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x 1761 if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.; 1744 1762 1745 y0 = fdNdxMM[i]; 1763 y0 = fdNdxMM[i]; 1746 yy1 = fdNdxMM[i+1]; 1764 yy1 = fdNdxMM[i+1]; 1747 //G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<" 1765 //G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1 1748 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4e 1766 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl; 1749 1767 1750 c = x1/x0; 1768 c = x1/x0; 1751 //G4cout<<" c = "<<c<< " yy1/y0= " << yy1/ 1769 //G4cout<<" c = "<<c<< " yy1/y0= " << yy1/y0 <<G4endl; 1752 a = log10(yy1/y0)/log10(c); 1770 a = log10(yy1/y0)/log10(c); 1753 << 1771 if(a > 10.0) return 0.; 1754 b = 0.0; << 1772 b = y0/pow(x0,a); 1755 if(a < 20.) b = y0/pow(x0,a); << 1756 1773 1757 a += 1.0; 1774 a += 1.0; 1758 if(a == 0) result = b*log(c); 1775 if(a == 0) result = b*log(c); 1759 else result = y0*(x1*pow(c,a-1) - x0 1776 else result = y0*(x1*pow(c,a-1) - x0)/a; 1760 a += 1.0; 1777 a += 1.0; 1761 1778 1762 if( a == 0 ) fIntegralMM[0] += b*log(c); 1779 if( a == 0 ) fIntegralMM[0] += b*log(c); 1763 else fIntegralMM[0] += y0*(x1*x1*p 1780 else fIntegralMM[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a; 1764 //G4cout<<"a = "<<a<<"; b = "<<b<<"; resul 1781 //G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl; 1765 return result; 1782 return result; 1766 1783 1767 } // end of SumOverInterMM 1784 } // end of SumOverInterMM 1768 1785 1769 ///////////////////////////////////////////// 1786 ////////////////////////////////////////////////////////////////////// 1770 // 1787 // 1771 // Calculation the PAI Plasmon integral cross 1788 // Calculation the PAI Plasmon integral cross-section inside 1772 // of interval of continuous values of photo- 1789 // of interval of continuous values of photo-ionisation Plasmon 1773 // cross-section. Parameter 'i' is the numbe 1790 // cross-section. Parameter 'i' is the number of interval. 1774 1791 1775 G4double G4PAIxSection::SumOverInterPlasmon( 1792 G4double G4PAIxSection::SumOverInterPlasmon( G4int i ) 1776 { 1793 { 1777 G4double x0,x1,y0,yy1,a,b,c,result; 1794 G4double x0,x1,y0,yy1,a,b,c,result; 1778 1795 1779 x0 = fSplineEnergy[i]; 1796 x0 = fSplineEnergy[i]; 1780 x1 = fSplineEnergy[i+1]; 1797 x1 = fSplineEnergy[i+1]; 1781 1798 1782 if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x 1799 if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.; 1783 1800 1784 y0 = fdNdxPlasmon[i]; 1801 y0 = fdNdxPlasmon[i]; 1785 yy1 = fdNdxPlasmon[i+1]; 1802 yy1 = fdNdxPlasmon[i+1]; 1786 c = x1/x0; << 1803 c =x1/x0; 1787 a = log10(yy1/y0)/log10(c); 1804 a = log10(yy1/y0)/log10(c); 1788 << 1805 if(a > 10.0) return 0.; 1789 b = 0.0; << 1806 // b = log10(y0) - a*log10(x0); 1790 if(a < 20.) b = y0/pow(x0,a); << 1807 b = y0/pow(x0,a); 1791 1808 1792 a += 1.0; 1809 a += 1.0; 1793 if(a == 0) result = b*log(x1/x0); 1810 if(a == 0) result = b*log(x1/x0); 1794 else result = y0*(x1*pow(c,a-1) - x0 1811 else result = y0*(x1*pow(c,a-1) - x0)/a; 1795 a += 1.0; 1812 a += 1.0; 1796 1813 1797 if( a == 0 ) fIntegralPlasmon[0] += b*log( 1814 if( a == 0 ) fIntegralPlasmon[0] += b*log(x1/x0); 1798 else fIntegralPlasmon[0] += y0*(x1 1815 else fIntegralPlasmon[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a; 1799 1816 1800 return result; 1817 return result; 1801 1818 1802 } // end of SumOverInterPlasmon 1819 } // end of SumOverInterPlasmon 1803 1820 1804 ///////////////////////////////////////////// 1821 ////////////////////////////////////////////////////////////////////// 1805 // 1822 // 1806 // Calculation the PAI resonance integral cro 1823 // Calculation the PAI resonance integral cross-section inside 1807 // of interval of continuous values of photo- 1824 // of interval of continuous values of photo-ionisation resonance 1808 // cross-section. Parameter 'i' is the numbe 1825 // cross-section. Parameter 'i' is the number of interval. 1809 1826 1810 G4double G4PAIxSection::SumOverInterResonance 1827 G4double G4PAIxSection::SumOverInterResonance( G4int i ) 1811 { 1828 { 1812 G4double x0,x1,y0,yy1,a,b,c,result; 1829 G4double x0,x1,y0,yy1,a,b,c,result; 1813 1830 1814 x0 = fSplineEnergy[i]; 1831 x0 = fSplineEnergy[i]; 1815 x1 = fSplineEnergy[i+1]; 1832 x1 = fSplineEnergy[i+1]; 1816 1833 1817 if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x 1834 if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.; 1818 1835 1819 y0 = fdNdxResonance[i]; 1836 y0 = fdNdxResonance[i]; 1820 yy1 = fdNdxResonance[i+1]; 1837 yy1 = fdNdxResonance[i+1]; 1821 c =x1/x0; 1838 c =x1/x0; 1822 a = log10(yy1/y0)/log10(c); 1839 a = log10(yy1/y0)/log10(c); 1823 << 1840 if(a > 10.0) return 0.; 1824 b = 0.0; << 1841 // b = log10(y0) - a*log10(x0); 1825 if(a < 20.) b = y0/pow(x0,a); << 1842 b = y0/pow(x0,a); 1826 1843 1827 a += 1.0; 1844 a += 1.0; 1828 if(a == 0) result = b*log(x1/x0); 1845 if(a == 0) result = b*log(x1/x0); 1829 else result = y0*(x1*pow(c,a-1) - x0 1846 else result = y0*(x1*pow(c,a-1) - x0)/a; 1830 a += 1.0; 1847 a += 1.0; 1831 1848 1832 if( a == 0 ) fIntegralResonance[0] += b*lo 1849 if( a == 0 ) fIntegralResonance[0] += b*log(x1/x0); 1833 else fIntegralResonance[0] += y0*( 1850 else fIntegralResonance[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a; 1834 1851 1835 return result; 1852 return result; 1836 1853 1837 } // end of SumOverInterResonance 1854 } // end of SumOverInterResonance 1838 1855 1839 ///////////////////////////////////////////// 1856 /////////////////////////////////////////////////////////////////////////////// 1840 // 1857 // 1841 // Integration of PAI cross-section for the c 1858 // Integration of PAI cross-section for the case of 1842 // passing across border between intervals 1859 // passing across border between intervals 1843 1860 1844 G4double G4PAIxSection::SumOverBorder( G4int 1861 G4double G4PAIxSection::SumOverBorder( G4int i , 1845 G4doub 1862 G4double en0 ) 1846 { 1863 { 1847 G4double x0,x1,y0,yy1,a,b,/*c,*/d,e0,result 1864 G4double x0,x1,y0,yy1,a,b,/*c,*/d,e0,result; 1848 1865 1849 e0 = en0; 1866 e0 = en0; 1850 x0 = fSplineEnergy[i]; 1867 x0 = fSplineEnergy[i]; 1851 x1 = fSplineEnergy[i+1]; 1868 x1 = fSplineEnergy[i+1]; 1852 y0 = fDifPAIxSection[i]; 1869 y0 = fDifPAIxSection[i]; 1853 yy1 = fDifPAIxSection[i+1]; 1870 yy1 = fDifPAIxSection[i+1]; 1854 1871 1855 //c = x1/x0; 1872 //c = x1/x0; 1856 d = e0/x0; 1873 d = e0/x0; 1857 a = log10(yy1/y0)/log10(x1/x0); 1874 a = log10(yy1/y0)/log10(x1/x0); >> 1875 if(a > 10.0) return 0.; 1858 1876 1859 if(fVerbose>0) G4cout<<"SumOverBorder, a = 1877 if(fVerbose>0) G4cout<<"SumOverBorder, a = "<<a<<G4endl; 1860 1878 1861 b = 0.0; << 1879 // b0 = log10(y0) - a*log10(x0); 1862 if(a < 20.) b = y0/pow(x0,a); << 1880 b = y0/pow(x0,a); // pow(10.,b); 1863 1881 1864 a += 1.; 1882 a += 1.; 1865 if( std::abs(a) < 1.e-6 ) << 1883 if( std::fabs(a) < 1.e-6 ) 1866 { 1884 { 1867 result = b*log(x0/e0); 1885 result = b*log(x0/e0); 1868 } 1886 } 1869 else 1887 else 1870 { 1888 { 1871 result = y0*(x0 - e0*pow(d,a-1))/a; 1889 result = y0*(x0 - e0*pow(d,a-1))/a; 1872 } 1890 } 1873 a += 1.; 1891 a += 1.; 1874 if( std::abs(a) < 1.e-6 ) << 1892 if( std::fabs(a) < 1.e-6 ) 1875 { 1893 { 1876 fIntegralPAIxSection[0] += b*log(x0/e0) 1894 fIntegralPAIxSection[0] += b*log(x0/e0); 1877 } 1895 } 1878 else 1896 else 1879 { 1897 { 1880 fIntegralPAIxSection[0] += y0*(x0*x0 - 1898 fIntegralPAIxSection[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a; 1881 } 1899 } 1882 x0 = fSplineEnergy[i - 1]; 1900 x0 = fSplineEnergy[i - 1]; 1883 x1 = fSplineEnergy[i - 2]; 1901 x1 = fSplineEnergy[i - 2]; 1884 y0 = fDifPAIxSection[i - 1]; 1902 y0 = fDifPAIxSection[i - 1]; 1885 yy1 = fDifPAIxSection[i - 2]; 1903 yy1 = fDifPAIxSection[i - 2]; 1886 1904 >> 1905 //c = x1/x0; 1887 d = e0/x0; 1906 d = e0/x0; 1888 a = log10(yy1/y0)/log10(x1/x0); 1907 a = log10(yy1/y0)/log10(x1/x0); 1889 << 1908 // b0 = log10(y0) - a*log10(x0); 1890 b = 0.0; << 1909 b = y0/pow(x0,a); 1891 if(a < 20.) b = y0/pow(x0,a); << 1892 << 1893 a += 1.; 1910 a += 1.; 1894 if( std::abs(a) < 1.e-6 ) << 1911 if( std::fabs(a) < 1.e-6 ) 1895 { 1912 { 1896 result += b*log(e0/x0); 1913 result += b*log(e0/x0); 1897 } 1914 } 1898 else 1915 else 1899 { 1916 { 1900 result += y0*(e0*pow(d,a-1) - x0)/a; 1917 result += y0*(e0*pow(d,a-1) - x0)/a; 1901 } 1918 } 1902 a += 1.; 1919 a += 1.; 1903 if( std::abs(a) < 1.e-6 ) << 1920 if( std::fabs(a) < 1.e-6 ) 1904 { 1921 { 1905 fIntegralPAIxSection[0] += b*log(e0/x0) 1922 fIntegralPAIxSection[0] += b*log(e0/x0); 1906 } 1923 } 1907 else 1924 else 1908 { 1925 { 1909 fIntegralPAIxSection[0] += y0*(e0*e0*po 1926 fIntegralPAIxSection[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a; 1910 } 1927 } 1911 return result; 1928 return result; 1912 1929 1913 } 1930 } 1914 1931 1915 ///////////////////////////////////////////// 1932 /////////////////////////////////////////////////////////////////////// 1916 1933 1917 G4double G4PAIxSection::SumOverBorderdEdx( G4 << 1934 G4double G4PAIxSection::SumOverBorderdEdx( G4int i , >> 1935 G4double en0 ) 1918 { 1936 { 1919 G4double x0,x1,y0,yy1,a,b,d,e0,result; << 1937 G4double x0,x1,y0,yy1,a,b,/*c,*/d,e0,result; 1920 1938 1921 e0 = en0; 1939 e0 = en0; 1922 x0 = fSplineEnergy[i]; 1940 x0 = fSplineEnergy[i]; 1923 x1 = fSplineEnergy[i+1]; 1941 x1 = fSplineEnergy[i+1]; 1924 y0 = fDifPAIxSection[i]; 1942 y0 = fDifPAIxSection[i]; 1925 yy1 = fDifPAIxSection[i+1]; 1943 yy1 = fDifPAIxSection[i+1]; 1926 1944 >> 1945 //c = x1/x0; 1927 d = e0/x0; 1946 d = e0/x0; 1928 a = log10(yy1/y0)/log10(x1/x0); 1947 a = log10(yy1/y0)/log10(x1/x0); 1929 << 1948 if(a > 10.0) return 0.; 1930 b = 0.0; << 1949 // b0 = log10(y0) - a*log10(x0); 1931 if(a < 20.) b = y0/pow(x0,a); << 1950 b = y0/pow(x0,a); // pow(10.,b); 1932 1951 1933 a += 2; 1952 a += 2; 1934 if(a == 0) 1953 if(a == 0) 1935 { 1954 { 1936 result = b*log(x0/e0); 1955 result = b*log(x0/e0); 1937 } 1956 } 1938 else 1957 else 1939 { 1958 { 1940 result = y0*(x0*x0 - e0*e0*pow(d,a-2))/ 1959 result = y0*(x0*x0 - e0*e0*pow(d,a-2))/a; 1941 } 1960 } 1942 x0 = fSplineEnergy[i - 1]; 1961 x0 = fSplineEnergy[i - 1]; 1943 x1 = fSplineEnergy[i - 2]; 1962 x1 = fSplineEnergy[i - 2]; 1944 y0 = fDifPAIxSection[i - 1]; 1963 y0 = fDifPAIxSection[i - 1]; 1945 yy1 = fDifPAIxSection[i - 2]; 1964 yy1 = fDifPAIxSection[i - 2]; 1946 1965 1947 // c = x1/x0; 1966 // c = x1/x0; 1948 d = e0/x0; 1967 d = e0/x0; 1949 a = log10(yy1/y0)/log10(x1/x0); 1968 a = log10(yy1/y0)/log10(x1/x0); 1950 << 1969 // b0 = log10(y0) - a*log10(x0); 1951 b = 0.0; << 1970 b = y0/pow(x0,a); 1952 if(a < 20.) b = y0/pow(x0,a); << 1953 << 1954 a += 2; 1971 a += 2; 1955 if(a == 0) 1972 if(a == 0) 1956 { 1973 { 1957 result += b*log(e0/x0); 1974 result += b*log(e0/x0); 1958 } 1975 } 1959 else 1976 else 1960 { 1977 { 1961 result += y0*(e0*e0*pow(d,a-2) - x0*x0) 1978 result += y0*(e0*e0*pow(d,a-2) - x0*x0)/a; 1962 } 1979 } 1963 return result; 1980 return result; 1964 1981 1965 } 1982 } 1966 1983 1967 ///////////////////////////////////////////// 1984 /////////////////////////////////////////////////////////////////////////////// 1968 // 1985 // 1969 // Integration of Cerenkov cross-section for 1986 // Integration of Cerenkov cross-section for the case of 1970 // passing across border between intervals 1987 // passing across border between intervals 1971 1988 1972 G4double G4PAIxSection::SumOverBordCerenkov( << 1989 G4double G4PAIxSection::SumOverBordCerenkov( G4int i , >> 1990 G4double en0 ) 1973 { 1991 { 1974 G4double x0,x1,y0,yy1,a,b,e0,c,d,result; 1992 G4double x0,x1,y0,yy1,a,b,e0,c,d,result; 1975 1993 1976 e0 = en0; 1994 e0 = en0; 1977 x0 = fSplineEnergy[i]; 1995 x0 = fSplineEnergy[i]; 1978 x1 = fSplineEnergy[i+1]; 1996 x1 = fSplineEnergy[i+1]; 1979 y0 = fdNdxCerenkov[i]; 1997 y0 = fdNdxCerenkov[i]; 1980 yy1 = fdNdxCerenkov[i+1]; 1998 yy1 = fdNdxCerenkov[i+1]; 1981 1999 1982 //G4cout<<"SumBordC, i = "<<i<<"; en0 = "< << 2000 // G4cout<<G4endl; 1983 //<<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl << 2001 // G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1 >> 2002 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl; 1984 c = x1/x0; 2003 c = x1/x0; 1985 d = e0/x0; 2004 d = e0/x0; 1986 a = log10(yy1/y0)/log10(c); 2005 a = log10(yy1/y0)/log10(c); 1987 //G4cout << " a= " << a << " c=" << c < << 2006 if(a > 10.0) return 0.; 1988 << 2007 // b0 = log10(y0) - a*log10(x0); 1989 b = 0.0; << 2008 b = y0/pow(x0,a); // pow(10.,b0); 1990 if(a < 20.) b = y0/pow(x0,a); << 1991 2009 1992 a += 1.0; 2010 a += 1.0; 1993 if( a == 0 ) result = b*log(x0/e0); 2011 if( a == 0 ) result = b*log(x0/e0); 1994 else result = y0*(x0 - e0*pow(d,a- 2012 else result = y0*(x0 - e0*pow(d,a-1))/a; 1995 a += 1.0; 2013 a += 1.0; 1996 2014 1997 if( a == 0 ) fIntegralCerenkov[0] += b*log 2015 if( a == 0 ) fIntegralCerenkov[0] += b*log(x0/e0); 1998 else fIntegralCerenkov[0] += y0*(x 2016 else fIntegralCerenkov[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a; 1999 2017 >> 2018 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "<<b<<"; result = "<<result<<G4endl; >> 2019 2000 x0 = fSplineEnergy[i - 1]; 2020 x0 = fSplineEnergy[i - 1]; 2001 x1 = fSplineEnergy[i - 2]; 2021 x1 = fSplineEnergy[i - 2]; 2002 y0 = fdNdxCerenkov[i - 1]; 2022 y0 = fdNdxCerenkov[i - 1]; 2003 yy1 = fdNdxCerenkov[i - 2]; 2023 yy1 = fdNdxCerenkov[i - 2]; 2004 2024 2005 // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1 2025 // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1 2006 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4 2026 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl; 2007 2027 2008 c = x1/x0; 2028 c = x1/x0; 2009 d = e0/x0; 2029 d = e0/x0; 2010 a = log10(yy1/y0)/log10(c); << 2030 a = log10(yy1/y0)/log10(x1/x0); 2011 << 2031 // b0 = log10(y0) - a*log10(x0); 2012 b = 0.0; << 2032 b = y0/pow(x0,a); // pow(10.,b0); 2013 if(a < 20.) b = y0/pow(x0,a); << 2014 2033 2015 a += 1.0; 2034 a += 1.0; 2016 if( a == 0 ) result += b*log(e0/x0); 2035 if( a == 0 ) result += b*log(e0/x0); 2017 else result += y0*(e0*pow(d,a-1) - 2036 else result += y0*(e0*pow(d,a-1) - x0 )/a; 2018 a += 1.0; 2037 a += 1.0; 2019 2038 2020 if( a == 0 ) fIntegralCerenkov[0] += b*log << 2039 if( a == 0 ) fIntegralCerenkov[0] += b*log(e0/x0); 2021 else fIntegralCerenkov[0] += y0*(e << 2040 else fIntegralCerenkov[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a; >> 2041 >> 2042 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = " >> 2043 // <<b<<"; result = "<<result<<G4endl; 2022 2044 2023 //G4cout<<" a="<< a <<" b="<< b <<" res << 2024 return result; 2045 return result; >> 2046 2025 } 2047 } 2026 2048 2027 ///////////////////////////////////////////// 2049 /////////////////////////////////////////////////////////////////////////////// 2028 // 2050 // 2029 // Integration of MM-Cerenkov cross-section f 2051 // Integration of MM-Cerenkov cross-section for the case of 2030 // passing across border between intervals 2052 // passing across border between intervals 2031 2053 2032 G4double G4PAIxSection::SumOverBordMM( G4int << 2054 G4double G4PAIxSection::SumOverBordMM( G4int i , >> 2055 G4double en0 ) 2033 { 2056 { 2034 G4double x0,x1,y0,yy1,a,b,e0,c,d,result; 2057 G4double x0,x1,y0,yy1,a,b,e0,c,d,result; 2035 2058 2036 e0 = en0; 2059 e0 = en0; 2037 x0 = fSplineEnergy[i]; 2060 x0 = fSplineEnergy[i]; 2038 x1 = fSplineEnergy[i+1]; 2061 x1 = fSplineEnergy[i+1]; 2039 y0 = fdNdxMM[i]; 2062 y0 = fdNdxMM[i]; 2040 yy1 = fdNdxMM[i+1]; 2063 yy1 = fdNdxMM[i+1]; 2041 2064 >> 2065 // G4cout<<G4endl; 2042 // G4cout<<"SumBordC, i = "<<i<<"; en0 = 2066 // G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1 2043 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G 2067 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl; 2044 c = x1/x0; 2068 c = x1/x0; 2045 d = e0/x0; 2069 d = e0/x0; 2046 a = log10(yy1/y0)/log10(c); 2070 a = log10(yy1/y0)/log10(c); 2047 << 2071 if(a > 10.0) return 0.; 2048 if(a > 20.0) b = 0.0; << 2072 // b0 = log10(y0) - a*log10(x0); 2049 else b = y0/pow(x0,a); << 2073 b = y0/pow(x0,a); // pow(10.,b0); 2050 2074 2051 a += 1.0; 2075 a += 1.0; 2052 if( a == 0 ) result = b*log(x0/e0); 2076 if( a == 0 ) result = b*log(x0/e0); 2053 else result = y0*(x0 - e0*pow(d,a- 2077 else result = y0*(x0 - e0*pow(d,a-1))/a; 2054 a += 1.0; 2078 a += 1.0; 2055 2079 2056 if( a == 0 ) fIntegralMM[0] += b*log(x0/e0 2080 if( a == 0 ) fIntegralMM[0] += b*log(x0/e0); 2057 else fIntegralMM[0] += y0*(x0*x0 - 2081 else fIntegralMM[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a; 2058 2082 2059 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b << 2083 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "<<b<<"; result = "<<result<<G4endl; 2060 2084 2061 x0 = fSplineEnergy[i - 1]; 2085 x0 = fSplineEnergy[i - 1]; 2062 x1 = fSplineEnergy[i - 2]; 2086 x1 = fSplineEnergy[i - 2]; 2063 y0 = fdNdxMM[i - 1]; 2087 y0 = fdNdxMM[i - 1]; 2064 yy1 = fdNdxMM[i - 2]; 2088 yy1 = fdNdxMM[i - 2]; 2065 2089 2066 // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1 2090 // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1 2067 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4 2091 // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl; 2068 2092 2069 c = x1/x0; 2093 c = x1/x0; 2070 d = e0/x0; 2094 d = e0/x0; 2071 a = log10(yy1/y0)/log10(x1/x0); 2095 a = log10(yy1/y0)/log10(x1/x0); 2072 << 2096 // b0 = log10(y0) - a*log10(x0); 2073 if(a > 20.0) b = 0.0; << 2097 b = y0/pow(x0,a); // pow(10.,b0); 2074 else b = y0/pow(x0,a); << 2075 2098 2076 a += 1.0; 2099 a += 1.0; 2077 if( a == 0 ) result += b*log(e0/x0); 2100 if( a == 0 ) result += b*log(e0/x0); 2078 else result += y0*(e0*pow(d,a-1) - 2101 else result += y0*(e0*pow(d,a-1) - x0 )/a; 2079 a += 1.0; 2102 a += 1.0; 2080 2103 2081 if( a == 0 ) fIntegralMM[0] += b*log(e0/x0 << 2104 if( a == 0 ) fIntegralMM[0] += b*log(e0/x0); 2082 else fIntegralMM[0] += y0*(e0*e0*p << 2105 else fIntegralMM[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a; 2083 2106 2084 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b 2107 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = " 2085 // <<b<<"; result = "<<result<<G4endl; 2108 // <<b<<"; result = "<<result<<G4endl; 2086 2109 2087 return result; 2110 return result; 2088 2111 2089 } 2112 } 2090 2113 2091 ///////////////////////////////////////////// 2114 /////////////////////////////////////////////////////////////////////////////// 2092 // 2115 // 2093 // Integration of Plasmon cross-section for t 2116 // Integration of Plasmon cross-section for the case of 2094 // passing across border between intervals 2117 // passing across border between intervals 2095 2118 2096 G4double G4PAIxSection::SumOverBordPlasmon( G 2119 G4double G4PAIxSection::SumOverBordPlasmon( G4int i , 2097 2120 G4double en0 ) 2098 { 2121 { 2099 G4double x0,x1,y0,yy1,a,b,c,d,e0,result; 2122 G4double x0,x1,y0,yy1,a,b,c,d,e0,result; 2100 2123 2101 e0 = en0; 2124 e0 = en0; 2102 x0 = fSplineEnergy[i]; 2125 x0 = fSplineEnergy[i]; 2103 x1 = fSplineEnergy[i+1]; 2126 x1 = fSplineEnergy[i+1]; 2104 y0 = fdNdxPlasmon[i]; 2127 y0 = fdNdxPlasmon[i]; 2105 yy1 = fdNdxPlasmon[i+1]; 2128 yy1 = fdNdxPlasmon[i+1]; 2106 2129 2107 c = x1/x0; 2130 c = x1/x0; 2108 d = e0/x0; 2131 d = e0/x0; 2109 a = log10(yy1/y0)/log10(c); 2132 a = log10(yy1/y0)/log10(c); 2110 << 2133 if(a > 10.0) return 0.; 2111 if(a > 20.0) b = 0.0; << 2134 // b0 = log10(y0) - a*log10(x0); 2112 else b = y0/pow(x0,a); << 2135 b = y0/pow(x0,a); //pow(10.,b); 2113 2136 2114 a += 1.0; 2137 a += 1.0; 2115 if( a == 0 ) result = b*log(x0/e0); 2138 if( a == 0 ) result = b*log(x0/e0); 2116 else result = y0*(x0 - e0*pow(d,a- 2139 else result = y0*(x0 - e0*pow(d,a-1))/a; 2117 a += 1.0; 2140 a += 1.0; 2118 2141 2119 if( a == 0 ) fIntegralPlasmon[0] += b*log( 2142 if( a == 0 ) fIntegralPlasmon[0] += b*log(x0/e0); 2120 else fIntegralPlasmon[0] += y0*(x0 2143 else fIntegralPlasmon[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a; 2121 2144 2122 x0 = fSplineEnergy[i - 1]; 2145 x0 = fSplineEnergy[i - 1]; 2123 x1 = fSplineEnergy[i - 2]; 2146 x1 = fSplineEnergy[i - 2]; 2124 y0 = fdNdxPlasmon[i - 1]; 2147 y0 = fdNdxPlasmon[i - 1]; 2125 yy1 = fdNdxPlasmon[i - 2]; 2148 yy1 = fdNdxPlasmon[i - 2]; 2126 2149 2127 c = x1/x0; 2150 c = x1/x0; 2128 d = e0/x0; 2151 d = e0/x0; 2129 a = log10(yy1/y0)/log10(c); 2152 a = log10(yy1/y0)/log10(c); 2130 << 2153 // b0 = log10(y0) - a*log10(x0); 2131 if(a > 20.0) b = 0.0; << 2154 b = y0/pow(x0,a);// pow(10.,b0); 2132 else b = y0/pow(x0,a); << 2133 2155 2134 a += 1.0; 2156 a += 1.0; 2135 if( a == 0 ) result += b*log(e0/x0); 2157 if( a == 0 ) result += b*log(e0/x0); 2136 else result += y0*(e0*pow(d,a-1) - 2158 else result += y0*(e0*pow(d,a-1) - x0)/a; 2137 a += 1.0; 2159 a += 1.0; 2138 2160 2139 if( a == 0 ) fIntegralPlasmon[0] += b*lo 2161 if( a == 0 ) fIntegralPlasmon[0] += b*log(e0/x0); 2140 else fIntegralPlasmon[0] += y0*( 2162 else fIntegralPlasmon[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a; 2141 2163 2142 return result; 2164 return result; >> 2165 2143 } 2166 } 2144 2167 2145 ///////////////////////////////////////////// 2168 /////////////////////////////////////////////////////////////////////////////// 2146 // 2169 // 2147 // Integration of resonance cross-section for 2170 // Integration of resonance cross-section for the case of 2148 // passing across border between intervals 2171 // passing across border between intervals 2149 2172 2150 G4double G4PAIxSection::SumOverBordResonance( 2173 G4double G4PAIxSection::SumOverBordResonance( G4int i , 2151 2174 G4double en0 ) 2152 { 2175 { 2153 G4double x0,x1,y0,yy1,a,b,c,d,e0,result; 2176 G4double x0,x1,y0,yy1,a,b,c,d,e0,result; 2154 2177 2155 e0 = en0; 2178 e0 = en0; 2156 x0 = fSplineEnergy[i]; 2179 x0 = fSplineEnergy[i]; 2157 x1 = fSplineEnergy[i+1]; 2180 x1 = fSplineEnergy[i+1]; 2158 y0 = fdNdxResonance[i]; 2181 y0 = fdNdxResonance[i]; 2159 yy1 = fdNdxResonance[i+1]; 2182 yy1 = fdNdxResonance[i+1]; 2160 2183 2161 c = x1/x0; 2184 c = x1/x0; 2162 d = e0/x0; 2185 d = e0/x0; 2163 a = log10(yy1/y0)/log10(c); 2186 a = log10(yy1/y0)/log10(c); 2164 << 2187 if(a > 10.0) return 0.; 2165 if(a > 20.0) b = 0.0; << 2188 // b0 = log10(y0) - a*log10(x0); 2166 else b = y0/pow(x0,a); << 2189 b = y0/pow(x0,a); //pow(10.,b); 2167 2190 2168 a += 1.0; 2191 a += 1.0; 2169 if( a == 0 ) result = b*log(x0/e0); 2192 if( a == 0 ) result = b*log(x0/e0); 2170 else result = y0*(x0 - e0*pow(d,a- 2193 else result = y0*(x0 - e0*pow(d,a-1))/a; 2171 a += 1.0; 2194 a += 1.0; 2172 2195 2173 if( a == 0 ) fIntegralResonance[0] += b*lo 2196 if( a == 0 ) fIntegralResonance[0] += b*log(x0/e0); 2174 else fIntegralResonance[0] += y0*( 2197 else fIntegralResonance[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a; 2175 2198 2176 x0 = fSplineEnergy[i - 1]; 2199 x0 = fSplineEnergy[i - 1]; 2177 x1 = fSplineEnergy[i - 2]; 2200 x1 = fSplineEnergy[i - 2]; 2178 y0 = fdNdxResonance[i - 1]; 2201 y0 = fdNdxResonance[i - 1]; 2179 yy1 = fdNdxResonance[i - 2]; 2202 yy1 = fdNdxResonance[i - 2]; 2180 2203 2181 c = x1/x0; 2204 c = x1/x0; 2182 d = e0/x0; 2205 d = e0/x0; 2183 a = log10(yy1/y0)/log10(c); 2206 a = log10(yy1/y0)/log10(c); 2184 << 2207 // b0 = log10(y0) - a*log10(x0); 2185 if(a > 20.0) b = 0.0; << 2208 b = y0/pow(x0,a);// pow(10.,b0); 2186 else b = y0/pow(x0,a); << 2187 2209 2188 a += 1.0; 2210 a += 1.0; 2189 if( a == 0 ) result += b*log(e0/x0); 2211 if( a == 0 ) result += b*log(e0/x0); 2190 else result += y0*(e0*pow(d,a-1) - 2212 else result += y0*(e0*pow(d,a-1) - x0)/a; 2191 a += 1.0; 2213 a += 1.0; 2192 2214 2193 if( a == 0 ) fIntegralResonance[0] += b*lo << 2215 if( a == 0 ) fIntegralResonance[0] += b*log(e0/x0); 2194 else fIntegralResonance[0] += y0*( << 2216 else fIntegralResonance[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a; 2195 2217 2196 return result; 2218 return result; 2197 2219 2198 } 2220 } 2199 2221 2200 ///////////////////////////////////////////// 2222 ///////////////////////////////////////////////////////////////////////// 2201 // 2223 // 2202 // Returns random PAI-total energy loss over 2224 // Returns random PAI-total energy loss over step 2203 2225 2204 G4double G4PAIxSection::GetStepEnergyLoss( G4 2226 G4double G4PAIxSection::GetStepEnergyLoss( G4double step ) 2205 { 2227 { 2206 G4long numOfCollisions; 2228 G4long numOfCollisions; 2207 G4double meanNumber, loss = 0.0; 2229 G4double meanNumber, loss = 0.0; 2208 2230 2209 // G4cout<<" G4PAIxSection::GetStepEnergyLo 2231 // G4cout<<" G4PAIxSection::GetStepEnergyLoss "<<G4endl; 2210 2232 2211 meanNumber = fIntegralPAIxSection[1]*step; 2233 meanNumber = fIntegralPAIxSection[1]*step; 2212 numOfCollisions = G4Poisson(meanNumber); 2234 numOfCollisions = G4Poisson(meanNumber); 2213 2235 2214 // G4cout<<"numOfCollisions = "<<numOfColli << 2236 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl; 2215 2237 2216 while(numOfCollisions) 2238 while(numOfCollisions) 2217 { 2239 { 2218 loss += GetEnergyTransfer(); << 2240 loss += GetEnergyTransfer(); 2219 numOfCollisions--; << 2241 numOfCollisions--; 2220 // Loop checking, 03-Aug-2015, Vladimir I << 2221 } 2242 } 2222 // G4cout<<"PAI energy loss = "<<loss/keV<< 2243 // G4cout<<"PAI energy loss = "<<loss/keV<<" keV"<<G4endl; 2223 2244 2224 return loss; 2245 return loss; 2225 } 2246 } 2226 2247 2227 ///////////////////////////////////////////// 2248 ///////////////////////////////////////////////////////////////////////// 2228 // 2249 // 2229 // Returns random PAI-total energy transfer i 2250 // Returns random PAI-total energy transfer in one collision 2230 2251 2231 G4double G4PAIxSection::GetEnergyTransfer() 2252 G4double G4PAIxSection::GetEnergyTransfer() 2232 { 2253 { 2233 G4int iTransfer ; 2254 G4int iTransfer ; 2234 2255 2235 G4double energyTransfer, position; 2256 G4double energyTransfer, position; 2236 2257 2237 position = fIntegralPAIxSection[1]*G4Unifor 2258 position = fIntegralPAIxSection[1]*G4UniformRand(); 2238 2259 2239 for( iTransfer = 1; iTransfer <= fSplineNum 2260 for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ ) 2240 { 2261 { 2241 if( position >= fIntegralPAIxSection[iTra << 2262 if( position >= fIntegralPAIxSection[iTransfer] ) break; 2242 } 2263 } 2243 if(iTransfer > fSplineNumber) iTransfer--; 2264 if(iTransfer > fSplineNumber) iTransfer--; 2244 2265 2245 energyTransfer = fSplineEnergy[iTransfer]; 2266 energyTransfer = fSplineEnergy[iTransfer]; 2246 2267 2247 if(iTransfer > 1) 2268 if(iTransfer > 1) 2248 { 2269 { 2249 energyTransfer -= (fSplineEnergy[iTransfe 2270 energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand(); 2250 } 2271 } 2251 return energyTransfer; 2272 return energyTransfer; 2252 } 2273 } 2253 2274 2254 ///////////////////////////////////////////// 2275 ///////////////////////////////////////////////////////////////////////// 2255 // 2276 // 2256 // Returns random Cerenkov energy loss over s 2277 // Returns random Cerenkov energy loss over step 2257 2278 2258 G4double G4PAIxSection::GetStepCerenkovLoss( 2279 G4double G4PAIxSection::GetStepCerenkovLoss( G4double step ) 2259 { 2280 { 2260 G4long numOfCollisions; 2281 G4long numOfCollisions; 2261 G4double meanNumber, loss = 0.0; 2282 G4double meanNumber, loss = 0.0; 2262 2283 2263 // G4cout<<" G4PAIxSection::GetStepCerenkov 2284 // G4cout<<" G4PAIxSection::GetStepCerenkovLoss "<<G4endl; 2264 2285 2265 meanNumber = fIntegralCerenkov[1]*step; 2286 meanNumber = fIntegralCerenkov[1]*step; 2266 numOfCollisions = G4Poisson(meanNumber); 2287 numOfCollisions = G4Poisson(meanNumber); 2267 2288 2268 // G4cout<<"numOfCollisions = "<<numOfCol 2289 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl; 2269 2290 2270 while(numOfCollisions) 2291 while(numOfCollisions) 2271 { 2292 { 2272 loss += GetCerenkovEnergyTransfer(); 2293 loss += GetCerenkovEnergyTransfer(); 2273 numOfCollisions--; 2294 numOfCollisions--; 2274 // Loop checking, 03-Aug-2015, Vladimir I << 2275 } 2295 } 2276 // G4cout<<"PAI Cerenkov loss = "<<loss/keV 2296 // G4cout<<"PAI Cerenkov loss = "<<loss/keV<<" keV"<<G4endl; 2277 2297 2278 return loss; 2298 return loss; 2279 } 2299 } 2280 2300 2281 ///////////////////////////////////////////// 2301 ///////////////////////////////////////////////////////////////////////// 2282 // 2302 // 2283 // Returns random MM-Cerenkov energy loss ove 2303 // Returns random MM-Cerenkov energy loss over step 2284 2304 2285 G4double G4PAIxSection::GetStepMMLoss( G4doub 2305 G4double G4PAIxSection::GetStepMMLoss( G4double step ) 2286 { 2306 { 2287 G4long numOfCollisions; 2307 G4long numOfCollisions; 2288 G4double meanNumber, loss = 0.0; 2308 G4double meanNumber, loss = 0.0; 2289 2309 2290 // G4cout<<" G4PAIxSection::GetStepMMLoss " 2310 // G4cout<<" G4PAIxSection::GetStepMMLoss "<<G4endl; 2291 2311 2292 meanNumber = fIntegralMM[1]*step; 2312 meanNumber = fIntegralMM[1]*step; 2293 numOfCollisions = G4Poisson(meanNumber); 2313 numOfCollisions = G4Poisson(meanNumber); 2294 2314 2295 // G4cout<<"numOfCollisions = "<<numOfCol 2315 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl; 2296 2316 2297 while(numOfCollisions) 2317 while(numOfCollisions) 2298 { 2318 { 2299 loss += GetMMEnergyTransfer(); 2319 loss += GetMMEnergyTransfer(); 2300 numOfCollisions--; 2320 numOfCollisions--; 2301 // Loop checking, 03-Aug-2015, Vladimir I << 2302 } 2321 } 2303 // G4cout<<"PAI MM-Cerenkov loss = "<<loss/ 2322 // G4cout<<"PAI MM-Cerenkov loss = "<<loss/keV<<" keV"<<G4endl; 2304 2323 2305 return loss; 2324 return loss; 2306 } 2325 } 2307 2326 2308 ///////////////////////////////////////////// 2327 ///////////////////////////////////////////////////////////////////////// 2309 // 2328 // 2310 // Returns Cerenkov energy transfer in one co 2329 // Returns Cerenkov energy transfer in one collision 2311 2330 2312 G4double G4PAIxSection::GetCerenkovEnergyTran 2331 G4double G4PAIxSection::GetCerenkovEnergyTransfer() 2313 { 2332 { 2314 G4int iTransfer ; 2333 G4int iTransfer ; 2315 2334 2316 G4double energyTransfer, position; 2335 G4double energyTransfer, position; 2317 2336 2318 position = fIntegralCerenkov[1]*G4UniformRa 2337 position = fIntegralCerenkov[1]*G4UniformRand(); 2319 2338 2320 for( iTransfer = 1; iTransfer <= fSplineNum 2339 for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ ) 2321 { 2340 { 2322 if( position >= fIntegralCerenkov[iTr 2341 if( position >= fIntegralCerenkov[iTransfer] ) break; 2323 } 2342 } 2324 if(iTransfer > fSplineNumber) iTransfer--; 2343 if(iTransfer > fSplineNumber) iTransfer--; 2325 2344 2326 energyTransfer = fSplineEnergy[iTransfer]; 2345 energyTransfer = fSplineEnergy[iTransfer]; 2327 2346 2328 if(iTransfer > 1) 2347 if(iTransfer > 1) 2329 { 2348 { 2330 energyTransfer -= (fSplineEnergy[iTransfe 2349 energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand(); 2331 } 2350 } 2332 return energyTransfer; 2351 return energyTransfer; 2333 } 2352 } 2334 2353 2335 ///////////////////////////////////////////// 2354 ///////////////////////////////////////////////////////////////////////// 2336 // 2355 // 2337 // Returns MM-Cerenkov energy transfer in one 2356 // Returns MM-Cerenkov energy transfer in one collision 2338 2357 2339 G4double G4PAIxSection::GetMMEnergyTransfer() 2358 G4double G4PAIxSection::GetMMEnergyTransfer() 2340 { 2359 { 2341 G4int iTransfer ; 2360 G4int iTransfer ; 2342 2361 2343 G4double energyTransfer, position; 2362 G4double energyTransfer, position; 2344 2363 2345 position = fIntegralMM[1]*G4UniformRand(); 2364 position = fIntegralMM[1]*G4UniformRand(); 2346 2365 2347 for( iTransfer = 1; iTransfer <= fSplineNum 2366 for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ ) 2348 { 2367 { 2349 if( position >= fIntegralMM[iTransfer] ) << 2368 if( position >= fIntegralMM[iTransfer] ) break; 2350 } 2369 } 2351 if(iTransfer > fSplineNumber) iTransfer--; 2370 if(iTransfer > fSplineNumber) iTransfer--; 2352 2371 2353 energyTransfer = fSplineEnergy[iTransfer]; 2372 energyTransfer = fSplineEnergy[iTransfer]; 2354 2373 2355 if(iTransfer > 1) 2374 if(iTransfer > 1) 2356 { 2375 { 2357 energyTransfer -= (fSplineEnergy[iTransfe 2376 energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand(); 2358 } 2377 } 2359 return energyTransfer; 2378 return energyTransfer; 2360 } 2379 } 2361 2380 2362 ///////////////////////////////////////////// 2381 ///////////////////////////////////////////////////////////////////////// 2363 // 2382 // 2364 // Returns random plasmon energy loss over st 2383 // Returns random plasmon energy loss over step 2365 2384 2366 G4double G4PAIxSection::GetStepPlasmonLoss( G 2385 G4double G4PAIxSection::GetStepPlasmonLoss( G4double step ) 2367 { 2386 { 2368 G4long numOfCollisions; 2387 G4long numOfCollisions; 2369 G4double meanNumber, loss = 0.0; 2388 G4double meanNumber, loss = 0.0; 2370 2389 2371 // G4cout<<" G4PAIxSection::GetStepPlasmonL 2390 // G4cout<<" G4PAIxSection::GetStepPlasmonLoss "<<G4endl; 2372 2391 2373 meanNumber = fIntegralPlasmon[1]*step; 2392 meanNumber = fIntegralPlasmon[1]*step; 2374 numOfCollisions = G4Poisson(meanNumber); 2393 numOfCollisions = G4Poisson(meanNumber); 2375 2394 2376 // G4cout<<"numOfCollisions = "<<numOfCol 2395 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl; 2377 2396 2378 while(numOfCollisions) 2397 while(numOfCollisions) 2379 { 2398 { 2380 loss += GetPlasmonEnergyTransfer(); 2399 loss += GetPlasmonEnergyTransfer(); 2381 numOfCollisions--; 2400 numOfCollisions--; 2382 // Loop checking, 03-Aug-2015, Vladimir I << 2383 } 2401 } 2384 // G4cout<<"PAI Plasmon loss = "<<loss/keV< 2402 // G4cout<<"PAI Plasmon loss = "<<loss/keV<<" keV"<<G4endl; 2385 2403 2386 return loss; 2404 return loss; 2387 } 2405 } 2388 2406 2389 ///////////////////////////////////////////// 2407 ///////////////////////////////////////////////////////////////////////// 2390 // 2408 // 2391 // Returns plasmon energy transfer in one col 2409 // Returns plasmon energy transfer in one collision 2392 2410 2393 G4double G4PAIxSection::GetPlasmonEnergyTrans 2411 G4double G4PAIxSection::GetPlasmonEnergyTransfer() 2394 { 2412 { 2395 G4int iTransfer ; 2413 G4int iTransfer ; 2396 2414 2397 G4double energyTransfer, position; 2415 G4double energyTransfer, position; 2398 2416 2399 position = fIntegralPlasmon[1]*G4UniformRan 2417 position = fIntegralPlasmon[1]*G4UniformRand(); 2400 2418 2401 for( iTransfer = 1; iTransfer <= fSplineNum 2419 for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ ) 2402 { 2420 { 2403 if( position >= fIntegralPlasmon[iTransfe << 2421 if( position >= fIntegralPlasmon[iTransfer] ) break; 2404 } 2422 } 2405 if(iTransfer > fSplineNumber) iTransfer--; 2423 if(iTransfer > fSplineNumber) iTransfer--; 2406 2424 2407 energyTransfer = fSplineEnergy[iTransfer]; 2425 energyTransfer = fSplineEnergy[iTransfer]; 2408 2426 2409 if(iTransfer > 1) 2427 if(iTransfer > 1) 2410 { 2428 { 2411 energyTransfer -= (fSplineEnergy[iTransfe 2429 energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand(); 2412 } 2430 } 2413 return energyTransfer; 2431 return energyTransfer; 2414 } 2432 } 2415 2433 2416 ///////////////////////////////////////////// 2434 ///////////////////////////////////////////////////////////////////////// 2417 // 2435 // 2418 // Returns random resonance energy loss over 2436 // Returns random resonance energy loss over step 2419 2437 2420 G4double G4PAIxSection::GetStepResonanceLoss( 2438 G4double G4PAIxSection::GetStepResonanceLoss( G4double step ) 2421 { 2439 { 2422 G4long numOfCollisions; 2440 G4long numOfCollisions; 2423 G4double meanNumber, loss = 0.0; 2441 G4double meanNumber, loss = 0.0; 2424 2442 2425 // G4cout<<" G4PAIxSection::GetStepCreLosnk 2443 // G4cout<<" G4PAIxSection::GetStepCreLosnkovs "<<G4endl; 2426 2444 2427 meanNumber = fIntegralResonance[1]*step; 2445 meanNumber = fIntegralResonance[1]*step; 2428 numOfCollisions = G4Poisson(meanNumber); 2446 numOfCollisions = G4Poisson(meanNumber); 2429 2447 2430 // G4cout<<"numOfCollisions = "<<numOfCol 2448 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl; 2431 2449 2432 while(numOfCollisions) 2450 while(numOfCollisions) 2433 { 2451 { 2434 loss += GetResonanceEnergyTransfer(); 2452 loss += GetResonanceEnergyTransfer(); 2435 numOfCollisions--; 2453 numOfCollisions--; 2436 // Loop checking, 03-Aug-2015, Vladimir I << 2437 } 2454 } 2438 // G4cout<<"PAI resonance loss = "<<loss/ke 2455 // G4cout<<"PAI resonance loss = "<<loss/keV<<" keV"<<G4endl; 2439 2456 2440 return loss; 2457 return loss; 2441 } 2458 } 2442 2459 2443 2460 2444 ///////////////////////////////////////////// 2461 ///////////////////////////////////////////////////////////////////////// 2445 // 2462 // 2446 // Returns resonance energy transfer in one c 2463 // Returns resonance energy transfer in one collision 2447 2464 2448 G4double G4PAIxSection::GetResonanceEnergyTra 2465 G4double G4PAIxSection::GetResonanceEnergyTransfer() 2449 { 2466 { 2450 G4int iTransfer ; 2467 G4int iTransfer ; 2451 2468 2452 G4double energyTransfer, position; 2469 G4double energyTransfer, position; 2453 2470 2454 position = fIntegralResonance[1]*G4UniformR 2471 position = fIntegralResonance[1]*G4UniformRand(); 2455 2472 2456 for( iTransfer = 1; iTransfer <= fSplineNum 2473 for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ ) 2457 { 2474 { 2458 if( position >= fIntegralResonance[iTrans << 2475 if( position >= fIntegralResonance[iTransfer] ) break; 2459 } 2476 } 2460 if(iTransfer > fSplineNumber) iTransfer--; 2477 if(iTransfer > fSplineNumber) iTransfer--; 2461 2478 2462 energyTransfer = fSplineEnergy[iTransfer]; 2479 energyTransfer = fSplineEnergy[iTransfer]; 2463 2480 2464 if(iTransfer > 1) 2481 if(iTransfer > 1) 2465 { 2482 { 2466 energyTransfer -= (fSplineEnergy[iTransfe 2483 energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand(); 2467 } 2484 } 2468 return energyTransfer; 2485 return energyTransfer; 2469 } 2486 } 2470 2487 2471 2488 2472 ///////////////////////////////////////////// 2489 ///////////////////////////////////////////////////////////////////////// 2473 // 2490 // 2474 // Returns Rutherford energy transfer in one 2491 // Returns Rutherford energy transfer in one collision 2475 2492 2476 G4double G4PAIxSection::GetRutherfordEnergyTr 2493 G4double G4PAIxSection::GetRutherfordEnergyTransfer() 2477 { 2494 { 2478 G4int iTransfer ; 2495 G4int iTransfer ; 2479 2496 2480 G4double energyTransfer, position; 2497 G4double energyTransfer, position; 2481 2498 2482 position = (fIntegralPlasmon[1]-fIntegralRe 2499 position = (fIntegralPlasmon[1]-fIntegralResonance[1])*G4UniformRand(); 2483 2500 2484 for( iTransfer = 1; iTransfer <= fSplineNum 2501 for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ ) 2485 { 2502 { 2486 if( position >= (fIntegralPlasmon[iTransf << 2503 if( position >= (fIntegralPlasmon[iTransfer]-fIntegralResonance[iTransfer]) ) break; 2487 } 2504 } 2488 if(iTransfer > fSplineNumber) iTransfer--; 2505 if(iTransfer > fSplineNumber) iTransfer--; 2489 2506 2490 energyTransfer = fSplineEnergy[iTransfer]; 2507 energyTransfer = fSplineEnergy[iTransfer]; 2491 2508 2492 if(iTransfer > 1) 2509 if(iTransfer > 1) 2493 { 2510 { 2494 energyTransfer -= (fSplineEnergy[iTransfe 2511 energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand(); 2495 } 2512 } 2496 return energyTransfer; 2513 return energyTransfer; 2497 } 2514 } 2498 2515 2499 ///////////////////////////////////////////// 2516 ///////////////////////////////////////////////////////////////////////////// 2500 // 2517 // 2501 2518 2502 void G4PAIxSection::CallError(G4int i, const 2519 void G4PAIxSection::CallError(G4int i, const G4String& methodName) const 2503 { 2520 { 2504 G4String head = "G4PAIxSection::" + methodN 2521 G4String head = "G4PAIxSection::" + methodName + "()"; 2505 G4ExceptionDescription ed; 2522 G4ExceptionDescription ed; 2506 ed << "Wrong index " << i << " fSplineNumbe 2523 ed << "Wrong index " << i << " fSplineNumber= " << fSplineNumber; 2507 G4Exception(head,"pai001",FatalException,ed 2524 G4Exception(head,"pai001",FatalException,ed); 2508 } 2525 } 2509 2526 2510 ///////////////////////////////////////////// 2527 ///////////////////////////////////////////////////////////////////////////// 2511 // 2528 // 2512 // Init array of Lorentz factors 2529 // Init array of Lorentz factors 2513 // 2530 // 2514 2531 2515 G4int G4PAIxSection::fNumberOfGammas = 111; 2532 G4int G4PAIxSection::fNumberOfGammas = 111; 2516 2533 2517 const G4double G4PAIxSection::fLorentzFactor[ 2534 const G4double G4PAIxSection::fLorentzFactor[112] = // fNumberOfGammas+1 2518 { 2535 { 2519 0.0, 2536 0.0, 2520 1.094989e+00, 1.107813e+00, 1.122369e+00, 1.1 2537 1.094989e+00, 1.107813e+00, 1.122369e+00, 1.138890e+00, 1.157642e+00, 2521 1.178925e+00, 1.203082e+00, 1.230500e+00, 1.2 2538 1.178925e+00, 1.203082e+00, 1.230500e+00, 1.261620e+00, 1.296942e+00, // 10 2522 1.337032e+00, 1.382535e+00, 1.434181e+00, 1.4 2539 1.337032e+00, 1.382535e+00, 1.434181e+00, 1.492800e+00, 1.559334e+00, 2523 1.634850e+00, 1.720562e+00, 1.817845e+00, 1.9 2540 1.634850e+00, 1.720562e+00, 1.817845e+00, 1.928263e+00, 2.053589e+00, // 20 2524 2.195835e+00, 2.357285e+00, 2.540533e+00, 2.7 2541 2.195835e+00, 2.357285e+00, 2.540533e+00, 2.748522e+00, 2.984591e+00, 2525 3.252533e+00, 3.556649e+00, 3.901824e+00, 4.2 2542 3.252533e+00, 3.556649e+00, 3.901824e+00, 4.293602e+00, 4.738274e+00, // 30 2526 5.242981e+00, 5.815829e+00, 6.466019e+00, 7.2 2543 5.242981e+00, 5.815829e+00, 6.466019e+00, 7.203990e+00, 8.041596e+00, 2527 8.992288e+00, 1.007133e+01, 1.129606e+01, 1.2 2544 8.992288e+00, 1.007133e+01, 1.129606e+01, 1.268614e+01, 1.426390e+01, // 40 2528 1.605467e+01, 1.808721e+01, 2.039417e+01, 2.3 2545 1.605467e+01, 1.808721e+01, 2.039417e+01, 2.301259e+01, 2.598453e+01, 2529 2.935771e+01, 3.318630e+01, 3.753180e+01, 4.2 2546 2.935771e+01, 3.318630e+01, 3.753180e+01, 4.246399e+01, 4.806208e+01, // 50 2530 5.441597e+01, 6.162770e+01, 6.981310e+01, 7.9 2547 5.441597e+01, 6.162770e+01, 6.981310e+01, 7.910361e+01, 8.964844e+01, 2531 1.016169e+02, 1.152013e+02, 1.306197e+02, 1.4 2548 1.016169e+02, 1.152013e+02, 1.306197e+02, 1.481198e+02, 1.679826e+02, // 60 2532 1.905270e+02, 2.161152e+02, 2.451581e+02, 2.7 2549 1.905270e+02, 2.161152e+02, 2.451581e+02, 2.781221e+02, 3.155365e+02, 2533 3.580024e+02, 4.062016e+02, 4.609081e+02, 5.2 2550 3.580024e+02, 4.062016e+02, 4.609081e+02, 5.230007e+02, 5.934765e+02, // 70 2534 6.734672e+02, 7.642575e+02, 8.673056e+02, 9.8 2551 6.734672e+02, 7.642575e+02, 8.673056e+02, 9.842662e+02, 1.117018e+03, 2535 1.267692e+03, 1.438709e+03, 1.632816e+03, 1.8 2552 1.267692e+03, 1.438709e+03, 1.632816e+03, 1.853128e+03, 2.103186e+03, // 80 2536 2.387004e+03, 2.709140e+03, 3.074768e+03, 3.4 2553 2.387004e+03, 2.709140e+03, 3.074768e+03, 3.489760e+03, 3.960780e+03, 2537 4.495394e+03, 5.102185e+03, 5.790900e+03, 6.5 2554 4.495394e+03, 5.102185e+03, 5.790900e+03, 6.572600e+03, 7.459837e+03, // 90 2538 8.466860e+03, 9.609843e+03, 1.090714e+04, 1.2 2555 8.466860e+03, 9.609843e+03, 1.090714e+04, 1.237959e+04, 1.405083e+04, 2539 1.594771e+04, 1.810069e+04, 2.054434e+04, 2.3 2556 1.594771e+04, 1.810069e+04, 2.054434e+04, 2.331792e+04, 2.646595e+04, // 100 2540 3.003901e+04, 3.409446e+04, 3.869745e+04, 4.3 2557 3.003901e+04, 3.409446e+04, 3.869745e+04, 4.392189e+04, 4.985168e+04, 2541 5.658206e+04, 6.422112e+04, 7.289153e+04, 8.2 2558 5.658206e+04, 6.422112e+04, 7.289153e+04, 8.273254e+04, 9.390219e+04, // 110 2542 1.065799e+05 2559 1.065799e+05 2543 }; 2560 }; 2544 2561 2545 ///////////////////////////////////////////// 2562 /////////////////////////////////////////////////////////////////////// 2546 // 2563 // 2547 // The number of gamma for creation of splin 2564 // The number of gamma for creation of spline (near ion-min , G ~ 4 ) 2548 // 2565 // 2549 2566 2550 const 2567 const 2551 G4int G4PAIxSection::fRefGammaNumber = 29; 2568 G4int G4PAIxSection::fRefGammaNumber = 29; 2552 2569 2553 2570 2554 // 2571 // 2555 // end of G4PAIxSection implementation file 2572 // end of G4PAIxSection implementation file 2556 // 2573 // 2557 ///////////////////////////////////////////// 2574 //////////////////////////////////////////////////////////////////////////// 2558 2575 2559 2576