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