]> 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
Updating macro to check the various methods of the AliMUONMCDataInterface
[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];
148 TMath::Sort(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;
222 Float_t cellY = (Float_t) (ktwobysqrt3*celldumY/10);
223 Float_t cellX = (Float_t) (celldumX/10 - (celldumY/2.)/10);
224
225 //
226 // Cell X centroid is back transformed
227 //
228 if (ismn < 12)
229 {
230 celldataX[ihit] = (Int_t) (cellX - (24-1) + cellY/2.);
231 }
232 else if (ismn >= 12 && ismn <= 23)
233 {
234 celldataX[ihit] = (Int_t) (cellX - (48-1) + cellY/2.);
235 }
236 celldataY[ihit] = (Int_t) cellY;
8c7250c5 237 }
238
239 pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY);
240 pmdcont->Add(pmdcl);
241 }
562718f9 242 fPMDclucont->Clear();
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
562718f9 384 AliPMDcludata *pmdcludata = 0;
385 TArrayI testncl;
386 TArrayI testindex;
8c7250c5 387 const Int_t kndim = 4500;
562718f9 388 Int_t i, j, k, i1, i2, id, icl, itest, ihld;
389 Int_t ig, nsupcl, clno, clX,clY;
390 Int_t clxy[15];
391 Int_t ncl[kndim], iord[kndim];
392 Float_t clusdata[6];
393 Double_t x1, y1, z1, x2, y2, z2, rr;
394 Double_t x[kndim], y[kndim], z[kndim];
395 Double_t xc[kndim], yc[kndim], zc[kndim], cells[kndim];
396 Double_t rcl[kndim], rcs[kndim];
397
398 for(Int_t kk = 0; kk < 15; kk++)
399 {
400 if( kk < 6 )clusdata[kk] = 0.;
401 }
402
8c7250c5 403 // nsupcl = # of superclusters; ncl[i]= # of cells in supercluster i
404 // x, y and z store (x,y) coordinates of and energy deposited in a cell
405 // xc, yc store (x,y) coordinates of the cluster center
562718f9 406 // zc stores the energy deposited in a cluster, rc is cluster radius
8c7250c5 407
562718f9 408 clno = -1;
8c7250c5 409 nsupcl = -1;
562718f9 410 for(i = 0; i < 4500; i++)
411 {
412 ncl[i] = -1;
8c7250c5 413 }
562718f9 414 for(i = 0; i < incr; i++)
415 {
416 if(fInfcl[0][i] != nsupcl)
417 {
418 nsupcl++;
419 }
420 if (nsupcl > 4500)
421 {
422 AliWarning("RefClust: Too many superclusters!");
423 nsupcl = 4500;
424 break;
425 }
426 ncl[nsupcl]++;
427 }
428
8c7250c5 429 AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
430 incr+1,nsupcl+1));
562718f9 431
432 id = -1;
433 icl = -1;
434 for(i = 0; i < nsupcl; i++)
435 {
436 if(ncl[i] == 0)
8c7250c5 437 {
562718f9 438 id++;
439 icl++;
440 // one cell super-clusters --> single cluster
441 // cluster center at the centyer of the cell
442 // cluster radius = half cell dimension
443 if (clno >= 5000)
8c7250c5 444 {
562718f9 445 AliWarning("RefClust: Too many clusters! more than 5000");
446 return;
8c7250c5 447 }
562718f9 448 clno++;
449 i1 = fInfcl[1][id];
450 i2 = fInfcl[2][id];
451 Int_t i12 = i1 + i2*kNDIMX;
452 clusdata[0] = fCoord[0][i1][i2];
453 clusdata[1] = fCoord[1][i1][i2];
454 clusdata[2] = edepcell[i12];
455 clusdata[3] = 1.;
456 clusdata[4] = 0.0;
457 clusdata[5] = 0.0;
458
459 //cell information
460 clX = (Int_t) fCoord[0][i1][i2]*10;
461 clY = (Int_t) fCoord[1][i1][i2]*10;
462 clxy[0] = clX*10000 + clY;
463 for(Int_t icltr = 1; icltr < 15; icltr++)
464 {
465 clxy[icltr] = -1;
466 }
467 pmdcludata = new AliPMDcludata(clusdata,clxy);
468 fPMDclucont->Add(pmdcludata);
469
470
8c7250c5 471 }
562718f9 472 else if(ncl[i] == 1)
8c7250c5 473 {
562718f9 474 // two cell super-cluster --> single cluster
475 // cluster center is at ener. dep.-weighted mean of two cells
476 // cluster radius == half cell dimension
477 id++;
478 icl++;
479 if (clno >= 5000)
8c7250c5 480 {
481 AliWarning("RefClust: Too many clusters! more than 5000");
482 return;
483 }
562718f9 484 clno++;
485 i1 = fInfcl[1][id];
486 i2 = fInfcl[2][id];
487 Int_t i12 = i1 + i2*kNDIMX;
488
489 x1 = fCoord[0][i1][i2];
490 y1 = fCoord[1][i1][i2];
491 z1 = edepcell[i12];
492
493 id++;
494 i1 = fInfcl[1][id];
495 i2 = fInfcl[2][id];
496 i12 = i1 + i2*kNDIMX;
497
498 x2 = fCoord[0][i1][i2];
499 y2 = fCoord[1][i1][i2];
500 z2 = edepcell[i12];
501
502 clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
503 clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
504 clusdata[2] = z1+z2;
505 clusdata[3] = 2.;
506 clusdata[4] = (TMath::Sqrt(z1*z2))/(z1+z2);
507 clusdata[5] = 0.0;
508
509 clX = (Int_t) x1*10;
510 clY = (Int_t) y1*10;
511 clxy[0] = clX*10000 + clY;
512
513 clX = (Int_t) x2*10;
514 clY = (Int_t) y2*10;
515 clxy[1] = clX*10000 + clY;
516
517 for(Int_t icltr = 2; icltr < 15; icltr++)
8c7250c5 518 {
562718f9 519 clxy[icltr] = -1;
8c7250c5 520 }
562718f9 521 pmdcludata = new AliPMDcludata(clusdata, clxy);
522 fPMDclucont->Add(pmdcludata);
8c7250c5 523 }
562718f9 524 else{
525 id++;
526 iord[0] = 0;
527 // super-cluster of more than two cells - broken up into smaller
528 // clusters gaussian centers computed. (peaks separated by > 1 cell)
529 // Begin from cell having largest energy deposited This is first
530 // cluster center
531 // *****************************************************************
532 // NOTE --- POSSIBLE MODIFICATION: ONE MAY NOT BREAKING SUPERCLUSTERS
533 // IF NO. OF CELLS IS NOT TOO LARGE ( SAY 5 OR 6 )
534 // SINCE WE EXPECT THE SUPERCLUSTER
535 // TO BE A SINGLE CLUSTER
536 //*******************************************************************
537
538 i1 = fInfcl[1][id];
539 i2 = fInfcl[2][id];
540 Int_t i12 = i1 + i2*kNDIMX;
541
542 x[0] = fCoord[0][i1][i2];
543 y[0] = fCoord[1][i1][i2];
544 z[0] = edepcell[i12];
545
546 iord[0] = 0;
547 for(j = 1; j <= ncl[i]; j++)
548 {
549
550 id++;
551 i1 = fInfcl[1][id];
552 i2 = fInfcl[2][id];
553 Int_t i12 = i1 + i2*kNDIMX;
554 iord[j] = j;
555 x[j] = fCoord[0][i1][i2];
556 y[j] = fCoord[1][i1][i2];
557 z[j] = edepcell[i12];
558 }
559
560 // arranging cells within supercluster in decreasing order
561 for(j = 1; j <= ncl[i];j++)
562 {
563 itest = 0;
564 ihld = iord[j];
565 for(i1 = 0; i1 < j; i1++)
566 {
567 if(itest == 0 && z[iord[i1]] < z[ihld])
568 {
569 itest = 1;
570 for(i2 = j-1;i2 >= i1;i2--)
571 {
572 iord[i2+1] = iord[i2];
573 }
574 iord[i1] = ihld;
575 }
576 }
577 }
578
579
580 // compute the number of clusters and their centers ( first
581 // guess )
582 // centers must be separated by cells having smaller ener. dep.
583 // neighbouring centers should be either strong or well-separated
584 ig = 0;
585 xc[ig] = x[iord[0]];
586 yc[ig] = y[iord[0]];
587 zc[ig] = z[iord[0]];
588 for(j = 1; j <= ncl[i]; j++)
589 {
590 itest = -1;
591 x1 = x[iord[j]];
592 y1 = y[iord[j]];
593 for(k = 0; k <= ig; k++)
594 {
595 x2 = xc[k];
596 y2 = yc[k];
597 rr = Distance(x1,y1,x2,y2);
598 //************************************************************
599 // finetuning cluster splitting
600 // the numbers zc/4 and zc/10 may need to be changed.
601 // Also one may need to add one more layer because our
602 // cells are smaller in absolute scale
603 //************************************************************
604
605
606 if( rr >= 1.1 && rr < 1.8 && z[iord[j]] > zc[k]/4.) itest++;
607 if( rr >= 1.8 && rr < 2.1 && z[iord[j]] > zc[k]/10.) itest++;
608 if( rr >= 2.1)itest++;
609 }
610
611 if(itest == ig)
612 {
613 ig++;
614 xc[ig] = x1;
615 yc[ig] = y1;
616 zc[ig] = z[iord[j]];
617 }
618 }
619 ClustDetails(ncl[i], ig, x[0], y[0] ,z[0], xc[0], yc[0], zc[0],
620 rcl[0], rcs[0], cells[0], testncl, testindex);
621
622 Int_t pp = 0;
623 for(j = 0; j <= ig; j++)
624 {
625 clno++;
626 if (clno >= 5000)
627 {
628 AliWarning("RefClust: Too many clusters! more than 5000");
629 return;
630 }
631 clusdata[0] = xc[j];
632 clusdata[1] = yc[j];
633 clusdata[2] = zc[j];
634 clusdata[4] = rcl[j];
635 clusdata[5] = rcs[j];
636 if(ig == 0)
637 {
638 clusdata[3] = ncl[i];
639 }
640 else
641 {
642 clusdata[3] = cells[j];
643 }
644 // cell information
645 Int_t ncellcls = testncl[j];
646 if( ncellcls < 15 )
647 {
648 for(Int_t kk = 1; kk <= ncellcls; kk++)
649 {
650 Int_t ll = testindex[pp];
651 clX = (Int_t) x[ll]*10;
652 clY = (Int_t) y[ll]*10;
653 clxy[kk-1] = (Int_t) clX*10000 + clY ;
654 pp++;
655 }
656 for(Int_t icltr = ncellcls ; icltr < 15; icltr++)
657 {
658 clxy[icltr] = -1;
659 }
660 }
661 pmdcludata = new AliPMDcludata(clusdata, clxy);
662 fPMDclucont->Add(pmdcludata);
663 }
664 testncl.Set(0);
665 testindex.Set(0);
666 }
8c7250c5 667 }
8c7250c5 668}
8c7250c5 669// ------------------------------------------------------------------------ //
562718f9 670void AliPMDClusteringV2::ClustDetails(Int_t ncell, Int_t nclust, Double_t &x,
671 Double_t &y, Double_t &z, Double_t &xc,
672 Double_t &yc, Double_t &zc,
673 Double_t &rcl, Double_t &rcs,
674 Double_t &cells, TArrayI &testncl,
675 TArrayI &testindex)
8c7250c5 676{
677 // function begins
678 //
8c7250c5 679
562718f9 680 const Int_t kndim1 = 2000;
681 const Int_t kndim2 = 10;
682 const Int_t kndim3 = 400;
8c7250c5 683
562718f9 684 Int_t i, j, k, i1, i2;
685 Int_t cluster[kndim1][kndim2], cell[kndim1][kndim3];
686
687 Double_t x1, y1, x2, y2, rr, b, c, r1, r2;
688 Double_t sumx, sumy, sumxy, sumxx, sum, sum1, sumyy;
689 Double_t xx[kndim1], yy[kndim1], zz[kndim1];
690 Double_t xxc[kndim1], yyc[kndim1],clustcell[kndim3][kndim1];
691 Double_t str[kndim1], str1[kndim1],xcl[kndim1], ycl[kndim1], cln[kndim1];
8c7250c5 692
562718f9 693 for(i = 0; i <= nclust; i++)
694 {
695 xxc[i] = *(&xc+i);
696 yyc[i] = *(&yc+i);
697 str[i] = 0.;
698 str1[i] = 0.;
8c7250c5 699 }
562718f9 700 for(i = 0; i <= ncell; i++)
701 {
702 xx[i] = *(&x+i);
703 yy[i] = *(&y+i);
704 zz[i] = *(&z+i);
705 }
706
707 // INITIALIZE
708
709 for(i = 0; i < kndim1; i++)
710 {
711 for(j = 0; j < kndim2; j++)
8c7250c5 712 {
562718f9 713 cluster[i][j] = 0;
714
715 }
716 for(j = 0; j < kndim3; j++)
717 {
718 cell[i][j] = 0;
719 clustcell[j][i] = 0.;
8c7250c5 720 }
562718f9 721 }
722
723
724 if(nclust > 0)
725 {
726 // more than one cluster
727 // checking cells shared between several clusters.
728 // First check if the cell is within
729 // one cell unit ( nearest neighbour). Else,
730 // if it is within 1.74 cell units ( next nearest )
731 // Else if it is upto 2 cell units etc.
732
733 for (i = 0; i <= ncell; i++)
8c7250c5 734 {
562718f9 735 x1 = xx[i];
736 y1 = yy[i];
737 cluster[i][0] = 0;
738 // distance <= 1 cell unit
739 for(j = 0; j <= nclust; j++)
8c7250c5 740 {
741 x2 = xxc[j];
742 y2 = yyc[j];
743 rr = Distance(x1, y1, x2, y2);
562718f9 744 if(rr <= 1.)
8c7250c5 745 {
746 cluster[i][0]++;
747 i1 = cluster[i][0];
748 cluster[i][i1] = j;
749 }
750 }
562718f9 751 // next nearest neighbour
752 if(cluster[i][0] == 0)
753 {
754 for(j=0; j<=nclust; j++)
755 {
756 x2 = xxc[j];
757 y2 = yyc[j];
758 rr = Distance(x1, y1, x2, y2);
759 if(rr <= TMath::Sqrt(3.))
760 {
761 cluster[i][0]++;
762 i1 = cluster[i][0];
763 cluster[i][i1] = j;
764 }
765 }
766 }
767 // next-to-next nearest neighbour
768 if(cluster[i][0] == 0)
769 {
770 for(j=0; j<=nclust; j++)
771 {
772 x2 = xxc[j];
773 y2 = yyc[j];
774 rr = Distance(x1, y1, x2, y2);
775 if(rr <= 2.)
776 {
777 cluster[i][0]++;
778 i1 = cluster[i][0];
779 cluster[i][i1] = j;
780 }
781 }
782 }
783 // one more
784 if(cluster[i][0] == 0)
785 {
786 for(j = 0; j <= nclust; j++)
787 {
788 x2 = xxc[j];
789 y2 = yyc[j];
790 rr = Distance(x1, y1, x2, y2);
791 if(rr <= 2.7)
792 {
793 cluster[i][0]++;
794 i1 = cluster[i][0];
795 cluster[i][i1] = j;
796 }
797 }
798 }
8c7250c5 799 }
562718f9 800
801 // computing cluster strength. Some cells are shared.
802 for(i = 0; i <= ncell; i++)
8c7250c5 803 {
562718f9 804 if(cluster[i][0] != 0)
8c7250c5 805 {
562718f9 806 i1 = cluster[i][0];
807 for(j = 1; j <= i1; j++)
8c7250c5 808 {
562718f9 809 i2 = cluster[i][j];
810 str[i2] += zz[i]/i1;
8c7250c5 811 }
812 }
813 }
562718f9 814
815 for(k = 0; k < 5; k++)
8c7250c5 816 {
562718f9 817 for(i = 0; i <= ncell; i++)
8c7250c5 818 {
562718f9 819 if(cluster[i][0] != 0)
8c7250c5 820 {
562718f9 821 i1=cluster[i][0];
822 sum=0.;
823 for(j = 1; j <= i1; j++)
824 {
825 sum += str[cluster[i][j]];
826 }
827
828 for(j = 1; j <= i1; j++)
829 {
830 i2 = cluster[i][j];
831 str1[i2] += zz[i]*str[i2]/sum;
832 clustcell[i2][i] = zz[i]*str[i2]/sum;
833 }
8c7250c5 834 }
835 }
562718f9 836
837
838 for(j = 0; j <= nclust; j++)
839 {
840 str[j] = str1[j];
841 str1[j] = 0.;
842 }
8c7250c5 843 }
562718f9 844
845 for(i = 0; i <= nclust; i++)
846 {
847 sumx = 0.;
848 sumy = 0.;
849 sum = 0.;
850 sum1 = 0.;
851 for(j = 0; j <= ncell; j++)
852 {
853 if(clustcell[i][j] != 0)
854 {
855 sumx += clustcell[i][j]*xx[j];
856 sumy += clustcell[i][j]*yy[j];
857 sum += clustcell[i][j];
858 sum1 += clustcell[i][j]/zz[j];
859 }
860 }
861 //** xcl and ycl are cluster centroid positions ( center of gravity )
862
863 xcl[i] = sumx/sum;
864 ycl[i] = sumy/sum;
865 cln[i] = sum1;
866 sumxx = 0.;
867 sumyy = 0.;
868 sumxy = 0.;
869 for(j = 0; j <= ncell; j++)
870 {
871 sumxx += clustcell[i][j]*(xx[j]-xcl[i])*(xx[j]-xcl[i])/sum;
872 sumyy += clustcell[i][j]*(yy[j]-ycl[i])*(yy[j]-ycl[i])/sum;
873 sumxy += clustcell[i][j]*(xx[j]-xcl[i])*(yy[j]-ycl[i])/sum;
874 }
875 b = sumxx+sumyy;
876 c = sumxx*sumyy-sumxy*sumxy;
877 // ******************r1 and r2 are major and minor axes ( r1 > r2 ).
878 r1 = b/2.+TMath::Sqrt(b*b/4.-c);
879 r2 = b/2.-TMath::Sqrt(b*b/4.-c);
880 // final assignments to proper external variables
881 *(&xc + i) = xcl[i];
882 *(&yc + i) = ycl[i];
883 *(&zc + i) = str[i];
884 *(&cells + i) = cln[i];
885 *(&rcl+i) = r1;
886 *(&rcs+i) = r2;
8c7250c5 887 }
562718f9 888
889 //To get the cell position in a cluster
890
891 for(Int_t ii=0; ii<= ncell; ii++)
892 {
893 Int_t jj = cluster[ii][0];
894 for(Int_t kk=1; kk<= jj; kk++)
8c7250c5 895 {
562718f9 896 Int_t ll = cluster[ii][kk];
897 cell[ll][0]++;
898 cell[ll][cell[ll][0]] = ii;
8c7250c5 899 }
562718f9 900 }
901
902 testncl.Set(nclust+1);
903 Int_t counter = 0;
904
905 for(Int_t ii=0; ii <= nclust; ii++)
906 {
907 testncl[ii] = cell[ii][0];
908 counter += testncl[ii];
909 }
910 testindex.Set(counter);
911 Int_t ll = 0;
912 for(Int_t ii=0; ii<= nclust; ii++)
913 {
914 for(Int_t jj = 1; jj<= testncl[ii]; jj++)
915 {
916 Int_t kk = cell[ii][jj];
917 testindex[ll] = kk;
918 ll++;
919 }
920 }
921
922 }
923 else
924 {
8c7250c5 925 sumx = 0.;
926 sumy = 0.;
927 sum = 0.;
928 sum1 = 0.;
562718f9 929 i = 0;
930 for(j = 0; j <= ncell; j++)
931 {
932 sumx += zz[j]*xx[j];
933 sumy += zz[j]*yy[j];
934 sum += zz[j];
935 sum1++;
8c7250c5 936 }
8c7250c5 937 xcl[i] = sumx/sum;
938 ycl[i] = sumy/sum;
939 cln[i] = sum1;
562718f9 940 sumxx = 0.;
941 sumyy = 0.;
942 sumxy = 0.;
943 for(j = 0; j <= ncell; j++)
944 {
945 sumxx += clustcell[i][j]*(xx[j]-xcl[i])*(xx[j]-xcl[i])/sum;
946 sumyy += clustcell[i][j]*(yy[j]-ycl[i])*(yy[j]-ycl[i])/sum;
947 sumxy += clustcell[i][j]*(xx[j]-xcl[i])*(yy[j]-ycl[i])/sum;
948 }
949 b = sumxx+sumyy;
950 c = sumxx*sumyy-sumxy*sumxy;
951 r1 = b/2.+ TMath::Sqrt(b*b/4.-c);
952 r2 = b/2.- TMath::Sqrt(b*b/4.-c);
953
954 // To get the cell position in a cluster
955 testncl.Set(nclust+1);
956 testindex.Set(ncell);
957 cell[0][0] = ncell;
958 testncl[0] = cell[0][0];
959 Int_t ll = 0;
960 for(Int_t ii=1; ii<=ncell; ii++)
961 {
962 cell[0][ii]=ii;
963 //clustcell[0][ii]=1.;
964 Int_t kk = cell[0][ii];
965 testindex[ll] = kk;
966 ll++;
967 }
968 // final assignments
969 *(&xc + i) = xcl[i];
970 *(&yc + i) = ycl[i];
971 *(&zc + i) = str[i];
8c7250c5 972 *(&cells + i) = cln[i];
562718f9 973 *(&rcl+i) = r1;
974 *(&rcs+i) = r2;
8c7250c5 975 }
8c7250c5 976}
977
978// ------------------------------------------------------------------------ //
979Double_t AliPMDClusteringV2::Distance(Double_t x1, Double_t y1,
980 Double_t x2, Double_t y2)
981{
562718f9 982 return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
8c7250c5 983}
984// ------------------------------------------------------------------------ //
985void AliPMDClusteringV2::SetEdepCut(Float_t decut)
986{
987 fCutoff = decut;
988}
989// ------------------------------------------------------------------------ //