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 "AliPMDisocell.h"
54 #include "AliPMDClustering.h"
55 #include "AliPMDClusteringV1.h"
59 ClassImp(AliPMDClusteringV1)
61 const Double_t AliPMDClusteringV1::fgkSqroot3by2=0.8660254; // sqrt(3.)/2.
63 AliPMDClusteringV1::AliPMDClusteringV1():
64 fPMDclucont(new TObjArray()),
67 for(Int_t i = 0; i < kNDIMX; i++)
69 for(Int_t j = 0; j < kNDIMY; j++)
71 fCoord[0][i][j] = i+j/2.;
72 fCoord[1][i][j] = fgkSqroot3by2*j;
76 // ------------------------------------------------------------------------ //
77 AliPMDClusteringV1::AliPMDClusteringV1(const AliPMDClusteringV1& pmdclv1):
78 AliPMDClustering(pmdclv1),
83 AliError("Copy constructor not allowed ");
86 // ------------------------------------------------------------------------ //
87 AliPMDClusteringV1 &AliPMDClusteringV1::operator=(const AliPMDClusteringV1& /*pmdclv1*/)
90 AliError("Assignment operator not allowed ");
93 // ------------------------------------------------------------------------ //
94 AliPMDClusteringV1::~AliPMDClusteringV1()
98 // ------------------------------------------------------------------------ //
99 void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn,
100 Int_t celltrack[48][96],
101 Int_t cellpid[48][96],
102 Double_t celladc[48][96],
103 TObjArray *pmdisocell, TObjArray *pmdcont)
105 // main function to call other necessary functions to do clustering
108 AliPMDcluster *pmdcl = 0;
110 const float ktwobysqrt3 = 1.1547; // 2./sqrt(3.)
111 const Int_t kNmaxCell = 19; // # of cells surrounding a cluster center
113 Int_t i, j, nmx1, incr, id, jd;
114 Int_t celldataX[kNmaxCell], celldataY[kNmaxCell];
115 Int_t celldataTr[kNmaxCell], celldataPid[kNmaxCell];
116 Float_t celldataAdc[kNmaxCell];
118 Double_t cutoff, ave;
119 Double_t edepcell[kNMX];
122 Double_t *cellenergy = new Double_t [11424];
125 // call the isolated cell search method
127 CalculateIsoCell(idet, ismn, celladc, pmdisocell);
129 // ndimXr and ndimYr are different because of different module size
139 else if (ismn >= 12 && ismn <= 23)
145 for (i =0; i < 11424; i++)
151 for (i = 0; i < kNDIMX; i++)
153 for (j = 0; j < kNDIMY; j++)
160 for (id = 0; id < ndimXr; id++)
162 for (jd = 0; jd < ndimYr; jd++)
165 i = id+(ndimYr/2-1)-(jd/2);
167 Int_t ij = i + j*kNDIMX;
171 cellenergy[ij] = celladc[jd][id];
173 else if (ismn >= 12 && ismn <= 23)
175 cellenergy[ij] = celladc[id][jd];
180 for (i = 0; i < kNMX; i++)
182 edepcell[i] = cellenergy[i];
185 delete [] cellenergy;
188 TMath::Sort((Int_t)kNMX,edepcell,iord1);// order the data
189 cutoff = fCutoff; // cutoff to discard cells
192 for(i = 0;i < kNMX; i++)
198 if(edepcell[i] > cutoff )
204 AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
206 if (nmx1 == 0) nmx1 = 1;
208 AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
211 incr = CrClust(ave, cutoff, nmx1,iord1, edepcell );
212 RefClust(incr,edepcell);
213 Int_t nentries1 = fPMDclucont->GetEntries();
214 AliDebug(1,Form("Detector Plane = %d Serial Module No = %d Number of clusters = %d",idet, ismn, nentries1));
215 AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
217 for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
219 AliPMDcludata *pmdcludata =
220 (AliPMDcludata*)fPMDclucont->UncheckedAt(ient1);
221 Float_t cluXC = pmdcludata->GetClusX();
222 Float_t cluYC = pmdcludata->GetClusY();
223 Float_t cluADC = pmdcludata->GetClusADC();
224 Float_t cluCELLS = pmdcludata->GetClusCells();
225 Float_t cluSIGX = pmdcludata->GetClusSigmaX();
226 Float_t cluSIGY = pmdcludata->GetClusSigmaY();
228 Float_t cluY0 = ktwobysqrt3*cluYC;
229 Float_t cluX0 = cluXC - cluY0/2.;
232 // Cluster X centroid is back transformed
237 clusdata[0] = cluX0 - (24-1) + cluY0/2.;
239 else if ( ismn >= 12 && ismn <= 23)
241 clusdata[0] = cluX0 - (48-1) + cluY0/2.;
245 clusdata[2] = cluADC;
246 clusdata[3] = cluCELLS;
247 clusdata[4] = cluSIGX;
248 clusdata[5] = cluSIGY;
251 // Cells associated with a cluster
254 for (Int_t ihit = 0; ihit < kNmaxCell; ihit++)
256 Int_t cellrow = pmdcludata->GetCellXY(ihit)/10000;
257 Int_t cellcol = pmdcludata->GetCellXY(ihit)%10000;
261 celldataX[ihit] = cellrow - (24-1) + int(cellcol/2.);
263 else if (ismn >= 12 && ismn <= 23)
265 celldataX[ihit] = cellrow - (48-1) + int(cellcol/2.);
268 celldataY[ihit] = cellcol;
270 Int_t irow = celldataX[ihit];
271 Int_t icol = celldataY[ihit];
273 if ((irow >= 0 && irow < 48) && (icol >= 0 && icol < 96))
275 celldataTr[ihit] = celltrack[irow][icol];
276 celldataPid[ihit] = cellpid[irow][icol];
277 celldataAdc[ihit] = (Float_t) celladc[irow][icol];
281 celldataTr[ihit] = -1;
282 celldataPid[ihit] = -1;
283 celldataAdc[ihit] = -1;
287 pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY,
288 celldataTr, celldataPid, celldataAdc);
292 fPMDclucont->Delete();
294 // ------------------------------------------------------------------------ //
295 Int_t AliPMDClusteringV1::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
296 Int_t iord1[], Double_t edepcell[])
298 // Does crude clustering
299 // Finds out only the big patch by just searching the
302 const Int_t kndim = 4609;
303 Int_t i,j,k,id1,id2,icl, numcell, clust[2][kndim];
304 Int_t jd1,jd2, icell, cellcount;
305 static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
307 AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
309 for (j = 0; j < kNDIMX; j++)
311 for(k = 0; k < kNDIMY; k++)
313 fInfocl[0][j][k] = 0;
314 fInfocl[1][j][k] = 0;
317 for(i=0; i < kNMX; i++)
325 if(edepcell[j] <= cutoff)
327 fInfocl[0][id1][id2] = -1;
331 // ---------------------------------------------------------------
332 // crude clustering begins. Start with cell having largest adc
333 // count and loop over the cells in descending order of adc count
334 // ---------------------------------------------------------------
339 for(icell = 0; icell <= nmx1; icell++)
345 if(fInfocl[0][id1][id2] == 0 )
350 fInfocl[0][id1][id2] = 1;
351 fInfocl[1][id1][id2] = icl;
352 fInfcl[0][cellcount] = icl;
353 fInfcl[1][cellcount] = id1;
354 fInfcl[2][cellcount] = id2;
356 clust[0][numcell] = id1;
357 clust[1][numcell] = id2;
359 for(i = 1; i < kndim; i++)
363 // ---------------------------------------------------------------
364 // check for adc count in neib. cells. If ne 0 put it in this clust
365 // ---------------------------------------------------------------
366 for(i = 0; i < 6; i++)
368 jd1 = id1 + neibx[i];
369 jd2 = id2 + neiby[i];
370 if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
371 fInfocl[0][jd1][jd2] == 0)
374 fInfocl[0][jd1][jd2] = 2;
375 fInfocl[1][jd1][jd2] = icl;
376 clust[0][numcell] = jd1;
377 clust[1][numcell] = jd2;
379 fInfcl[0][cellcount] = icl;
380 fInfcl[1][cellcount] = jd1;
381 fInfcl[2][cellcount] = jd2;
384 // ---------------------------------------------------------------
385 // check adc count for neighbour's neighbours recursively and
386 // if nonzero, add these to the cluster.
387 // ---------------------------------------------------------------
388 for(i = 1; i < kndim;i++)
394 for(j = 0; j < 6 ; j++)
396 jd1 = id1 + neibx[j];
397 jd2 = id2 + neiby[j];
398 if( (jd1 >= 0 && jd1 < kNDIMX) &&
399 (jd2 >= 0 && jd2 < kNDIMY) &&
400 fInfocl[0][jd1][jd2] == 0 )
402 fInfocl[0][jd1][jd2] = 2;
403 fInfocl[1][jd1][jd2] = icl;
405 clust[0][numcell] = jd1;
406 clust[1][numcell] = jd2;
408 fInfcl[0][cellcount] = icl;
409 fInfcl[1][cellcount] = jd1;
410 fInfcl[2][cellcount] = jd2;
419 // ------------------------------------------------------------------------ //
420 void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
422 // Does the refining of clusters
423 // Takes the big patch and does gaussian fitting and
424 // finds out the more refined clusters
427 AliPMDcludata *pmdcludata = 0;
429 const Int_t kNmaxCell = 19; // # of cells surrounding a cluster center
431 Int_t ndim = incr + 1;
436 Int_t i, j, k, i1, i2, id, icl, itest,ihld, ig, nsupcl,clno, t1, t2;
438 Double_t x1, y1, z1, x2, y2, z2, rr;
440 ncl = new Int_t [ndim];
441 clxy = new Int_t [kNmaxCell];
444 for(i = 0; i<ndim; i++)
447 if (i < 6) clusdata[i] = 0.;
448 if (i < kNmaxCell) clxy[i] = 0;
451 // clno counts the final clusters
452 // nsupcl = # of superclusters; ncl[i]= # of cells in supercluster i
453 // x, y and z store (x,y) coordinates of and energy deposited in a cell
454 // xc, yc store (x,y) coordinates of the cluster center
455 // zc stores the energy deposited in a cluster
456 // rc is cluster radius
461 for(i = 0; i <= incr; i++)
463 if(fInfcl[0][i] != nsupcl)
469 AliWarning("RefClust: Too many superclusters!");
476 AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
481 for(i = 0; i <= nsupcl; i++)
489 AliWarning("RefClust: Too many clusters! more than 4608");
496 i12 = i1 + i2*kNDIMX;
498 clusdata[0] = fCoord[0][i1][i2];
499 clusdata[1] = fCoord[1][i1][i2];
500 clusdata[2] = edepcell[i12];
505 clxy[0] = i1*10000 + i2;
507 for(Int_t icltr = 1; icltr < kNmaxCell; icltr++)
511 pmdcludata = new AliPMDcludata(clusdata,clxy);
512 fPMDclucont->Add(pmdcludata);
520 AliWarning("RefClust: Too many clusters! more than 4608");
526 i12 = i1 + i2*kNDIMX;
528 x1 = fCoord[0][i1][i2];
529 y1 = fCoord[1][i1][i2];
532 clxy[0] = i1*10000 + i2;
538 i22 = i1 + i2*kNDIMX;
539 x2 = fCoord[0][i1][i2];
540 y2 = fCoord[1][i1][i2];
543 clxy[1] = i1*10000 + i2;
546 for(Int_t icltr = 2; icltr < kNmaxCell; icltr++)
551 clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
552 clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
557 pmdcludata = new AliPMDcludata(clusdata,clxy);
558 fPMDclucont->Add(pmdcludata);
563 Int_t *iord, *tc, *t;
564 Double_t *x, *y, *z, *xc, *yc, *zc;
566 iord = new Int_t [ncl[i]+1];
567 tc = new Int_t [ncl[i]+1];
568 t = new Int_t [ncl[i]+1];
570 x = new Double_t [ncl[i]+1];
571 y = new Double_t [ncl[i]+1];
572 z = new Double_t [ncl[i]+1];
573 xc = new Double_t [ncl[i]+1];
574 yc = new Double_t [ncl[i]+1];
575 zc = new Double_t [ncl[i]+1];
577 for( k = 0; k < ncl[i]+1; k++)
590 // super-cluster of more than two cells - broken up into smaller
591 // clusters gaussian centers computed. (peaks separated by > 1 cell)
592 // Begin from cell having largest energy deposited This is first
596 i12 = i1 + i2*kNDIMX;
598 x[0] = fCoord[0][i1][i2];
599 y[0] = fCoord[1][i1][i2];
600 z[0] = edepcell[i12];
601 t[0] = i1*10000 + i2;
605 for(j = 1; j <= ncl[i]; j++)
610 i12 = i1 + i2*kNDIMX;
613 x[j] = fCoord[0][i1][i2];
614 y[j] = fCoord[1][i1][i2];
615 z[j] = edepcell[i12];
616 t[j] = i1*10000 + i2;
620 // arranging cells within supercluster in decreasing order
622 for(j = 1;j <= ncl[i]; j++)
626 for(i1 = 0; i1 < j; i1++)
628 if(itest == 0 && z[iord[i1]] < z[ihld])
631 for(i2 = j-1; i2 >= i1; i2--)
633 iord[i2+1] = iord[i2];
639 /* MODIFICATION PART STARTS (Tapan July 2008)
640 iord[0] is the cell with highest ADC in the crude-cluster
641 ig is the number of local maxima in the crude-cluster
642 For the higest peak we make ig=0 which means first local maximum.
643 Next we go down in terms of the ADC sequence and find out if any
644 more of the cells form local maxima. The definition of local
645 maxima is that all its neighbours are of less ADC compared to it.
652 Int_t ivalid = 0, icount = 0;
654 for(j=1;j<=ncl[i];j++)
660 rr=Distance(x1,y1,xc[ig],yc[ig]);
662 // Check the cells which are outside the neighbours (rr>1.2)
667 for(Int_t j1=1;j1<j;j1++)
670 Float_t rr1=Distance(x1,y1,x[iord[j1]],y[iord[j1]]);
671 if(rr1>1.2) ivalid++;
673 if(ivalid == icount && z1>0.5*zc[ig])
686 // We use simple Gaussian weighting. (Tapan Jan 2005)
687 // compute the number of cells belonging to each cluster.
688 // cell can be shared between several clusters
689 // in the ratio of cluster energy deposition
691 // (1) number of cells belonging to a cluster (ig) and
692 // (2) total ADC of the cluster (ig)
693 // (3) x and y positions of the cluster
700 Double_t *totaladc, *totaladc2, *ncell,*weight;
701 Double_t *xclust, *yclust, *sigxclust, *sigyclust;
702 Double_t *ax, *ay, *ax2, *ay2;
705 status = new Int_t [ncl[i]+1];
706 cellXY = new Int_t *[ncl[i]+1];
708 cellCount = new Int_t [ig+1];
709 totaladc = new Double_t [ig+1];
710 totaladc2 = new Double_t [ig+1];
711 ncell = new Double_t [ig+1];
712 weight = new Double_t [ig+1];
713 xclust = new Double_t [ig+1];
714 yclust = new Double_t [ig+1];
715 sigxclust = new Double_t [ig+1];
716 sigyclust = new Double_t [ig+1];
717 ax = new Double_t [ig+1];
718 ay = new Double_t [ig+1];
719 ax2 = new Double_t [ig+1];
720 ay2 = new Double_t [ig+1];
722 for(j = 0; j < ncl[i]+1; j++)
725 cellXY[j] = new Int_t[ig+1];
728 for(Int_t kcl = 0; kcl < ig+1; kcl++)
743 for(j = 0; j < ncl[i]+1; j++)
748 Double_t sumweight, gweight;
750 for(j = 0;j <= ncl[i]; j++)
757 for(Int_t kcl=0; kcl<=ig; kcl++)
761 rr = Distance(x1,y1,x2,y2);
768 totaladc2[kcl] = pow(z1,2);
778 for(j = 0; j <= ncl[i]; j++)
791 for(Int_t kcl = 0; kcl <= ig; kcl++)
795 rr = Distance(x1,y1,x2,y2);
796 gweight = exp(-(rr*rr)/(2*(1.2*1.2)));
797 weight[kcl] = zc[kcl] * gweight;
798 sumweight = sumweight + weight[kcl];
800 if(weight[kcl] > max)
807 cellXY[cellCount[maxweight]][maxweight] = iord[j];
809 cellCount[maxweight]++;
811 for(Int_t kcl = 0; kcl <= ig; kcl++)
815 rr = Distance(x1,y1,x2,y2);
817 // For calculating weighted mean:
818 // Weighted_mean = (\sum w_i x_i) / (\sum w_i)
822 totaladc[kcl] = totaladc[kcl] + z1*(weight[kcl]/sumweight);
823 ncell[kcl] = ncell[kcl] + (weight[kcl]/sumweight);
824 ax[kcl] = ax[kcl] + (x1 * z1 *weight[kcl]/sumweight);
825 ay[kcl] = ay[kcl] + (y1 * z1 *weight[kcl]/sumweight);
827 // For calculating weighted sigma:
828 // Wieghted_sigma= ((\sum w_i)/((\sum w_i)^2 - \sum (w_i^2))) \sum w_i (x_i - \Hat\mu)^2
829 // I assume here x1,y1 are cluster centers, otherwise life becomes too difficult
830 totaladc2[kcl] = totaladc2[kcl] + pow((z1*(weight[kcl]/sumweight)),2);
831 ax2[kcl] = ax2[kcl] + (z1 *weight[kcl]/sumweight) * pow((x1-x2),2);
832 ay2[kcl] = ay2[kcl] + (z1 *weight[kcl]/sumweight) * pow((y1-y2),2);
838 for(Int_t kcl = 0; kcl <= ig; kcl++)
842 if(totaladc[kcl]>0.)xclust[kcl] = (ax[kcl])/ totaladc[kcl];
843 if(totaladc[kcl]>0.)yclust[kcl] = (ay[kcl])/ totaladc[kcl];
845 if(totaladc[kcl]>0.)sigxclust[kcl] = (totaladc[kcl]/(pow(totaladc[kcl],2)-totaladc2[kcl]))*ax2[kcl];
846 if(totaladc[kcl]>0.)sigyclust[kcl] = (totaladc[kcl]/(pow(totaladc[kcl],2)-totaladc2[kcl]))*ay2[kcl];
849 for(j = 0; j < cellCount[kcl]; j++) clno++;
853 AliWarning("RefClust: Too many clusters! more than 4608");
856 clusdata[0] = xclust[kcl];
857 clusdata[1] = yclust[kcl];
858 clusdata[2] = totaladc[kcl];
859 clusdata[3] = ncell[kcl];
860 if(sigxclust[kcl] > sigyclust[kcl])
862 clusdata[4] = pow(sigxclust[kcl],0.5);
863 clusdata[5] = pow(sigyclust[kcl],0.5);
867 clusdata[4] = pow(sigyclust[kcl],0.5);
868 clusdata[5] = pow(sigxclust[kcl],0.5);
874 for (Int_t ii = 0; ii < cellCount[kcl]; ii++)
878 clxy[Ncell] = t[cellXY[ii][kcl]];
883 pmdcludata = new AliPMDcludata(clusdata,clxy);
884 fPMDclucont->Add(pmdcludata);
898 for(Int_t jj = 0; jj < ncl[i]+1; jj++) delete [] cellXY[jj];
918 // ------------------------------------------------------------------------ //
919 Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1,
920 Double_t x2, Double_t y2)
922 return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
924 // ------------------------------------------------------------------------ //
925 void AliPMDClusteringV1::CalculateIsoCell(Int_t idet, Int_t ismn, Double_t celladc[][96], TObjArray *pmdisocell)
927 // Does isolated cell search for offline calibration
929 AliPMDisocell *isocell = 0;
931 const Int_t kMaxRow = 48;
932 const Int_t kMaxCol = 96;
933 const Int_t kCellNeighbour = 6;
937 Int_t neibx[6] = {1,0,-1,-1,0,1};
938 Int_t neiby[6] = {0,1,1,0,-1,-1};
941 for(Int_t irow = 0; irow < kMaxRow; irow++)
943 for(Int_t icol = 0; icol < kMaxCol; icol++)
945 if(celladc[irow][icol] > 0)
948 for(Int_t ii = 0; ii < kCellNeighbour; ii++)
950 id1 = irow + neibx[ii];
951 jd1 = icol + neiby[ii];
952 Float_t adc = (Float_t) celladc[id1][jd1];
956 if(isocount == kCellNeighbour)
958 Float_t cadc = (Float_t) celladc[irow][icol];
960 isocell = new AliPMDisocell(idet,ismn,irow,icol,cadc);
961 pmdisocell->Add(isocell);
965 } // neigh cell cond.
972 // ------------------------------------------------------------------------ //
973 void AliPMDClusteringV1::SetEdepCut(Float_t decut)
977 // ------------------------------------------------------------------------ //