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 : PMDClusteringV1.cxx, Version 00 //
22 // Date : September 26 2002 //
24 // clustering code for alice pmd //
26 //-----------------------------------------------------//
28 /* --------------------------------------------------------------------
29 Code developed by S. C. Phatak, Institute of Physics,
30 Bhubaneswar 751 005 ( phatak@iopb.res.in ) Given the energy deposited
31 ( or ADC value ) in each cell of supermodule ( pmd or cpv ), the code
32 builds up superclusters and breaks them into clusters. The input is
33 in array edepcell[kNMX] and cluster information is in a
34 TObjarray. Integer clno gives total number of clusters in the
37 fClusters is the only global ( public ) variables.
38 Others are local ( private ) to the code.
39 At the moment, the data is read for whole detector ( all supermodules
40 and pmd as well as cpv. This will have to be modify later )
41 LAST UPDATE : October 23, 2002
42 -----------------------------------------------------------------------*/
44 #include <Riostream.h>
47 #include <TObjArray.h>
51 #include "AliPMDcludata.h"
52 #include "AliPMDcluster.h"
53 #include "AliPMDClustering.h"
54 #include "AliPMDClusteringV1.h"
58 ClassImp(AliPMDClusteringV1)
60 const Double_t AliPMDClusteringV1::fgkSqroot3by2=0.8660254; // sqrt(3.)/2.
62 AliPMDClusteringV1::AliPMDClusteringV1():
63 fPMDclucont(new TObjArray()),
66 for(Int_t i = 0; i < kNDIMX; i++)
68 for(Int_t j = 0; j < kNDIMY; j++)
70 fCoord[0][i][j] = i+j/2.;
71 fCoord[1][i][j] = fgkSqroot3by2*j;
75 // ------------------------------------------------------------------------ //
76 AliPMDClusteringV1::AliPMDClusteringV1(const AliPMDClusteringV1& pmdclv1):
77 AliPMDClustering(pmdclv1),
82 AliError("Copy constructor not allowed ");
85 // ------------------------------------------------------------------------ //
86 AliPMDClusteringV1 &AliPMDClusteringV1::operator=(const AliPMDClusteringV1& /*pmdclv1*/)
89 AliError("Assignment operator not allowed ");
92 // ------------------------------------------------------------------------ //
93 AliPMDClusteringV1::~AliPMDClusteringV1()
97 // ------------------------------------------------------------------------ //
98 void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn,
99 Int_t celltrack[48][96],
100 Int_t cellpid[48][96],
101 Double_t celladc[48][96],
104 // main function to call other necessary functions to do clustering
107 AliPMDcluster *pmdcl = 0;
109 const float ktwobysqrt3 = 1.1547; // 2./sqrt(3.)
110 const Int_t kNmaxCell = 19; // # of cells surrounding a cluster center
112 Int_t i, j, nmx1, incr, id, jd;
113 Int_t celldataX[kNmaxCell], celldataY[kNmaxCell];
114 Int_t celldataTr[kNmaxCell], celldataPid[kNmaxCell];
115 Float_t celldataAdc[kNmaxCell];
117 Double_t cutoff, ave;
118 Double_t edepcell[kNMX];
121 Double_t cellenergy[11424];
123 // ndimXr and ndimYr are different because of different module size
133 else if (ismn >= 12 && ismn <= 23)
139 for (i =0; i < 11424; i++)
145 for (i = 0; i < kNDIMX; i++)
147 for (j = 0; j < kNDIMY; j++)
154 for (id = 0; id < ndimXr; id++)
156 for (jd = 0; jd < ndimYr; jd++)
159 i = id+(ndimYr/2-1)-(jd/2);
161 Int_t ij = i + j*kNDIMX;
165 cellenergy[ij] = celladc[jd][id];
167 else if (ismn >= 12 && ismn <= 23)
169 cellenergy[ij] = celladc[id][jd];
174 for (i = 0; i < kNMX; i++)
176 edepcell[i] = cellenergy[i];
180 TMath::Sort((Int_t)kNMX,edepcell,iord1);// order the data
181 cutoff = fCutoff; // cutoff to discard cells
184 for(i = 0;i < kNMX; i++)
190 if(edepcell[i] > cutoff )
196 AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
198 if (nmx1 == 0) nmx1 = 1;
200 AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
203 incr = CrClust(ave, cutoff, nmx1,iord1, edepcell );
204 RefClust(incr,edepcell);
205 Int_t nentries1 = fPMDclucont->GetEntries();
206 AliDebug(1,Form("Detector Plane = %d Serial Module No = %d Number of clusters = %d",idet, ismn, nentries1));
207 AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
209 for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
211 AliPMDcludata *pmdcludata =
212 (AliPMDcludata*)fPMDclucont->UncheckedAt(ient1);
213 Float_t cluXC = pmdcludata->GetClusX();
214 Float_t cluYC = pmdcludata->GetClusY();
215 Float_t cluADC = pmdcludata->GetClusADC();
216 Float_t cluCELLS = pmdcludata->GetClusCells();
217 Float_t cluSIGX = pmdcludata->GetClusSigmaX();
218 Float_t cluSIGY = pmdcludata->GetClusSigmaY();
220 Float_t cluY0 = ktwobysqrt3*cluYC;
221 Float_t cluX0 = cluXC - cluY0/2.;
224 // Cluster X centroid is back transformed
229 clusdata[0] = cluX0 - (24-1) + cluY0/2.;
231 else if ( ismn >= 12 && ismn <= 23)
233 clusdata[0] = cluX0 - (48-1) + cluY0/2.;
237 clusdata[2] = cluADC;
238 clusdata[3] = cluCELLS;
239 clusdata[4] = cluSIGX;
240 clusdata[5] = cluSIGY;
243 // Cells associated with a cluster
246 for (Int_t ihit = 0; ihit < kNmaxCell; ihit++)
248 Int_t cellrow = pmdcludata->GetCellXY(ihit)/10000;
249 Int_t cellcol = pmdcludata->GetCellXY(ihit)%10000;
253 celldataX[ihit] = cellrow - (24-1) + int(cellcol/2.);
255 else if (ismn >= 12 && ismn <= 23)
257 celldataX[ihit] = cellrow - (48-1) + int(cellcol/2.);
260 celldataY[ihit] = cellcol;
262 Int_t irow = celldataX[ihit];
263 Int_t icol = celldataY[ihit];
267 if ((irow >= 0 && irow < 96) && (icol >= 0 && icol < 48))
269 celldataTr[ihit] = celltrack[icol][irow];
270 celldataPid[ihit] = cellpid[icol][irow];
271 celldataAdc[ihit] = (Float_t) celladc[icol][irow];
275 celldataTr[ihit] = -1;
276 celldataPid[ihit] = -1;
277 celldataAdc[ihit] = -1;
280 else if (ismn >= 12 && ismn < 24)
282 if ((irow >= 0 && irow < 48) && (icol >= 0 && icol < 96))
284 celldataTr[ihit] = celltrack[irow][icol];
285 celldataPid[ihit] = cellpid[irow][icol];
286 celldataAdc[ihit] = (Float_t) celladc[irow][icol];
291 celldataTr[ihit] = -1;
292 celldataPid[ihit] = -1;
293 celldataAdc[ihit] = -1;
299 pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY,
300 celldataTr, celldataPid, celldataAdc);
304 fPMDclucont->Delete();
306 // ------------------------------------------------------------------------ //
307 Int_t AliPMDClusteringV1::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
308 Int_t iord1[], Double_t edepcell[])
310 // Does crude clustering
311 // Finds out only the big patch by just searching the
314 const Int_t kndim = 4609;
315 Int_t i,j,k,id1,id2,icl, numcell, clust[2][kndim];
316 Int_t jd1,jd2, icell, cellcount;
317 static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
319 AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
321 for (j = 0; j < kNDIMX; j++)
323 for(k = 0; k < kNDIMY; k++)
325 fInfocl[0][j][k] = 0;
326 fInfocl[1][j][k] = 0;
329 for(i=0; i < kNMX; i++)
337 if(edepcell[j] <= cutoff)
339 fInfocl[0][id1][id2] = -1;
343 // ---------------------------------------------------------------
344 // crude clustering begins. Start with cell having largest adc
345 // count and loop over the cells in descending order of adc count
346 // ---------------------------------------------------------------
351 for(icell = 0; icell <= nmx1; icell++)
357 if(fInfocl[0][id1][id2] == 0 )
362 fInfocl[0][id1][id2] = 1;
363 fInfocl[1][id1][id2] = icl;
364 fInfcl[0][cellcount] = icl;
365 fInfcl[1][cellcount] = id1;
366 fInfcl[2][cellcount] = id2;
368 clust[0][numcell] = id1;
369 clust[1][numcell] = id2;
371 for(i = 1; i < kndim; i++)
375 // ---------------------------------------------------------------
376 // check for adc count in neib. cells. If ne 0 put it in this clust
377 // ---------------------------------------------------------------
378 for(i = 0; i < 6; i++)
380 jd1 = id1 + neibx[i];
381 jd2 = id2 + neiby[i];
382 if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
383 fInfocl[0][jd1][jd2] == 0)
386 fInfocl[0][jd1][jd2] = 2;
387 fInfocl[1][jd1][jd2] = icl;
388 clust[0][numcell] = jd1;
389 clust[1][numcell] = jd2;
391 fInfcl[0][cellcount] = icl;
392 fInfcl[1][cellcount] = jd1;
393 fInfcl[2][cellcount] = jd2;
396 // ---------------------------------------------------------------
397 // check adc count for neighbour's neighbours recursively and
398 // if nonzero, add these to the cluster.
399 // ---------------------------------------------------------------
400 for(i = 1; i < kndim;i++)
406 for(j = 0; j < 6 ; j++)
408 jd1 = id1 + neibx[j];
409 jd2 = id2 + neiby[j];
410 if( (jd1 >= 0 && jd1 < kNDIMX) &&
411 (jd2 >= 0 && jd2 < kNDIMY) &&
412 fInfocl[0][jd1][jd2] == 0 )
414 fInfocl[0][jd1][jd2] = 2;
415 fInfocl[1][jd1][jd2] = icl;
417 clust[0][numcell] = jd1;
418 clust[1][numcell] = jd2;
420 fInfcl[0][cellcount] = icl;
421 fInfcl[1][cellcount] = jd1;
422 fInfcl[2][cellcount] = jd2;
431 // ------------------------------------------------------------------------ //
432 void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
434 // Does the refining of clusters
435 // Takes the big patch and does gaussian fitting and
436 // finds out the more refined clusters
439 AliPMDcludata *pmdcludata = 0;
441 const Int_t kNmaxCell = 19; // # of cells surrounding a cluster center
443 Int_t ndim = incr + 1;
448 Int_t i, j, k, i1, i2, id, icl, itest,ihld, ig, nsupcl,clno, t1, t2;
450 Double_t x1, y1, z1, x2, y2, z2, rr;
452 ncl = new Int_t [ndim];
453 clxy = new Int_t [kNmaxCell];
456 for(i = 0; i<ndim; i++)
459 if (i < 6) clusdata[i] = 0.;
460 if (i < kNmaxCell) clxy[i] = 0;
463 // clno counts the final clusters
464 // nsupcl = # of superclusters; ncl[i]= # of cells in supercluster i
465 // x, y and z store (x,y) coordinates of and energy deposited in a cell
466 // xc, yc store (x,y) coordinates of the cluster center
467 // zc stores the energy deposited in a cluster
468 // rc is cluster radius
473 for(i = 0; i <= incr; i++)
475 if(fInfcl[0][i] != nsupcl)
481 AliWarning("RefClust: Too many superclusters!");
488 AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
493 for(i = 0; i <= nsupcl; i++)
501 AliWarning("RefClust: Too many clusters! more than 4608");
508 i12 = i1 + i2*kNDIMX;
510 clusdata[0] = fCoord[0][i1][i2];
511 clusdata[1] = fCoord[1][i1][i2];
512 clusdata[2] = edepcell[i12];
517 clxy[0] = i1*10000 + i2;
519 for(Int_t icltr = 1; icltr < kNmaxCell; icltr++)
524 pmdcludata = new AliPMDcludata(clusdata,clxy);
525 fPMDclucont->Add(pmdcludata);
533 AliWarning("RefClust: Too many clusters! more than 4608");
539 i12 = i1 + i2*kNDIMX;
541 x1 = fCoord[0][i1][i2];
542 y1 = fCoord[1][i1][i2];
545 clxy[0] = i1*10000 + i2;
551 i22 = i1 + i2*kNDIMX;
552 x2 = fCoord[0][i1][i2];
553 y2 = fCoord[1][i1][i2];
556 clxy[1] = i1*10000 + i2;
559 for(Int_t icltr = 2; icltr < kNmaxCell; icltr++)
564 clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
565 clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
570 pmdcludata = new AliPMDcludata(clusdata,clxy);
571 fPMDclucont->Add(pmdcludata);
576 Int_t *iord, *tc, *t;
577 Double_t *x, *y, *z, *xc, *yc, *zc;
579 iord = new Int_t [ncl[i]+1];
580 tc = new Int_t [ncl[i]+1];
581 t = new Int_t [ncl[i]+1];
583 x = new Double_t [ncl[i]+1];
584 y = new Double_t [ncl[i]+1];
585 z = new Double_t [ncl[i]+1];
586 xc = new Double_t [ncl[i]+1];
587 yc = new Double_t [ncl[i]+1];
588 zc = new Double_t [ncl[i]+1];
590 for( k = 0; k < ncl[i]+1; k++)
603 // super-cluster of more than two cells - broken up into smaller
604 // clusters gaussian centers computed. (peaks separated by > 1 cell)
605 // Begin from cell having largest energy deposited This is first
609 i12 = i1 + i2*kNDIMX;
611 x[0] = fCoord[0][i1][i2];
612 y[0] = fCoord[1][i1][i2];
613 z[0] = edepcell[i12];
614 t[0] = i1*10000 + i2;
618 for(j = 1; j <= ncl[i]; j++)
623 i12 = i1 + i2*kNDIMX;
626 x[j] = fCoord[0][i1][i2];
627 y[j] = fCoord[1][i1][i2];
628 z[j] = edepcell[i12];
629 t[j] = i1*10000 + i2;
633 // arranging cells within supercluster in decreasing order
635 for(j = 1;j <= ncl[i]; j++)
639 for(i1 = 0; i1 < j; i1++)
641 if(itest == 0 && z[iord[i1]] < z[ihld])
644 for(i2 = j-1; i2 >= i1; i2--)
646 iord[i2+1] = iord[i2];
652 /* MODIFICATION PART STARTS (Tapan July 2008)
653 iord[0] is the cell with highest ADC in the crude-cluster
654 ig is the number of local maxima in the crude-cluster
655 For the higest peak we make ig=0 which means first local maximum.
656 Next we go down in terms of the ADC sequence and find out if any
657 more of the cells form local maxima. The definition of local
658 maxima is that all its neighbours are of less ADC compared to it.
665 Int_t ivalid = 0, icount = 0;
667 for(j=1;j<=ncl[i];j++)
673 rr=Distance(x1,y1,xc[ig],yc[ig]);
675 // Check the cells which are outside the neighbours (rr>1.2)
680 for(Int_t j1=1;j1<j;j1++)
683 Float_t rr1=Distance(x1,y1,x[iord[j1]],y[iord[j1]]);
684 if(rr1>1.2) ivalid++;
686 if(ivalid == icount && z1>0.5*zc[ig])
699 // We use simple Gaussian weighting. (Tapan Jan 2005)
700 // compute the number of cells belonging to each cluster.
701 // cell can be shared between several clusters
702 // in the ratio of cluster energy deposition
704 // (1) number of cells belonging to a cluster (ig) and
705 // (2) total ADC of the cluster (ig)
706 // (3) x and y positions of the cluster
713 Double_t *totaladc, *totaladc2, *ncell,*weight;
714 Double_t *xclust, *yclust, *sigxclust, *sigyclust;
715 Double_t *ax, *ay, *ax2, *ay2;
718 status = new Int_t [ncl[i]+1];
719 cellXY = new Int_t *[ncl[i]+1];
721 cellCount = new Int_t [ig+1];
722 totaladc = new Double_t [ig+1];
723 totaladc2 = new Double_t [ig+1];
724 ncell = new Double_t [ig+1];
725 weight = new Double_t [ig+1];
726 xclust = new Double_t [ig+1];
727 yclust = new Double_t [ig+1];
728 sigxclust = new Double_t [ig+1];
729 sigyclust = new Double_t [ig+1];
730 ax = new Double_t [ig+1];
731 ay = new Double_t [ig+1];
732 ax2 = new Double_t [ig+1];
733 ay2 = new Double_t [ig+1];
735 for(j = 0; j < ncl[i]+1; j++)
738 cellXY[j] = new Int_t[ig+1];
741 for(Int_t kcl = 0; kcl < ig+1; kcl++)
756 for(j = 0; j < ncl[i]+1; j++)
761 Double_t sumweight, gweight;
763 for(j = 0;j <= ncl[i]; j++)
770 for(Int_t kcl=0; kcl<=ig; kcl++)
774 rr = Distance(x1,y1,x2,y2);
781 totaladc2[kcl] = z1*z1;
791 for(j = 0; j <= ncl[i]; j++)
804 for(Int_t kcl = 0; kcl <= ig; kcl++)
808 rr = Distance(x1,y1,x2,y2);
809 gweight = exp(-(rr*rr)/(2*(1.2*1.2)));
810 weight[kcl] = zc[kcl] * gweight;
811 sumweight = sumweight + weight[kcl];
813 if(weight[kcl] > max)
820 cellXY[cellCount[maxweight]][maxweight] = iord[j];
822 cellCount[maxweight]++;
826 totaladc[maxweight] += z1;
827 ax[maxweight] += x1*z1;
828 ay[maxweight] += y1*z1;
829 totaladc2[maxweight] += z1*z1;
830 ax2[maxweight] += z1*(x1-x2)*(x1-x2);
831 ay2[maxweight] += z1*(y1-y2)*(y1-y2);
837 for(Int_t kcl = 0; kcl <= ig; kcl++)
840 if(totaladc[kcl] > 0.)
842 xclust[kcl] = (ax[kcl])/ totaladc[kcl];
843 yclust[kcl] = (ay[kcl])/ totaladc[kcl];
846 Float_t sqtotadc = totaladc[kcl]*totaladc[kcl];
847 if(totaladc2[kcl] >= sqtotadc)
849 sigxclust[kcl] = 0.25;
850 sigyclust[kcl] = 0.25;
854 sigxclust[kcl] = (totaladc[kcl]/(sqtotadc-totaladc2[kcl]))*ax2[kcl];
855 sigyclust[kcl] = (totaladc[kcl]/(sqtotadc-totaladc2[kcl]))*ay2[kcl];
859 for(j = 0; j < cellCount[kcl]; j++) clno++;
863 AliWarning("RefClust: Too many clusters! more than 4608");
866 clusdata[0] = xclust[kcl];
867 clusdata[1] = yclust[kcl];
868 clusdata[2] = totaladc[kcl];
869 clusdata[3] = ncell[kcl];
872 if(sigxclust[kcl] > sigyclust[kcl])
874 clusdata[4] = TMath::Sqrt(sigxclust[kcl]);
875 clusdata[5] = TMath::Sqrt(sigyclust[kcl]);
879 clusdata[4] = TMath::Sqrt(sigyclust[kcl]);
880 clusdata[5] = TMath::Sqrt(sigxclust[kcl]);
886 for (Int_t ii = 0; ii < cellCount[kcl]; ii++)
890 clxy[Ncell] = t[cellXY[ii][kcl]];
895 pmdcludata = new AliPMDcludata(clusdata,clxy);
896 fPMDclucont->Add(pmdcludata);
910 for(Int_t jj = 0; jj < ncl[i]+1; jj++) delete [] cellXY[jj];
930 // ------------------------------------------------------------------------ //
931 Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1,
932 Double_t x2, Double_t y2)
934 return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
936 // ------------------------------------------------------------------------ //
937 void AliPMDClusteringV1::SetEdepCut(Float_t decut)
941 // ------------------------------------------------------------------------ //