Coding conventions (Ionut)
[u/mrichter/AliRoot.git] / TUHKMgen / UHKM / RandArrayFunction.cxx
CommitLineData
03896fc4 1//
2//
3// Nikolai Amelin, Ludmila Malinina, Timur Pocheptsov (C) JINR/Dubna
4// amelin@sunhe.jinr.ru, malinina@sunhe.jinr.ru, pocheptsov@sunhe.jinr.ru
5// November. 2, 2005
6//
7//
8// This class is taken from the GEANT4 tool kit and changed!!!!!
b1c2e580 9
03896fc4 10#include <TError.h>
b1c2e580 11#include "RandArrayFunction.h"
b1c2e580 12
13
786056a2 14RandArrayFunction::RandArrayFunction(const Double_t *aProbFunc, Int_t theProbSize, Int_t intType):
15 fIntegralPdf(),
16 fNBins(theProbSize),
17 fOneOverNbins(0),
18 fInterpolationType(intType)
b1c2e580 19{
20 PrepareTable(aProbFunc);
21}
22
786056a2 23RandArrayFunction::RandArrayFunction(Int_t theProbSize, Int_t intType):
24 fIntegralPdf(),
25 fNBins(theProbSize),
26 fOneOverNbins(0),
27 fInterpolationType(intType)
b1c2e580 28{}
29
30void RandArrayFunction::PrepareTable(const Double_t* aProbFunc) {
31 //Prepares fIntegralPdf.
32 if(fNBins < 1) {
33 Error("RandArrayFunction::PrepareTable",
34 "RandArrayFunction constructed with no bins - will use flat distribution.");
35 UseFlatDistribution();
36 return;
37 }
38
39 fIntegralPdf.resize(fNBins + 1);
40 fIntegralPdf[0] = 0;
41 Int_t ptn;
42 for (ptn = 0; ptn < fNBins; ++ptn ) {
43 Double_t weight = aProbFunc[ptn];
44 if (weight < 0.) {
45 // We can't stomach negative bin contents, they invalidate the
46 // search algorithm when the distribution is fired.
47 Warning("RandArrayFunction::PrepareTable",
48 "RandArrayFunction constructed with negative-weight bin %d == %f -- will substitute 0 weight",
49 ptn, weight);
50 weight = 0.;
51 }
52 fIntegralPdf[ptn + 1] = fIntegralPdf[ptn] + weight;
53 }
54
55 if (fIntegralPdf[fNBins] <= 0.) {
56 Warning("RandArrayFunction::PrepareTable",
57 "RandArrayFunction constructed with nothing in bins - will use flat distribution");
58 UseFlatDistribution();
59 return;
60 }
61
62 for (ptn = 0; ptn < fNBins + 1; ++ptn)
63 fIntegralPdf[ptn] /= fIntegralPdf[fNBins];
64
65 // And another useful variable is ...
66 fOneOverNbins = 1.0 / fNBins;
67 // One last chore:
68 if (fInterpolationType && fInterpolationType != 1) {
69 Info("RandArrayFunction::PrepareTable",
70 "RandArrayFunction does not recognize fInterpolationType %d \n"
71 "Will use type 0 (continuous linear interpolation)", fInterpolationType);
72 fInterpolationType = 0;
73 }
74}
75
76void RandArrayFunction::UseFlatDistribution() {
77 //Called only by PrepareTable in case of user error.
78 fNBins = 1;
79 fIntegralPdf.resize(2);
80 fIntegralPdf[0] = 0;
81 fIntegralPdf[1] = 1;
82 fOneOverNbins = 1.0;
83}
84
85Double_t RandArrayFunction::MapRandom(Double_t rand) const {
86 // Private method to take the random (however it is created) and map it
87 // according to the distribution.
88
89 Int_t nBelow = 0; // largest k such that I[k] is known to be <= rand
90 Int_t nAbove = fNBins; // largest k such that I[k] is known to be > rand
91 Int_t middle;
92
93 while (nAbove > nBelow+1) {
94 middle = (nAbove + nBelow+1)>>1;
95 rand >= fIntegralPdf[middle] ? nBelow = middle : nAbove = middle;
96 }// after this loop, nAbove is always nBelow+1 and they straddle rad:
97
98 /*assert ( nAbove = nBelow+1 );
99 assert ( fIntegralPdf[nBelow] <= rand );
100 assert ( fIntegralPdf[nAbove] >= rand );*/
101 // If a defective engine produces rand=1, that will
102 // still give sensible results so we relax the > rand assertion
103
104 if (fInterpolationType == 1) {
105 return nBelow * fOneOverNbins;
106 }
107 else {
108 Double_t binMeasure = fIntegralPdf[nAbove] - fIntegralPdf[nBelow];
109 // binMeasure is always aProbFunc[nBelow],
110 // but we don't have aProbFunc any more so we subtract.
111
112 if (!binMeasure) {
113 // rand lies right in a bin of measure 0. Simply return the center
114 // of the range of that bin. (Any value between k/N and (k+1)/N is
115 // equally good, in this rare case.)
116 return (nBelow + .5) * fOneOverNbins;
117 }
118
119 Double_t binFraction = (rand - fIntegralPdf[nBelow]) / binMeasure;
120
121 return (nBelow + binFraction) * fOneOverNbins;
122 }
123}
124
125void RandArrayFunction::FireArray(Int_t size, Double_t *vect) const {
126 for (Int_t i = 0; i < size; ++i)
127 vect[i] = Fire();
128}