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