Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/eRosita/physics/src/G4RDeBremsstrahlungSpectrum.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 // GEANT4 Class file
 30 //
 31 //
 32 // File name:     G4RDeBremsstrahlungSpectrum
 33 //
 34 // Author:        V.Ivanchenko (Vladimir.Ivanchenko@cern.ch)
 35 //
 36 // Creation date: 29 September 2001
 37 //
 38 // Modifications:
 39 // 10.10.01  MGP  Revision to improve code quality and consistency with design
 40 // 15.11.01  VI   Update spectrum model Bethe-Haitler spectrum at high energy
 41 // 30.05.02  VI   Update interpolation between 2 last energy points in the
 42 //                parametrisation
 43 // 21.02.03  V.Ivanchenko    Energy bins are defined in the constructor
 44 // 28.02.03  V.Ivanchenko    Filename is defined in the constructor
 45 //
 46 // -------------------------------------------------------------------
 47 
 48 #include "G4RDeBremsstrahlungSpectrum.hh"
 49 #include "G4RDBremsstrahlungParameters.hh"
 50 #include "Randomize.hh"
 51 #include "G4SystemOfUnits.hh"
 52 
 53 
 54 G4RDeBremsstrahlungSpectrum::G4RDeBremsstrahlungSpectrum(const G4DataVector& bins,
 55   const G4String& name):G4RDVEnergySpectrum(),
 56   lowestE(0.1*eV),
 57   xp(bins)
 58 {
 59   length = xp.size();
 60   theBRparam = new G4RDBremsstrahlungParameters(name,length+1);
 61 
 62   verbose = 0;
 63 }
 64 
 65 
 66 G4RDeBremsstrahlungSpectrum::~G4RDeBremsstrahlungSpectrum()
 67 {
 68   delete theBRparam;
 69 }
 70 
 71 
 72 G4double G4RDeBremsstrahlungSpectrum::Probability(G4int Z,
 73                                                 G4double tmin,
 74                                                 G4double tmax,
 75                                                 G4double e,
 76                                                 G4int,
 77                                        const G4ParticleDefinition*) const
 78 {
 79   G4double tm = std::min(tmax, e);
 80   G4double t0 = std::max(tmin, lowestE);
 81   if(t0 >= tm) return 0.0;
 82 
 83   t0 /= e;
 84   tm /= e;
 85 
 86   G4double z0 = lowestE/e;
 87   G4DataVector p;
 88 
 89   // Access parameters
 90   for (size_t i=0; i<=length; i++) {
 91     p.push_back(theBRparam->Parameter(i, Z, e));
 92   }
 93 
 94   G4double x  = IntSpectrum(t0, tm, p);
 95   G4double y  = IntSpectrum(z0, 1.0, p);
 96 
 97 
 98   if(1 < verbose) {
 99     G4cout << "tcut(MeV)= " << tmin/MeV
100            << "; tMax(MeV)= " << tmax/MeV
101            << "; t0= " << t0
102            << "; tm= " << tm
103            << "; xp[0]= " << xp[0]
104            << "; z= " << z0
105            << "; val= " << x
106            << "; nor= " << y
107            << G4endl;
108   }
109   p.clear();
110 
111   if(y > 0.0) x /= y;
112   else        x  = 0.0;
113   //  if(x < 0.0) x  = 0.0;
114 
115   return x;
116 }
117 
118 
119 G4double G4RDeBremsstrahlungSpectrum::AverageEnergy(G4int Z,
120                                                   G4double tmin,
121                                                   G4double tmax,
122                                                   G4double e,
123                                                   G4int,
124               const G4ParticleDefinition*) const
125 {
126   G4double tm = std::min(tmax, e);
127   G4double t0 = std::max(tmin, lowestE);
128   if(t0 >= tm) return 0.0;
129 
130   t0 /= e;
131   tm /= e;
132 
133   G4double z0 = lowestE/e;
134 
135   G4DataVector p;
136 
137   // Access parameters
138   for (size_t i=0; i<=length; i++) {
139       p.push_back(theBRparam->Parameter(i, Z, e));
140   }
141 
142   G4double x  = AverageValue(t0, tm, p);
143   G4double y  = IntSpectrum(z0, 1.0, p);
144 
145   // Add integrant over lowest energies
146   G4double zmin = tmin/e;
147   if(zmin < t0) {
148     G4double c  = std::sqrt(theBRparam->ParameterC(Z));
149     x += p[0]*(t0 - zmin - c*(std::atan(t0/c) - std::atan(zmin/c)));
150   }
151   x *= e;
152 
153   if(1 < verbose) {
154     G4cout << "tcut(MeV)= " << tmin/MeV
155            << "; tMax(MeV)= " << tmax/MeV
156            << "; e(MeV)= " << e/MeV
157            << "; t0= " << t0
158            << "; tm= " << tm
159            << "; y= " << y
160            << "; x= " << x
161            << G4endl;
162   }
163   p.clear();
164 
165   if(y > 0.0) x /= y;
166   else        x  = 0.0;
167   //  if(x < 0.0) x  = 0.0;
168 
169   return x;
170 }
171 
172 
173 G4double G4RDeBremsstrahlungSpectrum::SampleEnergy(G4int Z,
174                                                  G4double tmin,
175                                                  G4double tmax,
176                                                  G4double e,
177                                                  G4int,
178              const G4ParticleDefinition*) const
179 {
180   G4double tm = std::min(tmax, e);
181   G4double t0 = std::max(tmin, lowestE);
182   if(t0 >= tm) return 0.0;
183 
184   t0 /= e;
185   tm /= e;
186 
187   G4DataVector p;
188 
189   for (size_t i=0; i<=length; i++) {
190     p.push_back(theBRparam->Parameter(i, Z, e));
191   }
192   G4double amaj = std::max(p[length], 1. - (p[1] - p[0])*xp[0]/(xp[1] - xp[0]) );
193 
194   G4double amax = std::log(tm);
195   G4double amin = std::log(t0);
196   G4double tgam, q, fun;
197 
198   do {
199     G4double x = amin + G4UniformRand()*(amax - amin);
200     tgam = std::exp(x);
201     fun = Function(tgam, p);
202 
203     if(fun > amaj) {
204           G4cout << "WARNING in G4RDeBremsstrahlungSpectrum::SampleEnergy:"
205                  << " Majoranta " << amaj
206                  << " < " << fun
207                  << G4endl;
208     }
209 
210     q = amaj * G4UniformRand();
211   } while (q > fun);
212 
213   tgam *= e;
214 
215   p.clear();
216 
217   return tgam;
218 }
219 
220 G4double G4RDeBremsstrahlungSpectrum::IntSpectrum(G4double xMin,
221                                                 G4double xMax,
222             const G4DataVector& p) const
223 {
224   G4double x1 = std::min(xMin, xp[0]);
225   G4double x2 = std::min(xMax, xp[0]);
226   G4double sum = 0.0;
227 
228   if(x1 < x2) {
229     G4double k = (p[1] - p[0])/(xp[1] - xp[0]);
230     sum += (1. - k*xp[0])*std::log(x2/x1) + k*(x2 - x1);
231   }
232 
233   for (size_t i=0; i<length-1; i++) {
234     x1 = std::max(xMin, xp[i]);
235     x2 = std::min(xMax, xp[i+1]);
236     if(x1 < x2) {
237       G4double z1 = p[i];
238       G4double z2 = p[i+1];
239       sum += z2 - z1 + std::log(x2/x1)*(z1*x2 - z2*x1)/(x2 - x1);
240     }
241   }
242   if(sum < 0.0) sum = 0.0;
243   return sum;
244 }
245 
246 G4double G4RDeBremsstrahlungSpectrum::AverageValue(G4double xMin,
247                                                  G4double xMax,
248              const G4DataVector& p) const
249 {
250   G4double x1 = std::min(xMin, xp[0]);
251   G4double x2 = std::min(xMax, xp[0]);
252   G4double z1 = x1;
253   G4double z2 = x2;
254   G4double sum = 0.0;
255 
256   if(x1 < x2) {
257     G4double k = (p[1] - p[0])/(xp[1] - xp[0]);
258     sum += (z2 - z1)*(1. - k*xp[0]);
259     z1 *= x1;
260     z2 *= x2;
261     sum += 0.5*k*(z2 - z1);
262   }
263 
264   for (size_t i=0; i<length-1; i++) {
265     x1 = std::max(xMin, xp[i]);
266     x2 = std::min(xMax, xp[i+1]);
267     if(x1 < x2) {
268       z1 = p[i];
269       z2 = p[i+1];
270       sum += 0.5*(z2 - z1)*(x2 + x1) + z1*x2 - z2*x1;
271     }
272   }
273   if(sum < 0.0) sum = 0.0;
274   return sum;
275 }
276 
277 G4double G4RDeBremsstrahlungSpectrum::Function(G4double x,
278                const G4DataVector& p) const
279 {
280   G4double f = 0.0;
281 
282   if(x <= xp[0]) {
283     f = p[0] + (p[1] - p[0])*(x - xp[0])/(xp[1] - xp[0]);
284 
285   } else {
286 
287     for (size_t i=0; i<length-1; i++) {
288 
289       if(x <= xp[i+1]) {
290         f = p[i] + (p[i+1] - p[i])*(x - xp[i])/(xp[i+1] - xp[i]);
291         break;
292       }
293     }
294   }
295   return f;
296 }
297 
298 void G4RDeBremsstrahlungSpectrum::PrintData() const
299 { theBRparam->PrintData(); }
300 
301 G4double G4RDeBremsstrahlungSpectrum::Excitation(G4int , G4double ) const
302 {
303   return 0.0;
304 }
305 
306 G4double G4RDeBremsstrahlungSpectrum::MaxEnergyOfSecondaries(G4double kineticEnergy,
307                  G4int, // Z,
308                  const G4ParticleDefinition*) const
309 {
310   return kineticEnergy;
311 }
312