//=====================================================================
//	PndEmcWaveform.cxx
//
// 	Class to hold waveforms created from Emc Hits
//
//	 Software developed for the BaBar Detector at the SLAC B-Factory.
// Adapted for the PANDA experiment at GSI		
//
// 	P.D.Strother 	Imperial College
//      Naveen Gunawardane  Imperial College
// Dima Melnichuk - adaption for PANDA		
//-----------------------


#include "PndEmcWaveform.h"
#include "PndEmcHit.h"
#include "PndEmcCRRCPulseshape.h"
#include "PndEmcCR2RCPulseshape.h"
#include "PndEmcMapper.h"
#include "PndDetectorList.h" 

#include "TRandom.h"
#include "TClonesArray.h"

#include "TClass.h"
#include "TBuffer.h"
#include "Riostream.h"
//---------------
// C++ Headers --
//---------------

#include <iostream>
#include <algorithm>
#include "assert.h"

using std::cout;
using std::endl;


//----------------
// Constructors --
//----------------
Double_t PndEmcWaveform::BarrelOverlapTime = 700;//(120*12.5 - 800);ns
Double_t PndEmcWaveform::ForwardOverlapTime = 390;//(60*10 - 210);
Double_t PndEmcWaveform::ShashylikOverlapTime = 130;//(60*5.5 - 200);

PndEmcWaveform::PndEmcWaveform():
	fTrackId(-1)
	,fDetectorId(-1)
	,fWaveformLength(0)
	,fSignal(0,0.)
	,fSignalError(0,0.)
	,fHitIndex(-1)
	,fSampleRate(0.)
	, fBaselineValue(0.)
	,FairTimeStamp()
{}

PndEmcWaveform::PndEmcWaveform(int trackId, long detId, Double_t sampleRate, long waveform_length, Int_t hitIndex, Double_t time) :
	fTrackId(trackId)
	,fDetectorId(detId)
	,fWaveformLength(waveform_length)
	,fSignal(waveform_length,0.)
	,fSignalError(waveform_length,0.)
	,fHitIndex(hitIndex)
	,fSampleRate(sampleRate)
	, fBaselineValue(0.)
	,FairTimeStamp(time)
{
  if(hitIndex>=0) SetLink(FairLink("EmcHit", hitIndex));
}

PndEmcWaveform::PndEmcWaveform(int trackId, long detId, const std::vector<Double_t>& signal, Int_t hitIndex) : 
	fTrackId(trackId), 
	fDetectorId(detId),
	fWaveformLength(signal.size()),
	fSignal(signal),
	fHitIndex(hitIndex)
{
  if(hitIndex>=0) SetLink(FairLink("EmcHit", hitIndex));
}

PndEmcWaveform::PndEmcWaveform(long detId, const std::vector<Double_t>& signal, const FairMultiLinkedData& links) :
	fTrackId(-1),
	fDetectorId(detId),
	fWaveformLength(signal.size()),
	fSignal(signal),
	fHitIndex(-1) {
	SetLinks(links);
}

/* auto generated by compiler
PndEmcWaveform::PndEmcWaveform(const PndEmcWaveform& copy):
	fTrackId(copy.fTrackId)
	,fDetectorId ( copy.fDetectorId)
	,fWaveformLength(copy.fWaveformLength)
	,fHitIndex(copy.fHitIndex)
	,fSignal(copy.fSignal)
  ,fSignalError(copy.fSignalError)
	,fSampleRate(copy.fSampleRate)
  ,fBaselineValue(copy.fBaselineValue)
	,fEvt(copy.fEvt) 
  ,FairTimeStamp(copy)
{
}
*/

/*auto generated by compiler
PndEmcWaveform& PndEmcWaveform::operator=(const PndEmcWaveform &copy){
	if(this != &copy){
		fTrackId = copy.fTrackId;
		fDetectorId = copy.fDetectorId;
		fWaveformLength = copy.fWaveformLength;
		fHitIndex = copy.fHitIndex;
		fSignal = copy.fSignal;  
		fSignalError = copy.fSignalError;
		fEvt = copy.fEvt;
		fSampleRate = copy.fSampleRate;
		fBaselineValue = copy.fBaselineValue;
		FairTimeStamp::operator=(copy);
  }
  return *this;
}
*/

//--------------
// destructor --
//--------------

PndEmcWaveform::~PndEmcWaveform()
{
	fSignal.clear();
	fSignalError.clear();
	fEvt.clear();
}


//-------------
// selectors --
//-------------

Double_t 
PndEmcWaveform::GetScale(Double_t samplerate, PndEmcAbsPulseshape *pulseshape) const
{
	return GetNormalisation(samplerate, pulseshape);
}

Double_t 
PndEmcWaveform::GetNormalisation(Double_t samplerate, PndEmcAbsPulseshape *pulseshape) const
{
	// function to return the equivalent pulse height for a 1 gev pulse
	// used  in emchitstowaveform to determine electronics noise scale

	PndEmcWaveform newWaveform(*this);
	newWaveform.clearAndReset();

	PndEmcHit *gevhit=new PndEmcHit();
	gevhit->SetEnergy(1.0);
	gevhit->SetTime(0.);
	newWaveform.UpdateWaveform(gevhit, 0, false, 1., 0., samplerate, pulseshape);
	delete gevhit;

	Double_t maximum=newWaveform.Max();
	return maximum;
}


//-------------
// modifiers --
//-------------

void 
PndEmcWaveform::UpdateWaveform(PndEmcHit *hit
		, Double_t pePerMeV
		, Bool_t usePhotonStatistic
		, Double_t excessNoiseFactor
		, Double_t firstADCBinTime
		, Double_t sampleRate
		, PndEmcAbsPulseshape *pulseshape
		, Double_t EnergyError)
{
	Double_t energy=(Double_t)hit->GetEnergy();
	Double_t time= GetTimeStamp();//firstADCBinTime = 0.;

	MakeWaveform(energy, time, pePerMeV, usePhotonStatistic, excessNoiseFactor, firstADCBinTime, sampleRate, pulseshape, EnergyError);

}

void 
PndEmcWaveform::MakeWaveform(Double_t energy
		, Double_t time
		, Double_t pePerMeV
		, Bool_t usePhotonStatistic
		, Double_t excessNoiseFactor
		, Double_t firstADCBinTime
		, Double_t sampleRate
		, PndEmcAbsPulseshape *pulseshape
		, Double_t EnergyError)
{
	Double_t amplitude;
	Double_t photonStatFactor;
	if (usePhotonStatistic)
	{
		Double_t crystalPhotonsMeV = 1.0e3 * energy * pePerMeV;
		photonStatFactor = gRandom->Gaus(1,sqrt(excessNoiseFactor/crystalPhotonsMeV));
	}
	else
	{
		photonStatFactor=1.;
	}

	amplitude = energy*photonStatFactor;

	Double_t local_time(0.);
	Double_t time_offset=firstADCBinTime;//nano seconds

	Double_t sumCharge(0.);
	Double_t sumChargeErr(0.);

	for (Int_t i=0;i<fWaveformLength;i++)
	{
		local_time = i/sampleRate;//seconds
		fSignal[i]+=(pulseshape->value(local_time,amplitude,time_offset));
		sumCharge += fSignal[i];
		sumChargeErr += sqrt(fabs(fSignal[i]));
	}
	Double_t Coff = EnergyError*sumCharge/sumChargeErr;
	for(Int_t i=0;i<fWaveformLength;++i){
		fSignalError[i] = fSignal[i]*Coff;
	}

}


	void 
PndEmcWaveform::AddElecNoise(Double_t width)
{

	for (int i=0;i<fWaveformLength;i++)
	{
		Double_t ran_noise;
		ran_noise=gRandom->Gaus(0,1)*width;
		fSignal[i]+=ran_noise;
	}
}


	void 
PndEmcWaveform::Digitise(Double_t oneBitResolution)
{
	for (  int i=0;i<fWaveformLength;i++)
	{
		fSignal[i]=(Double_t) ( long (fSignal[i]/oneBitResolution+64)-64 ) * oneBitResolution;
	}
}


	void 
PndEmcWaveform::AddElecNoiseAndDigitise(Double_t noise_width,
		Double_t oneBitResolution, Double_t EnergyError)
{
	// Do both e_noise and digitisation.  

	for (  int i=0;i<fWaveformLength;i++)
	{
		Double_t ran_noise=gRandom->Gaus(0,1)*noise_width;
		fSignal[i]+=ran_noise;
		if (fSignal[i]<0) fSignal[i]=0;
		fSignal[i]=(Double_t) ( long (fSignal[i]/oneBitResolution+64)-64 ) * oneBitResolution;
	}    
	Double_t sumCharge(0.);
	Double_t sumChargeErr(0.);
	for (Int_t i=0;i<fWaveformLength;i++)
	{
		sumCharge += fSignal[i];
		sumChargeErr += sqrt(fSignal[i]);
	}
	Double_t Coff = EnergyError*sumCharge/sumChargeErr;
	for(Int_t i=0;i<fWaveformLength;++i){
		fSignalError[i] = fSignal[i]*Coff;
	}
  //calculate baseline
	static Int_t baseline_length =100;
	for(Int_t j=0;j<baseline_length;++j){
		fBaselineValue += gRandom->Gaus(0,1)*noise_width;
	}
	fBaselineValue/=baseline_length;
}

	void 
PndEmcWaveform::AddShapedElecNoiseAndDigitise(Double_t noise_width,
		Double_t oneBitResolution, PndEmcAbsPulseshape *pulseshape, Double_t firstADCBinTime, Double_t sampleRate, Double_t EnergyError )
{
	// Do both e_noise and digitisation.
	Double_t t;
	for (  int i=0;i<fWaveformLength;i++)
	{
		t = i/sampleRate;
		Double_t ran_noise=gRandom->Gaus(0,1)*noise_width;
		fSignal[i] += pulseshape->value(t,ran_noise,firstADCBinTime);
	}
	Double_t sumCharge(0.);
	Double_t sumChargeErr(0.);
	for (Int_t i=0;i<fWaveformLength;i++)
	{
		sumCharge += fSignal[i];
		sumChargeErr += sqrt(fSignal[i]);
	}
	Double_t Coff = EnergyError*sumCharge/sumChargeErr;
	for(Int_t i=0;i<fWaveformLength;++i){
		fSignalError[i] = fSignal[i]*Coff;
	}
	//calculate baseline
	static Int_t baseline_length =100;
	for(Int_t j=0;j<baseline_length;++j){
		t =  j/sampleRate;
		Double_t ran_noise=gRandom->Gaus(0,1)*noise_width;
		fBaselineValue += pulseshape->value(t,ran_noise,firstADCBinTime);
	}
	fBaselineValue/=baseline_length;

}

	Double_t 
PndEmcWaveform::Max()
{
	Double_t _max;
	_max=*max_element(fSignal.begin(),fSignal.end());
	return _max;
}


void 
PndEmcWaveform::clearAndReset(){
	// reset values of fSignal to 0
	fill(fSignal.begin(),fSignal.end(),0);
}

PndEmcTwoCoordIndex* PndEmcWaveform::GetTCI() const
{
	PndEmcMapper *emcMap=PndEmcMapper::Instance();
	PndEmcTwoCoordIndex* tci=emcMap->GetTCI(fDetectorId);
	return tci;
}
bool PndEmcWaveform::operator == (const PndEmcWaveform& otherWave) const
{
	if(fDetectorId != otherWave.fDetectorId) return false;
	return true;
}
bool PndEmcWaveform::operator < (const PndEmcWaveform& otherWave) const
{
	Double_t overlapped_time = 0.;
	if(GetModule() == 3) overlapped_time = PndEmcWaveform::ForwardOverlapTime;
	else if(GetModule() == 5) overlapped_time = PndEmcWaveform::ShashylikOverlapTime;
	else overlapped_time = PndEmcWaveform::BarrelOverlapTime;

	if(GetDetectorId() < otherWave.GetDetectorId())
		return true;
	else if	(GetDetectorId() == otherWave.GetDetectorId()){
		if(GetTimeStamp() < otherWave.GetTimeStamp())
			return otherWave.GetTimeStamp() >  GetActiveTime() - overlapped_time;
		else
			return GetTimeStamp() >  otherWave.GetActiveTime() - overlapped_time;
	}else
		return false;

}
bool PndEmcWaveform::operator != (const PndEmcWaveform& otherWave) const
{
	if(fDetectorId != otherWave.fDetectorId) return true;
	return false;
}
bool PndEmcWaveform::equal(FairTimeStamp* data)
{
	PndEmcWaveform* other = (PndEmcWaveform*)(data);
	if(GetDetectorId() == other->GetDetectorId()) return true;
	return false;
}

PndEmcWaveform& PndEmcWaveform::operator += (const PndEmcWaveform& otherWave)
{
	if(GetTimeStamp()> otherWave.GetTimeStamp())// current wave earlier
	{
		std::cerr<<"Please make sure the eariler waveform += the later waveform"<<std::endl;
		return *this;
	}
	//++ fPileupCount ;

	const std::vector<Int_t>& evtList = otherWave.GetEvtList();
	for(Int_t i=0; i< evtList.size();++i){
		AddEvt(evtList[i]);
	}

	Int_t k =0;
	Int_t IDX = 0;
	for(; (IDX < fWaveformLength) && (k < otherWave.fWaveformLength); ++ IDX){
		if((GetTimeStamp() + IDX/fSampleRate*1.0e9) < otherWave.GetTimeStamp()){
			continue;
		}
		fSignal[IDX] += otherWave.fSignal[k];
		fSignalError[IDX] = sqrt(fSignalError[IDX]*fSignalError[IDX] + otherWave.fSignalError[k]*otherWave.fSignalError[k]);
		++k;
	}
	if(k < otherWave.fWaveformLength){
		fWaveformLength += otherWave.fWaveformLength - k ;
		for(; IDX < fWaveformLength; ++IDX, ++k){
			fSignal.push_back(otherWave.fSignal[k]);
			fSignalError.push_back(otherWave.fSignalError[k]);
		}
	}

	return *this;
}

TGraphErrors* PndEmcWaveform::ToTGraph() const 
{

	//free this object outside 
	TGraphErrors* g = new TGraphErrors(fSignal.size());
	for(Int_t i = 0; i< fSignal.size(); ++i){
		g->SetPoint(i, GetTimeStamp()/1.e9 + Double_t(i)/fSampleRate, fSignal[i]);
		g->SetPointError(i, 0, fSignalError[i]);
	}
	return g;
}

void PndEmcWaveform::SetWaveform(std::vector<Double_t>&signal,Int_t length)
{
	fSignal = signal; 
	fWaveformLength=length;
}
Double_t PndEmcWaveform::Integral() const {
	Double_t sum(0.);
	for(Int_t i=0;i<fSignal.size();++i)
		sum += fSignal[i];
	return sum;
}

ClassImp(PndEmcWaveform)