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