#include "TRandom3.h"
#include "TH1.h"
#include "func.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TApplication.h"

double event_generator(TRandom3* EG, double npe = 1., double gain = 60., double sig = 6., double xtalk = 0.15, double ped = 100., double sig_ped = 7.)
{
  int npixel = EG->Poisson(npe);
  npixel += EG->Poisson((double)npixel*xtalk);
  double width = sqrt(sig_ped*sig_ped + (double)npixel*sig*sig);
  double tmpx = EG->Gaus(ped+gain*(double)npixel,width);
  return tmpx;
}

void gen()
{
  gStyle->SetOptFit(1);

//   double factor[4];
//   factor[0] = atof(gApplication->Argv(2));
//   factor[1] = atof(gApplication->Argv(3));
//   factor[2] = atof(gApplication->Argv(4));
//   factor[3] = atof(gApplication->Argv(5));
//   int xtalk_flag = atoi(gApplication->Argv(6));

//   cout<<"scaling factors are: "<<endl;
//   cout<<"channel 0: "<<factor[0]<<endl;
//   cout<<"channel 1: "<<factor[1]<<endl;
//   cout<<"channel 2: "<<factor[2]<<endl;
//   cout<<"channel 3: "<<factor[3]<<endl;
//   cout<<"xtalk flag: "<<xtalk_flag<<endl;

  int xtalk_flag = 0;

  TRandom3 *EG = new TRandom3(0);
  
  double npe = 0.5;
  double gain = 60;
  double sig = 6;
  double xtalk = 0;//0.15;
  double ped = 50;
  double sig_ped = 7;
  
  int ntot = 4;
  TH1F* hh[16];
  char hname[256], htitle[256];

  for(int i = 0; i<ntot; i++)
    {
      sprintf(hname,"histo_ch%02d",i);
      sprintf(htitle,"ADC Channel #%02d; ADC; Counts",i);
      hh[i] = new TH1F(hname,htitle,1000,0,1000);
    }
  TH1F* hh_sum = new TH1F("histo_total","Sum of all channels; ADC; Counts", 2000,0,2000);
  
  for(int i = 0; i<10000000; i++)
    {
      double adc_sum = 0;
      for(int j = 0; j<ntot; j++)
	{
	  //	  double adc_tmp = event_generator(EG, npe, gain*factor[j], sig*sqrt(factor[j]), xtalk, ped, sig_ped);
	  double adc_tmp = event_generator(EG, npe, gain, sig, xtalk, ped, sig_ped);
	  hh[j]->Fill(adc_tmp);
	  adc_sum += adc_tmp;
	}
      hh_sum->Fill(adc_sum);
      if((i+1)%500000 == 0) printf("%d events generated ...\n", i+1);
    }
  TCanvas* c1 = new TCanvas("c1","",800,1000);

  TPad* cp[16];
  cp[0] = new TPad("c1_1","c1_1",0.01,0.71,0.49,0.99);
  cp[1] = new TPad("c1_2","c1_2",0.51,0.71,0.99,0.99);
  cp[2] = new TPad("c1_3","c1_3",0.01,0.41,0.49,0.69);
  cp[3] = new TPad("c1_4","c1_4",0.51,0.41,0.99,0.69);
  cpt = new TPad("c1_t","c1_t",0.01,0.01,0.99,0.39);

  double par[7], par_err[7];
//   gStyle->SetStatW(0.25);
//   gStyle->SetStatH(0.25);
  

  for(int i = 0; i<ntot; i++)
    {
      c1->cd();
      cp[i]->Draw();
      cp[i]->cd();
      hh[i]->Draw();
//       par[5] = 50;
//       fit(hh[i],par,par_err,xtalk_flag);
      c1->Update();
    }

//   gStyle->SetStatW(0.3);
//   gStyle->SetStatH(0.4);
  c1->cd();
  cpt->Draw();
  cpt->cd();
  hh_sum->Draw();
  c1->Update();
  //  par[5] = 200;
  //  fit(hh_sum,par,par_err,xtalk_flag);
  
//   sprintf(hname,"plot_%03d_%03d_%03d_%03d.png",int(factor[0]*100.+0.5),int(factor[1]*100.+0.5),int(factor[2]*100.+0.5),int(factor[3]*100.+0.5));
//   c1->Print(hname);
// //   sprintf(hname,"plot_tot_%03d_%03d_%03d_%03d.png",int(factor[0]*100.),int(factor[1]*100.),int(factor[2]*100.),int(factor[3]*100.));
// //   c2->Print(hname);
//   sprintf(hname,"result_%03d_%03d_%03d_%03d.txt",int(factor[0]*100.+0.5),int(factor[1]*100.+0.5),int(factor[2]*100.+0.5),int(factor[3]*100.+0.5));
//   FILE* result = fopen(hname,"w");
//   fprintf(result,"Set Parameters: \n\n");
//   fprintf(result,"npe: \t%.6f\n",npe);
//   fprintf(result,"gain: \t%.6f\n",gain);
//   fprintf(result,"sigma: \t%.6f\n",sig);
//   fprintf(result,"xtalk: \t%.6f\n",xtalk);
//   fprintf(result,"ped: \t%.6f\n",ped);
//   fprintf(result,"sig_ped: \t%.6f\n",sig_ped);
//   for(int i = 0; i<ntot; i++)
//     fprintf(result,"rel_gain ch %d: \t %.6f\n",i,factor[0]);
//   fprintf(result,"\nFit Parameters: \n\n");
//   fprintf(result,"Chi-sqare: \t%.6f\n",hh_sum->GetFunction("pf")->GetChisquare());
//   fprintf(result,"ndf: \t%d\n",hh_sum->GetFunction("pf")->GetNDF());
//   fprintf(result,"npe: \t%.6f\t+-\t%.6f\n",par[0],par_err[0]);
//   fprintf(result,"gain: \t%.6f\t+-\t%.6f\n",par[3],par_err[3]);
//   fprintf(result,"sigma: \t%.6f\t+-\t%.6f\n",par[1],par_err[1]);
//   fprintf(result,"xtalk: \t%.6f\t+-\t%.6f\n",par[6],par_err[6]);
//   fprintf(result,"ped: \t%.6f\t+-\t%.6f\n",par[5],par_err[5]);
//   fprintf(result,"sig_ped: \t%.6f\t+-\t%.6f\n",par[4],par_err[4]);
//   fprintf(result,"amplitude: \t%.6f\t+-\t%.6f\n",par[2],par_err[2]);
}