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