clustering parameter is set to decide the clustering algo
[u/mrichter/AliRoot.git] / PMD / AliPMDClusteringV2.cxx
CommitLineData
8c7250c5 1/***************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15
090026bf 16/* $Id$ */
17
8c7250c5 18//-----------------------------------------------------//
19// //
20// Source File : PMDClusteringV2.cxx //
21// //
22// clustering code for alice pmd //
23// //
24//-----------------------------------------------------//
25
26/* --------------------------------------------------------------------
27 Code developed by S. C. Phatak, Institute of Physics,
28 Bhubaneswar 751 005 ( phatak@iopb.res.in ) Given the energy deposited
29 ( or ADC value ) in each cell of supermodule ( pmd or cpv ), the code
30 builds up superclusters and breaks them into clusters. The input is
562718f9 31 in TObjarray and cluster information is in TObjArray.
32 integer clno gives total number of clusters in the supermodule.
33 fClusters is the global ( public ) variables.
8c7250c5 34 Others are local ( private ) to the code.
35 At the moment, the data is read for whole detector ( all supermodules
36 and pmd as well as cpv. This will have to be modify later )
37 LAST UPDATE : October 23, 2002
38-----------------------------------------------------------------------*/
39
090026bf 40#include <Riostream.h>
41#include <TMath.h>
8c7250c5 42#include <TObjArray.h>
562718f9 43#include <TArrayI.h>
8c7250c5 44
562718f9 45#include "AliPMDcludata.h"
8c7250c5 46#include "AliPMDcluster.h"
47#include "AliPMDClustering.h"
48#include "AliPMDClusteringV2.h"
49#include "AliLog.h"
50
51ClassImp(AliPMDClusteringV2)
52
53const Double_t AliPMDClusteringV2::fgkSqroot3by2=0.8660254; // sqrt(3.)/2.
54
55AliPMDClusteringV2::AliPMDClusteringV2():
562718f9 56 fPMDclucont(new TObjArray()),
4755c3f1 57 fCutoff(0.0),
58 fClusParam(0)
8c7250c5 59{
60 for(int i = 0; i < kNDIMX; i++)
61 {
62 for(int j = 0; j < kNDIMY; j++)
63 {
64 fCoord[0][i][j] = i+j/2.;
65 fCoord[1][i][j] = fgkSqroot3by2*j;
8c7250c5 66 }
67 }
68}
69// ------------------------------------------------------------------------ //
562718f9 70
71
72AliPMDClusteringV2::AliPMDClusteringV2(const AliPMDClusteringV2& pmdclv2):
73 AliPMDClustering(pmdclv2),
74 fPMDclucont(0),
4755c3f1 75 fCutoff(0),
76 fClusParam(0)
562718f9 77{
78 // copy constructor
79 AliError("Copy constructor not allowed ");
80
81}
82// ------------------------------------------------------------------------ //
83AliPMDClusteringV2 &AliPMDClusteringV2::operator=(const AliPMDClusteringV2& /*pmdclv2*/)
84{
85 // copy constructor
86 AliError("Assignment operator not allowed ");
87 return *this;
88}
89// ------------------------------------------------------------------------ //
8c7250c5 90AliPMDClusteringV2::~AliPMDClusteringV2()
91{
562718f9 92 delete fPMDclucont;
8c7250c5 93}
94// ------------------------------------------------------------------------ //
562718f9 95
96void AliPMDClusteringV2::DoClust(Int_t idet, Int_t ismn,
920e13db 97 Int_t celltrack[48][96],
98 Int_t cellpid[48][96],
99 Double_t celladc[48][96],
22bd512d 100 TObjArray *pmdcont)
8c7250c5 101{
102 // main function to call other necessary functions to do clustering
103 //
104 AliPMDcluster *pmdcl = 0;
105
2c1131dd 106 const Float_t ktwobysqrt3 = 1.1547; // 2./sqrt(3.)
107 const Int_t kNmaxCell = 19; // # of cells surrounding a cluster center
562718f9 108 Int_t i, j, nmx1, incr, id, jd;
109 Int_t ndimXr = 0;
110 Int_t ndimYr = 0;
2c1131dd 111 Int_t celldataX[kNmaxCell], celldataY[kNmaxCell];
920e13db 112 Int_t celldataTr[kNmaxCell], celldataPid[kNmaxCell];
b0e4d1e1 113 Float_t celldataAdc[kNmaxCell];
562718f9 114 Float_t clusdata[6];
8c7250c5 115 Double_t cutoff, ave;
562718f9 116 Double_t edepcell[kNMX];
8c7250c5 117
8c7250c5 118
119 if (ismn < 12)
120 {
121 ndimXr = 96;
122 ndimYr = 48;
123 }
124 else if (ismn >= 12 && ismn <= 23)
125 {
126 ndimXr = 48;
127 ndimYr = 96;
128 }
562718f9 129
78fc1b96 130 for (i =0; i < kNMX; i++)
8c7250c5 131 {
562718f9 132 edepcell[i] = 0.;
8c7250c5 133 }
562718f9 134
8c7250c5 135 for (id = 0; id < ndimXr; id++)
136 {
137 for (jd = 0; jd < ndimYr; jd++)
138 {
562718f9 139 j = jd;
140 i = id + (ndimYr/2-1) - (jd/2);
141 Int_t ij = i + j*kNDIMX;
8c7250c5 142 if (ismn < 12)
143 {
562718f9 144 edepcell[ij] = celladc[jd][id];
8c7250c5 145 }
146 else if (ismn >= 12 && ismn <= 23)
147 {
562718f9 148 edepcell[ij] = celladc[id][jd];
8c7250c5 149 }
150
151 }
152 }
153
562718f9 154 Int_t iord1[kNMX];
df4e6759 155 TMath::Sort((Int_t)kNMX,edepcell,iord1);// order the data
8c7250c5 156 cutoff = fCutoff; // cutoff used to discard cells having ener. dep.
562718f9 157 ave = 0.;
158 nmx1 = -1;
8c7250c5 159
562718f9 160 for(i = 0;i < kNMX; i++)
8c7250c5 161 {
562718f9 162 if(edepcell[i] > 0.)
163 {
164 ave += edepcell[i];
165 }
166 if(edepcell[i] > cutoff )
167 {
168 nmx1++;
169 }
8c7250c5 170 }
562718f9 171
8c7250c5 172 AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
562718f9 173
174 if (nmx1 == 0)
175 {
176 nmx1 = 1;
177 }
178 ave = ave/nmx1;
179
8c7250c5 180 AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
181 kNMX,ave));
8c7250c5 182
562718f9 183 incr = CrClust(ave, cutoff, nmx1,iord1, edepcell);
184 RefClust(incr,edepcell );
185
186 Int_t nentries1 = fPMDclucont->GetEntries();
187 AliDebug(1,Form("Detector Plane = %d Serial Module No = %d Number of clusters = %d",idet, ismn, nentries1));
188 AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
189 for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
8c7250c5 190 {
562718f9 191 AliPMDcludata *pmdcludata =
192 (AliPMDcludata*)fPMDclucont->UncheckedAt(ient1);
193 Float_t cluXC = pmdcludata->GetClusX();
194 Float_t cluYC = pmdcludata->GetClusY();
195 Float_t cluADC = pmdcludata->GetClusADC();
196 Float_t cluCELLS = pmdcludata->GetClusCells();
197 Float_t cluSIGX = pmdcludata->GetClusSigmaX();
198 Float_t cluSIGY = pmdcludata->GetClusSigmaY();
199
8c7250c5 200 Float_t cluY0 = ktwobysqrt3*cluYC;
201 Float_t cluX0 = cluXC - cluY0/2.;
562718f9 202
8c7250c5 203 //
204 // Cluster X centroid is back transformed
205 //
206 if (ismn < 12)
207 {
208 clusdata[0] = cluX0 - (24-1) + cluY0/2.;
209 }
562718f9 210 else if (ismn >= 12 && ismn <= 23)
8c7250c5 211 {
212 clusdata[0] = cluX0 - (48-1) + cluY0/2.;
213 }
214
562718f9 215 clusdata[1] = cluY0;
216 clusdata[2] = cluADC;
217 clusdata[3] = cluCELLS;
218 clusdata[4] = cluSIGX;
219 clusdata[5] = cluSIGY;
8c7250c5 220 //
221 // Cells associated with a cluster
222 //
2c1131dd 223 for (Int_t ihit = 0; ihit < kNmaxCell; ihit++)
8c7250c5 224 {
562718f9 225 Int_t dummyXY = pmdcludata->GetCellXY(ihit);
226
227 Int_t celldumY = dummyXY%10000;
228 Int_t celldumX = dummyXY/10000;
c1339151 229 Float_t cellY = (Float_t) celldumY/10;
230 Float_t cellX = (Float_t) celldumX/10;
231
562718f9 232 //
233 // Cell X centroid is back transformed
234 //
235 if (ismn < 12)
236 {
c1339151 237 celldataX[ihit] = (Int_t) ((cellX - (24-1) + cellY/2.) + 0.5);
562718f9 238 }
239 else if (ismn >= 12 && ismn <= 23)
240 {
c1339151 241 celldataX[ihit] = (Int_t) ((cellX - (48-1) + cellY/2.) + 0.5 );
562718f9 242 }
b0e4d1e1 243 celldataY[ihit] = (Int_t) (cellY + 0.5);
b6d6a9b5 244
245 Int_t irow = celldataX[ihit];
246 Int_t icol = celldataY[ihit];
247
248 if ((irow >= 0 && irow < 48) && (icol >= 0 && icol < 96))
249 {
250 celldataTr[ihit] = celltrack[irow][icol];
251 celldataPid[ihit] = cellpid[irow][icol];
252 celldataAdc[ihit] = (Float_t) celladc[irow][icol];
253 }
254 else
255 {
256 celldataTr[ihit] = -1;
257 celldataPid[ihit] = -1;
258 celldataAdc[ihit] = -1;
259 }
260
8c7250c5 261 }
262
920e13db 263 pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY,
b0e4d1e1 264 celldataTr, celldataPid, celldataAdc);
8c7250c5 265 pmdcont->Add(pmdcl);
266 }
2c1131dd 267 fPMDclucont->Delete();
8c7250c5 268}
269// ------------------------------------------------------------------------ //
562718f9 270Int_t AliPMDClusteringV2::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
271 Int_t iord1[], Double_t edepcell[])
8c7250c5 272{
273 // Does crude clustering
274 // Finds out only the big patch by just searching the
275 // connected cells
276 //
277
562718f9 278 Int_t i,j,k,id1,id2,icl, numcell;
279 Int_t jd1,jd2, icell, cellcount;
280 Int_t clust[2][5000];
281 static Int_t neibx[6] = {1,0,-1,-1,0,1}, neiby[6] = {0,1,1,0,-1,-1};
8c7250c5 282
283 // neibx and neiby define ( incremental ) (i,j) for the neighbours of a
284 // cell. There are six neighbours.
285 // cellcount --- total number of cells having nonzero ener dep
286 // numcell --- number of cells in a given supercluster
562718f9 287
8c7250c5 288 AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
289
562718f9 290 for (j=0; j < kNDIMX; j++)
291 {
292 for(k=0; k < kNDIMY; k++)
293 {
294 fInfocl[0][j][k] = 0;
295 fInfocl[1][j][k] = 0;
296 }
297 }
298
299 for(i=0; i < kNMX; i++)
300 {
301 fInfcl[0][i] = -1;
302
303 j = iord1[i];
304 id2 = j/kNDIMX;
305 id1 = j-id2*kNDIMX;
306
307 if(edepcell[j] <= cutoff)
308 {
309 fInfocl[0][id1][id2] = -1;
310 }
8c7250c5 311 }
8c7250c5 312 // ---------------------------------------------------------------
313 // crude clustering begins. Start with cell having largest adc
314 // count and loop over the cells in descending order of adc count
315 // ---------------------------------------------------------------
562718f9 316 icl = -1;
317 cellcount = -1;
318 for(icell=0; icell <= nmx1; icell++)
319 {
320 j = iord1[icell];
321 id2 = j/kNDIMX;
322 id1 = j-id2*kNDIMX;
323 if(fInfocl[0][id1][id2] == 0 )
324 {
325 // ---------------------------------------------------------------
326 // icl -- cluster #, numcell -- # of cells in it, clust -- stores
327 // coordinates of the cells in a cluster, fInfocl[0][i1][i2] is 1 for
328 // primary and 2 for secondary cells,
329 // fInfocl[1][i1][i2] stores cluster #
330 // ---------------------------------------------------------------
331 icl++;
332 numcell = 0;
333 cellcount++;
334 fInfocl[0][id1][id2] = 1;
335 fInfocl[1][id1][id2] = icl;
336 fInfcl[0][cellcount] = icl;
337 fInfcl[1][cellcount] = id1;
338 fInfcl[2][cellcount] = id2;
339
340 clust[0][numcell] = id1;
341 clust[1][numcell] = id2;
342 for(i = 1; i < 5000; i++)
343 {
344 clust[0][i] = -1;
345 }
346 // ---------------------------------------------------------------
347 // check for adc count in neib. cells. If ne 0 put it in this clust
348 // ---------------------------------------------------------------
349 for(i = 0; i < 6; i++)
350 {
351 jd1 = id1 + neibx[i];
352 jd2 = id2 + neiby[i];
8c7250c5 353 if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
562718f9 354 fInfocl[0][jd1][jd2] == 0)
355 {
356 numcell++;
357 fInfocl[0][jd1][jd2] = 2;
358 fInfocl[1][jd1][jd2] = icl;
359 clust[0][numcell] = jd1;
360 clust[1][numcell] = jd2;
361 cellcount++;
362 fInfcl[0][cellcount] = icl;
363 fInfcl[1][cellcount] = jd1;
364 fInfcl[2][cellcount] = jd2;
365 }
366 }
367 // ---------------------------------------------------------------
368 // check adc count for neighbour's neighbours recursively and
369 // if nonzero, add these to the cluster.
370 // ---------------------------------------------------------------
371 for(i = 1;i < 5000; i++)
372 {
373 if(clust[0][i] != -1)
374 {
375 id1 = clust[0][i];
376 id2 = clust[1][i];
377 for(j = 0; j < 6 ; j++)
378 {
379 jd1 = id1 + neibx[j];
380 jd2 = id2 + neiby[j];
381 if( (jd1 >= 0 && jd1 < kNDIMX) &&
382 (jd2 >= 0 && jd2 < kNDIMY)
383 && fInfocl[0][jd1][jd2] == 0 )
384 {
385 fInfocl[0][jd1][jd2] = 2;
386 fInfocl[1][jd1][jd2] = icl;
387 numcell++;
388 clust[0][numcell] = jd1;
389 clust[1][numcell] = jd2;
390 cellcount++;
391 fInfcl[0][cellcount] = icl;
392 fInfcl[1][cellcount] = jd1;
393 fInfcl[2][cellcount] = jd2;
394 }
395 }
396 }
8c7250c5 397 }
8c7250c5 398 }
8c7250c5 399 }
8c7250c5 400 return cellcount;
401}
402// ------------------------------------------------------------------------ //
562718f9 403 void AliPMDClusteringV2::RefClust(Int_t incr, Double_t edepcell[])
8c7250c5 404{
405 // Does the refining of clusters
406 // Takes the big patch and does gaussian fitting and
407 // finds out the more refined clusters
8c7250c5 408
c1339151 409 const Float_t ktwobysqrt3 = 1.1547;
2c1131dd 410 const Int_t kNmaxCell = 19;
411
562718f9 412 AliPMDcludata *pmdcludata = 0;
c1339151 413
78fc1b96 414 Int_t i12;
562718f9 415 Int_t i, j, k, i1, i2, id, icl, itest, ihld;
416 Int_t ig, nsupcl, clno, clX,clY;
2c1131dd 417 Int_t clxy[kNmaxCell];
c1339151 418
562718f9 419 Float_t clusdata[6];
420 Double_t x1, y1, z1, x2, y2, z2, rr;
c1339151 421
422 Int_t kndim = incr + 1;
423
424 TArrayI testncl;
425 TArrayI testindex;
426
427 Int_t *ncl, *iord;
428
429 Double_t *x, *y, *z, *xc, *yc, *zc, *cells, *rcl, *rcs;
430
431 ncl = new Int_t [kndim];
432 iord = new Int_t [kndim];
433 x = new Double_t [kndim];
434 y = new Double_t [kndim];
435 z = new Double_t [kndim];
436 xc = new Double_t [kndim];
437 yc = new Double_t [kndim];
438 zc = new Double_t [kndim];
439 cells = new Double_t [kndim];
440 rcl = new Double_t [kndim];
441 rcs = new Double_t [kndim];
562718f9 442
443 for(Int_t kk = 0; kk < 15; kk++)
444 {
445 if( kk < 6 )clusdata[kk] = 0.;
446 }
447
8c7250c5 448 // nsupcl = # of superclusters; ncl[i]= # of cells in supercluster i
449 // x, y and z store (x,y) coordinates of and energy deposited in a cell
450 // xc, yc store (x,y) coordinates of the cluster center
562718f9 451 // zc stores the energy deposited in a cluster, rc is cluster radius
8c7250c5 452
562718f9 453 clno = -1;
8c7250c5 454 nsupcl = -1;
c1339151 455
456 for(i = 0; i < kndim; i++)
562718f9 457 {
458 ncl[i] = -1;
8c7250c5 459 }
c1339151 460 for(i = 0; i <= incr; i++)
562718f9 461 {
462 if(fInfcl[0][i] != nsupcl)
463 {
464 nsupcl++;
465 }
466 if (nsupcl > 4500)
467 {
468 AliWarning("RefClust: Too many superclusters!");
469 nsupcl = 4500;
470 break;
471 }
472 ncl[nsupcl]++;
473 }
474
8c7250c5 475 AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
476 incr+1,nsupcl+1));
562718f9 477
478 id = -1;
479 icl = -1;
c1339151 480 for(i = 0; i <= nsupcl; i++)
562718f9 481 {
482 if(ncl[i] == 0)
8c7250c5 483 {
562718f9 484 id++;
485 icl++;
486 // one cell super-clusters --> single cluster
487 // cluster center at the centyer of the cell
488 // cluster radius = half cell dimension
489 if (clno >= 5000)
8c7250c5 490 {
562718f9 491 AliWarning("RefClust: Too many clusters! more than 5000");
492 return;
8c7250c5 493 }
562718f9 494 clno++;
495 i1 = fInfcl[1][id];
496 i2 = fInfcl[2][id];
78fc1b96 497 i12 = i1 + i2*kNDIMX;
562718f9 498 clusdata[0] = fCoord[0][i1][i2];
499 clusdata[1] = fCoord[1][i1][i2];
500 clusdata[2] = edepcell[i12];
501 clusdata[3] = 1.;
502 clusdata[4] = 0.0;
503 clusdata[5] = 0.0;
504
505 //cell information
c1339151 506
507 clY = (Int_t)((ktwobysqrt3*fCoord[1][i1][i2])*10);
508 clX = (Int_t)((fCoord[0][i1][i2] - clY/20.)*10);
509 clxy[0] = clX*10000 + clY ;
510
2c1131dd 511 for(Int_t icltr = 1; icltr < kNmaxCell; icltr++)
562718f9 512 {
513 clxy[icltr] = -1;
514 }
515 pmdcludata = new AliPMDcludata(clusdata,clxy);
516 fPMDclucont->Add(pmdcludata);
517
518
8c7250c5 519 }
562718f9 520 else if(ncl[i] == 1)
8c7250c5 521 {
562718f9 522 // two cell super-cluster --> single cluster
523 // cluster center is at ener. dep.-weighted mean of two cells
524 // cluster radius == half cell dimension
525 id++;
526 icl++;
527 if (clno >= 5000)
8c7250c5 528 {
529 AliWarning("RefClust: Too many clusters! more than 5000");
530 return;
531 }
562718f9 532 clno++;
533 i1 = fInfcl[1][id];
534 i2 = fInfcl[2][id];
78fc1b96 535 i12 = i1 + i2*kNDIMX;
562718f9 536
537 x1 = fCoord[0][i1][i2];
538 y1 = fCoord[1][i1][i2];
539 z1 = edepcell[i12];
540
541 id++;
542 i1 = fInfcl[1][id];
543 i2 = fInfcl[2][id];
544 i12 = i1 + i2*kNDIMX;
545
546 x2 = fCoord[0][i1][i2];
547 y2 = fCoord[1][i1][i2];
548 z2 = edepcell[i12];
549
550 clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
551 clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
552 clusdata[2] = z1+z2;
553 clusdata[3] = 2.;
554 clusdata[4] = (TMath::Sqrt(z1*z2))/(z1+z2);
555 clusdata[5] = 0.0;
556
c1339151 557 clY = (Int_t)((ktwobysqrt3*y1)*10);
558 clX = (Int_t)((x1 - clY/20.)*10);
559 clxy[0] = clX*10000 + clY ;
560
c1339151 561 clY = (Int_t)((ktwobysqrt3*y2)*10);
562 clX = (Int_t)((x2 - clY/20.)*10);
563 clxy[1] = clX*10000 + clY ;
564
2c1131dd 565 for(Int_t icltr = 2; icltr < kNmaxCell; icltr++)
8c7250c5 566 {
562718f9 567 clxy[icltr] = -1;
8c7250c5 568 }
562718f9 569 pmdcludata = new AliPMDcludata(clusdata, clxy);
570 fPMDclucont->Add(pmdcludata);
8c7250c5 571 }
562718f9 572 else{
573 id++;
574 iord[0] = 0;
575 // super-cluster of more than two cells - broken up into smaller
576 // clusters gaussian centers computed. (peaks separated by > 1 cell)
577 // Begin from cell having largest energy deposited This is first
578 // cluster center
579 // *****************************************************************
580 // NOTE --- POSSIBLE MODIFICATION: ONE MAY NOT BREAKING SUPERCLUSTERS
581 // IF NO. OF CELLS IS NOT TOO LARGE ( SAY 5 OR 6 )
582 // SINCE WE EXPECT THE SUPERCLUSTER
583 // TO BE A SINGLE CLUSTER
584 //*******************************************************************
585
586 i1 = fInfcl[1][id];
587 i2 = fInfcl[2][id];
78fc1b96 588 i12 = i1 + i2*kNDIMX;
562718f9 589
590 x[0] = fCoord[0][i1][i2];
591 y[0] = fCoord[1][i1][i2];
592 z[0] = edepcell[i12];
593
594 iord[0] = 0;
595 for(j = 1; j <= ncl[i]; j++)
596 {
597
598 id++;
599 i1 = fInfcl[1][id];
600 i2 = fInfcl[2][id];
78fc1b96 601 i12 = i1 + i2*kNDIMX;
562718f9 602 iord[j] = j;
603 x[j] = fCoord[0][i1][i2];
604 y[j] = fCoord[1][i1][i2];
605 z[j] = edepcell[i12];
606 }
607
608 // arranging cells within supercluster in decreasing order
609 for(j = 1; j <= ncl[i];j++)
610 {
611 itest = 0;
612 ihld = iord[j];
613 for(i1 = 0; i1 < j; i1++)
614 {
615 if(itest == 0 && z[iord[i1]] < z[ihld])
616 {
617 itest = 1;
618 for(i2 = j-1;i2 >= i1;i2--)
619 {
620 iord[i2+1] = iord[i2];
621 }
622 iord[i1] = ihld;
623 }
624 }
625 }
626
627
628 // compute the number of clusters and their centers ( first
629 // guess )
630 // centers must be separated by cells having smaller ener. dep.
631 // neighbouring centers should be either strong or well-separated
632 ig = 0;
633 xc[ig] = x[iord[0]];
634 yc[ig] = y[iord[0]];
635 zc[ig] = z[iord[0]];
636 for(j = 1; j <= ncl[i]; j++)
637 {
638 itest = -1;
639 x1 = x[iord[j]];
640 y1 = y[iord[j]];
641 for(k = 0; k <= ig; k++)
642 {
643 x2 = xc[k];
644 y2 = yc[k];
645 rr = Distance(x1,y1,x2,y2);
646 //************************************************************
647 // finetuning cluster splitting
648 // the numbers zc/4 and zc/10 may need to be changed.
649 // Also one may need to add one more layer because our
650 // cells are smaller in absolute scale
651 //************************************************************
652
653
654 if( rr >= 1.1 && rr < 1.8 && z[iord[j]] > zc[k]/4.) itest++;
655 if( rr >= 1.8 && rr < 2.1 && z[iord[j]] > zc[k]/10.) itest++;
656 if( rr >= 2.1)itest++;
657 }
658
659 if(itest == ig)
660 {
661 ig++;
662 xc[ig] = x1;
663 yc[ig] = y1;
664 zc[ig] = z[iord[j]];
665 }
666 }
c1339151 667 ClustDetails(ncl[i], ig, x, y ,z, xc, yc, zc, rcl, rcs, cells,
668 testncl, testindex);
562718f9 669
670 Int_t pp = 0;
671 for(j = 0; j <= ig; j++)
672 {
673 clno++;
674 if (clno >= 5000)
675 {
676 AliWarning("RefClust: Too many clusters! more than 5000");
677 return;
678 }
679 clusdata[0] = xc[j];
680 clusdata[1] = yc[j];
681 clusdata[2] = zc[j];
682 clusdata[4] = rcl[j];
683 clusdata[5] = rcs[j];
684 if(ig == 0)
685 {
2c1131dd 686 clusdata[3] = ncl[i] + 1;
562718f9 687 }
688 else
689 {
690 clusdata[3] = cells[j];
691 }
692 // cell information
693 Int_t ncellcls = testncl[j];
2c1131dd 694 if( ncellcls < kNmaxCell )
562718f9 695 {
696 for(Int_t kk = 1; kk <= ncellcls; kk++)
697 {
698 Int_t ll = testindex[pp];
c1339151 699 clY = (Int_t)((ktwobysqrt3*y[ll])*10);
700 clX = (Int_t)((x[ll] - clY/20.)*10);
701 clxy[kk-1] = clX*10000 + clY ;
c1339151 702
562718f9 703 pp++;
704 }
2c1131dd 705 for(Int_t icltr = ncellcls ; icltr < kNmaxCell; icltr++)
562718f9 706 {
707 clxy[icltr] = -1;
708 }
709 }
710 pmdcludata = new AliPMDcludata(clusdata, clxy);
711 fPMDclucont->Add(pmdcludata);
712 }
713 testncl.Set(0);
714 testindex.Set(0);
715 }
8c7250c5 716 }
c1339151 717 delete [] ncl;
718 delete [] iord;
719 delete [] x;
720 delete [] y;
721 delete [] z;
722 delete [] xc;
723 delete [] yc;
724 delete [] zc;
725 delete [] cells;
726 delete [] rcl;
727 delete [] rcs;
8c7250c5 728}
8c7250c5 729// ------------------------------------------------------------------------ //
c1339151 730void AliPMDClusteringV2::ClustDetails(Int_t ncell, Int_t nclust, Double_t x[],
731 Double_t y[], Double_t z[],Double_t xc[],
732 Double_t yc[], Double_t zc[],
733 Double_t rcl[], Double_t rcs[],
734 Double_t cells[], TArrayI &testncl,
562718f9 735 TArrayI &testindex)
8c7250c5 736{
737 // function begins
738 //
8c7250c5 739
c1339151 740 Int_t kndim1 = ncell + 1;//ncell
741 Int_t kndim2 = 20;
742 Int_t kndim3 = nclust + 1;//nclust
78fc1b96 743
562718f9 744 Int_t i, j, k, i1, i2;
562718f9 745 Double_t x1, y1, x2, y2, rr, b, c, r1, r2;
746 Double_t sumx, sumy, sumxy, sumxx, sum, sum1, sumyy;
8c7250c5 747
c1339151 748 Double_t *str, *str1, *xcl, *ycl, *cln;
749 Int_t **cell;
750 Int_t ** cluster;
751 Double_t **clustcell;
752 str = new Double_t [kndim3];
753 str1 = new Double_t [kndim3];
754 xcl = new Double_t [kndim3];
755 ycl = new Double_t [kndim3];
756 cln = new Double_t [kndim3];
562718f9 757
c1339151 758 clustcell = new Double_t *[kndim3];
759 cell = new Int_t *[kndim3];
760 cluster = new Int_t *[kndim1];
562718f9 761 for(i = 0; i < kndim1; i++)
762 {
c1339151 763 cluster[i] = new Int_t [kndim2];
764 }
765
766 for(i = 0; i < kndim3; i++)
767 {
768 str[i] = 0;
769 str1[i] = 0;
770 xcl[i] = 0;
771 ycl[i] = 0;
772 cln[i] = 0;
773
774 cell[i] = new Int_t [kndim2];
775 clustcell[i] = new Double_t [kndim1];
776 for(j = 0; j < kndim1; j++)
8c7250c5 777 {
c1339151 778 clustcell[i][j] = 0;
562718f9 779 }
c1339151 780 for(j = 0; j < kndim2; j++)
562718f9 781 {
c1339151 782 cluster[i][j] = 0;
783 cell[i][j] = 0;
8c7250c5 784 }
562718f9 785 }
c1339151 786
562718f9 787 if(nclust > 0)
788 {
789 // more than one cluster
790 // checking cells shared between several clusters.
791 // First check if the cell is within
792 // one cell unit ( nearest neighbour). Else,
793 // if it is within 1.74 cell units ( next nearest )
794 // Else if it is upto 2 cell units etc.
795
796 for (i = 0; i <= ncell; i++)
8c7250c5 797 {
c1339151 798 x1 = x[i];
799 y1 = y[i];
562718f9 800 cluster[i][0] = 0;
2c1131dd 801
562718f9 802 // distance <= 1 cell unit
2c1131dd 803
562718f9 804 for(j = 0; j <= nclust; j++)
8c7250c5 805 {
c1339151 806 x2 = xc[j];
807 y2 = yc[j];
8c7250c5 808 rr = Distance(x1, y1, x2, y2);
562718f9 809 if(rr <= 1.)
8c7250c5 810 {
811 cluster[i][0]++;
812 i1 = cluster[i][0];
813 cluster[i][i1] = j;
814 }
815 }
562718f9 816 // next nearest neighbour
817 if(cluster[i][0] == 0)
818 {
819 for(j=0; j<=nclust; j++)
820 {
c1339151 821 x2 = xc[j];
822 y2 = yc[j];
562718f9 823 rr = Distance(x1, y1, x2, y2);
824 if(rr <= TMath::Sqrt(3.))
825 {
826 cluster[i][0]++;
827 i1 = cluster[i][0];
828 cluster[i][i1] = j;
829 }
830 }
831 }
832 // next-to-next nearest neighbour
833 if(cluster[i][0] == 0)
834 {
835 for(j=0; j<=nclust; j++)
836 {
c1339151 837 x2 = xc[j];
838 y2 = yc[j];
562718f9 839 rr = Distance(x1, y1, x2, y2);
840 if(rr <= 2.)
841 {
842 cluster[i][0]++;
843 i1 = cluster[i][0];
844 cluster[i][i1] = j;
845 }
846 }
847 }
848 // one more
849 if(cluster[i][0] == 0)
850 {
851 for(j = 0; j <= nclust; j++)
852 {
c1339151 853 x2 = xc[j];
854 y2 = yc[j];
562718f9 855 rr = Distance(x1, y1, x2, y2);
856 if(rr <= 2.7)
857 {
858 cluster[i][0]++;
859 i1 = cluster[i][0];
860 cluster[i][i1] = j;
861 }
862 }
863 }
8c7250c5 864 }
562718f9 865
866 // computing cluster strength. Some cells are shared.
867 for(i = 0; i <= ncell; i++)
8c7250c5 868 {
562718f9 869 if(cluster[i][0] != 0)
8c7250c5 870 {
562718f9 871 i1 = cluster[i][0];
872 for(j = 1; j <= i1; j++)
8c7250c5 873 {
562718f9 874 i2 = cluster[i][j];
c1339151 875 str[i2] += z[i]/i1;
8c7250c5 876 }
877 }
878 }
562718f9 879
880 for(k = 0; k < 5; k++)
8c7250c5 881 {
562718f9 882 for(i = 0; i <= ncell; i++)
8c7250c5 883 {
562718f9 884 if(cluster[i][0] != 0)
8c7250c5 885 {
562718f9 886 i1=cluster[i][0];
887 sum=0.;
888 for(j = 1; j <= i1; j++)
889 {
890 sum += str[cluster[i][j]];
891 }
892
893 for(j = 1; j <= i1; j++)
894 {
895 i2 = cluster[i][j];
c1339151 896 str1[i2] += z[i]*str[i2]/sum;
897 clustcell[i2][i] = z[i]*str[i2]/sum;
562718f9 898 }
8c7250c5 899 }
900 }
562718f9 901
902
903 for(j = 0; j <= nclust; j++)
904 {
905 str[j] = str1[j];
906 str1[j] = 0.;
907 }
8c7250c5 908 }
562718f9 909
910 for(i = 0; i <= nclust; i++)
911 {
912 sumx = 0.;
913 sumy = 0.;
914 sum = 0.;
915 sum1 = 0.;
916 for(j = 0; j <= ncell; j++)
917 {
918 if(clustcell[i][j] != 0)
919 {
c1339151 920 sumx += clustcell[i][j]*x[j];
921 sumy += clustcell[i][j]*y[j];
562718f9 922 sum += clustcell[i][j];
c1339151 923 sum1 += clustcell[i][j]/z[j];
562718f9 924 }
925 }
926 //** xcl and ycl are cluster centroid positions ( center of gravity )
927
928 xcl[i] = sumx/sum;
929 ycl[i] = sumy/sum;
930 cln[i] = sum1;
931 sumxx = 0.;
932 sumyy = 0.;
933 sumxy = 0.;
934 for(j = 0; j <= ncell; j++)
935 {
c1339151 936 sumxx += clustcell[i][j]*(x[j]-xcl[i])*(x[j]-xcl[i])/sum;
937 sumyy += clustcell[i][j]*(y[j]-ycl[i])*(y[j]-ycl[i])/sum;
938 sumxy += clustcell[i][j]*(x[j]-xcl[i])*(y[j]-ycl[i])/sum;
562718f9 939 }
940 b = sumxx+sumyy;
941 c = sumxx*sumyy-sumxy*sumxy;
942 // ******************r1 and r2 are major and minor axes ( r1 > r2 ).
943 r1 = b/2.+TMath::Sqrt(b*b/4.-c);
944 r2 = b/2.-TMath::Sqrt(b*b/4.-c);
945 // final assignments to proper external variables
c1339151 946 xc[i] = xcl[i];
947 yc[i] = ycl[i];
948 zc[i] = str[i];
949 cells[i] = cln[i];
950 rcl[i] = r1;
951 rcs[i] = r2;
952
8c7250c5 953 }
562718f9 954
955 //To get the cell position in a cluster
956
957 for(Int_t ii=0; ii<= ncell; ii++)
958 {
959 Int_t jj = cluster[ii][0];
960 for(Int_t kk=1; kk<= jj; kk++)
8c7250c5 961 {
562718f9 962 Int_t ll = cluster[ii][kk];
963 cell[ll][0]++;
964 cell[ll][cell[ll][0]] = ii;
8c7250c5 965 }
562718f9 966 }
967
968 testncl.Set(nclust+1);
969 Int_t counter = 0;
970
971 for(Int_t ii=0; ii <= nclust; ii++)
972 {
973 testncl[ii] = cell[ii][0];
974 counter += testncl[ii];
975 }
976 testindex.Set(counter);
977 Int_t ll = 0;
978 for(Int_t ii=0; ii<= nclust; ii++)
979 {
980 for(Int_t jj = 1; jj<= testncl[ii]; jj++)
981 {
982 Int_t kk = cell[ii][jj];
983 testindex[ll] = kk;
984 ll++;
985 }
986 }
987
988 }
c1339151 989 else if(nclust == 0)
562718f9 990 {
8c7250c5 991 sumx = 0.;
992 sumy = 0.;
993 sum = 0.;
994 sum1 = 0.;
562718f9 995 i = 0;
996 for(j = 0; j <= ncell; j++)
997 {
c1339151 998 sumx += z[j]*x[j];
999 sumy += z[j]*y[j];
1000 sum += z[j];
562718f9 1001 sum1++;
8c7250c5 1002 }
8c7250c5 1003 xcl[i] = sumx/sum;
1004 ycl[i] = sumy/sum;
1005 cln[i] = sum1;
562718f9 1006 sumxx = 0.;
1007 sumyy = 0.;
1008 sumxy = 0.;
1009 for(j = 0; j <= ncell; j++)
1010 {
c1339151 1011 sumxx += clustcell[i][j]*(x[j]-xcl[i])*(x[j]-xcl[i])/sum;
1012 sumyy += clustcell[i][j]*(y[j]-ycl[i])*(y[j]-ycl[i])/sum;
1013 sumxy += clustcell[i][j]*(x[j]-xcl[i])*(y[j]-ycl[i])/sum;
562718f9 1014 }
1015 b = sumxx+sumyy;
1016 c = sumxx*sumyy-sumxy*sumxy;
1017 r1 = b/2.+ TMath::Sqrt(b*b/4.-c);
1018 r2 = b/2.- TMath::Sqrt(b*b/4.-c);
1019
1020 // To get the cell position in a cluster
1021 testncl.Set(nclust+1);
c1339151 1022 testindex.Set(ncell+1);
2c1131dd 1023 cell[0][0] = ncell + 1;
562718f9 1024 testncl[0] = cell[0][0];
1025 Int_t ll = 0;
c1339151 1026 for(Int_t ii = 1; ii <= ncell; ii++)
562718f9 1027 {
1028 cell[0][ii]=ii;
562718f9 1029 Int_t kk = cell[0][ii];
1030 testindex[ll] = kk;
1031 ll++;
1032 }
1033 // final assignments
c1339151 1034 xc[i] = xcl[i];
1035 yc[i] = ycl[i];
c1339151 1036 zc[i] = sum;
1037 cells[i] = cln[i];
1038 rcl[i] = r1;
1039 rcs[i] = r2;
1040 }
1041 for(i = 0; i < kndim3; i++)
1042 {
1043 delete [] clustcell[i];
1044 delete [] cell[i];
1045 }
1046 delete [] clustcell;
1047 delete [] cell;
1048 for(i = 0; i <kndim1 ; i++)
1049 {
1050 delete [] cluster[i];
8c7250c5 1051 }
c1339151 1052 delete [] cluster;
1053 delete [] str;
1054 delete [] str1;
1055 delete [] xcl;
1056 delete [] ycl;
1057 delete [] cln;
8c7250c5 1058}
1059
1060// ------------------------------------------------------------------------ //
1061Double_t AliPMDClusteringV2::Distance(Double_t x1, Double_t y1,
1062 Double_t x2, Double_t y2)
1063{
562718f9 1064 return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
8c7250c5 1065}
1066// ------------------------------------------------------------------------ //
1067void AliPMDClusteringV2::SetEdepCut(Float_t decut)
1068{
1069 fCutoff = decut;
1070}
1071// ------------------------------------------------------------------------ //
4755c3f1 1072void AliPMDClusteringV2::SetClusteringParam(Int_t cluspar)
1073{
1074 fClusParam = cluspar;
1075}
1076// ------------------------------------------------------------------------ //