source: fact/tools/rootmacros/PulseTemplates/Sample.C@ 16858

#include "Sample.h"

Sample::Sample()
{
    InitMembers();
    mpRndNumList = new vector<int>;
}

Sample::Sample(int size)
{
    InitMembers();
    mpRndNumList = new vector<int>;
    mpRndNumList->resize(size);
}

Sample::Sample(vector<int> *RndNumList)
{
    InitMembers();
    mpRndNumList    =   RndNumList;
}

Sample::Sample(vector<int> *RndNumList, int size)
{
    InitMembers();
    mpRndNumList    =   RndNumList;
    mpRndNumList->resize(size);
}

Sample::~Sample(void)
{
    if (mpRndNumList != NULL)
    {
        delete mpRndNumList;
    }
}

void Sample::InitMembers()
{
    mpRndNumList    = NULL;
    mMinNumber      = 0;
    mMaxNumber      = 0;
    mSampleSize     = 0;
    mVerbosityLvl   = 0;
    mSeed           = 4357;
}

//============================================================================
//ACCESS
//============================================================================
void    Sample::SetMinNumber(int min){  mMinNumber  = min; return;}
void    Sample::SetMaxNumber(int max){  mMaxNumber  = max; return;}
void    Sample::SetSampleSize(int size){mSampleSize = size; return;}
void    Sample::SetSeed(int seed)
{
    mSeed = seed;
    mRandom.SetSeed(mSeed);
    return;
}
void    Sample::SetVerbosityLvl(int VerbosityLvl){ mVerbosityLvl = VerbosityLvl; return;}
int     Sample::GetMinNumber(){return   mMinNumber;}
int     Sample::GetMaxNumber(){return   mMaxNumber;}
int     Sample::GetSampleSize(){return  mSampleSize;}
int     Sample::GetSeed(){return  mSeed;}

void    Sample::GetRndListValues(vector<int> *RndNumList)
{
    vector<int>::iterator it;

    for ( it=mpRndNumList->begin() ; it < mpRndNumList->end(); it++ )
    {
        RndNumList->push_back(*it);
    }

    return;
}

vector<int>*    Sample::GetRndListPtr()
{
    return mpRndNumList;
}

//============================================================================
//OPERATIONS
//============================================================================
int
Sample::GenerateRandomInt(
        int     min,
        int     max)
{
//    int rndNum = min + mRandom.Integer(max - min);
    int rndNum = (int) mRandom.Uniform(min, max);
//    cout << rndNum << endl;
    return rndNum;
}

int
Sample::GenerateRandomInt()
{
    return GenerateRandomInt( mMinNumber, mMaxNumber );
//    return 0;
}

void
Sample::GenerateSample()
{
    GenerateSample( mpRndNumList, mMinNumber, mMaxNumber, mSampleSize);
    return;
}

void
Sample::GenerateSample(
        vector<int>*    sampleVector,
        int     min,
        int     max,
        int     size)
{
    if ( size > max - min)
    {
        cout << "sample size is larger than range of sample, will set size to range" << endl;
        size = max - min + 1;
    }
    //resize destination vector for pulled numbers
    sampleVector->resize(size);

    //calculate qunatity of numbers in range
    int qunatityOfNumbers = max - min + 1;

    //crate a vector list of ordered numbers in defined range
    vector<int> listOfNumbers;
    listOfNumbers.resize(qunatityOfNumbers);

    //fill a list of ordered numbers in defined range
    for (int i = min; i <= max; i++)
    {
            listOfNumbers.at(i)=i;
    }

    //container to fill in numbers with ordering
    set<int> sampleSet;

    //container for insert result
    bool result;
    for (int i = 0; i < size; i++)
    {
        int randomNumber = GenerateRandomInt( 0, qunatityOfNumbers-i);
        result = sampleSet.insert(listOfNumbers.at(randomNumber)).second;
        if (result)
        {
//        cout << "rndNR " << randomNumber << endl;
            listOfNumbers.erase(listOfNumbers.begin()+randomNumber);
        }
        else if (!result)
        {
            cout << " pulled number exists, pulling again" << endl;
            i--;
        }
    }

    set<int>::iterator it;

    int counter = 0;
    for ( it=sampleSet.begin() ; it != sampleSet.end(); it++ )
    {
        sampleVector->at(counter)=*it;
        counter++;
    }

     return;
}


int
Sample::BootstrapVector(vector<double> *inVector, vector<double> *outVector)
{
    //get size of sample to be bootstrapped
    int sampleSize = inVector->size();

    //Vector with positions in original sample
    vector<int> entryID (sampleSize,0);

    //calculate wich entries from inVector will be put into outVector
    BootstrapSample(&entryID, 0, sampleSize, sampleSize );

    vector<int>::iterator it;

    //Loop over entryID vector to fill content from inVector to outVector
    int counter = 0;
    for ( it=entryID.begin() ; it != entryID.end(); it++ )
    {
        outVector->at(counter) = inVector->at(*it);
        counter++;
    }

    return counter + 1;
}

int
Sample::BootstrapTH1(TH1* inputHisto, TH1*  outHisto)
{
    //reset outHisto in case it is reused
    outHisto->Reset();

    //compute number of bins for 1-d histogram h1
    int nbins = inputHisto->GetXaxis()->GetNbins();

    if (mVerbosityLvl > 1) cout << "nbins: " << nbins << endl;
    //we need to get the binning

    //vector to contain entries of TH1
    vector<double>    entries;

    //quantity of entries in bin
    int     quantity    = 0;        //number of entries in a bin
    double  value       = 0;        //value of a bin

    //Loop over bins to fill entries vektor
    for (int i=1;i<nbins;i++)
    {
        //value       = inputHisto->GetBinLowEdge(i);
        value       = inputHisto->GetBinCenter(i);
        quantity    = inputHisto->GetBinContent(i);
        //Loop over bin quantities
        for (int j = 0; j < quantity; j++)
        {
            // fill entries vektor with value of a bin
            // as many times as the bin has entries
            entries.push_back(value);
        }
    }

    //get size of sample to be bootstrapped
    int sampleSize = entries.size();

    if (mVerbosityLvl > 1) cout << "sampleSize: " << sampleSize << endl;

    //Vector with positions in original sample
    vector<int> entryID (sampleSize,0);

    //calculate a list with random int numbers between 0 ad sampleSize
    // and fill it into entryID vector
    BootstrapSample(&entryID, 0, sampleSize, sampleSize );

    //Loop over entryID vector to bootstrap the histogram
    int counter = 0;
    for ( unsigned int i = 0 ; i < entryID.size(); i++ )
    {
        //fill values pulled from entries vektor into outHisto
        outHisto->Fill(
                    entries.at( entryID.at(i) )
                    );
        counter++;
    }

    //return the nummber of filled entries
    return counter + 1;
}

void
Sample::BootstrapSample()
{
    BootstrapSample( mpRndNumList, mMinNumber, mMaxNumber, mSampleSize );

    return;
}

void
Sample::BootstrapSample(
        vector<int> *sampleVector,
        int   numMinEvent,
        int   numMaxEvent,
        int   size)
{
    if (size == 0){
        if (mVerbosityLvl > 1) cout << "Bootstrapping: size of vector = 0; nothing to do" << endl;
        return;
    }

    //resize the sample vector to size of boostrapped sample
    sampleVector->resize(size);

    //list of rndnumbers pulled with putting back
    multiset<int> sampleSet;

    //loop over samplesize to generate random numbers to fill into sampleset
    for (int i = 0; i < size; i++)
    {
        int randomNumber = GenerateRandomInt( numMinEvent, numMaxEvent);
        sampleSet.insert(randomNumber);
    }

    multiset<int>::iterator it;

//    set<int>::iterator it;

    int counter = 0;
    // loop over list of rndnumbers and fill their entries into sampleVector
    // entries are vector positions
    for ( it=sampleSet.begin() ; it != sampleSet.end(); it++ )
    {
        sampleVector->at(counter)=*it;
        counter++;
    }

    return;
}
//NOTE: crashes for some reason if size is smaller than max-min