]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PMD/AliPMDClusteringV2.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
AliTOFRecoParam now inheriting from AliDetectorRecoParam.
[u/mrichter/AliRoot.git] / PMD / AliPMDClusteringV2.cxx
... / ...
CommitLineData
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
16/* $Id$ */
17
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
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.
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
40#include <Riostream.h>
41#include <TMath.h>
42#include <TObjArray.h>
43#include <TArrayI.h>
44
45#include "AliPMDcludata.h"
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():
56 fPMDclucont(new TObjArray()),
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;
65 }
66 }
67}
68// ------------------------------------------------------------------------ //
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// ------------------------------------------------------------------------ //
88AliPMDClusteringV2::~AliPMDClusteringV2()
89{
90 delete fPMDclucont;
91}
92// ------------------------------------------------------------------------ //
93
94void AliPMDClusteringV2::DoClust(Int_t idet, Int_t ismn,
95 Double_t celladc[48][96], TObjArray *pmdcont)
96{
97 // main function to call other necessary functions to do clustering
98 //
99 AliPMDcluster *pmdcl = 0;
100
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;
106 Int_t celldataX[15], celldataY[15];
107 Float_t clusdata[6];
108 Double_t cutoff, ave;
109 Double_t edepcell[kNMX];
110
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 }
122
123 for (i =0; i < kNMX; i++)
124 {
125 edepcell[i] = 0.;
126 }
127
128 for (id = 0; id < ndimXr; id++)
129 {
130 for (jd = 0; jd < ndimYr; jd++)
131 {
132 j = jd;
133 i = id + (ndimYr/2-1) - (jd/2);
134 Int_t ij = i + j*kNDIMX;
135 if (ismn < 12)
136 {
137 edepcell[ij] = celladc[jd][id];
138 }
139 else if (ismn >= 12 && ismn <= 23)
140 {
141 edepcell[ij] = celladc[id][jd];
142 }
143
144 }
145 }
146
147 Int_t iord1[kNMX];
148 TMath::Sort((Int_t)kNMX,edepcell,iord1);// order the data
149 cutoff = fCutoff; // cutoff used to discard cells having ener. dep.
150 ave = 0.;
151 nmx1 = -1;
152
153 for(i = 0;i < kNMX; i++)
154 {
155 if(edepcell[i] > 0.)
156 {
157 ave += edepcell[i];
158 }
159 if(edepcell[i] > cutoff )
160 {
161 nmx1++;
162 }
163 }
164
165 AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
166
167 if (nmx1 == 0)
168 {
169 nmx1 = 1;
170 }
171 ave = ave/nmx1;
172
173 AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
174 kNMX,ave));
175
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++)
183 {
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
193 Float_t cluY0 = ktwobysqrt3*cluYC;
194 Float_t cluX0 = cluXC - cluY0/2.;
195
196 //
197 // Cluster X centroid is back transformed
198 //
199 if (ismn < 12)
200 {
201 clusdata[0] = cluX0 - (24-1) + cluY0/2.;
202 }
203 else if (ismn >= 12 && ismn <= 23)
204 {
205 clusdata[0] = cluX0 - (48-1) + cluY0/2.;
206 }
207
208 clusdata[1] = cluY0;
209 clusdata[2] = cluADC;
210 clusdata[3] = cluCELLS;
211 clusdata[4] = cluSIGX;
212 clusdata[5] = cluSIGY;
213 //
214 // Cells associated with a cluster
215 //
216 for (Int_t ihit = 0; ihit < 15; ihit++)
217 {
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 Float_t cellY = (Float_t) celldumY/10;
226 Float_t cellX = (Float_t) celldumX/10;
227
228
229 //
230 // Cell X centroid is back transformed
231 //
232 if (ismn < 12)
233 {
234 celldataX[ihit] = (Int_t) ((cellX - (24-1) + cellY/2.) + 0.5);
235 }
236 else if (ismn >= 12 && ismn <= 23)
237 {
238 celldataX[ihit] = (Int_t) ((cellX - (48-1) + cellY/2.) + 0.5 );
239 }
240 celldataY[ihit] = (Int_t) (cellY + 0.5);
241 }
242
243 pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY);
244 pmdcont->Add(pmdcl);
245 }
246 fPMDclucont->Clear();
247}
248// ------------------------------------------------------------------------ //
249Int_t AliPMDClusteringV2::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
250 Int_t iord1[], Double_t edepcell[])
251{
252 // Does crude clustering
253 // Finds out only the big patch by just searching the
254 // connected cells
255 //
256
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};
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
266
267 AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
268
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 }
290 }
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 // ---------------------------------------------------------------
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];
332 if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
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 }
376 }
377 }
378 }
379 return cellcount;
380}
381// ------------------------------------------------------------------------ //
382 void AliPMDClusteringV2::RefClust(Int_t incr, Double_t edepcell[])
383{
384 // Does the refining of clusters
385 // Takes the big patch and does gaussian fitting and
386 // finds out the more refined clusters
387
388 const Float_t ktwobysqrt3 = 1.1547;
389 AliPMDcludata *pmdcludata = 0;
390
391 Int_t i12;
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];
395
396 Float_t clusdata[6];
397 Double_t x1, y1, z1, x2, y2, z2, rr;
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];
419
420 for(Int_t kk = 0; kk < 15; kk++)
421 {
422 if( kk < 6 )clusdata[kk] = 0.;
423 }
424
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
428 // zc stores the energy deposited in a cluster, rc is cluster radius
429
430 clno = -1;
431 nsupcl = -1;
432
433 for(i = 0; i < kndim; i++)
434 {
435 ncl[i] = -1;
436 }
437 for(i = 0; i <= incr; i++)
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
452 AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
453 incr+1,nsupcl+1));
454
455 id = -1;
456 icl = -1;
457 for(i = 0; i <= nsupcl; i++)
458 {
459 if(ncl[i] == 0)
460 {
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)
467 {
468 AliWarning("RefClust: Too many clusters! more than 5000");
469 return;
470 }
471 clno++;
472 i1 = fInfcl[1][id];
473 i2 = fInfcl[2][id];
474 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
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
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
496 }
497 else if(ncl[i] == 1)
498 {
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)
505 {
506 AliWarning("RefClust: Too many clusters! more than 5000");
507 return;
508 }
509 clno++;
510 i1 = fInfcl[1][id];
511 i2 = fInfcl[2][id];
512 i12 = i1 + i2*kNDIMX;
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
534 clY = (Int_t)((ktwobysqrt3*y1)*10);
535 clX = (Int_t)((x1 - clY/20.)*10);
536 clxy[0] = clX*10000 + clY ;
537
538 clY = (Int_t)((ktwobysqrt3*y2)*10);
539 clX = (Int_t)((x2 - clY/20.)*10);
540 clxy[1] = clX*10000 + clY ;
541
542 for(Int_t icltr = 2; icltr < 15; icltr++)
543 {
544 clxy[icltr] = -1;
545 }
546 pmdcludata = new AliPMDcludata(clusdata, clxy);
547 fPMDclucont->Add(pmdcludata);
548 }
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];
565 i12 = i1 + i2*kNDIMX;
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];
578 i12 = i1 + i2*kNDIMX;
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 }
644 ClustDetails(ncl[i], ig, x, y ,z, xc, yc, zc, rcl, rcs, cells,
645 testncl, testindex);
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 {
663 clusdata[3] = ncl[i] + 1;//ajay
664 //clusdata[3] = ncl[i] ;
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];
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 ;
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 }
698 }
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;
710}
711// ------------------------------------------------------------------------ //
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,
717 TArrayI &testindex)
718{
719 // function begins
720 //
721
722 Int_t kndim1 = ncell + 1;//ncell
723 Int_t kndim2 = 20;
724 Int_t kndim3 = nclust + 1;//nclust
725
726 Int_t i, j, k, i1, i2;
727 Double_t x1, y1, x2, y2, rr, b, c, r1, r2;
728 Double_t sumx, sumy, sumxy, sumxx, sum, sum1, sumyy;
729
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];
739
740 clustcell = new Double_t *[kndim3];
741 cell = new Int_t *[kndim3];
742 cluster = new Int_t *[kndim1];
743 for(i = 0; i < kndim1; i++)
744 {
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++)
759 {
760 clustcell[i][j] = 0;
761 }
762 for(j = 0; j < kndim2; j++)
763 {
764 cluster[i][j] = 0;
765 cell[i][j] = 0;
766 }
767 }
768
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++)
779 {
780 x1 = x[i];
781 y1 = y[i];
782 cluster[i][0] = 0;
783 // distance <= 1 cell unit
784 for(j = 0; j <= nclust; j++)
785 {
786 x2 = xc[j];
787 y2 = yc[j];
788 rr = Distance(x1, y1, x2, y2);
789 if(rr <= 1.)
790 {
791 cluster[i][0]++;
792 i1 = cluster[i][0];
793 cluster[i][i1] = j;
794 }
795 }
796 // next nearest neighbour
797 if(cluster[i][0] == 0)
798 {
799 for(j=0; j<=nclust; j++)
800 {
801 x2 = xc[j];
802 y2 = yc[j];
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 {
817 x2 = xc[j];
818 y2 = yc[j];
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 {
833 x2 = xc[j];
834 y2 = yc[j];
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 }
844 }
845
846 // computing cluster strength. Some cells are shared.
847 for(i = 0; i <= ncell; i++)
848 {
849 if(cluster[i][0] != 0)
850 {
851 i1 = cluster[i][0];
852 for(j = 1; j <= i1; j++)
853 {
854 i2 = cluster[i][j];
855 str[i2] += z[i]/i1;
856 }
857 }
858 }
859
860 for(k = 0; k < 5; k++)
861 {
862 for(i = 0; i <= ncell; i++)
863 {
864 if(cluster[i][0] != 0)
865 {
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];
876 str1[i2] += z[i]*str[i2]/sum;
877 clustcell[i2][i] = z[i]*str[i2]/sum;
878 }
879 }
880 }
881
882
883 for(j = 0; j <= nclust; j++)
884 {
885 str[j] = str1[j];
886 str1[j] = 0.;
887 }
888 }
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 {
900 sumx += clustcell[i][j]*x[j];
901 sumy += clustcell[i][j]*y[j];
902 sum += clustcell[i][j];
903 sum1 += clustcell[i][j]/z[j];
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 {
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;
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
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
933 }
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++)
941 {
942 Int_t ll = cluster[ii][kk];
943 cell[ll][0]++;
944 cell[ll][cell[ll][0]] = ii;
945 }
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 }
969 else if(nclust == 0)
970 {
971 sumx = 0.;
972 sumy = 0.;
973 sum = 0.;
974 sum1 = 0.;
975 i = 0;
976 for(j = 0; j <= ncell; j++)
977 {
978 sumx += z[j]*x[j];
979 sumy += z[j]*y[j];
980 sum += z[j];
981 sum1++;
982 }
983 xcl[i] = sumx/sum;
984 ycl[i] = sumy/sum;
985 cln[i] = sum1;
986 sumxx = 0.;
987 sumyy = 0.;
988 sumxy = 0.;
989 for(j = 0; j <= ncell; j++)
990 {
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;
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);
1002 //testindex.Set(ncell);
1003 testindex.Set(ncell+1);
1004 cell[0][0] = ncell + 1;//ajay
1005 //cell[0][0] = ncell;
1006 testncl[0] = cell[0][0];
1007 Int_t ll = 0;
1008 for(Int_t ii = 1; ii <= ncell; ii++)
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
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];
1035 }
1036 delete [] cluster;
1037 delete [] str;
1038 delete [] str1;
1039 delete [] xcl;
1040 delete [] ycl;
1041 delete [] cln;
1042}
1043
1044// ------------------------------------------------------------------------ //
1045Double_t AliPMDClusteringV2::Distance(Double_t x1, Double_t y1,
1046 Double_t x2, Double_t y2)
1047{
1048 return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
1049}
1050// ------------------------------------------------------------------------ //
1051void AliPMDClusteringV2::SetEdepCut(Float_t decut)
1052{
1053 fCutoff = decut;
1054}
1055// ------------------------------------------------------------------------ //