#define Tracer_cxx
#include "Tracer.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>

TH1D* Tracer::Loop(int channel, int pulse, int fadc)
{
//   In a ROOT session, you can do:
//      Root > .L Tracer.C
//      Root > Tracer t
//      Root > t.GetEntry(12); // Fill t data members with entry number 12
//      Root > t.Show();       // Show values of entry 12
//      Root > t.Show(16);     // Read and show values of entry 16
//      Root > t.Loop();       // Loop on all entries
//

//     This is the loop skeleton where:
//    jentry is the global entry number in the chain
//    ientry is the entry number in the current Tree
//  Note that the argument to GetEntry must be:
//    jentry for TChain::GetEntry
//    ientry for TTree::GetEntry and TBranch::GetEntry
//
//       To read only selected branches, Insert statements like:
// METHOD1:
//    fChain->SetBranchStatus("*",0);  // disable all branches
//    fChain->SetBranchStatus("branchname",1);  // activate branchname
// METHOD2: replace line
//    fChain->GetEntry(jentry);       //read all branches
//by  b_branchname->GetEntry(ientry); //read only this branch
   if (fChain == 0) return 0;

   char title[80];
   sprintf(title,"pulse integral spectrum, channel %d", channel);
   TH1D *h1d = new TH1D("spec",title,500,-2.,498.);

   Long64_t nentries = fChain->GetEntriesFast();

   Long64_t nbytes = 0, nb = 0;
   for (Long64_t jentry=2; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);
      nbytes += nb;
      if (jentry/100000*100000 == jentry) {
         std::cout << jentry << std::endl;
      }

      int ped_window[] = {1, 12};
      double mean[2] = {0,0};
      double rms[2] = {0,0};
      int count = 0;
      for (int i=ped_window[0]; i < ped_window[1]; ++i) {
         double h0 = fadc_raw[fadc][channel][i-1];
         double h1 = fadc_raw[fadc][channel][fWindowWidth-i];
         mean[0] += h0;
         mean[1] += h1;
         rms[0] += h0*h0;
         rms[1] += h1*h1;
         ++count;
      }
      mean[0] /= count;
      mean[1] /= count;
      rms[0] /= count;
      rms[1] /= count;
      rms[0] = sqrt(rms[0] - mean[0]*mean[0] + 1e-10);
      rms[1] = sqrt(rms[1] - mean[1]*mean[1] + 1e-10);
      if (rms[0] > 0.03 * mean[0] || rms[1] > 0.03 * mean[1]) {
         continue;
      }
      double pedestal = mean[0];
      int pulse_window_width = 20;

      double pulse_max = 0;
      for (int i0=0; i0 < fWindowWidth; ++i0) {
         double pulse_sum = 0;
         for (int i=0; i < pulse_window_width; ++i) {
            if (i+i0 < fWindowWidth) {
               pulse_sum += fadc_raw[fadc][channel][i+i0];
            }
         }
         pulse_max = (pulse_sum > pulse_max)? pulse_sum : pulse_max;
      }
      pulse_max -= pedestal * pulse_window_width;
      h1d->Fill(pulse_max / pulse_window_width);
   }
   h1d->Draw();
   return h1d;
}

TH1D* Tracer::Spec(int channel, int pulse, int fadc)
{
   if (fChain == 0) return 0;

   char title[80];
   sprintf(title,"triggered pulse spectrum, channel %d", channel);
   if (channel < 6) {
      TH1D *h1d = new TH1D("spec",title,500,-2.,498.);
   }
   else {
      TH1D *h1d = new TH1D("spec",title,500,-2.,498.);
   }

   Long64_t nentries = fChain->GetEntriesFast();

   Long64_t nbytes = 0, nb = 0;
   for (Long64_t jentry=2; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);
      nbytes += nb;
      if (jentry/100000*100000 == jentry) {
         std::cout << jentry << std::endl;
      }

      int ped_window[] = {1, 12};
      double mean[2] = {0,0};
      double rms[2] = {0,0};
      int count = 0;
      for (int i=ped_window[0]; i < ped_window[1]; ++i) {
         double h0 = fadc_raw[fadc][channel][i-1];
         double h1 = fadc_raw[fadc][channel][fWindowWidth-i];
         mean[0] += h0;
         mean[1] += h1;
         rms[0] += h0*h0;
         rms[1] += h1*h1;
         ++count;
      }
      mean[0] /= count;
      mean[1] /= count;
      rms[0] /= count;
      rms[1] /= count;
      rms[0] = sqrt(rms[0] - mean[0]*mean[0] + 1e-10);
      rms[1] = sqrt(rms[1] - mean[1]*mean[1] + 1e-10);
      if (rms[0] > 0.03 * mean[0] || rms[1] > 0.03 * mean[1]) {
         continue;
      }
      double pedestal = mean[0];

      int pulse_window[2];
      if (channel > 5) {
         pulse_window[0] = 107;
         pulse_window[1] = 118;
      }
      else {
         pulse_window[0] = 89;
         pulse_window[1] = 100;
      }
#if TRIGGERED_BY_PULSER
      pulse_window[0] += 28;
      pulse_window[1] += 28;
#endif

      double pulse_sum = 0;
      for (int i=pulse_window[0]; i < pulse_window[1]; ++i) {
         pulse_sum += fadc_raw[fadc][channel][i];
      }
      h1d->Fill(pulse_sum/(pulse_window[1] - pulse_window[0]) - pedestal);
   }
   h1d->Draw();
   return h1d;
}

Double_t Tracer::PulseInt(TH1D* hist)
{
   int ped_window[] = {1, 12};
   int pulse_window_width = 20;
   double mean[2] = {0,0};
   double rms[2] = {0,0};
   int count = 0;
   for (int i=ped_window[0]; i < ped_window[1]; ++i) {
      double h0 = hist->GetBinContent(i);
      double h1 = hist->GetBinContent(hist->GetNbinsX()+1-i);
      mean[0] += h0;
      mean[1] += h1;
      rms[0] += h0*h0;
      rms[1] += h1*h1;
      ++count;
   }
   mean[0] /= count;
   mean[1] /= count;
   rms[0] /= count;
   rms[1] /= count;
   rms[0] = sqrt(rms[0] - mean[0]*mean[0] + 1e-10);
   rms[1] = sqrt(rms[1] - mean[1]*mean[1] + 1e-10);
   if (rms[0] > 0.03 * mean[0] || rms[1] > 0.03 * mean[1]) {
      std::cout << mean[0] << ", " << mean[1] << ", "
                << rms[0] << ", " << rms[1] << std::endl;
      return 0;
   }
   double pedestal = mean[0];

   double pulse_max = 0;
   for (int i0=1; i0 <= hist->GetNbinsX(); ++i0) {
      double pulse_sum = 0;
      for (int i=0; i < pulse_window_width; ++i) {
         if (i+i0 <= hist->GetNbinsX()) {
            pulse_sum += hist->GetBinContent(i+i0);
         }
      }
      pulse_max = (pulse_sum > pulse_max)? pulse_sum : pulse_max;
   }
   pulse_max -= pedestal * pulse_window_width;
   return pulse_max / pulse_window_width;
}