generate.h

Go to the documentation of this file.

/* emacs: this is -*- c++ -*- */
 
#ifndef GENERATE_H
#define GENERATE_H
 
 
#include <cmath>
#include <vector>
#include <iostream>
 
#include "TH1D.h"
 
#include "BasicRandom.h"
 
 
namespace Resplot::generate { 

double get_mean( const std::vector<double>& m );
double get_rms( const std::vector<double>& m );
 

TH1D* integrate( TH1D* hin );

void ScaleIntegral(TH1D* hin );

TH1D* PDF(TH1D* hin);

TH1D* smooth( TH1D* hin, bool sym=true );
 

class generator_base { 
public:
  virtual ~generator_base() { }
  virtual double generate() const = 0;
};
 
 

class breit_generator : public generator_base { 
 
public: 
  
  virtual double generate() const {
    double x = std::tan((r.uniform()-0.5)*M_PI);
    return x;
  }
  
  static inline thread_local BasicRandom r;
 
};
 
 
 

class hist_generator : public generator_base { 
 
public:
 
  hist_generator(TH1D* h, bool smooth_=true );
  
  hist_generator(const hist_generator &) = delete;
  hist_generator operator=(const hist_generator &) = delete;
  
  virtual ~hist_generator() { delete m_s; } // if ( m_random ) delete m_random; }

  double generate() const {
    return invert(s_random.uniform());
  }
 
  TH1D* histogram()       { return m_s; }
  TH1D* rawhistogram()    { return m_raw; }
  TH1D* smoothhistogram() { return m_smooth; }
 
private:
 
  int getbin( double y ) const {
    for ( unsigned i=0 ; i<m_y.size() ; i++ ) if ( y<=m_y[i] ) return i;
    return m_y.size()-1;
  } 
    
  double invert( double y ) const {
    int i = getbin(y);
    if ( i==0 ) return 0;
    else return (y-m_y[i-1])*m_dxdy[i-1]+m_x[i-1];
  }
 
 
private: 
 
  std::vector<double> m_y;
  std::vector<double> m_x;
  
 
  std::vector<double> m_dx;
  std::vector<double> m_dy;
  
  std::vector<double> m_dxdy;
 
  TH1D* m_s;   // accumated histogram
  TH1D* m_raw;    // raw histogram
  TH1D* m_smooth; // smoothed raw histogram
  
  //  static BasicRandom* m_random;
 
  static inline thread_local BasicRandom s_random = BasicRandom(true);
  
};
 
 
 

 
class experiment { 
 
public:
 
  experiment( TH1D* h, int Nexperiments=10, int fevents=0 );
  experiment(const experiment& ) = delete;
  experiment operator=(const experiment& ) = delete;
 
  double hmean() const { return m_hmean; }  
  double hrms()  const { return m_hrms; }  
 
  double mean()       const { return m_global_mean; }  
  double mean_error() const { return m_global_mean_error; }  
 
  double rms()       const { return m_global_rms; }  
  double rms_error() const { return m_global_rms_error; }  
 
  generator_base* gen() { return m_gen; }
 
  TH1D* rmshisto() { return m_THrms; } 
  
private:
 
  int m_N; 
 
  std::vector<double> m_mean;  
  std::vector<double> m_rms;   
 
  double m_hmean; 
  double m_hrms;  

  double m_global_mean;
  double m_global_mean_error;
 
  double m_global_rms;
  double m_global_rms_error;
    
  generator_base* m_gen;
 
  TH1D* m_THrms;
};
 
 
 
} // namespace generate
 
#endif // GENERATE_H