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[11424];
125 // call the isolated cell search method
127 FindIsoCell(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];
186 TMath::Sort((Int_t)kNMX,edepcell,iord1);// order the data
187 cutoff = fCutoff; // cutoff to discard cells
190 for(i = 0;i < kNMX; i++)
196 if(edepcell[i] > cutoff )
202 AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
204 if (nmx1 == 0) nmx1 = 1;
206 AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
209 incr = CrClust(ave, cutoff, nmx1,iord1, edepcell );
210 RefClust(incr,edepcell);
211 Int_t nentries1 = fPMDclucont->GetEntries();
212 AliDebug(1,Form("Detector Plane = %d Serial Module No = %d Number of clusters = %d",idet, ismn, nentries1));
213 AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
215 for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
217 AliPMDcludata *pmdcludata =
218 (AliPMDcludata*)fPMDclucont->UncheckedAt(ient1);
219 Float_t cluXC = pmdcludata->GetClusX();
220 Float_t cluYC = pmdcludata->GetClusY();
221 Float_t cluADC = pmdcludata->GetClusADC();
222 Float_t cluCELLS = pmdcludata->GetClusCells();
223 Float_t cluSIGX = pmdcludata->GetClusSigmaX();
224 Float_t cluSIGY = pmdcludata->GetClusSigmaY();
226 Float_t cluY0 = ktwobysqrt3*cluYC;
227 Float_t cluX0 = cluXC - cluY0/2.;
230 // Cluster X centroid is back transformed
235 clusdata[0] = cluX0 - (24-1) + cluY0/2.;
237 else if ( ismn >= 12 && ismn <= 23)
239 clusdata[0] = cluX0 - (48-1) + cluY0/2.;
243 clusdata[2] = cluADC;
244 clusdata[3] = cluCELLS;
245 clusdata[4] = cluSIGX;
246 clusdata[5] = cluSIGY;
249 // Cells associated with a cluster
252 for (Int_t ihit = 0; ihit < kNmaxCell; ihit++)
254 Int_t cellrow = pmdcludata->GetCellXY(ihit)/10000;
255 Int_t cellcol = pmdcludata->GetCellXY(ihit)%10000;
259 celldataX[ihit] = cellrow - (24-1) + int(cellcol/2.);
261 else if (ismn >= 12 && ismn <= 23)
263 celldataX[ihit] = cellrow - (48-1) + int(cellcol/2.);
266 celldataY[ihit] = cellcol;
268 Int_t irow = celldataX[ihit];
269 Int_t icol = celldataY[ihit];
273 if ((irow >= 0 && irow < 96) && (icol >= 0 && icol < 48))
275 celldataTr[ihit] = celltrack[icol][irow];
276 celldataPid[ihit] = cellpid[icol][irow];
277 celldataAdc[ihit] = (Float_t) celladc[icol][irow];
281 celldataTr[ihit] = -1;
282 celldataPid[ihit] = -1;
283 celldataAdc[ihit] = -1;
286 else if (ismn >= 12 && ismn < 24)
288 if ((irow >= 0 && irow < 48) && (icol >= 0 && icol < 96))
290 celldataTr[ihit] = celltrack[irow][icol];
291 celldataPid[ihit] = cellpid[irow][icol];
292 celldataAdc[ihit] = (Float_t) celladc[irow][icol];
297 celldataTr[ihit] = -1;
298 celldataPid[ihit] = -1;
299 celldataAdc[ihit] = -1;
305 pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY,
306 celldataTr, celldataPid, celldataAdc);
310 fPMDclucont->Delete();
312 // ------------------------------------------------------------------------ //
313 Int_t AliPMDClusteringV1::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
314 Int_t iord1[], Double_t edepcell[])
316 // Does crude clustering
317 // Finds out only the big patch by just searching the
320 const Int_t kndim = 4609;
321 Int_t i,j,k,id1,id2,icl, numcell, clust[2][kndim];
322 Int_t jd1,jd2, icell, cellcount;
323 static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
325 AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
327 for (j = 0; j < kNDIMX; j++)
329 for(k = 0; k < kNDIMY; k++)
331 fInfocl[0][j][k] = 0;
332 fInfocl[1][j][k] = 0;
335 for(i=0; i < kNMX; i++)
343 if(edepcell[j] <= cutoff)
345 fInfocl[0][id1][id2] = -1;
349 // ---------------------------------------------------------------
350 // crude clustering begins. Start with cell having largest adc
351 // count and loop over the cells in descending order of adc count
352 // ---------------------------------------------------------------
357 for(icell = 0; icell <= nmx1; icell++)
363 if(fInfocl[0][id1][id2] == 0 )
368 fInfocl[0][id1][id2] = 1;
369 fInfocl[1][id1][id2] = icl;
370 fInfcl[0][cellcount] = icl;
371 fInfcl[1][cellcount] = id1;
372 fInfcl[2][cellcount] = id2;
374 clust[0][numcell] = id1;
375 clust[1][numcell] = id2;
377 for(i = 1; i < kndim; i++)
381 // ---------------------------------------------------------------
382 // check for adc count in neib. cells. If ne 0 put it in this clust
383 // ---------------------------------------------------------------
384 for(i = 0; i < 6; i++)
386 jd1 = id1 + neibx[i];
387 jd2 = id2 + neiby[i];
388 if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
389 fInfocl[0][jd1][jd2] == 0)
392 fInfocl[0][jd1][jd2] = 2;
393 fInfocl[1][jd1][jd2] = icl;
394 clust[0][numcell] = jd1;
395 clust[1][numcell] = jd2;
397 fInfcl[0][cellcount] = icl;
398 fInfcl[1][cellcount] = jd1;
399 fInfcl[2][cellcount] = jd2;
402 // ---------------------------------------------------------------
403 // check adc count for neighbour's neighbours recursively and
404 // if nonzero, add these to the cluster.
405 // ---------------------------------------------------------------
406 for(i = 1; i < kndim;i++)
412 for(j = 0; j < 6 ; j++)
414 jd1 = id1 + neibx[j];
415 jd2 = id2 + neiby[j];
416 if( (jd1 >= 0 && jd1 < kNDIMX) &&
417 (jd2 >= 0 && jd2 < kNDIMY) &&
418 fInfocl[0][jd1][jd2] == 0 )
420 fInfocl[0][jd1][jd2] = 2;
421 fInfocl[1][jd1][jd2] = icl;
423 clust[0][numcell] = jd1;
424 clust[1][numcell] = jd2;
426 fInfcl[0][cellcount] = icl;
427 fInfcl[1][cellcount] = jd1;
428 fInfcl[2][cellcount] = jd2;
437 // ------------------------------------------------------------------------ //
438 void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
440 // Does the refining of clusters
441 // Takes the big patch and does gaussian fitting and
442 // finds out the more refined clusters
445 AliPMDcludata *pmdcludata = 0;
447 const Int_t kNmaxCell = 19; // # of cells surrounding a cluster center
449 Int_t ndim = incr + 1;
454 Int_t i, j, k, i1, i2, id, icl, itest,ihld, ig, nsupcl,clno, t1, t2;
456 Double_t x1, y1, z1, x2, y2, z2, rr;
458 ncl = new Int_t [ndim];
459 clxy = new Int_t [kNmaxCell];
462 for(i = 0; i<ndim; i++)
465 if (i < 6) clusdata[i] = 0.;
466 if (i < kNmaxCell) clxy[i] = 0;
469 // clno counts the final clusters
470 // nsupcl = # of superclusters; ncl[i]= # of cells in supercluster i
471 // x, y and z store (x,y) coordinates of and energy deposited in a cell
472 // xc, yc store (x,y) coordinates of the cluster center
473 // zc stores the energy deposited in a cluster
474 // rc is cluster radius
479 for(i = 0; i <= incr; i++)
481 if(fInfcl[0][i] != nsupcl)
487 AliWarning("RefClust: Too many superclusters!");
494 AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
499 for(i = 0; i <= nsupcl; i++)
507 AliWarning("RefClust: Too many clusters! more than 4608");
514 i12 = i1 + i2*kNDIMX;
516 clusdata[0] = fCoord[0][i1][i2];
517 clusdata[1] = fCoord[1][i1][i2];
518 clusdata[2] = edepcell[i12];
523 clxy[0] = i1*10000 + i2;
525 for(Int_t icltr = 1; icltr < kNmaxCell; icltr++)
530 pmdcludata = new AliPMDcludata(clusdata,clxy);
531 fPMDclucont->Add(pmdcludata);
539 AliWarning("RefClust: Too many clusters! more than 4608");
545 i12 = i1 + i2*kNDIMX;
547 x1 = fCoord[0][i1][i2];
548 y1 = fCoord[1][i1][i2];
551 clxy[0] = i1*10000 + i2;
557 i22 = i1 + i2*kNDIMX;
558 x2 = fCoord[0][i1][i2];
559 y2 = fCoord[1][i1][i2];
562 clxy[1] = i1*10000 + i2;
565 for(Int_t icltr = 2; icltr < kNmaxCell; icltr++)
570 clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
571 clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
576 pmdcludata = new AliPMDcludata(clusdata,clxy);
577 fPMDclucont->Add(pmdcludata);
582 Int_t *iord, *tc, *t;
583 Double_t *x, *y, *z, *xc, *yc, *zc;
585 iord = new Int_t [ncl[i]+1];
586 tc = new Int_t [ncl[i]+1];
587 t = new Int_t [ncl[i]+1];
589 x = new Double_t [ncl[i]+1];
590 y = new Double_t [ncl[i]+1];
591 z = new Double_t [ncl[i]+1];
592 xc = new Double_t [ncl[i]+1];
593 yc = new Double_t [ncl[i]+1];
594 zc = new Double_t [ncl[i]+1];
596 for( k = 0; k < ncl[i]+1; k++)
609 // super-cluster of more than two cells - broken up into smaller
610 // clusters gaussian centers computed. (peaks separated by > 1 cell)
611 // Begin from cell having largest energy deposited This is first
615 i12 = i1 + i2*kNDIMX;
617 x[0] = fCoord[0][i1][i2];
618 y[0] = fCoord[1][i1][i2];
619 z[0] = edepcell[i12];
620 t[0] = i1*10000 + i2;
624 for(j = 1; j <= ncl[i]; j++)
629 i12 = i1 + i2*kNDIMX;
632 x[j] = fCoord[0][i1][i2];
633 y[j] = fCoord[1][i1][i2];
634 z[j] = edepcell[i12];
635 t[j] = i1*10000 + i2;
639 // arranging cells within supercluster in decreasing order
641 for(j = 1;j <= ncl[i]; j++)
645 for(i1 = 0; i1 < j; i1++)
647 if(itest == 0 && z[iord[i1]] < z[ihld])
650 for(i2 = j-1; i2 >= i1; i2--)
652 iord[i2+1] = iord[i2];
658 /* MODIFICATION PART STARTS (Tapan July 2008)
659 iord[0] is the cell with highest ADC in the crude-cluster
660 ig is the number of local maxima in the crude-cluster
661 For the higest peak we make ig=0 which means first local maximum.
662 Next we go down in terms of the ADC sequence and find out if any
663 more of the cells form local maxima. The definition of local
664 maxima is that all its neighbours are of less ADC compared to it.
671 Int_t ivalid = 0, icount = 0;
673 for(j=1;j<=ncl[i];j++)
679 rr=Distance(x1,y1,xc[ig],yc[ig]);
681 // Check the cells which are outside the neighbours (rr>1.2)
686 for(Int_t j1=1;j1<j;j1++)
689 Float_t rr1=Distance(x1,y1,x[iord[j1]],y[iord[j1]]);
690 if(rr1>1.2) ivalid++;
692 if(ivalid == icount && z1>0.5*zc[ig])
705 // We use simple Gaussian weighting. (Tapan Jan 2005)
706 // compute the number of cells belonging to each cluster.
707 // cell can be shared between several clusters
708 // in the ratio of cluster energy deposition
710 // (1) number of cells belonging to a cluster (ig) and
711 // (2) total ADC of the cluster (ig)
712 // (3) x and y positions of the cluster
719 Double_t *totaladc, *totaladc2, *ncell,*weight;
720 Double_t *xclust, *yclust, *sigxclust, *sigyclust;
721 Double_t *ax, *ay, *ax2, *ay2;
724 status = new Int_t [ncl[i]+1];
725 cellXY = new Int_t *[ncl[i]+1];
727 cellCount = new Int_t [ig+1];
728 totaladc = new Double_t [ig+1];
729 totaladc2 = new Double_t [ig+1];
730 ncell = new Double_t [ig+1];
731 weight = new Double_t [ig+1];
732 xclust = new Double_t [ig+1];
733 yclust = new Double_t [ig+1];
734 sigxclust = new Double_t [ig+1];
735 sigyclust = new Double_t [ig+1];
736 ax = new Double_t [ig+1];
737 ay = new Double_t [ig+1];
738 ax2 = new Double_t [ig+1];
739 ay2 = new Double_t [ig+1];
741 for(j = 0; j < ncl[i]+1; j++)
744 cellXY[j] = new Int_t[ig+1];
747 for(Int_t kcl = 0; kcl < ig+1; kcl++)
762 for(j = 0; j < ncl[i]+1; j++)
767 Double_t sumweight, gweight;
769 for(j = 0;j <= ncl[i]; j++)
776 for(Int_t kcl=0; kcl<=ig; kcl++)
780 rr = Distance(x1,y1,x2,y2);
787 totaladc2[kcl] = z1*z1;
797 for(j = 0; j <= ncl[i]; j++)
810 for(Int_t kcl = 0; kcl <= ig; kcl++)
814 rr = Distance(x1,y1,x2,y2);
815 gweight = exp(-(rr*rr)/(2*(1.2*1.2)));
816 weight[kcl] = zc[kcl] * gweight;
817 sumweight = sumweight + weight[kcl];
819 if(weight[kcl] > max)
826 cellXY[cellCount[maxweight]][maxweight] = iord[j];
828 cellCount[maxweight]++;
832 totaladc[maxweight] += z1;
833 ax[maxweight] += x1*z1;
834 ay[maxweight] += y1*z1;
835 totaladc2[maxweight] += z1*z1;
836 ax2[maxweight] += z1*(x1-x2)*(x1-x2);
837 ay2[maxweight] += z1*(y1-y2)*(y1-y2);
843 for(Int_t kcl = 0; kcl <= ig; kcl++)
846 if(totaladc[kcl] > 0.)
848 xclust[kcl] = (ax[kcl])/ totaladc[kcl];
849 yclust[kcl] = (ay[kcl])/ totaladc[kcl];
852 Float_t sqtotadc = totaladc[kcl]*totaladc[kcl];
853 if(totaladc2[kcl] >= sqtotadc)
855 sigxclust[kcl] = 0.25;
856 sigyclust[kcl] = 0.25;
860 sigxclust[kcl] = (totaladc[kcl]/(sqtotadc-totaladc2[kcl]))*ax2[kcl];
861 sigyclust[kcl] = (totaladc[kcl]/(sqtotadc-totaladc2[kcl]))*ay2[kcl];
865 for(j = 0; j < cellCount[kcl]; j++) clno++;
869 AliWarning("RefClust: Too many clusters! more than 4608");
872 clusdata[0] = xclust[kcl];
873 clusdata[1] = yclust[kcl];
874 clusdata[2] = totaladc[kcl];
875 clusdata[3] = ncell[kcl];
878 if(sigxclust[kcl] > sigyclust[kcl])
880 clusdata[4] = TMath::Sqrt(sigxclust[kcl]);
881 clusdata[5] = TMath::Sqrt(sigyclust[kcl]);
885 clusdata[4] = TMath::Sqrt(sigyclust[kcl]);
886 clusdata[5] = TMath::Sqrt(sigxclust[kcl]);
892 for (Int_t ii = 0; ii < cellCount[kcl]; ii++)
896 clxy[Ncell] = t[cellXY[ii][kcl]];
901 pmdcludata = new AliPMDcludata(clusdata,clxy);
902 fPMDclucont->Add(pmdcludata);
916 for(Int_t jj = 0; jj < ncl[i]+1; jj++) delete [] cellXY[jj];
936 // ------------------------------------------------------------------------ //
937 Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1,
938 Double_t x2, Double_t y2)
940 return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
942 // ------------------------------------------------------------------------ //
943 void AliPMDClusteringV1::FindIsoCell(Int_t idet, Int_t ismn, Double_t celladc[][96], TObjArray *pmdisocell)
945 // Does isolated cell search for offline calibration
947 AliPMDisocell *isocell = 0;
949 const Int_t kMaxRow = 48;
950 const Int_t kMaxCol = 96;
951 const Int_t kCellNeighbour = 6;
955 Int_t neibx[6] = {1,0,-1,-1,0,1};
956 Int_t neiby[6] = {0,1,1,0,-1,-1};
959 for(Int_t irow = 0; irow < kMaxRow; irow++)
961 for(Int_t icol = 0; icol < kMaxCol; icol++)
963 if(celladc[irow][icol] > 0)
966 for(Int_t ii = 0; ii < kCellNeighbour; ii++)
968 id1 = irow + neibx[ii];
969 jd1 = icol + neiby[ii];
970 if (id1 < 0) id1 = 0;
971 if (id1 > kMaxRow-1) id1 = kMaxRow - 1;
972 if (jd1 < 0) jd1 = 0;
973 if (jd1 > kMaxCol-1) jd1 = kMaxCol - 1;
974 Float_t adc = (Float_t) celladc[id1][jd1];
978 if(isocount == kCellNeighbour)
980 Float_t cadc = (Float_t) celladc[irow][icol];
982 isocell = new AliPMDisocell(idet,ismn,irow,icol,cadc);
983 pmdisocell->Add(isocell);
987 } // neigh cell cond.
994 // ------------------------------------------------------------------------ //
995 void AliPMDClusteringV1::SetEdepCut(Float_t decut)
999 // ------------------------------------------------------------------------ //