#include "Rtypes.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "MultiWF.h"
#include <string.h>
#include "TCanvas.h"
#include "TF1.h"

ClassImp(MultiWF)

 MultiWF::MultiWF(Char_t *name, Float_t newt0,Int_t type)
{
  // constructor 
  //	 Char_t *name  ;  file name
  //	 Float_t newt0 ;  time shift to correct time scale 

 Int_t ok,i=0; 
 Float_t *buf=new Float_t[400];
 Date=TArrayF(5);
 if((in=fopen(name,"rb"))==NULL) {printf("Can't read file!!!!"); return;};
 

switch(type)
   {
   case 0:
  HeaderSize=44;
  ok=fread(buf,sizeof(Float_t),11,in);
  swooip(buf,11);
  Date[0]=buf[0]; Date[1]=buf[1]; Date[2]=buf[2]; Date[3]=buf[3]; Date[4]=buf[4];
  t0=buf[5]; deltat=buf[6]; Points=(Int_t)buf[7]; Voltage=buf[8]; Current=buf[9]; Temperature=buf[10];
  Voltage=TArrayF(1); 
  Temperature=TArrayF(1);
  Current=TArrayF(1);
  NumU=1;
  Voltage[0]=buf[8]; Current[0]=buf[9]; Temperature[0]=buf[10];
  break;
   case 1:
  HeaderSize=16;   
  ok=fread(buf,sizeof(Float_t),4,in); 
  swooip(buf,4); 
  for(i=0;i<5;i++) Date[i]=0;
  deltat=buf[1];
  Voltage=TArrayF(1);
  Temperature=TArrayF(1);
  Current=TArrayF(1);
  Voltage[0]=0; Current[0]=0; Temperature[0]=0;
  t0=buf[0]+deltat*buf[2];
  Points=(Int_t)(buf[3]-buf[2]); 
  NumU=1;
  break;
   case 2:
  ok=fread(buf,sizeof(Float_t),10,in);
  swooip(buf,11);
  Date[0]=buf[0]; Date[1]=buf[1]; Date[2]=buf[2]; Date[3]=buf[3]; Date[4]=buf[4];
  t0=buf[5]; deltat=buf[6]; Points=(Int_t)buf[7]; NumU=(Int_t) buf[9]; NumUAver=(Int_t)buf[8]; 
  HeaderSize=40+3*NumU*sizeof(Float_t);
  Voltage=TArrayF(NumU);
  Temperature=TArrayF(NumU);
  Current=TArrayF(NumU);
  //  Info();
  printf("kuracn");
  ok=fread(buf,sizeof(Float_t),NumU*3,in);
  printf("kurac ok=%dn",ok);
  swooip(buf,3*NumU);
  for(i=0;i<NumU;i++) {Voltage[i]=buf[i]; Current[i]=buf[i+NumU]; Temperature[i]=buf[i+NumU*2];}

   break;
   }
  
  if(newt0!=-1111) {t0=newt0*1e-9; cm_mint=t0; cm_maxt=0; printf("User start time    : %en",t0);}
  chis=new TH1F("Event View","Event View",Points,t0*1e9,(t0+Points*deltat)*1e9);
  delete buf;
  CheckSize();
  CCmode=0;
  cevent=1;
  printf("DATE= %4.0f %4.0f %4.0f %4.0f %4.0fn",Date[0],Date[1],Date[2],Date[3],Date[4]);
  printf("Points per event   : %dn",Points);
  printf("Time between points: %en",deltat);
  printf("Start time         : %en",t0);
  printf("Num. of U [%d each]: %dn",NumUAver,NumU);
    printf("Voltage            : %4.3f,",Voltage[0]);
  for(i=1;i<NumU;i++) printf("%4.3f,",Voltage[i]);
  printf("nCurrent            : %4.3e,",Current[0]);
  for(i=1;i<NumU;i++) printf("%4.3e,",Current[i]);
  printf("nTemperature        : %4.3f,",Temperature[0]);
  for(i=1;i<NumU;i++) printf("%4.3f,",Temperature[i]);
  printf("n");
  cmcut=1e6;
  Weight=1;
}


 MultiWF::~MultiWF()
{ 
  //destructor
if(in!=NULL) fclose(in);
if(pt!=NULL)  delete pt;
if(chis!=NULL) delete chis;
Clear();
}

 void  MultiWF::swoo(char *a, char *b) {
  // byte swaping (LABVIEW,HPUX g++)<->(LINUX g++, WINNT cl) 
  char c = *a;
  *a = *b; *b = c;
}

 void MultiWF::swooip(float *in, int s) {
 // byte swaping (LABVIEW,HPUX g++)<->(LINUX g++, WINNT cl) 
  char *sr, b;
  while(s--) {
    sr=(char *)in;
    swoo(&sr[0], &sr[3]);
    swoo(&sr[1], &sr[2]);
    in++;
  }
}

 TH1F *MultiWF::AverPlot(Int_t num)
{
  Int_t inc=Events/num;
  TH1F *his;
  TH1F *plot=new TH1F("aveplot","aveplot",num,0,Events);
  for(int i=1;i<=num;i++)
    {
      his=Average((i-1)*inc+1,i*inc);// printf("i=%d, %f",k, his->GetMaximumBin());
      his->SetAxisRange(his->GetBinCenter(1),his->GetBinCenter(1024),"X");
      plot->SetBinContent(i, his->GetMaximum());
      delete his;
    }
  return plot;
}


 Int_t MultiWF::CheckSize()
{
  // Retruns estimated number of events in the file
Int_t curpos=ftell(in);
fseek(in,0,SEEK_END);
Int_t filesize=ftell(in);
rewind(in);
fseek(in,curpos,SEEK_CUR);
Events=(filesize-HeaderSize)/(sizeof(float)*Points);
if((filesize-HeaderSize)%(sizeof(float)*Points)==0) 
  {printf("Number of events in file=%dn",Events); if(NumU==1) NumUAver=Events; return 0;} 
 else 
  {printf("Sucpicious file!!!!! Around %d events in the file!n",Events); return 1;}; 

}


 void MultiWF::GetNextEvent(Int_t CM)
{
  // Gets next event in the file 
  //	         Int_t CM   ; ignore=0  ,  correct common mode=1 (this is base line shift) 
  Int_t ok,i; 
  Float_t *buf=new Float_t[Points];
  ok=fread(buf,sizeof(Float_t),Points,in);
  swooip(buf,Points);
  for(i=1;i<=Points;i++){ chis->SetBinContent(i,Weight*buf[i-1]);}// printf("%f ",buf[i-1]);} printf("n");
  if(CM) CommonMode();
  cevent++;
  delete buf;
}

 void MultiWF::CommonMode()
{
  // Corrects base line shift
Int_t i=0,num=chis->GetNbinsX();
Int_t right,left;
Double_t cm;
right=chis->GetXaxis()->FindBin(cm_maxt);
left=chis->GetXaxis()->FindBin(cm_mint);
cmode=chis->Integral(left,right);
cmode=cmode/((Float_t )right-(Float_t )left);
//printf("left bin=%d , right bin=%d , common mode=%f n",left,right,cm);
for(i=1;i<num;i++) chis->SetBinContent(i,chis->GetBinContent(i)-cmode);  	

}

 Float_t MultiWF::Integral(Float_t mint,Float_t maxt)
{
  //Calculates integral of the current histogram (chis)!
  //		Float_t mint	; integration time interval
  //		Float_t maxt	;
Int_t left,right;
if(mint==-1111) left=1; else left=chis->GetXaxis()->FindBin(mint);
if(maxt==-1111) right=chis->GetNbinsX(); else right=chis->GetXaxis()->FindBin(maxt); 
//printf("left bin=%d , right bin=%dn",left,right);
 return(chis->Integral(left,right)*deltat*1e9);
}

 void MultiWF::EShaping(Elec *el,TH1F *his,Int_t what)
{
  //copies the current histogram in TH1F *his and shapes it with Elec *el 
  // Int_t what; 0-> preamplifier->CR-RC^2 shaper ,  1->deconvolution for the TCT transfer function
  
chis->Copy(*his);
Double_t xmax,xmin;
switch(what)
  {
case 0: xmax=his->GetXaxis()->GetXmax();
        xmin=his->GetXaxis()->GetXmin();
        his->GetXaxis()->SetLimits(xmin*1e-9,xmax*1e-9);
        el->preamp(his);
        el->CRshape(his);
        el->RCshape(his);
        el->RCshape(his); 
	his->GetXaxis()->SetLimits(xmin*1e9,xmax*1e9);
	break;
case 1: el->Revpreamp(his,1e9); break;
  }
}	

 void MultiWF::WF(Int_t start,Int_t end,Int_t step ,Float_t stime,Float_t etime,Int_t aver)
{
  // Display waveform dependence on time of the measurement
  //			Int_t start   ;   the first event to be displayed
  //			Int_t end     ;   the first event to be displayed
  //                    Int_t step    ;   step size between events
  //			Float_t stime ;   time window of display
  //			Float_t etime ;   time window of display
  //			Int_t aver    ;   number of events to average first event 
  //					  Average(start,start+aver) etc.
TH1F *his;
Int_t colori[]={1,2,3,4,5,6,7,13,28,30,34,38,40,31,46,49,1,2,3,4,5,6,7};
Int_t color,k=0;
  for(Int_t i=start;i<=end;i+=step)
    {
  his=Average(i,i+aver); 
  color=colori[k];
  his->SetLineColor(color);
  his->SetLineWidth(2);
  if(i==start) 
    {
    his->GetXaxis()->SetRange(his->GetXaxis()->FindBin(stime),his->GetXaxis()->FindBin(etime));
    his->SetTitle("Current shape dependence on time");
    his->SetXTitle("t[ns]");
    his->SetYTitle("I [V/50#Omega]");
    his->DrawCopy(); 
    Legend(his,start,end,step);
    } 
  else his->DrawCopy("SAME");
  k++;
    }
}

 void MultiWF::Legend(TH1F *ch,Int_t start, Int_t end,Int_t Step)
{
  //Draws legend
Float_t minx,miny,maxy,maxx,x1,x2,y1,y2;
TString title,utit="t=";
Char_t v[7]; 
Int_t color,cii=0;;
Int_t colori[]={1,2,3,4,5,6,7,13,28,30,34,38,40,31,46,49,1,2,3,4,5,6,7};
TText *text;

minx=ch->GetXaxis()->GetBinCenter(ch->GetXaxis()->GetFirst());
maxx=ch->GetXaxis()->GetBinCenter(ch->GetXaxis()->GetLast());
miny=ch->GetMinimum();
maxy=ch->GetMaximum();

x1=(maxx-minx)*0.6+minx; 
x2=(maxx-minx)*0.9+minx;
y2=(maxy-miny)*0.35+miny;
y1=(maxy-miny)*0.95+miny;

//printf("coords: x1=%f y1=%f x2=%f y2=%fn",x1,y1,x2,y2);


if(pt!=NULL) delete pt;
pt=new TPaveText(x1,y1,x2,y2);

for(Int_t i=start;i<=end;i+=Step)
  { 
   i2a(v,(Int_t) i,0); 
   //   color=i/7*40+i%7+1; 
   color=colori[cii];
   title=utit+v; title=title+" ";
   text=pt->AddText((const char *)title);
   text->SetTextColor((Color_t)color);
   text->SetTextSize(0.05);
   cii++;
  }
pt->Draw();

}



 TH1F *MultiWF::Time(Float_t time,Float_t low, Float_t high,Int_t rebin,Float_t stime)
{
  // Used with known trigger interval period
  // Each event charge (integral) is drawn vs. time of its occurence
  //			Float_t time ; 1/frequency of repetition
  //			Float_t low  ; integration time interval
  //			Float_t high ;
  //			Int_t rebin  ; rebining factor (average charge is shown)
  //			Float_t stime ; start time of the histogram (used with delayed DAQ)
Int_t i,j;
Float_t sum=0;
Float_t etime=time*Events+stime;
TH1F *th=new TH1F("th","CCE vs. Time",Events/rebin,stime,etime);
if(low==-1111) high=chis->GetBinCenter(Points-1);
for(i=0;i<Events/rebin;i++)
  {
for(j=0;j<rebin;j++) {GetNextEvent(CCmode); sum+=Integral(low,high);} 
     sum/=rebin;
     th->SetBinContent(i,sum);
     sum=0;
  }

 th->SetXTitle("time[s]");
 th->SetYTitle("Charge [arb.]");
 th->SetLabelSize(0.045,"X");
 th->SetLabelSize(0.045,"Y"); 
 rewind(in);
 fseek(in,11*sizeof(Float_t),SEEK_SET);
 cevent=1;
return th;
}

 TH1F *MultiWF::Specter(Float_t low, Float_t high,Float_t weight,Int_t num,Float_t first,Float_t last)
{
  // Specter of all recorded events 
  //			Float_t low   ; integration time interval
  //			Float_t high  ;
  //			Float_t weight; normalization
  //			Int_t num     ; number of bins
  //			Int_t first   ; histogram limits
  //			Int_t last    ; 
  
Int_t i;
TH1F *spectrum=new TH1F("Spectrum","Spectrum",num,first,last);
if(low==-1111) high=chis->GetBinCenter(Points-1);
for(i=0;i<Events;i++)
  {
    GetNextEvent(CCmode);
    spectrum->Fill(Integral(low,high)*weight);
  }
 spectrum->SetXTitle("Charge [arb.]");
 spectrum->SetYTitle("Entries");
 spectrum->SetLabelSize(0.045,"X");
 spectrum->SetLabelSize(0.045,"Y"); 
 rewind(in);
 fseek(in,11*sizeof(Float_t),SEEK_SET);
 cevent=1;
return spectrum;
}

 TH1F *MultiWF::CMHisto()
{
TH1F *cmodehis=new TH1F("CModeHis","CModeHis",200,-100e-3,100e-3);
Rewind();
for(Int_t i=0;i<Events;i++)
  {
    GetNextEvent(CCmode);
    cmodehis->Fill(cmode);    
  }
 return cmodehis;
}

 TH1F *MultiWF::PHSpecter(Float_t weight,Int_t num,Float_t first,Float_t last,Float_t time_lo,Float_t time_hi,Int_t mode)
{
  // Pulse height specter of all recorded events 
  //			Float_t weight; normalization
  //			Int_t num     ; number of bins
  //			Int_t first   ; histogram limits
  //			Int_t last    ; 

Int_t i,j;
Int_t timelo;
Int_t timehi;
Float_t val,max=-1e6;
Char_t xname[100];



TH1F *spectrum=new TH1F("Spectrum","Spectrum",num,first,last);

if(time_lo==-1e6) timelo=1; 
else timelo=chis->GetXaxis()->FindBin(time_lo);
if(time_lo==-1e6) timehi=chis->GetNbinsX(); else 
timehi=chis->GetXaxis()->FindBin(time_hi);
// TF1 *maxf=new TF1("maxf","pol2",time_lo,time_hi);
 TF1 *maxf=new TF1("maxf","pol0",time_lo,time_hi);
printf("Time window selected: %d %d (bins)n",timelo,timehi);
//		TCanvas *c1=new TCanvas("c1","c1",800,600);
for(i=0;i<Events;i++)
  {
    GetNextEvent(CCmode);
    if(TMath::Abs(cmode)<cmcut)
      {
	    switch(mode)
	      {
	      case 0:	
        	for(j=timelo;j<=timehi;j++) 
	        {
                val=chis->GetBinContent(j);
	        max=max>val?max:val;
		//		printf("%d %d %f %fn",i,cevent-1,max,cmode);
		}
	        break;
	      case 1:
		
		chis->Fit("maxf","RNQ");
		//			max=-TMath::Power(maxf->GetParameter(1),2)/(4*maxf->GetParameter(2))+maxf->GetParameter(0);
		max=maxf->GetParameter(0);
		//	
		//		printf("max=%fn",max);
		//		c1->Update();
		//		while(!getchar());
	        break;
	      }
	    //	    chis->SetAxisRange(0,100);
	   
	    //	    if(chis->GetMaximumBin()==timelo)
	    // if(max>0.11) printf("%f %dn",max,i);
	spectrum->Fill(max*weight);
      }  
    //    spectrum->Fill(chis->GetMaximum()*weight);
    max=-1e6;
  }
sprintf(xname,"Signal [%4.1f mV]",1000./weight); 
 spectrum->SetXTitle(xname);
 spectrum->SetYTitle("Entries");
 spectrum->SetLabelSize(0.045,"X");
 spectrum->SetLabelSize(0.045,"Y"); 
 rewind(in);
 fseek(in,11*sizeof(Float_t),SEEK_SET);
 cevent=1;
return spectrum;
}

 TH1F *MultiWF::WaveHisto(Float_t weight,Int_t num,Float_t first,Float_t last,Int_t startbin,Int_t endbin)
{
  // Pulse height specter of all recorded events 
  //			Float_t weight; normalization
  //			Int_t num     ; number of bins
  //			Int_t first   ; histogram limits
  //			Int_t last    ; 
  //                    Int_t startbin; start bin of the waveform histogram
  //                    Int_t endbin;   end bin of the waveform histogram
Int_t i,j;
Char_t xname[100];
TH1F *spectrum=new TH1F("Spectrum","Spectrum",num,first,last);
 if(startbin==-1) startbin=0;
 if(endbin==-1) endbin=Points;
for(i=0;i<Events;i++)
  {
    GetNextEvent(CCmode);
    for(j=startbin;j<endbin;j++) spectrum->Fill(chis->GetBinContent(j)*weight);
  }
sprintf(xname,"Signal [%4.1f mv]",1000./weight); 
 spectrum->SetXTitle(xname);
 spectrum->SetYTitle("Entries");
 spectrum->SetLabelSize(0.045,"X");
 spectrum->SetLabelSize(0.045,"Y"); 
 rewind(in);
 fseek(in,11*sizeof(Float_t),SEEK_SET);
 cevent=1;
return spectrum;
}


 void  MultiWF::SetNewFile(Char_t *name)
{
  if(in!=NULL) fclose(in);
  if((in=fopen(name,"rb"))==NULL) {printf("Can't read file!!!!"); return;};
  cevent=0;
}

 void MultiWF::Rewind(Int_t index)
{
  //Rewinds the file poiter to the desired event
  //		      Int_t index ; event number (defualt=0, begining of the file)
  if(cevent>index); {rewind(in);  fseek(in,HeaderSize,SEEK_SET); cevent=1;}
  if(cevent==index) return;
  while(cevent<index) GetNextEvent(CCmode);
}

 TH1F *MultiWF::Average(Int_t start, Int_t end)
{
  // Averaging events (Waveforms)
  //		Int_t start ; first event
  //		Int_t end   ; last event
  //		Function returns pointer to averaged histogram
Int_t i,ANum=0;
Rewind(start);
TH1F *ah=new TH1F("AverHis","Averaged Histograms",chis->GetNbinsX(),chis->GetXaxis()->GetXmin(),chis->GetXaxis()->GetXmax());
 for(i=start; i<=end; i++) 
   { 
     GetNextEvent(CCmode); 
     //To be written - filter function eliminating the 
     ah->Add(chis); 
     ANum++;
   }
ah->Scale(1/(Float_t) ANum);
return ah;
}

 TH1F *MultiWF::Add(Int_t start, Int_t end)
{
  // Averaging events (Waveforms)
  //		Int_t start ; first event
  //		Int_t end   ; last event
  //		Function returns pointer to averaged histogram
Int_t i;
if(end==0) end=Events-1; 
if(start==0) start=1;
Rewind(start);
TH1F *ah=new TH1F("sum","Cumulative histogram",chis->GetNbinsX(),chis->GetXaxis()->GetXmin(),chis->GetXaxis()->GetXmax());
 for(i=start; i<=end; i++) { GetNextEvent(CCmode); ah->Add(chis); }
 printf("Total entries in the histogram : %dn",ah->Integral());
return ah;
}



 void MultiWF::Info()
{
  Int_t i=0;
  // Prints information about measurement
  printf("DATE= %4.0f %4.0f %4.0f %4.0f %4.0fn",Date[0],Date[1],Date[2],Date[3],Date[4]);
  printf("Points per event   : %dn",Points);
  printf("Time between points: %en",deltat);
  printf("Start time         : %en",t0);
  printf("Num. of U [%d each]: %dn",NumUAver,NumU);
  printf("Voltage            : %4.3f,",Voltage[0]);
  for(i=1;i<NumU;i++) printf("%4.3f,",Voltage[i]);
  printf("nCurrent            : %4.3e,",Current[0]);
  for(i=1;i<NumU;i++) printf("%4.3e,",Current[i]);
  printf("nTemperature        : %4.3f,",Temperature[0]);
  for(i=1;i<NumU;i++) printf("%4.3f,",Temperature[i]);
  printf("n");
  printf("*********************************n");

  printf("*********************************n");
  printf("Current Event      : %dn",cevent);
  printf("Common mode time   : [%4.2f, %4.2f] n",cm_mint,cm_maxt);
}



 Int_t MultiWF::i2a(Char_t *v,Float_t vol,Int_t sign)
{
  // int to ascii with signs (same as in MeasureWF)
Int_t k=0;
Int_t dec=10000;
       if(sign) if(vol>0) v[k++]='+'; else {v[k++]='-'; vol=-vol;}
       if(!sign) vol=vol>0?vol:-vol; 
     for(Int_t j=0;j<=2;j++) { if((Int_t) vol/dec!=0) v[k++]=(Char_t)(((Int_t) vol)/dec)+48; dec=dec/10;}
       v[k++]=(Char_t)(((Int_t) vol%100)/10)+48;
       v[k++]=(Char_t)((Int_t) vol%10)+48; 
       v[k++]='0';
       return k-1;
}


 void MultiWF::FitShape(Int_t order,Float_t peaking_time)
{
  TF1 *fitf=new TF1("filter","[0]*((x-[1])*[2]/[3])^[2]*exp(-(x-[1])*[2]/[3])",0,1500);
  fitf->SetParameter(0,0.2);
  fitf->SetParameter(1,200);
  fitf->SetParameter(2,(Float_t )order);
  fitf->SetParLimits(2,-1,-1);
  fitf->SetParameter(3,peaking_time); 
  fitf->SetParLimits(3,-1,-1);
  chis->Fit("filter","R");
}


Double_t langaufun(Double_t *x, Double_t *par) {

   //Fit parameters:
   //par[0]=Width (scale) parameter of Landau density
   //par[1]=Most Probable (MP, location) parameter of Landau density
   //par[2]=Total area (integral -inf to inf, normalization constant)
   //par[3]=Width (sigma) of convoluted Gaussian function
   //
   //In the Landau distribution (represented by the CERNLIB approximation), 
   //the maximum is located at x=-0.22278298 with the location parameter=0.
   //This shift is corrected within this function, so that the actual
   //maximum is identical to the MP parameter.

      // Numeric constants
      Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
      Double_t mpshift  = -0.22278298;       // Landau maximum location

      // Control constants
      Double_t np = 300.0;      // number of convolution steps
      Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas

      // Variables
      Double_t xx;
      Double_t mpc;
      Double_t fland;
      Double_t sum = 0.0;
      Double_t xlow,xupp;
      Double_t step;
      Double_t i;


      // MP shift correction
      mpc = par[1] - mpshift * par[0]; 

      // Range of convolution integral
      xlow = x[0] - sc * par[3];
      xupp = x[0] + sc * par[3];

      step = (xupp-xlow) / np;

      // Convolution integral of Landau and Gaussian by sum
      for(i=1.0; i<=np/2; i++) {
         xx = xlow + (i-.5) * step;
         fland = TMath::Landau(xx,mpc,par[0]) / par[0];
         sum += fland * TMath::Gaus(x[0],xx,par[3]);

         xx = xupp - (i-.5) * step;
         fland = TMath::Landau(xx,mpc,par[0]) / par[0];
         sum += fland * TMath::Gaus(x[0],xx,par[3]);
      }

      return (par[2] * step * sum * invsq2pi / par[3]);
}




 MeasureWF *MultiWF::GetWF()
{
  Float_t Delta;
  Int_t i=0,ix;
  MeasureWF *MWF=new MeasureWF(NumU);
  TH1F *his;
    for(i=0;i<NumU;i++)
        {
	his=Average(i*NumUAver+1,(i+1)*NumUAver);
	MWF->AddHisto(i,GetU(i),his);
	}  
    return MWF;
}