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