#include <utils/yodaUtility.h>
#include "TrkFunctions.cpp"
#include <iostream>
#include <fstream>
#include <TPaveText.h>
#include <TLatex.h>
#include <TCanvas.h>
#include <TGraph.h>
#include <TStyle.h>
#include <TObjString.h>
#include <PscuHeader.h>
#include <EventHeader.h>
#include <CalibTrk1Event.h>
#include <CalibTrk2Event.h>

/**
 * fillpedsig  --> is a function to calculate the mean value for
 *                 sigma, variance of sigma, pedestal and variance of pedestal
 *                 and fill a root.file whit these values
 *
 * autor: D.Fedele
 * version 1.0
 * Parameters:
 * 	file - the path to the root file to analyse
 *      outdir - total path of output file 
 * 
 */

void fillpedsig(TString file,TString outdir)
{
  // 
  //   obtain information about the data file and select the output file
  const string filepath=file.Data(); 
  Int_t dwpos = filepath.find("DW_");
  Int_t dwpos1 = filepath.find(".root");
  TString fpath=(filepath.c_str());
  TString base,ffile,date1,date2,date;
  stringcopy(base,fpath,0,dwpos);
  stringcopy(ffile,fpath,dwpos,dwpos1);

  stringcopy(date1,fpath,dwpos+3,dwpos+9);
  stringcopy(date2,fpath,dwpos+10,dwpos1);

  TString out;
  stringstream outfile;
  if(outdir.Length()==0){
    out = base;
  }else{
    out = outdir;
  };
  outfile<<out.Data()<<"pedsig.root";
  
  //
  //  initialise structure where data will store 
  struct caltrk_def ctrk;
  Int_t nevents=0;
  Float_t OBT[2];
  
  OBT[0]=0;
  OBT[1]=0;
  ctrk.good0[0]=0;
  ctrk.good0[1]=0;
  for(Int_t i=0;i<12;i++){ 
    ctrk.daqmode[i]=0;
    ctrk.dspnum[i]=0;
    ctrk.calibnum[i]=0;
    ctrk.ncalev[i]=0;
    ctrk.calfl[i]=0;
    ctrk.ped1[i]=0;
    ctrk.ped2[i]=0;
    ctrk.ped3[i]=0;
    ctrk.sig1[i]=0;
    ctrk.sig2[i]=0;
    ctrk.sig3[i]=0;
    ctrk.nbad1[i]=0;
    ctrk.nbad2[i]=0;
    ctrk.nbad3[i]=0;
    ctrk.crc_hc[i]=0;
    ctrk.crc_c[i][0]=0;
    ctrk.crc_c[i][1]=0;
    ctrk.crc_c[i][2]=0;
    for(Int_t iii=0;iii<3072;iii++){
      ctrk.dspped[i][iii]=0;
      ctrk.dspsig[i][iii]=0;
      ctrk.dspbad[i][iii]=0;
    }
  }

  //
  // open data files
  pamela::EventHeader *eh1=0,*eh2=0;
  pamela::PscuHeader *ph1=0,*ph2=0;
  pamela::CalibTrk1Event *trk1 = 0;
  pamela::CalibTrk2Event *trk2 = 0;

  TFile *datafile = new TFile(file);
  if ( !datafile ){
    //  printf("No data file, exiting...\n");
    return;
  };  
  TTree *otr1,*otr2;
  
  otr1 = (TTree*)datafile->Get("CalibTrk1");
  otr1->SetBranchAddress("CalibTrk1", &trk1);
  otr1->SetBranchAddress("Header",&eh1);
  otr2 = (TTree*)datafile->Get("CalibTrk2");
  otr2->SetBranchAddress("CalibTrk2", &trk2);
  otr2->SetBranchAddress("Header",&eh2);
  

  if(otr1->GetEntries()==otr2->GetEntries())
    nevents = otr1->GetEntries();
  else{
    //    printf("WARNING: CalibTrk1 entries is different from CalibTrk2 entries");
    return;}

  if (nevents<=0) {
    datafile->Close();
    //    printf("No calibration packets found, exiting...\n");
    return;
  };

  //
  // Define variables
  Int_t ndsp =0,numbad[12][12];
  Float_t pedav[12][12],pedavtemp[12][12],sigav[12][12],sigavtemp[12][12];
  Float_t pedvar[12][12],pedvartemp[12][12],sigvar[12][12],sigvartemp[12][12];
  Long64_t datefile=0;


  //
  // open output files
  TFile *pedsig;
  TTree *pstree;
  ifstream pedsigtmp(outfile.str().c_str());
  if(pedsigtmp.fail()){
    pedsigtmp.close();
    pedsig = new TFile(outfile.str().c_str(),"new");
    pstree=new TTree("pstree","pedsig");
    pstree->Branch("pedav",pedav,"pedav[12][12]/F");
    pstree->Branch("pedvar",pedvar,"pedvar[12][12]/F");
    pstree->Branch("sigav",sigav,"sigav[12][12]/F");
    pstree->Branch("sigvar",sigvar,"sigvar[12][12]/F");
    pstree->Branch("datefile",&datefile,"datefile/L");
    pstree->Branch("numbad",numbad,"numbad[12][12]/I");
  }
  if(!pedsigtmp.fail()){
    pedsigtmp.close();
    pedsig = new TFile(outfile.str().c_str(),"update");
    pstree = (TTree*)pedsig->Get("pstree");
    pstree->SetBranchAddress("pedav",pedav);
    pstree->SetBranchAddress("pedvar",pedvar);
    pstree->SetBranchAddress("sigav",sigav);
    pstree->SetBranchAddress("sigvar",sigvar);
    pstree->SetBranchAddress("datefile",&datefile);
    pstree->SetBranchAddress("numbad",numbad);
  }

  datefile=atoll(date1+date2);
  //**********************************************************************
  //
  //   LOOP OVER EVENTS
  //
  //**********************************************************************

  for (Int_t i = 0; i < nevents; i++){     
      
    otr1->GetEntry(i); 
    otr2->GetEntry(i); 
      
    //
    // fill data storage structure
    ctrk.good0[0]=trk1->good0;
    ctrk.good0[1]=trk2->good0;
    for (Int_t m = 0; m < 6; m++){
      Int_t plane=m;
      ph1 = eh1->GetPscuHeader();
      OBT[0]= ph1->GetOrbitalTime();
      ctrk.daqmode[trk1->DSPnumber[m]-1]=trk1->DAQmode[m];
      ctrk.dspnum[trk1->DSPnumber[m]-1]=trk1->DSPnumber[m];
      ctrk.calibnum[trk1->DSPnumber[m]-1]=trk1->calibnumber[m];
      ctrk.ncalev[trk1->DSPnumber[m]-1]=trk1->ncalib_event[m];
      ctrk.ped1[trk1->DSPnumber[m]-1]=trk1->ped_l1[m];
      ctrk.ped2[trk1->DSPnumber[m]-1]=trk1->ped_l2[m];
      ctrk.ped3[trk1->DSPnumber[m]-1]=trk1->ped_l3[m];
      ctrk.sig1[trk1->DSPnumber[m]-1]=trk1->sig_l1[m];
      ctrk.sig2[trk1->DSPnumber[m]-1]=trk1->sig_l2[m];
      ctrk.sig3[trk1->DSPnumber[m]-1]=trk1->sig_l3[m];
      ctrk.nbad1[trk1->DSPnumber[m]-1]=trk1->nbad_l1[m];
      ctrk.nbad2[trk1->DSPnumber[m]-1]=trk1->nbad_l2[m];
      ctrk.nbad3[trk1->DSPnumber[m]-1]=trk1->nbad_l3[m];
      ctrk.calfl[trk1->DSPnumber[m]-1]=trk1->cal_flag[m];
      ctrk.crc_c[trk1->DSPnumber[m]-1][0]=trk1->crc_cal[m][0];
      ctrk.crc_c[trk1->DSPnumber[m]-1][1]=trk1->crc_cal[m][1];
      ctrk.crc_c[trk1->DSPnumber[m]-1][2]=trk1->crc_cal[m][2];
      ctrk.crc_hc[trk1->DSPnumber[m]-1]=trk1->crc_hcal[m];
      for (Int_t j = 0; j < 3072; j++){
	ctrk.dspped[trk1->DSPnumber[m]-1][j]=trk1->DSPped_par[m][j];
	ctrk.dspsig[trk1->DSPnumber[m]-1][j]=trk1->DSPsig_par[m][j];
	ctrk.dspbad[trk1->DSPnumber[m]-1][j]=trk1->DSPbad_par[m][j];
      }
      ph2 = eh2->GetPscuHeader();
      OBT[1]= ph2->GetOrbitalTime();
      ctrk.daqmode[trk2->DSPnumber[m]-1]=trk2->DAQmode[m];
      ctrk.dspnum[trk2->DSPnumber[m]-1]=trk2->DSPnumber[m];
      ctrk.calibnum[trk2->DSPnumber[m]-1]=trk2->calibnumber[m];
      ctrk.ncalev[trk2->DSPnumber[m]-1]=trk2->ncalib_event[m];
      ctrk.ped1[trk2->DSPnumber[m]-1]=trk2->ped_l1[m];
      ctrk.ped2[trk2->DSPnumber[m]-1]=trk2->ped_l2[m];
      ctrk.ped3[trk2->DSPnumber[m]-1]=trk2->ped_l3[m];
      ctrk.sig1[trk2->DSPnumber[m]-1]=trk2->sig_l1[m];
      ctrk.sig2[trk2->DSPnumber[m]-1]=trk2->sig_l2[m];
      ctrk.sig3[trk2->DSPnumber[m]-1]=trk2->sig_l3[m];
      ctrk.nbad1[trk2->DSPnumber[m]-1]=trk2->nbad_l1[m];
      ctrk.nbad2[trk2->DSPnumber[m]-1]=trk2->nbad_l2[m];
      ctrk.nbad3[trk2->DSPnumber[m]-1]=trk2->nbad_l3[m];
      ctrk.calfl[trk2->DSPnumber[m]-1]=trk2->cal_flag[m];
      ctrk.crc_c[trk1->DSPnumber[m]-1][0]=trk2->crc_cal[m][0];
      ctrk.crc_c[trk1->DSPnumber[m]-1][1]=trk2->crc_cal[m][1];
      ctrk.crc_c[trk1->DSPnumber[m]-1][2]=trk2->crc_cal[m][2];
      ctrk.crc_hc[trk1->DSPnumber[m]-1]=trk2->crc_hcal[m];
      for (Int_t j = 0; j < 3072; j++){
	ctrk.dspped[trk2->DSPnumber[m]-1][j]=trk2->DSPped_par[m][j];
	ctrk.dspsig[trk2->DSPnumber[m]-1][j]=trk2->DSPsig_par[m][j];
	ctrk.dspbad[trk2->DSPnumber[m]-1][j]=trk2->DSPbad_par[m][j];
      }
    }      
  
    //
    // inizialise variables
    for(Int_t n = 0; n<12; n++){
      for(Int_t nm = 0; nm<12; nm++){
	pedav[n][nm]=0;
	pedavtemp[n][nm]=0;
	sigav[n][nm]=0;
	sigavtemp[n][nm]=0;
	pedvar[n][nm]=0;
	pedvartemp[n][nm]=0;
	sigvar[n][nm]=0;
	sigvartemp[n][nm]=0;
	numbad[n][nm]=0;
      }
    }

    //
    //  calculate the mean value for sigma, variance of sigma, pedestal, variance of pedestal
    Int_t nn=0;
    for(Int_t n = 0; n<12; n++){
      
      ndsp = ctrk.dspnum[n];
      nn = ndsp-1;

      for(Int_t j = 0; j < 3072; j++){
	if(ctrk.dspbad[nn][j]==0){
	  sigavtemp[nn][j/256]+=ctrk.dspsig[nn][j];
	  pedavtemp[nn][j/256]+=ctrk.dspped[nn][j];
	}
	else 
	  numbad[nn][j/256]+=1;
      }
      for(Int_t ii=0;ii<12;ii++){
	pedav[nn][ii]=pedavtemp[nn][ii]/256;  
	sigav[nn][ii]=sigavtemp[nn][ii]/256;  
      }
      for(Int_t j = 0; j < 3072; j++){
	if(ctrk.dspbad[nn][j]==0){
	  sigvartemp[nn][j/256]+=pow((ctrk.dspsig[nn][j]-sigav[nn][j/256]),2);
	  pedvartemp[nn][j/256]+=pow((ctrk.dspped[nn][j]-pedav[nn][j/256]),2);
	};
      }
      for(Int_t ii=0;ii<12;ii++){
	pedvar[nn][ii]=sqrt(pedvartemp[nn][ii]/255); 
	sigvar[nn][ii]=sqrt(sigvartemp[nn][ii]/255); 
      }
    }

    //
    // Fill the output Tree
    pstree->Fill();    

  };//end loop on events

  //
  //close all files
  datafile->Close();  
  pedsig->Write();
  pedsig->Close();  
  
  return;
}