physics/src/G4RDeBremsstrahlungSpectrum.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 // GEANT4 Class file
0030 //
0031 //
0032 // File name:     G4RDeBremsstrahlungSpectrum
0033 //
0034 // Author:        V.Ivanchenko (Vladimir.Ivanchenko@cern.ch)
0035 //
0036 // Creation date: 29 September 2001
0037 //
0038 // Modifications:
0039 // 10.10.01  MGP  Revision to improve code quality and consistency with design
0040 // 15.11.01  VI   Update spectrum model Bethe-Haitler spectrum at high energy
0041 // 30.05.02  VI   Update interpolation between 2 last energy points in the
0042 //                parametrisation
0043 // 21.02.03  V.Ivanchenko    Energy bins are defined in the constructor
0044 // 28.02.03  V.Ivanchenko    Filename is defined in the constructor
0045 //
0046 // -------------------------------------------------------------------
0047
0048 #include "G4RDeBremsstrahlungSpectrum.hh"
0049 #include "G4RDBremsstrahlungParameters.hh"
0050 #include "Randomize.hh"
0051 #include "G4SystemOfUnits.hh"
0052
0053
0054 G4RDeBremsstrahlungSpectrum::G4RDeBremsstrahlungSpectrum(const G4DataVector& bins,
0055   const G4String& name):G4RDVEnergySpectrum(),
0056   lowestE(0.1*eV),
0057   xp(bins)
0058 {
0059   length = xp.size();
0060   theBRparam = new G4RDBremsstrahlungParameters(name,length+1);
0061
0062   verbose = 0;
0063 }
0064
0065
0066 G4RDeBremsstrahlungSpectrum::~G4RDeBremsstrahlungSpectrum()
0067 {
0068   delete theBRparam;
0069 }
0070
0071
0072 G4double G4RDeBremsstrahlungSpectrum::Probability(G4int Z,
0073                                                 G4double tmin,
0074                                                 G4double tmax,
0075                                                 G4double e,
0076                                                 G4int,
0077                                        const G4ParticleDefinition*) const
0078 {
0079   G4double tm = std::min(tmax, e);
0080   G4double t0 = std::max(tmin, lowestE);
0081   if(t0 >= tm) return 0.0;
0082
0083   t0 /= e;
0084   tm /= e;
0085
0086   G4double z0 = lowestE/e;
0087   G4DataVector p;
0088
0089   // Access parameters
0090   for (size_t i=0; i<=length; i++) {
0091     p.push_back(theBRparam->Parameter(i, Z, e));
0092   }
0093
0094   G4double x  = IntSpectrum(t0, tm, p);
0095   G4double y  = IntSpectrum(z0, 1.0, p);
0096
0097
0098   if(1 < verbose) {
0099     G4cout << "tcut(MeV)= " << tmin/MeV
0100            << "; tMax(MeV)= " << tmax/MeV
0101            << "; t0= " << t0
0102            << "; tm= " << tm
0103            << "; xp[0]= " << xp[0]
0104            << "; z= " << z0
0105            << "; val= " << x
0106            << "; nor= " << y
0107            << G4endl;
0108   }
0109   p.clear();
0110
0111   if(y > 0.0) x /= y;
0112   else        x  = 0.0;
0113   //  if(x < 0.0) x  = 0.0;
0114
0115   return x;
0116 }
0117
0118
0119 G4double G4RDeBremsstrahlungSpectrum::AverageEnergy(G4int Z,
0120                                                   G4double tmin,
0121                                                   G4double tmax,
0122                                                   G4double e,
0123                                                   G4int,
0124                           const G4ParticleDefinition*) const
0125 {
0126   G4double tm = std::min(tmax, e);
0127   G4double t0 = std::max(tmin, lowestE);
0128   if(t0 >= tm) return 0.0;
0129
0130   t0 /= e;
0131   tm /= e;
0132
0133   G4double z0 = lowestE/e;
0134
0135   G4DataVector p;
0136
0137   // Access parameters
0138   for (size_t i=0; i<=length; i++) {
0139       p.push_back(theBRparam->Parameter(i, Z, e));
0140   }
0141
0142   G4double x  = AverageValue(t0, tm, p);
0143   G4double y  = IntSpectrum(z0, 1.0, p);
0144
0145   // Add integrant over lowest energies
0146   G4double zmin = tmin/e;
0147   if(zmin < t0) {
0148     G4double c  = std::sqrt(theBRparam->ParameterC(Z));
0149     x += p[0]*(t0 - zmin - c*(std::atan(t0/c) - std::atan(zmin/c)));
0150   }
0151   x *= e;
0152
0153   if(1 < verbose) {
0154     G4cout << "tcut(MeV)= " << tmin/MeV
0155            << "; tMax(MeV)= " << tmax/MeV
0156            << "; e(MeV)= " << e/MeV
0157            << "; t0= " << t0
0158            << "; tm= " << tm
0159            << "; y= " << y
0160            << "; x= " << x
0161            << G4endl;
0162   }
0163   p.clear();
0164
0165   if(y > 0.0) x /= y;
0166   else        x  = 0.0;
0167   //  if(x < 0.0) x  = 0.0;
0168
0169   return x;
0170 }
0171
0172
0173 G4double G4RDeBremsstrahlungSpectrum::SampleEnergy(G4int Z,
0174                                                  G4double tmin,
0175                                                  G4double tmax,
0176                                                  G4double e,
0177                                                  G4int,
0178                          const G4ParticleDefinition*) const
0179 {
0180   G4double tm = std::min(tmax, e);
0181   G4double t0 = std::max(tmin, lowestE);
0182   if(t0 >= tm) return 0.0;
0183
0184   t0 /= e;
0185   tm /= e;
0186
0187   G4DataVector p;
0188
0189   for (size_t i=0; i<=length; i++) {
0190     p.push_back(theBRparam->Parameter(i, Z, e));
0191   }
0192   G4double amaj = std::max(p[length], 1. - (p[1] - p[0])*xp[0]/(xp[1] - xp[0]) );
0193
0194   G4double amax = std::log(tm);
0195   G4double amin = std::log(t0);
0196   G4double tgam, q, fun;
0197
0198   do {
0199     G4double x = amin + G4UniformRand()*(amax - amin);
0200     tgam = std::exp(x);
0201     fun = Function(tgam, p);
0202
0203     if(fun > amaj) {
0204           G4cout << "WARNING in G4RDeBremsstrahlungSpectrum::SampleEnergy:"
0205                  << " Majoranta " << amaj
0206                  << " < " << fun
0207                  << G4endl;
0208     }
0209
0210     q = amaj * G4UniformRand();
0211   } while (q > fun);
0212
0213   tgam *= e;
0214
0215   p.clear();
0216
0217   return tgam;
0218 }
0219
0220 G4double G4RDeBremsstrahlungSpectrum::IntSpectrum(G4double xMin,
0221                                                 G4double xMax,
0222                         const G4DataVector& p) const
0223 {
0224   G4double x1 = std::min(xMin, xp[0]);
0225   G4double x2 = std::min(xMax, xp[0]);
0226   G4double sum = 0.0;
0227
0228   if(x1 < x2) {
0229     G4double k = (p[1] - p[0])/(xp[1] - xp[0]);
0230     sum += (1. - k*xp[0])*std::log(x2/x1) + k*(x2 - x1);
0231   }
0232
0233   for (size_t i=0; i<length-1; i++) {
0234     x1 = std::max(xMin, xp[i]);
0235     x2 = std::min(xMax, xp[i+1]);
0236     if(x1 < x2) {
0237       G4double z1 = p[i];
0238       G4double z2 = p[i+1];
0239       sum += z2 - z1 + std::log(x2/x1)*(z1*x2 - z2*x1)/(x2 - x1);
0240     }
0241   }
0242   if(sum < 0.0) sum = 0.0;
0243   return sum;
0244 }
0245
0246 G4double G4RDeBremsstrahlungSpectrum::AverageValue(G4double xMin,
0247                                                  G4double xMax,
0248                          const G4DataVector& p) const
0249 {
0250   G4double x1 = std::min(xMin, xp[0]);
0251   G4double x2 = std::min(xMax, xp[0]);
0252   G4double z1 = x1;
0253   G4double z2 = x2;
0254   G4double sum = 0.0;
0255
0256   if(x1 < x2) {
0257     G4double k = (p[1] - p[0])/(xp[1] - xp[0]);
0258     sum += (z2 - z1)*(1. - k*xp[0]);
0259     z1 *= x1;
0260     z2 *= x2;
0261     sum += 0.5*k*(z2 - z1);
0262   }
0263
0264   for (size_t i=0; i<length-1; i++) {
0265     x1 = std::max(xMin, xp[i]);
0266     x2 = std::min(xMax, xp[i+1]);
0267     if(x1 < x2) {
0268       z1 = p[i];
0269       z2 = p[i+1];
0270       sum += 0.5*(z2 - z1)*(x2 + x1) + z1*x2 - z2*x1;
0271     }
0272   }
0273   if(sum < 0.0) sum = 0.0;
0274   return sum;
0275 }
0276
0277 G4double G4RDeBremsstrahlungSpectrum::Function(G4double x,
0278                          const G4DataVector& p) const
0279 {
0280   G4double f = 0.0;
0281
0282   if(x <= xp[0]) {
0283     f = p[0] + (p[1] - p[0])*(x - xp[0])/(xp[1] - xp[0]);
0284
0285   } else {
0286
0287     for (size_t i=0; i<length-1; i++) {
0288
0289       if(x <= xp[i+1]) {
0290         f = p[i] + (p[i+1] - p[i])*(x - xp[i])/(xp[i+1] - xp[i]);
0291         break;
0292       }
0293     }
0294   }
0295   return f;
0296 }
0297
0298 void G4RDeBremsstrahlungSpectrum::PrintData() const
0299 { theBRparam->PrintData(); }
0300
0301 G4double G4RDeBremsstrahlungSpectrum::Excitation(G4int , G4double ) const
0302 {
0303   return 0.0;
0304 }
0305
0306 G4double G4RDeBremsstrahlungSpectrum::MaxEnergyOfSecondaries(G4double kineticEnergy,
0307                                G4int, // Z,
0308                                const G4ParticleDefinition*) const
0309 {
0310   return kineticEnergy;
0311 }