#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"
#include "TProfile.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_sigma5(Int_t crate_select=37, Int_t slot_select=3, Int_t channel_select=0, Int_t sample_pick=0)
{

  // pedestal_studies script modified to compute max and min rms from data 

gStyle->SetPalette(1,0);
gStyle->SetOptStat(kTRUE);;
gStyle->SetOptStat(111111);
// 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
  
  TFile* f = new TFile("bcal_s4_3281.0.root");
  // TFile* f = new TFile("hd_rawdata_003620_000.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 = 10000;
  const int max = 10000;
  const int maxids = 768;

  /*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 = maxids*max > nsize? nsize : maxids*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;

   // open list output file

   ofstream outlist;
   outlist.open ("pedestal_studiesA_sigma5.list");
             

   TH1D *hped[maxids];
   TH1D *hpedrel[maxids];
   TH2D *h2ped = new TH2D("h2ped","R3281 RMS Pedestals all",48,0,48,16,0,16);
   TH2D *h2pedrel = new TH2D("h2pedrel","R3281 RMS Pedestal-Window Ave",48,0,48,16,0,16);
   TH1D *h1ped = new TH1D("h1ped","R3281 RMS Pedestals all",nbins,0,5);
   TH1D *h1pedrel = new TH1D("h1pedrel","R3281 RMS Pedestal-Window Ave",nbins,0,5);

   // Double_t Ped0=102;  // for R3620
   Double_t Ped0=98;  // for R3281
   Double_t delta=10; 

  for (int j=0; j<maxids; j++) {

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

       // define histogram to keep results

      sprintf (string,"hped_h%03d\n",j);
      // printf("string=%s",string);
      hped[j] = new TH1D(string,string,nbins,Ped0-delta,Ped0+delta);
      sprintf (string,"hpedrel_h%03d\n",j);
      // printf("string=%s",string);
      hpedrel[j] = new TH1D(string,string,nbins,Ped0-delta,Ped0+delta);
  }

  Int_t evnt_ndx=0;
  Double_t Window_sample[max];

  for (Int_t ievent=0; ievent<maxevents; ievent++){
    tr->GetEntry(ievent);
    if (ievent%1000==0)  cout << " event counter=" << ievent << endl;
    

    Int_t crateid=-1;
    if (rocid <39)     // use for south
      crateid = (rocid-37)*12 + slot -3;
    else
      crateid = (rocid-38)*12 + slot -3;
      if (slot > 10) crateid = crateid -2;

    /*if (rocid <33)       // use for north
      crateid = (rocid-31)*12 + slot -3;
    else
      crateid = (rocid-32)*12 + slot -3;
      if (slot > 10) crateid = crateid -2;*/

    crateid = crateid*16 + channel;

    // cout << " ievent=" << ievent << " rocid=" << rocid << " slot=" << slot << " channel=" << channel << " crateid=" << crateid << endl; 
    outlist  << " ievent=" << ievent << " rocid=" << rocid << " slot=" << slot << " channel=" << channel << " crateid=" << crateid << endl;

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

    // first compute mean for window

    Double_t pedmean = 0;

    for (Int_t jj=0; jj<_samp; jj++){
      
      // Get adc value
      int adc_value = (*waveform)[jj];
      pedmean += adc_value;
    }

    pedmean /= _samp;

    for (Int_t jj=0; jj<_samp; jj++){
      
      // Get adc value
      int adc_value = (*waveform)[jj];
      hped[crateid]->Fill(adc_value);
      hpedrel[crateid]->Fill(adc_value-pedmean + Ped0);
    }
   
    // outlist << "ievent=" << ievent << " crateid=" << crateid << " ievent/maxids=" << ievent/maxids << endl;
  } // end over event loop


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

   c0->SetGridx();
   c0->SetGridy();
   c0->Divide(4,4);

   for (Int_t j=0; j<16; j++) {
      c0->cd(j+1);
      Int_t cr=0;
      Int_t ndx = cr*16 + j;
      hped[ndx]->Draw();
   }

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

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

   for (Int_t j=0; j<16; j++) {
      c1->cd(j+1);
      Int_t cr=0;
      Int_t ndx = cr*16 + j;
      hpedrel[ndx]->Draw();
   }

   // find max and min values of pedestal

   for (Int_t j=0; j<maxids; j++) {
     Int_t cr=j/16;
     Int_t sl=j%16;
     Double_t rmsped=hped[j]->GetRMS();
     Double_t rmspedrel=hpedrel[j]->GetRMS();

     outlist << " j=" << j << " cr=" << cr << " sl=" << sl << " rmsped=" << rmsped << " rmspedrel=" << rmspedrel << endl;
     h2ped->Fill(cr,sl,rmsped);
     h2pedrel->Fill(cr,sl,rmspedrel);
     if (rmsped > 0) h1ped->Fill(rmsped);    // eliminate channels that have pedestal offsets
     h1pedrel->Fill(rmspedrel);
   }

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

   c2->cd(1);
   h2ped->Draw("colz");
   h2ped->GetZaxis()->SetRangeUser(0,delta/2);
   c2->cd(2);
   h2pedrel->Draw("colz");
   h2pedrel->GetZaxis()->SetRangeUser(0,delta/2);

   c2->cd(3);
   h1ped->Draw();
   c2->cd(4);
   h1pedrel->Draw();

   // 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_sigma5.pdf(");
   c0->SaveAs(string);
   sprintf (string,"pedestal_studiesA_sigma5.pdf");
   c1->SaveAs(string);
   sprintf (string,"pedestal_studiesA_sigma5.pdf)");
   c2->SaveAs(string);

   outlist.close();




}