#ifndef HParticleMetaMatcher_h
#define HParticleMetaMatcher_h

#include "hgeomtransform.h"
#include "hgeomvector.h"
#include "hreconstructor.h"

#include "hrpccluster.h"
#include "hemccluster.h"
#include "htofcluster.h"
#include "hparticlecand.h"

#include "hparticleRK.h"
#include "hgeomvector.h"
#include "htool.h"


#include "TVector3.h"
#include "TCollection.h"     // kIterForward + kIterBackward

#include <vector>
#include <map>
#include <iostream>
#include <iomanip>


#define EMCMAXCELL 163
#define RPCMAXCELL 31
#define iTOFMAXCELL 3

using namespace std;

// from hydra
class HCategory;
class HMetaMatchPar;
class HRpcGeomPar;
class HRpcGeomCellPar;
class HEmcGeomPar;
class HiTofGeomPar;
class HTofGeomPar;
class HRich700DigiPar;
class HParticleCand;
class HParticleCandSim;
class HSpecGeomPar;
class HMdcSizesCells;
class HMagnetPar;
class HMdcTrackGFieldPar;
class HParticleMetaMatcher;
class HGeantKine;
class HRichCal;



//##################################################################
//##################################################################

//----------------------------------------------------------------
// store infos of wires per candidate
class HParticleWireInfo : public TObject {
public:

    Int_t ar[4][6][2]; // mod,lay,2 wires/lay
    HParticleCand* c;
    Bool_t isGood;  // no error h
    Bool_t isUsed;  // kIsUsed
    Bool_t isSharedInner;// shares wire
    Bool_t isSharedOuter;// shares wire
    Int_t  nSharedInner; // number of shared wires inner MDC
    Int_t  nSharedOuter; // number of shared wires outer MDC
    void reset(){
	c      = 0;
	isGood = kFALSE;
	isUsed = kFALSE;
	isSharedInner = kFALSE;
	isSharedOuter = kFALSE;
	nSharedInner = 0;
	nSharedOuter = 0;


	for(Int_t m = 0; m < 4; m ++){
	    for(Int_t l = 0; l < 6; l ++){
		ar[m][l][0] =-1;
		ar[m][l][1] =-1;
	    }
	}
    }
    void print() {

        cout<<"-------------------------------------------------------"<<endl;
	cout<<"isGood "<<isGood
	    <<" isUsed "<<isUsed<<" sector "<<c->getSector()<<" iSeg "<<c->getInnerSegInd()<<" oSeg "<<c->getOuterSegInd()
	    <<" isSharedInner "<<isSharedInner
	    <<" isSharedOuter "<<isSharedOuter
	    <<" nSharedInner " <<setw(3)<<nSharedInner
	    <<" nSharedOuter " <<setw(3)<<nSharedOuter
	    <<endl;
	printWires();
	cout<<"-------------------------------------------------------"<<endl;
    }

    void printWires(){
	for(Int_t m = 0; m < 4; m ++){
	    for(Int_t l = 0; l < 6; l ++){
		cout<<"m "<<m<<" l "<<l<<" c1 "<<ar[m][l][0]<<" c2 "<<ar[m][l][1]<<endl;
	    }
	}
    }

    void printSharedWires(HParticleWireInfo& w,Int_t io=2,Int_t range=0){
	Int_t nWR = abs(range);
        Int_t ctW = 0;
	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 ++){
                if(ar[m][l][0] == -1 && ar[m][l][1] == -1) continue; // no wires in layer
		Int_t c1_1 = ar  [m][l][0];
		Int_t c2_1 = ar  [m][l][1];
		Int_t c1_2 = w.ar[m][l][0];
		Int_t c2_2 = w.ar[m][l][1];

		Int_t min,max;
                min = max = -1000;

		if(c1_1 > -1 ){
		    min = c1_1 - nWR;
		    max = c1_1 + nWR;
		}

		if(c2_1 > -1 ){
		    if(min == -1000 ) {  // not yet set
			min = c2_1 - nWR;
			max = c2_1 + nWR;
		    } else {
                        if(c2_1 > c1_1) max = c2_1 + nWR;
                        else            min = c2_1 - nWR;
		    }

		}

		//check boundaries
		if(min < 0)       min = 0;
		if(max > 220 - 1) max = 220 - 1;

		if(c1_2 <= max && c1_2 >= min ) { ctW++; cout<<setw(3)<<ctW<<" m "<<m<<" l "<<l<<" c intern  ("<<ar[m][l][0]<<","<<ar[m][l][1]<<"), (min,max) =  ("<<min<<","<<max<<")  :  c extern "<<c1_2<<endl;}
		if(c2_2 <= max && c2_2 >= min ) { ctW++; cout<<setw(3)<<ctW<<" m "<<m<<" l "<<l<<" c intern  ("<<ar[m][l][0]<<","<<ar[m][l][1]<<"), (min,max) =  ("<<min<<","<<max<<")  :  c extern "<<c2_2<<endl;}

	    }
	}
    }

    Int_t nSharedWires(HParticleWireInfo& w,Int_t io=2,Int_t range=2, Bool_t quitEarly=kFALSE)
    {
	// returns n shared Wires if inner/outer Segments (io = 0 ->inner, = 1 -> outer, =2 ->both)
	// share wires inside a range +- wire number.
        // if quitEarly == kTRUE (deafult kFALSE ) will return after the first match

	//----------------------------------------------------------------
	// check if any wire (+- wRange) of
	// cand has been used more than once.
	Int_t nWR = abs(range);
        Int_t ctW = 0;
	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 = ar  [m][l][0];
		Int_t c2_1 = ar  [m][l][1];
		Int_t c1_2 = w.ar[m][l][0];
		Int_t c2_2 = w.ar[m][l][1];

		if( c1_1 ==-1 && c2_1 ==-1 ) continue ; // no cell in layer

		Int_t min,max;
                min = max = -1000;
		if(c1_1 > -1 ){
		    min = c1_1 - nWR;
		    max = c1_1 + nWR;
		}

		if(c2_1 > -1 ){
		    if(min == -1000 ) {  // not yet set
			min = c2_1 - nWR;
			max = c2_1 + nWR;
		    } else {
                        if(c2_1 > c1_1) max = c2_1 + nWR;
                        else            min = c2_1 - nWR;
		    }

		}

		if(min == -1000) continue; // no cell in  ar[m][l] , nothing to do

		//check boundaries
		if(min < 0)       min = 0;
		if(max > 220 - 1) max = 220 - 1;

		if(c1_2 <= max && c1_2 >= min ) { ctW++; if (quitEarly) return 1; }
		if(c2_2 <= max && c2_2 >= min ) { ctW++; if (quitEarly) return 1; }
	    }
	}

	return ctW;

    }


    Int_t getNWiresSeg(Int_t io){

	// io = 0 : inner seg, 1 : outer seg : 2  both segs

	Int_t n[2] = {0};

        if(io < 0 || io > 2) return -1;

	for(Int_t m = 0; m < 4; m ++){
	    for(Int_t l = 0; l < 6; l ++){
		if(ar[m][l][0] > -1) n[m%2]++;
		if(ar[m][l][1] > -1) n[m%2]++;
	    }
	}

        if(io == 0) return n[0];
        else if(io == 1) return n[1];
        else             return n[0]+n[1];

    }

    Int_t getNLaySeg(Int_t io){

	// io = 0 : inner seg, 1 : outer seg : 2  both segs

	Int_t n[2] = {0};

        if(io < 0 || io > 2) return -1;

	for(Int_t m = 0; m < 4; m ++){
	    for(Int_t l = 0; l < 6; l ++){
		if(ar[m][l][0] > -1 || ar[m][l][1] > -1) n[m%2]++;
	    }
	}

        if(io == 0) return n[0];
        else if(io == 1) return n[1];
        else             return n[0]+n[1];

    }

    Int_t getNWiresMod(Int_t mod){

	Int_t n[4] = {0};

        if(mod < 0 || mod > 3) return -1;

	for(Int_t m = 0; m < 4; m ++){
	    for(Int_t l = 0; l < 6; l ++){
		if(ar[m][l][0] > -1) n[m]++;
		if(ar[m][l][1] > -1) n[m]++;
	    }
	}

        return n[mod];
    }

    Int_t getNLayMod(Int_t mod){

	Int_t n[4] = {0};

        if(mod < 0 || mod > 3) return -1;

	for(Int_t m = 0; m < 4; m ++){
	    for(Int_t l = 0; l < 6; l ++){
		if(ar[m][l][0] > -1 || ar[m][l][1] > -1) n[m]++;
	    }
	}

        return n[mod];
    }


    ClassDef(HParticleWireInfo,0)
};
//----------------------------------------------------------------

//##################################################################
//##################################################################

//----------------------------------------------------------------
// do some book keeping of all wires + candidates in event
class HParticleWireManager : public TObject {

#define MAXWIRE 220

    Int_t wires[6][4][6][MAXWIRE]; //! wire address -> ntimes      sec,mod,lay,cell
    Int_t pid ;
    Int_t nWRange;             // +- wire to check in flagging


    Int_t nWires;              // number of wire in event
    Int_t nWiresUsed;          // number of wires used in event incl. doubles
    vector<HParticleWireInfo> vWireInfo;// vector of wireinfo for each HParticleCand
    vector<Int_t> vsegind;      // inner segment ind for kIsUsed particles
    map   <HParticleCand*,HParticleWireInfo> mCandWireInfo; // map HParticleCand -> WireInfo

    HCategory* pParticleCandCat;        //! HParticleCand category
    HCategory* pSegCat;                 //! HMdcSeg category
    HMdcSizesCells* sizescells;         //!
    HParticleMetaMatcher* matcher; //!


    Int_t ctCand;
    Int_t ctCandUsed;
    Int_t ctevent;
    Int_t ctSecCand[6];
    Int_t ctSecCandUsed[6];
    Int_t ctSecSegDouble[6];
    Int_t ctSecWire[6];
    Int_t ctSecWireDouble[6];
    Bool_t fDebug;          //! switch print on/off

protected:
    void  fillWireInfo(HParticleWireInfo& wi,Int_t flag);
    void  flagDoubleUsedWires();
    void  fillWire    (Int_t s,Int_t m,Int_t l,Int_t c,Int_t flag);  // count up a single cell
    void  clear();    // at start of event


public:
    HParticleWireManager();
    void  setDebug(Bool_t deb=kTRUE){fDebug = deb;}
    void  setWireRange(Int_t r)     { nWRange = r;}       // blocked +- range arround used wire in layer (0-n)
    Int_t getNWires    ()           { return nWires;}     // number of wires in event
    Int_t getNWiresUsed()           { return nWiresUsed;}
    Int_t isUsedNtimes(Int_t s,Int_t m,Int_t l,Int_t c);  // how many time this cell has been used ?
    vector<HParticleWireInfo>&  getWireInfo(){ return vWireInfo;}
    Bool_t getWireInfo(HParticleCand* c ,HParticleWireInfo& w);  // HParticleCand -> WireInfo
    Bool_t getWireInfo(UInt_t i,HParticleWireInfo& w,HParticleCand* c = 0);  // index of HParticleCand -> WireInfo
    Bool_t sharedWires(HParticleWireInfo& w1,HParticleWireInfo& w2     ,Int_t io=0,Int_t range = -1); // compare 2 WireInfo objects if they shared wires
    Bool_t sharedWires(HParticleCand* c1,HParticleCand* c2,Int_t io=0,Int_t range = -1); // compare 2 HParticleCand objects if they shared wires
    void   getWireInfoDirect (HParticleCand* cand,HParticleWireInfo& w,Bool_t fillStat=kFALSE);
    void   setMatcher(HParticleMetaMatcher* m) {matcher = m;}
    void  init();
    Int_t execute();

    void  print();

    ClassDef(HParticleWireManager,0)

};

//##################################################################
//##################################################################



//----------------------------------------------------------------
class HParticleMetaPredict : public TObject {

public:

    HParticleMetaPredict(){
	reset();
    }

    HParticleCand* c;
    Int_t system;
    Int_t nTof;
    Int_t nRpc;
    Int_t nMeta;
    Int_t addrTof [2];
    Int_t addrRpc [2];
    Int_t addrMeta[2];
    Double_t pathTof[2];
    Double_t pathRpc[2];
    Int_t    segInd [2];
    HGeomVector hitTof[2];
    HGeomVector hitRpc[2];
    HGeomVector meta;
    HGeomVector metaSec;

    vector<UInt_t> vcandTof; // index in output vector of candidates which share TOF
    vector<UInt_t> vcandRpc; // index in output vector of candidates which share RPC

    vector<Int_t> vPredmatchReal;  // linaddr of predicted found in real meta
    vector<Int_t> vPrednomatchReal;// linaddr of predicted not found in real meta
    vector<Int_t> vRealmatchPred;  // linaddr of real found in predicted meta
    vector<Int_t> vRealnomatchPred;// linaddr of real not found in predicted meta

    void reset(){
	c = 0;
	system = -1;
	nTof = 0;
	nRpc = 0;
	nMeta = 0;
	for(Int_t i = 0 ; i < 2; i++){
	    addrTof[i] =-1;
	    addrRpc[i] =-1;
	    addrMeta[i]=-1;
	    segInd  [i]=-1;
            pathTof[i] = 0;
            pathRpc[i] = 0;
	    hitTof[i].setXYZ(0,0,0);
	    hitRpc[i].setXYZ(0,0,0);
	}
	meta.setXYZ(0,0,0);
	metaSec.setXYZ(0,0,0);
        vcandTof.clear();
	vcandRpc.clear();

	vPredmatchReal  .clear();
	vPrednomatchReal.clear();
	vRealmatchPred  .clear();
	vRealnomatchPred.clear();
    }

    void set(HParticleCand* ca,Int_t sys,Int_t hit,Int_t s,Int_t m,Int_t ce,HGeomVector& vhit,Double_t path){

	if(ca == 0) return;
	if(hit < 0 || hit > 1) return;
	if(s < 0 || s > 5)     return;
	if(m == -1)            return;

	c = ca;
        segInd[0] = c->getInnerSegInd();
        segInd[1] = c->getOuterSegInd();

        nMeta = c ->getNMetaHits();

	if(sys == 1){
	    addrTof[hit] = convertAddr(s,m,ce);
	    hitTof [hit] = vhit;
            pathTof[hit] = path;
	    if(nTof < 2 ) nTof++;
	} else {
	    addrRpc[hit] = convertAddr(s,m,ce);
	    hitRpc [hit] = vhit;
            pathRpc[hit] = path;
	    if(nRpc < 2 ) nRpc++;
	}

	if(system == -1) system = sys;
	if(system == 0 && sys == 1 ) system = 2;
	if(system == 1 && sys == 0 ) system = 2;

	return;
    }

    void convertAddr(Int_t linadd,Int_t& s,Int_t&m,Int_t& ce){
	if(linadd < 0) { s = m = ce = -1; return;}
	Int_t d4,d5;
	HTool::getDimIndex(linadd,s,m,ce,d4,d5,6,10,100);
        return;
    }

    Int_t convertAddr(Int_t& s,Int_t&m,Int_t& ce){
        return HTool::getLinearIndex(s,6,m,10,ce,100);
    }

    Bool_t metaMatchPredict(){

	if(c->isRpcClstUsed()){
	    for(Int_t i = 0; i < 2; i++){
		if(addrMeta [i] == -1 ) continue;
		for(Int_t j = 0; j < nRpc; j++){
		    if(addrRpc[j] == addrMeta [i]) return kTRUE;
		}
	    }
	    return kFALSE;
	} else if (c->isTofClstUsed() || c->isTofHitUsed()){
	    for(Int_t i = 0; i < 2; i++){
		if(addrMeta [i] == -1 ) continue;
		for(Int_t j = 0; j < nTof; j++){
		    if(addrTof[j] == addrTof [i]) return kTRUE;
		}
	    }
	    return kFALSE;
	}
	return kFALSE;
    }

    void metaMatchPredictAll(){

	if(c->isRpcClstUsed()){
            // compare meta to pred
	    for(Int_t i = 0; i < 2; i++){
		if(addrMeta [i] == -1 ) continue;
		Bool_t found = kFALSE;
		for(Int_t j = 0; j < nRpc; j++){
		    if(addrRpc[j] == addrMeta [i]) vRealmatchPred.push_back(addrMeta[i]);
		}
                if(!found) vRealnomatchPred.push_back(addrMeta[i]);
	    }

            // compare pred to meta
	    for(Int_t i = 0; i < 2; i++){
		if(addrRpc [i] == -1 ) continue;
		Bool_t found = kFALSE;
		for(Int_t j = 0; j < 2; j++){
		    if(addrMeta [j] == -1 ) continue;
		    if(addrRpc[i] == addrMeta [j]) vPredmatchReal.push_back(addrRpc[i]);
		}
                if(!found) vPrednomatchReal.push_back(addrRpc[i]);
	    }
	    return;

	} else if (c->isTofClstUsed() || c->isTofHitUsed()){


            // compare meta to pred
	    for(Int_t i = 0; i < 2; i++){
		if(addrMeta [i] == -1 ) continue;
		Bool_t found = kFALSE;
		for(Int_t j = 0; j < nTof; j++){
		    if(addrTof[j] == addrMeta [i]) vRealmatchPred.push_back(addrMeta[i]);
		}
		if(!found) vRealnomatchPred.push_back(addrMeta[i]);
	    }

	    // compare pred to meta
	    for(Int_t i = 0; i < 2; i++){
		if(addrTof [i] == -1 ) continue;
		Bool_t found = kFALSE;
		for(Int_t j = 0; j < 2; j++){
		    if(addrMeta [j] == -1 ) continue;
		    if(addrTof[i] == addrMeta [j]) vPredmatchReal.push_back(addrTof[i]);
		}
		if(!found) vPrednomatchReal.push_back(addrTof[i]);
	    }

	    return ;
	}
	return ;
    }

    ClassDef(HParticleMetaPredict,0)

};


typedef  struct {

 Int_t sec;
 Int_t row;
 Int_t col;
 void set(Int_t s,Int_t r,Int_t c) {sec = s; row = r; col =  c;}

} richaddress;




class HParticleMetaMatcher : public HReconstructor {
    friend class HParticleWireManager;

protected:

    HCategory*       fCatEmcCluster;  //! pointer to the Emc cluster category
    HCategory*       fCatTofHit;      //! pointer to the Tof hit category
    HCategory*       fCatTofCluster;  //! pointer to the Tof cluster category
    HCategory*       fCatRpcCluster;  //! pointer to the Rpc cluster category.
    HCategory*       fCatParticleCand; //! HParticleCand category

    HRpcGeomPar*     fRpcGeometry;    //! Rpc geometry
    HRpcGeomCellPar* pGeomCellPar;    //! Rpc cell info
    HiTofGeomPar*    fiTofGeometry;   //! iTof geometry
    HEmcGeomPar*     fEmcGeometry;    //! Emc geometry
    HTofGeomPar*     fTofGeometry;    //! TOF geometry
    HSpecGeomPar*    pSpecGeomPar;    //! sector lab trans
    HRich700DigiPar* fRich700digipar; //! inverse look up of Rich700
    HMdcSizesCells*  fSizesCells;     //! helper container
    HMagnetPar*      fMagnetPar;      //! magnet current + polarity
    HMdcTrackGFieldPar* fFieldPar;    //! fieldmap

    HGeomTransform   labSecTrans[6];       //! [sec]
    HGeomTransform   modSecTransRpc[6];    //! [sec]
    HGeomTransform   modSecTransEmc[6];    //! [sec]
    HGeomTransform   modSecTransTof[6][8]; //! [sec][mod]
    HGeomTransform   cellSecTransiTof[6][3]; //! [sec][cell]

    TVector3         normVecRpc[6][3];   //! [sec][plane col{0,2,4}, {1,3,5}, mid]  normal vector on the RPC module in the sector coordinate system
    TVector3         normVecEmc[6];      //! [sec] normal vector on the Emc module in the sector coordinate system
    TVector3         normVecTof[6][8];   //! [sec][mod] normal vector on each Tof module in the sector coordinate system
    TVector3         normVeciTof[6];     //! [sec] normal vector on each iTof cell in the sector coordinate system

    TVector3         centerVecRpc[6][3]; //! [sec][plane col{0,2,4}, {1,3,5}, mid] center vector on the RPC module in the sector coordinate system
    TVector3         centerVecEmc[6];    //! [sec] center vector on the Emc module in the sector coordinate system
    TVector3         centerVecTof[6][8]; //! [sec][mod] center vector on each Tof module in the sector coordinate system
    TVector3         centerVeciTof[6]  ; //! [sec] center vector on each iTof cell in the sector coordinate system

    TVector3         centerVecRK[6];  //! [sec] point on plane of RK track definition
    TVector3         normVecRK  [6];  //! [sec] norm vec for plane of RK track definition
    TVector3         centerVecMdc[6][4]; //! [sec][mod] point on MDC MOD plane
    TVector3         normVecMdc  [6][4]; //! [sec][mod] norm  on MDC MOD plane
    Int_t            mdcInstalled[6][4]; //! [sec][mod] installed  =1 , not installed =0
    HParticleRKPlane planeMdc    [6][5]; //! [sec][start + 4mod]

    HGeomVector cellRPC[6][6][RPCMAXCELL][4];  //! [sec][col][cell]  4points
    HGeomVector cellEMC[6][EMCMAXCELL][4];     //! [sec][cell]       4points
    HGeomVector cellTOF[6][8][8][4];           //! [sec][mod][cell]  4points
    HGeomVector celliTOF[6][3][4];             //! [sec][cell][cell]  4points




    Float_t DPlanesRpc;        //! RPC displacement for collumns

    //------------------------------------
    // return fake objects for recalcXXX() functions
    HRpcCluster rpccluster;   //!
    HTofCluster tofcluster;   //!
    HEmcCluster emccluster;   //!

    //------------------------------------
    // MDC hit on meta calculated by  recalcXXX() functions
    HGeomVector  trackMetaMod; //!
    HGeomVector  trackMetaSec; //!
    HGeomVector  trackMetaLab; //!

    Bool_t fDebug;          //! switch print on/off
    Bool_t frunWireManager; //! default true  : fill ManagerContent in execute()
    Bool_t frunFlagMETADouble; //! default false : flag double used META CELL by canddidates pointing to the cell (even if they have no match META hit)
    Bool_t fuseEMC;         //! default true
    Bool_t fuseRICH700;     //! default true
    Bool_t fuseiTOF;        //! default false
    Bool_t fuseRK;          //! default true : use RK to estimate MDC hits (otherwise sraight line) from HParticleCand
    Bool_t fuseSeg;         //! default false : use MDC segs from input if available to fill wire statistics
    Bool_t fInitOK;         //! check if reinit has been ok
    HParticleWireManager  wiremanager; // store used wires in event
    HParticleRK          rkpropagator; // propagate particles through the detector

    Bool_t traceToMeta(HParticleCand* cand,TVector3 &metaHit, const TVector3 &metaCenter, const TVector3 &metaNorm,Double_t& path);
public:

    HParticleMetaMatcher(const Text_t *name = "HParticleMetaMatcher",const Text_t *title ="HParticleMetaMatcher");

    virtual ~HParticleMetaMatcher() { ; }
    void     setDebug         (Bool_t debug=kTRUE) {fDebug = debug;}
    void     setRunWireManager(Bool_t run=kTRUE)   {frunWireManager = run;}
    void     setRunFlagMetaDoubles(Bool_t run=kTRUE)   {frunFlagMETADouble = run;}
    void     setUseEMC(Bool_t use=kTRUE)           {fuseEMC = use;}
    void     setUseiTof(Bool_t use=kTRUE)          {fuseiTOF = use;}
    void     setUseRich700(Bool_t use=kTRUE)       {fuseRICH700 = use;}
    void     setUseRK (Bool_t use=kTRUE)           {fuseRK = use;}
    void     setUseSeg(Bool_t use=kTRUE)           {fuseSeg = use;}
    HParticleWireManager&  getWireManager()  { return wiremanager;}
    HParticleRK&           getRKPropagator() { return rkpropagator;}


    const HRpcCluster*  recalcRpc(HParticleCand* cand);
    const HTofCluster*  recalcTof(HParticleCand* cand);
    const HEmcCluster*  recalcEmc(HParticleCand* cand);

    const HGeomVector&  getTrackMetaMod() { return trackMetaMod;}  // only valid direct following recalcXXX() call!
    const HGeomVector&  getTrackMetaSec() { return trackMetaSec;}  // only valid direct following recalcXXX() call!
    const HGeomVector&  getTrackMetaLab() { return trackMetaLab;}  // only valid direct following recalcXXX() call!

    Bool_t              predictRpcCell(HParticleCand* cand,HGeomVector& hit1,HGeomVector& hit2,Int_t& s,Int_t& col1,Int_t& cell1,Int_t& col2,Int_t& cell2,Double_t& path1,Double_t& path2,Int_t sys=1);
    Bool_t              predictTofCell(HParticleCand* cand,HGeomVector& hit0,HGeomVector& hit1,Int_t& s,Int_t& mod0,Int_t& cell0,Int_t& mod1,Int_t& cell1,Double_t& path1,Double_t& path2,Int_t sys=1);
    Bool_t              predictEmcCell(HParticleCand* cand,HGeomVector& hit,Int_t& s,Int_t& pos,Int_t& cell,Double_t& path,Int_t sys=1);
    Bool_t              predictiTofCell(HParticleCand* cand,HGeomVector& hit,Int_t& s, Int_t& cell, Double_t& path,Int_t sys=1);
    Bool_t              predictRich700Hit(HParticleCand* cand,Double_t zvertex, Int_t& s ,Int_t& col,Int_t& row, HGeomVector& hit, vector<richaddress>& vrichcaladd, Int_t range = 10);
    Bool_t              predictRich700Hit(HParticleCand* cand,Double_t zvertex, Int_t& s ,Int_t& col,Int_t& row, HGeomVector& hit, vector<HRichCal*>& vrichcal, Int_t range = 10);
    Bool_t              predictMdcHit   (HParticleCand* cand,HGeomVector& hit,Int_t m,Double_t& path);
    Bool_t              predictMdcHitSec(HParticleCand* cand,HGeomVector& hit,Int_t m,Double_t& path);
    Bool_t              predictMdcHitLab(HParticleCand* cand,HGeomVector& hit,Int_t m,Double_t& path);

    Bool_t              calcLineRpcCell(const Int_t s,const TVectorD& state,HGeomVector& hit1,HGeomVector& hit2,Int_t& col1,Int_t& cell1,Int_t& col2,Int_t& cell2,Double_t& path1,Double_t& path2,Int_t sys=1);
    Bool_t              calcLineTofCell(const Int_t s,const TVectorD& state,HGeomVector& hit0,HGeomVector& hit1,Int_t& mod0,Int_t& cell0,Int_t& mod1,Int_t& cell1,Double_t& path1,Double_t& path2,Int_t sys=1);
    Bool_t              calcLineEmcCell(const Int_t s,const TVectorD& state,HGeomVector& hit,Int_t& pos,Int_t& cell,Double_t& path,Int_t sys=1);
    Bool_t              calcLineMdcHit (const Int_t s,const TVectorD& state,HGeomVector& hit,const Int_t m,Double_t& path,Int_t sys=1);
    Bool_t              calcLineiTofCell(const Int_t s,const TVectorD& state,HGeomVector& hit,Int_t& cell,Double_t& path,Int_t sys=1);

    Bool_t              isInRpcCell(HGeomVector& vmod,Int_t s,Int_t col,Int_t c,Float_t xRes = 0.,Float_t yRes = 0.);
    Bool_t              isInTofCell(HGeomVector& vmod,Int_t s,Int_t col,Int_t c,Float_t xRes = 0.,Float_t yRes = 0.);
    Bool_t              isInEmcCell(HGeomVector& vmod,Int_t s,Int_t c,Float_t xRes = 0.,Float_t yRes = 0.);
    Bool_t              isIniTofCell(HGeomVector& vmod,Int_t s,Int_t c,Float_t xRes = 0.,Float_t yRes = 0.);
    void                distToEdgeTofCell(HGeomVector& vmod,Int_t s,Int_t m, Int_t c, Float_t& distLeft,Float_t& distRight,Float_t& distLow,Float_t& distUp);
    void                distToEdgeRpcCell(HGeomVector& vmod,Int_t s,Int_t m, Int_t c, Float_t& distLeft,Float_t& distRight,Float_t& distLow,Float_t& distUp);
    void                distToEdgeEmcCell(HGeomVector& vmod,Int_t s, Int_t c, Float_t& distLeft,Float_t& distRight,Float_t& distLow,Float_t& distUp);
    void                distToEdgeiTofCell(HGeomVector& vcell,Int_t s,Int_t c, Float_t& distLeft,Float_t& distRight,Float_t& distLow,Float_t& distUp);

    const HGeomVector*  getRpcCellGeom(Int_t s,Int_t col,Int_t cell);
    const HGeomVector*  getTofCellGeom(Int_t s,Int_t mod,Int_t cell);
    const HGeomVector*  getEmcCellGeom(Int_t s,Int_t cell);
    const HGeomVector*  getiTofCellGeom(Int_t s,Int_t cell);
    void                getEmcCellArray(Int_t s,Int_t cell,vector<Int_t>& vcells);
    void                getWireInfoDirect(HParticleCand* cand,HParticleWireInfo& w,Bool_t fillStat=kFALSE);

    Bool_t              getTrackStateRK(HParticleCand* cand,TVector3& pos,TVector3& momv,Double_t mom=0.,Int_t charge=0);
    Bool_t              getTrackStateRK(HParticleCand* cand,TVectorD& sv,Double_t mom=0.,Int_t charge=0);
    Bool_t              propagateToPlane(Int_t chrg,TVectorD& stateTo, const TVectorD& stateFrom,const HParticleRKPlane& planeTo,Double_t& pathL,Bool_t direction=kIterForward);
    Bool_t              propagateToPCA  (Int_t chrg,TVectorD& stateTo, const TVectorD& stateFrom,const TVector3& point,Double_t& pathL,Bool_t direction=kIterBackward);
    Bool_t              getMDCStateRK  (HParticleCand* cand,TVectorD& st0 ,TVectorD& st1 ,TVectorD& st2 ,TVectorD& st3 ,TVectorD& pathlength,Double_t mom=0.,Int_t charge=0);
    Bool_t              getMDCPosRK    (HParticleCand* cand,TVector3& hit0,TVector3& hit1,TVector3& hit2,TVector3& hit3,TVectorD& pathlength,Double_t mom=0.,Int_t charge=0);
    Bool_t              getMDCPosMomRK (HParticleCand* cand,TVector3& hit0,TVector3& hit1,TVector3& hit2,TVector3& hit3,TVector3& mom0,TVector3& mom1,TVector3& mom2,TVector3& mom3,TVectorD& pathlength,Double_t mom=0.,Int_t charge=0);

    HParticleRKPlane    getPlaneRK (Int_t s);
    HParticleRKPlane    getPlaneMDC(Int_t s,Int_t col);
    HParticleRKPlane    getPlaneTOF(Int_t s,Int_t m);
    HParticleRKPlane    getPlaneRPC(Int_t s,Int_t m);
    HParticleRKPlane    getPlaneEMC(Int_t s);
    HParticleRKPlane    getPlaneiTOF(Int_t s);
    static Bool_t       calcState  (Int_t& s,TVectorD& state, const TVector3& pos,const Double_t theta, const Double_t phi, const Double_t mom);
    static Bool_t       kineToParticleCand(HGeantKine* kine,HParticleCandSim& cand, Bool_t realmom=kFALSE);
    static Bool_t       kineToToState     (HGeantKine* kine,TVectorD& state,Int_t& charge,Bool_t realmom=kFALSE,const Int_t forceSector=-1);
    static Bool_t       stateLabToSecSys(const Int_t s, const TVectorD&    stateLab,TVectorD& stateSec,Bool_t to=kTRUE);
    static Bool_t       pointLabToSecSys(const Int_t s, const TVector3&    pointLab,TVector3& pointSec,Bool_t to=kTRUE);
    static Bool_t       pointLabToSecSys(const Int_t s, const HGeomVector& pointLab,TVector3& pointSec,Bool_t to=kTRUE);
    Bool_t              getTargetMiddlePoint(const Int_t s,HGeomVector& midpoint);
    Bool_t              getDistanceToTargetMiddlePoint(const Int_t s,const TVectorD& stateRK, Double_t& dist);
    Bool_t              propagateToVertexStateRK   (HParticleCand* cand,TVectorD& stStart ,TVectorD& stPCA  ,const TVector3& vertexLab,Double_t& pathlength,Double_t mom=0.,Int_t charge=0);
    Bool_t              propagateToVertexPosRK     (HParticleCand* cand,TVector3& posStart,TVector3& posPCA ,const TVector3& vertexLab,Double_t& pathlength,Double_t mom=0.,Int_t charge=0);
    Bool_t              propagateToVertexPosMomRK  (HParticleCand* cand,TVector3& posStart,TVector3& posPCA,TVector3& momStart,TVector3& momPCA,const TVector3& vertexLab,Double_t& pathlength,Double_t mom=0.,Int_t charge=0);
    Bool_t              propagateToVertexStateLine (HParticleCand* cand,TVectorD& stStart ,TVectorD& stPCA ,const TVector3& vertexLab,Double_t& pathlength);
    Bool_t              propagateToVertexPosLine   (HParticleCand* cand,TVector3& posStart,TVector3& posPCA,const TVector3& vertexLab,Double_t& pathlength);
    Bool_t              propagateToVertexPosMomLine(HParticleCand* cand,TVector3& posStart,TVector3& posPCA,TVector3& momStart,TVector3& momPCA,const TVector3& vertexLab,Double_t& pathlength);

    Bool_t              printRichCals(vector<HRichCal*>& vrichcal, Int_t col,Int_t row, Int_t range);
    Bool_t              loopFlagDoubleUsedMetaCells(vector<HParticleMetaPredict>& vpredict,UInt_t print = 0);
    HiTofGeomPar*       getiTofGeomPar()    { return fiTofGeometry; }
    HRich700DigiPar*    getRich700DigiPar() { return fRich700digipar; }

    Bool_t  init    ();
    Bool_t  reinit  ();
    Int_t   execute ();
    Bool_t  finalize() { return kTRUE; }

    ClassDef(HParticleMetaMatcher,0)
};

#endif // HParticleMetaMatcher_h