physics/src/G4LowEnergyBremsstrahlung.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 // --------------------------------------------------------------
0028 //
0029 // File name:     G4LowEnergyBremsstrahlung
0030 //
0031 // Author:        Alessandra Forti, Vladimir Ivanchenko
0032 //
0033 // Creation date: March 1999
0034 //
0035 // Modifications:
0036 // 18.04.2000 V.L.
0037 //  - First implementation of continuous energy loss.
0038 // 17.02.2000 Veronique Lefebure
0039 //  - correct bug : the gamma energy was not deposited when the gamma was
0040 //    not produced when its energy was < cutForLowEnergySecondaryPhotons
0041 //
0042 // Added Livermore data table construction methods A. Forti
0043 // Modified BuildMeanFreePath to read new data tables A. Forti
0044 // Modified PostStepDoIt to insert sampling with with EEDL data A. Forti
0045 // Added SelectRandomAtom A. Forti
0046 // Added map of the elements A. Forti
0047 // 20.09.00 update printout V.Ivanchenko
0048 // 24.04.01 V.Ivanchenko remove RogueWave
0049 // 29.09.2001 V.Ivanchenko: major revision based on design iteration
0050 // 10.10.2001 MGP Revision to improve code quality and consistency with design
0051 // 18.10.2001 MGP Revision to improve code quality
0052 // 28.10.2001 VI  Update printout
0053 // 29.11.2001 VI  New parametrisation
0054 // 30.07.2002 VI  Fix in restricted energy loss
0055 // 21.01.2003 VI  Cut per region
0056 // 21.02.2003 V.Ivanchenko    Energy bins for spectrum are defined here
0057 // 28.02.03 V.Ivanchenko    Filename is defined in the constructor
0058 // 24.03.2003 P.Rodrigues Changes to accommodate new angular generators
0059 // 20.05.2003 MGP  Removed memory leak related to angularDistribution
0060 // 06.11.2003 MGP  Improved user interface to select angular distribution model
0061 //
0062 // --------------------------------------------------------------
0063
0064 #include "G4LowEnergyBremsstrahlung.hh"
0065 #include "G4RDeBremsstrahlungSpectrum.hh"
0066 #include "G4RDBremsstrahlungCrossSectionHandler.hh"
0067 #include "G4RDVBremAngularDistribution.hh"
0068 #include "G4RDModifiedTsai.hh"
0069 #include "G4RDGenerator2BS.hh"
0070 #include "G4RDGenerator2BN.hh"
0071 #include "G4RDVDataSetAlgorithm.hh"
0072 #include "G4RDLogLogInterpolation.hh"
0073 #include "G4RDVEMDataSet.hh"
0074 #include "G4EnergyLossTables.hh"
0075 #include "G4PhysicalConstants.hh"
0076 #include "G4SystemOfUnits.hh"
0077 #include "G4UnitsTable.hh"
0078 #include "G4Electron.hh"
0079 #include "G4Gamma.hh"
0080 #include "G4ProductionCutsTable.hh"
0081
0082
0083 G4LowEnergyBremsstrahlung::G4LowEnergyBremsstrahlung(const G4String& nam)
0084   : G4eLowEnergyLoss(nam),
0085   crossSectionHandler(0),
0086   theMeanFreePath(0),
0087   energySpectrum(0)
0088 {
0089   cutForPhotons = 0.;
0090   verboseLevel = 0;
0091   generatorName = "tsai";
0092   angularDistribution = new G4RDModifiedTsai("TsaiGenerator"); // default generator
0093 //  angularDistribution->PrintGeneratorInformation();
0094   TsaiAngularDistribution = new G4RDModifiedTsai("TsaiGenerator");
0095 }
0096
0097 /*
0098 G4LowEnergyBremsstrahlung::G4LowEnergyBremsstrahlung(const G4String& nam, G4RDVBremAngularDistribution* distribution)
0099   : G4eLowEnergyLoss(nam),
0100     crossSectionHandler(0),
0101     theMeanFreePath(0),
0102     energySpectrum(0),
0103     angularDistribution(distribution)
0104 {
0105   cutForPhotons = 0.;
0106   verboseLevel = 0;
0107
0108   angularDistribution->PrintGeneratorInformation();
0109
0110   TsaiAngularDistribution = new G4RDModifiedTsai("TsaiGenerator");
0111 }
0112 */
0113
0114 G4LowEnergyBremsstrahlung::~G4LowEnergyBremsstrahlung()
0115 {
0116   if(crossSectionHandler) delete crossSectionHandler;
0117   if(energySpectrum) delete energySpectrum;
0118   if(theMeanFreePath) delete theMeanFreePath;
0119   delete angularDistribution;
0120   delete TsaiAngularDistribution;
0121   energyBins.clear();
0122 }
0123
0124
0125 void G4LowEnergyBremsstrahlung::BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
0126 {
0127   if(verboseLevel > 0) {
0128     G4cout << "G4LowEnergyBremsstrahlung::BuildPhysicsTable start"
0129            << G4endl;
0130       }
0131
0132   cutForSecondaryPhotons.clear();
0133
0134   // Create and fill BremsstrahlungParameters once
0135   if( energySpectrum != 0 ) delete energySpectrum;
0136   energyBins.clear();
0137   for(size_t i=0; i<15; i++) {
0138     G4double x = 0.1*((G4double)i);
0139     if(i == 0)  x = 0.01;
0140     if(i == 10) x = 0.95;
0141     if(i == 11) x = 0.97;
0142     if(i == 12) x = 0.99;
0143     if(i == 13) x = 0.995;
0144     if(i == 14) x = 1.0;
0145     energyBins.push_back(x);
0146   }
0147   const G4String dataName("/brem/br-sp.dat");
0148   energySpectrum = new G4RDeBremsstrahlungSpectrum(energyBins,dataName);
0149
0150   if(verboseLevel > 0) {
0151     G4cout << "G4LowEnergyBremsstrahlungSpectrum is initialized"
0152            << G4endl;
0153       }
0154
0155   // Create and fill G4RDCrossSectionHandler once
0156
0157   if( crossSectionHandler != 0 ) delete crossSectionHandler;
0158   G4RDVDataSetAlgorithm* interpolation = new G4RDLogLogInterpolation();
0159   G4double lowKineticEnergy  = GetLowerBoundEloss();
0160   G4double highKineticEnergy = GetUpperBoundEloss();
0161   G4int    totBin = GetNbinEloss();
0162   crossSectionHandler = new G4RDBremsstrahlungCrossSectionHandler(energySpectrum, interpolation);
0163   crossSectionHandler->Initialise(0,lowKineticEnergy, highKineticEnergy, totBin);
0164   crossSectionHandler->LoadShellData("brem/br-cs-");
0165
0166   if (verboseLevel > 0) {
0167     G4cout << GetProcessName()
0168            << " is created; Cross section data: "
0169            << G4endl;
0170     crossSectionHandler->PrintData();
0171     G4cout << "Parameters: "
0172            << G4endl;
0173     energySpectrum->PrintData();
0174   }
0175
0176   // Build loss table for Bremsstrahlung
0177
0178   BuildLossTable(aParticleType);
0179
0180   if(verboseLevel > 0) {
0181     G4cout << "The loss table is built"
0182            << G4endl;
0183       }
0184
0185   if (&aParticleType==G4Electron::Electron()) {
0186
0187     RecorderOfElectronProcess[CounterOfElectronProcess] = (*this).theLossTable;
0188     CounterOfElectronProcess++;
0189     PrintInfoDefinition();
0190
0191   } else {
0192
0193     RecorderOfPositronProcess[CounterOfPositronProcess] = (*this).theLossTable;
0194     CounterOfPositronProcess++;
0195   }
0196
0197   // Build mean free path data using cut values
0198
0199   if( theMeanFreePath != 0 ) delete theMeanFreePath;
0200   theMeanFreePath = crossSectionHandler->
0201                     BuildMeanFreePathForMaterials(&cutForSecondaryPhotons);
0202
0203   if(verboseLevel > 0) {
0204     G4cout << "The MeanFreePath table is built"
0205            << G4endl;
0206       }
0207
0208   // Build common DEDX table for all ionisation processes
0209
0210   BuildDEDXTable(aParticleType);
0211
0212   if(verboseLevel > 0) {
0213     G4cout << "G4LowEnergyBremsstrahlung::BuildPhysicsTable end"
0214            << G4endl;
0215       }
0216
0217 }
0218
0219
0220 void G4LowEnergyBremsstrahlung::BuildLossTable(const G4ParticleDefinition& )
0221 {
0222   // Build table for energy loss due to soft brems
0223   // the tables are built for *MATERIALS* binning is taken from LowEnergyLoss
0224
0225   G4double lowKineticEnergy  = GetLowerBoundEloss();
0226   G4double highKineticEnergy = GetUpperBoundEloss();
0227   size_t totBin = GetNbinEloss();
0228
0229   //  create table
0230
0231   if (theLossTable) {
0232       theLossTable->clearAndDestroy();
0233       delete theLossTable;
0234   }
0235   const G4ProductionCutsTable* theCoupleTable=
0236         G4ProductionCutsTable::GetProductionCutsTable();
0237   size_t numOfCouples = theCoupleTable->GetTableSize();
0238   theLossTable = new G4PhysicsTable(numOfCouples);
0239
0240   // Clean up the vector of cuts
0241   cutForSecondaryPhotons.clear();
0242
0243   // Loop for materials
0244
0245   for (size_t j=0; j<numOfCouples; j++) {
0246
0247     // create physics vector and fill it
0248     G4PhysicsLogVector* aVector = new G4PhysicsLogVector(lowKineticEnergy,
0249                                  highKineticEnergy,
0250                              totBin);
0251
0252     // get material parameters needed for the energy loss calculation
0253     const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
0254     const G4Material* material= couple->GetMaterial();
0255
0256     // the cut cannot be below lowest limit
0257     G4double tCut = (*(theCoupleTable->GetEnergyCutsVector(0)))[j];
0258     tCut = std::min(highKineticEnergy, tCut);
0259     cutForSecondaryPhotons.push_back(tCut);
0260
0261     const G4ElementVector* theElementVector = material->GetElementVector();
0262     size_t NumberOfElements = material->GetNumberOfElements() ;
0263     const G4double* theAtomicNumDensityVector =
0264       material->GetAtomicNumDensityVector();
0265     if(verboseLevel > 1) {
0266       G4cout << "Energy loss for material # " << j
0267              << " tCut(keV)= " << tCut/keV
0268              << G4endl;
0269       }
0270
0271     // now comes the loop for the kinetic energy values
0272     for (size_t i = 0; i<totBin; i++) {
0273
0274       G4double lowEdgeEnergy = aVector->GetLowEdgeEnergy(i);
0275       G4double ionloss = 0.;
0276
0277       // loop for elements in the material
0278       for (size_t iel=0; iel<NumberOfElements; iel++ ) {
0279         G4int Z = (G4int)((*theElementVector)[iel]->GetZ());
0280         G4double e = energySpectrum->AverageEnergy(Z, 0.0, tCut, lowEdgeEnergy);
0281         G4double cs= crossSectionHandler->FindValue(Z, lowEdgeEnergy);
0282         ionloss   += e * cs  * theAtomicNumDensityVector[iel];
0283         if(verboseLevel > 1) {
0284           G4cout << "Z= " << Z
0285                  << "; tCut(keV)= " << tCut/keV
0286                  << "; E(keV)= " << lowEdgeEnergy/keV
0287                  << "; Eav(keV)= " << e/keV
0288                  << "; cs= " << cs
0289          << "; loss= " << ionloss
0290                  << G4endl;
0291         }
0292       }
0293       aVector->PutValue(i,ionloss);
0294     }
0295     theLossTable->insert(aVector);
0296   }
0297 }
0298
0299
0300 G4VParticleChange* G4LowEnergyBremsstrahlung::PostStepDoIt(const G4Track& track,
0301                                const G4Step& step)
0302 {
0303   aParticleChange.Initialize(track);
0304
0305   const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
0306   G4double kineticEnergy = track.GetKineticEnergy();
0307   G4int index = couple->GetIndex();
0308   G4double tCut = cutForSecondaryPhotons[index];
0309
0310   // Control limits
0311   if(tCut >= kineticEnergy)
0312      return G4VContinuousDiscreteProcess::PostStepDoIt(track, step);
0313
0314   G4int Z = crossSectionHandler->SelectRandomAtom(couple, kineticEnergy);
0315
0316   G4double tGamma = energySpectrum->SampleEnergy(Z, tCut, kineticEnergy, kineticEnergy);
0317   G4double totalEnergy = kineticEnergy + electron_mass_c2;
0318   G4double finalEnergy = kineticEnergy - tGamma; // electron/positron final energy
0319   G4double theta = 0;
0320
0321   if((kineticEnergy < 1*MeV && kineticEnergy > 1*keV && generatorName == "2bn")){
0322       theta = angularDistribution->PolarAngle(kineticEnergy,finalEnergy,Z);
0323   }else{
0324       theta = TsaiAngularDistribution->PolarAngle(kineticEnergy,finalEnergy,Z);
0325   }
0326
0327   G4double phi   = twopi * G4UniformRand();
0328   G4double dirZ  = std::cos(theta);
0329   G4double sinTheta  = std::sqrt(1. - dirZ*dirZ);
0330   G4double dirX  = sinTheta*std::cos(phi);
0331   G4double dirY  = sinTheta*std::sin(phi);
0332
0333   G4ThreeVector gammaDirection (dirX, dirY, dirZ);
0334   G4ThreeVector electronDirection = track.GetMomentumDirection();
0335
0336   //
0337   // Update the incident particle
0338   //
0339   gammaDirection.rotateUz(electronDirection);
0340
0341   // Kinematic problem
0342   if (finalEnergy < 0.) {
0343     tGamma += finalEnergy;
0344     finalEnergy = 0.0;
0345   }
0346
0347   G4double momentum = std::sqrt((totalEnergy + electron_mass_c2)*kineticEnergy);
0348
0349   G4double finalX = momentum*electronDirection.x() - tGamma*gammaDirection.x();
0350   G4double finalY = momentum*electronDirection.y() - tGamma*gammaDirection.y();
0351   G4double finalZ = momentum*electronDirection.z() - tGamma*gammaDirection.z();
0352
0353   aParticleChange.SetNumberOfSecondaries(1);
0354   G4double norm = 1./std::sqrt(finalX*finalX + finalY*finalY + finalZ*finalZ);
0355   aParticleChange.ProposeMomentumDirection(finalX*norm, finalY*norm, finalZ*norm);
0356   aParticleChange.ProposeEnergy( finalEnergy );
0357
0358   // create G4DynamicParticle object for the gamma
0359   G4DynamicParticle* aGamma= new G4DynamicParticle (G4Gamma::Gamma(),
0360                             gammaDirection, tGamma);
0361   aParticleChange.AddSecondary(aGamma);
0362
0363   return G4VContinuousDiscreteProcess::PostStepDoIt(track, step);
0364 }
0365
0366
0367 void G4LowEnergyBremsstrahlung::PrintInfoDefinition()
0368 {
0369   G4String comments = "Total cross sections from EEDL database.";
0370   comments += "\n      Gamma energy sampled from a parameterised formula.";
0371   comments += "\n      Implementation of the continuous dE/dx part.";
0372   comments += "\n      At present it can be used for electrons ";
0373   comments += "in the energy range [250eV,100GeV].";
0374   comments += "\n      The process must work with G4LowEnergyIonisation.";
0375
0376   G4cout << G4endl << GetProcessName() << ":  " << comments << G4endl;
0377 }
0378
0379 G4bool G4LowEnergyBremsstrahlung::IsApplicable(const G4ParticleDefinition& particle)
0380 {
0381   return (  (&particle == G4Electron::Electron())  );
0382 }
0383
0384
0385 G4double G4LowEnergyBremsstrahlung::GetMeanFreePath(const G4Track& track,
0386                             G4double,
0387                             G4ForceCondition* cond)
0388 {
0389   *cond = NotForced;
0390   G4int index = (track.GetMaterialCutsCouple())->GetIndex();
0391   const G4RDVEMDataSet* data = theMeanFreePath->GetComponent(index);
0392   G4double meanFreePath = data->FindValue(track.GetKineticEnergy());
0393   return meanFreePath;
0394 }
0395
0396 void G4LowEnergyBremsstrahlung::SetCutForLowEnSecPhotons(G4double cut)
0397 {
0398   cutForPhotons = cut;
0399 }
0400
0401 void G4LowEnergyBremsstrahlung::SetAngularGenerator(G4RDVBremAngularDistribution* distribution)
0402 {
0403   angularDistribution = distribution;
0404   angularDistribution->PrintGeneratorInformation();
0405 }
0406
0407 void G4LowEnergyBremsstrahlung::SetAngularGenerator(const G4String& name)
0408 {
0409   if (name == "tsai")
0410     {
0411       delete angularDistribution;
0412       angularDistribution = new G4RDModifiedTsai("TsaiGenerator");
0413       generatorName = name;
0414     }
0415   else if (name == "2bn")
0416     {
0417       delete angularDistribution;
0418       angularDistribution = new G4RDGenerator2BN("2BNGenerator");
0419       generatorName = name;
0420     }
0421   else if (name == "2bs")
0422     {
0423        delete angularDistribution;
0424        angularDistribution = new G4RDGenerator2BS("2BSGenerator");
0425        generatorName = name;
0426     }
0427   else
0428     {
0429       G4Exception("G4LowEnergyBremsstrahlung::SetAngularGenerator()",
0430                   "InvalidSetup", FatalException, "Generator does not exist!");
0431     }
0432
0433   angularDistribution->PrintGeneratorInformation();
0434 }
0435