//_HADES_CLASS_DESCRIPTION
////////////////////////////////////////////////////////////////
//
//  HEmcCalibrater:
//
//  Calibrates all fired cells in EmcRaw category and fills
//  the EmcCal category
//
//  embedding = 0 no embedding
//            = 1 embedding
////////////////////////////////////////////////////////////////

#include "hemccalibrater.h"
#include "emcdef.h"
#include "hemcraw.h"
#include "hemccalqa.h"
#include "hemccalsim.h"
#include "hemcdetector.h"
#include "hemccalpar.h"
#include "hemccellstatuspar.h"
#include "hemccalibraterpar.h"
#include "hstartdef.h"
#include "hstart2hit.h"
#include "hades.h"
#include "hcategory.h"
#include "hdebug.h"
#include "hevent.h"
#include "hiterator.h"
#include "hruntimedb.h"
#include "hrun.h"
#include "hspectrometer.h"

#include <TMath.h>

#include <iomanip>
#include <iostream>
#include <stdlib.h>
#include <string.h>
using namespace std;

#define SELECT_FIRST_OR_LAST // Select the first or the last pair, if commented, the last pair will be stored

ClassImp(HEmcCalibrater)
HEmcCalibrater::HEmcCalibrater(void) {
  // default constructor
  pRawCat = NULL;
  pCalCat = NULL;
  iter    = NULL;
  pCalpar = NULL;
  embedding = 0;
}

HEmcCalibrater::HEmcCalibrater(const Text_t *name, const Text_t *title)
                      : HReconstructor(name, title) {
  // constructor
  pRawCat = NULL;
  pCalCat = NULL;
  iter    = NULL;
  pCalpar = NULL;
  embedding = 0;
}

HEmcCalibrater::~HEmcCalibrater(void) {
  //destructor deletes the iterator on the raw category
  if (NULL != iter) {
    delete iter;
    iter = NULL;
  }
}

Bool_t HEmcCalibrater::init(void) {
  // gets the calibration parameter container
  // gets the EmcRaw category and creates the EmcCal category
  // creates an iterator which loops over all fired cells in EmcRaw
  HEmcDetector* pDet = (HEmcDetector*)gHades->getSetup()->getDetector("Emc");
  if (!pDet) {
    Error("init", "No Emc found.");
    return kFALSE;
  }

  pCalpar = (HEmcCalPar*)gHades->getRuntimeDb()->getContainer("EmcCalPar");
  if (!pCalpar) return kFALSE;

  pStatuspar = (HEmcCellStatusPar*)gHades->getRuntimeDb()->getContainer("EmcCellStatusPar");
  if (!pStatuspar) return kFALSE;

  pCalibpar = (HEmcCalibraterPar*)gHades->getRuntimeDb()->getContainer("EmcCalibraterPar");
  if (!pCalibpar) return kFALSE;

  pRawCat = gHades->getCurrentEvent()->getCategory(catEmcRaw);
  if (!pRawCat) {
    Warning("init()", "HEmcRaw category not available!");
  }

  if(gHades->getEmbeddingMode()>0) embedding = 1;

  pCalCat = (HCategory*)(((HEvent*)(gHades->getCurrentEvent()))->getCategory(catEmcCal));
  if(pCalCat == NULL) {
    if (embedding == 0) pCalCat = pDet->buildMatrixCategory("HEmcCal",0.5); 
    else                pCalCat = pDet->buildMatrixCategory("HEmcCalSim",0.5);
    if(pCalCat == NULL) return kFALSE;
    gHades->getCurrentEvent()->addCategory(catEmcCal,pCalCat,"Emc");
  }

  pStartHitCat = (HCategory*)(((HEvent*)(gHades->getCurrentEvent()))->getCategory(catStart2Hit));
  if(pStartHitCat == NULL) {
    Info("HEmcCalibrater::init()", "No Start2Hit category, no Start time correction\n");
  }

#ifdef WITHCALQA
  pQACat = (HCategory*)(((HEvent*)(gHades->getCurrentEvent()))->getCategory(catEmcCalQA));
  if(pQACat == NULL) {
    pQACat = pDet->buildMatrixCategory("HEmcCalQA",0.5);
    if(pQACat == NULL) return kFALSE;
    gHades->getCurrentEvent()->addCategory(catEmcCalQA, pQACat,"Emc");
  }
#endif

  if(pRawCat) iter = (HIterator*)pRawCat->MakeIterator();
  else iter = 0;

  loc.set(2, 0, 0);
  fActive = kTRUE;

  return kTRUE;
}

Float_t twcCorrection(Float_t t, Float_t * twcPars)
{
  if (t > twcPars[2])
    return twcPars[0] + twcPars[1]/TMath::Sqrt(t - twcPars[2]);

  return 0.0;
}

Int_t HEmcCalibrater::execute(void) {
  // calibrates all fired cells

  if(!pRawCat) return 0;

  HEmcRaw *pRaw = NULL;
  HEmcCal *pCal = NULL;

  Int_t  sec   = -1;
  Int_t  cell  = 0;
  Int_t  pos   = 0;
  Char_t row   = -1;
  Char_t col   = -1;
  UInt_t nfasthits = 0, nslowhits=0;

  Float_t rawTimeLeading   = 0.;
  Float_t rawTimeTrailing  = 0.;
  Float_t twcTimeLeading   = 0.;
  Float_t rawTimeLeadingS  = 0.;
  Float_t rawTimeTrailingS = 0.;
  Float_t rawWidth         = 0.;
  Float_t calTimeLeading   = 0.;

  Float_t parCellData[8] = {1.F, 0.F, 0.F, 1.F, 0.F, 0.F, 0.F, 0.F};

  //--------------------------------------------------
  // start time extraction. For embedding no extra work
  // has to be done
  Float_t start_time = 0.0;
  if (pStartHitCat && pStartHitCat->getEntries()>0) {
    HStart2Hit *pStartHit = (HStart2Hit *) pStartHitCat->getObject(0);
    if (pStartHit && pStartHit->getFlag())
        start_time = pStartHit->getTime();
  }
  //--------------------------------------------------

  Float_t diff_min = pCalibpar->getMatchWindowMin();
  Float_t diff_max = pCalibpar->getMatchWindowMax();

  Float_t accept_min = pCalibpar->getAcceptWindowMin();
  Float_t accept_max = pCalibpar->getAcceptWindowMax();

  iter->Reset();
  while ((pRaw = (HEmcRaw *)iter->Next()) != 0) {
    Float_t rawLead   = 0.F;
    Float_t rawTot  = 0.F;
    Float_t calTime   = 0.F;
    Float_t calEnergy = 0.F;

    nfasthits = pRaw->getFastMultiplicity();
    nslowhits = pRaw->getSlowMultiplicity();

    if (nfasthits == 0 && nslowhits == 0)
        continue;

    pRaw->getAddress(sec, cell, row, col);

    if (sec < 0)
        continue;

    pos = HEmcDetector::getPositionFromCell(cell);

    if (pStatuspar->getCellStatus(sec, pos) == 0)
        continue;

    loc[0] = sec;
    loc[1] = cell;

    // get the calibration parameters
    HEmcCalParCell &pPar = (*pCalpar)[sec][cell];
    pPar.getData(parCellData);

    Int_t hits_cnt = 0;

#ifdef WITHCALQA
    HEmcCalQA *pQA = (HEmcCalQA*)pQACat->getObject(loc);
    if (!pQA) {
      pQA = (HEmcCalQA *)pQACat->getSlot(loc);
      if (pQA) {
        pQA = new(pQA) HEmcCalQA;
        pQA->setAddress((Char_t)sec, (UChar_t)cell, row, col);
      } else {
        Error("execute()", "Can't get slot for CalQA sector=%i, cell=%i", sec, cell);
        return -1;
      }
    } else {
      Error("execute()", "Slot already exists for CalQA sector=%i, cell=%i", sec, cell);
      return -1;
    }
#endif
    // loop over number of fast hits in raw level
    UInt_t j = 0; // fast trailing
    UInt_t k = 0; // slow leading
    UInt_t l = 0; // slow trailing
    for (UInt_t i = 0; i < nfasthits; ++i)
    {
      rawTimeLeading = pRaw->getFastTimeLeading(i);
      if (0.0 == rawTimeLeading) continue;

#ifdef WITHCALQA
      for (; j < nfasthits; ++j)
      {
        rawTimeTrailing = pRaw->getFastTimeTrailing(j);
        if (0.0 == rawTimeTrailing) continue;

        if (rawTimeTrailing < rawTimeLeading) continue;
        pQA->addFastHit(rawTimeLeading * parCellData[0] + parCellData[1], (rawTimeTrailing - rawTimeLeading) * parCellData[0]);
        break;
      }
#endif

      for (; k < nslowhits; ++k)
      {
        rawTimeLeadingS = pRaw->getSlowTimeLeading(k);
        if (0.0 == rawTimeLeadingS) continue;

        Float_t diff = rawTimeLeadingS - rawTimeLeading;
        if (diff < diff_min) continue;      // slow leading to fast
        if (diff > diff_max) break;         // slow leading to late

#ifdef WITHCALQA
        for (UInt_t _l = l; _l < nslowhits; ++_l)
        {
          rawTimeTrailingS = pRaw->getSlowTimeTrailing(_l);
          if (0.0 == rawTimeTrailingS) continue;

          if (rawTimeTrailingS < rawTimeLeadingS) continue;
          pQA->addSlowHit(rawTimeLeadingS * parCellData[0] + parCellData[1], (rawTimeTrailingS - rawTimeLeadingS) * parCellData[0]);
          break;
        }
#endif

        // search for matching slow trailing
        for (; l < nslowhits; ++l)
        {
          rawTimeTrailingS = pRaw->getSlowTimeTrailing(l);
          if (0.0 == rawTimeTrailingS) continue;

          if (rawTimeTrailingS < rawTimeLeadingS) continue;

          rawWidth = rawTimeTrailingS - rawTimeLeading;

          // calc Time Walk Correction
          Float_t twc = twcCorrection(rawWidth, (parCellData+5));
          twcTimeLeading = rawTimeLeading - twc;

          // correct the time and width
          calTimeLeading = parCellData[0] * twcTimeLeading + parCellData[1];
          Float_t calWidth = rawWidth * parCellData[0];

#ifdef WITHCALQA
          pQA->addMatchedHit(parCellData[0] * rawTimeLeading + parCellData[1], calWidth);
#endif

          Float_t corrTime = calTimeLeading - start_time;

          /* If the leading edge is outside accepting range, completely ignore this pair. */
          if (corrTime < accept_min || corrTime > accept_max) break;

          ++hits_cnt;
#ifdef SELECT_FIRST_OR_LAST
          if (hits_cnt > 1) break;
#endif

          rawLead = rawTimeLeading;
          rawTot = rawWidth;
          calTime = corrTime;
          calEnergy = parCellData[2]+exp( parCellData[3]+ calWidth*parCellData[4]);

          break;
        }
        break;
      }
    }

    // store the last found hit
    if (hits_cnt)
    {
      pCal = (HEmcCal*)pCalCat->getObject(loc);
      if (!pCal) {
        pCal = (HEmcCal *)pCalCat->getSlot(loc);
        if (pCal) {
        if(embedding == 0) pCal = new(pCal) HEmcCal;
        else               pCal = new(pCal) HEmcCalSim;
        pCal->setAddress((Char_t)sec, (UChar_t)cell, row, col);
        } else {
        Error("execute()", "Can't get slot sector=%i, cell=%i", sec, cell);
        return -1;
        }
      } else {
        Error("execute()", "Slot already exists for sector=%i, cell=%i", sec, cell);
        return -1;
      }

      // fill cal hits
      pCal->setLead(rawLead);
      pCal->setWidth(rawTot);
      pCal->setTime(calTime);
      pCal->setEnergy(calEnergy);

      pCal->setNHits(hits_cnt);
    }
  }
  return 0;
}