AliPMDClusteringV2::AliPMDClusteringV2():
fPMDclucont(new TObjArray()),
- fCutoff(0.0)
+ fCutoff(0.0),
+ fClusParam(0)
{
for(int i = 0; i < kNDIMX; i++)
{
AliPMDClusteringV2::AliPMDClusteringV2(const AliPMDClusteringV2& pmdclv2):
AliPMDClustering(pmdclv2),
fPMDclucont(0),
- fCutoff(0)
+ fCutoff(0),
+ fClusParam(0)
{
// copy constructor
AliError("Copy constructor not allowed ");
// ------------------------------------------------------------------------ //
void AliPMDClusteringV2::DoClust(Int_t idet, Int_t ismn,
- Double_t celladc[48][96], TObjArray *pmdcont)
+ Int_t celltrack[48][96],
+ Int_t cellpid[48][96],
+ Double_t celladc[48][96],
+ TObjArray *pmdcont)
{
// main function to call other necessary functions to do clustering
//
AliPMDcluster *pmdcl = 0;
- const float ktwobysqrt3 = 1.1547; // 2./sqrt(3.)
-
- Int_t i, j, nmx1, incr, id, jd;
+ const Float_t ktwobysqrt3 = 1.1547; // 2./sqrt(3.)
+ const Int_t kNmaxCell = 19; // # of cells surrounding a cluster center
+ Int_t i = 0, j = 0, nmx1 = 0;
+ Int_t incr = 0, id = 0, jd = 0;
Int_t ndimXr = 0;
Int_t ndimYr = 0;
- Int_t celldataX[15], celldataY[15];
- Float_t clusdata[6];
- Double_t cutoff, ave;
+ Int_t celldataX[kNmaxCell], celldataY[kNmaxCell];
+ Int_t celldataTr[kNmaxCell], celldataPid[kNmaxCell];
+ Float_t celldataAdc[kNmaxCell];
+ Float_t clusdata[6] = {0.,0.,0.,0.,0.,0.};
+ Double_t cutoff = 0., ave = 0.;
Double_t edepcell[kNMX];
ndimYr = 96;
}
- for (Int_t i =0; i < kNMX; i++)
+ for (i =0; i < kNMX; i++)
{
edepcell[i] = 0.;
}
}
}
- Int_t iord1[kNMX];
- TMath::Sort(kNMX,edepcell,iord1);// order the data
+ // the dimension of iord1 is increased twice
+ Int_t iord1[2*kNMX];
+ TMath::Sort((Int_t)kNMX,edepcell,iord1);// order the data
cutoff = fCutoff; // cutoff used to discard cells having ener. dep.
ave = 0.;
nmx1 = -1;
//
// Cells associated with a cluster
//
- for (Int_t ihit = 0; ihit < 15; ihit++)
+ for (Int_t ihit = 0; ihit < kNmaxCell; ihit++)
{
Int_t dummyXY = pmdcludata->GetCellXY(ihit);
Int_t celldumY = dummyXY%10000;
Int_t celldumX = dummyXY/10000;
- Float_t cellY = (Float_t) (ktwobysqrt3*celldumY/10);
- Float_t cellX = (Float_t) (celldumX/10 - (celldumY/2.)/10);
-
+ Float_t cellY = (Float_t) celldumY/10;
+ Float_t cellX = (Float_t) celldumX/10;
+
//
// Cell X centroid is back transformed
//
if (ismn < 12)
{
- celldataX[ihit] = (Int_t) (cellX - (24-1) + cellY/2.);
+ celldataX[ihit] = (Int_t) ((cellX - (24-1) + cellY/2.) + 0.5);
}
else if (ismn >= 12 && ismn <= 23)
{
- celldataX[ihit] = (Int_t) (cellX - (48-1) + cellY/2.);
+ celldataX[ihit] = (Int_t) ((cellX - (48-1) + cellY/2.) + 0.5 );
}
- celldataY[ihit] = (Int_t) cellY;
+ celldataY[ihit] = (Int_t) (cellY + 0.5);
+
+ Int_t irow = celldataX[ihit];
+ Int_t icol = celldataY[ihit];
+
+ if ((irow >= 0 && irow < 48) && (icol >= 0 && icol < 96))
+ {
+ celldataTr[ihit] = celltrack[irow][icol];
+ celldataPid[ihit] = cellpid[irow][icol];
+ celldataAdc[ihit] = (Float_t) celladc[irow][icol];
+ }
+ else
+ {
+ celldataTr[ihit] = -1;
+ celldataPid[ihit] = -1;
+ celldataAdc[ihit] = -1;
+ }
+
}
- pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY);
+ pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY,
+ celldataTr, celldataPid, celldataAdc);
pmdcont->Add(pmdcl);
}
- fPMDclucont->Clear();
+ fPMDclucont->Delete();
}
// ------------------------------------------------------------------------ //
Int_t AliPMDClusteringV2::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
// connected cells
//
- Int_t i,j,k,id1,id2,icl, numcell;
- Int_t jd1,jd2, icell, cellcount;
+ Int_t i = 0, j = 0, k = 0, id1 =0, id2 = 0, icl = 0, numcell = 0;
+ Int_t jd1 = 0, jd2 = 0, icell = 0, cellcount = 0;
Int_t clust[2][5000];
static Int_t neibx[6] = {1,0,-1,-1,0,1}, neiby[6] = {0,1,1,0,-1,-1};
// Takes the big patch and does gaussian fitting and
// finds out the more refined clusters
+ const Float_t ktwobysqrt3 = 1.1547;
+ const Int_t kNmaxCell = 19;
+
AliPMDcludata *pmdcludata = 0;
+
+ Int_t i12 = 0;
+ Int_t i = 0, j = 0, k = 0;
+ Int_t i1 = 0, i2 = 0, id = 0, icl = 0, itest = 0, ihld = 0;
+ Int_t ig = 0, nsupcl = 0, clno = 0, clX = 0, clY = 0;
+ Int_t clxy[kNmaxCell];
+
+ Float_t clusdata[6] = {0.,0.,0.,0.,0.,0.};
+ Double_t x1 = 0., y1 = 0., z1 = 0., x2 = 0., y2 = 0., z2 = 0., rr = 0.;
+
+ Int_t kndim = incr + 1;
+
TArrayI testncl;
TArrayI testindex;
- const Int_t kndim = 4500;
- Int_t i, j, k, i1, i2, id, icl, itest, ihld;
- Int_t ig, nsupcl, clno, clX,clY;
- Int_t clxy[15];
- Int_t ncl[kndim], iord[kndim];
- Float_t clusdata[6];
- Double_t x1, y1, z1, x2, y2, z2, rr;
- Double_t x[kndim], y[kndim], z[kndim];
- Double_t xc[kndim], yc[kndim], zc[kndim], cells[kndim];
- Double_t rcl[kndim], rcs[kndim];
+
+ Int_t *ncl, *iord;
+
+ Double_t *x, *y, *z, *xc, *yc, *zc, *cells, *rcl, *rcs;
+
+ ncl = new Int_t [kndim];
+ iord = new Int_t [kndim];
+ x = new Double_t [kndim];
+ y = new Double_t [kndim];
+ z = new Double_t [kndim];
+ xc = new Double_t [kndim];
+ yc = new Double_t [kndim];
+ zc = new Double_t [kndim];
+ cells = new Double_t [kndim];
+ rcl = new Double_t [kndim];
+ rcs = new Double_t [kndim];
for(Int_t kk = 0; kk < 15; kk++)
{
clno = -1;
nsupcl = -1;
- for(i = 0; i < 4500; i++)
+
+ for(i = 0; i < kndim; i++)
{
ncl[i] = -1;
}
- for(i = 0; i < incr; i++)
+ for(i = 0; i <= incr; i++)
{
if(fInfcl[0][i] != nsupcl)
{
id = -1;
icl = -1;
- for(i = 0; i < nsupcl; i++)
+ for(i = 0; i <= nsupcl; i++)
{
if(ncl[i] == 0)
{
clno++;
i1 = fInfcl[1][id];
i2 = fInfcl[2][id];
- Int_t i12 = i1 + i2*kNDIMX;
+ i12 = i1 + i2*kNDIMX;
clusdata[0] = fCoord[0][i1][i2];
clusdata[1] = fCoord[1][i1][i2];
clusdata[2] = edepcell[i12];
clusdata[5] = 0.0;
//cell information
- clX = (Int_t) fCoord[0][i1][i2]*10;
- clY = (Int_t) fCoord[1][i1][i2]*10;
- clxy[0] = clX*10000 + clY;
- for(Int_t icltr = 1; icltr < 15; icltr++)
+
+ clY = (Int_t)((ktwobysqrt3*fCoord[1][i1][i2])*10);
+ clX = (Int_t)((fCoord[0][i1][i2] - clY/20.)*10);
+ clxy[0] = clX*10000 + clY ;
+
+ for(Int_t icltr = 1; icltr < kNmaxCell; icltr++)
{
clxy[icltr] = -1;
}
clno++;
i1 = fInfcl[1][id];
i2 = fInfcl[2][id];
- Int_t i12 = i1 + i2*kNDIMX;
+ i12 = i1 + i2*kNDIMX;
x1 = fCoord[0][i1][i2];
y1 = fCoord[1][i1][i2];
clusdata[4] = (TMath::Sqrt(z1*z2))/(z1+z2);
clusdata[5] = 0.0;
- clX = (Int_t) x1*10;
- clY = (Int_t) y1*10;
- clxy[0] = clX*10000 + clY;
-
- clX = (Int_t) x2*10;
- clY = (Int_t) y2*10;
- clxy[1] = clX*10000 + clY;
-
- for(Int_t icltr = 2; icltr < 15; icltr++)
+ clY = (Int_t)((ktwobysqrt3*y1)*10);
+ clX = (Int_t)((x1 - clY/20.)*10);
+ clxy[0] = clX*10000 + clY ;
+
+ clY = (Int_t)((ktwobysqrt3*y2)*10);
+ clX = (Int_t)((x2 - clY/20.)*10);
+ clxy[1] = clX*10000 + clY ;
+
+ for(Int_t icltr = 2; icltr < kNmaxCell; icltr++)
{
clxy[icltr] = -1;
}
i1 = fInfcl[1][id];
i2 = fInfcl[2][id];
- Int_t i12 = i1 + i2*kNDIMX;
+ i12 = i1 + i2*kNDIMX;
x[0] = fCoord[0][i1][i2];
y[0] = fCoord[1][i1][i2];
id++;
i1 = fInfcl[1][id];
i2 = fInfcl[2][id];
- Int_t i12 = i1 + i2*kNDIMX;
+ i12 = i1 + i2*kNDIMX;
iord[j] = j;
x[j] = fCoord[0][i1][i2];
y[j] = fCoord[1][i1][i2];
zc[ig] = z[iord[j]];
}
}
- ClustDetails(ncl[i], ig, x[0], y[0] ,z[0], xc[0], yc[0], zc[0],
- rcl[0], rcs[0], cells[0], testncl, testindex);
+ ClustDetails(ncl[i], ig, x, y ,z, xc, yc, zc, rcl, rcs, cells,
+ testncl, testindex);
Int_t pp = 0;
for(j = 0; j <= ig; j++)
clusdata[5] = rcs[j];
if(ig == 0)
{
- clusdata[3] = ncl[i];
+ clusdata[3] = ncl[i] + 1;
}
else
{
}
// cell information
Int_t ncellcls = testncl[j];
- if( ncellcls < 15 )
+ if( ncellcls < kNmaxCell )
{
for(Int_t kk = 1; kk <= ncellcls; kk++)
{
Int_t ll = testindex[pp];
- clX = (Int_t) x[ll]*10;
- clY = (Int_t) y[ll]*10;
- clxy[kk-1] = (Int_t) clX*10000 + clY ;
+ clY = (Int_t)((ktwobysqrt3*y[ll])*10);
+ clX = (Int_t)((x[ll] - clY/20.)*10);
+ clxy[kk-1] = clX*10000 + clY ;
+
pp++;
}
- for(Int_t icltr = ncellcls ; icltr < 15; icltr++)
+ for(Int_t icltr = ncellcls ; icltr < kNmaxCell; icltr++)
{
clxy[icltr] = -1;
}
testindex.Set(0);
}
}
+ delete [] ncl;
+ delete [] iord;
+ delete [] x;
+ delete [] y;
+ delete [] z;
+ delete [] xc;
+ delete [] yc;
+ delete [] zc;
+ delete [] cells;
+ delete [] rcl;
+ delete [] rcs;
}
// ------------------------------------------------------------------------ //
-void AliPMDClusteringV2::ClustDetails(Int_t ncell, Int_t nclust, Double_t &x,
- Double_t &y, Double_t &z, Double_t &xc,
- Double_t &yc, Double_t &zc,
- Double_t &rcl, Double_t &rcs,
- Double_t &cells, TArrayI &testncl,
+void AliPMDClusteringV2::ClustDetails(Int_t ncell, Int_t nclust, Double_t x[],
+ Double_t y[], Double_t z[],Double_t xc[],
+ Double_t yc[], Double_t zc[],
+ Double_t rcl[], Double_t rcs[],
+ Double_t cells[], TArrayI &testncl,
TArrayI &testindex)
{
// function begins
//
- const Int_t kndim1 = 2000;
- const Int_t kndim2 = 10;
- const Int_t kndim3 = 400;
-
- Int_t i, j, k, i1, i2;
- Int_t cluster[kndim1][kndim2], cell[kndim1][kndim3];
-
- Double_t x1, y1, x2, y2, rr, b, c, r1, r2;
- Double_t sumx, sumy, sumxy, sumxx, sum, sum1, sumyy;
- Double_t xx[kndim1], yy[kndim1], zz[kndim1];
- Double_t xxc[kndim1], yyc[kndim1],clustcell[kndim3][kndim1];
- Double_t str[kndim1], str1[kndim1],xcl[kndim1], ycl[kndim1], cln[kndim1];
+ Int_t kndim1 = ncell + 1;//ncell
+ Int_t kndim2 = 20;
+ Int_t kndim3 = nclust + 1;//nclust
- for(i = 0; i <= nclust; i++)
- {
- xxc[i] = *(&xc+i);
- yyc[i] = *(&yc+i);
- str[i] = 0.;
- str1[i] = 0.;
- }
- for(i = 0; i <= ncell; i++)
- {
- xx[i] = *(&x+i);
- yy[i] = *(&y+i);
- zz[i] = *(&z+i);
- }
+ Int_t i = 0, j = 0, k = 0, i1 = 0, i2 = 0;
+ Double_t x1 = 0., y1 = 0., x2 = 0., y2 = 0.;
+ Double_t rr = 0., b = 0., c = 0., r1 = 0., r2 = 0.;
+ Double_t sumx = 0., sumy = 0., sumxy = 0.;
+ Double_t sumxx = 0., sum = 0., sum1 = 0., sumyy = 0.;
- // INITIALIZE
+ Double_t *str, *str1, *xcl, *ycl, *cln;
+ Int_t **cell;
+ Int_t ** cluster;
+ Double_t **clustcell;
+ str = new Double_t [kndim3];
+ str1 = new Double_t [kndim3];
+ xcl = new Double_t [kndim3];
+ ycl = new Double_t [kndim3];
+ cln = new Double_t [kndim3];
+ clustcell = new Double_t *[kndim3];
+ cell = new Int_t *[kndim3];
+ cluster = new Int_t *[kndim1];
for(i = 0; i < kndim1; i++)
{
- for(j = 0; j < kndim2; j++)
+ cluster[i] = new Int_t [kndim2];
+ }
+
+ for(i = 0; i < kndim3; i++)
+ {
+ str[i] = 0;
+ str1[i] = 0;
+ xcl[i] = 0;
+ ycl[i] = 0;
+ cln[i] = 0;
+
+ cell[i] = new Int_t [kndim2];
+ clustcell[i] = new Double_t [kndim1];
+ for(j = 0; j < kndim1; j++)
{
- cluster[i][j] = 0;
-
+ clustcell[i][j] = 0;
}
- for(j = 0; j < kndim3; j++)
+ for(j = 0; j < kndim2; j++)
{
- cell[i][j] = 0;
- clustcell[j][i] = 0.;
+ cluster[i][j] = 0;
+ cell[i][j] = 0;
}
}
-
-
+
if(nclust > 0)
{
// more than one cluster
for (i = 0; i <= ncell; i++)
{
- x1 = xx[i];
- y1 = yy[i];
+ x1 = x[i];
+ y1 = y[i];
cluster[i][0] = 0;
+
// distance <= 1 cell unit
+
for(j = 0; j <= nclust; j++)
{
- x2 = xxc[j];
- y2 = yyc[j];
+ x2 = xc[j];
+ y2 = yc[j];
rr = Distance(x1, y1, x2, y2);
if(rr <= 1.)
{
{
for(j=0; j<=nclust; j++)
{
- x2 = xxc[j];
- y2 = yyc[j];
+ x2 = xc[j];
+ y2 = yc[j];
rr = Distance(x1, y1, x2, y2);
if(rr <= TMath::Sqrt(3.))
{
{
for(j=0; j<=nclust; j++)
{
- x2 = xxc[j];
- y2 = yyc[j];
+ x2 = xc[j];
+ y2 = yc[j];
rr = Distance(x1, y1, x2, y2);
if(rr <= 2.)
{
{
for(j = 0; j <= nclust; j++)
{
- x2 = xxc[j];
- y2 = yyc[j];
+ x2 = xc[j];
+ y2 = yc[j];
rr = Distance(x1, y1, x2, y2);
if(rr <= 2.7)
{
for(j = 1; j <= i1; j++)
{
i2 = cluster[i][j];
- str[i2] += zz[i]/i1;
+ str[i2] += z[i]/i1;
}
}
}
for(j = 1; j <= i1; j++)
{
i2 = cluster[i][j];
- str1[i2] += zz[i]*str[i2]/sum;
- clustcell[i2][i] = zz[i]*str[i2]/sum;
+ str1[i2] += z[i]*str[i2]/sum;
+ clustcell[i2][i] = z[i]*str[i2]/sum;
}
}
}
{
if(clustcell[i][j] != 0)
{
- sumx += clustcell[i][j]*xx[j];
- sumy += clustcell[i][j]*yy[j];
+ sumx += clustcell[i][j]*x[j];
+ sumy += clustcell[i][j]*y[j];
sum += clustcell[i][j];
- sum1 += clustcell[i][j]/zz[j];
+ sum1 += clustcell[i][j]/z[j];
}
}
//** xcl and ycl are cluster centroid positions ( center of gravity )
sumxy = 0.;
for(j = 0; j <= ncell; j++)
{
- sumxx += clustcell[i][j]*(xx[j]-xcl[i])*(xx[j]-xcl[i])/sum;
- sumyy += clustcell[i][j]*(yy[j]-ycl[i])*(yy[j]-ycl[i])/sum;
- sumxy += clustcell[i][j]*(xx[j]-xcl[i])*(yy[j]-ycl[i])/sum;
+ sumxx += clustcell[i][j]*(x[j]-xcl[i])*(x[j]-xcl[i])/sum;
+ sumyy += clustcell[i][j]*(y[j]-ycl[i])*(y[j]-ycl[i])/sum;
+ sumxy += clustcell[i][j]*(x[j]-xcl[i])*(y[j]-ycl[i])/sum;
}
b = sumxx+sumyy;
c = sumxx*sumyy-sumxy*sumxy;
r1 = b/2.+TMath::Sqrt(b*b/4.-c);
r2 = b/2.-TMath::Sqrt(b*b/4.-c);
// final assignments to proper external variables
- *(&xc + i) = xcl[i];
- *(&yc + i) = ycl[i];
- *(&zc + i) = str[i];
- *(&cells + i) = cln[i];
- *(&rcl+i) = r1;
- *(&rcs+i) = r2;
+ xc[i] = xcl[i];
+ yc[i] = ycl[i];
+ zc[i] = str[i];
+ cells[i] = cln[i];
+ rcl[i] = r1;
+ rcs[i] = r2;
+
}
//To get the cell position in a cluster
}
}
- else
+ else if(nclust == 0)
{
sumx = 0.;
sumy = 0.;
i = 0;
for(j = 0; j <= ncell; j++)
{
- sumx += zz[j]*xx[j];
- sumy += zz[j]*yy[j];
- sum += zz[j];
+ sumx += z[j]*x[j];
+ sumy += z[j]*y[j];
+ sum += z[j];
sum1++;
}
xcl[i] = sumx/sum;
sumxy = 0.;
for(j = 0; j <= ncell; j++)
{
- sumxx += clustcell[i][j]*(xx[j]-xcl[i])*(xx[j]-xcl[i])/sum;
- sumyy += clustcell[i][j]*(yy[j]-ycl[i])*(yy[j]-ycl[i])/sum;
- sumxy += clustcell[i][j]*(xx[j]-xcl[i])*(yy[j]-ycl[i])/sum;
+ sumxx += clustcell[i][j]*(x[j]-xcl[i])*(x[j]-xcl[i])/sum;
+ sumyy += clustcell[i][j]*(y[j]-ycl[i])*(y[j]-ycl[i])/sum;
+ sumxy += clustcell[i][j]*(x[j]-xcl[i])*(y[j]-ycl[i])/sum;
}
b = sumxx+sumyy;
c = sumxx*sumyy-sumxy*sumxy;
// To get the cell position in a cluster
testncl.Set(nclust+1);
- testindex.Set(ncell);
- cell[0][0] = ncell;
+ testindex.Set(ncell+1);
+ cell[0][0] = ncell + 1;
testncl[0] = cell[0][0];
Int_t ll = 0;
- for(Int_t ii=1; ii<=ncell; ii++)
+ for(Int_t ii = 1; ii <= ncell; ii++)
{
cell[0][ii]=ii;
- //clustcell[0][ii]=1.;
Int_t kk = cell[0][ii];
testindex[ll] = kk;
ll++;
}
// final assignments
- *(&xc + i) = xcl[i];
- *(&yc + i) = ycl[i];
- *(&zc + i) = str[i];
- *(&cells + i) = cln[i];
- *(&rcl+i) = r1;
- *(&rcs+i) = r2;
+ xc[i] = xcl[i];
+ yc[i] = ycl[i];
+ zc[i] = sum;
+ cells[i] = cln[i];
+ rcl[i] = r1;
+ rcs[i] = r2;
}
+ for(i = 0; i < kndim3; i++)
+ {
+ delete [] clustcell[i];
+ delete [] cell[i];
+ }
+ delete [] clustcell;
+ delete [] cell;
+ for(i = 0; i <kndim1 ; i++)
+ {
+ delete [] cluster[i];
+ }
+ delete [] cluster;
+ delete [] str;
+ delete [] str1;
+ delete [] xcl;
+ delete [] ycl;
+ delete [] cln;
}
// ------------------------------------------------------------------------ //
fCutoff = decut;
}
// ------------------------------------------------------------------------ //
+void AliPMDClusteringV2::SetClusteringParam(Int_t cluspar)
+{
+ fClusParam = cluspar;
+}
+// ------------------------------------------------------------------------ //
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff