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