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"
57 ClassImp(AliPMDClusteringV1)
59 const Double_t AliPMDClusteringV1::fgkSqroot3by2=0.8660254; // sqrt(3.)/2.
61 AliPMDClusteringV1::AliPMDClusteringV1():
62 fPMDclucont(new TObjArray()),
65 for(Int_t i = 0; i < kNDIMX; i++)
67 for(Int_t j = 0; j < kNDIMY; j++)
69 fCoord[0][i][j] = i+j/2.;
70 fCoord[1][i][j] = fgkSqroot3by2*j;
74 // ------------------------------------------------------------------------ //
75 AliPMDClusteringV1::AliPMDClusteringV1(const AliPMDClusteringV1& pmdclv1):
76 AliPMDClustering(pmdclv1),
81 AliError("Copy constructor not allowed ");
84 // ------------------------------------------------------------------------ //
85 AliPMDClusteringV1 &AliPMDClusteringV1::operator=(const AliPMDClusteringV1& /*pmdclv1*/)
88 AliError("Assignment operator not allowed ");
91 // ------------------------------------------------------------------------ //
92 AliPMDClusteringV1::~AliPMDClusteringV1()
96 // ------------------------------------------------------------------------ //
97 void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn,
98 Double_t celladc[48][96], TObjArray *pmdcont)
100 // main function to call other necessary functions to do clustering
103 AliPMDcluster *pmdcl = 0;
105 Int_t i, j, nmx1, incr, id, jd;
106 Int_t celldataX[15], celldataY[15];
108 Double_t cutoff, ave;
109 Double_t edepcell[kNMX];
112 Double_t *cellenergy = new Double_t [11424];// Ajay
114 const float ktwobysqrt3 = 1.1547; // 2./sqrt(3.)
116 // ndimXr and ndimYr are different because of different module size
126 else if (ismn >= 12 && ismn <= 23)
132 for (i =0; i < 11424; i++)
139 for (i = 0; i < kNDIMX; i++)
141 for (j = 0; j < kNDIMY; j++)
143 fCellTrNo[i][j] = -1;
149 for (id = 0; id < ndimXr; id++)
151 for (jd = 0; jd < ndimYr; jd++)
154 i = id+(ndimYr/2-1)-(jd/2);
156 Int_t ij = i + j*kNDIMX;
160 cellenergy[ij] = celladc[jd][id];//Ajay
161 fCellTrNo[i][j] = jd*10000+id; // for association
163 else if (ismn >= 12 && ismn <= 23)
165 cellenergy[ij] = celladc[id][jd];//Ajay
166 fCellTrNo[i][j] = id*10000+jd; // for association
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
227 clusdata[0] = cluX0 - (24-1) + cluY0/2.;
229 else if ( ismn >= 12 && ismn <= 23)
231 clusdata[0] = cluX0 - (48-1) + cluY0/2.;
235 clusdata[2] = cluADC;
236 clusdata[3] = cluCELLS;
237 clusdata[4] = cluSIGX;
238 clusdata[5] = cluSIGY;
241 // Cells associated with a cluster
244 for (Int_t ihit = 0; ihit < 15; ihit++)
248 celldataX[ihit] = pmdcludata->GetCellXY(ihit)%10000;
249 celldataY[ihit] = pmdcludata->GetCellXY(ihit)/10000;
251 else if (ismn >= 12 && ismn <= 23)
253 celldataX[ihit] = pmdcludata->GetCellXY(ihit)/10000;
254 celldataY[ihit] = pmdcludata->GetCellXY(ihit)%10000;
257 pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY);
261 fPMDclucont->Delete();
264 // ------------------------------------------------------------------------ //
265 Int_t AliPMDClusteringV1::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
266 Int_t iord1[], Double_t edepcell[])
268 // Does crude clustering
269 // Finds out only the big patch by just searching the
272 Int_t i,j,k,id1,id2,icl, numcell, clust[2][5000];
273 Int_t jd1,jd2, icell, cellcount;
274 static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
276 AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
278 for (j = 0; j < kNDIMX; j++)
280 for(k = 0; k < kNDIMY; k++)
282 fInfocl[0][j][k] = 0;
283 fInfocl[1][j][k] = 0;
286 for(i=0; i < kNMX; i++)
294 if(edepcell[j] <= cutoff)
296 fInfocl[0][id1][id2] = -1;
300 // ---------------------------------------------------------------
301 // crude clustering begins. Start with cell having largest adc
302 // count and loop over the cells in descending order of adc count
303 // ---------------------------------------------------------------
308 for(icell = 0; icell <= nmx1; icell++)
314 if(fInfocl[0][id1][id2] == 0 )
319 fInfocl[0][id1][id2] = 1;
320 fInfocl[1][id1][id2] = icl;
321 fInfcl[0][cellcount] = icl;
322 fInfcl[1][cellcount] = id1;
323 fInfcl[2][cellcount] = id2;
325 clust[0][numcell] = id1;
326 clust[1][numcell] = id2;
328 for(i = 1; i < 5000; i++)
332 // ---------------------------------------------------------------
333 // check for adc count in neib. cells. If ne 0 put it in this clust
334 // ---------------------------------------------------------------
335 for(i = 0; i < 6; i++)
337 jd1 = id1 + neibx[i];
338 jd2 = id2 + neiby[i];
339 if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
340 fInfocl[0][jd1][jd2] == 0)
343 fInfocl[0][jd1][jd2] = 2;
344 fInfocl[1][jd1][jd2] = icl;
345 clust[0][numcell] = jd1;
346 clust[1][numcell] = jd2;
348 fInfcl[0][cellcount] = icl;
349 fInfcl[1][cellcount] = jd1;
350 fInfcl[2][cellcount] = jd2;
353 // ---------------------------------------------------------------
354 // check adc count for neighbour's neighbours recursively and
355 // if nonzero, add these to the cluster.
356 // ---------------------------------------------------------------
357 for(i = 1; i < 5000;i++)
363 for(j = 0; j < 6 ; j++)
365 jd1 = id1 + neibx[j];
366 jd2 = id2 + neiby[j];
367 if( (jd1 >= 0 && jd1 < kNDIMX) &&
368 (jd2 >= 0 && jd2 < kNDIMY) &&
369 fInfocl[0][jd1][jd2] == 0 )
371 fInfocl[0][jd1][jd2] = 2;
372 fInfocl[1][jd1][jd2] = icl;
374 clust[0][numcell] = jd1;
375 clust[1][numcell] = jd2;
377 fInfcl[0][cellcount] = icl;
378 fInfcl[1][cellcount] = jd1;
379 fInfcl[2][cellcount] = jd2;
388 // ------------------------------------------------------------------------ //
389 void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
391 // Does the refining of clusters
392 // Takes the big patch and does gaussian fitting and
393 // finds out the more refined clusters
398 AliPMDcludata *pmdcludata = 0;
400 Int_t *cellCount = 0x0;
401 Int_t **cellXY = 0x0;
402 const Int_t kdim = 4609;
405 Int_t i, j, k, i1, i2, id, icl, itest;
407 Int_t ig, nsupcl,clno;
409 Int_t ncl[kdim], iord[kdim], lev1[kdim], lev2[kdim];
412 Double_t x1, y1, z1, x2, y2, z2, dist,rr,sum;
413 Double_t x[kdim], y[kdim], z[kdim];
414 Double_t xc[kdim], yc[kdim], zc[kdim], cells[kdim], rc[kdim];
417 for(i = 0; i<kdim; i++)
421 if (i < 6) clusdata[i] = 0.;
422 if (i < 15) clxy[i] = 0;
425 // clno counts the final clusters
426 // nsupcl = # of superclusters; ncl[i]= # of cells in supercluster i
427 // x, y and z store (x,y) coordinates of and energy deposited in a cell
428 // xc, yc store (x,y) coordinates of the cluster center
429 // zc stores the energy deposited in a cluster
430 // rc is cluster radius
435 for(i = 0; i <= incr; i++)
437 if(fInfcl[0][i] != nsupcl)
443 AliWarning("RefClust: Too many superclusters!");
450 AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
455 for(i = 0; i <= nsupcl; i++)
463 AliWarning("RefClust: Too many clusters! more than 5000");
470 i12 = i1 + i2*kNDIMX;
472 clusdata[0] = fCoord[0][i1][i2];
473 clusdata[1] = fCoord[1][i1][i2];
474 clusdata[2] = edepcell[i12];
479 clxy[0] = fCellTrNo[i1][i2]; //association
480 for(Int_t icltr = 1; icltr < 15; icltr++)
484 pmdcludata = new AliPMDcludata(clusdata,clxy);
485 fPMDclucont->Add(pmdcludata);
493 AliWarning("RefClust: Too many clusters! more than 5000");
499 i12 = i1 + i2*kNDIMX;
501 x1 = fCoord[0][i1][i2];
502 y1 = fCoord[1][i1][i2];
504 clxy[0] = fCellTrNo[i1][i2]; //asso
509 Int_t i22 = i1 + i2*kNDIMX;
510 x2 = fCoord[0][i1][i2];
511 y2 = fCoord[1][i1][i2];
513 clxy[1] = fCellTrNo[i1][i2]; //asso
514 for(Int_t icltr = 2; icltr < 15; icltr++)
519 clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
520 clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
525 pmdcludata = new AliPMDcludata(clusdata,clxy);
526 fPMDclucont->Add(pmdcludata);
532 // super-cluster of more than two cells - broken up into smaller
533 // clusters gaussian centers computed. (peaks separated by > 1 cell)
534 // Begin from cell having largest energy deposited This is first
538 i12 = i1 + i2*kNDIMX;
540 x[0] = fCoord[0][i1][i2];
541 y[0] = fCoord[1][i1][i2];
542 z[0] = edepcell[i12];
544 t[0] = fCellTrNo[i1][i2]; //asso
547 for(j = 1; j <= ncl[i]; j++)
552 i12 = i1 + i2*kNDIMX;
555 x[j] = fCoord[0][i1][i2];
556 y[j] = fCoord[1][i1][i2];
557 z[j] = edepcell[i12];
559 t[j] = fCellTrNo[i1][i2]; //asso
562 // arranging cells within supercluster in decreasing order
565 for(j = 1;j <= ncl[i]; j++)
569 for(i1 = 0; i1 < j; i1++)
571 if(itest == 0 && z[iord[i1]] < z[ihld])
574 for(i2 = j-1; i2 >= i1; i2--)
576 iord[i2+1] = iord[i2];
584 Int_t imaxdim = ncl[i] + 1;
585 TMath::Sort(imaxdim,z,iord);// order the data
588 // compute the number of Gaussians and their centers ( first
590 // centers must be separated by cells having smaller ener. dep.
591 // neighbouring centers should be either strong or well-separated
596 for(j = 1; j <= ncl[i]; j++)
601 for(k = 0; k <= ig; k++)
605 rr = Distance(x1,y1,x2,y2);
606 if(rr >= 1.1 && rr < 1.8 && z[iord[j]] > zc[k]/4.)itest++;
607 if(rr >= 1.8 && rr < 2.1 && z[iord[j]] > zc[k]/10.)itest++;
608 if( rr >= 2.1)itest++;
618 GaussFit(ncl[i], ig, x[0], y[0] ,z[0], xc[0], yc[0], zc[0], rc[0]);
620 // compute the number of cells belonging to each cluster.
621 // cell is shared between several clusters ( if they are equidistant
622 // from it ) in the ratio of cluster energy deposition
625 cellCount = new Int_t [ig+1];
626 cellXY = new Int_t *[jj];
627 for(Int_t ij = 0; ij < 15; ij++) cellXY[ij] = new Int_t [ig+1];
629 for(j = 0; j <= ig; j++)
637 for(j = 0; j <= ncl[i]; j++)
641 for(k = 0; k <= ig; k++)
643 dist = Distance(x[j], y[j], xc[k], yc[k]);
644 if(dist < TMath::Sqrt(3.) )
647 if (cellCount[k] < 15)
649 cellXY[cellCount[k]][k] = t[j];
676 for(k = 1; k <= lev1[0]; k++)
680 for(k = 1; k <= lev1[0]; k++)
682 cells[lev1[k]] += zc[lev1[k]]/sum;
695 for( k = 1; k <= lev2[0]; k++)
699 for(k = 1; k <= lev2[0]; k++)
701 cells[lev2[k]] += zc[lev2[k]]/sum;
708 // zero rest of the cell array
710 for( k = 0; k <= ig; k++)
712 for(Int_t icltr = cellCount[k]; icltr < 15; icltr++)
714 cellXY[icltr][k] = -1;
719 for(j = 0; j <= ig; j++)
724 AliWarning("RefClust: Too many clusters! more than 5000");
734 clusdata[3] = ncl[i];
738 clusdata[3] = cells[j];
742 for (Int_t ii=0; ii < 15; ii++)
744 clxy[ii] = cellXY[ii][j];
746 pmdcludata = new AliPMDcludata(clusdata,clxy);
747 fPMDclucont->Add(pmdcludata);
750 for(jj = 0; jj < 15; jj++)delete [] cellXY[jj];
755 // ------------------------------------------------------------------------ //
756 void AliPMDClusteringV1::GaussFit(Int_t ncell, Int_t nclust, Double_t &x,
757 Double_t &y ,Double_t &z, Double_t &xc,
758 Double_t &yc, Double_t &zc, Double_t &rc)
760 // Does gaussian fitting
763 const Int_t kdim = 4609;
764 Int_t i, j, i1, i2, novar, idd, jj;
765 Int_t neib[kdim][50];
767 Double_t sum, dx, dy, str, str1, aint, sum1, rr, dum;
768 Double_t x1, x2, y1, y2;
769 Double_t xx[kdim], yy[kdim], zz[kdim], xxc[kdim], yyc[kdim];
770 Double_t a[kdim], b[kdim], c[kdim], d[kdim], ha[kdim], hb[kdim];
771 Double_t hc[kdim], hd[kdim], zzc[kdim], rrc[kdim];
781 for(i = 0; i <= ncell; i++)
788 for(i=0; i<=nclust; i++)
796 for(i = 0; i <= nclust; i++)
798 zzc[i] = str/str1*zzc[i];
807 for(i = 0; i <= ncell; i++)
812 for(j = 0; j <= nclust; j++)
816 if(Distance(x1,y1,x2,y2) <= 3.)
825 for(i1 = 0; i1 <= ncell; i1++)
829 for(i2 = 1; i2 <= idd; i2++)
832 dx = xx[i1] - xxc[jj];
833 dy = yy[i1] - yyc[jj];
834 dum = rrc[j]*rrc[jj] + rr*rr;
835 aint += exp(-(dx*dx+dy*dy)/dum)*zzc[idd]*rr*rr/dum;
837 sum += (aint - zz[i1])*(aint - zz[i1])/str;
841 for(i = 0; i <= nclust; i++)
843 a[i] = xxc[i] + 0.6*(rnd.Uniform() - 0.5);
844 b[i] = yyc[i] + 0.6*(rnd.Uniform() - 0.5);
845 c[i] = zzc[i]*(1.+ ( rnd.Uniform() - 0.5)*0.2);
847 d[i] = rrc[i]*(1.+ ( rnd.Uniform() - 0.5)*0.1);
854 for(i = 0; i <= nclust; i++)
856 c[i] = c[i]*str/str1;
860 for(i1 = 0; i1 <= ncell; i1++)
864 for(i2 = 1; i2 <= idd; i2++)
869 dum = d[jj]*d[jj]+rr*rr;
870 aint += exp(-(dx*dx+dy*dy)/dum)*c[i2]*rr*rr/dum;
872 sum1 += (aint - zz[i1])*(aint - zz[i1])/str;
877 for(i2 = 0; i2 <= nclust; i2++)
886 for(j = 0; j <= nclust; j++)
894 // ------------------------------------------------------------------------ //
895 Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1,
896 Double_t x2, Double_t y2)
898 return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
900 // ------------------------------------------------------------------------ //
901 void AliPMDClusteringV1::SetEdepCut(Float_t decut)
905 // ------------------------------------------------------------------------ //