using namespace std;
#include <TRandom.h>
#include <time.h>
#include "hwalldigitizer.h"
#include "walldef.h"
#include "hdebug.h"
#include "hades.h"
#include "hiterator.h"
#include "hruntimedb.h"
#include "hspectrometer.h"
#include "hwalldetector.h"
#include "hgeantwall.h"
#include "hwallrawsim.h"
#include "hevent.h"
#include "hcategory.h"
#include "hlocation.h"
//#include "hwalldigipar.h"
#include "hwallrawsimfilter.h"
#include <iostream> 
#include <iomanip>

//*-- Author : D.Vasiliev
//*-- Modified: 07/04/2005 F.Krizek
//*-- Modified: 27/10/2002 D.Zovinec
//*-- Modified: 13/02/2002 by D.Zovinec
//*-- Modified: 30/11/2000 by R.Holzmann
//*-- Modified: 16/12/99 by J.Bielcik
//*-- Modified: 9/12/99 by V.Charriere
//*-- Modified: 8/11/99 by R.Holzmann
//*-- Modified: 24/10/99 by D.Vasiliev

//_HADES_CLASS_DESCRIPTION 
/////////////////////////////////////////////////////////////////////
//
//  HTofDigitizer digitizes data, puts output values into
//  raw data category
//
//  Need to implement quenching?
//
/////////////////////////////////////////////////////////////////////

 void HWallDigitizer::initParContainer() {

  // fDigiPar=(HWallDigiPar *)gHades->getRuntimeDb()->getContainer("WallDigiPar");
}

 HWallDigitizer::HWallDigitizer(void) {

  fGeantCat=0;
  fRawCat=0;
  fDigiPar=0;
  fLoc.set(1,0); 
  iterGeant=0;
  iterWallRaw=0;
}

 HWallDigitizer::HWallDigitizer(Text_t *name,Text_t *title) :
               HReconstructor(name,title) {

  fGeantCat=0;
  fRawCat=0;
  fDigiPar=0;
  fLoc.set(1,0);
  iterGeant=0;
  iterWallRaw=0;
}

 HWallDigitizer::~HWallDigitizer(void) {
}

 Bool_t HWallDigitizer::init(void) {

  time_t curtime;
  initParContainer();
  fGeantCat = gHades->getCurrentEvent()->getCategory(catWallGeantRaw);
  if (!fGeantCat) {
    fGeantCat = gHades->getSetup()->getDetector("Wall")->buildCategory(catWallGeantRaw);
    if (!fGeantCat) return kFALSE;
    else gHades->getCurrentEvent()->addCategory(catWallGeantRaw,fGeantCat,"Wall");
  }

  fRawCat = gHades->getCurrentEvent()->getCategory(catWallRaw);
  if (!fRawCat) {
     HWallDetector* wall=(HWallDetector*)(gHades->getSetup()->getDetector("Wall"));
     //??????????Matrix or linear category
     fRawCat=wall->buildMatrixCategory("HWallRawSim",0.5F);
     if (!fRawCat) return kFALSE;
     else gHades->getCurrentEvent()->addCategory(catWallRaw,fRawCat,"Wall");
  }

  iterGeant = (HIterator *)fGeantCat->MakeIterator("native");
  iterWallRaw = (HIterator *)fRawCat->MakeIterator("native");

  time(&curtime);

  return kTRUE;
}

 Int_t HWallDigitizer::execute(void) {

  //const Float_t deg2rad = 0.017453293; // deg to rad conversion.
  const Float_t quantEff = 0.24;       // PMT quantum efficiency.
  const Float_t photYield = 11100.0;   // photon yield in scintillator (phot/MeV).
  const Float_t timeResZero = 0.18;     // time resolution in the centre of the strip [ns].
  const Float_t relAmplResol = 0.08;   // sigma of Gaus distribution.
  const Float_t minEnerRelease = 1.8;  // minimum energy release (MeV/cm)

  HGeantWall* geant = 0;
  HWallRawSim* raw = 0;
  HWallRawSimFilter fRawFilter;

  Float_t slopeTDC;
  Int_t thrCFD, thrADC;

  Int_t   numTrack, numTrack1 = -1, numTrack2 = -1;      //track numbers
  Float_t trackLen;
  Float_t time, charge;
  Int_t timeCh, chargeCh;
  Float_t prevTime, prevCharge;
  Float_t timeResol, amplResol, chargeRef;
 

  Float_t geaTof = 0.;
  Float_t geaTof1[383];
  Float_t geaTof2[383];
  Float_t geaEner = 0.;
  Float_t geaX = 0.;
  Float_t geaY = 0.;     // used by Tofino only
  Float_t geaMom = 0.;

  iterGeant->Reset();   // this works only in split mode=2
                   // (for 0 and 1 the iterator must be recreated)

  while ((geant=(HGeantWall *)iterGeant->Next())!=0) {
    
    fLoc[0] = geant->getCell();   

    raw = (HWallRawSim*) fRawCat->getObject(fLoc);   // test if cell in use
    if(raw) {
       raw->incNHit();  // increment counter for additional hits
       numTrack1 = raw->getNTrack1();
       numTrack2 = raw->getNTrack2();
       prevTime = raw->getTime();    
       prevCharge = raw->getCharge();       
    }
    else {
       prevTime = 100000.;
       prevCharge = 0.;
       raw = (HWallRawSim*) fRawCat->getNewSlot(fLoc);  // get new slot
       if(!raw) continue;
       raw = new(raw) HWallRawSim;
       raw->setNHit(1);
    }

    //HWallDigiParCell& cell=(*fDigiPar)[fLoc[0]];     
    // slopeTDC = cell.getTDCSlope();
    //thrCFD = cell.getCFDThreshold();
    //thrADC = cell.getADCThreshold();
       
      slopeTDC = 0.065;
      thrCFD = 1;
      thrADC = 1;

    geant->getHit(geaEner,geaX,geaY,geaTof,geaMom,trackLen);

    numTrack=geant->getTrack();

      {
	time = geaTof;
        
	timeResol = timeResZero;
	time = gRandom->Gaus(time,timeResol);            


        timeCh = (Int_t) (time/slopeTDC);
        if (timeCh < 0) timeCh = 0;
	if (timeCh > 4095) timeCh = 4095;
        time = ( Float_t ) (timeCh);

        if(raw->getNHit()>1){
          if(geaTof < geaTof1[fLoc[0]]){
            numTrack2=numTrack1;
            numTrack1=numTrack;
            geaTof2[fLoc[0]]=geaTof1[fLoc[0]];
            geaTof1[fLoc[0]]=geaTof;
          } else {
            if(geaTof < geaTof2[fLoc[0]]){
              numTrack2=numTrack;              
              geaTof2[fLoc[0]]=geaTof;
            }
          }
        }

        if(time > prevTime) time = prevTime;        

	//charge = geaEner*photYield*quantEff*0.5*(1 - cos(ar*deg2rad))*exp(-(hl-geaX)/al);
        charge = geaEner*photYield*quantEff;
	amplResol = charge*relAmplResol;
	charge = gRandom->Gaus(charge,amplResol);
	
	chargeRef = photYield*quantEff*minEnerRelease;

	//if(fLoc[1] <= 3) chargeRef *= 3.;
	//if(fLoc[1]>3 && fLoc[1]<=7) chargeRef *= 2.;

	chargeCh = (Int_t) ((charge/chargeRef)*256. + prevCharge);
        if (chargeCh < 0) chargeCh = 0;
	if (chargeCh > 4095)  chargeCh = 4095;
	charge = (Float_t)chargeCh;
      

      } //else {

	//time = 4095.;      
	//charge = 4095.;    
      // }
    raw->setTime(time);   
    raw->setCharge(charge);      
    raw->setCell((Int_t) fLoc[0]);     
    if(raw->getNHit()>1){
      raw->setNTrack1(numTrack1);
      raw->setNTrack2(numTrack2);
    } else {
      raw->setNTrack1(numTrack);
      raw->setNTrack2(-1);
      geaTof1[fLoc[0]]=geaTof;
      geaTof2[fLoc[0]]=100000.;
    }
  }

  // Exclusion of hits with charge less than ADC threshold.
  // Time of hits with charge less than CFD threshold is set
  // to zero. These hits are excluded later in the hitfinder.

  iterWallRaw->Reset();
  while ( (raw=(HWallRawSim *)iterWallRaw->Next())!=NULL) {  
    fLoc[0] = raw->getCell();
    time=raw->getTime();
    charge=raw->getCharge();
 

    // overflow suppression
    if(((Int_t)time)>=4095) raw->setTime(0.0);
  
    //  HWallDigiParCell& cell=(*fDigiPar)[fLoc[0]];
    thrCFD = 50;//cell.getCFDThreshold();
    thrADC = 1;// cell.getADCThreshold();  

    // CFD and ADC thresholds
    if(((Int_t)charge)<thrCFD){
      raw->setTime(0.0);
      if(((Int_t)charge)<thrADC){
        raw->setCharge(0.0);
      }

    }
    
  }

  fRawCat->filter(fRawFilter);

  return 0;
}

ClassImp(HWallDigitizer)