#include <inttypes.h>
#include <iostream>
#include <fstream>
#include <utility>
#include <vector>
#include <map>
#include <stdio.h>
#include <stdlib.h>
#include <string>
#include <cstring>
#include <sstream>

#include "TLatex.h"
#include "TPaveStats.h"
#include "TGraphPainter.h"
#include "TString.h"
#include "TCollection.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TH2D.h"
#include "TH1F.h"
#include "TF1.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TMinuit.h"
#include "TKey.h"
#include "TDatime.h"
#include "TAxis.h"
#include "TLine.h"
#include "TTree.h"
#include "TBranch.h"
#include "TStyle.h"

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

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

   Double_t func;

   func = sigma1/sqrt(x1);
   //  cout << "x1=" << x1 << " sigma1=" << sigma1 << " func=" << func << endl;

   return func;
}

void pedestal_studiesA_sigma(Int_t crate_select=37, Int_t slot_select=3, Int_t channel_select=0)
{

  // pedestal_studies script modified to compute rms values from data

gStyle->SetPalette(1,0);
gStyle->SetOptStat(kTRUE);
// gStyle->SetOptStat(kFALSE);
// gStyle->SetOptFit(kTRUE);
gStyle->SetOptFit(kFALSE);
// 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);
//

  // Read input file generated by DAQ and converted to ROOT using DAQTree plugin
    
  // TString fileName = "bcal_roc5_14.0.root";
  
  TFile* f = new TFile("bcal_s4_3281.0.root");
  //TTree* tr = (TTree*)f->Get("channeltree");
  TTree* tr = (TTree*)f->Get("Df250WindowRawData");
  
  uint32_t channelnum;         /// Arbitrary global channel number
  uint32_t eventnum;           /// Event number
  uint32_t rocid;              /// Crate number
  uint32_t slot;               /// Slot number in crate
  uint32_t channel;            /// Channel number in slot
  uint32_t itrigger; ///trigger number for cases when this hit was read in a multi-event block (DDAQAddress)
  uint32_t nsamples;           /// Number of samples in the event
  std::vector<uint32_t>* waveform = 0;   /// STL vector of samples in the waveform for the event
  uint32_t w_integral;         /// Sum of all samples in the waveform
  uint32_t w_min;              /// Minimum sample in the waveform
  uint32_t w_max;              /// Maximum sample in the waveform
  uint32_t w_samp1;            /// First sample in the waveform
  
  tr->SetBranchAddress("channelnum",&channelnum);
  tr->SetBranchAddress("eventnum",&eventnum);
  tr->SetBranchAddress("rocid",&rocid);
  tr->SetBranchAddress("slot",&slot);
  tr->SetBranchAddress("channel",&channel);
  tr->SetBranchAddress("itrigger",&itrigger);
  tr->SetBranchAddress("nsamples",&nsamples);
  tr->SetBranchAddress("waveform",&waveform);
  tr->SetBranchAddress("w_integral",&w_integral);
  tr->SetBranchAddress("w_min",&w_min);
  tr->SetBranchAddress("w_max",&w_max);
  tr->SetBranchAddress("w_samp1",&w_samp1);
  
// Combos is 1, 2, 5, 10
//  int combo[] = {1,2,5,10,50};
//  const int _crate = 6, _slot = 20, _ch = 16, _co = sizeof(combo)/sizeof(*combo), _samp = 100;
//  TH1F* histPed[_crate][_slot][_ch][_co];
//  TH1F* histSample[_crate][_slot][_ch][_samp];

  Int_t nbins=100;
  
  // Entries over events
  const int max = 100;

  /*Int_t crate_select=37;
  Int_t slot_select=3;
  Int_t channel_select=8;*/
  Int_t crate = crate_select;

  Int_t nsize = (Int_t) tr->GetEntries();
  Int_t maxevents = 768*max > nsize? nsize : 768*max;

  const int _samp = 100;

   char string[256];

   // Double_t pedestal=100;
   // Double_t ped_sigma=1.2;

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

  TH1I *hwindowI = new TH1I("hwindowI","Readout Window", nbins,xmin,xmax);     

  TH1D *hsumpedDI[9];
  Int_t nhist=0;

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

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

       // define histogram to keep results

      sprintf (string,"SumpedDI_j%03d_bins\n",j);
      printf("string=%s",string);
      hsumpedDI[nhist] = new TH1D(string,string,40,90,110);
      hsumpedDI[nhist]->Sumw2();
      nhist++;

  }

  for (Int_t ievent=0; ievent<maxevents; ievent++){
    tr->GetEntry(ievent);
    // File is looped as follows:
    // in order of event go through all channels in a slot and all slot in a crate and all crates
 
    if (ievent%1000==0)  cout << " event counter=" << ievent << endl;
    

    if (crate != rocid) continue;
    if (slot != slot_select) continue;
    if (channel != channel_select) continue;

    hwindowI->Reset();

    for (Int_t jj=0; jj<_samp; jj++){
      
      // Get adc value
      int adc_value = (*waveform)[jj];
      hwindowI->Fill(jj,adc_value);
    }
  
    /*TCanvas *c0 = new TCanvas("c0", "c0",200,10,700,700);
   c0->SetBorderMode(0);
   c0->SetFillColor(0);

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

   hwindowI->SetTitle("");
   hwindowI->GetYaxis()->SetTitle("Pedestal");
   hwindowI->GetXaxis()->SetTitle("FADC samples");
   hwindowI->SetLineColor(2);
   hwindowI->Draw();*/
    nhist = 0;

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

      cout << " j=" << j << " nhist=" << nhist << endl;
      Double_t sumpedDI=0;

      for (int jj=1; jj<j+1; jj++) {
         sumpedDI += hwindowI->GetBinContent(jj);
         // cout << "j=" << j << " jj=" << jj << " content=" << hwindowI->GetBinContent(jj) << endl;               
      }
      sumpedDI /= j;
      cout << " j=" << j << " sumpedDI=" << sumpedDI << endl;   
      hsumpedDI[nhist]->Fill(sumpedDI);
      nhist++;

  } 

  }  // end loop over events in tree
      

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

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

   for (Int_t jj=0; jj<nhist; jj++) {

     // cout << "jj=" << jj << " nhist=" << nhist << endl;    
      c1->cd(jj+1);
      hsumpedDI[jj]->GetYaxis()->SetTitle("Counts");
      hsumpedDI[jj]->GetXaxis()->SetTitle("Double Average IPed (counts)");
      hsumpedDI[jj]->SetLineColor(2);
      hsumpedDI[jj]->Draw("h");
      }


   // extract mean and sigmas from histograms


  Double_t avepedDI[9];

  Double_t aveped_errDI[9];
  Double_t sigmapedDI[9];

  Double_t sigmaped_errDI[9];

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

  nhist=0;

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

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

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

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

    nhist++;
    }

  TGraph *gavepedDI = new TGraph(7,ndxvec,avepedDI);
  TGraph *gsigmapedDI = new TGraph(7,ndxvec,sigmapedDI);

   // define 1/sqrt function


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

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

   gsigmapedDI->SetMarkerStyle(20);
   gsigmapedDI->SetMarkerColor(2);
   gsigmapedDI->GetYaxis()->SetRangeUser(0,1.3);
   gsigmapedDI->SetTitle("R3281, averages over samples in the same event");
   gsigmapedDI->GetYaxis()->SetTitle("Data RMS(counts)");
   gsigmapedDI->GetXaxis()->SetTitle("Number of entries in average");
   // gsigmapedDI->Fit("sqrt");
   gsigmapedDI->Draw("Ap");
   sqrt->SetLineColor(1);
   sqrt->SetLineWidth(2);
   sqrt->Draw("samel"); 

   c5->cd(2);

  Double_t ymin=95;
  Double_t ymax=101;

   gavepedDI->SetMarkerStyle(20);
   gavepedDI->SetMarkerColor(2);
   gavepedDI->GetYaxis()->SetRangeUser(ymin,ymax);
   gavepedDI->SetTitle("R3281, averages over samples in the same event");
   gavepedDI->GetYaxis()->SetTitle("Mean (counts)");
   gavepedDI->GetXaxis()->SetTitle("Number of entries in average");
   gavepedDI->Draw("Ap");

   // ouput legend

   /*TLegend *leg = new TLegend(0.4,0.7,0.7,0.9);
   leg->AddEntry(gavepedDI,"Real 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_studiesA_sigma.pdf(");
     c0->SaveAs(string);*/
   sprintf (string,"pedestal_studiesA_sigma.pdf(");
   c1->SaveAs(string);
   sprintf (string,"pedestal_studiesA_sigma.pdf)");
   c5->SaveAs(string);



}