using namespace std;

#include "hitofdigitizer.h"

#include "hades.h"
#include "hruntimedb.h"
#include "hspectrometer.h"
#include "hdetector.h"
#include "hitofdetector.h"

#include "hrecevent.h"
#include "hlocation.h"
#include "hpartialevent.h"
#include "hcategory.h"
#include "hlinearcategory.h"
#include "hmatrixcategory.h"


#include "hitofdef.h"
#include "hstartdef.h"
#include "hstart2hit.h"
#include "hgeantkine.h"
#include "hitofcalsim.h"


#include "htool.h"
#include "hiterator.h"
#include "hgeanttof.h"
#include "hevent.h"

#include "hitofdigitpar.h"

#include "TRandom.h"
#include "TMath.h"

#include <time.h>
#include <iostream>
#include <iomanip>
#include <algorithm>


//_HADES_CLASS_DESCRIPTION 
/////////////////////////////////////////////////////////////////////
//
//  HiTofDigitizer digitizes HGeant data, puts output values into
//  cal data category
//  The input data are read from the HGeantTof category. After calculating
//  TOF and Charge of the hit etc the output is stored in the HiTofCalSim
//  category.
//
//   ----------------------                                 -----------
//  |   HiTofCalibrater    |                               | HGeantTof |
//  |   (embedding mode)   | \                              -----------
//  |                      |  \                                 ||
//   ----------------------    \                                || input to digitizer
//                              \                               \/
//                           ---------------    read real     ------------------
//                          | HiTofCalSim   | ------------>  | HiTofDigitizer   |
//                           ---------------  <-----------   |                  |
//                                            write output    ------------------
//
//
//  In the case of TRACK EMBEDDING of simulated tracks into
//  experimental data the real data are written by the HiTofCalibrater into
//  HiTofCalSim category. The real hits are taken into
//  account by the digitizer (adding of charges, sorting by first hit in
//  photo multiplier). The embedding mode is recognized
//  automatically by analyzing the
//  gHades->getEmbeddingMode() > 0 flag.
//
//  The output object HiTofCalSim contains additional to the tof and
//  Eloss a list of up to 5 tracks hitting the itof cell. The track numbers
//  are stored by increasing tof.
//
// For the suppression of secondaries produced in the TOF itself there are
// several configuration possibilities:
// They can be switched via setStoreFirstTrack(Int_t flag):
//     flag ==  0 realistic (secondaries included)
//              1 primary particle is stored
//              2 (default) the first track number entering the tof in SAME SECTOR is stored
//              3 as 2 but SAME SECTOR && CELL
//
/////////////////////////////////////////////////////////////////////

ClassImp(HiTofDigitizer)

HiTofDigitizer*  HiTofDigitizer::piTofDigi = NULL;

HiTofDigitizer::HiTofDigitizer(const Text_t *name,const Text_t *title) :
    HReconstructor(name,title)
{

  fGeantCat      = NULL;
  fCalCat        = NULL;
  fStartHitCat   = NULL;
  piTofDigitPar  = NULL;
  piTofGeomPar   = NULL;
  pSpecGeomPar   = NULL;
  fLoc.set(3,0,0,0);

  iterGeant      = NULL;
  iterCal        = NULL;
  piTofDigi      = this;
  storeFirstTrack= 2;
  fstartTimeSmearing = 0;
}

HiTofDigitizer::~HiTofDigitizer(void)
{

    if(iterGeant) delete iterGeant;
    if(iterCal)   delete iterCal;
}

Bool_t HiTofDigitizer::sortTime(tofmeasurement& a,tofmeasurement& b)
{
    return (a.tdcTime < b.tdcTime);
}

Bool_t HiTofDigitizer::init(void)
{

    piTofDigitPar = (HiTofDigitPar *) gHades->getRuntimeDb()->getContainer("iTofDigitPar");
    if(piTofDigitPar == NULL)
    {
	Error("Error in HiTofDigitizer::init()","Container iTofDigitPar not initialized");
	return kFALSE;
    }

    piTofGeomPar = (HiTofGeomPar *) gHades->getRuntimeDb()->getContainer("iTofGeomPar");
    if(piTofGeomPar == NULL)
    {
	Error("Error in HiTofDigitizer::init()","Container iTofGeomPar not initialized");
	return kFALSE;
    }

    pSpecGeomPar = (HSpecGeomPar*) gHades->getRuntimeDb()->getContainer("SpecGeomPar");
    if(pSpecGeomPar == NULL)
    {
	Error("Error in HiTofDigitizer::init()","Container SpecGeomPar not initialized");
	return kFALSE;
    }

    fGeantCat = gHades->getCurrentEvent()->getCategory(catTofGeantRaw);
    if (!fGeantCat)
    {   Info("HiTofDigitizer::init()","No catTofGeantRaw in input!");

    } else iterGeant = (HIterator *)fGeantCat->MakeIterator("native");


    fGeantKineCat=(HLinearCategory*)gHades->getCurrentEvent()->getCategory(catGeantKine);
    if(!fGeantKineCat){
	Error("HiTofDigitizer::init()","No catGeantKine in input!");
	return kFALSE;
    }

    fCalCat = gHades->getCurrentEvent()->getCategory(catiTofCal);
    if (!fCalCat)
    {   HiTofDetector* itofdet = ((HiTofDetector*)(gHades->getSetup()->getDetector("iTof")));
	fCalCat = itofdet->buildMatrixCategory("HiTofCalSim");
	if (!fCalCat) {
	    Error("HiTofDigitizer::init()","Could not create catiTofCal!");
	    return kFALSE;
	}
	else gHades->getCurrentEvent()->addCategory(catiTofCal,fCalCat,"iTof");

    } else iterCal   = (HIterator *)fCalCat  ->MakeIterator("native");



    fStartHitCat = gHades->getCurrentEvent()->getCategory(catStart2Hit);
    if (!fStartHitCat){
	Info("HiTofDigitizer::init()","No HStart2HitCat in input!");
    }

    return kTRUE;
}

Bool_t HiTofDigitizer::reinit(void)
{

    if(piTofGeomPar && pSpecGeomPar) {
	piTofGeomPar ->initGeomHelpers(pSpecGeomPar);
    }

    return kTRUE;
}

Int_t HiTofDigitizer::execute(void)
{
    const Int_t embeddingmode = gHades->getEmbeddingMode(); // flag embedding 1=realistic 2:keep geant

    //-------------------------------------------------
    // loop over raw sim category and set the default values
    // for embedding mode the real data
    if(embeddingmode > 0){
	HiTofCalSim* pCal = NULL;
	iterCal->Reset();
	while ((pCal = (HiTofCalSim *)iterCal->Next()) != 0)
	{
	    pCal->addTrack(gHades->getEmbeddingRealTrackId(),pCal->getTime());  // ???? mean or fast time ?
	}
    }
    //-------------------------------------------------

    //------------------------------------------------------------------------
    // getting start time
    fstartTimeSmearing = 0.0;  // for simulation/embedding of sim tracks  pStartH->getResolution()
                               // startTime  = startTimeSmearing  for simulation
                               // startTime != startTimeSmearing  for embedding
                               // startTimeSmearing == -1000      for real data
    if (fStartHitCat && fStartHitCat->getEntries() > 0 ) {
	HStart2Hit* pStartH = 0;
	if((pStartH = (HStart2Hit *) fStartHitCat->getObject(0)) != NULL) {
	    if(pStartH->getResolution() != -1000) fstartTimeSmearing = pStartH->getResolution();
	}
    }
    //------------------------------------------------------------------------

    //-------------------------------------------------
    // fill temporary array from geant.
    fillArray();
    //-------------------------------------------------
    fillOutput();
    //-------------------------------------------------

    return 0;
}

Bool_t HiTofDigitizer::finalize(void)
{
    return kTRUE;
}

Float_t HiTofDigitizer::calcDistToPMT(Float_t distPMT[],Int_t s,Int_t c,Float_t geaX,Float_t geaY)
{
    // calculate dist to all PMT by geomtry + geant Hits in cell system
    // returns the shortest dist;
    HGeomVector p;

    p.setXYZ(geaX,geaY,0);
    piTofGeomPar->distPMT(distPMT,p,s,c);

    Float_t shortDist = 10000;
    for(Int_t i = 0 ; i < ITOFNPMT; i++){
       if(distPMT[i] < shortDist) shortDist = distPMT[i];
    }

    return shortDist;
}

void  HiTofDigitizer::fillArray()
{
    //-------------------------------------------------
    // fill temporary array from geant.
    // all sim data will be collected

    for(Int_t s = 0; s < 6; s++){
	for(Int_t cell = 0; cell < ITOFMAXCELL; cell++){
	    fcellhits[s][cell].vmeasurements.clear();
	}
    }
    //-------------------------------------------------
    // loop over geant category
    iterGeant->Reset();
    HGeantTof* pGeant = NULL;
    Float_t geaTof, geaEner, geaX, geaY, geaMom, trackLen;  // Geant output

    while ( (pGeant = (HGeantTof *)iterGeant->Next()) != NULL)
    {  // loop over HGeant hits

	fLoc[1] = pGeant->getModule();
	if(fLoc[1] < 22) continue;     // skip TOF hits

	fLoc[0] = pGeant->getSector();
	fLoc[2] = fLoc[1] - 22; // to start counting from 0

	fLoc[1] = 0;

	tofmeasurement mtof;

	//-------------------------------------------------
	// get GEANT values of the hit
	pGeant->getHit(geaEner,geaX,geaY,geaTof,geaMom,trackLen);
	Float_t shortDistPMT = calcDistToPMT(mtof.distPMTAll,fLoc[0],fLoc[2],geaX,geaY); // ????
	//-------------------------------------------------

        mtof.s        = fLoc[0];
        mtof.c        = fLoc[2];
        mtof.numTrack = pGeant->getTrack();
        mtof.geaTof   = geaTof;
        mtof.geaEloss = geaEner;
        mtof.geaX     = geaX;
        mtof.geaY     = geaY;
        mtof.distPMT  = shortDistPMT;

	//-------------------------------------------------
	// find the first track ID entering the TOF
	// depending on storeFirstTrack flag
	// if replacing of tof track number is used
	HGeantTof* pOld = pGeant;
	HGeantTof* pNew = pGeant;
	mtof.numFirstTrack = findFirstHitIniTof(pOld,&pNew);
	//-------------------------------------------------

	//-------------------------------------------------
	// calulation of measured times
	mtof.propTime     = piTofDigitPar->calcLightPropTime(fLoc[0],fLoc[2],mtof.distPMT);  //????
        Float_t totalTime = mtof.geaTof + mtof.propTime;
	Float_t timeResol = piTofDigitPar->calcTimeResolution(fLoc[0],fLoc[2],mtof.distPMT); //????
	totalTime         = gRandom->Gaus(totalTime,timeResol);

        mtof.resSmear     = totalTime - mtof.geaTof - mtof.propTime;
        mtof.startTime    = fstartTimeSmearing ;
        mtof.tdcTime      = totalTime - fstartTimeSmearing ;
	mtof.eloss        = piTofDigitPar->calcEnergyLoss(fLoc[0],fLoc[2],mtof.distPMT,mtof.geaEloss);
	//-------------------------------------------------

	fcellhits[fLoc[0]][fLoc[2]].vmeasurements.push_back(mtof);

    } // end loop over geant category

    // sort by fastest time
    for(Int_t s = 0; s < 6; s++){
	for(Int_t cell = 0; cell < ITOFMAXCELL; cell++){
            sort(fcellhits[s][cell].vmeasurements.begin(),fcellhits[s][cell].vmeasurements.end(),this->sortTime);
	}
    }
    /*
    for(Int_t s = 0; s < 6; s++){
	for(Int_t cell = 0; cell < ITOFMAXCELL; cell++){
	    if(fcellhits[s][cell].vmeasurements.size() > 0) {
		for(UInt_t i = 0; i < fcellhits[s][cell].vmeasurements.size(); i++){
		    fcellhits[s][cell].vmeasurements[i].print();
		}
	    }
	}
    } */
}

void HiTofDigitizer::fillOutput()
{
    // write to output category (and merge with real data)
    // fill list of tracks/geaTof other than the fastest hit

    const Int_t embeddingmode = gHades->getEmbeddingMode(); // flag embedding 1=realistic 2:keep geant

    for(Int_t s = 0; s < 6; s++){
	for(Int_t c = 0; c < ITOFMAXCELL; c++)
	{
            vector<tofmeasurement>& vm = fcellhits[s][c].vmeasurements;
	    if(vm.size() == 0) continue;

	    fLoc[0] = s;
	    fLoc[1] = 0;
	    fLoc[2] = c;

            Float_t fastTimeGeant  = vm[0].tdcTime;
            Float_t fastTrackGeant = vm[0].numFirstTrack;

            Float_t sumEloss    = 0;
            Float_t sumGEloss   = 0;
	    Float_t sumSumEloss = 0;

            Bool_t isGood = kFALSE;
	    for (UInt_t i = 0; i < vm.size(); i ++) {
		sumEloss    += vm[i].eloss;
		sumSumEloss += vm[i].eloss;
		sumGEloss   += vm[i].geaEloss;
	    }

            // smearing energy deposit with gausian resolution
            // ??? TODO: should we allow to go negative ?
            // For the moment those values would be cut by min energy cut
	    sumEloss    = gRandom->Gaus(sumEloss,piTofDigitPar->getSigmaEloss());
	    sumSumEloss = sumEloss;

            if(sumEloss > piTofDigitPar->getMinEnergyCut()) isGood = kTRUE;

	    //-------------------------------------------------
	    // embedding mode:
	    // get TofCalSim object from output category
	    // if the same adress object exists already sorting has
	    // to be done to find the earliest hit
	    HiTofCalSim* pCalOut = (HiTofCalSim*) fCalCat->getObject(fLoc);
	    if (pCalOut)
	    {
                HiTofCalSim pCalOld(*pCalOut);
		// test if this cell has already seen a hit
		// this should happen only in embedding mode
		// when the cell has been hit by real data before
		// TODO : 1.  compare times to set the fastest time
                //        2.  sum Eloss of hits
                //        3.  add track numbers = geaTof vals (depending on embeding mode)

		if(embeddingmode > 0)
		{

		    if(pCalOut->getEloss() > 0 )      sumEloss    += pCalOut->getEloss();
		    if(pCalOut->getElossSummed() > 0 )sumSumEloss += pCalOut->getElossSummed();

		    if(fastTimeGeant < pCalOld.getTime()) {
			pCalOut->setTime    (fastTimeGeant);  // fast time
                        pCalOut->setMeanTime(fastTimeGeant);  // ????

			if(embeddingmode != 2)// not keep geant mode
			{
			    pCalOut->fillTrackData(0,fastTrackGeant, fastTimeGeant);
			}
		    }

		    UInt_t start = 0;
                    if(embeddingmode != 2) start = 1; // fastest track already set

		    Bool_t doneReal = kFALSE;
		    for (UInt_t i = start; i < vm.size(); i ++) {

			if(embeddingmode != 2 && !doneReal && vm[i].tdcTime > pCalOld.getTime()){
			    doneReal = kTRUE;
			    pCalOut->addTrack(gHades->getEmbeddingRealTrackId(),pCalOld.getTime());
			}

			if( pCalOut->getNHit() < pCalOut->getMaxTrack()) { // stop writing track infos if array is full
			    pCalOut->addTrack(vm[i].numFirstTrack,vm[i].tdcTime);
			}
		    }

		    pCalOut->setEloss      (sumEloss);
                    pCalOut->setElossSummed(sumSumEloss);

		}
	    }
	    else
	    {
		// this cell has not yet been hit, ask for a new slot in cat
		// this case should be true in case we are running pure
		// simulation or cell has not not been hit by real data
		// in both cases we have to simply copy the sim data
		// from the temporay array to the output

		if(!isGood) continue;

		pCalOut = (HiTofCalSim*) fCalCat->getSlot(fLoc);  // get new slot
		if(!pCalOut)
		{
		    Error("execute()","Could not retrieve new Slot in Category catiTofCal!");
		    continue;
		}


		pCalOut = new(pCalOut) HiTofCalSim();
		pCalOut->setAddress     (s,0,c);
		pCalOut->setTime        (fastTimeGeant);
		pCalOut->setMeanTime    (fastTimeGeant);
		pCalOut->setMultiplicity(1);
		pCalOut->setEloss       (sumEloss);
		pCalOut->setElossSummed (sumSumEloss);

		for (UInt_t i = 0; i < vm.size(); i ++) {
		    if( pCalOut->getNHit() < pCalOut->getMaxTrack()) { // stop writing track infos if array is full
			pCalOut->addTrack(vm[i].numFirstTrack,vm[i].tdcTime);
		    }
		}


	    } // cell
	    //-------------------------------------------------
	}
    } // sec
    //-------------------------------------------------
}

Int_t HiTofDigitizer::findFirstHitIniTof(HGeantTof* poldTof,HGeantTof** pnewTof)
{
    //-------------------------------------------------
    // find the first track ID entering the iTOF
    // Used to suppress the secondaries created in the
    // iTOF itself.
    //        0 realistic (secondaries included)
    //        1 primary particle is stored
    //        2 (default) the first track number entering the tof in SAME SECTOR is stored
    //        3 as 2 but SAME SECTOR && CELL

    Int_t numTrack = poldTof->getTrack();

    if(numTrack <= 0) return numTrack; // nothing to do for real data
    HGeantKine* kine = 0;

    //--------------------------------------------------------
    // return the track number for
    // the selected option storeFirstTrack

    //--------------------------------------
    // case 0
    if(storeFirstTrack == 0) return numTrack;
    //--------------------------------------
    // case 1
    if(storeFirstTrack == 1)
    {   // return track number of primary particle
	// of the given track
	kine = (HGeantKine*)fGeantKineCat->getObject(numTrack - 1);
	kine = HGeantKine::getPrimary(numTrack,fGeantKineCat);
	return kine->getTrack();
    }

    //--------------------------------------
    // case 2 and 3
    kine = (HGeantKine*)fGeantKineCat->getObject(numTrack - 1);
    if(kine->getParentTrack() == 0){ return numTrack; } // nothing to do

    Int_t s,c;
    s = poldTof->getSector();
    c = poldTof->getModule() - 22;

    Int_t first = 0;
    Int_t tempTrack = numTrack;
    while((kine = kine->getParent(tempTrack,fGeantKineCat)) != 0)
    {
	first = kine->getFirstTofHit();
	if(first != -1)
	{
	    // track is still in TOF/iTOF

	    // now we have to check if it is in TOF or iTOF
	    HGeantTof* gtof = (HGeantTof*)fGeantCat->getObject(first);

	    Int_t s1,c1;
	    s1 = c1 = 0;
	    if(gtof->getModule() >= 22)
	    {
		// inside iTOF
		s1 = gtof->getSector();
		c1 = gtof->getModule() - 22;

		if(storeFirstTrack == 2 &&
		   s == s1)
		{   // case 2 :: check only sector
		    tempTrack  = kine->getTrack();
		    (*pnewTof) = gtof;
		}
		else if(storeFirstTrack == 3 &&
			s == s1 && c == c1)
		{   // case 3 :: check for same cell
		    tempTrack  = kine->getTrack();
		    (*pnewTof) = gtof;
		}
		else {
		    // track has no iTOF hit any longer
		    // which fulfills the condition
		    break;
		}

	    } else {
		// track is in TOF
		break;
	    }
	}
	else {
	    // track has no iTOF hit any longer,
	    // so the previous object was the one we
	    // searched  for
	    break;
	}
    }

    return tempTrack;
}