]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PMD/AliPMDClusteringV1.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
coverity defects fixed
[u/mrichter/AliRoot.git] / PMD / AliPMDClusteringV1.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 : PMDClusteringV1.cxx, Version 00 //
21// //
22// Date : September 26 2002 //
23// //
24// clustering code for alice pmd //
25// //
26//-----------------------------------------------------//
27
28/* --------------------------------------------------------------------
29 Code developed by S. C. Phatak, Institute of Physics,
30 Bhubaneswar 751 005 ( phatak@iopb.res.in ) Given the energy deposited
31 ( or ADC value ) in each cell of supermodule ( pmd or cpv ), the code
32 builds up superclusters and breaks them into clusters. The input is
33 in array edepcell[kNMX] and cluster information is in a
34 TObjarray. Integer clno gives total number of clusters in the
35 supermodule.
36
37 fClusters is the only global ( public ) variables.
38 Others are local ( private ) to the code.
39 At the moment, the data is read for whole detector ( all supermodules
40 and pmd as well as cpv. This will have to be modify later )
41 LAST UPDATE : October 23, 2002
42-----------------------------------------------------------------------*/
43
44#include <Riostream.h>
45#include <TMath.h>
46#include <TNtuple.h>
47#include <TObjArray.h>
48#include "TRandom.h"
49#include <stdio.h>
50
51#include "AliPMDcludata.h"
52#include "AliPMDcluster.h"
53#include "AliPMDClustering.h"
54#include "AliPMDClusteringV1.h"
55#include "AliLog.h"
56
57
58ClassImp(AliPMDClusteringV1)
59
60const Double_t AliPMDClusteringV1::fgkSqroot3by2=0.8660254; // sqrt(3.)/2.
61
62AliPMDClusteringV1::AliPMDClusteringV1():
63 fPMDclucont(new TObjArray()),
64 fCutoff(0.0),
65 fClusParam(0)
66{
67 for(Int_t i = 0; i < kNDIMX; i++)
68 {
69 for(Int_t j = 0; j < kNDIMY; j++)
70 {
71 fCoord[0][i][j] = i+j/2.;
72 fCoord[1][i][j] = fgkSqroot3by2*j;
73 }
74 }
75}
76// ------------------------------------------------------------------------ //
77AliPMDClusteringV1::AliPMDClusteringV1(const AliPMDClusteringV1& pmdclv1):
78 AliPMDClustering(pmdclv1),
79 fPMDclucont(0),
80 fCutoff(0),
81 fClusParam(0)
82{
83 // copy constructor
84 AliError("Copy constructor not allowed ");
85
86}
87// ------------------------------------------------------------------------ //
88AliPMDClusteringV1 &AliPMDClusteringV1::operator=(const AliPMDClusteringV1& /*pmdclv1*/)
89{
90 // copy constructor
91 AliError("Assignment operator not allowed ");
92 return *this;
93}
94// ------------------------------------------------------------------------ //
95AliPMDClusteringV1::~AliPMDClusteringV1()
96{
97 delete fPMDclucont;
98}
99// ------------------------------------------------------------------------ //
100void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn,
101 Int_t celltrack[48][96],
102 Int_t cellpid[48][96],
103 Double_t celladc[48][96],
104 TObjArray *pmdcont)
105{
106 // main function to call other necessary functions to do clustering
107 //
108
109 AliPMDcluster *pmdcl = 0;
110
111 const float ktwobysqrt3 = 1.1547; // 2./sqrt(3.)
112 const Int_t kNmaxCell = 19; // # of cells surrounding a cluster center
113
114 Int_t i = 0, j = 0, nmx1 = 0;
115 Int_t incr = 0, id = 0, jd = 0;
116 Int_t celldataX[kNmaxCell], celldataY[kNmaxCell];
117 Int_t celldataTr[kNmaxCell], celldataPid[kNmaxCell];
118 Float_t celldataAdc[kNmaxCell];
119 Float_t clusdata[6] = {0.,0.,0.,0.,0.,0.};
120 Double_t cutoff, ave;
121 Double_t edepcell[kNMX];
122 Double_t cellenergy[kNMX];
123
124 // ndimXr and ndimYr are different because of different module size
125
126 Int_t ndimXr = 0;
127 Int_t ndimYr = 0;
128
129 if (ismn < 12)
130 {
131 ndimXr = 96;
132 ndimYr = 48;
133 }
134 else if (ismn >= 12 && ismn <= 23)
135 {
136 ndimXr = 48;
137 ndimYr = 96;
138 }
139
140 for (i = 0; i < kNMX; i++)
141 {
142 cellenergy[i] = 0.;
143 }
144
145 Int_t kk = 0;
146 for (i = 0; i < kNDIMX; i++)
147 {
148 for (j = 0; j < kNDIMY; j++)
149 {
150 edepcell[kk] = 0.;
151 kk++;
152 }
153 }
154
155 for (id = 0; id < ndimXr; id++)
156 {
157 for (jd = 0; jd < ndimYr; jd++)
158 {
159 j = jd;
160 i = id+(ndimYr/2-1)-(jd/2);
161
162 Int_t ij = i + j*kNDIMX;
163
164 if (ismn < 12)
165 {
166 cellenergy[ij] = celladc[jd][id];
167 }
168 else if (ismn >= 12 && ismn <= 23)
169 {
170 cellenergy[ij] = celladc[id][jd];
171 }
172 }
173 }
174
175 for (i = 0; i < kNMX; i++)
176 {
177 edepcell[i] = cellenergy[i];
178 }
179
180 Bool_t jsort = true;
181 // the dimension of iord1 is increased twice
182 Int_t iord1[2*kNMX];
183 TMath::Sort((Int_t)kNMX,edepcell,iord1,jsort);// order the data
184 cutoff = fCutoff; // cutoff to discard cells
185 ave = 0.;
186 nmx1 = -1;
187 for(i = 0;i < kNMX; i++)
188 {
189 if(edepcell[i] > 0.)
190 {
191 ave += edepcell[i];
192 }
193 if(edepcell[i] > cutoff )
194 {
195 nmx1++;
196 }
197 }
198
199 AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
200
201 if (nmx1 == 0) nmx1 = 1;
202 ave = ave/nmx1;
203 AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
204 kNMX,ave));
205
206 incr = CrClust(ave, cutoff, nmx1,iord1, edepcell );
207 RefClust(incr,edepcell);
208 Int_t nentries1 = fPMDclucont->GetEntries();
209 AliDebug(1,Form("Detector Plane = %d Serial Module No = %d Number of clusters = %d",idet, ismn, nentries1));
210 AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
211
212 for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
213 {
214 AliPMDcludata *pmdcludata =
215 (AliPMDcludata*)fPMDclucont->UncheckedAt(ient1);
216 Float_t cluXC = pmdcludata->GetClusX();
217 Float_t cluYC = pmdcludata->GetClusY();
218 Float_t cluADC = pmdcludata->GetClusADC();
219 Float_t cluCELLS = pmdcludata->GetClusCells();
220 Float_t cluSIGX = pmdcludata->GetClusSigmaX();
221 Float_t cluSIGY = pmdcludata->GetClusSigmaY();
222
223 Float_t cluY0 = ktwobysqrt3*cluYC;
224 Float_t cluX0 = cluXC - cluY0/2.;
225
226 //
227 // Cluster X centroid is back transformed
228 //
229
230 if (ismn < 12)
231 {
232 clusdata[0] = cluX0 - (24-1) + cluY0/2.;
233 }
234 else if ( ismn >= 12 && ismn <= 23)
235 {
236 clusdata[0] = cluX0 - (48-1) + cluY0/2.;
237 }
238
239 clusdata[1] = cluY0;
240 clusdata[2] = cluADC;
241 clusdata[3] = cluCELLS;
242 clusdata[4] = cluSIGX;
243 clusdata[5] = cluSIGY;
244
245 //
246 // Cells associated with a cluster
247 //
248
249 for (Int_t ihit = 0; ihit < kNmaxCell; ihit++)
250 {
251 Int_t cellrow = pmdcludata->GetCellXY(ihit)/10000;
252 Int_t cellcol = pmdcludata->GetCellXY(ihit)%10000;
253
254 if (ismn < 12)
255 {
256 celldataX[ihit] = cellrow - (24-1) + int(cellcol/2.);
257 }
258 else if (ismn >= 12 && ismn <= 23)
259 {
260 celldataX[ihit] = cellrow - (48-1) + int(cellcol/2.);
261 }
262
263 celldataY[ihit] = cellcol;
264
265 Int_t irow = celldataX[ihit];
266 Int_t icol = celldataY[ihit];
267
268 if (ismn < 12)
269 {
270 if ((irow >= 0 && irow < 96) && (icol >= 0 && icol < 48))
271 {
272 celldataTr[ihit] = celltrack[icol][irow];
273 celldataPid[ihit] = cellpid[icol][irow];
274 celldataAdc[ihit] = (Float_t) celladc[icol][irow];
275 }
276 else
277 {
278 celldataTr[ihit] = -1;
279 celldataPid[ihit] = -1;
280 celldataAdc[ihit] = -1;
281 }
282 }
283 else if (ismn >= 12 && ismn < 24)
284 {
285 if ((irow >= 0 && irow < 48) && (icol >= 0 && icol < 96))
286 {
287 celldataTr[ihit] = celltrack[irow][icol];
288 celldataPid[ihit] = cellpid[irow][icol];
289 celldataAdc[ihit] = (Float_t) celladc[irow][icol];
290
291 }
292 else
293 {
294 celldataTr[ihit] = -1;
295 celldataPid[ihit] = -1;
296 celldataAdc[ihit] = -1;
297 }
298 }
299
300 }
301
302 pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY,
303 celldataTr, celldataPid, celldataAdc);
304 pmdcont->Add(pmdcl);
305 }
306
307 fPMDclucont->Delete();
308}
309// ------------------------------------------------------------------------ //
310Int_t AliPMDClusteringV1::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
311 Int_t iord1[], Double_t edepcell[])
312{
313 // Does crude clustering
314 // Finds out only the big patch by just searching the
315 // connected cells
316 //
317 const Int_t kndim = 4609;
318 Int_t i=0,j=0,k=0,id1=0,id2=0,icl=0, numcell=0;
319 Int_t jd1=0,jd2=0, icell=0, cellcount=0;
320 Int_t clust[2][kndim];
321 static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
322
323 AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
324
325 for (j = 0; j < kNDIMX; j++)
326 {
327 for(k = 0; k < kNDIMY; k++)
328 {
329 fInfocl[0][j][k] = 0;
330 fInfocl[1][j][k] = 0;
331 }
332 }
333 for(i=0; i < kNMX; i++)
334 {
335 fInfcl[0][i] = -1;
336
337 j = iord1[i];
338 id2 = j/kNDIMX;
339 id1 = j-id2*kNDIMX;
340
341 if(edepcell[j] <= cutoff)
342 {
343 fInfocl[0][id1][id2] = -1;
344 }
345 }
346
347 // ---------------------------------------------------------------
348 // crude clustering begins. Start with cell having largest adc
349 // count and loop over the cells in descending order of adc count
350 // ---------------------------------------------------------------
351
352 icl = -1;
353 cellcount = -1;
354
355 for(icell = 0; icell <= nmx1; icell++)
356 {
357 j = iord1[icell];
358 id2 = j/kNDIMX;
359 id1 = j-id2*kNDIMX;
360
361 if(fInfocl[0][id1][id2] == 0 )
362 {
363 icl++;
364 numcell = 0;
365 cellcount++;
366 fInfocl[0][id1][id2] = 1;
367 fInfocl[1][id1][id2] = icl;
368 fInfcl[0][cellcount] = icl;
369 fInfcl[1][cellcount] = id1;
370 fInfcl[2][cellcount] = id2;
371
372 clust[0][numcell] = id1;
373 clust[1][numcell] = id2;
374
375 for(i = 1; i < kndim; i++)
376 {
377 clust[0][i]=0;
378 }
379 // ---------------------------------------------------------------
380 // check for adc count in neib. cells. If ne 0 put it in this clust
381 // ---------------------------------------------------------------
382 for(i = 0; i < 6; i++)
383 {
384 jd1 = id1 + neibx[i];
385 jd2 = id2 + neiby[i];
386 if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
387 fInfocl[0][jd1][jd2] == 0)
388 {
389 numcell++;
390 fInfocl[0][jd1][jd2] = 2;
391 fInfocl[1][jd1][jd2] = icl;
392 clust[0][numcell] = jd1;
393 clust[1][numcell] = jd2;
394 cellcount++;
395 fInfcl[0][cellcount] = icl;
396 fInfcl[1][cellcount] = jd1;
397 fInfcl[2][cellcount] = jd2;
398 }
399 }
400 // ---------------------------------------------------------------
401 // check adc count for neighbour's neighbours recursively and
402 // if nonzero, add these to the cluster.
403 // ---------------------------------------------------------------
404 for(i = 1; i < kndim;i++)
405 {
406 if(clust[0][i] != 0)
407 {
408 id1 = clust[0][i];
409 id2 = clust[1][i];
410 for(j = 0; j < 6 ; j++)
411 {
412 jd1 = id1 + neibx[j];
413 jd2 = id2 + neiby[j];
414 if( (jd1 >= 0 && jd1 < kNDIMX) &&
415 (jd2 >= 0 && jd2 < kNDIMY) &&
416 fInfocl[0][jd1][jd2] == 0 )
417 {
418 fInfocl[0][jd1][jd2] = 2;
419 fInfocl[1][jd1][jd2] = icl;
420 numcell++;
421 clust[0][numcell] = jd1;
422 clust[1][numcell] = jd2;
423 cellcount++;
424 fInfcl[0][cellcount] = icl;
425 fInfcl[1][cellcount] = jd1;
426 fInfcl[2][cellcount] = jd2;
427 }
428 }
429 }
430 }
431 }
432 }
433 return cellcount;
434}
435// ------------------------------------------------------------------------ //
436void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
437{
438 // Does the refining of clusters
439 // Takes the big patch and does gaussian fitting and
440 // finds out the more refined clusters
441 //
442
443 AliPMDcludata *pmdcludata = 0;
444
445 const Int_t kNmaxCell = 19; // # of cells surrounding a cluster center
446
447 Int_t ndim = incr + 1;
448
449 Int_t *ncl = 0x0;
450 Int_t *clxy = 0x0;
451 Int_t i12 = 0, i22 = 0;
452 Int_t i = 0, j = 0, k = 0;
453 Int_t i1 = 0, i2 = 0, id = 0, icl = 0;
454 Int_t itest = 0, ihld = 0, ig = 0;
455 Int_t nsupcl = 0, clno = 0, t1 = 0, t2 = 0;
456 Float_t clusdata[6];
457 Double_t x1 = 0, y1 = 0, z1 = 0, x2 = 0, y2 = 0, z2 = 0, rr = 0;
458
459 ncl = new Int_t [ndim];
460 clxy = new Int_t [kNmaxCell];
461
462 // Initialisation
463 for(i = 0; i<ndim; i++)
464 {
465 ncl[i] = -1;
466 }
467 for(i = 0; i<6; i++)
468 {
469 clusdata[i] = 0.;
470 }
471 for(i = 0; i<19; i++)
472 {
473 clxy[i] = 0;
474 }
475
476 // clno counts the final clusters
477 // nsupcl = # of superclusters; ncl[i]= # of cells in supercluster i
478 // x, y and z store (x,y) coordinates of and energy deposited in a cell
479 // xc, yc store (x,y) coordinates of the cluster center
480 // zc stores the energy deposited in a cluster
481 // rc is cluster radius
482
483 clno = -1;
484 nsupcl = -1;
485
486 for(i = 0; i <= incr; i++)
487 {
488 if(fInfcl[0][i] != nsupcl)
489 {
490 nsupcl++;
491 }
492 if (nsupcl > ndim)
493 {
494 AliWarning("RefClust: Too many superclusters!");
495 nsupcl = ndim;
496 break;
497 }
498 ncl[nsupcl]++;
499 }
500
501 AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
502 incr+1,nsupcl+1));
503 id = -1;
504 icl = -1;
505
506 for(i = 0; i <= nsupcl; i++)
507 {
508 if(ncl[i] == 0)
509 {
510 id++;
511 icl++;
512 if (clno >= 4608)
513 {
514 AliWarning("RefClust: Too many clusters! more than 4608");
515 delete [] ncl;
516 delete [] clxy;
517 return;
518 }
519 clno++;
520 i1 = fInfcl[1][id];
521 i2 = fInfcl[2][id];
522
523 i12 = i1 + i2*kNDIMX;
524
525 clusdata[0] = fCoord[0][i1][i2];
526 clusdata[1] = fCoord[1][i1][i2];
527 clusdata[2] = edepcell[i12];
528 clusdata[3] = 1.;
529 clusdata[4] = 999.5;
530 clusdata[5] = 999.5;
531
532 clxy[0] = i1*10000 + i2;
533
534 for(Int_t icltr = 1; icltr < kNmaxCell; icltr++)
535 {
536 clxy[icltr] = -1;
537 }
538
539 pmdcludata = new AliPMDcludata(clusdata,clxy);
540 fPMDclucont->Add(pmdcludata);
541 }
542 else if(ncl[i] == 1)
543 {
544 id++;
545 icl++;
546 if (clno >= 4608)
547 {
548 AliWarning("RefClust: Too many clusters! more than 4608");
549 delete [] ncl;
550 delete [] clxy;
551
552 return;
553 }
554 clno++;
555 i1 = fInfcl[1][id];
556 i2 = fInfcl[2][id];
557 i12 = i1 + i2*kNDIMX;
558
559 x1 = fCoord[0][i1][i2];
560 y1 = fCoord[1][i1][i2];
561 z1 = edepcell[i12];
562
563 clxy[0] = i1*10000 + i2;
564
565 id++;
566 i1 = fInfcl[1][id];
567 i2 = fInfcl[2][id];
568
569 i22 = i1 + i2*kNDIMX;
570 x2 = fCoord[0][i1][i2];
571 y2 = fCoord[1][i1][i2];
572 z2 = edepcell[i22];
573
574 clxy[1] = i1*10000 + i2;
575
576
577 for(Int_t icltr = 2; icltr < kNmaxCell; icltr++)
578 {
579 clxy[icltr] = -1;
580 }
581
582 clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
583 clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
584 clusdata[2] = z1+z2;
585 clusdata[3] = 2.;
586 clusdata[4] = 0.5;
587 clusdata[5] = 0.0;
588 pmdcludata = new AliPMDcludata(clusdata,clxy);
589 fPMDclucont->Add(pmdcludata);
590 }
591 else
592 {
593 Int_t *iord, *tc, *t;
594 Double_t *x, *y, *z, *xc, *yc, *zc;
595
596 iord = new Int_t [ncl[i]+1];
597 tc = new Int_t [ncl[i]+1];
598 t = new Int_t [ncl[i]+1];
599 x = new Double_t [ncl[i]+1];
600 y = new Double_t [ncl[i]+1];
601 z = new Double_t [ncl[i]+1];
602 xc = new Double_t [ncl[i]+1];
603 yc = new Double_t [ncl[i]+1];
604 zc = new Double_t [ncl[i]+1];
605
606 for( k = 0; k < ncl[i]+1; k++)
607 {
608 iord[k] = -1;
609 t[k] = -1;
610 tc[k] = -1;
611 x[k] = -1;
612 y[k] = -1;
613 z[k] = -1;
614 xc[k] = -1;
615 yc[k] = -1;
616 zc[k] = -1;
617 }
618 id++;
619 // super-cluster of more than two cells - broken up into smaller
620 // clusters gaussian centers computed. (peaks separated by > 1 cell)
621 // Begin from cell having largest energy deposited This is first
622 // cluster center
623 i1 = fInfcl[1][id];
624 i2 = fInfcl[2][id];
625 i12 = i1 + i2*kNDIMX;
626
627 x[0] = fCoord[0][i1][i2];
628 y[0] = fCoord[1][i1][i2];
629 z[0] = edepcell[i12];
630 t[0] = i1*10000 + i2;
631
632
633 iord[0] = 0;
634 for(j = 1; j <= ncl[i]; j++)
635 {
636 id++;
637 i1 = fInfcl[1][id];
638 i2 = fInfcl[2][id];
639 i12 = i1 + i2*kNDIMX;
640 iord[j] = j;
641 x[j] = fCoord[0][i1][i2];
642 y[j] = fCoord[1][i1][i2];
643 z[j] = edepcell[i12];
644 t[j] = i1*10000 + i2;
645 }
646
647 // arranging cells within supercluster in decreasing order
648
649 for(j = 1;j <= ncl[i]; j++)
650 {
651 itest = 0;
652 ihld = iord[j];
653 for(i1 = 0; i1 < j; i1++)
654 {
655 if(itest == 0 && z[iord[i1]] < z[ihld])
656 {
657 itest = 1;
658 for(i2 = j-1; i2 >= i1; i2--)
659 {
660 iord[i2+1] = iord[i2];
661 }
662 iord[i1] = ihld;
663 }
664 }
665 }
666
667 if (fClusParam == 1)
668 {
669 // Clustering algorithm returns from here for PP collisions
670 // for pp, only the output of crude clusterng is taken
671 // sigx and sigy are not calculated at this moment
672
673 Double_t supx=0., supy=0., supz=0.;
674
675 for(j = 0;j <= ncl[i]; j++)
676 {
677 supx += x[iord[j]]*z[iord[j]];
678 supy += y[iord[j]]*z[iord[j]];
679 supz += z[iord[j]];
680 if(j < 19)
681 {
682 clxy[j] = t[iord[j]];
683 }
684 }
685
686 if( ncl[i] + 1 < 19)
687 {
688 for(Int_t ncel = ncl[i] + 1; ncel < kNmaxCell; ncel ++ )
689 {
690 clxy[ncel] = -1;
691 }
692 }
693 clusdata[0] = supx/supz;
694 clusdata[1] = supy/supz;
695 clusdata[2] = supz;
696 clusdata[3] = ncl[i]+1;
697 clusdata[4] = 0.5;
698 clusdata[5] = 0.0;
699 pmdcludata = new AliPMDcludata(clusdata,clxy);
700 fPMDclucont->Add(pmdcludata);
701 }
702
703 /* MODIFICATION PART STARTS (Tapan July 2008)
704 iord[0] is the cell with highest ADC in the crude-cluster
705 ig is the number of local maxima in the crude-cluster
706 For the higest peak we make ig=0 which means first local maximum.
707 Next we go down in terms of the ADC sequence and find out if any
708 more of the cells form local maxima. The definition of local
709 maxima is that all its neighbours are of less ADC compared to it.
710 */
711
712 if (fClusParam == 2)
713 {
714 // This part is to split the supercluster
715 //
716 ig = 0;
717 xc[ig] = x[iord[0]];
718 yc[ig] = y[iord[0]];
719 zc[ig] = z[iord[0]];
720 tc[ig] = t[iord[0]];
721 Int_t ivalid = 0, icount = 0;
722
723 for(j=1;j<=ncl[i];j++)
724 {
725 x1 = x[iord[j]];
726 y1 = y[iord[j]];
727 z1 = z[iord[j]];
728 t1 = t[iord[j]];
729 rr=Distance(x1,y1,xc[ig],yc[ig]);
730
731 // Check the cells which are outside the neighbours (rr>1.2)
732 if(rr>1.2 )
733 {
734 ivalid=0;
735 icount=0;
736 for(Int_t j1=1;j1<j;j1++)
737 {
738 icount++;
739 Float_t rr1=Distance(x1,y1,x[iord[j1]],y[iord[j1]]);
740 if(rr1>1.2) ivalid++;
741 }
742 if(ivalid == icount && z1>0.5*zc[ig])
743 {
744 ig++;
745 xc[ig]=x1;
746 yc[ig]=y1;
747 zc[ig]=z1;
748 tc[ig]=t1;
749 }
750 }
751 }
752
753 icl=icl+ig+1;
754
755 // We use simple Gaussian weighting. (Tapan Jan 2005)
756 // compute the number of cells belonging to each cluster.
757 // cell can be shared between several clusters
758 // in the ratio of cluster energy deposition
759 // To calculate:
760 // (1) number of cells belonging to a cluster (ig) and
761 // (2) total ADC of the cluster (ig)
762 // (3) x and y positions of the cluster
763
764
765 Int_t *cellCount;
766 Int_t **cellXY;
767
768 Int_t *status;
769 Double_t *totaladc, *totaladc2, *ncell,*weight;
770 Double_t *xclust, *yclust, *sigxclust, *sigyclust;
771 Double_t *ax, *ay, *ax2, *ay2;
772
773
774 status = new Int_t [ncl[i]+1];
775 cellXY = new Int_t *[ncl[i]+1];
776
777 cellCount = new Int_t [ig+1];
778 totaladc = new Double_t [ig+1];
779 totaladc2 = new Double_t [ig+1];
780 ncell = new Double_t [ig+1];
781 weight = new Double_t [ig+1];
782 xclust = new Double_t [ig+1];
783 yclust = new Double_t [ig+1];
784 sigxclust = new Double_t [ig+1];
785 sigyclust = new Double_t [ig+1];
786 ax = new Double_t [ig+1];
787 ay = new Double_t [ig+1];
788 ax2 = new Double_t [ig+1];
789 ay2 = new Double_t [ig+1];
790
791 for(j = 0; j < ncl[i]+1; j++)
792 {
793 status[j] = 0;
794 cellXY[j] = new Int_t[ig+1];
795 }
796 //initialization
797 for(Int_t kcl = 0; kcl < ig+1; kcl++)
798 {
799 cellCount[kcl] = 0;
800 totaladc[kcl] = 0.;
801 totaladc2[kcl] = 0.;
802 ncell[kcl] = 0.;
803 weight[kcl] = 0.;
804 xclust[kcl] = 0.;
805 yclust[kcl] = 0.;
806 sigxclust[kcl] = 0.;
807 sigyclust[kcl] = 0.;
808 ax[kcl] = 0.;
809 ay[kcl] = 0.;
810 ax2[kcl] = 0.;
811 ay2[kcl] = 0.;
812 for(j = 0; j < ncl[i]+1; j++)
813 {
814 cellXY[j][kcl] = 0;
815 }
816 }
817 Double_t sumweight, gweight;
818
819 for(j = 0;j <= ncl[i]; j++)
820 {
821 x1 = x[iord[j]];
822 y1 = y[iord[j]];
823 z1 = z[iord[j]];
824 t1 = t[iord[j]];
825
826 for(Int_t kcl=0; kcl<=ig; kcl++)
827 {
828 x2 = xc[kcl];
829 y2 = yc[kcl];
830 rr = Distance(x1,y1,x2,y2);
831 t2 = tc[kcl];
832
833 if(rr==0)
834 {
835 ncell[kcl] = 1.;
836 totaladc[kcl] = z1;
837 totaladc2[kcl] = z1*z1;
838 ax[kcl] = x1 * z1;
839 ay[kcl] = y1 * z1;
840 ax2[kcl] = 0.;
841 ay2[kcl] = 0.;
842 status[j] = 1;
843 }
844 }
845 }
846
847 for(j = 0; j <= ncl[i]; j++)
848 {
849 Int_t maxweight = 0;
850 Double_t max = 0.;
851
852 if(status[j] == 0)
853 {
854 x1 = x[iord[j]];
855 y1 = y[iord[j]];
856 z1 = z[iord[j]];
857 t1 = t[iord[j]];
858 sumweight = 0.;
859
860 for(Int_t kcl = 0; kcl <= ig; kcl++)
861 {
862 x2 = xc[kcl];
863 y2 = yc[kcl];
864 rr = Distance(x1,y1,x2,y2);
865 gweight = exp(-(rr*rr)/(2*(1.2*1.2)));
866 weight[kcl] = zc[kcl] * gweight;
867 sumweight = sumweight + weight[kcl];
868
869 if(weight[kcl] > max)
870 {
871 max = weight[kcl];
872 maxweight = kcl;
873 }
874 }
875
876 cellXY[cellCount[maxweight]][maxweight] = iord[j];
877
878 cellCount[maxweight]++;
879
880 x2 = xc[maxweight];
881 y2 = yc[maxweight];
882 totaladc[maxweight] += z1;
883 ax[maxweight] += x1*z1;
884 ay[maxweight] += y1*z1;
885 totaladc2[maxweight] += z1*z1;
886 ax2[maxweight] += z1*(x1-x2)*(x1-x2);
887 ay2[maxweight] += z1*(y1-y2)*(y1-y2);
888 ncell[maxweight]++;
889
890 }
891 }
892
893 for(Int_t kcl = 0; kcl <= ig; kcl++)
894 {
895 if(totaladc[kcl] > 0.)
896 {
897 xclust[kcl] = (ax[kcl])/ totaladc[kcl];
898 yclust[kcl] = (ay[kcl])/ totaladc[kcl];
899
900 //natasha
901 Float_t sqtotadc = totaladc[kcl]*totaladc[kcl];
902 if(totaladc2[kcl] >= sqtotadc)
903 {
904 sigxclust[kcl] = 0.25;
905 sigyclust[kcl] = 0.25;
906 }
907 else
908 {
909 sigxclust[kcl] = (totaladc[kcl]/(sqtotadc-totaladc2[kcl]))*ax2[kcl];
910 sigyclust[kcl] = (totaladc[kcl]/(sqtotadc-totaladc2[kcl]))*ay2[kcl];
911 }
912 }
913
914 for(j = 0; j < cellCount[kcl]; j++) clno++;
915
916 if (clno >= 4608)
917 {
918 AliWarning("RefClust: Too many clusters! more than 4608");
919
920 delete [] cellCount;
921 for(Int_t jj = 0; jj < ncl[i]+1; jj++) delete [] cellXY[jj];
922
923 delete [] status;
924 delete [] totaladc;
925 delete [] totaladc2;
926 delete [] ncell;
927 delete [] xclust;
928 delete [] yclust;
929 delete [] sigxclust;
930 delete [] sigyclust;
931 delete [] ax;
932 delete [] ay;
933 delete [] ax2;
934 delete [] ay2;
935 delete [] weight;
936
937 delete [] iord;
938 delete [] tc;
939 delete [] t;
940 delete [] x;
941 delete [] y;
942 delete [] z;
943 delete [] xc;
944 delete [] yc;
945 delete [] zc;
946
947
948 delete [] ncl;
949 delete [] clxy;
950
951 return;
952 }
953 clusdata[0] = xclust[kcl];
954 clusdata[1] = yclust[kcl];
955 clusdata[2] = totaladc[kcl];
956 clusdata[3] = ncell[kcl];
957
958 if(sigxclust[kcl] > sigyclust[kcl])
959 {
960 clusdata[4] = TMath::Sqrt(sigxclust[kcl]);
961 clusdata[5] = TMath::Sqrt(sigyclust[kcl]);
962 }
963 else
964 {
965 clusdata[4] = TMath::Sqrt(sigyclust[kcl]);
966 clusdata[5] = TMath::Sqrt(sigxclust[kcl]);
967 }
968
969 clxy[0] = tc[kcl];
970
971 Int_t Ncell=1;
972 for (Int_t ii = 0; ii < cellCount[kcl]; ii++)
973 {
974 if(ii<18)
975 {
976 clxy[Ncell] = t[cellXY[ii][kcl]];
977 Ncell++;
978 }
979 }
980
981 pmdcludata = new AliPMDcludata(clusdata,clxy);
982 fPMDclucont->Add(pmdcludata);
983 }
984 delete [] cellCount;
985 for(Int_t jj = 0; jj < ncl[i]+1; jj++) delete [] cellXY[jj];
986
987 delete [] status;
988 delete [] totaladc;
989 delete [] totaladc2;
990 delete [] ncell;
991 delete [] xclust;
992 delete [] yclust;
993 delete [] sigxclust;
994 delete [] sigyclust;
995 delete [] ax;
996 delete [] ay;
997 delete [] ax2;
998 delete [] ay2;
999 delete [] weight;
1000
1001 }
1002
1003 delete [] iord;
1004 delete [] tc;
1005 delete [] t;
1006 delete [] x;
1007 delete [] y;
1008 delete [] z;
1009 delete [] xc;
1010 delete [] yc;
1011 delete [] zc;
1012
1013
1014 }
1015 }
1016 delete [] ncl;
1017 delete [] clxy;
1018}
1019// ------------------------------------------------------------------------ //
1020Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1,
1021 Double_t x2, Double_t y2)
1022{
1023 return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
1024}
1025// ------------------------------------------------------------------------ //
1026void AliPMDClusteringV1::SetEdepCut(Float_t decut)
1027{
1028 fCutoff = decut;
1029}
1030// ------------------------------------------------------------------------ //
1031void AliPMDClusteringV1::SetClusteringParam(Int_t cluspar)
1032{
1033 fClusParam = cluspar;
1034}
1035// ------------------------------------------------------------------------ //