#include <vector>

void pedestal_studies_samples(Double_t Nsamples=4)
{

// histogram average values for a fixed number of samples

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;

   Double_t pedestal=99.4;
   Double_t ped_sigma=1.2;
   Int_t nbins=1024;

   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);
  
   TRandom3 *r1 = new TRandom3();
   UInt_t  iseed = 23424111513;
   r1->SetSeed(iseed);


  const int Npeds=11;
  TH1D *hsumpedD[Npeds];
  TH1D *hsumpedDI[Npeds];
  TH1I *hsumpedI[Npeds];

  for (int nhist=0; nhist<Npeds; nhist++) {

     pedestal += 0.1;

     // reset and then fill pedestal windows

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

      Int_t niter = nbins/Nsamples;
      cout << " niter=" << niter << " pedestal=" << pedestal << " nhist=" << nhist << endl;

       // define histogram to keep results

       sprintf (string,"SumpedD_ped=%04d\n",pedestal*10);
      // printf("string=%s",string);
      hsumpedD[nhist] = new TH1D(string,string,500,95,105);
      hsumpedD[nhist]->Sumw2();

      sprintf (string,"SumpedDI_ped=%04d\n",pedestal*10);
      // printf("string=%s",string);
      hsumpedDI[nhist] = new TH1D(string,string,500,95,105);
      hsumpedDI[nhist]->Sumw2();

      sprintf (string,"SumpedI_ped=%04d\n",pedestal*10);
      // 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<Nsamples;k++) {
                    Int_t ndx = jj*Nsamples + k;
                    sumpedD += hwindowD->GetBinContent(ndx+1);
                    sumpedDI += hwindowI->GetBinContent(ndx+1);
                    sumpedI += hwindowI->GetBinContent(ndx+1);
	}
              sumpedD /= Nsamples;
              sumpedDI /= Nsamples;
              sumpedI /= Nsamples;
              hsumpedD[nhist]->Fill(sumpedD);
              hsumpedDI[nhist]->Fill(sumpedDI);
              hsumpedI[nhist]->Fill(sumpedI);
       }
  }
      

   // extract mean and sigmas from histograms


  TH1D *Ped4D = new TH1D("Ped4D","Pedestal mean D for 4 Samples",500, 95,105);
  TH1D *Ped4DI = new TH1D("Ped4DI","Pedestal mean DI for 4 Samples",500, 95,105);
  TH1D *Ped4I = new TH1D("Ped4I","Pedestal mean I for 4 Samples",500, 95,105);

  TH1D *DeltaPed4D = new TH1D("DeltaPed4D","Delta Pedestal D for 4 Samples",500, -5,5);
  TH1D *DeltaPed4DI = new TH1D("DeltaPed4DI","Delta Pedestal DI for 4 Samples",500, -5,5);
  TH1D *DeltaPed4I = new TH1D("DeltaPed4I","Delta Pedestal I for 4 Samples",500, -5,5);

  TH1D *Sigma4D = new TH1D("Sigma4D","Sigma mean D for 4 Samples",200, 0,2);
  TH1D *Sigma4DI = new TH1D("Sigma4DI","Sigma mean DI for 4 Samples",200, 0,2);
  TH1D *Sigma4I = new TH1D("Sigma4I","Sigma mean I for 4 Samples",200, 0,2);

  pedestal = 99.4;

  for (int nhist=0; nhist<Npeds; nhist++) {

      pedestal += 0.1;

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

      Ped4D->Fill(avepedD);
      DeltaPed4D->Fill(avepedD-pedestal);
      Sigma4D->Fill(sigmapedD);

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

      Ped4DI->Fill(avepedDI);
      DeltaPed4DI->Fill(avepedDI-pedestal);
      Sigma4DI->Fill(sigmapedDI);

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

      Ped4I->Fill(avepedI);
      DeltaPed4I->Fill(avepedI-pedestal);
      Sigma4I->Fill(sigmapedI);

    }

  
   TCanvas *c1 = new TCanvas("c1", "c1",200,10,700,700);
   c1->Divide(2,2);

   c1->cd(1);
   c1_1->SetGridx();
   c1_1->SetGridy();

   Ped4D->SetTitle(""); 
   Ped4D->GetYaxis()->SetTitle("Counts");
   Ped4D->GetXaxis()->SetTitle("Mean Pedestal");
   Ped4D->SetMarkerStyle(20);
   Ped4D->SetMarkerColor(1);
   Ped4D->Draw("");  

   sprintf (string,"Nsamples=%d\n",Nsamples);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.92,string);
   t1->SetNDC();
   t1->SetTextSize(0.07);
   t1->Draw();


   c1->cd(2);
   c1_2->SetGridx();
   c1_2->SetGridy();

   Ped4DI->SetTitle(""); 
   Ped4DI->GetYaxis()->SetTitle("Counts");
   Ped4DI->GetXaxis()->SetTitle("Mean Pedestal");
   Ped4DI->SetMarkerStyle(20);
   Ped4DI->SetMarkerColor(1);
   Ped4DI->Draw("");  

   c1->cd(3);
   c1_3->SetGridx();
   c1_3->SetGridy();

   Ped4I->SetTitle(""); 
   Ped4I->GetYaxis()->SetTitle("Counts");
   Ped4I->GetXaxis()->SetTitle("Mean Pedestal");
   Ped4I->SetMarkerStyle(20);
   Ped4I->SetMarkerColor(1);
   Ped4I->Draw("");  


   TCanvas *c2 = new TCanvas("c2", "c2",200,10,700,700);
   c2->Divide(2,2);

   c2->cd(1);
   c2_1->SetGridx();
   c2_1->SetGridy();

   Sigma4D->SetTitle(""); 
   Sigma4D->GetYaxis()->SetTitle("Counts");
   Sigma4D->GetXaxis()->SetTitle("Mean Sigma");
   Sigma4D->SetMarkerStyle(20);
   Sigma4D->SetMarkerColor(1);
   Sigma4D->Draw("");   

   sprintf (string,"Nsamples=%d\n",Nsamples);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.92,string);
   t1->SetNDC();
   t1->SetTextSize(0.07);
   t1->Draw();

   c2->cd(2);
   c2_2->SetGridx();
   c2_2->SetGridy();

   Sigma4DI->SetTitle(""); 
   Sigma4DI->GetYaxis()->SetTitle("Counts");
   Sigma4DI->GetXaxis()->SetTitle("Mean Sigma");
   Sigma4DI->SetMarkerStyle(20);
   Sigma4DI->SetMarkerColor(1);
   Sigma4DI->Draw("");  

   c2->cd(3);
   c2_3->SetGridx();
   c2_3->SetGridy();

   Sigma4I->SetTitle(""); 
   Sigma4I->GetYaxis()->SetTitle("Counts");
   Sigma4I->GetXaxis()->SetTitle("Mean Sigma");
   Sigma4I->SetMarkerStyle(20);
   Sigma4I->SetMarkerColor(1);
   Sigma4I->Draw(""); 


   TCanvas *c3 = new TCanvas("c3", "c3",200,10,700,700);
   c3->Divide(2,2);

   c3->cd(1);
   c3_1->SetGridx();
   c3_1->SetGridy();

   DeltaPed4D->SetTitle(""); 
   DeltaPed4D->GetYaxis()->SetTitle("Counts");
   DeltaPed4D->GetXaxis()->SetTitle("Mean DeltaPedestal");
   DeltaPed4D->SetMarkerStyle(20);
   DeltaPed4D->SetMarkerColor(1);
   DeltaPed4D->Draw("");  

   sprintf (string,"Nsamples=%d\n",Nsamples);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.92,string);
   t1->SetNDC();
   t1->SetTextSize(0.07);
   t1->Draw();


   c3->cd(2);
   c3_2->SetGridx();
   c3_2->SetGridy();

   DeltaPed4DI->SetTitle(""); 
   DeltaPed4DI->GetYaxis()->SetTitle("Counts");
   DeltaPed4DI->GetXaxis()->SetTitle("Mean DeltaPedestal");
   DeltaPed4DI->SetMarkerStyle(20);
   DeltaPed4DI->SetMarkerColor(1);
   DeltaPed4DI->Draw("");  

   c3->cd(3);
   c3_3->SetGridx();
   c3_3->SetGridy();

   DeltaPed4I->SetTitle(""); 
   DeltaPed4I->GetYaxis()->SetTitle("Counts");
   DeltaPed4I->GetXaxis()->SetTitle("Mean DeltaPedestal");
   DeltaPed4I->SetMarkerStyle(20);
   DeltaPed4I->SetMarkerColor(1);
   DeltaPed4I->Draw("");  


   // 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,"pedestal_studies_%03dsamples.pdf(",Nsamples);
   c1->SaveAs(string);
   sprintf (string,"pedestal_studies_%03dsamples.pdf",Nsamples);
   c3->SaveAs(string);
   sprintf (string,"pedestal_studies_%03dsamples.pdf)",Nsamples);
   c2->SaveAs(string);


}