#define analyse2_cxx
#include "analyse2.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <string.h>
#include <stdio.h>
#include <iostream>
#include <vector>
#include <utility>
#include <TFile.h>
#include <TTree.h>
#include <TSystemDirectory.h>
#include <gsl/gsl_statistics.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_complex.h>

#include <TF1.h>
#include <TGraphErrors.h>
using namespace std;



void analyse2::Loop()
{
//   In a ROOT session, you can do:
//      Root > .L analyse2.C
//      Root > analyse2 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->ve types of modulesSetBranchStatus("*",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;
   
   Double_t totalcurrent = 0.;
   Double_t totaltime = 0.;
   const int length = 100;
   double array[length] = {0.};
   double arrayavg = 0.;
    TF1 * bkgfunc = new TF1("bkgfunc","[0]+[1]*TMath::Cos(x*[2]*(2*3.14159)+[3])");
    bkgfunc->SetParameters( bkgfitPar0[64], bkgfitPar1[64], bkgfitPar2[64], bkgfitPar3[64]);
    TF1 * gausfunc = new TF1("gausfunc","gaus(0)+[3]");
    TF1 * gausfunc2 = new TF1("gausfunc2","gaus(0)+[3]");

   // TF1 * gausfunc2 = new TF1("gausfunc2","[0]*exp(-((x-[1])/[2])^2) +[3]*exp(-((x-[1])/[4])^2) + [5]");
    TGraphErrors * gausgraph;
  
    int numoverthresh=0;
    double maxchannelamp = 0.;
    double threshold = 10.;
    int maxchannel = 0;
    int numtocalc = 0;
    int  sidenumtocalc= 0;
    Long64_t nentries = fChain->GetEntries();
    Long64_t ientry = LoadTree(1);
    fChain->GetEntry(1);  
    Double_t basetimeoffset = time;
    Double_t channelamp[128] = {0.};
    Double_t abschannelamp[128] = {0.};
    Double_t channelamp2[128] = {0.};
    Double_t abschannelamp2[128] = {0.};

    Bool_t beamonflag = false;
    Int_t waitcounter = 0;
    Int_t lastfit = 0;
    Int_t lastfit2 = 0;
    Double_t channellist_gsl[128];
    Double_t channellist_gsl2[128];

    Long64_t nbytes = 0, nb = 0;
    for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      // if (Cut(ientry) < 0) continue;
      time_1 = time-basetimeoffset;
      ic1_1 = ic1;
      ic2_1 = ic2;
      mw1_focusx_1 = mw1_focusx;
      mw1_focusy_1 =  mw1_focusy;
      mw2_focusx_1 = mw2_focusx;
      mw2_focusy_1 = mw2_focusy;
      mw1_posx_1 = mw1_posx;
      mw1_posy_1 = mw1_posy;
      mw2_posx_1 = mw2_posx;
      mw2_posy_1=   mw2_posy;
      eventid_1 = jentry;
      if (jentry % 1000 ==0) cout << jentry <<  " events completed..." << endl; 
      
      //calculate a rolling average of the signal
      arrayavg = 0;
      
      for (int i = 1; i<length;i++){
	array[i-1] = array[i];
      }
      if(ic1_1>-100) array[length-1] = ic1_1 ;
      
      for (int i = 0; i<length;i++){
	arrayavg += array[i]/double(length);
      }
      rollingavg = double(arrayavg);
      // cout << rollingavg << endl;
      
      //rerun the baseline fitter using the data in front of each spill
      if(rollingavg>0.5 && !beamonflag && jentry>5000) {
	beamonflag = true; beamon = 1;
	cout << "Beam On: " << jentry << endl;

	Baseline(true,1800,int(jentry-2000),true);
	//	bkgfunc->SetParameters( bkgfitPar0[128], bkgfitPar1[128], bkgfitPar2[128], bkgfitPar3[128]);
        
      }
      if(rollingavg <0.25 && beamonflag  && jentry>5000 )
	{
	beamonflag = false;
	beamon = 0;
	}



      if (ic1_1>0.01) {totaltime+=0.100;  totalcurrent+=ic1_1;}
      for (int ch = 4; ch < 128; ch++){
      	bkgfunc->SetParameters( bkgfitPar0[ch], bkgfitPar1[ch], bkgfitPar2[ch], bkgfitPar3[ch]); //set the parameters for the common mode subtraction
	channelamp[ch] = channels[ch]- baseline[ch] - bkgfunc->Eval(time_1); //subtract the baseline
	if (beamon==1&&ic1_1>0.01) th2d_beamSignal_channel->Fill(ch,channelamp[ch]);
	if (beamon==0&&ic1_1<=0.01) th2d_bkg_channel->Fill(ch,channelamp[ch]);

      }
      //      calculate mean and integral
      beamSignal_1 = 0.;
      numoverthresh=0;
      maxchannelamp = 0.;
      threshold = 20.;
      maxchannel = 0;
      //  cout << (time_1-basetimeoffset)<< endl;
      for (int ch = 4; ch < 64; ch++){
	//	bkgfunc->SetParameters( bkgfitPar0[ch], bkgfitPar1[ch], bkgfitPar2[ch], bkgfitPar3[ch]); //set the parameters for the common mode subtraction
	//	channelamp[ch] = channels[ch]- baseline[ch] - bkgfunc->Eval(time_1); //subtract the baseline
	//channelamp_smooth[ch] = channelamp[ch];//to be smoothed later
	if (channelamp[ch] < 0.) {abschannelamp[ch] = 0.;}
	else {abschannelamp[ch] = channelamp[ch];}

	if (ic1_1>0.01) th2d_beamSignal_channel->Fill(ch,channelamp[ch]);
	if (ic1_1<=0.01) th2d_bkg_channel->Fill(ch,channelamp[ch]);

	if (	channelamp[ch] >threshold) {
	  numoverthresh++;
	  if ( 	channelamp[ch] > maxchannelamp) {
	    maxchannelamp = 	channelamp[ch] ;
	    maxchannel = ch;
	  }
	}
	channellist_gsl[ch] = 0;
	channelamp_smooth[ch] = 0.;
      }
      numoverthresh+=4;
      numtocalc = 0;
      sidenumtocalc = 0;
      for (int ch = maxchannel-numoverthresh/2 ; ch < maxchannel + numoverthresh/2; ch++){
	if (ch>=4 && ch<=63){
	beamSignal_1 +=  channelamp[ch]; //integrate the signal around the peak channel
	channelamp_smooth[numtocalc] =channelamp[ch] ;
	channellist_gsl[numtocalc] = ch;
	numtocalc++;
	}
      }
      beamSidebandNoise_1 = 0.;
      for (int ch = 4; ch <  maxchannel-numoverthresh/2;  ch++){
	if (ch>=4 && ch<=63){
	  beamSidebandNoise_1 +=  channelamp[ch]; //integrate the noise outside the peak
	  sidenumtocalc++;
	}
      }
      for (int ch = maxchannel+numoverthresh/2; ch <  64;  ch++){
	if (ch>=4 && ch<=63){
	  beamSidebandNoise_1 +=  channelamp[ch]; //integrate the noise outside the peak
	  sidenumtocalc++;
	}
      }
      beamSidebandNoise_1 = beamSidebandNoise_1 /double(sidenumtocalc)*double(numtocalc);

      //  	if (ic1_1<=5 && ic2_1<=5 && jentry<1000) {	  cout << beamSignal_1 << " "<<(time_1)<<" "  << bkgfunc->Eval((time_1)) << endl;	}
      //  beamSignal_1 -= bkgfunc->Eval(time_1);

       //statistical analysis of the beam position
      // FFTsmoothing(channelamp_smooth);//Fourier transform smoothing



      beamPosX_1 = gsl_stats_wmean(channelamp_smooth,1,channellist_gsl,1,numtocalc); //calculate the weighted mean
      beamFocusX_1 = gsl_stats_wsd_with_fixed_mean(channelamp_smooth,1,channellist_gsl,1,numtocalc,beamPosX_1); //SD
      beamSkewX_1 = gsl_stats_wskew_m_sd(channelamp_smooth,1,channellist_gsl,1,numtocalc,beamPosX_1,beamFocusX_1); //skewness (symmetry)
      beamKurtX_1 = gsl_stats_wkurtosis_m_sd(channelamp_smooth,1,channellist_gsl,1,numtocalc,beamPosX_1,beamFocusX_1); //excess kurtosis (well behaved tails)
      beamNumX_1 = numtocalc;
      beamFocusX_1 *=2.3548;//SD-->FWHM


      if (ic1_1>0.01 && beamSignal_1>5)  th2d_mw1_beamPosX->Fill(mw1_posx_1,beamPosX_1);
      if (ic1_1>0.01 && beamSignal_1>5)  th2d_ic1_beamSignal->Fill(ic1_1,beamSignal_1);
      if (ic1_1>0.01 && beamSignal_1>5)  th1d_beamSignal_ic1ratio->Fill(beamSignal_1/ic1_1);


      //fit with a gaussian function;
      if (beamon==1&&ic1_1>0.05){
	if (numtocalc>10&&numtocalc<50){
	gausfunc->SetParameters(beamSignal_1/(sqrt(2)*beamFocusX_1/2.3548),beamPosX_1,beamFocusX_1/2.3548,0.);
	gausfunc2->SetParameters(beamSignal_1/(sqrt(2)*beamFocusX_1/2.3548),128-beamPosX_1,beamFocusX_1/2.3548,0.);

	}
	else{
	  gausfunc->SetParameters(500.,32.,10.,0.);
	  gausfunc2->SetParameters(500.,92.,10.,0.);
	}

	
	gausfunc->SetParLimits(0,0.,10000.);
	gausfunc->SetParLimits(1,10.,50.);
	gausfunc->SetParLimits(2,0.5,30.);
	gausfunc->SetParLimits(3,-150.,150.);

	gausfunc2->SetParLimits(0,0.,10000.);
	gausfunc2->SetParLimits(1,64,124);
	gausfunc2->SetParLimits(2,0.5,30.);
	gausfunc2->SetParLimits(3,-150.,150.);

	    gausgraph = new TGraphErrors(128,channellist,channelamp,errorx,errory);
	    lastfit = gausgraph->Fit(gausfunc,"QRN","",4,60); // single gaussian fit
	    lastfit2 = gausgraph->Fit(gausfunc2,"QRN","",68,124);// second peak

	    beamPosX_fit = gausfunc->GetParameter(1);
	    beamFocusX_fit =2.3548* gausfunc->GetParameter(2);
	    beamChi2_fit = gausfunc->GetChisquare()/gausfunc->GetNDF();
	    beamPeakX_fit = gausfunc->GetParameter(0);

	    beamPosX_fit2 = gausfunc2->GetParameter(1);
	    beamFocusX_fit2 =2.3548* gausfunc2->GetParameter(2);
	    beamChi2_fit2 = gausfunc2->GetChisquare()/gausfunc2->GetNDF();
	    beamPeakX_fit2 = gausfunc2->GetParameter(0);
	   
	    for (int ch = 4; ch < 64; ch++){
	      th2d_fitdiff_channel->Fill(ch, double(channelamp[ch]-gausfunc->Eval(ch)));
	    }
	    for (int ch = 64; ch < 128; ch++){
	      th2d_fit2diff_channel->Fill(ch, double(channelamp[ch]-gausfunc2->Eval(ch)));
	    }
	    
      }
      
      
      newdata->Fill();

    }
    TVectorD v(10);
    v[0] = totalcurrent;
    v[1] = totaltime;
    v.Write("icinfo");
   th2d_mw1_beamPosX->Write();
   th2d_ic1_beamSignal->Write();
   th1d_beamSignal_ic1ratio->Write();
   th2d_beamSignal_channel->Write();
   th2d_bkg_channel->Write();
   th2d_fitdiff_channel->Write();
   th2d_fit2diff_channel->Write();
     newdata->Write();
}


int main(int argc, char **argv)
{
  // Working directories
  const char *dirname = "/work/leverington/beamprofilemonitor/hitdata/HIT_01-12-2016/with_timestamp/";
  const char *pin_dirname = "/work/leverington/beamprofilemonitor/hitdata/HIT_01-12-2016/with_timestamp/pin/";
  const char *ext = ".root";
  TSystemDirectory pin_dir(pin_dirname, pin_dirname);

  if (true)
  {
    TSystemFile* file;
    TString fname = argv[1];
    

    
    // fname = file->GetName();

      // execute single PiN_run***.root
      if ( fname.EndsWith(ext) && !fname.BeginsWith("SAVE") && fname.BeginsWith("PiN_run"))
      {
        analyse2 *mon = new analyse2();
        printf("File name: %s \n", fname.Data());

        // Main part
        // Initialize(DIRName, FileName, baselineEvents, prelimEvents, beamLevel, firstFibermat, readOutFrequency in Hz, integrationTime in us)
        mon->Initialize(dirname, fname.Data(), 3200, 10000, 1., true, 1000., 0.);
	mon->Baseline(true,1000,1,true);
	mon->Loop(); //analysis loop

	//	cout << th2d_mw1_beamPosX->GetCorrelationFactor() << endl;
	mon->Save();//save output tree file
	//	mon->Close();//close output tree file

        delete mon;
      
      }
  
  }
  return 0;
}