Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 27 // ------------------------------------------- 28 // 29 // File name: G4LowEnergyBremsstrahlung 30 // 31 // Author: Alessandra Forti, Vladimir I 32 // 33 // Creation date: March 1999 34 // 35 // Modifications: 36 // 18.04.2000 V.L. 37 // - First implementation of continuous energ 38 // 17.02.2000 Veronique Lefebure 39 // - correct bug : the gamma energy was not d 40 // not produced when its energy was < cutFo 41 // 42 // Added Livermore data table construction met 43 // Modified BuildMeanFreePath to read new data 44 // Modified PostStepDoIt to insert sampling wi 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 bas 50 // 10.10.2001 MGP Revision to improve code qua 51 // 18.10.2001 MGP Revision to improve code qua 52 // 28.10.2001 VI Update printout 53 // 29.11.2001 VI New parametrisation 54 // 30.07.2002 VI Fix in restricted energy los 55 // 21.01.2003 VI Cut per region 56 // 21.02.2003 V.Ivanchenko Energy bins for 57 // 28.02.03 V.Ivanchenko Filename is define 58 // 24.03.2003 P.Rodrigues Changes to accommoda 59 // 20.05.2003 MGP Removed memory leak related 60 // 06.11.2003 MGP Improved user interface to 61 // 62 // ------------------------------------------- 63 64 #include "G4LowEnergyBremsstrahlung.hh" 65 #include "G4RDeBremsstrahlungSpectrum.hh" 66 #include "G4RDBremsstrahlungCrossSectionHandle 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::G4LowEnergyBremsstr 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(" 93 // angularDistribution->PrintGeneratorInforma 94 TsaiAngularDistribution = new G4RDModifiedTs 95 } 96 97 /* 98 G4LowEnergyBremsstrahlung::G4LowEnergyBremsstr 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->PrintGeneratorInformati 109 110 TsaiAngularDistribution = new G4RDModifiedTs 111 } 112 */ 113 114 G4LowEnergyBremsstrahlung::~G4LowEnergyBremsst 115 { 116 if(crossSectionHandler) delete crossSectionH 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::BuildPhysicsTa 126 { 127 if(verboseLevel > 0) { 128 G4cout << "G4LowEnergyBremsstrahlung::Buil 129 << G4endl; 130 } 131 132 cutForSecondaryPhotons.clear(); 133 134 // Create and fill BremsstrahlungParameters 135 if( energySpectrum != 0 ) delete energySpect 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 G4RDeBremsstrahlungSpec 149 150 if(verboseLevel > 0) { 151 G4cout << "G4LowEnergyBremsstrahlungSpectr 152 << G4endl; 153 } 154 155 // Create and fill G4RDCrossSectionHandler o 156 157 if( crossSectionHandler != 0 ) delete crossS 158 G4RDVDataSetAlgorithm* interpolation = new G 159 G4double lowKineticEnergy = GetLowerBoundEl 160 G4double highKineticEnergy = GetUpperBoundEl 161 G4int totBin = GetNbinEloss(); 162 crossSectionHandler = new G4RDBremsstrahlung 163 crossSectionHandler->Initialise(0,lowKinetic 164 crossSectionHandler->LoadShellData("brem/br- 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[CounterOfElectro 188 CounterOfElectronProcess++; 189 PrintInfoDefinition(); 190 191 } else { 192 193 RecorderOfPositronProcess[CounterOfPositro 194 CounterOfPositronProcess++; 195 } 196 197 // Build mean free path data using cut value 198 199 if( theMeanFreePath != 0 ) delete theMeanFre 200 theMeanFreePath = crossSectionHandler-> 201 BuildMeanFreePathForMateri 202 203 if(verboseLevel > 0) { 204 G4cout << "The MeanFreePath table is built 205 << G4endl; 206 } 207 208 // Build common DEDX table for all ionisatio 209 210 BuildDEDXTable(aParticleType); 211 212 if(verboseLevel > 0) { 213 G4cout << "G4LowEnergyBremsstrahlung::Buil 214 << G4endl; 215 } 216 217 } 218 219 220 void G4LowEnergyBremsstrahlung::BuildLossTable 221 { 222 // Build table for energy loss due to soft b 223 // the tables are built for *MATERIALS* binn 224 225 G4double lowKineticEnergy = GetLowerBoundEl 226 G4double highKineticEnergy = GetUpperBoundEl 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::GetProductionCu 237 size_t numOfCouples = theCoupleTable->GetTab 238 theLossTable = new G4PhysicsTable(numOfCoupl 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 G4Physic 249 highKineticEnergy, 250 totBin); 251 252 // get material parameters needed for the 253 const G4MaterialCutsCouple* couple = theCo 254 const G4Material* material= couple->GetMat 255 256 // the cut cannot be below lowest limit 257 G4double tCut = (*(theCoupleTable->GetEner 258 tCut = std::min(highKineticEnergy, tCut); 259 cutForSecondaryPhotons.push_back(tCut); 260 261 const G4ElementVector* theElementVector = 262 size_t NumberOfElements = material->GetNum 263 const G4double* theAtomicNumDensityVector 264 material->GetAtomicNumDensityVector(); 265 if(verboseLevel > 1) { 266 G4cout << "Energy loss for material # " 267 << " tCut(keV)= " << tCut/keV 268 << G4endl; 269 } 270 271 // now comes the loop for the kinetic ener 272 for (size_t i = 0; i<totBin; i++) { 273 274 G4double lowEdgeEnergy = aVector->GetLow 275 G4double ionloss = 0.; 276 277 // loop for elements in the material 278 for (size_t iel=0; iel<NumberOfElements; 279 G4int Z = (G4int)((*theElementVector)[ 280 G4double e = energySpectrum->AverageEn 281 G4double cs= crossSectionHandler->Find 282 ionloss += e * cs * theAtomicNumDen 283 if(verboseLevel > 1) { 284 G4cout << "Z= " << Z 285 << "; tCut(keV)= " << tCut/ke 286 << "; E(keV)= " << lowEdgeEne 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:: 301 const G4Step& step) 302 { 303 aParticleChange.Initialize(track); 304 305 const G4MaterialCutsCouple* couple = track.G 306 G4double kineticEnergy = track.GetKineticEne 307 G4int index = couple->GetIndex(); 308 G4double tCut = cutForSecondaryPhotons[index 309 310 // Control limits 311 if(tCut >= kineticEnergy) 312 return G4VContinuousDiscreteProcess::Post 313 314 G4int Z = crossSectionHandler->SelectRandomA 315 316 G4double tGamma = energySpectrum->SampleEner 317 G4double totalEnergy = kineticEnergy + elect 318 G4double finalEnergy = kineticEnergy - tGamm 319 G4double theta = 0; 320 321 if((kineticEnergy < 1*MeV && kineticEnergy > 322 theta = angularDistribution->PolarAngle( 323 }else{ 324 theta = TsaiAngularDistribution->PolarAn 325 } 326 327 G4double phi = twopi * G4UniformRand(); 328 G4double dirZ = std::cos(theta); 329 G4double sinTheta = std::sqrt(1. - dirZ*dir 330 G4double dirX = sinTheta*std::cos(phi); 331 G4double dirY = sinTheta*std::sin(phi); 332 333 G4ThreeVector gammaDirection (dirX, dirY, di 334 G4ThreeVector electronDirection = track.GetM 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 + 348 349 G4double finalX = momentum*electronDirection 350 G4double finalY = momentum*electronDirection 351 G4double finalZ = momentum*electronDirection 352 353 aParticleChange.SetNumberOfSecondaries(1); 354 G4double norm = 1./std::sqrt(finalX*finalX + 355 aParticleChange.ProposeMomentumDirection(fin 356 aParticleChange.ProposeEnergy( finalEnergy ) 357 358 // create G4DynamicParticle object for the g 359 G4DynamicParticle* aGamma= new G4DynamicPart 360 gammaDirection, tGamma); 361 aParticleChange.AddSecondary(aGamma); 362 363 return G4VContinuousDiscreteProcess::PostSte 364 } 365 366 367 void G4LowEnergyBremsstrahlung::PrintInfoDefin 368 { 369 G4String comments = "Total cross sections fr 370 comments += "\n Gamma energy sampled fr 371 comments += "\n Implementation of the c 372 comments += "\n At present it can be us 373 comments += "in the energy range [250eV,100G 374 comments += "\n The process must work w 375 376 G4cout << G4endl << GetProcessName() << ": 377 } 378 379 G4bool G4LowEnergyBremsstrahlung::IsApplicable 380 { 381 return ( (&particle == G4Electron::Electron 382 } 383 384 385 G4double G4LowEnergyBremsstrahlung::GetMeanFre 386 G4double, 387 G4ForceCondition* cond) 388 { 389 *cond = NotForced; 390 G4int index = (track.GetMaterialCutsCouple() 391 const G4RDVEMDataSet* data = theMeanFreePath 392 G4double meanFreePath = data->FindValue(trac 393 return meanFreePath; 394 } 395 396 void G4LowEnergyBremsstrahlung::SetCutForLowEn 397 { 398 cutForPhotons = cut; 399 } 400 401 void G4LowEnergyBremsstrahlung::SetAngularGene 402 { 403 angularDistribution = distribution; 404 angularDistribution->PrintGeneratorInformati 405 } 406 407 void G4LowEnergyBremsstrahlung::SetAngularGene 408 { 409 if (name == "tsai") 410 { 411 delete angularDistribution; 412 angularDistribution = new G4RDModifiedTs 413 generatorName = name; 414 } 415 else if (name == "2bn") 416 { 417 delete angularDistribution; 418 angularDistribution = new G4RDGenerator2 419 generatorName = name; 420 } 421 else if (name == "2bs") 422 { 423 delete angularDistribution; 424 angularDistribution = new G4RDGenerator 425 generatorName = name; 426 } 427 else 428 { 429 G4Exception("G4LowEnergyBremsstrahlung:: 430 "InvalidSetup", FatalExcepti 431 } 432 433 angularDistribution->PrintGeneratorInformati 434 } 435 436