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