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 // 28 // 29 // File name: G4LowEnergyBremsstrahlung 29 // File name: G4LowEnergyBremsstrahlung 30 // 30 // 31 // Author: Alessandra Forti, Vladimir I 31 // Author: Alessandra Forti, Vladimir Ivanchenko 32 // 32 // 33 // Creation date: March 1999 33 // Creation date: March 1999 34 // 34 // 35 // Modifications: 35 // Modifications: 36 // 18.04.2000 V.L. 36 // 18.04.2000 V.L. 37 // - First implementation of continuous energ 37 // - First implementation of continuous energy loss. 38 // 17.02.2000 Veronique Lefebure 38 // 17.02.2000 Veronique Lefebure 39 // - correct bug : the gamma energy was not d 39 // - correct bug : the gamma energy was not deposited when the gamma was 40 // not produced when its energy was < cutFo 40 // not produced when its energy was < cutForLowEnergySecondaryPhotons 41 // 41 // 42 // Added Livermore data table construction met 42 // Added Livermore data table construction methods A. Forti 43 // Modified BuildMeanFreePath to read new data 43 // Modified BuildMeanFreePath to read new data tables A. Forti 44 // Modified PostStepDoIt to insert sampling wi 44 // Modified PostStepDoIt to insert sampling with with EEDL data A. Forti 45 // Added SelectRandomAtom A. Forti 45 // Added SelectRandomAtom A. Forti 46 // Added map of the elements A. Forti 46 // Added map of the elements A. Forti 47 // 20.09.00 update printout V.Ivanchenko 47 // 20.09.00 update printout V.Ivanchenko 48 // 24.04.01 V.Ivanchenko remove RogueWave 48 // 24.04.01 V.Ivanchenko remove RogueWave 49 // 29.09.2001 V.Ivanchenko: major revision bas 49 // 29.09.2001 V.Ivanchenko: major revision based on design iteration 50 // 10.10.2001 MGP Revision to improve code qua 50 // 10.10.2001 MGP Revision to improve code quality and consistency with design 51 // 18.10.2001 MGP Revision to improve code qua 51 // 18.10.2001 MGP Revision to improve code quality 52 // 28.10.2001 VI Update printout 52 // 28.10.2001 VI Update printout 53 // 29.11.2001 VI New parametrisation 53 // 29.11.2001 VI New parametrisation 54 // 30.07.2002 VI Fix in restricted energy los 54 // 30.07.2002 VI Fix in restricted energy loss 55 // 21.01.2003 VI Cut per region 55 // 21.01.2003 VI Cut per region 56 // 21.02.2003 V.Ivanchenko Energy bins for 56 // 21.02.2003 V.Ivanchenko Energy bins for spectrum are defined here 57 // 28.02.03 V.Ivanchenko Filename is define 57 // 28.02.03 V.Ivanchenko Filename is defined in the constructor 58 // 24.03.2003 P.Rodrigues Changes to accommoda 58 // 24.03.2003 P.Rodrigues Changes to accommodate new angular generators 59 // 20.05.2003 MGP Removed memory leak related 59 // 20.05.2003 MGP Removed memory leak related to angularDistribution 60 // 06.11.2003 MGP Improved user interface to 60 // 06.11.2003 MGP Improved user interface to select angular distribution model 61 // 61 // 62 // ------------------------------------------- 62 // -------------------------------------------------------------- 63 63 64 #include "G4LowEnergyBremsstrahlung.hh" 64 #include "G4LowEnergyBremsstrahlung.hh" 65 #include "G4RDeBremsstrahlungSpectrum.hh" 65 #include "G4RDeBremsstrahlungSpectrum.hh" 66 #include "G4RDBremsstrahlungCrossSectionHandle 66 #include "G4RDBremsstrahlungCrossSectionHandler.hh" 67 #include "G4RDVBremAngularDistribution.hh" 67 #include "G4RDVBremAngularDistribution.hh" 68 #include "G4RDModifiedTsai.hh" 68 #include "G4RDModifiedTsai.hh" 69 #include "G4RDGenerator2BS.hh" 69 #include "G4RDGenerator2BS.hh" 70 #include "G4RDGenerator2BN.hh" 70 #include "G4RDGenerator2BN.hh" 71 #include "G4RDVDataSetAlgorithm.hh" 71 #include "G4RDVDataSetAlgorithm.hh" 72 #include "G4RDLogLogInterpolation.hh" 72 #include "G4RDLogLogInterpolation.hh" 73 #include "G4RDVEMDataSet.hh" 73 #include "G4RDVEMDataSet.hh" 74 #include "G4EnergyLossTables.hh" 74 #include "G4EnergyLossTables.hh" 75 #include "G4PhysicalConstants.hh" 75 #include "G4PhysicalConstants.hh" 76 #include "G4SystemOfUnits.hh" 76 #include "G4SystemOfUnits.hh" 77 #include "G4UnitsTable.hh" 77 #include "G4UnitsTable.hh" 78 #include "G4Electron.hh" 78 #include "G4Electron.hh" 79 #include "G4Gamma.hh" 79 #include "G4Gamma.hh" 80 #include "G4ProductionCutsTable.hh" 80 #include "G4ProductionCutsTable.hh" 81 81 82 82 83 G4LowEnergyBremsstrahlung::G4LowEnergyBremsstr 83 G4LowEnergyBremsstrahlung::G4LowEnergyBremsstrahlung(const G4String& nam) 84 : G4eLowEnergyLoss(nam), 84 : G4eLowEnergyLoss(nam), 85 crossSectionHandler(0), 85 crossSectionHandler(0), 86 theMeanFreePath(0), 86 theMeanFreePath(0), 87 energySpectrum(0) 87 energySpectrum(0) 88 { 88 { 89 cutForPhotons = 0.; 89 cutForPhotons = 0.; 90 verboseLevel = 0; 90 verboseLevel = 0; 91 generatorName = "tsai"; 91 generatorName = "tsai"; 92 angularDistribution = new G4RDModifiedTsai(" 92 angularDistribution = new G4RDModifiedTsai("TsaiGenerator"); // default generator 93 // angularDistribution->PrintGeneratorInforma 93 // angularDistribution->PrintGeneratorInformation(); 94 TsaiAngularDistribution = new G4RDModifiedTs 94 TsaiAngularDistribution = new G4RDModifiedTsai("TsaiGenerator"); 95 } 95 } 96 96 97 /* 97 /* 98 G4LowEnergyBremsstrahlung::G4LowEnergyBremsstr 98 G4LowEnergyBremsstrahlung::G4LowEnergyBremsstrahlung(const G4String& nam, G4RDVBremAngularDistribution* distribution) 99 : G4eLowEnergyLoss(nam), 99 : G4eLowEnergyLoss(nam), 100 crossSectionHandler(0), 100 crossSectionHandler(0), 101 theMeanFreePath(0), 101 theMeanFreePath(0), 102 energySpectrum(0), 102 energySpectrum(0), 103 angularDistribution(distribution) 103 angularDistribution(distribution) 104 { 104 { 105 cutForPhotons = 0.; 105 cutForPhotons = 0.; 106 verboseLevel = 0; 106 verboseLevel = 0; 107 107 108 angularDistribution->PrintGeneratorInformati 108 angularDistribution->PrintGeneratorInformation(); 109 109 110 TsaiAngularDistribution = new G4RDModifiedTs 110 TsaiAngularDistribution = new G4RDModifiedTsai("TsaiGenerator"); 111 } 111 } 112 */ 112 */ 113 113 114 G4LowEnergyBremsstrahlung::~G4LowEnergyBremsst 114 G4LowEnergyBremsstrahlung::~G4LowEnergyBremsstrahlung() 115 { 115 { 116 if(crossSectionHandler) delete crossSectionH 116 if(crossSectionHandler) delete crossSectionHandler; 117 if(energySpectrum) delete energySpectrum; 117 if(energySpectrum) delete energySpectrum; 118 if(theMeanFreePath) delete theMeanFreePath; 118 if(theMeanFreePath) delete theMeanFreePath; 119 delete angularDistribution; 119 delete angularDistribution; 120 delete TsaiAngularDistribution; 120 delete TsaiAngularDistribution; 121 energyBins.clear(); 121 energyBins.clear(); 122 } 122 } 123 123 124 124 125 void G4LowEnergyBremsstrahlung::BuildPhysicsTa 125 void G4LowEnergyBremsstrahlung::BuildPhysicsTable(const G4ParticleDefinition& aParticleType) 126 { 126 { 127 if(verboseLevel > 0) { 127 if(verboseLevel > 0) { 128 G4cout << "G4LowEnergyBremsstrahlung::Buil 128 G4cout << "G4LowEnergyBremsstrahlung::BuildPhysicsTable start" 129 << G4endl; 129 << G4endl; 130 } 130 } 131 131 132 cutForSecondaryPhotons.clear(); 132 cutForSecondaryPhotons.clear(); 133 133 134 // Create and fill BremsstrahlungParameters 134 // Create and fill BremsstrahlungParameters once 135 if( energySpectrum != 0 ) delete energySpect 135 if( energySpectrum != 0 ) delete energySpectrum; 136 energyBins.clear(); 136 energyBins.clear(); 137 for(size_t i=0; i<15; i++) { 137 for(size_t i=0; i<15; i++) { 138 G4double x = 0.1*((G4double)i); 138 G4double x = 0.1*((G4double)i); 139 if(i == 0) x = 0.01; 139 if(i == 0) x = 0.01; 140 if(i == 10) x = 0.95; 140 if(i == 10) x = 0.95; 141 if(i == 11) x = 0.97; 141 if(i == 11) x = 0.97; 142 if(i == 12) x = 0.99; 142 if(i == 12) x = 0.99; 143 if(i == 13) x = 0.995; 143 if(i == 13) x = 0.995; 144 if(i == 14) x = 1.0; 144 if(i == 14) x = 1.0; 145 energyBins.push_back(x); 145 energyBins.push_back(x); 146 } 146 } 147 const G4String dataName("/brem/br-sp.dat"); 147 const G4String dataName("/brem/br-sp.dat"); 148 energySpectrum = new G4RDeBremsstrahlungSpec 148 energySpectrum = new G4RDeBremsstrahlungSpectrum(energyBins,dataName); 149 149 150 if(verboseLevel > 0) { 150 if(verboseLevel > 0) { 151 G4cout << "G4LowEnergyBremsstrahlungSpectr 151 G4cout << "G4LowEnergyBremsstrahlungSpectrum is initialized" 152 << G4endl; 152 << G4endl; 153 } 153 } 154 154 155 // Create and fill G4RDCrossSectionHandler o 155 // Create and fill G4RDCrossSectionHandler once 156 156 157 if( crossSectionHandler != 0 ) delete crossS 157 if( crossSectionHandler != 0 ) delete crossSectionHandler; 158 G4RDVDataSetAlgorithm* interpolation = new G 158 G4RDVDataSetAlgorithm* interpolation = new G4RDLogLogInterpolation(); 159 G4double lowKineticEnergy = GetLowerBoundEl 159 G4double lowKineticEnergy = GetLowerBoundEloss(); 160 G4double highKineticEnergy = GetUpperBoundEl 160 G4double highKineticEnergy = GetUpperBoundEloss(); 161 G4int totBin = GetNbinEloss(); 161 G4int totBin = GetNbinEloss(); 162 crossSectionHandler = new G4RDBremsstrahlung 162 crossSectionHandler = new G4RDBremsstrahlungCrossSectionHandler(energySpectrum, interpolation); 163 crossSectionHandler->Initialise(0,lowKinetic 163 crossSectionHandler->Initialise(0,lowKineticEnergy, highKineticEnergy, totBin); 164 crossSectionHandler->LoadShellData("brem/br- 164 crossSectionHandler->LoadShellData("brem/br-cs-"); 165 165 166 if (verboseLevel > 0) { 166 if (verboseLevel > 0) { 167 G4cout << GetProcessName() 167 G4cout << GetProcessName() 168 << " is created; Cross section data 168 << " is created; Cross section data: " 169 << G4endl; 169 << G4endl; 170 crossSectionHandler->PrintData(); 170 crossSectionHandler->PrintData(); 171 G4cout << "Parameters: " 171 G4cout << "Parameters: " 172 << G4endl; 172 << G4endl; 173 energySpectrum->PrintData(); 173 energySpectrum->PrintData(); 174 } 174 } 175 175 176 // Build loss table for Bremsstrahlung 176 // Build loss table for Bremsstrahlung 177 177 178 BuildLossTable(aParticleType); 178 BuildLossTable(aParticleType); 179 179 180 if(verboseLevel > 0) { 180 if(verboseLevel > 0) { 181 G4cout << "The loss table is built" 181 G4cout << "The loss table is built" 182 << G4endl; 182 << G4endl; 183 } 183 } 184 184 185 if (&aParticleType==G4Electron::Electron()) 185 if (&aParticleType==G4Electron::Electron()) { 186 186 187 RecorderOfElectronProcess[CounterOfElectro 187 RecorderOfElectronProcess[CounterOfElectronProcess] = (*this).theLossTable; 188 CounterOfElectronProcess++; 188 CounterOfElectronProcess++; 189 PrintInfoDefinition(); 189 PrintInfoDefinition(); 190 190 191 } else { 191 } else { 192 192 193 RecorderOfPositronProcess[CounterOfPositro 193 RecorderOfPositronProcess[CounterOfPositronProcess] = (*this).theLossTable; 194 CounterOfPositronProcess++; 194 CounterOfPositronProcess++; 195 } 195 } 196 196 197 // Build mean free path data using cut value 197 // Build mean free path data using cut values 198 198 199 if( theMeanFreePath != 0 ) delete theMeanFre 199 if( theMeanFreePath != 0 ) delete theMeanFreePath; 200 theMeanFreePath = crossSectionHandler-> 200 theMeanFreePath = crossSectionHandler-> 201 BuildMeanFreePathForMateri 201 BuildMeanFreePathForMaterials(&cutForSecondaryPhotons); 202 202 203 if(verboseLevel > 0) { 203 if(verboseLevel > 0) { 204 G4cout << "The MeanFreePath table is built 204 G4cout << "The MeanFreePath table is built" 205 << G4endl; 205 << G4endl; 206 } 206 } 207 207 208 // Build common DEDX table for all ionisatio 208 // Build common DEDX table for all ionisation processes 209 209 210 BuildDEDXTable(aParticleType); 210 BuildDEDXTable(aParticleType); 211 211 212 if(verboseLevel > 0) { 212 if(verboseLevel > 0) { 213 G4cout << "G4LowEnergyBremsstrahlung::Buil 213 G4cout << "G4LowEnergyBremsstrahlung::BuildPhysicsTable end" 214 << G4endl; 214 << G4endl; 215 } 215 } 216 216 217 } 217 } 218 218 219 219 220 void G4LowEnergyBremsstrahlung::BuildLossTable 220 void G4LowEnergyBremsstrahlung::BuildLossTable(const G4ParticleDefinition& ) 221 { 221 { 222 // Build table for energy loss due to soft b 222 // Build table for energy loss due to soft brems 223 // the tables are built for *MATERIALS* binn 223 // the tables are built for *MATERIALS* binning is taken from LowEnergyLoss 224 224 225 G4double lowKineticEnergy = GetLowerBoundEl 225 G4double lowKineticEnergy = GetLowerBoundEloss(); 226 G4double highKineticEnergy = GetUpperBoundEl 226 G4double highKineticEnergy = GetUpperBoundEloss(); 227 size_t totBin = GetNbinEloss(); 227 size_t totBin = GetNbinEloss(); 228 228 229 // create table 229 // create table 230 230 231 if (theLossTable) { 231 if (theLossTable) { 232 theLossTable->clearAndDestroy(); 232 theLossTable->clearAndDestroy(); 233 delete theLossTable; 233 delete theLossTable; 234 } 234 } 235 const G4ProductionCutsTable* theCoupleTable= 235 const G4ProductionCutsTable* theCoupleTable= 236 G4ProductionCutsTable::GetProductionCu 236 G4ProductionCutsTable::GetProductionCutsTable(); 237 size_t numOfCouples = theCoupleTable->GetTab 237 size_t numOfCouples = theCoupleTable->GetTableSize(); 238 theLossTable = new G4PhysicsTable(numOfCoupl 238 theLossTable = new G4PhysicsTable(numOfCouples); 239 239 240 // Clean up the vector of cuts 240 // Clean up the vector of cuts 241 cutForSecondaryPhotons.clear(); 241 cutForSecondaryPhotons.clear(); 242 242 243 // Loop for materials 243 // Loop for materials 244 244 245 for (size_t j=0; j<numOfCouples; j++) { 245 for (size_t j=0; j<numOfCouples; j++) { 246 246 247 // create physics vector and fill it 247 // create physics vector and fill it 248 G4PhysicsLogVector* aVector = new G4Physic 248 G4PhysicsLogVector* aVector = new G4PhysicsLogVector(lowKineticEnergy, 249 highKineticEnergy, 249 highKineticEnergy, 250 totBin); 250 totBin); 251 251 252 // get material parameters needed for the 252 // get material parameters needed for the energy loss calculation 253 const G4MaterialCutsCouple* couple = theCo 253 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j); 254 const G4Material* material= couple->GetMat 254 const G4Material* material= couple->GetMaterial(); 255 255 256 // the cut cannot be below lowest limit 256 // the cut cannot be below lowest limit 257 G4double tCut = (*(theCoupleTable->GetEner 257 G4double tCut = (*(theCoupleTable->GetEnergyCutsVector(0)))[j]; 258 tCut = std::min(highKineticEnergy, tCut); 258 tCut = std::min(highKineticEnergy, tCut); 259 cutForSecondaryPhotons.push_back(tCut); 259 cutForSecondaryPhotons.push_back(tCut); 260 260 261 const G4ElementVector* theElementVector = 261 const G4ElementVector* theElementVector = material->GetElementVector(); 262 size_t NumberOfElements = material->GetNum 262 size_t NumberOfElements = material->GetNumberOfElements() ; 263 const G4double* theAtomicNumDensityVector 263 const G4double* theAtomicNumDensityVector = 264 material->GetAtomicNumDensityVector(); 264 material->GetAtomicNumDensityVector(); 265 if(verboseLevel > 1) { 265 if(verboseLevel > 1) { 266 G4cout << "Energy loss for material # " 266 G4cout << "Energy loss for material # " << j 267 << " tCut(keV)= " << tCut/keV 267 << " tCut(keV)= " << tCut/keV 268 << G4endl; 268 << G4endl; 269 } 269 } 270 270 271 // now comes the loop for the kinetic ener 271 // now comes the loop for the kinetic energy values 272 for (size_t i = 0; i<totBin; i++) { 272 for (size_t i = 0; i<totBin; i++) { 273 273 274 G4double lowEdgeEnergy = aVector->GetLow 274 G4double lowEdgeEnergy = aVector->GetLowEdgeEnergy(i); 275 G4double ionloss = 0.; 275 G4double ionloss = 0.; 276 276 277 // loop for elements in the material 277 // loop for elements in the material 278 for (size_t iel=0; iel<NumberOfElements; 278 for (size_t iel=0; iel<NumberOfElements; iel++ ) { 279 G4int Z = (G4int)((*theElementVector)[ 279 G4int Z = (G4int)((*theElementVector)[iel]->GetZ()); 280 G4double e = energySpectrum->AverageEn 280 G4double e = energySpectrum->AverageEnergy(Z, 0.0, tCut, lowEdgeEnergy); 281 G4double cs= crossSectionHandler->Find 281 G4double cs= crossSectionHandler->FindValue(Z, lowEdgeEnergy); 282 ionloss += e * cs * theAtomicNumDen 282 ionloss += e * cs * theAtomicNumDensityVector[iel]; 283 if(verboseLevel > 1) { 283 if(verboseLevel > 1) { 284 G4cout << "Z= " << Z 284 G4cout << "Z= " << Z 285 << "; tCut(keV)= " << tCut/ke 285 << "; tCut(keV)= " << tCut/keV 286 << "; E(keV)= " << lowEdgeEne 286 << "; E(keV)= " << lowEdgeEnergy/keV 287 << "; Eav(keV)= " << e/keV 287 << "; Eav(keV)= " << e/keV 288 << "; cs= " << cs 288 << "; cs= " << cs 289 << "; loss= " << ionloss 289 << "; loss= " << ionloss 290 << G4endl; 290 << G4endl; 291 } 291 } 292 } 292 } 293 aVector->PutValue(i,ionloss); 293 aVector->PutValue(i,ionloss); 294 } 294 } 295 theLossTable->insert(aVector); 295 theLossTable->insert(aVector); 296 } 296 } 297 } 297 } 298 298 299 299 300 G4VParticleChange* G4LowEnergyBremsstrahlung:: 300 G4VParticleChange* G4LowEnergyBremsstrahlung::PostStepDoIt(const G4Track& track, 301 const G4Step& step) 301 const G4Step& step) 302 { 302 { 303 aParticleChange.Initialize(track); 303 aParticleChange.Initialize(track); 304 304 305 const G4MaterialCutsCouple* couple = track.G 305 const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple(); 306 G4double kineticEnergy = track.GetKineticEne 306 G4double kineticEnergy = track.GetKineticEnergy(); 307 G4int index = couple->GetIndex(); 307 G4int index = couple->GetIndex(); 308 G4double tCut = cutForSecondaryPhotons[index 308 G4double tCut = cutForSecondaryPhotons[index]; 309 309 310 // Control limits 310 // Control limits 311 if(tCut >= kineticEnergy) 311 if(tCut >= kineticEnergy) 312 return G4VContinuousDiscreteProcess::Post 312 return G4VContinuousDiscreteProcess::PostStepDoIt(track, step); 313 313 314 G4int Z = crossSectionHandler->SelectRandomA 314 G4int Z = crossSectionHandler->SelectRandomAtom(couple, kineticEnergy); 315 315 316 G4double tGamma = energySpectrum->SampleEner 316 G4double tGamma = energySpectrum->SampleEnergy(Z, tCut, kineticEnergy, kineticEnergy); 317 G4double totalEnergy = kineticEnergy + elect 317 G4double totalEnergy = kineticEnergy + electron_mass_c2; 318 G4double finalEnergy = kineticEnergy - tGamm 318 G4double finalEnergy = kineticEnergy - tGamma; // electron/positron final energy 319 G4double theta = 0; 319 G4double theta = 0; 320 320 321 if((kineticEnergy < 1*MeV && kineticEnergy > 321 if((kineticEnergy < 1*MeV && kineticEnergy > 1*keV && generatorName == "2bn")){ 322 theta = angularDistribution->PolarAngle( 322 theta = angularDistribution->PolarAngle(kineticEnergy,finalEnergy,Z); 323 }else{ 323 }else{ 324 theta = TsaiAngularDistribution->PolarAn 324 theta = TsaiAngularDistribution->PolarAngle(kineticEnergy,finalEnergy,Z); 325 } 325 } 326 326 327 G4double phi = twopi * G4UniformRand(); 327 G4double phi = twopi * G4UniformRand(); 328 G4double dirZ = std::cos(theta); 328 G4double dirZ = std::cos(theta); 329 G4double sinTheta = std::sqrt(1. - dirZ*dir 329 G4double sinTheta = std::sqrt(1. - dirZ*dirZ); 330 G4double dirX = sinTheta*std::cos(phi); 330 G4double dirX = sinTheta*std::cos(phi); 331 G4double dirY = sinTheta*std::sin(phi); 331 G4double dirY = sinTheta*std::sin(phi); 332 332 333 G4ThreeVector gammaDirection (dirX, dirY, di 333 G4ThreeVector gammaDirection (dirX, dirY, dirZ); 334 G4ThreeVector electronDirection = track.GetM 334 G4ThreeVector electronDirection = track.GetMomentumDirection(); 335 335 336 // 336 // 337 // Update the incident particle 337 // Update the incident particle 338 // 338 // 339 gammaDirection.rotateUz(electronDirection); 339 gammaDirection.rotateUz(electronDirection); 340 340 341 // Kinematic problem 341 // Kinematic problem 342 if (finalEnergy < 0.) { 342 if (finalEnergy < 0.) { 343 tGamma += finalEnergy; 343 tGamma += finalEnergy; 344 finalEnergy = 0.0; 344 finalEnergy = 0.0; 345 } 345 } 346 346 347 G4double momentum = std::sqrt((totalEnergy + 347 G4double momentum = std::sqrt((totalEnergy + electron_mass_c2)*kineticEnergy); 348 348 349 G4double finalX = momentum*electronDirection 349 G4double finalX = momentum*electronDirection.x() - tGamma*gammaDirection.x(); 350 G4double finalY = momentum*electronDirection 350 G4double finalY = momentum*electronDirection.y() - tGamma*gammaDirection.y(); 351 G4double finalZ = momentum*electronDirection 351 G4double finalZ = momentum*electronDirection.z() - tGamma*gammaDirection.z(); 352 352 353 aParticleChange.SetNumberOfSecondaries(1); 353 aParticleChange.SetNumberOfSecondaries(1); 354 G4double norm = 1./std::sqrt(finalX*finalX + 354 G4double norm = 1./std::sqrt(finalX*finalX + finalY*finalY + finalZ*finalZ); 355 aParticleChange.ProposeMomentumDirection(fin 355 aParticleChange.ProposeMomentumDirection(finalX*norm, finalY*norm, finalZ*norm); 356 aParticleChange.ProposeEnergy( finalEnergy ) 356 aParticleChange.ProposeEnergy( finalEnergy ); 357 357 358 // create G4DynamicParticle object for the g 358 // create G4DynamicParticle object for the gamma 359 G4DynamicParticle* aGamma= new G4DynamicPart 359 G4DynamicParticle* aGamma= new G4DynamicParticle (G4Gamma::Gamma(), 360 gammaDirection, tGamma); 360 gammaDirection, tGamma); 361 aParticleChange.AddSecondary(aGamma); 361 aParticleChange.AddSecondary(aGamma); 362 362 363 return G4VContinuousDiscreteProcess::PostSte 363 return G4VContinuousDiscreteProcess::PostStepDoIt(track, step); 364 } 364 } 365 365 366 366 367 void G4LowEnergyBremsstrahlung::PrintInfoDefin 367 void G4LowEnergyBremsstrahlung::PrintInfoDefinition() 368 { 368 { 369 G4String comments = "Total cross sections fr 369 G4String comments = "Total cross sections from EEDL database."; 370 comments += "\n Gamma energy sampled fr 370 comments += "\n Gamma energy sampled from a parameterised formula."; 371 comments += "\n Implementation of the c 371 comments += "\n Implementation of the continuous dE/dx part."; 372 comments += "\n At present it can be us 372 comments += "\n At present it can be used for electrons "; 373 comments += "in the energy range [250eV,100G 373 comments += "in the energy range [250eV,100GeV]."; 374 comments += "\n The process must work w 374 comments += "\n The process must work with G4LowEnergyIonisation."; 375 375 376 G4cout << G4endl << GetProcessName() << ": 376 G4cout << G4endl << GetProcessName() << ": " << comments << G4endl; 377 } 377 } 378 378 379 G4bool G4LowEnergyBremsstrahlung::IsApplicable 379 G4bool G4LowEnergyBremsstrahlung::IsApplicable(const G4ParticleDefinition& particle) 380 { 380 { 381 return ( (&particle == G4Electron::Electron 381 return ( (&particle == G4Electron::Electron()) ); 382 } 382 } 383 383 384 384 385 G4double G4LowEnergyBremsstrahlung::GetMeanFre 385 G4double G4LowEnergyBremsstrahlung::GetMeanFreePath(const G4Track& track, 386 G4double, 386 G4double, 387 G4ForceCondition* cond) 387 G4ForceCondition* cond) 388 { 388 { 389 *cond = NotForced; 389 *cond = NotForced; 390 G4int index = (track.GetMaterialCutsCouple() 390 G4int index = (track.GetMaterialCutsCouple())->GetIndex(); 391 const G4RDVEMDataSet* data = theMeanFreePath 391 const G4RDVEMDataSet* data = theMeanFreePath->GetComponent(index); 392 G4double meanFreePath = data->FindValue(trac 392 G4double meanFreePath = data->FindValue(track.GetKineticEnergy()); 393 return meanFreePath; 393 return meanFreePath; 394 } 394 } 395 395 396 void G4LowEnergyBremsstrahlung::SetCutForLowEn 396 void G4LowEnergyBremsstrahlung::SetCutForLowEnSecPhotons(G4double cut) 397 { 397 { 398 cutForPhotons = cut; 398 cutForPhotons = cut; 399 } 399 } 400 400 401 void G4LowEnergyBremsstrahlung::SetAngularGene 401 void G4LowEnergyBremsstrahlung::SetAngularGenerator(G4RDVBremAngularDistribution* distribution) 402 { 402 { 403 angularDistribution = distribution; 403 angularDistribution = distribution; 404 angularDistribution->PrintGeneratorInformati 404 angularDistribution->PrintGeneratorInformation(); 405 } 405 } 406 406 407 void G4LowEnergyBremsstrahlung::SetAngularGene 407 void G4LowEnergyBremsstrahlung::SetAngularGenerator(const G4String& name) 408 { 408 { 409 if (name == "tsai") 409 if (name == "tsai") 410 { 410 { 411 delete angularDistribution; 411 delete angularDistribution; 412 angularDistribution = new G4RDModifiedTs 412 angularDistribution = new G4RDModifiedTsai("TsaiGenerator"); 413 generatorName = name; 413 generatorName = name; 414 } 414 } 415 else if (name == "2bn") 415 else if (name == "2bn") 416 { 416 { 417 delete angularDistribution; 417 delete angularDistribution; 418 angularDistribution = new G4RDGenerator2 418 angularDistribution = new G4RDGenerator2BN("2BNGenerator"); 419 generatorName = name; 419 generatorName = name; 420 } 420 } 421 else if (name == "2bs") 421 else if (name == "2bs") 422 { 422 { 423 delete angularDistribution; 423 delete angularDistribution; 424 angularDistribution = new G4RDGenerator 424 angularDistribution = new G4RDGenerator2BS("2BSGenerator"); 425 generatorName = name; 425 generatorName = name; 426 } 426 } 427 else 427 else 428 { 428 { 429 G4Exception("G4LowEnergyBremsstrahlung:: 429 G4Exception("G4LowEnergyBremsstrahlung::SetAngularGenerator()", 430 "InvalidSetup", FatalExcepti 430 "InvalidSetup", FatalException, "Generator does not exist!"); 431 } 431 } 432 432 433 angularDistribution->PrintGeneratorInformati 433 angularDistribution->PrintGeneratorInformation(); 434 } 434 } 435 435 436 436