#include "hrpcdigitizer.h"
#include "hades.h"
#include "hruntimedb.h"
#include "hphysicsconstants.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   : 08/06/2006 D. Gonzalez-Diaz
//*-- Modified : 12/06/2006 D. Gonzalez-Diaz
//*-- Modified : 13/12/2009 D. Gonzalez-Diaz
//*-- Modified : 24/11/2010 P. Cabanelas DigiPar container moved to
//*-- Modified :                         condition style.
//*-- Modified : 22/06/2012 A.Mangiarotti General code cleanup.
//*-- Modified :            Geometry parameters used in several places
//*-- Modified :            are now stored as local variables (see D,
//*-- Modified :            geaXLoc). Reading of the digitiser
//*-- Modified :            parameters, which are constant, moved out
//*-- Modified :            of the hit processing loop (see vprop,
//*-- Modified :            sigma_el, sigma_T, t_offset, Qmean, gap)!
//*-- Modified :            The inefficiency per gap for normally
//*-- Modified :            incident particles ineff_gap_n is
//*-- Modified :            calculated once for all outside the hit
//*-- Modified :            loop. Namespace TMath used consistently
//*-- Modified :            everywhere. The procedure has been
//*-- Modified :            simplified. First the ideal left and
//*-- Modified :            right times left->gtime and right->gtime
//*-- Modified :            are calculated and after the resolutions
//*-- Modified :            are added to obtain left ->time and
//*-- Modified :            right->time.
//*-- Modified : 20/09/2012 A.Mangiarotti Special treatment of
//*-- Modified :            particles produced in the gap.
//*-- Modified : 20/11/2012 A.Mangiarotti modified to include the new
//*-- Modified :            model II with parametrised dependences on
//*-- Modified :            the beta of the particles. Efficiency and
//*-- Modified :            time of flight resolution employ a
//*-- Modified :            sigmoid function. The charge distribution
//*-- Modified :            is now a Landau with parameters depending
//*-- Modified :            on beta with second degree polynomials.
//*-- Modified :            An attempt to correct on average for the
//*-- Modified :            number of active gaps is included. The
//*-- Modified :            old digitizer is now available as model
//*-- Modified :            I. If the old parameter container is
//*-- Modified :            found, model I is used automatically.
//*-- Modified : 27/11/2012 A.Mangiarotti new cell wise mode of
//*-- Modified :            operation controllable with a flag
//*-- Modified :            (fMode=1). If the new mode is selected,
//*-- Modified :            the gaps are first reorganised into
//*-- Modified :            cells before applying the response with
//*-- Modified :            model I or model II. Two new internal
//*-- Modified :            structs celltrack and celldat added to
//*-- Modified :            handle the reorganisation of the gaps
//*-- Modified :            into cells. If model II is used in cell
//*-- Modified :            wise mode, no further corrections are
//*-- Modified :            needed for the average number of active
//*-- Modified :            gaps.
/////////////////////////////////////////////////////////////////////
//
//  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();
    cellobjects.clear();
    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("HRpcDigitizer::init()","No RpcGeomCellPar Parameters");
        return kFALSE;
    }
    fDigiPar    =(HRpcDigiPar *)    gHades->getRuntimeDb()->getContainer("RpcDigiPar");
    if (!fDigiPar){
        Error("HRpcDigitizer::init()","No RpcDigiPar Parameters");
        return kFALSE;
    }
    return kTRUE;
}

Bool_t HRpcDigitizer::init(void) {

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


    if(!rpc){
        Error("HRpcDigitizer::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
    cellobjects.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;

    Float_t geaTof = 0., geaEner = 0., geaX = 0., geaY = 0., geaZ = 0., geaXLoc = 0.;
    Float_t geaMom = 0., geaLen = 0., geaLocLen = 0.;
    Float_t sigma0_T, sigma1_T, sigma2_T, sigma3_T, sigma_el, vprop, Pgap, t_offset;
    Float_t T_smearing, D, gap, Qmean0, Qmean1, Qmean2, Qwid0, Qwid1, Qwid2, Qgap;
    Float_t Eff0, Eff1, Eff2, Eff3;
    Float_t ineff_gap_n, ineff_gap, mass, beta, beta2;
    Float_t corr_ngap,sigma_T,X_smearing,Qmean,Qwid;
    Int_t   RefL=-1,trackNumber,gpid,mode,ngap;


    //Load the relevant parameters

    vprop       = fDigiPar->getVprop();                                  //[mm/ns]
    sigma_el    = (fDigiPar->getSigmaX()) / vprop;                       //[ns]
    sigma0_T    = (fDigiPar->getSigmaT())  / 1000.;                      //[ns]
    sigma1_T    = (fDigiPar->getSigmaT1()) / 1000.;                      //[ns]
    sigma2_T    = (fDigiPar->getSigmaT2());                              //pure number
    sigma3_T    = (fDigiPar->getSigmaT3());                              //pure number
    t_offset    = (fDigiPar->getToff()) / 1000.;                         //[ns]
    Qmean0      = (fDigiPar->getQmean());                                //[pC]
    Qmean1      = (fDigiPar->getQmean1());                               //[pC]
    Qmean2      = (fDigiPar->getQmean2());                               //[pC]
    Qwid0       = (fDigiPar->getQwid());                                 //[pC]
    Qwid1       = (fDigiPar->getQwid1());                                //[pC]
    Qwid2       = (fDigiPar->getQwid2());                                //[pC]
    Eff0        = (fDigiPar->getEff());                                  //pure number
    Eff1        = (fDigiPar->getEff1());                                 //pure number
    Eff2        = (fDigiPar->getEff2());                                 //pure number
    Eff3        = (fDigiPar->getEff3());                                 //pure number
    gap         = fDigiPar->getGap();                                    //[mm]
    mode        = fDigiPar->getMode();                                   //mode flag
    if((Eff1==0.)||(Eff2<=0.)||(Eff3<=0.)) {
      //Digitizing model I.
      ineff_gap_n = TMath::Sqrt(TMath::Sqrt(1. - fDigiPar->getEff()));
    } else {
	ineff_gap_n = 0.;
    }
    if((sigma1_T<=0.)||(sigma2_T<=0.)||(sigma3_T<=0.)) {
      //Digitizing model I.
      sigma_el *= TMath::Sqrt2();
    }


    //---------------------------------------------------------------
    // 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);
        // If the gap enetering was missed in Geant, assume the gap as local track length
        if (geaLocLen<0.) geaLocLen=gap;

        // Determine the beta of the particle.
        trackNumber = geantrpc->getTrack();
        HGeantKine* kine = (HGeantKine*)fKineCat->getObject(trackNumber-1);
        if(!kine) {
          Warning("HRpcDigitizer::execute()",
                  "missing kine entry for track %i, assuming beta=0",trackNumber);
          beta=0.;
        } else {
          gpid =  kine->getID();
          //From the Geant particle type find the mass (only charged
          //particles need to be considered).

	  mass=HPhysicsConstants::mass(gpid);
	  if(mass < 0.){
	      Warning("HRpcDigitizer::execute()","unknown particle id %i, assuming proton mass",gpid);
	      mass=HPhysicsConstants::mass(14);
          }
          beta=mass/geaMom;
          beta=1./TMath::Sqrt(1.+beta*beta);
        }

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


        if (mode==0) {
        /////////////////////////////////////////////////////////////
        /////////////////////////////////////////////////////////////
        //////////// BEGIN GAP WISE OPEREATION (MODE=0) /////////////
        /////////////////////////////////////////////////////////////
        /////////////////////////////////////////////////////////////


          if((Eff1==0.)||(Eff2<=0.)||(Eff3<=0.)) {
            //Digitizing model I.
            ineff_gap = TMath::Power(ineff_gap_n,(geaLocLen/gap));    //correction for angle of incidence
          } else {
            //Digitizing model II.
            ineff_gap = Eff0+Eff1/(1.+TMath::Exp(-Eff2*(beta-Eff3))); //efficeny of the cell
            ineff_gap = TMath::Sqrt(TMath::Sqrt(1. - ineff_gap));     //convert to inefficency of a gap
            ineff_gap = TMath::Power(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
          // objects are stored linear
          Int_t ind = fLoc[0] * maxCol * maxCell + fLoc[1] * maxCell + fLoc[2];

          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("HRpcDigitizer::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;
          }


          //Change to cell coordinate system (with origin at the cell center)
          //Note that fGeomCellPar->getX(column,cell) corresponds to the left-down corner
          D       = fGeomCellPar->getLength(fLoc[1],fLoc[2]);      // Cell length in mm
          geaXLoc = fGeomCellPar->getX(fLoc[1],fLoc[2]) - geaX;    // Position in cell frame

          //Calculate GEANT time left/right without smearing for sorting
          left->gtime  = geaTof + geaXLoc / vprop;         //Geant Time at left-end
          right->gtime = geaTof + (D-geaXLoc) / vprop;     //Geant Time at right-end


          //Common to left and right!. Intrinsic RPC fluctuations
          if((sigma1_T<=0.)||(sigma2_T<=0.)||(sigma3_T<=0.)) {
            //Digitizing model I.
            //Common to left and right!. Intrinsic RPC fluctuations
            T_smearing  = gRandom->Gaus(t_offset, sigma0_T);
            left ->time = left ->gtime + T_smearing + gRandom->Gaus(0.,sigma_el);
            right->time = right->gtime + T_smearing + gRandom->Gaus(0.,sigma_el);
          } else {
            //Digitizing model II.
            sigma_T = sigma0_T+sigma1_T/(1.+TMath::Exp(-sigma2_T*(beta-sigma3_T)));
            //cout << sigma0_T << " " << sigma1_T << " " << sigma2_T << " " << sigma3_T << " "
            //     << beta << " " << sigma_T << endl;
            //Attempt to approximately correct for the number of active gaps
            Float_t norm=0.;
            corr_ngap=0.;
            for(Int_t i = 1; i < 5; i ++) {
              Float_t w4i=TMath::Binomial(4,i)*TMath::Power(Pgap,i)*TMath::Power(ineff_gap,4-i);
              norm += w4i;
              corr_ngap += w4i/TMath::Sqrt((Float_t)i);
            }
            //The correction factor for the sigma of the Gaussian is the
	    //inverse of corr_ngap/norm, to have a multiplicative correction
            //it is calculated here once for all
            corr_ngap = norm/corr_ngap;
            //cout << beta << " " << Eff0+Eff1/(1.+TMath::Exp(-Eff2*(beta-Eff3))) 
            //     << " " << Pgap << " " << corr_ngap << " " << sigma_T << endl;
            T_smearing = gRandom->Gaus(t_offset, sigma_T*corr_ngap);
            //The position and time resolutions are decoupled.
            X_smearing = gRandom->Gaus(0.,sigma_el*corr_ngap);
            //X_smearing = 0.;
            left ->time = left ->gtime + T_smearing + X_smearing;
            right->time = right->gtime + T_smearing - X_smearing;
            //cout << fLoc[0] << " "  << fLoc[1] << " " << fLoc[2] << " " 
            //     << geantrpc->getGap() << " "
            //     << left ->gtime << " " << left ->time << " "
            //     << right->gtime << " " << right->gtime << " " << endl;
          }

          //Generated charge.
          if(Qwid0<=0.) {
            //Digitizing model I.
            //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. * Qmean0) / 4. * (geaLocLen / gap); //correction for angle of incidence
          } else {
            //Digitizing model II.
            //Use Landau distribution with a mean and width parametrised with
            //second order polynomials in beta.
            beta2 = beta*beta;
            Qmean = Qmean0+Qmean1*beta+Qmean2*beta2;
            Qwid  = Qwid0 +Qwid1*beta +Qwid2*beta2;
            //Attempt to approximately correct for the number of active gaps
            Float_t norm=0.;
            corr_ngap=0.;
            for(Int_t i = 1; i < 5; i ++) {
              Float_t w4i=TMath::Binomial(4,i)*TMath::Power(Pgap,i)*TMath::Power(ineff_gap,4-i);
              norm += w4i;
              corr_ngap += w4i*TMath::Sqrt((Float_t)i);
            }
            //The correction factor for the sigma of the Gaussian is the
	    //inverse of corr_ngap/norm, to have a multiplicative correction
            //it is calculated here once for all. Note that this is very
	    //approximate for a Landau distribution.
            corr_ngap = norm/corr_ngap;
            //per gap with correction for angle of incidence
            Qmean = Qmean/(4.*Pgap) * (geaLocLen / gap);
            //per gap with correction for angle of incidence
            Qwid  = Qwid*corr_ngap * (geaLocLen / gap);
            Qgap  = gRandom->Landau(Qmean,Qwid);
          }
          left ->charge = Qgap;
          right->charge = Qgap;

          //---------------------------------------------------------

        /////////////////////////////////////////////////////////////
        /////////////////////////////////////////////////////////////
        ///////////// END GAP WISE OPEREATION (MODE=0) //////////////
        /////////////////////////////////////////////////////////////
        /////////////////////////////////////////////////////////////


        } else {


        /////////////////////////////////////////////////////////////
        /////////////////////////////////////////////////////////////
        //////////// BEGIN CELL WISE OPEREATION (MODE=1) ////////////
        /////////////////////////////////////////////////////////////
        /////////////////////////////////////////////////////////////
          //---------------------------------------------------------
          // collect data
          // objects are stored linear
          Int_t ind = fLoc[0] * maxCol * maxCell + fLoc[1] * maxCell + fLoc[2];

          if(cellobjects[ind] == 0) {
            // first time we have to create the object
            cellobjects[ind] = new celldat;
            //rpcobjects[ind]->linIndex= ind;
            cellobjects[ind]->sec     = fLoc[0];
            cellobjects[ind]->col     = fLoc[1];
            cellobjects[ind]->cell    = fLoc[2];
          }
          celldat* cdat = cellobjects[ind];
          //---------------------------------------------------------
          // create new objects for the track crossing the cell
          // and add them to the list of tracks for this cell
          celltrack* celltr;
          // CASE I: no tacks in this cell yet
          if(cdat->celltr.size()==0) {
            celltr = new celltrack;
            celltr->reset();
            cdat->celltr.push_back(celltr);
          // CASE II: this cell has already tracks associated
          // with it
          } else {
            //Look if current track is already present in the list
            for(UInt_t i=0;i<cdat->celltr.size();i++){
              celltr = cdat->celltr[i];
              //CASE II a: the current track is already associated,
              //so select it
              if(celltr->track == trackNumber) break;
            }
            //CASE II b: the current track is not associated,
            //so add it
            if(celltr->track != trackNumber) {
              celltr = new celltrack;
              celltr->reset();
              cdat->celltr.push_back(celltr);
            }
          }
          //---------------------------------------------------------

          ngap = geantrpc->getGap();
          //The condition on the time takes care of possible multiple
          //rentries of the same track.
          if(celltr->gaptime[ngap]==0.) {
            celltr->gaptime[ngap] = geaTof;
            celltr->gappos[ngap]  = geaX;
            celltr->gapbeta[ngap] = beta;
            celltr->track         = trackNumber;
            celltr->geantrpcIndex = fGeantRpcCat->getIndex(geantrpc);
          } else if(celltr->gaptime[ngap]<geaTof) {
            celltr->gaptime[ngap] = geaTof;
            celltr->gappos[ngap]  = geaX;
            celltr->gapbeta[ngap] = beta;
            celltr->track         = trackNumber;
            celltr->geantrpcIndex = fGeantRpcCat->getIndex(geantrpc);
            Warning("HRpcDigitizer::execute()",
                    "reantry of track=%i in Sector=%i Module=%i Cell=%i",
                    trackNumber,fLoc[0],fLoc[1],fLoc[2]);
          }

        /////////////////////////////////////////////////////////////
        /////////////////////////////////////////////////////////////
        ///////////// END CELL WISE OPEREATION (MODE=1) /////////////
        /////////////////////////////////////////////////////////////
        /////////////////////////////////////////////////////////////

        }

        //end of algorithm

    } // end of HGeantRpc loop


    //---------------------------------------------------------------
    // If in cell mode, process cells

    if(mode>0) {
      for(UInt_t i = 0; i < cellobjects.size(); i ++) {
        celldat* cdat = cellobjects[i];
        if(cdat){
          fLoc[0] = cdat->sec;
          fLoc[1] = cdat->col;
          fLoc[2] = cdat->cell;

          //Change to cell coordinate system (with origin at the cell center)
          //Note that fGeomCellPar->getX(column,cell) corresponds to the left-down corner
          D       = fGeomCellPar->getLength(fLoc[1],fLoc[2]);           // Cell length in mm

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

          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];


          for(UInt_t j = 0; j < cdat->celltr.size(); j ++) {
            celltrack* celltr = cdat->celltr[j];
            celltr->calcMeans();

            //-------------------------------------------------------
            // 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   = celltr->geantrpcIndex;
            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("HRpcDigitizer::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;
            }

            //Change to cell coordinate system (with origin at the cell center)
            geaXLoc = fGeomCellPar->getX(fLoc[1],fLoc[2]) - celltr->pos;  // Position in cell frame
            //Calculate GEANT time left/right without smearing for sorting
            geaTof = celltr->time;                                        //Geant Time
            left->gtime  = geaTof + geaXLoc / vprop;                      //Geant Time at left-end
            right->gtime = geaTof + (D-geaXLoc) / vprop;                  //Geant Time at right-end
            beta = celltr->beta;

            //Common to left and right!. Intrinsic RPC fluctuations
            if((sigma1_T<=0.)||(sigma2_T<=0.)||(sigma3_T<=0.)) {
              //Digitizing model I.
              //Common to left and right!. Intrinsic RPC fluctuations
              T_smearing  = gRandom->Gaus(t_offset, sigma0_T);
              left ->time = left ->gtime + T_smearing + gRandom->Gaus(0.,sigma_el);
              right->time = right->gtime + T_smearing + gRandom->Gaus(0.,sigma_el);
            } else {
              //Digitizing model II.
              sigma_T = sigma0_T+sigma1_T/(1.+TMath::Exp(-sigma2_T*(beta-sigma3_T)));
              //cout << sigma0_T << " " << sigma1_T << " " << sigma2_T << " " << sigma3_T << " "
              //     << beta << " " << sigma_T << endl;
              T_smearing = gRandom->Gaus(t_offset, sigma_T);
              //The position and time resolutions are decoupled.
              X_smearing = gRandom->Gaus(0.,sigma_el);
              //X_smearing = 0.;
              left ->time = left ->gtime + T_smearing + X_smearing;
              right->time = right->gtime + T_smearing - X_smearing;
              //cout << fLoc[0] << " "  << fLoc[1] << " " << fLoc[2] << " " 
              //     << left ->gtime << " " << left ->time << " "
              //     << right->gtime << " " << right->gtime << " " << endl;
            }

            //Generated charge.
            if(Qwid0<=0.) {
              //Digitizing model I.
              //Charge distribution assumed to be flat per cell. 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. * Qmean0) / 4. * (geaLocLen / gap); //correction for angle of incidence
            } else {
              //Digitizing model II.
              //Use Landau distribution with a mean and width parametrised with
              //second order polynomials in beta.
              beta2 = beta*beta;
              Qmean = Qmean0+Qmean1*beta+Qmean2*beta2;
              //correction for angle of incidence
              Qmean = Qmean * (geaLocLen / gap);
              Qwid  = Qwid0 +Qwid1*beta +Qwid2*beta2;
              //correction for angle of incidence
              Qwid  = Qwid * (geaLocLen / gap);
              Qgap  = gRandom->Landau(Qmean,Qwid);
            }
            left ->charge = Qgap;
            right->charge = Qgap;

            //-------------------------------------------------------

          }  // end of celltrack loop

        }

      }  // end of cellobjects 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
            if(((mode==0)&&((sigma1_T<=0.)||(sigma2_T<=0.)||(sigma3_T<=0.)))
               ||(mode>0)) {
              //Digitizing model I and II (except gap-wise)
              //Take the fastest in the left and in the right.
              cal->setRightTime(dat->getSmallestTof(kTRUE )); // right
              cal->setLeftTime (dat->getSmallestTof(kFALSE)); // left
            } else {
              //Digitizing model II gap-wise.
              //Take the average in the left and in the right.
              cal->setRightTime  (dat->getMeanTof(kTRUE ));   // right
              cal->setLeftTime   (dat->getMeanTof(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;
                    // remember used track
                    tracklist.push_back(tr);
                    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;
                    // remember used track
                    tracklist.push_back(tr);
                    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;
        }
    }
    for(UInt_t i = 0; i < cellobjects.size(); i ++){
        if(cellobjects[i]){
            cellobjects[i]->reset();
            delete cellobjects[i];
            cellobjects[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)