void mc_integration()
{
  enum { ACCEPT, TRANSFORM };
  //int mode = ACCEPT;    
  int mode = TRANSFORM;    
  double cut  = 5;
  double xmax = 10;
  int N        = 100000;
  int n_accept = 0;
  int n_trans  = 0;
  double trueval = 2*(1-exp(-cut));

  TF1* f = new TF1( "f", func, 0, xmax, 1 );
  f->SetParameter( 0, cut );
  f->SetNpx( 1000 );
  f->Draw();
  
  // accept/reject 
  if( mode == ACCEPT ) {
    for( int i = 1; i <= N; ++i ) {
      
      double x = gRandom->Uniform(xmax);
      double y = gRandom->Uniform();
      if( f->Eval( x ) > y ) ++n_accept;
       
      // total area: xmax * 1
      if( i % 10000 == 0 ) {
	cout << "Accept/Reject:\t" <<  n_accept << "\t" << i << "\t" << xmax*double(n_accept)/double(i) << "\t" << trueval << endl;
      }
    }
  } 
  // transformation
  else if( mode == TRANSFORM ) {
    for( int i = 1; i <= N; ++i ) {

      // transformation method: -ln(y) and twice the integral to the cut
      double y = gRandom->Rndm();
      double trans = -log( y );
      if( trans < cut ) ++n_trans;
      if( i % 10000 == 0 ) {
	cout << "Transformation:\t" << n_trans << "\t" << i << "\t" << 2*double(n_trans)/double(i) << "\t" << trueval << endl; 
      }
    }
  }
}

// function with a sharp peak
// analytical integral: 2(1-exp(-par[0])), for par[0] = 5: 1.98652410600182905e+00
double func( double* x, double* par ) 
{
  double x0 = x[0] - par[0];
  if( x[0] > 5 ) { 
    return exp( -x0 );
  } else {
    return exp( x0 );
  }
}