TUHKMgen/UHKM/RandArrayFunction.cxx

   1 #include <TError.h>
   2
   3 #include "RandArrayFunction.h"
   4 /*
   5
   6         Nikolai Amelin, Ludmila Malinina, Timur Pocheptsov (C) JINR/Dubna
   7       amelin@sunhe.jinr.ru, malinina@sunhe.jinr.ru, pocheptsov@sunhe.jinr.ru
   8                            November. 2, 2005
   9
  10 */
  11 //This class is taken from the GEANT4 tool kit  and changed!!!!!
  12
  13
  14 RandArrayFunction::RandArrayFunction(const Double_t *aProbFunc, Int_t theProbSize, Int_t intType):
  15   fIntegralPdf(),
  16   fNBins(theProbSize),
  17   fOneOverNbins(0),
  18   fInterpolationType(intType)
  19 {
  20   PrepareTable(aProbFunc);
  21 }
  22
  23 RandArrayFunction::RandArrayFunction(Int_t theProbSize, Int_t intType):
  24   fIntegralPdf(),
  25   fNBins(theProbSize),
  26   fOneOverNbins(0),
  27   fInterpolationType(intType)
  28 {}
  29
  30 void 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
  76 void 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
  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.
  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
 125 void RandArrayFunction::FireArray(Int_t size, Double_t *vect) const {
 126   for (Int_t i = 0; i < size; ++i)
 127     vect[i] = Fire();
 128 }