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