void makedata(double *data, const int N, double mu, double sig){

  //makedata generates some Gaussian 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" data

  //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(7134554);

  for(int ievent=0;ievent<N;ievent++){
      data[ievent] = r.Gaus(mu,sig);
  }
}


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

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

}


void fitmodel(){
    //N is the number of data points
    const int N = 10000;
    const float pi = 3.1415;
    //array to store decay times
    double t[N],  ll;
    const int ntry = 50;
    
    double mu_[ntry], mumax = 0.0, llmax = -10000000.0, ll_[ntry];
    double s_[ntry], smax = 0.0;
    
    double mu_low = 120.0, mu_high = 130.0;
    double s_low = 0.5, s_high = 2.5;                // s => sig^2
    
    double dmu = (mu_high - mu_low)/double(ntry);
    double ds =  (s_high - s_low)/double(ntry);

    
  //this is like data taking of experiment
  //this call will fill the array t[] with "measured" Gaussian data
  double true_mu = 125.0, true_sig = 1.0;
  makedata(t, N, true_mu,true_sig);

  //histogram for data for visualization
    
    TCanvas *c1 = new TCanvas("c1","multipads",900,700); // root command for plotting many figures
    gStyle->SetOptStat(0);
    c1->Divide(2,2);
    
    
  TH1F *hdata = new TH1F("data", "gaussian sig", 100, true_mu - 5.*true_sig, true_mu + 5.*true_sig);
  TH2F * h2a = (TH2F *) new TH2F("h2a","LL map",ntry,mu_low,mu_high,ntry,s_low,s_high);


  //fill the histo
  histdata(t,hdata,N);
    
  //draw the histo
  c1->cd(1);
  hdata->Draw();


    
  int jj = 0;
  for (double mu = mu_low; mu <= mu_high; mu+=dmu){
      
      for (double s = s_low; s <= s_high; s+=ds){
          
        double ll = 0.0;
          
          for (int ii = 0; ii < N; ii++){
             ll += - 0.5*pow((t[ii]-mu), 2)/s  - 0.5*log(s)  - 0.5*log(2.0*pi) ;
              
              
          }

      if (ll > llmax) {llmax = ll; mumax = mu; smax=s;}
          
          h2a->Fill(mu,s,ll);


    jj++;
    cout<<jj<<" mu = "<<mu<<" sig^2 = "<<s <<" log likelihood = "<<ll<<endl;
      }
  }
  cout<<"max of log likelihood "<<llmax<<" occurs at mu = "<<mumax<<" sig = "<<sqrt(smax)<<endl;

  //graph drawing
  //create a graph to plot of ll in mu Vs sig^2
    
    c1->cd(2);
    h2a->Draw("lego");
    
    c1->cd(3);
    h2a->Draw("surf2z");
    
//    c1->cd(4);
//    h2a->Draw("colz");
    
    c1->cd(4);
    h2a->Draw("contz");
    
//    c1->cd(4);
//    h2a->Draw("cont1");
    
    h2a->SetTitle("log likelihood");
    h2a->GetXaxis()->SetTitle("mu");
    h2a->GetYaxis()->SetTitle("sig^2");

    
}