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 // 27 // 28 // G4InitXscPAI.cc -- class implementation fil 28 // G4InitXscPAI.cc -- class implementation file 29 // 29 // 30 // GEANT 4 class implementation file 30 // GEANT 4 class implementation file 31 // 31 // 32 // For information related to this code, pleas 32 // For information related to this code, please, contact 33 // the Geant4 Collaboration. 33 // the Geant4 Collaboration. 34 // 34 // 35 // R&D: Vladimir.Grichine@cern.ch 35 // R&D: Vladimir.Grichine@cern.ch 36 // 36 // 37 // History: 37 // History: 38 // 38 // 39 39 40 40 41 41 42 #include "G4InitXscPAI.hh" 42 #include "G4InitXscPAI.hh" 43 43 44 #include "globals.hh" 44 #include "globals.hh" 45 #include "G4PhysicalConstants.hh" 45 #include "G4PhysicalConstants.hh" 46 #include "G4SystemOfUnits.hh" 46 #include "G4SystemOfUnits.hh" 47 #include "G4ios.hh" 47 #include "G4ios.hh" 48 #include "G4Poisson.hh" 48 #include "G4Poisson.hh" 49 #include "G4Integrator.hh" 49 #include "G4Integrator.hh" 50 #include "G4Material.hh" 50 #include "G4Material.hh" 51 #include "G4MaterialCutsCouple.hh" 51 #include "G4MaterialCutsCouple.hh" 52 #include "G4SandiaTable.hh" 52 #include "G4SandiaTable.hh" 53 53 54 54 55 55 56 // Local class constants 56 // Local class constants 57 57 58 const G4double G4InitXscPAI::fDelta = 0 58 const G4double G4InitXscPAI::fDelta = 0.005 ; // energy shift from interval border 59 const G4int G4InitXscPAI::fPAIbin = 1 59 const G4int G4InitXscPAI::fPAIbin = 100 ; // size of energy transfer vectors 60 const G4double G4InitXscPAI::fSolidDensity = 0 60 const G4double G4InitXscPAI::fSolidDensity = 0.05*g/cm3 ; // ~gas-solid border 61 61 62 ////////////////////////////////////////////// 62 ////////////////////////////////////////////////////////////////// 63 // 63 // 64 // Constructor 64 // Constructor 65 // 65 // 66 66 67 using namespace std; 67 using namespace std; 68 68 69 G4InitXscPAI::G4InitXscPAI( const G4MaterialCu 69 G4InitXscPAI::G4InitXscPAI( const G4MaterialCutsCouple* matCC) 70 : fPAIxscVector(nullptr), 70 : fPAIxscVector(nullptr), 71 fPAIdEdxVector(nullptr), 71 fPAIdEdxVector(nullptr), 72 fPAIphotonVector(nullptr), 72 fPAIphotonVector(nullptr), 73 fPAIelectronVector(nullptr), 73 fPAIelectronVector(nullptr), 74 fChCosSqVector(nullptr), 74 fChCosSqVector(nullptr), 75 fChWidthVector(nullptr) 75 fChWidthVector(nullptr) 76 { 76 { 77 G4int i, j, matIndex; 77 G4int i, j, matIndex; 78 78 79 fDensity = matCC->GetMaterial()->Get 79 fDensity = matCC->GetMaterial()->GetDensity(); 80 fElectronDensity = matCC->GetMaterial()->Get 80 fElectronDensity = matCC->GetMaterial()->GetElectronDensity(); 81 matIndex = (G4int)matCC->GetMaterial << 81 matIndex = matCC->GetMaterial()->GetIndex(); 82 82 83 fSandia = new G4SandiaTable(matInde 83 fSandia = new G4SandiaTable(matIndex); 84 fIntervalNumber = fSandia->GetMaxInterval() 84 fIntervalNumber = fSandia->GetMaxInterval()-1; 85 85 86 fMatSandiaMatrix = new G4OrderedTable(); 86 fMatSandiaMatrix = new G4OrderedTable(); 87 87 88 for (i = 0; i < fIntervalNumber; ++i) << 88 for (i = 0; i < fIntervalNumber; i++) 89 { 89 { 90 fMatSandiaMatrix->push_back(new G4DataVect 90 fMatSandiaMatrix->push_back(new G4DataVector(5,0.)); 91 } 91 } 92 for (i = 0; i < fIntervalNumber; ++i) << 92 for (i = 0; i < fIntervalNumber; i++) 93 { 93 { 94 (*(*fMatSandiaMatrix)[i])[0] = fSandia->Ge 94 (*(*fMatSandiaMatrix)[i])[0] = fSandia->GetSandiaMatTable(i,0); 95 95 96 for(j = 1; j < 5 ; ++j) << 96 for(j = 1; j < 5 ; j++) 97 { 97 { 98 (*(*fMatSandiaMatrix)[i])[j] = fSandia-> 98 (*(*fMatSandiaMatrix)[i])[j] = fSandia->GetSandiaMatTable(i,j)*fDensity; 99 } 99 } 100 } 100 } 101 KillCloseIntervals(); 101 KillCloseIntervals(); 102 Normalisation(); 102 Normalisation(); 103 fBetaGammaSq = fTmax = 0.0; 103 fBetaGammaSq = fTmax = 0.0; 104 fIntervalTmax = fCurrentInterval = 0; 104 fIntervalTmax = fCurrentInterval = 0; 105 } 105 } 106 106 107 107 108 108 109 109 110 ////////////////////////////////////////////// 110 //////////////////////////////////////////////////////////////////////////// 111 // 111 // 112 // Destructor 112 // Destructor 113 113 114 G4InitXscPAI::~G4InitXscPAI() 114 G4InitXscPAI::~G4InitXscPAI() 115 { 115 { 116 delete fPAIxscVector; << 116 if(fPAIxscVector) delete fPAIxscVector; 117 delete fPAIdEdxVector; << 117 if(fPAIdEdxVector) delete fPAIdEdxVector; 118 delete fPAIphotonVector; << 118 if(fPAIphotonVector) delete fPAIphotonVector; 119 delete fPAIelectronVector; << 119 if(fPAIelectronVector) delete fPAIelectronVector; 120 delete fChCosSqVector; << 120 if(fChCosSqVector) delete fChCosSqVector; 121 delete fChWidthVector; << 121 if(fChWidthVector) delete fChWidthVector; 122 delete fSandia; 122 delete fSandia; 123 delete fMatSandiaMatrix; 123 delete fMatSandiaMatrix; 124 } 124 } 125 125 126 ////////////////////////////////////////////// 126 //////////////////////////////////////////////////////////////////////// 127 // 127 // 128 // Kill close intervals, recalculate fInterval 128 // Kill close intervals, recalculate fIntervalNumber 129 129 130 void G4InitXscPAI::KillCloseIntervals() 130 void G4InitXscPAI::KillCloseIntervals() 131 { 131 { 132 G4int i, j, k; 132 G4int i, j, k; 133 G4double energy1, energy2; 133 G4double energy1, energy2; 134 134 135 for( i = 0 ; i < fIntervalNumber - 1 ; i++ ) 135 for( i = 0 ; i < fIntervalNumber - 1 ; i++ ) 136 { 136 { 137 energy1 = (*(*fMatSandiaMatrix)[i])[0]; 137 energy1 = (*(*fMatSandiaMatrix)[i])[0]; 138 energy2 = (*(*fMatSandiaMatrix)[i+1])[0]; 138 energy2 = (*(*fMatSandiaMatrix)[i+1])[0]; 139 139 140 if( energy2 - energy1 > 1.5*fDelta*(energy 140 if( energy2 - energy1 > 1.5*fDelta*(energy1 + energy2) ) continue ; 141 else 141 else 142 { 142 { 143 for(j = i; j < fIntervalNumber-1; j++) 143 for(j = i; j < fIntervalNumber-1; j++) 144 { 144 { 145 for( k = 0; k < 5; k++ ) 145 for( k = 0; k < 5; k++ ) 146 { 146 { 147 (*(*fMatSandiaMatrix)[j])[k] = (*(*f 147 (*(*fMatSandiaMatrix)[j])[k] = (*(*fMatSandiaMatrix)[j+1])[k]; 148 } 148 } 149 } 149 } 150 fIntervalNumber-- ; 150 fIntervalNumber-- ; 151 i-- ; 151 i-- ; 152 } 152 } 153 } 153 } 154 154 155 } 155 } 156 156 157 ////////////////////////////////////////////// 157 //////////////////////////////////////////////////////////////////////// 158 // 158 // 159 // Kill close intervals, recalculate fInterval 159 // Kill close intervals, recalculate fIntervalNumber 160 160 161 void G4InitXscPAI::Normalisation() 161 void G4InitXscPAI::Normalisation() 162 { 162 { 163 G4int i, j; 163 G4int i, j; 164 G4double energy1, energy2, /*delta,*/ cof; / 164 G4double energy1, energy2, /*delta,*/ cof; // , shift; 165 165 166 energy1 = (*(*fMatSandiaMatrix)[fIntervalNum 166 energy1 = (*(*fMatSandiaMatrix)[fIntervalNumber-1])[0]; 167 energy2 = 2.*(*(*fMatSandiaMatrix)[fInterval 167 energy2 = 2.*(*(*fMatSandiaMatrix)[fIntervalNumber-1])[0]; 168 168 169 169 170 cof = RutherfordIntegral(fIntervalNumber-1,e 170 cof = RutherfordIntegral(fIntervalNumber-1,energy1,energy2); 171 171 172 for( i = fIntervalNumber-2; i >= 0; i-- ) 172 for( i = fIntervalNumber-2; i >= 0; i-- ) 173 { 173 { 174 energy1 = (*(*fMatSandiaMatrix)[i])[0]; 174 energy1 = (*(*fMatSandiaMatrix)[i])[0]; 175 energy2 = (*(*fMatSandiaMatrix)[i+1])[0]; 175 energy2 = (*(*fMatSandiaMatrix)[i+1])[0]; 176 176 177 cof += RutherfordIntegral(i,energy1,energy 177 cof += RutherfordIntegral(i,energy1,energy2); 178 // G4cout<<"norm. cof = "<<cof<<G4endl; 178 // G4cout<<"norm. cof = "<<cof<<G4endl; 179 } 179 } 180 fNormalizationCof = 2*pi*pi*hbarc*hbarc*fin 180 fNormalizationCof = 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2 ; 181 fNormalizationCof *= fElectronDensity; 181 fNormalizationCof *= fElectronDensity; 182 //delta = fNormalizationCof - cof; 182 //delta = fNormalizationCof - cof; 183 fNormalizationCof /= cof; 183 fNormalizationCof /= cof; 184 // G4cout<<"G4InitXscPAI::fNormalizationCof 184 // G4cout<<"G4InitXscPAI::fNormalizationCof/cof = "<<fNormalizationCof 185 // <<"; at delta ="<<delta<<G4endl ; 185 // <<"; at delta ="<<delta<<G4endl ; 186 186 187 for (i = 0; i < fIntervalNumber; i++) // ren 187 for (i = 0; i < fIntervalNumber; i++) // renormalisation on QM sum rule 188 { 188 { 189 for(j = 1; j < 5 ; j++) 189 for(j = 1; j < 5 ; j++) 190 { 190 { 191 (*(*fMatSandiaMatrix)[i])[j] *= fNormali 191 (*(*fMatSandiaMatrix)[i])[j] *= fNormalizationCof; 192 } 192 } 193 } 193 } 194 /* 194 /* 195 if(delta > 0) // shift the first energy inte 195 if(delta > 0) // shift the first energy interval 196 { 196 { 197 for(i=1;i<100;i++) 197 for(i=1;i<100;i++) 198 { 198 { 199 energy1 = (1.-i/100.)*(*(*fMatSandiaMatr 199 energy1 = (1.-i/100.)*(*(*fMatSandiaMatrix)[0])[0]; 200 energy2 = (*(*fMatSandiaMatrix)[0])[0]; 200 energy2 = (*(*fMatSandiaMatrix)[0])[0]; 201 shift = RutherfordIntegral(0,energy1,e 201 shift = RutherfordIntegral(0,energy1,energy2); 202 G4cout<<shift<<"\t"; 202 G4cout<<shift<<"\t"; 203 if(shift >= delta) break; 203 if(shift >= delta) break; 204 } 204 } 205 (*(*fMatSandiaMatrix)[0])[0] = energy1; 205 (*(*fMatSandiaMatrix)[0])[0] = energy1; 206 cof += shift; 206 cof += shift; 207 } 207 } 208 else if(delta < 0) 208 else if(delta < 0) 209 { 209 { 210 for(i=1;i<100;i++) 210 for(i=1;i<100;i++) 211 { 211 { 212 energy1 = (*(*fMatSandiaMatrix)[0])[0]; 212 energy1 = (*(*fMatSandiaMatrix)[0])[0]; 213 energy2 = (*(*fMatSandiaMatrix)[0])[0] + 213 energy2 = (*(*fMatSandiaMatrix)[0])[0] + 214 ( (*(*fMatSandiaMatrix)[0])[0] - 214 ( (*(*fMatSandiaMatrix)[0])[0] - (*(*fMatSandiaMatrix)[0])[0] )*i/100.; 215 shift = RutherfordIntegral(0,energy1,e 215 shift = RutherfordIntegral(0,energy1,energy2); 216 if( shift >= std::abs(delta) ) break; 216 if( shift >= std::abs(delta) ) break; 217 } 217 } 218 (*(*fMatSandiaMatrix)[0])[0] = energy2; 218 (*(*fMatSandiaMatrix)[0])[0] = energy2; 219 cof -= shift; 219 cof -= shift; 220 } 220 } 221 G4cout<<G4cout<<"G4InitXscPAI::fNormalizatio 221 G4cout<<G4cout<<"G4InitXscPAI::fNormalizationCof/cof = "<<fNormalizationCof/cof 222 <<"; at delta ="<<delta<<" and i = " 222 <<"; at delta ="<<delta<<" and i = "<<i<<G4endl ; 223 */ 223 */ 224 } 224 } 225 225 226 ////////////////////////////////////////////// 226 //////////////////////////////////////////////////////////////////// 227 // 227 // 228 // Integration over electrons that could be co 228 // Integration over electrons that could be considered 229 // quasi-free at energy transfer of interest 229 // quasi-free at energy transfer of interest 230 230 231 G4double G4InitXscPAI::RutherfordIntegral( G4i 231 G4double G4InitXscPAI::RutherfordIntegral( G4int k, 232 G4double x1, 232 G4double x1, 233 G4double x2 ) 233 G4double x2 ) 234 { 234 { 235 G4double c1, c2, c3, a1, a2, a3, a4 ; 235 G4double c1, c2, c3, a1, a2, a3, a4 ; 236 236 237 a1 = (*(*fMatSandiaMatrix)[k])[1]; 237 a1 = (*(*fMatSandiaMatrix)[k])[1]; 238 a2 = (*(*fMatSandiaMatrix)[k])[2]; 238 a2 = (*(*fMatSandiaMatrix)[k])[2]; 239 a3 = (*(*fMatSandiaMatrix)[k])[3]; 239 a3 = (*(*fMatSandiaMatrix)[k])[3]; 240 a4 = (*(*fMatSandiaMatrix)[k])[4]; 240 a4 = (*(*fMatSandiaMatrix)[k])[4]; 241 // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<< 241 // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl; 242 c1 = (x2 - x1)/x1/x2 ; 242 c1 = (x2 - x1)/x1/x2 ; 243 c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2 ; 243 c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2 ; 244 c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x 244 c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ; 245 // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "< 245 // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<<c2<<"; "<<"c3 = "<<c3<<G4endl; 246 246 247 return a1*log(x2/x1) + a2*c1 + a3*c2/2 + a 247 return a1*log(x2/x1) + a2*c1 + a3*c2/2 + a4*c3/3 ; 248 248 249 } // end of RutherfordIntegral 249 } // end of RutherfordIntegral 250 250 251 ////////////////////////////////////////////// 251 /////////////////////////////////////////////////////////////// 252 // 252 // 253 // Integrate photo-absorption cross-section f 253 // Integrate photo-absorption cross-section from I1 up to omega 254 254 255 G4double G4InitXscPAI::IntegralTerm(G4double o 255 G4double G4InitXscPAI::IntegralTerm(G4double omega) 256 { 256 { 257 G4int i; 257 G4int i; 258 G4double energy1, energy2, result = 0.; 258 G4double energy1, energy2, result = 0.; 259 259 260 for( i = 0; i <= fIntervalTmax; i++ ) 260 for( i = 0; i <= fIntervalTmax; i++ ) 261 { 261 { 262 if(i == fIntervalTmax) 262 if(i == fIntervalTmax) 263 { 263 { 264 energy1 = (*(*fMatSandiaMatrix)[i])[0]; 264 energy1 = (*(*fMatSandiaMatrix)[i])[0]; 265 result += RutherfordIntegral(i,energy1,o 265 result += RutherfordIntegral(i,energy1,omega); 266 } 266 } 267 else 267 else 268 { 268 { 269 if( omega <= (*(*fMatSandiaMatrix)[i+1]) 269 if( omega <= (*(*fMatSandiaMatrix)[i+1])[0]) 270 { 270 { 271 energy1 = (*(*fMatSandiaMatrix)[i])[0] 271 energy1 = (*(*fMatSandiaMatrix)[i])[0]; 272 result += RutherfordIntegral(i,energy1 272 result += RutherfordIntegral(i,energy1,omega); 273 break; 273 break; 274 } 274 } 275 else 275 else 276 { 276 { 277 energy1 = (*(*fMatSandiaMatrix)[i])[0] 277 energy1 = (*(*fMatSandiaMatrix)[i])[0]; 278 energy2 = (*(*fMatSandiaMatrix)[i+1])[ 278 energy2 = (*(*fMatSandiaMatrix)[i+1])[0]; 279 result += RutherfordIntegral(i,energy1 279 result += RutherfordIntegral(i,energy1,energy2); 280 } 280 } 281 } 281 } 282 // G4cout<<"IntegralTerm<<"("<<omega<<")"< 282 // G4cout<<"IntegralTerm<<"("<<omega<<")"<<" = "<<result<<G4endl; 283 } 283 } 284 return result; 284 return result; 285 } 285 } 286 286 287 287 288 ////////////////////////////////////////////// 288 //////////////////////////////////////////////////////////////// 289 // 289 // 290 // Imaginary part of dielectric constant 290 // Imaginary part of dielectric constant 291 // (G4int k - interval number, G4double en1 - 291 // (G4int k - interval number, G4double en1 - energy point) 292 292 293 G4double G4InitXscPAI::ImPartDielectricConst( 293 G4double G4InitXscPAI::ImPartDielectricConst( G4int k , 294 G4double energy1 294 G4double energy1 ) 295 { 295 { 296 G4double energy2,energy3,energy4,a1,a2,a3,a 296 G4double energy2,energy3,energy4,a1,a2,a3,a4,result; 297 297 298 a1 = (*(*fMatSandiaMatrix)[k])[1]; 298 a1 = (*(*fMatSandiaMatrix)[k])[1]; 299 a2 = (*(*fMatSandiaMatrix)[k])[2]; 299 a2 = (*(*fMatSandiaMatrix)[k])[2]; 300 a3 = (*(*fMatSandiaMatrix)[k])[3]; 300 a3 = (*(*fMatSandiaMatrix)[k])[3]; 301 a4 = (*(*fMatSandiaMatrix)[k])[4]; 301 a4 = (*(*fMatSandiaMatrix)[k])[4]; 302 302 303 energy2 = energy1*energy1; 303 energy2 = energy1*energy1; 304 energy3 = energy2*energy1; 304 energy3 = energy2*energy1; 305 energy4 = energy3*energy1; 305 energy4 = energy3*energy1; 306 306 307 result = a1/energy1+a2/energy2+a3/energy3+ 307 result = a1/energy1+a2/energy2+a3/energy3+a4/energy4 ; 308 result *= hbarc/energy1 ; 308 result *= hbarc/energy1 ; 309 309 310 return result ; 310 return result ; 311 311 312 } // end of ImPartDielectricConst 312 } // end of ImPartDielectricConst 313 313 314 ////////////////////////////////////////////// 314 //////////////////////////////////////////////////////////////// 315 // 315 // 316 // Modulus squared of dielectric constant 316 // Modulus squared of dielectric constant 317 // (G4int k - interval number, G4double omega 317 // (G4int k - interval number, G4double omega - energy point) 318 318 319 G4double G4InitXscPAI::ModuleSqDielectricConst 319 G4double G4InitXscPAI::ModuleSqDielectricConst( G4int k , 320 G4double omega ) 320 G4double omega ) 321 { 321 { 322 G4double eIm2, eRe2, result; 322 G4double eIm2, eRe2, result; 323 323 324 result = ImPartDielectricConst(k,omega); 324 result = ImPartDielectricConst(k,omega); 325 eIm2 = result*result; 325 eIm2 = result*result; 326 326 327 result = RePartDielectricConst(omega); 327 result = RePartDielectricConst(omega); 328 eRe2 = result*result; 328 eRe2 = result*result; 329 329 330 result = eIm2 + eRe2; 330 result = eIm2 + eRe2; 331 331 332 return result ; 332 return result ; 333 } 333 } 334 334 335 335 336 ////////////////////////////////////////////// 336 ////////////////////////////////////////////////////////////////////////////// 337 // 337 // 338 // Real part of dielectric constant minus unit 338 // Real part of dielectric constant minus unit: epsilon_1 - 1 339 // (G4double enb - energy point) 339 // (G4double enb - energy point) 340 // 340 // 341 341 342 G4double G4InitXscPAI::RePartDielectricConst(G 342 G4double G4InitXscPAI::RePartDielectricConst(G4double enb) 343 { 343 { 344 G4int i; 344 G4int i; 345 G4double x0, x02, x03, x04, x05, x1, x2, a1 345 G4double x0, x02, x03, x04, x05, x1, x2, a1,a2,a3,a4,xx1 ,xx2 , xx12, 346 c1, c2, c3, cof1, cof2, xln1, xln2 346 c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result ; 347 347 348 x0 = enb ; 348 x0 = enb ; 349 result = 0 ; 349 result = 0 ; 350 350 351 for( i = 0; i < fIntervalNumber-1; i++) 351 for( i = 0; i < fIntervalNumber-1; i++) 352 { 352 { 353 x1 = (*(*fMatSandiaMatrix)[i])[0]; 353 x1 = (*(*fMatSandiaMatrix)[i])[0]; 354 x2 = (*(*fMatSandiaMatrix)[i+1])[0] ; 354 x2 = (*(*fMatSandiaMatrix)[i+1])[0] ; 355 355 356 a1 = (*(*fMatSandiaMatrix)[i])[1]; 356 a1 = (*(*fMatSandiaMatrix)[i])[1]; 357 a2 = (*(*fMatSandiaMatrix)[i])[2]; 357 a2 = (*(*fMatSandiaMatrix)[i])[2]; 358 a3 = (*(*fMatSandiaMatrix)[i])[3]; 358 a3 = (*(*fMatSandiaMatrix)[i])[3]; 359 a4 = (*(*fMatSandiaMatrix)[i])[4]; 359 a4 = (*(*fMatSandiaMatrix)[i])[4]; 360 360 361 if( std::abs(x0-x1) < 0.5*(x0+x1)*fDelta 361 if( std::abs(x0-x1) < 0.5*(x0+x1)*fDelta ) 362 { 362 { 363 if(x0 >= x1) x0 = x1*(1+fDelta); 363 if(x0 >= x1) x0 = x1*(1+fDelta); 364 else x0 = x1*(1-fDelta); 364 else x0 = x1*(1-fDelta); 365 } 365 } 366 if( std::abs(x0-x2) < 0.5*(x0+x2)*fDelta 366 if( std::abs(x0-x2) < 0.5*(x0+x2)*fDelta ) 367 { 367 { 368 if(x0 >= x2) x0 = x2*(1+fDelta); 368 if(x0 >= x2) x0 = x2*(1+fDelta); 369 else x0 = x2*(1-fDelta); 369 else x0 = x2*(1-fDelta); 370 } 370 } 371 xx1 = x1 - x0 ; 371 xx1 = x1 - x0 ; 372 xx2 = x2 - x0 ; 372 xx2 = x2 - x0 ; 373 xx12 = xx2/xx1 ; 373 xx12 = xx2/xx1 ; 374 374 375 if( xx12 < 0 ) xx12 = -xx12; 375 if( xx12 < 0 ) xx12 = -xx12; 376 376 377 xln1 = log(x2/x1) ; 377 xln1 = log(x2/x1) ; 378 xln2 = log(xx12) ; 378 xln2 = log(xx12) ; 379 xln3 = log((x2 + x0)/(x1 + x0)) ; 379 xln3 = log((x2 + x0)/(x1 + x0)) ; 380 380 381 x02 = x0*x0 ; 381 x02 = x0*x0 ; 382 x03 = x02*x0 ; 382 x03 = x02*x0 ; 383 x04 = x03*x0 ; 383 x04 = x03*x0 ; 384 x05 = x04*x0; 384 x05 = x04*x0; 385 385 386 c1 = (x2 - x1)/x1/x2 ; 386 c1 = (x2 - x1)/x1/x2 ; 387 c2 = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2 ; 387 c2 = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2 ; 388 c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x 388 c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ; 389 389 390 result -= (a1/x02 + a3/x04)*xln1 ; 390 result -= (a1/x02 + a3/x04)*xln1 ; 391 result -= (a2/x02 + a4/x04)*c1 ; 391 result -= (a2/x02 + a4/x04)*c1 ; 392 result -= a3*c2/2/x02 ; 392 result -= a3*c2/2/x02 ; 393 result -= a4*c3/3/x02 ; 393 result -= a4*c3/3/x02 ; 394 394 395 cof1 = a1/x02 + a3/x04 ; 395 cof1 = a1/x02 + a3/x04 ; 396 cof2 = a2/x03 + a4/x05 ; 396 cof2 = a2/x03 + a4/x05 ; 397 397 398 result += 0.5*(cof1 +cof2)*xln2 ; 398 result += 0.5*(cof1 +cof2)*xln2 ; 399 result += 0.5*(cof1 - cof2)*xln3 ; 399 result += 0.5*(cof1 - cof2)*xln3 ; 400 } 400 } 401 result *= 2*hbarc/pi ; 401 result *= 2*hbarc/pi ; 402 402 403 return result ; 403 return result ; 404 404 405 } // end of RePartDielectricConst 405 } // end of RePartDielectricConst 406 406 407 ////////////////////////////////////////////// 407 ////////////////////////////////////////////////////////////////////// 408 // 408 // 409 // PAI differential cross-section in terms of 409 // PAI differential cross-section in terms of 410 // simplified Allison's equation 410 // simplified Allison's equation 411 // 411 // 412 412 413 G4double G4InitXscPAI::DifPAIxSection( G4doubl 413 G4double G4InitXscPAI::DifPAIxSection( G4double omega ) 414 { 414 { 415 G4int i = fCurrentInterval; 415 G4int i = fCurrentInterval; 416 G4double betaGammaSq = fBetaGammaSq; 416 G4double betaGammaSq = fBetaGammaSq; 417 G4double integralTerm = IntegralTerm(omega); 417 G4double integralTerm = IntegralTerm(omega); 418 G4double be2,cof,x1,x2,x3,x4,x5,x6,x7,x8,res 418 G4double be2,cof,x1,x2,x3,x4,x5,x6,x7,x8,result ; 419 G4double epsilonRe = RePartDielectricConst(o 419 G4double epsilonRe = RePartDielectricConst(omega); 420 G4double epsilonIm = ImPartDielectricConst(i 420 G4double epsilonIm = ImPartDielectricConst(i,omega); 421 G4double be4 ; 421 G4double be4 ; 422 static const G4double betaBohr2 = fine_struc 422 static const G4double betaBohr2 = fine_structure_const*fine_structure_const ; 423 static const G4double betaBohr4 = betaBohr2* 423 static const G4double betaBohr4 = betaBohr2*betaBohr2*4.0 ; 424 be2 = betaGammaSq/(1 + betaGammaSq) ; 424 be2 = betaGammaSq/(1 + betaGammaSq) ; 425 be4 = be2*be2 ; 425 be4 = be2*be2 ; 426 426 427 cof = 1 ; 427 cof = 1 ; 428 x1 = log(2*electron_mass_c2/omega) ; 428 x1 = log(2*electron_mass_c2/omega) ; 429 429 430 if( betaGammaSq < 0.01 ) x2 = log(be2) ; 430 if( betaGammaSq < 0.01 ) x2 = log(be2) ; 431 else 431 else 432 { 432 { 433 x2 = -log( (1/betaGammaSq - epsilonRe)* 433 x2 = -log( (1/betaGammaSq - epsilonRe)* 434 (1/betaGammaSq - epsilonRe) + 434 (1/betaGammaSq - epsilonRe) + 435 epsilonIm*epsilonIm )/2 ; 435 epsilonIm*epsilonIm )/2 ; 436 } 436 } 437 if( epsilonIm == 0.0 || betaGammaSq < 0.01 437 if( epsilonIm == 0.0 || betaGammaSq < 0.01 ) 438 { 438 { 439 x6=0 ; 439 x6=0 ; 440 } 440 } 441 else 441 else 442 { 442 { 443 x3 = -epsilonRe + 1/betaGammaSq ; 443 x3 = -epsilonRe + 1/betaGammaSq ; 444 x5 = -1 - epsilonRe + be2*((1 +epsilonRe) 444 x5 = -1 - epsilonRe + be2*((1 +epsilonRe)*(1 + epsilonRe) + 445 epsilonIm*epsilonIm) ; 445 epsilonIm*epsilonIm) ; 446 446 447 x7 = atan2(epsilonIm,x3) ; 447 x7 = atan2(epsilonIm,x3) ; 448 x6 = x5 * x7 ; 448 x6 = x5 * x7 ; 449 } 449 } 450 // if(fImPartDielectricConst[i] == 0) x6 = 450 // if(fImPartDielectricConst[i] == 0) x6 = 0 ; 451 451 452 x4 = ((x1 + x2)*epsilonIm + x6)/hbarc ; 452 x4 = ((x1 + x2)*epsilonIm + x6)/hbarc ; 453 // if( x4 < 0.0 ) x4 = 0.0 ; 453 // if( x4 < 0.0 ) x4 = 0.0 ; 454 x8 = (1 + epsilonRe)*(1 + epsilonRe) + 454 x8 = (1 + epsilonRe)*(1 + epsilonRe) + 455 epsilonIm*epsilonIm; 455 epsilonIm*epsilonIm; 456 456 457 result = (x4 + cof*integralTerm/omega/omega 457 result = (x4 + cof*integralTerm/omega/omega) ; 458 if(result < 1.0e-8) result = 1.0e-8 ; 458 if(result < 1.0e-8) result = 1.0e-8 ; 459 result *= fine_structure_const/be2/pi ; 459 result *= fine_structure_const/be2/pi ; 460 // result *= (1-exp(-beta/betaBohr))*(1-e 460 // result *= (1-exp(-beta/betaBohr))*(1-exp(-beta/betaBohr)) ; 461 // result *= (1-exp(-be2/betaBohr2)) ; 461 // result *= (1-exp(-be2/betaBohr2)) ; 462 result *= (1-exp(-be4/betaBohr4)) ; 462 result *= (1-exp(-be4/betaBohr4)) ; 463 if(fDensity >= fSolidDensity) 463 if(fDensity >= fSolidDensity) 464 { 464 { 465 result /= x8 ; 465 result /= x8 ; 466 } 466 } 467 return result ; 467 return result ; 468 468 469 } // end of DifPAIxSection 469 } // end of DifPAIxSection 470 470 471 ////////////////////////////////////////////// 471 ////////////////////////////////////////////////////////////////////// 472 // 472 // 473 // Differential PAI dEdx(omega)=omega*dNdx(ome 473 // Differential PAI dEdx(omega)=omega*dNdx(omega) 474 // 474 // 475 475 476 G4double G4InitXscPAI::DifPAIdEdx( G4double om 476 G4double G4InitXscPAI::DifPAIdEdx( G4double omega ) 477 { 477 { 478 G4double dEdx = omega*DifPAIxSection(omega); 478 G4double dEdx = omega*DifPAIxSection(omega); 479 return dEdx; 479 return dEdx; 480 } 480 } 481 481 482 ////////////////////////////////////////////// 482 ////////////////////////////////////////////////////////////////////////// 483 // 483 // 484 // Calculation od dN/dx of collisions with cre 484 // Calculation od dN/dx of collisions with creation of Cerenkov pseudo-photons 485 485 486 G4double G4InitXscPAI::PAIdNdxCherenkov( G4dou 486 G4double G4InitXscPAI::PAIdNdxCherenkov( G4double omega ) 487 { 487 { 488 G4int i = fCurrentInterval; 488 G4int i = fCurrentInterval; 489 G4double betaGammaSq = fBetaGammaSq; 489 G4double betaGammaSq = fBetaGammaSq; 490 G4double epsilonRe = RePartDielectricConst(o 490 G4double epsilonRe = RePartDielectricConst(omega); 491 G4double epsilonIm = ImPartDielectricConst(i 491 G4double epsilonIm = ImPartDielectricConst(i,omega); 492 492 493 G4double /*cof,*/ logarithm, x3, x5, argumen 493 G4double /*cof,*/ logarithm, x3, x5, argument, modul2, dNdxC ; 494 G4double be2, be4; 494 G4double be2, be4; 495 495 496 //cof = 1.0 ; 496 //cof = 1.0 ; 497 static const G4double cofBetaBohr = 4.0 ; 497 static const G4double cofBetaBohr = 4.0 ; 498 static const G4double betaBohr2 = fine_str 498 static const G4double betaBohr2 = fine_structure_const*fine_structure_const ; 499 static const G4double betaBohr4 = betaBohr 499 static const G4double betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr ; 500 500 501 be2 = betaGammaSq/(1 + betaGammaSq) ; 501 be2 = betaGammaSq/(1 + betaGammaSq) ; 502 be4 = be2*be2 ; 502 be4 = be2*be2 ; 503 503 504 if( betaGammaSq < 0.01 ) logarithm = log(1. 504 if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq) ; // 0.0 ; 505 else 505 else 506 { 506 { 507 logarithm = -log( (1/betaGammaSq - epsil 507 logarithm = -log( (1/betaGammaSq - epsilonRe)* 508 (1/betaGammaSq - epsilonRe) 508 (1/betaGammaSq - epsilonRe) + 509 epsilonIm*epsilonIm )*0.5 ; 509 epsilonIm*epsilonIm )*0.5 ; 510 logarithm += log(1+1.0/betaGammaSq) ; 510 logarithm += log(1+1.0/betaGammaSq) ; 511 } 511 } 512 512 513 if( epsilonIm == 0.0 || betaGammaSq < 0.01 513 if( epsilonIm == 0.0 || betaGammaSq < 0.01 ) 514 { 514 { 515 argument = 0.0 ; 515 argument = 0.0 ; 516 } 516 } 517 else 517 else 518 { 518 { 519 x3 = -epsilonRe + 1.0/betaGammaSq ; 519 x3 = -epsilonRe + 1.0/betaGammaSq ; 520 x5 = -1.0 - epsilonRe + 520 x5 = -1.0 - epsilonRe + 521 be2*((1.0 +epsilonRe)*(1.0 + epsilon 521 be2*((1.0 +epsilonRe)*(1.0 + epsilonRe) + 522 epsilonIm*epsilonIm) ; 522 epsilonIm*epsilonIm) ; 523 if( x3 == 0.0 ) argument = 0.5*pi; 523 if( x3 == 0.0 ) argument = 0.5*pi; 524 else argument = atan2(epsilonI 524 else argument = atan2(epsilonIm,x3) ; 525 argument *= x5 ; 525 argument *= x5 ; 526 } 526 } 527 dNdxC = ( logarithm*epsilonIm + argument )/ 527 dNdxC = ( logarithm*epsilonIm + argument )/hbarc ; 528 528 529 if(dNdxC < 1.0e-8) dNdxC = 1.0e-8 ; 529 if(dNdxC < 1.0e-8) dNdxC = 1.0e-8 ; 530 530 531 dNdxC *= fine_structure_const/be2/pi ; 531 dNdxC *= fine_structure_const/be2/pi ; 532 532 533 dNdxC *= (1-exp(-be4/betaBohr4)) ; 533 dNdxC *= (1-exp(-be4/betaBohr4)) ; 534 534 535 if(fDensity >= fSolidDensity) 535 if(fDensity >= fSolidDensity) 536 { 536 { 537 modul2 = (1.0 + epsilonRe)*(1.0 + epsilo 537 modul2 = (1.0 + epsilonRe)*(1.0 + epsilonRe) + 538 epsilonIm*epsilonIm; 538 epsilonIm*epsilonIm; 539 dNdxC /= modul2 ; 539 dNdxC /= modul2 ; 540 } 540 } 541 return dNdxC ; 541 return dNdxC ; 542 542 543 } // end of PAIdNdxCerenkov 543 } // end of PAIdNdxCerenkov 544 544 545 ////////////////////////////////////////////// 545 ////////////////////////////////////////////////////////////////////////// 546 // 546 // 547 // Calculation od dN/dx of collisions with cre 547 // Calculation od dN/dx of collisions with creation of longitudinal EM 548 // excitations (plasmons, delta-electrons) 548 // excitations (plasmons, delta-electrons) 549 549 550 G4double G4InitXscPAI::PAIdNdxPlasmon( G4doubl 550 G4double G4InitXscPAI::PAIdNdxPlasmon( G4double omega ) 551 { 551 { 552 G4int i = fCurrentInterval; 552 G4int i = fCurrentInterval; 553 G4double betaGammaSq = fBetaGammaSq; 553 G4double betaGammaSq = fBetaGammaSq; 554 G4double integralTerm = IntegralTerm(omega); 554 G4double integralTerm = IntegralTerm(omega); 555 G4double epsilonRe = RePartDielectricConst(o 555 G4double epsilonRe = RePartDielectricConst(omega); 556 G4double epsilonIm = ImPartDielectricConst(i 556 G4double epsilonIm = ImPartDielectricConst(i,omega); 557 557 558 G4double cof, resonance, modul2, dNdxP ; 558 G4double cof, resonance, modul2, dNdxP ; 559 G4double be2, be4; 559 G4double be2, be4; 560 560 561 cof = 1 ; 561 cof = 1 ; 562 static const G4double cofBetaBohr = 4.0 ; 562 static const G4double cofBetaBohr = 4.0 ; 563 static const G4double betaBohr2 = fine_st 563 static const G4double betaBohr2 = fine_structure_const*fine_structure_const ; 564 static const G4double betaBohr4 = betaBoh 564 static const G4double betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr ; 565 565 566 be2 = betaGammaSq/(1 + betaGammaSq) ; 566 be2 = betaGammaSq/(1 + betaGammaSq) ; 567 be4 = be2*be2 ; 567 be4 = be2*be2 ; 568 568 569 resonance = log(2*electron_mass_c2*be2/ome 569 resonance = log(2*electron_mass_c2*be2/omega) ; 570 resonance *= epsilonIm/hbarc ; 570 resonance *= epsilonIm/hbarc ; 571 571 572 572 573 dNdxP = ( resonance + cof*integralTerm/omeg 573 dNdxP = ( resonance + cof*integralTerm/omega/omega ) ; 574 574 575 if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8 ; 575 if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8 ; 576 576 577 dNdxP *= fine_structure_const/be2/pi ; 577 dNdxP *= fine_structure_const/be2/pi ; 578 dNdxP *= (1-exp(-be4/betaBohr4)) ; 578 dNdxP *= (1-exp(-be4/betaBohr4)) ; 579 579 580 if( fDensity >= fSolidDensity ) 580 if( fDensity >= fSolidDensity ) 581 { 581 { 582 modul2 = (1 + epsilonRe)*(1 + epsilonRe) 582 modul2 = (1 + epsilonRe)*(1 + epsilonRe) + 583 epsilonIm*epsilonIm; 583 epsilonIm*epsilonIm; 584 dNdxP /= modul2 ; 584 dNdxP /= modul2 ; 585 } 585 } 586 return dNdxP ; 586 return dNdxP ; 587 587 588 } // end of PAIdNdxPlasmon 588 } // end of PAIdNdxPlasmon 589 589 590 ////////////////////////////////////////////// 590 //////////////////////////////////////////////////////////////////////// 591 // 591 // 592 // Calculation of the PAI integral cross-secti 592 // Calculation of the PAI integral cross-section 593 // = specific primary ionisation, 1/cm 593 // = specific primary ionisation, 1/cm 594 // 594 // 595 595 596 void G4InitXscPAI::IntegralPAIxSection(G4doubl 596 void G4InitXscPAI::IntegralPAIxSection(G4double bg2, G4double Tmax) 597 { 597 { 598 G4int i,k,i1,i2; 598 G4int i,k,i1,i2; 599 G4double energy1, energy2, result = 0.; 599 G4double energy1, energy2, result = 0.; 600 600 601 fBetaGammaSq = bg2; 601 fBetaGammaSq = bg2; 602 fTmax = Tmax; 602 fTmax = Tmax; 603 603 604 delete fPAIxscVector; << 604 if(fPAIxscVector) delete fPAIxscVector; 605 605 606 fPAIxscVector = new G4PhysicsLogVector( (*(* 606 fPAIxscVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin); 607 fPAIxscVector->PutValue(fPAIbin-1,result); 607 fPAIxscVector->PutValue(fPAIbin-1,result); 608 608 609 for( i = fIntervalNumber - 1; i >= 0; i-- ) 609 for( i = fIntervalNumber - 1; i >= 0; i-- ) 610 { 610 { 611 if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) 611 if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break; 612 } 612 } 613 if (i < 0) i = 0; // Tmax should be more tha 613 if (i < 0) i = 0; // Tmax should be more than 614 // first ionisation potent 614 // first ionisation potential 615 fIntervalTmax = i; 615 fIntervalTmax = i; 616 616 617 G4Integrator<G4InitXscPAI,G4double(G4InitXsc 617 G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral; 618 618 619 for( k = fPAIbin - 2; k >= 0; k-- ) 619 for( k = fPAIbin - 2; k >= 0; k-- ) 620 { 620 { 621 energy1 = fPAIxscVector->GetLowEdgeEnergy( 621 energy1 = fPAIxscVector->GetLowEdgeEnergy(k); 622 energy2 = fPAIxscVector->GetLowEdgeEnergy( 622 energy2 = fPAIxscVector->GetLowEdgeEnergy(k+1); 623 623 624 for( i = fIntervalTmax; i >= 0; i-- ) 624 for( i = fIntervalTmax; i >= 0; i-- ) 625 { 625 { 626 if( energy2 > (*(*fMatSandiaMatrix)[i])[ 626 if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break; 627 } 627 } 628 if(i < 0) i = 0; 628 if(i < 0) i = 0; 629 i2 = i; 629 i2 = i; 630 630 631 for( i = fIntervalTmax; i >= 0; i-- ) 631 for( i = fIntervalTmax; i >= 0; i-- ) 632 { 632 { 633 if( energy1 > (*(*fMatSandiaMatrix)[i])[ 633 if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break; 634 } 634 } 635 if(i < 0) i = 0; 635 if(i < 0) i = 0; 636 i1 = i; 636 i1 = i; 637 637 638 if( i1 == i2 ) 638 if( i1 == i2 ) 639 { 639 { 640 fCurrentInterval = i1; 640 fCurrentInterval = i1; 641 result += integral.Legendre10(this,&G4In 641 result += integral.Legendre10(this,&G4InitXscPAI::DifPAIxSection, 642 energy1,en 642 energy1,energy2); 643 fPAIxscVector->PutValue(k,result); 643 fPAIxscVector->PutValue(k,result); 644 } 644 } 645 else 645 else 646 { 646 { 647 for( i = i2; i >= i1; i-- ) 647 for( i = i2; i >= i1; i-- ) 648 { 648 { 649 fCurrentInterval = i; 649 fCurrentInterval = i; 650 650 651 if( i==i2 ) result += integral. 651 if( i==i2 ) result += integral.Legendre10(this, 652 &G4InitXscPAI::DifP 652 &G4InitXscPAI::DifPAIxSection, 653 (*(*fMatSandiaMatri 653 (*(*fMatSandiaMatrix)[i])[0] ,energy2); 654 654 655 else if( i == i1 ) result += integral.Legend 655 else if( i == i1 ) result += integral.Legendre10(this, 656 &G4InitXscPAI::DifP 656 &G4InitXscPAI::DifPAIxSection,energy1, 657 (*(*fMatSandiaMatri 657 (*(*fMatSandiaMatrix)[i+1])[0]); 658 658 659 else result += integral. 659 else result += integral.Legendre10(this, 660 &G4InitXscPAI::DifP 660 &G4InitXscPAI::DifPAIxSection, 661 (*(*fMatSandiaMatrix)[i 661 (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]); 662 } 662 } 663 fPAIxscVector->PutValue(k,result); 663 fPAIxscVector->PutValue(k,result); 664 } 664 } 665 // G4cout<<k<<"\t"<<result<<G4endl; 665 // G4cout<<k<<"\t"<<result<<G4endl; 666 } 666 } 667 return ; 667 return ; 668 } 668 } 669 669 670 670 671 ////////////////////////////////////////////// 671 //////////////////////////////////////////////////////////////////////// 672 // 672 // 673 // Calculation of the PAI integral dEdx 673 // Calculation of the PAI integral dEdx 674 // = mean energy loss per unit length, keV/cm 674 // = mean energy loss per unit length, keV/cm 675 // 675 // 676 676 677 void G4InitXscPAI::IntegralPAIdEdx(G4double bg 677 void G4InitXscPAI::IntegralPAIdEdx(G4double bg2, G4double Tmax) 678 { 678 { 679 G4int i,k,i1,i2; 679 G4int i,k,i1,i2; 680 G4double energy1, energy2, result = 0.; 680 G4double energy1, energy2, result = 0.; 681 681 682 fBetaGammaSq = bg2; 682 fBetaGammaSq = bg2; 683 fTmax = Tmax; 683 fTmax = Tmax; 684 684 685 delete fPAIdEdxVector; << 685 if(fPAIdEdxVector) delete fPAIdEdxVector; 686 686 687 fPAIdEdxVector = new G4PhysicsLogVector( (*( 687 fPAIdEdxVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin); 688 fPAIdEdxVector->PutValue(fPAIbin-1,result); 688 fPAIdEdxVector->PutValue(fPAIbin-1,result); 689 689 690 for( i = fIntervalNumber - 1; i >= 0; i-- ) 690 for( i = fIntervalNumber - 1; i >= 0; i-- ) 691 { 691 { 692 if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) 692 if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break; 693 } 693 } 694 if (i < 0) i = 0; // Tmax should be more tha 694 if (i < 0) i = 0; // Tmax should be more than 695 // first ionisation potent 695 // first ionisation potential 696 fIntervalTmax = i; 696 fIntervalTmax = i; 697 697 698 G4Integrator<G4InitXscPAI,G4double(G4InitXsc 698 G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral; 699 699 700 for( k = fPAIbin - 2; k >= 0; k-- ) 700 for( k = fPAIbin - 2; k >= 0; k-- ) 701 { 701 { 702 energy1 = fPAIdEdxVector->GetLowEdgeEnergy 702 energy1 = fPAIdEdxVector->GetLowEdgeEnergy(k); 703 energy2 = fPAIdEdxVector->GetLowEdgeEnergy 703 energy2 = fPAIdEdxVector->GetLowEdgeEnergy(k+1); 704 704 705 for( i = fIntervalTmax; i >= 0; i-- ) 705 for( i = fIntervalTmax; i >= 0; i-- ) 706 { 706 { 707 if( energy2 > (*(*fMatSandiaMatrix)[i])[ 707 if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break; 708 } 708 } 709 if(i < 0) i = 0; 709 if(i < 0) i = 0; 710 i2 = i; 710 i2 = i; 711 711 712 for( i = fIntervalTmax; i >= 0; i-- ) 712 for( i = fIntervalTmax; i >= 0; i-- ) 713 { 713 { 714 if( energy1 > (*(*fMatSandiaMatrix)[i])[ 714 if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break; 715 } 715 } 716 if(i < 0) i = 0; 716 if(i < 0) i = 0; 717 i1 = i; 717 i1 = i; 718 718 719 if( i1 == i2 ) 719 if( i1 == i2 ) 720 { 720 { 721 fCurrentInterval = i1; 721 fCurrentInterval = i1; 722 result += integral.Legendre10(this,&G4In 722 result += integral.Legendre10(this,&G4InitXscPAI::DifPAIdEdx, 723 energy1,en 723 energy1,energy2); 724 fPAIdEdxVector->PutValue(k,result); 724 fPAIdEdxVector->PutValue(k,result); 725 } 725 } 726 else 726 else 727 { 727 { 728 for( i = i2; i >= i1; i-- ) 728 for( i = i2; i >= i1; i-- ) 729 { 729 { 730 fCurrentInterval = i; 730 fCurrentInterval = i; 731 731 732 if( i==i2 ) result += integral. 732 if( i==i2 ) result += integral.Legendre10(this, 733 &G4InitXscPAI::DifP 733 &G4InitXscPAI::DifPAIdEdx, 734 (*(*fMatSandiaMatri 734 (*(*fMatSandiaMatrix)[i])[0] ,energy2); 735 735 736 else if( i == i1 ) result += integral.Legend 736 else if( i == i1 ) result += integral.Legendre10(this, 737 &G4InitXscPAI::DifP 737 &G4InitXscPAI::DifPAIdEdx,energy1, 738 (*(*fMatSandiaMatri 738 (*(*fMatSandiaMatrix)[i+1])[0]); 739 739 740 else result += integral. 740 else result += integral.Legendre10(this, 741 &G4InitXscPAI::DifP 741 &G4InitXscPAI::DifPAIdEdx, 742 (*(*fMatSandiaMatrix)[i 742 (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]); 743 } 743 } 744 fPAIdEdxVector->PutValue(k,result); 744 fPAIdEdxVector->PutValue(k,result); 745 } 745 } 746 // G4cout<<k<<"\t"<<result<<G4endl; 746 // G4cout<<k<<"\t"<<result<<G4endl; 747 } 747 } 748 return ; 748 return ; 749 } 749 } 750 750 751 ////////////////////////////////////////////// 751 //////////////////////////////////////////////////////////////////////// 752 // 752 // 753 // Calculation of the PAI Cerenkov integral cr 753 // Calculation of the PAI Cerenkov integral cross-section 754 // fIntegralCrenkov[1] = specific Crenkov ioni 754 // fIntegralCrenkov[1] = specific Crenkov ionisation, 1/cm 755 // and fIntegralCerenkov[0] = mean Cerenkov lo 755 // and fIntegralCerenkov[0] = mean Cerenkov loss per cm in keV/cm 756 756 757 void G4InitXscPAI::IntegralCherenkov(G4double 757 void G4InitXscPAI::IntegralCherenkov(G4double bg2, G4double Tmax) 758 { 758 { 759 G4int i,k,i1,i2; 759 G4int i,k,i1,i2; 760 G4double energy1, energy2, beta2, module2, c 760 G4double energy1, energy2, beta2, module2, cos2, width, result = 0.; 761 761 762 fBetaGammaSq = bg2; 762 fBetaGammaSq = bg2; 763 fTmax = Tmax; 763 fTmax = Tmax; 764 beta2 = bg2/(1+bg2); 764 beta2 = bg2/(1+bg2); 765 765 766 delete fPAIphotonVector; << 766 if(fPAIphotonVector) delete fPAIphotonVector; 767 delete fChCosSqVector; << 767 if(fChCosSqVector) delete fChCosSqVector; 768 delete fChWidthVector; << 768 if(fChWidthVector) delete fChWidthVector; 769 769 770 fPAIphotonVector = new G4PhysicsLogVector( ( 770 fPAIphotonVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin); 771 fChCosSqVector = new G4PhysicsLogVector( (*( 771 fChCosSqVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin); 772 fChWidthVector = new G4PhysicsLogVector( (*( 772 fChWidthVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin); 773 773 774 fPAIphotonVector->PutValue(fPAIbin-1,result) 774 fPAIphotonVector->PutValue(fPAIbin-1,result); 775 fChCosSqVector->PutValue(fPAIbin-1,1.); 775 fChCosSqVector->PutValue(fPAIbin-1,1.); 776 fChWidthVector->PutValue(fPAIbin-1,1e-7); 776 fChWidthVector->PutValue(fPAIbin-1,1e-7); 777 777 778 for( i = fIntervalNumber - 1; i >= 0; i-- ) 778 for( i = fIntervalNumber - 1; i >= 0; i-- ) 779 { 779 { 780 if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) 780 if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break; 781 } 781 } 782 if (i < 0) i = 0; // Tmax should be more tha 782 if (i < 0) i = 0; // Tmax should be more than 783 // first ionisation potent 783 // first ionisation potential 784 fIntervalTmax = i; 784 fIntervalTmax = i; 785 785 786 G4Integrator<G4InitXscPAI,G4double(G4InitXsc 786 G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral; 787 787 788 for( k = fPAIbin - 2; k >= 0; k-- ) 788 for( k = fPAIbin - 2; k >= 0; k-- ) 789 { 789 { 790 energy1 = fPAIphotonVector->GetLowEdgeEner 790 energy1 = fPAIphotonVector->GetLowEdgeEnergy(k); 791 energy2 = fPAIphotonVector->GetLowEdgeEner 791 energy2 = fPAIphotonVector->GetLowEdgeEnergy(k+1); 792 792 793 for( i = fIntervalTmax; i >= 0; i-- ) 793 for( i = fIntervalTmax; i >= 0; i-- ) 794 { 794 { 795 if( energy2 > (*(*fMatSandiaMatrix)[i])[ 795 if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break; 796 } 796 } 797 if(i < 0) i = 0; 797 if(i < 0) i = 0; 798 i2 = i; 798 i2 = i; 799 799 800 for( i = fIntervalTmax; i >= 0; i-- ) 800 for( i = fIntervalTmax; i >= 0; i-- ) 801 { 801 { 802 if( energy1 > (*(*fMatSandiaMatrix)[i])[ 802 if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break; 803 } 803 } 804 if(i < 0) i = 0; 804 if(i < 0) i = 0; 805 i1 = i; 805 i1 = i; 806 806 807 module2 = ModuleSqDielectricConst(i1,energ 807 module2 = ModuleSqDielectricConst(i1,energy1); 808 cos2 = RePartDielectricConst(energy1)/m 808 cos2 = RePartDielectricConst(energy1)/module2/beta2; 809 width = ImPartDielectricConst(i1,energy1 809 width = ImPartDielectricConst(i1,energy1)/module2/beta2; 810 810 811 fChCosSqVector->PutValue(k,cos2); 811 fChCosSqVector->PutValue(k,cos2); 812 fChWidthVector->PutValue(k,width); 812 fChWidthVector->PutValue(k,width); 813 813 814 if( i1 == i2 ) 814 if( i1 == i2 ) 815 { 815 { 816 fCurrentInterval = i1; 816 fCurrentInterval = i1; 817 result += integral.Legendre10(this,&G4In 817 result += integral.Legendre10(this,&G4InitXscPAI::PAIdNdxCherenkov, 818 energy1,en 818 energy1,energy2); 819 fPAIphotonVector->PutValue(k,result); 819 fPAIphotonVector->PutValue(k,result); 820 820 821 } 821 } 822 else 822 else 823 { 823 { 824 for( i = i2; i >= i1; i-- ) 824 for( i = i2; i >= i1; i-- ) 825 { 825 { 826 fCurrentInterval = i; 826 fCurrentInterval = i; 827 827 828 if( i==i2 ) result += integral. 828 if( i==i2 ) result += integral.Legendre10(this, 829 &G4InitXscPAI::PAId 829 &G4InitXscPAI::PAIdNdxCherenkov, 830 (*(*fMatSandiaMatri 830 (*(*fMatSandiaMatrix)[i])[0] ,energy2); 831 831 832 else if( i == i1 ) result += integral.Legend 832 else if( i == i1 ) result += integral.Legendre10(this, 833 &G4InitXscPAI::PAId 833 &G4InitXscPAI::PAIdNdxCherenkov,energy1, 834 (*(*fMatSandiaMatri 834 (*(*fMatSandiaMatrix)[i+1])[0]); 835 835 836 else result += integral. 836 else result += integral.Legendre10(this, 837 &G4InitXscPAI::PAId 837 &G4InitXscPAI::PAIdNdxCherenkov, 838 (*(*fMatSandiaMatrix)[i 838 (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]); 839 } 839 } 840 fPAIphotonVector->PutValue(k,result); 840 fPAIphotonVector->PutValue(k,result); 841 } 841 } 842 // G4cout<<k<<"\t"<<result<<G4endl; 842 // G4cout<<k<<"\t"<<result<<G4endl; 843 } 843 } 844 return; 844 return; 845 } // end of IntegralCerenkov 845 } // end of IntegralCerenkov 846 846 847 ////////////////////////////////////////////// 847 //////////////////////////////////////////////////////////////////////// 848 // 848 // 849 // Calculation of the PAI Plasmon integral cro 849 // Calculation of the PAI Plasmon integral cross-section 850 // fIntegralPlasmon[1] = splasmon primary ioni 850 // fIntegralPlasmon[1] = splasmon primary ionisation, 1/cm 851 // and fIntegralPlasmon[0] = mean plasmon loss 851 // and fIntegralPlasmon[0] = mean plasmon loss per cm in keV/cm 852 852 853 void G4InitXscPAI::IntegralPlasmon(G4double bg 853 void G4InitXscPAI::IntegralPlasmon(G4double bg2, G4double Tmax) 854 { 854 { 855 G4int i,k,i1,i2; 855 G4int i,k,i1,i2; 856 G4double energy1, energy2, result = 0.; 856 G4double energy1, energy2, result = 0.; 857 857 858 fBetaGammaSq = bg2; 858 fBetaGammaSq = bg2; 859 fTmax = Tmax; 859 fTmax = Tmax; 860 860 861 delete fPAIelectronVector; << 861 if(fPAIelectronVector) delete fPAIelectronVector; 862 862 863 fPAIelectronVector = new G4PhysicsLogVector( 863 fPAIelectronVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin); 864 fPAIelectronVector->PutValue(fPAIbin-1,resul 864 fPAIelectronVector->PutValue(fPAIbin-1,result); 865 865 866 for( i = fIntervalNumber - 1; i >= 0; i-- ) 866 for( i = fIntervalNumber - 1; i >= 0; i-- ) 867 { 867 { 868 if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) 868 if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break; 869 } 869 } 870 if (i < 0) i = 0; // Tmax should be more tha 870 if (i < 0) i = 0; // Tmax should be more than 871 // first ionisation potent 871 // first ionisation potential 872 fIntervalTmax = i; 872 fIntervalTmax = i; 873 873 874 G4Integrator<G4InitXscPAI,G4double(G4InitXsc 874 G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral; 875 875 876 for( k = fPAIbin - 2; k >= 0; k-- ) 876 for( k = fPAIbin - 2; k >= 0; k-- ) 877 { 877 { 878 energy1 = fPAIelectronVector->GetLowEdgeEn 878 energy1 = fPAIelectronVector->GetLowEdgeEnergy(k); 879 energy2 = fPAIelectronVector->GetLowEdgeEn 879 energy2 = fPAIelectronVector->GetLowEdgeEnergy(k+1); 880 880 881 for( i = fIntervalTmax; i >= 0; i-- ) 881 for( i = fIntervalTmax; i >= 0; i-- ) 882 { 882 { 883 if( energy2 > (*(*fMatSandiaMatrix)[i])[ 883 if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break; 884 } 884 } 885 if(i < 0) i = 0; 885 if(i < 0) i = 0; 886 i2 = i; 886 i2 = i; 887 887 888 for( i = fIntervalTmax; i >= 0; i-- ) 888 for( i = fIntervalTmax; i >= 0; i-- ) 889 { 889 { 890 if( energy1 > (*(*fMatSandiaMatrix)[i])[ 890 if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break; 891 } 891 } 892 if(i < 0) i = 0; 892 if(i < 0) i = 0; 893 i1 = i; 893 i1 = i; 894 894 895 if( i1 == i2 ) 895 if( i1 == i2 ) 896 { 896 { 897 fCurrentInterval = i1; 897 fCurrentInterval = i1; 898 result += integral.Legendre10(this,&G4In 898 result += integral.Legendre10(this,&G4InitXscPAI::PAIdNdxPlasmon, 899 energy1,en 899 energy1,energy2); 900 fPAIelectronVector->PutValue(k,result); 900 fPAIelectronVector->PutValue(k,result); 901 } 901 } 902 else 902 else 903 { 903 { 904 for( i = i2; i >= i1; i-- ) 904 for( i = i2; i >= i1; i-- ) 905 { 905 { 906 fCurrentInterval = i; 906 fCurrentInterval = i; 907 907 908 if( i==i2 ) result += integral. 908 if( i==i2 ) result += integral.Legendre10(this, 909 &G4InitXscPAI::PAId 909 &G4InitXscPAI::PAIdNdxPlasmon, 910 (*(*fMatSandiaMatri 910 (*(*fMatSandiaMatrix)[i])[0] ,energy2); 911 911 912 else if( i == i1 ) result += integral.Legend 912 else if( i == i1 ) result += integral.Legendre10(this, 913 &G4InitXscPAI::PAId 913 &G4InitXscPAI::PAIdNdxPlasmon,energy1, 914 (*(*fMatSandiaMatri 914 (*(*fMatSandiaMatrix)[i+1])[0]); 915 915 916 else result += integral. 916 else result += integral.Legendre10(this, 917 &G4InitXscPAI::PAId 917 &G4InitXscPAI::PAIdNdxPlasmon, 918 (*(*fMatSandiaMatrix)[i 918 (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]); 919 } 919 } 920 fPAIelectronVector->PutValue(k,result); 920 fPAIelectronVector->PutValue(k,result); 921 } 921 } 922 // G4cout<<k<<"\t"<<result<<G4endl; 922 // G4cout<<k<<"\t"<<result<<G4endl; 923 } 923 } 924 return; 924 return; 925 } // end of IntegralPlasmon 925 } // end of IntegralPlasmon 926 926 927 927 928 ////////////////////////////////////////////// 928 ///////////////////////////////////////////////////////////////////////// 929 // 929 // 930 // 930 // 931 931 932 G4double G4InitXscPAI::GetPhotonLambda( G4doub 932 G4double G4InitXscPAI::GetPhotonLambda( G4double omega ) 933 { 933 { 934 G4int i ; 934 G4int i ; 935 G4double omega2, omega3, omega4, a1, a2, a3, 935 G4double omega2, omega3, omega4, a1, a2, a3, a4, lambda ; 936 936 937 omega2 = omega*omega ; 937 omega2 = omega*omega ; 938 omega3 = omega2*omega ; 938 omega3 = omega2*omega ; 939 omega4 = omega2*omega2 ; 939 omega4 = omega2*omega2 ; 940 940 941 for(i = 0; i < fIntervalNumber;i++) 941 for(i = 0; i < fIntervalNumber;i++) 942 { 942 { 943 if( omega < (*(*fMatSandiaMatrix)[i])[0] ) 943 if( omega < (*(*fMatSandiaMatrix)[i])[0] ) break ; 944 } 944 } 945 if( i == 0 ) 945 if( i == 0 ) 946 { 946 { 947 G4cout<<"Warning: energy in G4InitXscPAI:: 947 G4cout<<"Warning: energy in G4InitXscPAI::GetPhotonLambda < I1"<<G4endl; 948 } 948 } 949 else i-- ; 949 else i-- ; 950 950 951 a1 = (*(*fMatSandiaMatrix)[i])[1]; 951 a1 = (*(*fMatSandiaMatrix)[i])[1]; 952 a2 = (*(*fMatSandiaMatrix)[i])[2]; 952 a2 = (*(*fMatSandiaMatrix)[i])[2]; 953 a3 = (*(*fMatSandiaMatrix)[i])[3]; 953 a3 = (*(*fMatSandiaMatrix)[i])[3]; 954 a4 = (*(*fMatSandiaMatrix)[i])[4]; 954 a4 = (*(*fMatSandiaMatrix)[i])[4]; 955 955 956 lambda = 1./(a1/omega + a2/omega2 + a3/omega 956 lambda = 1./(a1/omega + a2/omega2 + a3/omega3 + a4/omega4); 957 957 958 return lambda ; 958 return lambda ; 959 } 959 } 960 960 961 ////////////////////////////////////////////// 961 ///////////////////////////////////////////////////////////////////////// 962 // 962 // 963 // 963 // 964 964 965 ////////////////////////////////////////////// 965 ///////////////////////////////////////////////////////////////////////// 966 // 966 // 967 // 967 // 968 968 969 G4double G4InitXscPAI::GetStepEnergyLoss( G4do 969 G4double G4InitXscPAI::GetStepEnergyLoss( G4double step ) 970 { 970 { 971 G4double loss = 0.0 ; 971 G4double loss = 0.0 ; 972 loss *= step; 972 loss *= step; 973 973 974 return loss ; 974 return loss ; 975 } 975 } 976 976 977 ////////////////////////////////////////////// 977 ///////////////////////////////////////////////////////////////////////// 978 // 978 // 979 // 979 // 980 980 981 G4double G4InitXscPAI::GetStepCerenkovLoss( G4 981 G4double G4InitXscPAI::GetStepCerenkovLoss( G4double step ) 982 { 982 { 983 G4double loss = 0.0 ; 983 G4double loss = 0.0 ; 984 loss *= step; 984 loss *= step; 985 985 986 return loss ; 986 return loss ; 987 } 987 } 988 988 989 ////////////////////////////////////////////// 989 ///////////////////////////////////////////////////////////////////////// 990 // 990 // 991 // 991 // 992 992 993 G4double G4InitXscPAI::GetStepPlasmonLoss( G4d 993 G4double G4InitXscPAI::GetStepPlasmonLoss( G4double step ) 994 { 994 { 995 995 996 996 997 G4double loss = 0.0 ; 997 G4double loss = 0.0 ; 998 loss *= step; 998 loss *= step; 999 return loss ; 999 return loss ; 1000 } 1000 } 1001 1001 1002 1002 1003 // 1003 // 1004 // end of G4InitXscPAI implementation file 1004 // end of G4InitXscPAI implementation file 1005 // 1005 // 1006 ///////////////////////////////////////////// 1006 //////////////////////////////////////////////////////////////////////////// 1007 1007 1008 1008