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

Last change on this file since 1241 was 1196, checked in by garnier, 16 years ago

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