//----------------------------------------------------------------------
// File and Version Information:
//      $Id: //
// Description:
//      Class PndEmcHitsToWaveform. Module to take the hit list for the 
//      calorimeter and make ADC waveforms from them.
// 
//	 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 "PndEmcHitsToWaveform.h"

#include "PndEmcHit.h"
#include "PndEmcWaveform.h"
#include "PndEmcMapper.h"
#include "PndEmcStructure.h"
#include "PndEmcDigiPar.h"		

#include "CbmRootManager.h"
#include "CbmRunAna.h"
#include "CbmRuntimeDb.h"

#include "TStopwatch.h"
#include "TROOT.h"
#include "TClonesArray.h"

#include <iostream>
#include <map>
#include <string>

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

PndEmcHitsToWaveform::PndEmcHitsToWaveform(Int_t verbose)
{
	fVerbose = verbose;
	fEnergyCutGeV=1e-5; //GeV (to sparsify waveform list)
}

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

PndEmcHitsToWaveform::~PndEmcHitsToWaveform()
{
}


InitStatus PndEmcHitsToWaveform::Init()
{
	// Get RootManager
	CbmRootManager* ioman = CbmRootManager::Instance();
	if ( ! ioman )
	{
		cout << "-E- PndEmcHitsToWaveform::Init: "
		<< "RootManager not instantiated!" << endl;
		return kFATAL;
	}
	
	// Get input array
	fHitArray = (TClonesArray*) ioman->GetObject("EmcHit");
	if ( ! fHitArray ) {
		cout << "-W- PndEmcHitsToWaveform::Init: "
		<< "No EmcHit array!" << endl;
		return kERROR;
	}
	
	// Create and register output array
	fWaveformArray = new TClonesArray("PndEmcWaveform");
	
	ioman->Register("EmcWaveform","Emc",fWaveformArray,kTRUE);
	
	cout << "-I- PndEmcHitsToWaveform: Intialization successfull" << endl;

	// Geometry loading
// 	TFile *infile = ioman->GetInFile();
// 	TGeoManager *geoMan = (TGeoManager*) infile->Get("CBMGeom");
	TGeoManager *geoMan = (TGeoManager*) gROOT->FindObject("CBMGeom");

	fNBits=fDigiPar->GetNBits();
	fDetectedPhotonsPerMeV=fDigiPar->GetDetectedPhotonsPerMeV();
	fEnergyRange=fDigiPar->GetEnergyRange(); //GeV
	fExcessNoiseFactor=fDigiPar->GetExcessNoiseFactor();
	fFirstSamplePhase=fDigiPar->GetFirstSamplePhase();
	fNumber_of_samples_in_waveform=fDigiPar->GetNumber_of_samples_in_waveform();
	fShaping_diff_time=fDigiPar->GetShaping_diff_time();     //s
	fShaping_int_time=fDigiPar->GetShaping_int_time();      //s
	fCrystal_time_constant=fDigiPar->GetCrystal_time_constant();  //s
	fIncoherent_elec_noise_width_GeV=fDigiPar->GetIncoherent_elec_noise_width_GeV(); //GeV
	fSampleRate=fDigiPar->GetSampleRate();
	fUse_shaped_noise=fDigiPar->GetUse_shaped_noise();
	fUse_photon_statistic=fDigiPar->GetUse_photon_statistic();
	
	// Test how parameters were read from DB.
	cout<<"nBits "<<fNBits<<endl;
	cout<<"detectedPhotonsPerMeV "<<fDetectedPhotonsPerMeV<<endl;
	cout<<"energyRange "<<fEnergyRange<<endl;
	cout<<"excessNoiseFactor "<<fExcessNoiseFactor<<endl;
	cout<<"firstSamplePhase "<<fFirstSamplePhase<<endl;
	cout<<"number_of_samples_in_waveform "<<fNumber_of_samples_in_waveform<<endl;
	cout<<"Shaping_diff_time "<<fShaping_diff_time<<endl;
	cout<<"Shaping_int_time "<<fShaping_int_time<<endl;
	cout<<"crystal_time_constant "<<fCrystal_time_constant<<endl;
	cout<<"incoherent_elec_noise_width_GeV "<<fIncoherent_elec_noise_width_GeV<<endl;
	cout<<"sampleRate "<<fSampleRate<<endl;
	cout<<"use_shaped_noise "<<fUse_shaped_noise<<endl;
	cout<<"use_photon_statistic "<<fUse_photon_statistic<<endl;

	fMapVersion=fDigiPar->GetMapperVersion();
	cout<<"fMapVersion: "<<fMapVersion<<endl;
	PndEmcMapper::Instance(fMapVersion);
	PndEmcStructure::Instance(geoMan);

	// Calculate 1 bit resolution (in units of FADC amplitude)
	PndEmcWaveform *tmpwaveform;
	tmpwaveform=new PndEmcWaveform(0,101010001, fShaping_diff_time, fShaping_int_time, fSampleRate, 0.0,0.0,fCrystal_time_constant,fNumber_of_samples_in_waveform,fExcessNoiseFactor,false);

	fGevPeakAnalogue = tmpwaveform->get_scale();
	
	fOneBitResolution=fEnergyRange/((double) (1<<fNBits))*fGevPeakAnalogue;
	
	fFirstADCBinTime=fFirstSamplePhase/fSampleRate;
	
	tciMap=PndEmcMapper::Instance(0)->GetTciMap();

	delete tmpwaveform;
	
	return kSUCCESS;
}

void PndEmcHitsToWaveform::Exec(Option_t* opt)
{
	TStopwatch timer;
	if (fVerbose>0){
		timer.Start();
	}
	// Reset output array
	if ( ! fWaveformArray ) Fatal("Exec", "No Waveform Array");
	fWaveformArray->Delete();
	// Variable declaration
	PndEmcHit* theHit = NULL;

	for(tciIter=tciMap.begin();tciIter!=tciMap.end();++tciIter)
	{	
		 PndEmcWaveform* tempWaveform = new PndEmcWaveform(0,tciIter->first,
				       fShaping_diff_time,fShaping_int_time,
				       fSampleRate,
				       fFirstADCBinTime,
				       fDetectedPhotonsPerMeV,
				       fCrystal_time_constant,
				       fNumber_of_samples_in_waveform,
						 fExcessNoiseFactor);

		wf_map[tciIter->first]=(tempWaveform);
	}


	// Loop over PndEmcHits to add them to correspondent wavefoorms
	Int_t nHits = fHitArray->GetEntriesFast();
	cout<<"Hit array contains "<<nHits<< " hits"<<endl;
	for (Int_t iHit=0; iHit<nHits; iHit++) {
		theHit = (PndEmcHit*) fHitArray->At(iHit);
		int detId=theHit->GetDetectorId();

		// Both endcaps are already implemented w PndEmcMapper (21.05.07)
	        int module = detId/100000000;

		PndEmcWaveform *theWaveform = wf_map[detId];
		if (theWaveform==0) Fatal("PndEmcHitsToWaveform::Exec","No Detector Id in map");
		theWaveform->update_waveform(theHit);
	}

///////////////////// Sparsify //////////////////////
////////////////////////////////////////////////////
	// Sparsify waveform list, i.e keep only those waveforms which contain hits or their neigbours
	int detId_tmp;

	for (wf_iter=wf_map.begin();wf_iter!=wf_map.end();wf_iter++)
	{
		PndEmcWaveform* theWaveform = wf_iter->second;
		int detId=theWaveform->GetDetectorId();
				
		if (theWaveform->GetNHits()>0)
		  {
		    if (theWaveform->max()>fEnergyCutGeV*theWaveform->get_scale())
		      {
			keep_index.push_back(detId);
			
			PndEmcTwoCoordIndex* theIndex=theWaveform->GetTCI();
			PndEmcNeighbourStore theNeighbourStore=theIndex->itsNeighbours();
			PndEmcNeighbourStore::iterator theNeighbourIterator; 
			PndEmcTwoCoordIndex *theNeighbourIndex;
			for (theNeighbourIterator = theNeighbourStore.begin(); theNeighbourIterator != theNeighbourStore.end(); ++theNeighbourIterator ) 
			  {
			    theNeighbourIndex = (PndEmcTwoCoordIndex*)(*theNeighbourIterator); 
			    detId_tmp =theNeighbourIndex->itsIndex();
			    keep_index.push_back(detId_tmp);
			  }
			
		      }
		    else
		      {
			std::cout<<"Waveform failed energy cut "<<endl;
			std::cout<<"theWaveform->max()= "<<theWaveform->max()<<" fEnergyCut*theWaveform->get_scale()= "<<fEnergyCutGeV*theWaveform->get_scale()<<endl;
		      }
		  }
	}
	

	// Copy waveforms which passed cryteria to wf_vec
	for (wf_iter=wf_map.begin();wf_iter!=wf_map.end();++wf_iter)
	{
		PndEmcWaveform* theWaveform = wf_iter->second;
		int detId=theWaveform->GetDetectorId();
		keep_index_iter = find(keep_index.begin(),keep_index.end(),detId);
		if (keep_index_iter != keep_index.end()) // i.e. detId is in keep_index
		{
			wf_vec.push_back(theWaveform);
		}
   }


	//cout<<"The number of waveform left after sparsification = "<<wf_vec.size()<<endl;

//////////////// end of sparsification section //////
/////////////////////////////////////////////////////

	// Noise is added after sparsifycation to save calculation time
	// There are two options how to add noise (before and after shaping) 

	if (fUse_shaped_noise==0)
	{
		Int_t iwf=0;
		for (wf_vec_iter=wf_vec.begin();wf_vec_iter!=wf_vec.end();++wf_vec_iter)
		{
			PndEmcWaveform* tempWaveform = *wf_vec_iter;
			tempWaveform->add_elec_noise_and_digitise(fIncoherent_elec_noise_width_GeV*fGevPeakAnalogue,fOneBitResolution);
			new((*fWaveformArray)[iwf]) PndEmcWaveform(*tempWaveform); 
			iwf++;
		}
	}
	else
	{
		Int_t iwf=0;
		for (wf_vec_iter=wf_vec.begin();wf_vec_iter!=wf_vec.end();++wf_vec_iter)
		{
			PndEmcWaveform* tempWaveform = *wf_vec_iter;
			tempWaveform->add_shaped_elec_noise_and_digitise(fIncoherent_elec_noise_width_GeV*fGevPeakAnalogue,fOneBitResolution);
			new((*fWaveformArray)[iwf]) PndEmcWaveform(*tempWaveform); 
			iwf++;
		}
	}



	for(tciIter=tciMap.begin();tciIter!=tciMap.end();++tciIter)
	{	
		delete wf_map[tciIter->first];
	}


	wf_vec.clear();
	keep_index.clear();

	if (fVerbose>0){
		timer.Stop();
		Double_t rtime = timer.RealTime();
		Double_t ctime = timer.CpuTime();
		cout << "PndEmcHitsToWaveform, Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
	}


}

void PndEmcHitsToWaveform::SetParContainers() {

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

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

  // Get Emc digitisation parameter container
  fDigiPar = (PndEmcDigiPar*) db->getContainer("PndEmcDigiPar");
 
}


ClassImp(PndEmcHitsToWaveform)