#include <iostream>
#include "MyRandom.h"
#include "TMath.h"

ClassImp(MyRandom)
  ///////////////////////////////////////////////////////
  // Class used to extend the number of functions managed
  // by the class TRandom3 in order to generate random numbers
  // Origin: M. Masera 17/10/02
  ////////////////////////////////////////////////////////
//_________________________________________________________________
  MyRandom::MyRandom(UInt_t seed) {
  // Default constructor
  SetName("MyRandom");
  SetTitle("Random number generator: Mersenne-Twistor");
  SetSeed(seed);
}
//_________________________________________________________________
MyRandom::~MyRandom() {
  // Default destructor
}
//_________________________________________________________________
Double_t MyRandom::Funct1(Double_t theta,Double_t alpha) {
  // Function funct1:  1/(sin**2 theta + alpha*cos**2 theta)
  // with positive alpha
  if(alpha<=0.){
    cerr << "MyRandom::Funct1: Invalid parameter alpha = "<<alpha<<endl;
    return 0;
  }
  Double_t tmp1 = TMath::Sin(theta);
  tmp1*=tmp1;
  Double_t tmp2 = TMath::Cos(theta);
  tmp2*=tmp2;
  tmp1=tmp1+alpha*tmp2;
  return 1./tmp1;
}
//_________________________________________________________________
Double_t MyRandom::Funct1RndmByInversion(Double_t alpha) {
  // returns a deviate distributed according to Funct1, using the 
  // inversion technique
 
  if(alpha<=0.){
    cerr << "MyRandom::Funct1RndmByInversion: ";
    cerr<<"Invalid parameter alpha = "<<alpha<<endl;
    return 0;
  }

  Double_t tmp1 = TMath::Sqrt(alpha);
  Double_t tmp2 = TMath::Tan(2.*TMath::Pi()*Rndm());
  Double_t resu = TMath::ATan(tmp1*tmp2);
  // TMath::ATan() return a value between -pi/2 and pi/2
  // a uniform random number is thrown in order to extend the range of
  // the sampled funtion over 2*pi.
 
  if(Rndm()<0.5){
    // pi/2 < theta < 3/2*pi
    resu =  TMath::Pi() - resu;
  }
  else {
    // 0 < theta < pi/2 and 3/2*pi < theta < 2*pi
    if(resu<0.)resu += 2.*TMath::Pi();
  }

  return resu;
}
//_________________________________________________________________
Double_t MyRandom::Funct1RndmByRejection(Double_t alpha) {
  // returns a deviate distributed according to Funct1, using the 
  // rejection technique
  if(alpha<=0.){
    cerr << "MyRandom::Funct1RndmByRejection: ";
    cerr<<"Invalid parameter alpha = "<<alpha<<endl;
    return 0;
  }
  Bool_t throwing = kTRUE;
  Double_t fbig=1./alpha;
  Double_t x,y;
  while(throwing){
    x=2.*Rndm()*TMath::Pi();
    y=Rndm()*fbig;
    if(Funct1(x,alpha)>=y)throwing = kFALSE;
  }
  return x;
}