#include "PndEmcFWEndcapDigi.h"

#include "PndEmcWaveformData.h"
#include "PndEmcMultiWaveform.h"
#include "PndEmcPSAFPGA/PndEmcPSAFPGAPileupAnalyser.h"
#include "PndEmcPSAFPGA/PndEmcHighLowPSA.h"
#include "PndEmcDigi.h"
#include "PndEmcFWEndcapDigiPar.h"		
#include "PndEmcGeoPar.h"		
#include "PndEmcRecoPar.h"				
#include "PndEmcMapper.h"
#include "PndEmcStructure.h"
#include "PndEmcSimCrystalCalibrator.h"

#include "FairRootManager.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"

#include "TClonesArray.h"
#include "TStopwatch.h"
#include "TF1.h"

#include <iostream>
#include <vector>
#include <utility>

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

PndEmcFWEndcapDigi::PndEmcFWEndcapDigi(Int_t verbose, Bool_t storedigis):
  fWaveformArray(NULL), fDigiArray(NULL), fEnergyDigiThreshold(0), fDigiPosMethod(""), fEmcDigiRescaleFactor(0.), fEmcDigiPositionDepthPWO(0), fEmcDigiPositionDepthShashlyk(0), fHighgainPSA(NULL), fLowgainPSA(NULL), fHighLowPSA(verbose), fCalibrator(NULL), fDigiPar(NULL), fRecoPar(NULL), fGeoPar(NULL), fVerbose(verbose), fStoreDigis(storedigis), fTimeOrderedDigi(kFALSE), fNrOfEvents(0), evtCounter(0), nFwDigiProg(0)
{
}

//--------------
// Destructor --
//--------------
PndEmcFWEndcapDigi::~PndEmcFWEndcapDigi()
{
	if (fDigiArray!= 0) delete fDigiArray;
}

/**
 * @brief Init Task
 * 
 * Prepares the TClonesArray of PndEmcMultiWaveform for reading and PndEmcDigi for writing.
 * Also reads the EMC parameters and prepares the pulseshapes
 * (PndEmcAbsPulseshape) and pulse shape analyser (PndEmcAbsPSA) as well as the
 * calibrator (PndEmcSimCrystalCalibrator).
 * 
 * @return InitStatus
 * @retval kSUCCESS success
 */
InitStatus PndEmcFWEndcapDigi::Init()
{
	// Get RootManager
	FairRootManager* ioman = FairRootManager::Instance();
	if ( ! ioman )
	{
		cout << "-E- PndEmcFWEndcapDigi::Init: "
		<< "RootManager not instantiated!" << endl;
		return kFATAL;
	}

	// Get nr of events	(divide by 100 for progress counter)
	fNrOfEvents = ((ioman->GetInTree())->GetEntriesFast())/100;
	
	// Get input array
	fWaveformArray = (TClonesArray*) ioman->GetObject("EmcMultiWaveform");
	if (!fWaveformArray) {
		//check if EmcWaveform contains MultiWaveforms
		fWaveformArray = (TClonesArray*) ioman->GetObject("EmcWaveform");
		if((!fWaveformArray) || (!fWaveformArray->GetClass()->InheritsFrom("PndEmcMultiWaveform"))){
		cout << "-W- PndEmcFWEndcapDigi::Init: "
		<< "No PndEmcWaveform array containing multi waveforms!" << endl;
		return kERROR;
		}
if (fWaveformArray) cout << "[INFO\t] Waveform array present!" << endl;

	}
	
	if(!fDigiPar) {
		cout << "-E- PndEmcFWEndcapTimebasedWaveforms::Init: "
			<< "no DigiPar containter found" << endl;
		return kFATAL;
	}
	
	// Create and register output array
	fDigiArray = ioman->Register("EmcDigi","PndEmcDigi", "Emc", fStoreDigis);

	//position methods
	fEmcDigiPositionDepthPWO=fRecoPar->GetEmcDigiPositionDepthPWO();
	fEmcDigiPositionDepthShashlyk=fRecoPar->GetEmcDigiPositionDepthShashlyk();
	
	cout<<"fEmcDigiPositionDepthPWO: "<<fEmcDigiPositionDepthPWO<<endl;
	cout<<"fEmcDigiPositionDepthShashlyk: "<<fEmcDigiPositionDepthShashlyk<<endl;

	//digi pos method				
	fDigiPosMethod = "depth";	//ATTN: hardcoded...do we need this at all?
	fEmcDigiRescaleFactor = 1.08;	
	if (!fDigiPosMethod.CompareTo("surface")) {
		PndEmcDigi::selectDigiPositionMethod( PndEmcDigi::surface, 0., 0., 1.0 );
	} else if (!fDigiPosMethod.CompareTo("depth")) {
	   PndEmcDigi::selectDigiPositionMethod(PndEmcDigi::depth, fEmcDigiPositionDepthPWO, fEmcDigiPositionDepthShashlyk, fEmcDigiRescaleFactor);
	} else {
		cout << "-W- PndEmcFWEndcapDigi::Init: " << "Unknown digi position method!" << endl;
		return kERROR;
	}

	// psa
	if(!fDigiPar->GetPsaTypeLow().IsNull()) {

		if(fLowgainPSA) {
			std::cout << "-W in PndEmcFWEndcapDigi::Init: Lowgain PSA already set. Skipping default initialization" << std::endl;
		} else if(!fDigiPar->GetPsaTypeLow().CompareTo("PSAFPGAPileupAnalyser")) {
			PndEmcPSAFPGAPileupAnalyser* psa = new PndEmcPSAFPGAPileupAnalyser();
			psa->SetVerbose(fVerbose);

			TObjArray* rparas = fDigiPar->GetRValueParLow().Tokenize(";");
			TF1 R_thres("R_thres_Low", rparas->GetEntriesFast()>0 ? rparas->UncheckedAt(0)->GetName() : "", 1, TMath::Power(2, fDigiPar->GetNBits()));				
			TF1 R_mean("R_mean_Low", rparas->GetEntriesFast()>1 ? rparas->UncheckedAt(1)->GetName() : "", 1, TMath::Power(2, fDigiPar->GetNBits()));				

			delete rparas;

			psa->Init(std::vector<Double_t>(fDigiPar->GetPsaParLow().GetArray(), fDigiPar->GetPsaParLow().GetArray()+fDigiPar->GetPsaParLow().GetSize()), R_thres.IsZombie() ? NULL : &R_thres, R_mean.IsZombie() ? NULL : & R_mean);

			fLowgainPSA = psa;

		} else {
			std::cerr << "-E in PndEmcFWEndcapDigi::Init could not find PSA of type: " << fDigiPar->GetPsaTypeLow() << std::endl;
		}
	}

	if(!fDigiPar->GetPsaTypeHigh().IsNull()) {

		if(fHighgainPSA) {
			std::cout << "-W in PndEmcFWEndcapDigi::Init: Highgain PSA already set. Skipping default initialization" << std::endl;
		} else if(!fDigiPar->GetPsaTypeHigh().CompareTo("PSAFPGAPileupAnalyser")) {
			PndEmcPSAFPGAPileupAnalyser* psa = new PndEmcPSAFPGAPileupAnalyser();
			psa->SetVerbose(fVerbose);

			TObjArray* rparas = fDigiPar->GetRValueParHigh().Tokenize(";");
			TF1 R_thres("R_thres_High", rparas->GetEntriesFast()>0 ? rparas->UncheckedAt(0)->GetName() : "", 1, TMath::Power(2, fDigiPar->GetNBits()));				
			TF1 R_mean("R_mean_High", rparas->GetEntriesFast()>1 ? rparas->UncheckedAt(1)->GetName() : "", 1, TMath::Power(2, fDigiPar->GetNBits()));				

			delete rparas;

			psa->Init(std::vector<Double_t>(fDigiPar->GetPsaParHigh().GetArray(), fDigiPar->GetPsaParHigh().GetArray()+fDigiPar->GetPsaParHigh().GetSize()), R_thres.IsZombie() ? NULL : &R_thres, R_mean.IsZombie() ? NULL : & R_mean);

			fHighgainPSA = psa;
		} else {
			std::cerr << "-E in PndEmcFWEndcapDigi::Init could not find PSA of type: " << fDigiPar->GetPsaTypeHigh() << std::endl;
		}
	}


	if(fHighgainPSA && fLowgainPSA) {
		fHighLowPSA.Init(fHighgainPSA, fLowgainPSA, fDigiPar->GetSignalOverflowHigh(), 0, 1);
	} else {
		std::cerr << "-E- PndEmcFWEndcapDigi::Init " << "No highgain and/or lowgain psa" << std::endl;
		return kFATAL;
	}

	//calibration
	if(fCalibrator != NULL) {
			std::cout << "-W in PndEmcFWEndcapDigi::Init: Calibrator already set. Skipping default initialization"  << std::endl;
	} else {
		fCalibrator = new PndEmcSimCrystalCalibrator();
		for(Int_t imod=1; imod<=5; ++imod) { //TODO: different calibration constants for different modules
			fCalibrator->SetCalibration(imod, fDigiPar->GetCalibHigh(), 0, fDigiPar->GetSignalOverflowHigh()); 	//idx:0 <-> highgain
			fCalibrator->SetCalibration(imod, fDigiPar->GetCalibLow(), 1, 0); 					//idx:1 <-> lowgain
		}
		fCalibrator->Init();
	}

	fEnergyDigiThreshold=fDigiPar->GetEnergyDigiThreshold();

	fGeoPar->InitEmcMapper();
	PndEmcStructure::Instance();
	
	cout << "-I- PndEmcFWEndcapDigi: Intialization successful" << endl;
	return kSUCCESS;
}

/**
 * @brief Runs the task.
 * 
 * The task loops over the waveforms and uses the pulse shape analyser (PndEmcAbsPSA) to
 * extract signal height and timing. The calibrator (PndEmcAbsCrystalCalibrator) is then
 * used to calculate the energy. If the energy is above the threshold (@ref fEnergyDigiThreshold),
 * a PndEmcDigi is created.
 * 
 * @param opt unused
 * @return void
 */
void PndEmcFWEndcapDigi::Exec(Option_t*)
{
	TStopwatch timer;
	if (fVerbose>2){
		timer.Start();
	}

	if (evtCounter-1==nFwDigiProg*fNrOfEvents && nFwDigiProg<100) {
		nFwDigiProg += 1;
		cout << "\r[INFO\t] PndEmcFWEndcapDigi: " << nFwDigiProg << "\% completed." << std::flush;	
	}
	evtCounter++;

	fDigiArray->Delete();

	Double_t MinDiff=0;
	Int_t MinDiffIdx=0;
	std::vector<Double_t> diffVector;
	FairRootManager* ioman = FairRootManager::Instance();

	Double_t energy;
	Double_t digi_time;
	Int_t hitIndex;
	Int_t nHits;
	Int_t detId;
	Int_t trackId;
	Int_t nWaveforms = fWaveformArray->GetEntriesFast();
	PndEmcWaveform* theWaveform;
	if (fVerbose>3) cout << "PndEmcFWEndcapDigi, #input waveforms: " << nWaveforms << " in Event: " << FairRootManager::Instance()->GetEntryNr() << endl;
	for (Int_t iWaveform=0; iWaveform<nWaveforms; iWaveform++) {

		theWaveform = (PndEmcWaveform*) fWaveformArray->At(iWaveform);
		hitIndex=theWaveform->GetHitIndex();
		detId=theWaveform->GetDetectorId();
		trackId=theWaveform->GetTrackId();

		//Double_t timeshift; // how maximum is shifted (not used?)

		nHits = fHighLowPSA.Process(theWaveform);	

		for(Int_t iHit=0; iHit<nHits; ++iHit) {
			fHighLowPSA.GetHit(iHit, energy, digi_time);	
			//PndEmcAbsCrystalCalibrator::CalibrationStatus_t CalS //[R.K. 03/2017] unused variable?tat;

			//cout << "#" << iHit << ".\tDetId: " << detId << "\tsample time: " << digi_time << "\traw energy: " << energy << "\t";
			//CalStat =  //[R.K. 03/2017] unused variable?
      fCalibrator->Calibrate(energy, detId, fHighLowPSA.GetWaveformIdx(iHit));

			Double_t sampleRate = theWaveform->GetSampleRate(); // 8e7 Hz

			digi_time/=sampleRate;//Hz
			digi_time*=1e9;		// back to ns

			if (energy>fEnergyDigiThreshold) {
				Double_t timestamp=theWaveform->GetTimeStamp() + digi_time;
				//Double_t timestamp=theWaveform->GetTimeStamp();
//cout.precision(20);
//cout << "WF time:"<< theWaveform->GetTimeStamp() << " digitime:"<<digi_time<<", energy:"<<energy<<endl;
				//std::cout << "PndEmcFWEndcapDigi::Exec Waveform TS: " << theWaveform->GetTimeStamp() << " digiTime: " << digi_time << std::endl;

				PndEmcDigi* myDigi = new((*fDigiArray)[fDigiArray->GetEntriesFast()]) PndEmcDigi(trackId, detId, energy, timestamp, hitIndex);
//cout << "// --- Added digi with energy: " << myDigi->GetEnergy() << " x: " << myDigi->where().x() << " y: " << myDigi->where().y() << " t: " << myDigi->GetTimeStamp() << endl;
				myDigi->ResetLinks();
				myDigi->AddLinks(theWaveform->GetLinksWithType(ioman->GetBranchId("EmcHit")));


		// select hit with energy closest to extracted digi energy to use in new timestamp
				/*FairMultiLinkedData digiLinks = myDigi->GetLinksWithType(ioman->GetBranchId("EmcHit"));
				diffVector.clear();
				//cout << "links:" << digiLinks.GetNLinks() << std::flush;
				for (int iLink=0; iLink<digiLinks.GetNLinks(); iLink++){
					PndEmcHit* theHit = (PndEmcHit*)ioman->GetCloneOfLinkData(digiLinks.GetLink(iLink));
					if(theHit) {
						diffVector.push_back( TMath::Abs( theHit->GetEnergy()-energy ) );	
						//diffVector.push_back( TMath::Abs( ioman->GetEventTime() + theHit->GetTime()*1.0e9-timestamp ) );	
//cout << "Hittime:" << ioman->GetEventTime() + theHit->GetTime()*1e9	<< ", timestamp:" << timestamp << ", energy:" << theHit->GetEnergy() << endl;
					} else {
						cout << "-E in PndEmcViewClusters::Exec FairLink: " << digiLinks.GetLink(iLink) << " to EmcHit delivers null" << std::endl;
					}
				}
				if (diffVector.size()>0) {
					MinDiff = diffVector[0];
					MinDiffIdx = 0;
					for (int e=0; e<diffVector.size(); e++) { // get smallest difference
						if (diffVector[e] < MinDiff) {
							MinDiff = diffVector[e];	
							MinDiffIdx = e;
						}
					}
//if(diffVector.size()==1)cout << MinDiff << " (" << diffVector.size() << ")" << endl;
					PndEmcHit* bestHit = (PndEmcHit*)ioman->GetCloneOfLinkData(digiLinks.GetLink(MinDiffIdx));
//cout << ioman->GetEventTime() + bestHit->GetTime()*1.0e9 << "-" << timestamp << "=\t" << (ioman->GetEventTime() + bestHit->GetTime()*1.0e9-timestamp) << endl;
					myDigi->SetTimeStamp( rand.Gaus( ioman->GetEventTime() + bestHit->GetTime()*1.0e9, GetEnergyDepTimeResolution(energy) ) ); // assign timestamp to digi using time of hit, smeared with energy-dependent Gaussian (only if the digi was not a pure noise hit, i.e. it came from an EmcHit. In that case, leave the timestamp as defined by the waveform)
				}*/
		// finished selecting hit

			}
				
		}

	}

	if (fVerbose>3) {
		cout << "PndEmcFWEndcapDigi, #output digis: " << fDigiArray->GetEntriesFast() << endl;
		timer.Stop();
		Double_t rtime = timer.RealTime();
		Double_t ctime = timer.CpuTime();
		cout << "PndEmcFWEndcapDigi, Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
	}
}

Double_t PndEmcFWEndcapDigi::GetEnergyDepTimeResolution(Double_t E)
{	
	// calculate energy-dependent time resolution. Formula taken from http://dx.doi.org/10.1016/j.nima.2011.06.044
	return 0.55+5.5*TMath::Exp(-(TMath::Log(2)/0.0250)*E);
}

void PndEmcFWEndcapDigi::SetParContainers() 
{
  // Get run and runtime database
  FairRun* run = FairRun::Instance();
  if ( ! run ) Fatal("SetParContainers", "No analysis run");

  FairRuntimeDb* db = run->GetRuntimeDb();
  if ( ! db ) Fatal("SetParContainers", "No runtime database");

  // Get Emc digitisation parameter container
  fDigiPar = dynamic_cast<PndEmcFWEndcapDigiPar*>(db->getContainer("PndEmcFWEndcapDigiPar"));
  
  // Get Emc reconstruction parameter container
  fRecoPar = (PndEmcRecoPar*) db->getContainer("PndEmcRecoPar");

  // Get Emc geometry parameter container
  fGeoPar = (PndEmcGeoPar*) db->getContainer("PndEmcGeoPar");
}

void PndEmcFWEndcapDigi::SetStorageOfData(Bool_t val)
{
  fStoreDigis = val;
  return;
}

ClassImp(PndEmcFWEndcapDigi);