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){

  //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;

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

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

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

}

void fitmodel(){
  const int N = 1000;
  double t[N];

  //this is like data taking of experiment
  //this call will fill the array t[] with "measured" decay times
  makedata(t, N);

  //histogram for data for visualization
  TH1F *hdata = new TH1F("data", "gaussian sig + exp bkg", 100, 0., 20.);
  //fill the histo 
  histdata(t,hdata,N);
  //draw the histo
  hdata->Draw();

  const int ntry = 500;
  double tau_[ntry], taumax = 0.0, llmax = -10000000.0, ll_[ntry];

  double tau_low = 2.0, tau_high = 7.0;
  double dtau = (tau_high - tau_low)/double(ntry);
  int jj = 0;
  for (double tau = tau_low; tau <= tau_high; tau+=dtau){
    double ll = 0.0;  
    for (int ii = 0; ii < N; ii++){
      ll += -log(tau) - (t[ii]/tau);
    } 

    if (ll>llmax) {llmax = ll; taumax = tau;} 
    ll_[jj] = ll;
    tau_[jj] = tau;

    jj++;
    cout<<jj<<" tau = "<<tau<<" log likelihood = "<<ll<<endl;
  }
  cout<<"max of log likelihood "<<llmax<<" occurs at "<<taumax<<endl;

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

 
}