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 array fEdepCell[kNDIMX][kNDIMY] and cluster information is in array
32 fClusters[5][5000]. integer fClno gives total number of clusters in the
35 fEdepCell, fClno and fClusters are the only global ( public ) variables.
36 Others are local ( private ) to the code.
37 At the moment, the data is read for whole detector ( all supermodules
38 and pmd as well as cpv. This will have to be modify later )
39 LAST UPDATE : October 23, 2002
40 -----------------------------------------------------------------------*/
42 #include <Riostream.h>
44 #include <TObjArray.h>
47 #include "AliPMDcluster.h"
48 #include "AliPMDClustering.h"
49 #include "AliPMDClusteringV2.h"
52 ClassImp(AliPMDClusteringV2)
54 const Double_t AliPMDClusteringV2::fgkSqroot3by2=0.8660254; // sqrt(3.)/2.
56 AliPMDClusteringV2::AliPMDClusteringV2():
60 for(int i = 0; i < kNDIMX; i++)
62 for(int j = 0; j < kNDIMY; j++)
64 fCoord[0][i][j] = i+j/2.;
65 fCoord[1][i][j] = fgkSqroot3by2*j;
70 // ------------------------------------------------------------------------ //
71 AliPMDClusteringV2::~AliPMDClusteringV2()
75 // ------------------------------------------------------------------------ //
76 void AliPMDClusteringV2::DoClust(Int_t idet, Int_t ismn, Double_t celladc[48][96], TObjArray *pmdcont)
78 // main function to call other necessary functions to do clustering
80 AliPMDcluster *pmdcl = 0;
82 Int_t i, i1, i2, j, nmx1, incr, id, jd;
83 Int_t celldataX[15], celldataY[15];
87 const float ktwobysqrt3 = 1.1547; // 2./sqrt(3.)
97 else if (ismn >= 12 && ismn <= 23)
103 for (Int_t i =0; i < kNDIMX; i++)
105 for (Int_t j =0; j < kNDIMY; j++)
112 for (id = 0; id < ndimXr; id++)
114 for (jd = 0; jd < ndimYr; jd++)
117 i=id+(ndimYr/2-1)-(jd/2);
121 fEdepCell[i][j] = celladc[jd][id];
123 else if (ismn >= 12 && ismn <= 23)
125 fEdepCell[i][j] = celladc[id][jd];
131 Order(); // order the data
132 cutoff = fCutoff; // cutoff used to discard cells having ener. dep.
140 if (fEdepCell[i1][i2] > 0.) {ave = ave + fEdepCell[i1][i2];}
141 if (fEdepCell[i1][i2] > cutoff ) nmx1 = nmx1 + 1;
143 // nmx1 --- number of cells having ener dep >= cutoff
145 AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
147 if (nmx1 == 0) nmx1 = 1;
150 AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
153 incr = CrClust(ave, cutoff, nmx1);
156 AliDebug(1,Form("Detector Plane = %d Serial Module No = %d Number of clusters = %d",idet, ismn, fClno));
158 for(i1=0; i1<=fClno; i1++)
160 Float_t cluXC = (Float_t) fClusters[0][i1];
161 Float_t cluYC = (Float_t) fClusters[1][i1];
162 Float_t cluADC = (Float_t) fClusters[2][i1];
163 Float_t cluCELLS = (Float_t) fClusters[3][i1];
164 Float_t sigmaX = (Float_t) fClusters[4][i1];
165 Float_t sigmaY = (Float_t) fClusters[5][i1];
166 Float_t cluY0 = ktwobysqrt3*cluYC;
167 Float_t cluX0 = cluXC - cluY0/2.;
169 // Cluster X centroid is back transformed
173 clusdata[0] = cluX0 - (24-1) + cluY0/2.;
175 else if (ismn >= 12 && ismn <= 23)
177 clusdata[0] = cluX0 - (48-1) + cluY0/2.;
181 clusdata[2] = cluADC;
182 clusdata[3] = cluCELLS;
183 clusdata[4] = sigmaX;
184 clusdata[5] = sigmaY;
187 // Cells associated with a cluster
189 for (Int_t ihit = 0; ihit < 15; ihit++)
191 celldataX[ihit] = 1; // dummy nos. -- will be changed
192 celldataY[ihit] = 1; // dummy nos. -- will be changed
195 pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY);
199 // ------------------------------------------------------------------------ //
200 void AliPMDClusteringV2::Order()
203 // sorts the ADC values from higher to lower
206 // matrix fEdepCell converted into
207 // one dimensional array dd. adum a place holder for double
208 int i, j, i1, i2, iord1[kNMX];
210 // ordering is stored in iord1, original array not ordered
212 // define arrays dd and iord1
213 for(i1=0; i1 < kNDIMX; i1++)
215 for(i2=0; i2 < kNDIMY; i2++)
219 dd[i] = fEdepCell[i1][i2];
222 // sort and store sorting information in iord1
224 TMath::Sort(kNMX,dd,iord1);
226 // store the sorted information in fIord for later use
227 for(i=0; i<kNMX; i++)
236 // ------------------------------------------------------------------------ //
237 Int_t AliPMDClusteringV2::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1)
239 // Does crude clustering
240 // Finds out only the big patch by just searching the
244 int i,j,k,id1,id2,icl, numcell;
245 int jd1,jd2, icell, cellcount;
247 static int neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
249 // neibx and neiby define ( incremental ) (i,j) for the neighbours of a
250 // cell. There are six neighbours.
251 // cellcount --- total number of cells having nonzero ener dep
252 // numcell --- number of cells in a given supercluster
253 // ofstream ofl0("cells_loc",ios::out);
254 // initialize fInfocl[2][kNDIMX][kNDIMY]
256 AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
258 for (j=0; j < kNDIMX; j++){
259 for(k=0; k < kNDIMY; k++){
260 fInfocl[0][j][k] = 0;
261 fInfocl[1][j][k] = 0;
264 for(i=0; i < kNMX; i++){
268 if(fEdepCell[id1][id2] <= cutoff){fInfocl[0][id1][id2]=-1;}
270 // ---------------------------------------------------------------
271 // crude clustering begins. Start with cell having largest adc
272 // count and loop over the cells in descending order of adc count
273 // ---------------------------------------------------------------
276 for(icell=0; icell <= nmx1; icell++){
279 if(fInfocl[0][id1][id2] == 0 ){
280 // ---------------------------------------------------------------
281 // icl -- cluster #, numcell -- # of cells in it, clust -- stores
282 // coordinates of the cells in a cluster, fInfocl[0][i1][i2] is 1 for
283 // primary and 2 for secondary cells,
284 // fInfocl[1][i1][i2] stores cluster #
285 // ---------------------------------------------------------------
288 cellcount = cellcount + 1;
289 fInfocl[0][id1][id2]=1;
290 fInfocl[1][id1][id2]=icl;
291 fInfcl[0][cellcount]=icl;
292 fInfcl[1][cellcount]=id1;
293 fInfcl[2][cellcount]=id2;
295 clust[0][numcell]=id1;
296 clust[1][numcell]=id2;
297 for(i=1; i<5000; i++)clust[0][i] = -1;
298 // ---------------------------------------------------------------
299 // check for adc count in neib. cells. If ne 0 put it in this clust
300 // ---------------------------------------------------------------
304 if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
305 fInfocl[0][jd1][jd2] == 0){
307 fInfocl[0][jd1][jd2]=2;
308 fInfocl[1][jd1][jd2]=icl;
309 clust[0][numcell]=jd1;
310 clust[1][numcell]=jd2;
311 cellcount=cellcount+1;
312 fInfcl[0][cellcount]=icl;
313 fInfcl[1][cellcount]=jd1;
314 fInfcl[2][cellcount]=jd2;
317 // ---------------------------------------------------------------
318 // check adc count for neighbour's neighbours recursively and
319 // if nonzero, add these to the cluster.
320 // ---------------------------------------------------------------
321 for(i=1;i < 5000;i++){
322 if(clust[0][i] != -1){
328 if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
329 fInfocl[0][jd1][jd2] == 0 ){
330 fInfocl[0][jd1][jd2] = 2;
331 fInfocl[1][jd1][jd2] = icl;
332 numcell = numcell + 1;
333 clust[0][numcell] = jd1;
334 clust[1][numcell] = jd2;
335 cellcount = cellcount+1;
336 fInfcl[0][cellcount] = icl;
337 fInfcl[1][cellcount] = jd1;
338 fInfcl[2][cellcount] = jd2;
345 // for(icell=0; icell<=cellcount; icell++){
346 // ofl0 << fInfcl[0][icell] << " " << fInfcl[1][icell] << " " <<
347 // fInfcl[2][icell] << endl;
351 // ------------------------------------------------------------------------ //
352 void AliPMDClusteringV2::RefClust(Int_t incr)
354 // Does the refining of clusters
355 // Takes the big patch and does gaussian fitting and
356 // finds out the more refined clusters
359 const Int_t kndim = 4500;
361 int i, j, k, i1, i2, id, icl, itest;
364 int ncl[kndim], iord[kndim];
366 double x1, y1, z1, x2, y2, z2;
369 double x[kndim], y[kndim], z[kndim];
370 double xc[kndim], yc[kndim], zc[kndim], cells[kndim];
371 double rcl[kndim], rcs[kndim];
373 // fClno counts the final clusters
374 // nsupcl = # of superclusters; ncl[i]= # of cells in supercluster i
375 // x, y and z store (x,y) coordinates of and energy deposited in a cell
376 // xc, yc store (x,y) coordinates of the cluster center
377 // zc stores the energy deposited in a cluster
378 // rc is cluster radius
379 // finally the cluster information is put in 2-dimensional array clusters
380 // ofstream ofl1("checking.5",ios::app);
384 for(i=0; i<4500; i++){ncl[i]=-1;}
385 for(i=0; i<incr; i++){
386 if(fInfcl[0][i] != nsupcl){ nsupcl=nsupcl+1; }
388 AliWarning("RefClust: Too many superclusters!");
392 ncl[nsupcl]=ncl[nsupcl]+1;
395 AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
400 for(i=0; i<nsupcl; i++){
404 // one cell super-clusters --> single cluster
405 // cluster center at the centyer of the cell
406 // cluster radius = half cell dimension
408 AliWarning("RefClust: Too many clusters! more than 5000");
414 fClusters[0][fClno] = fCoord[0][i1][i2];
415 fClusters[1][fClno] = fCoord[1][i1][i2];
416 fClusters[2][fClno] = fEdepCell[i1][i2];
417 fClusters[3][fClno] = 1.;
418 fClusters[4][fClno] = 0.0;
419 fClusters[5][fClno] = 0.0;
420 //ofl1 << icl << " " << fCoord[0][i1][i2] << " " << fCoord[1][i1][i2] <<
421 //" " << fEdepCell[i1][i2] << " " << fClusters[3][fClno] <<endl;
422 }else if(ncl[i] == 1){
423 // two cell super-cluster --> single cluster
424 // cluster center is at ener. dep.-weighted mean of two cells
425 // cluster radius == half cell dimension
429 AliWarning("RefClust: Too many clusters! more than 5000");
435 x1 = fCoord[0][i1][i2];
436 y1 = fCoord[1][i1][i2];
437 z1 = fEdepCell[i1][i2];
442 x2 = fCoord[0][i1][i2];
443 y2 = fCoord[1][i1][i2];
444 z2 = fEdepCell[i1][i2];
446 fClusters[0][fClno] = (x1*z1+x2*z2)/(z1+z2);
447 fClusters[1][fClno] = (y1*z1+y2*z2)/(z1+z2);
448 fClusters[2][fClno] = z1+z2;
449 fClusters[3][fClno] = 2.;
450 fClusters[4][fClno] = sqrt(z1*z2)/(z1+z2);
451 fClusters[5][fClno] = 0; // sigma large nonzero, sigma small zero
453 //ofl1 << icl << " " << fClusters[0][fClno] << " " << fClusters[1][fClno]
454 // << " " << fClusters[2][fClno] << " " <<fClusters[3][fClno] <<endl;
459 // super-cluster of more than two cells - broken up into smaller
460 // clusters gaussian centers computed. (peaks separated by > 1 cell)
461 // Begin from cell having largest energy deposited This is first
463 // *****************************************************************
464 // NOTE --- POSSIBLE MODIFICATION: ONE MAY NOT BREAKING SUPERCLUSTERS
465 // IF NO. OF CELLS IS NOT TOO LARGE ( SAY 5 OR 6 )
466 // SINCE WE EXPECT THE SUPERCLUSTER
467 // TO BE A SINGLE CLUSTER
468 //*******************************************************************
472 x[0] = fCoord[0][i1][i2];
473 y[0] = fCoord[1][i1][i2];
474 z[0] = fEdepCell[i1][i2];
476 for(j=1;j<=ncl[i];j++){
482 x[j] = fCoord[0][i1][i2];
483 y[j] = fCoord[1][i1][i2];
484 z[j] = fEdepCell[i1][i2];
486 // arranging cells within supercluster in decreasing order
487 for(j=1;j<=ncl[i];j++)
491 for(i1=0; i1<j; i1++)
493 if(itest == 0 && z[iord[i1]] < z[ihld])
496 for(i2=j-1;i2>=i1;i2--)
498 iord[i2+1] = iord[i2];
505 // compute the number of clusters and their centers ( first
507 // centers must be separated by cells having smaller ener. dep.
508 // neighbouring centers should be either strong or well-separated
513 for(j=1;j<=ncl[i];j++){
520 rr = Distance(x1,y1,x2,y2);
521 //***************************************************************
522 // finetuning cluster splitting
523 // the numbers zc/4 and zc/10 may need to be changed.
524 // Also one may need to add one more layer because our
525 // cells are smaller in absolute scale
526 //****************************************************************
529 if( rr >= 1.1 && rr < 1.8 && z[iord[j]] > zc[k]/4.)
531 if( rr >= 1.8 && rr < 2.1 && z[iord[j]] > zc[k]/10.)
533 if( rr >= 2.1)itest++;
543 ClustDetails(ncl[i], ig, x[0], y[0] ,z[0], xc[0], yc[0], zc[0],
544 rcl[0], rcs[0], cells[0]);
552 AliWarning("RefClust: Too many clusters! more than 5000");
556 fClusters[0][fClno] = xc[j];
557 fClusters[1][fClno] = yc[j];
558 fClusters[2][fClno] = zc[j];
559 fClusters[4][fClno] = rcl[j];
560 fClusters[5][fClno] = rcs[j];
563 fClusters[3][fClno] = ncl[i];
567 fClusters[3][fClno] = cells[j];
577 // ------------------------------------------------------------------------ //
579 void AliPMDClusteringV2::ClustDetails(Int_t ncell, Int_t nclust,
580 Double_t &x, Double_t &y, Double_t &z,
581 Double_t &xc, Double_t &yc, Double_t &zc,
582 Double_t &rcl, Double_t &rcs,
588 const Int_t kndim1 = 4500;
589 const Int_t kndim2 = 10;
590 const Int_t kndim3 = 100;
593 int cluster[kndim1][kndim2];
595 double x1, y1, x2, y2, rr;
596 double sumx, sumy, sumxy, sumxx;
597 double sum, sum1, sumyy;
600 double xx[kndim1], yy[kndim1], zz[kndim1];
601 double xxc[kndim1], yyc[kndim1];
606 double xcl[kndim1], ycl[kndim1], cln[kndim1];
607 double clustcell[kndim1][kndim3];
609 for(i=0; i<=nclust; i++){
615 for(i=0; i<=ncell; i++){
621 for(i=0; i<4500; i++){
622 for(j=0; j<100; j++){
636 // more than one cluster
637 // checking cells shared between several clusters.
638 // First check if the cell is within
639 // one cell unit ( nearest neighbour). Else,
640 // if it is within 1.74 cell units ( next nearest )
641 // Else if it is upto 2 cell units etc.
643 for (i=0; i<=ncell; i++){
647 // distance <= 1 cell unit
648 for(j=0; j<=nclust; j++)
652 rr = Distance(x1, y1, x2, y2);
660 // next nearest neighbour
661 if(cluster[i][0] == 0)
663 for(j=0; j<=nclust; j++)
667 rr = Distance(x1, y1, x2, y2);
676 // next-to-next nearest neighbour
677 if(cluster[i][0] == 0)
679 for(j=0; j<=nclust; j++)
683 rr = Distance(x1, y1, x2, y2);
693 if(cluster[i][0] == 0)
695 for(j=0; j<=nclust; j++)
699 rr = Distance(x1, y1, x2, y2);
711 // computing cluster strength. Some cells are shared.
712 for(i=0; i<=ncell; i++){
713 if(cluster[i][0] != 0){
715 for(j=1; j<=i1; j++){
717 str[i2] = str[i2]+zz[i]/i1;
724 for(i=0; i<=ncell; i++)
726 if(cluster[i][0] != 0)
732 sum=sum+str[cluster[i][j]];
738 str1[i2] = str1[i2] + zz[i]*str[i2]/sum;
739 clustcell[i2][i] = zz[i]*str[i2]/sum;
745 for(j=0; j<=nclust; j++)
752 for(i=0; i<=nclust; i++){
757 for(j=0; j<=ncell; j++){
758 if(clustcell[i][j] != 0){
759 sumx = sumx+clustcell[i][j]*xx[j];
760 sumy = sumy+clustcell[i][j]*yy[j];
761 sum = sum+clustcell[i][j];
762 sum1 = sum1+clustcell[i][j]/zz[j];
765 //***** xcl and ycl are cluster centroid positions ( center of gravity )
773 for(j=0; j<=ncell; j++){
774 sumxx = sumxx+clustcell[i][j]*(xx[j]-xcl[i])*(xx[j]-xcl[i])/sum;
775 sumyy = sumyy+clustcell[i][j]*(yy[j]-ycl[i])*(yy[j]-ycl[i])/sum;
776 sumxy = sumxy+clustcell[i][j]*(xx[j]-xcl[i])*(yy[j]-ycl[i])/sum;
779 c = sumxx*sumyy-sumxy*sumxy;
780 // ******************r1 and r2 are major and minor axes ( r1 > r2 ).
781 r1 = b/2.+sqrt(b*b/4.-c);
782 r2 = b/2.-sqrt(b*b/4.-c);
783 // final assignments to proper external variables
787 *(&cells + i) = cln[i];
797 for(j=0; j<=ncell; j++){
798 sumx = sumx+zz[j]*xx[j];
799 sumy = sumy+zz[j]*yy[j];
809 for(j=0; j<=ncell; j++){
810 sumxx = sumxx+clustcell[i][j]*(xx[j]-xcl[i])*(xx[j]-xcl[i])/sum;
811 sumyy = sumyy+clustcell[i][j]*(yy[j]-ycl[i])*(yy[j]-ycl[i])/sum;
812 sumxy = sumxy+clustcell[i][j]*(xx[j]-xcl[i])*(yy[j]-ycl[i])/sum;
815 c = sumxx*sumyy-sumxy*sumxy;
816 r1 = b/2.+sqrt(b*b/4.-c);
817 r2 = b/2.-sqrt(b*b/4.-c);
822 *(&cells + i) = cln[i];
828 // ------------------------------------------------------------------------ //
829 Double_t AliPMDClusteringV2::Distance(Double_t x1, Double_t y1,
830 Double_t x2, Double_t y2)
832 return sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
834 // ------------------------------------------------------------------------ //
835 void AliPMDClusteringV2::SetEdepCut(Float_t decut)
839 // ------------------------------------------------------------------------ //