#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"

void pedestal_studiesA(Int_t crate_select, Int_t slot_select, Int_t channel_select) {

gStyle->SetPalette(1,0);
// gStyle->SetOptStat(kTRUE);
gStyle->SetOptStat(kFALSE);
// gStyle->SetOptFit(kTRUE);
gStyle->SetOptFit(kFALSE);
// gStyle->SetOptFit(11111);
// gStyle->SetOptStat(111111);
gStyle->SetPadRightMargin(0.15);
gStyle->SetPadLeftMargin(0.15);
gStyle->SetPadBottomMargin(0.15);
// gStyle->SetFillColor(0);
  

   char fname[132];

  // 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;

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

  const int _crate = 1, _slot = 1, _ch = 16, _samp = 100;
  char histName[132];
  TH1F* histPed[_ch];


  // get tree

  tr->GetEntry(0);
  
  // Entries over events
  const int max = 100;

  // define arrays

  Double_t Sx2[max];
  Double_t Sx[max];
  Double_t S1[max];
  Double_t sigma[max];;
  Double_t Sxy[max][max];
  Double_t Sxy1[max][max];
  Int_t ndx=0;

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

  printf ("\nMumber of events to process=%d\n\n",maxevents);


    // initialize variables;
  for (Int_t j=0; j<_samp; j++) {
    Sx2[j] = 0;
    Sx[j] = 0;
    S1[j] = 0;
    for (Int_t k=0; k<_samp; k++) {
      Sxy[j][k] = 0;
      Sxy1[j][k] = 0;
    }
  }

  TH2D *Rhoxy = new TH2D("Rhoxy","Matrix of correllations",nbins,0,nbins,nbins,0,nbins);
  
  
  // Loop Over Tree
  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;

    int wsize = waveform->size(); // 100 samples
    // cout << "wsize=" << wsize << endl;
    
    // ******************
    // Loop Over Waveform
    // ******************
    // cout << endl;
    for (Int_t jj=0; jj<_samp; jj++){
      
      // Get adc value
      int adc_value = (*waveform)[jj];

      // cout << "ndx=" << ndx << " event=" << ievent << " sample=" << jj << " adc sample=" << adc_value << endl;
      
      // update arrays
      Sx2[jj] += adc_value*adc_value;
      Sx[jj] += adc_value;
      S1[jj] += 1;

      // second loop over waveform for covariance matrix
      for (Int_t k=0; k<_samp; k++){
      
	// Get adc value
	int adc_value2 = (*waveform)[k];

	// cout << "ndx=" << ndx << " event=" << ievent << " sample=" << k << " adc sample=" << adc_value << endl;
      
        Sxy[jj][k] += adc_value*adc_value2;
        Sxy1[jj][k] += 1;
        // Inner waveform loop
      }

    } // End Loop Over Waveform
    ndx++;
    
  } // End Loop Over Tree

  // compute sigmas

  for (Int_t k=0; k<_samp; k++) {
    Double_t sigma2 = Sx2[k]/S1[k] - (Sx[k]/S1[k])*(Sx[k]/S1[k]);
    sigma[k] = sigma2 > 0? sqrt(sigma2) : 0;
    // cout << " k=" << k << " Sx2=" << Sx2[k] << " Sx=" << Sx[k] << " S1=" << S1[k] << " sigma[k]=" << sigma[k] << endl;
  }

  for (Int_t j=0; j<_samp; j++) {
    for (Int_t k=0; k<_samp; k++) {
      Double_t rho = Sxy[j][k]/Sxy1[j][k] -(Sx[j]/S1[j])*(Sx[k]/S1[k]) ;
      rho /= sigma[j]*sigma[k];
      // cout << "j=" << j << " k=" << k << " rho=" << rho << " sigmaj=" << sigma[j] << " sigmak=" << sigma[k] << " Sxy1=" << Sxy1[j][k] <<   endl;
      Rhoxy->Fill(j,k,rho);
    }
  } 

  // compute average

  Double_t rho_sum=0;
  Double_t rho_1=0;
  
  for (Int_t j=0; j<_samp; j++) {
    for (Int_t k=j+1; k<_samp; k++) {
      Double_t rho = Sxy[j][k]/Sxy1[j][k] -(Sx[j]/S1[j])*(Sx[k]/S1[k]) ;
      rho /= sigma[j]*sigma[k];
      rho_sum += rho;
      rho_1 += 1;
    }
  } 

  Double_t rho_mean = rho_sum/rho_1;
  cout << " rho_mean=" << rho_mean << endl;
  
  // plot histograms

   TCanvas *c1 = new TCanvas("c1","c1 pedestal_studiesA",200,10,700,700);
   c1->SetBorderMode(0);
   c1->SetFillColor(0);
   //c1->Divide(2,2);
   // c1->SetGridx();
   // c1->SetGridy();
   // c1->SetLogz();

   Float_t zmin=-0.1;
   Float_t zmax=0.1;

   //c1->cd(1);

   sprintf(fname,"Rhoxy %02d %02d %02d, mean=%.4f\n",crate_select,slot_select,channel_select,rho_mean);
   Rhoxy->SetTitle(fname);
   Rhoxy->GetZaxis()->SetRangeUser(zmin,zmax);
   Rhoxy->GetYaxis()->SetTitle("Sample Number");
   Rhoxy->GetXaxis()->SetTitle("Sample Number");
   Rhoxy->Draw("colz");

     /*sprintf (string,"#sigma=%.2f\n",sigma);
     printf("string=%s",string);
     t1 = new TLatex(0.4,0.5,string);
     t1->SetNDC();
     t1->SetTextSize(0.07);
     t1->Draw();*/


   sprintf(fname,"pedestal_studiesA_%02d_%02d_%02d_c1.pdf",crate_select,slot_select,channel_select);
   c1->SaveAs(fname);
   sprintf(fname,"pedestal_studiesA_%02d_%02d_%02d_c1.png",crate_select,slot_select,channel_select);
   c1->SaveAs(fname);


}