new definitions
[u/mrichter/AliRoot.git] / PMD / AliPMDClusteringV1.cxx
CommitLineData
3edbbba2 1/***************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
090026bf 16/* $Id$ */
17
3edbbba2 18//-----------------------------------------------------//
19// //
01c4d84a 20// Source File : PMDClusteringV1.cxx, Version 00 //
3edbbba2 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
fd30f88e 33 in array fEdepCell[kNDIMX][kNDIMY] and cluster information is in a
34 TObjarray. Integer clno gives total number of clusters in the
3edbbba2 35 supermodule.
36
fd30f88e 37 fEdepCell and fClusters are the only global ( public ) variables.
3edbbba2 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
090026bf 44#include <Riostream.h>
45#include <TMath.h>
3edbbba2 46#include <TNtuple.h>
47#include <TObjArray.h>
48#include <stdio.h>
49
fd30f88e 50#include "AliPMDcludata.h"
3edbbba2 51#include "AliPMDcluster.h"
52#include "AliPMDClustering.h"
53#include "AliPMDClusteringV1.h"
54#include "AliLog.h"
fd30f88e 55#include "TRandom.h"
3edbbba2 56
57ClassImp(AliPMDClusteringV1)
58
59const Double_t AliPMDClusteringV1::fgkSqroot3by2=0.8660254; // sqrt(3.)/2.
60
61AliPMDClusteringV1::AliPMDClusteringV1():
fd30f88e 62 pmdclucont(new TObjArray()),
3edbbba2 63 fCutoff(0.0)
64{
fd30f88e 65 for(Int_t i = 0; i < kNDIMX; i++)
3edbbba2 66 {
fd30f88e 67 for(Int_t j = 0; j < kNDIMY; j++)
3edbbba2 68 {
69 fCoord[0][i][j] = i+j/2.;
70 fCoord[1][i][j] = fgkSqroot3by2*j;
71 fEdepCell[i][j] = 0;
72 }
73 }
74}
75// ------------------------------------------------------------------------ //
76AliPMDClusteringV1::~AliPMDClusteringV1()
77{
fd30f88e 78 delete pmdclucont;
3edbbba2 79}
80// ------------------------------------------------------------------------ //
fd30f88e 81void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn,
82 Double_t celladc[48][96], TObjArray *pmdcont)
3edbbba2 83{
84 // main function to call other necessary functions to do clustering
85 //
86 AliPMDcluster *pmdcl = 0;
87 /*
88 int id and jd defined to read the input data.
89 It is assumed that for data we have 0 <= id <= 48
90 and 0 <= jd <=96
91 */
01c4d84a 92
fd30f88e 93 Int_t i, i1, i2, j, nmx1, incr, id, jd;
3edbbba2 94 Int_t celldataX[15], celldataY[15];
0709b99f 95 Float_t clusdata[6];
3edbbba2 96
97 Double_t cutoff, ave;
98
99 const float ktwobysqrt3 = 1.1547; // 2./sqrt(3.)
100
01c4d84a 101 // ndimXr and ndimYr are different because of different module size
102
fd30f88e 103 Int_t ndimXr = 0;
104 Int_t ndimYr = 0;
01c4d84a 105
106 if (ismn < 12)
107 {
108 ndimXr = 96;
109 ndimYr = 48;
110 }
111 else if (ismn >= 12 && ismn <= 23)
112 {
113 ndimXr = 48;
114 ndimYr = 96;
115 }
116
fd30f88e 117 for (Int_t i = 0; i < kNDIMX; i++)
3edbbba2 118 {
fd30f88e 119 for (Int_t j = 0; j < kNDIMY; j++)
01c4d84a 120 {
121 fEdepCell[i][j] = 0;
122 fCellTrNo[i][j] = -1;
123 }
124 }
125
126 for (id = 0; id < ndimXr; id++)
127 {
128 for (jd = 0; jd < ndimYr; jd++)
3edbbba2 129 {
fd30f88e 130 j = jd;
131 i = id+(ndimYr/2-1)-(jd/2);
01c4d84a 132
133 if (ismn < 12)
134 {
135 fEdepCell[i][j] = celladc[jd][id];
136 fCellTrNo[i][j] = jd*10000+id; /* for association */
137 }
138 else if (ismn >= 12 && ismn <= 23)
139 {
140 fEdepCell[i][j] = celladc[id][jd];
141 fCellTrNo[i][j] = id*10000+jd; /* for association */
142 }
143
3edbbba2 144 }
145 }
fd30f88e 146 Order(); // order the data
3edbbba2 147 cutoff = fCutoff; // cutoff used to discard cells having ener. dep.
fd30f88e 148 ave = 0.;
149 nmx1 = -1;
3edbbba2 150
fd30f88e 151 for(j = 0;j < kNMX; j++)
3edbbba2 152 {
153 i1 = fIord[0][j];
154 i2 = fIord[1][j];
fd30f88e 155 if(fEdepCell[i1][i2] > 0.)
156 {
157 ave += fEdepCell[i1][i2];
158 }
159 if(fEdepCell[i1][i2] > cutoff )
160 {
161 nmx1++;
162 }
3edbbba2 163 }
3edbbba2 164
165 AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
166
3edbbba2 167 if (nmx1 == 0) nmx1 = 1;
fd30f88e 168 ave = ave/nmx1;
3edbbba2 169
170 AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
171 kNMX,ave));
172
173 incr = CrClust(ave, cutoff, nmx1);
174 RefClust(incr);
fd30f88e 175 Int_t nentries1 = pmdclucont->GetEntries();
176 AliDebug(1,Form("Detector Plane = %d Serial Module No = %d Number of clusters = %d",idet, ismn, nentries1));
177 AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
178 for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
3edbbba2 179 {
fd30f88e 180 AliPMDcludata *pmdcludata =
181 (AliPMDcludata*)pmdclucont->UncheckedAt(ient1);
182 Float_t cluXC = pmdcludata->GetClusX();
183 Float_t cluYC = pmdcludata->GetClusY();
184 Float_t cluADC = pmdcludata->GetClusADC();
185 Float_t cluCELLS = pmdcludata->GetClusCells();
186 Float_t cluSIGX = pmdcludata->GetClusSigmaX();
187 Float_t cluSIGY = pmdcludata->GetClusSigmaY();
188
3edbbba2 189 Float_t cluY0 = ktwobysqrt3*cluYC;
190 Float_t cluX0 = cluXC - cluY0/2.;
fd30f88e 191
3edbbba2 192 //
193 // Cluster X centroid is back transformed
194 //
01c4d84a 195 if (ismn < 12)
196 {
197 clusdata[0] = cluX0 - (24-1) + cluY0/2.;
198 }
fd30f88e 199 else if ( ismn >= 12 && ismn <= 23)
01c4d84a 200 {
201 clusdata[0] = cluX0 - (48-1) + cluY0/2.;
202 }
fd30f88e 203
01c4d84a 204 clusdata[1] = cluY0;
205 clusdata[2] = cluADC;
206 clusdata[3] = cluCELLS;
fd30f88e 207 clusdata[4] = cluSIGX;
208 clusdata[5] = cluSIGY;
209
3edbbba2 210 //
211 // Cells associated with a cluster
212 //
213 for (Int_t ihit = 0; ihit < 15; ihit++)
214 {
fd30f88e 215
01c4d84a 216 if (ismn < 12)
217 {
218 celldataX[ihit] = fClTr[ihit][i1]%10000;
219 celldataY[ihit] = fClTr[ihit][i1]/10000;
220 }
221 else if (ismn >= 12 && ismn <= 23)
222 {
223 celldataX[ihit] = fClTr[ihit][i1]/10000;
224 celldataY[ihit] = fClTr[ihit][i1]%10000;
225 }
3edbbba2 226 }
3edbbba2 227 pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY);
228 pmdcont->Add(pmdcl);
229 }
fd30f88e 230
231 pmdclucont->Clear();
232
3edbbba2 233}
234// ------------------------------------------------------------------------ //
235void AliPMDClusteringV1::Order()
236{
237 // Sorting algorithm
238 // sorts the ADC values from higher to lower
239 //
fd30f88e 240 Int_t i, j, i1, i2;
241 Int_t iord1[kNMX];
242 Double_t dd[kNMX];
243
3edbbba2 244 for(i1=0; i1 < kNDIMX; i1++)
245 {
246 for(i2=0; i2 < kNDIMY; i2++)
247 {
248 i = i1 + i2*kNDIMX;
249 iord1[i] = i;
250 dd[i] = fEdepCell[i1][i2];
251 }
252 }
253
254 TMath::Sort(kNMX,dd,iord1); //PH Using much better algorithm...
3edbbba2 255 for(i=0; i<kNMX; i++)
256 {
257 j = iord1[i];
258 i2 = j/kNDIMX;
259 i1 = j-i2*kNDIMX;
260 fIord[0][i]=i1;
261 fIord[1][i]=i2;
262 }
263}
264// ------------------------------------------------------------------------ //
fd30f88e 265Int_t AliPMDClusteringV1::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1)
3edbbba2 266{
fd30f88e 267 // Does crude clustering
3edbbba2 268 // Finds out only the big patch by just searching the
269 // connected cells
270 //
fd30f88e 271 Int_t i,j,k,id1,id2,icl, numcell, clust[2][5000];
272 Int_t jd1,jd2, icell, cellcount;
273 static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
274
3edbbba2 275 // neibx and neiby define ( incremental ) (i,j) for the neighbours of a
276 // cell. There are six neighbours.
277 // cellcount --- total number of cells having nonzero ener dep
278 // numcell --- number of cells in a given supercluster
fd30f88e 279
3edbbba2 280 AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
281
fd30f88e 282 for (j = 0; j < kNDIMX; j++)
283 {
284 for(k = 0; k < kNDIMY; k++)
285 {
286 fInfocl[0][j][k] = 0;
287 fInfocl[1][j][k] = 0;
288 }
289 }
290 for(i=0; i < kNMX; i++)
291 {
292 fInfcl[0][i] = -1;
293 id1 = fIord[0][i];
294 id2 = fIord[1][i];
295 if(fEdepCell[id1][id2] <= cutoff)
296 {
297 fInfocl[0][id1][id2] = -1;
298 }
3edbbba2 299 }
3edbbba2 300 // ---------------------------------------------------------------
301 // crude clustering begins. Start with cell having largest adc
302 // count and loop over the cells in descending order of adc count
303 // ---------------------------------------------------------------
fd30f88e 304 icl = -1;
305 cellcount = -1;
3edbbba2 306
fd30f88e 307 for(icell = 0; icell <= nmx1; icell++)
308 {
309 id1 = fIord[0][icell];
310 id2 = fIord[1][icell];
311 if(fInfocl[0][id1][id2] == 0 )
312 {
313 icl++;
314 numcell = 0;
315 cellcount++;
316 fInfocl[0][id1][id2] = 1;
317 fInfocl[1][id1][id2] = icl;
318 fInfcl[0][cellcount] = icl;
319 fInfcl[1][cellcount] = id1;
320 fInfcl[2][cellcount] = id2;
321
322 clust[0][numcell] = id1;
323 clust[1][numcell] = id2;
324
325 for(i = 1; i < 5000; i++)
326 {
327 clust[0][i]=0;
328 }
329 // ---------------------------------------------------------------
330 // check for adc count in neib. cells. If ne 0 put it in this clust
331 // ---------------------------------------------------------------
332 for(i = 0; i < 6; i++)
333 {
334 jd1 = id1 + neibx[i];
335 jd2 = id2 + neiby[i];
336 if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
337 fInfocl[0][jd1][jd2] == 0)
338 {
339 numcell++;
340 fInfocl[0][jd1][jd2] = 2;
341 fInfocl[1][jd1][jd2] = icl;
342 clust[0][numcell] = jd1;
343 clust[1][numcell] = jd2;
344 cellcount++;
345 fInfcl[0][cellcount] = icl;
346 fInfcl[1][cellcount] = jd1;
347 fInfcl[2][cellcount] = jd2;
348 }
349 }
350 // ---------------------------------------------------------------
351 // check adc count for neighbour's neighbours recursively and
352 // if nonzero, add these to the cluster.
353 // ---------------------------------------------------------------
354 for(i = 1; i < 5000;i++)
355 {
356 if(clust[0][i] != 0)
357 {
358 id1 = clust[0][i];
359 id2 = clust[1][i];
360 for(j = 0; j < 6 ; j++)
361 {
362 jd1 = id1 + neibx[j];
363 jd2 = id2 + neiby[j];
364 if( (jd1 >= 0 && jd1 < kNDIMX) &&
365 (jd2 >= 0 && jd2 < kNDIMY) &&
366 fInfocl[0][jd1][jd2] == 0 )
367 {
368 fInfocl[0][jd1][jd2] = 2;
369 fInfocl[1][jd1][jd2] = icl;
370 numcell++;
371 clust[0][numcell] = jd1;
372 clust[1][numcell] = jd2;
373 cellcount++;
374 fInfcl[0][cellcount] = icl;
375 fInfcl[1][cellcount] = jd1;
376 fInfcl[2][cellcount] = jd2;
377 }
378 }
379 }
3edbbba2 380 }
3edbbba2 381 }
3edbbba2 382 }
3edbbba2 383 return cellcount;
384}
385// ------------------------------------------------------------------------ //
fd30f88e 386void AliPMDClusteringV1::RefClust(Int_t incr)
3edbbba2 387{
388 // Does the refining of clusters
389 // Takes the big patch and does gaussian fitting and
390 // finds out the more refined clusters
391 //
0709b99f 392
fd30f88e 393 const Int_t kdim = 4500;
394 Int_t i, j, k, i1, i2, id, icl, itest,ihld, ig, nsupcl,clno;
395 Double_t x1, y1, z1, x2, y2, z2, dist,rr,sum;
0709b99f 396
fd30f88e 397 Int_t t[kdim],cellCount[kdim];
398 Int_t ncl[kdim], iord[kdim], lev1[20], lev2[20];
399 Double_t x[kdim], y[kdim], z[kdim];
400 Double_t xc[kdim], yc[kdim], zc[kdim], cells[kdim], rc[kdim];
401
402 Float_t clusdata[6];
403 for(Int_t kk = 0; kk < 6; kk++)
01c4d84a 404 {
fd30f88e 405 clusdata[kk] = 0.;
01c4d84a 406 }
0709b99f 407
fd30f88e 408 //asso
409 for(i = 0; i<kdim; i++)
410 {
411 t[i] = -1;
412 cellCount[i] = 0;
413 }
414 // clno counts the final clusters
3edbbba2 415 // nsupcl = # of superclusters; ncl[i]= # of cells in supercluster i
416 // x, y and z store (x,y) coordinates of and energy deposited in a cell
417 // xc, yc store (x,y) coordinates of the cluster center
418 // zc stores the energy deposited in a cluster
419 // rc is cluster radius
fd30f88e 420
421 clno = -1;
3edbbba2 422 nsupcl = -1;
fd30f88e 423 for(i = 0; i < kdim; i++)
424 {
425 ncl[i] = -1;
3edbbba2 426 }
fd30f88e 427 for(i = 0; i <= incr; i++)
428 {
429 if(fInfcl[0][i] != nsupcl)
01c4d84a 430 {
fd30f88e 431 nsupcl++;
01c4d84a 432 }
fd30f88e 433 if (nsupcl > kdim)
01c4d84a 434 {
fd30f88e 435 AliWarning("RefClust: Too many superclusters!");
436 nsupcl = kdim;
437 break;
01c4d84a 438 }
fd30f88e 439 ncl[nsupcl]++;
3edbbba2 440 }
fd30f88e 441
442 AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
443 incr+1,nsupcl+1));
444 id = -1;
445 icl = -1;
446 for(i = 0; i <= nsupcl; i++)
447 {
448 if(ncl[i] == 0)
01c4d84a 449 {
fd30f88e 450 id++;
451 icl++;
452 if (clno >= 5000)
453 {
454 AliWarning("RefClust: Too many clusters! more than 5000");
455 return;
456 }
457 clno++;
458 i1 = fInfcl[1][id];
459 i2 = fInfcl[2][id];
460
461 clusdata[0] = fCoord[0][i1][i2];
462 clusdata[1] = fCoord[1][i1][i2];
463 clusdata[2] = fEdepCell[i1][i2];
464 clusdata[3] = 1.;
465 clusdata[4] = 0.5;
466 clusdata[5] = 0.0;
467 pmdcludata = new AliPMDcludata(clusdata);
468 pmdclucont->Add(pmdcludata);
469
470 //association
471
472 fClTr[0][clno] = fCellTrNo[i1][i2];
473 for(Int_t icltr = 1; icltr < 14; icltr++)
474 {
475 fClTr[icltr][clno] = -1;
3edbbba2 476 }
3edbbba2 477 }
fd30f88e 478 else if(ncl[i] == 1)
479 {
480 id++;
481 icl++;
482 if (clno >= 5000)
483 {
484 AliWarning("RefClust: Too many clusters! more than 5000");
485 return;
486 }
487 clno++;
488 i1 = fInfcl[1][id];
489 i2 = fInfcl[2][id];
490 x1 = fCoord[0][i1][i2];
491 y1 = fCoord[1][i1][i2];
492 z1 = fEdepCell[i1][i2];
493
494 //asso
495 fClTr[0][clno] = fCellTrNo[i1][i2];
496 //
497
498 id = id+1;
499 i1 = fInfcl[1][id];
500 i2 = fInfcl[2][id];
501 x2 = fCoord[0][i1][i2];
502 y2 = fCoord[1][i1][i2];
503 z2 = fEdepCell[i1][i2];
504
505 //asso
506
507 fClTr[1][clno] = fCellTrNo[i1][i2];
508 for(Int_t icltr = 2; icltr < 14; icltr++)
509 {
510 fClTr[icltr][clno] = -1;
511 }
512 //
513
514 clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
515 clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
516 clusdata[2] = z1+z2;
517 clusdata[3] = 2.;
518 clusdata[4] = 0.5;
519 clusdata[5] = 0.0;
520 pmdcludata = new AliPMDcludata(clusdata);
521 pmdclucont->Add(pmdcludata);
3edbbba2 522 }
fd30f88e 523 else
524 {
525 //asso
526 for(Int_t icg = 0; icg < kdim; icg++)
527 {
528 cellCount[icg]=0;
529 }
530 //
531
532 id++;
533 iord[0] = 0;
534 // super-cluster of more than two cells - broken up into smaller
535 // clusters gaussian centers computed. (peaks separated by > 1 cell)
536 // Begin from cell having largest energy deposited This is first
537 // cluster center
538 i1 = fInfcl[1][id];
539 i2 = fInfcl[2][id];
540 x[0] = fCoord[0][i1][i2];
541 y[0] = fCoord[1][i1][i2];
542 z[0] = fEdepCell[i1][i2];
543
544 //asso
545 t[0] = fCellTrNo[i1][i2];
546 //
547
548 iord[0] = 0;
549 for(j = 1; j <= ncl[i]; j++)
550 {
551 id++;
552 i1 = fInfcl[1][id];
553 i2 = fInfcl[2][id];
554 iord[j] = j;
555 x[j] = fCoord[0][i1][i2];
556 y[j] = fCoord[1][i1][i2];
557 z[j] = fEdepCell[i1][i2];
01c4d84a 558 //asso
fd30f88e 559 t[j] = fCellTrNo[i1][i2];
01c4d84a 560 //
3edbbba2 561 }
fd30f88e 562
563
564 // arranging cells within supercluster in decreasing order
565
566 for(j = 1;j <= ncl[i]; j++)
567 {
568 itest = 0;
569 ihld = iord[j];
570 for(i1 = 0; i1 < j; i1++)
571 {
572 if(itest == 0 && z[iord[i1]] < z[ihld])
573 {
574 itest = 1;
575 for(i2 = j-1; i2 >= i1; i2--)
576 {
577 iord[i2+1] = iord[i2];
578 }
579 iord[i1] = ihld;
580 }
581 }
3edbbba2 582 }
fd30f88e 583 // compute the number of Gaussians and their centers ( first
584 // guess )
585 // centers must be separated by cells having smaller ener. dep.
586 // neighbouring centers should be either strong or well-separated
587 ig=0;
588 xc[ig] = x[iord[0]];
589 yc[ig] = y[iord[0]];
590 zc[ig] = z[iord[0]];
591 for(j = 1; j <= ncl[i]; j++)
592 {
593 itest = -1;
594 x1 = x[iord[j]];
595 y1 = y[iord[j]];
596 for(k = 0; k <= ig; k++)
597 {
598 x2 = xc[k];
599 y2 = yc[k];
600 rr = Distance(x1,y1,x2,y2);
601 if( rr >= 1.1 && rr < 1.8 && z[iord[j]] > zc[k]/4.)
602 {
603 itest++;
604 }
605 if( rr >= 1.8 && rr < 2.1 && z[iord[j]] > zc[k]/10.)
606 {
607 itest++;
608 }
609 if( rr >= 2.1)
610 {
611 itest++;
612 }
613 }
614 if(itest == ig)
615 {
616 ig++;
617 xc[ig] = x1;
618 yc[ig] = y1;
619 zc[ig] = z[iord[j]];
620 }
3edbbba2 621 }
fd30f88e 622
623 GaussFit(ncl[i], ig, x[0], y[0] ,z[0], xc[0], yc[0], zc[0], rc[0]);
624 icl += ig+1;
625 // compute the number of cells belonging to each cluster.
626 // cell is shared between several clusters ( if they are equidistant
627 // from it ) in the ratio of cluster energy deposition
628 for(j = 0; j <= ig; j++)
01c4d84a 629 {
fd30f88e 630 cells[j]=0.;
631 }
632 if(ig > 0)
633 {
634 for(j = 0; j <= ncl[i]; j++)
635 {
636 lev1[j] = 0;
637 lev2[j] = 0;
638 for(k = 0; k <= ig; k++)
639 {
640 dist = Distance(x[j], y[j], xc[k], yc[k]);
641 if(dist < TMath::Sqrt(3.) )
642 {
643 //asso
644 fClTr[cellCount[k]][clno+k+1] = t[j];
645 cellCount[k]++;
646 //
647 lev1[0]++;
648 i1 = lev1[0];
649 lev1[i1] = k;
650 }
651 else
652 {
653 if(dist < 2.1)
654 {
655 lev2[0]++;
656 i1 = lev2[0];
657 lev2[i1] = k;
658 }
659 }
660 }
661 if(lev1[0] != 0)
662 {
663 if(lev1[0] == 1)
664 {
665 cells[lev1[1]]++;
666 }
667 else
668 {
669 sum=0.;
670 for(k = 1; k <= lev1[0]; k++)
671 {
672 sum += zc[lev1[k]];
673 }
674 for(k = 1; k <= lev1[0]; k++)
675 {
676 cells[lev1[k]] += zc[lev1[k]]/sum;
677 }
678 }
679 }
680 else
681 {
682 if(lev2[0] == 0)
683 {
684 cells[lev2[1]]++;
685 }
686 else
687 {
688 sum=0.;
689 for( k = 1; k <= lev2[0]; k++)
690 {
691 sum += zc[lev2[k]];
692 }
693 for(k = 1; k <= lev2[0]; k++)
694 {
695 cells[lev2[k]] += zc[lev2[k]]/sum;
696 }
697 }
698 }
699 }
700 }
701
702 // zero rest of the cell array
703 //asso
704 for( k = 0; k <= ig; k++)
705 {
706 for(Int_t icltr = cellCount[k]; icltr < 14; icltr++)
707 {
708 fClTr[icltr][clno] = -1;
709 }
710 }
711 //
712
713 for(j = 0; j <= ig; j++)
714 {
715 clno++;
716 if (clno >= 5000)
717 {
718 AliWarning("RefClust: Too many clusters! more than 5000");
719 return;
720 }
721 clusdata[0] = xc[j];
722 clusdata[1] = yc[j];
723 clusdata[2] = zc[j];
724 clusdata[4] = rc[j];
725 clusdata[5] = 0.0;
726 if(ig == 0)
727 {
728 clusdata[3] = ncl[i];
729 }
730 else
731 {
732 clusdata[3] = cells[j];
733 }
734 pmdcludata = new AliPMDcludata(clusdata);
735 pmdclucont->Add(pmdcludata);
01c4d84a 736 }
737 }
3edbbba2 738 }
3edbbba2 739}
740// ------------------------------------------------------------------------ //
fd30f88e 741void AliPMDClusteringV1::GaussFit(Int_t ncell, Int_t nclust, Double_t &x,
742 Double_t &y ,Double_t &z, Double_t &xc,
743 Double_t &yc, Double_t &zc, Double_t &rc)
3edbbba2 744{
745 // Does gaussian fitting
746 //
fd30f88e 747 const Int_t kdim = 4500;
748 Int_t i, j, i1, i2, novar, idd, jj;
749 Double_t sum, dx, dy, str, str1, aint, sum1, rr, dum;
750 Double_t x1, x2, y1, y2;
751 Double_t xx[kdim], yy[kdim], zz[kdim], xxc[kdim], yyc[kdim];
752 Double_t a[kdim], b[kdim], c[kdim], d[kdim], ha[kdim], hb[kdim];
753 Double_t hc[kdim], hd[kdim], zzc[kdim], rrc[kdim];
754 Int_t neib[kdim][50];
755
756 TRandom rnd;
757
3edbbba2 758 str = 0.;
759 str1 = 0.;
760 rr = 0.3;
761 novar = 0;
fd30f88e 762 j = 0;
3edbbba2 763
fd30f88e 764 for(i = 0; i <= ncell; i++)
3edbbba2 765 {
766 xx[i] = *(&x+i);
767 yy[i] = *(&y+i);
768 zz[i] = *(&z+i);
fd30f88e 769 str += zz[i];
3edbbba2 770 }
771 for(i=0; i<=nclust; i++)
772 {
773 xxc[i] = *(&xc+i);
774 yyc[i] = *(&yc+i);
775 zzc[i] = *(&zc+i);
fd30f88e 776 str1 += zzc[i];
3edbbba2 777 rrc[i] = 0.5;
778 }
fd30f88e 779 for(i = 0; i <= nclust; i++)
3edbbba2 780 {
781 zzc[i] = str/str1*zzc[i];
782 ha[i] = xxc[i];
783 hb[i] = yyc[i];
784 hc[i] = zzc[i];
785 hd[i] = rrc[i];
786 x1 = xxc[i];
787 y1 = yyc[i];
788 }
fd30f88e 789 for(i = 0; i <= ncell; i++)
790 {
791 idd = 0;
792 x1 = xx[i];
793 y1 = yy[i];
794 for(j = 0; j <= nclust; j++)
795 {
796 x2 = xxc[j];
797 y2 = yyc[j];
798 if(Distance(x1,y1,x2,y2) <= 3.)
799 {
800 idd++;
801 neib[i][idd] = j;
802 }
803 }
804 neib[i][0] = idd;
3edbbba2 805 }
fd30f88e 806 sum = 0.;
807 for(i1 = 0; i1 <= ncell; i1++)
808 {
809 aint = 0.;
810 idd = neib[i1][0];
811 for(i2 = 1; i2 <= idd; i2++)
812 {
813 jj = neib[i1][i2];
814 dx = xx[i1] - xxc[jj];
815 dy = yy[i1] - yyc[jj];
816 dum = rrc[j]*rrc[jj] + rr*rr;
817 aint += exp(-(dx*dx+dy*dy)/dum)*zzc[idd]*rr*rr/dum;
818 }
819 sum += (aint - zz[i1])*(aint - zz[i1])/str;
820 }
821 str1 = 0.;
822 for(i = 0; i <= nclust; i++)
823 {
824 a[i] = xxc[i] + 0.6*(rnd.Uniform() - 0.5);
825 b[i] = yyc[i] + 0.6*(rnd.Uniform() - 0.5);
826 c[i] = zzc[i]*(1.+ ( rnd.Uniform() - 0.5)*0.2);
827 str1 += zzc[i];
828 d[i] = rrc[i]*(1.+ ( rnd.Uniform() - 0.5)*0.1);
829
830 if(d[i] < 0.25)
831 {
832 d[i]=0.25;
833 }
3edbbba2 834 }
fd30f88e 835 for(i = 0; i <= nclust; i++)
836 {
837 c[i] = c[i]*str/str1;
3edbbba2 838 }
fd30f88e 839 sum1=0.;
840 for(i1 = 0; i1 <= ncell; i1++)
841 {
842 aint = 0.;
843 idd = neib[i1][0];
844 for(i2 = 1; i2 <= idd; i2++)
845 {
846 jj = neib[i1][i2];
847 dx = xx[i1] - a[jj];
848 dy = yy[i1] - b[jj];
849 dum = d[jj]*d[jj]+rr*rr;
850 aint += exp(-(dx*dx+dy*dy)/dum)*c[i2]*rr*rr/dum;
851 }
852 sum1 += (aint - zz[i1])*(aint - zz[i1])/str;
3edbbba2 853 }
854
fd30f88e 855 if(sum1 < sum)
856 {
857 for(i2 = 0; i2 <= nclust; i2++)
858 {
859 xxc[i2] = a[i2];
860 yyc[i2] = b[i2];
861 zzc[i2] = c[i2];
862 rrc[i2] = d[i2];
863 sum = sum1;
864 }
865 }
866 for(j = 0; j <= nclust; j++)
867 {
868 *(&xc+j) = xxc[j];
869 *(&yc+j) = yyc[j];
870 *(&zc+j) = zzc[j];
871 *(&rc+j) = rrc[j];
3edbbba2 872 }
3edbbba2 873}
874// ------------------------------------------------------------------------ //
fd30f88e 875Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1,
876 Double_t x2, Double_t y2)
3edbbba2 877{
fd30f88e 878 return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
3edbbba2 879}
880// ------------------------------------------------------------------------ //
881void AliPMDClusteringV1::SetEdepCut(Float_t decut)
882{
883 fCutoff = decut;
884}
885// ------------------------------------------------------------------------ //