//=====================================================================
//
//  $Id: $
//
//	EmcWaveform.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 "EmcWaveform.h"
#include "CbmEmcHit.h"
#include "EmcCRRCPulseshape.h"
#include "EmcCR2RCPulseshape.h"
#include "EmcMapper.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 --
//----------------

EmcWaveform::EmcWaveform(){}



EmcWaveform::EmcWaveform(int trackId, long detId, double tau1, double tau2, double sampleRate,double firstADCTime, 
			 double pePerMeV, double tauCrystal, long waveform_length, double excessNoiseFactor,bool use_photon_statistic):
	fTrackId(trackId),
	fDetId(detId),
	fTauCrystal(tauCrystal),
	fTau1(tau1),
	fTau2(tau2),
	fSampleRate(sampleRate),
	fFirstADCBinTime ( firstADCTime),
	fPePerMeV(pePerMeV),
	fExcessNoiseFactor(excessNoiseFactor),
	fUsePhotonStatistic(use_photon_statistic),
	fWaveformLength(waveform_length),
	fSignal(waveform_length,0.),
	fPulseshape(tau1,tau2,tauCrystal)
{
	EmcMapper *fEmcMap=EmcMapper::Instance(1);
	fTCI=fEmcMap->GetTCI(detId);
}

EmcWaveform::EmcWaveform(const EmcWaveform &copy)
{
		fTrackId=copy.fTrackId;
		fDetId = copy.fDetId;
		fTCI=copy.fTCI;
		fTau1=copy.fTau1;
		fTau2=copy.fTau2;
		fSampleRate=copy.fSampleRate;
		fTauCrystal=copy.fTauCrystal;
		fFirstADCBinTime =copy.fFirstADCBinTime;
		fPePerMeV=copy.fPePerMeV;
		fExcessNoiseFactor=copy.fExcessNoiseFactor;
		fUsePhotonStatistic=copy.fUsePhotonStatistic;
		fWaveformLength=copy.fWaveformLength;
		fHitList=copy.fHitList;
		fSignal=copy.fSignal;
		fPulseshape=copy.fPulseshape;
}

EmcWaveform&
EmcWaveform::operator=(const EmcWaveform &copy){
  if (this != &copy){
		fTrackId=copy.fTrackId;
		fDetId = copy.fDetId;
		fTCI=copy.fTCI;
		fTau1=copy.fTau1;
		fTau2=copy.fTau2;
		fSampleRate=copy.fSampleRate;
		fTauCrystal=copy.fTauCrystal;
		fFirstADCBinTime =copy.fFirstADCBinTime;
		fPePerMeV=copy.fPePerMeV;
		fExcessNoiseFactor=copy.fExcessNoiseFactor;
		fUsePhotonStatistic=copy.fUsePhotonStatistic;
		fWaveformLength=copy.fWaveformLength;
		fHitList=copy.fHitList;
		fSignal=copy.fSignal;
		fPulseshape=copy.fPulseshape;
  }
  return *this;
}

//--------------
// Destructor --
//--------------

EmcWaveform::~EmcWaveform()
{
}

    
//-------------
// Selectors --
//-------------

double 
EmcWaveform::get_scale() const
{
	return getNormalisation();
}

double 
EmcWaveform::getNormalisation() const
{
  // Function to return the equivalent pulse height for a 1 GeV pulse
  // Used  in EmcHitsToWaveform to determine electronics noise scale

	EmcWaveform newWaveform(*this);
	newWaveform.clearAndReset();
	
	CbmEmcHit *gevHit=new CbmEmcHit();
	gevHit->SetEnergy(1.0);
	gevHit->SetTime(0.);
	newWaveform.update_waveform(gevHit);
 	delete gevHit;

	double max=newWaveform.max();
	return max;
}

    
//-------------
// Modifiers --
//-------------

void
EmcWaveform::addHitToList( CbmEmcHit* hit )
{
	fHitList.push_back(hit);
}  

void 
EmcWaveform::update_waveform(CbmEmcHit *hit)
{
	fHitList.push_back(hit);
	fTrackId=hit->GetTrackId(); // Take track Id from the first hit

	double energy=(double)hit->GetEnergy();
	double time=(double)hit->GetTime();
	
	make_waveform(energy, time);

}

void 
EmcWaveform::make_waveform(double energy, double time)
{
	double amplitude;
	double photonStatFactor;
	if (fUsePhotonStatistic)
	{
		double crystalPhotonsMeV = 1.0e3 * energy * fPePerMeV;
		photonStatFactor = gRandom->Gaus(1,sqrt(fExcessNoiseFactor/crystalPhotonsMeV));
   }
	else
	{
		photonStatFactor=1.;
	}

	amplitude = energy*photonStatFactor;
	
	double t;
	double time_offset=time;
	for (  int i=0;i<fWaveformLength;i++)
	{
		t= time+fFirstADCBinTime+i/fSampleRate;
		fSignal[i]+=(fPulseshape.value(t,amplitude,time_offset));
	}
}

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


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

void 
EmcWaveform::add_elec_noise_and_digitise(double noise_width,
					 double oneBitResolution)
{
	// Do both e_noise and digitisation.  
	
	for (  int i=0;i<fWaveformLength;i++)
	{
		double ran_noise=gRandom->Gaus(0,1)*noise_width;
		fSignal[i]+=ran_noise;
		if (fSignal[i]<0) fSignal[i]=0;
		fSignal[i]=(double) ( long (fSignal[i]/oneBitResolution+64)-64 ) * oneBitResolution;
	}    
}

void 
EmcWaveform::add_shaped_elec_noise_and_digitise(double noise_width,
						double oneBitResolution)
{
  // Do both e_noise and digitisation.

	// "0" corresponds to timing constant of signal
	// for noise can be assumed 0 (it is taken tauInt*1e-5 to avoid nan)
	EmcCRRCPulseshape* fPulseshape2= new EmcCRRCPulseshape(fTau1,fTau2,fTau1*1e-5);

	double t;
	for (  int i=0;i<fWaveformLength;i++)
	{
			t= fFirstADCBinTime+i/fSampleRate;
			double ran_noise=gRandom->Gaus(0,1)*noise_width;
			fSignal[i]+=fPulseshape2->value(t,ran_noise,0);
	}    
	
	for (  int i=0;i<fWaveformLength;i++)
	{
		fSignal[i]=(double) ( long (fSignal[i]/oneBitResolution+64)-64 ) * oneBitResolution;
	}
	delete fPulseshape2;
}

void EmcWaveform::fitPeak(double& ampl, double& pos,
			  int peakBin) const
{
// parabolic fit
	ampl=pos=-1.;
	if (peakBin>0 && peakBin<fWaveformLength-1)
	{
		long theBin(peakBin);
	
		double pValue = fSignal[peakBin];
		double pPosition = double(peakBin);
	
		double leftValue = fSignal[theBin-1];
		double rightValue = fSignal[theBin+1];
		if (leftValue<pValue && rightValue<pValue) {
			double d = 0.25*(rightValue-leftValue);
			double b = pValue-0.5*(leftValue+rightValue);
			pValue += d*d/b;
			pPosition += d/b;
		}
		ampl=pValue;
		pos=pPosition;
  }
}

void EmcWaveform::fitPeak(double& ampl, double& pos, int start, int end) const
{
	std::vector<double>::const_iterator p;
	p=max_element(fSignal.begin()+start,fSignal.begin()+end);
	int pPosition = distance(fSignal.begin(),p);
	fitPeak(ampl,pos,pPosition);
}
  
void EmcWaveform::fitPeak(double& ampl, double& pos) const
{
	std::vector<double>::const_iterator p;
	p=max_element(fSignal.begin(),fSignal.end());
	int pPosition = distance(fSignal.begin(),p);
	fitPeak(ampl,pos,pPosition);
}

double 
EmcWaveform::max()
{
	double _max;
	_max=*max_element(fSignal.begin(),fSignal.end());
	return _max;
}

void 
EmcWaveform::clearAndReset(){
	// reset values of fSignal to 0
	fill(fSignal.begin(),fSignal.end(),0);
	// clear CbmEmcHit list
	fHitList.clear();
}

int EmcWaveform::GetNHits() const
{
	return fHitList.size();
}

void EmcWaveform::Streamer(TBuffer &R__b)
{
   // Stream an object of class EmcWaveform.

   if (R__b.IsReading()) {
      EmcWaveform::Class()->ReadBuffer(R__b, this);
		EmcMapper *fEmcMap=EmcMapper::Instance(1);
		fTCI=fEmcMap->GetTCI(fDetId);

   } else {
      EmcWaveform::Class()->WriteBuffer(R__b, this);
   }
}

ClassImp(EmcWaveform)