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 "hmdcgeomobj.h"
#include "hmdcgetcontainers.h"
#include "hmdccal1sim.h"
#include "hmdcclussim.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 <stdlib.h>
//*-- AUTHOR : Pechenov Vladimir
//*-- Modified : 07/02/2003 by V. Pechenov
//*-- Modified : 05/06/2002 by V. Pechenov
//*-- Modified : 10/05/2001 by V. Pechenov
//*-- Modified : 07/03/2001 by V. Pechenov
//*-- Modified : 12/07/2000 by V.Pechenov
//_HADES_CLASS_DESCRIPTION
////////////////////////////////////////////////////////////////
// HMdc34ClFinder
//
// Track finder for outer segment MDC:
// Constructor:
//
////////////////////////////////////////////////////////////////
ClassImp(HMdcProjPlot)
ClassImp(HMdc34ClFinderLayer)
ClassImp(HMdc34ClFinderMod)
ClassImp(HMdc34ClFinderSec)
ClassImp(HMdc34ClFinder)
Char_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;
memset(plotBArSc,0,sizeBAr);
rootPlot=0;
rootPlotV2=0;
if(size>wtArrSize) {
if(weights) delete [] weights;
weights=new Char_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);
}
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=yBinsPos=0;
}
void HMdcProjPlot::setEdges(Float_t iyMin, Float_t ixMinD, Float_t ixMaxD,
Float_t iyMax, Float_t ixMin, Float_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;
Float_t aL23=(xMinD-xMin)/(yMin-yMax);
Float_t aL01=(xMaxD-xMax)/(yMin-yMax);
for(Int_t ny=0; ny<nBinY; ny++) {
Float_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]) {
if(ny==0) yMinL[nx]=1;
else yMinL[nx]=ny;
} else break;
}
}
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=%gn",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(Char_t* name, 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(Char_t* name, 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;
}
//----------Layer-----------------------------
HMdc34ClFinderLayer::HMdc34ClFinderLayer(Int_t sec, Int_t mod, Int_t layer) {
// Geting of pointers to par. cont.
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;
}
HMdc34ClFinderLayer::~HMdc34ClFinderLayer() {
if(yBin) {
delete [] yBin;
yBin=0;
}
if(xBin1) {
delete [] xBin1;
xBin1=0;
delete [] xBin2;
xBin2=0;
}
if(cellShUp) {
delete [] cellShUp;
delete [] cellShDown;
cellShUp=cellShDown=0;
}
yBin=0;
nBinX=0;
}
Bool_t HMdc34ClFinderLayer::createArrayBins(Short_t nBins) {
if( nBinX != nBins ) {
if(yBin) delete [] yBin;
nBinX=nBins;
yBin=new Short_t [nBins];
if(yBin==0) return kFALSE;
}
return kTRUE;
}
//------------Module----------------------------
HMdc34ClFinderMod::HMdc34ClFinderMod(Int_t sec, Int_t mod) {
// constructor creates an array of pointers of type HMdc34ClFinderLayer
//
array = new TObjArray(6);
for (Int_t layer = 0; layer < 6; layer++) {
(*array)[layer] = new HMdc34ClFinderLayer(sec, mod, layer);
}
}
HMdc34ClFinderMod::~HMdc34ClFinderMod() {
// destructor
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].getNCells();
return nHits;
}
void HMdc34ClFinderMod::clear(void) {
for(Int_t lay=0; lay<6; lay++) if((*array)[lay]) (*this)[lay].clear();
}
//----------Sector------------------------------------
HMdcClFnStack* HMdc34ClFinderSec::stack=0;
Short_t HMdc34ClFinderSec::xMinClLines[300];
Short_t HMdc34ClFinderSec::xMaxClLines[300];
HKickPlane2* HMdc34ClFinderSec::fkick=0;
HMdcList12GroupCells HMdc34ClFinderSec::cLCells2[500];
Bool_t HMdc34ClFinderSec::cStatus[500];
Short_t HMdc34ClFinderSec::cNMergedClus[500];
Double_t HMdc34ClFinderSec::cSumWt[500];
Int_t HMdc34ClFinderSec::cNBins[500];
Double_t HMdc34ClFinderSec::cX[500];
Double_t HMdc34ClFinderSec::cY[500];
Double_t HMdc34ClFinderSec::cMeanX[500];
Double_t HMdc34ClFinderSec::cMeanY[500];
Double_t HMdc34ClFinderSec::cMeanXX[500];
Double_t HMdc34ClFinderSec::cMeanYY[500];
Double_t HMdc34ClFinderSec::cMeanYX[500];
Double_t HMdc34ClFinderSec::cMeanXWt[500];
Double_t HMdc34ClFinderSec::cMeanYWt[500];
Double_t HMdc34ClFinderSec::cMeanXXWt[500];
Double_t HMdc34ClFinderSec::cMeanYYWt[500];
HMdc34ClFinderSec::HMdc34ClFinderSec(Int_t sec, Int_t inBinX, Int_t inBinY) {
// constructor creates an array of pointers of type HMdc34ClFinderMod
//
sector=sec;
HMdcDetector* fMdcDet=(HMdcDetector*)gHades->getSetup()->getDetector("Mdc");
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]) {
prPlotSeg2=new HMdcProjPlot(mSeg[1],inBinX,inBinY);
nBinX=prPlotSeg2->nBinX;
if(stack==0) stack=new HMdcClFnStack(1000);
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(3,sector,0,0);
radToDeg=180./acos(Double_t(-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];
}
HMdc34ClFinderSec::~HMdc34ClFinderSec() {
// destructor
if(array) {
array->Delete();
delete array;
}
if(prPlotSeg2) {
delete prPlotSeg2;
prPlotSeg2=0;
}
if(stack) delete stack;
stack=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) {
counter=0;
for(Int_t mod=2; mod<4; mod++) {
if((*array)[mod]) (*this)[mod].clear();
maxAmp[mod]=-1;
}
notEnoughWrs=kTRUE;
}
void HMdc34ClFinderSec::setMinBin(Int_t *mBin) {
for(int i=0;i<4;i++) minAmp[i]=mBin[i];
notEnoughWrs=testMaxAmp();
}
Bool_t HMdc34ClFinderSec::testMaxAmp(void) {
// Calculating checking max.amplitude
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.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* fClus,
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);
// Errors ????????????????????????????????????????????????
Double_t cosPhiS=cos(phi+TMath::Pi()*0.5);
Double_t sinPhiS=sin(phi+TMath::Pi()*0.5);
Double_t x0 = r0*cosPhiS;
Double_t y0 = r0*sinPhiS;
// let length=1.
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(TMath::Power(cosPhiS*errR0,2)+TMath::Power(y0*errPh,2));
Double_t errY0=sqrt(TMath::Power(sinPhiS*errR0,2)+TMath::Power(x0*errPh,2));
//errY0=errY0/10.; //!!!!!!!!!!!!!!!!!!!!!!!! ??????????
Double_t errDirX=sqrt(TMath::Power(sinPhi*sinTh*errPh,2)+
TMath::Power(cosPhi*cosTh*errTh,2));
Double_t errDirY=sqrt(TMath::Power(cosPhi*sinTh*errPh,2)+
TMath::Power(sinPhi*cosTh*errTh,2));
//errDirY=errDirY/10.; //!!!!!!!!!!!!!!!!!!!!!!!! ??????????
Double_t errDirZ=TMath::Abs(sinTh*errTh);
// segOnkick[4]:
// [0]---[1] ^ y
// | | |
// | | |
// [2]---[3] |
// X<-------------
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); // hit
dirSeg1[4].setXYZ(dirX,dirY,dirZ); // hit
indexPar=fClusCat->getIndex(fClus);
indexFCh=-1;
indexLCh=-2;
Int_t nClustCh=findClustersSeg2();
if(nClustCh>0) fClus->setIndexChilds(indexFCh,indexLCh);
return nClustCh;
}
Int_t HMdc34ClFinderSec::findClustersSeg2(HMdcClus* fClus, Int_t *mBin){
if(fClus->getIOSeg() != 0) return 0;
if(mBin!=0) setMinBin(mBin);
if(notEnoughWrs) return 0;
Double_t dirX=fClus->getX();
Double_t dirY=fClus->getY();
Double_t dirZ=fClus->getZ();
Double_t errDirX=fClus->getErrX();
Double_t errDirY=fClus->getErrY();
// Double_t errDirZ=fClus->getZ();
Double_t x0=fClus->getXTarg();
Double_t y0=fClus->getYTarg();
Double_t z0=fClus->getZTarg();
Double_t errX0=fClus->getErrXTarg();
Double_t errY0=fClus->getErrYTarg();
Double_t errZ0=fClus->getErrZTarg();
// segOnkick[4]:
// [0]---[1] ^ y
// | | |
// | | |
// [2]---[3] |
// X<-------------
seg1[0].setXYZ(x0-errX0,y0-errY0,z0-errZ0);
dirSeg1[0].setXYZ(dirX+errDirX,dirY+errDirY,
fClus->getZOnPrPlane(dirX+errDirX,dirY+errDirY));
seg1[1].setXYZ(x0+errX0,y0-errY0,z0-errZ0);
dirSeg1[1].setXYZ(dirX-errDirX,dirY+errDirY,
fClus->getZOnPrPlane(dirX-errDirX,dirY+errDirY));
seg1[2].setXYZ(x0-errX0,y0+errY0,z0-errZ0);
dirSeg1[2].setXYZ(dirX+errDirX,dirY-errDirY,
fClus->getZOnPrPlane(dirX+errDirX,dirY-errDirY));
seg1[3].setXYZ(x0+errX0,y0+errY0,z0-errZ0);
dirSeg1[3].setXYZ(dirX-errDirX,dirY-errDirY,
fClus->getZOnPrPlane(dirX-errDirX,dirY-errDirY));
seg1[4].setXYZ(x0,y0,z0); // hit
dirSeg1[4].setXYZ(dirX,dirY,dirZ); // hit
for(Int_t np=0; np<5; np++) dirSeg1[np]-=seg1[np];
indexPar=fClusCat->getIndex(fClus);
indexFCh=-1;
indexLCh=-2;
Int_t nClustCh=findClustersSeg2();
if(nClustCh>0) fClus->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(
(TMath::Abs(segOnKick[0](0)-segOnKick[1](0))+
TMath::Abs(segOnKick[2](0)-segOnKick[3](0)))*0.25,
(TMath::Abs(segOnKick[0](1)-segOnKick[2](1))+
TMath::Abs(segOnKick[1](1)-segOnKick[3](1)))*0.25,
(TMath::Abs(segOnKick[0](2)-segOnKick[2](2))+
TMath::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); // D-A*x
al[np]=tmp-segOnKick[np](2); // D-A*x-z
bl[np]=tmp*segOnKick[np](1); // (D-A*x)*y
cl[np]=segOnKick[np](2)+prPlotSeg2->B()*segOnKick[np](1); // z+B*y
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) {
// Filling proj.plot for one MDC
cPlMod=(mod==2) ? prPlotSeg2->plModF:prPlotSeg2->plModS;
if(cPlMod==0) return;
cFMod=&((*this)[mod]);
cXMinM=xCMin[mod];
cXMaxM=xCMax[mod];
Int_t maxB4Sc=0;
Int_t minB4Sc=prPlotSeg2->size;
Int_t minAm=minAmp[mod];
if(minAm<3) minAm=3;
Char_t* weights=prPlotSeg2->weights;
UChar_t* plotBArSc=prPlotSeg2->plotBArSc;
Int_t layList[6]={2,3,1,4,0,5}; // order of layers at the pr.plot filling
Int_t nFiredLay=0;
for(Int_t il=0; il<6; il++) {
Int_t lay=layList[il];
setLayerVar(lay);
if(!test) continue;
nFiredLay++;
if(maxAmp[mod]-nFiredLay+1>=minAm) { // determination min(max)Bin[y]
if(nFiredLay<minAm) makeLayProjV1(lay); // filling pr.plot
else { // ...+ amp.checking
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 { // filling in minBin[y]-maxBin[y] only
if(nFiredLay<minAm) makeLayProjV2(lay); // filling pr.plot
else { // ...+ amp.checking
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;
}
void HMdc34ClFinderSec::setLayerVar(Int_t lay) {
// Setting of layer var. for filling proj.plot
cFLay=&((*cFMod)[lay]);
cell=cFLay->getFirstCell();
if(cell<0) {
test=kFALSE;
return;
}
if(cell >= cFLay->nCells) {
Warning("makePlotSeg2","Cell %i >= num.cells(=%i)",cell,cFLay->nCells);
test=kFALSE;
return;
}
nBinYM2=prPlotSeg2->nBinY-2;
shUp=TMath::Max(calcYbin(0,0,cell+1),calcYbin(0,1,cell+1));
shDown=TMath::Min(calcYbin(1,0,cell),calcYbin(1,1,cell));
add=1<<lay;
cFLay->cellShUp[cell]=shUp;
cFLay->cellShDown[cell]=shDown;
nYLinesM1=cFLay->nYLines-1;
xBin1L=cFLay->xBin1;
xBin2L=cFLay->xBin2;
shUpN=0;
shDownN=0;
test=kTRUE;
}
Bool_t HMdc34ClFinderSec::setNextCell(void) {
cell=cFLay->next(cell);
if(cell>=0) {
shDownN=TMath::Min(calcYbin(1,0,cell),calcYbin(1,1,cell));
shUpN=TMath::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) {
// Determination and testing X region for pr.plot filling
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) {
// plot filling and filled region determination
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) {
// plot filling and filled region determination ofMDC4 in region determined
// in MDC3
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) {
// plot filling in filled regions only
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;
Char_t* weights=prPlotSeg2->weights;
UChar_t* plotBArSc=prPlotSeg2->plotBArSc;
UChar_t* plModF=prPlotSeg2->plModF;
Int_t layList[6]={2,3,1,4,0,5}; // order of layers at the pr.plot filling
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];
setLayerVar(lay);
if(!test) 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;
Char_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}; // order of layers at the pr.plot filling
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) ? kFALSE:kTRUE;
cPlMod=(!isMod4) ? prPlotSeg2->plModF:prPlotSeg2->plModS;
if(cPlMod==0) continue;
setLayerVar(lay);
if(!test) 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(Char_t* name, 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++) {
setLayerVar(lay);
if(!test) continue;
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) {
//upDo=0 - up rib, =1 down rib
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) {
Int_t nClusters=0;
nClsArr=0;
locClus[1]=1;
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;
while((n=HMdcBArray::nextAndUnset(b1,b2,n))>=0) {
calcClusterSeg2(nBin4+n);
nClsArr++;
if(nClsArr >= 500) break;
}
if(nClsArr >= 500) {
Warning("scanPlotSeg2"," Num. of clusters in sector %i > 500n",sector);
memset(prPlotSeg2->plotBArSc,0,prPlotSeg2->sizeBAr);
break;
}
}
Int_t nClus=nClsArr;
while(nClus>1) {
Bool_t nomerg=kTRUE;
for(Int_t cl1=0; cl1<nClsArr-1; cl1++) {
if(!cStatus[cl1]) continue;
HMdcList12GroupCells& cLCSeg2=cLCells2[cl1];
for(Int_t cl2=cl1+1; cl2<nClsArr; cl2++) {
if(!cStatus[cl2]) continue;
HMdcList12GroupCells* cLCSeg2s=&(cLCells2[cl2]);
Float_t dY=cY[cl1]-cY[cl2];
if(dY>100.) break;
if(fabs(dY) > 30.) continue; // 30. mm !???
if(fabs(cX[cl1]-cX[cl2]) > 100.) continue; // 100. mm !???
if(realTypeClFinder==2) { // realTypeClFinder==2 when typeClFinder==2 +...!
if(/*typeClFinder==2 &&*/ cLCSeg2.compare(cLCSeg2s,0,11)<6) continue;
} else {
Int_t lev=(minAmp[2]>4) ? 4:3; //???
if(minAmp[2]>0 && cLCSeg2.compare(cLCSeg2s,0, 5)<lev) continue;
lev=(minAmp[3]>4) ? 4:3; //???
if(minAmp[3]>0 && cLCSeg2.compare(cLCSeg2s,6,11)<lev) continue;
}
cLCSeg2.add(cLCSeg2s);
cStatus[cl2]=kFALSE;
cSumWt[cl1] += cSumWt[cl2];
cNBins[cl1] += cNBins[cl2];
cMeanX[cl1] += cMeanX[cl2];
cMeanY[cl1] += cMeanY[cl2];
cMeanXX[cl1] += cMeanXX[cl2];
cMeanYY[cl1] += cMeanYY[cl2];
cMeanYX[cl1] += cMeanYX[cl2];
cMeanXWt[cl1] += cMeanXWt[cl2];
cMeanYWt[cl1] += cMeanYWt[cl2];
cMeanXXWt[cl1] += cMeanXXWt[cl2];
cMeanYYWt[cl1] += cMeanYYWt[cl2];
cNMergedClus[cl1] += cNMergedClus[cl2];
cX[cl1]=cMeanXWt[cl1]/cSumWt[cl1];
cY[cl1]=cMeanYWt[cl1]/cSumWt[cl1];
nomerg=kFALSE;
nClus--;
}
}
if(nomerg || nClus==1) break;
}
// Filling containers:
for(Int_t cl=0; cl<nClsArr; cl++) {
if(!cStatus[cl]) continue;
HMdcList12GroupCells& cLCSeg2=cLCells2[cl];
Float_t errX=sqrt(cMeanXXWt[cl]/cSumWt[cl]-cX[cl]*cX[cl] +
prPlotSeg2->stX*prPlotSeg2->stX*0.25);
Float_t errY=sqrt(cMeanYYWt[cl]/cSumWt[cl]-cY[cl]*cY[cl] +
prPlotSeg2->stY*prPlotSeg2->stY*0.25);
// 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 |
Double_t sigma1,sigma2,alpha;
if(cNBins[cl]<2) sigma1=sigma2=alpha=0.;
else {
Double_t nBn=cNBins[cl];
Double_t eXX = cMeanXX[cl]*nBn- cMeanX[cl]*cMeanX[cl];
Double_t eYX = cMeanYX[cl]*nBn- cMeanY[cl]*cMeanX[cl];
Double_t eYY = cMeanYY[cl]*nBn- cMeanY[cl]*cMeanY[cl];
Double_t norm = (eXX+eYY)/2.;
Double_t eYX2 = eYX*eYX;
Double_t c=eXX*eYY-eYX2;
sigma1=norm+sqrt(norm*norm-c);
sigma2=norm-sqrt(norm*norm-c);
if(sigma2<0.) sigma2=0.;
Double_t e1=sigma1-eXX;
Double_t e2=sigma1-eYY;
sigma1=sqrt(sigma1)/nBn;
sigma2=sqrt(sigma2)/nBn;
alpha=atan2(sqrt(eYX2+e1*e1),sqrt(eYX2+e2*e2))*radToDeg;
if(eYX<0.) alpha=180.-alpha;
}
locClus[2]=counter;
HMdcClus* clus = (HMdcClus*)fClusCat->getSlot(locClus);
if(!clus) {
Warning("fillClusCut2","S.%i No slot HMdcClus available. Size of catMdcClus is %i !",sector+1,fClusCat->getEntries());
return cl-1;
}
if(isGeant) clus=(HMdcClus*)(new(clus) HMdcClusSim(cLCSeg2));
else clus=new(clus) HMdcClus(cLCSeg2);
indexLCh=fClusCat->getIndex(clus);
if(indexFCh<0) indexFCh=indexLCh;
counter++;
nClusters++;
Int_t nLayM1=clus->getNLayersMod(0);
Int_t nLayM2=clus->getNLayersMod(1);
clus->setAddress(locClus);
if(realTypeClFinder==2 && (nLayM1<minAmp[2] || nLayM2<minAmp[3])) {
clus->setMod( (nLayM1>=nLayM2) ? 2:3);
clus->setTypeClFinder(2);
} else {
clus->setMod((mSeg[1]==3) ? -2 : (mSeg[1]+1) ); // mSeg[1]+1=mSeg[1]-1+2
clus->setTypeClFinder(0);
}
clus->setMinCl(minAmp[2],minAmp[3]);
clus->setPrPlane(prPlotSeg2->A(),prPlotSeg2->B(),prPlotSeg2->D());
clus->setTarg(segOnKick[4]);
clus->setErrTarg(dirSeg1[4]);
clus->setSumWt(cSumWt[cl]);
clus->setNBins(cNBins[cl]);
clus->setXY(cX[cl],errX,cY[cl],errY);
clus->setNMergClust(cNMergedClus[cl]);
Int_t nDrTmM1=clus->getNDrTimes(0,5);
Int_t nDrTmM2=clus->getNDrTimes(6,11);
if(nDrTmM1>0) {
clus->setClusSizeM1(cNBins[cl]);
clus->setNDrTimesM1(nDrTmM1);
clus->setNMergClustM1(cNMergedClus[cl]);
clus->setShapeM1(sigma1,sigma2,alpha);
} else clus->clearMod1Par();
if(nDrTmM2>0) {
clus->setClusSizeM2(cNBins[cl]);
clus->setNDrTimesM2(nDrTmM2);
clus->setNMergClustM2(cNMergedClus[cl]);
clus->setShapeM2(sigma1,sigma2,alpha);
} else clus->clearMod2Par();
clus->setIndexParent(indexPar);
if(clus->isGeant()) {
((HMdcClusSim*)clus)->calcTrList();
((HMdcClusSim*)clus)->calcTrListMod(cLCSeg2,0);
((HMdcClusSim*)clus)->calcTrListMod(cLCSeg2,1);
}
}
return nClusters;
}
void HMdc34ClFinderSec::calcClusterSeg2(Int_t nBinF) {
Int_t nBin=nBinF;
stack->init();
stack->push(nBin);
Int_t minAmpS=minAmp[2]+minAmp[3]-1;
Int_t nLinesInCl=1;
Int_t nDownLine=nBin/prPlotSeg2->nBinX;
xMinClLines[0]=xMaxClLines[0]=nBin%prPlotSeg2->nBinX;
Int_t xMinCl=xMinClLines[0];
Int_t xMaxCl=xMinCl;
Int_t& cNBinsI=cNBins[nClsArr];
Double_t& cMeanXI=cMeanX[nClsArr];
Double_t& cMeanYI=cMeanY[nClsArr];
Double_t& cMeanXXI=cMeanXX[nClsArr];
Double_t& cMeanYYI=cMeanYY[nClsArr];
Double_t& cMeanYXI=cMeanYX[nClsArr];
Double_t& cSumWtI=cSumWt[nClsArr];
Double_t& cMeanXWtI=cMeanXWt[nClsArr];
Double_t& cMeanYWtI=cMeanYWt[nClsArr];
Double_t& cMeanXXWtI=cMeanXXWt[nClsArr];
Double_t& cMeanYYWtI=cMeanYYWt[nClsArr];
HMdcList12GroupCells& list=cLCells2[nClsArr];
cStatus[nClsArr]=kTRUE;
cNMergedClus[nClsArr]=1;
list.clear();
cNBinsI=0;
cMeanXI=cMeanYI=cMeanXXI=cMeanYYI=cMeanYXI=0.;
cSumWtI=cMeanXWtI=cMeanYWtI=cMeanXXWtI=cMeanYYWtI=0.;
while((nBin=stack->pop()) >= 0) {
Int_t nx=nBin%prPlotSeg2->nBinX;
Int_t ny=nBin/prPlotSeg2->nBinX;
Double_t wt=prPlotSeg2->weights[nBin]-minAmpS;
//region of cluster ???
Double_t x=prPlotSeg2->xBinsPos[nx];
Double_t y=prPlotSeg2->yBinsPos[ny];
cSumWtI += wt;
cMeanXWtI += x*wt;
cMeanYWtI += y*wt;
cMeanXXWtI += x*x*wt;
cMeanYYWtI += y*y*wt;
cNBinsI++;
cMeanXI += x;
cMeanYI += y;
cMeanXXI += x*x;
cMeanYYI += y*y;
cMeanYXI += y*x;
// At the scaning of the cluster y step(in bins) can't be > 1
Int_t nyL=ny-nDownLine;
if(nyL>=nLinesInCl) {
if(nyL<300) {
xMinClLines[nyL]=xMaxClLines[nyL]=nx;
nLinesInCl++;
} else {
if(!HMdc34ClFinder::getQuietMode())
{
Warning("calcClusterSeg2","Sec.%i The cluster size is very big!",
sector+1);
}
}
} else {
if(nx<xMinClLines[nyL]) xMinClLines[nyL]=nx;
else if(nx>xMaxClLines[nyL]) xMaxClLines[nyL]=nx;
}
if(nx<xMinCl) xMinCl=nx;
else if(nx>xMaxCl) xMaxCl=nx;
for(Int_t ib=0; ib<8; ib++) {
Int_t nBinTs=nBin+nearbyBins[ib];
if( HMdcBArray::testAndUnset(prPlotSeg2->plotBArSc,nBinTs) ) stack->push(nBinTs);
}
stack->checkSize(); // If stack overflowed one will increased on 50%:
}
// Filling of cluster:
Int_t nUpLine=nDownLine+nLinesInCl-1;
for(Int_t mod=2; mod<4; mod++) {
if(!lMods[mod]) continue;
HMdc34ClFinderMod& fMod=(*this)[mod];
Int_t startLay=(mod-2)*6;
for(Int_t lay=0; lay<6; lay++) {
Int_t layer=startLay+lay;
HMdc34ClFinderLayer& fLay=fMod[lay];
for(Int_t cell=fLay.getFirstCell(); cell>=0; cell=fLay.next(cell)) {
if( fLay.yBin[xMinCl]+fLay.cellShUp[cell]<nDownLine &&
fLay.yBin[xMaxCl]+fLay.cellShUp[cell]<nDownLine ) continue;
if( fLay.yBin[xMinCl]+fLay.cellShDown[cell]>nUpLine &&
fLay.yBin[xMaxCl]+fLay.cellShDown[cell]>nUpLine ) break;
for(Int_t line=0; line<nLinesInCl; line++) {
Int_t nYbin=line+nDownLine;
Int_t tYmin=fLay.yBin[xMinClLines[line]]-nYbin;
Int_t tYmax=fLay.yBin[xMaxClLines[line]]-nYbin;
Int_t dY1b=tYmin+fLay.cellShUp[cell];
Int_t dY2b=tYmax+fLay.cellShUp[cell];
//if (dY1b==0 || dY2b==0) - optimizaciya, proverit' na time !!!!???
if(dY1b<0 && dY2b<0) break;
if(dY1b) dY1b=(dY1b>0) ? 1:-1;
if(dY2b) dY2b=(dY2b>0) ? 1:-1;
// if dY1b*dY2b<=0 the line cross cell projection
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 the line cross the cell projection
// dY1a*dY1b<0 The line is in the cell projection
if(dY1a*dY2a>0 && dY1a*dY1b > 0) continue;
}
Int_t nDeleted=list.setTime(layer,cell,fLay.getTime(cell));
if(nDeleted&&!HMdc34ClFinder::getQuietMode())
{
Warning("calcClusterSeg2",
"%i cells was removed (size of cluster is very big)",
nDeleted);
}
break;
}
}
}
}
cX[nClsArr]=cMeanXWtI/cSumWtI;
cY[nClsArr]=cMeanYWtI/cSumWtI;
}
void HMdc34ClFinderSec::printClusters(void){
// clusters.print();
}
//---------------------------------------------------------
HMdc34ClFinder* HMdc34ClFinder::fMdc34ClFinder=0;
Bool_t HMdc34ClFinder::quietmode=kTRUE;
HMdc34ClFinder::HMdc34ClFinder(const char* name,const char* title,
const char* context)
: HParSet(name,title,context) {
// constructor creates an array of pointers of type HMdc34ClFinderSec
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 = HMdcGetContainers::getCatMdcClus(kTRUE);
}
Bool_t HMdc34ClFinder::initContainer(void) {
if( !fMdcDet || !fSizesCells->initContainer() ||
!HMdcGetContainers::isInited(fSpecGeomPar) ||
!HMdcGetContainers::isInited(fKickPlane) ||
!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();
HMdc34ClFinderSec::fkick=fKickPlane;
for (Int_t sec = 0; sec < 6; sec++) {
if(!fMdcDet->isSectorActive(sec)) continue;
if( !(*array)[sec] ) {
Int_t nLM2=HMdcTrackDSet::getTrFnNLayersInMod()[sec*4+2];
Int_t nLM3=HMdcTrackDSet::getTrFnNLayersInMod()[sec*4+3];
if((mods[sec*4+2]&&nLM2>0) || (mods[sec*4+3]&&nLM3>0)) { // for mod.3&4 now !!!
(*array)[sec] = new HMdc34ClFinderSec(sec,320,886); //????0,0???
(*this)[sec].fClusCat=fMdcClusCat;
} else continue;
}
// initialization of container ---
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() {
// destructor
if(array) {
array->Delete();
delete array;
}
fMdc34ClFinder=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) { // *!*
// clears the container
for(Int_t s=0;s<6;s++) if((*array)[s]) (*this)[s].clear();
}
void HMdc34ClFinder::printClusters(Int_t sec){
Int_t s1=0;
Int_t s2=5;
if( sec>=0 && sec<=5 ) s1=s2=sec;
for(Int_t s=s1;s<=s2;s++) {
if( (*array)[s] ) {
printf("-Sector %i ----Clusters:n",s+1);
HMdc34ClFinderSec& sec=(*this)[s];
sec.printClusters();
}
}
}
Bool_t HMdc34ClFinder::calcTarget(Int_t sec){
//Geting size of target
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];
fsec.target[0]=(fSpecGeomPar->getTarget(0)->getTransform()).getTransVector();
fsec.target[0].setZ( fsec.target[0].getZ() +
fSpecGeomPar->getTarget(0)->getPoint(0)->getZ());
fsec.target[1]=(fSpecGeomPar->getTarget(nT)->getTransform()).getTransVector();
fsec.target[1].setZ( fsec.target[1].getZ() +
fSpecGeomPar->getTarget(nT)->getPoint(2)->getZ());
const HGeomTransform* trans=(*fSizesCells)[sec].getLabTrans();
if(&trans == 0) return kFALSE;
fsec.target[0]=trans->transTo(fsec.target[0]);
fsec.target[1]=trans->transTo(fsec.target[1]);
return kTRUE;
}
Bool_t HMdc34ClFinder::calcProjPlaneSeg2(Int_t sec){
// Calculation of project plane:
//
// a plane between MDC3 and MDC4 - -1.<par<1.
// plane 1,2,3,4,5 or 6 MDC3 - par=-6,-5,-4,-3,-2 or -1
// plane 1,2,3,4,5 or 6 MDC2 - par=+1,+2,+3,+4,+5 or +6
// midplane MDC1 - proj.plane=-7
// midplane MDC2 - proj.plane=+7
// Automatical selection: par>=10. (for two MDC
// in sector = 0.15, for one MDC - midplane of this MDC)
// Float_t par=-7.; //In a future will be geted from Par.Container
if(!fSizesCells->hasChanged()) return kTRUE;
HMdc34ClFinderSec& fsec=(*this)[sec];
HMdcSizesCellsSec& fSCellsSec=(*fSizesCells)[sec];
if( !&fSCellsSec ) return kFALSE;
Float_t par=10.; //In a future will be geted from Par.Container
Int_t nL=(Int_t)par;
// (mdc3 don't exist) (mdc4 don't exist)
if(nL<-7||nL>7 || (nL<0 && fsec.mSeg[1]==2) || (nL>0 && fsec.mSeg[1]==1) ||
//(one mdc exist only)
(nL==0 && fsec.mSeg[1]!=3)) {
nL=0;
par=0.1;
}
Char_t *text="as projection plane";
HGeomTransform prPl;
if( (nL==0 && fsec.mSeg[1]<3) || nL==-7 || nL==7 ) {
// the project plane - the middle plane of MDC3 or MDC4
Int_t mod=(nL==0) ? fsec.mSeg[1]+1 : (nL+7)/14+2;
prPl=*(fSCellsSec[mod].getSecTrans());
printf("n===> Sec.%i Seg.2: Using middle plane of MDC%i %sn",
sec+1,mod+1,text);
} else if( nL==0 && fsec.mSeg[1]==3 ) {
// the project plane - between MDC3 & MDC4
const HGeomTransform* trLayer6=fSCellsSec[2][5].getSecTrans();
const HGeomTransform* trLayer1=fSCellsSec[3][0].getSecTrans();
// MDC3 & MDC4 are ~parallel:
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) { // proj.plane will parallel to MDC3
prPl.setTransform(*trLayer6);
mult=newDist/distToL6;
} else { // proj.plane will parallel to MDC4
prPl.setTransform(*trLayer1);
mult=newDist/distToL1;
}
HGeomVector prTr(prPl.getTransVector());
prTr*=mult;
prPl.setTransVector(prTr);
printf("n===> Sec.%i Seg.2: Using plane between MDC 3 & 4 (p=%g) %sn",
sec+1,par,text);
} else {
// the project plane - one of the layers
Int_t mod=(nL<0) ? 2:3;
if(nL<0) nL+=7;
prPl=*(fSCellsSec[mod][nL-1].getSecTrans());
printf("n===> Sec.%i Seg.2: Using MDC%i, layer %i %sn",
sec+1,mod+1,nL,text);
}
fsec.prPlotSeg2->setPlanePar(prPl);
return kTRUE;
}
Bool_t HMdc34ClFinder::calcSizePlotSeg2(Int_t sec){
// Calculation of plot's size:
if( !fKickPlane->hasChanged() && !fMdcGeomPar->hasChanged() &&
!fSizesCells->hasChanged()) return kTRUE;
HMdc34ClFinderSec& fsec=(*this)[sec];
if(!fsec.mSeg[0]) return kFALSE;
HMdcSizesCellsSec& fSCellsSec=(*fSizesCells)[sec];
if( !&fSCellsSec ) return kFALSE;
HGeomVector vect[2][4];
HGeomVector pKick[4];
const HGeomTransform* trans=0;
HGeomVector newP[4][2]; //[4]-num.lines, [2]-the firt&last points of the line
HMdcPlane plane;
for(Int_t mod=0; mod<2; mod++) {
if(!(fsec.lMods[mod])) 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=(fsec.mSeg[0]==3) ? mod : 0; // for one MDC2 in seg. only
vect[indx][point]=*fpoint;
vect[indx][point].setZ(0.);
vect[indx][point]=trans->transFrom(vect[indx][point]);
if(mod==1 && fsec.mSeg[0]==3) {
Int_t nTag=(point==0 || point==3) ? 0:1;
vect[0][point]-=fsec.target[nTag];
plane.calcIntersection(fsec.target[nTag],vect[0][point],vect[0][point]);
}
}
}
if(fsec.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]-fsec.target[0];
else dir=vect[0][point]-fsec.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;
}
}
// pKick[0-3] pointers on the kick-plane
Float_t yMin=1.e+10;
Float_t yMax=-1.e+10;
for(Int_t mod=2; mod<4; mod++) {
if(!(fsec.lMods[mod])) continue;
HGeomCompositeVolume *fVolMdc=fGetCont->getGeomCompositeVolume(mod);
Int_t indx=(fsec.mSeg[1]==3) ? mod-2 : 0; // for one MDC2 in seg. only
trans=fSCellsSec[mod].getSecTrans();
for(Int_t point=0; point<4; point++) {
vect[indx][point]=*(fVolMdc->getPoint(point));
vect[indx][point].setZ(0.); //midplane of MDC
vect[indx][point]=trans->transFrom(vect[indx][point]);
vect[indx][point]-=pKick[point];
fsec.prPlotSeg2->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);
}
}
Float_t xMaxD=xLine(vect[0][0],vect[0][1],yMin);
Float_t xMinD=xLine(vect[0][2],vect[0][3],yMin);
Float_t xMax=xLine(vect[0][0],vect[0][1],yMax);
Float_t xMin=xLine(vect[0][2],vect[0][3],yMax);
if(fsec.mSeg[1]==3) {
Float_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.prPlotSeg2->setEdges(yMin, xMinD, xMaxD, yMax, xMin, xMax);
return kTRUE;
}
void HMdc34ClFinder::calcCrossLines(HGeomVector& p1l1, HGeomVector& p2l1,
HGeomVector& p1l2, HGeomVector& p2l2, HGeomVector& cross) {
// Calculeting a cross of 2 lines (p1l1-p2l1, p1l2-p2l2) on the one (!)
// plane. Z seted to 0.
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.);
}
Float_t HMdc34ClFinder::xLine(HGeomVector& p1,HGeomVector& p2,Float_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];
if(!fsec.mSeg[0]) return kFALSE;
HMdcSizesCellsSec& fSCellsSec=(*fSizesCells)[sec];
if( !&fSCellsSec ) return kFALSE;
HGeomVector wire[2];
HGeomVector rib[4];
HGeomVector point[2];
Int_t mod=(fsec.mSeg[0]>1) ? 1:0;
const HGeomTransform* trans=fSCellsSec[mod].getSecTrans();
if(!trans) return kFALSE;
HGeomCompositeVolume *fVolMdc=fGetCont->getGeomCompositeVolume(mod);
if(!fVolMdc) return kFALSE;
HGeomVector tag=fsec.target[0]+fsec.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.lMods[mod]) continue;
fVolMdc=fGetCont->getGeomCompositeVolume(mod);
if(!fVolMdc) return kFALSE;
HMdcSizesCellsMod& fSCellsMod=fSCellsSec[mod];
if(!&fSCellsMod) return kFALSE;
HMdc34ClFinderMod& fMod=fsec[mod];
for(Int_t lay=0; lay<6; lay++) {
HMdcSizesCellsLayer& fSCellsLay=fSCellsMod[lay];
HGeomVolume* fVolLayer=fVolMdc->getComponent(lay);
if(!fVolLayer) {
Error("calcWiresProj",
"The HGeomVolume (S%i|M%i|L%i) doesn't exist!n",
sec+1,mod+1,lay+1);
return kFALSE;
}
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(); // y=a*x+ln[1]
// calc. of the wire's projection
Double_t bWire=(100.*pitch-wOfSet)/cosWiOr; // 100.-cell n.100
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;
fsec.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
Float_t al=(point[0](1)-point[1](1))/(point[0](0)-point[1](0)); // (y1-y2)/(x1-x2)
Float_t x=(fsec.prPlotSeg2->nBinX-1) * fsec.prPlotSeg2->stX +
fsec.prPlotSeg2->xMin;
Float_t yLeft=al*( x - point[1](0))+point[1](1); // (x-x2)*a+y2
Float_t yRight=al*(fsec.prPlotSeg2->xMin - point[1](0))+point[1](1);
Short_t nYLeft =Short_t((yLeft - fsec.prPlotSeg2->yMin) /
fsec.prPlotSeg2->stY);
Short_t nYRight=Short_t((yRight - fsec.prPlotSeg2->yMin) /
fsec.prPlotSeg2->stY);
Short_t yShift=TMath::Min(nYLeft,nYRight);
if(TMath::Abs(nYLeft-nYRight)>=fsec.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,fsec.prPlotSeg2->nBinY);
return kFALSE;
}
HMdc34ClFinderLayer& fLayer=fMod[lay];
if( !fLayer.createArrayBins(fsec.prPlotSeg2->nBinX) ) return kFALSE;
for(Int_t nx=0; nx<fLayer.nBinX; nx++) {
Float_t x=nx * fsec.prPlotSeg2->stX + fsec.prPlotSeg2->xMin;
Float_t y=al*(x-point[1](0))+point[1](1);
Short_t nY=Short_t((y - fsec.prPlotSeg2->yMin) / fsec.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;
}
}
// calc. of the position of the cells' ribs:
// Z(x)=A1*x+Y0[nCell]; Z(x)=A2*x+Z0[nCell]; Y0&Z0 - at x=0;
// Y0[nCell]=Ystep*nCell + Yshift; Z0[nCell]=Zstep*nCell + Zshift(=0!);
//--- cells 0 & 100 ---
for(Int_t cell=0; cell<101; cell+=100) {
bWire=(cell*pitch-wOfSet-pitch/2.)/cosWiOr;
// rib[0]-rib[1] - the down and left rib of cell
// rib[2]-rib[3] - the down and right rib of cell
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+bWire, 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.;
//Calc. the point(index) number in HMdc34ClFinderSec::segOnKick[4]:
if(rib[1](2)+fsec.prPlotSeg2->B()*rib[1](1)-fsec.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)+fsec.prPlotSeg2->B()*rib[3](1)-fsec.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::setCellsList(HMdcEvntListCells& event) {
for(Int_t s=0;s<6;s++) if((*array)[s]) (*this)[s].setCellsList(event[s]);
}
void HMdc34ClFinderSec::setCellsList(HMdcSecListCells& event) {
counter=0;
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]);
}
ROOT page - Class index - Class Hierarchy - Top of the page
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.