//----------------------------------------------------------------------
// File and Version Information:
//      $Id: //
// Description:
//      Class PndEmcMultiWaveformToCalibratedDigi. Module to take the ADC waveforms and produces digi.
// 
//	 Software developed for the BaBar Detector at the SLAC B-Factory.
// Adapted for the PANDA experiment at GSI		
//		
// Author List:
//      Phil Strother                  Original Author
// Dima Melnichuk - adaption for PANDA		
// Copyright Information:
//      Copyright (C) 1996             Imperial College
//
//----------------------------------------------------------------------

#include "PndEmcMultiWaveformToCalibratedDigi.h"

#include "PndEmcMultiWaveform.h"
#include "PndEmcDigi.h"
#include "PndEmcDigiPar.h"		
#include "PndEmcRecoPar.h"				
#include "PndEmcAsicPulseshape.h"
#include "PndEmcPSAParabolic.h"
#include "PndEmcPSAMatchedDigiFilter.h"
#include "PndEmcSimCrystalCalibrator.h"

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

#include "TClonesArray.h"
#include "TStopwatch.h"
		
#include <iostream>
#include <vector>
#include <utility>
//#include <map>

using std::cout;
using std::endl;
//using std::fstream;

PndEmcMultiWaveformToCalibratedDigi::PndEmcMultiWaveformToCalibratedDigi(Int_t verbose, Bool_t storedigis):
  fWaveformArray(new TClonesArray()), fDigiArray(new TClonesArray()), fSampleRate(0), fSampleRate_PMT(0), fEnergyDigiThreshold(0), fASIC_Shaping_int_time(0), fPMT_Shaping_int_time(0), fPMT_Shaping_diff_time(0), fCrystal_time_constant(0), fShashlyk_time_constant(0), fNumber_of_samples_in_waveform(0), fNumber_of_samples_in_waveform_pmt(0), fDigiPosMethod(0), fEmcDigiRescaleFactor(0), fEmcDigiPositionDepthPWO(0), fEmcDigiPositionDepthShashlyk(0), fPulseshape(0), fPulseshape_pmt(0), fpsaAlgorithm(0), fpsaAlgorithm_pmt(0), fDigiPar(new PndEmcDigiPar()), fRecoPar(new PndEmcRecoPar()), fVerbose(verbose), fStoreDigis(storedigis), fTimeOrderedDigi(kFALSE), fWaveform_Signal_MaxTimediff(0), fWaveform_Signal_Overflow(0)
{
  fDigiPosMethod="depth";// "surface" or "depth"
  fEmcDigiRescaleFactor=1.08;
  //fPndEmcDigiPositionDepth=6.2;
}

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

PndEmcMultiWaveformToCalibratedDigi::~PndEmcMultiWaveformToCalibratedDigi()
{
	if (fDataBuffer!= 0) delete fDataBuffer;
}


InitStatus PndEmcMultiWaveformToCalibratedDigi::Init()
{
	// Get RootManager
	FairRootManager* ioman = FairRootManager::Instance();
	if ( ! ioman )
	{
		cout << "-E- PndEmcMultiWaveformToCalibratedDigi::Init: "
		<< "RootManager not instantiated!" << endl;
		return kFATAL;
	}
	
	// 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- PndEmcMultiWaveformToCalibratedDigi::Init: "
		<< "No PndEmcWaveform array!" << endl;
		return kERROR;
		}

	}
	
	// Create and register output array
//	fDigiArray = new TClonesArray("PndEmcDigi");
//	ioman->Register("EmcDigi","Emc",fDigiArray,fStoreDigis);
	
	fDataBuffer = new PndEmcDigiWriteoutBuffer("EmcDigi", "Emc", fStoreDigis);
	fDataBuffer = (PndEmcDigiWriteoutBuffer*)ioman->RegisterWriteoutBuffer("EmcDigi", fDataBuffer);
//   fDataBuffer->ActivateBuffering(fTimeOrderedDigi);
   fDataBuffer->ActivateBuffering(kFALSE);


	fSampleRate=fDigiPar->GetSampleRate();
	fSampleRate_PMT=fDigiPar->GetSampleRate_PMT();
	fASIC_Shaping_int_time=fDigiPar->GetASIC_Shaping_int_time();      //s
	fPMT_Shaping_int_time=fDigiPar->GetPMT_Shaping_int_time();      //s
	fPMT_Shaping_diff_time=fDigiPar->GetPMT_Shaping_diff_time();      //s
	fCrystal_time_constant=fDigiPar->GetCrystal_time_constant();  //s
	fShashlyk_time_constant=fDigiPar->GetShashlyk_time_constant();  //s
	fNumber_of_samples_in_waveform=fDigiPar->GetNumber_of_samples_in_waveform();
	fNumber_of_samples_in_waveform_pmt=fDigiPar->GetNumber_of_samples_in_waveform_pmt();
	fEnergyDigiThreshold=fDigiPar->GetEnergyDigiThreshold();
	fWaveform_Signal_Overflow = fDigiPar->GetWaveform_Signal_Overflow();
	fWaveform_Signal_MaxTimediff = fDigiPar->GetWaveform_Signal_MaxTimediff();
	fEmcDigiPositionDepthPWO=fRecoPar->GetEmcDigiPositionDepthPWO();
	fEmcDigiPositionDepthShashlyk=fRecoPar->GetEmcDigiPositionDepthShashlyk();
	
	cout<<"fEmcDigiPositionDepthPWO: "<<fEmcDigiPositionDepthPWO<<endl;
	cout<<"fEmcDigiPositionDepthShashlyk: "<<fEmcDigiPositionDepthShashlyk<<endl;
				
	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- PndEmcMultiWaveformToCalibratedDigi::Init: "
		<< "Unknown digi position method!" << endl;
		return kERROR;
	}
	
	fPulseshape= new PndEmcAsicPulseshape(fASIC_Shaping_int_time,fCrystal_time_constant);
	fPulseshape_pmt= new PndEmcCRRCPulseshape(fPMT_Shaping_int_time,fPMT_Shaping_diff_time,fShashlyk_time_constant);

	if(fpsaAlgorithm == NULL){
		// Matched digital filter
		// Parameters of the filter are hardcoded at the moment
		// For different pulseshape in barrel and endcaps different filters should be implemented
		std::vector<Double_t> params;
		params.push_back(30); // width
		params.push_back(fSampleRate); // Sample rate
		fpsaAlgorithm = new PndEmcPSAMatchedDigiFilter(params,fPulseshape);
	}

	if(fpsaAlgorithm_pmt == NULL){
		// Pulse shape analysis algorithm.
		// Simple parabolic fit.
		fpsaAlgorithm_pmt = new PndEmcPSAParabolic();
	}

	if(fCalibrator == NULL){
		PndEmcSimCrystalCalibrator *SimCalibrator = new PndEmcSimCrystalCalibrator();
		// Determine normalisation constant for PndEmcWaveform
		PndEmcWaveform *tmpwaveform=new PndEmcWaveform(0,101010001, fNumber_of_samples_in_waveform);
		PndEmcWaveform *tmpwaveform2=new PndEmcWaveform(0,101010001, fNumber_of_samples_in_waveform_pmt);
		PndEmcHit *gevHit=new PndEmcHit();
		gevHit->SetEnergy(1.0);
		gevHit->SetTime(0.);
		tmpwaveform->UpdateWaveform(gevHit, 0, false, 1., 0., fSampleRate, fPulseshape);
		tmpwaveform2->UpdateWaveform(gevHit, 0, false, 1., 0., fSampleRate_PMT, fPulseshape_pmt);
		Double_t tmpPeakPosition;
		Double_t WfNormalisation;
		fpsaAlgorithm->Process(tmpwaveform,WfNormalisation,tmpPeakPosition);
		for(Int_t i = 1; i <5;i++){
			SimCalibrator->SetCalibration(i,WfNormalisation);
		}
		fpsaAlgorithm_pmt->Process(tmpwaveform2,WfNormalisation,tmpPeakPosition);
		SimCalibrator->SetCalibration(5,WfNormalisation);
		delete tmpwaveform;
		delete tmpwaveform2;
		fCalibrator=SimCalibrator;
	}
	fCalibrator->Init();
	cout << "-I- PndEmcMultiWaveformToCalibratedDigi: Intialization successfull" << endl;
	
	return kSUCCESS;
}

void PndEmcMultiWaveformToCalibratedDigi::Exec(Option_t* opt)
{
	TStopwatch timer;
	if (fVerbose>2){
		timer.Start();
	}
	
	Double_t EventTime = FairRootManager::Instance()->GetEventTime();

	Double_t peakPosition;
	Double_t energy;
	Double_t digi_time;
	Int_t hitIndex;
	std::vector< std::vector< std::pair<Double_t,Double_t > > > hits;
	Double_t DigiInfo[2];
	Int_t nHits;
	Int_t HitsInEvent;
	Int_t detId;
	Int_t trackId;
	Int_t module;
	Int_t nWaveforms = fWaveformArray->GetEntriesFast();
	PndEmcAbsPSA *thePSA;
	Double_t theEnergyNorm;
	Double_t theSampleRate;
	Int_t nSignal;
	PndEmcMultiWaveform* theWaveform;
//   cout << fDataBuffer->GetAllData().size() << " Entries in Buffer" << endl;
//   cout<<"PndEmcMultiWaveformToCalibratedDigi: "<<nWaveforms<<" waveforms to convert"<<endl;
	HitsInEvent=0;
	for (Int_t iWaveform=0; iWaveform<nWaveforms; iWaveform++) {
		theWaveform = (PndEmcMultiWaveform*) fWaveformArray->At(iWaveform);
		hitIndex=theWaveform->GetHitIndex();
		detId=theWaveform->GetDetectorId();
		trackId=theWaveform->GetTrackId();
		module=theWaveform->GetModule();
		nSignal=theWaveform->GetNumberOfWaveforms();
		Double_t timeshift; // how maximum is shifted
		if(module == 5){
			thePSA = fpsaAlgorithm_pmt;
			theSampleRate = fSampleRate_PMT;
		} else {
			thePSA = fpsaAlgorithm;
			theSampleRate = fSampleRate;
		}
		for(Int_t i = 0; i< hits.size(); i++){
			hits[i].clear();
		}
		hits.resize(nSignal);
		for(Int_t iSignal = 0; iSignal < nSignal; iSignal++){
			theWaveform->SetActiveWaveform(iSignal);
			nHits = thePSA->Process(theWaveform);
//         std::cout << "found " << nHits << "Hits in Signal " << iSignal << std::endl;
			for(Int_t jHit = 0; jHit<nHits; jHit++){
				thePSA->GetHit(jHit,energy,digi_time);
				hits[iSignal].push_back(std::make_pair(energy,digi_time));
			}
		}
//      std::cout << "hits has " << hits.size() << " Signals" << std::endl;
//      std::cout << "Signal [0] has " <<hits[0].size() << " hits" << std::endl;

	
		for(Int_t iHit = 0 ; iHit< hits[0].size(); iHit++){
			energy=hits[0][iHit].first;
			digi_time=hits[0][iHit].second;
			Bool_t badCal=kFALSE;
/*
			if(hits.size()>0 && energy > fWaveform_Signal_Overflow){
				badCal=kTRUE;
				for(Int_t iHit1 = 0; iHit1 < hits[1].size();iHit1++){
					if(TMath::Abs(digi_time - hits[1][iHit1].second)< fWaveform_Signal_MaxTimediff){
						energy=hits[1][iHit1].first;
						if(fCalibrator->Calibrate(energy,detId,1)!=PndEmcAbsCrystalCalibrator::kCALOK){
							badCal=kTRUE;
						}else{
							badCal=kFALSE;
							break;
						}
					}
				}
			}else{
				if(fCalibrator->Calibrate(energy,detId,0)!=PndEmcAbsCrystalCalibrator::kCALOK){
					badCal=kTRUE;
				}
			}
			if(badCal){
				continue;
			}
*/
			PndEmcAbsCrystalCalibrator::CalibrationStatus_t CalStat;
//         cout << "DetId: " << detId << " raw: " << energy;
			CalStat = fCalibrator->Calibrate(energy,detId,0);
//         cout << " Energy: " << energy;
			Bool_t goodcal=kFALSE;
			switch(CalStat){
				case PndEmcAbsCrystalCalibrator::kCALOK:
					goodcal=kTRUE;
//               cout << " good";
					break;
				case PndEmcAbsCrystalCalibrator::kCALOVERFLOW:
//               cout << " overflow time: " << digi_time;
					for(Int_t iHit1 = 0; iHit1 < hits[1].size();iHit1++){
//                  cout << " lowgain time: " << hits[1][iHit1].second;
						if(TMath::Abs(digi_time - hits[1][iHit1].second)< fWaveform_Signal_MaxTimediff){
							energy=hits[1][iHit1].first;
//                        cout << " Lowgain: " << energy;
							if(fCalibrator->Calibrate(energy,detId,1)!=PndEmcAbsCrystalCalibrator::kCALOK){
//                        cout << " bad";
								continue;
							}else{
//                        cout << " LowEnergy: " << energy;
								goodcal=kTRUE;
								break;
							}
						}
					}
					break;
				default:
					break;
			}
			if(!goodcal){
//            cout << " discarded" << endl;
				continue;
			}
//         cout << endl;


			digi_time/=fSampleRate;
			digi_time*=1e9;//ns
//         std::cout << "digi energy: " << energy << " digi time: " << digi_time << std::endl;
			if (energy>fEnergyDigiThreshold)
			{
				Double_t timestamp=EventTime+digi_time;	
				PndEmcDigi* myDigi = new PndEmcDigi(trackId,detId, energy, timestamp, hitIndex);
				myDigi->AddLink(FairLink("EmcWaveform", iWaveform));
				fDataBuffer->FillNewData(myDigi, 300, 300); // 300 ns
				HitsInEvent++;
				if (fVerbose>2){
					cout<<"timestamp="<<timestamp<<" EventTime="<<EventTime<<" digi_time="<<digi_time<<endl;
				}
			}
		}
	}
//   cout << fDataBuffer->GetAllData().size() << " Entries in Buffer" << endl;
//   cout << HitsInEvent << " Hits in this Event" << endl;
	if (fVerbose>2){
		timer.Stop();
		Double_t rtime = timer.RealTime();
		Double_t ctime = timer.CpuTime();
		cout << "PndEmcMultiWaveformToCalibratedDigi, Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
	}

}

void PndEmcMultiWaveformToCalibratedDigi::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 = (PndEmcDigiPar*) db->getContainer("PndEmcDigiPar");
  
  // Get Emc reconstruction parameter container
  fRecoPar = (PndEmcRecoPar*) db->getContainer("PndEmcRecoPar");
 
}

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

ClassImp(PndEmcMultiWaveformToCalibratedDigi);