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