using namespace std;
#include "hmdc34clfinder.h"
#include "hades.h"
#include "hmatrixcategory.h"
#include "hruntimedb.h"
#include "hspectrometer.h"
#include "hmdcdetector.h"
#include <iostream>
#include <iomanip>
#include "hmdctrackddef.h"
#include "hmdcgeomobj.h"
#include "hmdcgetcontainers.h"
#include "hmdccal1sim.h"
#include "hmdcclussim.h"
#include "hmdccluster.h"
#include "hmdcseg.h"
#include "TH2.h"
#include "TMath.h"
#include "hkickplane2.h"
#include "hgeomtransform.h"
#include "hmdcgeomstruct.h"
#include "hspecgeompar.h"
#include "hgeomvolume.h"
#include "hgeomcompositevolume.h"
#include "hmdclayergeompar.h"
#include "hmdcsizescells.h"
#include "hmdcgeompar.h"
#include "hmdcclfnstack.h"
#include "hmdctrackdset.h"
#include "hmdcgeanttrack.h"
#include <stdlib.h>
ClassImp(HMdcProjPlot)
ClassImp(HMdc34ClFinderLayer)
ClassImp(HMdc34ClFinderMod)
ClassImp(HMdc34ClFinderSec)
ClassImp(HMdc34ClFinder)
using namespace TMath;
UChar_t* HMdcProjPlot::weights=0;
Int_t HMdcProjPlot::wtArrSize=0;
HMdcProjPlot::HMdcProjPlot(UChar_t mSeg, Int_t inBinX, Int_t inBinY) {
nBinX = inBinX;
nBinY = inBinY;
size = nBinX*nBinY;
size = (size/32 + ((size%32 > 0) ? 1:0))*32;
xMinL = new Short_t [nBinY];
xMaxL = new Short_t [nBinY];
yMinL = new Short_t [nBinX];
if(mSeg&1) plModF=new UChar_t [size];
else plModF=0;
if(mSeg&2) plModS=new UChar_t [size];
else plModS=0;
clearArr();
sizeBAr = size/8;
plotBArSc = new UChar_t [sizeBAr];
plotBArSc4b = (Int_t *)plotBArSc;
clearBitArr();
rootPlot = 0;
rootPlotV2 = 0;
if(size>wtArrSize) {
if(weights) delete [] weights;
weights = new UChar_t [size];
wtArrSize = size;
}
xBinsPos = new Double_t [nBinX];
yBinsPos = new Double_t [nBinY];
}
void HMdcProjPlot::clearArr(void) {
if(plModF) memset(plModF,0,size);
if(plModS) memset(plModS,0,size);
}
void HMdcProjPlot::clearBitArr(void) {
memset(plotBArSc,0,sizeBAr);
}
HMdcProjPlot::~HMdcProjPlot() {
if(xMinL) delete [] xMinL;
if(xMaxL) delete [] xMaxL;
if(yMinL) delete [] yMinL;
xMinL = xMaxL = yMinL = 0;
if(plModF) delete [] plModF;
if(plModS) delete [] plModS;
plModF=plModS=0;
if(plotBArSc) delete [] plotBArSc;
plotBArSc=0;
if(rootPlot) {
delete rootPlot;
rootPlot = 0;
}
if(rootPlotV2) {
delete rootPlotV2;
rootPlotV2 = 0;
}
if(weights) {
delete [] weights;
weights = 0;
wtArrSize = 0;
}
if(xBinsPos) delete [] xBinsPos;
if(yBinsPos) delete [] yBinsPos;
xBinsPos = 0;
yBinsPos = 0;
}
void HMdcProjPlot::setEdges(Double_t iyMin, Double_t ixMinD, Double_t ixMaxD,
Double_t iyMax, Double_t ixMin, Double_t ixMax) {
stY = (iyMax-iyMin)/(nBinY-2);
stX = (ixMax-ixMin)/(nBinX-2);
yMin = iyMin-stY;
yMax = iyMax+stY;
xMin = ixMin-stX;
xMax = ixMax+stX;
xMinD = ixMinD;
xMaxD = ixMaxD;
Double_t aL23 = (xMinD-xMin)/(yMin-yMax);
Double_t aL01 = (xMaxD-xMax)/(yMin-yMax);
for(Int_t ny=0; ny<nBinY; ny++) {
Double_t yl = (ny+1)*stY+yMin;
xMinL[ny] = Short_t((((yl-yMax)*aL23+xMin)-xMin)/stX);
if(xMinL[ny] < 1) xMinL[ny]=1;
xMaxL[ny] = Short_t((((yl-yMax)*aL01+xMax)-xMin)/stX);
if(xMaxL[ny] > nBinX-2) xMaxL[ny]=nBinX-2;
}
for(Int_t nx=0; nx<nBinX; nx++) {
yMinL[nx]=nBinY;
for(Int_t ny=nBinY-2; ny>=0; ny--) {
if(nx<xMinL[ny] || nx>xMaxL[ny]) break;
yMinL[nx] = ny==0 ? 1 : ny;
}
}
xFirstBin = 0.5*stX+xMin;
yFirstBin = 0.5*stY+yMin;
for(Int_t n=0; n<nBinX; n++) xBinsPos[n]=(((Double_t)n)+0.5)*stX+xMin;
for(Int_t n=0; n<nBinY; n++) yBinsPos[n]=(((Double_t)n)+0.5)*stY+yMin;
}
void HMdcProjPlot::print(void) {
printf("----- Project plane: -----\n");
HMdcPlane::print();
printf(" Sizes: nBinX=%i, nBinY=%i, stepX=%g, stepY=%g\n",nBinX,nBinY,stX,stY);
printf(" Points (x:y) n. 0=(%.1f:%.1f) 1=(%.1f:%.1f)\n",xMaxD,yMin,xMax,yMax);
printf(" 2=(%.1f:%.1f) 3=(%.1f:%.1f)\n",xMin,yMax,xMinD,yMin);
}
TH2C* HMdcProjPlot::getPlot(const Char_t* name,const Char_t* title){
if(!rootPlot) rootPlot=new TH2C(name,title,nBinX,xMin,xMax,nBinY,yMin,yMax);
else {
rootPlot->Reset();
rootPlot->SetName(name);
rootPlot->SetTitle(title);
}
rootPlot->SetMaximum(12.);
rootPlot->SetMinimum(0.);
rootPlot->Fill(0.,0.,0);
for (Int_t nBin=0; nBin<size; nBin++) {
Int_t wt=0;
if(plModF && plModF[nBin]) wt+=HMdcBArray::getNSet(plModF[nBin]);
if(plModS && plModS[nBin]) wt+=HMdcBArray::getNSet(plModS[nBin]);
if(wt==0) continue;
rootPlot->Fill(xBinsPos[nBin%nBinX],yBinsPos[nBin/nBinX],(Stat_t)wt);
}
return rootPlot;
}
TH2C* HMdcProjPlot::getPlotV2(const Char_t* name,const Char_t* title){
if(!rootPlotV2) rootPlotV2=new TH2C(name,title,nBinX,xMin,xMax,nBinY,yMin,yMax);
else {
rootPlotV2->Reset();
rootPlotV2->SetName(name);
rootPlotV2->SetTitle(title);
}
rootPlotV2->SetMaximum(18.);
rootPlotV2->SetMinimum(0.);
rootPlotV2->Fill(0.,0.,0);
for (Int_t nBin=0; nBin<size; nBin++) {
Int_t wt=0;
if(plModF && plModF[nBin]) wt+=HMdcBArray::getNSet(plModF[nBin]);
if(plModS && plModS[nBin]) wt+=HMdcBArray::getNSet(plModS[nBin])+6;
if(wt==0) continue;
rootPlotV2->Fill(xBinsPos[nBin/nBinX],yBinsPos[nBin%nBinX],(Stat_t)wt);
}
return rootPlotV2;
}
HMdc34ClFinderLayer::HMdc34ClFinderLayer(Int_t sec, Int_t mod, Int_t layer) {
cellShUp=cellShDown=0;
HMdcGetContainers* fGetCont=HMdcGetContainers::getObject();
if( !fGetCont ) return;
HMdcGeomStruct* pMdc=fGetCont->getMdcGeomStruct();
if( !pMdc ) return;
nCells = ((*pMdc)[sec][mod])[layer];
cellShUp = new Short_t [nCells];
cellShDown = new Short_t [nCells];
yBin = 0;
nBinX = 0;
xBin1 = 0;
xBin2 = 0;
cells = 0;
}
HMdc34ClFinderLayer::~HMdc34ClFinderLayer() {
if(yBin) delete [] yBin;
if(xBin1) delete [] xBin1;
if(xBin2) delete [] xBin2;
if(cellShUp) delete [] cellShUp;
if(cellShDown) delete [] cellShDown;
yBin = 0;
xBin1 = 0;
xBin2 = 0;
cellShUp = 0;
cellShDown = 0;
nBinX = 0;
}
Bool_t HMdc34ClFinderLayer::createArrayBins(Short_t nBins) {
if(nBinX == nBins) return kTRUE;
if(yBin) delete [] yBin;
nBinX = nBins;
yBin = new Short_t [nBins];
if(yBin==0) return kFALSE;
return kTRUE;
}
HMdc34ClFinderMod::HMdc34ClFinderMod(Int_t sec, Int_t mod) {
array = new TObjArray(6);
for (Int_t layer = 0; layer < 6; layer++) {
(*array)[layer] = new HMdc34ClFinderLayer(sec, mod, layer);
}
}
HMdc34ClFinderMod::~HMdc34ClFinderMod() {
if(array) {
array->Delete();
delete array;
}
}
Int_t HMdc34ClFinderMod::getNCells(void) {
Int_t nHits = 0;
for(Int_t lay=0; lay<6; lay++) nHits += (*this)[lay].cells->getNCells();
return nHits;
}
HMdc34ClFinderSec::HMdc34ClFinderSec(Int_t sec, Int_t inBinX, Int_t inBinY) {
sector = sec;
HMdcDetector* fMdcDet = HMdcGetContainers::getObject()->getMdcDetector();
const Int_t* nLM = HMdcTrackDSet::getTrFnNLayersInMod()+sector*4;
for(Int_t mod=0;mod<4;mod++)
lMods[mod]=(fMdcDet->getModule(sector,mod)&&nLM[mod]>0) ? 1:0;
mSeg[0] = (lMods[0] ? 1:0) + (lMods[1] ? 2:0);
mSeg[1] = (lMods[2] ? 1:0) + (lMods[3] ? 2:0);
array = new TObjArray(4);
typeClFinder = 0;
for (Int_t mod = 2; mod < 4; mod++) {
if(lMods[mod]) (*array)[mod] = new HMdc34ClFinderMod(sector,mod);
}
prPlotSeg2 = 0;
if(mSeg[1] > 0) {
prPlotSeg2 = new HMdcProjPlot(mSeg[1],inBinX,inBinY);
nBinX = prPlotSeg2->nBinX;
nearbyBins[0] = -1;
nearbyBins[1] = +1;
nearbyBins[2] = -nBinX;
nearbyBins[3] = +nBinX;
nearbyBins[4] = -1-nBinX;
nearbyBins[5] = -1+nBinX;
nearbyBins[6] = 1-nBinX;
nearbyBins[7] = 1+nBinX;
}
isGeant = HMdcGetContainers::isGeant();
locClus.set(2,sector,1);
for(Int_t mod=2;mod<4;mod++) {
if(lMods[mod]) {
xCMin[mod] = new Int_t [prPlotSeg2->nBinY];
xCMax[mod] = new Int_t [prPlotSeg2->nBinY];
for(Int_t y=0;y<prPlotSeg2->nBinY;y++) {
xCMin[mod][y] = prPlotSeg2->size;
xCMax[mod][y] = -1;
}
} else {
xCMin[mod] = 0;
xCMax[mod] = 0;
}
}
xMinM2 = xCMin[2];
xMaxM2 = xCMax[2];
pClustersArrs = HMdcClustersArrs::getObject();
clusArr = pClustersArrs->getArray1();
if(clusArr == 0) clusArr = pClustersArrs->createArray1(500);
clusArrSize = pClustersArrs->getArray1Size();
mergeFlag = kTRUE;
}
HMdc34ClFinderSec::~HMdc34ClFinderSec() {
if(array) {
array->Delete();
delete array;
}
if(prPlotSeg2) delete prPlotSeg2;
prPlotSeg2 = 0;
for(Int_t mod=2; mod<4; mod++) {
if(xCMin[mod]) delete [] xCMin[mod];
xCMin[mod] = 0;
if(xCMax[mod]) delete [] xCMax[mod];
xCMax[mod] = 0;
}
}
void HMdc34ClFinderSec::clear(void) {
maxAmp[2] = -1;
maxAmp[3] = -1;
notEnoughWrs = kTRUE;
}
void HMdc34ClFinderSec::setMinBin(Int_t *mBin) {
for(Int_t i=0;i<4;i++) minAmp[i]=mBin[i];
notEnoughWrs = testMaxAmp();
}
Bool_t HMdc34ClFinderSec::testMaxAmp(void) {
for(Int_t mod=2;mod<4;mod++) {
maxAmp[mod]=0;
if(minAmp[mod]<=0) continue;
HMdc34ClFinderMod& fMod=(*this)[mod];
for(Int_t lay=0;lay<6;lay++) {
HMdc34ClFinderLayer& fLay=fMod[lay];
if(fLay.cells->getFirstCell()>=0) maxAmp[mod]++;
}
}
if(minAmp[2]>0 && minAmp[3]>0) {
if(minAmp[2]>maxAmp[2] || minAmp[3]>maxAmp[3]) return kTRUE;
} else {
if(minAmp[2]>maxAmp[2]) return kTRUE;
if(minAmp[3]>maxAmp[3]) return kTRUE;
}
return kFALSE;
}
Int_t HMdc34ClFinderSec::findClustersSeg2(HMdcSeg* fSeg, HMdcClus* pClus,
Int_t *mBin){
if(mBin!=0) setMinBin(mBin);
if(notEnoughWrs) return 0;
Float_t z0,errZ0, r0,errR0, theta,errTh, phi,errPh;
fSeg->getZR(z0,errZ0,r0,errR0);
fSeg->getThetaPhi(theta,errTh,phi,errPh);
Double_t cosPhiS = cos(phi+Pi()*0.5);
Double_t sinPhiS = sin(phi+Pi()*0.5);
Double_t x0 = r0*cosPhiS;
Double_t y0 = r0*sinPhiS;
Double_t cosPhi = cos(phi);
Double_t sinPhi = sin(phi);
Double_t cosTh = cos(theta);
Double_t sinTh = sin(theta);
Double_t dirX = cosPhi*sinTh;
Double_t dirY = sinPhi*sinTh;
Double_t dirZ = cosTh;
Double_t errX0 = sqrt(Power(cosPhiS*errR0,2)+Power(y0*errPh,2));
Double_t errY0 = sqrt(Power(sinPhiS*errR0,2)+Power(x0*errPh,2));
Double_t errDirX = sqrt(Power(sinPhi*sinTh*errPh,2)+
Power(cosPhi*cosTh*errTh,2));
Double_t errDirY = sqrt(Power(cosPhi*sinTh*errPh,2)+
Power(sinPhi*cosTh*errTh,2));
Double_t errDirZ = Abs(sinTh*errTh);
seg1[0].setXYZ(x0+errX0,y0+errY0,z0-errZ0);
dirSeg1[0].setXYZ(dirX+errDirX,dirY+errDirY,dirZ-errDirZ);
seg1[1].setXYZ(x0-errX0,y0+errY0,z0-errZ0);
dirSeg1[1].setXYZ(dirX-errDirX,dirY+errDirY,dirZ-errDirZ);
seg1[2].setXYZ(x0+errX0,y0-errY0,z0-errZ0);
dirSeg1[2].setXYZ(dirX+errDirX,dirY-errDirY,dirZ+errDirZ);
seg1[3].setXYZ(x0-errX0,y0-errY0,z0-errZ0);
dirSeg1[3].setXYZ(dirX-errDirX,dirY-errDirY,dirZ+errDirZ);
seg1[4].setXYZ(x0,y0,z0);
dirSeg1[4].setXYZ(dirX,dirY,dirZ);
indexPar = pClus->getOwnIndex();
indexFCh = -1;
indexLCh = -2;
Int_t nClustCh = findClustersSeg2();
if(nClustCh>0) pClus->setIndexChilds(indexFCh,indexLCh);
return nClustCh;
}
Int_t HMdc34ClFinderSec::findClustersSeg2(HMdcClus* pClus, Int_t *mBin){
if(pClus->getIOSeg() != 0) return 0;
if(mBin!=0) setMinBin(mBin);
if(notEnoughWrs) return 0;
Double_t dirX = pClus->getX();
Double_t dirY = pClus->getY();
Double_t dirZ = pClus->getZ();
Double_t errDirX = pClus->getErrX();
Double_t errDirY = pClus->getErrY();
Double_t x0 = pClus->getXTarg();
Double_t y0 = pClus->getYTarg();
Double_t z0 = pClus->getZTarg();
Double_t errX0 = pClus->getErrXTarg();
Double_t errY0 = pClus->getErrYTarg();
Double_t errZ0 = pClus->getErrZTarg();
seg1[0].setXYZ(x0-errX0,y0-errY0,z0-errZ0);
dirSeg1[0].setXYZ(dirX+errDirX,dirY+errDirY,
pClus->getZOnPrPlane(dirX+errDirX,dirY+errDirY));
seg1[1].setXYZ(x0+errX0,y0-errY0,z0-errZ0);
dirSeg1[1].setXYZ(dirX-errDirX,dirY+errDirY,
pClus->getZOnPrPlane(dirX-errDirX,dirY+errDirY));
seg1[2].setXYZ(x0-errX0,y0+errY0,z0-errZ0);
dirSeg1[2].setXYZ(dirX+errDirX,dirY-errDirY,
pClus->getZOnPrPlane(dirX+errDirX,dirY-errDirY));
seg1[3].setXYZ(x0+errX0,y0+errY0,z0-errZ0);
dirSeg1[3].setXYZ(dirX-errDirX,dirY-errDirY,
pClus->getZOnPrPlane(dirX-errDirX,dirY-errDirY));
seg1[4].setXYZ(x0,y0,z0);
dirSeg1[4].setXYZ(dirX,dirY,dirZ);
for(Int_t np=0; np<5; np++) dirSeg1[np] -= seg1[np];
indexPar = pClus->getOwnIndex();
indexFCh = -1;
indexLCh = -2;
Int_t nClustCh = findClustersSeg2();
if(nClustCh>0) pClus->setIndexChilds(indexFCh,indexLCh);
return nClustCh;
}
Int_t HMdc34ClFinderSec::findClustersSeg2(void){
for(Int_t np=0; np<5; np++) {
if( !(fkick->calcIntersection(seg1[np],dirSeg1[np],segOnKick[np]))) {
Error("findClustersSeg2",
"Sec.%i: The point %i of HMDCSeg hit on the kick plane can't be calculated.",
sector+1,np);
return 0;
}
if(np==4) {
dirSeg1[4].setXYZ( (Abs(segOnKick[0](0)-segOnKick[1](0))+
Abs(segOnKick[2](0)-segOnKick[3](0)))*0.25,
(Abs(segOnKick[0](1)-segOnKick[2](1))+
Abs(segOnKick[1](1)-segOnKick[3](1)))*0.25,
(Abs(segOnKick[0](2)-segOnKick[2](2))+
Abs(segOnKick[1](2)-segOnKick[3](2)))*0.25);
break;
}
if(segOnKick[np](0)<prPlotSeg2->xMin || segOnKick[np](0)>prPlotSeg2->xMax) {
Warning("findClustersSeg2","S%i: The X=%g of the point %i out of the region (%g - %g).",
sector+1, segOnKick[np](0), np, prPlotSeg2->xMin, prPlotSeg2->xMax);
return 0;
}
Double_t tmp = prPlotSeg2->D() - prPlotSeg2->A()*segOnKick[np](0);
al[np] = tmp-segOnKick[np](2);
bl[np] = tmp*segOnKick[np](1);
cl[np] = segOnKick[np](2)+prPlotSeg2->B()*segOnKick[np](1);
nbX[np] = Int_t((segOnKick[np](0)-prPlotSeg2->xMin)/prPlotSeg2->stX);
}
#if DEBUG_LEVEL>2
printf("---Sec.%i---- The region of HMdcSeg-hit on the kick plane:\n",sector+1);
for(Int_t np=0; np<4; np++) segOnKick[np].print();
printf("x-x(0,1)=%g x-x(2,3)=%g y-y(0,2)=%g y-y(1,3)=%g\n",
segOnKick[0](0)-segOnKick[1](0),segOnKick[2](0)-segOnKick[3](0),
segOnKick[0](1)-segOnKick[2](1),segOnKick[1](1)-segOnKick[3](1));
#endif
realTypeClFinder=0;
if(minAmp[2]>0 && minAmp[3]>0) {
if(typeClFinder==2) {
realTypeClFinder=2;
makeSeg2PlotTp2();
} else makeSeg2Plot();
}
else if(minAmp[2]>0) makeModS2Plot(2);
else if(minAmp[3]>0) makeModS2Plot(3);
else return 0;
Int_t nClusters = scanPlotSeg2();
if(realTypeClFinder==2) clearPrSegTp2();
else clearPrMod();
return nClusters;
}
void HMdc34ClFinderSec::clearPrMod(Int_t mod) {
Int_t m1=(mod==-1) ? 2:mod;
Int_t m2=(mod==-1) ? 3:mod;
if(m1<2 || m2>3) return;
for(Int_t m=m1;m<=m2;m++) {
UChar_t* hPlModM=(m==2) ? prPlotSeg2->plModF:prPlotSeg2->plModS;
if(!hPlModM) continue;
Int_t* xMaxM=xCMax[m];
Int_t* xMinM=xCMin[m];
for(Int_t y=0;y<prPlotSeg2->nBinY;y++) {
if(xMaxM[y]<0) continue;
memset(hPlModM+xMinM[y],0,xMaxM[y]-xMinM[y]+1);
xMinM[y]=prPlotSeg2->size;
xMaxM[y]=-1;
}
}
}
void HMdc34ClFinderSec::clearPrSegTp2(void) {
for(Int_t y=0;y<prPlotSeg2->nBinY;y++) {
if(xMaxM2[y]<0) continue;
memset(prPlotSeg2->plModF+xMinM2[y],0,xMaxM2[y]-xMinM2[y]+1);
memset(prPlotSeg2->plModS+xMinM2[y],0,xMaxM2[y]-xMinM2[y]+1);
xMinM2[y]=prPlotSeg2->size;
xMaxM2[y]=-1;
}
}
void HMdc34ClFinderSec::makeModS2Plot(Int_t mod) {
Int_t module = mod;
cPlMod=(module==2) ? prPlotSeg2->plModF:prPlotSeg2->plModS;
if(cPlMod==0) return;
cFMod=&((*this)[module]);
cXMinM = xCMin[module];
cXMaxM = xCMax[module];
Int_t maxB4Sc = 0;
Int_t minB4Sc = prPlotSeg2->size;
Int_t minAm = minAmp[module];
if(minAm<3) minAm=3;
UChar_t* weights = prPlotSeg2->weights;
UChar_t* plotBArSc = prPlotSeg2->plotBArSc;
Int_t layList[6] = {2,3,1,4,0,5};
Int_t nFiredLay = 0;
for(Int_t il=0; il<6; il++) {
Int_t lay=layList[il];
if(!calcLayerProjVar(lay)) continue;
nFiredLay++;
if(maxAmp[module]-nFiredLay+1>=minAm) {
if(nFiredLay<minAm) makeLayProjV1(lay);
else {
while(test) {
if(setNextCell()) {
setYLimits();
for(Int_t ny=ny1; ny<=ny2; ny++) {
if(setXLimits(ny,prPlotSeg2->xMinL[ny],
prPlotSeg2->xMaxL[ny])) continue;
if(nbF<cXMinM[ny]) cXMinM[ny]=nbF;
if(nbL>cXMaxM[ny]) cXMaxM[ny]=nbL;
UChar_t* bt=cPlMod+nbF;
for(Int_t nb=nbF; nb<=nbL; nb++) {
*bt |= add;
Char_t wt = HMdcBArray::getNSet(bt);
bt++;
if( wt<minAm ) continue;
weights[nb] = wt;
HMdcBArray::set(plotBArSc,nb);
if(nb<minB4Sc) minB4Sc=nb;
if(nb>maxB4Sc) maxB4Sc=nb;
}
}
shDown=shDownN;
}
shUp=shUpN;
}
}
} else {
if(nFiredLay<minAm) makeLayProjV2(lay);
else {
while(test) {
if(setNextCell()) {
setYLimits();
for(Int_t ny=ny1; ny<=ny2; ny++) {
if(cXMaxM[ny]<0) continue;
if(setXLimits(ny,cXMinM[ny]%nBinX,cXMaxM[ny]%nBinX)) continue;
UChar_t* bt=cPlMod+nbF;
for(Int_t nb=nbF; nb<=nbL; nb++) {
*bt |= add;
Char_t wt=HMdcBArray::getNSet(bt);
bt++;
if( wt<minAm ) continue;
weights[nb]=wt;
HMdcBArray::set(plotBArSc,nb);
if(nb<minB4Sc) minB4Sc=nb;
if(nb>maxB4Sc) maxB4Sc=nb;
}
}
shDown=shDownN;
}
shUp=shUpN;
}
}
}
}
prPlotSeg2->minBin4Sc=minB4Sc/32;
prPlotSeg2->maxBin4Sc=maxB4Sc/32;
}
Bool_t HMdc34ClFinderSec::calcLayerProjVar(Int_t lay) {
cFLay = &((*cFMod)[lay]);
cell = cFLay->cells->getFirstCell();
test = kFALSE;
if(cell<0) return kFALSE;
if(cell >= cFLay->nCells) {
Warning("calcLayerProjVar","Cell %i >= num.cells(=%i)",cell,cFLay->nCells);
cell = -1;
return kFALSE;
}
nBinYM2 = prPlotSeg2->nBinY-2;
shUp = Max(calcYbin(0,0,cell+1),calcYbin(0,1,cell+1));
shDown = Min(calcYbin(1,0,cell),calcYbin(1,1,cell));
add = 1<<lay;
nYLinesM1 = cFLay->nYLines-1;
xBin1L = cFLay->xBin1;
xBin2L = cFLay->xBin2;
shUpN = 0;
shDownN = 0;
test = kTRUE;
cFLay->cellShUp[cell] = shUp;
cFLay->cellShDown[cell] = shDown;
return kTRUE;
}
Bool_t HMdc34ClFinderSec::setNextCell(void) {
cell=cFLay->cells->next(cell);
if(cell>=0) {
shDownN = Min(calcYbin(1,0,cell),calcYbin(1,1,cell));
shUpN = Max(calcYbin(0,0,cell+1),calcYbin(0,1,cell+1));
cFLay->cellShUp[cell] = shUpN;
cFLay->cellShDown[cell] = shDownN;
if(shDownN-shUp<=1) return kFALSE;
} else {
test = kFALSE;
shDownN = shUp+1;
}
return kTRUE;
}
void HMdc34ClFinderSec::setYLimits(void) {
shDown+=cFLay->yFirst;
shUp+=cFLay->yFirst;
ny1=(shDown<=0) ? 1 : shDown;
ny2=(nYLinesM1+shUp>nBinYM2) ? nBinYM2 : nYLinesM1+shUp;
}
Bool_t HMdc34ClFinderSec::setXLimits(Int_t ny, Short_t nx1,Short_t nx2) {
Int_t iy1=ny-shUp;
Int_t iy2=ny-shDown;
if(cFLay->wOrType>0) {
if(iy1>=0 && xBin1L[iy1]>nx1) nx1=xBin1L[iy1];
if(iy2<=nYLinesM1 && xBin2L[iy2]<nx2) nx2=xBin2L[iy2];
} else {
if(iy2<=nYLinesM1 && xBin1L[iy2]>nx1) nx1=xBin1L[iy2];
if(iy1>=0 && xBin2L[iy1]<nx2) nx2=xBin2L[iy1];
}
if(nx2<nx1) return kTRUE;
Int_t shift=ny * nBinX;
nbF=nx1+shift;
nbL=nx2+shift;
return kFALSE;
}
void HMdc34ClFinderSec::makeLayProjV1(Int_t lay) {
while(test) {
if(setNextCell()) {
setYLimits();
for(Int_t ny=ny1; ny<=ny2; ny++) {
if(setXLimits(ny,prPlotSeg2->xMinL[ny],prPlotSeg2->xMaxL[ny])) continue;
if(nbF<cXMinM[ny]) cXMinM[ny] = nbF;
if(nbL>cXMaxM[ny]) cXMaxM[ny] = nbL;
UChar_t* bin = cPlMod+nbF;
UChar_t* binMax = cPlMod+nbL;
for(; bin<=binMax; bin++) *bin |= add;
}
shDown=shDownN;
}
shUp=shUpN;
}
}
void HMdc34ClFinderSec::makeLayProjV1b(Int_t lay) {
while(test) {
if(setNextCell()) {
setYLimits();
for(Int_t ny=ny1; ny<=ny2; ny++) {
if(xMaxM2[ny]<0) continue;
if(setXLimits(ny,xMinM2[ny]%nBinX,xMaxM2[ny]%nBinX)) continue;
if(nbF<cXMinM[ny]) cXMinM[ny]=nbF;
if(nbL>cXMaxM[ny]) cXMaxM[ny]=nbL;
UChar_t* bin =cPlMod+nbF;
UChar_t* binMax=cPlMod+nbL;
for(; bin<=binMax; bin++) *bin |= add;
}
shDown=shDownN;
}
shUp=shUpN;
}
}
void HMdc34ClFinderSec::makeLayProjV2(Int_t lay) {
while(test) {
if(setNextCell()) {
setYLimits();
for(Int_t ny=ny1; ny<=ny2; ny++) {
if(cXMaxM[ny]<0) continue;
if(setXLimits(ny,cXMinM[ny]%nBinX,cXMaxM[ny]%nBinX)) continue;
UChar_t* bin =cPlMod+nbF;
UChar_t* binMax=cPlMod+nbL;
for(; bin<=binMax; bin++) *bin |= add;
}
shDown=shDownN;
}
shUp=shUpN;
}
}
void HMdc34ClFinderSec::makeSeg2Plot(void) {
Int_t mBm3 = minAmp[2];
Int_t maxB4Sc = 0;
Int_t minB4Sc = prPlotSeg2->size;
UChar_t* weights = prPlotSeg2->weights;
UChar_t* plotBArSc = prPlotSeg2->plotBArSc;
UChar_t* plModF=prPlotSeg2->plModF;
Int_t layList[6]={2,3,1,4,0,5};
for(Int_t mod=2; mod<4; mod++) {
cFMod=&((*this)[mod]);
cXMinM=xCMin[mod];
cXMaxM=xCMax[mod];
Int_t minAm=minAmp[mod];
if(minAm<3) minAm=3;
cPlMod=(mod==2) ? prPlotSeg2->plModF:prPlotSeg2->plModS;
if(cPlMod==0) continue;
Int_t nFiredLay=0;
for(Int_t il=0; il<6; il++) {
Int_t lay=layList[il];
if(!calcLayerProjVar(lay)) continue;
nFiredLay++;
if(maxAmp[mod]-nFiredLay+1>=minAm) {
if(mod==2) makeLayProjV1(lay);
else if(nFiredLay<minAm) makeLayProjV1b(lay);
else {
while(test) {
if(setNextCell()) {
setYLimits();
for(Int_t ny=ny1; ny<=ny2; ny++) {
if(setXLimits(ny,prPlotSeg2->xMinL[ny],
prPlotSeg2->xMaxL[ny])) continue;
if(nbF<cXMinM[ny]) cXMinM[ny]=nbF;
if(nbL>cXMaxM[ny]) cXMaxM[ny]=nbL;
UChar_t* bt=cPlMod+nbF;
for(Int_t nb=nbF; nb<=nbL; nb++) {
*bt |= add;
Char_t wt=HMdcBArray::getNSet(bt);
bt++;
if( wt<minAm ) continue;
Char_t wt2=HMdcBArray::getNSet(plModF[nb]);
if(wt2<mBm3) continue;
weights[nb]=wt+wt2;
HMdcBArray::set(plotBArSc,nb);
if(nb<minB4Sc) minB4Sc=nb;
if(nb>maxB4Sc) maxB4Sc=nb;
}
}
shDown=shDownN;
}
shUp=shUpN;
}
}
} else {
if(mod==2 || nFiredLay<minAm) makeLayProjV2(lay);
else {
while(test) {
if(setNextCell()) {
setYLimits();
for(Int_t ny=ny1; ny<=ny2; ny++) {
if(cXMaxM[ny]<0) continue;
if(setXLimits(ny,cXMinM[ny]%nBinX,cXMaxM[ny]%nBinX)) continue;
UChar_t* bt=cPlMod+nbF;
for(Int_t nb=nbF; nb<=nbL; nb++) {
*bt |= add;
Char_t wt=HMdcBArray::getNSet(bt);
bt++;
if( wt<minAm ) continue;
Char_t wt2=HMdcBArray::getNSet(plModF[nb]);
if(wt2<mBm3) continue;
weights[nb]=wt+wt2;
HMdcBArray::set(plotBArSc,nb);
if(nb<minB4Sc) minB4Sc=nb;
if(nb>maxB4Sc) maxB4Sc=nb;
}
}
shDown=shDownN;
}
shUp=shUpN;
}
}
}
}
}
prPlotSeg2->minBin4Sc=minB4Sc/32;
prPlotSeg2->maxBin4Sc=maxB4Sc/32;
}
void HMdc34ClFinderSec::makeSeg2PlotTp2(void) {
Int_t minAm=minAmp[2]+minAmp[3];
if(minAm<3) minAm=3;
Int_t maxAm=maxAmp[2]+maxAmp[3];
Int_t maxB4Sc=0;
Int_t minB4Sc=prPlotSeg2->size;
UChar_t* weights = prPlotSeg2->weights;
UChar_t* plotBArSc = prPlotSeg2->plotBArSc;
UChar_t* plModF = prPlotSeg2->plModF;
UChar_t* plModS = prPlotSeg2->plModS;
Int_t nFiredLay=0;
cXMinM=xMinM2;
cXMaxM=xMaxM2;
Int_t layList[6]={2,3,1,4,0,5};
for(Int_t il=0; il<6; il++) {
Int_t lay=layList[il];
for(Int_t mod=2; mod<4; mod++) {
cFMod=&((*this)[mod]);
Bool_t isMod4 = mod!=2;
cPlMod=(!isMod4) ? prPlotSeg2->plModF:prPlotSeg2->plModS;
if(cPlMod==0) continue;
if(!calcLayerProjVar(lay)) continue;
nFiredLay++;
if(maxAm-nFiredLay+1>=minAm) {
if(nFiredLay<minAm) makeLayProjV1(lay);
else {
while(test) {
if(setNextCell()) {
setYLimits();
for(Int_t ny=ny1; ny<=ny2; ny++) {
if(setXLimits(ny,prPlotSeg2->xMinL[ny],
prPlotSeg2->xMaxL[ny])) continue;
if(nbF<cXMinM[ny]) cXMinM[ny]=nbF;
if(nbL>cXMaxM[ny]) cXMaxM[ny]=nbL;
UChar_t* bt=cPlMod+nbF;
for(Int_t nb=nbF; nb<=nbL; nb++) {
*bt |= add;
bt++;
Char_t wt=HMdcBArray::getNSet(plModF[nb])+
HMdcBArray::getNSet(plModS[nb]);
if(wt<minAm ) continue;
weights[nb]=wt;
HMdcBArray::set(plotBArSc,nb);
if(nb<minB4Sc) minB4Sc=nb;
if(nb>maxB4Sc) maxB4Sc=nb;
}
}
shDown=shDownN;
}
shUp=shUpN;
}
}
} else {
if(nFiredLay<minAm) makeLayProjV2(lay);
else {
while(test) {
if(setNextCell()) {
setYLimits();
for(Int_t ny=ny1; ny<=ny2; ny++) {
if(cXMaxM[ny]<0) continue;
if(setXLimits(ny,cXMinM[ny]%nBinX,cXMaxM[ny]%nBinX)) continue;
UChar_t* bt=cPlMod+nbF;
for(Int_t nb=nbF; nb<=nbL; nb++) {
*bt |= add;
bt++;
Char_t wt=HMdcBArray::getNSet(plModF[nb])+
HMdcBArray::getNSet(plModS[nb]);
if(wt<minAm ) continue;
weights[nb]=wt;
HMdcBArray::set(plotBArSc,nb);
if(nb<minB4Sc) minB4Sc=nb;
if(nb>maxB4Sc) maxB4Sc=nb;
}
}
shDown=shDownN;
}
shUp=shUpN;
}
}
}
}
}
prPlotSeg2->minBin4Sc=minB4Sc/32;
prPlotSeg2->maxBin4Sc=maxB4Sc/32;
}
TH2C* HMdc34ClFinderSec::getPlot(const Char_t* name,const Char_t* title,Int_t ver) {
for(Int_t mod=2; mod<4; mod++) {
cPlMod=(mod==2) ? prPlotSeg2->plModF:prPlotSeg2->plModS;
if(cPlMod==0) continue;
cFMod = &((*this)[mod]);
cXMinM = xCMin[mod];
cXMaxM = xCMax[mod];
for(Int_t lay=0; lay<6; lay++) if(calcLayerProjVar(lay)) makeLayProjV1(lay);
}
TH2C* plt;
if(ver==0) plt = prPlotSeg2->getPlot(name,title);
else plt = prPlotSeg2->getPlotV2(name,title);
clearPrMod();
return plt;
}
Int_t HMdc34ClFinderSec::calcYbin(Int_t upDo, Int_t leRi, Int_t cell) {
Int_t np=cFLay->nPSegOnKick[upDo][leRi];
Double_t xp=segOnKick[np](0);
Double_t yp=cFLay->tgY*xp+cell*cFLay->yStep+cFLay->y0[leRi];
Double_t zp=cFLay->tgZ*xp+cell*cFLay->zStep+cFLay->z0[leRi];
Double_t yPrPl=(al[np]*yp-bl[np]+segOnKick[np](1)*zp)/
(zp+prPlotSeg2->B()*yp-cl[np]);
return Int_t((yPrPl-prPlotSeg2->yMin)/prPlotSeg2->stY) - cFLay->yBin[nbX[np]];
}
Int_t HMdc34ClFinderSec::scanPlotSeg2(void) {
nClsArr=0;
Bool_t useFixedLevel = HMdc34ClFinder::getUseFixedLevelFlag();
Int_t maxB4Sc=prPlotSeg2->maxBin4Sc;
for(Int_t n4=prPlotSeg2->minBin4Sc; n4<=maxB4Sc; n4++) {
if(prPlotSeg2->plotBArSc4b[n4]==0) continue;
UChar_t *b1=prPlotSeg2->plotBArSc+n4*4;
UChar_t *b2=b1+3;
Int_t nBin4=n4*32;
Int_t n=-1;
if(useFixedLevel) while( (n=HMdcBArray::nextAndUnset(b1,b2,n))>=0 &&
calcClusterSeg2FixedLevel(nBin4+n) );
else while( (n=HMdcBArray::nextAndUnset(b1,b2,n))>=0 &&
calcClusterSeg2FloatLevel(nBin4+n) );
if(nClsArr >= clusArrSize) break;
}
mergeClustSeg2();
return fillClusCat();
}
Int_t HMdc34ClFinderSec::fillClusCat(void) {
Int_t nClusters = 0;
Bool_t isFirst = isGeant;
Bool_t useFixedLevel = HMdc34ClFinder::getUseFixedLevelFlag();
for(Int_t cl=0; cl<nClsArr; cl++) if(clusArr[cl].getStatus()) {
if(isFirst) {
HMdcGeantEvent* pGEvent = HMdcGeantEvent::getExObject();
if(pGEvent) pGEvent->clearOSegClus();
isFirst = kFALSE;
}
HMdcCluster& cluster = clusArr[cl];
cluster.calcClusParam();
HMdcClus* clus = (HMdcClus*)fClusCat->getNewSlot(locClus,&indexLCh);
if(!clus) {
Warning("fillClusCat","S.%i No slot HMdcClus available!",sector+1);
return cl-1;
}
if(isGeant) clus = (HMdcClus*)(new(clus) HMdcClusSim(cluster.getLCells2()));
else clus = new(clus) HMdcClus(cluster.getLCells2());
if(indexFCh<0) indexFCh=indexLCh;
nClusters++;
Int_t nLayM1=clus->getNLayersMod(0);
Int_t nLayM2=clus->getNLayersMod(1);
clus->setSecSegInd(sector,1,indexLCh);
if(realTypeClFinder==2 && (nLayM1<minAmp[2] || nLayM2<minAmp[3])) {
clus->setMod( (nLayM1>=nLayM2) ? 2 : 3);
if(useFixedLevel) clus->setTypeClFinder(2);
else clus->setTypeClFinder(128+2);
} else {
clus->setMod((mSeg[1]==3) ? -2 : (mSeg[1]+1) );
if(useFixedLevel) clus->setTypeClFinder(0);
else clus->setTypeClFinder(128);
}
clus->setMinCl(minAmp[2],minAmp[3]);
clus->setPrPlane(prPlotSeg2->A(),prPlotSeg2->B(),prPlotSeg2->D());
clus->setTarg(segOnKick[4]);
clus->setErrTarg(dirSeg1[4]);
clus->setIndexParent(indexPar);
cluster.fillClus(clus,1,isGeant);
}
return nClusters;
}
void HMdc34ClFinderSec::mergeClustSeg2(void) {
if(!mergeFlag) return;
Int_t nClus=nClsArr;
Int_t mergeLevM3 = minAmp[2]>4 ? 4:3;
Int_t mergeLevM4 = minAmp[3]>4 ? 4:3;
while(nClus>1) {
Bool_t nomerg=kTRUE;
for(Int_t cl1=0; cl1<nClsArr-1; cl1++) {
HMdcCluster& cls1=clusArr[cl1];
if(!cls1.getStatus()) continue;
HMdcList12GroupCells& cLCSeg2 = cls1.getLCells2();
for(Int_t cl2=cl1+1; cl2<nClsArr; cl2++) {
HMdcCluster& cls2=clusArr[cl2];
if(!cls2.getStatus()) continue;
HMdcList12GroupCells* cLCSeg2s = &(cls2.getLCells2());
Float_t dY=cls1.getY()-cls2.getY();
if(Abs(dY) > 30.) continue;
if(Abs(cls1.getX()-cls2.getX()) > 100.) continue;
if(realTypeClFinder==2) {
if( cLCSeg2.compare(cLCSeg2s,0,11)<6) continue;
} else {
if(minAmp[2]>0 && cLCSeg2.compare(cLCSeg2s,0, 5)<mergeLevM3) continue;
if(minAmp[3]>0 && cLCSeg2.compare(cLCSeg2s,6,11)<mergeLevM4) continue;
}
cls1.addClus(cls2);
nomerg=kFALSE;
nClus--;
}
}
if(nomerg || nClus==1) break;
}
}
void HMdc34ClFinderSec::initCluster(Int_t nBin,UChar_t amp) {
cluster = &(clusArr[nClsArr]);
cluster->init(1,prPlotSeg2->xFirstBin, prPlotSeg2->yFirstBin,
prPlotSeg2->stX, prPlotSeg2->stY);
Int_t ny = prPlotSeg2->ybin(nBin);
nLMinCl = nLMaxCl = ny;
xMinClLines[ny] = xMaxClLines[ny] = prPlotSeg2->xbin(nBin);
isClstrInited = kTRUE;
addBinInCluster(nBin,amp);
}
void HMdc34ClFinderSec::reinitCluster(Int_t nBin,UChar_t amp) {
cluster->clearBinStat();
addBinInCluster(nBin,amp);
Int_t ny = prPlotSeg2->ybin(nBin);
nLMinCl = nLMaxCl = ny;
xMinClLines[ny] = xMaxClLines[ny] = prPlotSeg2->xbin(nBin);
}
Bool_t HMdc34ClFinderSec::calcClusterSeg2FixedLevel(Int_t nBin) {
initCluster(nBin,prPlotSeg2->weights[nBin]);
stack->init();
do {
for(Int_t ib=0; ib<8; ib++) {
Int_t nb=nBin+nearbyBins[ib];
if(prPlotSeg2->testBin(nb)) {
stack->push(nb);
addBinInCluster(nb,prPlotSeg2->weights[nb]);
}
}
} while((nBin=stack->pop()) >= 0);
if( fillClusterSeg2() ) return kTRUE;
return kFALSE;
}
Bool_t HMdc34ClFinderSec::calcClusterSeg2FloatLevel(Int_t nBin) {
UChar_t amp = prPlotSeg2->weights[nBin];
initCluster(nBin,amp);
stacksArr->init(12);
while (kTRUE) {
for(Int_t ib=0; ib<8; ib++) {
Int_t nb=nBin+nearbyBins[ib];
if( !prPlotSeg2->testBin(nb) ) continue;
UChar_t ampN = prPlotSeg2->weights[nb];
stacksArr->push(nb,ampN);
if(ampN > amp) {
if(ib<7) stacksArr->push(nBin,amp);
amp = ampN;
if(isClstrInited) reinitCluster(nb,amp);
else initCluster(nb,amp);
break;
}
if(ampN==amp && isClstrInited) addBinInCluster(nb,amp);
}
Int_t ampP = stacksArr->pop(nBin);
if(ampP == amp) continue;
if(isClstrInited && !fillClusterSeg2()) return kFALSE;
if(ampP<0) break;
amp = ampP;
}
return kTRUE;
}
Bool_t HMdc34ClFinderSec::increaseClusterNum(void) {
if(nClsArr > 0) pClustersArrs->testCluster(0,nClsArr);
nClsArr++;
isClstrInited = kFALSE;
if(nClsArr < clusArrSize) return kTRUE;
nClsArr = clusArrSize;
Warning("increaseClusterNum"," Num.of clusters in sector %i > max\n",sector);
prPlotSeg2->clearBitArr();
return kFALSE;
}
void HMdc34ClFinderSec::addBinInCluster(Int_t nBin,UChar_t wt) {
Int_t nx = prPlotSeg2->xbin(nBin);
Int_t ny = prPlotSeg2->ybin(nBin);
cluster->addBin(nx, ny, wt);
if(ny<nLMinCl) {
nLMinCl = ny;
xMinClLines[ny] = xMaxClLines[ny] = nx;
} else if(ny>nLMaxCl) {
nLMaxCl = ny;
xMinClLines[ny] = xMaxClLines[ny] = nx;
}
else if(nx < xMinClLines[ny]) xMinClLines[ny] = nx;
else if(nx > xMaxClLines[ny]) xMaxClLines[ny] = nx;
}
Bool_t HMdc34ClFinderSec::fillClusterSeg2(void) {
cluster->calcXY();
Int_t xMinCl = MinElement(nLMaxCl-nLMinCl+1,xMinClLines+nLMinCl);
Int_t xMaxCl = MaxElement(nLMaxCl-nLMinCl+1,xMaxClLines+nLMinCl);
for(Int_t module=2; module<4; module++) {
if(!lMods[module]) continue;
HMdc34ClFinderMod& fMod=(*this)[module];
Int_t startLay = (module-2)*6;
for(Int_t lay=0; lay<6; lay++) {
Int_t layerSeg = startLay+lay;
HMdc34ClFinderLayer& fLay = fMod[lay];
HMdcLayListCells* cells = fLay.cells;
for(Int_t cell=cells->getFirstCell();cell>=0;cell=cells->next(cell)) {
if( fLay.yBin[xMinCl]+fLay.cellShUp[cell] < nLMinCl &&
fLay.yBin[xMaxCl]+fLay.cellShUp[cell] < nLMinCl ) continue;
if( fLay.yBin[xMinCl]+fLay.cellShDown[cell] > nLMaxCl &&
fLay.yBin[xMaxCl]+fLay.cellShDown[cell] > nLMaxCl ) break;
for(Int_t ny=nLMinCl; ny<=nLMaxCl; ny++) {
Int_t tYmin = fLay.yBin[xMinClLines[ny]]-ny;
Int_t tYmax = fLay.yBin[xMaxClLines[ny]]-ny;
Int_t dY1b = tYmin+fLay.cellShUp[cell];
Int_t dY2b = tYmax+fLay.cellShUp[cell];
if(dY1b<0 && dY2b<0) break;
if(dY1b) dY1b = dY1b>0 ? 1 : -1;
if(dY2b) dY2b = dY2b>0 ? 1 : -1;
if(dY1b*dY2b>0) {
Int_t dY1a = tYmin+fLay.cellShDown[cell];
Int_t dY2a = tYmax+fLay.cellShDown[cell];
if(dY1a) dY1a = dY1a>0 ? 1 : -1;
if(dY2a) dY2a = dY2a>0 ? 1 : -1;
if(dY1a*dY2a>0 && dY1a*dY1b > 0) continue;
}
Int_t nDeleted =
cluster->getLCells2().setTime(layerSeg,cell,cells->getTime(cell));
if(nDeleted && !HMdc34ClFinder::getQuietMode()) Warning(
"fillClusterSeg2","Too big cluster, %i cells was removed",nDeleted);
break;
}
}
}
}
return increaseClusterNum();
}
HMdc34ClFinder* HMdc34ClFinder::fMdc34ClFinder = 0;
Bool_t HMdc34ClFinder::quietMode = kTRUE;
Bool_t HMdc34ClFinder::useFixedLevel = kFALSE;
HMdc34ClFinder::HMdc34ClFinder(const Char_t* name,const Char_t* title,
const Char_t* context)
: HParSet(name,title,context) {
strcpy(detName,"Mdc");
array = new TObjArray(6);
fGetCont=HMdcGetContainers::getObject();
if( !fGetCont ) return;
fMdcDet = fGetCont->getMdcDetector();
fSpecGeomPar = fGetCont->getSpecGeomPar();
fSizesCells = HMdcSizesCells::getObject();
fKickPlane = fGetCont->getKickPlane();
fMdcGeomPar = fGetCont->getMdcGeomPar();
fMdcClusCat = fGetCont->getCatMdcClus(kTRUE);
xMinClLines = 0;
xMaxClLines = 0;
stack = 0;
stacksArr = 0 ;
mergeFlag = kTRUE;
nBinX = 320;
nBinY = 886;
}
Bool_t HMdc34ClFinder::initContainer(HMdcEvntListCells& event) {
if( !fMdcDet || !fSizesCells->initContainer() ) return kFALSE;
if( !HMdcGetContainers::isInited(fSpecGeomPar) ) return kFALSE;
if( !HMdcGetContainers::isInited(fKickPlane) ) return kFALSE;
if( !HMdcGetContainers::isInited(fMdcGeomPar) ) return kFALSE;
if( !status && (fSizesCells->hasChanged() || fSpecGeomPar->hasChanged() ||
fKickPlane->hasChanged()) || fMdcGeomPar->hasChanged()) {
changed=kTRUE;
if(!fMdcClusCat) return kFALSE;
Int_t *mods=fMdcDet->getModules();
if(!quietMode) printf("Project plot size of outer segment: %ix%i\n",
nBinX,nBinY);
for (Int_t sec = 0; sec < 6; sec++) {
if(!fMdcDet->isSectorActive(sec)) continue;
if( (*array)[sec] == 0) {
Int_t nLM2 = HMdcTrackDSet::getTrFnNLayersInMod()[sec*4+2];
Int_t nLM3 = HMdcTrackDSet::getTrFnNLayersInMod()[sec*4+3];
if((mods[sec*4+2]==0||nLM2==0) && (mods[sec*4+3]==0||nLM3==0)) continue;
HMdc34ClFinderSec* secF = new HMdc34ClFinderSec(sec,nBinX,nBinY);
(*array)[sec] = secF;
if(!mergeFlag) secF->donotMergeClusters();
secF->setClusCut(fMdcClusCat);
secF->setKickPlane(fKickPlane);
secF->setCellsList(event[sec]);
if(xMinClLines == 0) xMinClLines = new Short_t [nBinY];
if(xMaxClLines == 0) xMaxClLines = new Short_t [nBinY];
secF->setXMinClLines(xMinClLines);
secF->setXMaxClLines(xMaxClLines);
if(stacksArr==0) {
stacksArr = new HMdcClFnStacksArr(10000);
stack = stacksArr->getOneStack();
}
secF->setClFnStArr(stacksArr);
secF->setClFnStack(stack);
}
if(!calcTarget(sec)) return kFALSE;
if(!calcProjPlaneSeg2(sec)) return kFALSE;
if(!calcSizePlotSeg2(sec)) return kFALSE;
if(!calcWiresProj(sec)) return kFALSE;
#if DEBUG_LEVEL>2
(*this)[sec].prPlotSeg2->print();
#endif
}
if(versions[1]<0 || versions[2]<0) versions[1]=versions[2]=0;
else versions[2]++;
} else changed=kFALSE;
return kTRUE;
}
HMdc34ClFinder::~HMdc34ClFinder() {
if(array) {
array->Delete();
delete array;
}
fMdc34ClFinder=0;
if(xMinClLines) delete [] xMinClLines;
if(xMaxClLines) delete [] xMaxClLines;
xMinClLines = 0;
xMaxClLines = 0;
if(stacksArr) delete stacksArr;
stacksArr = 0;
}
HMdc34ClFinder* HMdc34ClFinder::getObject(void) {
if(!fMdc34ClFinder) fMdc34ClFinder=new HMdc34ClFinder();
return fMdc34ClFinder;
}
HMdc34ClFinder* HMdc34ClFinder::getExObject(void) {
return fMdc34ClFinder;
}
void HMdc34ClFinder::deleteCont(void) {
if(fMdc34ClFinder) delete fMdc34ClFinder;
}
void HMdc34ClFinder::clear(void) {
for(Int_t s=0;s<6;s++) if((*array)[s]) (*this)[s].clear();
}
Bool_t HMdc34ClFinder::calcTarget(Int_t sec){
if(!fSizesCells->hasChanged() && !fSpecGeomPar->hasChanged()) return kTRUE;
Int_t nT = fSpecGeomPar->getNumTargets()-1;
if( nT < 0 ) {
Error("calcTarget","Number of targets = %i!",nT+1);
return kFALSE;
}
HMdc34ClFinderSec& fsec=(*this)[sec];
HGeomVector* target = fsec.getTargetArr();
target[0] = (fSpecGeomPar->getTarget(0)->getTransform()).getTransVector();
target[0].setZ(target[0].getZ() +
fSpecGeomPar->getTarget(0)->getPoint(0)->getZ());
target[1] = (fSpecGeomPar->getTarget(nT)->getTransform()).getTransVector();
target[1].setZ(target[1].getZ() +
fSpecGeomPar->getTarget(nT)->getPoint(2)->getZ());
const HGeomTransform* trans=(*fSizesCells)[sec].getLabTrans();
if(&trans == 0) return kFALSE;
target[0] = trans->transTo(target[0]);
target[1] = trans->transTo(target[1]);
return kTRUE;
}
Bool_t HMdc34ClFinder::calcProjPlaneSeg2(Int_t sec){
if(!fSizesCells->hasChanged()) return kTRUE;
HMdc34ClFinderSec& fsec=(*this)[sec];
HMdcSizesCellsSec& fSCellsSec=(*fSizesCells)[sec];
if( !&fSCellsSec ) return kFALSE;
Double_t par=10.;
UChar_t* mSeg = fsec.getMSeg();
Int_t nL=(Int_t)par;
if(nL<-7||nL>7 || (nL<0 && mSeg[1]==2) || (nL>0 && mSeg[1]==1) ||
(nL==0 && mSeg[1]!=3)) {
nL = 0;
par = 0.1;
}
const Char_t *text="as projection plane";
HGeomTransform prPl;
if( (nL==0 && mSeg[1]<3) || nL==-7 || nL==7 ) {
Int_t mod = nL==0 ? mSeg[1]+1 : (nL+7)/14+2;
prPl=*(fSCellsSec[mod].getSecTrans());
if(!quietMode) printf("\n===> Sec.%i Seg.2: Using middle plane of MDC%i %s\n",
sec+1,mod+1,text);
} else if( nL==0 && mSeg[1]==3 ) {
const HGeomTransform* trLayer6=fSCellsSec[2][5].getSecTrans();
const HGeomTransform* trLayer1=fSCellsSec[3][0].getSecTrans();
Double_t distToL6=trLayer6->getTransVector().length();
Double_t distToL1=trLayer1->getTransVector().length();
Double_t newDist=distToL6+(distToL1-distToL6)*(1.+par)/2.;
Double_t mult=1.;
if(par<=0) {
prPl.setTransform(*trLayer6);
mult=newDist/distToL6;
} else {
prPl.setTransform(*trLayer1);
mult=newDist/distToL1;
}
HGeomVector prTr(prPl.getTransVector());
prTr*=mult;
prPl.setTransVector(prTr);
if(!quietMode) printf("\n===> Sec.%i Seg.2: Using plane between MDC 3 & 4 (p=%g) %s\n",
sec+1,par,text);
} else {
Int_t mod = nL<0 ? 2:3;
if(nL<0) nL+=7;
prPl=*(fSCellsSec[mod][nL-1].getSecTrans());
if(!quietMode) printf("\n===> Sec.%i Seg.2: Using MDC%i, layer %i %s\n",
sec+1,mod+1,nL,text);
}
fsec.getPlotSeg2()->setPlanePar(prPl);
return kTRUE;
}
Bool_t HMdc34ClFinder::calcSizePlotSeg2(Int_t sec){
if( !fKickPlane->hasChanged() && !fMdcGeomPar->hasChanged() &&
!fSizesCells->hasChanged()) return kTRUE;
HMdc34ClFinderSec& fsec=(*this)[sec];
UChar_t* mSeg = fsec.getMSeg();
if(mSeg[0] == 0) return kFALSE;
HMdcSizesCellsSec& fSCellsSec=(*fSizesCells)[sec];
if( !&fSCellsSec ) return kFALSE;
HGeomVector* target = fsec.getTargetArr();
HGeomVector vect[2][4];
HGeomVector pKick[4];
const HGeomTransform* trans=0;
HGeomVector newP[4][2];
HMdcPlane plane;
for(Int_t mod=0; mod<2; mod++) {
if(fsec.mdcFlag(mod) == 0) continue;
HGeomCompositeVolume *fVolMdc=fGetCont->getGeomCompositeVolume(mod);
if(!fVolMdc) return kFALSE;
HGeomVolume* fVolLayer=fVolMdc->getComponent(5);
if(!fVolLayer) return kFALSE;
trans=fSCellsSec[mod].getSecTrans();
plane.setPlanePar(*trans);
for(Int_t point=0; point<4; point++) {
HGeomVector *fpoint=fVolLayer->getPoint(point);
if(!fpoint) return kFALSE;
Int_t indx = mSeg[0]==3 ? mod : 0;
vect[indx][point] = *fpoint;
vect[indx][point].setZ(0.);
vect[indx][point] = trans->transFrom(vect[indx][point]);
if(mod==1 && mSeg[0]==3) {
Int_t nTag=(point==0 || point==3) ? 0:1;
vect[0][point] -= target[nTag];
plane.calcIntersection(target[nTag],vect[0][point],vect[0][point]);
}
}
}
if(mSeg[0]==3) {
for(Int_t mod=0; mod<2; mod++) {
for(Int_t point=0; point<4; point++)
vect[mod][point]=trans->transTo(vect[mod][point]);
}
newP[0][0] = vect[0][0](0)>vect[1][0](0) ? vect[0][0] : vect[1][0];
newP[0][1] = vect[0][1](0)>vect[1][1](0) ? vect[0][1] : vect[1][1];
newP[1][0] = vect[0][1](1)>vect[1][1](1) ? vect[0][1] : vect[1][1];
newP[1][1] = vect[0][2](1)>vect[1][2](1) ? vect[0][2] : vect[1][2];
newP[2][0] = vect[0][2](0)<vect[1][2](0) ? vect[0][2] : vect[1][2];
newP[2][1] = vect[0][3](0)<vect[1][3](0) ? vect[0][3] : vect[1][3];
newP[3][0] = vect[0][3](1)<vect[1][3](1) ? vect[0][3] : vect[1][3];
newP[3][1] = vect[0][0](1)<vect[1][0](1) ? vect[0][0] : vect[1][0];
calcCrossLines(newP[0][0],newP[0][1],newP[3][0],newP[3][1],vect[0][0]);
calcCrossLines(newP[0][0],newP[0][1],newP[1][0],newP[1][1],vect[0][1]);
calcCrossLines(newP[1][0],newP[1][1],newP[2][0],newP[2][1],vect[0][2]);
calcCrossLines(newP[2][0],newP[2][1],newP[3][0],newP[3][1],vect[0][3]);
for(Int_t point=0; point<4; point++)
vect[0][point] = trans->transFrom(vect[0][point]);
}
HGeomVector dir;
for(Int_t point=0; point<4; point++) {
if(point==0 || point==3) dir = vect[0][point]-target[0];
else dir = vect[0][point]-target[1];
if( !(fKickPlane->calcIntersection(vect[0][point],dir,pKick[point]))) {
Error("calcSizePlotSeg2",
"Sec.%i: The point %i on the kick plane can't be calculated.",
sec+1,point);
return kFALSE;
}
}
Double_t yMin = 1.e+10;
Double_t yMax = -1.e+10;
for(Int_t mod=2; mod<4; mod++) {
if(fsec.mdcFlag(mod) == 0) continue;
HGeomCompositeVolume *fVolMdc=fGetCont->getGeomCompositeVolume(mod);
Int_t indx = mSeg[1]==3 ? mod-2 : 0;
trans=fSCellsSec[mod].getSecTrans();
for(Int_t point=0; point<4; point++) {
vect[indx][point] = *(fVolMdc->getPoint(point));
vect[indx][point].setZ(0.);
vect[indx][point] = trans->transFrom(vect[indx][point]);
vect[indx][point] -= pKick[point];
fsec.getPlotSeg2()->calcIntersection(pKick[point],vect[indx][point],
vect[indx][point]);
if(vect[indx][point](1)<yMin) yMin = vect[indx][point](1);
if(vect[indx][point](1)>yMax) yMax = vect[indx][point](1);
}
}
Double_t xMaxD=xLine(vect[0][0],vect[0][1],yMin);
Double_t xMinD=xLine(vect[0][2],vect[0][3],yMin);
Double_t xMax=xLine(vect[0][0],vect[0][1],yMax);
Double_t xMin=xLine(vect[0][2],vect[0][3],yMax);
if(mSeg[1]==3) {
Double_t xnew=xLine(vect[1][0],vect[1][1],yMin);
if(xnew>xMaxD) xMaxD=xnew;
xnew=xLine(vect[1][2],vect[1][3],yMin);
if(xnew<xMinD) xMinD=xnew;
xnew=xLine(vect[1][0],vect[1][1],yMax);
if(xnew>xMax) xMax=xnew;
xnew=xLine(vect[1][2],vect[1][3],yMax);
if(xnew<xMin) xMin=xnew;
}
fsec.getPlotSeg2()->setEdges(yMin, xMinD, xMaxD, yMax, xMin, xMax);
return kTRUE;
}
void HMdc34ClFinder::calcCrossLines(HGeomVector& p1l1, HGeomVector& p2l1,
HGeomVector& p1l2, HGeomVector& p2l2, HGeomVector& cross) {
Double_t a1 = (p1l1(1)-p2l1(1))/(p1l1(0)-p2l1(0));
Double_t b1 = p1l1(1)-a1*p1l1(0);
Double_t a2 = (p1l2(1)-p2l2(1))/(p1l2(0)-p2l2(0));
Double_t b2 = p1l2(1)-a2*p1l2(0);
Double_t x = (b2-b1)/(a1-a2);
Double_t y = a2*x+b2;
cross.setXYZ(x,y,0.);
}
Double_t HMdc34ClFinder::xLine(HGeomVector& p1,HGeomVector& p2,Double_t yi) {
return (yi-p2(1))/(p1(1)-p2(1))*(p1(0)-p2(0))+p2(0);
}
Bool_t HMdc34ClFinder::calcWiresProj(Int_t sec) {
if( !fKickPlane->hasChanged() && !fMdcGeomPar->hasChanged() &&
!fSizesCells->hasChanged()) return kTRUE;
HMdc34ClFinderSec& fsec=(*this)[sec];
UChar_t* mSeg = fsec.getMSeg();
if(mSeg[0]==0) return kFALSE;
HMdcSizesCellsSec& fSCellsSec=(*fSizesCells)[sec];
if( !&fSCellsSec ) return kFALSE;
HMdcProjPlot* prPlotSeg2 = fsec.getPlotSeg2();
HGeomVector wire[2];
HGeomVector rib[4];
HGeomVector point[2];
Int_t mod = mSeg[0]>1 ? 1:0;
const HGeomTransform* trans=fSCellsSec[mod].getSecTrans();
if(!trans) return kFALSE;
HGeomVector* target = fsec.getTargetArr();
HGeomVector tag = target[0]+target[1];
tag /= 2.;
HGeomVector dir = trans->getTransVector()-tag;
dir /= dir.length();
HGeomVector pKick;
if(!(fKickPlane->calcIntersection(tag,dir,pKick))) {
Error("calcWiresProj",
"Sec.%i: The mid.point on the kick plane can't be calculated.",sec+1);
return kFALSE;
}
for(Int_t mod=2; mod<4; mod++) {
if(fsec.mdcFlag(mod) == 0) continue;
HMdcSizesCellsMod& fSCellsMod=fSCellsSec[mod];
if(!&fSCellsMod) return kFALSE;
HMdc34ClFinderMod& fMod=fsec[mod];
for(Int_t lay=0; lay<6; lay++) {
HMdcSizesCellsLayer& fSCellsLay=fSCellsMod[lay];
const HGeomTransform* secsys = fSCellsLay.getSecTrans();
Double_t catDist = fSCellsLay.getHalfCatDist();
Double_t pitch = fSCellsLay.getPitch();
Double_t cosWiOr = fSCellsLay.getCosWireOr();
Double_t wOfSet = fSCellsLay.getWireOffSet();
Double_t aWire = fSCellsLay.getTanWireOr();
Double_t bWire = (100.*pitch-wOfSet)/cosWiOr;
wire[0].setXYZ( 500., aWire*500.+bWire,0.);
wire[1].setXYZ(-500.,-aWire*500.+bWire,0.);
for(Int_t np=0; np<2; np++) {
wire[np] = secsys->transFrom(wire[np]);
wire[np] -= pKick;
prPlotSeg2->calcIntersection(pKick,wire[np],point[np]);
}
#if DEBUG_LEVEL>2
printf("M%i L%i C%i Size of wire on the proj.plane (poin1&point2):\n",
mod+1,lay+1,100+1);
point[0].print();
point[1].print();
#endif
Double_t al = (point[0](1)-point[1](1))/(point[0](0)-point[1](0));
Double_t x = (prPlotSeg2->nBinX-1) * prPlotSeg2->stX +
prPlotSeg2->xMin;
Double_t yLeft = al*(x - point[1](0))+point[1](1);
Double_t yRight = al*(prPlotSeg2->xMin - point[1](0))+point[1](1);
Short_t nYLeft = Short_t((yLeft - prPlotSeg2->yMin) /
prPlotSeg2->stY);
Short_t nYRight = Short_t((yRight - prPlotSeg2->yMin) /
prPlotSeg2->stY);
Short_t yShift = Min(nYLeft,nYRight);
if(Abs(nYLeft-nYRight) >= prPlotSeg2->nBinY) {
Error("calcWiresProj",
"S%i M%i L%i The region of Y bins of proj. wire >= nBinY(%i).",
sec+1,mod+1,lay+1,prPlotSeg2->nBinY);
return kFALSE;
}
HMdc34ClFinderLayer& fLayer=fMod[lay];
if( !fLayer.createArrayBins(prPlotSeg2->nBinX) ) return kFALSE;
for(Int_t nx=0; nx<fLayer.nBinX; nx++) {
Double_t x = nx * prPlotSeg2->stX + prPlotSeg2->xMin;
Double_t y = al*(x-point[1](0))+point[1](1);
Short_t nY = Short_t((y-prPlotSeg2->yMin) / prPlotSeg2->stY);
fLayer.yBin[nx]=nY-yShift;
}
Int_t nLines = fLayer.yBin[fLayer.nBinX - 1]-fLayer.yBin[0];
fLayer.nYLines = Abs(nLines)+1;
fLayer.wOrType = nLines>=0 ? 1 : -1;
fLayer.yFirst = nLines>=0 ? fLayer.yBin[0]:fLayer.yBin[fLayer.nBinX-1];
if(fLayer.xBin1) {
delete [] fLayer.xBin1;
delete [] fLayer.xBin2;
}
fLayer.xBin1 = new Short_t [fLayer.nYLines];
fLayer.xBin2 = new Short_t [fLayer.nYLines];
for(Int_t yb=0;yb<fLayer.nYLines;yb++) {
Int_t yl=yb+fLayer.yFirst;
fLayer.xBin1[yb] = -1;
fLayer.xBin2[yb] = -1;
for(Int_t xb=0;xb<fLayer.nBinX;xb++) {
Int_t dy=(fLayer.wOrType>=0) ? fLayer.yBin[xb]-yl:yl-fLayer.yBin[xb];
if(dy<0) continue;
if(dy>0 && fLayer.xBin2[yb]>=0) break;
if(fLayer.xBin1[yb]<0) fLayer.xBin1[yb]=xb;
fLayer.xBin2[yb] = dy==0 ? xb : xb-1;
}
}
for(Int_t cell=0; cell<101; cell+=100) {
Double_t bCell= (cell*pitch-wOfSet-pitch/2.)/cosWiOr;
for(Int_t n=0; n<4; n++) {
Double_t x = ((n&1) == 0) ? 500.: -500.;
Double_t z = (n < 2) ? -catDist : catDist;
rib[n].setXYZ(x, aWire*x+bCell, z);
rib[n]=secsys->transFrom(rib[n]);
}
Double_t dX = rib[0](0)-rib[1](0);
rib[1].setXYZ(0.,(rib[0](0)*rib[1](1)-rib[1](0)*rib[0](1))/dX,
(rib[0](0)*rib[1](2)-rib[1](0)*rib[0](2))/dX);
dX = rib[2](0)-rib[3](0);
rib[3].setXYZ(0.,(rib[2](0)*rib[3](1)-rib[3](0)*rib[2](1))/dX,
(rib[2](0)*rib[3](2)-rib[3](0)*rib[2](2))/dX);
if(cell==0) {
fLayer.y0[0] = rib[1](1);
fLayer.z0[0] = rib[1](2);
fLayer.y0[1] = rib[3](1);
fLayer.z0[1] = rib[3](2);
fLayer.tgY = (rib[0](1)-rib[1](1))/(rib[0](0)-rib[1](0));
fLayer.tgZ = (rib[0](2)-rib[1](2))/(rib[0](0)-rib[1](0));
} else {
fLayer.yStep = (rib[1](1)-fLayer.y0[0])/100.;
fLayer.zStep = (rib[1](2)-fLayer.z0[0])/100.;
if(rib[1](2)+prPlotSeg2->B()*rib[1](1)-prPlotSeg2->D()<0) {
fLayer.nPSegOnKick[0][0]=(fLayer.tgY>0) ? 2 : 3;
fLayer.nPSegOnKick[1][0]=(fLayer.tgY>0) ? 1 : 0;
} else {
fLayer.nPSegOnKick[0][0]=(fLayer.tgY>0) ? 1 : 0;
fLayer.nPSegOnKick[1][0]=(fLayer.tgY>0) ? 2 : 3;
}
if(rib[3](2)+prPlotSeg2->B()*rib[3](1)-prPlotSeg2->D()<0) {
fLayer.nPSegOnKick[0][1]=(fLayer.tgY>0) ? 2 : 3;
fLayer.nPSegOnKick[1][1]=(fLayer.tgY>0) ? 1 : 0;
} else {
fLayer.nPSegOnKick[0][1]=(fLayer.tgY>0) ? 1 : 0;
fLayer.nPSegOnKick[1][1]=(fLayer.tgY>0) ? 2 : 3;
}
}
}
}
}
return kTRUE;
}
void HMdc34ClFinder::setPrPlotSize(Int_t nx,Int_t ny) {
nBinX = nx;
nBinY = ny<=1 ? (nx*277)/100 : ny;
}
void HMdc34ClFinder::setCellsList(HMdcEvntListCells& event) {
for(Int_t s=0;s<6;s++) if((*array)[s]) (*this)[s].setCellsList(event[s]);
}
void HMdc34ClFinderSec::setCellsList(HMdcSecListCells& event) {
notEnoughWrs = kTRUE;
for(Int_t m=2;m<4;m++) {
if((*array)[m]) (*this)[m].setCellsList(event[m]);
maxAmp[m] = -1;
}
}
void HMdc34ClFinderMod::setCellsList(HMdcModListCells& event) {
for(Int_t l=0;l<6;l++) if((*array)[l]) (*this)[l].setCellsList(event[l]);
}
Last change: Sat May 22 12:59:33 2010
Last generated: 2010-05-22 12:59
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.