JRootTestkit.hh

Go to the documentation of this file.

#ifndef __JROOT__JROOTTESTKIT__
#define __JROOT__JROOTTESTKIT__
 
#include <cmath>
#include <cstddef>
#include <vector>
#include <algorithm>
 
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wall"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TRandom3.h"
#pragma GCC diagnostic pop
 
 
/**
 * \author mdejong
 */
 
namespace JROOT {}
namespace JPP { using namespace JROOT; }
 
namespace JROOT {
 
  /**
   * Fill 1D histogram according Poisson statistics with expectation values from given 1D function.
   *
   * \param  h1            histogram
   * \param  f1            function
   */
  template<class T>
  inline void FillRandom(TH1* h1, const T& f1)
  {
    for (Int_t ix = 1; ix <= h1->GetXaxis()->GetNbins(); ++ix){
 
      const double x = h1->GetXaxis()->GetBinCenter(ix);
 
      h1->SetBinContent(ix, gRandom->Poisson(f1(x)));
 
      if (h1->GetSumw2N() != 0) {
        h1->SetBinError(ix, sqrt(h1->GetBinContent(ix)));
      }
    }
  }
 
 
  /**
   * Fill 2D histogram according Poisson statistics with expectation values from given 2D function.
   *
   * \param  h2            histogram
   * \param  f2            function
   */
  template<class T>
  inline void FillRandom(TH2* h2, const T& f2)
  {
    for (Int_t ix = 1; ix <= h2->GetXaxis()->GetNbins(); ++ix){
      for (Int_t iy = 1; iy <= h2->GetYaxis()->GetNbins(); ++iy){
 
        const double x = h2->GetXaxis()->GetBinCenter(ix);
        const double y = h2->GetYaxis()->GetBinCenter(iy);
        
        h2->SetBinContent(ix, iy, gRandom->Poisson(f2(x,y)));
 
        if (h2->GetSumw2N() != 0) {
          h2->SetBinError(ix, iy, sqrt(h2->GetBinContent(ix,iy)));
        }
      }
    }
  }
 
 
  /**
   * Fill 3D histogram according Poisson statistics with expectation values from given 3D function.
   *
   * \param  h3            histogram
   * \param  f3            function
   */
  template<class T>
  inline void FillRandom(TH3* h3, const T& f3)
  {
    for (Int_t ix = 1; ix <= h3->GetXaxis()->GetNbins(); ++ix){
      for (Int_t iy = 1; iy <= h3->GetYaxis()->GetNbins(); ++iy){
        for (Int_t iz = 1; iz <= h3->GetZaxis()->GetNbins(); ++iz){
 
          const double x = h3->GetXaxis()->GetBinCenter(ix);
          const double y = h3->GetYaxis()->GetBinCenter(iy);
          const double z = h3->GetZaxis()->GetBinCenter(iz);
        
          h3->SetBinContent(ix, iy, iz, gRandom->Poisson(f3(x,y,z)));
 
          if (h3->GetSumw2N() != 0) {
            h3->SetBinError(ix, iy, iy, sqrt(h3->GetBinContent(ix,iy,iz)));
          }
        }
      }
    }
  }
 
 
  /**
   * Fill 1D histogram according PDF as given 1D function.
   *
   * \param  h1            histogram
   * \param  f1            function
   * \param  ns            number of entries
   */
  template<class T>
  inline void FillRandom(TH1* h1, const T& f1, const size_t ns)
  {
    using namespace std;
 
    struct tuple_t {
 
      Double_t x;
      double   W;
 
      inline bool operator<(const double W) const { return this->W < W; }
    };
 
    vector<tuple_t> buffer;
 
    double W = 0.0;
 
    for (Int_t ix = 1; ix <= h1->GetXaxis()->GetNbins(); ++ix){
 
      const Double_t x = h1->GetXaxis()->GetBinCenter(ix);
      
      W += f1(x);
 
      buffer.push_back({x, W});
    }
 
    for (size_t i = 0; i != ns; ++i) {
 
      typename vector<tuple_t>::const_iterator p = lower_bound(buffer.begin(), buffer.end(), W * gRandom->Rndm());
 
      h1->Fill(p->x);
    }
  }
 
 
  /**
   * Fill 2D histogram according PDF as given 2D function.
   *
   * \param  h2            histogram
   * \param  f2            function
   * \param  ns            number of entries
   */
  template<class T>
  inline void FillRandom(TH2* h2, const T& f2, const size_t ns)
  {
    using namespace std;
 
    struct tuple_t {
 
      Double_t x;
      Double_t y;
      double   W;
 
      inline bool operator<(const double W) const { return this->W < W; }
    };
 
    vector<tuple_t> buffer;
 
    double W = 0.0;
 
    for (Int_t ix = 1; ix <= h2->GetXaxis()->GetNbins(); ++ix){
      for (Int_t iy = 1; iy <= h2->GetYaxis()->GetNbins(); ++iy){
 
        const Double_t x = h2->GetXaxis()->GetBinCenter(ix);
        const Double_t y = h2->GetYaxis()->GetBinCenter(iy);
      
        W += f2(x,y);
 
        buffer.push_back({x, y, W});
      }
    }
 
    for (size_t i = 0; i != ns; ++i) {
 
      typename vector<tuple_t>::const_iterator p = lower_bound(buffer.begin(), buffer.end(), W * gRandom->Rndm());
 
      h2->Fill(p->x, p->y);
    }
  }
 
 
  /**
   * Fill 3D histogram according PDF as given 3D function.
   *
   * \param  h3            histogram
   * \param  f3            function
   * \param  ns            number of entries
   */
  template<class T>
  inline void FillRandom(TH3* h3, const T& f3, const size_t ns)
  {
    using namespace std;
 
    struct tuple_t {
 
      Double_t x;
      Double_t y;
      Double_t z;
      double   W;
 
      inline bool operator<(const double W) const { return this->W < W; }
    };
 
    vector<tuple_t> buffer;
 
    double W = 0.0;
 
    for (Int_t ix = 1; ix <= h3->GetXaxis()->GetNbins(); ++ix){
      for (Int_t iy = 1; iy <= h3->GetYaxis()->GetNbins(); ++iy){
        for (Int_t iz = 1; iz <= h3->GetYaxis()->GetNbins(); ++iz){
 
          const Double_t x = h3->GetXaxis()->GetBinCenter(ix);
          const Double_t y = h3->GetYaxis()->GetBinCenter(iy);
          const Double_t z = h3->GetZaxis()->GetBinCenter(iz);
      
          W += f3(x,y,z);
 
          buffer.push_back({x, y, z, W});
        }
      }
    }
 
    for (size_t i = 0; i != ns; ++i) {
 
      typename vector<tuple_t>::const_iterator p = lower_bound(buffer.begin(), buffer.end(), W * gRandom->Rndm());
 
      h3->Fill(p->x, p->y, p->z);
    }
  }
}
 
#endif