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()),
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),
84 AliError("Copy constructor not allowed ");
87 // ------------------------------------------------------------------------ //
88 AliPMDClusteringV1 &AliPMDClusteringV1::operator=(const AliPMDClusteringV1& /*pmdclv1*/)
91 AliError("Assignment operator not allowed ");
94 // ------------------------------------------------------------------------ //
95 AliPMDClusteringV1::~AliPMDClusteringV1()
99 // ------------------------------------------------------------------------ //
100 void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn,
101 Int_t celltrack[48][96],
102 Int_t cellpid[48][96],
103 Double_t celladc[48][96],
106 // main function to call other necessary functions to do clustering
109 AliPMDcluster *pmdcl = 0;
111 const float ktwobysqrt3 = 1.1547; // 2./sqrt(3.)
112 const Int_t kNmaxCell = 19; // # of cells surrounding a cluster center
114 Int_t i = 0, j = 0, nmx1 = 0;
115 Int_t incr = 0, id = 0, jd = 0;
116 Int_t celldataX[kNmaxCell], celldataY[kNmaxCell];
117 Int_t celldataTr[kNmaxCell], celldataPid[kNmaxCell];
118 Float_t celldataAdc[kNmaxCell];
119 Float_t clusdata[6] = {0.,0.,0.,0.,0.,0.};
120 Double_t cutoff, ave;
121 Double_t edepcell[kNMX];
122 Double_t cellenergy[kNMX];
124 // ndimXr and ndimYr are different because of different module size
134 else if (ismn >= 12 && ismn <= 23)
140 for (i = 0; i < kNMX; i++)
146 for (i = 0; i < kNDIMX; i++)
148 for (j = 0; j < kNDIMY; j++)
155 for (id = 0; id < ndimXr; id++)
157 for (jd = 0; jd < ndimYr; jd++)
160 i = id+(ndimYr/2-1)-(jd/2);
162 Int_t ij = i + j*kNDIMX;
166 cellenergy[ij] = celladc[jd][id];
168 else if (ismn >= 12 && ismn <= 23)
170 cellenergy[ij] = celladc[id][jd];
175 for (i = 0; i < kNMX; i++)
177 edepcell[i] = cellenergy[i];
181 // the dimension of iord1 is increased twice
183 TMath::Sort((Int_t)kNMX,edepcell,iord1,jsort);// order the data
184 cutoff = fCutoff; // cutoff to discard cells
187 for(i = 0;i < kNMX; i++)
193 if(edepcell[i] > cutoff )
199 AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
201 if (nmx1 == 0) nmx1 = 1;
203 AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
206 incr = CrClust(ave, cutoff, nmx1,iord1, edepcell );
207 RefClust(incr,edepcell);
208 Int_t nentries1 = fPMDclucont->GetEntries();
209 AliDebug(1,Form("Detector Plane = %d Serial Module No = %d Number of clusters = %d",idet, ismn, nentries1));
210 AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
212 for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
214 AliPMDcludata *pmdcludata =
215 (AliPMDcludata*)fPMDclucont->UncheckedAt(ient1);
216 Float_t cluXC = pmdcludata->GetClusX();
217 Float_t cluYC = pmdcludata->GetClusY();
218 Float_t cluADC = pmdcludata->GetClusADC();
219 Float_t cluCELLS = pmdcludata->GetClusCells();
220 Float_t cluSIGX = pmdcludata->GetClusSigmaX();
221 Float_t cluSIGY = pmdcludata->GetClusSigmaY();
223 Float_t cluY0 = ktwobysqrt3*cluYC;
224 Float_t cluX0 = cluXC - cluY0/2.;
227 // Cluster X centroid is back transformed
232 clusdata[0] = cluX0 - (24-1) + cluY0/2.;
234 else if ( ismn >= 12 && ismn <= 23)
236 clusdata[0] = cluX0 - (48-1) + cluY0/2.;
240 clusdata[2] = cluADC;
241 clusdata[3] = cluCELLS;
242 clusdata[4] = cluSIGX;
243 clusdata[5] = cluSIGY;
246 // Cells associated with a cluster
249 for (Int_t ihit = 0; ihit < kNmaxCell; ihit++)
251 Int_t cellrow = pmdcludata->GetCellXY(ihit)/10000;
252 Int_t cellcol = pmdcludata->GetCellXY(ihit)%10000;
256 celldataX[ihit] = cellrow - (24-1) + int(cellcol/2.);
258 else if (ismn >= 12 && ismn <= 23)
260 celldataX[ihit] = cellrow - (48-1) + int(cellcol/2.);
263 celldataY[ihit] = cellcol;
265 Int_t irow = celldataX[ihit];
266 Int_t icol = celldataY[ihit];
270 if ((irow >= 0 && irow < 96) && (icol >= 0 && icol < 48))
272 celldataTr[ihit] = celltrack[icol][irow];
273 celldataPid[ihit] = cellpid[icol][irow];
274 celldataAdc[ihit] = (Float_t) celladc[icol][irow];
278 celldataTr[ihit] = -1;
279 celldataPid[ihit] = -1;
280 celldataAdc[ihit] = -1;
283 else if (ismn >= 12 && ismn < 24)
285 if ((irow >= 0 && irow < 48) && (icol >= 0 && icol < 96))
287 celldataTr[ihit] = celltrack[irow][icol];
288 celldataPid[ihit] = cellpid[irow][icol];
289 celldataAdc[ihit] = (Float_t) celladc[irow][icol];
294 celldataTr[ihit] = -1;
295 celldataPid[ihit] = -1;
296 celldataAdc[ihit] = -1;
302 pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY,
303 celldataTr, celldataPid, celldataAdc);
307 fPMDclucont->Delete();
309 // ------------------------------------------------------------------------ //
310 Int_t AliPMDClusteringV1::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
311 Int_t iord1[], Double_t edepcell[])
313 // Does crude clustering
314 // Finds out only the big patch by just searching the
317 const Int_t kndim = 4609;
318 Int_t i=0,j=0,k=0,id1=0,id2=0,icl=0, numcell=0;
319 Int_t jd1=0,jd2=0, icell=0, cellcount=0;
320 Int_t clust[2][kndim];
321 static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
323 AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
325 for (j = 0; j < kNDIMX; j++)
327 for(k = 0; k < kNDIMY; k++)
329 fInfocl[0][j][k] = 0;
330 fInfocl[1][j][k] = 0;
333 for(i=0; i < kNMX; i++)
341 if(edepcell[j] <= cutoff)
343 fInfocl[0][id1][id2] = -1;
347 // ---------------------------------------------------------------
348 // crude clustering begins. Start with cell having largest adc
349 // count and loop over the cells in descending order of adc count
350 // ---------------------------------------------------------------
355 for(icell = 0; icell <= nmx1; icell++)
361 if(fInfocl[0][id1][id2] == 0 )
366 fInfocl[0][id1][id2] = 1;
367 fInfocl[1][id1][id2] = icl;
368 fInfcl[0][cellcount] = icl;
369 fInfcl[1][cellcount] = id1;
370 fInfcl[2][cellcount] = id2;
372 clust[0][numcell] = id1;
373 clust[1][numcell] = id2;
375 for(i = 1; i < kndim; i++)
379 // ---------------------------------------------------------------
380 // check for adc count in neib. cells. If ne 0 put it in this clust
381 // ---------------------------------------------------------------
382 for(i = 0; i < 6; i++)
384 jd1 = id1 + neibx[i];
385 jd2 = id2 + neiby[i];
386 if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
387 fInfocl[0][jd1][jd2] == 0)
390 fInfocl[0][jd1][jd2] = 2;
391 fInfocl[1][jd1][jd2] = icl;
392 clust[0][numcell] = jd1;
393 clust[1][numcell] = jd2;
395 fInfcl[0][cellcount] = icl;
396 fInfcl[1][cellcount] = jd1;
397 fInfcl[2][cellcount] = jd2;
400 // ---------------------------------------------------------------
401 // check adc count for neighbour's neighbours recursively and
402 // if nonzero, add these to the cluster.
403 // ---------------------------------------------------------------
404 for(i = 1; i < kndim;i++)
410 for(j = 0; j < 6 ; j++)
412 jd1 = id1 + neibx[j];
413 jd2 = id2 + neiby[j];
414 if( (jd1 >= 0 && jd1 < kNDIMX) &&
415 (jd2 >= 0 && jd2 < kNDIMY) &&
416 fInfocl[0][jd1][jd2] == 0 )
418 fInfocl[0][jd1][jd2] = 2;
419 fInfocl[1][jd1][jd2] = icl;
421 clust[0][numcell] = jd1;
422 clust[1][numcell] = jd2;
424 fInfcl[0][cellcount] = icl;
425 fInfcl[1][cellcount] = jd1;
426 fInfcl[2][cellcount] = jd2;
435 // ------------------------------------------------------------------------ //
436 void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
438 // Does the refining of clusters
439 // Takes the big patch and does gaussian fitting and
440 // finds out the more refined clusters
443 AliPMDcludata *pmdcludata = 0;
445 const Int_t kNmaxCell = 19; // # of cells surrounding a cluster center
447 Int_t ndim = incr + 1;
451 Int_t i12 = 0, i22 = 0;
452 Int_t i = 0, j = 0, k = 0;
453 Int_t i1 = 0, i2 = 0, id = 0, icl = 0;
454 Int_t itest = 0, ihld = 0, ig = 0;
455 Int_t nsupcl = 0, clno = 0, t1 = 0, t2 = 0;
457 Double_t x1 = 0, y1 = 0, z1 = 0, x2 = 0, y2 = 0, z2 = 0, rr = 0;
459 ncl = new Int_t [ndim];
460 clxy = new Int_t [kNmaxCell];
463 for(i = 0; i<ndim; i++)
471 for(i = 0; i<19; i++)
476 // clno counts the final clusters
477 // nsupcl = # of superclusters; ncl[i]= # of cells in supercluster i
478 // x, y and z store (x,y) coordinates of and energy deposited in a cell
479 // xc, yc store (x,y) coordinates of the cluster center
480 // zc stores the energy deposited in a cluster
481 // rc is cluster radius
486 for(i = 0; i <= incr; i++)
488 if(fInfcl[0][i] != nsupcl)
494 AliWarning("RefClust: Too many superclusters!");
501 AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
506 for(i = 0; i <= nsupcl; i++)
514 AliWarning("RefClust: Too many clusters! more than 4608");
523 i12 = i1 + i2*kNDIMX;
525 clusdata[0] = fCoord[0][i1][i2];
526 clusdata[1] = fCoord[1][i1][i2];
527 clusdata[2] = edepcell[i12];
532 clxy[0] = i1*10000 + i2;
534 for(Int_t icltr = 1; icltr < kNmaxCell; icltr++)
539 pmdcludata = new AliPMDcludata(clusdata,clxy);
540 fPMDclucont->Add(pmdcludata);
548 AliWarning("RefClust: Too many clusters! more than 4608");
557 i12 = i1 + i2*kNDIMX;
559 x1 = fCoord[0][i1][i2];
560 y1 = fCoord[1][i1][i2];
563 clxy[0] = i1*10000 + i2;
569 i22 = i1 + i2*kNDIMX;
570 x2 = fCoord[0][i1][i2];
571 y2 = fCoord[1][i1][i2];
574 clxy[1] = i1*10000 + i2;
577 for(Int_t icltr = 2; icltr < kNmaxCell; icltr++)
582 clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
583 clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
588 pmdcludata = new AliPMDcludata(clusdata,clxy);
589 fPMDclucont->Add(pmdcludata);
593 Int_t *iord, *tc, *t;
594 Double_t *x, *y, *z, *xc, *yc, *zc;
596 iord = new Int_t [ncl[i]+1];
597 tc = new Int_t [ncl[i]+1];
598 t = new Int_t [ncl[i]+1];
599 x = new Double_t [ncl[i]+1];
600 y = new Double_t [ncl[i]+1];
601 z = new Double_t [ncl[i]+1];
602 xc = new Double_t [ncl[i]+1];
603 yc = new Double_t [ncl[i]+1];
604 zc = new Double_t [ncl[i]+1];
606 for( k = 0; k < ncl[i]+1; k++)
619 // super-cluster of more than two cells - broken up into smaller
620 // clusters gaussian centers computed. (peaks separated by > 1 cell)
621 // Begin from cell having largest energy deposited This is first
625 i12 = i1 + i2*kNDIMX;
627 x[0] = fCoord[0][i1][i2];
628 y[0] = fCoord[1][i1][i2];
629 z[0] = edepcell[i12];
630 t[0] = i1*10000 + i2;
634 for(j = 1; j <= ncl[i]; j++)
639 i12 = i1 + i2*kNDIMX;
641 x[j] = fCoord[0][i1][i2];
642 y[j] = fCoord[1][i1][i2];
643 z[j] = edepcell[i12];
644 t[j] = i1*10000 + i2;
647 // arranging cells within supercluster in decreasing order
649 for(j = 1;j <= ncl[i]; j++)
653 for(i1 = 0; i1 < j; i1++)
655 if(itest == 0 && z[iord[i1]] < z[ihld])
658 for(i2 = j-1; i2 >= i1; i2--)
660 iord[i2+1] = iord[i2];
669 // Clustering algorithm returns from here for PP collisions
670 // for pp, only the output of crude clusterng is taken
671 // sigx and sigy are not calculated at this moment
673 Double_t supx=0., supy=0., supz=0.;
675 for(j = 0;j <= ncl[i]; j++)
677 supx += x[iord[j]]*z[iord[j]];
678 supy += y[iord[j]]*z[iord[j]];
682 clxy[j] = t[iord[j]];
688 for(Int_t ncel = ncl[i] + 1; ncel < kNmaxCell; ncel ++ )
693 clusdata[0] = supx/supz;
694 clusdata[1] = supy/supz;
696 clusdata[3] = ncl[i]+1;
699 pmdcludata = new AliPMDcludata(clusdata,clxy);
700 fPMDclucont->Add(pmdcludata);
703 /* MODIFICATION PART STARTS (Tapan July 2008)
704 iord[0] is the cell with highest ADC in the crude-cluster
705 ig is the number of local maxima in the crude-cluster
706 For the higest peak we make ig=0 which means first local maximum.
707 Next we go down in terms of the ADC sequence and find out if any
708 more of the cells form local maxima. The definition of local
709 maxima is that all its neighbours are of less ADC compared to it.
714 // This part is to split the supercluster
721 Int_t ivalid = 0, icount = 0;
723 for(j=1;j<=ncl[i];j++)
729 rr=Distance(x1,y1,xc[ig],yc[ig]);
731 // Check the cells which are outside the neighbours (rr>1.2)
736 for(Int_t j1=1;j1<j;j1++)
739 Float_t rr1=Distance(x1,y1,x[iord[j1]],y[iord[j1]]);
740 if(rr1>1.2) ivalid++;
742 if(ivalid == icount && z1>0.5*zc[ig])
755 // We use simple Gaussian weighting. (Tapan Jan 2005)
756 // compute the number of cells belonging to each cluster.
757 // cell can be shared between several clusters
758 // in the ratio of cluster energy deposition
760 // (1) number of cells belonging to a cluster (ig) and
761 // (2) total ADC of the cluster (ig)
762 // (3) x and y positions of the cluster
769 Double_t *totaladc, *totaladc2, *ncell,*weight;
770 Double_t *xclust, *yclust, *sigxclust, *sigyclust;
771 Double_t *ax, *ay, *ax2, *ay2;
774 status = new Int_t [ncl[i]+1];
775 cellXY = new Int_t *[ncl[i]+1];
777 cellCount = new Int_t [ig+1];
778 totaladc = new Double_t [ig+1];
779 totaladc2 = new Double_t [ig+1];
780 ncell = new Double_t [ig+1];
781 weight = new Double_t [ig+1];
782 xclust = new Double_t [ig+1];
783 yclust = new Double_t [ig+1];
784 sigxclust = new Double_t [ig+1];
785 sigyclust = new Double_t [ig+1];
786 ax = new Double_t [ig+1];
787 ay = new Double_t [ig+1];
788 ax2 = new Double_t [ig+1];
789 ay2 = new Double_t [ig+1];
791 for(j = 0; j < ncl[i]+1; j++)
794 cellXY[j] = new Int_t[ig+1];
797 for(Int_t kcl = 0; kcl < ig+1; kcl++)
812 for(j = 0; j < ncl[i]+1; j++)
817 Double_t sumweight, gweight;
819 for(j = 0;j <= ncl[i]; j++)
826 for(Int_t kcl=0; kcl<=ig; kcl++)
830 rr = Distance(x1,y1,x2,y2);
837 totaladc2[kcl] = z1*z1;
847 for(j = 0; j <= ncl[i]; j++)
860 for(Int_t kcl = 0; kcl <= ig; kcl++)
864 rr = Distance(x1,y1,x2,y2);
865 gweight = exp(-(rr*rr)/(2*(1.2*1.2)));
866 weight[kcl] = zc[kcl] * gweight;
867 sumweight = sumweight + weight[kcl];
869 if(weight[kcl] > max)
876 cellXY[cellCount[maxweight]][maxweight] = iord[j];
878 cellCount[maxweight]++;
882 totaladc[maxweight] += z1;
883 ax[maxweight] += x1*z1;
884 ay[maxweight] += y1*z1;
885 totaladc2[maxweight] += z1*z1;
886 ax2[maxweight] += z1*(x1-x2)*(x1-x2);
887 ay2[maxweight] += z1*(y1-y2)*(y1-y2);
893 for(Int_t kcl = 0; kcl <= ig; kcl++)
895 if(totaladc[kcl] > 0.)
897 xclust[kcl] = (ax[kcl])/ totaladc[kcl];
898 yclust[kcl] = (ay[kcl])/ totaladc[kcl];
901 Float_t sqtotadc = totaladc[kcl]*totaladc[kcl];
902 if(totaladc2[kcl] >= sqtotadc)
904 sigxclust[kcl] = 0.25;
905 sigyclust[kcl] = 0.25;
909 sigxclust[kcl] = (totaladc[kcl]/(sqtotadc-totaladc2[kcl]))*ax2[kcl];
910 sigyclust[kcl] = (totaladc[kcl]/(sqtotadc-totaladc2[kcl]))*ay2[kcl];
914 for(j = 0; j < cellCount[kcl]; j++) clno++;
918 AliWarning("RefClust: Too many clusters! more than 4608");
921 for(Int_t jj = 0; jj < ncl[i]+1; jj++) delete [] cellXY[jj];
954 clusdata[0] = xclust[kcl];
955 clusdata[1] = yclust[kcl];
956 clusdata[2] = totaladc[kcl];
957 clusdata[3] = ncell[kcl];
959 if(sigxclust[kcl] > sigyclust[kcl])
961 clusdata[4] = TMath::Sqrt(sigxclust[kcl]);
962 clusdata[5] = TMath::Sqrt(sigyclust[kcl]);
966 clusdata[4] = TMath::Sqrt(sigyclust[kcl]);
967 clusdata[5] = TMath::Sqrt(sigxclust[kcl]);
973 for (Int_t ii = 0; ii < cellCount[kcl]; ii++)
977 clxy[Ncell] = t[cellXY[ii][kcl]];
982 pmdcludata = new AliPMDcludata(clusdata,clxy);
983 fPMDclucont->Add(pmdcludata);
986 for(Int_t jj = 0; jj < ncl[i]+1; jj++) delete [] cellXY[jj];
1021 // ------------------------------------------------------------------------ //
1022 Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1,
1023 Double_t x2, Double_t y2)
1025 return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
1027 // ------------------------------------------------------------------------ //
1028 void AliPMDClusteringV1::SetEdepCut(Float_t decut)
1032 // ------------------------------------------------------------------------ //
1033 void AliPMDClusteringV1::SetClusteringParam(Int_t cluspar)
1035 fClusParam = cluspar;
1037 // ------------------------------------------------------------------------ //