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, j, nmx1, incr, id, jd;
115 Int_t celldataX[kNmaxCell], celldataY[kNmaxCell];
116 Int_t celldataTr[kNmaxCell], celldataPid[kNmaxCell];
117 Float_t celldataAdc[kNmaxCell];
119 Double_t cutoff, ave;
120 Double_t edepcell[kNMX];
123 Double_t cellenergy[11424];
125 // ndimXr and ndimYr are different because of different module size
135 else if (ismn >= 12 && ismn <= 23)
141 for (i =0; i < 11424; i++)
147 for (i = 0; i < kNDIMX; i++)
149 for (j = 0; j < kNDIMY; j++)
156 for (id = 0; id < ndimXr; id++)
158 for (jd = 0; jd < ndimYr; jd++)
161 i = id+(ndimYr/2-1)-(jd/2);
163 Int_t ij = i + j*kNDIMX;
167 cellenergy[ij] = celladc[jd][id];
169 else if (ismn >= 12 && ismn <= 23)
171 cellenergy[ij] = celladc[id][jd];
176 for (i = 0; i < kNMX; i++)
178 edepcell[i] = cellenergy[i];
182 TMath::Sort((Int_t)kNMX,edepcell,iord1);// order the data
183 cutoff = fCutoff; // cutoff to discard cells
186 for(i = 0;i < kNMX; i++)
192 if(edepcell[i] > cutoff )
198 AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
200 if (nmx1 == 0) nmx1 = 1;
202 AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
205 incr = CrClust(ave, cutoff, nmx1,iord1, edepcell );
206 RefClust(incr,edepcell);
207 Int_t nentries1 = fPMDclucont->GetEntries();
208 AliDebug(1,Form("Detector Plane = %d Serial Module No = %d Number of clusters = %d",idet, ismn, nentries1));
209 AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
211 for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
213 AliPMDcludata *pmdcludata =
214 (AliPMDcludata*)fPMDclucont->UncheckedAt(ient1);
215 Float_t cluXC = pmdcludata->GetClusX();
216 Float_t cluYC = pmdcludata->GetClusY();
217 Float_t cluADC = pmdcludata->GetClusADC();
218 Float_t cluCELLS = pmdcludata->GetClusCells();
219 Float_t cluSIGX = pmdcludata->GetClusSigmaX();
220 Float_t cluSIGY = pmdcludata->GetClusSigmaY();
222 Float_t cluY0 = ktwobysqrt3*cluYC;
223 Float_t cluX0 = cluXC - cluY0/2.;
226 // Cluster X centroid is back transformed
231 clusdata[0] = cluX0 - (24-1) + cluY0/2.;
233 else if ( ismn >= 12 && ismn <= 23)
235 clusdata[0] = cluX0 - (48-1) + cluY0/2.;
239 clusdata[2] = cluADC;
240 clusdata[3] = cluCELLS;
241 clusdata[4] = cluSIGX;
242 clusdata[5] = cluSIGY;
245 // Cells associated with a cluster
248 for (Int_t ihit = 0; ihit < kNmaxCell; ihit++)
250 Int_t cellrow = pmdcludata->GetCellXY(ihit)/10000;
251 Int_t cellcol = pmdcludata->GetCellXY(ihit)%10000;
255 celldataX[ihit] = cellrow - (24-1) + int(cellcol/2.);
257 else if (ismn >= 12 && ismn <= 23)
259 celldataX[ihit] = cellrow - (48-1) + int(cellcol/2.);
262 celldataY[ihit] = cellcol;
264 Int_t irow = celldataX[ihit];
265 Int_t icol = celldataY[ihit];
269 if ((irow >= 0 && irow < 96) && (icol >= 0 && icol < 48))
271 celldataTr[ihit] = celltrack[icol][irow];
272 celldataPid[ihit] = cellpid[icol][irow];
273 celldataAdc[ihit] = (Float_t) celladc[icol][irow];
277 celldataTr[ihit] = -1;
278 celldataPid[ihit] = -1;
279 celldataAdc[ihit] = -1;
282 else if (ismn >= 12 && ismn < 24)
284 if ((irow >= 0 && irow < 48) && (icol >= 0 && icol < 96))
286 celldataTr[ihit] = celltrack[irow][icol];
287 celldataPid[ihit] = cellpid[irow][icol];
288 celldataAdc[ihit] = (Float_t) celladc[irow][icol];
293 celldataTr[ihit] = -1;
294 celldataPid[ihit] = -1;
295 celldataAdc[ihit] = -1;
301 pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY,
302 celldataTr, celldataPid, celldataAdc);
306 fPMDclucont->Delete();
308 // ------------------------------------------------------------------------ //
309 Int_t AliPMDClusteringV1::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
310 Int_t iord1[], Double_t edepcell[])
312 // Does crude clustering
313 // Finds out only the big patch by just searching the
316 const Int_t kndim = 4609;
317 Int_t i,j,k,id1,id2,icl, numcell, clust[2][kndim];
318 Int_t jd1,jd2, icell, cellcount;
319 static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
321 AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
323 for (j = 0; j < kNDIMX; j++)
325 for(k = 0; k < kNDIMY; k++)
327 fInfocl[0][j][k] = 0;
328 fInfocl[1][j][k] = 0;
331 for(i=0; i < kNMX; i++)
339 if(edepcell[j] <= cutoff)
341 fInfocl[0][id1][id2] = -1;
345 // ---------------------------------------------------------------
346 // crude clustering begins. Start with cell having largest adc
347 // count and loop over the cells in descending order of adc count
348 // ---------------------------------------------------------------
353 for(icell = 0; icell <= nmx1; icell++)
359 if(fInfocl[0][id1][id2] == 0 )
364 fInfocl[0][id1][id2] = 1;
365 fInfocl[1][id1][id2] = icl;
366 fInfcl[0][cellcount] = icl;
367 fInfcl[1][cellcount] = id1;
368 fInfcl[2][cellcount] = id2;
370 clust[0][numcell] = id1;
371 clust[1][numcell] = id2;
373 for(i = 1; i < kndim; i++)
377 // ---------------------------------------------------------------
378 // check for adc count in neib. cells. If ne 0 put it in this clust
379 // ---------------------------------------------------------------
380 for(i = 0; i < 6; i++)
382 jd1 = id1 + neibx[i];
383 jd2 = id2 + neiby[i];
384 if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
385 fInfocl[0][jd1][jd2] == 0)
388 fInfocl[0][jd1][jd2] = 2;
389 fInfocl[1][jd1][jd2] = icl;
390 clust[0][numcell] = jd1;
391 clust[1][numcell] = jd2;
393 fInfcl[0][cellcount] = icl;
394 fInfcl[1][cellcount] = jd1;
395 fInfcl[2][cellcount] = jd2;
398 // ---------------------------------------------------------------
399 // check adc count for neighbour's neighbours recursively and
400 // if nonzero, add these to the cluster.
401 // ---------------------------------------------------------------
402 for(i = 1; i < kndim;i++)
408 for(j = 0; j < 6 ; j++)
410 jd1 = id1 + neibx[j];
411 jd2 = id2 + neiby[j];
412 if( (jd1 >= 0 && jd1 < kNDIMX) &&
413 (jd2 >= 0 && jd2 < kNDIMY) &&
414 fInfocl[0][jd1][jd2] == 0 )
416 fInfocl[0][jd1][jd2] = 2;
417 fInfocl[1][jd1][jd2] = icl;
419 clust[0][numcell] = jd1;
420 clust[1][numcell] = jd2;
422 fInfcl[0][cellcount] = icl;
423 fInfcl[1][cellcount] = jd1;
424 fInfcl[2][cellcount] = jd2;
433 // ------------------------------------------------------------------------ //
434 void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
436 // Does the refining of clusters
437 // Takes the big patch and does gaussian fitting and
438 // finds out the more refined clusters
441 AliPMDcludata *pmdcludata = 0;
443 const Int_t kNmaxCell = 19; // # of cells surrounding a cluster center
445 Int_t ndim = incr + 1;
450 Int_t i, j, k, i1, i2, id, icl, itest,ihld, ig, nsupcl,clno, t1, t2;
452 Double_t x1, y1, z1, x2, y2, z2, rr;
454 ncl = new Int_t [ndim];
455 clxy = new Int_t [kNmaxCell];
458 for(i = 0; i<ndim; i++)
461 if (i < 6) clusdata[i] = 0.;
462 if (i < kNmaxCell) clxy[i] = 0;
465 // clno counts the final clusters
466 // nsupcl = # of superclusters; ncl[i]= # of cells in supercluster i
467 // x, y and z store (x,y) coordinates of and energy deposited in a cell
468 // xc, yc store (x,y) coordinates of the cluster center
469 // zc stores the energy deposited in a cluster
470 // rc is cluster radius
475 for(i = 0; i <= incr; i++)
477 if(fInfcl[0][i] != nsupcl)
483 AliWarning("RefClust: Too many superclusters!");
490 AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
495 for(i = 0; i <= nsupcl; i++)
503 AliWarning("RefClust: Too many clusters! more than 4608");
510 i12 = i1 + i2*kNDIMX;
512 clusdata[0] = fCoord[0][i1][i2];
513 clusdata[1] = fCoord[1][i1][i2];
514 clusdata[2] = edepcell[i12];
519 clxy[0] = i1*10000 + i2;
521 for(Int_t icltr = 1; icltr < kNmaxCell; icltr++)
526 pmdcludata = new AliPMDcludata(clusdata,clxy);
527 fPMDclucont->Add(pmdcludata);
535 AliWarning("RefClust: Too many clusters! more than 4608");
541 i12 = i1 + i2*kNDIMX;
543 x1 = fCoord[0][i1][i2];
544 y1 = fCoord[1][i1][i2];
547 clxy[0] = i1*10000 + i2;
553 i22 = i1 + i2*kNDIMX;
554 x2 = fCoord[0][i1][i2];
555 y2 = fCoord[1][i1][i2];
558 clxy[1] = i1*10000 + i2;
561 for(Int_t icltr = 2; icltr < kNmaxCell; icltr++)
566 clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
567 clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
572 pmdcludata = new AliPMDcludata(clusdata,clxy);
573 fPMDclucont->Add(pmdcludata);
578 Int_t *iord, *tc, *t;
579 Double_t *x, *y, *z, *xc, *yc, *zc;
581 iord = new Int_t [ncl[i]+1];
582 tc = new Int_t [ncl[i]+1];
583 t = new Int_t [ncl[i]+1];
585 x = new Double_t [ncl[i]+1];
586 y = new Double_t [ncl[i]+1];
587 z = new Double_t [ncl[i]+1];
588 xc = new Double_t [ncl[i]+1];
589 yc = new Double_t [ncl[i]+1];
590 zc = new Double_t [ncl[i]+1];
592 for( k = 0; k < ncl[i]+1; k++)
605 // super-cluster of more than two cells - broken up into smaller
606 // clusters gaussian centers computed. (peaks separated by > 1 cell)
607 // Begin from cell having largest energy deposited This is first
611 i12 = i1 + i2*kNDIMX;
613 x[0] = fCoord[0][i1][i2];
614 y[0] = fCoord[1][i1][i2];
615 z[0] = edepcell[i12];
616 t[0] = i1*10000 + i2;
620 for(j = 1; j <= ncl[i]; j++)
625 i12 = i1 + i2*kNDIMX;
628 x[j] = fCoord[0][i1][i2];
629 y[j] = fCoord[1][i1][i2];
630 z[j] = edepcell[i12];
631 t[j] = i1*10000 + i2;
635 // arranging cells within supercluster in decreasing order
637 for(j = 1;j <= ncl[i]; j++)
641 for(i1 = 0; i1 < j; i1++)
643 if(itest == 0 && z[iord[i1]] < z[ihld])
646 for(i2 = j-1; i2 >= i1; i2--)
648 iord[i2+1] = iord[i2];
657 // Clustering algorithm returns from here for PP collisions
658 // for pp, only the output of crude clusterng is taken
659 // sigx and sigy are not calculated at this moment
661 Double_t supx=0., supy=0., supz=0.;
663 for(j = 0;j <= ncl[i]; j++)
665 supx += x[iord[j]]*z[iord[j]];
666 supy += y[iord[j]]*z[iord[j]];
670 clxy[j] = t[iord[j]];
676 for(Int_t ncel = ncl[i] + 1; ncel < kNmaxCell; ncel ++ )
681 clusdata[0] = supx/supz;
682 clusdata[1] = supy/supz;
684 clusdata[3] = ncl[i]+1;
687 pmdcludata = new AliPMDcludata(clusdata,clxy);
688 fPMDclucont->Add(pmdcludata);
691 /* MODIFICATION PART STARTS (Tapan July 2008)
692 iord[0] is the cell with highest ADC in the crude-cluster
693 ig is the number of local maxima in the crude-cluster
694 For the higest peak we make ig=0 which means first local maximum.
695 Next we go down in terms of the ADC sequence and find out if any
696 more of the cells form local maxima. The definition of local
697 maxima is that all its neighbours are of less ADC compared to it.
702 // This part is to split the supercluster
709 Int_t ivalid = 0, icount = 0;
711 for(j=1;j<=ncl[i];j++)
717 rr=Distance(x1,y1,xc[ig],yc[ig]);
719 // Check the cells which are outside the neighbours (rr>1.2)
724 for(Int_t j1=1;j1<j;j1++)
727 Float_t rr1=Distance(x1,y1,x[iord[j1]],y[iord[j1]]);
728 if(rr1>1.2) ivalid++;
730 if(ivalid == icount && z1>0.5*zc[ig])
743 // We use simple Gaussian weighting. (Tapan Jan 2005)
744 // compute the number of cells belonging to each cluster.
745 // cell can be shared between several clusters
746 // in the ratio of cluster energy deposition
748 // (1) number of cells belonging to a cluster (ig) and
749 // (2) total ADC of the cluster (ig)
750 // (3) x and y positions of the cluster
757 Double_t *totaladc, *totaladc2, *ncell,*weight;
758 Double_t *xclust, *yclust, *sigxclust, *sigyclust;
759 Double_t *ax, *ay, *ax2, *ay2;
762 status = new Int_t [ncl[i]+1];
763 cellXY = new Int_t *[ncl[i]+1];
765 cellCount = new Int_t [ig+1];
766 totaladc = new Double_t [ig+1];
767 totaladc2 = new Double_t [ig+1];
768 ncell = new Double_t [ig+1];
769 weight = new Double_t [ig+1];
770 xclust = new Double_t [ig+1];
771 yclust = new Double_t [ig+1];
772 sigxclust = new Double_t [ig+1];
773 sigyclust = new Double_t [ig+1];
774 ax = new Double_t [ig+1];
775 ay = new Double_t [ig+1];
776 ax2 = new Double_t [ig+1];
777 ay2 = new Double_t [ig+1];
779 for(j = 0; j < ncl[i]+1; j++)
782 cellXY[j] = new Int_t[ig+1];
785 for(Int_t kcl = 0; kcl < ig+1; kcl++)
800 for(j = 0; j < ncl[i]+1; j++)
805 Double_t sumweight, gweight;
807 for(j = 0;j <= ncl[i]; j++)
814 for(Int_t kcl=0; kcl<=ig; kcl++)
818 rr = Distance(x1,y1,x2,y2);
825 totaladc2[kcl] = z1*z1;
835 for(j = 0; j <= ncl[i]; j++)
848 for(Int_t kcl = 0; kcl <= ig; kcl++)
852 rr = Distance(x1,y1,x2,y2);
853 gweight = exp(-(rr*rr)/(2*(1.2*1.2)));
854 weight[kcl] = zc[kcl] * gweight;
855 sumweight = sumweight + weight[kcl];
857 if(weight[kcl] > max)
864 cellXY[cellCount[maxweight]][maxweight] = iord[j];
866 cellCount[maxweight]++;
870 totaladc[maxweight] += z1;
871 ax[maxweight] += x1*z1;
872 ay[maxweight] += y1*z1;
873 totaladc2[maxweight] += z1*z1;
874 ax2[maxweight] += z1*(x1-x2)*(x1-x2);
875 ay2[maxweight] += z1*(y1-y2)*(y1-y2);
881 for(Int_t kcl = 0; kcl <= ig; kcl++)
883 if(totaladc[kcl] > 0.)
885 xclust[kcl] = (ax[kcl])/ totaladc[kcl];
886 yclust[kcl] = (ay[kcl])/ totaladc[kcl];
889 Float_t sqtotadc = totaladc[kcl]*totaladc[kcl];
890 if(totaladc2[kcl] >= sqtotadc)
892 sigxclust[kcl] = 0.25;
893 sigyclust[kcl] = 0.25;
897 sigxclust[kcl] = (totaladc[kcl]/(sqtotadc-totaladc2[kcl]))*ax2[kcl];
898 sigyclust[kcl] = (totaladc[kcl]/(sqtotadc-totaladc2[kcl]))*ay2[kcl];
902 for(j = 0; j < cellCount[kcl]; j++) clno++;
906 AliWarning("RefClust: Too many clusters! more than 4608");
909 clusdata[0] = xclust[kcl];
910 clusdata[1] = yclust[kcl];
911 clusdata[2] = totaladc[kcl];
912 clusdata[3] = ncell[kcl];
914 if(sigxclust[kcl] > sigyclust[kcl])
916 clusdata[4] = TMath::Sqrt(sigxclust[kcl]);
917 clusdata[5] = TMath::Sqrt(sigyclust[kcl]);
921 clusdata[4] = TMath::Sqrt(sigyclust[kcl]);
922 clusdata[5] = TMath::Sqrt(sigxclust[kcl]);
928 for (Int_t ii = 0; ii < cellCount[kcl]; ii++)
932 clxy[Ncell] = t[cellXY[ii][kcl]];
937 pmdcludata = new AliPMDcludata(clusdata,clxy);
938 fPMDclucont->Add(pmdcludata);
941 for(Int_t jj = 0; jj < ncl[i]+1; jj++) delete [] cellXY[jj];
975 // ------------------------------------------------------------------------ //
976 Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1,
977 Double_t x2, Double_t y2)
979 return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
981 // ------------------------------------------------------------------------ //
982 void AliPMDClusteringV1::SetEdepCut(Float_t decut)
986 // ------------------------------------------------------------------------ //
987 void AliPMDClusteringV1::SetClusteringParam(Int_t cluspar)
989 fClusParam = cluspar;
991 // ------------------------------------------------------------------------ //