//*-- Author   : Pablo Cabanelas
//*-- Created  : 30/08/07
//*-- Modified : 18/12/07 Diego Gonzalez-Diaz
//*-- Modified : 13/12/09 Diego Gonzalez-Diaz 
//(removed hard-coded parameters)
//
//_HADES_CLASS_DESCRIPTION 
////////////////////////////////////////////////////////////////////
//
//  HRpcHitF: RPC detector hit finder
//
//  Gets calibrated data from RpcCal(Sim) and writes to RpcHit(Sim)
//
///////////////////////////////////////////////////////////////////

using namespace std;
#include "TRandom.h"
#include <time.h>
#include "hrpcdetector.h"
#include "hrpchitf.h"
#include "rpcdef.h"
#include "hrpccal.h"
#include "hrpccalsim.h"
#include "hrpchitsim.h"
#include "hgeantrpc.h"
#include "hrpcgeomcellpar.h"
#include "hrpcgeompar.h"
#include "hrpccellstatuspar.h"
#include "hrpcdigipar.h"
#include "hdebug.h"
#include "hades.h"
#include "hiterator.h"
#include "hruntimedb.h"
#include "hspectrometer.h"
#include "hspecgeompar.h"
#include "hdetector.h"
#include "hevent.h"
#include "hgeantkine.h"
#include "hcategory.h"
#include "hstart2detector.h"
#include "hstart2hit.h"
#include "hstartdef.h"
#include "hgeomvector.h"
#include "hgeomvector2.h"
#include "hgeomvolume.h"
#include "hgeomtransform.h"
#include "hrpcwtoqpar.h"
#include "hrpchitfpar.h"

#include <iostream>
#include <iomanip>

void HRpcHitF::initParContainer() {
    pRpcGeometry   = (HRpcGeomPar*)(gHades->getRuntimeDb()->getContainer("RpcGeomPar"));
    pGeomCellPar   = (HRpcGeomCellPar*)(gHades->getRuntimeDb()->getContainer("RpcGeomCellPar"));
    pCellStatusPar = (HRpcCellStatusPar*)(gHades->getRuntimeDb()->getContainer("RpcCellStatusPar"));
    pDigiPar       = (HRpcDigiPar*)(gHades->getRuntimeDb()->getContainer("RpcDigiPar"));
    pWtoQPar       = (HRpcWtoQPar*)(gHades->getRuntimeDb()->getContainer("RpcWtoQPar"));
    pHitFPar       = (HRpcHitFPar*)(gHades->getRuntimeDb()->getContainer("RpcHitFPar"));
    pSpecGeomPar   = (HSpecGeomPar*)(gHades->getRuntimeDb()->getContainer("SpecGeomPar"));
}

HRpcHitF::HRpcHitF(void)
{
    pCalCat        = NULL;
    pHitCat        = NULL;
    iter           = NULL;
    pRpcGeometry   = NULL;
    pGeomCellPar   = NULL;
    pCellStatusPar = NULL;
    pDigiPar       = NULL;
    pWtoQPar       = NULL;
    pHitFPar       = NULL;
    pStartHitCat   = NULL;
}

HRpcHitF::HRpcHitF(const Text_t *name, const Text_t *title) :
HReconstructor(name,title)
{
    pCalCat        = NULL;
    pHitCat        = NULL;
    iter           = NULL;
    pRpcGeometry   = NULL;
    pGeomCellPar   = NULL;
    pCellStatusPar = NULL;
    pDigiPar       = NULL;
    pWtoQPar       = NULL;
    pHitFPar       = NULL;
    pStartHitCat   = NULL;
}

HRpcHitF::~HRpcHitF(void)
{
    if (iter) delete iter;
    iter=NULL;
}

Bool_t HRpcHitF::init(void)
{

    initParContainer();

    if (!pRpcGeometry){
	Error("init","No RPC Geometry!!");
	return kFALSE;
    }
    if (!pGeomCellPar){
	Error("init","No RpcGeomCellPar Parameters");
	return kFALSE;
    }
    if (!pCellStatusPar){
	Error("init","No RpcCellStatusPar Parameters");
	return kFALSE;
    }
    if (!pDigiPar){
	Error("init","No RpcDigiPar Parameters");
	return kFALSE;
    }
    if (!pWtoQPar){
	Error("init","No RpcWtoQPar Parameters");
	return kFALSE;
    }
    if (!pHitFPar){
	Error("init","No RpcHitFPar Parameters");
	return kFALSE;
    }
    if (!pSpecGeomPar){
	Error("init","No SpecGeomPar Parameters");
	return kFALSE;
    }

    HRpcDetector *rpc;
    rpc=(HRpcDetector *)gHades->getSetup()->getDetector("Rpc");

    if(!rpc){
	Error("init","No Rpc Detector found");
	return kFALSE;
    }

    // Create category for Hit level of Start detector
    HStart2Detector *start;
    start=(HStart2Detector *)gHades->getSetup()->getDetector("Start");

    if (start) {
	pStartHitCat=gHades->getCurrentEvent()->getCategory(catStart2Hit);
	if (!pStartHitCat) Warning("init","Start hit level not defined; setting start time to 0");
    } else {
	Warning("init","Start detector not found; setting start time to 0");
	pStartHitCat=0;
    }

    // Create category for calibrated data
    pCalCat=gHades->getCurrentEvent()->getCategory(catRpcCal);
    if (!pCalCat) {
	Error("init","No RpcCal Category");
	return kFALSE;
    }

    // Decide whether we are running for simulation or real data
    HCategory* catKin=gHades->getCurrentEvent()->getCategory(catGeantKine);

    if(catKin) simulation=kTRUE;
    else       simulation=kFALSE;

    // Create category for hit data
    pHitCat=gHades->getCurrentEvent()->getCategory(catRpcHit);

    if (!pHitCat) {

	if (simulation) pHitCat=rpc->buildMatrixCategory("HRpcHitSim",0.5);
	else            pHitCat=rpc->buildMatrixCategory("HRpcHit",0.5);

	gHades->getCurrentEvent()->addCategory(catRpcHit,pHitCat,"Rpc");
    }

    if(simulation) {
	pGeantRpcCat = gHades->getCurrentEvent()->getCategory(catRpcGeantRaw);
	if (!pGeantRpcCat) {
	    Error("HRpcClusterF::init()","No HGeant RPC Category");
	    return kFALSE;
	}
    }

    iter=(HIterator *)pCalCat->MakeIterator();
    loc.set(3,0,0,0);

    return kTRUE;
}

Int_t HRpcHitF::execute(void)
{
    Float_t tof, charge, xCell, yCell, xSec, ySec, zSec, xMod, yMod, zMod, xLab, yLab, zLab;//[ns], [pC],[mm]
    Float_t startTime; //[ns]
    Float_t rad2deg = 180./TMath::Pi();
    Float_t theta, phi; //[degrees]
    Float_t sigma_X, sigma_Y, sigma_Z, sigma_T;
    Float_t vprop, DPlanes; //[mm/ns], [mm], [mm]
    Int_t trackL[10], trackR[10], trackLDgtr[10], trackRDgtr[10];
    Bool_t isLMotherAtBox[10], isRMotherAtBox[10];

    HRpcCalSim *pCal    = NULL; //Only reading, use the same class for HRpcCalSim and HRpcCal objects
    HRpcHitSim *pHitSim = NULL; //Reading and STORING in category, use 2 different classes.
    HRpcHit    *pHit    = NULL; //Reading and STORING in category, use 2 different classes.
    HStart2Hit *pStartH = NULL; //Start Hit object

    HGeomVector  rRpcInMod;
    HGeomVector2 rRpcInLab,rRpcInSec;
    //HGeomVector rRpcInLab;

    Float_t fQcut = pHitFPar->getQcut();

    //Getting start time. Only for real data!
    if(!simulation)
    {
	if (pStartHitCat) {
	    if((pStartH=(HStart2Hit *)pStartHitCat->getObject(0))!=NULL) {
		startTime =pStartH->getTime();
	    } else {
		startTime = 0.0;
	    }
	} else {
	    startTime = 0.0;
	}
    } else {
	startTime=0.0;
    }

    iter->Reset();

    while ((pCal = (HRpcCalSim *)iter->Next())!= NULL) {

	loc[0] = pCal->getSector();
	loc[1] = pCal->getColumn();
	loc[2] = pCal->getCell();
		
	//Skip deactivated cells
        if (pCellStatusPar->getCellStatus(loc[0],loc[1],loc[2])!=1) continue;

	//Getting cell geometry params for cell dimensions and operations inside the module.

	DPlanes = pGeomCellPar->getDPlanes();
	vprop   = pDigiPar->getVprop();

	//Find object RpcHit

	if (simulation) {
	    pHitSim = (HRpcHitSim*)pHitCat->getSlot(loc);
	    if ( !pHitSim ) {
		Error ("execute()", "Can't get slot in hit finder: sec %i, col %i, cell %i",loc[0],loc[1],loc[2]);
		return EXIT_FAILURE;
	    }
	    pHitSim = new(pHitSim) HRpcHitSim;
	}
	else {
	    pHit = (HRpcHit*)pHitCat->getSlot(loc);
	    if ( !pHit ) {
		Error ("execute()", "Can't get slot in hit finder: sec %i, col %i, cell %i",loc[0],loc[1],loc[2]);
		return EXIT_FAILURE;
	    }
	    pHit = new(pHit) HRpcHit;
	}

	//Calculate the variables for filling the object RpcHit

	charge = 0.5*( pCal->getRightCharge() + pCal->getLeftCharge() );

	//conversion from tdc width to fC
	if(!simulation) charge = pWtoQPar->getPar0() + pWtoQPar->getPar1()*charge + pWtoQPar->getPar2()*charge*charge +
	    pWtoQPar->getPar3()*pow(charge,3) + pWtoQPar->getPar4()*pow(charge,4) + pWtoQPar->getPar5()*pow(charge,5);

	tof    = 0.5*( pCal->getRightTime()   + pCal->getLeftTime()   ) - pGeomCellPar->getLength(loc[1],loc[2])/(2*vprop);
	xCell  = 0.5*( pCal->getLeftTime()    - pCal->getRightTime()  )*vprop;

	if(!simulation) {
	    xCell  = pCal->getRightTime()*vprop;  // FIXME in analysis: temporary!!!!!!
	    tof    = pCal->getLeftTime() - pGeomCellPar->getLength(loc[1],loc[2])/(2*vprop);  // FIXME in analysis: temporary!!!!!!
	    tof   -= startTime;
	}
	yCell  = 0.5*( pGeomCellPar->getWidth(loc[1],loc[2]) );

	//cout << "Tof antes: " << tof << endl;
	//Applying Slewing Correction
	if(!simulation) {
	    if(charge<fQcut) {
		const Float_t *params = pHitFPar->getPol5Pars(loc[0],loc[1],loc[2]);
		tof = tof - (params[0] + params[1]*charge + params[2]*pow(charge,2) + params[3]*pow(charge,3) + params[4]*pow(charge,4) + params[5]*pow(charge,5));
	    } else {
		const Float_t *params = pHitFPar->getPol1Pars(loc[0],loc[1],loc[2]);
		tof = tof - (params[0] + params[1]*charge);
	    }
	}

	if(!simulation)xMod   = pGeomCellPar->getX(loc[1],loc[2]) - pGeomCellPar->getLength(loc[1],loc[2])/2 - xCell;
	else           xMod   = pGeomCellPar->getX(loc[1],loc[2]) + pGeomCellPar->getLength(loc[1],loc[2])/2 + xCell;  // mod
	yMod   = pGeomCellPar->getY(loc[1],loc[2]) + yCell;
	zMod   = DPlanes*(loc[1]%2-0.5);

	rRpcInMod.setX(xMod);
	rRpcInMod.setY(yMod);
	rRpcInMod.setZ(zMod);

	//Flag for cell outliers
	Bool_t isInsideCell = true;
	Float_t extra    = pGeomCellPar->getWidth(loc[1],loc[2]) * TMath::Tan(24.5/rad2deg);
	Float_t minRange = pGeomCellPar->getX(loc[1],loc[2])-pGeomCellPar->getLength(loc[1],loc[2])-extra;
	Float_t maxRange = pGeomCellPar->getX(loc[1],loc[2])+extra;
	if(xMod<minRange || xMod>maxRange) isInsideCell = false;


	//Lab system: all sectors in Lab
	HGeomTransform& TransRpc2Lab = pRpcGeometry->getModule(loc[0],0)->getLabTransform();
	rRpcInLab = TransRpc2Lab.transFrom(rRpcInMod);

	xLab = rRpcInLab.getX();
	yLab = rRpcInLab.getY();
	zLab = rRpcInLab.getZ();

	//Sec system: all sectors referred as sector 0
	HGeomTransform &TransRpc2Sec = pSpecGeomPar->getSector(loc[0])->getTransform();
	rRpcInSec = TransRpc2Sec.transTo(rRpcInLab);

	xSec = rRpcInSec.getX();
	ySec = rRpcInSec.getY();
	zSec = rRpcInSec.getZ();

	theta = (zSec>0.) ? (rad2deg * TMath::ATan(ySec/zSec)) : 0.;
	phi = rad2deg * TMath::ATan2(yLab,xLab);
	if (phi < 0.) phi += 360.;

	//Obtain the resolutions from both the digitizer info and the geometrical one
	sigma_X = pDigiPar->getSigmaX();
	sigma_T = pDigiPar->getSigmaT()/1000;
	sigma_Z = (pGeomCellPar->getDeltaZ())/sqrt(12.);
	sigma_Y = (pGeomCellPar->getWidth(loc[1],loc[2]))/sqrt(12.);

	//Fill the Hit

	if(simulation) {
	    pHitSim->setAddress(pCal->getAddress());
	    pHitSim->setHit(tof,charge,xMod,yMod,zMod);
	    pHitSim->setXSec(xSec);
	    pHitSim->setYSec(ySec);
	    pHitSim->setZSec(zSec);
	    pHitSim->setXYZLab(xLab,yLab,zLab);
	    pHitSim->setRMS(sigma_T, sigma_X, sigma_Y, sigma_Z);

	    pHitSim->setRefL(pCal->getRefL());
	    pHitSim->setRefR(pCal->getRefR());
	    pHitSim->setRefLDgtr(pCal->getRefLDgtr());
	    pHitSim->setRefRDgtr(pCal->getRefRDgtr());

	    pCal->getLisAtBoxArray(isLMotherAtBox);
	    pCal->getRisAtBoxArray(isRMotherAtBox);
	    pCal->getTrackLArray(trackL);
	    pCal->getTrackRArray(trackR);
	    pCal->getTrackLDgtrArray(trackLDgtr);
	    pCal->getTrackRDgtrArray(trackRDgtr);

	    pHitSim->setLisAtBoxArray(isLMotherAtBox);
	    pHitSim->setRisAtBoxArray(isRMotherAtBox);

	    pHitSim->setTrackLArray(trackL);
	    pHitSim->setTrackRArray(trackR);
	    pHitSim->setTrackLDgtrArray(trackLDgtr);
	    pHitSim->setTrackRDgtrArray(trackRDgtr);

	}
	else {
	    pHit->setAddress(pCal->getAddress());
	    pHit->setLogBit(pCal->getLogBit());
	    pHit->setHit(tof,charge,xMod,yMod,zMod);
	    pHit->setXSec(xSec);
	    pHit->setYSec(ySec);
	    pHit->setZSec(zSec);
	    pHit->setXYZLab(xLab,yLab,zLab);
	    pHit->setRMS(sigma_T, sigma_X, sigma_Y, sigma_Z);
	    pHit->setTheta(theta);
	    pHit->setPhi(phi);
	    pHit->setInsideCellFlag(isInsideCell);
	}
    }

    return EXIT_SUCCESS;
}

ClassImp(HRpcHitF)