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