PMD DA informations : satya
[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
be8b7039 108 Int_t i = 0, j = 0, nmx1 = 0;
109 Int_t incr = 0, id = 0, jd = 0;
562718f9 110 Int_t ndimXr = 0;
111 Int_t ndimYr = 0;
2c1131dd 112 Int_t celldataX[kNmaxCell], celldataY[kNmaxCell];
920e13db 113 Int_t celldataTr[kNmaxCell], celldataPid[kNmaxCell];
b0e4d1e1 114 Float_t celldataAdc[kNmaxCell];
be8b7039 115 Float_t clusdata[6] = {0.,0.,0.,0.,0.,0.};
116 Double_t cutoff = 0., ave = 0.;
562718f9 117 Double_t edepcell[kNMX];
8c7250c5 118
8c7250c5 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 }
562718f9 130
78fc1b96 131 for (i =0; i < kNMX; i++)
8c7250c5 132 {
562718f9 133 edepcell[i] = 0.;
8c7250c5 134 }
562718f9 135
8c7250c5 136 for (id = 0; id < ndimXr; id++)
137 {
138 for (jd = 0; jd < ndimYr; jd++)
139 {
562718f9 140 j = jd;
141 i = id + (ndimYr/2-1) - (jd/2);
142 Int_t ij = i + j*kNDIMX;
8c7250c5 143 if (ismn < 12)
144 {
562718f9 145 edepcell[ij] = celladc[jd][id];
8c7250c5 146 }
147 else if (ismn >= 12 && ismn <= 23)
148 {
562718f9 149 edepcell[ij] = celladc[id][jd];
8c7250c5 150 }
151
152 }
153 }
154
5406008b 155 // the dimension of iord1 is increased twice
156 Int_t iord1[2*kNMX];
df4e6759 157 TMath::Sort((Int_t)kNMX,edepcell,iord1);// order the data
8c7250c5 158 cutoff = fCutoff; // cutoff used to discard cells having ener. dep.
562718f9 159 ave = 0.;
160 nmx1 = -1;
8c7250c5 161
562718f9 162 for(i = 0;i < kNMX; i++)
8c7250c5 163 {
562718f9 164 if(edepcell[i] > 0.)
165 {
166 ave += edepcell[i];
167 }
168 if(edepcell[i] > cutoff )
169 {
170 nmx1++;
171 }
8c7250c5 172 }
562718f9 173
8c7250c5 174 AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
562718f9 175
176 if (nmx1 == 0)
177 {
178 nmx1 = 1;
179 }
180 ave = ave/nmx1;
181
8c7250c5 182 AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
183 kNMX,ave));
8c7250c5 184
562718f9 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++)
8c7250c5 192 {
562718f9 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
8c7250c5 202 Float_t cluY0 = ktwobysqrt3*cluYC;
203 Float_t cluX0 = cluXC - cluY0/2.;
562718f9 204
8c7250c5 205 //
206 // Cluster X centroid is back transformed
207 //
208 if (ismn < 12)
209 {
210 clusdata[0] = cluX0 - (24-1) + cluY0/2.;
211 }
562718f9 212 else if (ismn >= 12 && ismn <= 23)
8c7250c5 213 {
214 clusdata[0] = cluX0 - (48-1) + cluY0/2.;
215 }
216
562718f9 217 clusdata[1] = cluY0;
218 clusdata[2] = cluADC;
219 clusdata[3] = cluCELLS;
220 clusdata[4] = cluSIGX;
221 clusdata[5] = cluSIGY;
8c7250c5 222 //
223 // Cells associated with a cluster
224 //
2c1131dd 225 for (Int_t ihit = 0; ihit < kNmaxCell; ihit++)
8c7250c5 226 {
562718f9 227 Int_t dummyXY = pmdcludata->GetCellXY(ihit);
228
229 Int_t celldumY = dummyXY%10000;
230 Int_t celldumX = dummyXY/10000;
c1339151 231 Float_t cellY = (Float_t) celldumY/10;
232 Float_t cellX = (Float_t) celldumX/10;
233
562718f9 234 //
235 // Cell X centroid is back transformed
236 //
237 if (ismn < 12)
238 {
c1339151 239 celldataX[ihit] = (Int_t) ((cellX - (24-1) + cellY/2.) + 0.5);
562718f9 240 }
241 else if (ismn >= 12 && ismn <= 23)
242 {
c1339151 243 celldataX[ihit] = (Int_t) ((cellX - (48-1) + cellY/2.) + 0.5 );
562718f9 244 }
b0e4d1e1 245 celldataY[ihit] = (Int_t) (cellY + 0.5);
b6d6a9b5 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
8c7250c5 263 }
264
920e13db 265 pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY,
b0e4d1e1 266 celldataTr, celldataPid, celldataAdc);
8c7250c5 267 pmdcont->Add(pmdcl);
268 }
2c1131dd 269 fPMDclucont->Delete();
8c7250c5 270}
271// ------------------------------------------------------------------------ //
562718f9 272Int_t AliPMDClusteringV2::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
273 Int_t iord1[], Double_t edepcell[])
8c7250c5 274{
275 // Does crude clustering
276 // Finds out only the big patch by just searching the
277 // connected cells
278 //
279
be8b7039 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;
562718f9 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};
8c7250c5 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
562718f9 289
8c7250c5 290 AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
291
562718f9 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 }
8c7250c5 313 }
8c7250c5 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 // ---------------------------------------------------------------
562718f9 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];
8c7250c5 355 if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
562718f9 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 }
8c7250c5 399 }
8c7250c5 400 }
8c7250c5 401 }
8c7250c5 402 return cellcount;
403}
404// ------------------------------------------------------------------------ //
562718f9 405 void AliPMDClusteringV2::RefClust(Int_t incr, Double_t edepcell[])
8c7250c5 406{
407 // Does the refining of clusters
408 // Takes the big patch and does gaussian fitting and
409 // finds out the more refined clusters
8c7250c5 410
c1339151 411 const Float_t ktwobysqrt3 = 1.1547;
2c1131dd 412 const Int_t kNmaxCell = 19;
413
562718f9 414 AliPMDcludata *pmdcludata = 0;
c1339151 415
be8b7039 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];
c1339151 421
be8b7039 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.;
c1339151 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];
562718f9 445
446 for(Int_t kk = 0; kk < 15; kk++)
447 {
448 if( kk < 6 )clusdata[kk] = 0.;
449 }
450
8c7250c5 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
562718f9 454 // zc stores the energy deposited in a cluster, rc is cluster radius
8c7250c5 455
562718f9 456 clno = -1;
8c7250c5 457 nsupcl = -1;
c1339151 458
459 for(i = 0; i < kndim; i++)
562718f9 460 {
461 ncl[i] = -1;
8c7250c5 462 }
c1339151 463 for(i = 0; i <= incr; i++)
562718f9 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
8c7250c5 478 AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
479 incr+1,nsupcl+1));
562718f9 480
481 id = -1;
482 icl = -1;
c1339151 483 for(i = 0; i <= nsupcl; i++)
562718f9 484 {
485 if(ncl[i] == 0)
8c7250c5 486 {
562718f9 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)
8c7250c5 493 {
562718f9 494 AliWarning("RefClust: Too many clusters! more than 5000");
495 return;
8c7250c5 496 }
562718f9 497 clno++;
498 i1 = fInfcl[1][id];
499 i2 = fInfcl[2][id];
78fc1b96 500 i12 = i1 + i2*kNDIMX;
562718f9 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
c1339151 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
2c1131dd 514 for(Int_t icltr = 1; icltr < kNmaxCell; icltr++)
562718f9 515 {
516 clxy[icltr] = -1;
517 }
518 pmdcludata = new AliPMDcludata(clusdata,clxy);
519 fPMDclucont->Add(pmdcludata);
520
521
8c7250c5 522 }
562718f9 523 else if(ncl[i] == 1)
8c7250c5 524 {
562718f9 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)
8c7250c5 531 {
532 AliWarning("RefClust: Too many clusters! more than 5000");
533 return;
534 }
562718f9 535 clno++;
536 i1 = fInfcl[1][id];
537 i2 = fInfcl[2][id];
78fc1b96 538 i12 = i1 + i2*kNDIMX;
562718f9 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
c1339151 560 clY = (Int_t)((ktwobysqrt3*y1)*10);
561 clX = (Int_t)((x1 - clY/20.)*10);
562 clxy[0] = clX*10000 + clY ;
563
c1339151 564 clY = (Int_t)((ktwobysqrt3*y2)*10);
565 clX = (Int_t)((x2 - clY/20.)*10);
566 clxy[1] = clX*10000 + clY ;
567
2c1131dd 568 for(Int_t icltr = 2; icltr < kNmaxCell; icltr++)
8c7250c5 569 {
562718f9 570 clxy[icltr] = -1;
8c7250c5 571 }
562718f9 572 pmdcludata = new AliPMDcludata(clusdata, clxy);
573 fPMDclucont->Add(pmdcludata);
8c7250c5 574 }
562718f9 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];
78fc1b96 591 i12 = i1 + i2*kNDIMX;
562718f9 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];
78fc1b96 604 i12 = i1 + i2*kNDIMX;
562718f9 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 }
c1339151 670 ClustDetails(ncl[i], ig, x, y ,z, xc, yc, zc, rcl, rcs, cells,
671 testncl, testindex);
562718f9 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 {
2c1131dd 689 clusdata[3] = ncl[i] + 1;
562718f9 690 }
691 else
692 {
693 clusdata[3] = cells[j];
694 }
695 // cell information
696 Int_t ncellcls = testncl[j];
2c1131dd 697 if( ncellcls < kNmaxCell )
562718f9 698 {
699 for(Int_t kk = 1; kk <= ncellcls; kk++)
700 {
701 Int_t ll = testindex[pp];
c1339151 702 clY = (Int_t)((ktwobysqrt3*y[ll])*10);
703 clX = (Int_t)((x[ll] - clY/20.)*10);
704 clxy[kk-1] = clX*10000 + clY ;
c1339151 705
562718f9 706 pp++;
707 }
2c1131dd 708 for(Int_t icltr = ncellcls ; icltr < kNmaxCell; icltr++)
562718f9 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 }
8c7250c5 719 }
c1339151 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;
8c7250c5 731}
8c7250c5 732// ------------------------------------------------------------------------ //
c1339151 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,
562718f9 738 TArrayI &testindex)
8c7250c5 739{
740 // function begins
741 //
8c7250c5 742
c1339151 743 Int_t kndim1 = ncell + 1;//ncell
744 Int_t kndim2 = 20;
745 Int_t kndim3 = nclust + 1;//nclust
78fc1b96 746
be8b7039 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.;
8c7250c5 752
c1339151 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];
562718f9 762
c1339151 763 clustcell = new Double_t *[kndim3];
764 cell = new Int_t *[kndim3];
765 cluster = new Int_t *[kndim1];
562718f9 766 for(i = 0; i < kndim1; i++)
767 {
c1339151 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++)
8c7250c5 782 {
c1339151 783 clustcell[i][j] = 0;
562718f9 784 }
c1339151 785 for(j = 0; j < kndim2; j++)
562718f9 786 {
c1339151 787 cluster[i][j] = 0;
788 cell[i][j] = 0;
8c7250c5 789 }
562718f9 790 }
c1339151 791
562718f9 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++)
8c7250c5 802 {
c1339151 803 x1 = x[i];
804 y1 = y[i];
562718f9 805 cluster[i][0] = 0;
2c1131dd 806
562718f9 807 // distance <= 1 cell unit
2c1131dd 808
562718f9 809 for(j = 0; j <= nclust; j++)
8c7250c5 810 {
c1339151 811 x2 = xc[j];
812 y2 = yc[j];
8c7250c5 813 rr = Distance(x1, y1, x2, y2);
562718f9 814 if(rr <= 1.)
8c7250c5 815 {
816 cluster[i][0]++;
817 i1 = cluster[i][0];
818 cluster[i][i1] = j;
819 }
820 }
562718f9 821 // next nearest neighbour
822 if(cluster[i][0] == 0)
823 {
824 for(j=0; j<=nclust; j++)
825 {
c1339151 826 x2 = xc[j];
827 y2 = yc[j];
562718f9 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 {
c1339151 842 x2 = xc[j];
843 y2 = yc[j];
562718f9 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 {
c1339151 858 x2 = xc[j];
859 y2 = yc[j];
562718f9 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 }
8c7250c5 869 }
562718f9 870
871 // computing cluster strength. Some cells are shared.
872 for(i = 0; i <= ncell; i++)
8c7250c5 873 {
562718f9 874 if(cluster[i][0] != 0)
8c7250c5 875 {
562718f9 876 i1 = cluster[i][0];
877 for(j = 1; j <= i1; j++)
8c7250c5 878 {
562718f9 879 i2 = cluster[i][j];
c1339151 880 str[i2] += z[i]/i1;
8c7250c5 881 }
882 }
883 }
562718f9 884
885 for(k = 0; k < 5; k++)
8c7250c5 886 {
562718f9 887 for(i = 0; i <= ncell; i++)
8c7250c5 888 {
562718f9 889 if(cluster[i][0] != 0)
8c7250c5 890 {
562718f9 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];
c1339151 901 str1[i2] += z[i]*str[i2]/sum;
902 clustcell[i2][i] = z[i]*str[i2]/sum;
562718f9 903 }
8c7250c5 904 }
905 }
562718f9 906
907
908 for(j = 0; j <= nclust; j++)
909 {
910 str[j] = str1[j];
911 str1[j] = 0.;
912 }
8c7250c5 913 }
562718f9 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 {
c1339151 925 sumx += clustcell[i][j]*x[j];
926 sumy += clustcell[i][j]*y[j];
562718f9 927 sum += clustcell[i][j];
c1339151 928 sum1 += clustcell[i][j]/z[j];
562718f9 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 {
c1339151 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;
562718f9 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
c1339151 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
8c7250c5 958 }
562718f9 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++)
8c7250c5 966 {
562718f9 967 Int_t ll = cluster[ii][kk];
968 cell[ll][0]++;
969 cell[ll][cell[ll][0]] = ii;
8c7250c5 970 }
562718f9 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 }
c1339151 994 else if(nclust == 0)
562718f9 995 {
8c7250c5 996 sumx = 0.;
997 sumy = 0.;
998 sum = 0.;
999 sum1 = 0.;
562718f9 1000 i = 0;
1001 for(j = 0; j <= ncell; j++)
1002 {
c1339151 1003 sumx += z[j]*x[j];
1004 sumy += z[j]*y[j];
1005 sum += z[j];
562718f9 1006 sum1++;
8c7250c5 1007 }
8c7250c5 1008 xcl[i] = sumx/sum;
1009 ycl[i] = sumy/sum;
1010 cln[i] = sum1;
562718f9 1011 sumxx = 0.;
1012 sumyy = 0.;
1013 sumxy = 0.;
1014 for(j = 0; j <= ncell; j++)
1015 {
c1339151 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;
562718f9 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);
c1339151 1027 testindex.Set(ncell+1);
2c1131dd 1028 cell[0][0] = ncell + 1;
562718f9 1029 testncl[0] = cell[0][0];
1030 Int_t ll = 0;
c1339151 1031 for(Int_t ii = 1; ii <= ncell; ii++)
562718f9 1032 {
1033 cell[0][ii]=ii;
562718f9 1034 Int_t kk = cell[0][ii];
1035 testindex[ll] = kk;
1036 ll++;
1037 }
1038 // final assignments
c1339151 1039 xc[i] = xcl[i];
1040 yc[i] = ycl[i];
c1339151 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];
8c7250c5 1056 }
c1339151 1057 delete [] cluster;
1058 delete [] str;
1059 delete [] str1;
1060 delete [] xcl;
1061 delete [] ycl;
1062 delete [] cln;
8c7250c5 1063}
1064
1065// ------------------------------------------------------------------------ //
1066Double_t AliPMDClusteringV2::Distance(Double_t x1, Double_t y1,
1067 Double_t x2, Double_t y2)
1068{
562718f9 1069 return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
8c7250c5 1070}
1071// ------------------------------------------------------------------------ //
1072void AliPMDClusteringV2::SetEdepCut(Float_t decut)
1073{
1074 fCutoff = decut;
1075}
1076// ------------------------------------------------------------------------ //
4755c3f1 1077void AliPMDClusteringV2::SetClusteringParam(Int_t cluspar)
1078{
1079 fClusParam = cluspar;
1080}
1081// ------------------------------------------------------------------------ //