3 // Nikolai Amelin, Ludmila Malinina, Timur Pocheptsov (C) JINR/Dubna
4 // amelin@sunhe.jinr.ru, malinina@sunhe.jinr.ru, pocheptsov@sunhe.jinr.ru
8 // This class is taken from the GEANT4 tool kit and changed!!!!!
11 #include "RandArrayFunction.h"
14 RandArrayFunction::RandArrayFunction(const Double_t *aProbFunc, Int_t theProbSize, Int_t intType):
18 fInterpolationType(intType)
20 PrepareTable(aProbFunc);
23 RandArrayFunction::RandArrayFunction(Int_t theProbSize, Int_t intType):
27 fInterpolationType(intType)
30 void RandArrayFunction::PrepareTable(const Double_t* aProbFunc) {
31 //Prepares fIntegralPdf.
33 Error("RandArrayFunction::PrepareTable",
34 "RandArrayFunction constructed with no bins - will use flat distribution.");
35 UseFlatDistribution();
39 fIntegralPdf.resize(fNBins + 1);
42 for (ptn = 0; ptn < fNBins; ++ptn ) {
43 Double_t weight = aProbFunc[ptn];
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",
52 fIntegralPdf[ptn + 1] = fIntegralPdf[ptn] + weight;
55 if (fIntegralPdf[fNBins] <= 0.) {
56 Warning("RandArrayFunction::PrepareTable",
57 "RandArrayFunction constructed with nothing in bins - will use flat distribution");
58 UseFlatDistribution();
62 for (ptn = 0; ptn < fNBins + 1; ++ptn)
63 fIntegralPdf[ptn] /= fIntegralPdf[fNBins];
65 // And another useful variable is ...
66 fOneOverNbins = 1.0 / fNBins;
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;
76 void RandArrayFunction::UseFlatDistribution() {
77 //Called only by PrepareTable in case of user error.
79 fIntegralPdf.resize(2);
85 Double_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.
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
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:
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
104 if (fInterpolationType == 1) {
105 return nBelow * fOneOverNbins;
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.
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;
119 Double_t binFraction = (rand - fIntegralPdf[nBelow]) / binMeasure;
121 return (nBelow + binFraction) * fOneOverNbins;
125 void RandArrayFunction::FireArray(Int_t size, Double_t *vect) const {
126 for (Int_t i = 0; i < size; ++i)