1 /***************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 //-----------------------------------------------------//
20 // Source File : PMDClusteringV2.cxx //
22 // clustering code for alice pmd //
24 //-----------------------------------------------------//
26 /* --------------------------------------------------------------------
27 Code developed by S. C. Phatak, Institute of Physics,
28 Bhubaneswar 751 005 ( phatak@iopb.res.in ) Given the energy deposited
29 ( or ADC value ) in each cell of supermodule ( pmd or cpv ), the code
30 builds up superclusters and breaks them into clusters. The input is
31 in TObjarray and cluster information is in TObjArray.
32 integer clno gives total number of clusters in the supermodule.
33 fClusters is the global ( public ) variables.
34 Others are local ( private ) to the code.
35 At the moment, the data is read for whole detector ( all supermodules
36 and pmd as well as cpv. This will have to be modify later )
37 LAST UPDATE : October 23, 2002
38 -----------------------------------------------------------------------*/
40 #include <Riostream.h>
42 #include <TObjArray.h>
45 #include "AliPMDcludata.h"
46 #include "AliPMDcluster.h"
47 #include "AliPMDClustering.h"
48 #include "AliPMDClusteringV2.h"
51 ClassImp(AliPMDClusteringV2)
53 const Double_t AliPMDClusteringV2::fgkSqroot3by2=0.8660254; // sqrt(3.)/2.
55 AliPMDClusteringV2::AliPMDClusteringV2():
56 fPMDclucont(new TObjArray()),
59 for(int i = 0; i < kNDIMX; i++)
61 for(int j = 0; j < kNDIMY; j++)
63 fCoord[0][i][j] = i+j/2.;
64 fCoord[1][i][j] = fgkSqroot3by2*j;
68 // ------------------------------------------------------------------------ //
71 AliPMDClusteringV2::AliPMDClusteringV2(const AliPMDClusteringV2& pmdclv2):
72 AliPMDClustering(pmdclv2),
77 AliError("Copy constructor not allowed ");
80 // ------------------------------------------------------------------------ //
81 AliPMDClusteringV2 &AliPMDClusteringV2::operator=(const AliPMDClusteringV2& /*pmdclv2*/)
84 AliError("Assignment operator not allowed ");
87 // ------------------------------------------------------------------------ //
88 AliPMDClusteringV2::~AliPMDClusteringV2()
92 // ------------------------------------------------------------------------ //
94 void AliPMDClusteringV2::DoClust(Int_t idet, Int_t ismn,
95 Double_t celladc[48][96], TObjArray *pmdcont)
97 // main function to call other necessary functions to do clustering
99 AliPMDcluster *pmdcl = 0;
101 const float ktwobysqrt3 = 1.1547; // 2./sqrt(3.)
103 Int_t i, j, nmx1, incr, id, jd;
106 Int_t celldataX[15], celldataY[15];
108 Double_t cutoff, ave;
109 Double_t edepcell[kNMX];
117 else if (ismn >= 12 && ismn <= 23)
123 for (Int_t i =0; i < kNMX; i++)
128 for (id = 0; id < ndimXr; id++)
130 for (jd = 0; jd < ndimYr; jd++)
133 i = id + (ndimYr/2-1) - (jd/2);
134 Int_t ij = i + j*kNDIMX;
137 edepcell[ij] = celladc[jd][id];
139 else if (ismn >= 12 && ismn <= 23)
141 edepcell[ij] = celladc[id][jd];
148 TMath::Sort(kNMX,edepcell,iord1);// order the data
149 cutoff = fCutoff; // cutoff used to discard cells having ener. dep.
153 for(i = 0;i < kNMX; i++)
159 if(edepcell[i] > cutoff )
165 AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
173 AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
176 incr = CrClust(ave, cutoff, nmx1,iord1, edepcell);
177 RefClust(incr,edepcell );
179 Int_t nentries1 = fPMDclucont->GetEntries();
180 AliDebug(1,Form("Detector Plane = %d Serial Module No = %d Number of clusters = %d",idet, ismn, nentries1));
181 AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
182 for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
184 AliPMDcludata *pmdcludata =
185 (AliPMDcludata*)fPMDclucont->UncheckedAt(ient1);
186 Float_t cluXC = pmdcludata->GetClusX();
187 Float_t cluYC = pmdcludata->GetClusY();
188 Float_t cluADC = pmdcludata->GetClusADC();
189 Float_t cluCELLS = pmdcludata->GetClusCells();
190 Float_t cluSIGX = pmdcludata->GetClusSigmaX();
191 Float_t cluSIGY = pmdcludata->GetClusSigmaY();
193 Float_t cluY0 = ktwobysqrt3*cluYC;
194 Float_t cluX0 = cluXC - cluY0/2.;
197 // Cluster X centroid is back transformed
201 clusdata[0] = cluX0 - (24-1) + cluY0/2.;
203 else if (ismn >= 12 && ismn <= 23)
205 clusdata[0] = cluX0 - (48-1) + cluY0/2.;
209 clusdata[2] = cluADC;
210 clusdata[3] = cluCELLS;
211 clusdata[4] = cluSIGX;
212 clusdata[5] = cluSIGY;
214 // Cells associated with a cluster
216 for (Int_t ihit = 0; ihit < 15; ihit++)
218 Int_t dummyXY = pmdcludata->GetCellXY(ihit);
220 Int_t celldumY = dummyXY%10000;
221 Int_t celldumX = dummyXY/10000;
222 Float_t cellY = (Float_t) (ktwobysqrt3*celldumY/10);
223 Float_t cellX = (Float_t) (celldumX/10 - (celldumY/2.)/10);
226 // Cell X centroid is back transformed
230 celldataX[ihit] = (Int_t) (cellX - (24-1) + cellY/2.);
232 else if (ismn >= 12 && ismn <= 23)
234 celldataX[ihit] = (Int_t) (cellX - (48-1) + cellY/2.);
236 celldataY[ihit] = (Int_t) cellY;
239 pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY);
242 fPMDclucont->Clear();
244 // ------------------------------------------------------------------------ //
245 Int_t AliPMDClusteringV2::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
246 Int_t iord1[], Double_t edepcell[])
248 // Does crude clustering
249 // Finds out only the big patch by just searching the
253 Int_t i,j,k,id1,id2,icl, numcell;
254 Int_t jd1,jd2, icell, cellcount;
255 Int_t clust[2][5000];
256 static Int_t neibx[6] = {1,0,-1,-1,0,1}, neiby[6] = {0,1,1,0,-1,-1};
258 // neibx and neiby define ( incremental ) (i,j) for the neighbours of a
259 // cell. There are six neighbours.
260 // cellcount --- total number of cells having nonzero ener dep
261 // numcell --- number of cells in a given supercluster
263 AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
265 for (j=0; j < kNDIMX; j++)
267 for(k=0; k < kNDIMY; k++)
269 fInfocl[0][j][k] = 0;
270 fInfocl[1][j][k] = 0;
274 for(i=0; i < kNMX; i++)
282 if(edepcell[j] <= cutoff)
284 fInfocl[0][id1][id2] = -1;
287 // ---------------------------------------------------------------
288 // crude clustering begins. Start with cell having largest adc
289 // count and loop over the cells in descending order of adc count
290 // ---------------------------------------------------------------
293 for(icell=0; icell <= nmx1; icell++)
298 if(fInfocl[0][id1][id2] == 0 )
300 // ---------------------------------------------------------------
301 // icl -- cluster #, numcell -- # of cells in it, clust -- stores
302 // coordinates of the cells in a cluster, fInfocl[0][i1][i2] is 1 for
303 // primary and 2 for secondary cells,
304 // fInfocl[1][i1][i2] stores cluster #
305 // ---------------------------------------------------------------
309 fInfocl[0][id1][id2] = 1;
310 fInfocl[1][id1][id2] = icl;
311 fInfcl[0][cellcount] = icl;
312 fInfcl[1][cellcount] = id1;
313 fInfcl[2][cellcount] = id2;
315 clust[0][numcell] = id1;
316 clust[1][numcell] = id2;
317 for(i = 1; i < 5000; i++)
321 // ---------------------------------------------------------------
322 // check for adc count in neib. cells. If ne 0 put it in this clust
323 // ---------------------------------------------------------------
324 for(i = 0; i < 6; i++)
326 jd1 = id1 + neibx[i];
327 jd2 = id2 + neiby[i];
328 if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
329 fInfocl[0][jd1][jd2] == 0)
332 fInfocl[0][jd1][jd2] = 2;
333 fInfocl[1][jd1][jd2] = icl;
334 clust[0][numcell] = jd1;
335 clust[1][numcell] = jd2;
337 fInfcl[0][cellcount] = icl;
338 fInfcl[1][cellcount] = jd1;
339 fInfcl[2][cellcount] = jd2;
342 // ---------------------------------------------------------------
343 // check adc count for neighbour's neighbours recursively and
344 // if nonzero, add these to the cluster.
345 // ---------------------------------------------------------------
346 for(i = 1;i < 5000; i++)
348 if(clust[0][i] != -1)
352 for(j = 0; j < 6 ; j++)
354 jd1 = id1 + neibx[j];
355 jd2 = id2 + neiby[j];
356 if( (jd1 >= 0 && jd1 < kNDIMX) &&
357 (jd2 >= 0 && jd2 < kNDIMY)
358 && fInfocl[0][jd1][jd2] == 0 )
360 fInfocl[0][jd1][jd2] = 2;
361 fInfocl[1][jd1][jd2] = icl;
363 clust[0][numcell] = jd1;
364 clust[1][numcell] = jd2;
366 fInfcl[0][cellcount] = icl;
367 fInfcl[1][cellcount] = jd1;
368 fInfcl[2][cellcount] = jd2;
377 // ------------------------------------------------------------------------ //
378 void AliPMDClusteringV2::RefClust(Int_t incr, Double_t edepcell[])
380 // Does the refining of clusters
381 // Takes the big patch and does gaussian fitting and
382 // finds out the more refined clusters
384 AliPMDcludata *pmdcludata = 0;
387 const Int_t kndim = 4500;
388 Int_t i, j, k, i1, i2, id, icl, itest, ihld;
389 Int_t ig, nsupcl, clno, clX,clY;
391 Int_t ncl[kndim], iord[kndim];
393 Double_t x1, y1, z1, x2, y2, z2, rr;
394 Double_t x[kndim], y[kndim], z[kndim];
395 Double_t xc[kndim], yc[kndim], zc[kndim], cells[kndim];
396 Double_t rcl[kndim], rcs[kndim];
398 for(Int_t kk = 0; kk < 15; kk++)
400 if( kk < 6 )clusdata[kk] = 0.;
403 // nsupcl = # of superclusters; ncl[i]= # of cells in supercluster i
404 // x, y and z store (x,y) coordinates of and energy deposited in a cell
405 // xc, yc store (x,y) coordinates of the cluster center
406 // zc stores the energy deposited in a cluster, rc is cluster radius
410 for(i = 0; i < 4500; i++)
414 for(i = 0; i < incr; i++)
416 if(fInfcl[0][i] != nsupcl)
422 AliWarning("RefClust: Too many superclusters!");
429 AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
434 for(i = 0; i < nsupcl; i++)
440 // one cell super-clusters --> single cluster
441 // cluster center at the centyer of the cell
442 // cluster radius = half cell dimension
445 AliWarning("RefClust: Too many clusters! more than 5000");
451 Int_t i12 = i1 + i2*kNDIMX;
452 clusdata[0] = fCoord[0][i1][i2];
453 clusdata[1] = fCoord[1][i1][i2];
454 clusdata[2] = edepcell[i12];
460 clX = (Int_t) fCoord[0][i1][i2]*10;
461 clY = (Int_t) fCoord[1][i1][i2]*10;
462 clxy[0] = clX*10000 + clY;
463 for(Int_t icltr = 1; icltr < 15; icltr++)
467 pmdcludata = new AliPMDcludata(clusdata,clxy);
468 fPMDclucont->Add(pmdcludata);
474 // two cell super-cluster --> single cluster
475 // cluster center is at ener. dep.-weighted mean of two cells
476 // cluster radius == half cell dimension
481 AliWarning("RefClust: Too many clusters! more than 5000");
487 Int_t i12 = i1 + i2*kNDIMX;
489 x1 = fCoord[0][i1][i2];
490 y1 = fCoord[1][i1][i2];
496 i12 = i1 + i2*kNDIMX;
498 x2 = fCoord[0][i1][i2];
499 y2 = fCoord[1][i1][i2];
502 clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
503 clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
506 clusdata[4] = (TMath::Sqrt(z1*z2))/(z1+z2);
511 clxy[0] = clX*10000 + clY;
515 clxy[1] = clX*10000 + clY;
517 for(Int_t icltr = 2; icltr < 15; icltr++)
521 pmdcludata = new AliPMDcludata(clusdata, clxy);
522 fPMDclucont->Add(pmdcludata);
527 // super-cluster of more than two cells - broken up into smaller
528 // clusters gaussian centers computed. (peaks separated by > 1 cell)
529 // Begin from cell having largest energy deposited This is first
531 // *****************************************************************
532 // NOTE --- POSSIBLE MODIFICATION: ONE MAY NOT BREAKING SUPERCLUSTERS
533 // IF NO. OF CELLS IS NOT TOO LARGE ( SAY 5 OR 6 )
534 // SINCE WE EXPECT THE SUPERCLUSTER
535 // TO BE A SINGLE CLUSTER
536 //*******************************************************************
540 Int_t i12 = i1 + i2*kNDIMX;
542 x[0] = fCoord[0][i1][i2];
543 y[0] = fCoord[1][i1][i2];
544 z[0] = edepcell[i12];
547 for(j = 1; j <= ncl[i]; j++)
553 Int_t i12 = i1 + i2*kNDIMX;
555 x[j] = fCoord[0][i1][i2];
556 y[j] = fCoord[1][i1][i2];
557 z[j] = edepcell[i12];
560 // arranging cells within supercluster in decreasing order
561 for(j = 1; j <= ncl[i];j++)
565 for(i1 = 0; i1 < j; i1++)
567 if(itest == 0 && z[iord[i1]] < z[ihld])
570 for(i2 = j-1;i2 >= i1;i2--)
572 iord[i2+1] = iord[i2];
580 // compute the number of clusters and their centers ( first
582 // centers must be separated by cells having smaller ener. dep.
583 // neighbouring centers should be either strong or well-separated
588 for(j = 1; j <= ncl[i]; j++)
593 for(k = 0; k <= ig; k++)
597 rr = Distance(x1,y1,x2,y2);
598 //************************************************************
599 // finetuning cluster splitting
600 // the numbers zc/4 and zc/10 may need to be changed.
601 // Also one may need to add one more layer because our
602 // cells are smaller in absolute scale
603 //************************************************************
606 if( rr >= 1.1 && rr < 1.8 && z[iord[j]] > zc[k]/4.) itest++;
607 if( rr >= 1.8 && rr < 2.1 && z[iord[j]] > zc[k]/10.) itest++;
608 if( rr >= 2.1)itest++;
619 ClustDetails(ncl[i], ig, x[0], y[0] ,z[0], xc[0], yc[0], zc[0],
620 rcl[0], rcs[0], cells[0], testncl, testindex);
623 for(j = 0; j <= ig; j++)
628 AliWarning("RefClust: Too many clusters! more than 5000");
634 clusdata[4] = rcl[j];
635 clusdata[5] = rcs[j];
638 clusdata[3] = ncl[i];
642 clusdata[3] = cells[j];
645 Int_t ncellcls = testncl[j];
648 for(Int_t kk = 1; kk <= ncellcls; kk++)
650 Int_t ll = testindex[pp];
651 clX = (Int_t) x[ll]*10;
652 clY = (Int_t) y[ll]*10;
653 clxy[kk-1] = (Int_t) clX*10000 + clY ;
656 for(Int_t icltr = ncellcls ; icltr < 15; icltr++)
661 pmdcludata = new AliPMDcludata(clusdata, clxy);
662 fPMDclucont->Add(pmdcludata);
669 // ------------------------------------------------------------------------ //
670 void AliPMDClusteringV2::ClustDetails(Int_t ncell, Int_t nclust, Double_t &x,
671 Double_t &y, Double_t &z, Double_t &xc,
672 Double_t &yc, Double_t &zc,
673 Double_t &rcl, Double_t &rcs,
674 Double_t &cells, TArrayI &testncl,
680 const Int_t kndim1 = 2000;
681 const Int_t kndim2 = 10;
682 const Int_t kndim3 = 400;
684 Int_t i, j, k, i1, i2;
685 Int_t cluster[kndim1][kndim2], cell[kndim1][kndim3];
687 Double_t x1, y1, x2, y2, rr, b, c, r1, r2;
688 Double_t sumx, sumy, sumxy, sumxx, sum, sum1, sumyy;
689 Double_t xx[kndim1], yy[kndim1], zz[kndim1];
690 Double_t xxc[kndim1], yyc[kndim1],clustcell[kndim3][kndim1];
691 Double_t str[kndim1], str1[kndim1],xcl[kndim1], ycl[kndim1], cln[kndim1];
693 for(i = 0; i <= nclust; i++)
700 for(i = 0; i <= ncell; i++)
709 for(i = 0; i < kndim1; i++)
711 for(j = 0; j < kndim2; j++)
716 for(j = 0; j < kndim3; j++)
719 clustcell[j][i] = 0.;
726 // more than one cluster
727 // checking cells shared between several clusters.
728 // First check if the cell is within
729 // one cell unit ( nearest neighbour). Else,
730 // if it is within 1.74 cell units ( next nearest )
731 // Else if it is upto 2 cell units etc.
733 for (i = 0; i <= ncell; i++)
738 // distance <= 1 cell unit
739 for(j = 0; j <= nclust; j++)
743 rr = Distance(x1, y1, x2, y2);
751 // next nearest neighbour
752 if(cluster[i][0] == 0)
754 for(j=0; j<=nclust; j++)
758 rr = Distance(x1, y1, x2, y2);
759 if(rr <= TMath::Sqrt(3.))
767 // next-to-next nearest neighbour
768 if(cluster[i][0] == 0)
770 for(j=0; j<=nclust; j++)
774 rr = Distance(x1, y1, x2, y2);
784 if(cluster[i][0] == 0)
786 for(j = 0; j <= nclust; j++)
790 rr = Distance(x1, y1, x2, y2);
801 // computing cluster strength. Some cells are shared.
802 for(i = 0; i <= ncell; i++)
804 if(cluster[i][0] != 0)
807 for(j = 1; j <= i1; j++)
815 for(k = 0; k < 5; k++)
817 for(i = 0; i <= ncell; i++)
819 if(cluster[i][0] != 0)
823 for(j = 1; j <= i1; j++)
825 sum += str[cluster[i][j]];
828 for(j = 1; j <= i1; j++)
831 str1[i2] += zz[i]*str[i2]/sum;
832 clustcell[i2][i] = zz[i]*str[i2]/sum;
838 for(j = 0; j <= nclust; j++)
845 for(i = 0; i <= nclust; i++)
851 for(j = 0; j <= ncell; j++)
853 if(clustcell[i][j] != 0)
855 sumx += clustcell[i][j]*xx[j];
856 sumy += clustcell[i][j]*yy[j];
857 sum += clustcell[i][j];
858 sum1 += clustcell[i][j]/zz[j];
861 //** xcl and ycl are cluster centroid positions ( center of gravity )
869 for(j = 0; j <= ncell; j++)
871 sumxx += clustcell[i][j]*(xx[j]-xcl[i])*(xx[j]-xcl[i])/sum;
872 sumyy += clustcell[i][j]*(yy[j]-ycl[i])*(yy[j]-ycl[i])/sum;
873 sumxy += clustcell[i][j]*(xx[j]-xcl[i])*(yy[j]-ycl[i])/sum;
876 c = sumxx*sumyy-sumxy*sumxy;
877 // ******************r1 and r2 are major and minor axes ( r1 > r2 ).
878 r1 = b/2.+TMath::Sqrt(b*b/4.-c);
879 r2 = b/2.-TMath::Sqrt(b*b/4.-c);
880 // final assignments to proper external variables
884 *(&cells + i) = cln[i];
889 //To get the cell position in a cluster
891 for(Int_t ii=0; ii<= ncell; ii++)
893 Int_t jj = cluster[ii][0];
894 for(Int_t kk=1; kk<= jj; kk++)
896 Int_t ll = cluster[ii][kk];
898 cell[ll][cell[ll][0]] = ii;
902 testncl.Set(nclust+1);
905 for(Int_t ii=0; ii <= nclust; ii++)
907 testncl[ii] = cell[ii][0];
908 counter += testncl[ii];
910 testindex.Set(counter);
912 for(Int_t ii=0; ii<= nclust; ii++)
914 for(Int_t jj = 1; jj<= testncl[ii]; jj++)
916 Int_t kk = cell[ii][jj];
930 for(j = 0; j <= ncell; j++)
943 for(j = 0; j <= ncell; j++)
945 sumxx += clustcell[i][j]*(xx[j]-xcl[i])*(xx[j]-xcl[i])/sum;
946 sumyy += clustcell[i][j]*(yy[j]-ycl[i])*(yy[j]-ycl[i])/sum;
947 sumxy += clustcell[i][j]*(xx[j]-xcl[i])*(yy[j]-ycl[i])/sum;
950 c = sumxx*sumyy-sumxy*sumxy;
951 r1 = b/2.+ TMath::Sqrt(b*b/4.-c);
952 r2 = b/2.- TMath::Sqrt(b*b/4.-c);
954 // To get the cell position in a cluster
955 testncl.Set(nclust+1);
956 testindex.Set(ncell);
958 testncl[0] = cell[0][0];
960 for(Int_t ii=1; ii<=ncell; ii++)
963 //clustcell[0][ii]=1.;
964 Int_t kk = cell[0][ii];
972 *(&cells + i) = cln[i];
978 // ------------------------------------------------------------------------ //
979 Double_t AliPMDClusteringV2::Distance(Double_t x1, Double_t y1,
980 Double_t x2, Double_t y2)
982 return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
984 // ------------------------------------------------------------------------ //
985 void AliPMDClusteringV2::SetEdepCut(Float_t decut)
989 // ------------------------------------------------------------------------ //