Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/event/src/G4SPSEneDistribution.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 ]

Diff markup

Differences between /event/src/G4SPSEneDistribution.cc (Version 11.3.0) and /event/src/G4SPSEneDistribution.cc (Version 8.0)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                    <<   3 // * DISCLAIMER                                                       *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th <<   5 // * The following disclaimer summarizes all the specific disclaimers *
  6 // * the Geant4 Collaboration.  It is provided <<   6 // * of contributors to this software. The specific disclaimers,which *
  7 // * conditions of the Geant4 Software License <<   7 // * govern, are listed with their locations in:                      *
  8 // * LICENSE and available at  http://cern.ch/ <<   8 // *   http://cern.ch/geant4/license                                  *
  9 // * include a list of copyright holders.      << 
 10 // *                                                9 // *                                                                  *
 11 // * Neither the authors of this software syst     10 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     11 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     12 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     13 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  <<  14 // * use.                                                             *
 16 // * for the full disclaimer and the limitatio << 
 17 // *                                               15 // *                                                                  *
 18 // * This  code  implementation is the result  <<  16 // * This  code  implementation is the  intellectual property  of the *
 19 // * technical work of the GEANT4 collaboratio <<  17 // * GEANT4 collaboration.                                            *
 20 // * By using,  copying,  modifying or  distri <<  18 // * By copying,  distributing  or modifying the Program (or any work *
 21 // * any work based  on the software)  you  ag <<  19 // * based  on  the Program)  you indicate  your  acceptance of  this *
 22 // * use  in  resulting  scientific  publicati <<  20 // * statement, and all its terms.                                    *
 23 // * acceptance of all terms of the Geant4 Sof << 
 24 // *******************************************     21 // ********************************************************************
 25 //                                                 22 //
 26 // G4SPSEneDistribution class implementation   <<  23 ///////////////////////////////////////////////////////////////////////////////
                                                   >>  24 //
                                                   >>  25 // MODULE:        G4SPSEneDistribution.cc
                                                   >>  26 //
                                                   >>  27 // Version:      1.0
                                                   >>  28 // Date:         5/02/04
                                                   >>  29 // Author:       Fan Lei 
                                                   >>  30 // Organisation: QinetiQ ltd.
                                                   >>  31 // Customer:     ESA/ESTEC
                                                   >>  32 //
                                                   >>  33 ///////////////////////////////////////////////////////////////////////////////
                                                   >>  34 //
                                                   >>  35 // CHANGE HISTORY
                                                   >>  36 // --------------
                                                   >>  37 //
                                                   >>  38 //
                                                   >>  39 // Version 1.0, 05/02/2004, Fan Lei, Created.
                                                   >>  40 //    Based on the G4GeneralParticleSource class in Geant4 v6.0
                                                   >>  41 //
                                                   >>  42 ///////////////////////////////////////////////////////////////////////////////
 27 //                                                 43 //
 28 // Author: Fan Lei, QinetiQ ltd - 05/02/2004   << 
 29 // Customer: ESA/ESTEC                         << 
 30 // Revisions: Andrew Green, Andrea Dotti       << 
 31 // ------------------------------------------- << 
 32 #include "G4SPSEneDistribution.hh"             << 
 33                                                << 
 34 #include "G4Exp.hh"                            << 
 35 #include "G4SystemOfUnits.hh"                  << 
 36 #include "G4UnitsTable.hh"                     << 
 37 #include "Randomize.hh"                            44 #include "Randomize.hh"
 38 #include "G4AutoLock.hh"                       <<  45 //#include <cmath>
 39 #include "G4Threading.hh"                      <<  46 
                                                   >>  47 #include "G4SPSEneDistribution.hh"
 40                                                    48 
 41 G4SPSEneDistribution::G4SPSEneDistribution()       49 G4SPSEneDistribution::G4SPSEneDistribution()
 42 {                                                  50 {
 43   G4MUTEXINIT(mutex);                          <<  51   //
 44                                                << 
 45   // Initialise all variables                      52   // Initialise all variables
                                                   >>  53   particle_energy = 1.0*MeV;
 46                                                    54 
 47   particle_energy = 1.0 * MeV;                 << 
 48   EnergyDisType = "Mono";                          55   EnergyDisType = "Mono";
 49   weight=1.;                                   <<  56   MonoEnergy = 1*MeV;
 50   MonoEnergy = 1 * MeV;                        << 
 51   Emin = 0.;                                       57   Emin = 0.;
 52   Emax = 1.e30;                                    58   Emax = 1.e30;
 53   alpha = 0.;                                      59   alpha = 0.;
 54   biasalpha = 0.;                              << 
 55   prob_norm = 1.0;                             << 
 56   Ezero = 0.;                                      60   Ezero = 0.;
 57   SE = 0.;                                         61   SE = 0.;
 58   Temp = 0.;                                       62   Temp = 0.;
 59   grad = 0.;                                       63   grad = 0.;
 60   cept = 0.;                                       64   cept = 0.;
                                                   >>  65   EnergySpec = true; // true - energy spectra, false - momentum spectra
                                                   >>  66   DiffSpec = true;  // true - differential spec, false integral spec
 61   IntType = "NULL"; // Interpolation type          67   IntType = "NULL"; // Interpolation type
                                                   >>  68   IPDFEnergyExist = false;
                                                   >>  69   IPDFArbExist = false;
 62                                                    70 
 63   ArbEmin = 0.;                                    71   ArbEmin = 0.;
 64   ArbEmax = 1.e30;                                 72   ArbEmax = 1.e30;
 65                                                    73 
 66   verbosityLevel = 0;                          <<  74   verbosityLevel = 0 ;
 67                                                <<  75 
 68   threadLocal_t& data = threadLocalData.Get(); << 
 69   data.Emax = Emax;                            << 
 70   data.Emin = Emin;                            << 
 71   data.alpha =alpha;                           << 
 72   data.cept = cept;                            << 
 73   data.Ezero = Ezero;                          << 
 74   data.grad = grad;                            << 
 75   data.particle_energy = 0.;                   << 
 76   data.particle_definition = nullptr;          << 
 77   data.weight = weight;                        << 
 78 }                                                  76 }
 79                                                    77 
 80 G4SPSEneDistribution::~G4SPSEneDistribution()      78 G4SPSEneDistribution::~G4SPSEneDistribution()
 81 {                                              <<  79 {}
 82   G4MUTEXDESTROY(mutex);                       << 
 83   if(Arb_grad_cept_flag)                       << 
 84   {                                            << 
 85     delete [] Arb_grad;                        << 
 86     delete [] Arb_cept;                        << 
 87   }                                            << 
 88                                                << 
 89   if(Arb_alpha_Const_flag)                     << 
 90   {                                            << 
 91     delete [] Arb_alpha;                       << 
 92     delete [] Arb_Const;                       << 
 93   }                                            << 
 94                                                << 
 95   if(Arb_ezero_flag)                           << 
 96   {                                            << 
 97     delete [] Arb_ezero;                       << 
 98   }                                            << 
 99   delete Bbody_x;                              << 
100   delete BBHist;                               << 
101   delete CP_x;                                 << 
102   delete CPHist;                               << 
103   for (auto & it : SplineInt)                  << 
104   {                                            << 
105     delete it;                                 << 
106     it = nullptr;                              << 
107   }                                            << 
108   SplineInt.clear();                           << 
109 }                                              << 
110                                                    80 
111 void G4SPSEneDistribution::SetEnergyDisType(co <<  81 void G4SPSEneDistribution::SetEnergyDisType(G4String DisType)
112 {                                                  82 {
113   G4AutoLock l(&mutex);                        << 
114   EnergyDisType = DisType;                         83   EnergyDisType = DisType;
115   if (EnergyDisType == "User")                 <<  84   if (EnergyDisType == "User"){
116   {                                            <<  85     UDefEnergyH = IPDFEnergyH = ZeroPhysVector ;
117     UDefEnergyH = IPDFEnergyH = ZeroPhysVector <<  86     IPDFEnergyExist = false ;
118     IPDFEnergyExist = false;                   <<  87   } else if ( EnergyDisType == "Arb"){
119   }                                            <<  88     ArbEnergyH =IPDFArbEnergyH = ZeroPhysVector ;
120   else if (EnergyDisType == "Arb")             << 
121   {                                            << 
122     ArbEnergyH = IPDFArbEnergyH = ZeroPhysVect << 
123     IPDFArbExist = false;                          89     IPDFArbExist = false;
124   }                                            <<  90   } else if (EnergyDisType == "Epn"){
125   else if (EnergyDisType == "Epn")             <<  91     UDefEnergyH = IPDFEnergyH = ZeroPhysVector ;
126   {                                            <<  92     IPDFEnergyExist = false ;
127     UDefEnergyH = IPDFEnergyH = ZeroPhysVector <<  93     EpnEnergyH = ZeroPhysVector ;
128     IPDFEnergyExist = false;                   << 
129     EpnEnergyH = ZeroPhysVector;               << 
130   }                                                94   }
131 }                                                  95 }
132                                                    96 
133 const G4String& G4SPSEneDistribution::GetEnerg << 
134 {                                              << 
135   G4AutoLock l(&mutex);                        << 
136   return EnergyDisType;                        << 
137 }                                              << 
138                                                << 
139 void G4SPSEneDistribution::SetEmin(G4double em     97 void G4SPSEneDistribution::SetEmin(G4double emi)
140 {                                                  98 {
141   G4AutoLock l(&mutex);                        << 
142   Emin = emi;                                      99   Emin = emi;
143   threadLocalData.Get().Emin = Emin;           << 
144 }                                              << 
145                                                << 
146 G4double G4SPSEneDistribution::GetEmin() const << 
147 {                                              << 
148   return threadLocalData.Get().Emin;           << 
149 }                                              << 
150                                                << 
151 G4double G4SPSEneDistribution::GetArbEmin()    << 
152 {                                              << 
153   G4AutoLock l(&mutex);                        << 
154   return ArbEmin;                              << 
155 }                                              << 
156                                                << 
157 G4double G4SPSEneDistribution::GetArbEmax()    << 
158 {                                              << 
159   G4AutoLock l(&mutex);                        << 
160   return ArbEmax;                              << 
161 }                                                 100 }
162                                                   101 
163 void G4SPSEneDistribution::SetEmax(G4double em    102 void G4SPSEneDistribution::SetEmax(G4double ema)
164 {                                                 103 {
165   G4AutoLock l(&mutex);                        << 
166   Emax = ema;                                     104   Emax = ema;
167   threadLocalData.Get().Emax = Emax;           << 
168 }                                              << 
169                                                << 
170 G4double G4SPSEneDistribution::GetEmax() const << 
171 {                                              << 
172   return threadLocalData.Get().Emax;           << 
173 }                                                 105 }
174                                                   106 
175 void G4SPSEneDistribution::SetMonoEnergy(G4dou    107 void G4SPSEneDistribution::SetMonoEnergy(G4double menergy)
176 {                                                 108 {
177   G4AutoLock l(&mutex);                        << 
178   MonoEnergy = menergy;                           109   MonoEnergy = menergy;
179 }                                                 110 }
180                                                   111 
181 void G4SPSEneDistribution::SetBeamSigmaInE(G4d    112 void G4SPSEneDistribution::SetBeamSigmaInE(G4double e)
182 {                                                 113 {
183   G4AutoLock l(&mutex);                        << 
184   SE = e;                                         114   SE = e;
185 }                                                 115 }
186 void G4SPSEneDistribution::SetAlpha(G4double a    116 void G4SPSEneDistribution::SetAlpha(G4double alp)
187 {                                                 117 {
188   G4AutoLock l(&mutex);                        << 
189   alpha = alp;                                    118   alpha = alp;
190   threadLocalData.Get().alpha = alpha;         << 
191 }                                              << 
192                                                << 
193 void G4SPSEneDistribution::SetBiasAlpha(G4doub << 
194 {                                              << 
195   G4AutoLock l(&mutex);                        << 
196   biasalpha = alp;                             << 
197   Biased = true;                               << 
198 }                                                 119 }
199                                                   120 
200 void G4SPSEneDistribution::SetTemp(G4double te    121 void G4SPSEneDistribution::SetTemp(G4double tem)
201 {                                                 122 {
202   G4AutoLock l(&mutex);                        << 
203   Temp = tem;                                     123   Temp = tem;
204 }                                                 124 }
205                                                   125 
206 void G4SPSEneDistribution::SetEzero(G4double e    126 void G4SPSEneDistribution::SetEzero(G4double eze)
207 {                                                 127 {
208   G4AutoLock l(&mutex);                        << 
209   Ezero = eze;                                    128   Ezero = eze;
210   threadLocalData.Get().Ezero = Ezero;         << 
211 }                                                 129 }
212                                                   130 
213 void G4SPSEneDistribution::SetGradient(G4doubl    131 void G4SPSEneDistribution::SetGradient(G4double gr)
214 {                                                 132 {
215   G4AutoLock l(&mutex);                        << 
216   grad = gr;                                      133   grad = gr;
217   threadLocalData.Get().grad = grad;           << 
218 }                                                 134 }
219                                                   135 
220 void G4SPSEneDistribution::SetInterCept(G4doub    136 void G4SPSEneDistribution::SetInterCept(G4double c)
221 {                                                 137 {
222   G4AutoLock l(&mutex);                        << 
223   cept = c;                                       138   cept = c;
224   threadLocalData.Get().cept = cept;           << 
225 }                                              << 
226                                                << 
227 const G4String& G4SPSEneDistribution::GetIntTy << 
228 {                                              << 
229   G4AutoLock l(&mutex);                        << 
230   return IntType;                              << 
231 }                                              << 
232                                                << 
233 void G4SPSEneDistribution::SetBiasRndm(G4SPSRa << 
234 {                                              << 
235   G4AutoLock l(&mutex);                        << 
236   eneRndm = a;                                 << 
237 }                                              << 
238                                                << 
239 void G4SPSEneDistribution::SetVerbosity(G4int  << 
240 {                                              << 
241   G4AutoLock l(&mutex);                        << 
242   verbosityLevel = a;                          << 
243 }                                              << 
244                                                << 
245 G4double G4SPSEneDistribution::GetWeight() con << 
246 {                                              << 
247   return threadLocalData.Get().weight;         << 
248 }                                              << 
249                                                << 
250 G4double G4SPSEneDistribution::GetMonoEnergy() << 
251 {                                              << 
252   G4AutoLock l(&mutex);                        << 
253   return MonoEnergy;                           << 
254 }                                              << 
255                                                << 
256 G4double G4SPSEneDistribution::GetSE()         << 
257 {                                              << 
258   G4AutoLock l(&mutex);                        << 
259   return SE;                                   << 
260 }                                              << 
261                                                << 
262 G4double G4SPSEneDistribution::Getalpha() cons << 
263 {                                              << 
264   return threadLocalData.Get().alpha;          << 
265 }                                              << 
266                                                << 
267 G4double G4SPSEneDistribution::GetEzero() cons << 
268 {                                              << 
269   return threadLocalData.Get().Ezero;          << 
270 }                                              << 
271                                                << 
272 G4double G4SPSEneDistribution::GetTemp()       << 
273 {                                              << 
274   G4AutoLock l(&mutex);                        << 
275   return Temp;                                 << 
276 }                                              << 
277                                                << 
278 G4double G4SPSEneDistribution::Getgrad() const << 
279 {                                              << 
280   return threadLocalData.Get().grad;           << 
281 }                                                 139 }
282                                                   140 
283 G4double G4SPSEneDistribution::Getcept() const << 141 void G4SPSEneDistribution::UserEnergyHisto(G4ThreeVector input)
284 {                                                 142 {
285   return threadLocalData.Get().cept;           << 143   G4double ehi, val;
286 }                                              << 144   ehi = input.x();
287                                                << 145   val = input.y();
288 G4PhysicsFreeVector G4SPSEneDistribution::GetU << 146   if(verbosityLevel > 1) {
289 {                                              << 
290   G4AutoLock l(&mutex);                        << 
291   return UDefEnergyH;                          << 
292 }                                              << 
293                                                << 
294 G4PhysicsFreeVector G4SPSEneDistribution::GetA << 
295 {                                              << 
296   G4AutoLock l(&mutex);                        << 
297   return ArbEnergyH;                           << 
298 }                                              << 
299                                                << 
300 void G4SPSEneDistribution::UserEnergyHisto(con << 
301 {                                              << 
302   G4AutoLock l(&mutex);                        << 
303   G4double ehi = input.x(),                    << 
304            val = input.y();                    << 
305   if (verbosityLevel > 1)                      << 
306   {                                            << 
307     G4cout << "In UserEnergyHisto" << G4endl;     147     G4cout << "In UserEnergyHisto" << G4endl;
308     G4cout << " " << ehi << " " << val << G4en    148     G4cout << " " << ehi << " " << val << G4endl;
309   }                                               149   }
310   UDefEnergyH.InsertValues(ehi, val);             150   UDefEnergyH.InsertValues(ehi, val);
311   Emax = ehi;                                     151   Emax = ehi;
312   threadLocalData.Get().Emax = Emax;           << 
313 }                                                 152 }
314                                                   153 
315 void G4SPSEneDistribution::ArbEnergyHisto(cons << 154 void G4SPSEneDistribution::ArbEnergyHisto(G4ThreeVector input)
316 {                                                 155 {
317   G4AutoLock l(&mutex);                        << 156   G4double ehi, val;
318   G4double ehi = input.x(),                    << 157   ehi = input.x();
319            val = input.y();                    << 158   val = input.y();
320   if (verbosityLevel > 1)                      << 159   if(verbosityLevel >1 ) {
321   {                                            << 
322     G4cout << "In ArbEnergyHisto" << G4endl;      160     G4cout << "In ArbEnergyHisto" << G4endl;
323     G4cout << " " << ehi << " " << val << G4en    161     G4cout << " " << ehi << " " << val << G4endl;
324   }                                               162   }
325   ArbEnergyH.InsertValues(ehi, val);              163   ArbEnergyH.InsertValues(ehi, val);
326 }                                                 164 }
327                                                   165 
328 void G4SPSEneDistribution::ArbEnergyHistoFile( << 166 void G4SPSEneDistribution::EpnEnergyHisto(G4ThreeVector input)
329 {                                                 167 {
330   G4AutoLock l(&mutex);                        << 
331   std::ifstream infile(filename, std::ios::in) << 
332   if (!infile)                                 << 
333   {                                            << 
334     G4Exception("G4SPSEneDistribution::ArbEner << 
335                 FatalException, "Unable to ope << 
336   }                                            << 
337   G4double ehi, val;                              168   G4double ehi, val;
338   while (infile >> ehi >> val)                 << 169   ehi = input.x();
339   {                                            << 170   val = input.y();
340     ArbEnergyH.InsertValues(ehi, val);         << 171   if(verbosityLevel > 1) {
341   }                                            << 
342 }                                              << 
343                                                << 
344 void G4SPSEneDistribution::EpnEnergyHisto(cons << 
345 {                                              << 
346   G4AutoLock l(&mutex);                        << 
347   G4double ehi = input.x(),                    << 
348            val = input.y();                    << 
349   if (verbosityLevel > 1)                      << 
350   {                                            << 
351     G4cout << "In EpnEnergyHisto" << G4endl;      172     G4cout << "In EpnEnergyHisto" << G4endl;
352     G4cout << " " << ehi << " " << val << G4en    173     G4cout << " " << ehi << " " << val << G4endl;
353   }                                               174   }
354   EpnEnergyH.InsertValues(ehi, val);              175   EpnEnergyH.InsertValues(ehi, val);
355   Emax = ehi;                                     176   Emax = ehi;
356   threadLocalData.Get().Emax = Emax;           << 
357   Epnflag = true;                                 177   Epnflag = true;
358 }                                                 178 }
359                                                   179 
360 void G4SPSEneDistribution::Calculate()            180 void G4SPSEneDistribution::Calculate()
361 {                                                 181 {
362   G4AutoLock l(&mutex);                        << 182   if(EnergyDisType == "Cdg")
363   if (EnergyDisType == "Cdg")                  << 
364   {                                            << 
365     CalculateCdgSpectrum();                       183     CalculateCdgSpectrum();
366   }                                            << 184   else if(EnergyDisType == "Bbody")
367   else if (EnergyDisType == "Bbody")           << 
368   {                                            << 
369     if(!BBhistInit)                            << 
370     {                                          << 
371       BBInitHists();                           << 
372     }                                          << 
373     CalculateBbodySpectrum();                     185     CalculateBbodySpectrum();
374   }                                            << 
375   else if (EnergyDisType == "CPow")            << 
376   {                                            << 
377     if(!CPhistInit)                            << 
378     {                                          << 
379       CPInitHists();                           << 
380     }                                          << 
381     CalculateCPowSpectrum();                   << 
382   }                                            << 
383 }                                                 186 }
384                                                   187 
385 void G4SPSEneDistribution::BBInitHists()  // M << 188 void G4SPSEneDistribution::CalculateCdgSpectrum()
386 {                                                 189 {
387   BBHist = new std::vector<G4double>(10001, 0. << 190   // This uses the spectrum from The INTEGRAL Mass Model (TIMM)
388   Bbody_x = new std::vector<G4double>(10001, 0 << 
389   BBhistInit = true;                           << 
390 }                                              << 
391                                                << 
392 void G4SPSEneDistribution::CPInitHists()  // M << 
393 {                                              << 
394   CPHist = new std::vector<G4double>(10001, 0. << 
395   CP_x = new std::vector<G4double>(10001, 0.0) << 
396   CPhistInit = true;                           << 
397 }                                              << 
398                                                << 
399 void G4SPSEneDistribution::CalculateCdgSpectru << 
400 {                                              << 
401   // This uses the spectrum from the INTEGRAL  << 
402   // to generate a Cosmic Diffuse X/gamma ray     191   // to generate a Cosmic Diffuse X/gamma ray spectrum.
403                                                << 192   G4double pfact[2] = {8.5, 112};
404   G4double pfact[2] = { 8.5, 112 };            << 193   G4double spind[2] = {1.4, 2.3};
405   G4double spind[2] = { 1.4, 2.3 };            << 194   G4double ene_line[3] = {1.*keV, 18.*keV, 1E6*keV};
406   G4double ene_line[3] = { 1. * keV, 18. * keV << 
407   G4int n_par;                                    195   G4int n_par;
408                                                   196 
409   ene_line[0] = threadLocalData.Get().Emin;    << 197   ene_line[0] = Emin;
410   if (threadLocalData.Get().Emin < 18 * keV)   << 198   if(Emin < 18*keV)
411   {                                            << 
412     n_par = 2;                                 << 
413     ene_line[2] = threadLocalData.Get().Emax;  << 
414     if (threadLocalData.Get().Emax < 18 * keV) << 
415     {                                             199     {
416       n_par = 1;                               << 200       n_par = 2;
417       ene_line[1] = threadLocalData.Get().Emax << 201       ene_line[2] = Emax;
                                                   >> 202       if(Emax < 18*keV)
                                                   >> 203   {
                                                   >> 204     n_par = 1;
                                                   >> 205     ene_line[1] = Emax;
                                                   >> 206   }
418     }                                             207     }
419   }                                            << 
420   else                                            208   else
421   {                                            << 209     {
422     n_par = 1;                                 << 210       n_par = 1;
423     pfact[0] = 112.;                           << 211       pfact[0] = 112.;
424     spind[0] = 2.3;                            << 212       spind[0] = 2.3;
425     ene_line[1] = threadLocalData.Get().Emax;  << 213       ene_line[1] = Emax;
426   }                                            << 214     }
427                                                << 215   
428   // Create a cumulative histogram             << 216   // Create a cumulative histogram.
429   //                                           << 
430   CDGhist[0] = 0.;                                217   CDGhist[0] = 0.;
431   G4double omalpha;                               218   G4double omalpha;
432   G4int i = 0;                                    219   G4int i = 0;
433   while (i < n_par)                            << 
434   {                                            << 
435     omalpha = 1. - spind[i];                   << 
436     CDGhist[i + 1] = CDGhist[i] + (pfact[i] /  << 
437                                 * (std::pow(en << 
438                                 - std::pow(ene << 
439     ++i;                                       << 
440   }                                            << 
441                                                   220 
442   // Normalise histo and divide by 1000 to mak << 221   while(i < n_par)
443   //                                           << 222     {
                                                   >> 223       omalpha = 1. - spind[i];
                                                   >> 224       CDGhist[i+1] = CDGhist[i] + (pfact[i]/omalpha)*
                                                   >> 225   (std::pow(ene_line[i+1],omalpha)-std::pow(ene_line[i],omalpha));
                                                   >> 226       i++;
                                                   >> 227     }
                                                   >> 228   
                                                   >> 229   // Normalise histo and divide by 1000 to make MeV.
444   i = 0;                                          230   i = 0;
445   while (i < n_par)                            << 231   while(i < n_par)
446   {                                            << 232     {
447     CDGhist[i + 1] = CDGhist[i + 1] / CDGhist[ << 233       CDGhist[i+1] = CDGhist[i+1]/CDGhist[n_par];
448     ++i;                                       << 234       //      G4cout << CDGhist[i] << CDGhist[n_par] << G4endl;
449   }                                            << 235       i++;
                                                   >> 236     }
450 }                                                 237 }
451                                                   238 
452 void G4SPSEneDistribution::CalculateBbodySpect << 239 void G4SPSEneDistribution::CalculateBbodySpectrum()
453 {                                                 240 {
454   // Create bbody spectrum                     << 241   // create bbody spectrum
455   // Proved very hard to integrate indefinitel << 242   // Proved very hard to integrate indefinitely, so different
456   // User inputs emin, emax and T. These are u << 243   // method. User inputs emin, emax and T. These are used to
457   // bin histogram.                            << 244   // create a 10,000 bin histogram.
458   // Use photon density spectrum = 2 nu**2/c**    245   // Use photon density spectrum = 2 nu**2/c**2 * (std::exp(h nu/kT)-1)
459   // = 2 E**2/h**2c**2 times the exponential      246   // = 2 E**2/h**2c**2 times the exponential
460                                                << 247   G4double erange = Emax - Emin;
461   G4double erange = threadLocalData.Get().Emax << 248   G4double steps = erange/10000.;
462   G4double steps = erange / 10000.;            << 249   G4double Bbody_y[10000];
463                                                << 250   G4double k = 8.6181e-11; //Boltzmann const in MeV/K
464   const G4double k = 8.6181e-11; //Boltzmann c << 251   G4double h = 4.1362e-21; // Plancks const in MeV s
465   const G4double h = 4.1362e-21; // Plancks co << 252   G4double c = 3e8; // Speed of light
466   const G4double c = 3e8; // Speed of light    << 253   G4double h2 = h*h;
467   const G4double h2 = h * h;                   << 254   G4double c2 = c*c;
468   const G4double c2 = c * c;                   << 
469   G4int count = 0;                             << 
470   G4double sum = 0.;                           << 
471   BBHist->at(0) = 0.;                          << 
472                                                << 
473   while (count < 10000)                        << 
474   {                                            << 
475     Bbody_x->at(count) = threadLocalData.Get() << 
476     G4double Bbody_y = (2. * std::pow(Bbody_x- << 
477                      / (h2*c2*(std::exp(Bbody_ << 
478     sum = sum + Bbody_y;                       << 
479     BBHist->at(count + 1) = BBHist->at(count)  << 
480     ++count;                                   << 
481   }                                            << 
482                                                << 
483   Bbody_x->at(10000) = threadLocalData.Get().E << 
484                                                << 
485   // Normalise cumulative histo                << 
486   //                                           << 
487   count = 0;                                   << 
488   while (count < 10001)                        << 
489   {                                            << 
490     BBHist->at(count) = BBHist->at(count) / su << 
491     ++count;                                   << 
492   }                                            << 
493 }                                              << 
494                                                << 
495 void G4SPSEneDistribution::CalculateCPowSpectr << 
496 {                                              << 
497   // Create cutoff power-law spectrum, x^a exp << 
498   // The integral of this function is an incom << 
499   // is only available in the Boost library.   << 
500   //                                           << 
501   // User inputs are emin, emax and alpha and  << 
502   // create a 10,000 bin histogram.            << 
503                                                << 
504   G4double erange = threadLocalData.Get().Emax << 
505   G4double steps = erange / 10000.;            << 
506   alpha = threadLocalData.Get().alpha ;        << 
507   Ezero = threadLocalData.Get().Ezero ;        << 
508                                                << 
509   G4int count = 0;                                255   G4int count = 0;
510   G4double sum = 0.;                              256   G4double sum = 0.;
511   CPHist->at(0) = 0.;                          << 257   BBHist[0] = 0.;
512                                                << 258   while(count < 10000)
513   while (count < 10000)                        << 259     {
514   {                                            << 260       Bbody_x[count] = Emin + G4double(count*steps);
515     CP_x->at(count) = threadLocalData.Get().Em << 261       Bbody_y[count] = (2.*std::pow(Bbody_x[count],2.))/
516     G4double CP_y = std::pow(CP_x->at(count),  << 262         (h2*c2*(std::exp(Bbody_x[count]/(k*Temp)) - 1.));
517                   * std::exp(-CP_x->at(count)  << 263       sum = sum + Bbody_y[count];
518     sum = sum + CP_y;                          << 264       BBHist[count+1] = BBHist[count] + Bbody_y[count];
519     CPHist->at(count + 1) = CPHist->at(count)  << 265       count++;
520     ++count;                                   << 266     }
521   }                                            << 
522                                                << 
523   CP_x->at(10000) = threadLocalData.Get().Emax << 
524                                                   267 
525   // Normalise cumulative histo                << 268   Bbody_x[10000] = Emax;
526   //                                           << 269   // Normalise cumulative histo.
527   count = 0;                                      270   count = 0;
528   while (count < 10001)                        << 271   while(count<10001)
529   {                                            << 272     {
530     CPHist->at(count) = CPHist->at(count) / su << 273       BBHist[count] = BBHist[count]/sum;
531     ++count;                                   << 274       count++;
532   }                                            << 275     }
533 }                                                 276 }
534                                                   277 
535 void G4SPSEneDistribution::InputEnergySpectra(    278 void G4SPSEneDistribution::InputEnergySpectra(G4bool value)
536 {                                                 279 {
537   G4AutoLock l(&mutex);                        << 280   // Allows user to specifiy spectrum is momentum
538                                                << 
539   // Allows user to specify spectrum is moment << 
540   //                                           << 
541   EnergySpec = value; // false if momentum        281   EnergySpec = value; // false if momentum
542   if (verbosityLevel > 1)                      << 282   if(verbosityLevel > 1)
543   {                                            << 
544     G4cout << "EnergySpec has value " << Energ    283     G4cout << "EnergySpec has value " << EnergySpec << G4endl;
545   }                                            << 
546 }                                                 284 }
547                                                   285 
548 void G4SPSEneDistribution::InputDifferentialSp    286 void G4SPSEneDistribution::InputDifferentialSpectra(G4bool value)
549 {                                                 287 {
550   G4AutoLock l(&mutex);                        << 
551                                                << 
552   // Allows user to specify integral or differ    288   // Allows user to specify integral or differential spectra
553   //                                           << 
554   DiffSpec = value; // true = differential, fa    289   DiffSpec = value; // true = differential, false = integral
555   if (verbosityLevel > 1)                      << 290   if(verbosityLevel > 1)
556   {                                            << 
557     G4cout << "Diffspec has value " << DiffSpe    291     G4cout << "Diffspec has value " << DiffSpec << G4endl;
558   }                                            << 
559 }                                                 292 }
560                                                   293 
561 void G4SPSEneDistribution::ArbInterpolate(cons << 294 void G4SPSEneDistribution::ArbInterpolate(G4String IType)
562 {                                                 295 {
563   G4AutoLock l(&mutex);                        << 296   if(EnergyDisType != "Arb")
564                                                << 297     G4cout << "Error: this is for arbitrary distributions" << G4endl;
565   IntType = IType;                                298   IntType = IType;
566   ArbEmax = ArbEnergyH.GetMaxEnergy();         << 299   ArbEmax = Emax;
567   ArbEmin = ArbEnergyH.Energy(0);              << 300   ArbEmin = Emin;
568                                                   301 
569   // Now interpolate points                       302   // Now interpolate points
570                                                << 303   if(IntType == "Lin")
571   if (IntType == "Lin") LinearInterpolation(); << 304     LinearInterpolation();
572   if (IntType == "Log") LogInterpolation();    << 305   if(IntType == "Log")
573   if (IntType == "Exp") ExpInterpolation();    << 306     LogInterpolation();
574   if (IntType == "Spline") SplineInterpolation << 307   if(IntType == "Exp")
                                                   >> 308     ExpInterpolation();
                                                   >> 309   if(IntType == "Spline")
                                                   >> 310     SplineInterpolation();  
575 }                                                 311 }
576                                                   312 
577 void G4SPSEneDistribution::LinearInterpolation << 313 void G4SPSEneDistribution::LinearInterpolation()
578 {                                                 314 {
579   // Method to do linear interpolation on the     315   // Method to do linear interpolation on the Arb points
580   // Calculate equation of each line segment,     316   // Calculate equation of each line segment, max 1024.
581   // Calculate Area under each segment            317   // Calculate Area under each segment
582   // Create a cumulative array which is then n    318   // Create a cumulative array which is then normalised Arb_Cum_Area
583                                                   319 
584   G4double Area_seg[1024]; // Stores area unde    320   G4double Area_seg[1024]; // Stores area under each segment
585   G4double sum = 0., Arb_x[1024]={0.}, Arb_y[1 << 321   G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
586   std::size_t i, count;                        << 322   G4int i, count;
587   std::size_t maxi = ArbEnergyH.GetVectorLengt << 323   G4int maxi = ArbEnergyH.GetVectorLength();
588   for (i = 0; i < maxi; ++i)                   << 324   for(i=0;i<maxi;i++) {
589   {                                            << 325     Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
590     Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(i); << 326     Arb_y[i] = ArbEnergyH(size_t(i));
591     Arb_y[i] = ArbEnergyH(i);                  << 
592   }                                               327   }
593                                                << 
594   // Points are now in x,y arrays. If the spec    328   // Points are now in x,y arrays. If the spectrum is integral it has to be
595   // made differential and if momentum it has  << 329   // made differential and if momentum it has to be made energy.
596                                                << 330   if(DiffSpec == false) {
597   if (!DiffSpec)                               << 
598   {                                            << 
599     // Converts integral point-wise spectra to    331     // Converts integral point-wise spectra to Differential
600     //                                         << 332     for( count=0;count < maxi-1;count++) {
601     for (count = 0; count < maxi-1; ++count)   << 333       Arb_y[count] = (Arb_y[count] - Arb_y[count+1])/(Arb_x[count+1]-Arb_x[count]);
602     {                                          << 
603       Arb_y[count] = (Arb_y[count] - Arb_y[cou << 
604                    / (Arb_x[count + 1] - Arb_x << 
605     }                                             334     }
606     --maxi;                                    << 335     maxi--;
607   }                                               336   }
608                                                << 337   //
609   if (!EnergySpec)                             << 338   if(EnergySpec == false) {
610   {                                            << 
611     // change currently stored values (emin et    339     // change currently stored values (emin etc) which are actually momenta
612     // to energies                             << 340     // to energies.
613     //                                         << 341     if(particle_definition == NULL)
614     G4ParticleDefinition* pdef = threadLocalDa << 342       G4cout << "Error: particle not defined" << G4endl;
615     if (pdef == nullptr)                       << 343     else {
616     {                                          << 
617       G4Exception("G4SPSEneDistribution::Linea << 
618                   "Event0302", FatalException, << 
619                   "Error: particle not defined << 
620     }                                          << 
621     else                                       << 
622     {                                          << 
623       // Apply Energy**2 = p**2c**2 + m0**2c**    344       // Apply Energy**2 = p**2c**2 + m0**2c**4
624       // p should be entered as E/c i.e. witho    345       // p should be entered as E/c i.e. without the division by c
625       // being done - energy equivalent        << 346       // being done - energy equivalent.
626                                                << 347       G4double mass = particle_definition->GetPDGMass();
627       G4double mass = pdef->GetPDGMass();      << 348       // convert point to energy unit and its value to per energy unit
628                                                << 
629       // Convert point to energy unit and its  << 
630       //                                       << 
631       G4double total_energy;                      349       G4double total_energy;
632       for (count = 0; count < maxi; ++count)   << 350       for(count=0;count<maxi;count++) {
633       {                                        << 351   total_energy = std::sqrt((Arb_x[count]*Arb_x[count]) 
634         total_energy = std::sqrt((Arb_x[count] << 352           + (mass*mass)); // total energy
635                      + (mass * mass)); // tota << 353   
636         Arb_y[count] = Arb_y[count] * Arb_x[co << 354   Arb_y[count] = Arb_y[count] * Arb_x[count]/total_energy; 
637         Arb_x[count] = total_energy - mass; // << 355   Arb_x[count] = total_energy - mass ;  // kinetic energy
638       }                                           356       }
639     }                                             357     }
640   }                                               358   }
641                                                << 359   //
642   i = 1;                                       << 360   i=1;
643                                                << 
644   Arb_grad = new G4double [1024];              << 
645   Arb_cept = new G4double [1024];              << 
646   Arb_grad_cept_flag = true;                   << 
647                                                << 
648   Arb_grad[0] = 0.;                               361   Arb_grad[0] = 0.;
649   Arb_cept[0] = 0.;                               362   Arb_cept[0] = 0.;
650   Area_seg[0] = 0.;                               363   Area_seg[0] = 0.;
651   Arb_Cum_Area[0] = 0.;                           364   Arb_Cum_Area[0] = 0.;
652   while (i < maxi)                             << 365   while(i < maxi)
653   {                                            << 
654     // calculate gradient and intercept for ea << 
655     //                                         << 
656     Arb_grad[i] = (Arb_y[i] - Arb_y[i - 1]) /  << 
657     if (verbosityLevel == 2)                   << 
658     {                                          << 
659       G4cout << Arb_grad[i] << G4endl;         << 
660     }                                          << 
661     if (Arb_grad[i] > 0.)                      << 
662     {                                          << 
663       if (verbosityLevel == 2)                 << 
664       {                                        << 
665          G4cout << "Arb_grad is positive" << G << 
666       }                                        << 
667       Arb_cept[i] = Arb_y[i] - (Arb_grad[i] *  << 
668     }                                          << 
669     else if (Arb_grad[i] < 0.)                 << 
670     {                                             366     {
671       if (verbosityLevel == 2)                 << 367       // calc gradient and intercept for each segment
672       {                                        << 368       Arb_grad[i] = (Arb_y[i] - Arb_y[i-1]) / (Arb_x[i] - Arb_x[i-1]);
673         G4cout << "Arb_grad is negative" << G4 << 369       if(verbosityLevel == 2)
674       }                                        << 370   G4cout << Arb_grad[i] << G4endl;
675       Arb_cept[i] = Arb_y[i] + (-Arb_grad[i] * << 371       if(Arb_grad[i] > 0.)
                                                   >> 372   {
                                                   >> 373     if(verbosityLevel == 2)
                                                   >> 374       G4cout << "Arb_grad is positive" << G4endl;
                                                   >> 375     Arb_cept[i] = Arb_y[i] - (Arb_grad[i] * Arb_x[i]);
                                                   >> 376   }
                                                   >> 377       else if(Arb_grad[i] < 0.)
                                                   >> 378   {
                                                   >> 379     if(verbosityLevel == 2)
                                                   >> 380       G4cout << "Arb_grad is negative" << G4endl;
                                                   >> 381     Arb_cept[i] = Arb_y[i] + (-Arb_grad[i] * Arb_x[i]);
                                                   >> 382   }
                                                   >> 383       else
                                                   >> 384   {
                                                   >> 385     if(verbosityLevel == 2)
                                                   >> 386       G4cout << "Arb_grad is 0." << G4endl;
                                                   >> 387     Arb_cept[i] = Arb_y[i];
                                                   >> 388   }
                                                   >> 389 
                                                   >> 390       Area_seg[i] = ((Arb_grad[i]/2)*(Arb_x[i]*Arb_x[i] - Arb_x[i-1]*Arb_x[i-1]) + Arb_cept[i]*(Arb_x[i] - Arb_x[i-1]));
                                                   >> 391       Arb_Cum_Area[i] = Arb_Cum_Area[i-1] + Area_seg[i];
                                                   >> 392       sum = sum + Area_seg[i];
                                                   >> 393       if(verbosityLevel == 2)
                                                   >> 394   G4cout << Arb_x[i] << Arb_y[i] << Area_seg[i] << sum << Arb_grad[i] << G4endl;
                                                   >> 395       i++;
676     }                                             396     }
677     else                                       << 397 
                                                   >> 398   i=0;
                                                   >> 399   while(i < maxi)
678     {                                             400     {
679       if (verbosityLevel == 2)                 << 401       Arb_Cum_Area[i] = Arb_Cum_Area[i]/sum; // normalisation
680       {                                        << 402       IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
681         G4cout << "Arb_grad is 0." << G4endl;  << 403       i++;
682       }                                        << 
683       Arb_cept[i] = Arb_y[i];                  << 
684     }                                             404     }
685                                                   405 
686     Area_seg[i] = ((Arb_grad[i] / 2)           << 406   if(verbosityLevel >= 1)
687                 * (Arb_x[i] * Arb_x[i] - Arb_x << 
688                 + Arb_cept[i] * (Arb_x[i] - Ar << 
689     Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Ar << 
690     sum = sum + Area_seg[i];                   << 
691     if (verbosityLevel == 2)                   << 
692     {                                             407     {
693       G4cout << Arb_x[i] << Arb_y[i] << Area_s << 408       G4cout << "Leaving LinearInterpolation" << G4endl;
694              << Arb_grad[i] << G4endl;         << 409       ArbEnergyH.DumpValues();
                                                   >> 410       IPDFArbEnergyH.DumpValues();
695     }                                             411     }
696     ++i;                                       << 
697   }                                            << 
698                                                << 
699   i = 0;                                       << 
700   while (i < maxi)                             << 
701   {                                            << 
702     Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum; / << 
703     IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_ << 
704     ++i;                                       << 
705   }                                            << 
706                                                << 
707   // now scale the ArbEnergyH, needed by Proba << 
708   //                                           << 
709   ArbEnergyH.ScaleVector(1., 1./sum);          << 
710                                                << 
711   if (verbosityLevel >= 1)                     << 
712   {                                            << 
713     G4cout << "Leaving LinearInterpolation" << << 
714     ArbEnergyH.DumpValues();                   << 
715     IPDFArbEnergyH.DumpValues();               << 
716   }                                            << 
717 }                                                 412 }
718                                                   413 
719 void G4SPSEneDistribution::LogInterpolation()  << 414 void G4SPSEneDistribution::LogInterpolation()
720 {                                                 415 {
721   // Interpolation based on Logarithmic equati    416   // Interpolation based on Logarithmic equations
722   // Generate equations of line segments          417   // Generate equations of line segments
723   // y = Ax**alpha => log y = alpha*logx + log    418   // y = Ax**alpha => log y = alpha*logx + logA
724   // Find area under line segments                419   // Find area under line segments
725   // Create normalised, cumulative array Arb_C << 420   // create normalised, cumulative array Arb_Cum_Area
726                                                << 
727   G4double Area_seg[1024]; // Stores area unde    421   G4double Area_seg[1024]; // Stores area under each segment
728   G4double sum = 0., Arb_x[1024]={0.}, Arb_y[1 << 422   G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
729   std::size_t i, count;                        << 423   G4int i, count;
730   std::size_t maxi = ArbEnergyH.GetVectorLengt << 424   G4int maxi = ArbEnergyH.GetVectorLength();
731   for (i = 0; i < maxi; ++i)                   << 425   for(i=0;i<maxi;i++) {
732   {                                            << 426     Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
733     Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(i); << 427     Arb_y[i] = ArbEnergyH(size_t(i));
734     Arb_y[i] = ArbEnergyH(i);                  << 
735   }                                               428   }
736                                                << 
737   // Points are now in x,y arrays. If the spec    429   // Points are now in x,y arrays. If the spectrum is integral it has to be
738   // made differential and if momentum it has  << 430   // made differential and if momentum it has to be made energy.
739                                                << 431   if(DiffSpec == false) {
740   if (!DiffSpec)                               << 
741   {                                            << 
742     // Converts integral point-wise spectra to    432     // Converts integral point-wise spectra to Differential
743     //                                         << 433     for( count=0;count<maxi-1;count++) {
744     for (count = 0; count < maxi-1; ++count)   << 434       Arb_y[count] = (Arb_y[count] - Arb_y[count+1])/(Arb_x[count+1]-Arb_x[count]);
745     {                                          << 
746       Arb_y[count] = (Arb_y[count] - Arb_y[cou << 
747                    / (Arb_x[count + 1] - Arb_x << 
748     }                                             435     }
749     --maxi;                                    << 436     maxi--;
750   }                                               437   }
751                                                << 438   //
752   if (!EnergySpec)                             << 439   if(EnergySpec == false) {
753   {                                            << 440     // change currently stored values (emin etc) which are actually momenta
754     // Change currently stored values (emin et << 441     // to energies.
755     // to energies                             << 442     if(particle_definition == NULL)
756                                                << 443       G4cout << "Error: particle not defined" << G4endl;
757     G4ParticleDefinition* pdef = threadLocalDa << 444     else {
758     if (pdef == nullptr)                       << 
759     {                                          << 
760       G4Exception("G4SPSEneDistribution::LogIn << 
761                   "Event0302", FatalException, << 
762                   "Error: particle not defined << 
763     }                                          << 
764     else                                       << 
765     {                                          << 
766       // Apply Energy**2 = p**2c**2 + m0**2c**    445       // Apply Energy**2 = p**2c**2 + m0**2c**4
767       // p should be entered as E/c i.e. witho    446       // p should be entered as E/c i.e. without the division by c
768       // being done - energy equivalent        << 447       // being done - energy equivalent.
769                                                << 448       G4double mass = particle_definition->GetPDGMass();
770       G4double mass = pdef->GetPDGMass();      << 449       // convert point to energy unit and its value to per energy unit
771                                                << 
772       // Convert point to energy unit and its  << 
773       //                                       << 
774       G4double total_energy;                      450       G4double total_energy;
775       for (count = 0; count < maxi; ++count)   << 451       for(count=0;count<maxi;count++) {
776       {                                        << 452   total_energy = std::sqrt((Arb_x[count]*Arb_x[count]) 
777         total_energy = std::sqrt((Arb_x[count] << 453           + (mass*mass)); // total energy
778         Arb_y[count] = Arb_y[count] * Arb_x[co << 454   
779         Arb_x[count] = total_energy - mass; // << 455   Arb_y[count] = Arb_y[count] * Arb_x[count]/total_energy; 
                                                   >> 456   Arb_x[count] = total_energy - mass ;  // kinetic energy
780       }                                           457       }
781     }                                             458     }
782   }                                               459   }
783                                                << 460   //
784   i = 1;                                       << 461   i=1;
785                                                << 
786   if ( Arb_ezero != nullptr ) { delete [] Arb_ << 
787   if ( Arb_Const != nullptr ) { delete [] Arb_ << 
788   Arb_alpha = new G4double [1024];             << 
789   Arb_Const = new G4double [1024];             << 
790   Arb_alpha_Const_flag = true;                 << 
791                                                << 
792   Arb_alpha[0] = 0.;                              462   Arb_alpha[0] = 0.;
793   Arb_Const[0] = 0.;                              463   Arb_Const[0] = 0.;
794   Area_seg[0] = 0.;                               464   Area_seg[0] = 0.;
795   Arb_Cum_Area[0] = 0.;                        << 465   if(Arb_x[0] <= 0. || Arb_y[0] <= 0.)
796   if (Arb_x[0] <= 0. || Arb_y[0] <= 0.)        << 
797   {                                            << 
798     G4cout << "You should not use log interpol << 
799            << G4endl;                          << 
800     G4cout << "These will be changed to 1e-20, << 
801            << G4endl;                          << 
802     if (Arb_x[0] <= 0.)                        << 
803     {                                          << 
804       Arb_x[0] = 1e-20;                        << 
805     }                                          << 
806     if (Arb_y[0] <= 0.)                        << 
807     {                                             466     {
808       Arb_y[0] = 1e-20;                        << 467       G4cout << "You should not use log interpolation with points <= 0." << G4endl;
                                                   >> 468       G4cout << "These will be changed to 1e-20, which may cause problems" << G4endl;
                                                   >> 469       if(Arb_x[0] <= 0.)
                                                   >> 470   Arb_x[0] = 1e-20;
                                                   >> 471       if(Arb_y[0] <= 0.)
                                                   >> 472   Arb_y[0] = 1e-20;
                                                   >> 473     }
                                                   >> 474 
                                                   >> 475   G4double alp;  
                                                   >> 476   while(i <maxi)
                                                   >> 477     {
                                                   >> 478       // Incase points are negative or zero
                                                   >> 479       if(Arb_x[i] <= 0. || Arb_y[i] <= 0.)
                                                   >> 480   {
                                                   >> 481     G4cout << "You should not use log interpolation with points <= 0." << G4endl;
                                                   >> 482     G4cout << "These will be changed to 1e-20, which may cause problems" << G4endl;
                                                   >> 483     if(Arb_x[i] <= 0.)
                                                   >> 484       Arb_x[i] = 1e-20;
                                                   >> 485     if(Arb_y[i] <= 0.)
                                                   >> 486       Arb_y[i] = 1e-20;
                                                   >> 487   }
                                                   >> 488 
                                                   >> 489       Arb_alpha[i] = (std::log10(Arb_y[i])-std::log10(Arb_y[i-1]))/(std::log10(Arb_x[i])-std::log10(Arb_x[i-1]));
                                                   >> 490       Arb_Const[i] = Arb_y[i]/(std::pow(Arb_x[i],Arb_alpha[i]));
                                                   >> 491       alp = Arb_alpha[i] + 1;
                                                   >> 492       Area_seg[i] = (Arb_Const[i]/alp) * (std::pow(Arb_x[i],alp) - std::pow(Arb_x[i-1],alp));
                                                   >> 493       sum = sum + Area_seg[i];
                                                   >> 494       Arb_Cum_Area[i] = Arb_Cum_Area[i-1] + Area_seg[i];
                                                   >> 495       if(verbosityLevel == 2)
                                                   >> 496   G4cout << Arb_alpha[i] << Arb_Const[i] << Area_seg[i] << G4endl;
                                                   >> 497       i++;
809     }                                             498     }
810   }                                            << 499   
811                                                << 500   i=0;
812   G4double alp;                                << 501   while(i<maxi)
813   while (i < maxi)                             << 
814   {                                            << 
815     // In case points are negative or zero     << 
816     //                                         << 
817     if (Arb_x[i] <= 0. || Arb_y[i] <= 0.)      << 
818     {                                          << 
819       G4cout << "You should not use log interp << 
820              << G4endl;                        << 
821       G4cout << "These will be changed to 1e-2 << 
822              << G4endl;                        << 
823       if (Arb_x[i] <= 0.)                      << 
824       {                                        << 
825         Arb_x[i] = 1e-20;                      << 
826       }                                        << 
827       if (Arb_y[i] <= 0.)                      << 
828       {                                        << 
829         Arb_y[i] = 1e-20;                      << 
830       }                                        << 
831     }                                          << 
832                                                << 
833     Arb_alpha[i] = (std::log10(Arb_y[i]) - std << 
834                  / (std::log10(Arb_x[i]) - std << 
835     Arb_Const[i] = Arb_y[i] / (std::pow(Arb_x[ << 
836     alp = Arb_alpha[i] + 1;                    << 
837     if (alp == 0.)                             << 
838     {                                          << 
839       Area_seg[i] = Arb_Const[i]               << 
840                   * (std::log(Arb_x[i]) - std: << 
841     }                                          << 
842     else                                       << 
843     {                                          << 
844       Area_seg[i] = (Arb_Const[i] / alp)       << 
845                   * (std::pow(Arb_x[i], alp) - << 
846     }                                          << 
847     sum = sum + Area_seg[i];                   << 
848     Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Ar << 
849     if (verbosityLevel == 2)                   << 
850     {                                             502     {
851       G4cout << Arb_alpha[i] << Arb_Const[i] < << 503       Arb_Cum_Area[i] = Arb_Cum_Area[i]/sum;
                                                   >> 504       IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
                                                   >> 505       i++;
852     }                                             506     }
853     ++i;                                       << 507   if(verbosityLevel >= 1)
854   }                                            << 
855                                                << 
856   i = 0;                                       << 
857   while (i < maxi)                             << 
858   {                                            << 
859     Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum;   << 
860     IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_ << 
861     ++i;                                       << 
862   }                                            << 
863                                                << 
864   // Now scale the ArbEnergyH, needed by Proba << 
865   //                                           << 
866   ArbEnergyH.ScaleVector(1., 1./sum);          << 
867                                                << 
868   if (verbosityLevel >= 1)                     << 
869   {                                            << 
870     G4cout << "Leaving LogInterpolation " << G    508     G4cout << "Leaving LogInterpolation " << G4endl;
871   }                                            << 
872 }                                                 509 }
873                                                   510 
874 void G4SPSEneDistribution::ExpInterpolation()  << 511 void G4SPSEneDistribution::ExpInterpolation()
875 {                                                 512 {
876   // Interpolation based on Exponential equati    513   // Interpolation based on Exponential equations
877   // Generate equations of line segments          514   // Generate equations of line segments
878   // y = Ae**-(x/e0) => ln y = -x/e0 + lnA        515   // y = Ae**-(x/e0) => ln y = -x/e0 + lnA
879   // Find area under line segments                516   // Find area under line segments
880   // Create normalised, cumulative array Arb_C << 517   // create normalised, cumulative array Arb_Cum_Area
881                                                << 
882   G4double Area_seg[1024]; // Stores area unde    518   G4double Area_seg[1024]; // Stores area under each segment
883   G4double sum = 0., Arb_x[1024]={0.}, Arb_y[1 << 519   G4double sum = 0., Arb_x[1024], Arb_y[1024], Arb_Cum_Area[1024];
884   std::size_t i, count;                        << 520   G4int i, count;
885   std::size_t maxi = ArbEnergyH.GetVectorLengt << 521   G4int maxi = ArbEnergyH.GetVectorLength();
886   for (i = 0; i < maxi; ++i)                   << 522   for(i=0;i<maxi;i++) {
887   {                                            << 523     Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
888     Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(i); << 524     Arb_y[i] = ArbEnergyH(size_t(i));
889     Arb_y[i] = ArbEnergyH(i);                  << 
890   }                                               525   }
891                                                << 
892   // Points are now in x,y arrays. If the spec    526   // Points are now in x,y arrays. If the spectrum is integral it has to be
893   // made differential and if momentum it has  << 527   // made differential and if momentum it has to be made energy.
894                                                << 528   if(DiffSpec == false) {
895   if (!DiffSpec)                               << 
896   {                                            << 
897     // Converts integral point-wise spectra to    529     // Converts integral point-wise spectra to Differential
898     //                                         << 530     for( count=0;count< maxi-1;count++) {
899     for (count = 0; count < maxi - 1; ++count) << 531       Arb_y[count] = (Arb_y[count] - Arb_y[count+1])/(Arb_x[count+1]-Arb_x[count]);
900     {                                          << 
901       Arb_y[count] = (Arb_y[count] - Arb_y[cou << 
902                    / (Arb_x[count + 1] - Arb_x << 
903     }                                             532     }
904     --maxi;                                    << 533     maxi--;
905   }                                               534   }
906                                                << 535   //
907   if (!EnergySpec)                             << 536   if(EnergySpec == false) {
908   {                                            << 537     // change currently stored values (emin etc) which are actually momenta
909     // Change currently stored values (emin et << 538     // to energies.
910     // to energies                             << 539     if(particle_definition == NULL)
911     //                                         << 540       G4cout << "Error: particle not defined" << G4endl;
912     G4ParticleDefinition* pdef = threadLocalDa << 541     else {
913     if (pdef == nullptr)                       << 
914     {                                          << 
915       G4Exception("G4SPSEneDistribution::ExpIn << 
916                   "Event0302", FatalException, << 
917                   "Error: particle not defined << 
918     }                                          << 
919     else                                       << 
920     {                                          << 
921       // Apply Energy**2 = p**2c**2 + m0**2c**    542       // Apply Energy**2 = p**2c**2 + m0**2c**4
922       // p should be entered as E/c i.e. witho    543       // p should be entered as E/c i.e. without the division by c
923       // being done - energy equivalent        << 544       // being done - energy equivalent.
924                                                << 545       G4double mass = particle_definition->GetPDGMass();
925       G4double mass = pdef->GetPDGMass();      << 546       // convert point to energy unit and its value to per energy unit
926                                                << 
927       // Convert point to energy unit and its  << 
928       //                                       << 
929       G4double total_energy;                      547       G4double total_energy;
930       for (count = 0; count < maxi; ++count)   << 548       for(count=0;count<maxi;count++) {
931       {                                        << 549   total_energy = std::sqrt((Arb_x[count]*Arb_x[count]) 
932         total_energy = std::sqrt((Arb_x[count] << 550           + (mass*mass)); // total energy
933         Arb_y[count] = Arb_y[count] * Arb_x[co << 551   
934         Arb_x[count] = total_energy - mass; // << 552   Arb_y[count] = Arb_y[count] * Arb_x[count]/total_energy; 
                                                   >> 553   Arb_x[count] = total_energy - mass ;  // kinetic energy
935       }                                           554       }
936     }                                             555     }
937   }                                               556   }
938                                                << 557   //
939   i = 1;                                       << 558   i=1;
940                                                << 
941   if ( Arb_ezero != nullptr ) { delete[] Arb_e << 
942   if ( Arb_Const != nullptr ) { delete[] Arb_C << 
943   Arb_ezero = new G4double [1024];             << 
944   Arb_Const = new G4double [1024];             << 
945   Arb_ezero_flag = true;                       << 
946                                                << 
947   Arb_ezero[0] = 0.;                              559   Arb_ezero[0] = 0.;
948   Arb_Const[0] = 0.;                              560   Arb_Const[0] = 0.;
949   Area_seg[0] = 0.;                               561   Area_seg[0] = 0.;
950   Arb_Cum_Area[0] = 0.;                           562   Arb_Cum_Area[0] = 0.;
951   while (i < maxi)                             << 563   while(i < maxi)
952   {                                            << 
953     G4double test = std::log(Arb_y[i]) - std:: << 
954     if (test > 0. || test < 0.)                << 
955     {                                          << 
956       Arb_ezero[i] = -(Arb_x[i] - Arb_x[i - 1] << 
957                    / (std::log(Arb_y[i]) - std << 
958       Arb_Const[i] = Arb_y[i] / (std::exp(-Arb << 
959       Area_seg[i] = -(Arb_Const[i] * Arb_ezero << 
960                   * (std::exp(-Arb_x[i] / Arb_ << 
961                    - std::exp(-Arb_x[i - 1] /  << 
962     }                                          << 
963     else                                       << 
964     {                                             564     {
965       G4Exception("G4SPSEneDistribution::ExpIn << 565       G4double test = std::log(Arb_y[i]) - std::log(Arb_y[i-1]);
966                   "Event0302", JustWarning,    << 566       if(test > 0. || test < 0.)
967                   "Flat line segment: problem, << 567   {
968       G4cout << "Flat line segment: problem" < << 568     Arb_ezero[i] = -(Arb_x[i] - Arb_x[i-1])/(std::log(Arb_y[i]) - std::log(Arb_y[i-1]));
969       Arb_ezero[i] = 0.;                       << 569     Arb_Const[i] = Arb_y[i]/(std::exp(-Arb_x[i]/Arb_ezero[i]));
970       Arb_Const[i] = 0.;                       << 570     Area_seg[i]=-(Arb_Const[i]*Arb_ezero[i])*(std::exp(-Arb_x[i]/Arb_ezero[i])
971       Area_seg[i] = 0.;                        << 571                 -std::exp(-Arb_x[i-1]/Arb_ezero[i]));
972     }                                          << 572   }
973     sum = sum + Area_seg[i];                   << 573       else 
974     Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + Ar << 574   {
975     if (verbosityLevel == 2)                   << 575     G4cout << "Flat line segment: problem" << G4endl;
                                                   >> 576     Arb_ezero[i] = 0.;
                                                   >> 577     Arb_Const[i] = 0.;
                                                   >> 578     Area_seg[i] = 0.;
                                                   >> 579   }
                                                   >> 580       sum = sum + Area_seg[i];
                                                   >> 581       Arb_Cum_Area[i] = Arb_Cum_Area[i-1] + Area_seg[i];
                                                   >> 582       if(verbosityLevel == 2)
                                                   >> 583   G4cout << Arb_ezero[i] << Arb_Const[i] << Area_seg[i] << G4endl;
                                                   >> 584       i++;
                                                   >> 585     }
                                                   >> 586   
                                                   >> 587   i=0;
                                                   >> 588   while(i<maxi)
976     {                                             589     {
977       G4cout << Arb_ezero[i] << Arb_Const[i] < << 590       Arb_Cum_Area[i] = Arb_Cum_Area[i]/sum;
                                                   >> 591       IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]);
                                                   >> 592       i++;
978     }                                             593     }
979     ++i;                                       << 594   if(verbosityLevel >= 1)
980   }                                            << 
981                                                << 
982   i = 0;                                       << 
983   while (i < maxi)                             << 
984   {                                            << 
985     Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum;   << 
986     IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_ << 
987     ++i;                                       << 
988   }                                            << 
989                                                << 
990   // Now scale the ArbEnergyH, needed by Proba << 
991   //                                           << 
992   ArbEnergyH.ScaleVector(1., 1./sum);          << 
993                                                << 
994   if (verbosityLevel >= 1)                     << 
995   {                                            << 
996     G4cout << "Leaving ExpInterpolation " << G    595     G4cout << "Leaving ExpInterpolation " << G4endl;
997   }                                            << 
998 }                                                 596 }
999                                                   597 
1000 void G4SPSEneDistribution::SplineInterpolatio << 598 void G4SPSEneDistribution::SplineInterpolation()
1001 {                                                599 {
1002   // Interpolation using Splines.                600   // Interpolation using Splines.
1003   // Create Normalised arrays, make x 0->1 an << 601   // Create Normalised arrays, make x 0->1 and y hold
1004   //                                          << 602   // the function (Energy)
1005   // Current method based on the above will n << 603   G4double  Arb_x[1024], Arb_y[1024];
1006   // New method is implemented below.         << 604   G4int i, count;
1007                                               << 605   G4int maxi = ArbEnergyH.GetVectorLength();
1008   G4double sum, Arb_x[1024]={0.}, Arb_y[1024] << 606   for(i=0;i<maxi;i++) {
1009   std::size_t i, count;                       << 607     Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i));
1010   std::size_t maxi = ArbEnergyH.GetVectorLeng << 608     Arb_y[i] = ArbEnergyH(size_t(i));
1011                                               << 
1012   for (i = 0; i < maxi; ++i)                  << 
1013   {                                           << 
1014     Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(i) << 
1015     Arb_y[i] = ArbEnergyH(i);                 << 
1016   }                                              609   }
1017                                               << 
1018   // Points are now in x,y arrays. If the spe    610   // Points are now in x,y arrays. If the spectrum is integral it has to be
1019   // made differential and if momentum it has << 611   // made differential and if momentum it has to be made energy.
1020                                               << 612   if(DiffSpec == false) {
1021   if (!DiffSpec)                              << 
1022   {                                           << 
1023     // Converts integral point-wise spectra t    613     // Converts integral point-wise spectra to Differential
1024     //                                        << 614     for( count=0;count< maxi-1;count++) {
1025     for (count = 0; count < maxi - 1; ++count << 615       Arb_y[count] = (Arb_y[count] - Arb_y[count+1])/(Arb_x[count+1]-Arb_x[count]);
1026     {                                         << 
1027       Arb_y[count] = (Arb_y[count] - Arb_y[co << 
1028                    / (Arb_x[count + 1] - Arb_ << 
1029     }                                            616     }
1030     --maxi;                                   << 617     maxi--;
1031   }                                              618   }
1032                                               << 619   //
1033   if (!EnergySpec)                            << 620   if(EnergySpec == false) {
1034   {                                           << 621     // change currently stored values (emin etc) which are actually momenta
1035     // Change currently stored values (emin e << 622     // to energies.
1036     // to energies                            << 623     if(particle_definition == NULL)
1037     //                                        << 624       G4cout << "Error: particle not defined" << G4endl;
1038     G4ParticleDefinition* pdef = threadLocalD << 625     else {
1039     if (pdef == nullptr)                      << 
1040     {                                         << 
1041       G4Exception("G4SPSEneDistribution::Spli << 
1042                   "Event0302", FatalException << 
1043                   "Error: particle not define << 
1044     }                                         << 
1045     else                                      << 
1046     {                                         << 
1047       // Apply Energy**2 = p**2c**2 + m0**2c*    626       // Apply Energy**2 = p**2c**2 + m0**2c**4
1048       // p should be entered as E/c i.e. with    627       // p should be entered as E/c i.e. without the division by c
1049       // being done - energy equivalent       << 628       // being done - energy equivalent.
1050                                               << 629       G4double mass = particle_definition->GetPDGMass();
1051       G4double mass = pdef->GetPDGMass();     << 630       // convert point to energy unit and its value to per energy unit
1052                                               << 
1053       // Convert point to energy unit and its << 
1054       //                                      << 
1055       G4double total_energy;                     631       G4double total_energy;
1056       for (count = 0; count < maxi; ++count)  << 632       for(count=0;count<maxi;count++) {
1057       {                                       << 633   total_energy = std::sqrt((Arb_x[count]*Arb_x[count]) 
1058         total_energy = std::sqrt((Arb_x[count << 634           + (mass*mass)); // total energy
1059         Arb_y[count] = Arb_y[count] * Arb_x[c << 635   
1060         Arb_x[count] = total_energy - mass; / << 636   Arb_y[count] = Arb_y[count] * Arb_x[count]/total_energy; 
                                                   >> 637   Arb_x[count] = total_energy - mass ;  // kinetic energy
1061       }                                          638       }
1062     }                                            639     }
1063   }                                              640   }
1064                                               << 641   //
1065   i = 1;                                      << 642   for(i=1;i<maxi;i++)
1066   Arb_Cum_Area[0] = 0.;                       << 643     Arb_y[i] += Arb_y[i-1];
1067   sum = 0.;                                   << 644 
1068   Splinetemp = new G4DataInterpolation(Arb_x, << 645   for(i=0;i<maxi;i++)
1069   G4double ei[101], prob[101];                << 646     Arb_y[i] /= Arb_y[maxi-1];
1070   for (auto & it : SplineInt)                 << 647   // now Arb_y is accumulated normalised probabilities
1071   {                                           << 648   /*  for(i=0; i<maxi;i++) {
1072     delete it;                                << 649       if(verbosityLevel >1) 
1073     it = nullptr;                             << 650        G4cout << i <<" "<< Arb_x[i] << " " << Arb_y[i] << G4endl;
1074   }                                           << 651        IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_y[i]);    
1075   SplineInt.clear();                          << 652    }
1076   SplineInt.resize(1024,nullptr);             << 653    Emax = IPDFArbEnergyH.GetLowEdgeEnergy(IPDFArbEnergyH.GetVectorLength()-1);
1077   while (i < maxi)                            << 654    Emin = IPDFArbEnergyH.GetLowEdgeEnergy(0);
1078   {                                           << 655   */
1079     // 100 step per segment for the integrati << 656   // Should now have normalised cumulative probabilities in Arb_y
1080                                               << 657   // and energy values in Arb_x.
1081     G4double de = (Arb_x[i] - Arb_x[i - 1])/1 << 658   // maxi = maxi + 1;
1082     G4double area = 0.;                       << 659   // Put y into x and x into y. The spline interpolation will then
1083                                               << 660   // go through x-axis to find where to interpolate (cum probability)
1084     for (count = 0; count < 101; ++count)     << 661   // then generate a y (which will now be energy).
1085     {                                         << 662   SplineInt = new G4DataInterpolation(Arb_y,Arb_x,maxi,1e30,1e30);
1086       ei[count] = Arb_x[i - 1] + de*count ;   << 663   if(verbosityLevel >1 )
1087       prob[count] =  Splinetemp->CubicSplineI << 
1088       if (prob[count] < 0.)                   << 
1089       {                                       << 
1090         G4ExceptionDescription ED;            << 
1091         ED << "Warning: G4DataInterpolation r << 
1092            << " " << ei[count] << G4endl;     << 
1093         G4Exception("G4SPSEneDistribution::Sp << 
1094                     FatalException, ED);      << 
1095       }                                       << 
1096       area += prob[count]*de;                 << 
1097     }                                         << 
1098     Arb_Cum_Area[i] = Arb_Cum_Area[i - 1] + a << 
1099     sum += area;                              << 
1100                                               << 
1101     prob[0] = prob[0]/(area/de);              << 
1102     for (count = 1; count < 100; ++count)     << 
1103     {                                            664     {
1104       prob[count] = prob[count-1] + prob[coun << 665       G4cout << SplineInt << G4endl;
                                                   >> 666       G4cout << SplineInt->LocateArgument(1.0) << G4endl;
1105     }                                            667     }
1106                                               << 668   if(verbosityLevel > 0 )
1107     SplineInt[i] = new G4DataInterpolation(pr << 
1108                                               << 
1109     // NOTE: i starts from 1!                 << 
1110     //                                        << 
1111     ++i;                                      << 
1112   }                                           << 
1113                                               << 
1114   i = 0;                                      << 
1115   while (i < maxi)                            << 
1116   {                                           << 
1117     Arb_Cum_Area[i] = Arb_Cum_Area[i] / sum;  << 
1118     IPDFArbEnergyH.InsertValues(Arb_x[i], Arb << 
1119     ++i;                                      << 
1120   }                                           << 
1121                                               << 
1122   // Now scale the ArbEnergyH, needed by Prob << 
1123   //                                          << 
1124   ArbEnergyH.ScaleVector(1., 1./sum);         << 
1125                                               << 
1126   if (verbosityLevel > 0)                     << 
1127   {                                           << 
1128     G4cout << "Leaving SplineInterpolation "     669     G4cout << "Leaving SplineInterpolation " << G4endl;
1129   }                                           << 
1130 }                                                670 }
1131                                                  671 
1132 void G4SPSEneDistribution::GenerateMonoEnerge    672 void G4SPSEneDistribution::GenerateMonoEnergetic()
1133 {                                                673 {
1134   // Method to generate MonoEnergetic particl << 674   // Method to generate MonoEnergetic particles.
1135                                               << 675   particle_energy = MonoEnergy;
1136   threadLocalData.Get().particle_energy = Mon << 
1137 }                                                676 }
1138                                                  677 
1139 void G4SPSEneDistribution::GenerateGaussEnerg    678 void G4SPSEneDistribution::GenerateGaussEnergies()
1140 {                                                679 {
1141   // Method to generate Gaussian particles    << 680   // Method to generate Gaussian particles.
1142                                               << 681   particle_energy = G4RandGauss::shoot(MonoEnergy,SE);
1143   G4double ene = G4RandGauss::shoot(MonoEnerg << 682   if (particle_energy < 0) particle_energy = 0.;
1144   if (ene < 0) ene = 0.;                      << 
1145   threadLocalData.Get().particle_energy = ene << 
1146 }                                                683 }
1147                                                  684 
1148 void G4SPSEneDistribution::GenerateLinearEner    685 void G4SPSEneDistribution::GenerateLinearEnergies(G4bool bArb = false)
1149 {                                                686 {
1150   G4double rndm;                                 687   G4double rndm;
1151   threadLocal_t& params = threadLocalData.Get << 688   G4double emaxsq = std::pow(Emax,2.); //Emax squared
1152   G4double emaxsq = std::pow(params.Emax, 2.) << 689   G4double eminsq = std::pow(Emin,2.); //Emin squared
1153   G4double eminsq = std::pow(params.Emin, 2.) << 690   G4double intersq = std::pow(cept,2.); //cept squared
1154   G4double intersq = std::pow(params.cept, 2. << 
1155                                                  691 
1156   if (bArb) rndm = G4UniformRand();              692   if (bArb) rndm = G4UniformRand();
1157   else      rndm = eneRndm->GenRandEnergy();  << 693   else rndm = eneRndm->GenRandEnergy();
1158                                               << 694   
1159   G4double bracket = ((params.grad / 2.)      << 695   G4double bracket = ((grad/2.)*(emaxsq - eminsq) + cept*(Emax-Emin));
1160                    * (emaxsq - eminsq)        << 
1161                    + params.cept * (params.Em << 
1162   bracket = bracket * rndm;                      696   bracket = bracket * rndm;
1163   bracket = bracket + (params.grad / 2.) * em << 697   bracket = bracket + (grad/2.)*eminsq + cept*Emin;
1164                                               << 
1165   // Now have a quad of form m/2 E**2 + cE -     698   // Now have a quad of form m/2 E**2 + cE - bracket = 0
1166   //                                          << 
1167   bracket = -bracket;                            699   bracket = -bracket;
1168                                               << 700   //  G4cout << "BRACKET" << bracket << G4endl;
1169   if (params.grad != 0.)                      << 701   if(grad != 0.)
1170   {                                           << 
1171     G4double sqbrack = (intersq - 4 * (params << 
1172     sqbrack = std::sqrt(sqbrack);             << 
1173     G4double root1 = -params.cept + sqbrack;  << 
1174     root1 = root1 / (2. * (params.grad / 2.)) << 
1175                                               << 
1176     G4double root2 = -params.cept - sqbrack;  << 
1177     root2 = root2 / (2. * (params.grad / 2.)) << 
1178                                               << 
1179     if (root1 > params.Emin && root1 < params << 
1180     {                                            702     {
1181       params.particle_energy = root1;         << 703       G4double sqbrack = (intersq - 4*(grad/2.)*(bracket));
                                                   >> 704       //      G4cout << "SQBRACK" << sqbrack << G4endl;
                                                   >> 705       sqbrack = std::sqrt(sqbrack);
                                                   >> 706       G4double root1 = -cept + sqbrack; 
                                                   >> 707       root1 = root1/(2.*(grad/2.));
                                                   >> 708       
                                                   >> 709       G4double root2 = -cept - sqbrack;
                                                   >> 710       root2 = root2/(2.*(grad/2.));
                                                   >> 711 
                                                   >> 712       //      G4cout << root1 << " roots " << root2 << G4endl;
                                                   >> 713 
                                                   >> 714       if(root1 > Emin && root1 < Emax)
                                                   >> 715   particle_energy = root1;
                                                   >> 716       if(root2 > Emin && root2 < Emax)
                                                   >> 717   particle_energy = root2;
1182     }                                            718     }
1183     if (root2 > params.Emin && root2 < params << 719   else if(grad == 0.)
1184     {                                         << 
1185       params.particle_energy = root2;         << 
1186     }                                         << 
1187   }                                           << 
1188   else if (params.grad == 0.)                 << 
1189   {                                           << 
1190     // have equation of form cE - bracket =0     720     // have equation of form cE - bracket =0
1191     //                                        << 721     particle_energy = bracket/cept;
1192     params.particle_energy = bracket / params << 
1193   }                                           << 
1194                                               << 
1195   if (params.particle_energy < 0.)            << 
1196   {                                           << 
1197     params.particle_energy = -params.particle << 
1198   }                                           << 
1199                                                  722 
1200   if (verbosityLevel >= 1)                    << 723   if(particle_energy < 0.)
1201   {                                           << 724     particle_energy = -particle_energy;
1202     G4cout << "Energy is " << params.particle << 725   
1203   }                                           << 726   if(verbosityLevel >= 1)
                                                   >> 727     G4cout << "Energy is " << particle_energy << G4endl;
1204 }                                                728 }
1205                                                  729 
1206 void G4SPSEneDistribution::GeneratePowEnergie    730 void G4SPSEneDistribution::GeneratePowEnergies(G4bool bArb = false)
1207 {                                                731 {
1208   // Method to generate particle energies dis << 732   // Method to generate particle energies distributed as
                                                   >> 733   // a powerlaw
1209                                                  734 
1210   G4double rndm;                                 735   G4double rndm;
1211   G4double emina, emaxa;                         736   G4double emina, emaxa;
1212                                               << 
1213   threadLocal_t& params = threadLocalData.Get << 
1214                                               << 
1215   emina = std::pow(params.Emin, params.alpha  << 
1216   emaxa = std::pow(params.Emax, params.alpha  << 
1217                                                  737 
1218   if (bArb) rndm = G4UniformRand();           << 738   emina = std::pow(Emin,alpha+1);
1219   else      rndm = eneRndm->GenRandEnergy();  << 739   emaxa = std::pow(Emax,alpha+1);
1220                                                  740 
1221   if (params.alpha != -1.)                    << 741   if (bArb) rndm = G4UniformRand();
1222   {                                           << 742   else rndm = eneRndm->GenRandEnergy();
1223     G4double ene = ((rndm * (emaxa - emina))  << 743   
1224     ene = std::pow(ene, (1. / (params.alpha + << 744   if(alpha != -1.)
1225     params.particle_energy = ene;             << 
1226   }                                           << 
1227   else                                        << 
1228   {                                           << 
1229     G4double ene = (std::log(params.Emin)     << 
1230                  + rndm * (std::log(params.Em << 
1231     params.particle_energy = std::exp(ene);   << 
1232   }                                           << 
1233   if (verbosityLevel >= 1)                    << 
1234   {                                           << 
1235     G4cout << "Energy is " << params.particle << 
1236   }                                           << 
1237 }                                             << 
1238                                               << 
1239 void G4SPSEneDistribution::GenerateCPowEnergi << 
1240 {                                             << 
1241   // Method to generate particle energies dis << 
1242   // cutoff power-law distribution            << 
1243   //                                          << 
1244   // CP_x holds Energies, and CPHist holds th << 
1245   // binary search to find correct bin then l << 
1246   // Use the earlier defined histogram + Rand << 
1247   // random numbers following the histos dist << 
1248                                               << 
1249   G4double rndm = eneRndm->GenRandEnergy();   << 
1250   G4int nabove = 10001, nbelow = 0, middle;   << 
1251                                               << 
1252   G4AutoLock l(&mutex);                       << 
1253   G4bool done = CPhistCalcd;                  << 
1254   l.unlock();                                 << 
1255                                               << 
1256   if(!done)                                   << 
1257   {                                           << 
1258     Calculate(); //This is has a lock inside, << 
1259     l.lock();                                 << 
1260     CPhistCalcd = true;                       << 
1261     l.unlock();                               << 
1262   }                                           << 
1263                                               << 
1264   // Binary search to find bin that rndm is i << 
1265   //                                          << 
1266   while (nabove - nbelow > 1)                 << 
1267   {                                           << 
1268     middle = (nabove + nbelow) / 2;           << 
1269     if (rndm == CPHist->at(middle))           << 
1270     {                                         << 
1271       break;                                  << 
1272     }                                         << 
1273     if (rndm < CPHist->at(middle))            << 
1274     {                                            745     {
1275       nabove = middle;                        << 746       particle_energy = ((rndm*(emaxa - emina)) + emina);
                                                   >> 747       particle_energy = std::pow(particle_energy,(1./(alpha+1.)));
1276     }                                            748     }
1277     else                                      << 749   else if(alpha == -1.)
1278     {                                            750     {
1279       nbelow = middle;                        << 751       particle_energy = (std::log(Emin) + rndm*(std::log(Emax) - std::log(Emin)));
                                                   >> 752       particle_energy = std::exp(particle_energy);
1280     }                                            753     }
1281   }                                           << 754   if(verbosityLevel >= 1)
1282                                               << 755     G4cout << "Energy is " << particle_energy << G4endl;
1283   // Now interpolate in that bin to find the  << 
1284   //                                          << 
1285   G4double x1, x2, y1, y2, t, q;              << 
1286   x1 = CP_x->at(nbelow);                      << 
1287   if(nbelow+1 == static_cast<G4int>(CP_x->siz << 
1288   {                                           << 
1289     x2 = CP_x->back();                        << 
1290   }                                           << 
1291   else                                        << 
1292   {                                           << 
1293     x2 = CP_x->at(nbelow + 1);                << 
1294   }                                           << 
1295   y1 = CPHist->at(nbelow);                    << 
1296   if(nbelow+1 == static_cast<G4int>(CPHist->s << 
1297   {                                           << 
1298     G4cout << CPHist->back() << G4endl;       << 
1299     y2 = CPHist->back();                      << 
1300   }                                           << 
1301   else                                        << 
1302   {                                           << 
1303     y2 = CPHist->at(nbelow + 1);              << 
1304   }                                           << 
1305   t = (y2 - y1) / (x2 - x1);                  << 
1306   q = y1 - t * x1;                            << 
1307                                               << 
1308   threadLocalData.Get().particle_energy = (rn << 
1309                                               << 
1310   if (verbosityLevel >= 1)                    << 
1311   {                                           << 
1312     G4cout << "Energy is " << threadLocalData << 
1313   }                                           << 
1314 }                                             << 
1315                                               << 
1316 void G4SPSEneDistribution::GenerateBiasPowEne << 
1317 {                                             << 
1318   // Method to generate particle energies dis << 
1319   // in biased power-law and calculate its we << 
1320                                               << 
1321   threadLocal_t& params = threadLocalData.Get << 
1322                                               << 
1323   G4double rndm;                              << 
1324   G4double emina, emaxa, emin, emax;          << 
1325                                               << 
1326   G4double normal = 1.;                       << 
1327                                               << 
1328   emin = params.Emin;                         << 
1329   emax = params.Emax;                         << 
1330                                               << 
1331   rndm = eneRndm->GenRandEnergy();            << 
1332                                               << 
1333   if (biasalpha != -1.)                       << 
1334   {                                           << 
1335     emina = std::pow(emin, biasalpha + 1);    << 
1336     emaxa = std::pow(emax, biasalpha + 1);    << 
1337     G4double ee = ((rndm * (emaxa - emina)) + << 
1338     params.particle_energy = std::pow(ee, (1. << 
1339     normal = 1./(1+biasalpha) * (emaxa - emin << 
1340   }                                           << 
1341   else                                        << 
1342   {                                           << 
1343     G4double ee = (std::log(emin) + rndm * (s << 
1344     params.particle_energy = std::exp(ee);    << 
1345     normal = std::log(emax) - std::log(emin); << 
1346   }                                           << 
1347   params.weight = GetProbability(params.parti << 
1348                 / (std::pow(params.particle_e << 
1349                                               << 
1350   if (verbosityLevel >= 1)                    << 
1351   {                                           << 
1352     G4cout << "Energy is " << params.particle << 
1353   }                                           << 
1354 }                                                756 }
1355                                                  757 
1356 void G4SPSEneDistribution::GenerateExpEnergie    758 void G4SPSEneDistribution::GenerateExpEnergies(G4bool bArb = false)
1357 {                                                759 {
1358   // Method to generate particle energies dis    760   // Method to generate particle energies distributed according
1359   // to an exponential curve                  << 761   // to an exponential curve.
1360                                               << 
1361   G4double rndm;                                 762   G4double rndm;
1362                                               << 763  
1363   if (bArb) rndm = G4UniformRand();              764   if (bArb) rndm = G4UniformRand();
1364   else      rndm = eneRndm->GenRandEnergy();  << 765   else rndm = eneRndm->GenRandEnergy();
1365                                                  766 
1366   threadLocal_t& params = threadLocalData.Get << 767   particle_energy = -Ezero*(std::log(rndm*(std::exp(-Emax/Ezero) - std::exp(-Emin/Ezero)) + 
1367   params.particle_energy = -params.Ezero      << 768         std::exp(-Emin/Ezero)));
1368                          * (std::log(rndm * ( << 769   if(verbosityLevel >= 1)
1369                                               << 770     G4cout << "Energy is " << particle_energy << G4endl;
1370                                            -  << 
1371                                               << 
1372                                    + std::exp << 
1373   if (verbosityLevel >= 1)                    << 
1374   {                                           << 
1375     G4cout << "Energy is " << params.particle << 
1376   }                                           << 
1377 }                                                771 }
1378                                                  772 
1379 void G4SPSEneDistribution::GenerateBremEnergi    773 void G4SPSEneDistribution::GenerateBremEnergies()
1380 {                                                774 {
1381   // Method to generate particle energies dis << 775   // Method to generate particle energies distributed according 
1382   // to a Bremstrahlung equation of the form  << 776   // to a Bremstrahlung equation of 
1383   // I = const*((kT)**1/2)*E*(e**(-E/kT))     << 777   // form I = const*((kT)**1/2)*E*(e**(-E/kT))
1384                                               << 778   
1385   G4double rndm = eneRndm->GenRandEnergy();   << 779   G4double rndm;
                                                   >> 780   rndm = eneRndm->GenRandEnergy();
1386   G4double expmax, expmin, k;                    781   G4double expmax, expmin, k;
1387                                               << 
1388   k = 8.6181e-11; // Boltzmann's const in MeV << 
1389   G4double ksq = std::pow(k, 2.); // k square << 
1390   G4double Tsq = std::pow(Temp, 2.); // Temp  << 
1391                                                  782 
1392   threadLocal_t& params = threadLocalData.Get << 783   k = 8.6181e-11; // Boltzmann's const in MeV/K
                                                   >> 784   G4double ksq = std::pow(k,2.); // k squared
                                                   >> 785   G4double Tsq = std::pow(Temp,2.); // Temp squared
1393                                                  786 
1394   expmax = std::exp(-params.Emax / (k * Temp) << 787   expmax = std::exp(-Emax/(k*Temp));
1395   expmin = std::exp(-params.Emin / (k * Temp) << 788   expmin = std::exp(-Emin/(k*Temp));
1396                                                  789 
1397   // If either expmax or expmin are zero then    790   // If either expmax or expmin are zero then this will cause problems
1398   // Most probably this will be because T is     791   // Most probably this will be because T is too low or E is too high
1399                                                  792 
1400   if (expmax == 0.)                           << 793   if(expmax == 0.)
1401   {                                           << 794     G4cout << "*****EXPMAX=0. Choose different E's or Temp" << G4endl;
1402     G4Exception("G4SPSEneDistribution::Genera << 795   if(expmin == 0.)
1403                 "Event0302", FatalException,  << 796     G4cout << "*****EXPMIN=0. Choose different E's or Temp" << G4endl;
1404                 "*****EXPMAX=0. Choose differ << 
1405   }                                           << 
1406   if (expmin == 0.)                           << 
1407   {                                           << 
1408     G4Exception("G4SPSEneDistribution::Genera << 
1409                 "Event0302", FatalException,  << 
1410                 "*****EXPMIN=0. Choose differ << 
1411   }                                           << 
1412                                                  797 
1413   G4double tempvar = rndm * ((-k) * Temp * (p << 798   G4double tempvar = rndm *((-k)*Temp*(Emax*expmax - Emin*expmin) -
1414                                           - p << 799     (ksq*Tsq*(expmax-expmin)));
1415                           - (ksq * Tsq * (exp << 
1416                                                  800 
1417   G4double bigc = (tempvar - k * Temp * param << 801   G4double bigc = (tempvar - k*Temp*Emin*expmin - ksq*Tsq*expmin)/(-k*Temp);
1418                  - ksq * Tsq * expmin) / (-k  << 
1419                                                  802 
1420   // This gives an equation of form: Ee(-E/kT    803   // This gives an equation of form: Ee(-E/kT) + kTe(-E/kT) - C =0
1421   // Solve this iteratively, step from Emin t    804   // Solve this iteratively, step from Emin to Emax in 1000 steps
1422   // and take the best solution.                 805   // and take the best solution.
1423                                                  806 
1424   G4double erange = params.Emax - params.Emin << 807   G4double erange = Emax - Emin;
1425   G4double steps = erange / 1000.;            << 808   G4double steps = erange/1000.;
1426   G4int i;                                       809   G4int i;
1427   G4double etest, diff, err = 100000.;        << 810   G4double etest, diff, err;
1428                                               << 811   
1429   for (i = 1; i < 1000; ++i)                  << 812   err = 100000.;
1430   {                                           << 
1431     etest = params.Emin + (i - 1) * steps;    << 
1432     diff = etest * (std::exp(-etest / (k * Te << 
1433          + k * Temp * (std::exp(-etest / (k * << 
1434                                               << 
1435     if (diff < 0.)                            << 
1436     {                                         << 
1437       diff = -diff;                           << 
1438     }                                         << 
1439                                                  813 
1440     if (diff < err)                           << 814   for(i=1; i<1000; i++)
1441     {                                            815     {
1442       err = diff;                             << 816       etest = Emin + (i-1)*steps;
1443       params.particle_energy = etest;         << 817       
1444     }                                         << 818       diff = etest*(std::exp(-etest/(k*Temp))) + k*Temp*(std::exp(-etest/(k*Temp))) - bigc;
1445   }                                           << 819 
1446   if (verbosityLevel >= 1)                    << 820       if(diff < 0.)
1447   {                                           << 821   diff = -diff;
1448     G4cout << "Energy is " << params.particle << 822 
1449   }                                           << 823       if(diff < err)
                                                   >> 824   {
                                                   >> 825     err = diff;
                                                   >> 826     particle_energy = etest;
                                                   >> 827   }
                                                   >> 828     }  
                                                   >> 829   if(verbosityLevel >= 1)
                                                   >> 830     G4cout << "Energy is " << particle_energy << G4endl;
1450 }                                                831 }
1451                                                  832 
1452 void G4SPSEneDistribution::GenerateBbodyEnerg    833 void G4SPSEneDistribution::GenerateBbodyEnergies()
1453 {                                                834 {
1454   // BBody_x holds Energies, and BBHist holds    835   // BBody_x holds Energies, and BBHist holds the cumulative histo.
1455   // Binary search to find correct bin then l << 836   // binary search to find correct bin then lin interpolation.
1456   // Use the earlier defined histogram + Rand    837   // Use the earlier defined histogram + RandGeneral method to generate
1457   // random numbers following the histos dist << 838   // random numbers following the histos distribution.
1458                                               << 839   G4double rndm;
1459   G4double rndm = eneRndm->GenRandEnergy();   << 840   G4int nabove, nbelow = 0, middle;
1460   G4int nabove = 10001, nbelow = 0, middle;   << 841   nabove = 10001;
1461                                               << 842   rndm = eneRndm->GenRandEnergy();
1462   G4AutoLock l(&mutex);                       << 
1463   G4bool done = BBhistCalcd;                  << 
1464   l.unlock();                                 << 
1465                                               << 
1466   if(!done)                                   << 
1467   {                                           << 
1468     Calculate(); //This is has a lock inside, << 
1469     l.lock();                                 << 
1470     BBhistCalcd = true;                       << 
1471     l.unlock();                               << 
1472   }                                           << 
1473                                                  843 
1474   // Binary search to find bin that rndm is i    844   // Binary search to find bin that rndm is in
1475   //                                          << 845   while(nabove-nbelow > 1)
1476   while (nabove - nbelow > 1)                 << 
1477   {                                           << 
1478     middle = (nabove + nbelow) / 2;           << 
1479     if (rndm == BBHist->at(middle))           << 
1480     {                                         << 
1481       break;                                  << 
1482     }                                         << 
1483     if (rndm < BBHist->at(middle))            << 
1484     {                                            846     {
1485       nabove = middle;                        << 847       middle = (nabove + nbelow)/2;
                                                   >> 848       if(rndm == BBHist[middle]) break;
                                                   >> 849       if(rndm < BBHist[middle]) nabove = middle;
                                                   >> 850       else nbelow = middle;
1486     }                                            851     }
1487     else                                      << 852   
                                                   >> 853   // Now interpolate in that bin to find the correct output value.
                                                   >> 854   G4double x1, x2, y1, y2, m, q;
                                                   >> 855   x1 = Bbody_x[nbelow];
                                                   >> 856   x2 = Bbody_x[nbelow+1];
                                                   >> 857   y1 = BBHist[nbelow];
                                                   >> 858   y2 = BBHist[nbelow+1];
                                                   >> 859   m = (y2-y1)/(x2-x1);
                                                   >> 860   q = y1 - m*x1;
                                                   >> 861   
                                                   >> 862   particle_energy = (rndm - q)/m;
                                                   >> 863 
                                                   >> 864   if(verbosityLevel >= 1)
1488     {                                            865     {
1489       nbelow = middle;                        << 866       G4cout << "Energy is " << particle_energy << G4endl;
1490     }                                            867     }
1491   }                                           << 
1492                                               << 
1493   // Now interpolate in that bin to find the  << 
1494   //                                          << 
1495   G4double x1, x2, y1, y2, t, q;              << 
1496   x1 = Bbody_x->at(nbelow);                   << 
1497                                               << 
1498   if(nbelow+1 == static_cast<G4int>(Bbody_x-> << 
1499   {                                           << 
1500     x2 = Bbody_x->back();                     << 
1501   }                                           << 
1502   else                                        << 
1503   {                                           << 
1504     x2 = Bbody_x->at(nbelow + 1);             << 
1505   }                                           << 
1506   y1 = BBHist->at(nbelow);                    << 
1507   if(nbelow+1 == static_cast<G4int>(BBHist->s << 
1508   {                                           << 
1509     G4cout << BBHist->back() << G4endl;       << 
1510     y2 = BBHist->back();                      << 
1511   }                                           << 
1512   else                                        << 
1513   {                                           << 
1514     y2 = BBHist->at(nbelow + 1);              << 
1515   }                                           << 
1516   t = (y2 - y1) / (x2 - x1);                  << 
1517   q = y1 - t * x1;                            << 
1518                                               << 
1519   threadLocalData.Get().particle_energy = (rn << 
1520                                               << 
1521   if (verbosityLevel >= 1)                    << 
1522   {                                           << 
1523     G4cout << "Energy is " << threadLocalData << 
1524   }                                           << 
1525 }                                                868 }
1526                                                  869 
1527 void G4SPSEneDistribution::GenerateCdgEnergie    870 void G4SPSEneDistribution::GenerateCdgEnergies()
1528 {                                                871 {
1529   // Generate random numbers, compare with va << 872   // Gen random numbers, compare with values in cumhist
1530   // to find appropriate part of spectrum and << 873   // to find appropriate part of spectrum and then 
1531   // generate energy in the usual inversion w << 874   // generate energy in the usual inversion way.
1532                                               << 875   //  G4double pfact[2] = {8.5, 112};
                                                   >> 876   // G4double spind[2] = {1.4, 2.3};
                                                   >> 877   // G4double ene_line[3] = {1., 18., 1E6};
1533   G4double rndm, rndm2;                          878   G4double rndm, rndm2;
1534   G4double ene_line[3]={0,0,0};               << 879   G4double ene_line[3];
1535   G4double omalpha[2]={0,0};                  << 880   G4double omalpha[2];
1536   threadLocal_t& params = threadLocalData.Get << 881   if(Emin < 18*keV && Emax < 18*keV)
1537   if (params.Emin < 18 * keV && params.Emax < << 882     {
1538   {                                           << 883       omalpha[0] = 1. - 1.4;
1539     omalpha[0] = 1. - 1.4;                    << 884       ene_line[0] = Emin;
1540     ene_line[0] = params.Emin;                << 885       ene_line[1] = Emax;
1541     ene_line[1] = params.Emax;                << 886     }
1542   }                                           << 887   if(Emin < 18*keV && Emax > 18*keV)
1543   if (params.Emin < 18 * keV && params.Emax > << 888     {
1544   {                                           << 889       omalpha[0] = 1. - 1.4;
1545     omalpha[0] = 1. - 1.4;                    << 890       omalpha[1] = 1. - 2.3;
1546     omalpha[1] = 1. - 2.3;                    << 891       ene_line[0] = Emin;
1547     ene_line[0] = params.Emin;                << 892       ene_line[1] = 18.;
1548     ene_line[1] = 18. * keV;                  << 893       ene_line[2] = Emax;
1549     ene_line[2] = params.Emax;                << 894     }
1550   }                                           << 895   if(Emin > 18*keV)
1551   if (params.Emin > 18 * keV)                 << 896     {
1552   {                                           << 897       omalpha[0] = 1. - 2.3;
1553     omalpha[0] = 1. - 2.3;                    << 898       ene_line[0] = Emin;
1554     ene_line[0] = params.Emin;                << 899       ene_line[1] = Emax;
1555     ene_line[1] = params.Emax;                << 900     }
1556   }                                           << 
1557   rndm = eneRndm->GenRandEnergy();               901   rndm = eneRndm->GenRandEnergy();
1558   rndm2 = eneRndm->GenRandEnergy();              902   rndm2 = eneRndm->GenRandEnergy();
1559                                                  903 
1560   G4int i = 0;                                   904   G4int i = 0;
1561   while (rndm >= CDGhist[i])                  << 905   while( rndm >= CDGhist[i])
1562   {                                           << 906     {
1563     ++i;                                      << 907       i++;
1564   }                                           << 908     }
                                                   >> 909   // Generate final energy.
                                                   >> 910   particle_energy = (std::pow(ene_line[i-1],omalpha[i-1]) + (std::pow(ene_line[i],omalpha[i-1])
                                                   >> 911           - std::pow(ene_line[i-1],omalpha[i-1]))*rndm2);
                                                   >> 912   particle_energy = std::pow(particle_energy,(1./omalpha[i-1]));
1565                                                  913 
1566   // Generate final energy                    << 914   if(verbosityLevel >= 1)
1567   //                                          << 915     G4cout << "Energy is " << particle_energy << G4endl;
1568   G4double ene = (std::pow(ene_line[i - 1], o << 
1569                + (std::pow(ene_line[i], omalp << 
1570                 - std::pow(ene_line[i - 1], o << 
1571   params.particle_energy = std::pow(ene, (1.  << 
1572                                               << 
1573   if (verbosityLevel >= 1)                    << 
1574   {                                           << 
1575     G4cout << "Energy is " << params.particle << 
1576   }                                           << 
1577 }                                                916 }
1578                                                  917 
1579 void G4SPSEneDistribution::GenUserHistEnergie    918 void G4SPSEneDistribution::GenUserHistEnergies()
1580 {                                                919 {
1581   // Histograms are DIFFERENTIAL              << 920   // Histograms are DIFFERENTIAL.
1582                                               << 921   //  G4cout << "In GenUserHistEnergies " << G4endl;
1583   G4AutoLock l(&mutex);                       << 922   if(IPDFEnergyExist == false)
1584                                               << 923     {
1585   if (!IPDFEnergyExist)                       << 924       G4int ii;
1586   {                                           << 925       G4int maxbin = G4int(UDefEnergyH.GetVectorLength());
1587     std::size_t ii;                           << 926       G4double bins[1024], vals[1024], sum;
1588     std::size_t maxbin = UDefEnergyH.GetVecto << 927       sum=0.;
1589     G4double bins[1024], vals[1024], sum;     << 928 
1590     for ( ii = 0 ; ii<1024 ; ++ii ) { bins[ii << 929       if((EnergySpec == false) && (particle_definition == NULL))
1591     sum = 0.;                                 << 930   G4cout << "Error: particle definition is NULL" << G4endl;
1592                                               << 931       
1593     if ( (!EnergySpec)                        << 932       if(maxbin > 1024)
1594       && (threadLocalData.Get().particle_defi << 933   {
1595     {                                         << 934     G4cout << "Maxbin > 1024" << G4endl;
1596       G4Exception("G4SPSEneDistribution::GenU << 935     G4cout << "Setting maxbin to 1024, other bins are lost" << G4endl;
1597                   "Event0302", FatalException << 936   }
1598                   "Error: particle definition << 
1599     }                                         << 
1600                                               << 
1601     if (maxbin > 1024)                        << 
1602     {                                         << 
1603       G4Exception("G4SPSEneDistribution::GenU << 
1604                   "Event0302", JustWarning,   << 
1605                  "Maxbin>1024\n Setting maxbi << 
1606       maxbin = 1024;                          << 
1607     }                                         << 
1608                                               << 
1609     if (!DiffSpec)                            << 
1610     {                                         << 
1611       G4cout << "Histograms are Differential! << 
1612     }                                         << 
1613     else                                      << 
1614     {                                         << 
1615       bins[0] = UDefEnergyH.GetLowEdgeEnergy( << 
1616       vals[0] = UDefEnergyH(0);               << 
1617       sum = vals[0];                          << 
1618       for (ii = 1; ii < maxbin; ++ii)         << 
1619       {                                       << 
1620         bins[ii] = UDefEnergyH.GetLowEdgeEner << 
1621         vals[ii] = UDefEnergyH(ii) + vals[ii  << 
1622         sum = sum + UDefEnergyH(ii);          << 
1623       }                                       << 
1624     }                                         << 
1625                                               << 
1626     if (!EnergySpec)                          << 
1627     {                                         << 
1628       G4double mass = threadLocalData.Get().p << 
1629                                               << 
1630       // Multiply the function (vals) up by t << 
1631       // to make the function counts/s (i.e.  << 
1632                                               << 
1633       for (ii = 1; ii < maxbin; ++ii)         << 
1634       {                                       << 
1635         vals[ii] = vals[ii] * (bins[ii] - bin << 
1636       }                                       << 
1637                                               << 
1638       // Put energy bins into new histo, plus << 
1639       // to make evals counts/s/energy        << 
1640       //                                      << 
1641       for (ii = 0; ii < maxbin; ++ii)         << 
1642       {                                       << 
1643         // kinetic energy                     << 
1644         //                                    << 
1645         bins[ii] = std::sqrt((bins[ii]*bins[i << 
1646       }                                       << 
1647       for (ii = 1; ii < maxbin; ++ii)         << 
1648       {                                       << 
1649         vals[ii] = vals[ii] / (bins[ii] - bin << 
1650       }                                       << 
1651       sum = vals[maxbin - 1];                 << 
1652       vals[0] = 0.;                           << 
1653     }                                         << 
1654     for (ii = 0; ii < maxbin; ++ii)           << 
1655     {                                         << 
1656       vals[ii] = vals[ii] / sum;              << 
1657       IPDFEnergyH.InsertValues(bins[ii], vals << 
1658     }                                         << 
1659                                                  937 
1660     IPDFEnergyExist = true;                   << 938       if(DiffSpec == false)
1661     if (verbosityLevel > 1)                   << 939   G4cout << "Histograms are Differential!!! " << G4endl;
1662     {                                         << 940       else
1663       IPDFEnergyH.DumpValues();               << 941   {
                                                   >> 942     bins[0] = UDefEnergyH.GetLowEdgeEnergy(size_t(0));
                                                   >> 943     vals[0] = UDefEnergyH(size_t(0));
                                                   >> 944     sum = vals[0];
                                                   >> 945     for(ii=1;ii<maxbin;ii++)
                                                   >> 946       {
                                                   >> 947         bins[ii] = UDefEnergyH.GetLowEdgeEnergy(size_t(ii));
                                                   >> 948         vals[ii] = UDefEnergyH(size_t(ii)) + vals[ii-1];
                                                   >> 949         sum = sum + UDefEnergyH(size_t(ii));
                                                   >> 950       }
                                                   >> 951   }
                                                   >> 952 
                                                   >> 953       if(EnergySpec == false)
                                                   >> 954   {
                                                   >> 955     G4double mass = particle_definition->GetPDGMass();
                                                   >> 956     // multiply the function (vals) up by the bin width
                                                   >> 957     // to make the function counts/s (i.e. get rid of momentum
                                                   >> 958     // dependence).
                                                   >> 959     for(ii=1;ii<maxbin;ii++)
                                                   >> 960       {
                                                   >> 961         vals[ii] = vals[ii] * (bins[ii] - bins[ii-1]);
                                                   >> 962       }
                                                   >> 963     // Put energy bins into new histo, plus divide by energy bin width
                                                   >> 964     // to make evals counts/s/energy
                                                   >> 965     for(ii=0;ii<maxbin;ii++)
                                                   >> 966       {
                                                   >> 967         bins[ii] = std::sqrt((bins[ii]*bins[ii]) + (mass*mass)) - mass; //kinetic energy
                                                   >> 968       }
                                                   >> 969     for(ii=1;ii<maxbin;ii++)
                                                   >> 970       {
                                                   >> 971         vals[ii] = vals[ii]/(bins[ii] - bins[ii-1]);
                                                   >> 972       }
                                                   >> 973     sum = vals[maxbin-1];
                                                   >> 974     vals[0] = 0.;
                                                   >> 975   }
                                                   >> 976       for(ii=0;ii<maxbin;ii++)
                                                   >> 977   {
                                                   >> 978     vals[ii] = vals[ii]/sum;
                                                   >> 979     IPDFEnergyH.InsertValues(bins[ii], vals[ii]);
                                                   >> 980   }
                                                   >> 981       
                                                   >> 982       // Make IPDFEnergyExist = true
                                                   >> 983       IPDFEnergyExist = true;
                                                   >> 984       if(verbosityLevel > 1)
                                                   >> 985   IPDFEnergyH.DumpValues();    
1664     }                                            986     }
1665   }                                           << 987   
1666   l.unlock();                                 << 
1667                                               << 
1668   // IPDF has been create so carry on            988   // IPDF has been create so carry on
1669   //                                          << 
1670   G4double rndm = eneRndm->GenRandEnergy();      989   G4double rndm = eneRndm->GenRandEnergy();
1671   threadLocalData.Get().particle_energy= IPDF << 990   particle_energy = IPDFEnergyH.GetEnergy(rndm);
1672                                               << 991   
1673   if (verbosityLevel >= 1)                    << 992   if(verbosityLevel >= 1)
1674   {                                           << 
1675     G4cout << "Energy is " << particle_energy    993     G4cout << "Energy is " << particle_energy << G4endl;
1676   }                                           << 
1677 }                                             << 
1678                                               << 
1679 G4double G4SPSEneDistribution::GetArbEneWeigh << 
1680 {                                             << 
1681   auto nbelow = IPDFArbEnergyH.FindBin(ene,(I << 
1682   G4double wei = 0.;                          << 
1683   if(IntType=="Lin")                          << 
1684   {                                           << 
1685     // note: grad[i] and cept[i] are calculat << 
1686     auto gr = Arb_grad[nbelow + 1];           << 
1687     auto ce = Arb_cept[nbelow + 1];           << 
1688     wei = ene*gr + ce;                        << 
1689   }                                           << 
1690   else if(IntType=="Log")                     << 
1691   {                                           << 
1692     auto alp = Arb_alpha[nbelow + 1];         << 
1693     auto cns = Arb_Const[nbelow + 1];         << 
1694     wei = cns * std::pow(ene,alp);            << 
1695   }                                           << 
1696   else if(IntType=="Exp")                     << 
1697   {                                           << 
1698     auto e0 = Arb_ezero[nbelow + 1];          << 
1699     auto cns = Arb_Const[nbelow + 1];         << 
1700     wei = cns * std::exp(-ene/e0);            << 
1701   }                                           << 
1702   else if(IntType=="Spline")                  << 
1703   {                                           << 
1704     wei = SplineInt[nbelow+1]->CubicSplineInt << 
1705   }                                           << 
1706   return wei;                                 << 
1707 }                                                994 }
1708                                                  995 
1709 void G4SPSEneDistribution::GenArbPointEnergie    996 void G4SPSEneDistribution::GenArbPointEnergies()
1710 {                                                997 {
1711   if (verbosityLevel > 0)                     << 998   if(verbosityLevel > 0)
1712   {                                           << 
1713     G4cout << "In GenArbPointEnergies" << G4e    999     G4cout << "In GenArbPointEnergies" << G4endl;
1714   }                                           << 1000   G4double rndm;
1715                                               << 1001   rndm = eneRndm->GenRandEnergy();
1716   G4double rndm = eneRndm->GenRandEnergy();   << 1002   if(IntType != "Spline")
1717                                               << 
1718   // Find the Bin, have x, y, no of points, a << 
1719   //                                          << 
1720   std::size_t nabove = IPDFArbEnergyH.GetVect << 
1721                                               << 
1722   // Binary search to find bin that rndm is i << 
1723   //                                          << 
1724   while (nabove - nbelow > 1)                 << 
1725   {                                           << 
1726     middle = (nabove + nbelow) / 2;           << 
1727     if (rndm == IPDFArbEnergyH(middle))       << 
1728     {                                         << 
1729       break;                                  << 
1730     }                                         << 
1731     if (rndm < IPDFArbEnergyH(middle))        << 
1732     {                                         << 
1733       nabove = middle;                        << 
1734     }                                         << 
1735     else                                      << 
1736     {                                         << 
1737       nbelow = middle;                        << 
1738     }                                         << 
1739   }                                           << 
1740   threadLocal_t& params = threadLocalData.Get << 
1741   if (IntType == "Lin")                       << 
1742   {                                           << 
1743     // Update thread-local copy of parameters << 
1744     //                                        << 
1745     params.Emax = IPDFArbEnergyH.GetLowEdgeEn << 
1746     params.Emin = IPDFArbEnergyH.GetLowEdgeEn << 
1747     params.grad = Arb_grad[nbelow + 1];       << 
1748     params.cept = Arb_cept[nbelow + 1];       << 
1749     GenerateLinearEnergies(true);             << 
1750   }                                           << 
1751   else if (IntType == "Log")                  << 
1752   {                                           << 
1753     params.Emax = IPDFArbEnergyH.GetLowEdgeEn << 
1754     params.Emin = IPDFArbEnergyH.GetLowEdgeEn << 
1755     params.alpha = Arb_alpha[nbelow + 1];     << 
1756     GeneratePowEnergies(true);                << 
1757   }                                           << 
1758   else if (IntType == "Exp")                  << 
1759   {                                           << 
1760     params.Emax = IPDFArbEnergyH.GetLowEdgeEn << 
1761     params.Emin = IPDFArbEnergyH.GetLowEdgeEn << 
1762     params.Ezero = Arb_ezero[nbelow + 1];     << 
1763     GenerateExpEnergies(true);                << 
1764   }                                           << 
1765   else if (IntType == "Spline")               << 
1766   {                                           << 
1767     params.Emax = IPDFArbEnergyH.GetLowEdgeEn << 
1768     params.Emin = IPDFArbEnergyH.GetLowEdgeEn << 
1769     params.particle_energy = -1e100;          << 
1770     rndm = eneRndm->GenRandEnergy();          << 
1771     while (params.particle_energy < params.Em << 
1772         || params.particle_energy > params.Em << 
1773     {                                         << 
1774       params.particle_energy =                << 
1775         SplineInt[nbelow+1]->CubicSplineInter << 
1776       rndm = eneRndm->GenRandEnergy();        << 
1777     }                                         << 
1778     if (verbosityLevel >= 1)                  << 
1779     {                                            1003     {
1780       G4cout << "Energy is " << params.partic << 1004       //      IPDFArbEnergyH.DumpValues();
                                                   >> 1005       // Find the Bin
                                                   >> 1006       // have x, y, no of points, and cumulative area distribution
                                                   >> 1007       G4int nabove, nbelow = 0, middle;
                                                   >> 1008       nabove = IPDFArbEnergyH.GetVectorLength();
                                                   >> 1009       //      G4cout << nabove << G4endl;
                                                   >> 1010       // Binary search to find bin that rndm is in
                                                   >> 1011       while(nabove-nbelow > 1)
                                                   >> 1012   {
                                                   >> 1013     middle = (nabove + nbelow)/2;
                                                   >> 1014     if(rndm == IPDFArbEnergyH(size_t(middle))) break;
                                                   >> 1015     if(rndm < IPDFArbEnergyH(size_t(middle))) nabove = middle;
                                                   >> 1016     else nbelow = middle;
                                                   >> 1017   }
                                                   >> 1018       if(IntType == "Lin")
                                                   >> 1019   {
                                                   >> 1020     Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow+1));
                                                   >> 1021     Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
                                                   >> 1022     grad = Arb_grad[nbelow+1];
                                                   >> 1023     cept = Arb_cept[nbelow+1];
                                                   >> 1024     //    G4cout << rndm << " " << Emax << " " << Emin << " " << grad << " " << cept << G4endl;
                                                   >> 1025     GenerateLinearEnergies(true);
                                                   >> 1026   }
                                                   >> 1027       else if(IntType == "Log")
                                                   >> 1028   {
                                                   >> 1029     Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow+1));
                                                   >> 1030     Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
                                                   >> 1031     alpha = Arb_alpha[nbelow+1];
                                                   >> 1032     //    G4cout << rndm << " " << Emax << " " << Emin << " " << alpha << G4endl;
                                                   >> 1033     GeneratePowEnergies(true);
                                                   >> 1034   }
                                                   >> 1035       else if(IntType == "Exp")
                                                   >> 1036   {
                                                   >> 1037     Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow+1));
                                                   >> 1038     Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow));
                                                   >> 1039     Ezero = Arb_ezero[nbelow+1];
                                                   >> 1040     //    G4cout << rndm << " " << Emax << " " << Emin << " " << Ezero << G4endl;
                                                   >> 1041     GenerateExpEnergies(true);
                                                   >> 1042   }
                                                   >> 1043     }
                                                   >> 1044   else if(IntType == "Spline")
                                                   >> 1045     {
                                                   >> 1046       if(verbosityLevel > 1)
                                                   >> 1047   G4cout << "IntType = Spline " << rndm << G4endl;
                                                   >> 1048       // in SplineInterpolation created SplineInt
                                                   >> 1049       // Now generate a random number put it into CubicSplineInterpolation
                                                   >> 1050       // and you should get out an energy!?!
                                                   >> 1051       particle_energy = -1e100;
                                                   >> 1052       while (particle_energy < Emin || particle_energy > Emax ) {
                                                   >> 1053   particle_energy = SplineInt->CubicSplineInterpolation(rndm);
                                                   >> 1054   rndm = eneRndm->GenRandEnergy();
                                                   >> 1055       } 
                                                   >> 1056       if(verbosityLevel >= 1)
                                                   >> 1057   G4cout << "Energy is " << particle_energy << G4endl;
1781     }                                            1058     }
1782   }                                           << 
1783   else                                           1059   else
1784   {                                           << 1060     G4cout << "Error: IntType unknown type" << G4endl;
1785     G4Exception("G4SPSEneDistribution::GenArb << 
1786                 FatalException, "Error: IntTy << 
1787   }                                           << 
1788 }                                                1061 }
1789                                                  1062 
1790 void G4SPSEneDistribution::GenEpnHistEnergies    1063 void G4SPSEneDistribution::GenEpnHistEnergies()
1791 {                                                1064 {
1792   // Firstly convert to energy if not already << 1065   //  G4cout << "In GenEpnHistEnergies " << Epnflag << G4endl;
1793                                               << 1066   
1794   G4AutoLock l(&mutex);                       << 1067   // Firstly convert to energy if not already done.
1795                                               << 1068   if(Epnflag == true)
1796   if (Epnflag)  // true means spectrum is epn << 1069     // epnflag = true means spectrum is epn, false means e.
1797   {                                           << 1070     {
1798     // Convert to energy by multiplying by A  << 1071       // convert to energy by multiplying by A number
1799     //                                        << 1072       ConvertEPNToEnergy();
1800     ConvertEPNToEnergy();                     << 1073       // EpnEnergyH will be replace by UDefEnergyH.
1801   }                                           << 1074       //      UDefEnergyH.DumpValues();
1802   if (!IPDFEnergyExist)                       << 1075     }
1803   {                                           << 1076 
1804     // IPDF has not been created, so create i << 1077   //  G4cout << "Creating IPDFEnergy if not already done so" << G4endl;
1805     //                                        << 1078   if(IPDFEnergyExist == false)
1806     G4double bins[1024], vals[1024], sum;     << 1079     {
1807     std::size_t ii;                           << 1080       // IPDF has not been created, so create it
1808     std::size_t maxbin = UDefEnergyH.GetVecto << 1081       G4double bins[1024],vals[1024], sum;
1809     bins[0] = UDefEnergyH.GetLowEdgeEnergy(0) << 1082       G4int ii;
1810     vals[0] = UDefEnergyH(0);                 << 1083       G4int maxbin = G4int(UDefEnergyH.GetVectorLength());
1811     sum = vals[0];                            << 1084       bins[0] = UDefEnergyH.GetLowEdgeEnergy(size_t(0));
1812     for (ii = 1; ii < maxbin; ++ii)           << 1085       vals[0] = UDefEnergyH(size_t(0));
1813     {                                         << 1086       sum = vals[0];
1814       bins[ii] = UDefEnergyH.GetLowEdgeEnergy << 1087       for(ii=1;ii<maxbin;ii++)
1815       vals[ii] = UDefEnergyH(ii) + vals[ii -  << 1088   {
1816       sum = sum + UDefEnergyH(ii);            << 1089     bins[ii] = UDefEnergyH.GetLowEdgeEnergy(size_t(ii));
1817     }                                         << 1090     vals[ii] = UDefEnergyH(size_t(ii)) + vals[ii-1];
1818                                               << 1091     sum = sum + UDefEnergyH(size_t(ii));
1819     l.lock();                                 << 1092   }
1820     for (ii = 0; ii < maxbin; ++ii)           << 1093       
1821     {                                         << 1094       for(ii=0;ii<maxbin;ii++)
1822       vals[ii] = vals[ii] / sum;              << 1095   {
1823       IPDFEnergyH.InsertValues(bins[ii], vals << 1096     vals[ii] = vals[ii]/sum;
                                                   >> 1097     IPDFEnergyH.InsertValues(bins[ii], vals[ii]);
                                                   >> 1098   }
                                                   >> 1099       // Make IPDFEpnExist = true
                                                   >> 1100       IPDFEnergyExist = true;
1824     }                                            1101     }
1825     IPDFEnergyExist = true;                   << 1102   //  IPDFEnergyH.DumpValues();
1826                                               << 
1827   }                                           << 
1828   l.unlock();                                 << 
1829                                               << 
1830   // IPDF has been create so carry on            1103   // IPDF has been create so carry on
1831   //                                          << 
1832   G4double rndm = eneRndm->GenRandEnergy();      1104   G4double rndm = eneRndm->GenRandEnergy();
1833   threadLocalData.Get().particle_energy = IPD << 1105   particle_energy = IPDFEnergyH.GetEnergy(rndm);
1834                                                  1106 
1835   if (verbosityLevel >= 1)                    << 1107   if(verbosityLevel >= 1)
1836   {                                           << 1108     G4cout << "Energy is " << particle_energy << G4endl;
1837     G4cout << "Energy is " << threadLocalData << 
1838   }                                           << 
1839 }                                                1109 }
1840                                                  1110 
1841 void G4SPSEneDistribution::ConvertEPNToEnergy << 1111 void G4SPSEneDistribution::ConvertEPNToEnergy()
1842 {                                                1112 {
1843   // Use this before particle generation to c << 1113   // Use this before particle generation to convert  the
1844   // currently stored histogram from energy/n << 1114   // currently stored histogram from energy/nucleon 
1845                                               << 1115   // to energy.
1846   threadLocal_t& params = threadLocalData.Get << 1116   //  G4cout << "In ConvertEpntoEnergy " << G4endl;
1847   if (params.particle_definition == nullptr)  << 1117   if(particle_definition==NULL)
1848   {                                           << 
1849     G4cout << "Error: particle not defined" <    1118     G4cout << "Error: particle not defined" << G4endl;
1850   }                                           << 
1851   else                                           1119   else
1852   {                                           << 
1853     // Need to multiply histogram by the numb << 
1854     // Baryon Number looks to hold the no. of << 
1855     //                                        << 
1856     G4int Bary = params.particle_definition-> << 
1857                                               << 
1858     // Change values in histogram, Read it ou << 
1859     //                                        << 
1860     std::size_t count, maxcount;              << 
1861     maxcount = EpnEnergyH.GetVectorLength();  << 
1862     G4double ebins[1024], evals[1024];        << 
1863     if (maxcount > 1024)                      << 
1864     {                                         << 
1865       G4Exception("G4SPSEneDistribution::Conv << 
1866                   "gps001", JustWarning,      << 
1867                   "Histogram contains more th << 
1868                    Those above 1024 will be i << 
1869       maxcount = 1024;                        << 
1870     }                                         << 
1871     if (maxcount < 1)                         << 
1872     {                                            1120     {
1873       G4Exception("G4SPSEneDistribution::Conv << 1121       // Need to multiply histogram by the number of nucleons.
1874                  "gps001", FatalException,    << 1122       // Baryon Number looks to hold the no. of nucleons.
1875                  "Histogram contains less tha << 1123       G4int Bary = particle_definition->GetBaryonNumber();
1876       return;                                 << 1124       //      G4cout << "Baryon No. " << Bary << G4endl;
1877     }                                         << 1125       // Change values in histogram, Read it out, delete it, re-create it
1878     for (count = 0; count < maxcount; ++count << 1126       G4int count, maxcount;
1879     {                                         << 1127       maxcount = G4int(EpnEnergyH.GetVectorLength());
1880       // Read out                             << 1128       //      G4cout << maxcount << G4endl;
1881       ebins[count] = EpnEnergyH.GetLowEdgeEne << 1129       G4double ebins[1024],evals[1024];
1882       evals[count] = EpnEnergyH(count);       << 1130       if(maxcount > 1024)
1883     }                                         << 1131   {
1884                                               << 1132     G4cout << "Histogram contains more than 1024 bins!" << G4endl;
1885     // Multiply the channels by the nucleon n << 1133     G4cout << "Those above 1024 will be ignored" << G4endl;
1886     //                                        << 1134     maxcount = 1024;
1887     for (count = 0; count < maxcount; ++count << 1135   }
1888     {                                         << 1136       for(count=0;count<maxcount;count++)
1889       ebins[count] = ebins[count] * Bary;     << 1137   {
1890     }                                         << 1138     // Read out
1891                                               << 1139     ebins[count] = EpnEnergyH.GetLowEdgeEnergy(size_t(count));
1892     // Set Emin and Emax                      << 1140     evals[count] = EpnEnergyH(size_t(count));
1893     //                                        << 1141   }
1894     params.Emin = ebins[0];                   << 1142 
1895     if (maxcount > 1)                         << 1143       // Multiply the channels by the nucleon number to give energies
1896     {                                         << 1144       for(count=0;count<maxcount;count++)
1897       params.Emax = ebins[maxcount - 1];      << 1145   {
1898     }                                         << 1146     ebins[count] = ebins[count] * Bary;
1899     else                                      << 1147   }
1900     {                                         << 1148 
1901       params.Emax = ebins[0];                 << 1149       // Set Emin and Emax
1902     }                                         << 1150       Emin = ebins[0];
1903                                               << 1151       Emax = ebins[maxcount-1];
1904     // Put energy bins into new histogram - U << 1152 
1905     //                                        << 1153       // Put energy bins into new histogram - UDefEnergyH.
1906     for (count = 0; count < maxcount; ++count << 1154       for(count=0;count<maxcount;count++)
1907     {                                         << 1155   {
1908       UDefEnergyH.InsertValues(ebins[count],  << 1156     UDefEnergyH.InsertValues(ebins[count], evals[count]);
1909     }                                         << 1157   }
1910     Epnflag = false; // so that you dont repe << 1158       Epnflag = false; //so that you dont repeat this method.
1911   }                                           << 1159     }  
1912 }                                                1160 }
1913                                                  1161 
1914 void G4SPSEneDistribution::ReSetHist(const G4 << 1162 //
                                                   >> 1163 void G4SPSEneDistribution::ReSetHist(G4String atype)
1915 {                                                1164 {
1916   G4AutoLock l(&mutex);                       << 1165   if (atype == "energy"){
1917   if (atype == "energy")                      << 1166     UDefEnergyH = IPDFEnergyH = ZeroPhysVector ;
1918   {                                           << 1167     IPDFEnergyExist = false ;
1919     UDefEnergyH = IPDFEnergyH = ZeroPhysVecto << 
1920     IPDFEnergyExist = false;                  << 
1921     Emin = 0.;                                   1168     Emin = 0.;
1922     Emax = 1e30;                              << 1169     Emax = 1e30;}  
1923   }                                           << 1170   else if ( atype == "arb"){
1924   else if (atype == "arb")                    << 1171     ArbEnergyH =IPDFArbEnergyH = ZeroPhysVector ;
1925   {                                           << 1172     IPDFArbExist = false;} 
1926     ArbEnergyH = IPDFArbEnergyH = ZeroPhysVec << 1173   else if ( atype == "epn"){
1927     IPDFArbExist = false;                     << 1174     UDefEnergyH = IPDFEnergyH = ZeroPhysVector ;
1928   }                                           << 1175     IPDFEnergyExist = false ;
1929   else if (atype == "epn")                    << 1176     EpnEnergyH = ZeroPhysVector ;}
1930   {                                           << 1177   else {
1931     UDefEnergyH = IPDFEnergyH = ZeroPhysVecto << 
1932     IPDFEnergyExist = false;                  << 
1933     EpnEnergyH = ZeroPhysVector;              << 
1934   }                                           << 
1935   else                                        << 
1936   {                                           << 
1937     G4cout << "Error, histtype not accepted "    1178     G4cout << "Error, histtype not accepted " << G4endl;
1938   }                                              1179   }
1939 }                                                1180 }
1940                                                  1181 
1941 G4double G4SPSEneDistribution::GenerateOne(G4    1182 G4double G4SPSEneDistribution::GenerateOne(G4ParticleDefinition* a)
1942 {                                                1183 {
1943   // Copy global shared status to thread-loca << 1184   particle_definition = a;
1944   //                                          << 1185   particle_energy = -1.;
1945   threadLocal_t& params = threadLocalData.Get << 1186   while ( (EnergyDisType == "Arb")? (particle_energy < ArbEmin || particle_energy > ArbEmax) 
1946   params.particle_definition=a;               << 1187     : (particle_energy < Emin || particle_energy > Emax) ) { 
1947   params.particle_energy=-1;                  << 1188     if(EnergyDisType == "Mono")
1948   if(applyEvergyWeight)                       << 1189       GenerateMonoEnergetic();
1949   {                                           << 1190     else if(EnergyDisType == "Lin")
1950     params.Emax = ArbEmax;                    << 1191       GenerateLinearEnergies();
1951     params.Emin = ArbEmin;                    << 1192     else if(EnergyDisType == "Pow")
1952   }                                           << 1193       GeneratePowEnergies();
1953   else                                        << 1194     else if(EnergyDisType == "Exp")
1954   {                                           << 1195       GenerateExpEnergies();
1955     params.Emax = Emax;                       << 1196     else if(EnergyDisType == "Gauss")
1956     params.Emin = Emin;                       << 1197       GenerateGaussEnergies();
1957   }                                           << 1198     else if(EnergyDisType == "Brem")
1958   params.alpha = alpha;                       << 1199       GenerateBremEnergies();
1959   params.Ezero = Ezero;                       << 1200     else if(EnergyDisType == "Bbody")
1960   params.grad = grad;                         << 1201       GenerateBbodyEnergies();
1961   params.cept = cept;                         << 1202     else if(EnergyDisType == "Cdg")
1962   params.weight = weight;                     << 1203       GenerateCdgEnergies();
1963   // particle_energy = -1.;                   << 1204     else if(EnergyDisType == "User")
1964                                               << 1205       GenUserHistEnergies();
1965   if((EnergyDisType == "Mono") && ((MonoEnerg << 1206     else if(EnergyDisType == "Arb")
1966   {                                           << 1207       GenArbPointEnergies();
1967     G4ExceptionDescription ed;                << 1208     else if(EnergyDisType == "Epn")
1968     ed << "MonoEnergy " << G4BestUnit(MonoEne << 1209       GenEpnHistEnergies();
1969        << " is outside of [Emin,Emax] = ["    << 
1970        << G4BestUnit(Emin,"Energy") << ", "   << 
1971        << G4BestUnit(Emax,"Energy") << ". Mon << 
1972     G4Exception("G4SPSEneDistribution::Genera << 
1973                 "GPS0001", JustWarning, ed);  << 
1974     params.particle_energy=MonoEnergy;        << 
1975     return params.particle_energy;            << 
1976   }                                           << 
1977   while ( (EnergyDisType == "Arb")            << 
1978         ?   (params.particle_energy < ArbEmin << 
1979           || params.particle_energy > ArbEmax << 
1980         :   (params.particle_energy < params. << 
1981           || params.particle_energy > params. << 
1982   {                                           << 
1983     if (Biased)                               << 
1984     {                                         << 
1985       GenerateBiasPowEnergies();              << 
1986     }                                         << 
1987     else                                         1210     else
1988     {                                         << 1211       G4cout << "Error: EnergyDisType has unusual value" << G4endl;
1989       if (EnergyDisType == "Mono")            << 
1990       {                                       << 
1991         GenerateMonoEnergetic();              << 
1992       }                                       << 
1993       else if (EnergyDisType == "Lin")        << 
1994       {                                       << 
1995         GenerateLinearEnergies(false);        << 
1996       }                                       << 
1997       else if (EnergyDisType == "Pow")        << 
1998       {                                       << 
1999         GeneratePowEnergies(false);           << 
2000       }                                       << 
2001       else if (EnergyDisType == "CPow")       << 
2002       {                                       << 
2003         GenerateCPowEnergies();               << 
2004       }                                       << 
2005       else if (EnergyDisType == "Exp")        << 
2006       {                                       << 
2007         GenerateExpEnergies(false);           << 
2008       }                                       << 
2009       else if (EnergyDisType == "Gauss")      << 
2010       {                                       << 
2011         GenerateGaussEnergies();              << 
2012       }                                       << 
2013       else if (EnergyDisType == "Brem")       << 
2014       {                                       << 
2015         GenerateBremEnergies();               << 
2016       }                                       << 
2017       else if (EnergyDisType == "Bbody")      << 
2018       {                                       << 
2019         GenerateBbodyEnergies();              << 
2020       }                                       << 
2021       else if (EnergyDisType == "Cdg")        << 
2022       {                                       << 
2023         GenerateCdgEnergies();                << 
2024       }                                       << 
2025       else if (EnergyDisType == "User")       << 
2026       {                                       << 
2027         GenUserHistEnergies();                << 
2028       }                                       << 
2029       else if (EnergyDisType == "Arb")        << 
2030       {                                       << 
2031         GenArbPointEnergies();                << 
2032       }                                       << 
2033       else if (EnergyDisType == "Epn")        << 
2034       {                                       << 
2035         GenEpnHistEnergies();                 << 
2036       }                                       << 
2037       else                                    << 
2038       {                                       << 
2039         G4cout << "Error: EnergyDisType has u << 
2040       }                                       << 
2041     }                                         << 
2042   }                                              1212   }
2043    return params.particle_energy;             << 1213   return particle_energy;
2044 }                                                1214 }
2045                                                  1215 
2046 G4double G4SPSEneDistribution::GetProbability << 
2047 {                                             << 
2048   G4double prob = 1.;                         << 
2049                                                  1216 
2050   threadLocal_t& params = threadLocalData.Get << 
2051   if (EnergyDisType == "Lin")                 << 
2052   {                                           << 
2053     if (prob_norm == 1.)                      << 
2054     {                                         << 
2055       prob_norm = 0.5*params.grad*params.Emax << 
2056                 + params.cept*params.Emax     << 
2057                 - 0.5*params.grad*params.Emin << 
2058                 - params.cept*params.Emin;    << 
2059     }                                         << 
2060     prob = params.cept + params.grad * ene;   << 
2061     prob /= prob_norm;                        << 
2062   }                                           << 
2063   else if (EnergyDisType == "Pow")            << 
2064   {                                           << 
2065     if (prob_norm == 1.)                      << 
2066     {                                         << 
2067       if (alpha != -1.)                       << 
2068       {                                       << 
2069         G4double emina = std::pow(params.Emin << 
2070         G4double emaxa = std::pow(params.Emax << 
2071         prob_norm = 1./(1.+alpha) * (emaxa -  << 
2072       }                                       << 
2073       else                                    << 
2074       {                                       << 
2075         prob_norm = std::log(params.Emax) - s << 
2076       }                                       << 
2077     }                                         << 
2078     prob = std::pow(ene, params.alpha)/prob_n << 
2079   }                                           << 
2080   else if (EnergyDisType == "Exp")            << 
2081   {                                           << 
2082     if (prob_norm == 1.)                      << 
2083     {                                         << 
2084       prob_norm = -params.Ezero*(std::exp(-pa << 
2085                                - std::exp(par << 
2086     }                                         << 
2087     prob = std::exp(-ene / params.Ezero);     << 
2088     prob /= prob_norm;                        << 
2089   }                                           << 
2090   else if (EnergyDisType == "Arb")            << 
2091   {                                           << 
2092     prob = ArbEnergyH.Value(ene);             << 
2093                                               << 
2094     if (prob <= 0.)                           << 
2095     {                                         << 
2096       G4cout << " Warning:G4SPSEneDistributio << 
2097              << prob << " " << ene << G4endl; << 
2098       prob = 1e-30;                           << 
2099     }                                         << 
2100   }                                           << 
2101   else                                        << 
2102   {                                           << 
2103     G4cout << "Error: EnergyDisType not suppo << 
2104   }                                           << 
2105                                                  1217 
2106   return prob;                                << 1218 
2107 }                                             << 1219 
                                                   >> 1220 
                                                   >> 1221 
                                                   >> 1222 
                                                   >> 1223 
2108                                                  1224