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