ROOT logo
//_HADES_CLASS_DESCRIPTION
////////////////////////////////////////////////////////////////////////////////
// Calculate running mean of a distribution of integer values e.g. multiplicities.
// calculate the sigma (RMS) of the distribution.
// A cutoff at a maximum value is defined as mean value plus sigma * a scaling factor.
// A minimum value can be specified as well (default = 1).
// If the mean value gets below 1, the cutoff value is calculated by (minimum value + 1) + scaling factor*(sigma=1)
// instead of mean value + sigma* scaling factor.
// This allows for recovery from large sets with input values of 0 (e.g. sparks).
// It accounts as well for the fact that input values are integer (0 or 1 or 2 ...), thus cutoff values must be > 1.
// The array size used for averaging can be specified in initParam (default = 2000).
// The values of the distribution have to be >= 0
// add an offset if the original distibution extends to negative values
// W. Koenig April 2015
////////////////////////////////////////////////////////////////////////////////

#include "hparticlerunningmean.h"
#include "TMath.h"
#include "TRandom.h"
#include <iostream>
using namespace std;

ClassImp(HParticleRunningMeanI)

Float_t HParticleRunningMeanI::calcMean(Int_t val)
{
  // calculate new mean and new sigma by replacing old  point by actual one.

    if(val < fmin) {
	fvalid = kFALSE;
	return fmean;
    }

    if(val <= fmax || fn < fminEvts) {
	fmax     = Int_t(max(fmean +fscaleFacSigma*fsigma + 0.5F,fmaxMin));
	fvalid   = kTRUE;
	fSum     += val - fnPoints[findex];
 	fSum2    += val*val - fnPoints[findex]*fnPoints[findex];
        fnPoints[findex] = val;
	if(++findex >= fnMax) findex = 0;
	fmean = Float_t(fSum)/Float_t(fn);
	fsigma = sqrt(float(fSum2)/float(fn) - fmean*fmean);
	if(fn < fnMax) ++fn;
    } else {
	fvalid = kFALSE;
    }
    return fmean;
};

Int_t HParticleRunningMeanI::initParam(const Int_t maxEvents, const Int_t minEvents,
				       const Float_t scaleFacSigma, const Float_t initMean, const Float_t initSigma, const Int_t min)
{
    // Reset moving average, clear multiplicity array,
    // Set minimum number of events to generate a valid multiplicity mean value
    // Set default scaling factor (mean + factor * sigma) to maximum multiplicity (in order to reject monsters from calculating mean)


    fnMax = maxEvents; // number of selected events which are used to calculate mean multiplicity
    if(fnPoints != NULL) delete[] fnPoints;
    fnPoints = new Int_t[fnMax];
    reset(kFALSE);
    if(minEvents <= fnMax) fminEvts = minEvents; //minimum events for valid mean value
    else                   fminEvts = fnMax;
    fscaleFacSigma         = scaleFacSigma; // scaling from mean multiplicity to max multiplicity (excluding outliers)
    fmin                   = max(min,0);
    fmaxMin                = (fmin + 1.0F) + fscaleFacSigma + 0.5F;
    fvalid                 = kFALSE;
    fmax                   = fmaxMin;
    initialMean            = initMean;
    initialSigma           = initSigma;
    setInitialMean(initialMean, initialSigma) ;
    return fnMax; // returns array size used for running mean
};

Int_t HParticleRunningMeanI::setInitialMean( const Float_t initMean, const Float_t initSigma)
{
    // create a set of initial mean values filling first 'fminEvts'
    // multiplicity array members. returns number of initialized multiplicity
    // array members (= fminEvts set in 'initParam' before
    if(initMean==-999999 || initSigma==-999999) return 0;
    initialMean  = initMean;
    initialSigma = initSigma;
    Float_t initialWidth = initSigma*sqrt(12.0F); // width of a box-like distribution
    if(initialWidth < 1.0F) initialWidth = 1.0F; // minimum needed to get correct mean for integers
    for (findex = 0; findex < fminEvts; ++findex) {
	fnPoints[findex] = Int_t(initMean + (gRandom->Uniform()-0.5F)*initialWidth + 0.5F);
	fSum += fnPoints[findex];
	fSum2 += fnPoints[findex]*fnPoints[findex]; // quadratic sum needed for sigma distribution
    };
    fn = fminEvts;
    fmean = Float_t(fSum)/Float_t(fn);
    fsigma=  sqrt(Float_t(fSum2)/Float_t(fn) - fmean*fmean);
    fmax = Int_t (max(fmean + fscaleFacSigma * initSigma + 0.5F,fmaxMin));

    return fn;
};

void HParticleRunningMeanI::reset(Bool_t full)
{
    //clear multiplicity sum, multiplicity array, reset findex counter
    for (Int_t i = 0; i < fnMax; ++i) fnPoints[i] = 0;
    fSum    = 0;
    fSum2   = 0;
    fn      = 1; //starts with 1, runs up while multiplicityPoints array is filled and stays constant after fnMax is reached
    findex  = 0; //running findex pointing to the last entry in the list of multiplicity points
    fmean   = 0.0F;
    fsigma  = 0.0F;
    if(full)setInitialMean(initialMean,initialSigma) ;
};

Bool_t HParticleRunningMeanI::isValid(void)
{
    // returns true if event no > fminEvts (reliable mean) and
    // last val < fmax (no outlier)
    return (fn >= fminEvts) && fvalid;
};
 hparticlerunningmean.cc:1
 hparticlerunningmean.cc:2
 hparticlerunningmean.cc:3
 hparticlerunningmean.cc:4
 hparticlerunningmean.cc:5
 hparticlerunningmean.cc:6
 hparticlerunningmean.cc:7
 hparticlerunningmean.cc:8
 hparticlerunningmean.cc:9
 hparticlerunningmean.cc:10
 hparticlerunningmean.cc:11
 hparticlerunningmean.cc:12
 hparticlerunningmean.cc:13
 hparticlerunningmean.cc:14
 hparticlerunningmean.cc:15
 hparticlerunningmean.cc:16
 hparticlerunningmean.cc:17
 hparticlerunningmean.cc:18
 hparticlerunningmean.cc:19
 hparticlerunningmean.cc:20
 hparticlerunningmean.cc:21
 hparticlerunningmean.cc:22
 hparticlerunningmean.cc:23
 hparticlerunningmean.cc:24
 hparticlerunningmean.cc:25
 hparticlerunningmean.cc:26
 hparticlerunningmean.cc:27
 hparticlerunningmean.cc:28
 hparticlerunningmean.cc:29
 hparticlerunningmean.cc:30
 hparticlerunningmean.cc:31
 hparticlerunningmean.cc:32
 hparticlerunningmean.cc:33
 hparticlerunningmean.cc:34
 hparticlerunningmean.cc:35
 hparticlerunningmean.cc:36
 hparticlerunningmean.cc:37
 hparticlerunningmean.cc:38
 hparticlerunningmean.cc:39
 hparticlerunningmean.cc:40
 hparticlerunningmean.cc:41
 hparticlerunningmean.cc:42
 hparticlerunningmean.cc:43
 hparticlerunningmean.cc:44
 hparticlerunningmean.cc:45
 hparticlerunningmean.cc:46
 hparticlerunningmean.cc:47
 hparticlerunningmean.cc:48
 hparticlerunningmean.cc:49
 hparticlerunningmean.cc:50
 hparticlerunningmean.cc:51
 hparticlerunningmean.cc:52
 hparticlerunningmean.cc:53
 hparticlerunningmean.cc:54
 hparticlerunningmean.cc:55
 hparticlerunningmean.cc:56
 hparticlerunningmean.cc:57
 hparticlerunningmean.cc:58
 hparticlerunningmean.cc:59
 hparticlerunningmean.cc:60
 hparticlerunningmean.cc:61
 hparticlerunningmean.cc:62
 hparticlerunningmean.cc:63
 hparticlerunningmean.cc:64
 hparticlerunningmean.cc:65
 hparticlerunningmean.cc:66
 hparticlerunningmean.cc:67
 hparticlerunningmean.cc:68
 hparticlerunningmean.cc:69
 hparticlerunningmean.cc:70
 hparticlerunningmean.cc:71
 hparticlerunningmean.cc:72
 hparticlerunningmean.cc:73
 hparticlerunningmean.cc:74
 hparticlerunningmean.cc:75
 hparticlerunningmean.cc:76
 hparticlerunningmean.cc:77
 hparticlerunningmean.cc:78
 hparticlerunningmean.cc:79
 hparticlerunningmean.cc:80
 hparticlerunningmean.cc:81
 hparticlerunningmean.cc:82
 hparticlerunningmean.cc:83
 hparticlerunningmean.cc:84
 hparticlerunningmean.cc:85
 hparticlerunningmean.cc:86
 hparticlerunningmean.cc:87
 hparticlerunningmean.cc:88
 hparticlerunningmean.cc:89
 hparticlerunningmean.cc:90
 hparticlerunningmean.cc:91
 hparticlerunningmean.cc:92
 hparticlerunningmean.cc:93
 hparticlerunningmean.cc:94
 hparticlerunningmean.cc:95
 hparticlerunningmean.cc:96
 hparticlerunningmean.cc:97
 hparticlerunningmean.cc:98
 hparticlerunningmean.cc:99
 hparticlerunningmean.cc:100
 hparticlerunningmean.cc:101
 hparticlerunningmean.cc:102
 hparticlerunningmean.cc:103
 hparticlerunningmean.cc:104
 hparticlerunningmean.cc:105
 hparticlerunningmean.cc:106
 hparticlerunningmean.cc:107
 hparticlerunningmean.cc:108
 hparticlerunningmean.cc:109
 hparticlerunningmean.cc:110
 hparticlerunningmean.cc:111
 hparticlerunningmean.cc:112
 hparticlerunningmean.cc:113
 hparticlerunningmean.cc:114
 hparticlerunningmean.cc:115
 hparticlerunningmean.cc:116