#include "hrpcdigitizer.h"
#include "hades.h"
#include "hruntimedb.h"
#include "hdebug.h"
#include "hspectrometer.h"
#include "hrpcdetector.h"
#include "hevent.h"
#include "hcategory.h"
#include "hiterator.h"
#include "hlocation.h"
#include "rpcdef.h"
#include "hgeantrpc.h"
#include "hgeantkine.h"
#include "hrpccalsim.h"
#include "hrpccal.h"
#include "hrpcdigipar.h"
#include "hrpcgeomcellpar.h"

#include "TRandom.h"

#include <iostream>
#include <iomanip>
using namespace std;

//*-- Author   : D. Gonzalez-Diaz 08/06/2006
//*-- Modified : D. Gonzalez-Diaz 12/06/2006
//*-- Modified : D. Gonzalez-Diaz 13/12/2009
//*-- Modified : P. Cabanelas     24/11/2010 (DigiPar container moved to condition style)
/////////////////////////////////////////////////////////////////////
//
//  HRpcDigitizer digitizes HGeantRpc data, puts output values into
//  HRpcCalSim data category.
//
//
/////////////////////////////////////////////////////////////////////

HRpcDigitizer::HRpcDigitizer(void) {
    initVars();
}

HRpcDigitizer::HRpcDigitizer(const Text_t *name, const Text_t *title) :
HReconstructor(name,title) {
    initVars();
}

HRpcDigitizer::~HRpcDigitizer(void) {

    if(iterGeantRpc) delete iterGeantRpc;
    if(iterRpcCal)   delete iterRpcCal;
    clearObjects();
    rpcobjects.clear();
}
void HRpcDigitizer::initVars(){

    fGeantRpcCat = 0;
    fCalCat      = 0;
    fKineCat     = 0;
    fGeomCellPar = 0;
    fDigiPar     = 0;
    fLoc.set(3,0,0,0);
    iterGeantRpc = 0;
    iterRpcCal   = 0;
}

Bool_t HRpcDigitizer::initParContainer() {
    fGeomCellPar=(HRpcGeomCellPar *)gHades->getRuntimeDb()->getContainer("RpcGeomCellPar");
    if (!fGeomCellPar){
	Error("init","No RpcGeomCellPar Parameters");
	return kFALSE;
    }
    fDigiPar    =(HRpcDigiPar *)    gHades->getRuntimeDb()->getContainer("RpcDigiPar");
    if (!fDigiPar){
	Error("init","No RpcDigiPar Parameters");
	return kFALSE;
    }
    return kTRUE;
}

Bool_t HRpcDigitizer::init(void) {

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


    if(!rpc){
	Error("init","No Rpc Detector found");
	return kFALSE;
    }
    //----------------------------------------------
    // working arrays
    maxCol  = rpc->getMaxColumns();
    maxCell = rpc->getMaxCells();
    rpcobjects.resize(6*maxCol*maxCell, 0 );  // size is constant over run time
    //----------------------------------------------

    //----------------------------------------------
    // parameters
    if(!initParContainer()) return kFALSE;
    //----------------------------------------------

    //----------------------------------------------
    // data
    fKineCat = gHades->getCurrentEvent()->getCategory(catGeantKine);
    if (!fKineCat) {
	Error("HRpcDigitizer::init()","HGeant kine input missing");
	return kFALSE;
    }

    fGeantRpcCat = gHades->getCurrentEvent()->getCategory(catRpcGeantRaw);
    if (!fGeantRpcCat) {
	Error("HRpcDigitizer::init()","HGeant RPC input missing");
	return kFALSE;
    }

    //Build the Calibration category

    fCalCat = gHades->getCurrentEvent()->getCategory(catRpcCal);
    if (!fCalCat) {
	fCalCat=rpc->buildMatrixCategory("HRpcCalSim",0.5);
	gHades->getCurrentEvent()->addCategory(catRpcCal,fCalCat,"Rpc");
    }

    iterGeantRpc = (HIterator *)fGeantRpcCat->MakeIterator("native");
    iterRpcCal   = (HIterator *)fCalCat     ->MakeIterator("native");
    //----------------------------------------------

    return kTRUE;
}

Int_t HRpcDigitizer::execute(void) {

    HGeantRpc*  geantrpc = 0;
    HRpcCalSim* cal      = 0;

    Int_t   RefL=-1;
    Float_t geaTof = 0., geaEner = 0.,   geaX = 0.,   geaXLoc = 0., geaY = 0.,
	      geaZ = 0.,  geaMom = 0., geaLen = 0., geaLocLen = 0., 
                x  = 0.,       y = 0.,      z = 0.,         p = 0.,   dE = 0., t = 0.;
    Float_t sigma_T, sigma_el, vprop, Pgap, t_offset, T_smearing, D, gap, Qmean, Qgap = 0.;


    //----------------------------------------------
    // loop over the HGeantrpc objects, i.e, over the gaps
    // and fill temporary working objects
    clearObjects();
    iterGeantRpc->Reset();
    while ((geantrpc=(HGeantRpc *)iterGeantRpc->Next()) != 0) {

	if(geantrpc->getDetectorID() < 0) continue; //Hit in virtual volume EBOX

	geantrpc->getHit(geaX, geaY, geaZ, geaTof, geaMom, geaEner); //[mm], [ns], [MeV] in module ref system
	geantrpc->getTLength(geaLen, geaLocLen);

	fLoc[0] = geantrpc->getSector();
	fLoc[1] = geantrpc->getColumn();
	fLoc[2] = geantrpc->getCell();

	//Load the relevant parameters

	//Note that fGeomCellPar->getX(column,cell) corresponds to the left-down corner

	vprop    = fDigiPar->getVprop();                         //[mm/ns]
	sigma_el = (fDigiPar->getSigmaX()) * 2 / vprop / sqrt(2.);//[ns]
	sigma_T  = (fDigiPar->getSigmaT()) / 1000;                //[ns]
	t_offset = (fDigiPar->getToff()) / 1000;                  //[ns]
	Qmean    = (fDigiPar->getQmean());			  //[pC]
	gap      = fDigiPar->getGap();				  //[mm]

	Float_t ineff_gap;
	ineff_gap = sqrt(sqrt(1 - fDigiPar->getEff()));
	ineff_gap = pow(ineff_gap,(geaLocLen/gap));             //correction for angle of incidence

	Pgap     = 1 - ineff_gap;        //Efficiency per gap


	//This may depend a bit on the digitization algorithm
	if(gRandom->Uniform(1) > Pgap) continue;

        //-------------------------------------------------------------
        // collect data
	Int_t ind = fLoc[0] * maxCol * maxCell + fLoc[1] * maxCell + fLoc[2]; // objects are stored linear

	if(rpcobjects[ind] == 0) { // first time we have to create the object
	    rpcobjects[ind] = new rpcdat;
	    rpcobjects[ind]->linIndex= ind;
	    rpcobjects[ind]->sec     = fLoc[0];
	    rpcobjects[ind]->col     = fLoc[1];
	    rpcobjects[ind]->cell    = fLoc[2];
	}
        rpcdat* dat = rpcobjects[ind];
        //------------------------------------------------
        // create new objects for left/right of the gap
        // and add them to the list of gaps for this cell
	gaptrack* left  = new gaptrack;
        gaptrack* right = new gaptrack;

	left ->reset();
	right->reset();

        dat->left .push_back(left);
        dat->right.push_back(right);
        //------------------------------------------------

        left ->refDgtr   = fGeantRpcCat->getIndex(geantrpc);
        left ->trackDgtr = ((HGeantRpc*)fGeantRpcCat->getObject(left->refDgtr))->getTrack();

        // right is the same as left
	right->refDgtr   = left->refDgtr;
	right->trackDgtr = left->trackDgtr;


	RefL  = findMother(left->refDgtr); // right is the same as left , do it only once

	if(RefL < 0){ //Mother not found. Very seldom. Keep info of daughter.
	    left ->ref     = left ->refDgtr;
            left ->track   = left ->trackDgtr;
            left ->isAtBox = kFALSE;
	    right->ref     = left ->refDgtr;
            right->track   = left ->trackDgtr;
	    right->isAtBox = kFALSE;
	  
	    Warning("execute()","mother of track not found in RPC box! will use RefLDgtr=%i, TrackLDgtr=%i",left->refDgtr,left->trackDgtr);
	}
	else {        //Mother found
	    left ->track   = ((HGeantRpc*)fGeantRpcCat->getObject(RefL))->getTrack();
            left ->ref     = RefL;
	    left ->isAtBox = kTRUE;

	    right->track   = left->track;
            right->ref     = left->ref;
	    right->isAtBox = kTRUE;
	}

        // calculate GEANT time left/right without smearing for sorting
	geantrpc->getHit(x, y, z, t, p, dE);
	right->gtime = t + (fGeomCellPar->getLength(fLoc[1],fLoc[2]) - (x - fGeomCellPar->getX(fLoc[1],fLoc[2]))) / fDigiPar->getVprop(); //Geant Time at right-end
	left ->gtime = t + (x - fGeomCellPar->getX(fLoc[1],fLoc[2])) / fDigiPar->getVprop();                      //Geant Time at left-end


	//Change to cell coordinate system (with origin at the cell center)
	D       =  fGeomCellPar->getLength(fLoc[1],fLoc[2]);                   // Cell length in mm
	geaXLoc =  geaX - fGeomCellPar->getX(fLoc[1],fLoc[2]);

	T_smearing = gRandom->Gaus(t_offset, sigma_T); //Common to left and right!. Intrinsic RPC fluctuations

	left ->time = geaTof + T_smearing + geaXLoc / vprop + gRandom->Gaus(0,sigma_el);
	right->time = geaTof + T_smearing + (D - geaXLoc) / vprop + gRandom->Gaus(0,sigma_el);

	//Digitizing model I. Take the fastest in the left and in the right. Charge distribution
	//assumed to be flat per gap. At present slewing correction is not included. Also charges in
	//left and in right are *exactly* the same by definition. A smearing in charge to account
	//for the resolution of the Q2W method can be easily included.


	Qgap = (gRandom->Uniform(1)) * ( 2 * Qmean) / 4 * (geaLocLen / gap); //correction for angle of incidence;
        left ->charge = Qgap;
	right->charge = Qgap;
	//-------------------------------------------------------------


	//end of algorithm

    } // end of HGeantRpc loop

    //------------------------------------------------------------
    // fill output

    Int_t  Track    [10];
    Int_t  TrackDgtr[10];
    Bool_t isAtBox  [10];
    vector <Int_t> tracklist;
    tracklist.reserve(20); // minimum size

    for(UInt_t i = 0; i < rpcobjects.size(); i ++) {
	rpcdat* dat = rpcobjects[i];
	if(dat){
	    fLoc[0] = dat->sec;
            fLoc[1] = dat->col;
            fLoc[2] = dat->cell;

            //-----------------------------------------------------------
            // test if cell in use. Needed in embedding mode
	    cal = (HRpcCalSim*) fCalCat->getObject(fLoc);   
	    if(cal) {
                // add track from real data for sorting procedure
		gaptrack* left  = new gaptrack;
		gaptrack* right = new gaptrack;

		left ->reset();
		right->reset();

		dat->left .push_back(left);
		dat->right.push_back(right);

		left ->time      = cal->getLeftTime();
		left ->gtime     = cal->getLeftTime();
                left ->charge    = cal->getLeftCharge();
                left ->track     = gHades->getEmbeddingRealTrackId();
		left ->trackDgtr = gHades->getEmbeddingRealTrackId();
                left ->isAtBox   = kTRUE;
		right->time      = cal->getRightTime();
		right->gtime     = cal->getRightTime();
                right->charge    = cal->getRightCharge();
                right->track     = gHades->getEmbeddingRealTrackId();
		right->trackDgtr = gHades->getEmbeddingRealTrackId();
                right->isAtBox   = kTRUE;
	    } else {
		cal = (HRpcCalSim*) fCalCat->getSlot(fLoc);
		if(!cal) Error("execute()","Error: could not allocate a new slot in HRpcCalSim!");
		cal = new(cal) HRpcCalSim;
		cal->setAddress(fLoc[0], fLoc[1], fLoc[2]);
	    }
            //-----------------------------------------------------------
            // sort by GEANT time and select the fastest hit
	    dat->sortTime(kTRUE);  // right
	    dat->sortTime(kFALSE); // left
	    cal->setRightTime  (dat->getSmallestTof(kTRUE )); // right
	    cal->setLeftTime   (dat->getSmallestTof(kFALSE)); // left
	    cal->setRightCharge(dat->getSumCharge  (kTRUE )); // right
	    cal->setLeftCharge (dat->getSumCharge  (kFALSE)); // left

	    cal->setRefR    (dat->right[0]->ref);
	    cal->setRefRDgtr(dat->right[0]->refDgtr);
            cal->setRefL    (dat->left [0]->ref);
	    cal->setRefLDgtr(dat->left [0]->refDgtr);
            //-----------------------------------------------------------


	    //------------------------------------------------------------
	    // fill list of tracks
            // remove doubles

	    //------------------------------------------------------------
            // right side
	    tracklist.clear();
	    Int_t ct = 0;

	    for(UInt_t j = 0; j < 10; j ++) {
		Track[j] = TrackDgtr[j] = -999;
		isAtBox[j] = kFALSE;
	    }

	    for(UInt_t j = 0; j < dat->right.size() && ct < 10; j ++ ){
		Int_t tr = dat->right[j]->track;
		if(find(tracklist.begin(),tracklist.end(),tr ) == tracklist.end()){ // track was not stored before
		    Track    [ct] = tr;
                    TrackDgtr[ct] = dat->right[j]->trackDgtr;
                    isAtBox  [ct] = dat->right[j]->isAtBox;
		    tracklist.push_back(tr);  // remember used track
                    ct ++;
		}
	    }

	    cal->setTrackRArray    (Track);
	    cal->setTrackRDgtrArray(TrackDgtr);
	    cal->setRisAtBoxArray  (isAtBox);
            cal->setNTracksR(ct);
	    //------------------------------------------------------------

	    //------------------------------------------------------------
	    // left side
	    tracklist.clear();
	    ct = 0;

	    for(UInt_t j = 0; j < 10; j ++) {
		Track[j] = TrackDgtr[j] = -999;
		isAtBox[j] = kFALSE;
	    }

	    for(UInt_t j = 0; j < dat->left.size() && ct < 10 ; j ++ ){
                Int_t tr = dat->left[j]->track;
		if(find(tracklist.begin(),tracklist.end(),tr) == tracklist.end()){ // track was not stored before
		    Track    [ct] = tr;
                    TrackDgtr[ct] = dat->left[j]->trackDgtr;
                    isAtBox  [ct] = dat->left[j]->isAtBox;
		    tracklist.push_back(tr); // remember used track
                    ct ++;
		}
	    }

	    cal->setTrackLArray    (Track);
	    cal->setTrackLDgtrArray(TrackDgtr);
	    cal->setLisAtBoxArray  (isAtBox);
            cal->setNTracksL(ct);
	    //------------------------------------------------------------

	    //------------------------------------------------------------
	}
    }
    //------------------------------------------------------------

    return 0;
}

void  HRpcDigitizer::clearObjects(){
    // delete objects in working array and
    // set pointer to 0. vector is still
    // not cleared.
    for(UInt_t i = 0; i < rpcobjects.size(); i ++){
	if(rpcobjects[i]){
	    rpcobjects[i]->reset();
	    delete rpcobjects[i];
	    rpcobjects[i] = 0;
	}
    }
}

Int_t HRpcDigitizer::findMother(Int_t Ref_initial) {

    // Algorithm for finding mother. When not found, -1 is returned
    // in order to tag the track. It very rarely happens and is
    // always associated to neutral particles at the box.

    Int_t detectorID =  0, track_mother = -1, Ref_final = -1;

    HGeantRpc*  geantrpc = 0;
    HGeantKine* kine     = 0;

    track_mother = ((HGeantRpc*)fGeantRpcCat->getObject(Ref_initial))->getTrack();
    kine         = (HGeantKine*)fKineCat    ->getObject(track_mother - 1);

    while(kine && detectorID >= 0){

	if(track_mother == 0)  { //End of loop. Mother not found
	    Ref_final = -999;
	    break;                     
	}

	kine = (HGeantKine*)fKineCat->getObject(track_mother - 1);

	if(kine && kine->getNRpcHits() == 0) { //Mother has no hit in RPC cells
	    track_mother = kine->getParentTrack(); 
	    continue;
	}

	Ref_final    = kine->getFirstRpcHit(); //Look for the first hit (time-wise) of the track
	geantrpc     = (HGeantRpc*)fGeantRpcCat->getObject(Ref_final);
	detectorID   = geantrpc->getDetectorID();
	track_mother = kine->getParentTrack();

    }

    return Ref_final;
}


ClassImp(HRpcDigitizer)