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 Double_t celladc[48][96], TObjArray *pmdcont)
101 // main function to call other necessary functions to do clustering
104 AliPMDcluster *pmdcl = 0;
106 const float ktwobysqrt3 = 1.1547; // 2./sqrt(3.)
107 const Int_t kNmaxCell = 19; // # of cells surrounding a cluster center
109 Int_t i, j, nmx1, incr, id, jd;
110 Int_t celldataX[kNmaxCell], celldataY[kNmaxCell];
112 Double_t cutoff, ave;
113 Double_t edepcell[kNMX];
116 Double_t *cellenergy = new Double_t [11424];
119 // ndimXr and ndimYr are different because of different module size
129 else if (ismn >= 12 && ismn <= 23)
135 for (i =0; i < 11424; i++)
142 for (i = 0; i < kNDIMX; i++)
144 for (j = 0; j < kNDIMY; j++)
151 for (id = 0; id < ndimXr; id++)
153 for (jd = 0; jd < ndimYr; jd++)
156 i = id+(ndimYr/2-1)-(jd/2);
158 Int_t ij = i + j*kNDIMX;
162 cellenergy[ij] = celladc[jd][id];
164 else if (ismn >= 12 && ismn <= 23)
166 cellenergy[ij] = celladc[id][jd];
171 for (i = 0; i < kNMX; i++)
173 edepcell[i] = cellenergy[i];
176 delete [] cellenergy;
179 TMath::Sort((Int_t)kNMX,edepcell,iord1);// order the data
180 cutoff = fCutoff; // cutoff to discard cells
183 for(i = 0;i < kNMX; i++)
189 if(edepcell[i] > cutoff )
195 AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
197 if (nmx1 == 0) nmx1 = 1;
199 AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
202 incr = CrClust(ave, cutoff, nmx1,iord1, edepcell );
203 RefClust(incr,edepcell);
204 Int_t nentries1 = fPMDclucont->GetEntries();
205 AliDebug(1,Form("Detector Plane = %d Serial Module No = %d Number of clusters = %d",idet, ismn, nentries1));
206 AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
208 for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
210 AliPMDcludata *pmdcludata =
211 (AliPMDcludata*)fPMDclucont->UncheckedAt(ient1);
212 Float_t cluXC = pmdcludata->GetClusX();
213 Float_t cluYC = pmdcludata->GetClusY();
214 Float_t cluADC = pmdcludata->GetClusADC();
215 Float_t cluCELLS = pmdcludata->GetClusCells();
216 Float_t cluSIGX = pmdcludata->GetClusSigmaX();
217 Float_t cluSIGY = pmdcludata->GetClusSigmaY();
219 Float_t cluY0 = ktwobysqrt3*cluYC;
220 Float_t cluX0 = cluXC - cluY0/2.;
223 // Cluster X centroid is back transformed
228 clusdata[0] = cluX0 - (24-1) + cluY0/2.;
230 else if ( ismn >= 12 && ismn <= 23)
232 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;
264 pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY);
268 fPMDclucont->Delete();
270 // ------------------------------------------------------------------------ //
271 Int_t AliPMDClusteringV1::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
272 Int_t iord1[], Double_t edepcell[])
274 // Does crude clustering
275 // Finds out only the big patch by just searching the
278 const Int_t kndim = 4609;
279 Int_t i,j,k,id1,id2,icl, numcell, clust[2][kndim];
280 Int_t jd1,jd2, icell, cellcount;
281 static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
283 AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
285 for (j = 0; j < kNDIMX; j++)
287 for(k = 0; k < kNDIMY; k++)
289 fInfocl[0][j][k] = 0;
290 fInfocl[1][j][k] = 0;
293 for(i=0; i < kNMX; i++)
301 if(edepcell[j] <= cutoff)
303 fInfocl[0][id1][id2] = -1;
307 // ---------------------------------------------------------------
308 // crude clustering begins. Start with cell having largest adc
309 // count and loop over the cells in descending order of adc count
310 // ---------------------------------------------------------------
315 for(icell = 0; icell <= nmx1; icell++)
321 if(fInfocl[0][id1][id2] == 0 )
326 fInfocl[0][id1][id2] = 1;
327 fInfocl[1][id1][id2] = icl;
328 fInfcl[0][cellcount] = icl;
329 fInfcl[1][cellcount] = id1;
330 fInfcl[2][cellcount] = id2;
332 clust[0][numcell] = id1;
333 clust[1][numcell] = id2;
335 for(i = 1; i < kndim; i++)
339 // ---------------------------------------------------------------
340 // check for adc count in neib. cells. If ne 0 put it in this clust
341 // ---------------------------------------------------------------
342 for(i = 0; i < 6; i++)
344 jd1 = id1 + neibx[i];
345 jd2 = id2 + neiby[i];
346 if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
347 fInfocl[0][jd1][jd2] == 0)
350 fInfocl[0][jd1][jd2] = 2;
351 fInfocl[1][jd1][jd2] = icl;
352 clust[0][numcell] = jd1;
353 clust[1][numcell] = jd2;
355 fInfcl[0][cellcount] = icl;
356 fInfcl[1][cellcount] = jd1;
357 fInfcl[2][cellcount] = jd2;
360 // ---------------------------------------------------------------
361 // check adc count for neighbour's neighbours recursively and
362 // if nonzero, add these to the cluster.
363 // ---------------------------------------------------------------
364 for(i = 1; i < kndim;i++)
370 for(j = 0; j < 6 ; j++)
372 jd1 = id1 + neibx[j];
373 jd2 = id2 + neiby[j];
374 if( (jd1 >= 0 && jd1 < kNDIMX) &&
375 (jd2 >= 0 && jd2 < kNDIMY) &&
376 fInfocl[0][jd1][jd2] == 0 )
378 fInfocl[0][jd1][jd2] = 2;
379 fInfocl[1][jd1][jd2] = icl;
381 clust[0][numcell] = jd1;
382 clust[1][numcell] = jd2;
384 fInfcl[0][cellcount] = icl;
385 fInfcl[1][cellcount] = jd1;
386 fInfcl[2][cellcount] = jd2;
395 // ------------------------------------------------------------------------ //
396 void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
398 // Does the refining of clusters
399 // Takes the big patch and does gaussian fitting and
400 // finds out the more refined clusters
403 AliPMDcludata *pmdcludata = 0;
405 const Int_t kNmaxCell = 19; // # of cells surrounding a cluster center
407 Int_t ndim = incr + 1;
412 Int_t i, j, k, i1, i2, id, icl, itest,ihld, ig, nsupcl,clno, t1, t2;
414 Double_t x1, y1, z1, x2, y2, z2, rr;
416 ncl = new Int_t [ndim];
417 clxy = new Int_t [kNmaxCell];
420 for(i = 0; i<ndim; i++)
423 if (i < 6) clusdata[i] = 0.;
424 if (i < kNmaxCell) clxy[i] = 0;
427 // clno counts the final clusters
428 // nsupcl = # of superclusters; ncl[i]= # of cells in supercluster i
429 // x, y and z store (x,y) coordinates of and energy deposited in a cell
430 // xc, yc store (x,y) coordinates of the cluster center
431 // zc stores the energy deposited in a cluster
432 // rc is cluster radius
437 for(i = 0; i <= incr; i++)
439 if(fInfcl[0][i] != nsupcl)
445 AliWarning("RefClust: Too many superclusters!");
452 AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
457 for(i = 0; i <= nsupcl; i++)
465 AliWarning("RefClust: Too many clusters! more than 4608");
472 i12 = i1 + i2*kNDIMX;
474 clusdata[0] = fCoord[0][i1][i2];
475 clusdata[1] = fCoord[1][i1][i2];
476 clusdata[2] = edepcell[i12];
481 clxy[0] = i1*10000 + i2;
483 for(Int_t icltr = 1; icltr < kNmaxCell; icltr++)
487 pmdcludata = new AliPMDcludata(clusdata,clxy);
488 fPMDclucont->Add(pmdcludata);
496 AliWarning("RefClust: Too many clusters! more than 4608");
502 i12 = i1 + i2*kNDIMX;
504 x1 = fCoord[0][i1][i2];
505 y1 = fCoord[1][i1][i2];
508 clxy[0] = i1*10000 + i2;
514 i22 = i1 + i2*kNDIMX;
515 x2 = fCoord[0][i1][i2];
516 y2 = fCoord[1][i1][i2];
519 clxy[1] = i1*10000 + i2;
522 for(Int_t icltr = 2; icltr < kNmaxCell; icltr++)
527 clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
528 clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
533 pmdcludata = new AliPMDcludata(clusdata,clxy);
534 fPMDclucont->Add(pmdcludata);
539 Int_t *iord, *tc, *t;
540 Double_t *x, *y, *z, *xc, *yc, *zc;
542 iord = new Int_t [ncl[i]+1];
543 tc = new Int_t [ncl[i]+1];
544 t = new Int_t [ncl[i]+1];
546 x = new Double_t [ncl[i]+1];
547 y = new Double_t [ncl[i]+1];
548 z = new Double_t [ncl[i]+1];
549 xc = new Double_t [ncl[i]+1];
550 yc = new Double_t [ncl[i]+1];
551 zc = new Double_t [ncl[i]+1];
553 for( k = 0; k < ncl[i]+1; k++)
566 // super-cluster of more than two cells - broken up into smaller
567 // clusters gaussian centers computed. (peaks separated by > 1 cell)
568 // Begin from cell having largest energy deposited This is first
572 i12 = i1 + i2*kNDIMX;
574 x[0] = fCoord[0][i1][i2];
575 y[0] = fCoord[1][i1][i2];
576 z[0] = edepcell[i12];
577 t[0] = i1*10000 + i2;
581 for(j = 1; j <= ncl[i]; j++)
586 i12 = i1 + i2*kNDIMX;
589 x[j] = fCoord[0][i1][i2];
590 y[j] = fCoord[1][i1][i2];
591 z[j] = edepcell[i12];
592 t[j] = i1*10000 + i2;
596 // arranging cells within supercluster in decreasing order
598 for(j = 1;j <= ncl[i]; j++)
602 for(i1 = 0; i1 < j; i1++)
604 if(itest == 0 && z[iord[i1]] < z[ihld])
607 for(i2 = j-1; i2 >= i1; i2--)
609 iord[i2+1] = iord[i2];
615 /* MODIFICATION PART STARTS (Tapan July 2008)
616 iord[0] is the cell with highest ADC in the crude-cluster
617 ig is the number of local maxima in the crude-cluster
618 For the higest peak we make ig=0 which means first local maximum.
619 Next we go down in terms of the ADC sequence and find out if any
620 more of the cells form local maxima. The definition of local
621 maxima is that all its neighbours are of less ADC compared to it.
628 Int_t ivalid = 0, icount = 0;
630 for(j=1;j<=ncl[i];j++)
636 rr=Distance(x1,y1,xc[ig],yc[ig]);
638 // Check the cells which are outside the neighbours (rr>1.2)
643 for(Int_t j1=1;j1<j;j1++)
646 Float_t rr1=Distance(x1,y1,x[iord[j1]],y[iord[j1]]);
647 if(rr1>1.2) ivalid++;
649 if(ivalid == icount && z1>0.5*zc[ig])
662 // We use simple Gaussian weighting. (Tapan Jan 2005)
663 // compute the number of cells belonging to each cluster.
664 // cell can be shared between several clusters
665 // in the ratio of cluster energy deposition
667 // (1) number of cells belonging to a cluster (ig) and
668 // (2) total ADC of the cluster (ig)
669 // (3) x and y positions of the cluster
676 Double_t *totaladc, *totaladc2, *ncell,*weight;
677 Double_t *xclust, *yclust, *sigxclust, *sigyclust;
678 Double_t *ax, *ay, *ax2, *ay2;
681 status = new Int_t [ncl[i]+1];
682 cellXY = new Int_t *[ncl[i]+1];
684 cellCount = new Int_t [ig+1];
685 totaladc = new Double_t [ig+1];
686 totaladc2 = new Double_t [ig+1];
687 ncell = new Double_t [ig+1];
688 weight = new Double_t [ig+1];
689 xclust = new Double_t [ig+1];
690 yclust = new Double_t [ig+1];
691 sigxclust = new Double_t [ig+1];
692 sigyclust = new Double_t [ig+1];
693 ax = new Double_t [ig+1];
694 ay = new Double_t [ig+1];
695 ax2 = new Double_t [ig+1];
696 ay2 = new Double_t [ig+1];
698 for(j = 0; j < ncl[i]+1; j++)
701 cellXY[j] = new Int_t[ig+1];
704 for(Int_t kcl = 0; kcl < ig+1; kcl++)
719 for(j = 0; j < ncl[i]+1; j++)
724 Double_t sumweight, gweight;
726 for(j = 0;j <= ncl[i]; j++)
733 for(Int_t kcl=0; kcl<=ig; kcl++)
737 rr = Distance(x1,y1,x2,y2);
744 totaladc2[kcl] = pow(z1,2);
754 for(j = 0; j <= ncl[i]; j++)
767 for(Int_t kcl = 0; kcl <= ig; kcl++)
771 rr = Distance(x1,y1,x2,y2);
772 gweight = exp(-(rr*rr)/(2*(1.2*1.2)));
773 weight[kcl] = zc[kcl] * gweight;
774 sumweight = sumweight + weight[kcl];
776 if(weight[kcl] > max)
783 cellXY[cellCount[maxweight]][maxweight] = iord[j];
785 cellCount[maxweight]++;
787 for(Int_t kcl = 0; kcl <= ig; kcl++)
791 rr = Distance(x1,y1,x2,y2);
793 // For calculating weighted mean:
794 // Weighted_mean = (\sum w_i x_i) / (\sum w_i)
798 totaladc[kcl] = totaladc[kcl] + z1*(weight[kcl]/sumweight);
799 ncell[kcl] = ncell[kcl] + (weight[kcl]/sumweight);
800 ax[kcl] = ax[kcl] + (x1 * z1 *weight[kcl]/sumweight);
801 ay[kcl] = ay[kcl] + (y1 * z1 *weight[kcl]/sumweight);
803 // For calculating weighted sigma:
804 // Wieghted_sigma= ((\sum w_i)/((\sum w_i)^2 - \sum (w_i^2))) \sum w_i (x_i - \Hat\mu)^2
805 // I assume here x1,y1 are cluster centers, otherwise life becomes too difficult
806 totaladc2[kcl] = totaladc2[kcl] + pow((z1*(weight[kcl]/sumweight)),2);
807 ax2[kcl] = ax2[kcl] + (z1 *weight[kcl]/sumweight) * pow((x1-x2),2);
808 ay2[kcl] = ay2[kcl] + (z1 *weight[kcl]/sumweight) * pow((y1-y2),2);
814 for(Int_t kcl = 0; kcl <= ig; kcl++)
818 if(totaladc[kcl]>0.)xclust[kcl] = (ax[kcl])/ totaladc[kcl];
819 if(totaladc[kcl]>0.)yclust[kcl] = (ay[kcl])/ totaladc[kcl];
821 if(totaladc[kcl]>0.)sigxclust[kcl] = (totaladc[kcl]/(pow(totaladc[kcl],2)-totaladc2[kcl]))*ax2[kcl];
822 if(totaladc[kcl]>0.)sigyclust[kcl] = (totaladc[kcl]/(pow(totaladc[kcl],2)-totaladc2[kcl]))*ay2[kcl];
825 for(j = 0; j < cellCount[kcl]; j++) clno++;
829 AliWarning("RefClust: Too many clusters! more than 4608");
832 clusdata[0] = xclust[kcl];
833 clusdata[1] = yclust[kcl];
834 clusdata[2] = totaladc[kcl];
835 clusdata[3] = ncell[kcl];
836 if(sigxclust[kcl] > sigyclust[kcl])
838 clusdata[4] = pow(sigxclust[kcl],0.5);
839 clusdata[5] = pow(sigyclust[kcl],0.5);
843 clusdata[4] = pow(sigyclust[kcl],0.5);
844 clusdata[5] = pow(sigxclust[kcl],0.5);
850 for (Int_t ii = 0; ii < cellCount[kcl]; ii++)
852 Float_t xx = x[cellXY[ii][kcl]];
853 Float_t yy = y[cellXY[ii][kcl]];
855 rr=Distance(xclust[kcl],yclust[kcl],xx,yy);
857 // We store only the nearest and nest-nearest neighbours
860 clxy[Ncell] = t[cellXY[ii][kcl]];
865 pmdcludata = new AliPMDcludata(clusdata,clxy);
866 fPMDclucont->Add(pmdcludata);
880 for(Int_t jj = 0; jj < ncl[i]+1; jj++) delete [] cellXY[jj];
900 // ------------------------------------------------------------------------ //
901 Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1,
902 Double_t x2, Double_t y2)
904 return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
906 // ------------------------------------------------------------------------ //
907 void AliPMDClusteringV1::SetEdepCut(Float_t decut)
911 // ------------------------------------------------------------------------ //