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