#include "hparticlemetamatcher.h"
#include "hades.h"
#include "hspectrometer.h"
#include "hruntimedb.h"
#include "hcategory.h"
#include "hcategorymanager.h"
#include "hgeomvolume.h"
#include "hgeomcompositevolume.h"
#include "hgeomtransform.h"

#include "hemcdetector.h"

#include "hrpcgeompar.h"
#include "hrpcgeomcellpar.h"
#include "hemcgeompar.h"
#include "htofgeompar.h"
#include "hspecgeompar.h"

#include "hparticlecand.h"

#include "rpcdef.h"
#include "emcdef.h"
#include "tofdef.h"

#include "TMath.h"


//_HADES_CLASS_DESCRIPTION
/////////////////////////////////////////////////////////////
//
//  HParticleMetamatcher
//
//
// usage :
// Int_t loopDST_task(
//		   TString infileList="/lustre/hebe/hades/dst/mar19/gen2/ag158ag/3200A/083/root/be1908311041501.hld_dst_mar19.root",
//		   TString outfile="test.root",Int_t nEvents=50)
//{
//    HLoop loop(kTRUE);
//
//    //-------------------------------------------------
//    if(1){
//	TString beamtime="mar19";
//
//	Int_t mdcMods[6][4]=
//	{ {1,1,1,1},
//	{1,1,1,1},
//	{1,1,1,1},
//	{1,1,1,1},
//	{1,1,1,1},
//	{1,1,1,1} };
//	TString asciiParFile     = "";
//	TString rootParFile      = "/cvmfs/hades.gsi.de/param/real/mar19/allParam_Mar19_gen2_18102019.root";
//	TString paramSource      = "root"; // root, ascii, oracle
//	TString paramrelease     = "MAR19_dst_gen3";  // now, APR12_gen2_dst APR12_gen5_dst
//	HDst::setupSpectrometer(beamtime,mdcMods,"rich,mdc,tof,rpc,emc,wall,start,tbox");
//	HDst::setupParameterSources(paramSource,asciiParFile,rootParFile,paramrelease);  // now, APR12_gen2_dst
//    }
//    //-------------------------------------------------
//
//    Bool_t ret =kFALSE;
//    if(infileList.Contains(",")){
//	ret = loop.addMultFiles(infileList);      // file1,file2,file3
//    } else{
//	ret = loop.addFiles(infileList); // myroot*.root
//    }
//
//    if(ret == 0) {
//	cout<<"READBACK: ERROR : cannot find inputfiles : "<<infileList.Data()<<endl;
//	return 1;
//    }
//
//    if(!loop.setInput("")) {   // all input categories
//	cout<<"READBACK: ERROR : cannot read input !"<<endl;
//	exit(1);
//    } // read all categories
//    loop.printCategories();
//    loop.printChain();
//
//    HCategory* candCat    = (HCategory*)HCategoryManager::getCategory(catParticleCand);
//    HCategory* evtInfoCat = (HCategory*)HCategoryManager::getCategory(catParticleEvtInfo);
//
//
//    //-------------------------------------------------
//    // added for the task
//    HTaskSet *masterTaskSet = gHades->getTaskSet("all");
//    HParticleMetaMatcher* matcher = new HParticleMetaMatcher();
//    matcher->setDebug();  // print infos to screen
//    masterTaskSet->add(matcher);
//    //-------------------------------------------------
//
//    Int_t entries = loop.getEntries();
//    if(nEvents < entries && nEvents >= 0 ) entries = nEvents;
//
//    for (Int_t i = 0; i < entries; i++) {
//	Int_t nbytes =  loop.nextEvent(i);             // get next event. categories will be cleared before
//	if(nbytes <= 0) { cout<<nbytes<<endl; break; } // last event reached
//	if(i%1000 == 0) cout<<"event "<<i<<endl;
//
//	HParticleEvtInfo* evtInfo=0;
//	evtInfo = HCategoryManager::getObject(evtInfo,evtInfoCat,0 );
//
//	if(evtInfo&&!evtInfo->isGoodEvent(Particle::kGoodTRIGGER|          // standard event selection apr12
//					  Particle::kGoodVertexClust|
//					  Particle::kGoodVertexCand|
//					  Particle::kGoodSTART|
//					  Particle::kNoPileUpSTART|
//					  Particle::kNoVETO|
//					  Particle::kGoodSTARTVETO|
//					  Particle::kGoodSTARTMETA
//					 )) continue;
//
//	if(candCat){
//
//	    Int_t size = candCat->getEntries();
//	    HParticleCand* cand=0;
//
//	    for(Int_t j = 0; j < size; j++) {
//		cand = HCategoryManager::getObject(cand,candCat,j);
//		if(cand) {
//		    if(!cand->isFlagBit(kIsUsed)) { continue; }           // skip particles rejected by HParticleTrackSorter
//
//		    if(cand->isRpcClstUsed()){
//		        const HRpcCluster* rpc = matcher->recalcRpc(cand); // reconstruct the rpchit from tracking +rk dx dy
//
//			Int_t s,col0,col1,cell0,cell1;
//                      HGeomVector hit0,hit1;
//			matcher->predictRpcCell(cand,hit0,hit1,s,col0,cell0,col1,cell1); // which rpc cell should this track fire?
//		    }
//		    if(0&&(cand->isTofClstUsed() || cand->isTofHitUsed()) ){
//			const HTofCluster* tof = matcher->recalcTof(cand);
//			Int_t s,mod0,mod1,cell0,cell1;
//                      HGeomVector hit0,hit1;
//			matcher->predictTofCell(cand,hit0,hit1,s,mod0,cell0,mod1,cell1); // which Tof cell should this track fire?
//		    }
//		    if(0&&cand->getEmcInd()>-1){
//       	        const HEmcCluster* emc = matcher->recalcEmc(cand);
//			Int_t s,pos,cell;
//			HGeomVector hit;
//                      vector<Int_t> vcells;
//			matcher->predictEmcCell(cand,hit,s,pos,cell); // which emc cell should this track fire?
//		    }
//		}
//	    } // end cand loop
//	} // end cand cat
//    } // end eventloop
//
//    //-------------------------------------------------
//    // added for the task
//    if(gHades)gHades->finalizeTasks();
//    //-------------------------------------------------
//
//    delete gHades;
//    return 0;
//}
//
/////////////////////////////////////////////////////////////

ClassImp(HParticleMetaMatcher)

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

    fCatRpcCluster  = NULL;
    fCatEmcCluster  = NULL;
    fCatTofCluster  = NULL;
    fCatTofHit      = NULL;

    fRpcGeometry    = NULL;
    fEmcGeometry    = NULL;
    fTofGeometry    = NULL;
    pSpecGeomPar    = NULL;

    DPlanesRpc = 0;
    fDebug = kFALSE;
}


const HGeomVector* HParticleMetaMatcher::getRpcCellGeom(Int_t s,Int_t col,Int_t cell)
{
    // returns null pointer if s,col,cell are invalid
    // returns HGeomVector* to first element of array of 4 HGeomVectors
    // for the shape of the cell in mod sys
    if(s < 0 || s > 5) return 0;
    if(col < 0 || col > 5) return 0;
    if(cell < 0 || cell > RPCMAXCELL-1) return 0;
    return &cellRPC[s][col][cell][0];
}

const HGeomVector* HParticleMetaMatcher::getTofCellGeom(Int_t s,Int_t mod,Int_t cell)
{
    // returns null pointer if s,mod,cell are invalid
    // returns HGeomVector* to first element of array of 4 HGeomVectors
    // for the shape of the cell in mod sys
    if(s < 0 || s > 5) return 0;
    if(mod < 0 || mod > 8) return 0;
    if(cell < 0 || cell > 8) return 0;
    return &cellRPC[s][mod][cell][0];
}

const HGeomVector* HParticleMetaMatcher::getEmcCellGeom(Int_t s,Int_t cell)
{
    // returns null pointer if s,cell are invalid
    // returns HGeomVector* to first element of array of 4 HGeomVectors
    // for the shape of the cell in mod sys
    if(s < 0 || s > 5) return 0;
    if(cell < 0 || cell > EMCMAXCELL-1) return 0;
    return &cellEMC[s][cell][0];
}


void HParticleMetaMatcher::calcSegPoints(TVector3& p1,TVector3& p2,Double_t phiseg,Double_t thetaseg,Double_t roseg,Double_t zseg)
{
    // in sector sys : calc p1 p2 points (not MDC Hit!) on MdcSeg

    Double_t Xseg[2],Yseg[2],Zseg[2];

    Xseg[0] = roseg*cos(phiseg+TMath::Pi()/2);
    Yseg[0] = roseg*sin(phiseg+TMath::Pi()/2);
    Zseg[0] = zseg;
    Xseg[1] = Xseg[0]+cos(phiseg)*sin(thetaseg);
    Yseg[1] = Yseg[0]+sin(phiseg)*sin(thetaseg);
    Zseg[1] = Zseg[0]+cos(thetaseg);

    p1.SetXYZ(Xseg[0],Yseg[0],Zseg[0]);
    p2.SetXYZ(Xseg[1],Yseg[1],Zseg[1]);
}

Bool_t HParticleMetaMatcher::findIntersectionLinePlane(TVector3 &pointIntersect,
					       const TVector3 &pos, const TVector3 &dir,
					       const TVector3 &planeCenter, const  TVector3 &planeNormal) {
    // Finds the intersection point of a straight line with any plane.
    // pointIntersect: the intersection point (return value).
    // pos:            a point on the straight line.
    // dir:            direction of the straight line.
    // planeCenter:    a point on the plane.
    // planeNormal:    normal vector of the plane.

    Double_t denom = planeNormal.Dot(dir);
    if (denom != 0.0) {
	Double_t t = ((planeCenter.x() - pos.x()) * planeNormal.x() +
		      (planeCenter.y() - pos.y()) * planeNormal.y() +
		      (planeCenter.z() - pos.z()) * planeNormal.z()) / denom;

	pointIntersect = pos + (t * dir);
	return kTRUE;
    } else {
	cout<<"Warning \"findIntersection()\" : , No intersection point found : (plane || track)"<<endl;
        return kFALSE;
    }
    return kFALSE;
}


Bool_t HParticleMetaMatcher::traceToMeta(const HParticleCand* cand,TVector3 &metaHit, const TVector3 &metaCenter, const TVector3 &metaNorm)
{
    // calculate intersection of outer MdcSeg (RK) straightline with META plane
    // in sector sys

    TVector3 p1;
    TVector3 p2;
    TVector3 dir;

    Double_t phiLab   = cand->getPhi2();    // outer MdcSeg (RK) 0-360 deg  lab sys

    Double_t phi = fmod(phiLab,60.F) + 60;  // sec sys

    Double_t theta = cand->getTheta2();  // deg
    Double_t r     = cand->getR2();      // mm
    Double_t z     = cand->getZ2();      // mm

    calcSegPoints(p1,p2,phi*TMath::DegToRad(),theta*TMath::DegToRad(),r,z);  // 2 points from segment in sec sys

    dir = p2 - p1;

    return  findIntersectionLinePlane(metaHit,p1,dir,metaCenter,metaNorm); // in sec system
}

Bool_t HParticleMetaMatcher::isInRpcCell(HGeomVector& vmod, Int_t s,Int_t col,Int_t c)
{
    //  check if vmodule (mod system) is inside cell (2D) Trapez
    //
    //       p1---------------p2
    //      b\                 /
    //        \               /
    //        p0-------------p3
    //              a
    //
    // checks if HGeomVector* vmodule


    HGeomVector& p0 = cellRPC[s][col][c][0];
    HGeomVector& p1 = cellRPC[s][col][c][1];

    if (p0.Y() > vmod.Y()) return kFALSE;   // first check y range
    if (p1.Y() < vmod.Y()) return kFALSE;

    // left side                            // x cut as function of y
    Double_t a = p1.X() - p0.X();
    Double_t b = p1.Y() - p0.Y();
    Double_t x = p0.X() + a * (( p0.Y() - vmod.Y())/b);

    if (x < vmod.X()) return kFALSE;

    HGeomVector& p2 = cellRPC[s][col][c][2];
    HGeomVector& p3 = cellRPC[s][col][c][3];

    // right side
    a = p2.X() - p3.X();
    b = p2.Y() - p3.Y();
    Double_t x1 = p3.X() + a * (( p3.Y() - vmod.Y())/b);

    if (x1 > vmod.X()) return kFALSE;

    if(fDebug){
	cout<<"isInRpcCell  s "<<s<<" col "<<col<<" c "<<setw(2)<<c
	    <<" p1 "<<setw(10)<<  p0.X() <<" , "<<setw(10)<< p0.Y()
	    <<" p2 "<<setw(10)<<  p1.X() <<" , "<<setw(10)<< p1.Y()
	    <<" p3 "<<setw(10)<<  p2.X() <<" , "<<setw(10)<< p2.Y()
	    <<" p4 "<<setw(10)<<  p3.X() <<" , "<<setw(10)<< p3.Y()
	    <<" hit "<<setw(10)<<vmod.X()<<" , "<<setw(10)<<vmod.Y()<<" calc x "<<setw(10)<<x<<" "<<x1
	    <<endl;
    }

    return kTRUE;
}
Bool_t HParticleMetaMatcher::isInEmcCell(HGeomVector& vmod,Int_t s,Int_t c)
{

    //  check if vmodule (mod system) is inside cell (2D) square
    //
    //        p1----p2
    //        |     |
    //        |     |
    //        p0----p3
    //
    //

    HGeomVector& p0 = cellEMC[s][c][0];
    HGeomVector& p2 = cellEMC[s][c][2];

    if (p0.Y() > vmod.Y()) return kFALSE;   // first check y range
    if (p2.Y() < vmod.Y()) return kFALSE;

    if (p0.X() < vmod.X()) return kFALSE;
    if (p2.X() > vmod.X()) return kFALSE;

    if(fDebug){
	HGeomVector& p1 = cellEMC[s][c][1];
	HGeomVector& p3 = cellEMC[s][c][3];
	Int_t cell = HEmcDetector::getCellFromPosition(c);
	cout<<"isInEmcCell  s "<<s<<" c "<<setw(3)<<cell<<" pos "<<setw(3)<<c
	    <<" p1 "<<setw(10)<<  p0.X() <<" , "<<setw(10)<< p0.Y()
	    <<" p2 "<<setw(10)<<  p1.X() <<" , "<<setw(10)<< p1.Y()
	    <<" p3 "<<setw(10)<<  p2.X() <<" , "<<setw(10)<< p2.Y()
	    <<" p4 "<<setw(10)<<  p3.X() <<" , "<<setw(10)<< p3.Y()
	    <<" hit "<<setw(10)<<vmod.X()<<" , "<<setw(10)<<vmod.Y()
	    <<endl;
        vector<Int_t> vcells;
	getEmcCellArray(s,cell,vcells);
	cout<<"neighboughrs: "<<flush;
	for(UInt_t k = 0; k < vcells.size(); k++){
	    if(k == 4)cout<<"|"<<vcells[k]<<"| "<<flush;
	    else      cout<<vcells[k]<<" "<<flush;
	}
        cout<<endl;
    }

    return kTRUE;
}

Bool_t HParticleMetaMatcher::isInTofCell(HGeomVector& vmod,Int_t s,Int_t m, Int_t c)
{

    //  check if vmodule (mod system) is inside cell (2D) square
    //
    //        p1----p2
    //        |     |
    //        |     |
    //        p0----p3
    //
    //

    HGeomVector& p0 = cellTOF[s][m][c][0];
    HGeomVector& p2 = cellTOF[s][m][c][2];

    if (p0.Y() > vmod.Y()) return kFALSE;   // first check y range
    if (p2.Y() < vmod.Y()) return kFALSE;

    if (p0.X() < vmod.X()) return kFALSE;
    if (p2.X() > vmod.X()) return kFALSE;

    if(fDebug){
	HGeomVector& p1 = cellTOF[s][m][c][1];
	HGeomVector& p3 = cellTOF[s][m][c][3];
	cout<<"isInTofCell  s "<<s<<" mod "<<m<<" c "<<setw(3)<<c
	    <<" p1 "<<setw(10)<<  p0.X() <<" , "<<setw(10)<< p0.Y()
	    <<" p2 "<<setw(10)<<  p1.X() <<" , "<<setw(10)<< p1.Y()
	    <<" p3 "<<setw(10)<<  p2.X() <<" , "<<setw(10)<< p2.Y()
	    <<" p4 "<<setw(10)<<  p3.X() <<" , "<<setw(10)<< p3.Y()
	    <<" hit "<<setw(10)<<vmod.X()<<" , "<<setw(10)<<vmod.Y()
	    <<endl;
    }

    return kTRUE;
}

void HParticleMetaMatcher::getEmcCellArray(Int_t s,Int_t cell,vector<Int_t>& vcells)
{
    // returns neighbour cell indices for a given cell
    // in the order
    //
    //  0   1    2
    //  3   4(cell) 5
    //  6   7    8
    //  into vcells (-1 if cell does not exist)
    vcells.clear();

    for(Int_t j = -18; j < 17; j += 17){ // 3 rows of 17 cells each
	for(Int_t i = 0; i < 3; i ++){   // 3 neigbouring cells in each row
	    if(HEmcDetector::getPositionFromCell(cell + j + i) == -1 ) vcells.push_back(-1);
	    else                                                       vcells.push_back(cell + j + i);
	}
    }
}

Bool_t HParticleMetaMatcher::init()
{
    // Get categories and coordinate transformations.
    // Call this after HADES has been initialized.

    if(!gHades) {
        Error("init()", "HADES not initialized.");
        return kFALSE;
    }

    HRuntimeDb *rtdb=gHades->getRuntimeDb();

    if(!rtdb) {
        Error("init()", "HADES runtime database not found.");
        return kFALSE;
    }

    fCatRpcCluster  = HCategoryManager::getCategory(catRpcCluster);
    fCatEmcCluster  = HCategoryManager::getCategory(catEmcCluster);
    fCatTofCluster  = HCategoryManager::getCategory(catTofCluster);
    fCatTofHit      = HCategoryManager::getCategory(catTofHit);

    if(!fCatEmcCluster) Warning("init()", "Category catEmcCluster not found.");
    if(!fCatRpcCluster) Warning("init()", "Category catRpcCluster not found.");
    if(!fCatTofCluster) Warning("init()", "Category catTofCluster not found.");
    if(!fCatTofHit)     Warning("init()", "Category catTofHit not found.");


    fRpcGeometry = (HRpcGeomPar*)    rtdb->getContainer("RpcGeomPar");
    pGeomCellPar = (HRpcGeomCellPar*)rtdb->getContainer("RpcGeomCellPar");
    pSpecGeomPar = (HSpecGeomPar*)   rtdb->getContainer("SpecGeomPar");
    fEmcGeometry = (HEmcGeomPar *)   rtdb->getContainer("EmcGeomPar");
    fTofGeometry = (HTofGeomPar *)   rtdb->getContainer("TofGeomPar");

    return kTRUE;
}

Bool_t HParticleMetaMatcher::reinit() {

    if(pGeomCellPar){
	DPlanesRpc  = pGeomCellPar->getDPlanes();
    }

    // Geometry transformation from lab to sector coordinate system.
    for(Int_t i = 0; i < 6; i++) {
        labSecTrans[i] = pSpecGeomPar->getSector(i)->getTransform();
    }

    //-------------------------------------------------------------
    // get RPC cell points
    {
	for(Int_t s = 0; s < 6; s++) {
            Int_t col = 0;
            Int_t cell = 0;
	    
	    HModGeomPar* fmodgeom = fRpcGeometry->getModule(s);
	    HGeomCompositeVolume* fMod    = fmodgeom->getRefVolume();
	    for(Int_t c = 0; c < fMod->getNumComponents(); c++) {

		HGeomVolume* fVol = fMod->getComponent(c);
		if(cell == 31) {
		    col++;
		    cell = 0;
                    continue;
		}
		if(0&& s == 0){

		cout<<" col "<<5-col<<" cell " <<cell<<"----------------------------------"<<endl;
		HGeomVector v = fVol->getTransform().getTransVector();
		cout<<c<<" mod  (x,y,z) "<<v.X()<<" "<<v.Y()<<" "<<v.Z()<<" "<<fVol->GetName()<<" "<<fVol->getShape()<<" "<<fVol->getNumPoints()
		    << " p0 "<<fVol->getPoint(0)->getX()<< " , "<<fVol->getPoint(0)->getY()
		    << " p1 "<<fVol->getPoint(1)->getX()<< " , "<<fVol->getPoint(1)->getY()
		    << " p2 "<<fVol->getPoint(2)->getX()<< " , "<<fVol->getPoint(2)->getY()
		    << " p3 "<<fVol->getPoint(3)->getX()<< " , "<<fVol->getPoint(3)->getY()
		    <<endl;
		}
		
		for(Int_t i = 0; i < 4; i++) {  // copy 4 points only
		    cellRPC[s][5-col][cell][i] = *fVol->getPoint(i);
		}
		cell++;
	    }
	}
    }

    //-------------------------------------------------------------
    // get TOF cell points
    {
	for(Int_t s = 0; s < 6; s ++) {
	    for(Int_t m = 0; m < 8; m ++) {

		HModGeomPar* fmodgeom = fTofGeometry->getModule(s,m);
		HGeomCompositeVolume* fMod    = fmodgeom->getRefVolume();
		for(Int_t c = 0; c < fMod->getNumComponents(); c++) {

		    HGeomVolume* fVol=fMod->getComponent(c);
                    HGeomVector v = fVol->getTransform().getTransVector();

		    for(Int_t i = 0; i < 4; i++) {  // copy 4 points only
		        cellTOF[s][m][c][i] = v + *fVol->getPoint(i);
		    }

		    if(0&& s == 0){
			cout<<" mod "<<m<<" c "<<c<<" mod  (x,y,z) "<<setw(8)<<v.X()<<" "<<setw(8)<<v.Y()<<" "<<setw(8)<<v.Z()<<" "<<fVol->GetName()<<" "<<fVol->getShape()<<" "<<fVol->getNumPoints()
			    << " p0 "<<setw(8)<<cellTOF[s][m][c][0].X()<< " , "<<setw(8)<<cellTOF[s][m][c][0].Y()
			    << " p1 "<<setw(8)<<cellTOF[s][m][c][1].X()<< " , "<<setw(8)<<cellTOF[s][m][c][1].Y()
			    << " p2 "<<setw(8)<<cellTOF[s][m][c][2].X()<< " , "<<setw(8)<<cellTOF[s][m][c][2].Y()
			    << " p3 "<<setw(8)<<cellTOF[s][m][c][3].X()<< " , "<<setw(8)<<cellTOF[s][m][c][3].Y()
			    <<endl;
		    }
		}
	    }
	}
    }
    //-------------------------------------------------------------
    // get EMC cell points
    {
	for(Int_t s = 0; s < 6; s++) {
	    HModGeomPar* fmodgeom = fEmcGeometry->getModule(s);
	    HGeomCompositeVolume* fMod = fmodgeom->getRefVolume();

            Int_t n = 0;
	    for(Int_t c = 0; c < fMod->getNumComponents(); c++) {
		HGeomVolume* fVol=fMod->getComponent(c);
		if(fVol == NULL || fVol->getNumPoints() != 8) {
		    continue;
		}
		HGeomVector v = fVol->getTransform().getTransVector();
		if(0&& s == 0){
		    cout<<c<<" "<<n<<" mod  (x,y,z) "<<v.X()<<" "<<v.Y()<<" "<<v.Z()
			<< " p0 "<<fVol->getPoint(0)->getX()<< " , "<<fVol->getPoint(0)->getY()
			<< " p1 "<<fVol->getPoint(1)->getX()<< " , "<<fVol->getPoint(1)->getY()
			<< " p2 "<<fVol->getPoint(2)->getX()<< " , "<<fVol->getPoint(2)->getY()
			<< " p3 "<<fVol->getPoint(3)->getX()<< " , "<<fVol->getPoint(3)->getY()
			<<endl;
		}
		for(Int_t i = 0; i < 4; i++) {  // copy 4 points only
		    cellEMC[s][n][i] = v + *fVol->getPoint(i);
		}
                if(0&& s == 0){
		    cout<<c<<" "<<n<<" mod  (x,y,z) "<<v.X()<<" "<<v.Y()<<" "<<v.Z()
			<< " p0 "<<cellEMC[s][n][0].getX()<< " , "<<cellEMC[s][n][0].getY()
			<< " p1 "<<cellEMC[s][n][1].getX()<< " , "<<cellEMC[s][n][1].getY()
			<< " p2 "<<cellEMC[s][n][2].getX()<< " , "<<cellEMC[s][n][2].getY()
			<< " p3 "<<cellEMC[s][n][3].getX()<< " , "<<cellEMC[s][n][3].getY()
			<<endl;
		}
                n++;
	    }
	}
    }

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

    for(Int_t iSec = 0; iSec < 6; iSec++) {
	if(fRpcGeometry) {
	    HModGeomPar *module = fRpcGeometry->getModule(iSec,0);
	    if(module) {
		// Transformation module->lab.
                HGeomTransform modTrans(module->getLabTransform());
                // Combine with transformation lab->sector
                // to get the transformation from module to sector.
                modTrans.transTo(labSecTrans[iSec]);

                modSecTransRpc[iSec] = modTrans;

                // Centre point plane 0
                HGeomVector r0_mod0(0.0, 0.0, -0.5*DPlanesRpc);  // plane col 0,2,4
                // Normal vector.
                HGeomVector rz_mod(0.0, 0.0, 1.0);
                HGeomVector nRpc = modTrans.transFrom(rz_mod) - modTrans.transFrom(r0_mod0);
                normVecRpc[iSec][0].SetXYZ(nRpc.getX(), nRpc.getY(), nRpc.getZ());
                HGeomVector cRpc = modTrans.transFrom(r0_mod0);
                centerVecRpc[iSec][0].SetXYZ(cRpc.getX(), cRpc.getY(), cRpc.getZ());

                // Centre point plane 1
                HGeomVector r0_mod1(0.0, 0.0, 0.5*DPlanesRpc);   // plane col 1,3,5

		nRpc = modTrans.transFrom(rz_mod) - modTrans.transFrom(r0_mod1);
                normVecRpc[iSec][1].SetXYZ(nRpc.getX(), nRpc.getY(), nRpc.getZ());
                cRpc = modTrans.transFrom(r0_mod1);
                centerVecRpc[iSec][1].SetXYZ(cRpc.getX(), cRpc.getY(), cRpc.getZ());

                HGeomVector r0_mod2(0.0, 0.0, 0);                 // plane between cols for cluster type 2

		nRpc = modTrans.transFrom(rz_mod) - modTrans.transFrom(r0_mod2);
                normVecRpc[iSec][2].SetXYZ(nRpc.getX(), nRpc.getY(), nRpc.getZ());
                cRpc = modTrans.transFrom(r0_mod2);
                centerVecRpc[iSec][2].SetXYZ(cRpc.getX(), cRpc.getY(), cRpc.getZ());

	    }
        }

	if (fEmcGeometry) {
	    HModGeomPar *pmodgeom = fEmcGeometry->getModule(iSec);
	    HGeomTransform modTrans(pmodgeom->getLabTransform());

	    modTrans.transTo(labSecTrans[iSec]);

	    modSecTransEmc[iSec] = modTrans;

	    HGeomVector r0_mod(0.0, 0.0, 0.);
	    HGeomVector rz_mod(0.0, 0.0, 1.);
	    HGeomVector nEmc = modTrans.transFrom(rz_mod) - modTrans.transFrom(r0_mod);
	    normVecEmc[iSec].SetXYZ(nEmc.getX(), nEmc.getY(), nEmc.getZ()) ;
	    HGeomVector cEmc = modTrans.transFrom(r0_mod);
	    centerVecEmc[iSec].SetXYZ(cEmc.getX(), cEmc.getY(), cEmc.getZ()) ;
	}


	if (fTofGeometry) {
	    for(Int_t iTofMod = 0; iTofMod < 8; iTofMod++) {
		HModGeomPar *module = fTofGeometry->getModule(iSec,iTofMod);
		if(module) {
		    HGeomTransform modTrans(module->getLabTransform());
		    modTrans.transTo(labSecTrans[iSec]);

		    modSecTransTof[iSec][iTofMod] = modTrans;

                    HGeomVector r0_mod(0.0, 0.0, 0.0);
                    HGeomVector rz_mod(0.0, 0.0, 1.0);
                    HGeomVector nTof = modTrans.transFrom(rz_mod) - modTrans.transFrom(r0_mod);
                    normVecTof[iSec][iTofMod].SetXYZ(nTof.getX(), nTof.getY(), nTof.getZ());
                    HGeomVector cTof = modTrans.transFrom(r0_mod);
                    centerVecTof[iSec][iTofMod].SetXYZ(cTof.getX(), cTof.getY(), cTof.getZ());
                }
            }
        }
    }
    return kTRUE;
}


const HRpcCluster* HParticleMetaMatcher::recalcRpc(HParticleCand* cand)
{
    // returns HRpcCluster pointer to an object owned by HParticleMetaMatcher
    // which has been reconstructed from an intersection of HParticleCand outer
    // segment of Runge Kutta with the META detector plane. The original META
    // hit position can be restorded using RK dx, dy from HParticleCand.
    // Works for candidates which have cand->isRpcClstUsed() set.
    // The hit point of RK on the META plane can be retrieved by
    // calling
    // const HGeomVector&  getTrackMetaMod()
    // const HGeomVector&  getTrackMetaSec()
    // const HGeomVector&  getTrackMetaLab()
    // directly after calling this function.


    if(!cand) return 0;
    if(cand->getOuterSegInd() == -1) return 0;
    if(!cand->isRpcClstUsed()) return 0;


    if(fDebug){ cout <<"matchRpc--------------------------------"<<endl;}
    Int_t sec   = cand->getSector();
    Int_t mod0  = cand->getMetaModule(0);
    Int_t mod1  = cand->getMetaModule(1);
    Int_t cell0 = cand->getMetaCell(0);
    Int_t cell1 = cand->getMetaCell(1);

    if(mod0 == -1) return 0;

    Short_t type = 0;
    if(mod0 > -1) type = 1;
    if(mod1 > -1) type = 2;

    Float_t dx = cand->getRkMetaDx();
    Float_t dy = cand->getRkMetaDy();

    if(fDebug){
	HRpcCluster* pRpc = 0 ;

	if(fCatRpcCluster) {
	    pRpc = (HRpcCluster*) fCatRpcCluster ->getObject(cand->getRpcInd());
	}

	if(pRpc){

	    Float_t x,y,z;
	    pRpc->getXYZLab(x,y,z);
	    cout<<"rpc sec "<<pRpc->getSector()<<" col "<<pRpc->getColumn1()<<" cell "<<setw(2)<<pRpc->getCell1()<<" type "<<type
		<<" ,pos mod "<<setw(10)<<pRpc->getXMod()
		<<" "<<setw(10)<<pRpc->getYMod()
		<<" "<<setw(10)<<pRpc->getZMod()
		<<" ,pos sec "<<setw(10)<<pRpc->getXSec()
		<<" "<<setw(10)<<pRpc->getYSec()
		<<" "<<setw(10)<<pRpc->getZSec()
		<<" ,pos lab "<<setw(10)<<x
		<<" "<<setw(10)<<y
		<<" "<<setw(10)<<z
		<<" ,theta "<<setw(10)<<pRpc->getTheta()
		<<" phi "<<setw(10)<<pRpc->getPhi()
		<<" dx "<<dx
		<<" dy "<<dy
                <<" type "<<pRpc->getClusterType()
		<<endl;
	}
    }

    TVector3 metaHit;
    Int_t ind = mod0%2; // finding the correct plane
    if(mod1 > -1) ind = 2;  // both cells, plane in middle

    traceToMeta(cand,metaHit,centerVecRpc[sec][ind],normVecRpc[sec][ind]);
    HGeomVector vsec(metaHit.X(),metaHit.Y(),metaHit.Z());    // hit point mdc on META sec sys

    trackMetaSec = vsec;

    HGeomTransform& TransSec = pSpecGeomPar->getSector(sec)->getTransform();
    HGeomVector vlab  = TransSec.transFrom(vsec); // lab sys

    trackMetaLab = vlab;

    HGeomTransform& modToLab = fRpcGeometry->getModule(sec,0)->getLabTransform();
    HGeomVector vmod = modToLab.transTo(vlab);  // mod sys

    trackMetaMod = vmod;

    vmod.setX(vmod.X()+dx);   // reconstruct the meta hit
    vmod.setY(vmod.Y()+dy);   // reconstruct the meta hit

    vlab = modToLab.transFrom(vmod); // lab sys
    vsec = TransSec.transTo(vlab);   // sec sys

    Double_t theta = TMath::RadToDeg() * TMath::ATan2(TMath::Sqrt(vlab.X()*vlab.X()+vlab.Y()*vlab.Y()),vlab.Z());
    Double_t phi   = TMath::RadToDeg() * TMath::ATan2(vlab.Y(),vlab.X());
    if (phi < 0.) phi += 360.;

    if(fDebug){
	cout<<"mdc sec "<<sec<<" col "<<mod0<<" cell "<<setw(2)<<cell0<<" type "<<type
	    <<" ,pos mod "<<setw(10)<<vmod.X()
	    <<" "<<setw(10)<<vmod.Y()
	    <<" "<<setw(10)<<vmod.Z()
	    <<" ,pos sec "<<setw(10)<<vsec.X()
	    <<" "<<setw(10)<<vsec.Y()
	    <<" "<<setw(10)<<vsec.Z()
	    <<" ,pos lab "<<setw(10)<<vlab.X()
	    <<" "<<setw(10)<<vlab.Y()
	    <<" "<<setw(10)<<vlab.Z()
	    <<" ,theta "<<setw(10)<<theta
	    <<" phi "<<setw(10)<<phi
            <<" type "<<type
	    <<endl;
    }
    rpccluster.setDetID1(sec,mod0,cell0);
    rpccluster.setDetID2(sec,mod1,cell1);
    rpccluster.setTheta(theta);
    rpccluster.setPhi(phi);
    rpccluster.setXYZLab(vlab.X(),vlab.Y(),vlab.Z());
    rpccluster.setRSEC(vsec.X(),vsec.Y(),vsec.Z());
    rpccluster.setCluster(cand->getTof(),1.,vmod.X(),vmod.Y(),vmod.Z());
    rpccluster.setClusterType(type);

    return &rpccluster;

}

const HTofCluster* HParticleMetaMatcher::recalcTof(HParticleCand* cand)
{
    // returns HTofCluster pointer to an object owned by HParticleMetaMatcher
    // which has been reconstructed from an intersection of HParticleCand outer
    // segment of Runge Kutta with the META detector plane. The original META
    // hit position can be restorded using RK dx, dy from HParticleCand.
    // Works for candidates which have cand->isTofClstUsed() or cand->isTofHitUsed() set.
    // The hit point of RK on the META plane can be retrieved by
    // calling
    // const HGeomVector&  getTrackMetaMod()
    // const HGeomVector&  getTrackMetaSec()
    // const HGeomVector&  getTrackMetaLab()
    // directly after calling this function.


    if(!cand) return 0;
    if(cand->getOuterSegInd() == -1) return 0;
    if(! (cand->isTofClstUsed() || cand->isTofHitUsed()) ) return 0;


    if(fDebug){ cout <<"matchTof--------------------------------"<<endl;}
    Int_t sec   = cand->getSector();
    Int_t mod0  = cand->getMetaModule(0);
    Int_t mod1  = cand->getMetaModule(1);
    Int_t cell0 = cand->getMetaCell(0);

    if(mod0 == -1) return 0;

    Short_t type = 0;
    if(mod0 > -1) type = 1;
    if(mod1 > -1) type = 2;

    Float_t dx = cand->getRkMetaDx();
    Float_t dy = cand->getRkMetaDy();

    if(fDebug){
	HTofCluster* pTof = 0 ;

	if(cand->isTofClstUsed() && fCatTofCluster) {
	    pTof = (HTofCluster*) fCatTofCluster ->getObject(cand->getTofClstInd());
	}
	if(cand->isTofHitUsed() && fCatTofHit) {
	    pTof = (HTofCluster*) fCatTofHit ->getObject(cand->getTofHitInd());
	}

	if(pTof){

	    Float_t x,y,z;
	    Float_t ph,th;
            pTof->getTheta(th);
            pTof->getPhi(ph);
	    pTof->getXYZLab(x,y,z);
	    cout<<"tof sec "<<(Int_t)pTof->getSector()<<" mod "<<(Int_t)pTof->getModule()<<" cell "<<setw(2)<<(Int_t)pTof->getCell()<<" type "<<type
		<<" ,pos mod "<<setw(10)<<pTof->getXpos()
		<<" ,pos lab "<<setw(10)<<x
		<<" "<<setw(10)<<y
		<<" "<<setw(10)<<z
		<<" ,theta "<<setw(10)<<th
		<<" phi "<<setw(10)<<ph
		<<" dx "<<dx
		<<" dy "<<dy
		<<endl;
	}
    }

    TVector3 metaHit;

    traceToMeta(cand,metaHit,centerVecTof[sec][mod0],normVecTof[sec][mod0]);
    HGeomVector vsec(metaHit.X(),metaHit.Y(),metaHit.Z());    // hit point mdc on META sec sys

    trackMetaSec = vsec;

    HGeomTransform& TransSec = pSpecGeomPar->getSector(sec)->getTransform();
    HGeomVector vlab  = TransSec.transFrom(vsec); // lab sys

    trackMetaLab = vlab;

    HGeomTransform& modToLab = fTofGeometry->getModule(sec,mod0)->getLabTransform();
    HGeomVector vmod = modToLab.transTo(vlab);  // mod sys

    trackMetaMod = vmod;

    vmod.setX(vmod.X()+dx);   // reconstruct the meta hit
    vmod.setY(vmod.Y()+dy);   // reconstruct the meta hit

    vlab = modToLab.transFrom(vmod); // lab sys
    vsec = TransSec.transTo(vlab);   // sec sys

    Double_t theta = TMath::RadToDeg() * TMath::ATan2(TMath::Sqrt(vlab.X()*vlab.X()+vlab.Y()*vlab.Y()),vlab.Z());
    Double_t phi   = TMath::RadToDeg() * TMath::ATan2(vlab.Y(),vlab.X());
    if (phi < 0.) phi += 360.;

    if(fDebug){
	cout<<"mdc sec "<<sec<<" mod "<<mod0<<" cell "<<setw(2)<<cell0<<" type "<<type
	    <<" ,pos mod "<<setw(10)<<vmod.X()
	    <<" ,pos lab "<<setw(10)<<vlab.X()
	    <<" "<<setw(10)<<vlab.Y()
	    <<" "<<setw(10)<<vlab.Z()
	    <<" ,theta "<<setw(10)<<theta
	    <<" phi "<<setw(10)<<phi
	    <<endl;
    }
    tofcluster.setSector(sec);
    tofcluster.setModule(mod0);
    tofcluster.setCell(cell0);
    tofcluster.setTheta(theta);
    tofcluster.setPhi(phi);
    tofcluster.setXYZLab(vlab.X(),vlab.Y(),vlab.Z());
    tofcluster.setXpos(vmod.X());
    tofcluster.setTof(cand->getTof());
    tofcluster.setEdep(cand->getTofdEdx());
    tofcluster.setClusterSize(type);

    return &tofcluster;

}

const HEmcCluster* HParticleMetaMatcher::recalcEmc(HParticleCand* cand)
{
    // returns HEmcCluster pointer to an object owned by HParticleMetaMatcher
    // which has been reconstructed from an intersection of HParticleCand outer
    // segment of Runge Kutta with the META detector plane. The original META
    // hit position can be restorded using RK dx dy from HParticleCand.
    // Works for candidates which have cand->getEmcInd()!=-1 set.
    // The hit point of RK on the META plane can be retrieved by
    // calling
    // const HGeomVector&  getTrackMetaMod()
    // const HGeomVector&  getTrackMetaSec()
    // const HGeomVector&  getTrackMetaLab()
    // directly after calling this function.
    //


    if(!cand) return 0;
    if(cand->getOuterSegInd() == -1) return 0;
    if(cand->getEmcInd() == -1) return 0;


    if(fDebug){ cout <<"matchEmc--------------------------------"<<endl;}
    Int_t sec   = cand->getSector();

    Short_t type = 0;

    Float_t dx = cand->getRkMetaDxEmc();
    Float_t dy = cand->getRkMetaDyEmc();

    if(fDebug){
	HEmcCluster* pEmc = 0 ;

	if(cand->getEmcInd() > -1 && fCatEmcCluster) {
	    pEmc = (HEmcCluster*) fCatEmcCluster ->getObject(cand->getEmcInd());
	}

	if(pEmc){

	    Float_t x,y,z;
	    pEmc->getXYZLab(x,y,z);
	    cout<<"emc sec "<<(Int_t)pEmc->getSector()<<" cell "<<setw(3)<<(Int_t)pEmc->getCell()<<" type "<<(Int_t)pEmc->getNCells()
		<<" ,pos mod "<<setw(10)<<pEmc->getXMod()
		<<" "<<setw(10)<<pEmc->getYMod()
		<<" ,pos lab "<<setw(10)<<x
		<<" "<<setw(10)<<y
		<<" "<<setw(10)<<z
		<<" ,theta "<<setw(10)<<pEmc->getTheta()
		<<" phi "<<setw(10)<<pEmc->getPhi()
		<<" dx "<<dx
		<<" dy "<<dy
		<<endl;
	}
    }

    TVector3 metaHit;

    traceToMeta(cand,metaHit,centerVecEmc[sec],normVecEmc[sec]);
    HGeomVector vsec(metaHit.X(),metaHit.Y(),metaHit.Z());    // hit point mdc on META sec sys

    trackMetaSec = vsec;

    HGeomTransform& TransSec = pSpecGeomPar->getSector(sec)->getTransform();
    HGeomVector vlab  = TransSec.transFrom(vsec); // lab sys

    trackMetaLab = vlab;

    HGeomTransform& modToLab = fEmcGeometry->getModule(sec,0)->getLabTransform();
    HGeomVector vmod = modToLab.transTo(vlab);  // mod sys

    trackMetaMod = vmod;

    vmod.setX(vmod.X()+dx);   // reconstruct the meta hit
    vmod.setY(vmod.Y()+dy);   // reconstruct the meta hit

    vlab = modToLab.transFrom(vmod); // lab sys
    vsec = TransSec.transTo(vlab);   // sec sys

    Double_t theta = TMath::RadToDeg() * TMath::ATan2(TMath::Sqrt(vlab.X()*vlab.X()+vlab.Y()*vlab.Y()),vlab.Z());
    Double_t phi   = TMath::RadToDeg() * TMath::ATan2(vlab.Y(),vlab.X());
    if (phi < 0.) phi += 360.;

    if(fDebug){
	cout<<"mdc sec "<<sec<<" cell "<<setw(3)<<0<<" type "<<type
	    <<" ,pos mod "<<setw(10)<<vmod.X()
	    <<" "<<setw(10)<<vmod.Y()
	    <<" ,pos lab "<<setw(10)<<vlab.X()
	    <<" "<<setw(10)<<vlab.Y()
	    <<" "<<setw(10)<<vlab.Z()
	    <<" ,theta "<<setw(10)<<theta
	    <<" phi "<<setw(10)<<phi
	    <<endl;
    }
    emccluster.setSector(sec);
    emccluster.setTheta(theta);
    emccluster.setPhi(phi);
    emccluster.setXYZLab(vlab.X(),vlab.Y(),vlab.Z());
    emccluster.setXYMod(vmod.X(),vmod.Y());
    emccluster.setTime(cand->getEmcTime());
    emccluster.setEnergy(cand->getEmcEnergy());

    return &emccluster;

}

Bool_t HParticleMetaMatcher::predictRpcCell(HParticleCand* cand,HGeomVector& hit0,HGeomVector& hit1,Int_t& s,Int_t& col0,Int_t& cell0,Int_t& col1,Int_t& cell1)
{
    // Reconstructs an intersection of HParticleCand outer
    // segment of Runge Kutta with the META detector plane.
    // 2 hits (hit0,hit2) on the Rpc planes of column (0,2,4)
    // and (1,3,5). Possible hit RPC cells are returned
    // by coordinates sector, (col0,cell0) (col1,cell1).
    // colX and cellX are -1 if no cell has been hit on this plane.
    // Works with all candidates which provide outer segments.


    if(!cand) return 0;
    if(cand->getOuterSegInd() == -1) return 0;

    if(fDebug){ cout <<"predictRpcCell--------------------------------"<<endl;}
    Int_t sec   = cand->getSector();


    if(fDebug && cand->getRpcInd() != -1){
	Int_t mod0  = cand->getMetaModule(0);
	Int_t mod1  = cand->getMetaModule(1);
	Int_t c0    = cand->getMetaCell(0);
	Int_t c1    = cand->getMetaCell(1);

	if(mod0 == -1) return 0;

	Short_t type = 0;
	if(mod0 > -1) type = 1;
	if(mod1 > -1) type = 2;


	HRpcCluster* pRpc = 0 ;

	if(fCatRpcCluster) {
	    pRpc = (HRpcCluster*) fCatRpcCluster ->getObject(cand->getRpcInd());
	}

	if(pRpc){

	    cout<<"rpc sec "<<pRpc->getSector()
		<<" col0 "<<mod0<<" cell0 "<<setw(2)<<c0
		<<" col1 "<<mod1<<" cell1 "<<setw(2)<<c1
		<<" type "<<type
		<<" ,pos mod "<<setw(10)<<pRpc->getXMod()
		<<" "<<setw(10)<<pRpc->getYMod()
		<<" "<<setw(10)<<pRpc->getZMod()
		<<endl;
	}
    }

    TVector3 metaHit;
    HGeomVector vmod[2];

    HGeomTransform& TransSec = pSpecGeomPar->getSector(sec)->getTransform();
    HGeomTransform& modToLab = fRpcGeometry->getModule(sec,0)->getLabTransform();

    for(Int_t i = 0; i < 2; i++){
	traceToMeta(cand,metaHit,centerVecRpc[sec][i],normVecRpc[sec][i]);
	HGeomVector vsec(metaHit.X(),metaHit.Y(),metaHit.Z());    // hit point mdc on META sec sys
	HGeomVector vlab  = TransSec.transFrom(vsec); // lab sys
	vmod[i] = modToLab.transTo(vlab);  // mod sys
    }

    hit0 = vmod[0];
    hit1 = vmod[1];
    s    = sec;
    col0 = col1 = -1;
    cell0= cell1= -1;

    Int_t ct = 0;
    for(Int_t col = 0; col < 6; col++){
	for(Int_t c = 0; c < RPCMAXCELL; c++){
	    HGeomVector v = vmod[col%2];
	    if(!isInRpcCell(v,sec,col,c)) continue;

	    if     (ct == 0) { col0 = col; cell0 = c;}
	    else if(ct == 1) { col1 = col; cell1 = c;}
	    else  {
		cout <<"predictRpcCell() : more than 2 Rpc cells found!"<<endl;
	    }
	    ct++;

	}
    }

    return ct > 0 ? kTRUE : kFALSE;

}

Bool_t HParticleMetaMatcher::predictTofCell(HParticleCand* cand,HGeomVector& hit0,HGeomVector& hit1,Int_t& s,Int_t& mod0,Int_t& cell0,Int_t& mod1,Int_t& cell1)
{
    // Reconstructs an intersection of HParticleCand outer
    // segment of Runge Kutta with the META detector plane.
    // 2 hits (hit0,hit2) on the TOF module planes cand be
    // returned. Possible hit TOF cells are returned
    // by coordinates sector, (mod0,cell0) (mod1,cell1).
    // colX and cellX are -1 if no cell has been hit on this plane.
    // The hit points will be 0,0,0 in this case.
    // Works with all candidates which provide outer segments.

    if(!cand) return 0;
    if(cand->getOuterSegInd() == -1) return 0;

    if(fDebug){ cout <<"predictTofCell--------------------------------"<<endl;}
    Int_t sec   = cand->getSector();



    if(fDebug){
	HTofCluster* pTof    = 0 ;

	if(cand->getTofClstInd() != -1 && fCatTofCluster) {
	    pTof = (HTofCluster*) fCatTofCluster ->getObject(cand->getTofClstInd());
	}
	if(cand->getTofHitInd() != -1 && fCatTofHit) {
	    pTof = (HTofCluster*) fCatTofHit ->getObject(cand->getTofHitInd());
	}

	if(pTof){
	    Int_t mod0  = cand->getMetaModule(0);
	    Int_t mod1  = cand->getMetaModule(1);
	    Int_t c0    = cand->getMetaCell(0);
	    Int_t c1    = cand->getMetaCell(1);

	    Short_t type = 0;
	    if(mod0 > -1) type = 1;
	    if(mod1 > -1) type = 2;

	    cout<<"tof sec "<<(Int_t)pTof->getSector()
		<<" mod0 "<<mod0<<" cell0 "<<setw(2)<<c0
		<<" mod1 "<<mod1<<" cell1 "<<setw(2)<<c1
		<<" type "<<type
		<<" ,pos mod "<<setw(10)<<pTof->getXpos()
		<<endl;
	}
    }


    TVector3 metaHit;
    HGeomVector vmod[8];

    HGeomTransform& TransSec = pSpecGeomPar->getSector(sec)->getTransform();

    for(Int_t i = 0; i < 8; i++){ // modules
	HGeomTransform& modToLab=fTofGeometry->getModule(sec,i)->getLabTransform();
	traceToMeta(cand,metaHit,centerVecTof[sec][i],normVecTof[sec][i]);
	HGeomVector vsec(metaHit.X(),metaHit.Y(),metaHit.Z());    // hit point mdc on META sec sys
	HGeomVector vlab  = TransSec.transFrom(vsec); // lab sys
	vmod[i] = modToLab.transTo(vlab);  // mod sys
    }

    s    = sec;
    mod0 = mod1 = -1;
    cell0 = cell1= -1;

    hit0.setXYZ(0,0,0);
    hit1.setXYZ(0,0,0);


    Int_t ct = 0;
    for(Int_t m = 0; m < 8; m++){
	for(Int_t c = 0; c < 8; c++){
	    if(!isInTofCell(vmod[m],sec,m,c)) continue;

	    if     (ct == 0) { mod0 = m; cell0 = c; hit0 = vmod[m];}
	    else if(ct == 1) { mod1 = m; cell1 = c; hit1 = vmod[m];}
	    else  {
		cout <<"predictTofCell() : more than 2 Tof rods found!"<<endl;
	    }
	    ct++;

	}
    }

    return ct > 0 ? kTRUE : kFALSE;

}

Bool_t HParticleMetaMatcher::predictEmcCell(HParticleCand* cand,HGeomVector& hit,Int_t& s,Int_t& pos,Int_t& cell)
{
    // Reconstructs an intersection of HParticleCand outer
    // segment of Runge Kutta with the META detector plane.
    // 1 hit on the EMC plane can be returned.
    // The possible hit cell returned  by coordinates sector, pos (0-162) , cell (0,254).
    // pos and cell are -1 if no cell has been hit on this plane.
    // Works with all candidates which provide outer segments.

    if(!cand) return 0;
    if(cand->getOuterSegInd() == -1) return 0;

    if(fDebug){ cout <<"predictEmcCell--------------------------------"<<endl;}
    Int_t sec   = cand->getSector();


    if(fDebug && cand->getEmcInd() != -1){

	HEmcCluster* pEmc = 0 ;

	if(fCatEmcCluster) {
	    pEmc = (HEmcCluster*) fCatEmcCluster ->getObject(cand->getEmcInd());
	}

	if(pEmc){

	    cout<<"emc sec "<<(Int_t)pEmc->getSector()<<" cell "<<setw(3)<<(Int_t)pEmc->getCell()<<" pos "<<setw(3)<<HEmcDetector::getPositionFromCell(pEmc->getCell())<<" type "<<(Int_t)pEmc->getNCells()
		<<" ,pos mod "<<setw(10)<<pEmc->getXMod()
		<<" "<<setw(10)<<pEmc->getYMod()
                <<" radius "<<cand->getMetaMatchRadiusEmc()
		<<endl;
	}
    }

    TVector3 metaHit;
    HGeomVector vmod;

    HGeomTransform& TransSec = pSpecGeomPar->getSector(sec)->getTransform();
    HGeomTransform& modToLab = fRpcGeometry->getModule(sec,0)->getLabTransform();

    traceToMeta(cand,metaHit,centerVecEmc[sec],normVecEmc[sec]);
    HGeomVector vsec(metaHit.X(),metaHit.Y(),metaHit.Z());    // hit point mdc on META sec sys
    HGeomVector vlab  = TransSec.transFrom(vsec); // lab sys



    vmod = modToLab.transTo(vlab);  // mod sys

    hit  = vmod;
    s    = sec;
    pos  = -1;
    cell = -1;

    for(Int_t c = 0; c < EMCMAXCELL; c++){
	if(!isInEmcCell(vmod,sec,c)) continue;
	pos = c;
        cell = HEmcDetector::getCellFromPosition(pos);
	break;
    }
    return pos > -1 ? kTRUE : kFALSE;

}