#include <vector>

void pedestal_studies(Double_t pedestal=100)
{

gROOT->Reset();
//TTree *Bfield = (TTree *) gROOT->FindObject("Bfield");
gStyle->SetPalette(1,0);
gStyle->SetOptStat(kTRUE);
gStyle->SetOptFit(kTRUE);
gStyle->SetOptFit(1111);
gStyle->SetPadRightMargin(0.15);
gStyle->SetPadLeftMargin(0.15);
gStyle->SetPadBottomMargin(0.15); 

gStyle->SetLabelSize(0.05,"XYZ");
gStyle->SetTitleSize(0.05,"XYZ");
gStyle->SetTitleOffset(1.5,"Y");
gStyle->SetNdivisions(505,"XYZ");

// gStyle->SetFillColor(0);
//
   char string[256];
   Int_t j,jj;

   Int_t exponent=10;
   Int_t nbins=pow(2,exponent);
   cout << "exponent=" << exponent << " nbins=" << nbins << endl;
   // Double_t pedestal=100;
   Double_t ped_sigma=1.2;

   Double_t xmin=0;
   Double_t xmax=nbins;
          
   // define histogram to hold pulses

  TH1F *hwindowD = new TH1F("hwindowD","Pedestal window (double)", nbins,xmin,xmax);
  TH1I *hwindowI = new TH1I("hwindowI","Readout Window (integer)", nbins,xmin,xmax);
  
   TCanvas *c0 = new TCanvas("c0", "c0",200,10,700,700);
   c0->SetBorderMode(0);
   c0->SetFillColor(0);

   c0->SetGridx();
   c0->SetGridy();

   c0->Divide(1,2);
   c0->cd(1);
   c0_1->SetBorderMode(0);
   c0_1->SetFillColor(0);

   hwindowD->SetTitle("");  
   hwindowD->GetXaxis()->SetLabelSize(0.05);
   hwindowD->GetXaxis()->SetTitleSize(0.05);
   hwindowD->GetYaxis()->SetLabelSize(0.05);
   hwindowD->GetYaxis()->SetTitleSize(0.05);
   hwindowD->GetYaxis()->SetTitleOffset(1.5);
   hwindowD->GetYaxis()->SetTitle("Pedestal (double)");
   hwindowD->GetXaxis()->SetTitle("FADC samples");
   hwindowD->GetXaxis()->SetNdivisions(505);
   hwindowD->SetLineColor(2);
   hwindowD->Draw();  

   c0->cd(2);
   c0_2->SetBorderMode(0);
   c0_2->SetFillColor(0);

   hwindowI->SetTitle("");
   hwindowI->GetYaxis()->SetTitle("Pedestal (integer)");
   hwindowI->GetXaxis()->SetTitle("FADC samples");
   hwindowI->SetLineColor(2);
   hwindowI->Draw();     

   TRandom3 *r1 = new TRandom3();
   UInt_t  iseed = 23424111513;
   r1->SetSeed(iseed);

   // fill pedestal windows

  for (int j=0; j<nbins; j++) {        
      Double_t z1 = r1->Gaus(pedestal,ped_sigma);
      hwindowD->Fill(j,z1);
      hwindowI->Fill(j,(int) z1);
  }

  TH1D *hsumpedD[9];
  TH1D *hsumpedDI[9];
  TH1I *hsumpedI[9];
  Int_t nhist=0;

  for (int j=1; j<nbins; j*=2) {

      Int_t niter = nbins/j;
      cout << " j=" << j << " niter=" << niter << " nhist=" << nhist << endl;

       // define histogram to keep results

       sprintf (string,"SumpedD_j%03d_bins\n",j);
      printf("string=%s",string);
      hsumpedD[nhist] = new TH1D(string,string,500,95,105);
      hsumpedD[nhist]->Sumw2();

      sprintf (string,"SumpedDI_j%03d_bins\n",j);
      printf("string=%s",string);
      hsumpedDI[nhist] = new TH1D(string,string,500,95,105);
      hsumpedDI[nhist]->Sumw2();

      sprintf (string,"SumpedI_j%03d_bins\n",j);
      printf("string=%s",string);
      hsumpedI[nhist] = new TH1I(string,string,500,95,105);
      hsumpedI[nhist]->Sumw2();

      for (int jj=0; jj<niter; jj++) {
             Double_t sumpedD=0;
             Double_t sumpedDI=0;
             Int_t sumpedI=0;

             for (int k=0;k<j;k++) {
                    Int_t ndx = jj*j + k;
                    sumpedD += hwindowD->GetBinContent(ndx+1);
                    sumpedDI += hwindowI->GetBinContent(ndx+1);
                    sumpedI += hwindowI->GetBinContent(ndx+1);
                    // cout << "j=" << j << " jj=" << jj << " k=" << k << " ndx=" << ndx << " content=" << hwindowD->GetBinContent(ndx+1) << endl;
	}
              sumpedD /= j;
              sumpedDI /= j;
              sumpedI /= j;
              hsumpedD[nhist]->Fill(sumpedD);
              hsumpedDI[nhist]->Fill(sumpedDI);
              hsumpedI[nhist]->Fill(sumpedI);
              // cout << " sumpedD=" << sumpedD  << " sumpedDI=" << sumpedDI << " sumpedI=" << sumpedI << endl;                  
              }
    nhist++;

    }
      

   TCanvas *c1 = new TCanvas("c1", "c1",200,10,1000,700);
   c1->SetBorderMode(0);
   c1->SetFillColor(0);

   c1->SetGridx();
   c1->SetGridy();
   c1->Divide(3,2);

   c1->cd(1);
   c1_1->SetBorderMode(0);
   c1_1->SetFillColor(0);

   // hsumpedI[0].SetTitle("");  
   hsumpedD[0]->GetYaxis()->SetTitle("Counts");
   hsumpedD[0]->GetXaxis()->SetTitle("Double Average RPed (counts)");
   hsumpedD[0]->SetLineColor(2);
   hsumpedD[0]->Draw("h");  

   c1->cd(2);
   c1_2->SetBorderMode(0);
   c1_2->SetFillColor(0);

   // hsumpedI[0]->SetTitle("");  
   hsumpedDI[0]->GetYaxis()->SetTitle("Counts");
   hsumpedDI[0]->GetXaxis()->SetTitle("Double Average IPed (counts)");
   hsumpedDI[0]->SetLineColor(2);
   hsumpedDI[0]->Draw("h");  

   c1->cd(3);
   c1_3->SetBorderMode(0);
   c1_3->SetFillColor(0);

   // hsumpedI[0]->SetTitle("");  
   hsumpedI[0]->GetYaxis()->SetTitle("Counts");
   hsumpedI[0]->GetXaxis()->SetTitle("Integer Average IPed(counts)");
   hsumpedI[0]->SetLineColor(2);
   hsumpedI[0]->Draw("h"); 

   c1->cd(4);
   c1_4->SetBorderMode(0);
   c1_4->SetFillColor(0);

   // hsumpedD[1]->SetTitle("");  
   hsumpedD[1]->GetYaxis()->SetTitle("Counts");
   hsumpedD[1]->GetXaxis()->SetTitle("Double Average RPed (counts)");
   hsumpedD[1]->SetLineColor(2);
   hsumpedD[1]->Draw("h");  

   c1->cd(5);
   c1_5->SetBorderMode(0);
   c1_5->SetFillColor(0);

   // hsumpedDII1]->SetTitle("");  
   hsumpedDI[1]->GetYaxis()->SetTitle("Counts");
   hsumpedDI[1]->GetXaxis()->SetTitle("Double Average IPed(counts)");
   hsumpedDI[1]->SetLineColor(2);
   hsumpedDI[1]->Draw("h");  

   c1->cd(6);
   c1_6->SetBorderMode(0);
   c1_6->SetFillColor(0);

   // hsumpedDI1]->SetTitle("");  
   hsumpedI[1]->GetYaxis()->SetTitle("Counts");
   hsumpedI[1]->GetXaxis()->SetTitle("Integer Average IPed (counts)");
   hsumpedI[1]->SetLineColor(2);
   hsumpedI[1]->Draw("h");  


   TCanvas *c2 = new TCanvas("c2", "c2",200,10,1000,700);
   c2->SetBorderMode(0);
   c2->SetFillColor(0);

   c2->SetGridx();
   c2->SetGridy();
   c2->Divide(3,2);

   c2->cd(1);
   c2_1->SetBorderMode(0);
   c2_1->SetFillColor(0);

   // hsumpedI[2].SetTitle("");  
   hsumpedD[2]->GetYaxis()->SetTitle("Counts");
   hsumpedD[2]->GetXaxis()->SetTitle("Double Average RPed (counts)");
   hsumpedD[2]->SetLineColor(2);
   hsumpedD[2]->Draw("h");  

   c2->cd(2);
   c2_2->SetBorderMode(0);
   c2_2->SetFillColor(0);

   // hsumpedDI[2]->SetTitle("");  
   hsumpedDI[2]->GetYaxis()->SetTitle("Counts");
   hsumpedDI[2]->GetXaxis()->SetTitle("Double Average IPed (counts)");
   hsumpedDI[2]->SetLineColor(2);
   hsumpedDI[2]->Draw("h");  

   c2->cd(3);
   c2_3->SetBorderMode(0);
   c2_3->SetFillColor(0);

   // hsumpedI[2]->SetTitle("");  
   hsumpedI[2]->GetYaxis()->SetTitle("Counts");
   hsumpedI[2]->GetXaxis()->SetTitle("Integer Average IPed (counts)");
   hsumpedI[2]->SetLineColor(2);
   hsumpedI[2]->Draw("h"); 

   c2->cd(4);
   c2_4->SetBorderMode(0);
   c2_4->SetFillColor(0);

   // hsumpedD[3]->SetTitle("");  
   hsumpedD[3]->GetYaxis()->SetTitle("Counts");
   hsumpedD[3]->GetXaxis()->SetTitle("Double Average RPed (counts)");
   hsumpedD[3]->SetLineColor(2);
   hsumpedD[3]->Draw("h");  

   c2->cd(5);
   c2_5->SetBorderMode(0);
   c2_5->SetFillColor(0);

   // hsumpedDI[3]->SetTitle("");  
   hsumpedDI[3]->GetYaxis()->SetTitle("Counts");
   hsumpedDI[3]->GetXaxis()->SetTitle("Double Average IPed (counts)");
   hsumpedDI[3]->SetLineColor(2);
   hsumpedDI[3]->Draw("h");  

   c2->cd(6);
   c2_6->SetBorderMode(0);
   c2_6->SetFillColor(0);

   // hsumpedDI1]->SetTitle("");  
   hsumpedI[3]->GetYaxis()->SetTitle("Counts");
   hsumpedI[3]->GetXaxis()->SetTitle("Integer Average IPed (counts)");
   hsumpedI[3]->SetLineColor(2);
   hsumpedI[3]->Draw("h");  

   TCanvas *c3 = new TCanvas("c3", "c3",200,10,1000,700);
   c3->SetBorderMode(0);
   c3->SetFillColor(0);

   c3->SetGridx();
   c3->SetGridy();
   c3->Divide(3,2);

   c3->cd(1);
   c3_1->SetBorderMode(0);
   c3_1->SetFillColor(0);

   // hsumpedI[4].SetTitle("");  
   hsumpedD[4]->GetYaxis()->SetTitle("Counts");
   hsumpedD[4]->GetXaxis()->SetTitle("Double Average RPed (counts)");
   hsumpedD[4]->SetLineColor(2);
   hsumpedD[4]->Draw("h");  

   c3->cd(2);
   c3_2->SetBorderMode(0);
   c3_2->SetFillColor(0);

   // hsumpedDI[4]->SetTitle("");  
   hsumpedDI[4]->GetYaxis()->SetTitle("Counts");
   hsumpedDI[4]->GetXaxis()->SetTitle("Double Average IPed (counts)");
   hsumpedDI[4]->SetLineColor(2);
   hsumpedDI[4]->Draw("h");  

   c3->cd(3);
   c3_3->SetBorderMode(0);
   c3_3->SetFillColor(0);

   // hsumpedI[4]->SetTitle("");  
   hsumpedI[4]->GetYaxis()->SetTitle("Counts");
   hsumpedI[4]->GetXaxis()->SetTitle("Integer Average IPed (counts)");
   hsumpedI[4]->SetLineColor(2);
   hsumpedI[4]->Draw("h"); 

   c3->cd(4);
   c3_4->SetBorderMode(0);
   c3_4->SetFillColor(0);

   // hsumpedD[5]->SetTitle("");  
   hsumpedD[5]->GetYaxis()->SetTitle("Counts");
   hsumpedD[5]->GetXaxis()->SetTitle("Double Average RPed (counts)");
   hsumpedD[5]->SetLineColor(2);
   hsumpedD[5]->Draw("h");  

   c3->cd(5);
   c3_5->SetBorderMode(0);
   c3_5->SetFillColor(0);

   // hsumpedDI[5]->SetTitle("");  
   hsumpedDI[5]->GetYaxis()->SetTitle("Counts");
   hsumpedDI[5]->GetXaxis()->SetTitle("Double Average IPed (counts)");
   hsumpedDI[5]->SetLineColor(2);
   hsumpedDI[5]->Draw("h");  

   c3->cd(6);
   c3_6->SetBorderMode(0);
   c3_6->SetFillColor(0);

   // hsumpedDI1]->SetTitle("");  
   hsumpedI[5]->GetYaxis()->SetTitle("Counts");
   hsumpedI[5]->GetXaxis()->SetTitle("Integer Average IPed  (counts)");
   hsumpedI[5]->SetLineColor(2);
   hsumpedI[5]->Draw("h");  
 

   TCanvas *c4 = new TCanvas("c4", "c4",200,10,1000,700);
   c4->SetBorderMode(0);
   c4->SetFillColor(0);

   c4->SetGridx();
   c4->SetGridy();
   c4->Divide(3,2);

   c4->cd(1);
   c4_1->SetBorderMode(0);
   c4_1->SetFillColor(0);

   // hsumpedI[6].SetTitle("");  
   hsumpedD[6]->GetYaxis()->SetTitle("Counts");
   hsumpedD[6]->GetXaxis()->SetTitle("Double Average RPed (counts)");
   hsumpedD[6]->SetLineColor(2);
   hsumpedD[6]->Draw("h");  

   c4->cd(2);
   c4_2->SetBorderMode(0);
   c4_2->SetFillColor(0);

   // hsumpedDI[6]->SetTitle("");  
   hsumpedDI[6]->GetYaxis()->SetTitle("Counts");
   hsumpedDI[6]->GetXaxis()->SetTitle("Double Average IPed (counts)");
   hsumpedDI[6]->SetLineColor(2);
   hsumpedDI[6]->Draw("h"); 

   c4->cd(3);
   c4_3->SetBorderMode(0);
   c4_3->SetFillColor(0);

   // hsumpedI[6]->SetTitle("");  
   hsumpedI[6]->GetYaxis()->SetTitle("Counts");
   hsumpedI[6]->GetXaxis()->SetTitle("Integer Average IPed (counts)");
   hsumpedI[6]->SetLineColor(2);
   hsumpedI[6]->Draw("h");   

   c4->cd(4);
   c4_4->SetBorderMode(0);
   c4_4->SetFillColor(0);

   // hsumpedD[7]->SetTitle("");  
   hsumpedD[7]->GetYaxis()->SetTitle("Counts");
   hsumpedD[7]->GetXaxis()->SetTitle("Double Average RPed (counts)");
   hsumpedD[7]->SetLineColor(2);
   hsumpedD[7]->Draw("h");  

   c4->cd(5);
   c4_5->SetBorderMode(0);
   c4_5->SetFillColor(0);

   // hsumpedDI1]->SetTitle("");  
   hsumpedDI[7]->GetYaxis()->SetTitle("Counts");
   hsumpedDI[7]->GetXaxis()->SetTitle("Double Average IPed (counts)");
   hsumpedDI[7]->SetLineColor(2);
   hsumpedDI[7]->Draw("h");   

   c4->cd(6);
   c4_6->SetBorderMode(0);
   c4_6->SetFillColor(0);

   // hsumpedDI1]->SetTitle("");  
   hsumpedI[7]->GetYaxis()->SetTitle("Counts");
   hsumpedI[7]->GetXaxis()->SetTitle("Integer Average IPed (counts)");
   hsumpedI[7]->SetLineColor(2);
   hsumpedI[7]->Draw("h");  


   // extract mean and sigmas from histograms


  Double_t avepedD[9];
  Double_t avepedDI[9];
  Double_t avepedI[9];

  Double_t aveped_errD[9];
  Double_t aveped_errDI[9];
  Double_t aveped_errI[9];

  Double_t sigmapedD[9];
  Double_t sigmapedDI[9];
  Double_t sigmapedI[9];

  Double_t sigmaped_errD[9];
  Double_t sigmaped_errDI[9];
  Double_t sigmaped_errI[9];

  Double_t ndxvec[9];
  Double_t ndxvec_err[9];

  Int_t nhist=0;

  for (int j=1; j<nbins; j*=2) {

      // cout << " j=" << j << " nhist=" << nhist << endl;

     ndxvec[nhist] = j;
     ndxvec_err[nhist] = 0.01;

      avepedD[nhist] = hsumpedD[nhist]->GetMean();
      aveped_errD[nhist] = hsumpedD[nhist]->GetMeanError();
      sigmapedD[nhist] = hsumpedD[nhist]->GetRMS();
      sigmaped_errD[nhist] = hsumpedD[nhist]->GetRMSError();

      avepedDI[nhist] = hsumpedDI[nhist]->GetMean();
      aveped_errDI[nhist] = hsumpedDI[nhist]->GetMeanError();
      sigmapedDI[nhist] = hsumpedDI[nhist]->GetRMS();
      sigmaped_errDI[nhist] = hsumpedDI[nhist]->GetRMSError();

      avepedI[nhist] = hsumpedI[nhist]->GetMean();
      aveped_errI[nhist] = hsumpedI[nhist]->GetMeanError();
      sigmapedI[nhist] = hsumpedI[nhist]->GetRMS();
      sigmaped_errI[nhist] = hsumpedI[nhist]->GetRMSError();

    nhist++;
    }

  // TGraphErrors *gavepedD = new TGraphErrors(9,ndxvec,avepedD,ndxvec_err,aveped_errD);
  TGraph *gavepedD = new TGraph(9,ndxvec,avepedD);
  TGraph *gsigmapedD = new TGraph(9,ndxvec,sigmapedD);
  TGraph *gavepedDI = new TGraph(9,ndxvec,avepedDI);
  TGraph *gsigmapedDI = new TGraph(9,ndxvec,sigmapedDI);
  TGraph *gavepedI = new TGraph(9,ndxvec,avepedI);
  TGraph *gsigmapedI = new TGraph(9,ndxvec,sigmapedI);

   // define 1/sqrt function


  TF1 *sqrt = new TF1("sqrt",sqrt_func,1,1000,1);
  sqrt->SetParName(0,"sigma1");
  sqrt->SetParameter(0,avepedD[0]);
  
   TCanvas *c5 = new TCanvas("c5", "c5",200,10,700,700);
   c5->Divide(1,2);

   c5->cd(1);
   c5_1->SetGridx();
   c5_1->SetGridy();
   c5_1->SetLogx();

   gsigmapedD->SetTitle(""); 
   gsigmapedD->GetYaxis()->SetTitle("Sigma");
   gsigmapedD->GetXaxis()->SetTitle("Entries in sum");
   gsigmapedD->SetMarkerStyle(20);
   gsigmapedD->SetMarkerColor(1);
   sqrt->SetParameter(0,avepedD[0]);
   sqrt->SetLineColor(1);
   gsigmapedD->Fit("sqrt");
   gsigmapedD->Draw("Ap");  

   gsigmapedDI->SetMarkerStyle(20);
   gsigmapedDI->SetMarkerColor(2);
   sqrt->SetParameter(0,avepedDI[0]);
   sqrt->SetLineColor(2);
   // gsigmapedDI->Fit("sqrt");
   gsigmapedDI->Draw("samep");  

   gsigmapedI->SetMarkerStyle(20);
   gsigmapedI->SetMarkerColor(4);
   sqrt->SetParameter(0,avepedDI[0]);
   sqrt->SetLineColor(4);
   // gsigmapedI->Fit("sqrt");
   gsigmapedI->Draw("samep");  

   c5->cd(2);
   c5_2->SetGridx();
   c5_2->SetGridy();
   c5_2->SetLogx();
   c5_2->SetBorderMode(0);
   c5_2->SetFillColor(0);

  Double_t ymin=98;
  Double_t ymax=102;

   gavepedD->SetTitle(""); 
   gavepedD->GetYaxis()->SetTitle("Mean");
   gavepedD->GetXaxis()->SetTitle("Entries in sum");
   gavepedD->GetYaxis()->SetRangeUser(ymin,ymax);
   gavepedD->SetMarkerStyle(20);
   gavepedD->SetMarkerColor(1);
   gavepedD->Draw("Ap");  

   gavepedDI->SetMarkerStyle(20);
   gavepedDI->SetMarkerColor(2);
   gavepedDI->Draw("samep");   

   gavepedI->SetMarkerStyle(20);
   gavepedI->SetMarkerColor(4);
   gavepedI->Draw("samep");  

   // ouput legend

  TLegend *leg = new TLegend(0.4,0.7,0.7,0.9);
   leg->AddEntry(gavepedD,"Real Ave, Real Ped","p");
   leg->AddEntry(gavepedDI,"Real Ave, Int Ped","p");
   leg->AddEntry(gavepedI,"Int Ave, Int Ped","p");
   leg->Draw();


   sprintf (string,"Mean Pedestal=%.2f\n",pedestal);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.92,string);
   t1->SetNDC();
   t1->SetTextSize(0.07);
   t1->Draw();

   sprintf (string,"pedestal_studies_%04d.pdf(",pedestal*10);
   c0->SaveAs(string);
   sprintf (string,"pedestal_studies_%04d.pdf",pedestal*10);
   c1->SaveAs(string);
   sprintf (string,"pedestal_studies_%04d.pdf",pedestal*10);
   c2->SaveAs(string);
   sprintf (string,"pedestal_studies_%04d.pdf",pedestal*10);
   c3->SaveAs(string);
   sprintf (string,"pedestal_studies_%04d.pdf",pedestal*10);
   c4->SaveAs(string);
   sprintf (string,"pedestal_studies_%04d.pdf)",pedestal*10);
   c5->SaveAs(string);
   sprintf (string,"pedestal_studies_%04d.png",pedestal*10);
   c5->SaveAs(string);



}
Double_t sqrt_func (Double_t *x, Double_t *par) {

   Double_t sigma1=par[0];
   Double_t x1=x[0];

   char string[256];

   Double_t func;

   func = sigma1/sqrt(x1);

   return func;
}