#include "hmdclookuptb.h"
#include "hspectrometer.h"
#include "hmdcdetector.h"
#include "hpario.h"
#include "hmdcgetcontainers.h"
#include <iostream.h>
#include "hmdcgeomstruct.h"
#include "hspecgeompar.h"
#include "hmdcgeompar.h"
#include "hmdcsizescells.h"
#include "hgeomvolume.h"
#include "hgeomcompositevolume.h"
#include "hmdclayergeompar.h"
#include "TObjArray.h"
#include "TH2.h"
#include "hmdcclussim.h"

//*-- AUTHOR : Pechenov Vladimir
//*-- Modified : 04/06/2002 by V.Pechenov
//*-- Modified : 09/05/2001 by V.Pechenov
//*-- Modified : 12/07/2000 by V.Pechenov
//*-- Modified : 23/05/2000 by V.Pechenov
//*-- Modified : 07/03/2000 by R. Holzmann
//*-- Modified : 02/12/99 by V.Pechenov
//*-- Modified : 26/10/99 by V.Pechenov
//*-- Modified : 20/05/99

////////////////////////////////////////////////////////////////
// HMdcLookUpTb
//
// Trackfinder for MDC1&2 and MDC3&4 if magnet off
//
////////////////////////////////////////////////////////////////

ClassImp(HMdcLookUpTbCell)
ClassImp(HMdcLookUpTbLayer)
ClassImp(HMdcLookUpTbMod)
ClassImp(HMdcLookUpTbSec)
ClassImp(HMdcLookUpTb)
ClassImp(HMdcCluster)

HMdcLookUpTbCell::HMdcLookUpTbCell() {
  nLines=0;
  xBinMin=xBinMax=0;
}

void HMdcLookUpTbCell::init(Int_t yBinMinT, Int_t yBinMaxT) {
  Int_t newNLines=yBinMaxT-yBinMinT+1;
  if(newNLines>nLines) {
    if(xBinMin) delete [] xBinMin;
    if(xBinMax) delete [] xBinMax;
    nLines=newNLines;
    xBinMin=new Short_t [nLines];
    xBinMax=new Short_t [nLines];
  }
  yBinMin=yBinMinT;
  yBinMax=yBinMaxT;
  line=0;
}

HMdcLookUpTbCell::~HMdcLookUpTbCell() {
  if(xBinMin) delete [] xBinMin;
  if(xBinMax) delete [] xBinMax;
}

Bool_t HMdcLookUpTbCell::addLine(Short_t nc1, Short_t nc2){
  if(line>=0 && line<nLines) {
    xBinMin[line]=nc1;
    xBinMax[line]=nc2;
    line++;
    return kTRUE;
  }
  else Error("addLine","line number > max. (%i > %i)",line,nLines);
  return kFALSE;
}

//----------Layer-----------------------------

HMdcLookUpTbLayer::HMdcLookUpTbLayer(Int_t sec, Int_t mod, Int_t layer) {
  // Geting of pointers to par. cont.
  HMdcGetContainers* parCont=HMdcGetContainers::getObject();
  if( !parCont ) return;
  HMdcGeomStruct* fMdcGeomStruct=parCont->getMdcGeomStruct();
  if( !fMdcGeomStruct ) return;
  Int_t nCells=((*fMdcGeomStruct)[sec][mod])[layer];
  array = new TObjArray(nCells);
  for(Int_t cell=0; cell<nCells; cell++) (*array)[cell]=new HMdcLookUpTbCell();
}

HMdcLookUpTbLayer::~HMdcLookUpTbLayer() {
  if(array) {
    array->Delete();
    delete array;
  }
}

Int_t HMdcLookUpTbLayer::getSize(void) {
  // return the size of the pointer array
  if(array) return array->GetEntries();
  else return -1;
}

//------------Module----------------------------
HMdcLookUpTbMod::HMdcLookUpTbMod(Int_t sec, Int_t mod) {
  // constructor creates an array of pointers of type HMdcLookUpTbMod
  array = new TObjArray(6);
  for (Int_t layer = 0; layer < 6; layer++)
    (*array)[layer] = new HMdcLookUpTbLayer(sec, mod, layer);
  nLayers=6;
}

HMdcLookUpTbMod::~HMdcLookUpTbMod() {
  // destructor
  if(array) {
    array->Delete();
    delete array;
  }
}

Int_t HMdcLookUpTbMod::getSize(void) {
  // returns the size of the pointer array
  if(array) return array->GetEntries();
  else return -1;
}

//----------HMdcCluster-------------------------------
HMdcCluster::HMdcCluster() {
  radToDeg=180./acos(Double_t(-1.));
}

void HMdcCluster::clear(void) {
  status=kTRUE;
  clusMerg=0;
  nMergedClus=1;
  lCells1.clear();
  lCells2.clear();
  nBins=0;
  meanX=meanY=meanXX=meanYY=meanYX=0.;
  sumWt=meanXWt=meanYWt=meanXXWt=meanYYWt=0.;
  clusMod1=clusMod2=0;
}

void HMdcCluster::addClus(HMdcCluster& clst2, Bool_t flag) {
  lCells1.add(&(clst2.lCells1));
  if(flag) lCells2.add(&(clst2.lCells2));
  clst2.status=kFALSE;
  clst2.clusMerg=this;
  sumWt += clst2.sumWt;
  nBins += clst2.nBins;
  meanX += clst2.meanX;
  meanY += clst2.meanY;
  meanXX += clst2.meanXX;
  meanYY += clst2.meanYY;
  meanYX += clst2.meanYX;
  meanXWt += clst2.meanXWt;
  meanYWt += clst2.meanYWt;
  meanXXWt += clst2.meanXXWt;
  meanYYWt += clst2.meanYYWt;
  nMergedClus += clst2.nMergedClus;
  calcXY();
}

void HMdcCluster::calcClusParam(Double_t xSt2, Double_t ySt2) {
  errX=sqrt(meanXXWt/sumWt-x*x + xSt2);
  errY=sqrt(meanYYWt/sumWt-y*y + ySt2);
  // Cluster shape:-------------------------------------------------
  //  eXX,eYY,eXY=eYX - covariance matrix componets
  //
  //    | eXX-sigma   aXY       |           | E1 |
  // A= |                       |         B=|    |
  //    | eYX         eYY-sigma |           | E2 |
  //
  // A^2=0;  ==> sigma1 & sigma2  - (sigma1 > sigma2)
  //
  //    | 0 |        | E1 |
  // AB=|   |;   ==> |    | - direction of main axis
  //    | 0 |        | E2 |
  if(nBins<2) {
    sigma1=sigma2=0.;
    alpha=0.;
  } else {
    Double_t nBn=nBins;
    Double_t eXX = meanXX*nBn- meanX*meanX;
    Double_t eYX = meanYX*nBn- meanY*meanX;
    Double_t eYY = meanYY*nBn- meanY*meanY;
    Double_t norm = (eXX+eYY)/2.;
    Double_t eYX2 = eYX*eYX;
    Double_t c=eXX*eYY-eYX2;
    Double_t sigma1D=norm+sqrt(norm*norm-c);
    Double_t sigma2D=norm-sqrt(norm*norm-c);
    if(sigma2<0.) sigma2=0.;
    Double_t e1=sigma1D-eXX;
    Double_t e2=sigma1D-eYY;
    sigma1=sqrt(sigma1D)/nBn;
    if(sigma2D<10e-5) sigma2=0.;
    else sigma2=sqrt(sigma2D)/nBn;
    alpha=atan2(sqrt(eYX2+e1*e1),sqrt(eYX2+e2*e2))*radToDeg;
    if(eYX<0.) alpha=180.-alpha;
  }
}

void HMdcCluster::fillClus(HMdcClus* fClus, Int_t nLst) {
  fClus->setSumWt(sumWt);
  fClus->setNBins(nBins);
  fClus->setXY(x,errX,y,errY);
  fClus->setNMergClust(nMergedClus);
  if(clusMod1==0) {
    Int_t nDrTm=(nLst==0) ? lCells1.getNDrTimes(0,5):lCells2.getNDrTimes(0,5);
    if(nDrTm>0) {
      fClus->setClusSizeM1(nBins);
      fClus->setNDrTimesM1(nDrTm);
      fClus->setNMergClustM1(nMergedClus);
      fClus->setShapeM1(sigma1,sigma2,alpha);
    } else fClus->clearMod1Par();
  } else {
    fClus->setClusSizeM1(clusMod1->nBins);
    fClus->setNDrTimesM1(clusMod1->lCells1.getNDrTimes(0,5));
    fClus->setNMergClustM1(clusMod1->nMergedClus);
    fClus->setShapeM1(clusMod1->sigma1,clusMod1->sigma2,clusMod1->alpha);
  }
  if(clusMod2==0) {
    Int_t nDrTm=(nLst==0) ? lCells1.getNDrTimes(6,11):lCells2.getNDrTimes(6,11);
    if(nDrTm>0) {
      fClus->setClusSizeM2(nBins);
      fClus->setNDrTimesM2(nDrTm);
      fClus->setNMergClustM2(nMergedClus);
      fClus->setShapeM2(sigma1,sigma2,alpha);
    } else fClus->clearMod2Par();
  } else {
    fClus->setClusSizeM2(clusMod2->nBins);
    fClus->setNDrTimesM2(clusMod2->lCells1.getNDrTimes(6,11));
    fClus->setNMergClustM2(clusMod2->nMergedClus);
    fClus->setShapeM2(clusMod2->sigma1,clusMod2->sigma2,clusMod2->alpha);
  }
}

//----------Sector------------------------------------

HMdcClFnStack* HMdcLookUpTbSec::stack=0;
Int_t HMdcLookUpTbSec::hPlModsSize=0;
UChar_t* HMdcLookUpTbSec::hPlMod[4]={0,0,0,0};
Int_t HMdcLookUpTbSec::sizeBArSt=0;
UChar_t* HMdcLookUpTbSec::plotBArSc=0;
UChar_t* HMdcLookUpTbSec::plotBArM[4]={0,0,0,0};
Short_t* HMdcLookUpTbSec::clusIndM1=0;
Int_t HMdcLookUpTbSec::clIndArrSzM1=0;
Short_t* HMdcLookUpTbSec::clusIndM2=0;
Int_t HMdcLookUpTbSec::clIndArrSzM2=0;
Int_t HMdcLookUpTbSec::clusArrSize=0;
HMdcCluster* HMdcLookUpTbSec::clusArr=0;
HMdcCluster* HMdcLookUpTbSec::clusArrM1=0;
HMdcCluster* HMdcLookUpTbSec::clusArrM2=0;

HMdcLookUpTbSec::HMdcLookUpTbSec(Int_t sec, Int_t nSegs,
                                 Int_t inBinX, Int_t inBinY) {
  // constructor creates an array of pointers of type HMdcLookUpTbMod
  sector=sec;
  nSegments=nSegs;
  array = new TObjArray(4);
  HMdcDetector* fMdcDet= HMdcGetContainers::getObject()->getMdcDetector();
  for (Int_t mod = 0; mod < nSegs*2; mod++)
    if(fMdcDet->getModule(sec,mod)) (*array)[mod]=new HMdcLookUpTbMod(sec,mod);
  hist=0;
  nBinX=(inBinX%2 == 0) ? inBinX:inBinX+1;
  nBinY=(inBinY%2 == 0) ? inBinY:inBinY+1;
  xBinsPos=new Double_t [nBinX];
  yBinsPos=new Double_t [nBinY];
  size=nBinX*nBinY;
  size=(size/32 + ((size%32 > 0) ? 1:0))*32;
  
  if(size>hPlModsSize) {
    if(hPlMod[0]) delete [] hPlMod[0];
    if(hPlMod[1]) delete [] hPlMod[1];
    hPlMod[0]=new UChar_t [size];
    hPlMod[1]=new UChar_t [size];
    if(nSegs==2 || hPlMod[2] || hPlMod[3]) {
      if(hPlMod[2]) delete [] hPlMod[2];
      if(hPlMod[3]) delete [] hPlMod[3];
      hPlMod[2]=new UChar_t [size];
      hPlMod[3]=new UChar_t [size];
    }
    hPlModsSize=size;
    clearPrArrs();
  }

  if(plotBArSc==0 || sizeBArSt<size/8) {
    if(plotBArSc) delete [] plotBArSc;
    for(Int_t mod=0;mod<4;mod++) if(plotBArM[mod]) delete [] plotBArM[mod];
    sizeBArSt=size/8;
    plotBArSc=new UChar_t [sizeBArSt];
    memset(plotBArSc,0,sizeBArSt);
    for(Int_t mod=0;mod<4;mod++) {
      plotBArM[mod]=new UChar_t [sizeBArSt];
      memset(plotBArM[mod],0,sizeBArSt);
    }
  }
  sizeBAr=size/8;
  typeClFinder=0;
  setLUpTb=kFALSE;
  neighbBins[0]=-1;
  neighbBins[1]=+1;
  neighbBins[2]=-nBinX;
  neighbBins[3]=+nBinX;
  neighbBins[4]=-1-nBinX;
  neighbBins[5]=-1+nBinX;
  neighbBins[6]=1-nBinX;
  neighbBins[7]=1+nBinX;
  if(stack==0) stack=new HMdcClFnStack(1000);
  isGeant = HMdcGetContainers::isGeant();
  locClus.set(3,sec,0,0);
  if(clusArrSize==0) {
    clusArrSize=500;    // Clusters array size!
    clusArr=new HMdcCluster [clusArrSize];
    clusArrM1=new HMdcCluster [clusArrSize]; 
    clusArrM2=new HMdcCluster [clusArrSize];
  }
  if(size>clIndArrSzM1) {
    if(clusIndM1) delete [] clusIndM1;
    clusIndM1=new Short_t [size];
    clIndArrSzM1=size;
  }
  if(size>clIndArrSzM2) {
    if(clusIndM2) delete [] clusIndM2;
    clusIndM2=new Short_t [size];
    clIndArrSzM2=size;
  }
  for(Int_t mod=0;mod<4;mod++) {
    xMin[mod]=new Int_t [nBinY];
    xMax[mod]=new Int_t [nBinY];
    for(Int_t y=0;y<nBinY;y++) {
      xMin[mod][y]=size;
      xMax[mod][y]=-1;;
    }
  }
  xMaxCl=new Short_t [nBinY];
  xMinCl=new Short_t [nBinY];
}

HMdcLookUpTbSec::~HMdcLookUpTbSec() {
  // destructor
  if(array) {
    array->Delete();
    delete array;
  }
  for(Int_t mod=0; mod<4; mod++) {
    if(hPlMod[mod]) delete [] hPlMod[mod];
    hPlMod[mod]=0;
  }
  hPlModsSize=0;
  if(stack) delete stack;
  stack=0;
  if(hist && hist->IsOnHeap()) delete hist;
  hist=0;
  if(plotBArSc) {
    delete [] plotBArSc;
    plotBArSc=0;
    sizeBArSt=0;
  }
  for(Int_t mod=0;mod<4;mod++) {
    if(plotBArM[mod]) {
      delete [] plotBArM[mod];
      plotBArM[mod]=0;
    }
  }
  delete [] xBinsPos;
  delete [] yBinsPos;
  for(Int_t mod=0;mod<4;mod++) {
    delete [] xMin[mod];
    delete [] xMax[mod];
  }
    
  if(clusArr)   delete [] clusArr;
  if(clusArrM1) delete [] clusArrM1; 
  if(clusArrM2) delete [] clusArrM2;
  clusArrSize=0;
  clusArr=clusArrM1=clusArrM2=0;

  if(clusIndM1) {
    delete [] clusIndM1;
    clusIndM1=0;
    clIndArrSzM1=0;
  }
  if(clusIndM2) {
    delete [] clusIndM2;
    clusIndM2=0;
    clIndArrSzM2=0;
  }
  delete [] xMaxCl;
  delete [] xMinCl;
}

void HMdcLookUpTbSec::clearPrArrs(void) {
  memset(hPlMod[0],0,size);
  memset(hPlMod[1],0,size);
  if(hPlMod[2]) memset(hPlMod[2],0,size);
  if(hPlMod[3]) memset(hPlMod[3],0,size);
}

void HMdcLookUpTbSec::clearPrMod(Int_t mod) {
  Int_t* hPlModM4=(Int_t*)(hPlMod[mod]);
  if(!hPlModM4) return;
  Int_t* xMaxM=xMax[mod];
  Int_t* xMinM=xMin[mod];
  Int_t* xMinMe=xMinM+nBinY;
  for(; xMinM<xMinMe; xMinM++) {
    if(*xMinM<=*xMaxM) {
      Int_t* eb=hPlModM4+((*xMaxM)>>2);
      for(Int_t* bt=hPlModM4+((*xMinM)>>2); bt<=eb; bt++) *bt=0;
      *xMinM=size;
      *xMaxM=-1;
    }
    xMaxM++;
  }  
}

Int_t HMdcLookUpTbSec::getSize(void) {
  // return the size of the pointer array
  if(array) return array->GetEntries();
  else return -1;
}

//May be numhits[i]=0 ... in another sub. ?????????????????
void HMdcLookUpTbSec::clearwk(void) {
  for(Int_t i=0; i<24; i++) nHits[i]=0;
  nHitsTot=0;
  setLUpTb=kFALSE;
  counter[0]=counter[1]=0;
}


void HMdcLookUpTbSec::setCell(Short_t mod, Short_t layer,
                              Short_t cell, Short_t nTimes){
  HMdcLookUpTbCell& fLUpTbLayer=(*this)[mod][layer][cell];
  if( &fLUpTbLayer && fLUpTbLayer.line>0 ) {
    Int_t ind1=mod*6+layer;
    Int_t ind2=nHits[ind1];
    if (ind2<250) {
      hits[ind1][ind2]=cell;
      hitsNTime[ind1][ind2]=nTimes;
      hitsDel[ind1][ind2]=0;
      nHits[ind1]++;
      nHitsTot++;
    }
    else Error("setCell","Too many hits in layer.");
  }
}

void HMdcLookUpTbSec::deleteHit(Short_t mod, Short_t layer,
                                Short_t cell, Short_t nTime){

 if(mod<2) {
   Short_t ind1=mod*6+layer;
   Short_t ind2=-1;
   for(Int_t i=0; i<nHits[ind1]; i++) {
     if(hits[ind1][i]==cell) { ind2=i; break;}
   }
   if (ind2>=0) {
     //hitsNTime[ind1][ind2]=1 - time1; =2 - tim2; =3 - time1&tim2
     if(hitsNTime[ind1][ind2]==1 || hitsNTime[ind1][ind2]==2 || nTime==0 )
       hitsDel[ind1][ind2]=-1;
     else if(hitsDel[ind1][ind2]==0) hitsDel[ind1][ind2]=3-nTime;
     else hitsDel[ind1][ind2]=-1;
   }
   else Error("deleteHit","Hit %iM%iL%iC is absent.\n",mod,layer,cell);
 }
/*
 HMdcLookUpTbLayer& fLUpTbLayer=(*this)[mod][layer];
 HMdcLookUpTbCell& fCell=fLUpTbLayer[cell];
 UChar_t nBin=(UChar_t) nCell;
 const UChar_t one=1;
 for(Int_t l=0; l<fCell.line; l++) {
   for(Int_t nb=fCell.xBinMin[l]; nb<=fCell.xBinMax[l]; nb++) {
     UChar_t bit=(one<<(nBin-fLUpTbLayer.nfcell[nb]));
     if(fLUpTbLayer.ncellbin[nb]&bit != 0) {
       fLUpTbLayer.ncellbin[nb]=fLUpTbLayer.ncellbin[nb]&(~bit);
       if(fLUpTbLayer.ncellbin[nb] == 0) hitPlS1[nb]--;
     }
   }
 }
*/
}

void HMdcLookUpTbSec::calcMaxAmp(void) {
  for(Int_t mod=0;mod<4;mod++) {
    maxAmp[mod]=0;
    for(Int_t layer=0;layer<6;layer++) if(nHits[mod*6+layer]>0) maxAmp[mod]++;
  }
}

void HMdcLookUpTbSec::makeSPlot(void) {
  maxBinBAr4Sc=0;
  minBinBAr4Sc=size;
  UChar_t addm=1;  // !!!
  Int_t nmods=-1;
  if(minAmp[3]>0) nmods=3;
  else if(minAmp[2]>0) nmods=2;
  else if(minAmp[1]>0) nmods=1;
  else if(minAmp[0]>0) nmods=0;
  Char_t minAm0=minAmp[0];
  Char_t minAm1=minAmp[1];
  Char_t minAm2=minAmp[2];
  if(nmods<0) return;
  for(Int_t mod=0; mod<=nmods; mod++) {
    cFMod=&((*this)[mod]);
    if( !cFMod ) continue;
    if(minAmp[mod]==0) continue;
    cHPlModM=hPlMod[mod];
    if(!cHPlModM) continue;
    cXMinM=xMin[mod];
    cXMaxM=xMax[mod];
    
    Int_t iLay=mod*6;
    Int_t nLay=(mod<nmods) ? 6:minAmp[mod]-1;
    // ??????????????????? eto poka staraya versiya
    for(Int_t layer=0; layer<nLay; layer++) {
      Int_t indLay=iLay+layer;
      makeLayProjV1(layer,indLay);
    }

    Char_t minAm=minAmp[mod];
    for(Int_t layer=nLay; layer<6; layer++) {
      Int_t i=iLay+layer;
      if(nHits[i]==0) continue;
      UChar_t add=addm<<layer;
      HMdcLookUpTbLayer& fLayer=(*cFMod)[layer];
      for(Int_t n=0; n<nHits[i]; n++) {
        if(hitsDel[i][n]<0) continue;
        HMdcLookUpTbCell& fCell=fLayer[hits[i][n]];
        for(Int_t ln=0; ln<fCell.line; ln++) {
          Int_t y=fCell.yBinMin+ln;
          Int_t shift=y * nBinX;
          Int_t nbL=fCell.xBinMax[ln]+shift;
          Int_t nbF=fCell.xBinMin[ln]+shift;
          if(nbF<cXMinM[y]) cXMinM[y]=nbF;
          if(nbL>cXMaxM[y]) cXMaxM[y]=nbL;
          UChar_t* bt=cHPlModM+nbF;
          for(Int_t nb=nbF; nb<=nbL; nb++) {
            *bt |= add;
            Char_t wt=HMdcBArray::getNSet(bt);
            bt++;
            if( wt<minAm ) continue;
            if(mod>0) {
              if(HMdcBArray::getNSet(hPlMod[0][nb])<minAm0) continue;
              if(mod>1) {
                if(HMdcBArray::getNSet(hPlMod[1][nb])<minAm1) continue;
                if(mod>2 && HMdcBArray::getNSet(hPlMod[2][nb])<minAm2) continue;
              }
            }
            HMdcBArray::set(plotBArSc,nb);
            if(nb<minBinBAr4Sc) minBinBAr4Sc=nb;
            if(nb>maxBinBAr4Sc) maxBinBAr4Sc=nb;
          }
        }
      }
    }
  }
  maxBinBAr4Sc/=32;
  minBinBAr4Sc/=32;
}

void HMdcLookUpTbSec::makeModPlot(Int_t mod,Int_t minAm) {
  if(minAm<2) return;
  cFMod=&((*this)[mod]);
  if( !cFMod ) return;
  cHPlModM=hPlMod[mod];
  cXMinM=xMin[mod];
  cXMaxM=xMax[mod];

  Int_t iLay=mod*6;
  Int_t indMeth=minAm-2;  
  Int_t lay[6]={2,3,1,4,0,5};   // order of layers at the pr.plot filling
  
  Int_t meth1[5]={1,2,3,2,1};   // if layer<meth1[minAm-2]      - makeLayProjV1
  Int_t meth2[5]={5,4,0,0,0};   // else if layer<meth2[minAm-2] - firs {}
  Int_t meth3[5]={5,4,3,4,5};   // else if layer<meth3[minAm-2] - makeLayProjV2
                                // else                         - second {}
  
//  Int_t mth1=meth1[indMeth];
//  Int_t mth2=meth2[indMeth];
//  Int_t mth3=meth3[indMeth];
  
//  Int_t* nHt=nHits+iLay;
  
//   if(minAm<6) {
//     if(minAm==5 && maxAmp[mod]==5/*???*/) {
//       // one layer can have 0 wiers only
//       if(nHt[lay[0]]>0) mth1=1;
//       else if(nHt[lay[5]]>0) mth3=5;
//     } else if(minAm==4) {
//       if(maxAmp[mod]==5) {
//         // one layer can have 0 wiers only
//         if(nHt[lay[3]]==0 || nHt[lay[4]]==0 || nHt[lay[5]]==0) meth1=2;
//         else mth3=4;
//       } else if(maxAmp[mod]==4) {
//       }
//     }
//     
//   }
    
  UChar_t* cPlotBAr=plotBArM[mod];
  Int_t& minBin=minBinBAr4M[mod];
  Int_t& maxBin=maxBinBAr4M[mod];
  maxBin=0;
  minBin=size;
  for(Int_t indLine=0; indLine<6; indLine++) {
    Int_t layer=lay[indLine];
    Int_t indLay=iLay+layer;
    if(nHits[indLay]==0) continue;
    if(indLine<meth1[indMeth]) makeLayProjV1(layer,indLay);
    else if(indLine<meth2[indMeth]) {
      // filling, min, max and bits seting
      UChar_t add=1<<layer;
      HMdcLookUpTbLayer& fLayer=(*cFMod)[layer];
      for(Int_t n=0; n<nHits[indLay]; n++) {
        if(hitsDel[indLay][n]<0) continue;
        HMdcLookUpTbCell& fCell=fLayer[hits[indLay][n]];
        for(Int_t ln=0; ln<fCell.line; ln++) {
          Int_t y=fCell.yBinMin+ln;
          if(cXMaxM[y]<0) continue;
          Int_t shift=y * nBinX;
          Int_t nbL=fCell.xBinMax[ln]+shift;
          Int_t nbF=fCell.xBinMin[ln]+shift;
          if(nbF<cXMinM[y]) cXMinM[y]=nbF;
          if(nbL>cXMaxM[y]) cXMaxM[y]=nbL;
          UChar_t* bt=cHPlModM+nbF; 
          for(Int_t nb=nbF; nb<=nbL; nb++) {
            *bt |= add;
            Char_t wt=HMdcBArray::getNSet(bt);
            bt++;
            if(wt<minAm) continue;
            HMdcBArray::set(cPlotBAr,nb);
            if(nb<minBin) minBin=nb;
            if(nb>maxBin) maxBin=nb;
          }
        }
      }      
    } else if(indLine<meth3[indMeth]) makeLayProjV2(layer,indLay);
    else {
      UChar_t add=1<<layer;
      HMdcLookUpTbLayer& fLayer=(*cFMod)[layer];
      for(Int_t n=0; n<nHits[indLay]; n++) {
        if(hitsDel[indLay][n]<0) continue;
        HMdcLookUpTbCell& fCell=fLayer[hits[indLay][n]];
        for(Int_t ln=0; ln<fCell.line; ln++) {
          Int_t y=fCell.yBinMin+ln;
          if(cXMaxM[y]<0) continue;
          Int_t shift=y * nBinX;
          Int_t nbL=fCell.xBinMax[ln]+shift;
          Int_t nbF=fCell.xBinMin[ln]+shift;
          if(nbL>cXMaxM[y]) nbL=cXMaxM[y];
          if(nbF<cXMinM[y]) nbF=cXMinM[y];
          UChar_t* bt=cHPlModM+nbF; 
          for(Int_t nb=nbF; nb<=nbL; nb++) {
            *bt |= add;
            Char_t wt=HMdcBArray::getNSet(bt);
            bt++;
            if(wt<minAm) continue;
            HMdcBArray::set(cPlotBAr,nb);
            if(nb<minBin) minBin=nb;
            if(nb>maxBin) maxBin=nb;
          }
        }
      }
    }
  }
  maxBin/=32;
  minBin/=32;
}

void HMdcLookUpTbSec::makeLayProjV1(Int_t layer, Int_t indLay) {
  // plot filling and filled region determination
  UChar_t add=1<<layer;
  HMdcLookUpTbLayer& fLayer=(*cFMod)[layer];
  for(Int_t n=0; n<nHits[indLay]; n++) {
    if(hitsDel[indLay][n]<0) continue;
    HMdcLookUpTbCell& fCell=fLayer[hits[indLay][n]];
    for(Int_t ln=0; ln<fCell.line; ln++) {
      Int_t y=fCell.yBinMin+ln;
      Int_t shift=y * nBinX;
      Int_t nbL=fCell.xBinMax[ln]+shift;
      Int_t nbF=fCell.xBinMin[ln]+shift;
      if(nbF<cXMinM[y]) cXMinM[y]=nbF;
      if(nbL>cXMaxM[y]) cXMaxM[y]=nbL;
      UChar_t* hPlModE=cHPlModM+nbL;
      for(UChar_t* bt=cHPlModM+nbF; bt<=hPlModE; bt++) *bt |= add;
    }
  }
}

void HMdcLookUpTbSec::makeLayProjV2(Int_t layer, Int_t indLay) {
  // plot filling in filled regions only
  UChar_t add=1<<layer;
  HMdcLookUpTbLayer& fLayer=(*cFMod)[layer];
  for(Int_t n=0; n<nHits[indLay]; n++) {
    if(hitsDel[indLay][n]<0) continue;
    HMdcLookUpTbCell& fCell=fLayer[hits[indLay][n]];
    for(Int_t ln=0; ln<fCell.line; ln++) {
      Int_t y=fCell.yBinMin+ln;
      if(cXMaxM[y]<0) continue;
      Int_t shift=y * nBinX;
      Int_t nbL=fCell.xBinMax[ln]+shift;
      Int_t nbF=fCell.xBinMin[ln]+shift;
      if(nbL>cXMaxM[y]) nbL=cXMaxM[y];
      if(nbF<cXMinM[y]) nbF=cXMinM[y];
      UChar_t* hPlModE=cHPlModM+nbL;
      for(UChar_t* bt=cHPlModM+nbF; bt<=hPlModE; bt++) *bt |= add;
    }
  }
}

Int_t HMdcLookUpTbSec::findClusters(Int_t *imax){
  nClusters=0;
  if(nHitsTot==0) return 0;
  calcMaxAmp();
  //isCoilOff - !!!??? ---------------------------------------
  nMods=0;
  for(Int_t mod=0; mod<4; mod++) {
    minAmp[mod]=(mod<nSegments*2) ? imax[mod]:0;
    if(minAmp[mod]>0) nMods++;
  }
  if(typeClFinder==1) {
    // chamber clust. finding
    for(Int_t m=0; m<4; m++) if(minAmp[m]>0) findClusInMod(m,minAmp[m]);
  } else {
    // combyned clust. finding
    nModSeg[0]=nModSeg[1]=0;
    for(Int_t m=0;m<4;m++) if(minAmp[m]>0) nModSeg[m/2] |= m%2+1;
    if(nModSeg[0]>0 && nModSeg[1]>0)
        findClusInSec(minAmp[0]+minAmp[1]+minAmp[2]+minAmp[3]);
    else {
      if(nModSeg[0]==3) findClusInSeg1();
        // findClusInSeg(0,minAmp[0]+minAmp[1]);
      else if(nModSeg[0]>0) findClusInMod(nModSeg[0]-1,minAmp[nModSeg[0]-1]);
      if(nModSeg[1]==3) findClusInSeg(1,minAmp[2]+minAmp[3]);
      else if(nModSeg[1]>0) findClusInMod(nModSeg[1]+1,minAmp[nModSeg[1]+1]);
    }
  }
  return nClusters;
}

void HMdcLookUpTbSec::getClusterSlot(Int_t seg, HMdcList12GroupCells& list) {
  locClus[1]=seg;
  locClus[2]=counter[seg];
  fClus = (HMdcClus*)fClusCat->getSlot(locClus);
  if(!fClus) {
    Warning("getClusterSlot","Sec.%i Segment %i No slot HMdcClus available",
        sector+1,seg+1);
    return;
  }
  if(isGeant) fClus=(HMdcClus*)(new(fClus) HMdcClusSim(list));
  else fClus=new(fClus) HMdcClus(list);
  fClus->setAddress(locClus);
  counter[seg]++;
}

void HMdcLookUpTbSec::fillWiresList(Int_t mod, HMdcList12GroupCells& list) {
  HMdcLookUpTbMod& fMod=(*this)[mod];
  if( !(&fMod) ) return;
  Int_t layAdd=(mod%2)*6;
  for(Int_t layer=0; layer<6; layer++) {
    Int_t iLayer=mod*6+layer;
    if(nHits[iLayer]<=0) continue;
    Int_t clLay=layAdd+layer;
    HMdcLookUpTbLayer& fLayer=fMod[layer];
    for(Int_t i=0; i<nHits[iLayer]; i++) {
      if(hitsDel[iLayer][i]<0) continue;
      Int_t cell=hits[iLayer][i];
      HMdcLookUpTbCell& fCell=fLayer[cell];
      Short_t yBinMax=fCell.yBinMax;
      if( nLMinCl > yBinMax ) continue;
      Short_t yBinMin=fCell.yBinMin;
      if( nLMaxCl < yBinMin ) continue;
//        if( nLMaxCl < yBinMin ) break;            Hits not sorted !
      Short_t lMax=(yBinMax < nLMaxCl) ? yBinMax : nLMaxCl;
      Short_t lMin=(yBinMin > nLMinCl) ? yBinMin : nLMinCl;
      for (Short_t nl=lMin; nl<=lMax; nl++) {
        if( fCell.xBinMin[nl-yBinMin] > xMaxCl[nl] ||
            fCell.xBinMax[nl-yBinMin] < xMinCl[nl] ) continue;
        list.setTime(clLay,cell,hitsNTime[iLayer][i]);
	break;
      }
    }
  }
}

void HMdcLookUpTbSec::findClusInSeg(Int_t seg, Short_t minAm){
  if(maxAmp[seg*2]<minAmp[seg*2] || maxAmp[seg*2+1]<minAmp[seg*2+1]) return;
  if(minAm<3) return;
  nClsArr=0;
  makeSPlot();
  UChar_t *hPlModF=(seg==0) ? hPlMod[0]:hPlMod[2];
  UChar_t *hPlModS=(seg==0) ? hPlMod[1]:hPlMod[3];
  Char_t maxClS=minAm-1;
  Int_t nmo=(seg+1)*2;
  cPlotBAr=plotBArSc;
  Int_t* plotBAr4b=(Int_t*) cPlotBAr;
  //---Cluster finding-------------------------------------
  for(Int_t n4=maxBinBAr4Sc; n4>=minBinBAr4Sc; n4--) {
    if(plotBAr4b[n4]==0) continue;
    UChar_t *b1=cPlotBAr+n4*4;
    UChar_t *b2=b1+3;
    Int_t nBin4=n4*32;
    Int_t nLBin=33;
    while((nLBin=HMdcBArray::prevAndUnset(b1,b2,nLBin))>=0) {
      //---New cluster---------------------------------------
      clus=&clusArr[nClsArr];
      Int_t nBin=nBin4+nLBin;   
      initCluster(nBin);
      while ((nBin=stack->pop()) >= 0) addBinInCluster(nBin,
            HMdcBArray::getNSet(hPlModF[nBin])+
            HMdcBArray::getNSet(hPlModS[nBin])-maxClS);
      //-Filling of cluster---------------------------------------------------
      clus->calcXY();
      for(Int_t mod=seg*2; mod<nmo; mod++) fillWiresList(mod,clus->lCells1);
      nClsArr++;
      if(nClsArr >= clusArrSize) break;
    }
    if(nClsArr >= clusArrSize) {
      Warning("findClusInSeg"," Num. of clusters in sector %i > max\n",sector);
      memset(cPlotBAr,0,sizeBAr);
      break;
    }
  }

  if(nClsArr>0) {
    if(nClsArr>1) mergeClusInSeg(seg,nClsArr,clusArr);
    fillClusCat(-2,seg);
  }
  for(Int_t mod=seg*2;mod<(seg+1)*2;mod++) if(minAmp[mod]>0) clearPrMod(mod);
}

void HMdcLookUpTbSec::findClusInSeg1(void){
  if(maxAmp[0]<minAmp[0] || maxAmp[1]<minAmp[1]) return;
  if(minAmp[0]+minAmp[1]<3) return;
  makeModPlot(0,minAmp[0]);
  makeModPlot(1,minAmp[1]);

  minBinBAr4Sc=(minBinBAr4M[0]>=minBinBAr4M[1])?minBinBAr4M[0]:minBinBAr4M[1];
  maxBinBAr4Sc=(maxBinBAr4M[0]<=maxBinBAr4M[1])?maxBinBAr4M[0]:maxBinBAr4M[1];
  Int_t* plotBAr4bM1=(Int_t *)(plotBArM[0]);
  Int_t* plotBAr4bM2=(Int_t *)(plotBArM[1]);
  Int_t* plotBArSc4b=(Int_t *)plotBArSc;
  for(Int_t n4=minBinBAr4Sc; n4<=maxBinBAr4Sc; n4++)
      plotBArSc4b[n4]=plotBAr4bM1[n4] & plotBAr4bM2[n4];

  scanPlotInMod(0,clusArrM1,nClsArrM1);
  if(nClsArrM1>1) mergeClusInMod(0,nClsArrM1,clusArrM1);

  scanPlotInMod(1,clusArrM2,nClsArrM2);
  if(nClsArrM2>1) mergeClusInMod(1,nClsArrM2,clusArrM2);

  scanPlotInSeg1(0,plotBArSc,clusArr,nClsArr);

  for(Int_t nCl=0; nCl<nClsArrM1; nCl++)
    if(clusArrM1[nCl].status) clusArrM1[nCl].calcClusParam(xHStep2,yHStep2);
  for(Int_t nCl=0; nCl<nClsArrM2; nCl++)
    if(clusArrM2[nCl].status) clusArrM2[nCl].calcClusParam(xHStep2,yHStep2);

  if(nClsArr>0) {
    if(nClsArr>1) mergeClusInSeg(0,nClsArr,clusArr);
    fillClusCat(-2,0);
  }
  clearPrMod(0);
  clearPrMod(1);
}

void HMdcLookUpTbSec::scanPlotInMod(Int_t mod,
    HMdcCluster* clArr,Int_t& nClArr) {
  Char_t maxClS=minAmp[mod]-1;
  nClArr=0;
  UChar_t* hPlModM=hPlMod[mod];
  cPlotBAr=plotBArM[mod];
  Int_t* plotBAr4b=(Int_t*) cPlotBAr;
  Int_t minBin=minBinBAr4M[mod];
  Short_t* clusInd=(mod==0) ? clusIndM1:clusIndM2;
  for(Int_t n4=maxBinBAr4M[mod]; n4>=minBin; n4--) {
    if(plotBAr4b[n4]==0) continue;
    UChar_t *b1=cPlotBAr+n4*4;
    UChar_t *b2=b1+3;
    Int_t nBin4=n4*32;
    Int_t nLBin=33;
    while((nLBin=HMdcBArray::prevAndUnset(b1,b2,nLBin))>=0) {
      //---New cluster---------------------------------------
      clus=&clArr[nClArr];
      Int_t nBin=nBin4+nLBin;     
      initCluster(nBin);
      while ((nBin=stack->pop()) >= 0) {
        clusInd[nBin]=nClArr;
        addBinInCluster(nBin, HMdcBArray::getNSet(hPlModM[nBin])-maxClS);
      }
      //-Filling of cluster---------------------------------------------------
      clus->calcXY();
      fillWiresList(mod,clus->lCells1);
      nClArr++;
      if(nClArr >= clusArrSize) break;
    }
    if(nClArr >= clusArrSize) {
      Warning("scanClusInMod"," Num. of clusters in sector %i > max\n",sector);
      memset(cPlotBAr,0,sizeBAr);
      break;
    }
  }
}

void HMdcLookUpTbSec::scanPlotInSeg1(Int_t seg, UChar_t* plotBAr,
    HMdcCluster* clArr, Int_t& nClArr) {
  Int_t m1=seg*2;
  Int_t m2=m1+1;
  Char_t maxClS=minAmp[m1]+minAmp[m2]-1;
  nClArr=0;
  UChar_t *hPlModF=(seg==0) ? hPlMod[0]:hPlMod[2];
  UChar_t *hPlModS=(seg==0) ? hPlMod[1]:hPlMod[3];
  cPlotBAr=plotBAr;
  Int_t* plotBAr4b=(Int_t*) cPlotBAr;
  for(Int_t n4=maxBinBAr4Sc; n4>=minBinBAr4Sc; n4--) {
    if(plotBAr4b[n4]==0) continue;
    UChar_t *b1=cPlotBAr+n4*4;
    UChar_t *b2=b1+3;
    Int_t nBin4=n4*32;
    Int_t nLBin=33;
    while((nLBin=HMdcBArray::prevAndUnset(b1,b2,nLBin))>=0) {
      //---New cluster---------------------------------------
      clus=&clArr[nClArr];
      Int_t nBin=nBin4+nLBin;
      initClusterT2(nBin);
      while ((nBin=stack->pop()) >= 0) addBinInCluster(nBin,
          HMdcBArray::getNSet(hPlModF[nBin])+
          HMdcBArray::getNSet(hPlModS[nBin])-maxClS);
      //-Filling of cluster--------------------------------------------------- 
      clus->calcXY();
      for(Int_t mod=m1;mod<=m2;mod++) fillWiresList(mod,clus->lCells1);
      nClArr++;
      if(nClArr >= clusArrSize) break;
    }
    if(nClArr >= clusArrSize) {
      Warning("scanClusInSeg"," Num. of clusters in sector %i > max\n",sector);
      memset(cPlotBAr,0,sizeBAr);
      break;
    }
  }
}

void HMdcLookUpTbSec::findClusInSec(Short_t minAm){
  for(Int_t mod=0;mod<4;mod++) 
      if(minAmp[mod]>0 && maxAmp[mod]<minAmp[mod]) return;
  if(minAm<3) return;
  makeSPlot();
  Char_t maxClS=minAm-1;
  nClsArr=0;
  cPlotBAr=plotBArSc;
  Int_t* plotBAr4b=(Int_t*) cPlotBAr;
  //---Cluster finding-------------------------------------
  for(Int_t n4=maxBinBAr4Sc; n4>=minBinBAr4Sc; n4--) {
    if(plotBAr4b[n4]==0) continue;
    UChar_t *b1=cPlotBAr+n4*4;
    UChar_t *b2=b1+3;
    Int_t nBin4=n4*32;
    Int_t nLBin=33;
    while((nLBin=HMdcBArray::prevAndUnset(b1,b2,nLBin))>=0) {
      //---New cluster---------------------------------------
      clus=&clusArr[nClsArr];
      Int_t nBin=nBin4+nLBin;    
      initCluster(nBin);
      while ((nBin=stack->pop()) >= 0) addBinInCluster(nBin,
          HMdcBArray::getNSet(hPlMod[0][nBin])+
          HMdcBArray::getNSet(hPlMod[1][nBin])+
          HMdcBArray::getNSet(hPlMod[2][nBin])+
          HMdcBArray::getNSet(hPlMod[3][nBin])-maxClS);
      //-End cluster------------------------------------------------------
      clus->calcXY();
      for(Int_t mod=0;mod<4;mod++) {
        if(minAmp[mod]==0) continue;
        if(mod<2) fillWiresList(mod,clus->lCells1);
        else  fillWiresList(mod,clus->lCells2);
      }
      nClsArr++;
      if(nClsArr >= clusArrSize) break;
    }
    if(nClsArr >= clusArrSize) {
      Warning("findClusInSec"," Num. of clusters in sector %i > max\n",sector);
      memset(cPlotBAr,0,sizeBAr);
      break;
    }
  }
  if(nClsArr>0) {
    if(nClsArr>1) mergeClusInSec(nClsArr, clusArr);
    fillClusCat(-nMods,-1);
  }
  for(Int_t mod=0;mod<4;mod++) if(minAmp[mod]>0) clearPrMod(mod);
}

void HMdcLookUpTbSec::mergeClusInMod(Int_t mod, Int_t nClArr,
    HMdcCluster* clArr) {
  Int_t nClus=nClArr;
  Int_t fLay=(mod==0 || mod==2) ? 0:6;
  Int_t lLay=fLay+5;
  while(nClus>1) {
    Bool_t nomerg=kTRUE;
    for(Int_t cl1=0; cl1<nClArr-1; cl1++) {
      HMdcCluster& cls1=clArr[cl1];
      if(!cls1.status) continue;
      HMdcList12GroupCells& listCells1=cls1.lCells1;
      for(Int_t cl2=cl1+1; cl2<nClArr; cl2++) {
        HMdcCluster& cls2=clArr[cl2];
        if(!cls2.status) continue;
        HMdcList12GroupCells& listCells2=cls2.lCells1;
        Float_t dY=cls1.y-cls2.y;
        if(dY>100.) break;  //????????? biylo zakomentirovano ???
        if(fabs(dY) > 30.) continue;                   // 10. mm !???
        if(fabs(cls1.x-cls2.x) > 100.) continue;        // 40. mm !???
        if(listCells1.compare(&listCells2,fLay,lLay)<4) continue;
        cls1.addClus(cls2,kFALSE);
        nomerg=kFALSE;
        nClus--;
      }
    }
    if(nomerg || nClus==1) break;
  }
}

void HMdcLookUpTbSec::mergeClusInSeg(Int_t seg, Int_t nClArr,
    HMdcCluster* clArr) {
  Int_t nClus=nClArr;
  while(nClus>1) {
    Bool_t nomerg=kTRUE;
    for(Int_t cl1=0; cl1<nClArr-1; cl1++) {
      HMdcCluster& cls1=clArr[cl1];
      if(!cls1.status) continue;
      HMdcList12GroupCells& listCells1=cls1.lCells1;
      for(Int_t cl2=cl1+1; cl2<nClArr; cl2++) {
        HMdcCluster& cls2=clArr[cl2];
        if(!cls2.status) continue;

        if(cls1.clusMod1 != cls2.clusMod1) continue;
        if(cls1.clusMod2 != cls2.clusMod2) continue;

        HMdcList12GroupCells& listCells2=cls2.lCells1;
        Float_t dY=cls1.y-cls2.y;
        if(dY>100.) break;  //????????? biylo zakomentirovano ???
        if(fabs(dY) > 30.) continue;                   // 10. mm !???
        if(fabs(cls1.x-cls2.x) > 100.) continue;       // 40. mm !???
        if(listCells1.compare(&listCells2,0, 5)<4) continue;
        if(listCells1.compare(&listCells2,6,11)<4) continue;
        cls1.addClus(cls2,kFALSE);
        nomerg=kFALSE;
        nClus--;
      }
    }
    if(nomerg || nClus==1) break;
  }
}

void HMdcLookUpTbSec::mergeClusInSec(Int_t nClArr, HMdcCluster* clArr) {
  Int_t nClus=nClArr;
  Int_t m1=minAmp[0];
  Int_t m2=minAmp[1];
  Int_t m3=minAmp[2];
  Int_t m4=minAmp[3];
  while(nClus>1) {
    Bool_t nomerg=kTRUE;
    for(Int_t cl1=0; cl1<nClArr-1; cl1++) {
      HMdcCluster& cls1=clArr[cl1];
      if(!cls1.status) continue;
      HMdcList12GroupCells& listCells1=cls1.lCells1;
      HMdcList12GroupCells& listCells2=cls1.lCells2;
      for(Int_t cl2=cl1+1; cl2<nClArr; cl2++) {
        HMdcCluster& cls2=clArr[cl2];
        if(!cls2.status) continue;
        HMdcList12GroupCells* listCells1s=&(cls2.lCells1);
        HMdcList12GroupCells* listCells2s=&(cls2.lCells2);
        Float_t dY=cls1.y-cls2.y;
        if(dY>100.) break;  //????????? biylo zakomentirovano ???
        if(fabs(dY) > 30.) continue;                   // 10. mm !???
        if(fabs(cls1.x-cls2.x) > 100.) continue;       // 40. mm !???
        if(m1>0 && listCells1.compare(listCells1s,0, 5)<4) continue;
        if(m2>0 && listCells1.compare(listCells1s,6,11)<4) continue;
        if(m3>0 && listCells2.compare(listCells2s,0, 5)<4) continue;
        if(m4>0 && listCells2.compare(listCells2s,6,11)<4) continue;
        cls1.addClus(cls2,kTRUE);
        nomerg=kFALSE;
        nClus--;
      }
    }
    if(nomerg || nClus==1) break;
  }
}

void HMdcLookUpTbSec::fillClusCat(Int_t mod, Int_t segp) {
  Int_t seg=(segp>=0) ? segp : 0;
  Int_t m1=minAmp[seg*2];
  Int_t m2=minAmp[seg*2+1];
  Int_t m3=0;
  Int_t m4=0;
  if(mod>=0) {
    if(mod==seg*2) m2=0;
    else if(mod==seg*2+1) m1=0;
  } else if(segp<0) {
    m3=minAmp[2];
    m4=minAmp[3];
  }
  // Filling containers:
  for(Int_t cl=0; cl<nClsArr; cl++) {
    HMdcCluster& cls=clusArr[cl];
    if(!cls.status) continue;
    getClusterSlot(seg,cls.lCells1);
    if(!fClus) return;
    cls.calcClusParam(xHStep2,yHStep2);
    nClusters++;
    fClus->setMod(mod);
    fClus->setTypeClFinder(typeClFinder);
    fClus->setMinCl(m1,m2);
    fClus->setPrPlane(prPlane.A(),prPlane.B(),prPlane.D());
    fClus->setTarget(target[0],eTarg[0],target[1],eTarg[1],target[2],eTarg[2]);
    cls.fillClus(fClus,0);
    if(fClus->isGeant()) ((HMdcClusSim*)fClus)->calcTrList();
    if(segp<0) {
      HMdcClus* fClus0=fClus;
      getClusterSlot(1,cls.lCells2);
      if(!fClus) return;
      fClus->setMod(mod);
      fClus->setTypeClFinder(typeClFinder);
      fClus->setMinCl(m3,m4);
      fClus->setPrPlane(prPlane.A(),prPlane.B(),prPlane.D());
      fClus->setTarget(target[0],eTarg[0],target[1],eTarg[1],
          target[2],eTarg[2]);
      cls.fillClus(fClus,1);
      fClus->setIndexParent(fClus0->getIndex());
      Int_t ind=fClus->getIndex();
      fClus0->setIndexChilds(ind,ind);
      if(fClus->isGeant()) ((HMdcClusSim*)fClus)->getNTracks();
    }
  }
}

void HMdcLookUpTbSec::findClusInMod(Int_t mod, Short_t minAm){
  if(maxAmp[mod]<minAmp[mod]) return;
  if(minAm<3) return;
  makeModPlot(mod,minAm);
  Char_t maxClS=minAmp[mod]-1;
  nClsArr=0;
  UChar_t* hPlModM=hPlMod[mod];
  cPlotBAr=plotBArM[mod];
  Int_t* plotBAr4b=(Int_t*) cPlotBAr;
  Int_t minBin=minBinBAr4M[mod];
  for(Int_t n4=maxBinBAr4M[mod]; n4>=minBin; n4--) {
    if(plotBAr4b[n4]==0) continue;
    UChar_t *b1=cPlotBAr+n4*4;
    UChar_t *b2=b1+3;
    Int_t nBin4=n4*32;
    Int_t nLBin=33;
    while((nLBin=HMdcBArray::prevAndUnset(b1,b2,nLBin))>=0) {
      //---New cluster---------------------------------------
      clus=&clusArr[nClsArr];
      Int_t nBin=nBin4+nLBin;    
      initCluster(nBin);
      while ((nBin=stack->pop()) >= 0) addBinInCluster(nBin,
            HMdcBArray::getNSet(hPlModM[nBin])-maxClS);
      //-Filling of cluster---------------------------------------------------
      clus->calcXY();
      fillWiresList(mod,clus->lCells1);
      nClsArr++;
      if(nClsArr >= clusArrSize) break;
    }
    if(nClsArr >= clusArrSize) {
      Warning("scanClusInMod"," Num. of clusters in sector %i > max\n",sector);
      memset(cPlotBAr,0,sizeBAr);
      break;
    }
  }
  if(nClsArr>0) {
    if(nClsArr>1) mergeClusInMod(mod,nClsArr,clusArr);
    fillClusCat(mod,mod/2);
  }
  clearPrMod(mod);
}

//---------------------------------------------------------
HMdcLookUpTb* HMdcLookUpTb::fMdcLookUpTb=0;

HMdcLookUpTb::HMdcLookUpTb(const char* name,const char* title,
                           const char* context)
            : HParSet(name,title,context) {
  // constructor creates an array of pointers of type HMdcLookUpTbSec
  strcpy(detName,"Mdc");
  fGetCont=HMdcGetContainers::getObject();
  if( !fGetCont ) return;
  array = new TObjArray(6);
  fMdcDet       = fGetCont->getMdcDetector();
  fMdcGeomPar   = fGetCont->getMdcGeomPar();
  fSpecGeomPar  = fGetCont->getSpecGeomPar();
  fLayerGeomPar = fGetCont->getMdcLayerGeomPar();
  fSizesCells   =  HMdcSizesCells::getObject();
  setPar(319, kFALSE);   // Project plot size is the same for all sectors!!!
  fMdcClusCat = HMdcGetContainers::getCatMdcClus(kTRUE);
  targLenInc[0]=targLenInc[1]=0.;
}

void HMdcLookUpTb::setPar(Int_t inBinX, Bool_t isCOff) {
  nBinX=(inBinX%2 == 0) ? inBinX:inBinX+1;
  if(nBinX<100) nBinX=320;
  nBinY=(nBinX*278)/100;
  if(nBinY%2 == 0) nBinY++;
  isCoilOff=isCOff;
}

Bool_t HMdcLookUpTb::initContainer() {
  // It is called from "reinit" of reconstractor!
  if( !fMdcDet || !HMdcGetContainers::isInited(fMdcGeomPar) ||
      !HMdcGetContainers::isInited(fSpecGeomPar) ||
      !HMdcGetContainers::isInited(fLayerGeomPar) ||
      !fSizesCells->initContainer()) return kFALSE;

  if( !status && (fSizesCells->hasChanged() || fSpecGeomPar->hasChanged() ||
      fMdcGeomPar->hasChanged() || fLayerGeomPar->hasChanged()) ) {
    changed=kTRUE;
    if(!fMdcClusCat) return kFALSE;
    for (Int_t sec = 0; sec < 6; sec++) {
      if(!fMdcDet->isSectorActive(sec))  continue;
      if( !(*array)[sec] ) {
        Int_t nSegs=(fMdcDet->getModule(sec,0)||fMdcDet->getModule(sec,1)) ?1:0;
        if(isCoilOff)
          nSegs=(fMdcDet->getModule(sec,2) || fMdcDet->getModule(sec,3)) ?2:1;
        else if(nSegs==0) continue;
        (*array)[sec] = new HMdcLookUpTbSec(sec,nSegs,nBinX,nBinY);
        (*this)[sec].fClusCat=fMdcClusCat;
      }
      // initialization of container ---
      if(!calcPrPlane(sec)) return kFALSE;
      if(!calcTarget(sec)) return kFALSE;
      if(!calcPlotSize(sec)) return kFALSE;
      if(!calcLookUpTb(sec)) return kFALSE;
      // --------------------------------
    }
    if(versions[1]<0 || versions[2]<0) versions[1]=versions[2]=0;
    else versions[2]++;
  } else changed=kFALSE;
  return kTRUE;
}

HMdcLookUpTb* HMdcLookUpTb::getExObject() {
  return fMdcLookUpTb;
}

HMdcLookUpTb* HMdcLookUpTb::getObject() {
  if(!fMdcLookUpTb) fMdcLookUpTb=new HMdcLookUpTb();
  return fMdcLookUpTb;
}

void HMdcLookUpTb::deleteCont() {
  if(fMdcLookUpTb) {
    delete fMdcLookUpTb;
    fMdcLookUpTb=0;
  }
}

HMdcLookUpTb::~HMdcLookUpTb() {
  // destructor
  HMdcGetContainers::deleteCont();
  HMdcSizesCells::deleteCont();
  if(array) {
    array->Delete();
    delete array;
  }
}

Int_t HMdcLookUpTb::getSize(void) {
  // return the size of the pointer array
  return array->GetEntries();
}

void HMdcLookUpTb::clear() {
  for(Int_t s=0;s<6;s++) {
   if( (*array)[s] ) {
      HMdcLookUpTbSec& sec=(*this)[s];
      sec.clearwk();
    }
  }
}

Int_t  HMdcLookUpTb::findClusters(Int_t *imax){
  Int_t ntot=0;
  for(Int_t sec=0;sec<6;sec++) {
    if( (*array)[sec] ) {
      HMdcLookUpTbSec& fSec=(*this)[sec];
      ntot+=fSec.findClusters(imax+4*sec);
    }
  }
  return ntot;
}

TH2C* HMdcLookUpTbSec::fillTH2C(Char_t* name, Char_t* title, Int_t type,
    Int_t bining){
  // Proj.plot filling without cluster finding.
  for(Int_t mod=0; mod<4; mod++) {
    cFMod=&((*this)[mod]);
    if( !cFMod ) continue;
    cHPlModM=hPlMod[mod];
    if(!hPlMod) continue;
    cXMinM=xMin[mod];
    cXMaxM=xMax[mod];
    Int_t iLay=mod*6;
    for(Int_t layer=0;layer<6;layer++) makeLayProjV1(layer,iLay+layer);
  }
  
  if(!hist) {
    plBining=(bining) ? 2:1;
    hist=new TH2C(name,title,nBinX/plBining,xLow,xUp,nBinY/plBining,yLow,yUp);
  } else {
    hist->Reset();
    hist->SetName(name);
    hist->SetTitle(title);
  }
  hist->SetMinimum(0.);
  hist->Fill(0.,0.,0);
  Float_t maximum=6.;
  if(type==2) maximum=4.5;
  else if(type==0) maximum=(hPlMod[2]) ? 24.:12.;
  else if(type==1) maximum=(hPlMod[2]) ? 60.:18.;
  hist->SetMaximum(maximum);
  if(nHitsTot==0) return hist;
  for (Int_t nx=0; nx<nBinX; nx+=plBining) {
    for (Int_t ny=0; ny<nBinY; ny+=plBining) {
      Int_t nBin=ny*nBinX+nx;
      Int_t nBin2=nBin+nBinX;
      UChar_t binM1=hPlMod[0][nBin];
      UChar_t binM2=hPlMod[1][nBin];
      if(plBining==2) {
        binM1 |= hPlMod[0][nBin+1] | hPlMod[0][nBin2] | hPlMod[0][nBin2+1];
        binM2 |= hPlMod[1][nBin+1] | hPlMod[1][nBin2] | hPlMod[1][nBin2+1];
      }
      UChar_t binM3=0;
      UChar_t binM4=0;
      if(hPlMod[2]) {
        binM3=hPlMod[2][nBin];
        binM4=hPlMod[3][nBin];
        if(plBining==2) {
          binM3|=hPlMod[2][nBin+1] | hPlMod[2][nBin2] | hPlMod[2][nBin2+1];
          binM4|=hPlMod[3][nBin+1] | hPlMod[3][nBin2] | hPlMod[3][nBin2+1];
        }
      }
      if( binM1==0 && binM2==0 && binM3==0 && binM4==0) continue;
      Int_t color=0;
      if(type==1) {
        if(binM4) color=HMdcBArray::getNSet(binM4)+18;
        else if(binM3) color+=HMdcBArray::getNSet(binM3)+12;
        else if(binM2) color+=HMdcBArray::getNSet(binM2)+6;
        else if(binM1) color+=HMdcBArray::getNSet(binM1);
      } else if(type==2) {
        if(binM4) color=4;
        else if(binM3) color=3;
        else if(binM2) color=2;
        else if(binM1) color=1;
      } else {
        color=HMdcBArray::getNSet(binM1)+HMdcBArray::getNSet(binM2)+
              HMdcBArray::getNSet(binM3)+HMdcBArray::getNSet(binM4);
      }
      hist->Fill(xBinsPos[nx],yBinsPos[ny], color);
    }
  }
  for(Int_t mod=0;mod<4;mod++) clearPrMod(mod);
  return hist;
}

Bool_t HMdcLookUpTb::calcLookUpTb(Int_t sec) {
  HMdcLookUpTbSec& fsec=(*this)[sec];
  HMdcSizesCellsSec& fSCSec = (*fSizesCells)[sec];
  HMdcTrap cellSize;
  //Modules:
  for(Int_t mod=0; mod<4; mod++) {
    if( !&(fsec[mod]) ) continue;
    for(Int_t layer=0; layer<6; layer++) {
      HMdcSizesCellsLayer& fSCLayer = fSCSec[mod][layer];
      Int_t nCells=fSCLayer.getSize();
      if(nCells<1) continue;
      HMdcLookUpTbLayer& fLUpTbLayer=fsec[mod][layer];
      //Cells:
      for(Int_t cell=0; cell<nCells; cell++) {
        HMdcLookUpTbCell& fLUpTbCell=fLUpTbLayer[cell];
        if( !&fLUpTbCell ) return kFALSE;
        fLUpTbCell.clear();
        HMdcSizesCellsCell& fSizesCellsCell = fSCLayer[cell];
        if( !&fSizesCellsCell || !fSizesCellsCell.getStatus()) continue;
        if(!fSizesCells->getCellVol(sec,mod,layer,cell,cellSize)) continue;
        HMdcTrapPlane cellProj;
        calcCellProj(sec,cellSize,cellProj);
        //Lines 0,1,2,3 - lines from points 0-1, 1-2, 2-3, 3-0 respectively
        Int_t yBinMinT=0;
        Int_t yBinMaxT=0;
        for(Int_t i=1; i<4; i++) {
          if(cellProj[yBinMinT].gety() > cellProj[i].gety()) yBinMinT=i;
	  if(cellProj[yBinMaxT].gety() < cellProj[i].gety()) yBinMaxT=i;
        }
        yBinMinT=(Int_t)((cellProj[yBinMinT].gety()-fsec.yLow)/fsec.yStep);
        yBinMaxT=(Int_t)((cellProj[yBinMaxT].gety()-fsec.yLow)/fsec.yStep);
        if(yBinMinT<1 || yBinMinT>fsec.nBinY-2) {
          Warning("calcLookUpTb","S%iM%iL%iC%i yBinMinT=%i isn't in 1 - %i !",
              sec+1,mod+1,layer+1,cell+1,yBinMinT,fsec.nBinY-2);
          yBinMinT=(yBinMinT<1) ? 1:fsec.nBinY-2;
        }
        if(yBinMaxT<1 || yBinMaxT>fsec.nBinY-1) {
          Warning("calcLookUpTb","S%iM%iL%iC%i  yBinMaxT=%i isn't in 1 - %i !",
              sec+1,mod+1,layer+1,cell+1,yBinMaxT,fsec.nBinY-2);
          yBinMaxT=(yBinMaxT<1) ? 1:fsec.nBinY-2;
        }
        fLUpTbCell.init(yBinMinT,yBinMaxT);
        Int_t nx1t,nx2t;
        nx1t=+1000000000;
        nx2t=-1000000000;
        for(Int_t ny=yBinMinT; ny<=yBinMaxT; ny++) {
          Double_t y=(ny+1)*fsec.yStep+fsec.yLow;
          Double_t x1,x2;
          calcX(cellProj, y, x1, x2);	
          Int_t nx1i=(Int_t)((x1-fsec.xLow)/fsec.xStep);
	  Int_t nx2i=(Int_t)((x2-fsec.xLow)/fsec.xStep);
	  if(nx1i < nx1t) nx1t=nx1i;
	  if(nx2i > nx2t) nx2t=nx2i;
          if(nx1i<1 || nx1i>fsec.nBinX-2) {
            Warning("calcLookUpTb","S%iM%iL%iC%i  nx1i=%i isn't in 1 - %i !",
              sec+1,mod+1,layer+1,cell+1,nx1i,fsec.nBinX-2);
            nx1i=(nx1i<1) ? 1:fsec.nBinX-2;
          }
          if(nx2i<1 || nx2i>fsec.nBinX-2) {
            Warning("calcLookUpTb","S%iM%iL%iC%i  nx2i=%i isn't in 1 - %i !",
              sec+1,mod+1,layer+1,cell+1,nx2i,fsec.nBinX-2);
            nx2i=(nx2i<1) ? 1:fsec.nBinX-2;
          }
          //Fill look up table:
          fLUpTbCell.addLine(nx1t,nx2t);
	  if(ny==yBinMaxT-1 && yBinMinT!=yBinMaxT) {
            if(!fLUpTbCell.addLine(nx1t,nx2t)) return kFALSE;
            ny++;
          }
	  nx1t=nx1i;
	  nx2t=nx2i;	
        }
      }
    }
  }
  return kTRUE;
}

void HMdcLookUpTb::calcX(HMdcTrapPlane& pr, Double_t y,
                           Double_t &xLow, Double_t &xUp) {
  xLow=1.0e+30;
  xUp=-1.0e+30;
  for(Int_t i=0; i<4; i++){
    Int_t u1=i;
    Int_t u2=(u1<3) ? u1+1:0;
    if(fabs(pr[u1].gety()-pr[u2].gety())<1.e-4) continue;
    if(pr[u1].getx()>pr[u2].getx()) {
      Int_t u=u1;  u1=u2;  u2=u;
    }
    Double_t x=(y-pr[u1].gety())/(pr[u1].gety()-pr[u2].gety())*
              (pr[u1].getx()-pr[u2].getx())+pr[u1].getx();
    if(x<pr[u1].getx() || x>pr[u2].getx()) continue;
    if(x<xLow) xLow=x;
    if(x>xUp) xUp=x;
  }
  //Attention pr[0]&pr[1] mast be >= pr[2]&pr[3] (0-3; 1-2)!!!!
  if(xLow<-2000) xLow=(pr[2].getx() < pr[3].getx()) ? pr[2].getx():pr[3].getx();
  if(xUp> 2000) xUp=(pr[0].getx() > pr[1].getx()) ? pr[0].getx():pr[1].getx();
}


Bool_t HMdcLookUpTb::calcPrPlane(Int_t sec){
  //Sectors:
  //  Calc. of projection plane:
  //  It's the plane between MDC0 and MDC1
  //  par=0.0 -> plane=layer 6, MDC1
  //  par=1.0 -> plane=layer 1, MDC2

  HGeomTransform coordSys;
  Double_t par=-100.;
  //  Double_t par=1.0;   // <<<===================== use MDC1
  Int_t nmod=0;
  Int_t firstMod=-1;
  Int_t lastMod=0;
  Int_t nMaxMods=(isCoilOff) ? 4:2;
  for(Int_t m=0; m<nMaxMods; m++) {
    if( fMdcDet->getModule(sec,m) ) {
      nmod++;
      if(firstMod<0) firstMod=m;
      lastMod=m;
    }
  }
  if(nmod==0) return kFALSE;
  Int_t iPar=-1;
  if(nmod==1) iPar=firstMod;
  else if(firstMod==0) par=0.575;   // determ. of par !!!
  else iPar=firstMod;                 // for beem time nov2000 with mdc 3
  Double_t transNew[3];
  if(iPar<0) {
    Double_t tanTheta1 = 0.0;
    Double_t tanTheta2=0.0;
    Double_t zMdc1=0.0;
    Double_t zMdc2=0.0;
    Double_t tmp;
    coordSys=*((*fSizesCells)[sec][firstMod][5].getSecTrans());
    tmp=(coordSys.getRotMatrix())(8);
    tmp=1./(tmp*tmp);
    tanTheta1=sqrt(tmp-1.);
    zMdc1=(coordSys.getTransVector()).getZ()*tmp;

    coordSys=*((*fSizesCells)[sec][lastMod][0].getSecTrans());
    tmp=(coordSys.getRotMatrix())(8);  tmp=1./(tmp*tmp);
    tanTheta2=sqrt(tmp-1.);
    zMdc2=(coordSys.getTransVector()).getZ()*tmp;

    if(zMdc1==0.0 && zMdc2==0.0) {
      Error("calcMdc12Proj","No MDC volume found.");
      return kFALSE;
    }
    printf(
    "\n===> Sec.%i: Using plane between MDC%i and MDC%i (p=%f) as projection plane\n",
        sec+1,firstMod+1,lastMod+1,par);

    Double_t tanTheta=tanTheta1*(1.-par)+tanTheta2*par;
    Double_t cosTheta=1/sqrt(tanTheta*tanTheta+1.);
    Double_t sinTheta=tanTheta*cosTheta;
    Double_t newRot[9];
    newRot[0]=1;
    newRot[1]=newRot[2]=newRot[3]=newRot[6]=0;
    newRot[4]=newRot[8]=cosTheta;                           // cos(Theta)
    newRot[5]=sinTheta;                                     // sin(Theta)
    newRot[7]=-sinTheta;                                    // -sin(Theta)
    HGeomRotation rotMatProj(newRot);
    transNew[0]=0;
    transNew[2]=(zMdc1*(1.-par)+zMdc2*par)*cosTheta*cosTheta; // mm
    transNew[1]=transNew[2]*tanTheta;

    coordSys.setRotMatrix(newRot);
    coordSys.setTransVector(transNew);
  }
  else {
    coordSys=*((*fSizesCells)[sec][iPar].getSecTrans());
    printf("\n===> Sec.%i: Using middle plane of MDC%i as projection plane\n",
                              sec+1,iPar+1);
  }
  ((*this)[sec]).prPlane.setPrPlane(coordSys);
  return kTRUE;
}

Bool_t HMdcLookUpTb::calcTarget(Int_t sec){
  //Geting target parameters
  Int_t nT=fSpecGeomPar->getNumTargets()-1;
  if( nT < 0 ) {
    Error("calcTarget","Number of targets = %i!",nT+1);
    return kFALSE;
  }
  HMdcLookUpTbSec& fsec=(*this)[sec];
  fsec.targVc[0]=(fSpecGeomPar->getTarget( 0)->getTransform()).getTransVector();
  fsec.targVc[0].setZ( fsec.targVc[0].getZ() +
      fSpecGeomPar->getTarget( 0)->getPoint(0)->getZ() - targLenInc[0]);
  fsec.targVc[1]=(fSpecGeomPar->getTarget(nT)->getTransform()).getTransVector();
  fsec.targVc[1].setZ( fsec.targVc[1].getZ() +
      fSpecGeomPar->getTarget(nT)->getPoint(2)->getZ() + targLenInc[1]);
  const HGeomTransform* trans=(*fSizesCells)[sec].getLabTrans();
  if(&trans == 0) return kFALSE;
  fsec.targVc[0]=trans->transTo(fsec.targVc[0]);
  fsec.targVc[1]=trans->transTo(fsec.targVc[1]);
  HGeomVector tar(fsec.targVc[0]+fsec.targVc[1]);
  tar*=0.5;
  HGeomVector errTar(fsec.targVc[0]-fsec.targVc[1]);
  errTar*=0.5;
  errTar=errTar.abs();
  for(Int_t i=0;i<3;i++) {
    fsec.target[i]=tar(i);
    fsec.eTarg[i]=errTar(i);
  }
  return kTRUE;
}

Bool_t HMdcLookUpTb::calcPlotSize(Int_t sec){
  HMdcLookUpTbSec& fsec=(*this)[sec];
  HGeomVector senVol;
  HGeomVector cross;
  fsec.xLow=fsec.yLow=1.0e+30;
  fsec.xUp=fsec.yUp=-1.0e+30;
  for(Int_t mod=0; mod<4; mod++) {
    if( &(fsec[mod]) == 0 ) continue;
    HGeomCompositeVolume* fComVol=fGetCont->getGeomCompositeVolume(mod);
    if(!fComVol) {
      Error("calcPlotSize","GeomCompositeVolume for MDC%i Sec.%i is absent.",
          mod+1,sec+1);
      return kFALSE;
    }

    for(Int_t layer=0; layer<6; layer+=5) {
      HGeomVolume* fGeomVolLay=fComVol->getComponent(layer);
      if(!fGeomVolLay || fGeomVolLay->getNumPoints()!=8) return kFALSE;
      HMdcSizesCellsLayer& fSizesCellsLay=(*fSizesCells)[sec][mod][layer];
      if( !&fSizesCells ) return kFALSE;
      HMdcLayerGeomParLay& fLayerGeomParLay=(*fLayerGeomPar)[sec][mod][layer];
      if( !&fLayerGeomParLay )return kFALSE;
      Double_t dstCathPl=(layer/2-1)*fLayerGeomParLay.getCatDist()*0.5;
      for(Int_t point=0; point<4; point++) {
        senVol=*(fGeomVolLay->getPoint(point)); // mm!
        senVol.setZ(dstCathPl);
        senVol=fSizesCellsLay.getSecTrans()->transFrom(senVol);
        for(Int_t targ=0; targ<2; targ++) {
          fsec.prPlane.calcIntersection(fsec.targVc[targ],
                       senVol-fsec.targVc[targ],cross);
          if(cross(0)<fsec.xLow) fsec.xLow=cross(0);
          else if(cross(0)>fsec.xUp) fsec.xUp=cross(0);
          if(cross(1)<fsec.yLow) fsec.yLow=cross(1);
          else if(cross(1)>fsec.yUp) fsec.yUp=cross(1);
        }
      }
    }
  }
  Double_t del=(fsec.xUp-fsec.xLow)*0.02;
  fsec.xLow -= del;
  fsec.xUp  += del;
  del=(fsec.yUp-fsec.yLow)*0.01;
  fsec.yLow -= del;
  fsec.yUp  += del;
  fsec.xStep=(fsec.xUp-fsec.xLow)/(Double_t)(fsec.nBinX-3);
  fsec.yStep=(fsec.yUp-fsec.yLow)/(Double_t)(fsec.nBinY-3);
  fsec.yLow -= fsec.yStep*1.5;            //must be bin=0 at limits of plot
  fsec.yUp += fsec.yStep*1.5;
  fsec.xLow -= fsec.xStep*1.5;
  fsec.xUp += fsec.xStep*1.5;
  fsec.xStep=(fsec.xUp-fsec.xLow)/(Double_t)(fsec.nBinX);
  fsec.yStep=(fsec.yUp-fsec.yLow)/(Double_t)(fsec.nBinY);
  Double_t xHStep=fsec.xStep/2;
  Double_t yHStep=fsec.yStep/2;
  fsec.xHStep2=xHStep*xHStep;
  fsec.yHStep2=yHStep*yHStep;
  for(Int_t n=0; n<fsec.nBinX; n++)
      fsec.xBinsPos[n]=((Double_t)n)*fsec.xStep+fsec.xLow+xHStep;
  for(Int_t n=0; n<fsec.nBinY; n++)
      fsec.yBinsPos[n]=((Double_t)n)*fsec.yStep+fsec.yLow+yHStep;
  return kTRUE;
}

Bool_t HMdcLookUpTb::calcCellProj(Int_t sec,HMdcTrap& cellSize,
    HMdcTrapPlane& cellProj) {
  // Calculation of cell projection on the proj. plane.
  // (For sector coordinat system only !!!)
  HGeomVector pProj[16];
  HMdcLookUpTbSec& fsec=(*this)[sec];
  for (Int_t i=0; i<8; i++) {
    for (Int_t j=0; j<2; j++) fsec.prPlane.calcIntersection(fsec.targVc[j],
      cellSize[i]-fsec.targVc[j],pProj[j*8+i]);
  }
  Double_t dist[2][8];
  Double_t sig=0.001;
  //Lines 1,5 for minmum:
  Bool_t indb=(fabs(pProj[0](0)-pProj[3](0)) > sig) ? kTRUE:kFALSE;
  for(Int_t i=0; i<2; i++) {
    for(Int_t j=0; j<2; j++) {
      dist[i][j*4]=pProj[i*8+j*4](1);
      if(indb) dist[i][j*4]=calcDist(pProj[i*8+j*4],pProj[i*8+j*4+3]);	
    }
  }
  //Lines 1,5 for maximum:
  Bool_t indu=(fabs(pProj[1](0)-pProj[2](0)) > sig) ? kTRUE:kFALSE;
  for(Int_t i=0; i<2; i++) {
    for(Int_t j=0; j<2; j++) {
      dist[i][j*4+1]=pProj[i*8+j*4+1](1);
      if(indu) dist[i][j*4+1]=calcDist(pProj[i*8+j*4+1],pProj[i*8+j*4+2]);
    }
  }
  // Minimum y:
  Int_t tb=0; Int_t pb=0;
  if(dist[1][0] < dist[tb][pb])  tb=1;
  if(dist[1][4] < dist[tb][pb]) {tb=1; pb=4;}
  if(dist[0][4] < dist[tb][pb]) {tb=0; pb=4;}
  //Maximum y:	
  Int_t tu=0; Int_t pu=1;
  if(dist[1][1] > dist[tu][pu])  tu=1;
  if(dist[1][5] > dist[tu][pu]) {tu=1; pu=5;}
  if(dist[0][5] > dist[tu][pu]) {tu=0; pu=5;}
  for(Int_t i=0; i<2; i++) {
   dist[i][0]=fabs(calcDist(pProj[i*8+0],pProj[i*8+1])); //line 0-1;
   dist[i][4]=fabs(calcDist(pProj[i*8+4],pProj[i*8+5])); //line 4-5;
   dist[i][2]=fabs(calcDist(pProj[i*8+2],pProj[i*8+3])); //line 2-3;
   dist[i][6]=fabs(calcDist(pProj[i*8+6],pProj[i*8+7])); //line 2-3
  }
  // Minimum-minmum for lines 0-1;4-5;  2-3;6-7:
  Int_t tl=0; Int_t pl=0;
  if(dist[1][0] > dist[tl][pl])  tl=1;
  if(dist[1][4] > dist[tl][pl]) {tl=1; pl=4;}
  if(dist[0][4] > dist[tl][pl]) {tl=0; pl=4;}
  Int_t tr=0; Int_t pr=2;
  if(dist[1][2] > dist[tr][pr])  tr=1;
  if(dist[1][6] > dist[tr][pr]) {tr=1; pr=6;}
  if(dist[0][6] > dist[tr][pr]) {tr=0; pr=6;}
  if(indb) {
    calcPoint(cellProj[0],pProj[tb*8+pb],pProj[tb*8+pb+3],
	                  pProj[tl*8+pl],pProj[tl*8+pl+1]);
    calcPoint(cellProj[3],pProj[tb*8+pb],pProj[tb*8+pb+3],
	                  pProj[tr*8+pr],pProj[tr*8+pr+1]);
  }
  if(indu) {
    calcPoint(cellProj[1],pProj[tu*8+pu],pProj[tu*8+pu+1],
	                  pProj[tl*8+pl],pProj[tl*8+pl+1]);
    calcPoint(cellProj[2],pProj[tu*8+pu],pProj[tu*8+pu+1],
	                  pProj[tr*8+pr],pProj[tr*8+pr+1]);
  }
  Double_t x,y;
  if(!indb) {
    x=(pProj[tb*8+pb](0)+pProj[tb*8+pb+3](0))*0.5;
    y=(pProj[tb*8+pb](1)+pProj[tb*8+pb+3](1))*0.5;
    if(y<0. || y>10000.) {
      Error("calcCellProj","Sec.%i p0,3  x=%f y=%f(!)",sec+1,x,y);
      return kFALSE;
    }
    cellProj[0].set(x,y);
    cellProj[3].set(x,y);
  }
  if(!indu) {
    x=(pProj[tb*8+pu](0)+pProj[tb*8+pu+1](0))*0.5;
    y=(pProj[tb*8+pu](1)+pProj[tb*8+pu+1](1))*0.5;
    if(y<0. || y>10000.) {
      Error("calcCellProj","Sec.%i p1,2  x=%f y=%f(!)",sec+1,x,y);
      return kFALSE;
    }
    cellProj[1].set(x,y);
    cellProj[2].set(x,y);
  }
  return kTRUE;
}

Double_t HMdcLookUpTb::calcDist(HGeomVector& p1, HGeomVector& p2) {
  return p1(1)-p1(0)*(p1(1)-p2(1))/(p1(0)-p2(0));
}

void HMdcLookUpTb::calcPoint(HMdcPointPlane& proj,
  HGeomVector& p1l1, HGeomVector& p2l1, HGeomVector& p1l2, HGeomVector& p2l2) {
  Double_t al1=(p1l1(1)-p2l1(1))/(p1l1(0)-p2l1(0));
  Double_t bl1=p1l1(1)-al1*p1l1(0);
  Double_t al2=(p1l2(1)-p2l2(1))/(p1l2(0)-p2l2(0));
  Double_t bl2=p1l2(1)-al2*p1l2(0);
  Double_t x=(bl2-bl1)/(al1-al2);
  Double_t y=al2*x+bl2;
  proj.set(x,y);
}