#include "GSignals.h"

ClassImp(GSignals)


 GSignals::GSignals()
{
Num=10;
Noise=16;
MIPSignal=200;
PeakTimes=TArrayF(Num);
Cuts=TArrayF(Num);	
sn =new TClonesArray("TH1F",Num);	
csignal =new TClonesArray("TH1F",Num);
csignalsn =new TClonesArray("TH1F",Num);
csize =new TClonesArray("TH1F",Num);
}



 GSignals::GSignals(Int_t num,Float_t noise, Float_t mipsignal, Float_t *peaktimes, Float_t *cuts) 
{

  TString snname="sn",csignalname="Signal_histogram_", csignalsnname="signal_SN_", csizename="Cluster_size_";
  TString name1,name2,name3,name4;
   Num=num;
   Noise=noise;
   MIPSignal=mipsignal;
   PeakTimes=TArrayF(Num,peaktimes);
   Cuts=TArrayF(3,cuts);	
   sn =new TClonesArray("TH1F",Num);	
   csignal =new TClonesArray("TH1F",Num);
   csignalsn =new TClonesArray("TH1F",Num);
   csize =new TClonesArray("TH1F",Num);

   TClonesArray &sncopy = *sn;
   TClonesArray &csignalcopy = *csignal;	
   TClonesArray &csignalsncopy = *csignalsn;	
   TClonesArray &csizecopy = *csize;

   rnd=TRandom(44);

   for(Int_t i=0;i<Num;i++)  
     { 
       char v[2]; v[0]=(Char_t) i+48; v[1]='\0';
       //printf("character %sn",v);
      name1=snname+v; name2=csignalname+v; name3=csignalsnname+v; name4=csizename+v;
  
       new(sncopy[i]) TH1F((const char *)name1,"SN",100,0,50);
       new(csignalcopy[i]) TH1F((const char *)name2,"Signal histogram",200,0,1e3);
       new(csignalsncopy[i]) TH1F((const char *)name3,"Signal after SN-cuts",160,0,80);
       new(csizecopy[i]) TH1F((const char *)name4,"Cluster size",10,0,10);	
     }

}


 GSignals::~GSignals()
{
sn->Delete();
csignal->Delete();
csignalsn->Delete();
csize->Delete();
}

 void GSignals::DrawHistogram(Int_t i,char *option,char *doption)
{
  if(! strcmp(option,"s")) ((TH1F *)csignal->At(i))->DrawCopy(doption);
  if(! strcmp(option,"sn")) ((TH1F *)sn->At(i))->DrawCopy(doption); 
  if(! strcmp(option,"ssc")) ((TH1F *)csignalsn->At(i))->DrawCopy(doption);
  if(! strcmp(option,"csize")) ((TH1F *)csize->At(i))->DrawCopy(doption);
}

 void GSignals::GetHistogram(Int_t i,TH1F *his,char *option)
{
  if(! strcmp(option,"s")) ((TH1F *)csignal->At(i))->Copy(*his);
  if(! strcmp(option,"sn")) ((TH1F *)sn->At(i))->Copy(*his);
  if(! strcmp(option,"ssc")) ((TH1F *)csignalsn->At(i))->Copy(*his);
  if(! strcmp(option,"csize")) ((TH1F *)csize->At(i))->Copy(*his);
}

 void GSignals::AddEvent(Int_t timenum, Int_t clsiz, Float_t *sig)
{
Float_t sumsig=0; //sum of all signals(over hit strips) in the event after addition of random noise
Float_t sumsn=0;  //sum of all SN(over hit strips) in the event after addition of random noise
Int_t cls=0;        // cluster size in the event
Int_t j=0;

TArrayF signal=TArrayF(clsiz,sig);
((TH1F *)csignal->At(timenum))->Fill(signal.GetSum()/1e4);

Float_t *SN=new Float_t[clsiz]; // SN
Float_t *SIG=new Float_t[clsiz]; // Signal
Float_t maxsn=0;
Int_t maxsni;
Int_t cent=0;

 for(j=0;j<clsiz;j++) 
   {
     SN[j]=(signal[j]/1e4+rnd.Gaus(0,Noise))/Noise;  //addition of random noise
     SIG[j]=(signal[j]/1e4+rnd.Gaus(0,Noise)); //addition of random noise
     if(SN[j]>maxsn) {maxsn=SN[j]; maxsni=j;} //max SN-strip  in event search
     //printf("SN[%d]=%f Signal=%fn",j,SN[j],signal[j]/1e4);
     //printf("event=%d index=%dn",ev,maxsni);
   }
 for(j=0;j<clsiz;j++) 
   {
     if(j!=maxsni) {if(SN[j]>Cuts[2]) {sumsn+=SN[j]; sumsig+=SIG[j]; cls++;}} else 
     if(SN[j]>Cuts[1]) {sumsn+=SN[j]; sumsig+=SIG[j]; cent=1; cls++;}
   }

if(cent==1 && sumsn>Cuts[0]) 
  { 
    ((TH1F *)sn->At(timenum))->Fill(sumsn);
    ((TH1F *)csize->At(timenum))->Fill(cls);
    ((TH1F *)csignalsn->At(timenum))->Fill(sumsig/10);
  }
 sumsn=0; sumsig=0; cls=0;
 delete [] SN;
 delete [] SIG;

}

 void GSignals::Clean()
{
  for(Int_t i=0;i<Num;i++)  
     { 
     ((TH1F *)sn->At(i))->Reset();
     ((TH1F *)csize->At(i))->Reset();
     ((TH1F *)csignal->At(i))->Reset();
     ((TH1F *)csignalsn->At(i))->Reset();
     }
}

 Double_t GSignals::FitLG(Int_t timenum,char *option, Double_t *par, Double_t *parerr)
{
Double_t GausLandConv(Double_t *,Double_t *);
Double_t param[4];
Double_t max=0;
TH1F *his=new TH1F();
GetHistogram(timenum,his,option);
param[0]=(Double_t) GetMax(his,0);
param[3]=(Double_t) GetMax(his,1);
param[2]=(Double_t) 3;//his->GetRMS(1);
param[1]=(Double_t) 1;
TF1 *f = new TF1("f",GausLandConv,0,70,4);
f->SetParameters(param);
his->Fit("f","RN");
max= FindLandauPeak(f);
 for(Int_t i=0;i<4;i++) {parerr[i]=f->GetParError(i); par[i]=f->GetParameter(i);}

delete his;
delete f;
return max;

}

 Double_t GSignals::FindLandauPeak(TF1 *land)
{
Double_t eps=0.01;
Double_t sigma=land->GetParameter(1);
Double_t mean=land->GetParameter(0);
Double_t x=mean-sigma;
while(land->Eval(x) > land->Eval(x-eps)) {x=x+eps;}  
//printf("Max Y=%f, Max X=%f, sigma=%f, mean=%fn",land->Eval(x-eps),x-eps,sigma,mean);
return x-eps;
}
Double_t GausLandConv(Double_t *x, Double_t *par) 
{
// Crude convolution of a Gaussian (mean= 0.) and a Landau
// distribution at x[0].
  Double_t HiBound = 80; //par[6]; // hi bound for integration
  Double_t LoBound = 0; //par[5]; // lo bound for integration
  Double_t NumIntBins = 320;// par[4]; // number of Riemann integration bins
Double_t YNorm = par[3]; // Normalization
Double_t GausSigma = par[2]; // sigma width of gaussian
Double_t XNorm = par[1]; // x axis normalization
Double_t LandMostProb = par[0]; // most probable of landau
Double_t StepSize = (HiBound - LoBound) / NumIntBins;


// Integrate over the region LoBound --> HiBound using 
// NumIntBins Riemannian bins

Double_t gaus, result = 0;
for ( Double_t tau = LoBound + StepSize/2; tau < HiBound; tau += StepSize ) {

// evaluate a gaussian curve at (x[0]-tau)
gaus = TMath::Exp( -1*(*x-tau)*(*x-tau)/(2*GausSigma*GausSigma) )/( TMath::Sqrt(2.*3.141592653)*GausSigma );

// add the area of this Riemann bin to the total integral result
result+=YNorm * TMath::Landau(tau,LandMostProb,XNorm) * gaus;
}
//result/=NumIntBins;
return result;
}


 Float_t GSignals::GetMax(TH1F *histo,Int_t opt) //X=0,Y=else
{
Int_t i,maxi=0;
Float_t max=-1e20;

 for(i=1;i<histo->GetNbinsX();i++) {//printf("%d %f %dn",i,max,maxi);
   if(max<histo->GetBinContent(i)) {max=histo->GetBinContent(i); maxi=i;}}
   if(opt==0) return(histo->GetBinCenter(maxi)); else return(max);
}