]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PMD/AliPMDClusteringV1.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
added prediction processor framework and implementations; moved MUON preprocessors
[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
2c1131dd 57
3edbbba2 58ClassImp(AliPMDClusteringV1)
59
60const Double_t AliPMDClusteringV1::fgkSqroot3by2=0.8660254; // sqrt(3.)/2.
61
62AliPMDClusteringV1::AliPMDClusteringV1():
562718f9 63 fPMDclucont(new TObjArray()),
3edbbba2 64 fCutoff(0.0)
65{
fd30f88e 66 for(Int_t i = 0; i < kNDIMX; i++)
3edbbba2 67 {
fd30f88e 68 for(Int_t j = 0; j < kNDIMY; j++)
3edbbba2 69 {
70 fCoord[0][i][j] = i+j/2.;
71 fCoord[1][i][j] = fgkSqroot3by2*j;
3edbbba2 72 }
73 }
74}
75// ------------------------------------------------------------------------ //
562718f9 76AliPMDClusteringV1::AliPMDClusteringV1(const AliPMDClusteringV1& pmdclv1):
77 AliPMDClustering(pmdclv1),
78 fPMDclucont(0),
79 fCutoff(0)
80{
81 // copy constructor
82 AliError("Copy constructor not allowed ");
83
84}
85// ------------------------------------------------------------------------ //
86AliPMDClusteringV1 &AliPMDClusteringV1::operator=(const AliPMDClusteringV1& /*pmdclv1*/)
87{
88 // copy constructor
89 AliError("Assignment operator not allowed ");
90 return *this;
91}
92// ------------------------------------------------------------------------ //
3edbbba2 93AliPMDClusteringV1::~AliPMDClusteringV1()
94{
562718f9 95 delete fPMDclucont;
3edbbba2 96}
97// ------------------------------------------------------------------------ //
fd30f88e 98void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn,
99 Double_t celladc[48][96], TObjArray *pmdcont)
3edbbba2 100{
101 // main function to call other necessary functions to do clustering
102 //
01c4d84a 103
562718f9 104 AliPMDcluster *pmdcl = 0;
3edbbba2 105
2c1131dd 106 const float ktwobysqrt3 = 1.1547; // 2./sqrt(3.)
107 const Int_t kNmaxCell = 19; // # of cells surrounding a cluster center
108
562718f9 109 Int_t i, j, nmx1, incr, id, jd;
2c1131dd 110 Int_t celldataX[kNmaxCell], celldataY[kNmaxCell];
562718f9 111 Float_t clusdata[6];
112 Double_t cutoff, ave;
113 Double_t edepcell[kNMX];
9c294f00 114
115
2c1131dd 116 Double_t *cellenergy = new Double_t [11424];
9c294f00 117
3edbbba2 118
01c4d84a 119 // ndimXr and ndimYr are different because of different module size
120
fd30f88e 121 Int_t ndimXr = 0;
122 Int_t ndimYr = 0;
01c4d84a 123
124 if (ismn < 12)
125 {
126 ndimXr = 96;
127 ndimYr = 48;
128 }
129 else if (ismn >= 12 && ismn <= 23)
130 {
131 ndimXr = 48;
132 ndimYr = 96;
133 }
134
78fc1b96 135 for (i =0; i < 11424; i++)
9c294f00 136 {
137 cellenergy[i] = 0.;
138 }
139
140
562718f9 141 Int_t kk = 0;
78fc1b96 142 for (i = 0; i < kNDIMX; i++)
3edbbba2 143 {
78fc1b96 144 for (j = 0; j < kNDIMY; j++)
01c4d84a 145 {
562718f9 146 edepcell[kk] = 0.;
147 kk++;
01c4d84a 148 }
149 }
150
151 for (id = 0; id < ndimXr; id++)
152 {
153 for (jd = 0; jd < ndimYr; jd++)
3edbbba2 154 {
fd30f88e 155 j = jd;
156 i = id+(ndimYr/2-1)-(jd/2);
01c4d84a 157
9c294f00 158 Int_t ij = i + j*kNDIMX;
9c294f00 159
01c4d84a 160 if (ismn < 12)
161 {
2c1131dd 162 cellenergy[ij] = celladc[jd][id];
01c4d84a 163 }
164 else if (ismn >= 12 && ismn <= 23)
165 {
2c1131dd 166 cellenergy[ij] = celladc[id][jd];
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
2c1131dd 180 cutoff = fCutoff; // cutoff to discard cells
562718f9 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 //
2c1131dd 225
01c4d84a 226 if (ismn < 12)
227 {
228 clusdata[0] = cluX0 - (24-1) + cluY0/2.;
229 }
fd30f88e 230 else if ( ismn >= 12 && ismn <= 23)
01c4d84a 231 {
232 clusdata[0] = cluX0 - (48-1) + cluY0/2.;
233 }
fd30f88e 234
2c1131dd 235
01c4d84a 236 clusdata[1] = cluY0;
237 clusdata[2] = cluADC;
238 clusdata[3] = cluCELLS;
fd30f88e 239 clusdata[4] = cluSIGX;
240 clusdata[5] = cluSIGY;
241
3edbbba2 242 //
243 // Cells associated with a cluster
244 //
562718f9 245
2c1131dd 246 for (Int_t ihit = 0; ihit < kNmaxCell; ihit++)
3edbbba2 247 {
2c1131dd 248 Int_t cellrow = pmdcludata->GetCellXY(ihit)/10000;
249 Int_t cellcol = pmdcludata->GetCellXY(ihit)%10000;
250
01c4d84a 251 if (ismn < 12)
252 {
2c1131dd 253 celldataX[ihit] = cellrow - (24-1) + int(cellcol/2.);
01c4d84a 254 }
255 else if (ismn >= 12 && ismn <= 23)
256 {
2c1131dd 257 celldataX[ihit] = cellrow - (48-1) + int(cellcol/2.);
01c4d84a 258 }
2c1131dd 259
260 celldataY[ihit] = cellcol;
3edbbba2 261 }
2c1131dd 262
263
3edbbba2 264 pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY);
265 pmdcont->Add(pmdcl);
266 }
fd30f88e 267
91e6e2a0 268 fPMDclucont->Delete();
3edbbba2 269}
270// ------------------------------------------------------------------------ //
562718f9 271Int_t AliPMDClusteringV1::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
272 Int_t iord1[], Double_t edepcell[])
3edbbba2 273{
fd30f88e 274 // Does crude clustering
3edbbba2 275 // Finds out only the big patch by just searching the
276 // connected cells
277 //
2c1131dd 278 const Int_t kndim = 4609;
279 Int_t i,j,k,id1,id2,icl, numcell, clust[2][kndim];
fd30f88e 280 Int_t jd1,jd2, icell, cellcount;
281 static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
282
3edbbba2 283 AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
284
fd30f88e 285 for (j = 0; j < kNDIMX; j++)
286 {
287 for(k = 0; k < kNDIMY; k++)
288 {
289 fInfocl[0][j][k] = 0;
290 fInfocl[1][j][k] = 0;
291 }
292 }
293 for(i=0; i < kNMX; i++)
294 {
295 fInfcl[0][i] = -1;
562718f9 296
297 j = iord1[i];
298 id2 = j/kNDIMX;
299 id1 = j-id2*kNDIMX;
300
301 if(edepcell[j] <= cutoff)
fd30f88e 302 {
303 fInfocl[0][id1][id2] = -1;
304 }
3edbbba2 305 }
562718f9 306
3edbbba2 307 // ---------------------------------------------------------------
308 // crude clustering begins. Start with cell having largest adc
309 // count and loop over the cells in descending order of adc count
310 // ---------------------------------------------------------------
562718f9 311
fd30f88e 312 icl = -1;
313 cellcount = -1;
3edbbba2 314
fd30f88e 315 for(icell = 0; icell <= nmx1; icell++)
316 {
562718f9 317 j = iord1[icell];
318 id2 = j/kNDIMX;
319 id1 = j-id2*kNDIMX;
320
fd30f88e 321 if(fInfocl[0][id1][id2] == 0 )
322 {
323 icl++;
324 numcell = 0;
325 cellcount++;
326 fInfocl[0][id1][id2] = 1;
327 fInfocl[1][id1][id2] = icl;
328 fInfcl[0][cellcount] = icl;
329 fInfcl[1][cellcount] = id1;
330 fInfcl[2][cellcount] = id2;
331
332 clust[0][numcell] = id1;
333 clust[1][numcell] = id2;
334
2c1131dd 335 for(i = 1; i < kndim; i++)
fd30f88e 336 {
337 clust[0][i]=0;
338 }
339 // ---------------------------------------------------------------
340 // check for adc count in neib. cells. If ne 0 put it in this clust
341 // ---------------------------------------------------------------
342 for(i = 0; i < 6; i++)
343 {
344 jd1 = id1 + neibx[i];
345 jd2 = id2 + neiby[i];
346 if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
347 fInfocl[0][jd1][jd2] == 0)
348 {
349 numcell++;
350 fInfocl[0][jd1][jd2] = 2;
351 fInfocl[1][jd1][jd2] = icl;
352 clust[0][numcell] = jd1;
353 clust[1][numcell] = jd2;
354 cellcount++;
355 fInfcl[0][cellcount] = icl;
356 fInfcl[1][cellcount] = jd1;
357 fInfcl[2][cellcount] = jd2;
358 }
359 }
360 // ---------------------------------------------------------------
361 // check adc count for neighbour's neighbours recursively and
362 // if nonzero, add these to the cluster.
363 // ---------------------------------------------------------------
2c1131dd 364 for(i = 1; i < kndim;i++)
fd30f88e 365 {
366 if(clust[0][i] != 0)
367 {
368 id1 = clust[0][i];
369 id2 = clust[1][i];
370 for(j = 0; j < 6 ; j++)
371 {
372 jd1 = id1 + neibx[j];
373 jd2 = id2 + neiby[j];
374 if( (jd1 >= 0 && jd1 < kNDIMX) &&
375 (jd2 >= 0 && jd2 < kNDIMY) &&
376 fInfocl[0][jd1][jd2] == 0 )
377 {
378 fInfocl[0][jd1][jd2] = 2;
379 fInfocl[1][jd1][jd2] = icl;
380 numcell++;
381 clust[0][numcell] = jd1;
382 clust[1][numcell] = jd2;
383 cellcount++;
384 fInfcl[0][cellcount] = icl;
385 fInfcl[1][cellcount] = jd1;
386 fInfcl[2][cellcount] = jd2;
387 }
388 }
389 }
3edbbba2 390 }
3edbbba2 391 }
3edbbba2 392 }
3edbbba2 393 return cellcount;
394}
395// ------------------------------------------------------------------------ //
562718f9 396void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
3edbbba2 397{
398 // Does the refining of clusters
399 // Takes the big patch and does gaussian fitting and
400 // finds out the more refined clusters
401 //
22d01768 402
22d01768 403 AliPMDcludata *pmdcludata = 0;
2c1131dd 404
405 const Int_t kNmaxCell = 19; // # of cells surrounding a cluster center
406
407 Int_t ndim = incr + 1;
408
409 Int_t *ncl = 0x0;
410 Int_t *clxy = 0x0;
411 Int_t i12, i22;
412 Int_t i, j, k, i1, i2, id, icl, itest,ihld, ig, nsupcl,clno, t1, t2;
562718f9 413 Float_t clusdata[6];
2c1131dd 414 Double_t x1, y1, z1, x2, y2, z2, rr;
415
416 ncl = new Int_t [ndim];
417 clxy = new Int_t [kNmaxCell];
418
562718f9 419 // Initialisation
2c1131dd 420 for(i = 0; i<ndim; i++)
421 {
422 ncl[i] = -1;
562718f9 423 if (i < 6) clusdata[i] = 0.;
2c1131dd 424 if (i < kNmaxCell) clxy[i] = 0;
fd30f88e 425 }
562718f9 426
fd30f88e 427 // clno counts the final clusters
3edbbba2 428 // nsupcl = # of superclusters; ncl[i]= # of cells in supercluster i
429 // x, y and z store (x,y) coordinates of and energy deposited in a cell
430 // xc, yc store (x,y) coordinates of the cluster center
431 // zc stores the energy deposited in a cluster
432 // rc is cluster radius
fd30f88e 433
434 clno = -1;
3edbbba2 435 nsupcl = -1;
562718f9 436
fd30f88e 437 for(i = 0; i <= incr; i++)
438 {
439 if(fInfcl[0][i] != nsupcl)
01c4d84a 440 {
fd30f88e 441 nsupcl++;
01c4d84a 442 }
2c1131dd 443 if (nsupcl > ndim)
01c4d84a 444 {
fd30f88e 445 AliWarning("RefClust: Too many superclusters!");
2c1131dd 446 nsupcl = ndim;
fd30f88e 447 break;
01c4d84a 448 }
fd30f88e 449 ncl[nsupcl]++;
3edbbba2 450 }
fd30f88e 451
452 AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
453 incr+1,nsupcl+1));
454 id = -1;
455 icl = -1;
562718f9 456
fd30f88e 457 for(i = 0; i <= nsupcl; i++)
458 {
459 if(ncl[i] == 0)
01c4d84a 460 {
fd30f88e 461 id++;
462 icl++;
2c1131dd 463 if (clno >= 4608)
fd30f88e 464 {
2c1131dd 465 AliWarning("RefClust: Too many clusters! more than 4608");
fd30f88e 466 return;
467 }
468 clno++;
469 i1 = fInfcl[1][id];
470 i2 = fInfcl[2][id];
471
78fc1b96 472 i12 = i1 + i2*kNDIMX;
562718f9 473
fd30f88e 474 clusdata[0] = fCoord[0][i1][i2];
475 clusdata[1] = fCoord[1][i1][i2];
562718f9 476 clusdata[2] = edepcell[i12];
fd30f88e 477 clusdata[3] = 1.;
478 clusdata[4] = 0.5;
479 clusdata[5] = 0.0;
562718f9 480
2c1131dd 481 clxy[0] = i1*10000 + i2;
482
483 for(Int_t icltr = 1; icltr < kNmaxCell; icltr++)
fd30f88e 484 {
562718f9 485 clxy[icltr] = -1;
3edbbba2 486 }
562718f9 487 pmdcludata = new AliPMDcludata(clusdata,clxy);
488 fPMDclucont->Add(pmdcludata);
3edbbba2 489 }
fd30f88e 490 else if(ncl[i] == 1)
491 {
492 id++;
493 icl++;
2c1131dd 494 if (clno >= 4608)
fd30f88e 495 {
2c1131dd 496 AliWarning("RefClust: Too many clusters! more than 4608");
fd30f88e 497 return;
498 }
499 clno++;
500 i1 = fInfcl[1][id];
501 i2 = fInfcl[2][id];
78fc1b96 502 i12 = i1 + i2*kNDIMX;
562718f9 503
fd30f88e 504 x1 = fCoord[0][i1][i2];
505 y1 = fCoord[1][i1][i2];
562718f9 506 z1 = edepcell[i12];
2c1131dd 507
508 clxy[0] = i1*10000 + i2;
509
562718f9 510 id++;
fd30f88e 511 i1 = fInfcl[1][id];
512 i2 = fInfcl[2][id];
562718f9 513
2c1131dd 514 i22 = i1 + i2*kNDIMX;
fd30f88e 515 x2 = fCoord[0][i1][i2];
516 y2 = fCoord[1][i1][i2];
562718f9 517 z2 = edepcell[i22];
2c1131dd 518
519 clxy[1] = i1*10000 + i2;
520
521
522 for(Int_t icltr = 2; icltr < kNmaxCell; icltr++)
fd30f88e 523 {
562718f9 524 clxy[icltr] = -1;
fd30f88e 525 }
fd30f88e 526
527 clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
528 clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
529 clusdata[2] = z1+z2;
530 clusdata[3] = 2.;
531 clusdata[4] = 0.5;
532 clusdata[5] = 0.0;
562718f9 533 pmdcludata = new AliPMDcludata(clusdata,clxy);
534 fPMDclucont->Add(pmdcludata);
3edbbba2 535 }
fd30f88e 536 else
537 {
2c1131dd 538
539 Int_t *iord, *tc, *t;
540 Double_t *x, *y, *z, *xc, *yc, *zc;
541
542 iord = new Int_t [ncl[i]+1];
543 tc = new Int_t [ncl[i]+1];
544 t = new Int_t [ncl[i]+1];
545
546 x = new Double_t [ncl[i]+1];
547 y = new Double_t [ncl[i]+1];
548 z = new Double_t [ncl[i]+1];
549 xc = new Double_t [ncl[i]+1];
550 yc = new Double_t [ncl[i]+1];
551 zc = new Double_t [ncl[i]+1];
552
553 for( k = 0; k < ncl[i]+1; k++)
554 {
555 iord[k] = -1;
556 t[k] = -1;
557 tc[k] = -1;
558 x[k] = -1;
559 y[k] = -1;
560 z[k] = -1;
561 xc[k] = -1;
562 yc[k] = -1;
563 zc[k] = -1;
564 }
fd30f88e 565 id++;
fd30f88e 566 // super-cluster of more than two cells - broken up into smaller
567 // clusters gaussian centers computed. (peaks separated by > 1 cell)
568 // Begin from cell having largest energy deposited This is first
569 // cluster center
570 i1 = fInfcl[1][id];
571 i2 = fInfcl[2][id];
78fc1b96 572 i12 = i1 + i2*kNDIMX;
562718f9 573
fd30f88e 574 x[0] = fCoord[0][i1][i2];
575 y[0] = fCoord[1][i1][i2];
562718f9 576 z[0] = edepcell[i12];
2c1131dd 577 t[0] = i1*10000 + i2;
fd30f88e 578
2c1131dd 579
fd30f88e 580 iord[0] = 0;
581 for(j = 1; j <= ncl[i]; j++)
582 {
583 id++;
584 i1 = fInfcl[1][id];
585 i2 = fInfcl[2][id];
78fc1b96 586 i12 = i1 + i2*kNDIMX;
562718f9 587
fd30f88e 588 iord[j] = j;
589 x[j] = fCoord[0][i1][i2];
590 y[j] = fCoord[1][i1][i2];
562718f9 591 z[j] = edepcell[i12];
2c1131dd 592 t[j] = i1*10000 + i2;
562718f9 593
3edbbba2 594 }
fd30f88e 595
fd30f88e 596 // arranging cells within supercluster in decreasing order
2c1131dd 597
fd30f88e 598 for(j = 1;j <= ncl[i]; j++)
599 {
600 itest = 0;
601 ihld = iord[j];
602 for(i1 = 0; i1 < j; i1++)
603 {
604 if(itest == 0 && z[iord[i1]] < z[ihld])
605 {
606 itest = 1;
607 for(i2 = j-1; i2 >= i1; i2--)
608 {
609 iord[i2+1] = iord[i2];
610 }
611 iord[i1] = ihld;
612 }
613 }
3edbbba2 614 }
2c1131dd 615 /* MODIFICATION PART STARTS (Tapan July 2008)
616 iord[0] is the cell with highest ADC in the crude-cluster
617 ig is the number of local maxima in the crude-cluster
618 For the higest peak we make ig=0 which means first local maximum.
619 Next we go down in terms of the ADC sequence and find out if any
620 more of the cells form local maxima. The definition of local
621 maxima is that all its neighbours are of less ADC compared to it.
622 */
562718f9 623 ig = 0;
fd30f88e 624 xc[ig] = x[iord[0]];
625 yc[ig] = y[iord[0]];
626 zc[ig] = z[iord[0]];
2c1131dd 627 tc[ig] = t[iord[0]];
628 Int_t ivalid = 0, icount = 0;
629
630 for(j=1;j<=ncl[i];j++)
fd30f88e 631 {
2c1131dd 632 x1 = x[iord[j]];
633 y1 = y[iord[j]];
634 z1 = z[iord[j]];
635 t1 = t[iord[j]];
636 rr=Distance(x1,y1,xc[ig],yc[ig]);
637
638 // Check the cells which are outside the neighbours (rr>1.2)
639 if(rr>1.2 )
fd30f88e 640 {
2c1131dd 641 ivalid=0;
642 icount=0;
643 for(Int_t j1=1;j1<j;j1++)
644 {
645 icount++;
646 Float_t rr1=Distance(x1,y1,x[iord[j1]],y[iord[j1]]);
647 if(rr1>1.2) ivalid++;
648 }
649 if(ivalid == icount && z1>0.5*zc[ig])
650 {
651 ig++;
652 xc[ig]=x1;
653 yc[ig]=y1;
654 zc[ig]=z1;
655 tc[ig]=t1;
656 }
657 }
3edbbba2 658 }
2c1131dd 659
660 icl=icl+ig+1;
661
662 // We use simple Gaussian weighting. (Tapan Jan 2005)
fd30f88e 663 // compute the number of cells belonging to each cluster.
2c1131dd 664 // cell can be shared between several clusters
665 // in the ratio of cluster energy deposition
666 // To calculate:
667 // (1) number of cells belonging to a cluster (ig) and
668 // (2) total ADC of the cluster (ig)
669 // (3) x and y positions of the cluster
670
671
672 Int_t *cellCount;
673 Int_t **cellXY;
674
675 Int_t *status;
676 Double_t *totaladc, *totaladc2, *ncell,*weight;
677 Double_t *xclust, *yclust, *sigxclust, *sigyclust;
678 Double_t *ax, *ay, *ax2, *ay2;
679
680
681 status = new Int_t [ncl[i]+1];
682 cellXY = new Int_t *[ncl[i]+1];
683
562718f9 684 cellCount = new Int_t [ig+1];
2c1131dd 685 totaladc = new Double_t [ig+1];
686 totaladc2 = new Double_t [ig+1];
687 ncell = new Double_t [ig+1];
688 weight = new Double_t [ig+1];
689 xclust = new Double_t [ig+1];
690 yclust = new Double_t [ig+1];
691 sigxclust = new Double_t [ig+1];
692 sigyclust = new Double_t [ig+1];
693 ax = new Double_t [ig+1];
694 ay = new Double_t [ig+1];
695 ax2 = new Double_t [ig+1];
696 ay2 = new Double_t [ig+1];
697
698 for(j = 0; j < ncl[i]+1; j++)
01c4d84a 699 {
2c1131dd 700 status[j] = 0;
701 cellXY[j] = new Int_t[ig+1];
fd30f88e 702 }
2c1131dd 703 //initialization
704 for(Int_t kcl = 0; kcl < ig+1; kcl++)
705 {
706 cellCount[kcl] = 0;
707 totaladc[kcl] = 0.;
708 totaladc2[kcl] = 0.;
709 ncell[kcl] = 0.;
710 weight[kcl] = 0.;
711 xclust[kcl] = 0.;
712 yclust[kcl] = 0.;
713 sigxclust[kcl] = 0.;
714 sigyclust[kcl] = 0.;
715 ax[kcl] = 0.;
716 ay[kcl] = 0.;
717 ax2[kcl] = 0.;
718 ay2[kcl] = 0.;
719 for(j = 0; j < ncl[i]+1; j++)
720 {
721 cellXY[j][kcl] = 0;
722 }
723 }
724 Double_t sumweight, gweight;
562718f9 725
2c1131dd 726 for(j = 0;j <= ncl[i]; j++)
fd30f88e 727 {
2c1131dd 728 x1 = x[iord[j]];
729 y1 = y[iord[j]];
730 z1 = z[iord[j]];
731 t1 = t[iord[j]];
732
733 for(Int_t kcl=0; kcl<=ig; kcl++)
fd30f88e 734 {
2c1131dd 735 x2 = xc[kcl];
736 y2 = yc[kcl];
737 rr = Distance(x1,y1,x2,y2);
738 t2 = tc[kcl];
739
740 if(rr==0)
fd30f88e 741 {
2c1131dd 742 ncell[kcl] = 1.;
743 totaladc[kcl] = z1;
744 totaladc2[kcl] = pow(z1,2);
745 ax[kcl] = x1 * z1;
746 ay[kcl] = y1 * z1;
747 ax2[kcl]= 0.;
748 ay2[kcl]= 0.;
749 status[j] = 1;
fd30f88e 750 }
2c1131dd 751 }
752 }
753
754 for(j = 0; j <= ncl[i]; j++)
755 {
756 Int_t maxweight = 0;
757 Double_t max = 0.;
758
759 if(status[j] == 0)
760 {
761 x1 = x[iord[j]];
762 y1 = y[iord[j]];
763 z1 = z[iord[j]];
764 t1 = t[iord[j]];
765 sumweight = 0.;
766
767 for(Int_t kcl = 0; kcl <= ig; kcl++)
fd30f88e 768 {
2c1131dd 769 x2 = xc[kcl];
770 y2 = yc[kcl];
771 rr = Distance(x1,y1,x2,y2);
772 gweight = exp(-(rr*rr)/(2*(1.2*1.2)));
773 weight[kcl] = zc[kcl] * gweight;
774 sumweight = sumweight + weight[kcl];
775
776 if(weight[kcl] > max)
fd30f88e 777 {
2c1131dd 778 max = weight[kcl];
779 maxweight = kcl;
fd30f88e 780 }
781 }
2c1131dd 782
783 cellXY[cellCount[maxweight]][maxweight] = iord[j];
784
785 cellCount[maxweight]++;
786
787 for(Int_t kcl = 0; kcl <= ig; kcl++)
fd30f88e 788 {
2c1131dd 789 x2 = xc[kcl];
790 y2 = yc[kcl];
791 rr = Distance(x1,y1,x2,y2);
792 t2 = tc[kcl];
793 // For calculating weighted mean:
794 // Weighted_mean = (\sum w_i x_i) / (\sum w_i)
795
796 if(sumweight>0.)
fd30f88e 797 {
2c1131dd 798 totaladc[kcl] = totaladc[kcl] + z1*(weight[kcl]/sumweight);
799 ncell[kcl] = ncell[kcl] + (weight[kcl]/sumweight);
800 ax[kcl] = ax[kcl] + (x1 * z1 *weight[kcl]/sumweight);
801 ay[kcl] = ay[kcl] + (y1 * z1 *weight[kcl]/sumweight);
802
803 // For calculating weighted sigma:
804 // Wieghted_sigma= ((\sum w_i)/((\sum w_i)^2 - \sum (w_i^2))) \sum w_i (x_i - \Hat\mu)^2
805 // I assume here x1,y1 are cluster centers, otherwise life becomes too difficult
806 totaladc2[kcl] = totaladc2[kcl] + pow((z1*(weight[kcl]/sumweight)),2);
807 ax2[kcl] = ax2[kcl] + (z1 *weight[kcl]/sumweight) * pow((x1-x2),2);
808 ay2[kcl] = ay2[kcl] + (z1 *weight[kcl]/sumweight) * pow((y1-y2),2);
fd30f88e 809 }
810 }
811 }
812 }
813
2c1131dd 814 for(Int_t kcl = 0; kcl <= ig; kcl++)
fd30f88e 815 {
2c1131dd 816
817 if(totaladc[kcl]>0){
818 if(totaladc[kcl]>0.)xclust[kcl] = (ax[kcl])/ totaladc[kcl];
819 if(totaladc[kcl]>0.)yclust[kcl] = (ay[kcl])/ totaladc[kcl];
820
821 if(totaladc[kcl]>0.)sigxclust[kcl] = (totaladc[kcl]/(pow(totaladc[kcl],2)-totaladc2[kcl]))*ax2[kcl];
822 if(totaladc[kcl]>0.)sigyclust[kcl] = (totaladc[kcl]/(pow(totaladc[kcl],2)-totaladc2[kcl]))*ay2[kcl];
823 }
824
825 for(j = 0; j < cellCount[kcl]; j++) clno++;
826
827 if (clno >= 4608)
fd30f88e 828 {
2c1131dd 829 AliWarning("RefClust: Too many clusters! more than 4608");
fd30f88e 830 return;
831 }
2c1131dd 832 clusdata[0] = xclust[kcl];
833 clusdata[1] = yclust[kcl];
834 clusdata[2] = totaladc[kcl];
835 clusdata[3] = ncell[kcl];
836 if(sigxclust[kcl] > sigyclust[kcl])
fd30f88e 837 {
2c1131dd 838 clusdata[4] = pow(sigxclust[kcl],0.5);
839 clusdata[5] = pow(sigyclust[kcl],0.5);
fd30f88e 840 }
841 else
842 {
2c1131dd 843 clusdata[4] = pow(sigyclust[kcl],0.5);
844 clusdata[5] = pow(sigxclust[kcl],0.5);
fd30f88e 845 }
562718f9 846
2c1131dd 847 clxy[0] = tc[kcl];
562718f9 848
2c1131dd 849 Int_t Ncell=1;
850 for (Int_t ii = 0; ii < cellCount[kcl]; ii++)
562718f9 851 {
2c1131dd 852 Float_t xx = x[cellXY[ii][kcl]];
853 Float_t yy = y[cellXY[ii][kcl]];
854
855 rr=Distance(xclust[kcl],yclust[kcl],xx,yy);
856 // Natasha
857 // We store only the nearest and nest-nearest neighbours
858 if(rr<2.2)
859 {
860 clxy[Ncell] = t[cellXY[ii][kcl]];
861 Ncell++;
862 }
863 }
864
562718f9 865 pmdcludata = new AliPMDcludata(clusdata,clxy);
866 fPMDclucont->Add(pmdcludata);
01c4d84a 867 }
2c1131dd 868
869 delete [] iord;
870 delete [] tc;
871 delete [] t;
872 delete [] x;
873 delete [] y;
874 delete [] z;
875 delete [] xc;
876 delete [] yc;
877 delete [] zc;
878
562718f9 879 delete [] cellCount;
2c1131dd 880 for(Int_t jj = 0; jj < ncl[i]+1; jj++) delete [] cellXY[jj];
881
882 delete [] status;
883 delete [] totaladc;
884 delete [] totaladc2;
885 delete [] ncell;
886 delete [] xclust;
887 delete [] yclust;
888 delete [] sigxclust;
889 delete [] sigyclust;
890 delete [] ax;
891 delete [] ay;
892 delete [] ax2;
893 delete [] ay2;
894 delete [] weight;
fd30f88e 895 }
3edbbba2 896 }
2c1131dd 897 delete [] ncl;
898 delete [] clxy;
3edbbba2 899}
900// ------------------------------------------------------------------------ //
fd30f88e 901Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1,
902 Double_t x2, Double_t y2)
3edbbba2 903{
fd30f88e 904 return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
3edbbba2 905}
906// ------------------------------------------------------------------------ //
907void AliPMDClusteringV1::SetEdepCut(Float_t decut)
908{
909 fCutoff = decut;
910}
911// ------------------------------------------------------------------------ //