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