//*--- Author: Jochen Markert

// from ROOT
#include "TF2.h"
#include "TRandom.h"

// form hydra
#include "hcategory.h"
#include "hgeantkine.h"
#include "hgeantmdc.h"
#include "hlocation.h"
#include "hmdccal1.h"
#include "hmdccal2par.h"
#include "hmdccal2parsim.h"
#include "hmdcgetcontainers.h"
#include "hmagnetpar.h"
#include "hmdcseg.h"
#include "hmdcidealclasses.h"
#include "hmdcsizescells.h"
#include "hmdctrackgfieldpar.h"
#include "hmdctrkcand.h"
#include "hmdcidealclasses.h"
#include "hparoraio.h"
#include "hrootsource.h"
#include "hruntimedb.h"
#include "hmdcdetector.h"
#include "htofdetector.h"
#include "hshowerdetector.h"
#include "hspectrometer.h"

#include "hkaldef.h"
#include "hkaldetcradle.h"
#include "hkalinput.h"
#include "hkalmdchit.h"


#include <iostream>
using namespace std;

ClassImp(HKalInput)

//_HADES_CLASS_DESCRIPTION
///////////////////////////////////////////////////////////////////////////////
//    HKalMdcImput is a helper class to fill a TObjArray of hits in sector
//    coordinate system to be used in the Kalman filter. Input can be provided
//    from HMdcTrkCand (real tracking) or HMdcTrkCandIdeal (ideal tracking from
//    from GEANT hit points). It will fill 4 (2x2) hits from inner/outer segment.
//    Errors for real tracking are taken from HMdcSeg, for ideal tracking 200/100 mu
//    in x/y. The errors can be scaled by a factor. The geometrical transformations
//    are performed via HMdcSizesCells, thus the container has to initialized
//    before usage.
//
//-----------------------------------------------------------------------
//    HOWTO use the HKalMdcInput class::
//    before Event loop:
//
//    // assume Spectrometer setup + Hades creation is performed ....
//
//    HKalDetCradle *pCradleHades = new HKalDetCradle();
//    HKalMdcInput *input = new HKalMdcInput(refID, pCradleHades);   // fill the input hits
//    input->setPrint(0);                       // debug print (off == 0)
//
//    TObjArray allhits;                        // array of hits to be fed into kalman
//    Int_t measDim = 2;                        // dimension of hit vector (default = 2)
//    Double_t scaleErr = 1.;                   // scaling factor for errors (default = 1.)
//
//    Inside eventloop:
//
//        //-----------------------------------------------------------------------
//        // Run Kalman filter on hits.
//        Int_t nrHits = 0;
//        input->resetIterMdcTrkCand();
//	  	  while( (nrHits = input->nextMdcTrackCand(allhits, measDim, scaleErr)) )
//        {
//	       // ..... kalman
//
//	      }
//-----------------------------------------------------------------------
///////////////////////////////////////////////////////////////////////////////

//  -----------------------------------
//  Ctors and Dtor
//  -----------------------------------

HKalInput::HKalInput() : TObject() {
    // It is recommended to use the constructor
    // HKalInput(Int_t refID, HKalDetCradle *detCradle)
    // instead where the object will be initialized correctly.
    // This empty constructor sets all pointers to NULL. To use the object
    // the following still has to be done afterwards:
    // Do spectrometer setup + Hades creation ....
    // To set the pointer to detector cradle:
    // setDetectorCradle(detCradle)  // set pointer to detector cradle
    // To get pointer to categories, HMdcSizesCells, magnet and B-field
    // parameters:
    // init(refID)                   // gets categories

    pCradleHades = NULL;
    fSizesCells  = NULL;
    pField       = NULL;
    pMagnet      = NULL;
    setPrint(kFALSE);
}

HKalInput::HKalInput(Int_t refID, HKalDetCradle *detCradle) : TObject() {
    // Assumes spectrometer setup + Hades creation have been performed ....
    // This constructor sets all needed class members.
    // refID:     runtime ID
    // detCradle: pointer to Hades detector cradle

    pCradleHades = NULL;
    fSizesCells  = NULL;
    pField       = NULL;
    pMagnet      = NULL;
    setDetectorCradle(detCradle);
    init(refID);
    setPrint(kFALSE);
}

HKalInput::~HKalInput() {
}

//  -----------------------------------
//  Implementation of public methods
//  -----------------------------------

void HKalInput::calcSegPoints(const HMdcSeg *seg, Double_t& x1,  Double_t& y1, Double_t& z1, Double_t& x2, Double_t& y2, Double_t& z2) {
    // calculates 2 coordinate points from segment.

    Double_t teta = seg->getTheta();
    Double_t phi  = seg->getPhi();
    Double_t z    = seg->getZ();
    Double_t r    = seg->getR();
    Double_t pi2  = acos(-1.)/2.;

    Double_t X = r * cos(phi + pi2);
    Double_t Y = r * sin(phi + pi2);
    Double_t Z = z;

    x1 = X;
    y1 = Y;
    z1 = Z;
    x2 = X + cos(phi) * sin(teta);
    y2 = Y + sin(phi) * sin(teta);
    z2 = Z + cos(teta);
}


void HKalInput::clearCategories() {
    // Clear the categories (for macros in MakeCode style).

    if(kineCat)            kineCat           ->Clear();
    if(mdcTrkCandCat)      mdcTrkCandCat     ->Clear();
    if(mdcTrkCandCat)      mdcTrkCandCat     ->Clear();
    if(mdcSegCat)          mdcSegCat         ->Clear();
    if(mdcHitCat)          mdcHitCat         ->Clear();
    if(mdcTrkCandIdealCat) mdcTrkCandIdealCat->Clear();
    if(mdcSegIdealCat)     mdcSegIdealCat    ->Clear();
    if(mdcHitIdealCat)     mdcHitIdealCat    ->Clear();
}

void HKalInput::getCategories() {
    // Retries the pointer to needed input categories. Creates
    // iterators.

    geantMdcCat         = (HCategory*)gHades->getCurrentEvent()->getCategory(catMdcGeantRaw);
    kineCat             = (HCategory*)gHades->getCurrentEvent()->getCategory(catGeantKine);

    mdcTrkCandCat       = (HCategory*)gHades->getCurrentEvent()->getCategory(catMdcTrkCand);
    mdcSegCat           = (HCategory*)gHades->getCurrentEvent()->getCategory(catMdcSeg);
    mdcHitCat           = (HCategory*)gHades->getCurrentEvent()->getCategory(catMdcHit);

    mdcTrkCandIdealCat  = (HCategory*)gHades->getCurrentEvent()->getCategory(catMdcTrkCandIdeal);
    mdcSegIdealCat      = (HCategory*)gHades->getCurrentEvent()->getCategory(catMdcSegIdeal);
    mdcHitIdealCat      = (HCategory*)gHades->getCurrentEvent()->getCategory(catMdcHitIdeal);

    mdcCal1Cat          = (HCategory*)gHades->getCurrentEvent()->getCategory(catMdcCal1);

    if(kineCat) iterKine = (HIterator *) kineCat->MakeIterator("native");
    else        iterKine = NULL;

    if(mdcTrkCandCat) iterMdcTrkCand = (HIterator *) mdcTrkCandCat->MakeIterator("native");
    else              iterMdcTrkCand = NULL;

    if(mdcTrkCandIdealCat) iterMdcTrkCandIdeal = (HIterator *) mdcTrkCandIdealCat->MakeIterator("native");
    else                   iterMdcTrkCandIdeal = NULL;

}

Bool_t HKalInput::init(Int_t refID) {
    Bool_t bSuccess;
    if(!gHades) {
        Error("init()", "Hades object does not exist!");
        exit(1);
    }

    getCategories();                   // get Pointer to categories
    bSuccess  = initSizesCells(refID);             // init HMdcSizesCells
    if(!bSuccess) {
        Warning("init()", "Could not initialize SizesCells.");
    }
    bSuccess &= initBField(refID);                 // init Bfield parameters
    if(!bSuccess) {
        Warning("init()", "Could not initialize B-field.");
    }
    return bSuccess;
}


Bool_t HKalInput::initBField(Int_t refID) {
    // Initialize HMdcTrackGFieldPar and HMagnetPar for run refID.
    // If the container already exists it will not
    // be created a second time.

    HRuntimeDb* rtdb = gHades->getRuntimeDb();
    if(!pField) {
        pField =  (HMdcTrackGFieldPar*)(rtdb->getContainer("MdcTrackGFieldPar"));
    }
    if(!pMagnet) {
        pMagnet = (HMagnetPar*)( rtdb->getContainer("MagnetPar"));
    }
    return (pField && pMagnet);
}


Bool_t HKalInput::initSizesCells(Int_t refID) {
    // Initialize HMdcSizesCells for run refID.
    // If the container already exists it will not
    // be created a second time.

    fSizesCells = HMdcSizesCells::getExObject();  // check if is already there

    if(!fSizesCells) {
        fSizesCells = HMdcSizesCells::getObject();

        return fSizesCells->initContainer();
    } else {
        return kTRUE;
    }
}

HMdcSizesCells* HKalInput::getSizesCells() {
    if(!fSizesCells) {
        Error("getSizesCells()","HMdcSizesCells is NULL!");
        exit(1);
    }
    return fSizesCells;
}



HMdcTrkCand* HKalInput::nextMdcTrackCand(TObjArray &allhits, Int_t measDim, Double_t scaleErr) {
    // Fills 4 hits in sector coordinates system to allhits array.
    // Loops over HMdcTrkCand Category for the next for hits. The
    // iterator for the category has to be reset by hand at the begin
    // of each event. Return 0 if no next candidate has been found.
    // The errors of the hits will be taken from HMdcSeg and can be scaled
    // by scaleErr.

    if(fSizesCells == NULL){
        Error("nextMdcTrackCand()","HMdcSizeCells is NULL!");
        exit(1);
    }
    if(pCradleHades == NULL) {
        Error("nextMdcTrackCand()","DetectorCadle is NULL!");
        exit(1);
    }
    Int_t nrhits = 0;

    if(!iterMdcTrkCand) return NULL;

    HMdcTrkCand *mdcTrkCand;

    // This is the SEGment class of MDC. This data container
    // holds the information about a straight track let between
    // a pair of MDCs.

    Int_t   indexSeg    [2];   // index in catMdcSeg for inner/outer seg
    HMdcSeg *mdcSeg     [2];   // segment pointer for inner/outer seg
    Int_t   mdcSegIndHit[4];
    HMdcHit *hit        [4];   // hit pointer for both hits in segment

    Double_t coordHit   [3];
    Double_t coordErrHit[3];
    Double_t coordErrLay[3];

    // Collecting track candidates.
    for(Int_t ihit = 0; ihit < 4; ihit ++) mdcSegIndHit[ihit] = -1;

    allhits.Delete();
    while ((mdcTrkCand = (HMdcTrkCand*)iterMdcTrkCand->Next()) != 0) {
        // number of track candidates for outer segment
        Short_t nOuterSeg = mdcTrkCand->getNCandForSeg1();
        // no candidates in outer segments
        if(nOuterSeg == 0){
            continue;
        } else {
            if(nOuterSeg == -1)  continue;
            // loop through track candidates
            Short_t nextObjInd = 0;
            while(nextObjInd >= 0) {
                if(bPrint > 0) cout<<"get index sec "<<(Int_t)mdcTrkCand->getSec()<<endl;
                // Next HMdcTrkCand and segment getting:

                for(Int_t ihit = 0; ihit < 4; ihit ++) mdcSegIndHit[ihit] = -1;

                nextObjInd = mdcTrkCand->getNextCandInd();
                if(nextObjInd >= 0) {
                    mdcTrkCand = (HMdcTrkCand*)mdcTrkCandCat->getObject(nextObjInd);
                    for(Int_t io = 0; io < 2; io ++) { // loop inner/outer segment
                        indexSeg[io] =  mdcTrkCand->getSegInd(io);
                        mdcSeg  [io] = (HMdcSeg*) mdcSegCat->getObject(indexSeg[io]);
                        if(!mdcSeg  [io]){ cout<<"zero segment!"; continue;}
                        for(Int_t nhits = 0; nhits < 2; nhits ++) { // loop hits in segment
                            mdcSegIndHit[io * 2 + nhits] = mdcSeg[io]->getHitInd(nhits);
                            if(bPrint > 0) cout<<io * 2 + nhits<<" "<< mdcSegIndHit[io * 2 + nhits]<<endl;
                        }
                    }
                }


                //-----------------------------------------------------------------------
                // fill hit array for Kalman filter

                // If a Hit doesn't exist -1 is returned.
                if(mdcSegIndHit[0] != -1 && mdcSegIndHit[1] != -1 &&
                        mdcSegIndHit[2] != -1 && mdcSegIndHit[3] != -1 ) {
                    if(bPrint > 0) cout<<"filling sec "<<(Int_t)mdcSeg[0]->getSec()<<endl;
                    Int_t sec,mod;
                    //-------inner segment----------------------------
                    HMdcHit* locHit; // temp hit pointer
                    for(Int_t io = 0; io < 2; io ++) {              // loop inner/outer segment
                        for(Int_t nhits = 0; nhits < 2; nhits ++) { // loop hits in segment
                            locHit = hit[io * 2 + nhits]; //?
                            locHit = (HMdcHit*)mdcHitCat->getObject(mdcSegIndHit[io * 2 + nhits]);
                            if(!locHit){ cout<<"zero hit!"; continue;}
                            locHit->getSecMod(sec, mod);
                            const HKalMdcMeasLayer* hitdmeas = pCradleHades->getMdcLayerAt(mod, sec, 0);

                            // Get HMdcHit coordinates and errors in local Mdc coordinate system.
                            coordHit   [0] = locHit->getX();
                            coordHit   [1] = locHit->getY();
                            coordHit   [2] = 0;

                            (*fSizesCells)[sec][mod].transFrom(coordHit[0], coordHit[1], coordHit[2]);

                            coordErrLay[0] = locHit->getErrX() * scaleErr;
                            coordErrLay[1] = locHit->getErrY() * scaleErr;
                            coordErrLay[2] = 0.03;

                            (*fSizesCells)[sec][mod].transFrom(coordErrHit[0], coordErrHit[1], coordErrHit[2]);
                            (*fSizesCells)[sec][mod].transFrom(coordErrLay[0], coordErrLay[1], coordErrLay[2]);
                            for(Int_t i = 0; i < 3; i++) {
                                coordErrHit[i] = TMath::Abs(coordHit[i] - coordErrLay[i]);
                            }
                            allhits.Add(new HKalMdcHit(&coordHit[0], &coordErrHit[0], *hitdmeas, measDim));
                            nrhits++;

                        }
                    }
                    return mdcTrkCand;
                } // end filling
            } // end while outer cand
        } // going ncandfor outer > 0
    }

    return mdcTrkCand;
}

HMdcTrkCandIdeal* HKalInput::nextMdcTrackCandIdeal(TObjArray &allhits, Int_t measDim,
        Double_t errX, Double_t errY, Double_t scaleHit, Double_t scaleErr) {
    // Fills 4 hits in sector coordinates system to allhits array.
    // Loops over HMdcTrkCandIdeal Category for the next four hits. The
    // iterator for the category has to be reset by hand at the begin
    // of each event. Return 0 if no next candidate has been found.
    // The hit coordinates will be smeared by as Gaussian function
    // with mean values of 200/100 mu (x:y). The smearing can be scaled
    // by scaleHit. Pass scaleHit=0 to do no coordinates smearing.
    // The errors of the hits will be taken constant 200/100 mu (x:y)
    // and can be scaled by scaleErr.

    if(fSizesCells == NULL) {
        Error("nextMdcTrackCand()","HMdcSizeCells is NULL!");
        exit(1);
    }
    if(pCradleHades == NULL) {
        Error("nextMdcTrackCandIdeal()","DetectorCradle is NULL!");
        exit(1);
    }


    if(!iterMdcTrkCandIdeal) return NULL;

    HMdcTrkCandIdeal *mdcTrkCand;
    HMdcTrkCandIdeal *mdcTrkCandOut = 0;

    // This is the SEGment class of MDC. This data container
    // holds the information about a straight track let between
    // a pair of MDCs.

    Int_t       indexSeg[2];   // index in catMdcSeg for inner/outer seg
    HMdcSegIdeal *mdcSeg[2];   // segment pointer for inner/outer seg
    Int_t   mdcSegIndHit[4];
    HMdcHitIdeal *hit   [4];   // hit pointer for both hits in segment

    Double_t coordHit   [3];
    Double_t coordGeant [3];
    Double_t coordErrHit[3];
    Double_t coordErrLay[3];

    // Collecting track candidates.
    for(Int_t ihit = 0; ihit < 4; ihit ++) mdcSegIndHit[ihit] = -1;

    allhits.Delete();
    while ((mdcTrkCand = (HMdcTrkCandIdeal*)iterMdcTrkCandIdeal->Next()) != 0) {
        // number of track candidates for outer segment
        Short_t nOuterSeg = mdcTrkCand->getNCandForSeg1();
        // no candidates in outer segments
        if(nOuterSeg == 0){
            continue;
        } else {
            if(nOuterSeg == -1)  continue;
            // loop through track candidates
            Short_t nextObjInd = 0;
            while(nextObjInd >= 0) {
                if(bPrint > 0) cout<<"get index sec "<<(Int_t)mdcTrkCand->getSec()<<endl;
                // Next HMdcTrkCand and segment getting:

                for(Int_t ihit = 0; ihit < 4; ihit ++) mdcSegIndHit[ihit] = -1;

                nextObjInd = mdcTrkCand->getNextCandInd();
                if(nextObjInd >= 0) {
                    mdcTrkCand = (HMdcTrkCandIdeal*)mdcTrkCandIdealCat->getObject(nextObjInd);
                    mdcTrkCandOut = mdcTrkCand;
                    for(Int_t io = 0; io < 2; io ++) { // loop inner/outer segment
                        indexSeg[io] =  mdcTrkCand->getSegInd(io);
                        mdcSeg  [io] = (HMdcSegIdeal*) mdcSegIdealCat->getObject(indexSeg[io]);
                        if(!mdcSeg  [io]){ cout<<"zero segment!"; continue;}
                        for(Int_t nhits = 0; nhits < 2; nhits ++) { // loop hits in segment
                            mdcSegIndHit[io * 2 + nhits] = mdcSeg[io]->getHitInd(nhits);
                            if(bPrint > 0) cout<<io * 2 + nhits<<" "<< mdcSegIndHit[io * 2 + nhits]<<endl;
                        }
                    }
                }


                //-----------------------------------------------------------------------
                // fill hit array for  Kalman filter

                // If a Hit doesn't exist -1 is returned.
                if(mdcSegIndHit[0] != -1 && mdcSegIndHit[1] != -1 &&
                        mdcSegIndHit[2] != -1 && mdcSegIndHit[3] != -1 ) {
                    if(bPrint > 0) cout<<"filling sec "<<(Int_t)mdcSeg[0]->getSec()<<endl;
                    Int_t sec,mod;
                    //-------inner segment----------------------------
                    HMdcHit* locHit; // temp hit pointer
                    for(Int_t io = 0; io < 2; io ++) {              // loop inner/outer segment
                        for(Int_t nhits = 0; nhits < 2; nhits ++) { // loop hits in segment
                            locHit = hit[io * 2 + nhits]; //?
                            locHit = (HMdcHitIdeal*)mdcHitIdealCat->getObject(mdcSegIndHit[io * 2 + nhits]);
                            if(!locHit){ cout<<"zero hit!"; continue;}
                            locHit->getSecMod(sec, mod);
                            const HKalMdcMeasLayer* hitdmeas = pCradleHades->getMdcLayerAt(mod, sec, 0);
                            // Errors in Mdc coordinate system in mm.
                            // No errors from segment available
                            Double_t errX = 0.2;
                            Double_t errY = 0.1;
                            Double_t errZ = 0.003;
                            // Get HMdcHit coordinates and errors in local Mdc coordinate system.
                            coordHit   [0] = locHit->getX() + gRandom->Gaus() * errX * scaleHit;
                            coordHit   [1] = locHit->getY() + gRandom->Gaus() * errY * scaleHit;
                            coordHit   [2] = 0;
                            coordGeant [0] = locHit->getX();
                            coordGeant [1] = locHit->getY();
                            coordGeant [2] = 0.;

                            if(measDim > 2) {

                                // Create error hit point.
                                coordErrLay[0] = errX * scaleErr;
                                coordErrLay[1] = errY * scaleErr;
                                coordErrLay[2] = errZ;

                                // Transform hit and error hit to sector coordinates.
                                (*fSizesCells)[sec][mod].transFrom(coordHit[0], coordHit[1], coordHit[2]);
                                (*fSizesCells)[sec][mod].transFrom(coordErrLay[0], coordErrLay[1], coordErrLay[2]);
                                // The difference of hit and error hit points yield the error bars in
                                // sector coordinates.
                                for(Int_t i = 0; i < 3; i++) {
                                    coordErrHit[i] = TMath::Abs(coordHit[i] - coordErrLay[i]);
                                }
                            } else {
                                // Transform hit to sector coordinates. Don't transform errors.
                                (*fSizesCells)[sec][mod].transFrom(coordHit[0], coordHit[1], coordHit[2]);
                                coordErrHit[0] = errX * scaleErr;
                                coordErrHit[1] = errY * scaleErr;
                                coordErrHit[2] = errZ;
                            }
                            allhits.Add(new HKalMdcHit(&coordHit[0], &coordErrHit[0], *hitdmeas, measDim));
                        }
                    }
                    return mdcTrkCandOut;
                } // end filling
            } // end while outer cand
        } // going ncand for outer > 0
    }
    return mdcTrkCandOut;
}

HGeantKine* HKalInput::getGeantKine(const HMdcTrkCand *cand) {
    if(!cand)          { Error("getGeantKine()","HMdcTrkCand == NULL!"); return NULL; }
    if(!mdcSegCat)     { Error("getGeantKine()","mdcSegCat == NULL!"  ); return NULL; }
    HMdcSegSim* mdcSeg = (HMdcSegSim*) mdcSegCat->getObject(cand->getSegInd(0));
    Int_t ind = mdcSeg->getTrack(0);
    if(ind < 0)        { Error("getGeantKine()","ind < 0!"                 ); return NULL; }
    if(!kineCat)       { Error("getGeantKine()","kineCat == NULL!"         ); return NULL; }

    return (HGeantKine* )kineCat->getObject(ind-1);
}

HGeantKine* HKalInput::getGeantKine(const HMdcTrkCandIdeal *cand) {
    if(!cand)          {Error("getGeantKine()","HMdcTrkCandIdeal == NULL!"); return NULL;}
    if(!mdcSegIdealCat){Error("getGeantKine()","mdcSegIdealCat == NULL!"  ); return NULL;}
    HMdcSegIdeal* mdcSeg = (HMdcSegIdeal*) mdcSegIdealCat->getObject(cand->getSegInd(0));
    Int_t ind = mdcSeg->getTrack(0);
    if(ind < 0)        {Error("getGeantKine()","ind < 0!"                 ); return NULL;}
    if(!kineCat)       {Error("getGeantKine()","kineCat == NULL!"         ); return NULL;}

    return (HGeantKine* )kineCat->getObject(ind-1);
}

void HKalInput::getGeantMdcRawpoints(TObjArray &allhits, HGeantKine *kine, Bool_t midplane) {
    // Return TObjArray with pointer to HGeantMdc objects for the track.
    // if midplane == kTRUE  fill only points for layer == 6.
    // if midplane == KFALSE fill only points for layer < 6.

    if(!kine)         {Error("getGeantMdcRawpoints()","HGeantKine == NULL!"); return ;}

    allhits.Clear();

    HGeantMdc* geantMdc;
    kine->resetMdcIter();
    while( (geantMdc = (HGeantMdc*)kine->nextMdcHit()) != 0) {
        if( midplane && geantMdc->getLayer() != 6) continue;  // only hits for modules
        if(!midplane && geantMdc->getLayer() == 6) continue;  // only hits for layers
        allhits.Add(geantMdc);
    }
}

void HKalInput::getPosAndDir(TVector3 &pos, TVector3 &dir, HGeantMdc *geantMdc, Bool_t sectorCoord) {
    // Return pos and dir from HGeantMdc.
    // if sectorCoord == kTRUE the points will be transformed
    // to sector coordinate system.

    if(!geantMdc) { Error("getPosAndDir()","HGeantMdc == NULL!"); return ;}

    Double_t point   [3];
    Double_t pointDir[3];

    Float_t x,y,tof,p,th,ph;
    geantMdc->getHit(x,y,tof,p);

    point[0] = x;
    point[1] = y;
    point[2] = 0;

    geantMdc->getIncidence(th,ph);
    Double_t theta = th * TMath::DegToRad();    // 0 -360 deg
    Double_t phi   = ph * TMath::DegToRad();    // 0 -180 deg
    Double_t sinTh = sin(theta);
    pointDir[0]    = sinTh * cos(phi)         * 1000; // just for precision * 1000
    pointDir[1]    = sinTh * sin(phi)         * 1000;
    pointDir[2]    = sqrt(1. - sinTh * sinTh) * 1000;

    pointDir[0]   += point[0];
    pointDir[1]   += point[1];
    pointDir[2]   += point[2];

    if(sectorCoord){
        Int_t s = geantMdc->getSector();
        Int_t m = geantMdc->getModule();
        Int_t l = geantMdc->getLayer();

        if(l == 6) {
            ((*fSizesCells)[s][m]).transFrom(point   [0],point   [1],point   [2]); // transform to sector coordinate system
            ((*fSizesCells)[s][m]).transFrom(pointDir[0],pointDir[1],pointDir[2]); // transform to sector coordinate system
        } else {
            ((*fSizesCells)[s][m][l]).transFrom(point   [0],point   [1],point   [2]); // transform to sector coordinate system
            ((*fSizesCells)[s][m][l]).transFrom(pointDir[0],pointDir[1],pointDir[2]); // transform to sector coordinate system
        }
    }

    pos.SetX(point[0]);
    pos.SetY(point[1]);
    pos.SetZ(point[2]);

    dir.SetX(pointDir[0]);
    dir.SetY(pointDir[1]);
    dir.SetZ(pointDir[2]);

    dir -= pos ;
}

void HKalInput::getMdcSegs(const HMdcTrkCand *cand, HMdcSeg **segs) {
    if(!cand)      {Error("getGeantKine()","HMdcTrkCand == NULL!"); return ;}
    if(!mdcSegCat) {Error("getGeantKine()","mdcSegCat == NULL!"  ); return ;}

    segs[0]     = (HMdcSeg*) mdcSegCat->getObject(cand->getSegInd(0));
    if(cand->getSegInd(1)!=-1){
        segs[1] = (HMdcSeg*) mdcSegCat->getObject(cand->getSegInd(1));
    } else {
        segs[1] = NULL;
    }
}

void HKalInput::getMdcSegs(const HMdcTrkCandIdeal *cand, HMdcSegIdeal **segs) {
    if(!cand)          {Error("getGeantKine()","HMdcTrkCand == NULL!");      return ;}
    if(!mdcSegIdealCat){Error("getGeantKine()","mdcSegIdealCat == NULL!"  ); return ;}

    segs[0]     = (HMdcSegIdeal*) mdcSegIdealCat->getObject(cand->getSegInd(0));
    if(cand->getSegInd(1)!=-1){
        segs[1] = (HMdcSegIdeal*) mdcSegIdealCat->getObject(cand->getSegInd(1));
    } else {
        segs[1] = NULL;
    }
}


void HKalInput::getWire(HMdcSeg* const seg) {
    cout<<"call getWire"<<endl;

    HLocation loccal;
    loccal.set(4,0,0,0,0);

    HMdcCal2ParSim *cal2sim = (HMdcCal2ParSim*)gHades->getRuntimeDb()->getContainer("MdcCal2ParSim");
    HMdcCal2Par    *cal2    = (HMdcCal2Par   *)gHades->getRuntimeDb()->getContainer("MdcCal2Par");

    for(Int_t l = 0; l < 12; l ++) { // segment = 2 modules = 2x6 layers

        // input selection
        Int_t nCell = seg->getNCells(l); // number of cells per layer
        cout<<"number cells: "<<nCell<<endl;
        loccal[0]   = seg->getSec();
        Int_t ioseg = seg->getIOSeg(); // 0=inner,1=outer

        for(Int_t i = 0; i < nCell; i ++) {
            if(ioseg == 0)  {
                (l < 6) ? loccal[1] = 0 : loccal[1] = 1;
            } else if(ioseg == 1) {
                (l < 6) ? loccal[1] = 2 : loccal[1] = 3;
            }
            (l < 6) ? loccal[2] = l : loccal[2] = l - 6;

            // input selection
            loccal[3] = seg->getCell(l,i);  // get cell number

            Double_t x1, y1, z1, x2, y2, z2 = 0.;
            calcSegPoints(seg, x1, y1, z1, x2, y2, z2);
            Double_t alpha, mindist;
            (*fSizesCells)[loccal[0]][loccal[1]][loccal[2]].getImpact(x1, y1, z1, x2, y2, z2, loccal[3], alpha, mindist);

            HMdcCal1 *cal1 = (HMdcCal1*)mdcCal1Cat->getObject(loccal);
            if(cal1 != 0) {
                Float_t t1  = cal1->getTime1();
                Float_t t2  = cal1->getTime2();

                Double_t timeSimErr = cal2sim->calcTimeErr(loccal[0] ,loccal[1], alpha, mindist);
                Double_t timeSim    = cal2sim->calcTime(loccal[0] ,loccal[1], alpha, mindist);

                Double_t dist      = cal2->calcDistance(loccal[0], loccal[1], alpha, t1);
                Double_t distErr   = cal2->calcDistanceErr(loccal[0] ,loccal[1], alpha, t1);

                cout<<"------"<<endl;
                cout<<"measured time1: "<<t1<<endl;
                cout<<"measured time2: "<<t2<<endl;
                cout<<"angle: "<<alpha<<", min. distance: "<<mindist<<endl;
                cout<<"timeSim: "<<timeSim<<", timeSim Error: "<<timeSimErr<<endl;
                cout<<"calc. distance: "<<dist<<" with error "<<distErr<<endl;
            }
        }
    }
}