void drawgraph(double *tau_, double *ll_, const int ntry){   
   TCanvas *c2 = new TCanvas("c2","log likelihood example",200,10,700,500);
   c2->SetFillColor(42);
   c2->SetGrid();


   TGraph *gr = new TGraph(ntry,tau_,ll_);
   gr->SetLineColor(2);
   gr->SetLineWidth(4);
   gr->SetMarkerColor(4);
   gr->SetMarkerStyle(21);
   gr->SetTitle("ll vs tau");
   gr->GetXaxis()->SetTitle("tau");
   gr->GetYaxis()->SetTitle("log likelihood");
   gr->Draw("ACP");

   // TCanvas::Update() draws the frame, after which one can change it
   c2->Update();
   c2->GetFrame()->SetFillColor(21);
   c2->GetFrame()->SetBorderSize(12);
   c2->Modified();
}

void makedata(double *data, const int N, double tau, int seed){

  //makedata generates some exponentially distributes
  //numbers. Think that these are measured time intervals
  //between decays measured in an experiment
  //We are doing a pseudo-experiment (also called toy 
  //Monte Carlo experiment)

  //data is an array of double, to store the "measured"
  //decay times. 

  //TRandom2 is the random number generator class in ROOT
  //gRandom is a pointer to a TRandom2 object, that 
  //ROOT has by default whenever you enter a ROOT session
  //But we can make our own TRandom2 objects 
  //If you want to know about classes and objects in C++
  //you can follow one of the links given in the course homepage
  //

  TRandom2 r;

  r.SetSeed(seed);
  for(int ievent=0;ievent<N;ievent++){
      data[ievent] = r.Exp(tau);
  }
}

void histdata(double *data, TH1F *histo, const int N){

  for (int j = 0; j<N; j++){
    histo->Fill(data[j]);
  }

}

void  LLRdistribution(double *t, const int N, double tau, double hypo_tau0, double hypo_tau1, const int ntrial, TH1F *h0){
  //generate the distribution of LLR under a given hypothesis
  //THIS IS TO BE DONE BEFORE THE REAL DATA COMES
  //Note: if you could somehow analytically calculate the distribution
  //      of the LLR, this toy MC step will not be required 

  for (int itrial = 0; itrial < ntrial; itrial++){
    //this is like data taking of experiment
    //this call will fill the array t[] with "measured" decay times
    int seed = itrial + 54549;
    double ll1, ll0;

    makedata(t,N,tau,seed);
    
    //calculate llr
    ll0=ll1=0;  
    for (int ii = 0; ii < N; ii++){
      ll1 += -log(hypo_tau1) - (t[ii]/hypo_tau1);
      ll0 += -log(hypo_tau0) - (t[ii]/hypo_tau0);
    }
    
    double t_tau = ll1-ll0; //this is the log likelihood ratio (LLR) test statis

    //print ll0, ll1 every 1000 times
    if (!(itrial%100)) cout<<"logL for hypo0 "<<ll0<<" and for hypo1 "<<ll1<<endl;
    //fill the histo 
    h0->Fill(t_tau);
  }
}

void testmodel(){
  const int N = 100, ntrial = 10000;
  double t[N];
  
  double hypo_tau0 = 2.0, hypo_tau1=2.25, true_tau=2.2;
  
  //histogram for P(t_tau | H0)
  TH1F *h0 = new TH1F("h0", "P(t|H0)", 100, -5., 5.);
  TH1F *h1 = new TH1F("h1", "P(t|H1)", 100, -5., 5.);

  //generate the distribution of llr under hypothesis 0 (null hypothesis)
  LLRdistribution(t, N, hypo_tau0, hypo_tau0, hypo_tau1, ntrial, h0);

  //draw the histo
  h0->Draw();
    
  //generate the distribution of llr under hypothesis 1 (alt hypothesis)
  LLRdistribution(t, N, hypo_tau1, hypo_tau0, hypo_tau1, ntrial, h1);


  //draw the histo
  h1->SetLineColor(6);
  h1->Draw("same");
    

  //graph drawing
  //drawgraph(tau_,ll_, ntry); 
}