source: trunk/source/event/src/G4SPSEneDistribution.cc @ 1347

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago

geant4 tag 9.4
File size: 38.8 KB
Line 
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// MODULE:        G4SPSEneDistribution.cc
29//
30// Version:      1.0
31// Date:         5/02/04
32// Author:       Fan Lei
33// Organisation: QinetiQ ltd.
34// Customer:     ESA/ESTEC
35//
36///////////////////////////////////////////////////////////////////////////////
37//
38// CHANGE HISTORY
39// --------------
40//
41//
42// Version 1.0, 05/02/2004, Fan Lei, Created.
43//    Based on the G4GeneralParticleSource class in Geant4 v6.0
44//
45///////////////////////////////////////////////////////////////////////////////
46//
47#include "Randomize.hh"
48//#include <cmath>
49
50#include "G4SPSEneDistribution.hh"
51
52G4SPSEneDistribution::G4SPSEneDistribution() {
53        //
54        // Initialise all variables
55        particle_energy = 1.0 * MeV;
56
57        EnergyDisType = "Mono";
58        weight = 1.;
59        MonoEnergy = 1 * MeV;
60        Emin = 0.;
61        Emax = 1.e30;
62        alpha = 0.;
63        biasalpha = 0.;
64        prob_norm = 1.0;
65        Ezero = 0.;
66        SE = 0.;
67        Temp = 0.;
68        grad = 0.;
69        cept = 0.;
70        Biased = false; // not biased
71        EnergySpec = true; // true - energy spectra, false - momentum spectra
72        DiffSpec = true; // true - differential spec, false integral spec
73        IntType = "NULL"; // Interpolation type
74        IPDFEnergyExist = false;
75        IPDFArbExist = false;
76
77        ArbEmin = 0.;
78        ArbEmax = 1.e30;
79
80        verbosityLevel = 0;
81
82}
83
84G4SPSEneDistribution::~G4SPSEneDistribution() {
85}
86
87void G4SPSEneDistribution::SetEnergyDisType(G4String DisType) {
88        EnergyDisType = DisType;
89        if (EnergyDisType == "User") {
90                UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
91                IPDFEnergyExist = false;
92        } else if (EnergyDisType == "Arb") {
93                ArbEnergyH = IPDFArbEnergyH = ZeroPhysVector;
94                IPDFArbExist = false;
95        } else if (EnergyDisType == "Epn") {
96                UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
97                IPDFEnergyExist = false;
98                EpnEnergyH = ZeroPhysVector;
99        }
100}
101
102void G4SPSEneDistribution::SetEmin(G4double emi) {
103        Emin = emi;
104}
105
106void G4SPSEneDistribution::SetEmax(G4double ema) {
107        Emax = ema;
108}
109
110void G4SPSEneDistribution::SetMonoEnergy(G4double menergy) {
111        MonoEnergy = menergy;
112}
113
114void G4SPSEneDistribution::SetBeamSigmaInE(G4double e) {
115        SE = e;
116}
117void G4SPSEneDistribution::SetAlpha(G4double alp) {
118        alpha = alp;
119}
120
121void G4SPSEneDistribution::SetBiasAlpha(G4double alp) {
122        biasalpha = alp;
123        Biased = true;
124}
125
126void G4SPSEneDistribution::SetTemp(G4double tem) {
127        Temp = tem;
128}
129
130void G4SPSEneDistribution::SetEzero(G4double eze) {
131        Ezero = eze;
132}
133
134void G4SPSEneDistribution::SetGradient(G4double gr) {
135        grad = gr;
136}
137
138void G4SPSEneDistribution::SetInterCept(G4double c) {
139        cept = c;
140}
141
142void G4SPSEneDistribution::UserEnergyHisto(G4ThreeVector input) {
143        G4double ehi, val;
144        ehi = input.x();
145        val = input.y();
146        if (verbosityLevel > 1) {
147                G4cout << "In UserEnergyHisto" << G4endl;
148                G4cout << " " << ehi << " " << val << G4endl;
149        }
150        UDefEnergyH.InsertValues(ehi, val);
151        Emax = ehi;
152}
153
154void G4SPSEneDistribution::ArbEnergyHisto(G4ThreeVector input) {
155        G4double ehi, val;
156        ehi = input.x();
157        val = input.y();
158        if (verbosityLevel > 1) {
159                G4cout << "In ArbEnergyHisto" << G4endl;
160                G4cout << " " << ehi << " " << val << G4endl;
161        }
162        ArbEnergyH.InsertValues(ehi, val);
163}
164
165void G4SPSEneDistribution::ArbEnergyHistoFile(G4String filename) {
166        std::ifstream infile(filename, std::ios::in);
167        if (!infile)
168                G4Exception("Unable to open the histo ASCII file");
169        G4double ehi, val;
170        while (infile >> ehi >> val) {
171                ArbEnergyH.InsertValues(ehi, val);
172        }
173}
174
175void G4SPSEneDistribution::EpnEnergyHisto(G4ThreeVector input) {
176        G4double ehi, val;
177        ehi = input.x();
178        val = input.y();
179        if (verbosityLevel > 1) {
180                G4cout << "In EpnEnergyHisto" << G4endl;
181                G4cout << " " << ehi << " " << val << G4endl;
182        }
183        EpnEnergyH.InsertValues(ehi, val);
184        Emax = ehi;
185        Epnflag = true;
186}
187
188void G4SPSEneDistribution::Calculate() {
189        if (EnergyDisType == "Cdg")
190                CalculateCdgSpectrum();
191        else if (EnergyDisType == "Bbody")
192                CalculateBbodySpectrum();
193}
194
195void G4SPSEneDistribution::CalculateCdgSpectrum() {
196        // This uses the spectrum from The INTEGRAL Mass Model (TIMM)
197        // to generate a Cosmic Diffuse X/gamma ray spectrum.
198        G4double pfact[2] = { 8.5, 112 };
199        G4double spind[2] = { 1.4, 2.3 };
200        G4double ene_line[3] = { 1. * keV, 18. * keV, 1E6 * keV };
201        G4int n_par;
202
203        ene_line[0] = Emin;
204        if (Emin < 18 * keV) {
205                n_par = 2;
206                ene_line[2] = Emax;
207                if (Emax < 18 * keV) {
208                        n_par = 1;
209                        ene_line[1] = Emax;
210                }
211        } else {
212                n_par = 1;
213                pfact[0] = 112.;
214                spind[0] = 2.3;
215                ene_line[1] = Emax;
216        }
217
218        // Create a cumulative histogram.
219        CDGhist[0] = 0.;
220        G4double omalpha;
221        G4int i = 0;
222
223        while (i < n_par) {
224                omalpha = 1. - spind[i];
225                CDGhist[i + 1] = CDGhist[i] + (pfact[i] / omalpha) * (std::pow(
226                                ene_line[i + 1] / keV, omalpha) - std::pow(ene_line[i] / keV,
227                                omalpha));
228                i++;
229        }
230
231        // Normalise histo and divide by 1000 to make MeV.
232        i = 0;
233        while (i < n_par) {
234                CDGhist[i + 1] = CDGhist[i + 1] / CDGhist[n_par];
235                //      G4cout << CDGhist[i] << CDGhist[n_par] << G4endl;
236                i++;
237        }
238}
239
240void G4SPSEneDistribution::CalculateBbodySpectrum() {
241        // create bbody spectrum
242        // Proved very hard to integrate indefinitely, so different
243        // method. User inputs emin, emax and T. These are used to
244        // create a 10,000 bin histogram.
245        // Use photon density spectrum = 2 nu**2/c**2 * (std::exp(h nu/kT)-1)
246        // = 2 E**2/h**2c**2 times the exponential
247        G4double erange = Emax - Emin;
248        G4double steps = erange / 10000.;
249        G4double Bbody_y[10000];
250        G4double k = 8.6181e-11; //Boltzmann const in MeV/K
251        G4double h = 4.1362e-21; // Plancks const in MeV s
252        G4double c = 3e8; // Speed of light
253        G4double h2 = h * h;
254        G4double c2 = c * c;
255        G4int count = 0;
256        G4double sum = 0.;
257        BBHist[0] = 0.;
258        while (count < 10000) {
259                Bbody_x[count] = Emin + G4double(count * steps);
260                Bbody_y[count] = (2. * std::pow(Bbody_x[count], 2.)) / (h2 * c2
261                                * (std::exp(Bbody_x[count] / (k * Temp)) - 1.));
262                sum = sum + Bbody_y[count];
263                BBHist[count + 1] = BBHist[count] + Bbody_y[count];
264                count++;
265        }
266
267        Bbody_x[10000] = Emax;
268        // Normalise cumulative histo.
269        count = 0;
270        while (count < 10001) {
271                BBHist[count] = BBHist[count] / sum;
272                count++;
273        }
274}
275
276void G4SPSEneDistribution::InputEnergySpectra(G4bool value) {
277        // Allows user to specifiy spectrum is momentum
278        EnergySpec = value; // false if momentum
279        if (verbosityLevel > 1)
280                G4cout << "EnergySpec has value " << EnergySpec << G4endl;
281}
282
283void G4SPSEneDistribution::InputDifferentialSpectra(G4bool value) {
284        // Allows user to specify integral or differential spectra
285        DiffSpec = value; // true = differential, false = integral
286        if (verbosityLevel > 1)
287                G4cout << "Diffspec has value " << DiffSpec << G4endl;
288}
289
290void G4SPSEneDistribution::ArbInterpolate(G4String IType) {
291        if (EnergyDisType != "Arb")
292                G4cout << "Error: this is for arbitrary distributions" << G4endl;
293        IntType = IType;
294        ArbEmax = ArbEnergyH.GetMaxLowEdgeEnergy();
295        ArbEmin = ArbEnergyH.GetMinLowEdgeEnergy();
296
297        // Now interpolate points
298        if (IntType == "Lin")
299                LinearInterpolation();
300        if (IntType == "Log")
301                LogInterpolation();
302        if (IntType == "Exp")
303                ExpInterpolation();
304        if (IntType == "Spline")
305                SplineInterpolation();
306}
307
308void G4SPSEneDistribution::LinearInterpolation() {
309        // Method to do linear interpolation on the Arb points
310        // Calculate equation of each line segment, max 1024.
311        // Calculate Area under each segment
312        // Create a cumulative array which is then normalised Arb_Cum_Area
313
314        G4double Area_seg[1024]; // Stores area under each segment
315        G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
316        G4int i, count;
317        G4int maxi = ArbEnergyH.GetVectorLength();
318        for (i = 0; i < maxi; i++) {
319                Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
320                Arb_y[i] = ArbEnergyH(size_t(i));
321        }
322        // Points are now in x,y arrays. If the spectrum is integral it has to be
323        // made differential and if momentum it has to be made energy.
324        if (DiffSpec == false) {
325                // Converts integral point-wise spectra to Differential
326                for (count = 0; count < maxi - 1; count++) {
327                        Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
328                                        / (Arb_x[count + 1] - Arb_x[count]);
329                }
330                maxi--;
331        }
332        //
333        if (EnergySpec == false) {
334                // change currently stored values (emin etc) which are actually momenta
335                // to energies.
336                if (particle_definition == NULL)
337                        G4cout << "Error: particle not defined" << G4endl;
338                else {
339                        // Apply Energy**2 = p**2c**2 + m0**2c**4
340                        // p should be entered as E/c i.e. without the division by c
341                        // being done - energy equivalent.
342                        G4double mass = particle_definition->GetPDGMass();
343                        // convert point to energy unit and its value to per energy unit
344                        G4double total_energy;
345                        for (count = 0; count < maxi; count++) {
346                                total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
347                                                * mass)); // total energy
348
349                                Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
350                                Arb_x[count] = total_energy - mass; // kinetic energy
351                        }
352                }
353        }
354        //
355        i = 1;
356        Arb_grad[0] = 0.;
357        Arb_cept[0] = 0.;
358        Area_seg[0] = 0.;
359        Arb_Cum_Area[0] = 0.;
360        while (i < maxi) {
361                // calc gradient and intercept for each segment
362                Arb_grad[i] = (Arb_y[i] - Arb_y[i - 1]) / (Arb_x[i] - Arb_x[i - 1]);
363                if (verbosityLevel == 2)
364                        G4cout << Arb_grad[i] << G4endl;
365                if (Arb_grad[i] > 0.) {
366                        if (verbosityLevel == 2)
367                                G4cout << "Arb_grad is positive" << G4endl;
368                        Arb_cept[i] = Arb_y[i] - (Arb_grad[i] * Arb_x[i]);
369                } else if (Arb_grad[i] < 0.) {
370                        if (verbosityLevel == 2)
371                                G4cout << "Arb_grad is negative" << G4endl;
372                        Arb_cept[i] = Arb_y[i] + (-Arb_grad[i] * Arb_x[i]);
373                } else {
374                        if (verbosityLevel == 2)
375                                G4cout << "Arb_grad is 0." << G4endl;
376                        Arb_cept[i] = Arb_y[i];
377                }
378
379                Area_seg[i] = ((Arb_grad[i] / 2) * (Arb_x[i] * Arb_x[i] - Arb_x[i - 1]
380                                * Arb_x[i - 1]) + Arb_cept[i] * (Arb_x[i] - Arb_x[i - 1]));
381                Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
382                sum = sum + Area_seg[i];
383                if (verbosityLevel == 2)
384                        G4cout << Arb_x[i] << Arb_y[i] << Area_seg[i] << sum << Arb_grad[i]
385                                        << G4endl;
386                i++;
387        }
388
389        i = 0;
390        while (i < maxi) {
391                Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum; // normalisation
392                IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
393                i++;
394        }
395
396        // now scale the ArbEnergyH, needed by Probability()
397        ArbEnergyH.ScaleVector(1., 1./sum);
398
399        if (verbosityLevel >= 1) {
400                G4cout << "Leaving LinearInterpolation" << G4endl;
401                ArbEnergyH.DumpValues();
402                IPDFArbEnergyH.DumpValues();
403        }
404}
405
406void G4SPSEneDistribution::LogInterpolation() {
407        // Interpolation based on Logarithmic equations
408        // Generate equations of line segments
409        // y = Ax**alpha => log y = alpha*logx + logA
410        // Find area under line segments
411        // create normalised, cumulative array Arb_Cum_Area
412        G4double Area_seg[1024]; // Stores area under each segment
413        G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
414        G4int i, count;
415        G4int maxi = ArbEnergyH.GetVectorLength();
416        for (i = 0; i < maxi; i++) {
417                Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
418                Arb_y[i] = ArbEnergyH(size_t(i));
419        }
420        // Points are now in x,y arrays. If the spectrum is integral it has to be
421        // made differential and if momentum it has to be made energy.
422        if (DiffSpec == false) {
423                // Converts integral point-wise spectra to Differential
424                for (count = 0; count < maxi - 1; count++) {
425                        Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
426                                        / (Arb_x[count + 1] - Arb_x[count]);
427                }
428                maxi--;
429        }
430        //
431        if (EnergySpec == false) {
432                // change currently stored values (emin etc) which are actually momenta
433                // to energies.
434                if (particle_definition == NULL)
435                        G4cout << "Error: particle not defined" << G4endl;
436                else {
437                        // Apply Energy**2 = p**2c**2 + m0**2c**4
438                        // p should be entered as E/c i.e. without the division by c
439                        // being done - energy equivalent.
440                        G4double mass = particle_definition->GetPDGMass();
441                        // convert point to energy unit and its value to per energy unit
442                        G4double total_energy;
443                        for (count = 0; count < maxi; count++) {
444                                total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
445                                                * mass)); // total energy
446
447                                Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
448                                Arb_x[count] = total_energy - mass; // kinetic energy
449                        }
450                }
451        }
452        //
453        i = 1;
454        Arb_alpha[0] = 0.;
455        Arb_Const[0] = 0.;
456        Area_seg[0] = 0.;
457        Arb_Cum_Area[0] = 0.;
458        if (Arb_x[0] <= 0. || Arb_y[0] <= 0.) {
459                G4cout << "You should not use log interpolation with points <= 0."
460                                << G4endl;
461                G4cout << "These will be changed to 1e-20, which may cause problems"
462                                << G4endl;
463                if (Arb_x[0] <= 0.)
464                        Arb_x[0] = 1e-20;
465                if (Arb_y[0] <= 0.)
466                        Arb_y[0] = 1e-20;
467        }
468
469        G4double alp;
470        while (i < maxi) {
471                // Incase points are negative or zero
472                if (Arb_x[i] <= 0. || Arb_y[i] <= 0.) {
473                        G4cout << "You should not use log interpolation with points <= 0."
474                                        << G4endl;
475                        G4cout
476                                        << "These will be changed to 1e-20, which may cause problems"
477                                        << G4endl;
478                        if (Arb_x[i] <= 0.)
479                                Arb_x[i] = 1e-20;
480                        if (Arb_y[i] <= 0.)
481                                Arb_y[i] = 1e-20;
482                }
483
484                Arb_alpha[i] = (std::log10(Arb_y[i]) - std::log10(Arb_y[i - 1]))
485                                / (std::log10(Arb_x[i]) - std::log10(Arb_x[i - 1]));
486                Arb_Const[i] = Arb_y[i] / (std::pow(Arb_x[i], Arb_alpha[i]));
487                alp = Arb_alpha[i] + 1;
488                if (alp == 0.) {
489                  Area_seg[i] = Arb_Const[i] * (std::log(Arb_x[i]) - std::log(Arb_x[i - 1])); 
490                } else {
491                  Area_seg[i] = (Arb_Const[i] / alp) * (std::pow(Arb_x[i], alp)
492                                - std::pow(Arb_x[i - 1], alp));
493                }
494                sum = sum + Area_seg[i];
495                Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
496                if (verbosityLevel == 2)
497                        G4cout << Arb_alpha[i] << Arb_Const[i] << Area_seg[i] << G4endl;
498                i++;
499        }
500
501        i = 0;
502        while (i < maxi) {
503                Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum;
504                IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
505                i++;
506        }
507
508        // now scale the ArbEnergyH, needed by Probability()
509        ArbEnergyH.ScaleVector(1., 1./sum);
510
511        if (verbosityLevel >= 1)
512                G4cout << "Leaving LogInterpolation " << G4endl;
513}
514
515void G4SPSEneDistribution::ExpInterpolation() {
516        // Interpolation based on Exponential equations
517        // Generate equations of line segments
518        // y = Ae**-(x/e0) => ln y = -x/e0 + lnA
519        // Find area under line segments
520        // create normalised, cumulative array Arb_Cum_Area
521        G4double Area_seg[1024]; // Stores area under each segment
522        G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
523        G4int i, count;
524        G4int maxi = ArbEnergyH.GetVectorLength();
525        for (i = 0; i < maxi; i++) {
526                Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
527                Arb_y[i] = ArbEnergyH(size_t(i));
528        }
529        // Points are now in x,y arrays. If the spectrum is integral it has to be
530        // made differential and if momentum it has to be made energy.
531        if (DiffSpec == false) {
532                // Converts integral point-wise spectra to Differential
533                for (count = 0; count < maxi - 1; count++) {
534                        Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
535                                        / (Arb_x[count + 1] - Arb_x[count]);
536                }
537                maxi--;
538        }
539        //
540        if (EnergySpec == false) {
541                // change currently stored values (emin etc) which are actually momenta
542                // to energies.
543                if (particle_definition == NULL)
544                        G4cout << "Error: particle not defined" << G4endl;
545                else {
546                        // Apply Energy**2 = p**2c**2 + m0**2c**4
547                        // p should be entered as E/c i.e. without the division by c
548                        // being done - energy equivalent.
549                        G4double mass = particle_definition->GetPDGMass();
550                        // convert point to energy unit and its value to per energy unit
551                        G4double total_energy;
552                        for (count = 0; count < maxi; count++) {
553                                total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
554                                                * mass)); // total energy
555
556                                Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
557                                Arb_x[count] = total_energy - mass; // kinetic energy
558                        }
559                }
560        }
561        //
562        i = 1;
563        Arb_ezero[0] = 0.;
564        Arb_Const[0] = 0.;
565        Area_seg[0] = 0.;
566        Arb_Cum_Area[0] = 0.;
567        while (i < maxi) {
568                G4double test = std::log(Arb_y[i]) - std::log(Arb_y[i - 1]);
569                if (test > 0. || test < 0.) {
570                        Arb_ezero[i] = -(Arb_x[i] - Arb_x[i - 1]) / (std::log(Arb_y[i])
571                                        - std::log(Arb_y[i - 1]));
572                        Arb_Const[i] = Arb_y[i] / (std::exp(-Arb_x[i] / Arb_ezero[i]));
573                        Area_seg[i] = -(Arb_Const[i] * Arb_ezero[i]) * (std::exp(-Arb_x[i]
574                                        / Arb_ezero[i]) - std::exp(-Arb_x[i - 1] / Arb_ezero[i]));
575                } else {
576                        G4cout << "Flat line segment: problem" << G4endl;
577                        Arb_ezero[i] = 0.;
578                        Arb_Const[i] = 0.;
579                        Area_seg[i] = 0.;
580                }
581                sum = sum + Area_seg[i];
582                Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Area_seg[i];
583                if (verbosityLevel == 2)
584                        G4cout << Arb_ezero[i] << Arb_Const[i] << Area_seg[i] << G4endl;
585                i++;
586        }
587
588        i = 0;
589        while (i < maxi) {
590                Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum;
591                IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
592                i++;
593        }
594
595        // now scale the ArbEnergyH, needed by Probability()
596        ArbEnergyH.ScaleVector(1., 1./sum);
597
598        if (verbosityLevel >= 1)
599                G4cout << "Leaving ExpInterpolation " << G4endl;
600}
601
602void G4SPSEneDistribution::SplineInterpolation() {
603        // Interpolation using Splines.
604        // Create Normalised arrays, make x 0->1 and y hold
605        // the function (Energy)
606        //
607        // Current method based on the above will not work in all cases.
608        // new method is implemented below.
609 
610        G4double sum, Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
611        G4int i, count;
612
613        G4int maxi = ArbEnergyH.GetVectorLength();
614        for (i = 0; i < maxi; i++) {
615                Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
616                Arb_y[i] = ArbEnergyH(size_t(i));
617        }
618        // Points are now in x,y arrays. If the spectrum is integral it has to be
619        // made differential and if momentum it has to be made energy.
620        if (DiffSpec == false) {
621                // Converts integral point-wise spectra to Differential
622                for (count = 0; count < maxi - 1; count++) {
623                        Arb_y[count] = (Arb_y[count] - Arb_y[count + 1])
624                                        / (Arb_x[count + 1] - Arb_x[count]);
625                }
626                maxi--;
627        }
628        //
629        if (EnergySpec == false) {
630                // change currently stored values (emin etc) which are actually momenta
631                // to energies.
632                if (particle_definition == NULL)
633                        G4cout << "Error: particle not defined" << G4endl;
634                else {
635                        // Apply Energy**2 = p**2c**2 + m0**2c**4
636                        // p should be entered as E/c i.e. without the division by c
637                        // being done - energy equivalent.
638                        G4double mass = particle_definition->GetPDGMass();
639                        // convert point to energy unit and its value to per energy unit
640                        G4double total_energy;
641                        for (count = 0; count < maxi; count++) {
642                                total_energy = std::sqrt((Arb_x[count] * Arb_x[count]) + (mass
643                                                * mass)); // total energy
644
645                                Arb_y[count] = Arb_y[count] * Arb_x[count] / total_energy;
646                                Arb_x[count] = total_energy - mass; // kinetic energy
647                        }
648                }
649        }
650
651        //
652        i = 1;
653        Arb_Cum_Area[0] = 0.;
654        sum = 0.;
655        Splinetemp = new G4DataInterpolation(Arb_x, Arb_y, maxi, 0., 0.);
656        G4double ei[101],prob[101];
657        while (i < maxi) {
658          // 100 step per segment for the integration of area
659          G4double de = (Arb_x[i] - Arb_x[i - 1])/100.;
660          G4double area = 0.;
661
662          for (count = 0; count < 101; count++) {
663            ei[count] = Arb_x[i - 1] + de*count ;
664            prob[count] =  Splinetemp->CubicSplineInterpolation(ei[count]);
665            if (prob[count] < 0.) { 
666              G4cout <<   "Warning: G4DataInterpolation returns value < 0  " << prob[count] <<" "<<ei[count]<< G4endl;
667              G4Exception("         Please use an alternative method, e.g. Lin, for interpolation");
668            }
669            area += prob[count]*de;
670          }
671          Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + area;
672          sum += area; 
673
674          prob[0] = prob[0]/(area/de);
675          for (count = 1; count < 100; count++)
676            prob[count] = prob[count-1] + prob[count]/(area/de);
677
678          SplineInt[i] = new G4DataInterpolation(prob, ei, 101, 0., 0.);
679          // note i start from 1!
680          i++;
681        }
682        i = 0;
683        while (i < maxi) {
684                Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum; // normalisation
685                IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
686                i++;
687        }
688        // now scale the ArbEnergyH, needed by Probability()
689        ArbEnergyH.ScaleVector(1., 1./sum);
690
691        if (verbosityLevel > 0)
692          G4cout << "Leaving SplineInterpolation " << G4endl;
693}
694
695void G4SPSEneDistribution::GenerateMonoEnergetic() {
696        // Method to generate MonoEnergetic particles.
697        particle_energy = MonoEnergy;
698}
699
700void G4SPSEneDistribution::GenerateGaussEnergies() {
701        // Method to generate Gaussian particles.
702        particle_energy = G4RandGauss::shoot(MonoEnergy,SE);
703        if (particle_energy < 0) particle_energy = 0.;
704}
705
706void G4SPSEneDistribution::GenerateLinearEnergies(G4bool bArb = false) {
707        G4double rndm;
708        G4double emaxsq = std::pow(Emax, 2.); //Emax squared
709        G4double eminsq = std::pow(Emin, 2.); //Emin squared
710        G4double intersq = std::pow(cept, 2.); //cept squared
711
712        if (bArb)
713                rndm = G4UniformRand();
714        else
715                rndm = eneRndm->GenRandEnergy();
716
717        G4double bracket = ((grad / 2.) * (emaxsq - eminsq) + cept * (Emax - Emin));
718        bracket = bracket * rndm;
719        bracket = bracket + (grad / 2.) * eminsq + cept * Emin;
720        // Now have a quad of form m/2 E**2 + cE - bracket = 0
721        bracket = -bracket;
722        //  G4cout << "BRACKET" << bracket << G4endl;
723        if (grad != 0.) {
724                G4double sqbrack = (intersq - 4 * (grad / 2.) * (bracket));
725                //      G4cout << "SQBRACK" << sqbrack << G4endl;
726                sqbrack = std::sqrt(sqbrack);
727                G4double root1 = -cept + sqbrack;
728                root1 = root1 / (2. * (grad / 2.));
729
730                G4double root2 = -cept - sqbrack;
731                root2 = root2 / (2. * (grad / 2.));
732
733                //      G4cout << root1 << " roots " << root2 << G4endl;
734
735                if (root1 > Emin && root1 < Emax)
736                        particle_energy = root1;
737                if (root2 > Emin && root2 < Emax)
738                        particle_energy = root2;
739        } else if (grad == 0.)
740                // have equation of form cE - bracket =0
741                particle_energy = bracket / cept;
742
743        if (particle_energy < 0.)
744                particle_energy = -particle_energy;
745
746        if (verbosityLevel >= 1)
747                G4cout << "Energy is " << particle_energy << G4endl;
748}
749
750void G4SPSEneDistribution::GeneratePowEnergies(G4bool bArb = false) {
751        // Method to generate particle energies distributed as
752        // a power-law
753
754        G4double rndm;
755        G4double emina, emaxa;
756
757        emina = std::pow(Emin, alpha + 1);
758        emaxa = std::pow(Emax, alpha + 1);
759
760        if (bArb)
761                rndm = G4UniformRand();
762        else
763                rndm = eneRndm->GenRandEnergy();
764
765        if (alpha != -1.) {
766                particle_energy = ((rndm * (emaxa - emina)) + emina);
767                particle_energy = std::pow(particle_energy, (1. / (alpha + 1.)));
768        } else {
769                particle_energy = (std::log(Emin) + rndm * (std::log(Emax) - std::log(
770                                Emin)));
771                particle_energy = std::exp(particle_energy);
772        }
773        if (verbosityLevel >= 1)
774                G4cout << "Energy is " << particle_energy << G4endl;
775}
776
777void G4SPSEneDistribution::GenerateBiasPowEnergies() {
778        // Method to generate particle energies distributed as
779        // in biased power-law and calculate its weight
780
781        G4double rndm;
782        G4double emina, emaxa, emin, emax;
783
784        G4double normal = 1. ;
785
786        emin = Emin;
787        emax = Emax;
788        //      if (EnergyDisType == "Arb") {
789        //  emin = ArbEmin;
790        //  emax = ArbEmax;
791        //}
792
793        rndm = eneRndm->GenRandEnergy();
794
795        if (biasalpha != -1.) {
796                emina = std::pow(emin, biasalpha + 1);
797                emaxa = std::pow(emax, biasalpha + 1);
798                particle_energy = ((rndm * (emaxa - emina)) + emina);
799                particle_energy = std::pow(particle_energy, (1. / (biasalpha + 1.)));
800                normal = 1./(1+biasalpha) * (emaxa - emina);
801        } else {
802                particle_energy = (std::log(emin) + rndm * (std::log(emax) - std::log(
803                                emin)));
804                particle_energy = std::exp(particle_energy);
805                normal = std::log(emax) - std::log(emin) ;
806        }
807        weight = GetProbability(particle_energy) / (std::pow(particle_energy,biasalpha)/normal);
808
809        if (verbosityLevel >= 1)
810                G4cout << "Energy is " << particle_energy << G4endl;
811}
812
813void G4SPSEneDistribution::GenerateExpEnergies(G4bool bArb = false) {
814        // Method to generate particle energies distributed according
815        // to an exponential curve.
816        G4double rndm;
817
818        if (bArb)
819                rndm = G4UniformRand();
820        else
821                rndm = eneRndm->GenRandEnergy();
822
823        particle_energy = -Ezero * (std::log(rndm * (std::exp(-Emax / Ezero)
824                        - std::exp(-Emin / Ezero)) + std::exp(-Emin / Ezero)));
825        if (verbosityLevel >= 1)
826                G4cout << "Energy is " << particle_energy << G4endl;
827}
828
829void G4SPSEneDistribution::GenerateBremEnergies() {
830        // Method to generate particle energies distributed according
831        // to a Bremstrahlung equation of
832        // form I = const*((kT)**1/2)*E*(e**(-E/kT))
833
834        G4double rndm;
835        rndm = eneRndm->GenRandEnergy();
836        G4double expmax, expmin, k;
837
838        k = 8.6181e-11; // Boltzmann's const in MeV/K
839        G4double ksq = std::pow(k, 2.); // k squared
840        G4double Tsq = std::pow(Temp, 2.); // Temp squared
841
842        expmax = std::exp(-Emax / (k * Temp));
843        expmin = std::exp(-Emin / (k * Temp));
844
845        // If either expmax or expmin are zero then this will cause problems
846        // Most probably this will be because T is too low or E is too high
847
848        if (expmax == 0.)
849                G4cout << "*****EXPMAX=0. Choose different E's or Temp" << G4endl;
850        if (expmin == 0.)
851                G4cout << "*****EXPMIN=0. Choose different E's or Temp" << G4endl;
852
853        G4double tempvar = rndm * ((-k) * Temp * (Emax * expmax - Emin * expmin)
854                        - (ksq * Tsq * (expmax - expmin)));
855
856        G4double bigc = (tempvar - k * Temp * Emin * expmin - ksq * Tsq * expmin)
857                        / (-k * Temp);
858
859        // This gives an equation of form: Ee(-E/kT) + kTe(-E/kT) - C =0
860        // Solve this iteratively, step from Emin to Emax in 1000 steps
861        // and take the best solution.
862
863        G4double erange = Emax - Emin;
864        G4double steps = erange / 1000.;
865        G4int i;
866        G4double etest, diff, err;
867
868        err = 100000.;
869
870        for (i = 1; i < 1000; i++) {
871                etest = Emin + (i - 1) * steps;
872
873                diff = etest * (std::exp(-etest / (k * Temp))) + k * Temp * (std::exp(
874                                -etest / (k * Temp))) - bigc;
875
876                if (diff < 0.)
877                        diff = -diff;
878
879                if (diff < err) {
880                        err = diff;
881                        particle_energy = etest;
882                }
883        }
884        if (verbosityLevel >= 1)
885                G4cout << "Energy is " << particle_energy << G4endl;
886}
887
888void G4SPSEneDistribution::GenerateBbodyEnergies() {
889        // BBody_x holds Energies, and BBHist holds the cumulative histo.
890        // binary search to find correct bin then lin interpolation.
891        // Use the earlier defined histogram + RandGeneral method to generate
892        // random numbers following the histos distribution.
893        G4double rndm;
894        G4int nabove, nbelow = 0, middle;
895        nabove = 10001;
896        rndm = eneRndm->GenRandEnergy();
897
898        // Binary search to find bin that rndm is in
899        while (nabove - nbelow > 1) {
900                middle = (nabove + nbelow) / 2;
901                if (rndm == BBHist[middle])
902                        break;
903                if (rndm < BBHist[middle])
904                        nabove = middle;
905                else
906                        nbelow = middle;
907        }
908
909        // Now interpolate in that bin to find the correct output value.
910        G4double x1, x2, y1, y2, m, q;
911        x1 = Bbody_x[nbelow];
912        x2 = Bbody_x[nbelow + 1];
913        y1 = BBHist[nbelow];
914        y2 = BBHist[nbelow + 1];
915        m = (y2 - y1) / (x2 - x1);
916        q = y1 - m * x1;
917
918        particle_energy = (rndm - q) / m;
919
920        if (verbosityLevel >= 1) {
921                G4cout << "Energy is " << particle_energy << G4endl;
922        }
923}
924
925void G4SPSEneDistribution::GenerateCdgEnergies() {
926        // Gen random numbers, compare with values in cumhist
927        // to find appropriate part of spectrum and then
928        // generate energy in the usual inversion way.
929        //  G4double pfact[2] = {8.5, 112};
930        // G4double spind[2] = {1.4, 2.3};
931        // G4double ene_line[3] = {1., 18., 1E6};
932        G4double rndm, rndm2;
933        G4double ene_line[3];
934        G4double omalpha[2];
935        if (Emin < 18 * keV && Emax < 18 * keV) {
936                omalpha[0] = 1. - 1.4;
937                ene_line[0] = Emin;
938                ene_line[1] = Emax;
939        }
940        if (Emin < 18 * keV && Emax > 18 * keV) {
941                omalpha[0] = 1. - 1.4;
942                omalpha[1] = 1. - 2.3;
943                ene_line[0] = Emin;
944                ene_line[1] = 18. * keV;
945                ene_line[2] = Emax;
946        }
947        if (Emin > 18 * keV) {
948                omalpha[0] = 1. - 2.3;
949                ene_line[0] = Emin;
950                ene_line[1] = Emax;
951        }
952        rndm = eneRndm->GenRandEnergy();
953        rndm2 = eneRndm->GenRandEnergy();
954
955        G4int i = 0;
956        while (rndm >= CDGhist[i]) {
957                i++;
958        }
959        // Generate final energy.
960        particle_energy = (std::pow(ene_line[i - 1], omalpha[i - 1]) + (std::pow(
961                        ene_line[i], omalpha[i - 1]) - std::pow(ene_line[i - 1], omalpha[i
962                        - 1])) * rndm2);
963        particle_energy = std::pow(particle_energy, (1. / omalpha[i - 1]));
964
965        if (verbosityLevel >= 1)
966                G4cout << "Energy is " << particle_energy << G4endl;
967}
968
969void G4SPSEneDistribution::GenUserHistEnergies() {
970        // Histograms are DIFFERENTIAL.
971        //  G4cout << "In GenUserHistEnergies " << G4endl;
972        if (IPDFEnergyExist == false) {
973                G4int ii;
974                G4int maxbin = G4int(UDefEnergyH.GetVectorLength());
975                G4double bins[1024], vals[1024], sum;
976                sum = 0.;
977
978                if ((EnergySpec == false) && (particle_definition == NULL))
979                        G4cout << "Error: particle definition is NULL" << G4endl;
980
981                if (maxbin > 1024) {
982                        G4cout << "Maxbin > 1024" << G4endl;
983                        G4cout << "Setting maxbin to 1024, other bins are lost" << G4endl;
984                }
985
986                if (DiffSpec == false)
987                        G4cout << "Histograms are Differential!!! " << G4endl;
988                else {
989                        bins[0] = UDefEnergyH.GetLowEdgeEnergy(size_t(0));
990                        vals[0] = UDefEnergyH(size_t(0));
991                        sum = vals[0];
992                        for (ii = 1; ii < maxbin; ii++) {
993                                bins[ii] = UDefEnergyH.GetLowEdgeEnergy(size_t(ii));
994                                vals[ii] = UDefEnergyH(size_t(ii)) + vals[ii - 1];
995                                sum = sum + UDefEnergyH(size_t(ii));
996                        }
997                }
998
999                if (EnergySpec == false) {
1000                        G4double mass = particle_definition->GetPDGMass();
1001                        // multiply the function (vals) up by the bin width
1002                        // to make the function counts/s (i.e. get rid of momentum
1003                        // dependence).
1004                        for (ii = 1; ii < maxbin; ii++) {
1005                                vals[ii] = vals[ii] * (bins[ii] - bins[ii - 1]);
1006                        }
1007                        // Put energy bins into new histo, plus divide by energy bin width
1008                        // to make evals counts/s/energy
1009                        for (ii = 0; ii < maxbin; ii++) {
1010                                bins[ii] = std::sqrt((bins[ii] * bins[ii]) + (mass * mass))
1011                                                - mass; //kinetic energy
1012                        }
1013                        for (ii = 1; ii < maxbin; ii++) {
1014                                vals[ii] = vals[ii] / (bins[ii] - bins[ii - 1]);
1015                        }
1016                        sum = vals[maxbin - 1];
1017                        vals[0] = 0.;
1018                }
1019                for (ii = 0; ii < maxbin; ii++) {
1020                        vals[ii] = vals[ii] / sum;
1021                        IPDFEnergyH.InsertValues(bins[ii], vals[ii]);
1022                }
1023
1024                // Make IPDFEnergyExist = true
1025                IPDFEnergyExist = true;
1026                if (verbosityLevel > 1)
1027                        IPDFEnergyH.DumpValues();
1028        }
1029
1030        // IPDF has been create so carry on
1031        G4double rndm = eneRndm->GenRandEnergy();
1032        particle_energy = IPDFEnergyH.GetEnergy(rndm);
1033
1034        if (verbosityLevel >= 1)
1035                G4cout << "Energy is " << particle_energy << G4endl;
1036}
1037
1038void G4SPSEneDistribution::GenArbPointEnergies() {
1039        if (verbosityLevel > 0)
1040                G4cout << "In GenArbPointEnergies" << G4endl;
1041        G4double rndm;
1042        rndm = eneRndm->GenRandEnergy();
1043        //      IPDFArbEnergyH.DumpValues();
1044        // Find the Bin
1045        // have x, y, no of points, and cumulative area distribution
1046        G4int nabove, nbelow = 0, middle;
1047        nabove = IPDFArbEnergyH.GetVectorLength();
1048        //      G4cout << nabove << G4endl;
1049        // Binary search to find bin that rndm is in
1050        while (nabove - nbelow > 1) {
1051          middle = (nabove + nbelow) / 2;
1052          if (rndm == IPDFArbEnergyH(size_t(middle)))
1053            break;
1054          if (rndm < IPDFArbEnergyH(size_t(middle)))
1055            nabove = middle;
1056          else
1057            nbelow = middle;
1058        }
1059        if (IntType == "Lin") {
1060          Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1));
1061          Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
1062          grad = Arb_grad[nbelow + 1];
1063          cept = Arb_cept[nbelow + 1];
1064          //      G4cout << rndm << " " << Emax << " " << Emin << " " << grad << " " << cept << G4endl;
1065          GenerateLinearEnergies(true);
1066        } else if (IntType == "Log") {
1067          Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1));
1068          Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
1069          alpha = Arb_alpha[nbelow + 1];
1070          //      G4cout << rndm << " " << Emax << " " << Emin << " " << alpha << G4endl;
1071          GeneratePowEnergies(true);
1072        } else if (IntType == "Exp") {
1073          Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1));
1074          Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
1075          Ezero = Arb_ezero[nbelow + 1];
1076          //      G4cout << rndm << " " << Emax << " " << Emin << " " << Ezero << G4endl;
1077          GenerateExpEnergies(true);
1078        } else if (IntType == "Spline") {
1079          Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow + 1));
1080          Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
1081          particle_energy = -1e100;
1082          rndm = eneRndm->GenRandEnergy();
1083          while (particle_energy < Emin || particle_energy > Emax) {
1084            particle_energy = SplineInt[nbelow+1]->CubicSplineInterpolation(rndm);
1085            rndm = eneRndm->GenRandEnergy();
1086          }
1087          if (verbosityLevel >= 1)
1088            G4cout << "Energy is " << particle_energy << G4endl;
1089        } else
1090                G4cout << "Error: IntType unknown type" << G4endl;
1091}
1092
1093void G4SPSEneDistribution::GenEpnHistEnergies() {
1094        //  G4cout << "In GenEpnHistEnergies " << Epnflag << G4endl;
1095
1096        // Firstly convert to energy if not already done.
1097        if (Epnflag == true)
1098        // epnflag = true means spectrum is epn, false means e.
1099        {
1100                // convert to energy by multiplying by A number
1101                ConvertEPNToEnergy();
1102                // EpnEnergyH will be replace by UDefEnergyH.
1103                //      UDefEnergyH.DumpValues();
1104        }
1105
1106        //  G4cout << "Creating IPDFEnergy if not already done so" << G4endl;
1107        if (IPDFEnergyExist == false) {
1108                // IPDF has not been created, so create it
1109                G4double bins[1024], vals[1024], sum;
1110                G4int ii;
1111                G4int maxbin = G4int(UDefEnergyH.GetVectorLength());
1112                bins[0] = UDefEnergyH.GetLowEdgeEnergy(size_t(0));
1113                vals[0] = UDefEnergyH(size_t(0));
1114                sum = vals[0];
1115                for (ii = 1; ii < maxbin; ii++) {
1116                        bins[ii] = UDefEnergyH.GetLowEdgeEnergy(size_t(ii));
1117                        vals[ii] = UDefEnergyH(size_t(ii)) + vals[ii - 1];
1118                        sum = sum + UDefEnergyH(size_t(ii));
1119                }
1120
1121                for (ii = 0; ii < maxbin; ii++) {
1122                        vals[ii] = vals[ii] / sum;
1123                        IPDFEnergyH.InsertValues(bins[ii], vals[ii]);
1124                }
1125                // Make IPDFEpnExist = true
1126                IPDFEnergyExist = true;
1127        }
1128        //  IPDFEnergyH.DumpValues();
1129        // IPDF has been create so carry on
1130        G4double rndm = eneRndm->GenRandEnergy();
1131        particle_energy = IPDFEnergyH.GetEnergy(rndm);
1132
1133        if (verbosityLevel >= 1)
1134                G4cout << "Energy is " << particle_energy << G4endl;
1135}
1136
1137void G4SPSEneDistribution::ConvertEPNToEnergy() {
1138        // Use this before particle generation to convert  the
1139        // currently stored histogram from energy/nucleon
1140        // to energy.
1141        //  G4cout << "In ConvertEpntoEnergy " << G4endl;
1142        if (particle_definition == NULL)
1143                G4cout << "Error: particle not defined" << G4endl;
1144        else {
1145                // Need to multiply histogram by the number of nucleons.
1146                // Baryon Number looks to hold the no. of nucleons.
1147                G4int Bary = particle_definition->GetBaryonNumber();
1148                //      G4cout << "Baryon No. " << Bary << G4endl;
1149                // Change values in histogram, Read it out, delete it, re-create it
1150                G4int count, maxcount;
1151                maxcount = G4int(EpnEnergyH.GetVectorLength());
1152                //      G4cout << maxcount << G4endl;
1153                G4double ebins[1024], evals[1024];
1154                if (maxcount > 1024) {
1155                        G4cout << "Histogram contains more than 1024 bins!" << G4endl;
1156                        G4cout << "Those above 1024 will be ignored" << G4endl;
1157                        maxcount = 1024;
1158                }
1159                if (maxcount < 1) {
1160                        G4cout << "Histogram contains less than 1 bin!" << G4endl;
1161                        G4cout << "Redefine the histogram" << G4endl;
1162                        return;
1163                }
1164                for (count = 0; count < maxcount; count++) {
1165                        // Read out
1166                        ebins[count] = EpnEnergyH.GetLowEdgeEnergy(size_t(count));
1167                        evals[count] = EpnEnergyH(size_t(count));
1168                }
1169
1170                // Multiply the channels by the nucleon number to give energies
1171                for (count = 0; count < maxcount; count++) {
1172                        ebins[count] = ebins[count] * Bary;
1173                }
1174
1175                // Set Emin and Emax
1176                Emin = ebins[0];
1177                if (maxcount > 1)
1178                        Emax = ebins[maxcount - 1];
1179                else
1180                        Emax = ebins[0];
1181                // Put energy bins into new histogram - UDefEnergyH.
1182                for (count = 0; count < maxcount; count++) {
1183                        UDefEnergyH.InsertValues(ebins[count], evals[count]);
1184                }
1185                Epnflag = false; //so that you dont repeat this method.
1186        }
1187}
1188
1189//
1190void G4SPSEneDistribution::ReSetHist(G4String atype) {
1191        if (atype == "energy") {
1192                UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
1193                IPDFEnergyExist = false;
1194                Emin = 0.;
1195                Emax = 1e30;
1196        } else if (atype == "arb") {
1197                ArbEnergyH = IPDFArbEnergyH = ZeroPhysVector;
1198                IPDFArbExist = false;
1199        } else if (atype == "epn") {
1200                UDefEnergyH = IPDFEnergyH = ZeroPhysVector;
1201                IPDFEnergyExist = false;
1202                EpnEnergyH = ZeroPhysVector;
1203        } else {
1204                G4cout << "Error, histtype not accepted " << G4endl;
1205        }
1206}
1207
1208G4double G4SPSEneDistribution::GenerateOne(G4ParticleDefinition* a) {
1209        particle_definition = a;
1210        particle_energy = -1.;
1211
1212        while ((EnergyDisType == "Arb") ? (particle_energy < ArbEmin
1213                        || particle_energy > ArbEmax) : (particle_energy < Emin
1214                        || particle_energy > Emax)) {
1215                if (Biased) {
1216                        GenerateBiasPowEnergies();
1217                } else {
1218                        if (EnergyDisType == "Mono")
1219                                GenerateMonoEnergetic();
1220                        else if (EnergyDisType == "Lin")
1221                                GenerateLinearEnergies();
1222                        else if (EnergyDisType == "Pow")
1223                                GeneratePowEnergies();
1224                        else if (EnergyDisType == "Exp")
1225                                GenerateExpEnergies();
1226                        else if (EnergyDisType == "Gauss")
1227                                GenerateGaussEnergies();
1228                        else if (EnergyDisType == "Brem")
1229                                GenerateBremEnergies();
1230                        else if (EnergyDisType == "Bbody")
1231                                GenerateBbodyEnergies();
1232                        else if (EnergyDisType == "Cdg")
1233                                GenerateCdgEnergies();
1234                        else if (EnergyDisType == "User")
1235                                GenUserHistEnergies();
1236                        else if (EnergyDisType == "Arb")
1237                                GenArbPointEnergies();
1238                        else if (EnergyDisType == "Epn")
1239                                GenEpnHistEnergies();
1240                        else
1241                                G4cout << "Error: EnergyDisType has unusual value" << G4endl;
1242                }
1243        }
1244        return particle_energy;
1245}
1246
1247G4double G4SPSEneDistribution::GetProbability(G4double ene) {
1248        G4double prob = 1.;
1249
1250        if (EnergyDisType == "Lin") {
1251          if (prob_norm == 1.) {
1252            prob_norm = 0.5*grad*Emax*Emax + cept*Emax - 0.5*grad*Emin*Emin - cept*Emin;
1253          }
1254          prob = cept + grad * ene;
1255          prob /= prob_norm;
1256        }
1257        else if (EnergyDisType == "Pow") {
1258          if (prob_norm == 1.) {
1259            if (alpha != -1.) {
1260              G4double emina = std::pow(Emin, alpha + 1);
1261              G4double emaxa = std::pow(Emax, alpha + 1);
1262              prob_norm = 1./(1.+alpha) * (emaxa - emina);
1263            } else {
1264              prob_norm = std::log(Emax) - std::log(Emin) ;
1265            }
1266          }
1267          prob = std::pow(ene, alpha)/prob_norm;
1268        }
1269        else if (EnergyDisType == "Exp"){
1270          if (prob_norm == 1.) {
1271            prob_norm = -Ezero*(std::exp(-Emax/Ezero) - std::exp(Emin/Ezero));
1272          } 
1273          prob = std::exp(-ene / Ezero);
1274          prob /= prob_norm;
1275        }
1276        else if (EnergyDisType == "Arb") {
1277          prob = ArbEnergyH.Value(ene);
1278          //  prob = ArbEInt->CubicSplineInterpolation(ene);
1279          //G4double deltaY;
1280          //prob = ArbEInt->PolynomInterpolation(ene, deltaY);
1281          if (prob <= 0.) {
1282            //G4cout << " Warning:G4SPSEneDistribution::GetProbability: prob<= 0. "<<prob <<" "<<ene << " " <<deltaY<< G4endl;
1283            G4cout << " Warning:G4SPSEneDistribution::GetProbability: prob<= 0. "<<prob <<" "<<ene << G4endl;
1284            prob = 1e-30;
1285          }
1286          // already normalised
1287        }
1288        else
1289                G4cout << "Error: EnergyDisType not supported" << G4endl;
1290       
1291        return prob;
1292}
Note: See TracBrowser for help on using the repository browser.