// compile with:
// g++ -o plotchipzsuppeaks.exe plotchipzsuppeaks.cpp `root-config --cflags --glibs`
//
#include "TH1F.h"
#include "TProfile.h"
#include "TMath.h"
#include "TF1.h"
#include "TLegend.h"
#include "TCanvas.h"
 //
#include "TROOT.h"
//#include "TStyle.h"
#include "TTree.h"
#include "TFile.h"
#include "TGraphErrors.h"
#include "TStyle.h"
#include "TLine.h"
#include <string>
#include <cmath>
#include <iostream>
#include <fstream>
#include <vector>
#include <algorithm>
//
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>

using namespace std;

int plotchipzsuppeaks(int mcm = 0,const int iEvent = 1,int chip = 0, int saltronb = 0, const char *dir = "",const char *filename = "", const char *title = "test")
{
  char name[100];

  int peakthr = 300;
  int peaktype = 0;
  int peaktime;
  unsigned int ErrorCode[16] = {0};
  unsigned int WarningCode[16] = {0};
  float ratioLow;
  float ratioHigh;
  int nbpeakstype[3];
  TH1F* hPeaktype0[16];
  TH1F* hPeaktype1[16];
  TH1F* hPeaktype2[16];
  

  Char_t iChan; // Channel number
  Char_t iChip; // Channel number
  UShort_t iClLength; //number of samples
  UShort_t iClTime; // Time of sample
  UShort_t iMcm; // mcm dcl address
  UInt_t iEventnb;
  Int_t iEntry;
  Int_t timeSamples[1024];
  Int_t iSamples[1024]; // Samples
  TH1F* hSignal[16];
  TH1F* hClLength[16];
  TH1F* hClusters;
  TH1F* hRatioLow;
  TH1F* hRatioHigh;

  
  FILE *limitfile;
  sprintf(name, "%s.txt","plotchipzsuppeaks");
  limitfile = fopen(name,"r");
  if (limitfile == NULL) {
    return 4;
  }
  
  int item = 0;
  float peak_min_mean[3];
  float peak_max_mean[3];
  float peak_max_rms[3];
  float cluster_min_mean;
  float cluster_max_mean;
  float cluster_max_rms;
  int clusters_min;
  int clusters_max;
  
    
  while (fgets(name,sizeof(name),limitfile) != NULL) {
    if (name[0] != '#') {
      switch(item) {
      case 0: sscanf(name,"%f %f %f",&peak_min_mean[0],&peak_max_mean[0],&peak_max_rms[0]); break;
      case 1: sscanf(name,"%f %f %f",&peak_min_mean[1],&peak_max_mean[1],&peak_max_rms[1]); break;
      case 2: sscanf(name,"%f %f %f",&peak_min_mean[2],&peak_max_mean[2],&peak_max_rms[2]); break;
      case 3: sscanf(name,"%f %f %f %d %d",&cluster_min_mean,&cluster_max_mean,&cluster_max_rms,&clusters_min,&clusters_max); break;
      }
      item++;
    }
  }
  fclose(limitfile);
  if (item == 4) {
    printf("%f %f %f\n",peak_min_mean[0],peak_max_mean[0],peak_max_rms[0]);
    printf("%f %f %f\n",peak_min_mean[1],peak_max_mean[1],peak_max_rms[1]);
    printf("%f %f %f\n",peak_min_mean[2],peak_max_mean[2],peak_max_rms[2]);
    printf("%f %f %f %d %d\n",cluster_min_mean,cluster_max_mean,cluster_max_rms,clusters_min,clusters_max);
  }
  else return 5;
  
  sprintf(name, "%s/%s.root",dir,filename);
  TFile *inFile = new TFile(name);
  TTree *tree = (TTree*) inFile->Get("tree");
  if(!tree) {
    cerr << "Tree SALTRO events was not found: " << name << "\n";
    return 2;
  }
  sprintf(name, "Event_%d", iEvent);
  TCanvas *c1 = new TCanvas("c1",name,800,800);

  tree->SetBranchAddress("iEventnb",&iEventnb);
  tree->SetBranchAddress("iSamples",&iSamples);
  tree->SetBranchAddress("iChan",&iChan);
  tree->SetBranchAddress("iChip",&iChip);
  tree->SetBranchAddress("iClLength",&iClLength);
  tree->SetBranchAddress("iClTime",&iClTime);

  iEventnb = 0;
  iEntry = 0;
  for (Int_t j = 0; j < 16; j++) {
    sprintf(name, "hSignal_%d", j);
    hSignal[j] = ((TH1F *)(gROOT->FindObject(name)));
    if (hSignal[j]) delete hSignal[j];
    hSignal[j] = new TH1F(name, "",1024,-0.5,1023.5);
    hSignal[j]->SetMinimum(0);
    hSignal[j]->SetStats(0);
  }
  
  for (Int_t j = 0; j < 16; j++) {
    sprintf(name, "hClLegth_%d", j);
    hClLength[j] = ((TH1F *)(gROOT->FindObject(name)));
    if (hClLength[j]) delete hClLength[j];
    hClLength[j] = new TH1F(name, "",1024,-0.5,1023.5);
    hClLength[j]->SetMinimum(0);
  }

  
  sprintf(name, "hClusters");
  hClusters = ((TH1F *)(gROOT->FindObject(name)));
  if (hClusters) delete hClusters;
  hClusters = new TH1F(name, "",16,-0.5,15.5);
  hClusters->GetXaxis()->SetTitle("Clusters");
  hClusters->GetYaxis()->SetTitle("Channel");
  hClusters->SetMinimum(0);
  hClusters->SetStats(0);

  sprintf(name, "hRatioLow");
  hRatioLow = ((TH1F *)(gROOT->FindObject(name)));
  if (hRatioLow) delete hRatioLow;
  hRatioLow = new TH1F(name, "",120,0.4,1.6);
  hRatioLow->SetTitle("RatioLow");
  hRatioLow->GetXaxis()->SetTitle("Entries");
  hRatioLow->GetYaxis()->SetTitle("RatioLow");
  hRatioLow->SetMinimum(0);
  hRatioLow->SetStats(11110);
  sprintf(name, "hRatioHigh");
  hRatioHigh = ((TH1F *)(gROOT->FindObject(name)));
  if (hRatioHigh) delete hRatioHigh;
  hRatioHigh = new TH1F(name, "",120,0.4,1.6);
  hRatioHigh->SetTitle("RatioHigh");
  hRatioHigh->GetXaxis()->SetTitle("Entries");
  hRatioHigh->GetYaxis()->SetTitle("RatioHigh");
  hRatioHigh->SetMinimum(0);
  hRatioHigh->SetStats(11110);
  

  for (Int_t j = 0; j < 16; j++) {
    sprintf(name, "hPeaktype0_%d", j);
    hPeaktype0[j] = ((TH1F *)(gROOT->FindObject(name)));
    if (hPeaktype0[j]) delete hPeaktype0[j];
    hPeaktype0[j] = new TH1F(name, "",1024,-0.5,1023.5);
    hPeaktype0[j]->SetMinimum(0);
  }
  for (Int_t j = 0; j < 16; j++) {
    sprintf(name, "hPeaktype1_%d", j);
    hPeaktype1[j] = ((TH1F *)(gROOT->FindObject(name)));
    if (hPeaktype1[j]) delete hPeaktype1[j];
    hPeaktype1[j] = new TH1F(name, "",1024,-0.5,1023.5);
    hPeaktype1[j]->SetMinimum(0);
  }
  for (Int_t j = 0; j < 16; j++) {
    sprintf(name, "hPeaktype2_%d", j);
    hPeaktype2[j] = ((TH1F *)(gROOT->FindObject(name)));
    if (hPeaktype2[j]) delete hPeaktype2[j];
    hPeaktype2[j] = new TH1F(name, "",1024,-0.5,1023.5);
    hPeaktype2[j]->SetMinimum(0);
  }

  
  Int_t nentries = tree->GetEntries(); 
  printf("ENTRIES %d\n",nentries);

  
  
  gStyle->SetOptStat(1110);
  int maxClLength = 0;
  int oldChan = -1;
  int nbClusters = 0;
  nbpeakstype[0] =  nbpeakstype[1] = nbpeakstype[2] = 0;
  for (Int_t j = 0; j < nentries; j++) {
    tree->GetEntry(j);
    if (iChan != oldChan) {
      if (oldChan != -1) {
	peaktime = 253;
	if (nbClusters == 1) nbClusters = 12;
	if (nbClusters == 12) {
	for(Int_t l = 0; l < nbClusters; l++) {
	  if (l == 0) {
	    peaktype = 0;
	    if (timeSamples[peaktime] > peakthr)  {
	      ratioLow = float(timeSamples[peaktime-1])/float(timeSamples[peaktime]);
	      ratioHigh = float(timeSamples[peaktime+1])/float(timeSamples[peaktime]);
	      hRatioLow->Fill(ratioLow);
	      hRatioHigh->Fill(ratioHigh);
	      if ((ratioHigh < 0.9) && (ratioLow > 0.6)) peaktype = 1;
	      else if ((ratioLow <= 0.6) || (ratioHigh >= 0.9))  peaktype = 2;
	    }
	    nbpeakstype[peaktype]++;
	  }
	  if (peaktype == 0) hPeaktype0[oldChan]->Fill(timeSamples[peaktime]);
	  else if (peaktype == 1) hPeaktype1[oldChan]->Fill(timeSamples[peaktime]);
	  else if (peaktype == 2) hPeaktype2[oldChan]->Fill(timeSamples[peaktime]);
	  peaktime += 60;
	}
	}
      }
      memset(timeSamples,0,sizeof(timeSamples));
      oldChan = iChan;
      nbClusters = 0;
    }
    nbClusters++;
    for(Int_t l = 0; l < iClLength; l++) {
      timeSamples[iClTime-l] =  iSamples[l];
    }
    if (iEventnb == iEvent) {
      for(Int_t l = 0; l < iClLength; l++) {
	hSignal[iChan]->AddBinContent(iClTime-l, iSamples[l]);
      }
    }
    hClLength[iChan]->Fill(iClLength);
    if (iClLength > maxClLength) maxClLength = iClLength;
    hClusters->Fill(iChan);    
  }

  Int_t oldLevel = gErrorIgnoreLevel;
  gErrorIgnoreLevel = kWarning;

  for (Int_t j = 0; j < 16; j++) {
    c1->cd(iChan+1);
    sprintf(name, "Event %d Channel %d",iEvent,j);
    hSignal[j]->SetTitle(name);
    hSignal[j]->GetXaxis()->SetTitle("Sample #");
    hSignal[j]->GetYaxis()->SetTitle("ADC value");
    hSignal[j]->SetAxisRange(0,1024, "Y");
    hSignal[j]->SetAxisRange(650,830, "X");
    hSignal[j]->Draw();
    TLine *line1 = new TLine(252+60*7,0,252+60*7,1024);
    line1->SetLineColor(kBlue);
    line1->Draw();
    TLine *line2 = new TLine(252+60*8,0,252+60*8,1024);
    line2->SetLineColor(kRed);
    line2->Draw();
    TLine *line3 = new TLine(252+60*9,0,252+60*9,1024);
    line3->SetLineColor(kBlue);
    line3->Draw();
    if (j == 0) sprintf(name,"%s/chip%d/saltro%d-mcm%d-c%d-zerosup-%s.pdf(",dir,saltronb,saltronb,mcm,chip,title);
    else  sprintf(name,"%s/chip%d/saltro%d-mcm%d-c%d-zerosup-%s.pdf",dir,saltronb,saltronb,mcm,chip,title);
    c1->SaveAs(name,"pdf");
  }

  sprintf(name,"%s/chip%d/saltro%d-mcm%d-c%d-zerosup-%s.pdf",dir,saltronb,saltronb,mcm,chip,title);  
  for (Int_t j = 0; j < 16; j++) {
    c1->cd(iChan+1);
    sprintf(name, "Channel %d Cluster Length",j);
    hClLength[j]->SetTitle(name);
    hClLength[j]->GetXaxis()->SetTitle("Entries");
    hClLength[j]->GetYaxis()->SetTitle("Length");
    hClLength[j]->SetAxisRange(0, maxClLength+5, "X");
    hClLength[j]->Draw();
    sprintf(name,"%s/chip%d/saltro%d-mcm%d-c%d-zerosup-%s.pdf",dir,saltronb,saltronb,mcm,chip,title);
    c1->SaveAs(name,"pdf");
    if ((hClLength[j]->GetMean() < cluster_min_mean) || (hClLength[j]->GetMean() > cluster_max_mean)) ErrorCode[j] |= 0x0100;
    if (hClLength[j]->GetRMS() > cluster_max_rms) ErrorCode[j] |= 0x0080;
    if ((hClLength[j]->GetEntries() < clusters_min) || (hClLength[j]->GetEntries() > clusters_max)) ErrorCode[j] |= 0x0200;
  }

  hRatioLow->Draw();
  c1->SaveAs(name,"pdf");
  hRatioHigh->Draw();
  c1->SaveAs(name,"pdf");
  hClusters->Draw();
  sprintf(name,"%s/chip%d/saltro%d-mcm%d-c%d-zerosup-%s.pdf)",dir,saltronb,saltronb,mcm,chip,title);  
  c1->SaveAs(name,"pdf");
  
  printf("Peaktype0\n");
  for (Int_t j = 0; j < 16; j++) {
      c1->cd(j+1);
      sprintf(name, "Pulstrain peak0 ch%d %s",j ,title);
      hPeaktype0[j]->SetTitle(name);
      hPeaktype0[j]->GetYaxis()->SetTitle("Entries");
      hPeaktype0[j]->GetXaxis()->SetTitle("ADC value");
      //    hPeaktype1[j]->SetAxisRange(0, 1023, "X");
      //hPeaktype1[j]->SetAxisRange(0, 400, "X");
      hPeaktype0[j]->Draw();
      if (j == 0) sprintf(name,"%s/chip%d/saltro%d-mcm%d-c%d-zerosup-peaktype0-%s.pdf(",dir,saltronb,saltronb,mcm,chip,title);
      else if ( j == 15) sprintf(name,"%s/chip%d/saltro%d-mcm%d-c%d-zerosup-peaktype0-%s.pdf)",dir,saltronb,saltronb,mcm,chip,title);
      else  sprintf(name,"%s/chip%d/saltro%d-mcm%d-c%d-zerosup-peaktype0-%s.pdf",dir,saltronb,saltronb,mcm,chip,title);
      c1->SaveAs(name,"pdf");
      if ((hPeaktype0[j]->GetMean() < peak_min_mean[0]) || (hPeaktype0[j]->GetMean() > peak_max_mean[0])) ErrorCode[j] |= 0x0001;
      if (hPeaktype0[j]->GetRMS() > peak_max_rms[0]) ErrorCode[j] |= 0x0002;
      printf("%2d\t%6.2f\t%5.3f\t%6.1f\n",j+iChip*16,hPeaktype0[j]->GetMean(),hPeaktype0[j]->GetRMS(),hPeaktype0[j]->GetEntries());
  }
  printf("Peaktype1\n");
  for (Int_t j = 0; j < 16; j++) {
      c1->cd(j+1);
      sprintf(name, "Pulstrain peak1 ch%d %s",j ,title);
      hPeaktype1[j]->SetTitle(name);
      hPeaktype1[j]->GetYaxis()->SetTitle("Entries");
      hPeaktype1[j]->GetXaxis()->SetTitle("ADC value");
      //    hPeaktype1[j]->SetAxisRange(0, 1023, "X");
      //hPeaktype1[j]->SetAxisRange(0, 400, "X");
      hPeaktype1[j]->Draw();
      if (j == 0) sprintf(name,"%s/chip%d/saltro%d-mcm%d-c%d-zerosup-peaktype1-%s.pdf(",dir,saltronb,saltronb,mcm,chip,title);
      else if ( j == 15) sprintf(name,"%s/chip%d/saltro%d-mcm%d-c%d-zerosup-peaktype1-%s.pdf)",dir,saltronb,saltronb,mcm,chip,title);
      else  sprintf(name,"%s/chip%d/saltro%d-mcm%d-c%d-zerosup-peaktype1-%s.pdf",dir,saltronb,saltronb,mcm,chip,title);
      c1->SaveAs(name,"pdf");
      if ((hPeaktype1[j]->GetMean() < peak_min_mean[1]) || (hPeaktype1[j]->GetMean() > peak_max_mean[1])) {
	ErrorCode[j] |= 0x0004;
	//	printf("ERR %f %f %f\n",hPeaktype1[j]->GetMean(),peak_min_mean[1],peak_max_mean[1]);
      }
      if (hPeaktype1[j]->GetRMS() > peak_max_rms[1]) ErrorCode[j] |= 0x0008;
      printf("%2d\t%6.2f\t%5.3f\t%6.1f\n",j+iChip*16,hPeaktype1[j]->GetMean(),hPeaktype1[j]->GetRMS(),hPeaktype1[j]->GetEntries());
  }
  printf("Peaktype2\n");
  for (Int_t j = 0; j < 16; j++) {
      c1->cd(j+1);
      sprintf(name, "Pulstrain peak2 ch%d %s",j ,title);
      hPeaktype2[j]->SetTitle(name);
      hPeaktype2[j]->GetYaxis()->SetTitle("Entries");
      hPeaktype2[j]->GetXaxis()->SetTitle("ADC value");
      //    hPeaktype1[j]->SetAxisRange(0, 1023, "X");
      //hPeaktype1[j]->SetAxisRange(0, 400, "X");
      hPeaktype2[j]->Draw();
      if (j == 0) sprintf(name,"%s/chip%d/saltro%d-mcm%d-c%d-zerosup-peaktype2-%s.pdf(",dir,saltronb,saltronb,mcm,chip,title);
      else if ( j == 15) sprintf(name,"%s/chip%d/saltro%d-mcm%d-c%d-zerosup-peaktype2-%s.pdf)",dir,saltronb,saltronb,mcm,chip,title);
      else  sprintf(name,"%s/chip%d/saltro%d-mcm%d-c%d-zerosup-peaktype2-%s.pdf",dir,saltronb,saltronb,mcm,chip,title);
      c1->SaveAs(name,"pdf");
      if ((hPeaktype2[j]->GetMean() < peak_min_mean[2]) || (hPeaktype2[j]->GetMean() > peak_max_mean[2])) ErrorCode[j] |= 0x0010;
      if (hPeaktype2[j]->GetRMS() > peak_max_rms[2]) ErrorCode[j] |= 0x0020;
      printf("%2d\t%6.2f\t%5.3f\t%6.1f\n",j+iChip*16,hPeaktype2[j]->GetMean(),hPeaktype2[j]->GetRMS(),hPeaktype2[j]->GetEntries());
  }

  gErrorIgnoreLevel = oldLevel;   
  printf("PEAKTYPE %d %d %d\n",nbpeakstype[0],nbpeakstype[1],nbpeakstype[2]);

  // Store result in csvfile
  struct stat   fileBuffer;
  int csvexist = 0;
  if (stat (filename, &fileBuffer) == 0) csvexist = 1;

  FILE *csvfile;
  sprintf(name,"%s/chip%d/saltro%d-mcm%d-c%d-zerosup-%s.csv",dir,saltronb,saltronb,mcm,chip,title);
  //  sprintf(name, "%s/%s.csv",dir,filename);
  csvfile = fopen(name,"a");
  if (csvfile == NULL) {
    return 3;
  }
  if (csvexist == 0) {
    fprintf(csvfile,"Channel,Peak0Mean,Peak0RMS,Peak0Entries,");
    fprintf(csvfile,"Peak1Mean,Peak1RMS,Peak1Entries,");
    fprintf(csvfile,"Peak2Mean,Peak2RMS,Peak2Entries,");
    fprintf(csvfile,"ClLengthMean,CllengthRMS,ClLengthEntries,");
    fprintf(csvfile,"Error,Warning\n");   
  }
 
  for (int j = 0; j < 16; j++) {
    fprintf(csvfile,"%2d,%6.2f,%6.3f,%6.1f,",j+iChip*16,hPeaktype0[j]->GetMean(),hPeaktype0[j]->GetRMS(),hPeaktype0[j]->GetEntries());
    fprintf(csvfile,"%6.2f,%6.3f,%6.1f,",hPeaktype1[j]->GetMean(),hPeaktype1[j]->GetRMS(),hPeaktype1[j]->GetEntries());
    fprintf(csvfile,"%6.2f,%6.3f,%6.1f,",hPeaktype2[j]->GetMean(),hPeaktype2[j]->GetRMS(),hPeaktype2[j]->GetEntries());
    fprintf(csvfile,"%6.2f,%6.3f,%6.1f,",hClLength[j]->GetMean(),hClLength[j]->GetRMS(),hClLength[j]->GetEntries());
    fprintf(csvfile,"%08X,%08X\n",ErrorCode[j],WarningCode[j]);
  }
  fclose(csvfile);
  return 0;
}


int main (int argc, char *argv[])
{
  int iret = 0;
  if (argc != 8) {
    printf("Wrong number of arguments: argc = %d\nSyntax should be:\n./plotchipzsuppeaks.exe eventnb saltro mcm chip dir filename title\n",argc);
    return 1;
  }
  int eventnb = atoi(argv[1]);
  int saltronb = atoi(argv[2]);
  int mcm = atoi(argv[3]);
  int chip = atoi(argv[4]);
  iret = plotchipzsuppeaks(mcm,eventnb,chip,saltronb,argv[5],argv[6],argv[7]);
  return iret;
}