Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // >> 26 // $Id: G4PenelopeIonisationXSHandler.cc 99415 2016-09-21 09:05:43Z gcosmo $ 26 // 27 // 27 // Author: Luciano Pandola 28 // Author: Luciano Pandola 28 // 29 // 29 // History: 30 // History: 30 // -------- 31 // -------- 31 // 08 Mar 2012 L Pandola First complete i 32 // 08 Mar 2012 L Pandola First complete implementation 32 // 33 // 33 34 34 #include "G4PenelopeIonisationXSHandler.hh" 35 #include "G4PenelopeIonisationXSHandler.hh" 35 #include "G4PhysicalConstants.hh" 36 #include "G4PhysicalConstants.hh" 36 #include "G4SystemOfUnits.hh" 37 #include "G4SystemOfUnits.hh" 37 #include "G4ParticleDefinition.hh" 38 #include "G4ParticleDefinition.hh" 38 #include "G4Electron.hh" 39 #include "G4Electron.hh" 39 #include "G4Positron.hh" 40 #include "G4Positron.hh" 40 #include "G4PenelopeOscillatorManager.hh" 41 #include "G4PenelopeOscillatorManager.hh" 41 #include "G4PenelopeOscillator.hh" 42 #include "G4PenelopeOscillator.hh" 42 #include "G4PenelopeCrossSection.hh" 43 #include "G4PenelopeCrossSection.hh" 43 #include "G4PhysicsFreeVector.hh" 44 #include "G4PhysicsFreeVector.hh" 44 #include "G4PhysicsLogVector.hh" 45 #include "G4PhysicsLogVector.hh" 45 46 46 G4PenelopeIonisationXSHandler::G4PenelopeIonis 47 G4PenelopeIonisationXSHandler::G4PenelopeIonisationXSHandler(size_t nb) 47 :fXSTableElectron(nullptr),fXSTablePositron( << 48 :XSTableElectron(0),XSTablePositron(0), 48 fDeltaTable(nullptr),fEnergyGrid(nullptr) << 49 theDeltaTable(0),energyGrid(0) 49 { 50 { 50 fNBins = nb; << 51 nBins = nb; 51 G4double LowEnergyLimit = 100.0*eV; 52 G4double LowEnergyLimit = 100.0*eV; 52 G4double HighEnergyLimit = 100.0*GeV; 53 G4double HighEnergyLimit = 100.0*GeV; 53 fOscManager = G4PenelopeOscillatorManager::G << 54 oscManager = G4PenelopeOscillatorManager::GetOscillatorManager(); 54 fXSTableElectron = new << 55 XSTableElectron = new 55 std::map< std::pair<const G4Material*,G4do 56 std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>; 56 fXSTablePositron = new << 57 XSTablePositron = new 57 std::map< std::pair<const G4Material*,G4do 58 std::map< std::pair<const G4Material*,G4double>, G4PenelopeCrossSection*>; 58 59 59 fDeltaTable = new std::map<const G4Material* << 60 theDeltaTable = new std::map<const G4Material*,G4PhysicsFreeVector*>; 60 fEnergyGrid = new G4PhysicsLogVector(LowEner << 61 energyGrid = new G4PhysicsLogVector(LowEnergyLimit, 61 HighEnergyLimit, 62 HighEnergyLimit, 62 fNBins-1); //one hidden bin is a << 63 nBins-1); //one hidden bin is added 63 fVerboseLevel = 0; << 64 >> 65 verboseLevel = 0; 64 } 66 } 65 67 66 //....oooOO0OOooo........oooOO0OOooo........oo 68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 67 69 68 G4PenelopeIonisationXSHandler::~G4PenelopeIoni 70 G4PenelopeIonisationXSHandler::~G4PenelopeIonisationXSHandler() 69 { 71 { 70 if (fXSTableElectron) << 72 if (XSTableElectron) 71 { 73 { 72 for (auto& item : (*fXSTableElectron)) << 74 for (auto& item : (*XSTableElectron)) 73 { 75 { 74 //G4PenelopeCrossSection* tab = i->second; 76 //G4PenelopeCrossSection* tab = i->second; 75 delete item.second; 77 delete item.second; 76 } 78 } 77 delete fXSTableElectron; << 79 delete XSTableElectron; 78 fXSTableElectron = nullptr; << 80 XSTableElectron = nullptr; 79 } 81 } 80 82 81 if (fXSTablePositron) << 83 if (XSTablePositron) 82 { 84 { 83 for (auto& item : (*fXSTablePositron)) << 85 for (auto& item : (*XSTablePositron)) 84 { 86 { 85 //G4PenelopeCrossSection* tab = i->second; 87 //G4PenelopeCrossSection* tab = i->second; 86 delete item.second; 88 delete item.second; 87 } 89 } 88 delete fXSTablePositron; << 90 delete XSTablePositron; 89 fXSTablePositron = nullptr; << 91 XSTablePositron = nullptr; 90 } 92 } 91 if (fDeltaTable) << 93 if (theDeltaTable) 92 { 94 { 93 for (auto& item : (*fDeltaTable)) << 95 for (auto& item : (*theDeltaTable)) 94 delete item.second; 96 delete item.second; 95 delete fDeltaTable; << 97 delete theDeltaTable; 96 fDeltaTable = nullptr; << 98 theDeltaTable = nullptr; 97 } 99 } 98 if (fEnergyGrid) << 100 if (energyGrid) 99 delete fEnergyGrid; << 101 delete energyGrid; 100 102 101 if (fVerboseLevel > 2) << 103 if (verboseLevel > 2) 102 G4cout << "G4PenelopeIonisationXSHandler. 104 G4cout << "G4PenelopeIonisationXSHandler. Tables have been cleared" 103 << G4endl; 105 << G4endl; 104 } 106 } 105 107 106 //....oooOO0OOooo........oooOO0OOooo........oo 108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 107 109 108 const G4PenelopeCrossSection* 110 const G4PenelopeCrossSection* 109 G4PenelopeIonisationXSHandler::GetCrossSection 111 G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple(const G4ParticleDefinition* part, 110 const G4Material* mat, 112 const G4Material* mat, 111 const G4double cut) const 113 const G4double cut) const 112 { 114 { 113 if (part != G4Electron::Electron() && part ! 115 if (part != G4Electron::Electron() && part != G4Positron::Positron()) 114 { 116 { 115 G4ExceptionDescription ed; 117 G4ExceptionDescription ed; 116 ed << "Invalid particle: " << part->GetP 118 ed << "Invalid particle: " << part->GetParticleName() << G4endl; 117 G4Exception("G4PenelopeIonisationXSHandl 119 G4Exception("G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple()", 118 "em0001",FatalException,ed); 120 "em0001",FatalException,ed); 119 return nullptr; 121 return nullptr; 120 } 122 } 121 123 122 if (part == G4Electron::Electron()) 124 if (part == G4Electron::Electron()) 123 { 125 { 124 if (!fXSTableElectron) << 126 if (!XSTableElectron) 125 { 127 { 126 G4Exception("G4PenelopeIonisationXSHandler 128 G4Exception("G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple()", 127 "em0028",FatalException, 129 "em0028",FatalException, 128 "The Cross Section Table for e- was 130 "The Cross Section Table for e- was not initialized correctly!"); 129 return nullptr; 131 return nullptr; 130 } 132 } 131 std::pair<const G4Material*,G4double> th 133 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut); 132 if (fXSTableElectron->count(theKey)) //t << 134 if (XSTableElectron->count(theKey)) //table already built 133 return fXSTableElectron->find(theKey)->secon << 135 return XSTableElectron->find(theKey)->second; 134 else 136 else 135 return nullptr; 137 return nullptr; 136 } 138 } 137 139 138 if (part == G4Positron::Positron()) 140 if (part == G4Positron::Positron()) 139 { 141 { 140 if (!fXSTablePositron) << 142 if (!XSTablePositron) 141 { 143 { 142 G4Exception("G4PenelopeIonisationXSHandler 144 G4Exception("G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple()", 143 "em0028",FatalException, 145 "em0028",FatalException, 144 "The Cross Section Table for e+ was 146 "The Cross Section Table for e+ was not initialized correctly!"); 145 return nullptr; 147 return nullptr; 146 } 148 } 147 std::pair<const G4Material*,G4double> th 149 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut); 148 if (fXSTablePositron->count(theKey)) //t << 150 if (XSTablePositron->count(theKey)) //table already built 149 return fXSTablePositron->find(theKey)->secon << 151 return XSTablePositron->find(theKey)->second; 150 else 152 else 151 return nullptr; 153 return nullptr; 152 } 154 } 153 return nullptr; 155 return nullptr; 154 } 156 } 155 157 156 //....oooOO0OOooo........oooOO0OOooo........oo 158 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 157 159 158 void G4PenelopeIonisationXSHandler::BuildXSTab 160 void G4PenelopeIonisationXSHandler::BuildXSTable(const G4Material* mat,G4double cut, 159 const G4ParticleDefinition* part, 161 const G4ParticleDefinition* part, 160 G4bool isMaster) 162 G4bool isMaster) 161 { 163 { 162 //Just to check 164 //Just to check 163 if (!isMaster) 165 if (!isMaster) 164 G4Exception("G4PenelopeIonisationXSHandler 166 G4Exception("G4PenelopeIonisationXSHandler::BuildXSTable()", 165 "em0100",FatalException,"Worke 167 "em0100",FatalException,"Worker thread in this method"); 166 168 167 // 169 // 168 //This method fills the G4PenelopeCrossSecti 170 //This method fills the G4PenelopeCrossSection containers for electrons or positrons 169 //and for the given material/cut couple. The 171 //and for the given material/cut couple. The calculation is done as sum over the 170 //individual shells. 172 //individual shells. 171 //Equivalent of subroutines EINaT and PINaT 173 //Equivalent of subroutines EINaT and PINaT of Penelope 172 // 174 // 173 if (fVerboseLevel > 2) << 175 if (verboseLevel > 2) 174 { 176 { 175 G4cout << "G4PenelopeIonisationXSHandler 177 G4cout << "G4PenelopeIonisationXSHandler: going to build cross section table " << G4endl; 176 G4cout << "for " << part->GetParticleNam 178 G4cout << "for " << part->GetParticleName() << " in " << mat->GetName() << G4endl; 177 G4cout << "Cut= " << cut/keV << " keV" < 179 G4cout << "Cut= " << cut/keV << " keV" << G4endl; 178 } 180 } 179 181 180 std::pair<const G4Material*,G4double> theKey 182 std::pair<const G4Material*,G4double> theKey = std::make_pair(mat,cut); 181 //Check if the table already exists 183 //Check if the table already exists 182 if (part == G4Electron::Electron()) 184 if (part == G4Electron::Electron()) 183 { 185 { 184 if (fXSTableElectron->count(theKey)) //t << 186 if (XSTableElectron->count(theKey)) //table already built 185 return; 187 return; 186 } 188 } 187 if (part == G4Positron::Positron()) 189 if (part == G4Positron::Positron()) 188 { 190 { 189 if (fXSTablePositron->count(theKey)) //t << 191 if (XSTablePositron->count(theKey)) //table already built 190 return; 192 return; 191 } 193 } 192 194 193 //check if the material has been built 195 //check if the material has been built 194 if (!(fDeltaTable->count(mat))) << 196 if (!(theDeltaTable->count(mat))) 195 BuildDeltaTable(mat); 197 BuildDeltaTable(mat); 196 198 197 199 198 //Tables have been already created (checked 200 //Tables have been already created (checked by GetCrossSectionTableForCouple) 199 G4PenelopeOscillatorTable* theTable = fOscMa << 201 G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableIonisation(mat); 200 size_t numberOfOscillators = theTable->size( 202 size_t numberOfOscillators = theTable->size(); 201 203 202 if (fEnergyGrid->GetVectorLength() != fNBins << 204 if (energyGrid->GetVectorLength() != nBins) 203 { 205 { 204 G4ExceptionDescription ed; 206 G4ExceptionDescription ed; 205 ed << "Energy Grid looks not initialized 207 ed << "Energy Grid looks not initialized" << G4endl; 206 ed << fNBins << " " << fEnergyGrid->GetV << 208 ed << nBins << " " << energyGrid->GetVectorLength() << G4endl; 207 G4Exception("G4PenelopeIonisationXSHandl 209 G4Exception("G4PenelopeIonisationXSHandler::BuildXSTable()", 208 "em2030",FatalException,ed); 210 "em2030",FatalException,ed); 209 } 211 } 210 212 211 G4PenelopeCrossSection* XSEntry = new G4Pene << 213 G4PenelopeCrossSection* XSEntry = new G4PenelopeCrossSection(nBins,numberOfOscillators); 212 214 213 //loop on the energy grid 215 //loop on the energy grid 214 for (size_t bin=0;bin<fNBins;bin++) << 216 for (size_t bin=0;bin<nBins;bin++) 215 { 217 { 216 G4double energy = fEnergyGrid->GetLowEd << 218 G4double energy = energyGrid->GetLowEdgeEnergy(bin); 217 G4double XH0=0, XH1=0, XH2=0; 219 G4double XH0=0, XH1=0, XH2=0; 218 G4double XS0=0, XS1=0, XS2=0; 220 G4double XS0=0, XS1=0, XS2=0; 219 221 220 //oscillator loop 222 //oscillator loop 221 for (size_t iosc=0;iosc<numberOfOscilla 223 for (size_t iosc=0;iosc<numberOfOscillators;iosc++) 222 { 224 { 223 G4DataVector* tempStorage = nullptr; << 225 G4DataVector* tempStorage = 0; 224 226 225 G4PenelopeOscillator* theOsc = (*theTable 227 G4PenelopeOscillator* theOsc = (*theTable)[iosc]; 226 G4double delta = GetDensityCorrection(mat 228 G4double delta = GetDensityCorrection(mat,energy); 227 if (part == G4Electron::Electron()) 229 if (part == G4Electron::Electron()) 228 tempStorage = ComputeShellCrossSections 230 tempStorage = ComputeShellCrossSectionsElectron(theOsc,energy,cut,delta); 229 else if (part == G4Positron::Positron()) 231 else if (part == G4Positron::Positron()) 230 tempStorage = ComputeShellCrossSections 232 tempStorage = ComputeShellCrossSectionsPositron(theOsc,energy,cut,delta); 231 //check results are all right 233 //check results are all right 232 if (!tempStorage) 234 if (!tempStorage) 233 { 235 { 234 G4ExceptionDescription ed; 236 G4ExceptionDescription ed; 235 ed << "Problem in calculating the she 237 ed << "Problem in calculating the shell XS for shell # " 236 << iosc << G4endl; 238 << iosc << G4endl; 237 G4Exception("G4PenelopeIonisationXSHa 239 G4Exception("G4PenelopeIonisationXSHandler::BuildXSTable()", 238 "em2031",FatalException,ed); 240 "em2031",FatalException,ed); 239 delete XSEntry; 241 delete XSEntry; 240 return; 242 return; 241 } 243 } 242 if (tempStorage->size() != 6) 244 if (tempStorage->size() != 6) 243 { 245 { 244 G4ExceptionDescription ed; 246 G4ExceptionDescription ed; 245 ed << "Problem in calculating the she 247 ed << "Problem in calculating the shell XS " << G4endl; 246 ed << "Result has dimension " << temp 248 ed << "Result has dimension " << tempStorage->size() << " instead of 6" << G4endl; 247 G4Exception("G4PenelopeIonisationXSHa 249 G4Exception("G4PenelopeIonisationXSHandler::BuildXSTable()", 248 "em2031",FatalException,ed); 250 "em2031",FatalException,ed); 249 } 251 } 250 G4double stre = theOsc->GetOscillatorStre 252 G4double stre = theOsc->GetOscillatorStrength(); 251 253 252 XH0 += stre*(*tempStorage)[0]; 254 XH0 += stre*(*tempStorage)[0]; 253 XH1 += stre*(*tempStorage)[1]; 255 XH1 += stre*(*tempStorage)[1]; 254 XH2 += stre*(*tempStorage)[2]; 256 XH2 += stre*(*tempStorage)[2]; 255 XS0 += stre*(*tempStorage)[3]; 257 XS0 += stre*(*tempStorage)[3]; 256 XS1 += stre*(*tempStorage)[4]; 258 XS1 += stre*(*tempStorage)[4]; 257 XS2 += stre*(*tempStorage)[5]; 259 XS2 += stre*(*tempStorage)[5]; 258 XSEntry->AddShellCrossSectionPoint(bin,io 260 XSEntry->AddShellCrossSectionPoint(bin,iosc,energy,stre*(*tempStorage)[0]); 259 if (tempStorage) 261 if (tempStorage) 260 { 262 { 261 delete tempStorage; 263 delete tempStorage; 262 tempStorage = nullptr; << 264 tempStorage = 0; 263 } 265 } 264 } 266 } 265 XSEntry->AddCrossSectionPoint(bin,energ 267 XSEntry->AddCrossSectionPoint(bin,energy,XH0,XH1,XH2,XS0,XS1,XS2); 266 } 268 } 267 //Do (only once) the final normalization 269 //Do (only once) the final normalization 268 XSEntry->NormalizeShellCrossSections(); 270 XSEntry->NormalizeShellCrossSections(); 269 271 270 //Insert in the appropriate table 272 //Insert in the appropriate table 271 if (part == G4Electron::Electron()) 273 if (part == G4Electron::Electron()) 272 fXSTableElectron->insert(std::make_pair(th << 274 XSTableElectron->insert(std::make_pair(theKey,XSEntry)); 273 else if (part == G4Positron::Positron()) 275 else if (part == G4Positron::Positron()) 274 fXSTablePositron->insert(std::make_pair(th << 276 XSTablePositron->insert(std::make_pair(theKey,XSEntry)); 275 else 277 else 276 delete XSEntry; 278 delete XSEntry; 277 279 278 return; 280 return; 279 } 281 } 280 282 281 283 282 //....oooOO0OOooo........oooOO0OOooo........oo 284 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 283 285 284 G4double G4PenelopeIonisationXSHandler::GetDen 286 G4double G4PenelopeIonisationXSHandler::GetDensityCorrection(const G4Material* mat, 285 const G4double energy) cons 287 const G4double energy) const 286 { 288 { 287 G4double result = 0; 289 G4double result = 0; 288 if (!fDeltaTable) << 290 if (!theDeltaTable) 289 { 291 { 290 G4Exception("G4PenelopeIonisationXSHandl 292 G4Exception("G4PenelopeIonisationXSHandler::GetDensityCorrection()", 291 "em2032",FatalException, 293 "em2032",FatalException, 292 "Delta Table not initialized. Was Initia 294 "Delta Table not initialized. Was Initialise() run?"); 293 return 0; 295 return 0; 294 } 296 } 295 if (energy <= 0*eV) 297 if (energy <= 0*eV) 296 { 298 { 297 G4cout << "G4PenelopeIonisationXSHandler 299 G4cout << "G4PenelopeIonisationXSHandler::GetDensityCorrection()" << G4endl; 298 G4cout << "Invalid energy " << energy/eV 300 G4cout << "Invalid energy " << energy/eV << " eV " << G4endl; 299 return 0; 301 return 0; 300 } 302 } 301 G4double logene = G4Log(energy); << 303 G4double logene = std::log(energy); 302 304 303 if (fDeltaTable->count(mat)) << 305 if (theDeltaTable->count(mat)) 304 { 306 { 305 const G4PhysicsFreeVector* vec = fDeltaT << 307 const G4PhysicsFreeVector* vec = theDeltaTable->find(mat)->second; 306 result = vec->Value(logene); //the table 308 result = vec->Value(logene); //the table has delta vs. ln(E) 307 } 309 } 308 else 310 else 309 { 311 { 310 G4ExceptionDescription ed; 312 G4ExceptionDescription ed; 311 ed << "Unable to build table for " << ma 313 ed << "Unable to build table for " << mat->GetName() << G4endl; 312 G4Exception("G4PenelopeIonisationXSHandl 314 G4Exception("G4PenelopeIonisationXSHandler::GetDensityCorrection()", 313 "em2033",FatalException,ed); 315 "em2033",FatalException,ed); 314 } 316 } 315 317 316 return result; 318 return result; 317 } 319 } 318 320 319 //....oooOO0OOooo........oooOO0OOooo........oo 321 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 320 322 321 void G4PenelopeIonisationXSHandler::BuildDelta 323 void G4PenelopeIonisationXSHandler::BuildDeltaTable(const G4Material* mat) 322 { 324 { 323 G4PenelopeOscillatorTable* theTable = fOscMa << 325 G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableIonisation(mat); 324 G4double plasmaSq = fOscManager->GetPlasmaEn << 326 G4double plasmaSq = oscManager->GetPlasmaEnergySquared(mat); 325 G4double totalZ = fOscManager->GetTotalZ(mat << 327 G4double totalZ = oscManager->GetTotalZ(mat); 326 size_t numberOfOscillators = theTable->size( 328 size_t numberOfOscillators = theTable->size(); 327 329 328 if (fEnergyGrid->GetVectorLength() != fNBins << 330 if (energyGrid->GetVectorLength() != nBins) 329 { 331 { 330 G4ExceptionDescription ed; 332 G4ExceptionDescription ed; 331 ed << "Energy Grid for Delta table looks 333 ed << "Energy Grid for Delta table looks not initialized" << G4endl; 332 ed << fNBins << " " << fEnergyGrid->GetV << 334 ed << nBins << " " << energyGrid->GetVectorLength() << G4endl; 333 G4Exception("G4PenelopeIonisationXSHandl 335 G4Exception("G4PenelopeIonisationXSHandler::BuildDeltaTable()", 334 "em2030",FatalException,ed); 336 "em2030",FatalException,ed); 335 } 337 } 336 338 337 G4PhysicsFreeVector* theVector = new G4Physi << 339 G4PhysicsFreeVector* theVector = new G4PhysicsFreeVector(nBins); 338 340 339 //loop on the energy grid 341 //loop on the energy grid 340 for (size_t bin=0;bin<fNBins;bin++) << 342 for (size_t bin=0;bin<nBins;bin++) 341 { 343 { 342 G4double delta = 0.; 344 G4double delta = 0.; 343 G4double energy = fEnergyGrid->GetLowEdg << 345 G4double energy = energyGrid->GetLowEdgeEnergy(bin); 344 346 345 //Here calculate delta 347 //Here calculate delta 346 G4double gam = 1.0+(energy/electron_mass 348 G4double gam = 1.0+(energy/electron_mass_c2); 347 G4double gamSq = gam*gam; 349 G4double gamSq = gam*gam; 348 350 349 G4double TST = totalZ/(gamSq*plasmaSq); 351 G4double TST = totalZ/(gamSq*plasmaSq); 350 G4double wl2 = 0; 352 G4double wl2 = 0; 351 G4double fdel = 0; 353 G4double fdel = 0; 352 354 353 //loop on oscillators 355 //loop on oscillators 354 for (size_t i=0;i<numberOfOscillators;i+ 356 for (size_t i=0;i<numberOfOscillators;i++) 355 { 357 { 356 G4PenelopeOscillator* theOsc = (*theTable) 358 G4PenelopeOscillator* theOsc = (*theTable)[i]; 357 G4double wri = theOsc->GetResonanceEnergy( 359 G4double wri = theOsc->GetResonanceEnergy(); 358 fdel += theOsc->GetOscillatorStrength()/(w 360 fdel += theOsc->GetOscillatorStrength()/(wri*wri+wl2); 359 } 361 } 360 if (fdel >= TST) //if fdel < TST, delta 362 if (fdel >= TST) //if fdel < TST, delta = 0 361 { 363 { 362 //get last oscillator 364 //get last oscillator 363 G4PenelopeOscillator* theOsc = (*theTable) 365 G4PenelopeOscillator* theOsc = (*theTable)[numberOfOscillators-1]; 364 wl2 = theOsc->GetResonanceEnergy()*theOsc- 366 wl2 = theOsc->GetResonanceEnergy()*theOsc->GetResonanceEnergy(); 365 367 366 //First iteration 368 //First iteration 367 G4bool loopAgain = false; 369 G4bool loopAgain = false; 368 do 370 do 369 { 371 { 370 loopAgain = false; 372 loopAgain = false; 371 wl2 += wl2; 373 wl2 += wl2; 372 fdel = 0.; 374 fdel = 0.; 373 for (size_t i=0;i<numberOfOscillators; 375 for (size_t i=0;i<numberOfOscillators;i++) 374 { 376 { 375 G4PenelopeOscillator* theOscLocal1 = (*t 377 G4PenelopeOscillator* theOscLocal1 = (*theTable)[i]; 376 G4double wri = theOscLocal1->GetResonanc 378 G4double wri = theOscLocal1->GetResonanceEnergy(); 377 fdel += theOscLocal1->GetOscillatorStren 379 fdel += theOscLocal1->GetOscillatorStrength()/(wri*wri+wl2); 378 } 380 } 379 if (fdel > TST) 381 if (fdel > TST) 380 loopAgain = true; 382 loopAgain = true; 381 }while(loopAgain); 383 }while(loopAgain); 382 384 383 G4double wl2l = 0; 385 G4double wl2l = 0; 384 G4double wl2u = wl2; 386 G4double wl2u = wl2; 385 //second iteration 387 //second iteration 386 do 388 do 387 { 389 { 388 loopAgain = false; 390 loopAgain = false; 389 wl2 = 0.5*(wl2l+wl2u); 391 wl2 = 0.5*(wl2l+wl2u); 390 fdel = 0; 392 fdel = 0; 391 for (size_t i=0;i<numberOfOscillators; 393 for (size_t i=0;i<numberOfOscillators;i++) 392 { 394 { 393 G4PenelopeOscillator* theOscLocal2 = (*t 395 G4PenelopeOscillator* theOscLocal2 = (*theTable)[i]; 394 G4double wri = theOscLocal2->GetResonanc 396 G4double wri = theOscLocal2->GetResonanceEnergy(); 395 fdel += theOscLocal2->GetOscillatorStren 397 fdel += theOscLocal2->GetOscillatorStrength()/(wri*wri+wl2); 396 } 398 } 397 if (fdel > TST) 399 if (fdel > TST) 398 wl2l = wl2; 400 wl2l = wl2; 399 else 401 else 400 wl2u = wl2; 402 wl2u = wl2; 401 if ((wl2u-wl2l)>1e-12*wl2) 403 if ((wl2u-wl2l)>1e-12*wl2) 402 loopAgain = true; 404 loopAgain = true; 403 }while(loopAgain); 405 }while(loopAgain); 404 406 405 //Eventually get density correction 407 //Eventually get density correction 406 delta = 0.; 408 delta = 0.; 407 for (size_t i=0;i<numberOfOscillators;i++) 409 for (size_t i=0;i<numberOfOscillators;i++) 408 { 410 { 409 G4PenelopeOscillator* theOscLocal3 = ( 411 G4PenelopeOscillator* theOscLocal3 = (*theTable)[i]; 410 G4double wri = theOscLocal3->GetResona 412 G4double wri = theOscLocal3->GetResonanceEnergy(); 411 delta += theOscLocal3->GetOscillatorSt 413 delta += theOscLocal3->GetOscillatorStrength()* 412 G4Log(1.0+(wl2/(wri*wri))); << 414 std::log(1.0+(wl2/(wri*wri))); 413 } 415 } 414 delta = (delta/totalZ)-wl2/(gamSq*plasmaSq 416 delta = (delta/totalZ)-wl2/(gamSq*plasmaSq); 415 } 417 } 416 energy = std::max(1e-9*eV,energy); //pre 418 energy = std::max(1e-9*eV,energy); //prevents log(0) 417 theVector->PutValue(bin,G4Log(energy),de << 419 theVector->PutValue(bin,std::log(energy),delta); 418 } 420 } 419 fDeltaTable->insert(std::make_pair(mat,theVe << 421 theDeltaTable->insert(std::make_pair(mat,theVector)); 420 return; 422 return; 421 } 423 } 422 424 423 //....oooOO0OOooo........oooOO0OOooo........oo 425 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 424 G4DataVector* G4PenelopeIonisationXSHandler::C 426 G4DataVector* G4PenelopeIonisationXSHandler::ComputeShellCrossSectionsElectron(G4PenelopeOscillator* theOsc, 425 G4double energy, 427 G4double energy, 426 G4double cut, 428 G4double cut, 427 G4double delta) 429 G4double delta) 428 { 430 { 429 // 431 // 430 //This method calculates the hard and soft c 432 //This method calculates the hard and soft cross sections (H0-H1-H2-S0-S1-S2) for 431 //the given oscillator/cut and at the given 433 //the given oscillator/cut and at the given energy. 432 //It returns a G4DataVector* with 6 entries 434 //It returns a G4DataVector* with 6 entries (H0-H1-H2-S0-S1-S2) 433 //Equivalent of subroutines EINaT1 of Penelo 435 //Equivalent of subroutines EINaT1 of Penelope 434 // 436 // 435 // Results are _per target electron_ 437 // Results are _per target electron_ 436 // 438 // 437 G4DataVector* result = new G4DataVector(); 439 G4DataVector* result = new G4DataVector(); 438 for (size_t i=0;i<6;i++) 440 for (size_t i=0;i<6;i++) 439 result->push_back(0.); 441 result->push_back(0.); 440 G4double ionEnergy = theOsc->GetIonisationEn 442 G4double ionEnergy = theOsc->GetIonisationEnergy(); 441 443 442 //return a set of zero's if the energy it to 444 //return a set of zero's if the energy it too low to excite the current oscillator 443 if (energy < ionEnergy) 445 if (energy < ionEnergy) 444 return result; 446 return result; 445 447 446 G4double H0=0.,H1=0.,H2=0.; 448 G4double H0=0.,H1=0.,H2=0.; 447 G4double S0=0.,S1=0.,S2=0.; 449 G4double S0=0.,S1=0.,S2=0.; 448 450 449 //Define useful constants to be used in the 451 //Define useful constants to be used in the calculation 450 G4double gamma = 1.0+energy/electron_mass_c2 452 G4double gamma = 1.0+energy/electron_mass_c2; 451 G4double gammaSq = gamma*gamma; 453 G4double gammaSq = gamma*gamma; 452 G4double beta = (gammaSq-1.0)/gammaSq; 454 G4double beta = (gammaSq-1.0)/gammaSq; 453 G4double pielr2 = pi*classic_electr_radius*c 455 G4double pielr2 = pi*classic_electr_radius*classic_electr_radius; //pi*re^2 454 G4double constant = pielr2*2.0*electron_mass 456 G4double constant = pielr2*2.0*electron_mass_c2/beta; 455 G4double XHDT0 = G4Log(gammaSq)-beta; << 457 G4double XHDT0 = std::log(gammaSq)-beta; 456 458 457 G4double cpSq = energy*(energy+2.0*electron_ 459 G4double cpSq = energy*(energy+2.0*electron_mass_c2); 458 G4double cp = std::sqrt(cpSq); 460 G4double cp = std::sqrt(cpSq); 459 G4double amol = (energy/(energy+electron_mas 461 G4double amol = (energy/(energy+electron_mass_c2))*(energy/(energy+electron_mass_c2)); 460 462 461 // 463 // 462 // Distant interactions 464 // Distant interactions 463 // 465 // 464 G4double resEne = theOsc->GetResonanceEnergy 466 G4double resEne = theOsc->GetResonanceEnergy(); 465 G4double cutoffEne = theOsc->GetCutoffRecoil 467 G4double cutoffEne = theOsc->GetCutoffRecoilResonantEnergy(); 466 if (energy > resEne) 468 if (energy > resEne) 467 { 469 { 468 G4double cp1Sq = (energy-resEne)*(energy 470 G4double cp1Sq = (energy-resEne)*(energy-resEne+2.0*electron_mass_c2); 469 G4double cp1 = std::sqrt(cp1Sq); 471 G4double cp1 = std::sqrt(cp1Sq); 470 472 471 //Distant longitudinal interactions 473 //Distant longitudinal interactions 472 G4double QM = 0; 474 G4double QM = 0; 473 if (resEne > 1e-6*energy) 475 if (resEne > 1e-6*energy) 474 QM = std::sqrt((cp-cp1)*(cp-cp1)+electron_ma 476 QM = std::sqrt((cp-cp1)*(cp-cp1)+electron_mass_c2*electron_mass_c2)-electron_mass_c2; 475 else 477 else 476 { 478 { 477 QM = resEne*resEne/(beta*2.0*electron_mass 479 QM = resEne*resEne/(beta*2.0*electron_mass_c2); 478 QM = QM*(1.0-0.5*QM/electron_mass_c2); 480 QM = QM*(1.0-0.5*QM/electron_mass_c2); 479 } 481 } 480 G4double SDL1 = 0; 482 G4double SDL1 = 0; 481 if (QM < cutoffEne) 483 if (QM < cutoffEne) 482 SDL1 = G4Log(cutoffEne*(QM+2.0*electron_mass << 484 SDL1 = std::log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2))); 483 485 484 //Distant transverse interactions 486 //Distant transverse interactions 485 if (SDL1) 487 if (SDL1) 486 { 488 { 487 G4double SDT1 = std::max(XHDT0-delta,0.0); 489 G4double SDT1 = std::max(XHDT0-delta,0.0); 488 G4double SD1 = SDL1+SDT1; 490 G4double SD1 = SDL1+SDT1; 489 if (cut > resEne) 491 if (cut > resEne) 490 { 492 { 491 S1 = SD1; //XS1 493 S1 = SD1; //XS1 492 S0 = SD1/resEne; //XS0 494 S0 = SD1/resEne; //XS0 493 S2 = SD1*resEne; //XS2 495 S2 = SD1*resEne; //XS2 494 } 496 } 495 else 497 else 496 { 498 { 497 H1 = SD1; //XH1 499 H1 = SD1; //XH1 498 H0 = SD1/resEne; //XH0 500 H0 = SD1/resEne; //XH0 499 H2 = SD1*resEne; //XH2 501 H2 = SD1*resEne; //XH2 500 } 502 } 501 } 503 } 502 } 504 } 503 // 505 // 504 // Close collisions (Moller's cross section) 506 // Close collisions (Moller's cross section) 505 // 507 // 506 G4double wl = std::max(cut,cutoffEne); 508 G4double wl = std::max(cut,cutoffEne); 507 G4double ee = energy + ionEnergy; 509 G4double ee = energy + ionEnergy; 508 G4double wu = 0.5*ee; 510 G4double wu = 0.5*ee; 509 if (wl < wu-(1e-5*eV)) 511 if (wl < wu-(1e-5*eV)) 510 { 512 { 511 H0 += (1.0/(ee-wu)) - (1.0/(ee-wl)) - (1 513 H0 += (1.0/(ee-wu)) - (1.0/(ee-wl)) - (1.0/wu) + (1.0/wl) + 512 (1.0-amol)*G4Log(((ee-wu)*wl)/((ee-wl)*wu))/ << 514 (1.0-amol)*std::log(((ee-wu)*wl)/((ee-wl)*wu))/ee + 513 amol*(wu-wl)/(ee*ee); 515 amol*(wu-wl)/(ee*ee); 514 H1 += G4Log(wu/wl)+(ee/(ee-wu))-(ee/(ee- << 516 H1 += std::log(wu/wl)+(ee/(ee-wu))-(ee/(ee-wl)) + 515 (2.0-amol)*G4Log((ee-wu)/(ee-wl)) + << 517 (2.0-amol)*std::log((ee-wu)/(ee-wl)) + 516 amol*(wu*wu-wl*wl)/(2.0*ee*ee); 518 amol*(wu*wu-wl*wl)/(2.0*ee*ee); 517 H2 += (2.0-amol)*(wu-wl)+(wu*(2.0*ee-wu) 519 H2 += (2.0-amol)*(wu-wl)+(wu*(2.0*ee-wu)/(ee-wu)) - 518 (wl*(2.0*ee-wl)/(ee-wl)) + 520 (wl*(2.0*ee-wl)/(ee-wl)) + 519 (3.0-amol)*ee*G4Log((ee-wu)/(ee-wl)) + << 521 (3.0-amol)*ee*std::log((ee-wu)/(ee-wl)) + 520 amol*(wu*wu*wu-wl*wl*wl)/(3.0*ee*ee); 522 amol*(wu*wu*wu-wl*wl*wl)/(3.0*ee*ee); 521 wu = wl; 523 wu = wl; 522 } 524 } 523 wl = cutoffEne; 525 wl = cutoffEne; 524 526 525 if (wl > wu-(1e-5*eV)) 527 if (wl > wu-(1e-5*eV)) 526 { 528 { 527 (*result)[0] = constant*H0; 529 (*result)[0] = constant*H0; 528 (*result)[1] = constant*H1; 530 (*result)[1] = constant*H1; 529 (*result)[2] = constant*H2; 531 (*result)[2] = constant*H2; 530 (*result)[3] = constant*S0; 532 (*result)[3] = constant*S0; 531 (*result)[4] = constant*S1; 533 (*result)[4] = constant*S1; 532 (*result)[5] = constant*S2; 534 (*result)[5] = constant*S2; 533 return result; 535 return result; 534 } 536 } 535 537 536 S0 += (1.0/(ee-wu))-(1.0/(ee-wl)) - (1.0/wu) 538 S0 += (1.0/(ee-wu))-(1.0/(ee-wl)) - (1.0/wu) + (1.0/wl) + 537 (1.0-amol)*G4Log(((ee-wu)*wl)/((ee-wl)*wu) << 539 (1.0-amol)*std::log(((ee-wu)*wl)/((ee-wl)*wu))/ee + 538 amol*(wu-wl)/(ee*ee); 540 amol*(wu-wl)/(ee*ee); 539 S1 += G4Log(wu/wl)+(ee/(ee-wu))-(ee/(ee-wl)) << 541 S1 += std::log(wu/wl)+(ee/(ee-wu))-(ee/(ee-wl)) + 540 (2.0-amol)*G4Log((ee-wu)/(ee-wl)) + << 542 (2.0-amol)*std::log((ee-wu)/(ee-wl)) + 541 amol*(wu*wu-wl*wl)/(2.0*ee*ee); 543 amol*(wu*wu-wl*wl)/(2.0*ee*ee); 542 S2 += (2.0-amol)*(wu-wl)+(wu*(2.0*ee-wu)/(ee 544 S2 += (2.0-amol)*(wu-wl)+(wu*(2.0*ee-wu)/(ee-wu)) - 543 (wl*(2.0*ee-wl)/(ee-wl)) + 545 (wl*(2.0*ee-wl)/(ee-wl)) + 544 (3.0-amol)*ee*G4Log((ee-wu)/(ee-wl)) + << 546 (3.0-amol)*ee*std::log((ee-wu)/(ee-wl)) + 545 amol*(wu*wu*wu-wl*wl*wl)/(3.0*ee*ee); 547 amol*(wu*wu*wu-wl*wl*wl)/(3.0*ee*ee); 546 548 547 (*result)[0] = constant*H0; 549 (*result)[0] = constant*H0; 548 (*result)[1] = constant*H1; 550 (*result)[1] = constant*H1; 549 (*result)[2] = constant*H2; 551 (*result)[2] = constant*H2; 550 (*result)[3] = constant*S0; 552 (*result)[3] = constant*S0; 551 (*result)[4] = constant*S1; 553 (*result)[4] = constant*S1; 552 (*result)[5] = constant*S2; 554 (*result)[5] = constant*S2; 553 return result; 555 return result; 554 } 556 } 555 //....oooOO0OOooo........oooOO0OOooo........oo 557 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 556 G4DataVector* G4PenelopeIonisationXSHandler::C 558 G4DataVector* G4PenelopeIonisationXSHandler::ComputeShellCrossSectionsPositron(G4PenelopeOscillator* theOsc, 557 G4double energy, 559 G4double energy, 558 G4double cut, 560 G4double cut, 559 G4double delta) 561 G4double delta) 560 { 562 { 561 // 563 // 562 //This method calculates the hard and soft c 564 //This method calculates the hard and soft cross sections (H0-H1-H2-S0-S1-S2) for 563 //the given oscillator/cut and at the given 565 //the given oscillator/cut and at the given energy. 564 //It returns a G4DataVector* with 6 entries 566 //It returns a G4DataVector* with 6 entries (H0-H1-H2-S0-S1-S2) 565 //Equivalent of subroutines PINaT1 of Penelo 567 //Equivalent of subroutines PINaT1 of Penelope 566 // 568 // 567 // Results are _per target electron_ 569 // Results are _per target electron_ 568 // 570 // 569 G4DataVector* result = new G4DataVector(); 571 G4DataVector* result = new G4DataVector(); 570 for (size_t i=0;i<6;i++) 572 for (size_t i=0;i<6;i++) 571 result->push_back(0.); 573 result->push_back(0.); 572 G4double ionEnergy = theOsc->GetIonisationEn 574 G4double ionEnergy = theOsc->GetIonisationEnergy(); 573 575 574 //return a set of zero's if the energy it to 576 //return a set of zero's if the energy it too low to excite the current oscillator 575 if (energy < ionEnergy) 577 if (energy < ionEnergy) 576 return result; 578 return result; 577 579 578 G4double H0=0.,H1=0.,H2=0.; 580 G4double H0=0.,H1=0.,H2=0.; 579 G4double S0=0.,S1=0.,S2=0.; 581 G4double S0=0.,S1=0.,S2=0.; 580 582 581 //Define useful constants to be used in the 583 //Define useful constants to be used in the calculation 582 G4double gamma = 1.0+energy/electron_mass_c2 584 G4double gamma = 1.0+energy/electron_mass_c2; 583 G4double gammaSq = gamma*gamma; 585 G4double gammaSq = gamma*gamma; 584 G4double beta = (gammaSq-1.0)/gammaSq; 586 G4double beta = (gammaSq-1.0)/gammaSq; 585 G4double pielr2 = pi*classic_electr_radius*c 587 G4double pielr2 = pi*classic_electr_radius*classic_electr_radius; //pi*re^2 586 G4double constant = pielr2*2.0*electron_mass 588 G4double constant = pielr2*2.0*electron_mass_c2/beta; 587 G4double XHDT0 = G4Log(gammaSq)-beta; << 589 G4double XHDT0 = std::log(gammaSq)-beta; 588 590 589 G4double cpSq = energy*(energy+2.0*electron_ 591 G4double cpSq = energy*(energy+2.0*electron_mass_c2); 590 G4double cp = std::sqrt(cpSq); 592 G4double cp = std::sqrt(cpSq); 591 G4double amol = (energy/(energy+electron_mas 593 G4double amol = (energy/(energy+electron_mass_c2))*(energy/(energy+electron_mass_c2)); 592 G4double g12 = (gamma+1.0)*(gamma+1.0); 594 G4double g12 = (gamma+1.0)*(gamma+1.0); 593 //Bhabha coefficients 595 //Bhabha coefficients 594 G4double bha1 = amol*(2.0*g12-1.0)/(gammaSq- 596 G4double bha1 = amol*(2.0*g12-1.0)/(gammaSq-1.0); 595 G4double bha2 = amol*(3.0+1.0/g12); 597 G4double bha2 = amol*(3.0+1.0/g12); 596 G4double bha3 = amol*2.0*gamma*(gamma-1.0)/g 598 G4double bha3 = amol*2.0*gamma*(gamma-1.0)/g12; 597 G4double bha4 = amol*(gamma-1.0)*(gamma-1.0) 599 G4double bha4 = amol*(gamma-1.0)*(gamma-1.0)/g12; 598 600 599 // 601 // 600 // Distant interactions 602 // Distant interactions 601 // 603 // 602 G4double resEne = theOsc->GetResonanceEnergy 604 G4double resEne = theOsc->GetResonanceEnergy(); 603 G4double cutoffEne = theOsc->GetCutoffRecoil 605 G4double cutoffEne = theOsc->GetCutoffRecoilResonantEnergy(); 604 if (energy > resEne) 606 if (energy > resEne) 605 { 607 { 606 G4double cp1Sq = (energy-resEne)*(energy 608 G4double cp1Sq = (energy-resEne)*(energy-resEne+2.0*electron_mass_c2); 607 G4double cp1 = std::sqrt(cp1Sq); 609 G4double cp1 = std::sqrt(cp1Sq); 608 610 609 //Distant longitudinal interactions 611 //Distant longitudinal interactions 610 G4double QM = 0; 612 G4double QM = 0; 611 if (resEne > 1e-6*energy) 613 if (resEne > 1e-6*energy) 612 QM = std::sqrt((cp-cp1)*(cp-cp1)+electron_ma 614 QM = std::sqrt((cp-cp1)*(cp-cp1)+electron_mass_c2*electron_mass_c2)-electron_mass_c2; 613 else 615 else 614 { 616 { 615 QM = resEne*resEne/(beta*2.0*electron_mass 617 QM = resEne*resEne/(beta*2.0*electron_mass_c2); 616 QM = QM*(1.0-0.5*QM/electron_mass_c2); 618 QM = QM*(1.0-0.5*QM/electron_mass_c2); 617 } 619 } 618 G4double SDL1 = 0; 620 G4double SDL1 = 0; 619 if (QM < cutoffEne) 621 if (QM < cutoffEne) 620 SDL1 = G4Log(cutoffEne*(QM+2.0*electron_mass << 622 SDL1 = std::log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2))); 621 623 622 //Distant transverse interactions 624 //Distant transverse interactions 623 if (SDL1) 625 if (SDL1) 624 { 626 { 625 G4double SDT1 = std::max(XHDT0-delta,0.0); 627 G4double SDT1 = std::max(XHDT0-delta,0.0); 626 G4double SD1 = SDL1+SDT1; 628 G4double SD1 = SDL1+SDT1; 627 if (cut > resEne) 629 if (cut > resEne) 628 { 630 { 629 S1 = SD1; //XS1 631 S1 = SD1; //XS1 630 S0 = SD1/resEne; //XS0 632 S0 = SD1/resEne; //XS0 631 S2 = SD1*resEne; //XS2 633 S2 = SD1*resEne; //XS2 632 } 634 } 633 else 635 else 634 { 636 { 635 H1 = SD1; //XH1 637 H1 = SD1; //XH1 636 H0 = SD1/resEne; //XH0 638 H0 = SD1/resEne; //XH0 637 H2 = SD1*resEne; //XH2 639 H2 = SD1*resEne; //XH2 638 } 640 } 639 } 641 } 640 } 642 } >> 643 641 // 644 // 642 // Close collisions (Bhabha's cross section) 645 // Close collisions (Bhabha's cross section) 643 // 646 // 644 G4double wl = std::max(cut,cutoffEne); 647 G4double wl = std::max(cut,cutoffEne); 645 G4double wu = energy; 648 G4double wu = energy; 646 G4double energySq = energy*energy; 649 G4double energySq = energy*energy; 647 if (wl < wu-(1e-5*eV)) 650 if (wl < wu-(1e-5*eV)) 648 { 651 { 649 G4double wlSq = wl*wl; 652 G4double wlSq = wl*wl; 650 G4double wuSq = wu*wu; 653 G4double wuSq = wu*wu; 651 H0 += (1.0/wl) - (1.0/wu)- bha1*G4Log(wu << 654 H0 += (1.0/wl) - (1.0/wu)- bha1*std::log(wu/wl)/energy 652 + bha2*(wu-wl)/energySq 655 + bha2*(wu-wl)/energySq 653 - bha3*(wuSq-wlSq)/(2.0*energySq*energy) 656 - bha3*(wuSq-wlSq)/(2.0*energySq*energy) 654 + bha4*(wuSq*wu-wlSq*wl)/(3.0*energySq*energ 657 + bha4*(wuSq*wu-wlSq*wl)/(3.0*energySq*energySq); 655 H1 += G4Log(wu/wl) - bha1*(wu-wl)/energy << 658 H1 += std::log(wu/wl) - bha1*(wu-wl)/energy 656 + bha2*(wuSq-wlSq)/(2.0*energySq) 659 + bha2*(wuSq-wlSq)/(2.0*energySq) 657 - bha3*(wuSq*wu-wlSq*wl)/(3.0*energySq*energ 660 - bha3*(wuSq*wu-wlSq*wl)/(3.0*energySq*energy) 658 + bha4*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*e 661 + bha4*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energySq); 659 H2 += wu - wl - bha1*(wuSq-wlSq)/(2.0*en 662 H2 += wu - wl - bha1*(wuSq-wlSq)/(2.0*energy) 660 + bha2*(wuSq*wu-wlSq*wl)/(3.0*energySq) 663 + bha2*(wuSq*wu-wlSq*wl)/(3.0*energySq) 661 - bha3*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*e 664 - bha3*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energy) 662 + bha4*(wuSq*wuSq*wu-wlSq*wlSq*wl)/(5.0*ener 665 + bha4*(wuSq*wuSq*wu-wlSq*wlSq*wl)/(5.0*energySq*energySq); 663 wu = wl; 666 wu = wl; 664 } 667 } 665 wl = cutoffEne; 668 wl = cutoffEne; 666 669 667 if (wl > wu-(1e-5*eV)) 670 if (wl > wu-(1e-5*eV)) 668 { 671 { 669 (*result)[0] = constant*H0; 672 (*result)[0] = constant*H0; 670 (*result)[1] = constant*H1; 673 (*result)[1] = constant*H1; 671 (*result)[2] = constant*H2; 674 (*result)[2] = constant*H2; 672 (*result)[3] = constant*S0; 675 (*result)[3] = constant*S0; 673 (*result)[4] = constant*S1; 676 (*result)[4] = constant*S1; 674 (*result)[5] = constant*S2; 677 (*result)[5] = constant*S2; 675 return result; 678 return result; 676 } 679 } 677 680 678 G4double wlSq = wl*wl; 681 G4double wlSq = wl*wl; 679 G4double wuSq = wu*wu; 682 G4double wuSq = wu*wu; 680 683 681 S0 += (1.0/wl) - (1.0/wu) - bha1*G4Log(wu/wl << 684 S0 += (1.0/wl) - (1.0/wu) - bha1*std::log(wu/wl)/energy 682 + bha2*(wu-wl)/energySq 685 + bha2*(wu-wl)/energySq 683 - bha3*(wuSq-wlSq)/(2.0*energySq*energy) 686 - bha3*(wuSq-wlSq)/(2.0*energySq*energy) 684 + bha4*(wuSq*wu-wlSq*wl)/(3.0*energySq*ene 687 + bha4*(wuSq*wu-wlSq*wl)/(3.0*energySq*energySq); 685 688 686 S1 += G4Log(wu/wl) - bha1*(wu-wl)/energy << 689 S1 += std::log(wu/wl) - bha1*(wu-wl)/energy 687 + bha2*(wuSq-wlSq)/(2.0*energySq) 690 + bha2*(wuSq-wlSq)/(2.0*energySq) 688 - bha3*(wuSq*wu-wlSq*wl)/(3.0*energySq*ene 691 - bha3*(wuSq*wu-wlSq*wl)/(3.0*energySq*energy) 689 + bha4*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq 692 + bha4*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energySq); 690 693 691 S2 += wu - wl - bha1*(wuSq-wlSq)/(2.0*energy 694 S2 += wu - wl - bha1*(wuSq-wlSq)/(2.0*energy) 692 + bha2*(wuSq*wu-wlSq*wl)/(3.0*energySq) 695 + bha2*(wuSq*wu-wlSq*wl)/(3.0*energySq) 693 - bha3*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq 696 - bha3*(wuSq*wuSq-wlSq*wlSq)/(4.0*energySq*energy) 694 + bha4*(wuSq*wuSq*wu-wlSq*wlSq*wl)/(5.0*en 697 + bha4*(wuSq*wuSq*wu-wlSq*wlSq*wl)/(5.0*energySq*energySq); 695 698 696 (*result)[0] = constant*H0; 699 (*result)[0] = constant*H0; 697 (*result)[1] = constant*H1; 700 (*result)[1] = constant*H1; 698 (*result)[2] = constant*H2; 701 (*result)[2] = constant*H2; 699 (*result)[3] = constant*S0; 702 (*result)[3] = constant*S0; 700 (*result)[4] = constant*S1; 703 (*result)[4] = constant*S1; 701 (*result)[5] = constant*S2; 704 (*result)[5] = constant*S2; 702 705 703 return result; 706 return result; 704 } 707 } 705 708