#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 "hitofgeompar.h"
#include "htofgeompar.h"
#include "hrich700data.h"
#include "hrich700digipar.h"
#include "hspecgeompar.h"
#include "hmdcsizescells.h"
#include "hmagnetpar.h"
#include "hmdctrackgfieldpar.h"
#include "hkalgeomtools.h"

#include "hparticlecand.h"
#include "hparticlecandsim.h"
#include "hgeantkine.h"
#include "hparticletool.h"
#include "hmdcseg.h"
#include "hrichcal.h"

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

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

#include <algorithm>
#include <math.h> 

//_HADES_CLASS_DESCRIPTION
/////////////////////////////////////////////////////////////
//
//  HParticleMetamatcher
//
//  Helper class for analysis. Derived from HReconstructer, should be added to tasklist. Recalculation
//  of detector hits etc needs proper params and correct init(), reinit() of the task.
//
//  Typical applications:
//
//  Main configuration:
//  void     setRunWireManager(Bool_t run=kTRUE)   {frunWireManager = run;} : kTRUE (default)  : fill wire information for full event
//  void     setUseEMC(Bool_t use=kTRUE)           {fuseEMC = use;}         : kTRUE (default)  : init EMC detector and params. set KFALSE for apr12 beamtime (no EMC available)
//  void     setUseiTof(Bool_t use=kTRUE)          {fuseiTOF = use;}        : kFALSE (default) : init iTof detector and params. set kTRUE for example feb22
//  void     setUseRK (Bool_t use=kTRUE)           {fuseRK = use;}          : kTRUE (default)  : use RK for track propagation. kFALSE : use straight line propagation from store HParticleCand inner/outer segs
//                                                                            This setting has effect on predictXXX() functions and filling of wireinfo,
//                                                                            recalcXXX() functions
//  void     setUseSeg(Bool_t use=kTRUE)           {fuseSeg = use;}         : kFALSE (default) : do not try to read HMdcSeg cateory to fill wire infos
//
//
//
//  1. recalculation of HTofCluster, HRpcCluster, HEmcCluster (if they are not stored in root file)
//     from HParticleCand using the store meta qualtity, tof etc and meta cell,mod
//     using RK propagation (setUseRK(kTRUE)(default) or straight line approach setUseRK(kFALSE)
//     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!
//
//  2. prediction of detector hits using Rungekutta propagation through the field (default) or RK segment info stored
//     in HParticleCand. Options for different coord. sys of the hits (MOD,SEC,LAB)
//     setUseRK (Bool_t use=kTRUE)  kTRUE (default) : use RK propagation.works for all dsts since only inner RK segments are needed
//                                                    and the outer hits are recalculated
//                                  kFALSE          : use RK segment info stored in HParticleCand (will not work in older dst where outer seg info was not stored)
//
//     The functions provide path length from the definition plane of RK (500 mm parallel in front of MDCI or 650 mm for MDCII)
//     to the hit if RK is used. In case  mom of cand < 0 (no outer seg ) straight lines are used for calculation of hit points and
//     path = 0.
//
//     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  predictMdcHit (HParticleCand* cand,HGeomVector& hit,Int_t m,Double_t& path);
//     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);
//     Bool_t  predictRich700Hit(HParticleCand* cand,Double_t zvertex,Int_t& s, Int_t& col,Int_t& row, HGeomVector& hit, vector<HRichCal*>&   vrichcaladd,Int_t range);
//     .....
//
//  3. propagation through detector
//
//     done in sec. coord. sys.  particle track state  (x,y,z,px,py,pz  [mm, MeV/c]) using RK propagation (not dependend on the values setUseRK (Bool_t use=kTRUE) )
//     is calculated from HParticleCand. The optional arguments for momentum and charge can be used
//     to propagated candidates without mom or chrg, when out segment was missing or one wants to check influence
//     of mom and charge to the reconstructed hits. mom and charge replace the values in HParticleCand when not equal 0
//     The functions provide path length from the definition plane of RK (500 mm parallel in front of MDCI or 650 mm for MDCII)
//     to the hit.
//
//     Bool_t getTrackStateRK(HParticleCand* cand,TVectorD& sv,Double_t mom=0.,Int_t charge=0);
//     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);
//
//     For neutral particles / field off the hit points on detectors can be efficiently calculated using
//
//     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);
//
//
//  4. propagation to point of closest approach to vertex using Rungekutta or straight line approach
//     done in sec. coord. sys.  particle track state  (x,y,z,px,py,pz  [mm, MeV/c]), path length is provided too.
//     For convenience the vertex input is in LAB sys (as it is stored in EventHeader) and transformed internally to SEC sys
//     The functions provide path length from the definition plane of RK (500 mm parallel in front of MDCI or 650 mm for MDCII)
//     to the point of closest approach (PCA) to the provided vertex.
//
//     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);
//
//     to transform hits, states between LAB <-> SEC sys  use (default to =kTRUE LAB->SEC, kFALSE : SEC->LAB):
//     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);
//
//     Helper functions to work with trackstates, HParticleCand and HGeantKine
//
//    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);
//
//
//  5. retrieving MDC wire info for an HParticleCand.
//     a. REAL WIRES        : from HMdcSeg if the category is available in the input and setUseSeg(kTRUE) is set.
//     b. THEORETICAL WIRES : using RK propagation of HParticleCand to estimation hit points of MDCs and calculation
//                            theoretical fired wires by the define straight line of the 2 segment points
//     HParticleWireInfo objects keep the infos about fired wires of a candidate.
//     The object provied functions to get Number of wires and layer in wireinfo, as well a print functions
//     The wireinfo for each candidate can be retrieved in 2 ways:
//     a. filling for all candidates the full event in the execute() if  setRunWireManager(kTRUE)  : kTRUE (default)
//        In this case a statistic is filled for all candidates and information about double used wires
//        in the event is filled too.
//        HParticleWireManager& wiremanager = matcher->getWireManager();
//        manager.sharedWires(HParticleCand* c1,HParticleCand* c2,Int_t io=0,Int_t range = -1); // compare 2 HParticleCand objects if they shared wires
//        wiremanager.getWireInfo(HParticleCand* c ,HParticleWireInfo& w);
//        Int_t wiremanager.getNWires    ()  // number of wires in event (double used wires are counted only once)
//        Int_t wiremanager.getNWiresUsed()  // number of wires in event (double used wires are counted multiple times)
//        Int_t wiremanager.isUsedNtimes(Int_t s,Int_t m,Int_t l,Int_t c);  // how many time this cell has been used ?
//
//     b. filling  per track candidate if information about double used wires is not needed and calculation time
//        should be saved. setRunWireManager(kFALSE)
//        matcher->getWireInfoDirect(HParticleCand* cand,HParticleWireInfo& w,Bool_t fillStat=kFALSE);
//        Int_t  w.nSharedWires(HParticleWireInfo& w2     ,Int_t io=0,Int_t range = -1); // get number of shared wires
//        wiremanager.sharedWires(HParticleCand* c1,HParticleCand* c2,Int_t io,Int_t range) // Bool_t only, faster than above
//        w.printSharedWires(HParticleWireInfo& w2);
//
//
//  6. flagging HParticleCand for use of the same META cells using the predictXXXCell() functions.
//     c->setFlagBit(Particle::kIsMetaDoubleUse); in HLoop macro matcher->loopFlagDoubleUsedMetaCells(vector<HParticleMetaPredict>& vpredict)
//     can be called directly or adding HParticleMetaMatcher to the tasklist (setting setRunFlagMetaDoubles(kTRUE)). This function
//     overwrite the flag in HParticleCand.
//-----------------------------------------------------------
// 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
//     //matcher->setUseEMC(kFALSE);         // default kTRUE : before mar19 EMC was not available
//     //matcher->setUseiTof(kFALSE);        // default kFALSE : available for feb22
//     //matcher->setRunFlagMetaDoubles(kTRUE); // if flag function should be called automaticcally in the tasklist
//     //matcher->setRunWireManager(kFALSE); // default kTRUE : do not fill wire info for full event
//     //matcher->setUseSeg(kTRUE); // default kFALSE : use Segment from Input to fill wirefinfo
//     //matcher->setUseRK(kTRUE);  // default kTRUE : use RK propagator to calculate hits points in MDC and META, other wise straight lines from HParticelCand seg
//     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)
// 	{
//          //---------------------------------------------------------------------------
//          // flagging HParticleCand for use of the same META cells (using the prredictXXXCell() functions)
//          // When used here matcher->setRunFlagMetaDoubles(kFALSE) should be set to not exectute it from the tasklist
//          // if matcher has been added. The output vector can be used for further investigation
//          // c->isFlagBit(Particle::kIsMetaDoubleUse);
//          vector<HParticleMetaPredict> vpredict;
//          matcher->loopFlagDoubleUsedMetaCells(vpredict);
//          //---------------------------------------------------------------------------
//
//
// 	    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;
//                      Double_t path1,path2;
// 			HGeomVector hit0,hit1;
// 			matcher->predictRpcCell(cand,hit0,hit1,s,col0,cell0,col1,cell1,path1,path2); // 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;
//                      Double_t path1,path2;
// 			HGeomVector hit0,hit1;
// 			matcher->predictTofCell(cand,hit0,hit1,s,mod0,cell0,mod1,cell1,path1,path2); // which Tof cell should this track fire?
// 		    }
// 		    if(0&&cand->getEmcInd()>-1){
// 			const HEmcCluster* emc = matcher->recalcEmc(cand);
// 			Int_t s,pos,cell;
//                      Double_t path;
// 			HGeomVector hit;
// 			vector<Int_t> vcells;
// 			matcher->predictEmcCell(cand,hit,s,pos,cell,path); // 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){
//                      Double_t path = 0;
// 			HGeomVector hitmdc;
// 			HGeomVector hitmdcsec;
// 			HGeomVector hitmdclab;
// 			Int_t module = 0; // 0-3
// 			matcher->predictMdcHit   (cand,hitmdc   ,module,path);  // MDC coord system
// 			matcher->predictMdcHitSec(cand,hitmdcsec,module,path);  // sector coord system
// 			matcher->predictMdcHitLab(cand,hitmdclab,module,path);  // 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(HParticleMetaPredict)
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;
    vsegind.clear();       // 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)
	{
            HParticleWireInfo w;
	    getWireInfoDirect(c,w,kTRUE);
	}
    }
    flagDoubleUsedWires();  // set isGood flag per candidate
    //-------------------------------------------------------------------------


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

    return 0;
}


void   HParticleWireManager::getWireInfoDirect(HParticleCand* c,HParticleWireInfo& w,Bool_t fillStat)
{
    // Fill  HParticleWireInfo from HParticleCand
    // fillStat = kTRUE            :  fill internal wire statistics of HParticleWireManager  (used in execute())
    // fuseSeg  = kTRUE  (default kFALSE) wire info will be filled from HMdcSeg Input
    //          = kFALSE                  wire info will be filled from recalculated RK Hit points

    w.reset();

    if(matcher->fuseSeg&&!pSegCat) {
	Error("getWireInfoSeg()","catMdcSeg not available !");
        return;
    }

    if(!c) {
	Error("getWireInfoSeg()","HParticleCand*  is NULL !");
        return;
    }
    sizescells          = HMdcSizesCells::getExObject();
    if(!sizescells) return;
    //---------------------------------------------
    // calculate expected hits

    HMdcSeg* pSeg[2] = { 0 ,0};
    TVector3 hit0,hit1,hit2,hit3;
    if(matcher->fuseSeg){
	pSeg[0] = HParticleTool::getMdcSeg(c->getInnerSegInd());
	pSeg[1] = HParticleTool::getMdcSeg(c->getOuterSegInd());
    } else {
        TVectorD pathlength;
	if(!matcher->getMDCPosRK(c,hit0,hit1,hit2,hit3,pathlength)) return ;
    }

    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

    }

    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(matcher->fuseSeg){
		//---------------------------------------
		// 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<2) good =(*sizescells)[sec][mod][lay].calcCrossedCells(hit0.X(),hit0.Y(),hit0.Z(),hit1.X(),hit1.Y(),hit1.Z(),cell1,cell2);
		else      good =(*sizescells)[sec][mod][lay].calcCrossedCells(hit2.X(),hit2.Y(),hit2.Z(),hit3.X(),hit3.Y(),hit3.Z(),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;
		if(fillStat)fillWire(sec,mod,lay,c1,flag);
	    }
	    if(c2 > -1 && c2 != c1) {
		w.ar[mod][lay][1] = c2;
		if(fillStat)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;
    if(fillStat) {
	fillWireInfo(w,flag);
	Int_t ind = vWireInfo.size() -1;
	mCandWireInfo[c] = vWireInfo[ind];
    }
    return;
}

Bool_t HParticleWireManager::getWireInfo(HParticleCand* c ,HParticleWireInfo& w)
{
    // get wireinfo object by candidate pointer
    if(!c) return kFALSE;

    map<HParticleCand*,HParticleWireInfo>::iterator it;
    it = mCandWireInfo.find(c);
    if (it != mCandWireInfo.end()) w = it->second ;
    else return kFALSE;

    return kTRUE;
}

Bool_t HParticleWireManager::getWireInfo(UInt_t i,HParticleWireInfo& w,HParticleCand* c)
{
    // get wireinfo object by index (checking if it match candidate)

    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 ,

    Int_t nWR = nWRange;
    if(range !=-1) { nWR = abs(range);}

    Int_t nSh = w1.nSharedWires(w2,io,nWR,kTRUE); // quit Early by first match

    if(nSh > 0) return kTRUE;
    else        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;
    frunWireManager = kTRUE;
    frunFlagMETADouble = kFALSE;
    fuseEMC         = kTRUE;
    fuseiTOF        = kFALSE;
    fuseRK          = kTRUE;
    fuseSeg         = kFALSE;
    fInitOK         = kFALSE;
    wiremanager.setMatcher(this);

    memset(mdcInstalled,0,sizeof(Int_t)*24);

}


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

const HGeomVector* HParticleMetaMatcher::getiTofCellGeom(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 > iTOFMAXCELL-1) return 0;
    return &celliTOF[s][cell][0];
}

Bool_t HParticleMetaMatcher::traceToMeta(HParticleCand* cand,TVector3 &metaHit, const TVector3 &metaCenter, const TVector3 &metaNorm, Double_t& path)
{
    // calculate intersection of outer MdcSeg straightline with META plane
    // or RK propagation (fuseRK). RK provides path length too
    // in sector sys

    path = 0.;

    if(!fuseRK || cand->getMomentum() < 0){
	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 HKalGeomTools::findIntersectionLinePlane(metaHit,p1,dir,metaCenter,metaNorm);

    } else { // use RK propagation

	TVectorD stateFrom(6);
	if(!getTrackStateRK(cand,stateFrom)) return kFALSE;

	Int_t chrg = cand->getCharge();

	TVectorD stateTo(6);  // transported state at the destination layer
        HParticleRKPlane planeMeta(metaCenter,metaNorm);

	rkpropagator.propagateToPlane(chrg,stateTo,stateFrom,planeMeta,path, kIterForward);

	metaHit.SetXYZ(stateTo(0),stateTo(1),stateTo(2));

        return kTRUE;
    }
}

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 * (( vmod.Y() - p0.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 * (( vmod.Y() - p3.Y())/b);

    if (x1 > vmod.X() + xRes) {
	if(fDebug){
	    cout<<"isInRpcCell  s "<<s<<" col "<<col<<" c "<<setw(2)<<c
		<<" 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;
}

void HParticleMetaMatcher::distToEdgeRpcCell(HGeomVector& vmod,Int_t s,Int_t col, Int_t c, Float_t& distLeft,Float_t& distRight,Float_t& distLow,Float_t& distUp)
{
    //  check if vmodule (mod system) is inside cell (2D) Trapez
    //
    //       p1---------------p2
    //      b\                 /
    //        \               /
    //        p0-------------p3
    //       a
    //
    //  calulates distance to edge of cell
    //  negative values -> outside cell

    HGeomVector& p0 = cellRPC[s][col][c][0];
    HGeomVector& p1 = cellRPC[s][col][c][1];
    HGeomVector& p2 = cellRPC[s][col][c][2];
    HGeomVector& p3 = cellRPC[s][col][c][3];

    distLeft = distRight = distLow = distUp = 0;


    distUp    = p2.Y()   - vmod.Y();
    distLow   = vmod.Y() - p0.Y();

    distLeft  = p0.X()   - vmod.X();
    distRight = vmod.X() - p2.X();

    // 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 * (( vmod.Y() - p0.Y())/b);

    distLeft   = fabs(x - vmod.X());
    if(vmod.X() > x) distLeft *=-1;

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

    distRight   = fabs(x1 - vmod.X());
    if(vmod.X() < x1) distRight *=-1;

    return;
}


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)

    if(!fuseEMC || !fEmcGeometry) return kFALSE;
    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;
}

void HParticleMetaMatcher::distToEdgeEmcCell(HGeomVector& vmod,Int_t s, Int_t c, Float_t& distLeft,Float_t& distRight,Float_t& distLow,Float_t& distUp)
{

    //  distance (mod system) to  cell (2D) square
    //
    //        p1----p2
    //        |     |
    //        |     |
    //        p0----p3
    //
    // negative values -> outside cell

    if(!fuseEMC || !fEmcGeometry) return ;
    Int_t pos = HEmcDetector::getPositionFromCell(c);

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

    distLeft = distRight = distLow = distUp = 0;

    distLeft  = fabs(p0.X()   - vmod.X());
    if(vmod.X() > p0.X()  ) distLeft *= -1;
    distRight = fabs(vmod.X() - p2.X());
    if(vmod.X() < p2.X()  ) distRight *= -1;

    distUp    = fabs(p2.Y()   - vmod.Y());
    if(vmod.Y() > p2.Y()  ) distUp *= -1;
    distLow   = fabs(vmod.Y() - p0.Y());
    if(vmod.Y() < p0.Y()  ) distLow *= -1;

    return;
}

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::distToEdgeTofCell(HGeomVector& vmod,Int_t s,Int_t m, Int_t c, Float_t& distLeft,Float_t& distRight,Float_t& distLow,Float_t& distUp)
{

    //  distance (mod system) to  cell (2D) square
    //
    //        p1----p2
    //        |     |
    //        |     |
    //        p0----p3
    //
    // negative values -> outside cell

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

    distLeft = distRight = distLow = distUp = 0;

    distLeft  = fabs(p0.X()   - vmod.X());
    if(vmod.X() > p0.X()) distLeft *= -1;
    distRight = fabs(vmod.X() - p2.X());
    if(vmod.X() < p2.X()) distRight *= -1;

    distUp    = fabs(p2.Y()   - vmod.Y());
    if(vmod.Y() > p2.Y()) distUp *= -1;
    distLow   = fabs(vmod.Y() - p0.Y());
    if(vmod.Y() < p0.Y()) distLow *= -1;

    return;
}

Bool_t HParticleMetaMatcher::isIniTofCell(HGeomVector& vcell,Int_t s,Int_t c, Float_t xRes, Float_t yRes)
{
    //  check if vcell (cell system) is inside cell (2D) Trapez
    //
    //       p1---------------p2
    //      b\                 /
    //        \               /
    //        p0-------------p3
    //              a
    //
    // checks if HGeomVector* vmodule
    if(!fiTofGeometry) return kFALSE;

    if(s < 0 || s > 5 || c < 0 || c > iTOFMAXCELL - 1){
	Error("isInCell()","address out of bounds ! s = %d, c = %d !",s,c);
        return kFALSE;
    }


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

    if (p0.Y() > vcell.Y() + yRes){
	if(fDebug){
	    cout<<"isInCell  s "<<" 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)<<vcell.X()<<" , "<<setw(10)<<vcell.Y()<<" miss : hit below cell"<<endl;
	}
	return kFALSE;   // first check y range
    }

    if (p1.Y() < vcell.Y() - yRes) {
	if(fDebug){
	    cout<<"isInCell  s "<<s<<" 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)<<vcell.X()<<" , "<<setw(10)<<vcell.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 * (( vcell.Y() - p0.Y())/b);

    if (x < vcell.X() - xRes) {
	if(fDebug){
	    cout<<"isInCell  s "<<s<<" 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)<<vcell.X()<<" , "<<setw(10)<<vcell.Y()<<" x cut "<<x<<" miss : x hit out of cell left side"<<endl;
	}
	return kFALSE;
    }

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

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

    if (x1 > vcell.X() + xRes) {
	if(fDebug){
	    cout<<"isInCell  s "<<s<<" 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)<<vcell.X()<<" , "<<setw(10)<<vcell.Y()<<" x cut "<<x1<<", miss : x hit out of cell right side"<<endl;
	}
	return kFALSE;
    }
    if(fDebug){
	cout<<"isInCell  s "<<s<<" 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)<<vcell.X()<<" , "<<setw(10)<<vcell.Y()<<" calc x "<<setw(10)<<x<<" "<<x1<<" hit matched"
	    <<endl;
    }

    return kTRUE;
}

void HParticleMetaMatcher::distToEdgeiTofCell(HGeomVector& vcell,Int_t s, Int_t c, Float_t& distLeft,Float_t& distRight,Float_t& distLow,Float_t& distUp)
{
    //  check if vcell (cell system) is inside cell (2D) Trapez
    //
    //       p1---------------p2
    //      b\                 /
    //        \               /
    //        p0-------------p3
    //       a
    //
    //  calulates distance to edge of cell
    //  negative values -> outside cell

    HGeomVector& p0 = celliTOF[s][c][0];
    HGeomVector& p1 = celliTOF[s][c][1];
    HGeomVector& p2 = celliTOF[s][c][2];
    HGeomVector& p3 = celliTOF[s][c][3];

    distLeft = distRight = distLow = distUp = 0;


    distUp    = p2.Y()   - vcell.Y();
    distLow   = vcell.Y() - p0.Y();

    distLeft  = p0.X()   - vcell.X();
    distRight = vcell.X() - p2.X();

    // 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 * (( vcell.Y() - p0.Y())/b);

    distLeft   = fabs(x - vcell.X());
    if(vcell.X() > x) distLeft *=-1;

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

    distRight   = fabs(x1 - vcell.X());
    if(vcell.X() < x1) distRight *=-1;

    return;
}

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();
    if(!fuseEMC || !fEmcGeometry) return;

    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);
    fCatParticleCand= HCategoryManager::getCategory(catParticleCand);

    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.");
    if(!fCatParticleCand)Warning("init()", "Category catParticleCand not found.");


    fRpcGeometry = (HRpcGeomPar*)    rtdb->getContainer("RpcGeomPar");
    pGeomCellPar = (HRpcGeomCellPar*)rtdb->getContainer("RpcGeomCellPar");
    pSpecGeomPar = (HSpecGeomPar*)   rtdb->getContainer("SpecGeomPar");

    if (gHades->getSetup()->getDetector("Emc")) {
	if(fuseEMC) { fEmcGeometry = (HEmcGeomPar *) rtdb->getContainer("EmcGeomPar"); }
	else        { fEmcGeometry = 0; fuseEMC = 0; }
    } else  { fEmcGeometry = 0; fuseEMC = 0; }

    if (gHades->getSetup()->getDetector("iTof")) {
	if(fuseiTOF) { fiTofGeometry = (HiTofGeomPar *) rtdb->getContainer("iTofGeomPar"); }
	else         { fiTofGeometry = 0; fuseiTOF = 0; }
    } else  { fiTofGeometry = 0; fuseiTOF = 0; }

    if (gHades->getSetup()->getDetector("Rich")) {
	if(fuseRICH700) { fRich700digipar = (HRich700DigiPar *) rtdb->getContainer("Rich700DigiPar"); }
	else            { fRich700digipar = 0; fuseRICH700 = 0; }
    } else  { fRich700digipar = 0; fuseRICH700 = 0; }

    fTofGeometry = (HTofGeomPar *)       rtdb->getContainer("TofGeomPar");
    fMagnetPar   = (HMagnetPar *)        rtdb->getContainer("MagnetPar");
    fFieldPar    = (HMdcTrackGFieldPar*) rtdb->getContainer("MdcTrackGFieldPar");

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

    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
    if(fuseEMC && fEmcGeometry)
    {
	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++;
	    }
	}
    }

    //-------------------------------------------------------------
    // get iTof Geom helper
    if(fuseiTOF && fiTofGeometry)
    {
	fiTofGeometry->initGeomHelpers(pSpecGeomPar);
	for(Int_t s = 0; s < 6; s++)
	{
	    TVector3 norm,center;
	    fiTofGeometry->getCenterAndNormVecMod (s, center,norm);
	    centerVeciTof[s] = center;
	    normVeciTof  [s] = norm;

	    for(Int_t c = 0; c < 3; c++)
	    {
		cellSecTransiTof[s][c] = fiTofGeometry->getSecCellTrans(s, c);
                vector<HGeomVector> points;
		fiTofGeometry->getGeomVolumePointsCell( s,c, points);
		for(Int_t i = 0; i < 4; i++)
		{
                    points[i].setZ(0);
		    celliTOF[s][c][i] = points[i];
		}
	    }
	}
    }
    //-------------------------------------------------------------

    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 (fuseEMC&&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();

    //--------------------------------------------
    // RungeKutta  stuff

    HGeomTransform trans;
    Int_t done[6] = {0};
    HGeomVector v,vs,norm;
    for(Int_t s = 0; s < 6; s ++){
	for(Int_t m = 0; m < 4; m ++){
	    if(fSizesCells->modStatus(s, m)){
		mdcInstalled[s][m] = 1;
		HMdcSizesCellsMod& fSCMod = (*fSizesCells)[s][m];
		trans = *(fSCMod.getSecTrans());

		v.setXYZ(0.0, 0.0, 0.0);
		vs = trans.transFrom(v);

		HGeomVector r0_mod(0.0, 0.0, 0.0);
		HGeomVector rz_mod(0.0, 0.0, 1.0);
		norm = trans.transFrom(rz_mod) - trans.transFrom(r0_mod);

		centerVecMdc[s][m].SetXYZ(vs  .getX(),vs  .getY(),vs  .getZ());
		normVecMdc  [s][m].SetXYZ(norm.getX(),norm.getY(),norm.getZ());

		if(done[s] == 0 && m < 2){ // RK definition plane per sector (checking if mod 0 or 1 has to be used)
		    if(m == 0) v.setXYZ(0.0, 0.0, -500.0);
		    else       v.setXYZ(0.0, 0.0, -650.0);
		    vs = trans.transFrom(v);

		    centerVecRK[s].SetXYZ(vs  .getX(),vs  .getY(),vs  .getZ());
		    normVecRK  [s].SetXYZ(norm.getX(),norm.getY(),norm.getZ());
		    done[s] = 1;
		}
	    }
	}
    }

    Float_t fieldFactor = fMagnetPar->getScalingFactor();
    if(fMagnetPar->getPolarity() < 0) {
	fieldFactor = - fieldFactor;
    }

    rkpropagator.setFieldMap(fFieldPar->getPointer(), fieldFactor);

    for(Int_t s = 0; s < 6 ; s ++){
	planeMdc[s][0].setPlane(centerVecRK[s],normVecRK[s]);  // RK track is defined on this plane
	for(Int_t m = 0; m < 4; m++){
	    if(mdcInstalled[s][m]){
		planeMdc[s][m+1].setPlane(centerVecMdc[s][m],normVecMdc[s][m]);
	    }
	}
    }
    //--------------------------------------------
    fInitOK = kTRUE;
    return retval;
}

Int_t HParticleMetaMatcher::execute()
{
    if(frunWireManager) {
	wiremanager.execute();
	wiremanager.print();
    }
    
    if(frunFlagMETADouble){
	vector<HParticleMetaPredict> vpredict;
	loopFlagDoubleUsedMetaCells(vpredict);
    }

    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(!fInitOK) { Error("recalcRpc()","You have to init HParticleMetaMatcher via the tasklist!"); return 0;}
    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
    Double_t path = 0;
    traceToMeta(cand,metaHit,centerVecRpc[sec][ind],normVecRpc[sec][ind],path);
    HGeomVector vsec(metaHit.X(),metaHit.Y(),metaHit.Z());    // hit point mdc on META sec sys

    trackMetaSec = vsec;

    HGeomVector vlab  = labSecTrans[sec].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 = labSecTrans[sec].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(!fInitOK) { Error("recalcTof()","You have to init HParticleMetaMatcher via the tasklist!"); return 0;}
    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;
    Double_t path=0;
    traceToMeta(cand,metaHit,centerVecTof[sec][mod0],normVecTof[sec][mod0],path);
    HGeomVector vsec(metaHit.X(),metaHit.Y(),metaHit.Z());    // hit point mdc on META sec sys

    trackMetaSec = vsec;

    HGeomVector vlab  = labSecTrans[sec].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 = labSecTrans[sec].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(!fInitOK) { Error("recalcEmc()","You have to init HParticleMetaMatcher via the tasklist!"); return 0;}
    if(!fuseEMC || !fEmcGeometry) return NULL;

    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;

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

    trackMetaSec = vsec;

    HGeomVector vlab  = labSecTrans[sec].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 = labSecTrans[sec].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,
                                            Double_t& path1, Double_t& path2,
					    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(!fInitOK) { Error("predictRpcCell()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}

    path1 = path2 = 0;
    s = col0 = col1 = cell0 = cell1 = -1;
    hit0.setXYZ(0,0,0);
    hit1.setXYZ(0,0,0);

    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& modToLab = fRpcGeometry->getModule(sec,0)->getLabTransform();
    for(Int_t i = 0; i < 2; i++){
	if(i==0)traceToMeta(cand,metaHit,centerVecRpc[sec][i],normVecRpc[sec][i],path1);
	else    traceToMeta(cand,metaHit,centerVecRpc[sec][i],normVecRpc[sec][i],path2);

	HGeomVector vsecL(metaHit.X(),metaHit.Y(),metaHit.Z());    // hit point mdc on META sec sys
	vsec[i] = vsecL;
	HGeomVector vlabL  = labSecTrans[sec].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,
					    Double_t& path1, Double_t& path2,
					    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(!fInitOK) { Error("predictTofCell()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}
    path1 = path2 = 0;
    s = mod0 = mod1 = cell0 = cell1 = -1;
    hit0.setXYZ(0,0,0);
    hit1.setXYZ(0,0,0);

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

    Double_t pathLoc[8];
    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],pathLoc[i]);
	HGeomVector vsecL(metaHit.X(),metaHit.Y(),metaHit.Z());    // hit point mdc on META sec sys
        vsec[i] = vsecL;
	HGeomVector vlabL  = labSecTrans[sec].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;
                path1= pathLoc[m];
		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;
                path2= pathLoc[m];
		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! ct = "<<ct<<endl;
	    }
	    ct++;

	}
    }

    if(ct>2) cand->print();

    return ct > 0 ? kTRUE : kFALSE;

}

Bool_t HParticleMetaMatcher::predictEmcCell(HParticleCand* cand,HGeomVector& hit,Int_t& s,Int_t& pos,Int_t& cell,Double_t& path,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(!fInitOK) { Error("predictEmcCell()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}
    path = 0;
    s = pos = cell = -1;
    hit.setXYZ(0,0,0);

    if(!fuseEMC || !fEmcGeometry) return kFALSE;

    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& modToLab = fEmcGeometry->getModule(sec,0)->getLabTransform();
    traceToMeta(cand,metaHit,centerVecEmc[sec],normVecEmc[sec],path);
    HGeomVector vsec(metaHit.X(),metaHit.Y(),metaHit.Z());    // hit point mdc on META sec sys
    if(sys == 2) hit = vsec;

    HGeomVector vlab  = labSecTrans[s].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::predictiTofCell(HParticleCand* cand,HGeomVector& hit,Int_t& s, Int_t& cell, Double_t& path,
					     Int_t sys)
{
    // Reconstructs an intersection of HParticleCand inner
    // segment of Runge Kutta with the iTOF detector plane.
    // Possible hit iTOF cells are returned
    // by coordinates sector, cell. (-1 if no hit, hit 0,0,0).
    // The hit can be returned in different coordinate systems:
    // sys  1 = cell (default), 2 = sector, 3 = lab

    if(!fInitOK) { Error("predictiTofCell()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}
    path = 0;
    s = cell = -1;
    hit.setXYZ(0,0,0);
    path =-1;
    if(!cand) return 0;

    if(fDebug){ cout <<"predictiTofCell--------------------------------"<<endl;}
    s = cand->getSector();

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

    if(!fuseRK || cand->getMomentum() < 0){

	Int_t s    = cand->getSector();

       	TVector3 p1;
	TVector3 p2;
	TVector3 dir;

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

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

	Double_t theta = cand->getTheta();  // deg
	Double_t r     = cand->getR();      // mm
	Double_t z     = cand->getZ();      // mm

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

	dir = p2 - p1;

        TVector3 pointIntersect;
	HKalGeomTools::findIntersectionLinePlane(pointIntersect,p1,dir,centerVeciTof[s],normVeciTof[s]);

        HGeomVector vsec;
        HGeomVector vlab;
        HGeomVector vcell;

	vsec = pointIntersect;
        vlab = vsec;
	vlab = labSecTrans[s].transFrom(vlab);             // LAB

        Bool_t found = kFALSE;
	for(Int_t c = 0; c < ITOFMAXCELL; c ++){ // find which cell has been hit
	    vcell = vsec;
	    vcell = cellSecTransiTof[s][c].transTo(vcell);  // CELL

	    if(isIniTofCell(vcell,s,c)){
		cell = c;
                found = kTRUE;
		break;
	    }
	}

	if(found){
	    if(sys == 1) hit = vcell;
	    if(sys == 2) hit = vsec ;
	    if(sys == 3) hit = vlab ;
	}

	return kTRUE;

    } else { // use RK propagation

	TVectorD stateFrom(6);
	if(!getTrackStateRK(cand,stateFrom)) return kFALSE;
        s    = cand->getSector();
	Int_t chrg = cand->getCharge();

	TVectorD stateTo(6);  // transported state at the destination layer
        HParticleRKPlane planeiTof(centerVeciTof[s],normVeciTof[s]);

	rkpropagator.propagateToPlane(chrg,stateTo,stateFrom,planeiTof,path, kIterForward);

	HGeomVector vsec(stateTo(0),stateTo(1),stateTo(2));  // SEC
	HGeomVector vlab  = vsec;
	vlab = labSecTrans[s].transFrom(vlab);             // LAB
        HGeomVector vcell;

        Bool_t found = kFALSE;
	for(Int_t c = 0; c < ITOFMAXCELL; c ++){ // find which cell has been hit
	    vcell = vsec;
	    vcell = cellSecTransiTof[s][c].transTo(vcell);  // CELL

	    if(isIniTofCell(vcell,s,c)){
		cell = c;
                found = kTRUE;
		break;
	    }
	}

	if(found){
	    if(sys == 1) hit = vcell;
	    if(sys == 2) hit = vsec ;
	    if(sys == 3) hit = vlab ;
	}
        return kTRUE;
    }

    return  kFALSE;
}

Bool_t HParticleMetaMatcher::predictRich700Hit(HParticleCand* cand,Double_t zvertex,
					       Int_t& s, Int_t& col,Int_t& row, HGeomVector& hit, vector<richaddress>& vrichcaladd,
					       Int_t range
					      )
{
    // Reconstructs an intersection of HParticleCand inner
    // segment with RICH700 pmt plane (using reverse mapping). LAB SYS
    // Possible RICH hit PMT hits are returned
    // by coordinates (if no hit, hit 0,0,0).
    // s,col,row address of the hit.
    // vector vrichcal  struct richaddress {Int_t s,Int_t row,Int_t col}
    // contains all rich cal addresses in range +-range col and row arround
    // the hit point (region of interest)

    vrichcaladd.clear();
    s = col = row = -1;
    hit.setXYZ(0,0,0);

    if(cand == 0)            return kFALSE;
    if(fRich700digipar == 0) return kFALSE;

    Float_t theta = cand->getTheta();
    Float_t phi   = cand->getPhi();


    Int_t sC = cand->getSector();
    TVector3 p1v, p2v,poutv,dir;
    HGeomVector p1,p2,pout;


    //----------------------------------------------------------------------
    // calculate segment points for hit point on rich mirror
    // this hit point will be used for estimate of theta, phi
    // to suppress a part of the multiple scattering on the
    // rich shell+mirror
    HParticleTool::calcSegPoints(cand,p1v,p2v,kFALSE,kTRUE); //inner seg SEC
    dir = p2v - p1v;

    HParticleRKPlane pl;
    pl.setPlane(centerVecMdc[sC][0],normVecMdc[sC][0]);
    pl.findIntersection(poutv,p1v,dir); // SEC sys

    p1v = poutv;
    p2v = p1v + dir;

    pointLabToSecSys(sC,p1v,poutv,kFALSE); // LAB
    p1 = poutv;
    pointLabToSecSys(sC,p2v,poutv,kFALSE); // LAB
    p2 = poutv;

    // calculate the hit of segment on rich mirror
    HParticleTool::calcRichMirrorHit(p1,p2,pout,-1000.); // LAB , take the predefined zpos for RICH

    //----------------------------------------------------------------------
    // take into account emmission from vertex (z only)
    poutv.SetXYZ(pout.X(), pout.Y(), pout.Z() - zvertex);
    //----------------------------------------------------------------------

    //----------------------------------------------------------------------
    // theta, phi in candidate LAB
    // allow multiple scattering corr max +-5 deg
    theta = poutv.Theta() * TMath::RadToDeg();
    phi   = poutv.Phi()   * TMath::RadToDeg();
    if(phi < 0 ) phi += 360;
    if(fabs(theta-cand->getTheta()) > 5. || fabs(phi-cand->getPhi()) > 5.) {
            theta = cand->getTheta();
            phi   = cand->getPhi();
    }
    //----------------------------------------------------------------------
    //----------------------------------------------------------------------


    //----------------------------------------------------------------------
    // This is inverse transformation Theta, Phi -> XY Pmt
    // later calulate sec,col,row rich cal adressesn in window
    // +- range arround the prdeicted ring center
    Float_t x,y;
    Float_t thetaCor,phiCor;

    thetaCor = theta;
    phiCor   = phi;

    // correct theta, phi for emission from vertex
    fRich700digipar->getAlignedThetaPhiInv(theta,phi, thetaCor, phiCor);
    fRich700digipar->getRingCenterXY(thetaCor,phiCor, zvertex, x, y);  // x,y mm on plane

    Int_t pmtid = fRich700digipar->getPMTId(x,y);
    if(pmtid < 0) {
        Warning ("predictRich700Hit()","No pmtid for x = %8.3f, y = %8.3f, corr theta = %8.3f, phi = %8.3f, cand theta = %8.3f, phi = %8.3f !",x,y,thetaCor,phiCor,cand->getTheta(),cand->getPhi());
	return kFALSE;
    }

    HRich700PmtData* rdata = fRich700digipar->getPMTData(pmtid);

    if(!rdata) {
        Warning ("predictRich700Hit()","No pmtid data for id  %d !",pmtid);
	return kFALSE;
    }
    hit.setXYZ(x,y,rdata->fZ);  // LAB

    //----------------------------------------------------------------------
    // calculate sec, row, col of ring center
    Int_t fNofPixelsInRow = fRich700digipar->getNPixelInRow();
    Double_t pmtsize      = fRich700digipar->getPmtSize();

    Double_t dx   = x - (rdata->fX - pmtsize/2.) ;
    Double_t dy   = y - (rdata->fY - pmtsize/2.) ;
    Double_t step = (pmtsize)/fNofPixelsInRow;
    Double_t halfPmtSensSize = pmtsize / 2.;

    row = rdata->fIndY * fNofPixelsInRow + Int_t(dy/step);
    col = rdata->fIndX * fNofPixelsInRow + Int_t(dx/step);
    s   = fRich700digipar->getSector(x,y);
    //----------------------------------------------------------------------

    //----------------------------------------------------------------------
    // scan the region of interest +- range arround
    // hit center and store the adresses in a vector
    Int_t colL = col - range;
    if(colL < 0) colL = 0;
    Int_t colU = col + range;
    if(colU > RICH700_MAX_COLS) colU = RICH700_MAX_COLS;

    Int_t rowL = row - range;
    if(rowL < 0) rowL = 0;
    Int_t rowU = row + range;
    if(rowU > RICH700_MAX_ROWS) rowU = RICH700_MAX_ROWS;

    for(Int_t c = colL; c <= colU; c++){
	for(Int_t r = rowL; r <= rowU; r++){
	    Int_t id = fRich700digipar->getPMTId(c, r, kTRUE);
	    if(id < 0) continue;

	    rdata = fRich700digipar->getPMTData(id);

            UInt_t pix_x = c % fNofPixelsInRow;
            UInt_t pix_y = r % fNofPixelsInRow;

	    Double_t xLoc = (pix_x + 0.5) * step - halfPmtSensSize;
	    Double_t yLoc = (pix_y + 0.5) * step - halfPmtSensSize;

	    x = xLoc + rdata->fX;
	    y = yLoc + rdata->fY;

	    Int_t sec = fRich700digipar->getSector(x,y);

	    richaddress add;
	    add.set(sec,r,c);
            vrichcaladd.push_back(add);

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

    return kTRUE;

}


Bool_t HParticleMetaMatcher::predictRich700Hit(HParticleCand* cand,Double_t zvertex,
					       Int_t& s, Int_t& col,Int_t& row, HGeomVector& hit, vector<HRichCal*>& vrichcal,
					       Int_t range
					      )
{
    // Reconstructs an intersection of HParticleCand inner
    // segment with RICH700 pmt plane (using reverse mapping). LAB SYS
    // Possible RICH hit PMT hits are returned
    // by coordinates (if no hit, hit 0,0,0).
    // s,col,row address of the hit.
    // vector vrichcal contains all rich cal found in range +-range col
    // and row arround the hit point (region of interest)

    vrichcal.clear();
    vector<richaddress> vrichcaladd;
    if(! predictRich700Hit(cand,zvertex, s,col,row,hit,vrichcaladd, range)) return kFALSE;

    HCategory* richCalCat = (HCategory*)HCategoryManager::getCategory(catRichCal);
    if(!richCalCat) return kFALSE;
    HLocation loc;
    loc.set(3,0,0,0);

    HRichCal* cal;
    for(UInt_t k = 0; k < vrichcaladd.size(); k++ ){
	richaddress& ad = vrichcaladd[k];
        loc[0] = ad.sec;
        loc[1] = ad.row;
        loc[2] = ad.col;

	cal = (HRichCal*) richCalCat->getObject(loc);
	if(!cal) continue;
        vrichcal.push_back(cal);
    }

    return kTRUE;

}

Bool_t HParticleMetaMatcher::printRichCals(vector<HRichCal*>& vrichcal,Int_t col,Int_t row,Int_t range)
{
    // print richcals region of interest ring center (col,row) +- range
    // works in conunction with
    // Bool_t HParticleMetaMatcher::predictRich700Hit(HParticleCand* cand,Double_t zvertex,
    //                                               Int_t& s, Int_t& col,Int_t& row, HGeomVector& hit, vector<HRichCal*>& vrichcal,
    //                   			     Int_t range)

    if(vrichcal.size() == 0) return kTRUE;
    if(col < 0 || row < 0)   return kFALSE;
    cout<<"#########################################"<<endl;
    cout<<"center col "<<setw(3)<<col<<", row "<<setw(3)<<row<<", ncals "<< vrichcal.size()
	<<" range "<<setw(3)<<range
	<<" col  ["<<setw(3)<<col-range<<","<<setw(3)<<col+range<<"], row ["
	<<setw(3)<<row-range<<","<<setw(3)<<row+range<<"]"<<endl;

    TString line = std::string(2*range+1, '.');
    vector<TString > box;
    for(Int_t i = 0 ; i < 2*range+1; i++) box.push_back(line);
    HRichCal* richcal;
    Bool_t ok = kTRUE;
    for(UInt_t i = 0; i < vrichcal.size(); i++){
	richcal = vrichcal[i];
	Int_t r = richcal->getRow() - row + range;
	Int_t c = richcal->getCol() - col + range;
	if(r < 0 || r > 2*range || c< 0 || c > 2*range  ) {
	    cout<< "out of box -> r "<<setw(3)<<r<<" , c "<<setw(3)<<endl;
	    ok = kFALSE;
	    continue;
	}
	box[r][c] = 'x';

    }

    for(Int_t i = 0 ; i < 2*range+1; i++) cout<<box[i]<<endl;
    cout<<"#########################################"<<endl;

    return ok;
}

Bool_t HParticleMetaMatcher::predictMdcHit(HParticleCand* cand,HGeomVector& hit,Int_t m,Double_t& path)
{
    // returns  MDC hit [x,y] in  MDC system (z=0) for MDC module m [0-3]
    // recalculated from r,z,phi,theta  (RK or straight Line depending on fuseRK)
    // Different use Cases :
    // return fuseRK : kFALSE if no outer MDC was reconstructed and m>=2 or MDC Modules is
    //                        not in setup, no straight line reconstruction possible
    //        fuseRK : kTRUE  works if outer SegInd =-1 and MDC Module is in setup.
    //                        Hits of Missing modules will be  0,0,0

    if(!fInitOK) { Error("predictMdcHit()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}
    hit.setXYZ(0.,0.,0.);
    path = 0;
    if(!cand)                               return kFALSE;
    if(!mdcInstalled[cand->getSector()][m]) return kFALSE;

    HMdcSizesCellsMod& sizesMod= (*fSizesCells)[cand->getSector()][m];

    if(fuseRK){

	TVector3 hit0,hit1,hit2,hit3;
        TVectorD pathlength;
	if(!getMDCPosRK(cand,hit0,hit1,hit2,hit3,pathlength)) return kFALSE; // sec coord. sys
	path = pathlength(m);

        const HGeomTransform* trans = sizesMod.getSecTrans(); // sec <-> mod sys

	if(m == 0)      { hit = hit0; trans->transTo(hit);}  // sec -> mod
	else if(m == 1) { hit = hit1; trans->transTo(hit);}
	else if(m == 2) { hit = hit2; trans->transTo(hit);}
	else if(m == 3) { hit = hit3; trans->transTo(hit);}

    } else {
        if((cand->getOuterSegInd() < 0 && m >= 2)) return kFALSE;
        if(cand->getPhi2() == -1       && m >= 2)  return kFALSE; // in old dsts the our seg info is not filled!
	if(mdcInstalled[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);

	} else {
	    return kFALSE;
	}
    }
    return kTRUE;
}

Bool_t HParticleMetaMatcher::predictMdcHitSec(HParticleCand* cand,HGeomVector& hit,Int_t m,Double_t& path)
{
    // returns  MDC hit [x,y,z] in sec coord. system for MDC module m [0-3]
    // recalculeated from r,z,phi,theta
    // return fuseRK : kFALSE if no outer MDC was reconstructed and m>=2 or MDC Modules is
    //                        not in setup, no straight line reconstruction possible
    //        fuseRK : kTRUE  works if outer SegInd =-1 and MDC Module is in setup.
    //                        Hits of Missing modules will be  0,0,0
    if(!fInitOK) { Error("predictMdcHitSec()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}

    if(!predictMdcHit(cand,hit,m,path)) 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,Double_t& path)
{
    // returns  MDC hit [x,y,z] in Lab coord. system for MDC module m [0-3]
    // recalculeated from r,z,phi,theta
    // return fuseRK : kFALSE if no outer MDC was reconstructed and m>=2 or MDC Modules is
    //                        not in setup, no straight line reconstruction possible
    //        fuseRK : kTRUE  works if outer SegInd =-1 and MDC Module is in setup.
    //                        Hits of Missing modules will be  0,0,0
    if(!fInitOK) { Error("predictMdcHitLab()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}

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

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

}

Bool_t HParticleMetaMatcher::calcLineRpcCell(const Int_t s, const TVectorD& state,HGeomVector& hit0,HGeomVector& hit1,
					     Int_t& col0,Int_t& cell0,Int_t& col1,Int_t& cell1,
					     Double_t& path1, Double_t& path2,
					     Int_t sys)
{
    // Reconstructs an intersection of state (x,y,z,px,py,pz) [mm,MeV/c] SEC sys
    // in sector  s straight line 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.
    // The hit can be returned in different coordinate systems:
    // sys  1 = mod (default), 2 = sector, 3 = lab


    if(!fInitOK) { Error("calcLineRpcCell()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}

    path1 = path2 = 0;
    col0 = col1 = cell0 = cell1 = -1;
    hit0.setXYZ(0,0,0);
    hit1.setXYZ(0,0,0);


    if(fDebug){ cout <<"calcLineRpcCell--------------------------------"<<endl;}

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

    TVector3 pos(state(0),state(1),state(2));  // SEC sys
    TVector3 dir(state(3),state(4),state(5));

    HGeomVector vmod[2];
    HGeomVector vsec[2];
    HGeomVector vlab[2];

    HParticleRKPlane pl;
    TVector3 pointIntersect;

    HGeomTransform& modToLab = fRpcGeometry->getModule(s,0)->getLabTransform();
    HGeomVector diff;

    for(Int_t i = 0; i < 2; i++){
	pl.setPlane(centerVecRpc[s][i],normVecRpc[s][i]);
	pl.findIntersection(pointIntersect,pos,dir);  // SEC sys
	vsec[i] = pointIntersect;
        vlab[i] = labSecTrans[s].transFrom(vsec[i]); // LAB sys
        vmod[i] = modToLab.transTo(vlab[i]);   // MOD sys

	diff = vsec[i] - pos;

        if(i == 0) path1 = diff.length();
        else       path2 = diff.length();

    }

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

    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,s,col,c)) continue;

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

	}
    }

    return ct > 0 ? kTRUE : kFALSE;

}

Bool_t HParticleMetaMatcher::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)
{
    // Reconstructs an intersection of state (x,y,z,px,py,pz) [mm,Mev/c] SEC sys
    // straigth line 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.
    // The hit can be returned in different coordinate systems:
    // sys  1 = mod (default), 2 = sector, 3 = lab

    if(!fInitOK) { Error("calcLineTofCell()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}
    path1 = path2 = 0;
    mod0 = mod1 = cell0 = cell1 = -1;
    hit0.setXYZ(0,0,0);
    hit1.setXYZ(0,0,0);


    if(fDebug){ cout <<"calcLineTofCell--------------------------------"<<endl;}

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

    TVector3 pos(state(0),state(1),state(2)); // SEC sys
    TVector3 dir(state(3),state(4),state(5));

    HGeomVector vmod[8];
    HGeomVector vsec[8];
    HGeomVector vlab[8];
    Double_t pathLoc[8];

    HParticleRKPlane pl;
    TVector3 pointIntersect;

    HGeomVector diff;


    for(Int_t i = 0; i < 8; i++){ // modules
	HGeomTransform& modToLab = fTofGeometry->getModule(s,i)->getLabTransform();

	pl.setPlane(centerVecTof[s][i],normVecTof[s][i]);
	pl.findIntersection(pointIntersect,pos,dir); // SEC sys
	vsec[i] = pointIntersect;
	vlab[i] = labSecTrans[s].transFrom(vsec[i]); // LAB sys
	vmod[i] = modToLab.transTo(vlab[i]);   // MOD sys

	diff = vsec[i] - pos;
	pathLoc[i] = diff.length();
    }

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

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

	    if     (ct == 0) {
		mod0 = m; cell0 = c;
		path1= pathLoc[m];
		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;
		path2= pathLoc[m];
		if(sys == 1)hit1 = vmod[m];
		if(sys == 2)hit1 = vsec[m];
		if(sys == 3)hit1 = vlab[m];
	    } else  {
		cout <<"calcLineTofCell() : more than 2 Tof rods found!"<<endl;
	    }
	    ct++;

	}
    }

    return ct > 0 ? kTRUE : kFALSE;

}

Bool_t HParticleMetaMatcher::calcLineEmcCell(const Int_t s, const TVectorD& state,HGeomVector& hit,Int_t& pos,Int_t& cell,Double_t& path,Int_t sys)
{
    // Reconstructs an intersection of state (x,y,z,px,py,pz) [mm,Mev/c] SEC sys
    // straight line 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.
    // The hit can be returned in different coordinate systems:
    // sys  1 = mod (default), 2 = sector, 3 = lab

    if(!fInitOK) { Error("calLineEmcCell()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}
    path = 0;
    pos = cell = -1;
    hit.setXYZ(0,0,0);

    if(!fuseEMC || !fEmcGeometry) return kFALSE;

    if(fDebug){ cout <<"calcLineEmcCell--------------------------------"<<endl;}

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

    TVector3 vpos(state(0),state(1),state(2)); // SEC sys
    TVector3 vdir(state(3),state(4),state(5));

    HGeomVector vmod;
    HGeomVector vsec;
    HGeomVector vlab;

    HGeomTransform& modToLab = fEmcGeometry->getModule(s,0)->getLabTransform();

    HParticleRKPlane pl(centerVecEmc[s],normVecEmc[s]);
    TVector3 pointIntersect;
    HGeomVector diff;

    pl.findIntersection(pointIntersect,vpos,vdir);  // SEC sys

    vsec = pointIntersect;
    vlab = labSecTrans[s].transFrom(vsec); // LAB sys
    vmod = modToLab.transTo(vlab);   // MOD sys
    diff = vsec - pos;
    path = diff.length();

    if     (sys == 1) hit = vmod;
    else if(sys == 2) hit = vsec;
    else if(sys == 3) hit = vlab;

    pos  = -1;
    cell = -1;

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

}

Bool_t HParticleMetaMatcher::calcLineiTofCell(const Int_t s,const TVectorD& state,HGeomVector& hit,
					     Int_t& cell,Double_t& path, Int_t sys)
{
    // Reconstructs an intersection of state (x,y,z,px,py,pz) [mm,Mev/c] SEC sys
    // straigth line with the iTOF detector plane.
    // Possible hit iTOF cells are returned
    // by coordinates sector, cell.
    // cell is -1 if no cell has been hit on this plane.
    // The hit points will be 0,0,0 in this case.
    // The hit can be returned in different coordinate systems:
    // sys  1 = cell (default), 2 = sector, 3 = lab

    if(!fInitOK) { Error("calcLineiTofCell()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}
    path = 0;
    cell = -1;
    hit.setXYZ(0,0,0);


    if(fDebug){ cout <<"calcLineiTofCell--------------------------------"<<endl;}

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

    TVector3 pos(state(0),state(1),state(2)); // SEC sys
    TVector3 dir(state(3),state(4),state(5));

    HGeomVector vcell[iTOFMAXCELL];
    HGeomVector vsec[iTOFMAXCELL];
    HGeomVector vlab[iTOFMAXCELL];
    Double_t pathLoc[iTOFMAXCELL];

    TVector3 pointIntersect;

    HParticleRKPlane pl;
    pl.setPlane(centerVeciTof[s],normVeciTof[s]);
    pl.findIntersection(pointIntersect,pos,dir); // SEC sys


    HGeomVector diff;

    for(Int_t i = 0; i < iTOFMAXCELL; i++){ // modules

	vlab[i]  = pointIntersect;
	vlab[i]  = labSecTrans[s].transFrom(vlab[i]); // LAB sys
        vsec[i]  = pointIntersect;
	vcell[i] = pointIntersect;
	diff = vsec[i] - pos;     // SEC
	pathLoc[i] = diff.length();
	vcell[i] = cellSecTransiTof[s][i].transTo(vsec[i]);   // CELL sys
    }

    cell = -1;

    for(Int_t c = 0; c < iTOFMAXCELL; c++){
	if(isIniTofCell(vcell[c],s,c)) {

	    cell = c;
	    path = pathLoc[c];
	    if(sys == 1) hit = vcell[c];
	    if(sys == 2) hit = vsec[c];
	    if(sys == 3) hit = vlab[c];
            return kTRUE;
	}
    }

    return  kFALSE;

}

Bool_t HParticleMetaMatcher::calcLineMdcHit(const Int_t s,const TVectorD& state,HGeomVector& hit,const Int_t m,Double_t& path,Int_t sys)
{
    // returns  MDC hit [x,y,z] for MDC module m [0-3]
    // recalculated from state (x,y,z,px,py) [mm,Mev/c] SEC sys using straight Line
    // The hit can be returned in different coordinate systems:
    // sys  1 = mod (default), 2 = sector, 3 = lab

    if(!fInitOK) { Error("calcLineMdcHit()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}
    hit.setXYZ(0.,0.,0.);
    path = 0;
    if(m < 0 || m > 3)  return kFALSE;

    if(!mdcInstalled[s][m]) return kFALSE;

    TVector3 pos(state(0),state(1),state(2));
    TVector3 dir(state(3),state(4),state(5));
    TVector3 p2 = pos + dir;

    HGeomVector vmod;
    HGeomVector vsec;
    HGeomVector vlab;

    Double_t x,y,z = 0;
    HMdcSizesCellsMod& sizesMod = (*fSizesCells)[s][m];
    sizesMod.calcInterTrMod(pos.X(),pos.Y(),pos.Z(),p2.X(),p2.Y(),p2.Z(), x,y); // x,y on plane z=0
    vmod.setXYZ(x,y,z); // MOD sys

    if(sys == 1) {
	hit = vmod;
    }

    (*fSizesCells)[s][m].transFromZ0(x,y,z);
    vsec.setXYZ(x,y,z); // SEC sys

    HGeomVector diff = vsec - pos;
    path = diff.length();

    if(sys == 1 || sys == 2) {
        if(sys == 2)  hit = vsec;
	return kTRUE;
    }

    const HGeomTransform* transLab = (*fSizesCells)[s].getLabTrans();
    vlab = transLab->transFrom(vsec); // LAB sys
    hit = vlab;

    return kTRUE;
}

void  HParticleMetaMatcher::getWireInfoDirect(HParticleCand* cand,HParticleWireInfo& w,Bool_t fillStat)
{
    if(!fInitOK) { Error("getWireInfoDirect()","You have to init HParticleMetaMatcher via the tasklist!"); return ;}
    wiremanager.getWireInfoDirect (cand,w,fillStat);
}

Bool_t HParticleMetaMatcher::calcState(Int_t& s, TVectorD& state,const TVector3& pos, const Double_t theta, const Double_t phi, const Double_t mom)
{
    // calculate track state (x,y,z,px,py,pz) [mm,Mev/c] SEC sys in sector s
    // from theta,phi, pos(x,y,z) LAB sys [DEG [ 0 - 360 ]] and mom [Mev/c]

    if(state.GetNrows() != 6) state.ResizeTo(6);
    if(phi < 0. || phi > 360.) { ::Error("calcState()","Phi out of range [0-360][DEG] phi = %lf",phi); return kFALSE;}
    if(mom < 0.)               { ::Error("calcState()","momentum negative  mom = %lf"           ,mom); return kFALSE;}
    if(phi < 60.) s = 5;
    else          s = ((Int_t) (phi/60.)) - 1;

    TVectorD stateLab(6);
    // First set state in LAB sys
    state(0) = pos.X();
    state(1) = pos.Y();
    state(2) = pos.Z();

    TVector3 dir;
    dir.SetMagThetaPhi(mom,theta * TMath::DegToRad(),phi * TMath::DegToRad());

    state(3) = dir.X();
    state(4) = dir.Y();
    state(5) = dir.Z();

    return stateLabToSecSys(s,stateLab,state,kTRUE); // LAB -> SEC

}

Bool_t HParticleMetaMatcher::kineToParticleCand(HGeantKine* kine,HParticleCandSim& cand, Bool_t realmom)
{
    // fill candidate  sector, phi,theta,p , lorentz vector, pid from HGeantKine , track num etc.
    // sector is calculated from the direction of the particle.
    // if realmom (default kFALSE) is kTRUE, particles with charge > 1 are treated
    // correctly. In default mode kFALSE the momenta/charge are treated to be
    // compatible with RK reconstruction (charge +- 1) and momentum is divided
    // by charge. In this way the candidates can be used with the propagateXXX() functions.

    if(!kine) return kFALSE;
    Int_t id = kine->getID();

    if(HPhysicsConstants::mass(id < 0)) return kFALSE;

    HGeomVector mom,pos;
    kine->getMomentum(mom);
    kine->getVertex  (pos);
    HGeomVector p2 = pos + mom; // LAB

    Double_t zm, r0,theta,phi;

    Double_t phiLab   = kine->getPhiDeg  (kTRUE); // LAB sys 0-360 DEG
    Double_t thetaLab = kine->getThetaDeg();    // DEG
    Int_t s           = kine->getSector();      // sector from phiLab

    HParticleTool::calcMdcSeg(pos,p2,zm,r0,theta,phi);

    Int_t    chrg = HPhysicsConstants::charge(id);
    Double_t pout = mom.length();
    Double_t beta = HParticleTool::beta(id,pout);
    if(realmom) pout /= abs(chrg);

    if(!realmom) {
	if       (chrg < 0) {
	    chrg = -1;
	} else if(chrg > 0) {
	    chrg =  1;
	}
    }

    cand.setSector     (s);
    cand.setPhi        (phiLab);
    cand.setTheta      (thetaLab);
    cand.setMomentum   (pout);
    cand.setMomentumOrg(pout);
    cand.setBeta       (beta);
    cand.setBetaOrg    (beta);
    cand.setCharge     (chrg);
    cand.setPID        (kine->getID());
    cand.calc4vectorProperties(HPhysicsConstants::mass(id));  // sync TLorentzVector

    cand.setGeantPID           (id);
    cand.setGeantTrack         (kine->getTrack());
    cand.setGeantxMom          (mom.getX());
    cand.setGeantyMom          (mom.getY());
    cand.setGeantzMom          (mom.getZ());
    cand.setGeantxVertex       (pos.getX());
    cand.setGeantyVertex       (pos.getY());
    cand.setGeantzVertex       (pos.getZ());

    HLinearCategory* kineCat =(HLinearCategory*) HCategoryManager::getCategory(catGeantKine,0,"HGeantKine");
    if(kineCat){
        Int_t parentTrack = kine->getParentTrack();
	if(parentTrack >= 0){
	    HGeantKine* parent = HGeantKine::getParent(parentTrack,kineCat); // -1 inside
	    if(parent){
                cand.setGeantParentPID     (parent->getID());
                cand.setGeantParentTrackNum(parentTrack);
                Int_t grandParentTrack = parent->getParentTrack();
		if(grandParentTrack >= 0){
		    HGeantKine* grandparent = HGeantKine::getParent(grandParentTrack,kineCat);
                    if(grandparent){
			cand.setGeantGrandParentPID     (grandparent->getID());
			cand.setGeantGrandParentTrackNum(grandParentTrack);
		    }
		}
	    }
	}
    }

    cand.setGeantCreationMechanism(kine->getMechanism());
    cand.setGeantMediumNumber     (kine->getMedium());
    cand.setGeantGeninfo          (kine->getGeneratorInfo());
    cand.setGeantGeninfo1         (kine->getGeneratorInfo1());
    cand.setGeantGeninfo2         (kine->getGeneratorInfo2());
    cand.setGeantGenweight        (kine->getGeneratorWeight());
    return kTRUE;
}

Bool_t HParticleMetaMatcher::kineToToState(HGeantKine* kine,TVectorD& state,Int_t& charge,Bool_t realmom,const Int_t forceSector)
{
    // fill track state (x,y,z,px,py,pz)[mm,Mev/c, SEC sys]
    // sector from LAB ->SEC tranformation is calculated from the direction of the particle.
    // if realmom (default kFALSE) is kTRUE, particles with charge > 1 are treated
    // correctly. In default mode kFALSE the momenta/charge are treated to be
    // compatible with RK reconstruction (charge +- 1) and momentum is divided
    // by charge. In this way the candidates can be used with the propagateXXX() functions.
    // If forceSector !=1 (default -1) the state will transformed to the forced sector instead
    // of retrieving the sector from the position. This option is usefull for
    // the calcLineEmcCell() etc. where a match in a given sector should be checked.

    charge = 0 ;

    if(!kine) return kFALSE;
    Int_t id = kine->getID();

    if(HPhysicsConstants::mass(id < 0)) return kFALSE;

    if(forceSector !=-1 && (forceSector < 0 || forceSector > 5)){
	::Error("kineToState()","forceSector has illegal value %i",forceSector);
	return kFALSE;
    }


    HGeomVector mom,pos;
    kine->getMomentum(mom);
    kine->getVertex  (pos);
    HGeomVector p2Lab = pos + mom; // LAB
    Int_t s =-1;

    if(forceSector == -1) {
        s = kine->getSector();
    } else {
        s = forceSector;
    }

    TVector3 posSec,momSec;
    pointLabToSecSys(s,pos,posSec,kTRUE);  // LAB ->SEC
    pointLabToSecSys(s,mom,momSec,kTRUE);  // LAB ->SEC

    Int_t chrg = HPhysicsConstants::charge(id);

    Double_t scale = 1./abs(chrg);
    if(!realmom)  momSec *= scale;

    if(!realmom) {
	if       (chrg < 0) {
	    chrg = -1;
	} else if(chrg > 0) {
	    chrg =  1;
	}
    }

    charge = chrg;
    state(0) = posSec.X();
    state(1) = posSec.Y();
    state(2) = posSec.Z();
    state(3) = momSec.X();
    state(4) = momSec.Y();
    state(5) = momSec.Z();

    return kTRUE;
}

Bool_t  HParticleMetaMatcher::stateLabToSecSys(const Int_t s, const TVectorD& stateLab,TVectorD& stateSec,Bool_t to)
{
    // returns state (x,y,z,px,py,px) in sector coord. sys
    // to = kTRUE (default) Lab Sys-> Sec Sys
    // to = kFALSE          Sec Sys-> Lab Sys
    stateSec.ResizeTo(6);
    if(s < 0 || s > 5) return kFALSE;

    TVector3 posLab(stateLab(0),stateLab(1),stateLab(2));
    TVector3 momLab(stateLab(3),stateLab(4),stateLab(5));

    TVector3 posSec,momSec;
    pointLabToSecSys(s,posLab,posSec,to);
    pointLabToSecSys(s,momLab,momSec,to);

    stateSec(0) = posSec.X();
    stateSec(1) = posSec.Y();
    stateSec(2) = posSec.Z();
    stateSec(3) = momSec.X();
    stateSec(4) = momSec.Y();
    stateSec(5) = momSec.Z();

    return kTRUE;
}

Bool_t  HParticleMetaMatcher::pointLabToSecSys(const Int_t s, const TVector3& pointLab,TVector3& pointSec,Bool_t to)
{
    // returns point in sector coord. sys
    // to = kTRUE (default) Lab Sys-> Sec Sys
    // to = kFALSE          Sec Sys-> Lab Sys
    pointSec.SetXYZ(0,0,0);
    if(s < 0 || s > 5) return kFALSE;

    // get the rotation angle to Lab from sector
    Double_t  ph = 0;
    if(s < 5)  ph =  60.0f * s;
    else       ph = -60.0f;
    ph*=TMath::DegToRad();

    pointSec = pointLab;

    if(to) pointSec.RotateZ(-ph);
    else   pointSec.RotateZ( ph);

    return kTRUE;
}


Bool_t  HParticleMetaMatcher::pointLabToSecSys(const Int_t s, const HGeomVector& pointLab,TVector3& pointSec,Bool_t to)
{
    // returns point in sector coord. sys
    // to = kTRUE (default) Lab Sys-> Sec Sys
    // to = kFALSE          Sec Sys-> Lab Sys
    pointSec.SetXYZ(0,0,0);
    if(s < 0 || s > 5) return kFALSE;

    pointSec.SetXYZ(pointLab.getX(),pointLab.getY(),pointLab.getZ());

    // get the rotation angle to Lab from sector
    Double_t  ph = 0;
    if(s < 5)  ph =  60.0f * s;
    else       ph = -60.0f;
    ph*=TMath::DegToRad();

    if(to) pointSec.RotateZ(-ph);
    else   pointSec.RotateZ( ph);

    return kTRUE;
}

Bool_t   HParticleMetaMatcher::getTargetMiddlePoint(const Int_t s,HGeomVector& midpoint)
{
    // get target middle point (SEC Sys) used for Track reconstruction
    // to calculate the final path length
    if(!fInitOK) { Error("getTargetMiddlePoint()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}
    if(s<0||s>5) return kFALSE;
    if(!fSizesCells) return kFALSE;
    HMdcSizesCellsSec& fSCSec=(*fSizesCells)[s];
    midpoint = fSCSec.getTargetMiddlePoint();

    return kTRUE;
}

Bool_t   HParticleMetaMatcher::getDistanceToTargetMiddlePoint(const Int_t s,const TVectorD& stateRK, Double_t& dist)
{
    // get distance of RK state (SEC sys) at definition plane to target middle point
    // used for Track reconstruction to calculate the final path length
    dist = 0;
    if(!fInitOK) { Error("getDistanceToTargetMiddlePoint()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}
    if(s<0||s>5) return kFALSE;

    HMdcSizesCellsSec& fSCSec=(*fSizesCells)[s];
    const HGeomVector& midpoint = fSCSec.getTargetMiddlePoint();
    HGeomVector pos(stateRK(0),stateRK(1),stateRK(2));
    HGeomVector diff = pos - midpoint;
    dist = diff.length();
    return kTRUE;
}


Bool_t   HParticleMetaMatcher::getTrackStateRK(HParticleCand* cand,TVector3& pos,TVector3& momv,Double_t momentum,Int_t charge)
{
    pos .SetXYZ(0,0,0);
    momv.SetXYZ(0,0,0);
    if(!fInitOK) { Error("getTrackStateRK()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}
    TVectorD sv(6);
    Bool_t ret = getTrackStateRK(cand,sv,momentum,charge);
    if(!ret) return kFALSE;
    pos .SetXYZ(sv(0),sv(1),sv(2));
    momv.SetXYZ(sv(3),sv(4),sv(5));
    return kTRUE;
}


Bool_t  HParticleMetaMatcher::getTrackStateRK(HParticleCand* cand, TVectorD& sv,Double_t mom,Int_t charge)
{
    // calulates the trackstate from HParticleCand at the definition plane
    // in front of the inner MDCs (SEC SYS).
    // if the optional arguments charge and mom
    // are not 0 they will replace the   charge and momentum from candidate.
    // This option can be used to estimate outer MDC hits for candidates without
    // reconstructed hits

    if(!fInitOK) { Error("getTrackStateRK()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}

    if(!cand) return kFALSE;
    if(cand->getCharge() == 0 && mom == 0 && charge == 0) return kFALSE;  // needs momentum and charge to work
    Int_t s = cand->getSector();

    Int_t    chrg     = charge;
    Double_t momentum = fabs(mom);

    if(chrg     == 0) chrg     = cand->getCharge();
    if(momentum == 0) momentum = cand->getMomentum();

    // calc pos,dir from straight line
    TVector3 pos,dir,posAt;
    TVector3 p1,p2;
    HParticleTool::calcSegPoints(cand,p1,p2,kFALSE,kTRUE); // inner seg, sec coord. sys   2 points

    pos = p1;           //sec coord. sys
    dir = p2 - p1;      //sec coord. sys

    planeMdc[s][0].findIntersection(posAt,pos,dir);

    dir  = dir.Unit();
    dir *= momentum;

    sv(0) = posAt.X();
    sv(1) = posAt.Y();
    sv(2) = posAt.Z();
    sv(3) = dir.X();
    sv(4) = dir.Y();
    sv(5) = dir.Z();

    return kTRUE;
}

Bool_t  HParticleMetaMatcher::propagateToPlane(Int_t chrg,TVectorD& stateTo, const TVectorD& stateFrom,const HParticleRKPlane& planeTo,Double_t& pathL,Bool_t direction)
{
    // propagates track state stateFrom (x,y,px,py,pz) [mm,Mev/c, SEC sys] to
    // plane planeTo. The track state at the destination layer as well as the
    // path length is returned. The propagation is set kIterForward ( down beam)
    // or kIterBackward (beam up)
    if(!fInitOK) { Error("propagateToPlane()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}
    if(stateFrom.GetNrows() != 6) return kFALSE;
    if(stateTo  .GetNrows() != 6) stateTo.ResizeTo(6);

    return rkpropagator.propagateToPlane(chrg,stateTo,stateFrom,planeTo,pathL,direction);

}

Bool_t  HParticleMetaMatcher::propagateToPCA(Int_t chrg,TVectorD& stateTo,const TVectorD& stateFrom, const TVector3& point,Double_t& pathlength,Bool_t direction)
{
    // propagates track state stateFrom (x,y,px,py,pz) [mm,Mev/c, SEC sys] to
    // to point of clostest approach to point. The track state at the destination pca as well as the
    // path length is returned. The propagation is set kIterForward ( down beam)
    // or kIterBackward (beam up)
    if(!fInitOK) { Error("propagateToPCA()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}
    if(stateFrom.GetNrows() != 6) return kFALSE;
    if(stateTo  .GetNrows() != 6) stateTo.ResizeTo(6);
    return rkpropagator.propagateToPCA(chrg,stateTo,stateFrom,point,pathlength,direction);
}

Bool_t  HParticleMetaMatcher::getMDCStateRK(HParticleCand* cand,
					    TVectorD& st0,TVectorD& st1,TVectorD& st2,TVectorD& st3,
					    TVectorD& pathlength,Double_t mom,Int_t charge)
{
    // MDC hits recalculated using RK. sec coordinate system
    // if the optional arguments charge and mom
    // are not 0 they will replace the   charge and momentum from candidate.
    // This option can be used to estimate outer MDC hits for candidates without
    // reconstructed hits (sec coord. sys)
    // The track state at the 4 MDC are returned aswell as the path length to the MDC planes.
    // If MDC modules are not in the setup the cosponding hits a at 0,0,0 and the
    // tracklength for the module pathlength[m] = 0;
    if(!fInitOK) { Error("getMDCStateRK()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}

    TVectorD stateFrom(6);     // state at the origin layer
    TVectorD stateTo  (6);     // transported state at the destination layer
    if(!getTrackStateRK(cand,stateFrom,mom,charge)) return kFALSE;


    const Int_t s = cand->getSector();
    Int_t chrg;
    if(charge == 0) chrg = cand->getCharge();
    else            chrg = charge;

    st0.ResizeTo(6);   // output sec. coord. sys
    st1.ResizeTo(6);   // output sec. coord. sys
    st2.ResizeTo(6);   // output sec. coord. sys
    st3.ResizeTo(6);   // output sec. coord. sys

    pathlength.ResizeTo(4);  // output sec. coord. sys

    Double_t     L = 0;
    Double_t pathL = 0;
    for(Int_t m = 0; m < 4; m++){
	pathlength(m) = 0;
	if(mdcInstalled[s][m])
	{
	    HParticleRKPlane planeTo   = planeMdc[s][m+1];

	    if(!rkpropagator.propagateToPlane(chrg,stateTo,stateFrom,planeTo,pathL,kIterForward)){

		cout<<"broken RK start "<<centerVecRK[s].X()<<" "<<centerVecRK[s].Y()<<" "<<centerVecRK[s].Z()<<endl;
		cout<<"stateFrom m ="<<m<<" in "<<stateFrom(0)<<" , "<<stateFrom(1)<<" , "<<stateFrom(2)
		    <<" out "<<stateTo(0)<<" , "<<stateTo(1)<<" , "<<stateTo(2)
		    <<" s "<<s<<" phi "<<cand->getPhi()<<" theta "<<cand->getTheta()<<endl;
		return kFALSE;
	    }

	    L += pathL;

	    pathlength(m) = L;          // sec coord. sys

//#define debugRK

#ifdef debugRK

	    if(m<2 || cand->getPhi2()!=-1)
	    {
		Double_t x1,y1,z1, x2,y2,z2, x,y;
		if(m<2)HParticleTool::calcSegPoints(cand,x1,y1,z1,x2,y2,z2,kFALSE,kTRUE);
		else   HParticleTool::calcSegPoints(cand,x1,y1,z1,x2,y2,z2,kTRUE,kTRUE);
		(*fSizesCells)[s][m].calcInterTrMod(x1,y1,z1,x2,y2,z2,x,y);
		HGeomVector p(x,y,0);
		HGeomTransform trans =  *((*fSizesCells)[s][m].getSecTrans());
		HGeomVector pSec = trans.transFrom(p);   // sec coord. sys.


		cout<<"stateFrom m = "<<m<<" "<<setw(10)<<stateFrom(0)<<" , "<<setw(10)<<stateFrom(1)<<" , "<<setw(10)<<stateFrom(1)
		    <<" stateTo "<<setw(10)<<stateTo(0)<<" , "<<setw(10)<<stateTo(1)<<" , "<<setw(10)<<stateTo(2)
		    <<" psec "<<setw(10)<<pSec.getX()<<" , "<<setw(10)<<pSec.getY()<<" , "<<setw(10)<<pSec.getZ()
		    <<" s "<<s<<" phi "<<cand->getPhi()<<" theta "<<cand->getTheta()
		    <<" phi2 "<<cand->getPhi2()<<" theta2 "<<cand->getTheta2()<<endl;
	    }
#endif

	    if(m == 0) {
		st0 = stateTo;  // sec coord. sys.
	    } else if ( m == 1){
		st1 = stateTo;  // sec coord. sys.
	    } else if ( m == 2){
		st2 = stateTo;  // sec coord. sys.
	    } else if ( m == 3){
		st3 = stateTo;  // sec coord. sys.
	    }
	    stateFrom = stateTo;   // track or  sec coord. sys.
	}
    }

    return kTRUE;
}

Bool_t  HParticleMetaMatcher::getMDCPosRK(HParticleCand* cand,
					  TVector3& hit0,TVector3& hit1,TVector3& hit2,TVector3& hit3,
					  TVectorD& pathlength,Double_t mom,Int_t charge)
{
    // MDC hits recalculated using RK. sec coordinate system
    // if the optional arguments charge and mom
    // are not 0 they will replace the   charge and momentum from candidate.
    // This option can be used to estimate outer MDC hits for candidates without
    // reconstructed hits (sec coord. sys)
    // The hit points at the 4 MDC are returned aswell as the path length to the hits.
    // If MDC modules are not in the setup the cosponding hits a at 0,0,0 and the
    // tracklength for the module pathlength[m] = 0;

    hit0.SetXYZ(0.,0.,0);
    hit1.SetXYZ(0.,0.,0);
    hit2.SetXYZ(0.,0.,0);
    hit3.SetXYZ(0.,0.,0);
    if(!fInitOK) { Error("getMDCPOsRK()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}

    TVectorD st0(6);
    TVectorD st1(6);
    TVectorD st2(6);
    TVectorD st3(6);

    TVectorD zpos(4);

    if(!getMDCStateRK(cand,st0,st1,st2,st3,pathlength,mom,charge)) return kFALSE;   // sec coord. sys

    Int_t s = cand->getSector();

    if(mdcInstalled[s][0]){ hit0.SetXYZ(st0(0),st0(1),st0(2));}
    if(mdcInstalled[s][1]){ hit1.SetXYZ(st1(0),st1(1),st1(2));}
    if(mdcInstalled[s][2]){ hit2.SetXYZ(st2(0),st2(1),st2(2));}
    if(mdcInstalled[s][3]){ hit3.SetXYZ(st3(0),st3(1),st3(2));}

    return kTRUE;
}

Bool_t  HParticleMetaMatcher::getMDCPosMomRK(HParticleCand* cand,
					     TVector3& hit0,TVector3& hit1,TVector3& hit2,TVector3& hit3,
					     TVector3& mom0,TVector3& mom1,TVector3& mom2,TVector3& mom3,
					     TVectorD& pathlength,Double_t mom,Int_t charge)
{
    // MDC hits and dir recalculated using RK. sec coordinate system
    // if the optional arguments charge and mom
    // are not 0 they will replace the   charge and momentum from candidate.
    // This option can be used to estimate outer MDC hits for candidates without
    // reconstructed hits (sec coord. sys)
    // The hit points at the 4 MDC are returned aswell as the path length to the hits.
    // If MDC modules are not in the setup the cosponding hits a at 0,0,0 and the
    // tracklength for the module pathlength[m] = 0;

    hit0.SetXYZ(0.,0.,0);
    hit1.SetXYZ(0.,0.,0);
    hit2.SetXYZ(0.,0.,0);
    hit3.SetXYZ(0.,0.,0);
    mom0.SetXYZ(0.,0.,0);
    mom1.SetXYZ(0.,0.,0);
    mom2.SetXYZ(0.,0.,0);
    mom3.SetXYZ(0.,0.,0);
    if(!fInitOK) { Error("getMDCPosMom()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}

    TVectorD st0(6);
    TVectorD st1(6);
    TVectorD st2(6);
    TVectorD st3(6);

    if(!getMDCStateRK(cand,st0,st1,st2,st3,pathlength,mom,charge)) return kFALSE; // sec coord. sys

    Int_t s = cand->getSector();

    if(mdcInstalled[s][0]){ hit0.SetXYZ(st0(0),st0(1),st0(2)); mom0.SetXYZ(st0(3),st0(4),st0(5));}
    if(mdcInstalled[s][1]){ hit1.SetXYZ(st1(0),st1(1),st1(2)); mom1.SetXYZ(st1(3),st1(4),st1(5));}
    if(mdcInstalled[s][2]){ hit2.SetXYZ(st2(0),st2(1),st2(2)); mom2.SetXYZ(st2(3),st2(4),st2(5));}
    if(mdcInstalled[s][3]){ hit3.SetXYZ(st3(0),st3(1),st3(2)); mom3.SetXYZ(st3(3),st3(4),st3(5));}

    return kTRUE;
}

HParticleRKPlane  HParticleMetaMatcher::getPlaneMDC(Int_t s,Int_t m)
{
    // returns plane object for propagation with RK
    // In case sector or module is out of range an empth dummy object
    // is returned
    HParticleRKPlane pl;
    if(!fInitOK) { Error("getPlaneMDC()","You have to init HParticleMetaMatcher via the tasklist!"); return pl;}

    if(s<6&&s>=0 && m<4&&m>=0)
    {
           pl.setPlane(centerVecMdc[s][m],normVecMdc[s][m]);
	   return pl;

    } else {
        Error("getPlaneMDC()","sector or mod out of range, dummy object for s = %i, m=%i",s,m);
	return pl;
    }
}

HParticleRKPlane  HParticleMetaMatcher::getPlaneRK(Int_t s)
{
    // returns plane object for propagation with RK
    // in front of the first MDC
    // In case sector or module is out of range an empth dummy object
    // is returned
    HParticleRKPlane pl;
    if(!fInitOK) { Error("getPlaneRK()","You have to init HParticleMetaMatcher via the tasklist!"); return pl;}

    if(s<6&&s>=0)
    {
           pl.setPlane(centerVecRK[s],normVecRK[s]);
	   return pl;

    } else {
        Error("getPlaneRK()","sector of range, dummy object for s = %i",s);
	return pl;
    }
}

HParticleRKPlane HParticleMetaMatcher::getPlaneTOF(Int_t s,Int_t m)
{
    // returns plane object for propagation with RK
    // In case sector or module is out of range an empty dummy object is returned
    HParticleRKPlane pl;
    if(!fInitOK) { Error("getPlaneTOF()","You have to init HParticleMetaMatcher via the tasklist!"); return pl;}
    if(s<6&&s>=0 && m<8&&m>=0) {
        pl.setPlane(centerVecTof[s][m],normVecTof[s][m]);
	return pl;
    } else {
        Error("getPlaneTOF()","sector or mod out of range, dummy object for s = %i, m=%i",s,m);
	return pl;
    }
}

HParticleRKPlane HParticleMetaMatcher::getPlaneRPC(Int_t s,Int_t lay)
{
    // returns plane object for propagation with RK
    // lay   0 : col{0,2,4}, 1 : {1,3,5} 2: mid layer (for clusters using 1 cell in each layer)
    // In case sector or lay is out of range an empty dummy object is returned

    HParticleRKPlane pl;
    if(!fInitOK) { Error("getPlaneRPC()","You have to init HParticleMetaMatcher via the tasklist!"); return pl;}
    if(s<6&&s>=0 && lay<3&&lay>=0) {
        pl.setPlane(centerVecRpc[s][lay],normVecRpc[s][lay]);
	return pl;
    } else {
        Error("getPlaneRPC()","sector or lay out of range, dummy object for s = %i, l=%i",s,lay);
	return pl;
    }
}

HParticleRKPlane HParticleMetaMatcher::getPlaneEMC(Int_t s)
{
    // returns plane object for propagation with RK
    // In case sector is out of range an empty dummy object is returned

    HParticleRKPlane pl;
    if(!fInitOK) { Error("getPlaneEMC()","You have to init HParticleMetaMatcher via the tasklist!"); return pl;}
    if(s<6&&s>=0) {
        pl.setPlane(centerVecEmc[s],normVecEmc[s]);
	return pl;
    } else {
        Error("getPlaneEMC()","sector out of range, dummy object for s = %i",s);
	return pl;
    }
}

HParticleRKPlane HParticleMetaMatcher::getPlaneiTOF(Int_t s)
{
    // returns plane object for propagation with RK (SEC SYS)
    // In case sector is out of range an empty dummy object is returned
    HParticleRKPlane pl;
    if(!fInitOK) { Error("getPlaneiTOF()","You have to init HParticleMetaMatcher via the tasklist!"); return pl;}
    if(s<6&&s>=0) {
        pl.setPlane(centerVeciTof[s],normVeciTof[s]);
	return pl;
    } else {
        Error("getPlaneiTOF()","sector or mod out of range, dummy object for s = %i",s);
	return pl;
    }
}

Bool_t HParticleMetaMatcher::propagateToVertexStateRK(HParticleCand* cand,
						      TVectorD& stStart ,TVectorD& stPCA ,
						      const TVector3& vertexLab,Double_t& pathlength,Double_t mom,Int_t charge)
{
    // calculates track states (x,y,z,px,py,pz [MeV/c, mm, sec . coord. sys])
    // at RK startPlane and propagates it using rungekutta to the point of
    // closest approach to vertexLab (in Lab coord. sys. will be tranformed to sec. sys. internally).
    // if the optional arguments charge and mom
    // are not 0 they will replace the   charge and momentum from candidate.
    // This option can be used to estimate candidates without
    // reconstructed outer segment (charge) and momentum.
    stStart.ResizeTo(6);     // state at the origin layer
    stPCA  .ResizeTo(6);     // transported state at point of closest approach (PCA) to vertex
    pathlength = 0;
    if(!fInitOK) { Error("propagateToVertexStateRK()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}
    if(!getTrackStateRK(cand,stStart,mom,charge)) return kFALSE;

    Int_t chrg;
    if(charge == 0) chrg = cand->getCharge();
    else            chrg = charge;

    TVector3 vertexSec;
    pointLabToSecSys(cand->getSector(),vertexLab,vertexSec);
    return rkpropagator.propagateToPCA(chrg,stStart,stPCA,vertexSec,pathlength,kIterBackward);

}

Bool_t HParticleMetaMatcher::propagateToVertexPosRK(HParticleCand* cand,
						    TVector3& posStart ,TVector3& posPCA ,
						    const TVector3& vertexLab,Double_t& pathlength,Double_t mom,Int_t charge)
{
    // calculates track position (x,y,z [ mm, sec . coord. sys])
    // at RK startPlane and propagates it using rungekutta to the point of
    // closest approach to vertexLab (in Lab coord. sys. will be tranformed to sec. sys. internally) posPCA.
    // if the optional arguments charge and mom
    // are not 0 they will replace the   charge and momentum from candidate.
    // This option can be used to estimate candidates without
    // reconstructed outer segment (charge) and momentum.
    posStart.SetXYZ(0,0,0);
    posPCA  .SetXYZ(0,0,0);

    TVectorD stStart(6);     // state at the origin layer
    TVectorD stPCA(6);       // transported state at point of closest approach (PCA) to vertex
    pathlength = 0;
    if(!fInitOK) { Error("propagateToVertexRK()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}
    if(!getTrackStateRK(cand,stStart,mom,charge)) return kFALSE;

    Int_t chrg;
    if(charge == 0) chrg = cand->getCharge();
    else            chrg = charge;

    TVector3 vertexSec;
    pointLabToSecSys(cand->getSector(),vertexLab,vertexSec);
    if(!rkpropagator.propagateToPCA(chrg,stPCA,stStart,vertexSec,pathlength,kIterBackward)) return kFALSE;

    posStart.SetXYZ(stStart(0),stStart(1),stStart(2));
    posPCA  .SetXYZ(stPCA  (0),stPCA  (1),stPCA  (2));

    return kTRUE;
}

Bool_t HParticleMetaMatcher::propagateToVertexPosMomRK(HParticleCand* cand,
						       TVector3& posStart,TVector3& posPCA,TVector3& momStart,TVector3& momPCA,
						       const TVector3& vertexLab,Double_t& pathlength,Double_t mom,Int_t charge)
{
    // calculates track position (x,y,z [ mm, sec . coord. sys])  and momentum (px,py,pz [MeV/c])
    // at RK startPlane and propagates it using rungekutta to the point of
    // closest approach to vertexLab (in Lab coord. sys. will be tranformed to sec. sys. internally) posPCA.
    // if the optional arguments charge and mom
    // are not 0 they will replace the   charge and momentum from candidate.
    // This option can be used to estimate candidates without
    // reconstructed outer segment (charge) and momentum.
    posStart.SetXYZ(0,0,0);
    posPCA  .SetXYZ(0,0,0);
    momStart.SetXYZ(0,0,0);
    momPCA  .SetXYZ(0,0,0);

    TVectorD stStart(6);     // state at the origin layer
    TVectorD stPCA  (6);     // transported state at point of closest approach (PCA) to vertex
    pathlength = 0;
    if(!fInitOK) { Error("propagateToVertexPosMomRK()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}
    if(!getTrackStateRK(cand,stStart,mom,charge)) return kFALSE;

    Int_t chrg;
    if(charge == 0) chrg = cand->getCharge();
    else            chrg = charge;

    TVector3 vertexSec;
    pointLabToSecSys(cand->getSector(),vertexLab,vertexSec);
    if(!rkpropagator.propagateToPCA(chrg,stPCA,stStart,vertexSec,pathlength,kIterBackward)) return kFALSE;

    posStart.SetXYZ(stStart(0),stStart(1),stStart(2));
    posPCA  .SetXYZ(stPCA  (0),stPCA  (1),stPCA  (2));

    momStart.SetXYZ(stStart(3),stStart(4),stStart(5));
    momPCA  .SetXYZ(stPCA  (3),stPCA  (4),stPCA  (5));

    return kTRUE;
}




Bool_t HParticleMetaMatcher::propagateToVertexStateLine(HParticleCand* cand,
							TVectorD& stStart ,TVectorD& stPCA ,
							const TVector3& vertexLab,Double_t& pathlength)
{
    // calculates track states (x,y,z,px,py,pz [MeV/c, mm, sec . coord. sys])
    // at RK startPlane and propagates it using straight line to the point of
    // closest approach to vertexLab (in Lab coord. sys. will be tranformed to sec. sys. internally).
    stStart.ResizeTo(6);     // state at the origin layer
    stPCA  .ResizeTo(6);     // transported state at point of closest approach (PCA) to vertex
    pathlength = 0;
    if(!fInitOK) { Error("propagateToVertexStateLine()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}

    Double_t mom = cand->getMomentum();
    Int_t chrg   = cand->getCharge();
    Bool_t fake = kFALSE;
    if(mom < 0)   { mom  = 1000.;  fake = kTRUE;}
    if(chrg == 0) { chrg = 1;      fake = kTRUE;}

    if(!getTrackStateRK(cand,stStart,mom,chrg)) return kFALSE;

    TVector3 p1,p2,dir;
    HParticleTool::calcSegPoints(cand,p1,p2,kFALSE,kTRUE);

    if(fake){ // no mom or outer seg was used, define the dir by seg points
	TVector3 dir = p1-p2;
	stStart(3) = dir.X();
	stStart(4) = dir.Y();
	stStart(5) = dir.Z();
    }

    TVector3 pca,pos;
    Int_t    pcaFlag = 0;

    pos.SetXYZ(stStart(0),stStart(1),stStart(2)); // start pos of RK seg

    TVector3 vertexSec;
    pointLabToSecSys(cand->getSector(),vertexLab,vertexSec);
    HKalGeomTools::distancePointToLine(pca, pcaFlag, vertexSec,p1, p2);
    TVector3 L = pca - pos;
    pathlength = L.Mag();

    stPCA(0) = pca.X();
    stPCA(1) = pca.Y();
    stPCA(2) = pca.Z();
    stPCA(3) = stStart(3);
    stPCA(4) = stStart(4);
    stPCA(5) = stStart(5);

    return kTRUE;

}

Bool_t HParticleMetaMatcher::propagateToVertexPosLine(HParticleCand* cand,
						      TVector3& posStart ,TVector3& posPCA ,
						      const TVector3& vertexLab,Double_t& pathlength)
{
    // calculates track position (x,y,z [ mm, sec . coord. sys])
    // at RK startPlane and propagates it using straight line to the point of
    // closest approach to vertexLab (in Lab coord. sys. will be tranformed to sec. sys. internally) posPCA.
    TVectorD stStart(6);     // state at the origin layer
    TVectorD stPCA  (6);     // transported state at point of closest approach (PCA) to vertex
    pathlength = 0;
    posStart.SetXYZ(0,0,0);
    posPCA  .SetXYZ(0,0,0);
    if(!fInitOK) { Error("propagateToVertexPosLine()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}

    TVector3 vertexSec;
    pointLabToSecSys(cand->getSector(),vertexLab,vertexSec);
    if(!propagateToVertexStateLine(cand,stStart,stPCA,vertexSec,pathlength)) return kFALSE;

    posStart.SetXYZ(stStart(0),stStart(1),stStart(2));
    posPCA  .SetXYZ(stPCA  (0),stPCA  (1),stPCA  (2));

    return kTRUE;

}

Bool_t HParticleMetaMatcher::propagateToVertexPosMomLine(HParticleCand* cand,
							 TVector3& posStart ,TVector3& posPCA , TVector3& momStart ,TVector3& momPCA ,
							 const TVector3& vertexLab,Double_t& pathlength)
{
    // calculates track position (x,y,z [ mm, sec . coord. sys]) and mom (px,py,pz [MeV/c])
    // at RK startPlane and propagates it using straight line to the point of
    // closest approach to vertexLab (in Lab coord. sys. will be tranformed to sec. sys. internally) posPCA.
    TVectorD stStart(6);     // state at the origin layer
    TVectorD stPCA  (6);     // transported state at point of closest approach (PCA) to vertex
    pathlength = 0;
    posStart.SetXYZ(0,0,0);
    posPCA  .SetXYZ(0,0,0);
    momStart.SetXYZ(0,0,0);
    momPCA  .SetXYZ(0,0,0);
    if(!fInitOK) { Error("propagateToVertexPosMomLine()","You have to init HParticleMetaMatcher via the tasklist!"); return kFALSE;}

    TVector3 vertexSec;
    pointLabToSecSys(cand->getSector(),vertexLab,vertexSec);
    if(!propagateToVertexStateLine(cand,stStart,stPCA,vertexSec,pathlength)) return kFALSE;

    posStart.SetXYZ(stStart(0),stStart(1),stStart(2));
    posPCA  .SetXYZ(stPCA  (0),stPCA  (1),stPCA  (2));
    momStart.SetXYZ(stStart(3),stStart(4),stStart(5));
    momPCA  .SetXYZ(stPCA  (3),stPCA  (4),stPCA  (5));

    return kTRUE;

}

Bool_t HParticleMetaMatcher::loopFlagDoubleUsedMetaCells(vector<HParticleMetaPredict>& vpredict,UInt_t print)
{
    // loops over all HParticleCand and flags candidates by c->setFlagBit(Particle::kIsMetaDoubleUse)
    // if more than 1 candidate points to the same META cell (TOF or RPC) (even if the candidate
    // had no META hit hatched to it). Matching with META is done via out MDC Seg using a straight
    // line approach or RK extrapolation from inner Mdc Seg (setting by setUseRK()).
    // print option  bitwise  1  prediction           dec  1
    //                        2  real cand            dec  2
    //                        3  skip flag            dec  4
    //                        4  list flag            dec  8
    //                        5  missmatch            dec 16
    //                        6  stats per event      dec 32
    //                        7  stats cumulative     dec 64
    //                        print = 1+2+4+8+16+32+64;  // all prints
    vpredict.clear();
    if(fCatParticleCand)
    {

	UInt_t n = fCatParticleCand->getEntries();
	HParticleCand* cand = 0;
	Int_t mod0,cell0,mod1,cell1;
	Double_t path0, path1;
        HGeomVector hit0,hit1;
        Int_t sys = 1; // mod sys

        map<UInt_t, vector<UInt_t> > mTofUse;    // Prediction tof lin addr -> vector of ind of vector vpredict
        map<UInt_t, vector<UInt_t> > mRpcUse;    // Prediction rpc lin addr -> vector of ind of vector vpredict
        map<UInt_t, vector<UInt_t> > mTofMetaUse;// Real tof lin addr -> vector of ind of vector vpredict
        map<UInt_t, vector<UInt_t> > mRpcMetaUse;// Real rpc lin addr -> vector of ind of vector vpredict

	map<UInt_t, vector<UInt_t> >::iterator it;

	for(UInt_t i = 0; i < n; i++)
	{
	    cand = HCategoryManager::getObject(cand,fCatParticleCand,i);
	    if(cand){
		if(cand->getOuterSegInd() == -1) continue;
		cand->unsetFlagBit(Particle::kIsMetaDoubleUse);  // clean flag from before
		Int_t s = cand->getSector();

		HParticleMetaPredict p;
                Int_t ad = -1;
		if(predictTofCell(cand,hit0,hit1,   s,mod0,cell0,mod1,cell1,  path0,path1,  sys))
		{  // theoretical hit found on TOF
		    p.set(cand,1,0 ,s,mod0,cell0,hit0,path0);
                    p.set(cand,1,1 ,s,mod1,cell1,hit1,path1);

		    if( (print >> 0) & 0x01) {
		    cout<<"TOF pred seg ind "<<setw(6)<<cand->getOuterSegInd()<<" vec ind "<<setw(6)<<vpredict.size()<<" a1 "<<setw(6)<<p.addrTof[0]<<" a2 "<<setw(6)<<p.addrTof[1]
		        <<" m0 "<<setw(6)<<mod0<<" c0 "<<setw(6)<<cell0<<" m1 "<<setw(6)<<mod1<<" c1 "<<setw(6)<<cell1<<endl;
		    }

		    for(Int_t a = 0; a < 2; a++){
			ad = p.addrTof[a];
			if(ad != -1){
			    it = mTofUse.find(ad);
			    if(it == mTofUse.end()){ // new TOF address
				vector<UInt_t> v;
				v.push_back(vpredict.size());
				mTofUse[ad] = v;
			    } else {
				it->second.push_back(vpredict.size());
			    }
			}
		    }
		}
                
                if(predictRpcCell(cand,hit0,hit1,   s,mod0,cell0,mod1,cell1,  path0,path1,  sys))
		{  // theoretical hit found on RPC
		    p.set(cand,0,0 ,s,mod0,cell0,hit0,path0);
                    p.set(cand,0,1 ,s,mod1,cell1,hit1,path1);

		    if( (print >> 0) & 0x01) {
			cout<<"RPC pred seg ind "<<setw(6)<<cand->getOuterSegInd()<<" vec ind "<<setw(6)<<vpredict.size()<<" a1 "<<setw(6)<<p.addrRpc[0]<<" a2 "<<setw(6)<<p.addrRpc[1]
			    <<" m0 "<<setw(6)<<mod0<<" c0 "<<setw(6)<<cell0<<" m1 "<<setw(6)<<mod1<<" c1 "<<setw(6)<<cell1<<endl;
		    }

		    for(Int_t a = 0; a < 2; a++){
			ad = p.addrRpc[a];
			if(ad != -1){
			    it = mRpcUse.find(ad);
			    if(it == mRpcUse.end()){ // new RPC address
				vector<UInt_t> v;
				v.push_back(vpredict.size());
				mRpcUse[ad] = v;
			    } else {
				it->second.push_back(vpredict.size());
			    }
			}
		    }
		}

		if(p.c){

		    //-------------------------------------------------------------------
                    // check match between asigned META cells and predicted META cells
		    if(p.c->isTofClstUsed() || p.c->isTofHitUsed() ) {

			for(Int_t a = 0; a < 2; a++){
			    ad = p.addrMeta[a];
			    if(ad != -1){
				it = mTofMetaUse.find(ad);
				if(it == mTofMetaUse.end()){ // new TOF address
				    vector<UInt_t> v;
				    v.push_back(vpredict.size());
				    mTofMetaUse[ad] = v;
				} else {
				    it->second.push_back(vpredict.size());
				}
			    }
			}

			recalcTof(p.c);
			p.meta    = getTrackMetaMod();
			p.metaSec = getTrackMetaSec();
                        if(p.c->getMetaModule(0) != -1) p.addrMeta[0] = HTool::getLinearIndex(s,6,p.c->getMetaModule(0),10,p.c->getMetaCell(0),100);
                        if(p.c->getMetaModule(1) != -1) p.addrMeta[1] = HTool::getLinearIndex(s,6,p.c->getMetaModule(1),10,p.c->getMetaCell(1),100);

			if( (print >> 1) & 0x01) {
			    cout<<"TOF real vec ind "<<setw(6)<<vpredict.size()
				<<" a1 "<<setw(6)<<p.addrMeta[0]<<" a2 "<<setw(6)<<p.addrMeta[1]
				<<" pred a1 "<<setw(6)<<p.addrTof[0]<<" a2 "<<setw(6)<<p.addrTof[1]<<endl;
			}


		    }  else if(p.c->isRpcClstUsed()) {

			for(Int_t a = 0; a < 2; a++){
			    ad = p.addrMeta[a];
			    if(ad != -1){
				it = mRpcMetaUse.find(ad);
				if(it == mRpcMetaUse.end()){ // new RPC address
				    vector<UInt_t> v;
				    v.push_back(vpredict.size());
				    mRpcMetaUse[ad] = v;
				} else {
				    it->second.push_back(vpredict.size());
				}
			    }
			}

 			recalcRpc(p.c);
                        p.meta    = getTrackMetaMod();
                        p.metaSec = getTrackMetaSec();
                        if(p.c->getMetaModule(0) != -1) p.addrMeta[0] = HTool::getLinearIndex(s,6,p.c->getMetaModule(0),10,p.c->getMetaCell(0),100);
                        if(p.c->getMetaModule(1) != -1) p.addrMeta[1] = HTool::getLinearIndex(s,6,p.c->getMetaModule(1),10,p.c->getMetaCell(1),100);

			if( (print >> 1) & 0x01) {
			    cout<<"RPC real vec ind "<<setw(6)<<vpredict.size()
				<<" a1 "<<setw(6)<<p.addrMeta[0]<<" a2 "<<setw(6)<<p.addrMeta[1]
				<<" pred a1 "<<setw(6)<<p.addrRpc[0]<<" a2 "<<setw(6)<<p.addrRpc[1]<<endl;
			}
		    }
                    p.metaMatchPredictAll();
		    vpredict.push_back(p);
		}

	    } // end cand
	} // end cand loop

        //-------------------------------------------------------------------
	// finally flag the candidates which point to
        // same META cell
        static Int_t flagTofS   = 0;
        static Int_t flagRpcS   = 0;
        static Int_t matchTofS  = 0;
        static Int_t matchRpcS  = 0;
        static Int_t matchTofOneOffS = 0;
        static Int_t matchRpcOneOffS = 0;
        static Int_t nCandS     = 0;
        static Int_t nCandMetaS = 0;
        Int_t nCandMeta = 0;
        Int_t flagTof   = 0;
        Int_t flagRpc   = 0;
        Int_t matchTof  = 0;
        Int_t matchRpc  = 0;
        Int_t matchTofOneOff  = 0;
        Int_t matchRpcOneOff  = 0;

	map<Int_t,Int_t> mCt;
	map<Int_t,Int_t>::iterator itCt;

	for (it = mTofUse.begin(); it != mTofUse.end(); ++it) {
	    vector<UInt_t>& v = it->second;
	    if(v.size() > 1){ // same TOF Cell hit more than once

		mCt.clear();
		for(UInt_t i = 0 ; i < v.size(); i ++){ // find unique out seg matches
		    HParticleMetaPredict& p = vpredict[v[i]];
		    itCt = mCt.find(p.segInd[1]);
		    if(itCt == mCt.end()){
                       mCt[p.segInd[1]] = 1;
		    } else {
                       mCt[p.segInd[1]] ++;
		    }
		}

		if(mCt.size() == 1)  {
                    itCt = mCt.begin();
                    if( (print >> 2) & 0x01) cout<<"TOF skip outer seg "<<setw(6)<<itCt->first<<" count "<<setw(6)<<itCt->second<<endl;
		    continue; // only cands with same outer seg
		}
		// more than 1 unique outer seg : flag all candidates
		for(UInt_t i = 0 ; i < v.size(); i ++){
		    HParticleMetaPredict& p = vpredict[v[i]];
                    if( (print >> 3) & 0x01) cout<<"TOF list addr "<<setw(6)<<it->first<<" vec ind "<<setw(6)<<v[i]<<" total "<<setw(6)<<v.size()<<" i "<<setw(6)<<i<<endl;
		    p.c->setFlagBit(Particle::kIsMetaDoubleUse);
		    flagTof++;
		}
	    }
	}

        //-------------------------------------------------------------------
	for (it = mRpcUse.begin(); it != mRpcUse.end(); ++it) {
	    vector<UInt_t>& v = it->second;
	    if(v.size() > 1){ // same RPC Cell hit more than once
		mCt.clear();
		for(UInt_t i = 0 ; i < v.size(); i ++){ // find unique out seg matches
		    HParticleMetaPredict& p = vpredict[v[i]];
		    itCt = mCt.find(p.segInd[1]);
		    if(itCt == mCt.end()){
                       mCt[p.segInd[1]] = 1;
		    } else {
                       mCt[p.segInd[1]] ++;
		    }
		}

		if(mCt.size() == 1)  {
                    itCt = mCt.begin();
                    if( (print >> 2) & 0x01) cout<<"RPC skip flag outer seg "<<setw(6)<<itCt->first<<" count "<<setw(6)<<itCt->second<<endl;
		    continue; // only cands with same outer seg
		}

		// more than 1 unique outer seg : flag all candidates
		for(UInt_t i = 0 ; i < v.size(); i ++){
		    HParticleMetaPredict& p = vpredict[v[i]];
                    if( (print >> 3) & 0x01) cout<<"RPC flag addr "<<setw(6)<<it->first<<" seg ind "<<setw(6)<<p.segInd[1]<<" vec ind "<<setw(6)<<v[i]<<" total "<<setw(6)<<v.size()<<" i "<<setw(6)<<i<<endl;
                    p.c->setFlagBit(Particle::kIsMetaDoubleUse);
                    flagRpc++;
		}
	    }
	}


        //-------------------------------------------------------------------
        // copy indices of other matches to META cells
	Int_t s [2],m [2],c [2];
	Int_t sM[2],mM[2],cM[2];

	for (UInt_t i = 0 ; i < vpredict.size(); i++){
	    HParticleMetaPredict& p = vpredict[i];
	    if(p.c->getBeta() > 0) {
                nCandMeta++;


		if(p.metaMatchPredict()){
		} else {
		    if(p.c->isTofClstUsed() || p.c->isTofHitUsed() ) {
                        matchTof++;
			p.convertAddr(p.addrTof [0], s[0], m[0], c[0]);
			p.convertAddr(p.addrTof [1], s[1], m[1], c[1]);
			p.convertAddr(p.addrMeta[0],sM[0],mM[0],cM[0]);
			p.convertAddr(p.addrMeta[1],sM[1],mM[1],cM[1]);

			Bool_t found = kFALSE;
			for(UInt_t j = 0; j < 2; j++){
                            if(found) break;
			    if(p.addrTof [j] == -1 ) continue;
			    for(UInt_t k = 0; k < 2; k++){
				if(p.addrMeta[k] == -1) continue;
                                Int_t linRod  = m [j]*8 + c [j];
				Int_t linRodM = mM[j]*8 + cM[j];
				if(abs(linRod - linRodM) == 1) { matchTofOneOff++; found = kTRUE; break;}  // rod +-1
			    }
			}

			if( (print >> 4) & 0x01) {
			    Float_t distLeft,distRight,distLow,distUp;
			    distToEdgeTofCell(p.meta,sM[0],mM[0],cM[0], distLeft,distRight,distLow,distUp);

			    cout<<"TOF missmatch real/pred : seg ind "<<setw(4)<<p.c->getOuterSegInd()
				<<" vec ind "<<setw(4)<<i
				<<" a1 "     <<setw(4)<<p.addrTof[0] <<" a2 "<<setw(4)<<p.addrTof[1]
				<<" RPC a1 " <<setw(4)<<p.addrRpc[0] <<" a2 "<<setw(4)<<p.addrRpc[1]
				<<" Real a1 "<<setw(4)<<p.addrMeta[0]<<" a2 "<<setw(4)<<p.addrMeta[1]
				<<" m0,c0 ("  <<setw(2)<<m[0]<<","<<setw(2)<<c[0]<<") m1,c1 ("<<setw(2)<<m[1]<<","<<setw(2)<<c[1]<<") "
				<<" (x,y) 0: "<<setw(10)<<p.hitTof[0].X()<<","<<setw(10)<<p.hitTof[0].Y()
				<<" 1: "      <<setw(10)<<p.hitTof[1].X()<<","<<setw(10)<<p.hitTof[1].Y()
				<<" real: "   <<setw(10)<<p.meta.X()<<","<<setw(10)<<p.meta.Y()
				<<" m0,c0 ("  <<setw(2)<<mM[0]<<","<<setw(2)<<cM[0]<<") m1,c1 ("<<setw(2)<<mM[1]<<","<<setw(2)<<cM[1]<<")"
                                <<" one Off "<<(Int_t)found
				<<endl;
			    cout<<" dist Left "<<setw(10)<<distLeft<<" right "<<setw(10)<<distRight<<" Low "<<setw(10)<<distLow<<" up "<<setw(10)<<distUp<<endl;
			}
		    } else if(p.c->isRpcClstUsed()){
			matchRpc++;
			p.convertAddr(p.addrRpc [0], s[0], m[0], c[0]);
			p.convertAddr(p.addrRpc [1], s[1], m[1], c[1]);
			p.convertAddr(p.addrMeta[0],sM[0],mM[0],cM[0]);
			p.convertAddr(p.addrMeta[1],sM[1],mM[1],cM[1]);

			Bool_t found = kFALSE;
			for(UInt_t j = 0; j < 2; j++){
                            if(found) break;
			    if(p.addrRpc [j] == -1 ) continue;
			    for(UInt_t k = 0; k < 2; k++){
				if(p.addrMeta[k] == -1) continue;
                                Int_t col  = (m [j] + 2 - ((m [j]%2)*1))/3;    // col 0,1,2 for both planes front/back
                                Int_t colM = (mM[j] + 2 - ((mM[j]%2)*1))/3;

				if(abs(p.addrRpc[j] - p.addrMeta[k]) == 1)  { matchRpcOneOff++; found = kTRUE;  break;}    // cell +- 1
				if( col == colM  && c[j] == cM[k])          { matchRpcOneOff++; found = kTRUE;  break;}    // col front/back aligned  +  same cell
				if( col == colM  && abs(c[j] - cM[k]) == 1) { matchRpcOneOff++; found = kTRUE;  break;}    // col front/back aligned  +  cell +-1
			    }
			}

			if( (print >> 4) & 0x01) {
			    Float_t distLeft,distRight,distLow,distUp;
			    distToEdgeRpcCell(p.meta,sM[0],mM[0],cM[0], distLeft,distRight,distLow,distUp);

			    cout<<"RPC missmatch real/pred : seg ind "<<setw(4)<<p.c->getOuterSegInd()
				<<" vec ind "<<setw(4)<<i
				<<" a1 "     <<setw(4)<<p.addrRpc[0]<<" a2 "<<setw(4)<<p.addrRpc[1]
				<<" TOF a1 " <<setw(4)<<p.addrTof[0]<<" a2 "<<setw(4)<<p.addrTof[1]
				<<" Real a1 "<<setw(4)<<p.addrMeta[0]<<" a2 "<<setw(4)<<p.addrMeta[1]
				<<" m0,c0 ("  <<setw(2)<<m[0]<<","<<setw(2)<<c[0]<<") m1,c1 ("<<setw(2)<<m[1]<<","<<setw(2)<<c[1]<<") "
				<<" (x,y) 0: "<<setw(10)<<p.hitRpc[0].X()<<","<<setw(10)<<p.hitRpc[0].Y()
				<<" 1: "      <<setw(10)<<p.hitRpc[1].X()<<","<<setw(10)<<p.hitRpc[1].Y()
				<<" real: "   <<setw(10)<<p.meta.X()<<","<<setw(10)<<p.meta.Y()
				<<" m0,c0 ("  <<setw(2)<<mM[0]<<","<<setw(2)<<cM[0]<<") m1,c1 ("<<setw(2)<<mM[1]<<","<<setw(2)<<cM[1]<<")"
                                <<" one Off "<<(Int_t)found
				<<endl;
			    cout<<" dist Left "<<setw(10)<<distLeft<<" right "<<setw(10)<<distRight<<" Low "<<setw(10)<<distLow<<" up "<<setw(10)<<distUp<<endl;
			}
		    }
		} // mismatched

	    } //beta

	    for(Int_t j = 0; j < p.nTof; j++)
	    {
		Int_t ad = p.addrTof [j];
                if(ad == -1) continue;
		it = mTofUse.find(ad);
		if(it != mTofUse.end()){
		    vector<UInt_t>& v = it->second;
		    for(UInt_t k = 0; k < v.size(); k++){
                      if(v[k] != i) p.vcandTof.push_back(v[k]);
		    }
		}
	    }

	    for(Int_t j = 0; j < p.nRpc; j++)
	    {
		Int_t ad = p.addrRpc [j];
		if(ad == -1) continue;
		it = mRpcUse.find(ad);
		if(it != mRpcUse.end()){
		    vector<UInt_t>& v = it->second;
		    for(UInt_t k = 0; k < v.size(); k++){
                      if(v[k] != i) p.vcandRpc.push_back(v[k]);
		    }
		}
	    }

	} // end copy

        //-------------------------------------------------------------------
        // remove double entries
	for (UInt_t i = 0 ; i < vpredict.size(); i++){
	    HParticleMetaPredict& p = vpredict[i];

	    UInt_t nT = p.vcandTof.size();
	    UInt_t nR = p.vcandRpc.size();
            if(nT < 2 && nR < 2)   continue; // nothing to do
	    sort(p.vcandTof.begin(),p.vcandTof.end());
	    vector<UInt_t>::iterator it2 = unique(p.vcandTof.begin(),p.vcandTof.end());
            p.vcandTof.erase(it2, p.vcandTof.end());

	    sort(p.vcandRpc.begin(),p.vcandRpc.end());
	    it2 = unique(p.vcandRpc.begin(),p.vcandRpc.end());
            p.vcandRpc.erase(it2, p.vcandRpc.end());

	}


        flagTofS   += flagTof;
        flagRpcS   += flagRpc;
        matchTofS  += matchTof;
        matchRpcS  += matchRpc;
        matchTofOneOffS += matchTofOneOff;
        matchRpcOneOffS += matchRpcOneOff;
        nCandS     += vpredict.size();
        nCandMetaS += nCandMeta;

	if( (print >> 5) & 0x01) {
	    cout<<"--------------------------------------------------------------"<<endl;
	    cout<<"HParticleMetaMatcher::loopFlagDoubleUsedMetaCells()  : PER EVENT   : candidates inspected : "
		<<setw(6)<< vpredict.size()<<" nCandMeta "<<nCandMeta
		<<" mismatchedTof : "<<setw(6)<<matchTof<<" one Off Tof : "<<setw(6)<<matchTofOneOff
		<<" mismatchedRpc : "<<setw(6)<<matchRpc<<" one Off Rpc : "<<setw(6)<<matchRpcOneOff
		<<" : flagged TOF "<<setw(6)<<flagTof<<" flagged RPC "<<setw(6)<<flagRpc<<endl;
	    cout<<"--------------------------------------------------------------"<<endl;
	}

	if( (print >> 6) & 0x01) {
	    cout<<"--------------------------------------------------------------"<<endl;
	    cout<<"HParticleMetaMatcher::loopFlagDoubleUsedMetaCells()  : ACCUMULATED : candidates inspected : "
		<<setw(6)<< nCandS<<" nCandMeta "<<nCandMetaS
		<<" mismatchedTof : "<<setw(6)<<matchTofS<<" one Off Tof : "<<setw(6)<<matchTofOneOffS
		<<" mismatchedRpc : "<<setw(6)<<matchRpcS<<" one Off Rpc : "<<setw(6)<<matchRpcOneOffS
		<<" : flagged TOF "<<setw(6)<<flagTofS<<" flagged RPC "<<setw(6)<<flagRpcS<<endl;
	    cout<<"--------------------------------------------------------------"<<endl;
	}
	//-------------------------------------------------------------------




    } // if catParticleCand

    return kTRUE;
}