#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 "hmdcsizescells.h"

#include "hparticlecand.h"
#include "hparticletool.h"
#include "hmdcseg.h"

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

#include "TMath.h"
#include "TString.h"


//_HADES_CLASS_DESCRIPTION
/////////////////////////////////////////////////////////////
//
//  HParticleMetamatcher
//
//
// usage :
// 
// 
// 
// 
// 
// #include "hades.h"
// #include "hloop.h"
// #include "hdst.h"
// #include "hcategory.h"
// #include "hcategorymanager.h"
// #include "htaskset.h"
// #include "hgeomvector.h"
//
//
// #include "hparticlemetamatcher.h"
//
// #include "hparticlecand.h"
// #include "hparticleevtinfo.h"
// #include "htofcluster.h"
// #include "hrpccluster.h"
// #include "hemccluster.h"
//
// #include "hparticledef.h"
//
// #include "TString.h"
//
// #include <iostream>
// #include <iomanip>
//
// using namespace std;
//
// Int_t loopDST(
// 	      TString infileList="/lustre/hades/dst/mar19/gen5/ag158ag/3200A/083/root/be1908301523101.hld_dst_mar19_accepted.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_gen5_13042021.root";
// 	TString paramSource      = "root"; // root, ascii, oracle
// 	TString paramrelease     = "MAR19_dst_gen5";  // 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(0&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?
// 		    }
// 		    //------------------------------------------------------------------------
// 		    // MDC + tracking
// 		    HParticleWireManager& wM = matcher->getWireManager(); // object to retrive infos about number of wires in event,
// 		    // wires used by candidate, shared wires between candidates etc.
// 		    //wM.setDebug(); // print options of wire usage
// 		    //wM.setWireRange(2)  // blocked +- range arround used wire in layer (0-n)(default value for sharedWires(....) function
// 		    //Int_t nw     = wM.getNWires    ()           { return nWires;}     // number of wires in event
// 		    //Int_t nwused = wM.getNWiresUsed()           { return nWiresUsed;}
// 		    //Int_t n      = wM.isUsedNtimes(Int_t s,Int_t m,Int_t l,Int_t c);  // how many time this cell has been used ? size 6:4:6:220
// 		    //vector<HParticleWireInfo>& vwinfo =  wM.getWireInfo();
// 		    //Bool_t wM.getWireInfo(UInt_t i,HParticleWireInfo& w,HParticleCand* c = 0);  // index of HParticleCand -> WireInfo
// 		    //Bool_t wM.sharedWires(HParticleWireInfo& w1     ,HParticleWireInfo& w2     ,Int_t io=0,Int_t range = -1); // compare 2 WireInfo objects if they shared wires
// 		    //Bool_t wM.sharedWires(HParticleCand* c1,HParticleCand* c2,Int_t io=0,Int_t range = -1); // compare 2 HParticleCand objects if they shared wires
//
// 		    if(1){
// 			HGeomVector hitmdc;
// 			HGeomVector hitmdcsec;
// 			HGeomVector hitmdclab;
// 			Int_t module = 0; // 0-3
// 			matcher->predictMdcHit   (cand,hitmdc   ,module);  // MDC coord system
// 			matcher->predictMdcHitSec(cand,hitmdcsec,module);  // sector coord system
// 			matcher->predictMdcHitLab(cand,hitmdclab,module);  // lab coord system
// 		    }
//
// 		    //------------------------------------------------------------------------
// 		} // if cand
// 	    } // end cand loop
// 	} // if cand cat
//     } // end eventloop
//
//     //-------------------------------------------------
//     // added for the task
//     if(gHades)gHades->finalizeTasks();
//     //-------------------------------------------------
//
//     delete gHades;
//     return 0;
// }
/////////////////////////////////////////////////////////////

ClassImp(HParticleWireInfo)
ClassImp(HParticleWireManager)
ClassImp(HParticleMetaMatcher)


//----------------------------------------------------------------
HParticleWireManager::HParticleWireManager()
{
    nWRange = 2;
    ctevent = 0;
    ctCand     = 0;
    ctCandUsed = 0;
    nWires    = 0;
    nWiresUsed= 0;

    pParticleCandCat = 0;
    pSegCat = 0;
    fDebug  = kFALSE;
};

void HParticleWireManager::clear()
{

    // clear the count table per wire
    ctCand = 0;
    ctCandUsed = 0;

    for(Int_t s = 0; s < 6; s ++){
	ctSecCand[s] = 0;
	ctSecCandUsed[s] = 0;
	ctSecWire[s] = 0;
	ctSecWireDouble[s] = 0;
	ctSecSegDouble[s] = 0;
	for(Int_t m = 0; m < 4; m ++){
	    for(Int_t l = 0; l < 6; l ++){
		for(Int_t c = 0; c < MAXWIRE; c ++){
		    wires[s][m][l][c] = 0;
		}
	    }
	}
    }
    nWires = 0;
    nWiresUsed = 0;
    vWireInfo.clear();
    mCandWireInfo.clear();

}

void HParticleWireManager::init()
{
    pParticleCandCat    = HCategoryManager::getCategory(catParticleCand   ,kTRUE,"catParticleCand");
    pSegCat             = HCategoryManager::getCategory(catMdcSeg         ,kTRUE,"catMdcSeg");
    
}

void HParticleWireManager::fillWire(Int_t s,Int_t m,Int_t l,Int_t c,Int_t flag)
{
    // count up a single cell
    // flag = 0 normal case
    // flag = 1 inner segment has been reused
    if(flag == 1) return;
    if(c < 0) return;
    wires[s][m][l][c] ++ ;
    if(wires[s][m][l][c] == 1) {
	nWires++;
        ctSecWire[s]++;
    } else {
	ctSecWireDouble[s]++;
    }
    nWiresUsed++;
}  

void HParticleWireManager::fillWireInfo(HParticleWireInfo& wi,Int_t flag)
{
    // fill the wireinfo
    vWireInfo.push_back(wi);
    if(flag == 1){
	ctSecSegDouble[wi.c->getSector()]++;
    }
}


Int_t HParticleWireManager::execute()
{
    ctevent++;
    clear();
    sizescells          = HMdcSizesCells::getExObject();
    if(!sizescells)       { return 0;}
    if(!pParticleCandCat) { return 0;}
    // -------------------------------------------------------------------------


    HParticleCand*    c;
    vector<Int_t> vsegind;       // inner segment ind for kIsUsed particles

    //-------------------------------------------------------------------------
    // pre calculate the list of cells for each candidate
    //-------------------------------------------------------------------------

    Int_t n = pParticleCandCat->getEntries();
    for(Int_t i = 0; i < n; i ++)
    {
	c = HCategoryManager::getObject(c,pParticleCandCat,i);
	ctSecCand[c->getSector()]++;
	if(c->isFlagBit(kIsUsed)) {
	    vsegind.push_back(c->getInnerSegInd());
	    ctSecCandUsed[c->getSector()]++;
	}
    }

    for(Int_t i = 0; i < n; i ++)
    {
	c = HCategoryManager::getObject(c,pParticleCandCat,i);
	if(c)
	{
	    //---------------------------------------------
	    // calculate expected hits

	    HMdcSeg* pSeg[2] = { 0 ,0};
	    Double_t x1,y1,z1,x2,y2,z2;
	    Double_t x3,y3,z3,x4,y4,z4;

	    if(pSegCat){
		pSeg[0] = HParticleTool::getMdcSeg(c->getInnerSegInd());
		pSeg[1] = HParticleTool::getMdcSeg(c->getOuterSegInd());
	    } else {
		HParticleTool::calcSegPoints(c,x1,y1,z1,x2,y2,z2,kFALSE); // r,z, phi in Lab [deg], theta [deg]  -> 2 space points in sec coord system
		HParticleTool::calcSegPoints(c,x3,y3,z3,x4,y4,z4,kTRUE ); // r,z, phi in Lab [deg], theta [deg]  -> 2 space points in sec coord system
	    }

	    Int_t flag = 0;
	    if(!c->isFlagBit(kIsUsed))
	    {
		if(find(vsegind.begin(),vsegind.end(),c->getInnerSegInd()) != vsegind.end()) flag = 1;  // reused one Segment
		else    vsegind.push_back(c->getInnerSegInd()); // new segment

	    }

	    HParticleWireInfo w;
	    w.reset();
	    Bool_t goodAll = kTRUE;
	    Int_t sec = c->getSector();
	    Float_t cell1,cell2;
	    Int_t c1,c2;
            c1=c2=-1;
	    for(Int_t mod = 0; mod < 4; mod ++) { // inner + outer MDC
		for(Int_t lay = 0; lay < 6; lay ++) { // 6 layers
		    Bool_t good = kFALSE;

		    if(pSegCat){
			//---------------------------------------
			// HMdcSeg  in input. get wire info
			good = kTRUE;
			Int_t io = 0;         // 0 = inner seg, 1 = outer seg
			if(mod > 1) io = 1;

			Int_t l = lay;        // segment stores layers of 2 MDC
			if(mod == 1 || mod == 3) l += 6; // recalc to seg

			if(!pSeg[io]) {  continue; } // seg missing

			Int_t n = pSeg[io]->getNCells(l);  // number of cells in layer
			if(n > 0){
			    cell1 = pSeg[io]->getCell(l,0);
			    if(n == 2) cell2 = pSeg[io]->getCell(l,1);
			    else cell2 = cell1;

			    c1 = (Int_t)cell1;
			    c2 = (Int_t)cell2;
			}
			//---------------------------------------

		    } else {

			//---------------------------------------
			// HMdcSeg not in input. Recalculate fired cells
			// from intersection (prediction):
			//
			// calculation of cells crossed by line x1,y1,z1-x2,y2,z2
			// x1,y1,z1, x2,y2,z2 - sector coordinate system
			// return kFALSE if failed.
			// The line can cross 2 cell in a layer at max. If
			// only 1 cell is hit  cell1 and cell2 are equal

			if(mod<3) good =(*sizescells)[sec][mod][lay].calcCrossedCells(x1,y1,z1,x2,y2,z2,cell1,cell2);
			else      good =(*sizescells)[sec][mod][lay].calcCrossedCells(x3,y3,z3,x4,y4,z4,cell1,cell2);

			if(good){  // cells in range 0 - ncell
			    c1 = (Int_t)cell1;
			    c2 = (Int_t)cell2;

			    if(c1 > c2) { // sort
				Int_t tmp = c2;
				c2 = c1;
				c2 = tmp;
			    }
			} else { // both cells < 0 or > ncell
			    c1 =-1;
                            c2 =-1;
			}

		    }
		    //---------------------------------------
		    // filling wireinfo
		    if(c1 > -1) {
			w.ar[mod][lay][0] = c1;
                        fillWire(sec,mod,lay,c1,flag);
		    }
		    if(c2 > -1 && c2 != c1) {
			w.ar[mod][lay][1] = c2;
			fillWire(sec,mod,lay,c1,flag);
		    }

		    if(!good)  goodAll = kFALSE;
		    //---------------------------------------
		} // end loop lay
	    } // end loop of segment cells

	    w.c = c;
	    if(goodAll) w.isGood = kTRUE;
	    if(c->isFlagBit(kIsUsed)) w.isUsed = kTRUE;
	    fillWireInfo(w,flag);
	    Int_t ind = vWireInfo.size() -1;
	    mCandWireInfo[c] = vWireInfo[ind];

	}
    }
    flagDoubleUsedWires();  // set isGood flag per candidate
    //-------------------------------------------------------------------------


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

    return 0;
}


Bool_t HParticleWireManager::getWireInfo(UInt_t i,HParticleWireInfo& w,HParticleCand* c)
{
    if(i < vWireInfo.size()){

	w = vWireInfo[i];
	if(c != 0){
	    if(c != w.c) cout<<"getWireInfo(): Candidate does not match WireInfo!";

	}
	return kTRUE;

    } else {
	return kFALSE;
    }
}

Bool_t HParticleWireManager::sharedWires(HParticleWireInfo& w1,HParticleWireInfo& w2,Int_t io,Int_t range)
{
    // returns kTRUE if inner/outer Segments (io = 0 ->inner, = 1 -> outer, =2 ->both)
    // share wires inside a range +- wire number. The intern range nWRange is used
    // unless the provide range != -1 ,


    //----------------------------------------------------------------
    // check if any wire (+- wRange) of
    // cand has been used more than once.
    Int_t nWR = nWRange;
    if(range !=-1) { nWR = abs(range);}

    Int_t minMod = 0;
    Int_t maxMod = 4;

    if     (io == 0)  maxMod = 2;
    else if(io == 1)  minMod = 2;

    for(Int_t m = minMod ; m < maxMod; m ++){
	for(Int_t l = 0 ; l < 6; l ++){
	    Int_t c1_1 = w1.ar[m][l][0];
	    Int_t c2_1 = w1.ar[m][l][1];
	    Int_t c1_2 = w2.ar[m][l][0];
	    Int_t c2_2 = w2.ar[m][l][1];
	    Int_t min = c1_1 - nWR;
	    Int_t max = c2_1 + nWR;
	    if(min < 0)   min = 0;
	    if(max > MAXWIRE-1) max = MAXWIRE -1;
	    if(c1_2 <= max && c1_2 >= min ) return kTRUE;
	    if(c2_2 <= max && c2_2 >= min ) return kTRUE;
	}
    }
    return kFALSE;

}

Bool_t HParticleWireManager::sharedWires(HParticleCand* c1,HParticleCand* c2,Int_t io,Int_t range)
{
    // returns kTRUE if inner/outer Segments (io = 0 ->inner, = 1 -> outer, =2 ->both)
    // Segments share wires inside a range +- wire number. The intern range nWRange is used
    // unless the provide range != -1 ,

    if(mCandWireInfo.find(c1) == mCandWireInfo.end()) {
	cout<<"sharedWires() : HParticelCand not found in map!"<<endl;
	return kFALSE;
    }
    if(mCandWireInfo.find(c2) == mCandWireInfo.end()) {
	cout<<"sharedWires() : HParticelCand not found in map!"<<endl;
	return kFALSE;
    }
    HParticleWireInfo& w1 = mCandWireInfo[c1];
    HParticleWireInfo& w2 = mCandWireInfo[c2];
    return sharedWires(w1,w2, range);
}

Int_t HParticleWireManager::isUsedNtimes(Int_t s,Int_t m,Int_t l,Int_t c)
{
    // return how often a wire has been used
    // in the event (by HParticleCandidate)
    // fillWire() should be called before.
    // returns usage count of the wire
    // ( 0 if nothing filled,
    //   >0 used by HParticleCand (or MdcSeg if available)
    // )
    return wires[s][m][l][c];
}

void HParticleWireManager::flagDoubleUsedWires()
{
    // loop over all wireinfo objects
    // flag multiple used wires. The
    // comparison is exact (no range
    // nWRange is used)

    ctCand     = vWireInfo.size();
    ctCandUsed = 0;

    for(UInt_t i = 0 ; i < vWireInfo.size(); i++){
	HParticleWireInfo& wi = vWireInfo[i];
	Int_t s = wi.c->getSector();

	//----------------------------------------------------------------
	// check if any wire of
        // cand has been used more than once.
	for(Int_t m = 0 ; m < 4; m ++){
	    for(Int_t l = 0 ; l < 6; l ++){
		Int_t c1 = wi.ar[m][l][0];
		Int_t c2 = wi.ar[m][l][1];
		Int_t min = c1 ;
		Int_t max = c2 ;

		for(Int_t i = min; i <= max ; i++){
		    if (wires[s][m][l][i] > 1){
			wi.isGood = kFALSE;
			if(m<2) { wi.isSharedInner = kTRUE; wi.nSharedInner++;}
			else    { wi.isSharedOuter = kTRUE; wi.nSharedOuter++;}
		    }
		} // cell
	    } // layer
	} // module
	//----------------------------------------------------------------
    }
}

void  HParticleWireManager::print ()
{
    if(!fDebug) return;
    cout<<"---------------------------------------------"<<endl;
    cout<<setw(6)<<ctevent
	<<setw(5)<<" nWires "    << nWires
	<<setw(5)<<" nWiresUsed "<< nWiresUsed
	<<setw(3)<<" nCand "     <<ctCand
	<<setw(3)<<" nCandUsed " <<ctCandUsed
	<<endl;
    cout<<"cand        s "<<setw(3)<<ctSecCand[0]      <<" "<<setw(3)<<ctSecCand[1]      <<" "<<setw(3)<<ctSecCand[2]      <<" "<<setw(3)<<ctSecCand[3]      <<" "<<setw(3)<<ctSecCand[4]      <<" "<<setw(3)<<ctSecCand[5]<<endl;
    cout<<"cand used   s "<<setw(3)<<ctSecCandUsed[0]  <<" "<<setw(3)<<ctSecCandUsed[1]  <<" "<<setw(3)<<ctSecCandUsed[2]  <<" "<<setw(3)<<ctSecCandUsed[3]  <<" "<<setw(3)<<ctSecCandUsed[4]  <<" "<<setw(3)<<ctSecCandUsed[5]<<endl;
    cout<<"seg double  s "<<setw(3)<<ctSecSegDouble[0] <<" "<<setw(3)<<ctSecSegDouble[1] <<" "<<setw(3)<<ctSecSegDouble[2] <<" "<<setw(3)<<ctSecSegDouble[3] <<" "<<setw(3)<<ctSecSegDouble[4] <<" "<<setw(3)<<ctSecSegDouble[5]<<endl;
    cout<<"wire        s "<<setw(3)<<ctSecWire[0]      <<" "<<setw(3)<<ctSecWire[1]      <<" "<<setw(3)<<ctSecWire[2]      <<" "<<setw(3)<<ctSecWire[3]      <<" "<<setw(3)<<ctSecWire[4]      <<" "<<setw(3)<<ctSecWire[5]<<endl;
    cout<<"wire double s "<<setw(3)<<ctSecWireDouble[0]<<" "<<setw(3)<<ctSecWireDouble[1]<<" "<<setw(3)<<ctSecWireDouble[2]<<" "<<setw(3)<<ctSecWireDouble[3]<<" "<<setw(3)<<ctSecWireDouble[4]<<" "<<setw(3)<<ctSecWireDouble[5]<<endl;
    cout<<"---------------------------------------------"<<endl;

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





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){
	cout<<"HParticleMetaMatcher::getRpcCellGeom() : cell "<<cell<<" out of bounds"<<endl;
	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 &cellTOF[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];
}

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

    HParticleTool::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, Float_t xRes, Float_t yRes)
{
    //  check if vmodule (mod system) is inside cell (2D) Trapez
    //
    //       p1---------------p2
    //      b\                 /
    //        \               /
    //        p0-------------p3
    //              a
    //
    // checks if HGeomVector* vmodule
    // xRes,yRes : matching tolerance in x and y (default 0)

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

    if (p0.Y() > vmod.Y() + yRes){
	if(fDebug){
	    cout<<"isInRpcCell  s "<<s<<" col "<<col<<" c "<<setw(2)<<c
		<<" p0 "<<setw(10)<<  p0.X() <<" , "<<setw(10)<< p0.Y()
		<<" p1 "<<setw(10)<<  p1.X() <<" , "<<setw(10)<< p1.Y()
		<<" hit "<<setw(10)<<vmod.X()<<" , "<<setw(10)<<vmod.Y()<<" miss : hit below cell"<<endl;
	}
	return kFALSE;   // first check y range
    }
    if (p1.Y() < vmod.Y() - yRes) {
	if(fDebug){
	    cout<<"isInRpcCell  s "<<s<<" col "<<col<<" c "<<setw(2)<<c
		<<" p0 "<<setw(10)<<  p0.X() <<" , "<<setw(10)<< p0.Y()
		<<" p1 "<<setw(10)<<  p1.X() <<" , "<<setw(10)<< p1.Y()
		<<" hit "<<setw(10)<<vmod.X()<<" , "<<setw(10)<<vmod.Y()<<" miss : hit above cell"<<endl;
	}
	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() - xRes ) {
	if(fDebug){
	    cout<<"isInRpcCell  s "<<s<<" col "<<col<<" c "<<setw(2)<<c
		<<" p0 "<<setw(10)<<  p0.X() <<" , "<<setw(10)<< p0.Y()
		<<" p1 "<<setw(10)<<  p1.X() <<" , "<<setw(10)<< p1.Y()
		<<" hit "<<setw(10)<<vmod.X()<<" , "<<setw(10)<<vmod.Y()<<" x cut "<<x<<" miss : x hit out of cell left side"<<endl;
	}
	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() + xRes) {
	if(fDebug){
	    cout<<"isInRpcCell  s "<<s<<" col "<<col<<" c "<<setw(2)<<c
		<<" p0 "<<setw(10)<<  p0.X() <<" , "<<setw(10)<< p0.Y()
		<<" p1 "<<setw(10)<<  p1.X() <<" , "<<setw(10)<< p1.Y()
		<<" p2 "<<setw(10)<<  p2.X() <<" , "<<setw(10)<< p2.Y()
		<<" p3 "<<setw(10)<<  p3.X() <<" , "<<setw(10)<< p3.Y()
		<<" hit "<<setw(10)<<vmod.X()<<" , "<<setw(10)<<vmod.Y()<<" x cut "<<x1<<", miss : x hit out of cell right side"<<endl;
	}
	return kFALSE;
    }
    if(fDebug){
	cout<<"isInRpcCell  s "<<s<<" col "<<col<<" c "<<setw(2)<<c
	    <<" p0 "<<setw(10)<<  p0.X() <<" , "<<setw(10)<< p0.Y()
	    <<" p1 "<<setw(10)<<  p1.X() <<" , "<<setw(10)<< p1.Y()
	    <<" p2 "<<setw(10)<<  p2.X() <<" , "<<setw(10)<< p2.Y()
	    <<" p3 "<<setw(10)<<  p3.X() <<" , "<<setw(10)<< p3.Y()
	    <<" hit "<<setw(10)<<vmod.X()<<" , "<<setw(10)<<vmod.Y()<<" calc x "<<setw(10)<<x<<" "<<x1<<" hit matched"
	    <<endl;
    }

    return kTRUE;
}


Bool_t HParticleMetaMatcher::isInEmcCell(HGeomVector& vmod,Int_t s,Int_t c, Float_t xRes, Float_t yRes)
{

    //  check if vmodule (mod system) is inside cell (2D) square
    //
    //        p1----p2
    //        |     |
    //        |     |
    //        p0----p3
    //
    //
    // xRes,yRes : matching tolerance in x and y (default 0)

    Int_t pos = HEmcDetector::getPositionFromCell(c);

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

    if (p0.Y() > vmod.Y() + yRes) {
	if(fDebug){
	    cout<<"isInEmcCell  s "<<s<<" c "<<setw(2)<<c
		<<" p0 "<<setw(10)<<  p0.X() <<" , "<<setw(10)<< p0.Y()
		<<" p2 "<<setw(10)<<  p2.X() <<" , "<<setw(10)<< p2.Y()
		<<" hit "<<setw(10)<<vmod.X()<<" , "<<setw(10)<<vmod.Y()<<" miss : hit below cell"<<endl;
	}
	return kFALSE;   // first check y range
    }

    if (p2.Y() < vmod.Y() - yRes) {
	if(fDebug){
	    cout<<"isInEmcCell  s "<<s<<" c "<<setw(2)<<c
		<<" p0 "<<setw(10)<<  p0.X() <<" , "<<setw(10)<< p0.Y()
		<<" p2 "<<setw(10)<<  p2.X() <<" , "<<setw(10)<< p2.Y()
		<<" hit "<<setw(10)<<vmod.X()<<" , "<<setw(10)<<vmod.Y()<<" miss : hit above cell"<<endl;
	}
	return kFALSE;
    }

    if (p0.X() < vmod.X() - xRes) {
	if(fDebug){
	    cout<<"isInEmcCell  s "<<s<<" c "<<setw(2)<<c
		<<" p0 "<<setw(10)<<  p0.X() <<" , "<<setw(10)<< p0.Y()
		<<" p2 "<<setw(10)<<  p2.X() <<" , "<<setw(10)<< p2.Y()
		<<" hit "<<setw(10)<<vmod.X()<<" , "<<setw(10)<<vmod.Y()<<" miss : x hit out of cell left side"<<endl;
	}
	return kFALSE;
    }

    if (p2.X() > vmod.X() + xRes) {
	if(fDebug){
	    cout<<"isInEmcCell  s "<<s<<" c "<<setw(2)<<c
		<<" p0 "<<setw(10)<<  p0.X() <<" , "<<setw(10)<< p0.Y()
		<<" p2 "<<setw(10)<<  p2.X() <<" , "<<setw(10)<< p2.Y()
		<<" hit "<<setw(10)<<vmod.X()<<" , "<<setw(10)<<vmod.Y()<<" miss : x hit out of cell right side"<<endl;
	}
	return kFALSE;
    }

    if(fDebug){
	cout<<"isInEmcCell  s "<<s<<" c "<<setw(3)<<c<<" pos "<<setw(3)<<pos
	    <<" p0 "<<setw(10)<<  p0.X() <<" , "<<setw(10)<< p0.Y()
	    <<" p2 "<<setw(10)<<  p2.X() <<" , "<<setw(10)<< p2.Y()
	    <<" hit "<<setw(10)<<vmod.X()<<" , "<<setw(10)<<vmod.Y()<<" hit matched"
	    <<endl;
        vector<Int_t> vcells;
	getEmcCellArray(s,c,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, Float_t xRes, Float_t yRes)
{

    //  check if vmodule (mod system) is inside cell (2D) square
    //
    //        p1----p2
    //        |     |
    //        |     |
    //        p0----p3
    //
    //
    // xRes,yRes : matching tolerance in x and y (default 0)

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

    if (p0.Y() > vmod.Y() + yRes) {
	if(fDebug){
	    cout<<"isInTofCell  s "<<s<<" m "<<m<<" c "<<setw(2)<<c
		<<" p0 "<<setw(10)<<  p0.X() <<" , "<<setw(10)<< p0.Y()
		<<" p2 "<<setw(10)<<  p2.X() <<" , "<<setw(10)<< p2.Y()
		<<" hit "<<setw(10)<<vmod.X()<<" , "<<setw(10)<<vmod.Y()<<" miss : hit below cell"<<endl;
	}
	return kFALSE;   // first check y range
    }

    if (p2.Y() < vmod.Y() - yRes) {
	if(fDebug){
	    cout<<"isInTofCell  s "<<s<<" m "<<m<<" c "<<setw(2)<<c
		<<" p0 "<<setw(10)<<  p0.X() <<" , "<<setw(10)<< p0.Y()
		<<" p2 "<<setw(10)<<  p2.X() <<" , "<<setw(10)<< p2.Y()
		<<" hit "<<setw(10)<<vmod.X()<<" , "<<setw(10)<<vmod.Y()<<" miss : hit above cell"<<endl;
	}
	return kFALSE;
    }

    if (p0.X() < vmod.X() - xRes) {
	if(fDebug){
	    cout<<"isInTofCell  s "<<s<<" m "<<m<<" c "<<setw(2)<<c
		<<" p0 "<<setw(10)<<  p0.X() <<" , "<<setw(10)<< p0.Y()
		<<" p2 "<<setw(10)<<  p2.X() <<" , "<<setw(10)<< p2.Y()
		<<" hit "<<setw(10)<<vmod.X()<<" , "<<setw(10)<<vmod.Y()<<" miss : x hit out of cell left side"<<endl;
	}
	return kFALSE;
    }

    if (p2.X() > vmod.X() + xRes) {
	if(fDebug){
	    cout<<"isInTofCell  s "<<s<<" m "<<m<<" c "<<setw(2)<<c
		<<" p0 "<<setw(10)<<  p0.X() <<" , "<<setw(10)<< p0.Y()
		<<" p2 "<<setw(10)<<  p2.X() <<" , "<<setw(10)<< p2.Y()
		<<" hit "<<setw(10)<<vmod.X()<<" , "<<setw(10)<<vmod.Y()<<" miss : x hit out of cell right side"<<endl;
	}
	return kFALSE;
    }

    if(fDebug){
	cout<<"isInTofCell  s "<<s<<" mod "<<m<<" c "<<setw(3)<<c
	    <<" p0 "<<setw(10)<<  p0.X() <<" , "<<setw(10)<< p0.Y()
	    <<" p2 "<<setw(10)<<  p2.X() <<" , "<<setw(10)<< p2.Y()
	    <<" hit "<<setw(10)<<vmod.X()<<" , "<<setw(10)<<vmod.Y()<<" hit matched"
	    <<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");

    rtdb->getContainer("MdcTrackGFieldPar");
    rtdb->getContainer("MagnetPar");

    rtdb->getContainer("MdcRawStruct");
    rtdb->getContainer("MdcGeomStruct");
    rtdb->getContainer("MdcLookupRaw");
    rtdb->getContainer("MdcLookupGeom");
    rtdb->getContainer("MdcLayerCorrPar");
    rtdb->getContainer("MdcLayerGeomPar");
    rtdb->getContainer("MdcGeomPar");
    rtdb->getContainer("SpecGeomPar");

    wiremanager.init();

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

		TString name=fVol->GetName();

                if(name.Length()!=5) continue;

		col =(Int_t)(name[2]-'0')-1;

		if((Int_t)(name[3]) > 57){
		    cell=(Int_t)(name[3]-'A')-1 +10;  // A-V
		} else cell=(Int_t)(name[3]-'0')-1;   // 0-9

		if(cell > RPCMAXCELL-1) continue;

		Int_t colMap = col;    // take care about remaping as in hgeant
                if(col == 0) colMap=4;
                if(col == 1) colMap=5;
                if(col == 4) colMap=0;
                if(col == 5) colMap=1;

		if(0&&s == 0){

		    cout<<" col "<<col<<" "<<colMap<<" 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()<<" , "<<fVol->getPoint(0)->getZ()
			<< " p1 "<<fVol->getPoint(1)->getX()<< " , "<<fVol->getPoint(1)->getY()<<" , "<<fVol->getPoint(1)->getZ()
			<< " p2 "<<fVol->getPoint(2)->getX()<< " , "<<fVol->getPoint(2)->getY()<<" , "<<fVol->getPoint(2)->getZ()
			<< " p3 "<<fVol->getPoint(3)->getX()<< " , "<<fVol->getPoint(3)->getY()<<" , "<<fVol->getPoint(3)->getZ()
			<<endl;
		}

		for(Int_t i = 0; i < 4; i++) {  // copy 4 points only
		    cellRPC[s][colMap][cell][i] = *fVol->getPoint(i);
		}

	    }
	}
    }

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

    fSizesCells = HMdcSizesCells::getObject();
    Bool_t retval = fSizesCells->initContainer();

    return retval;
}

Int_t HParticleMetaMatcher::execute()
{
    wiremanager.execute();
    wiremanager.print();
    return 0;
}


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

    Int_t   isInCell = 0     ;   // Flag for cell outliers (xMod out of geometric cell bounds)
	                         // If cluster type == 2
	                         //   0 if both hits in the cluster are out of bounds
	                         //   1 if both are in
	                         //   2 if only first hit is in
	                         //   3 if only second hit is in
	                         // If cluster type ==1
	                         //   0 if out of bounds
	                         //   1 otherwise

    if(type == 1){ // 1 cell
         if(isInRpcCell(vmod,sec,mod0,cell0)) isInCell = 1;
    } else if(type == 2){ // 2 cell
	if(isInRpcCell(vmod,sec,mod0,cell0)){
            if(isInRpcCell(vmod,sec,mod1,cell1)) isInCell = 1;
            else                                 isInCell = 2;
	}
	if(isInRpcCell(vmod,sec,mod1,cell1)){
            if(!isInRpcCell(vmod,sec,mod0,cell0)) isInCell = 3;
	}
    }

    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);
    rpccluster.setInsideCellFlag(isInCell);


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

    Bool_t isInside = kFALSE;
    isInside = isInTofCell(vmod,sec,mod0,cell0);
    if(type == 2) isInside = isInTofCell(vmod,sec,mod1,cell1);

    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);
    if(isInside) tofcluster.setInsideCell();
    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,Int_t sys)
{
    // 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.
    // The hit can be returned in different coordinate systems:
    // sys  1 = mod (default), 2 = sector, 3 = lab


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

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

    if(sys < 0 || sys > 3) sys = 1;

    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];
    HGeomVector vsec[2];
    HGeomVector vlab[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 vsecL(metaHit.X(),metaHit.Y(),metaHit.Z());    // hit point mdc on META sec sys
	vsec[i] = vsecL;
	HGeomVector vlabL  = TransSec.transFrom(vsecL); // lab sys
	vlab[i] = vlabL;
	vmod[i] = modToLab.transTo(vlabL);  // mod sys
    }
    if(sys == 1) {
	hit0 = vmod[0];
	hit1 = vmod[1];
    } else if(sys == 2){
	hit0 = vsec[0];
	hit1 = vsec[1];
    } else if(sys == 3){
 	hit0 = vlab[0];
	hit1 = vlab[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,Int_t sys)
{
    // 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.
    // The hit can be returned in different coordinate systems:
    // sys  1 = mod (default), 2 = sector, 3 = lab

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

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

    if(sys < 0 || sys > 3) sys = 1;

    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];
    HGeomVector vsec[8];
    HGeomVector vlab[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 vsecL(metaHit.X(),metaHit.Y(),metaHit.Z());    // hit point mdc on META sec sys
        vsec[i] = vsecL;
	HGeomVector vlabL  = TransSec.transFrom(vsecL); // lab sys
	vlab[i] = vlabL;
	vmod[i] = modToLab.transTo(vlabL);  // 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;
		if(sys == 1)hit0 = vmod[m];
		if(sys == 2)hit0 = vsec[m];
		if(sys == 3)hit0 = vlab[m];
	    }
	    else if(ct == 1) {
		mod1 = m; cell1 = c;

		if(sys == 1)hit1 = vmod[m];
		if(sys == 2)hit1 = vsec[m];
		if(sys == 3)hit1 = vlab[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,Int_t sys)
{
    // 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.
    // The hit can be returned in different coordinate systems:
    // sys  1 = mod (default), 2 = sector, 3 = lab

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

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

    if(sys < 0 || sys > 3) sys = 1;

    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 = fEmcGeometry->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
    if(sys == 2) hit = vsec;

    HGeomVector vlab  = TransSec.transFrom(vsec); // lab sys
    if(sys == 3) hit = vlab;



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

    if(sys == 1) hit  = vmod;



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

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

}


Bool_t HParticleMetaMatcher::predictMdcHit(HParticleCand* cand,HGeomVector& hit,Int_t m)
{
    // returns  MDC hit [x,y] in  MDC system (z=0) for MDC module m [0-3]
    // recalculeated from r,z,phi,theta
    // return kFALSE if no outer MDC was reconstructed and m>=2

    if(cand->getOuterSegInd() < 0 && m > 2) return kFALSE;
    HMdcSizesCellsMod& sizesMod= (*fSizesCells)[cand->getSector()][m];

    Double_t x1,y1,z1,x2,y2,z2,x,y,z;
    if     (m >= 0 && m<2) HParticleTool::calcSegPoints(cand,x1,y1,z1,x2,y2,z2,kFALSE);
    else if(m >1   && m<4) HParticleTool::calcSegPoints(cand,x1,y1,z1,x2,y2,z2,kTRUE);

    z=0;
    sizesMod.calcInterTrMod(x1,y1,z1,x2,y2,z2, x,y); // x,y on plane z=0
    hit.setXYZ(x,y,z);
    return kTRUE;
}

Bool_t HParticleMetaMatcher::predictMdcHitSec(HParticleCand* cand,HGeomVector& hit,Int_t m)
{
    // returns  MDC hit [x,y,z] in sec coord. system for MDC module m [0-3]
    // recalculeated from r,z,phi,theta
    // return kFALSE if no outer MDC was reconstructed and m>=2

    if(!predictMdcHit(cand,hit,m)) return kFALSE;
    Double_t x,y,z;
    hit.getXYZ(x,y,z);
    (*fSizesCells)[cand->getSector()][m].transFromZ0(x,y,z);
    hit.setXYZ(x,y,z);
    return kTRUE;

}

Bool_t HParticleMetaMatcher::predictMdcHitLab(HParticleCand* cand,HGeomVector& hit,Int_t m)
{
    // returns  MDC hit [x,y,z] in Lab coord. system for MDC module m [0-3]
    // recalculeated from r,z,phi,theta
    // return kFALSE if no outer MDC was reconstructed and m>=2

    if(!predictMdcHitSec(cand,hit,m)) return kFALSE;

    const HGeomTransform* transLab = (*fSizesCells)[cand->getSector()].getLabTrans();
    hit = transLab->transFrom(hit);
    return kTRUE;

}