Adding TestShuttle temp and log folder's getter and setter
[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
562718f9 147 Int_t ij = i + j*kNDIMX;
148
01c4d84a 149 if (ismn < 12)
150 {
562718f9 151 edepcell[ij] = celladc[jd][id];
152 fCellTrNo[i][j] = jd*10000+id; // for association
01c4d84a 153 }
154 else if (ismn >= 12 && ismn <= 23)
155 {
562718f9 156 edepcell[ij] = celladc[id][jd];
157 fCellTrNo[i][j] = id*10000+jd; // for association
01c4d84a 158 }
159
3edbbba2 160 }
161 }
562718f9 162
163 Int_t iord1[kNMX];
164 TMath::Sort(kNMX,edepcell,iord1);// order the data
165 cutoff = fCutoff; // cutoff to discard cells
166 ave = 0.;
fd30f88e 167 nmx1 = -1;
562718f9 168 for(i = 0;i < kNMX; i++)
3edbbba2 169 {
562718f9 170 if(edepcell[i] > 0.)
fd30f88e 171 {
562718f9 172 ave += edepcell[i];
fd30f88e 173 }
562718f9 174 if(edepcell[i] > cutoff )
fd30f88e 175 {
176 nmx1++;
177 }
3edbbba2 178 }
3edbbba2 179
180 AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
181
3edbbba2 182 if (nmx1 == 0) nmx1 = 1;
fd30f88e 183 ave = ave/nmx1;
3edbbba2 184
185 AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
186 kNMX,ave));
562718f9 187
188 incr = CrClust(ave, cutoff, nmx1,iord1, edepcell );
189
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 //
562718f9 382
383 AliPMDcludata *pmdcludata = 0;
384
fd30f88e 385 const Int_t kdim = 4500;
562718f9 386
387 Int_t i, j, k, i1, i2, id, icl, itest,ihld, ig, nsupcl,clno;
388 Int_t lev1[20], lev2[20];
389 Int_t ncl[kdim], iord[kdim], t[kdim];
390 Int_t clxy[15];
391
392 Int_t *cellCount;
393 Int_t **cellXY;
394
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];
399
562718f9 400
401 // Initialisation
fd30f88e 402 for(i = 0; i<kdim; i++)
403 {
404 t[i] = -1;
562718f9 405 ncl[i] = -1;
406 if (i < 6) clusdata[i] = 0.;
407 if (i < 15) clxy[i] = 0;
fd30f88e 408 }
562718f9 409
fd30f88e 410 // clno counts the final clusters
3edbbba2 411 // nsupcl = # of superclusters; ncl[i]= # of cells in supercluster i
412 // x, y and z store (x,y) coordinates of and energy deposited in a cell
413 // xc, yc store (x,y) coordinates of the cluster center
414 // zc stores the energy deposited in a cluster
415 // rc is cluster radius
fd30f88e 416
417 clno = -1;
3edbbba2 418 nsupcl = -1;
562718f9 419
fd30f88e 420 for(i = 0; i <= incr; i++)
421 {
422 if(fInfcl[0][i] != nsupcl)
01c4d84a 423 {
fd30f88e 424 nsupcl++;
01c4d84a 425 }
fd30f88e 426 if (nsupcl > kdim)
01c4d84a 427 {
fd30f88e 428 AliWarning("RefClust: Too many superclusters!");
429 nsupcl = kdim;
430 break;
01c4d84a 431 }
fd30f88e 432 ncl[nsupcl]++;
3edbbba2 433 }
fd30f88e 434
435 AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
436 incr+1,nsupcl+1));
437 id = -1;
438 icl = -1;
562718f9 439
fd30f88e 440 for(i = 0; i <= nsupcl; i++)
441 {
442 if(ncl[i] == 0)
01c4d84a 443 {
fd30f88e 444 id++;
445 icl++;
446 if (clno >= 5000)
447 {
448 AliWarning("RefClust: Too many clusters! more than 5000");
449 return;
450 }
451 clno++;
452 i1 = fInfcl[1][id];
453 i2 = fInfcl[2][id];
454
562718f9 455 Int_t i12 = i1 + i2*kNDIMX;
456
fd30f88e 457 clusdata[0] = fCoord[0][i1][i2];
458 clusdata[1] = fCoord[1][i1][i2];
562718f9 459 clusdata[2] = edepcell[i12];
fd30f88e 460 clusdata[3] = 1.;
461 clusdata[4] = 0.5;
462 clusdata[5] = 0.0;
562718f9 463
464 clxy[0] = fCellTrNo[i1][i2]; //association
465 for(Int_t icltr = 1; icltr < 15; icltr++)
fd30f88e 466 {
562718f9 467 clxy[icltr] = -1;
3edbbba2 468 }
562718f9 469 pmdcludata = new AliPMDcludata(clusdata,clxy);
470 fPMDclucont->Add(pmdcludata);
3edbbba2 471 }
fd30f88e 472 else if(ncl[i] == 1)
473 {
474 id++;
475 icl++;
476 if (clno >= 5000)
477 {
478 AliWarning("RefClust: Too many clusters! more than 5000");
479 return;
480 }
481 clno++;
482 i1 = fInfcl[1][id];
483 i2 = fInfcl[2][id];
562718f9 484 Int_t i12 = i1 + i2*kNDIMX;
485
fd30f88e 486 x1 = fCoord[0][i1][i2];
487 y1 = fCoord[1][i1][i2];
562718f9 488 z1 = edepcell[i12];
489 clxy[0] = fCellTrNo[i1][i2]; //asso
490 id++;
fd30f88e 491 i1 = fInfcl[1][id];
492 i2 = fInfcl[2][id];
562718f9 493
494 Int_t i22 = i1 + i2*kNDIMX;
fd30f88e 495 x2 = fCoord[0][i1][i2];
496 y2 = fCoord[1][i1][i2];
562718f9 497 z2 = edepcell[i22];
498 clxy[1] = fCellTrNo[i1][i2]; //asso
499 for(Int_t icltr = 2; icltr < 15; icltr++)
fd30f88e 500 {
562718f9 501 clxy[icltr] = -1;
fd30f88e 502 }
fd30f88e 503
504 clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
505 clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
506 clusdata[2] = z1+z2;
507 clusdata[3] = 2.;
508 clusdata[4] = 0.5;
509 clusdata[5] = 0.0;
562718f9 510 pmdcludata = new AliPMDcludata(clusdata,clxy);
511 fPMDclucont->Add(pmdcludata);
3edbbba2 512 }
fd30f88e 513 else
514 {
fd30f88e 515 id++;
516 iord[0] = 0;
517 // super-cluster of more than two cells - broken up into smaller
518 // clusters gaussian centers computed. (peaks separated by > 1 cell)
519 // Begin from cell having largest energy deposited This is first
520 // cluster center
521 i1 = fInfcl[1][id];
522 i2 = fInfcl[2][id];
562718f9 523 Int_t i12 = i1 + i2*kNDIMX;
524
fd30f88e 525 x[0] = fCoord[0][i1][i2];
526 y[0] = fCoord[1][i1][i2];
562718f9 527 z[0] = edepcell[i12];
528
529 t[0] = fCellTrNo[i1][i2]; //asso
fd30f88e 530
531 iord[0] = 0;
532 for(j = 1; j <= ncl[i]; j++)
533 {
534 id++;
535 i1 = fInfcl[1][id];
536 i2 = fInfcl[2][id];
562718f9 537 Int_t i12 = i1 + i2*kNDIMX;
538
fd30f88e 539 iord[j] = j;
540 x[j] = fCoord[0][i1][i2];
541 y[j] = fCoord[1][i1][i2];
562718f9 542 z[j] = edepcell[i12];
543
544 t[j] = fCellTrNo[i1][i2]; //asso
3edbbba2 545 }
fd30f88e 546
fd30f88e 547 // arranging cells within supercluster in decreasing order
548
549 for(j = 1;j <= ncl[i]; j++)
550 {
551 itest = 0;
552 ihld = iord[j];
553 for(i1 = 0; i1 < j; i1++)
554 {
555 if(itest == 0 && z[iord[i1]] < z[ihld])
556 {
557 itest = 1;
558 for(i2 = j-1; i2 >= i1; i2--)
559 {
560 iord[i2+1] = iord[i2];
561 }
562 iord[i1] = ihld;
563 }
564 }
3edbbba2 565 }
fd30f88e 566 // compute the number of Gaussians and their centers ( first
567 // guess )
568 // centers must be separated by cells having smaller ener. dep.
569 // neighbouring centers should be either strong or well-separated
562718f9 570 ig = 0;
fd30f88e 571 xc[ig] = x[iord[0]];
572 yc[ig] = y[iord[0]];
573 zc[ig] = z[iord[0]];
574 for(j = 1; j <= ncl[i]; j++)
575 {
576 itest = -1;
577 x1 = x[iord[j]];
578 y1 = y[iord[j]];
579 for(k = 0; k <= ig; k++)
580 {
581 x2 = xc[k];
582 y2 = yc[k];
583 rr = Distance(x1,y1,x2,y2);
562718f9 584 if(rr >= 1.1 && rr < 1.8 && z[iord[j]] > zc[k]/4.)itest++;
585 if(rr >= 1.8 && rr < 2.1 && z[iord[j]] > zc[k]/10.)itest++;
586 if( rr >= 2.1)itest++;
fd30f88e 587 }
588 if(itest == ig)
589 {
590 ig++;
591 xc[ig] = x1;
592 yc[ig] = y1;
593 zc[ig] = z[iord[j]];
594 }
3edbbba2 595 }
fd30f88e 596 GaussFit(ncl[i], ig, x[0], y[0] ,z[0], xc[0], yc[0], zc[0], rc[0]);
597 icl += ig+1;
598 // compute the number of cells belonging to each cluster.
599 // cell is shared between several clusters ( if they are equidistant
600 // from it ) in the ratio of cluster energy deposition
562718f9 601
602 Int_t jj = 15;
603 cellCount = new Int_t [ig+1];
604 cellXY = new Int_t *[jj];
605 for(Int_t ij = 0; ij < 15; ij++) cellXY[ij] = new Int_t [ig+1];
606
fd30f88e 607 for(j = 0; j <= ig; j++)
01c4d84a 608 {
562718f9 609 cellCount[j] = 0;
610 cells[j] = 0.;
fd30f88e 611 }
562718f9 612
fd30f88e 613 if(ig > 0)
614 {
615 for(j = 0; j <= ncl[i]; j++)
616 {
617 lev1[j] = 0;
618 lev2[j] = 0;
619 for(k = 0; k <= ig; k++)
620 {
621 dist = Distance(x[j], y[j], xc[k], yc[k]);
622 if(dist < TMath::Sqrt(3.) )
623 {
624 //asso
562718f9 625 if (cellCount[k] < 15)
626 {
627 cellXY[cellCount[k]][k] = t[j];
628 }
fd30f88e 629 cellCount[k]++;
630 //
631 lev1[0]++;
632 i1 = lev1[0];
633 lev1[i1] = k;
634 }
635 else
636 {
637 if(dist < 2.1)
638 {
639 lev2[0]++;
640 i1 = lev2[0];
641 lev2[i1] = k;
642 }
643 }
644 }
645 if(lev1[0] != 0)
646 {
647 if(lev1[0] == 1)
648 {
649 cells[lev1[1]]++;
650 }
651 else
652 {
653 sum=0.;
654 for(k = 1; k <= lev1[0]; k++)
655 {
656 sum += zc[lev1[k]];
657 }
658 for(k = 1; k <= lev1[0]; k++)
659 {
660 cells[lev1[k]] += zc[lev1[k]]/sum;
661 }
662 }
663 }
664 else
665 {
666 if(lev2[0] == 0)
667 {
668 cells[lev2[1]]++;
669 }
670 else
671 {
672 sum=0.;
673 for( k = 1; k <= lev2[0]; k++)
674 {
675 sum += zc[lev2[k]];
676 }
677 for(k = 1; k <= lev2[0]; k++)
678 {
679 cells[lev2[k]] += zc[lev2[k]]/sum;
680 }
681 }
682 }
683 }
684 }
685
686 // zero rest of the cell array
687 //asso
688 for( k = 0; k <= ig; k++)
689 {
562718f9 690 for(Int_t icltr = cellCount[k]; icltr < 15; icltr++)
fd30f88e 691 {
562718f9 692 cellXY[icltr][k] = -1;
fd30f88e 693 }
694 }
695 //
696
697 for(j = 0; j <= ig; j++)
698 {
699 clno++;
700 if (clno >= 5000)
701 {
702 AliWarning("RefClust: Too many clusters! more than 5000");
703 return;
704 }
705 clusdata[0] = xc[j];
706 clusdata[1] = yc[j];
707 clusdata[2] = zc[j];
708 clusdata[4] = rc[j];
709 clusdata[5] = 0.0;
710 if(ig == 0)
711 {
712 clusdata[3] = ncl[i];
713 }
714 else
715 {
716 clusdata[3] = cells[j];
717 }
562718f9 718
719
720 for (Int_t ii=0; ii < 15; ii++)
721 {
722 clxy[ii] = cellXY[ii][j];
723 }
724 pmdcludata = new AliPMDcludata(clusdata,clxy);
725 fPMDclucont->Add(pmdcludata);
01c4d84a 726 }
562718f9 727 delete [] cellCount;
728 for(Int_t jj = 0; jj < 15; jj++)delete [] cellXY[jj];
729 delete [] cellXY;
01c4d84a 730 }
3edbbba2 731 }
3edbbba2 732}
733// ------------------------------------------------------------------------ //
fd30f88e 734void AliPMDClusteringV1::GaussFit(Int_t ncell, Int_t nclust, Double_t &x,
735 Double_t &y ,Double_t &z, Double_t &xc,
736 Double_t &yc, Double_t &zc, Double_t &rc)
3edbbba2 737{
738 // Does gaussian fitting
739 //
562718f9 740
fd30f88e 741 const Int_t kdim = 4500;
742 Int_t i, j, i1, i2, novar, idd, jj;
562718f9 743 Int_t neib[kdim][50];
744
fd30f88e 745 Double_t sum, dx, dy, str, str1, aint, sum1, rr, dum;
746 Double_t x1, x2, y1, y2;
747 Double_t xx[kdim], yy[kdim], zz[kdim], xxc[kdim], yyc[kdim];
748 Double_t a[kdim], b[kdim], c[kdim], d[kdim], ha[kdim], hb[kdim];
749 Double_t hc[kdim], hd[kdim], zzc[kdim], rrc[kdim];
fd30f88e 750
751 TRandom rnd;
752
3edbbba2 753 str = 0.;
754 str1 = 0.;
755 rr = 0.3;
756 novar = 0;
562718f9 757 j = 0;
3edbbba2 758
fd30f88e 759 for(i = 0; i <= ncell; i++)
3edbbba2 760 {
761 xx[i] = *(&x+i);
762 yy[i] = *(&y+i);
763 zz[i] = *(&z+i);
fd30f88e 764 str += zz[i];
3edbbba2 765 }
766 for(i=0; i<=nclust; i++)
767 {
768 xxc[i] = *(&xc+i);
769 yyc[i] = *(&yc+i);
770 zzc[i] = *(&zc+i);
fd30f88e 771 str1 += zzc[i];
3edbbba2 772 rrc[i] = 0.5;
773 }
fd30f88e 774 for(i = 0; i <= nclust; i++)
3edbbba2 775 {
776 zzc[i] = str/str1*zzc[i];
777 ha[i] = xxc[i];
778 hb[i] = yyc[i];
779 hc[i] = zzc[i];
780 hd[i] = rrc[i];
781 x1 = xxc[i];
782 y1 = yyc[i];
783 }
562718f9 784
fd30f88e 785 for(i = 0; i <= ncell; i++)
786 {
787 idd = 0;
788 x1 = xx[i];
789 y1 = yy[i];
790 for(j = 0; j <= nclust; j++)
791 {
792 x2 = xxc[j];
793 y2 = yyc[j];
794 if(Distance(x1,y1,x2,y2) <= 3.)
795 {
796 idd++;
797 neib[i][idd] = j;
798 }
799 }
800 neib[i][0] = idd;
3edbbba2 801 }
fd30f88e 802 sum = 0.;
803 for(i1 = 0; i1 <= ncell; i1++)
804 {
805 aint = 0.;
806 idd = neib[i1][0];
807 for(i2 = 1; i2 <= idd; i2++)
808 {
809 jj = neib[i1][i2];
810 dx = xx[i1] - xxc[jj];
811 dy = yy[i1] - yyc[jj];
812 dum = rrc[j]*rrc[jj] + rr*rr;
813 aint += exp(-(dx*dx+dy*dy)/dum)*zzc[idd]*rr*rr/dum;
814 }
815 sum += (aint - zz[i1])*(aint - zz[i1])/str;
816 }
817 str1 = 0.;
562718f9 818
fd30f88e 819 for(i = 0; i <= nclust; i++)
820 {
821 a[i] = xxc[i] + 0.6*(rnd.Uniform() - 0.5);
822 b[i] = yyc[i] + 0.6*(rnd.Uniform() - 0.5);
823 c[i] = zzc[i]*(1.+ ( rnd.Uniform() - 0.5)*0.2);
824 str1 += zzc[i];
825 d[i] = rrc[i]*(1.+ ( rnd.Uniform() - 0.5)*0.1);
826
827 if(d[i] < 0.25)
828 {
829 d[i]=0.25;
830 }
3edbbba2 831 }
fd30f88e 832 for(i = 0; i <= nclust; i++)
833 {
834 c[i] = c[i]*str/str1;
3edbbba2 835 }
fd30f88e 836 sum1=0.;
562718f9 837
fd30f88e 838 for(i1 = 0; i1 <= ncell; i1++)
839 {
840 aint = 0.;
841 idd = neib[i1][0];
842 for(i2 = 1; i2 <= idd; i2++)
843 {
844 jj = neib[i1][i2];
845 dx = xx[i1] - a[jj];
846 dy = yy[i1] - b[jj];
847 dum = d[jj]*d[jj]+rr*rr;
848 aint += exp(-(dx*dx+dy*dy)/dum)*c[i2]*rr*rr/dum;
849 }
850 sum1 += (aint - zz[i1])*(aint - zz[i1])/str;
3edbbba2 851 }
852
fd30f88e 853 if(sum1 < sum)
854 {
855 for(i2 = 0; i2 <= nclust; i2++)
856 {
857 xxc[i2] = a[i2];
858 yyc[i2] = b[i2];
859 zzc[i2] = c[i2];
860 rrc[i2] = d[i2];
861 sum = sum1;
862 }
863 }
864 for(j = 0; j <= nclust; j++)
865 {
866 *(&xc+j) = xxc[j];
867 *(&yc+j) = yyc[j];
868 *(&zc+j) = zzc[j];
869 *(&rc+j) = rrc[j];
3edbbba2 870 }
3edbbba2 871}
872// ------------------------------------------------------------------------ //
fd30f88e 873Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1,
874 Double_t x2, Double_t y2)
3edbbba2 875{
fd30f88e 876 return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
3edbbba2 877}
878// ------------------------------------------------------------------------ //
879void AliPMDClusteringV1::SetEdepCut(Float_t decut)
880{
881 fCutoff = decut;
882}
883// ------------------------------------------------------------------------ //