HistogramBase.cpp

Go to the documentation of this file.

#include "Mint/HistogramBase.h"

#include "Mint/StatisticsFinder.h"


HistogramBase::HistogramBase(int nBins) :
  _nBins(nBins),
  _binContents(nBins+1,0.0),
  _sumW2      (nBins+1,0.0),
  _min(-999.9),
  _max(-999.9),
  _minDensity(-999.9),
  _maxDensity(-999.9)
{
  WELCOME_LOG << "Good day from the HistogramBase() Constructor"; 
}

void HistogramBase::resetBinContents(int nBins){
  _nBins = nBins;
  _binContents = std::vector<double>(nBins+1,0.0);
  _sumW2       = std::vector<double>(nBins+1,0.0);
}

void HistogramBase::clear(){
  for (unsigned i = 0; i < _binContents.size(); i++){
    _binContents.at(i) = 0.0;
    _sumW2      .at(i) = 0.0;
  }
}



void HistogramBase::reserveCapacity(int nElements){
  _binContents.reserve(nElements+1);
  _sumW2      .reserve(nElements+1);
}


void HistogramBase::merge( const HistogramBase& other ){
  
  double overflowCont = _binContents.at(_nBins);
  overflowCont += other._binContents.at(other._nBins);
  
  double overflowSumW2 = _sumW2.at(_nBins);
  overflowSumW2 += other._sumW2.at(other._nBins);  
  
  int capacityNeeded = _nBins + other._nBins + 1;

  reserveCapacity(capacityNeeded);

  _binContents.at(_nBins) = other._binContents.at(0);
  _sumW2      .at(_nBins) = other._sumW2      .at(0);
  
  for (int i = 1; i < other._nBins; i++){
    _binContents.push_back( other._binContents.at(i) );
    _sumW2      .push_back( other._sumW2      .at(i) );    
  }

  _binContents.push_back( overflowCont  );
  _sumW2      .push_back( overflowSumW2 );

  _nBins += other._nBins;

}



void HistogramBase::saveBase(){

  TTree tree("HistogramBase", "HistogramBase");

  int    binNumber  = -1;
  double binContent = 0.0;
  double sumW2      = 0.0;

  tree.Branch("binNumber" , &binNumber );
  tree.Branch("binContent", &binContent);
  tree.Branch("sumW2"     , &sumW2     );

  for(unsigned int bin = 0; bin < _binContents.size(); bin++ ){
    binNumber  = bin;
    binContent = _binContents.at(bin);
    sumW2      = _sumW2      .at(bin);
    tree.Fill();
  }
  
  tree.Write();

}

void HistogramBase::saveBase(TString filename){

  TFile file(filename, "RECREATE");

  saveBase();

  file.Write();
  file.Close();

}

void HistogramBase::loadBase(TString filename){

  TFile file(filename, "READ");
  TTree* tree = (TTree*)file.Get("HistogramBase");

  int    binNumber  = -1;
  double binContent = 0.0;
  double sumW2      = 0.0;
  
  int nEntries = tree->GetEntries();

  tree->SetBranchAddress("binNumber" , &binNumber );
  tree->SetBranchAddress("binContent", &binContent);
  tree->SetBranchAddress("sumW2"     , &sumW2     );
  
  this->resetBinContents(nEntries - 1);


  for(int ent = 0; ent < nEntries; ent++){
    tree->GetEntry(ent);
    
    VERBOSE_LOG << "Bin = " << binNumber << "       Content = " << binContent << "        SumW2 = " << sumW2;

    _binContents.at((int)binNumber) = binContent;
    _sumW2      .at((int)binNumber) = sumW2;
  }
  
  file.Close();

}

HistogramBase::~HistogramBase(){
  GOODBYE_LOG << "Goodbye from the HistogramBase() Constructor"; 
}

void HistogramBase::fillBase(int binNum, double weight){
  binNum = checkBinNumber(binNum);
  _binContents[binNum] += weight;
  _sumW2[binNum]       += weight*weight;
}

int HistogramBase::checkBinNumber(int bin) const{
  if(bin == -1){
    //This is also a valid way to add to the underflow/overflow bin with no error message
    return _nBins;
  }
  if(bin > _nBins || bin < 0) {
    ERROR_LOG << "Bin " << bin << " does not exist! Adding to underflow/overflow bin " << _nBins;
    return _nBins;
  }
  return bin;
}

double HistogramBase::getMin() const{
  if (_min != -999.9) return _min;
  MinMaxFinder stats;
  for (int i = 0; i < _nBins; i++){
    stats.add( _binContents[i] );
  }
  return stats.getMin();
}

double HistogramBase::getMax() const{
  if (_max != -999.9) return _max;
  MinMaxFinder stats;
  for (int i = 0; i < _nBins; i++){
    stats.add( _binContents[i] );
  }
  return stats.getMax();
}

double HistogramBase::getMinDensity() const{
  if (_minDensity != -999.9) return _minDensity;
  MinMaxFinder stats;
  for (int i = 0; i < _nBins; i++){
    stats.add( getFrequencyDensity(i) );
  }
  return stats.getMin();
}

double HistogramBase::getMaxDensity() const{
  if (_maxDensity != -999.9) return _maxDensity;
  MinMaxFinder stats;
  for (int i = 0; i < _nBins; i++){
    stats.add( getFrequencyDensity(i) );
  }
  return stats.getMax();
}

void HistogramBase::divide(const HistogramBase& other){
 
  if (other._binContents.size() != _binContents.size()){
    ERROR_LOG << "Trying to divide histograms with different numbers of bins. Doing nothing.";
    return;
  }

  for(unsigned int i = 0; i < _binContents.size(); i++){

    double binCont1 =       _binContents.at(i);
    double binCont2 = other._binContents.at(i);
    double var1     =       _sumW2.at(i);
    double var2     = other._sumW2.at(i);
    double frac1sq  = var1/(binCont1*binCont1);
    double frac2sq  = var2/(binCont2*binCont2);
    
    double content = binCont1 / binCont2;
    double var     = content*content*(frac1sq + frac2sq);
    
    if (binCont2 == 0.0) content = 0.0;

    _binContents.at(i) = content;
    _sumW2.at(i)       = var;

  }
  
  //reset min and max to default values
  _min = -999.9;
  _max = -999.9;
}

void HistogramBase::multiply(const HistogramBase& other){

  if (other._binContents.size() != _binContents.size()){
    ERROR_LOG << "Trying to divide histograms with different numbers of bins. Doing nothing.";
    return;
  }

  for(unsigned int i = 0; i < _binContents.size(); i++){

    double binCont1 =       _binContents.at(i);
    double binCont2 = other._binContents.at(i);
    double var1     =       _sumW2.at(i);
    double var2     = other._sumW2.at(i);
    double frac1sq  = var1/(binCont1*binCont1);
    double frac2sq  = var2/(binCont2*binCont2);
    
    double content = binCont1 * binCont2;
    double var     = content*content*(frac1sq + frac2sq);

    _binContents.at(i) = content;
    _sumW2.at(i)       = var;

  }

  //reset min and max to default values
  _min = -999.9;
  _max = -999.9;

}

void HistogramBase::add(const HistogramBase& other){

  if (other._binContents.size() != _binContents.size()){
    ERROR_LOG << "Trying to divide histograms with different numbers of bins. Doing nothing.";
    return;
  }

  for(unsigned int i = 0; i < _binContents.size(); i++){

    double binCont1 =       _binContents.at(i);
    double binCont2 = other._binContents.at(i);
    double var1     =       _sumW2.at(i);
    double var2     = other._sumW2.at(i);
    
    double content = binCont1 + binCont2;
    double var     = var1 + var2;

    _binContents.at(i) = content;
    _sumW2.at(i)       = var;

  }

  //reset min and max to default values
  _min = -999.9;
  _max = -999.9;

}

void HistogramBase::minus(const HistogramBase& other){

  if (other._binContents.size() != _binContents.size()){
    ERROR_LOG << "Trying to divide histograms with different numbers of bins. Doing nothing.";
    return;
  }

  for(unsigned int i = 0; i < _binContents.size(); i++){

    double binCont1 =       _binContents.at(i);
    double binCont2 = other._binContents.at(i);
    double var1     =       _sumW2.at(i);
    double var2     = other._sumW2.at(i);
    
    double content = binCont1 - binCont2;
    double var     = var1 + var2;

    _binContents.at(i) = content;
    _sumW2.at(i)       = var;

  }

  //reset min and max to default values
  _min = -999.9;
  _max = -999.9;

}

void HistogramBase::pulls(const HistogramBase& other){

  if (other._binContents.size() != _binContents.size()){
    ERROR_LOG << "Trying to divide histograms with different numbers of bins. Doing nothing.";
    return;
  }

  for( int i = 0; i < getNBins(); i++){

    double binCont1 =       _binContents.at(i);
    double binCont2 = other._binContents.at(i);
    double var1     =       _sumW2.at(i);
    double var2     = other._sumW2.at(i);

    double content = binCont1 - binCont2;
    double var     = var1 + var2;

    _binContents.at(i) = content/sqrt(var);
    _sumW2.at(i)       = 0.0;

  }
  
  //reset min and max to default values
  _min = -999.9;
  _max = -999.9;  
}


void HistogramBase::pulls(const HistogramBase& other1, const HistogramBase& other2){

  if (other1._binContents.size() != _binContents.size() || other2._binContents.size() != _binContents.size()){
    ERROR_LOG << "Trying to divide histograms with different numbers of bins. Doing nothing.";
    return;
  }

  for( int i = 0; i < getNBins(); i++){

    double binCont1 = other1._binContents.at(i);
    double binCont2 = other2._binContents.at(i);
    double var1     = other1._sumW2.at(i);
    double var2     = other2._sumW2.at(i);

    double content = binCont1 - binCont2;
    double var     = var1 + var2;

    _binContents.at(i) = content/sqrt(var);
    _sumW2.at(i)       = 0.0;

  }
  
  //reset min and max to default values
  _min = -999.9;
  _max = -999.9;  
}


void HistogramBase::asymmetry(const HistogramBase& other){

  if (other._binContents.size() != _binContents.size()){
    ERROR_LOG << "Trying to divide histograms with different numbers of bins. Doing nothing.";
    return;
  }

  for( int i = 0; i < getNBins(); i++){

    double binCont1 =       _binContents.at(i);
    double binCont2 = other._binContents.at(i);
    double var1     =       _sumW2.at(i);
    double var2     = other._sumW2.at(i);

    double content = (binCont1 - binCont2)/(binCont1 + binCont2);
    
    double dzdx  = (2.0*binCont2)/((binCont1 + binCont2)*(binCont1 + binCont2));
    double dzdy  = (2.0*binCont1)/((binCont1 + binCont2)*(binCont1 + binCont2));
    double varsq = dzdx*dzdx*var1*var1 + dzdy*dzdy*var2*var2;

    _binContents.at(i) = content;
    _sumW2.at(i)       = varsq;

  }
  
  //reset min and max to default values
  _min = -999.9;
  _max = -999.9;  

}

void HistogramBase::asymmetry(const HistogramBase& other1, const HistogramBase& other2){

  if (other1._binContents.size() != _binContents.size() || other2._binContents.size() != _binContents.size()){
    ERROR_LOG << "Trying to divide histograms with different numbers of bins. Doing nothing.";
    return;
  }

  for( int i = 0; i < getNBins(); i++){

    double binCont1 = other1._binContents.at(i);
    double binCont2 = other2._binContents.at(i);
    double var1     = other1._sumW2.at(i);
    double var2     = other2._sumW2.at(i);

    double content = (binCont1 - binCont2)/(binCont1 + binCont2);
    
    double dzdx  = (2.0*binCont2)/((binCont1 + binCont2)*(binCont1 + binCont2));
    double dzdy  = (2.0*binCont1)/((binCont1 + binCont2)*(binCont1 + binCont2));
    double varsq = dzdx*dzdx*var1*var1 + dzdy*dzdy*var2*var2;

    _binContents.at(i) = content;
    _sumW2.at(i)       = varsq;

  }
  
  //reset min and max to default values
  _min = -999.9;
  _max = -999.9;  

}



void HistogramBase::randomiseWithinErrors(int seed){

  TRandom3 random(seed);

  for( int i = 0; i < getNBins(); i++){

    double mean  = _binContents.at(i);
    double sigma = sqrt(_sumW2.at(i));

    double newContent = random.Gaus(mean, sigma);
    if (newContent < 0.0) newContent = 0.0;
    _binContents.at(i) = newContent;
  }

}

double HistogramBase::pvalue(const HistogramBase& other, int ndof) const{
  
  if (ndof < 0.0) ndof = getNBins();
  double chi2 = this->chi2(other);

  return TMath::Prob(chi2, ndof);

}

double HistogramBase::chi2sig(const HistogramBase& other, int ndof) const{
  
  double pval = pvalue(other, ndof);
  double sigma = TMath::NormQuantile( pval/2.0 );
  return fabs(sigma);

}


double HistogramBase::chi2(const HistogramBase& other) const{

  double chi2 = 0.0;

  for( int i = 0; i < getNBins(); i++){

    double binCont1 =       _binContents.at(i);
    double binCont2 = other._binContents.at(i);
    double var1     =       _sumW2.at(i);
    double var2     = other._sumW2.at(i);
    
    double content = binCont1 - binCont2;
    double var     = var1 + var2;

    double thisTerm = (content * content) / var;

    chi2 += thisTerm;

  }

  return chi2;

}

void HistogramBase::drawPullHistogram(const HistogramBase& other, TString name, int nBins, double pmLimits) const{

  TH1D pulls ("pulls", "pulls", nBins, -pmLimits, pmLimits);
  pulls.GetXaxis()->SetTitle("pull");
  pulls.GetYaxis()->SetTitle("frequency");
  for( int i = 0; i < getNBins(); i++){

    double binCont1 =       _binContents.at(i);
    double binCont2 = other._binContents.at(i);
    double var1     =       _sumW2.at(i);
    double var2     = other._sumW2.at(i);
    
    double content = binCont1 - binCont2;
    double var     = var1 + var2;

    double pull = (content) / sqrt(var);
    pulls.Fill(pull);

  }

  RootPlotter1D plotter(&pulls);
  plotter.plot(name);

}

double HistogramBase::integral() const{
  double integral = 0.0;
  for(int i = 0; i < getNBins() ; i++){
    integral += _binContents.at(i);
  }
  return integral;
}

double HistogramBase::integralError() const{
  double var = 0.0;
  for(int i = 0; i < getNBins() ; i++){
    var += _sumW2.at(i);
  }
  return sqrt(var);
}

void HistogramBase::normalise(double area){

  double divisor  = integral()/area;
  double divisor2 = divisor*divisor;

  for( int i = 0; i < getNBins(); i++){
    _binContents.at(i) = _binContents.at(i)/divisor;
    _sumW2      .at(i) = _sumW2      .at(i)/divisor2;
  }

}

void HistogramBase::setBinContent(int bin, double val){
  bin = checkBinNumber(bin);
  _binContents[bin] = val;
}

void HistogramBase::setBinError  (int bin, double val){
  bin = checkBinNumber(bin);
  _sumW2[bin] = val*val;
}

double HistogramBase::getBinContent(int bin) const{
  bin = checkBinNumber(bin); 
  return _binContents[bin];
}

double HistogramBase::getBinError  (int bin) const{
  bin = checkBinNumber(bin);
  return sqrt(_sumW2[bin]);
}

double HistogramBase::getBinVolume(int bin) const{
  bin = checkBinNumber(bin);
  return 1.0;
}

double HistogramBase::getFrequencyDensity(int bin) const{

  return getBinContent(bin)/getBinVolume(bin);

}

void HistogramBase::makeFrequencyDensity(){
  
  for(int i = 0; i < _nBins; i++){
    double binVolume = this->getBinVolume(i);
    _binContents[i] = _binContents[i] / binVolume;
    _sumW2[i] = _sumW2[i] / (binVolume*binVolume);
  }

}

void HistogramBase::print(){

  for(int i = 0; i < _nBins; i++){
    INFO_LOG << "Bin Content " << i << ": " << _binContents[i] << "      SumW2: " << _sumW2[i];
  }
  INFO_LOG << "Overflow: " << _binContents[_nBins];

}