//----------------------------------------------------------------------
// File and Version Information:
//      $Id: //
// Description:
//      Class EmcHitsToWaveform. 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 "EmcHitsToWaveform.h"
#include "CbmRootManager.h"
#include "CbmRunAna.h"
#include "CbmRuntimeDb.h"

#include "CbmEmcHit.h"
#include "EmcWaveform.h"
#include "EmcMapper.h"
#include "EmcStructure.h"
#include "EmcDigiPar.h"		
#include <iostream>
#include <map>
#include <string>

//#include "TStopwatch.h"

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

EmcHitsToWaveform::EmcHitsToWaveform(string fileGeo)
{
	fFileGeo=fileGeo;
	fEnergyCutGeV=1e-5; //GeV (to sparsify waveform list)
}

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

EmcHitsToWaveform::~EmcHitsToWaveform()
{
}


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

	nBits=fDigiPar->GetNBits();
	detectedPhotonsPerMeV=fDigiPar->GetDetectedPhotonsPerMeV();
	energyRange=fDigiPar->GetEnergyRange(); //GeV
	excessNoiseFactor=fDigiPar->GetExcessNoiseFactor();
	firstSamplePhase=fDigiPar->GetFirstSamplePhase ();
	number_of_samples_in_waveform=fDigiPar->GetNumber_of_samples_in_waveform ();
	Shaping_diff_time=fDigiPar->GetShaping_diff_time();     //s
	Shaping_int_time=fDigiPar->GetShaping_int_time();      //s
	crystal_time_constant=fDigiPar->GetCrystal_time_constant ();  //s
	incoherent_elec_noise_width_GeV=fDigiPar->GetIncoherent_elec_noise_width_GeV(); //GeV
	sampleRate=fDigiPar->GetSampleRate();
	use_shaped_noise=fDigiPar->GetUse_shaped_noise();
	use_photon_statistic=fDigiPar->GetUse_photon_statistic();
	
	// Test how parameters were read from DB.
	cout<<"nBits "<<nBits<<endl;
	cout<<"detectedPhotonsPerMeV "<<detectedPhotonsPerMeV<<endl;
	cout<<"energyRange "<<energyRange<<endl;
	cout<<"excessNoiseFactor "<<excessNoiseFactor<<endl;
	cout<<"firstSamplePhase "<<firstSamplePhase<<endl;
	cout<<"number_of_samples_in_waveform "<<number_of_samples_in_waveform<<endl;
	cout<<"Shaping_diff_time "<<Shaping_diff_time<<endl;
	cout<<"Shaping_int_time "<<Shaping_int_time<<endl;
	cout<<"crystal_time_constant "<<crystal_time_constant<<endl;
	cout<<"incoherent_elec_noise_width_GeV "<<incoherent_elec_noise_width_GeV<<endl;
	cout<<"sampleRate "<<sampleRate<<endl;
	cout<<"use_shaped_noise "<<use_shaped_noise<<endl;
	cout<<"use_photon_statistic "<<use_photon_statistic<<endl;

	// map version 1 corresponds to geometry from emc_module12345.dat
	// other options are not implemented yet
	int map_ver=0;	
	if (!fFileGeo.compare("emc_module12345.dat"))
		map_ver=1;
	else 
		Fatal("EmcHitsToWaveform::Init()", "Mapper for the given geometry file is not defined");
  
	EmcMapper::Instance(map_ver);
	EmcStructure::Instance(fFileGeo);


	// Calculate 1 bit resolution (in units of FADC amplitude)
	EmcWaveform *tmpwaveform;
	tmpwaveform=new EmcWaveform(0,101010001, Shaping_diff_time, Shaping_int_time, sampleRate, 0.0,0.0,crystal_time_constant,number_of_samples_in_waveform,excessNoiseFactor,false);

	fGevPeakAnalogue = tmpwaveform->get_scale();
	
	fOneBitResolution=energyRange/((double) (1<<nBits))*fGevPeakAnalogue;
	
	firstADCBinTime=firstSamplePhase/sampleRate;
	
	tciMap=EmcMapper::Instance(0)->GetTciMap();

	delete tmpwaveform;
	
	return kSUCCESS;
}

void EmcHitsToWaveform::Exec(Option_t* opt)
{
	// Reset output array
	if ( ! fWaveformArray ) Fatal("Exec", "No Waveform Array");
	fWaveformArray->Delete();
	// Variable declaration
	CbmEmcHit* theHit = NULL;

	for(tciIter=tciMap.begin();tciIter!=tciMap.end();++tciIter)
	{	
		 EmcWaveform* tempWaveform = new EmcWaveform(0,tciIter->first,
				       Shaping_diff_time,Shaping_int_time,
				       sampleRate,
				       firstADCBinTime,
				       detectedPhotonsPerMeV,
				       crystal_time_constant,
				       number_of_samples_in_waveform,
						 excessNoiseFactor);

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


	// Loop over CbmEmcHits 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 = (CbmEmcHit*) fHitArray->At(iHit);
		int detId=theHit->GetDetectorId();

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

		EmcWaveform *theWaveform = wf_map[detId];
		if (theWaveform==0) Fatal("EmcHitsToWaveform::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++)
	{
		EmcWaveform* theWaveform = wf_iter->second;
		int detId=theWaveform->GetDetectorId();
				
		if (theWaveform->GetNHits()>0)
		  {
		    if (theWaveform->max()>fEnergyCutGeV*theWaveform->get_scale())
		      {
			keep_index.push_back(detId);
			
			TwoCoordIndex* theIndex=theWaveform->GetTCI();
			NeighbourStore theNeighbourStore=theIndex->itsNeighbours();
			NeighbourStore::iterator theNeighbourIterator; 
			TwoCoordIndex *theNeighbourIndex;
			for (theNeighbourIterator = theNeighbourStore.begin(); theNeighbourIterator != theNeighbourStore.end(); ++theNeighbourIterator ) 
			  {
			    theNeighbourIndex = (TwoCoordIndex*)(*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)
	{
		EmcWaveform* 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 (!use_shaped_noise)
	{
		Int_t iwf=0;
		for (wf_vec_iter=wf_vec.begin();wf_vec_iter!=wf_vec.end();++wf_vec_iter)
		{
			EmcWaveform* tempWaveform = *wf_vec_iter;
			tempWaveform->add_elec_noise_and_digitise(incoherent_elec_noise_width_GeV*fGevPeakAnalogue,fOneBitResolution);
			new((*fWaveformArray)[iwf]) EmcWaveform(*tempWaveform); 
			iwf++;
		}
	}
	else
	{
		Int_t iwf=0;
		for (wf_vec_iter=wf_vec.begin();wf_vec_iter!=wf_vec.end();++wf_vec_iter)
		{
			EmcWaveform* tempWaveform = *wf_vec_iter;
			tempWaveform->add_shaped_elec_noise_and_digitise(incoherent_elec_noise_width_GeV*fGevPeakAnalogue,fOneBitResolution);
			new((*fWaveformArray)[iwf]) EmcWaveform(*tempWaveform); 
			iwf++;
		}
	}



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


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



}

void EmcHitsToWaveform::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 = (EmcDigiPar*) db->getContainer("EmcDigiPar");
 
}


ClassImp(EmcHitsToWaveform)