]> 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
QA update by Sylwester
[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
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 //
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
78fc1b96 391 Int_t i12;
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];
78fc1b96 474 i12 = i1 + i2*kNDIMX;
562718f9 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
484 clY = (Int_t)((ktwobysqrt3*fCoord[1][i1][i2])*10);
485 clX = (Int_t)((fCoord[0][i1][i2] - clY/20.)*10);
486 clxy[0] = clX*10000 + clY ;
487
562718f9 488 for(Int_t icltr = 1; icltr < 15; icltr++)
489 {
490 clxy[icltr] = -1;
491 }
492 pmdcludata = new AliPMDcludata(clusdata,clxy);
493 fPMDclucont->Add(pmdcludata);
494
495
8c7250c5 496 }
562718f9 497 else if(ncl[i] == 1)
8c7250c5 498 {
562718f9 499 // two cell super-cluster --> single cluster
500 // cluster center is at ener. dep.-weighted mean of two cells
501 // cluster radius == half cell dimension
502 id++;
503 icl++;
504 if (clno >= 5000)
8c7250c5 505 {
506 AliWarning("RefClust: Too many clusters! more than 5000");
507 return;
508 }
562718f9 509 clno++;
510 i1 = fInfcl[1][id];
511 i2 = fInfcl[2][id];
78fc1b96 512 i12 = i1 + i2*kNDIMX;
562718f9 513
514 x1 = fCoord[0][i1][i2];
515 y1 = fCoord[1][i1][i2];
516 z1 = edepcell[i12];
517
518 id++;
519 i1 = fInfcl[1][id];
520 i2 = fInfcl[2][id];
521 i12 = i1 + i2*kNDIMX;
522
523 x2 = fCoord[0][i1][i2];
524 y2 = fCoord[1][i1][i2];
525 z2 = edepcell[i12];
526
527 clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
528 clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
529 clusdata[2] = z1+z2;
530 clusdata[3] = 2.;
531 clusdata[4] = (TMath::Sqrt(z1*z2))/(z1+z2);
532 clusdata[5] = 0.0;
533
c1339151 534 clY = (Int_t)((ktwobysqrt3*y1)*10);
535 clX = (Int_t)((x1 - clY/20.)*10);
536 clxy[0] = clX*10000 + clY ;
537
c1339151 538 clY = (Int_t)((ktwobysqrt3*y2)*10);
539 clX = (Int_t)((x2 - clY/20.)*10);
540 clxy[1] = clX*10000 + clY ;
541
562718f9 542 for(Int_t icltr = 2; icltr < 15; icltr++)
8c7250c5 543 {
562718f9 544 clxy[icltr] = -1;
8c7250c5 545 }
562718f9 546 pmdcludata = new AliPMDcludata(clusdata, clxy);
547 fPMDclucont->Add(pmdcludata);
8c7250c5 548 }
562718f9 549 else{
550 id++;
551 iord[0] = 0;
552 // super-cluster of more than two cells - broken up into smaller
553 // clusters gaussian centers computed. (peaks separated by > 1 cell)
554 // Begin from cell having largest energy deposited This is first
555 // cluster center
556 // *****************************************************************
557 // NOTE --- POSSIBLE MODIFICATION: ONE MAY NOT BREAKING SUPERCLUSTERS
558 // IF NO. OF CELLS IS NOT TOO LARGE ( SAY 5 OR 6 )
559 // SINCE WE EXPECT THE SUPERCLUSTER
560 // TO BE A SINGLE CLUSTER
561 //*******************************************************************
562
563 i1 = fInfcl[1][id];
564 i2 = fInfcl[2][id];
78fc1b96 565 i12 = i1 + i2*kNDIMX;
562718f9 566
567 x[0] = fCoord[0][i1][i2];
568 y[0] = fCoord[1][i1][i2];
569 z[0] = edepcell[i12];
570
571 iord[0] = 0;
572 for(j = 1; j <= ncl[i]; j++)
573 {
574
575 id++;
576 i1 = fInfcl[1][id];
577 i2 = fInfcl[2][id];
78fc1b96 578 i12 = i1 + i2*kNDIMX;
562718f9 579 iord[j] = j;
580 x[j] = fCoord[0][i1][i2];
581 y[j] = fCoord[1][i1][i2];
582 z[j] = edepcell[i12];
583 }
584
585 // arranging cells within supercluster in decreasing order
586 for(j = 1; j <= ncl[i];j++)
587 {
588 itest = 0;
589 ihld = iord[j];
590 for(i1 = 0; i1 < j; i1++)
591 {
592 if(itest == 0 && z[iord[i1]] < z[ihld])
593 {
594 itest = 1;
595 for(i2 = j-1;i2 >= i1;i2--)
596 {
597 iord[i2+1] = iord[i2];
598 }
599 iord[i1] = ihld;
600 }
601 }
602 }
603
604
605 // compute the number of clusters and their centers ( first
606 // guess )
607 // centers must be separated by cells having smaller ener. dep.
608 // neighbouring centers should be either strong or well-separated
609 ig = 0;
610 xc[ig] = x[iord[0]];
611 yc[ig] = y[iord[0]];
612 zc[ig] = z[iord[0]];
613 for(j = 1; j <= ncl[i]; j++)
614 {
615 itest = -1;
616 x1 = x[iord[j]];
617 y1 = y[iord[j]];
618 for(k = 0; k <= ig; k++)
619 {
620 x2 = xc[k];
621 y2 = yc[k];
622 rr = Distance(x1,y1,x2,y2);
623 //************************************************************
624 // finetuning cluster splitting
625 // the numbers zc/4 and zc/10 may need to be changed.
626 // Also one may need to add one more layer because our
627 // cells are smaller in absolute scale
628 //************************************************************
629
630
631 if( rr >= 1.1 && rr < 1.8 && z[iord[j]] > zc[k]/4.) itest++;
632 if( rr >= 1.8 && rr < 2.1 && z[iord[j]] > zc[k]/10.) itest++;
633 if( rr >= 2.1)itest++;
634 }
635
636 if(itest == ig)
637 {
638 ig++;
639 xc[ig] = x1;
640 yc[ig] = y1;
641 zc[ig] = z[iord[j]];
642 }
643 }
c1339151 644 ClustDetails(ncl[i], ig, x, y ,z, xc, yc, zc, rcl, rcs, cells,
645 testncl, testindex);
562718f9 646
647 Int_t pp = 0;
648 for(j = 0; j <= ig; j++)
649 {
650 clno++;
651 if (clno >= 5000)
652 {
653 AliWarning("RefClust: Too many clusters! more than 5000");
654 return;
655 }
656 clusdata[0] = xc[j];
657 clusdata[1] = yc[j];
658 clusdata[2] = zc[j];
659 clusdata[4] = rcl[j];
660 clusdata[5] = rcs[j];
661 if(ig == 0)
662 {
c1339151 663 clusdata[3] = ncl[i] + 1;//ajay
664 //clusdata[3] = ncl[i] ;
562718f9 665 }
666 else
667 {
668 clusdata[3] = cells[j];
669 }
670 // cell information
671 Int_t ncellcls = testncl[j];
672 if( ncellcls < 15 )
673 {
674 for(Int_t kk = 1; kk <= ncellcls; kk++)
675 {
676 Int_t ll = testindex[pp];
c1339151 677 clY = (Int_t)((ktwobysqrt3*y[ll])*10);
678 clX = (Int_t)((x[ll] - clY/20.)*10);
679 clxy[kk-1] = clX*10000 + clY ;
680
681
682 //clX = (Int_t) x[ll]*10;
683 //clY = (Int_t) y[ll]*10;
684 //clxy[kk-1] = (Int_t) clX*10000 + clY ;
562718f9 685 pp++;
686 }
687 for(Int_t icltr = ncellcls ; icltr < 15; icltr++)
688 {
689 clxy[icltr] = -1;
690 }
691 }
692 pmdcludata = new AliPMDcludata(clusdata, clxy);
693 fPMDclucont->Add(pmdcludata);
694 }
695 testncl.Set(0);
696 testindex.Set(0);
697 }
8c7250c5 698 }
c1339151 699 delete [] ncl;
700 delete [] iord;
701 delete [] x;
702 delete [] y;
703 delete [] z;
704 delete [] xc;
705 delete [] yc;
706 delete [] zc;
707 delete [] cells;
708 delete [] rcl;
709 delete [] rcs;
8c7250c5 710}
8c7250c5 711// ------------------------------------------------------------------------ //
c1339151 712void AliPMDClusteringV2::ClustDetails(Int_t ncell, Int_t nclust, Double_t x[],
713 Double_t y[], Double_t z[],Double_t xc[],
714 Double_t yc[], Double_t zc[],
715 Double_t rcl[], Double_t rcs[],
716 Double_t cells[], TArrayI &testncl,
562718f9 717 TArrayI &testindex)
8c7250c5 718{
719 // function begins
720 //
8c7250c5 721
c1339151 722 Int_t kndim1 = ncell + 1;//ncell
723 Int_t kndim2 = 20;
724 Int_t kndim3 = nclust + 1;//nclust
78fc1b96 725
562718f9 726 Int_t i, j, k, i1, i2;
562718f9 727 Double_t x1, y1, x2, y2, rr, b, c, r1, r2;
728 Double_t sumx, sumy, sumxy, sumxx, sum, sum1, sumyy;
8c7250c5 729
c1339151 730 Double_t *str, *str1, *xcl, *ycl, *cln;
731 Int_t **cell;
732 Int_t ** cluster;
733 Double_t **clustcell;
734 str = new Double_t [kndim3];
735 str1 = new Double_t [kndim3];
736 xcl = new Double_t [kndim3];
737 ycl = new Double_t [kndim3];
738 cln = new Double_t [kndim3];
562718f9 739
c1339151 740 clustcell = new Double_t *[kndim3];
741 cell = new Int_t *[kndim3];
742 cluster = new Int_t *[kndim1];
562718f9 743 for(i = 0; i < kndim1; i++)
744 {
c1339151 745 cluster[i] = new Int_t [kndim2];
746 }
747
748 for(i = 0; i < kndim3; i++)
749 {
750 str[i] = 0;
751 str1[i] = 0;
752 xcl[i] = 0;
753 ycl[i] = 0;
754 cln[i] = 0;
755
756 cell[i] = new Int_t [kndim2];
757 clustcell[i] = new Double_t [kndim1];
758 for(j = 0; j < kndim1; j++)
8c7250c5 759 {
c1339151 760 clustcell[i][j] = 0;
562718f9 761 }
c1339151 762 for(j = 0; j < kndim2; j++)
562718f9 763 {
c1339151 764 cluster[i][j] = 0;
765 cell[i][j] = 0;
8c7250c5 766 }
562718f9 767 }
c1339151 768
562718f9 769 if(nclust > 0)
770 {
771 // more than one cluster
772 // checking cells shared between several clusters.
773 // First check if the cell is within
774 // one cell unit ( nearest neighbour). Else,
775 // if it is within 1.74 cell units ( next nearest )
776 // Else if it is upto 2 cell units etc.
777
778 for (i = 0; i <= ncell; i++)
8c7250c5 779 {
c1339151 780 x1 = x[i];
781 y1 = y[i];
562718f9 782 cluster[i][0] = 0;
783 // distance <= 1 cell unit
784 for(j = 0; j <= nclust; j++)
8c7250c5 785 {
c1339151 786 x2 = xc[j];
787 y2 = yc[j];
8c7250c5 788 rr = Distance(x1, y1, x2, y2);
562718f9 789 if(rr <= 1.)
8c7250c5 790 {
791 cluster[i][0]++;
792 i1 = cluster[i][0];
793 cluster[i][i1] = j;
794 }
795 }
562718f9 796 // next nearest neighbour
797 if(cluster[i][0] == 0)
798 {
799 for(j=0; j<=nclust; j++)
800 {
c1339151 801 x2 = xc[j];
802 y2 = yc[j];
562718f9 803 rr = Distance(x1, y1, x2, y2);
804 if(rr <= TMath::Sqrt(3.))
805 {
806 cluster[i][0]++;
807 i1 = cluster[i][0];
808 cluster[i][i1] = j;
809 }
810 }
811 }
812 // next-to-next nearest neighbour
813 if(cluster[i][0] == 0)
814 {
815 for(j=0; j<=nclust; j++)
816 {
c1339151 817 x2 = xc[j];
818 y2 = yc[j];
562718f9 819 rr = Distance(x1, y1, x2, y2);
820 if(rr <= 2.)
821 {
822 cluster[i][0]++;
823 i1 = cluster[i][0];
824 cluster[i][i1] = j;
825 }
826 }
827 }
828 // one more
829 if(cluster[i][0] == 0)
830 {
831 for(j = 0; j <= nclust; j++)
832 {
c1339151 833 x2 = xc[j];
834 y2 = yc[j];
562718f9 835 rr = Distance(x1, y1, x2, y2);
836 if(rr <= 2.7)
837 {
838 cluster[i][0]++;
839 i1 = cluster[i][0];
840 cluster[i][i1] = j;
841 }
842 }
843 }
8c7250c5 844 }
562718f9 845
846 // computing cluster strength. Some cells are shared.
847 for(i = 0; i <= ncell; i++)
8c7250c5 848 {
562718f9 849 if(cluster[i][0] != 0)
8c7250c5 850 {
562718f9 851 i1 = cluster[i][0];
852 for(j = 1; j <= i1; j++)
8c7250c5 853 {
562718f9 854 i2 = cluster[i][j];
c1339151 855 str[i2] += z[i]/i1;
8c7250c5 856 }
857 }
858 }
562718f9 859
860 for(k = 0; k < 5; k++)
8c7250c5 861 {
562718f9 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 sum=0.;
868 for(j = 1; j <= i1; j++)
869 {
870 sum += str[cluster[i][j]];
871 }
872
873 for(j = 1; j <= i1; j++)
874 {
875 i2 = cluster[i][j];
c1339151 876 str1[i2] += z[i]*str[i2]/sum;
877 clustcell[i2][i] = z[i]*str[i2]/sum;
562718f9 878 }
8c7250c5 879 }
880 }
562718f9 881
882
883 for(j = 0; j <= nclust; j++)
884 {
885 str[j] = str1[j];
886 str1[j] = 0.;
887 }
8c7250c5 888 }
562718f9 889
890 for(i = 0; i <= nclust; i++)
891 {
892 sumx = 0.;
893 sumy = 0.;
894 sum = 0.;
895 sum1 = 0.;
896 for(j = 0; j <= ncell; j++)
897 {
898 if(clustcell[i][j] != 0)
899 {
c1339151 900 sumx += clustcell[i][j]*x[j];
901 sumy += clustcell[i][j]*y[j];
562718f9 902 sum += clustcell[i][j];
c1339151 903 sum1 += clustcell[i][j]/z[j];
562718f9 904 }
905 }
906 //** xcl and ycl are cluster centroid positions ( center of gravity )
907
908 xcl[i] = sumx/sum;
909 ycl[i] = sumy/sum;
910 cln[i] = sum1;
911 sumxx = 0.;
912 sumyy = 0.;
913 sumxy = 0.;
914 for(j = 0; j <= ncell; j++)
915 {
c1339151 916 sumxx += clustcell[i][j]*(x[j]-xcl[i])*(x[j]-xcl[i])/sum;
917 sumyy += clustcell[i][j]*(y[j]-ycl[i])*(y[j]-ycl[i])/sum;
918 sumxy += clustcell[i][j]*(x[j]-xcl[i])*(y[j]-ycl[i])/sum;
562718f9 919 }
920 b = sumxx+sumyy;
921 c = sumxx*sumyy-sumxy*sumxy;
922 // ******************r1 and r2 are major and minor axes ( r1 > r2 ).
923 r1 = b/2.+TMath::Sqrt(b*b/4.-c);
924 r2 = b/2.-TMath::Sqrt(b*b/4.-c);
925 // final assignments to proper external variables
c1339151 926 xc[i] = xcl[i];
927 yc[i] = ycl[i];
928 zc[i] = str[i];
929 cells[i] = cln[i];
930 rcl[i] = r1;
931 rcs[i] = r2;
932
8c7250c5 933 }
562718f9 934
935 //To get the cell position in a cluster
936
937 for(Int_t ii=0; ii<= ncell; ii++)
938 {
939 Int_t jj = cluster[ii][0];
940 for(Int_t kk=1; kk<= jj; kk++)
8c7250c5 941 {
562718f9 942 Int_t ll = cluster[ii][kk];
943 cell[ll][0]++;
944 cell[ll][cell[ll][0]] = ii;
8c7250c5 945 }
562718f9 946 }
947
948 testncl.Set(nclust+1);
949 Int_t counter = 0;
950
951 for(Int_t ii=0; ii <= nclust; ii++)
952 {
953 testncl[ii] = cell[ii][0];
954 counter += testncl[ii];
955 }
956 testindex.Set(counter);
957 Int_t ll = 0;
958 for(Int_t ii=0; ii<= nclust; ii++)
959 {
960 for(Int_t jj = 1; jj<= testncl[ii]; jj++)
961 {
962 Int_t kk = cell[ii][jj];
963 testindex[ll] = kk;
964 ll++;
965 }
966 }
967
968 }
c1339151 969 else if(nclust == 0)
562718f9 970 {
8c7250c5 971 sumx = 0.;
972 sumy = 0.;
973 sum = 0.;
974 sum1 = 0.;
562718f9 975 i = 0;
976 for(j = 0; j <= ncell; j++)
977 {
c1339151 978 sumx += z[j]*x[j];
979 sumy += z[j]*y[j];
980 sum += z[j];
562718f9 981 sum1++;
8c7250c5 982 }
8c7250c5 983 xcl[i] = sumx/sum;
984 ycl[i] = sumy/sum;
985 cln[i] = sum1;
562718f9 986 sumxx = 0.;
987 sumyy = 0.;
988 sumxy = 0.;
989 for(j = 0; j <= ncell; j++)
990 {
c1339151 991 sumxx += clustcell[i][j]*(x[j]-xcl[i])*(x[j]-xcl[i])/sum;
992 sumyy += clustcell[i][j]*(y[j]-ycl[i])*(y[j]-ycl[i])/sum;
993 sumxy += clustcell[i][j]*(x[j]-xcl[i])*(y[j]-ycl[i])/sum;
562718f9 994 }
995 b = sumxx+sumyy;
996 c = sumxx*sumyy-sumxy*sumxy;
997 r1 = b/2.+ TMath::Sqrt(b*b/4.-c);
998 r2 = b/2.- TMath::Sqrt(b*b/4.-c);
999
1000 // To get the cell position in a cluster
1001 testncl.Set(nclust+1);
c1339151 1002 //testindex.Set(ncell);
1003 testindex.Set(ncell+1);
1004 cell[0][0] = ncell + 1;//ajay
1005 //cell[0][0] = ncell;
562718f9 1006 testncl[0] = cell[0][0];
1007 Int_t ll = 0;
c1339151 1008 for(Int_t ii = 1; ii <= ncell; ii++)
562718f9 1009 {
1010 cell[0][ii]=ii;
1011 //clustcell[0][ii]=1.;
1012 Int_t kk = cell[0][ii];
1013 testindex[ll] = kk;
1014 ll++;
1015 }
1016 // final assignments
c1339151 1017 xc[i] = xcl[i];
1018 yc[i] = ycl[i];
1019 //zc[i] = str[i];//ajay
1020 zc[i] = sum;
1021 cells[i] = cln[i];
1022 rcl[i] = r1;
1023 rcs[i] = r2;
1024 }
1025 for(i = 0; i < kndim3; i++)
1026 {
1027 delete [] clustcell[i];
1028 delete [] cell[i];
1029 }
1030 delete [] clustcell;
1031 delete [] cell;
1032 for(i = 0; i <kndim1 ; i++)
1033 {
1034 delete [] cluster[i];
8c7250c5 1035 }
c1339151 1036 delete [] cluster;
1037 delete [] str;
1038 delete [] str1;
1039 delete [] xcl;
1040 delete [] ycl;
1041 delete [] cln;
8c7250c5 1042}
1043
1044// ------------------------------------------------------------------------ //
1045Double_t AliPMDClusteringV2::Distance(Double_t x1, Double_t y1,
1046 Double_t x2, Double_t y2)
1047{
562718f9 1048 return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
8c7250c5 1049}
1050// ------------------------------------------------------------------------ //
1051void AliPMDClusteringV2::SetEdepCut(Float_t decut)
1052{
1053 fCutoff = decut;
1054}
1055// ------------------------------------------------------------------------ //