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 "hwalldigipar.h"
#include "hwallgeompar.h"
#include "hgeantwall.h"
#include "hwallrawsim.h"
#include "hevent.h"
#include "hcategory.h"
#include "hlocation.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

// Modified by M.Golubeva 01.11.2006


//_HADES_CLASS_DESCRIPTION 

/////////////////////////////////////////////////////////////////////

//

//  HTofDigitizer digitizes data, puts output values into

//  raw data category

//

//  Need to implement quenching?

//

/////////////////////////////////////////////////////////////////////


void HWallDigitizer::initParContainer() {

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

HWallDigitizer::HWallDigitizer(void) {

  fGeantCat=0;
  fRawCat=0;
  pWallDigiPar=0;
  pWallGeomPar = NULL;
  fLoc.set(1,0); 
  iterGeant=0;
  iterWallRaw=0;
}

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

  fGeantCat=0;
  fRawCat=0;
  pWallDigiPar=0;
  pWallGeomPar = NULL;
  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"));
     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 timeResWall[3] = 
                 {0.42, 0.52, 0.67};  //values based on apr07 analysis

               //{0.35, 0.42, 0.60};  // time resolution in the centres of the 

                                       //small(40x40mm), middle (80x80mm) and

                                       //large (160x160mm) cells [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;

  TString cellName,str;
  Float_t slopeTDC;
  Int_t thrCFD, thrADC;
  Int_t cell=-1;
  fLoc.set(1,-1);

  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);
    }

    cell=fLoc[0];
    //HWallDigiParCell &pPar=(*pWallDigiPar)[cell];

    //slopeTDC = pPar.getTDC_Slope();

    slopeTDC = pWallDigiPar->getTDC_Slope(cell);
    //slopeTDC = 0.065;

    
    geant->getHit(geaEner,geaX,geaY,geaTof,geaMom,trackLen);
    
    numTrack=geant->getTrack();
    
    cellName=pWallGeomPar->getCellName(fLoc[0]);
    //cout <<fLoc[0] <<" " <<cellName <<endl;


    str=cellName.Replace(0,3,"",0);
    str.Replace(1,3,"",0);
    
    time = geaTof;
    
    //timeResol = timeResZero;

    timeResol = timeResWall[atoi(str)-1];
    //cout <<cell <<" " <<str <<" " <<atoi(str) <<" " <<timeResol  <<endl;

    time = gRandom->Gaus(time,timeResol);            
    
    
    timeCh = (Int_t) (time/slopeTDC);
    if (timeCh < 0) timeCh = 0;
    //if (timeCh > 40000) timeCh = 40000;

    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 > 40000)  chargeCh = 40000;

    if (chargeCh > 4095)  chargeCh = 4095;
    charge = (Float_t)chargeCh;
    
    
    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);
    
    cell=fLoc[0];
    //HWallDigiParCell &pPar=(*pWallDigiPar)[cell];

    //thrCFD = (Int_t)pPar->getCFD_Threshold();

    //thrADC = (Int_t)pPar->getADC_Threshold();  

    //thrCFD = 50;

    //thrADC = 1;


    thrCFD = (Int_t)pWallDigiPar->getCFD_Threshold(cell);
    thrADC = (Int_t)pWallDigiPar->getADC_Threshold(cell);  
    // 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)