#include <vector>

void pedestal_studies_correlations(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 arrays to store distribution parameters

  const int Nsamples=10;

  Double_t avepedD[Nsamples];
  Double_t aveped_errD[Nsamples];

  Double_t sigmapedD[Nsamples];
  Double_t sigmaped_errD[Nsamples];

  Double_t ndxvec[Nsamples];
  Double_t ndxvec_err[Nsamples];

  const int maxbins=1024;

  Double_t rho[maxbins][maxbins];

  Double_t sum=0;
       
   const int Ncorr=4;
   Double_t delta_cor=0.05;

   TGraph *gsigmapedD[Ncorr];

  for (int jc=0; jc<Ncorr; jc++)  {

  Double_t correlation = jc*delta_cor;
  Int_t ndx=0;

  for (int j=1; j<nbins; j*=2) {
 
      // averaging over "j" samples:

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

      // fill correlation matrix for that size

      sum = 0;

     for (int k=0; k<j; k++) {
            for (int kk=0; kk<j; kk++) {

                   if (k == kk) {
                       rho[k][kk] = ped_sigma*ped_sigma;
                       }
                   else {
                      rho[k][kk] = ped_sigma*ped_sigma * correlation;
                   }
                   sum +=   rho[k][kk];    
              }

    }
    sigmapedD[ndx] = sum > 0? sqrt( sum) /j : 0;
     ndxvec[ndx] = j;
     ndxvec_err[ndx] = 0.01;

    cout << " jc=" << jc << " correlation=" << correlation << " ndx=" << ndx << " sigmapedD=" << sigmapedD[ndx] << endl;
    ndx++;
}
    gsigmapedD[jc] = new TGraph(Nsamples,ndxvec,sigmapedD);

}

   // define 1/sqrt function


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

   // ouput legend

   TLegend *leg = new TLegend(0.25,0.7,0.55,0.9);

   for (jc=0; jc<Ncorr; jc++) {

   Double_t correlation = jc*delta_cor;
   gsigmapedD[jc]->SetTitle(""); 
   gsigmapedD[jc]->GetYaxis()->SetTitle("Sigma");
   gsigmapedD[jc]->GetXaxis()->SetTitle("Entries in sum");
   gsigmapedD[jc]->SetMarkerStyle(20+jc);
   gsigmapedD[jc]->SetMarkerColor(4);

   sprintf (string,"Correlation=%.2f\n",correlation);
   printf("string=%s",string);
   leg->AddEntry(gsigmapedD[jc],string,"p");

   if (jc == 0 ) {
        gsigmapedD[jc]->Fit("sqrt");
        sqrt->SetParameter(0,sigmapedD[0]);
        sqrt->SetLineColor(1);
        gsigmapedD[jc]->Draw("Ap");  
        }
  else {
        gsigmapedD[jc]->Draw("psame"); 
        }

}    // end loop over correlations


   leg->Draw();


   sprintf (string,"Effect of correlations on average \n");
   printf("string=%s",string);
   t1 = new TLatex(0.15,0.92,string);
   t1->SetNDC();
   t1->SetTextSize(0.03);
   t1->Draw();

   sprintf (string,"pedestal_studies_correlations_%04d.pdf)",pedestal*10);
   c5->SaveAs(string);
   sprintf (string,"pedestal_studies_correlations_%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;
}