Corrected identifier for tpcCalibRaw
[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
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++)
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;
e6ba3040 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
d270ca46 273 if (ismn < 12)
920e13db 274 {
d270ca46 275 if ((irow >= 0 && irow < 96) && (icol >= 0 && icol < 48))
276 {
277 celldataTr[ihit] = celltrack[icol][irow];
278 celldataPid[ihit] = cellpid[icol][irow];
279 celldataAdc[ihit] = (Float_t) celladc[icol][irow];
280 }
281 else
282 {
283 celldataTr[ihit] = -1;
284 celldataPid[ihit] = -1;
285 celldataAdc[ihit] = -1;
286 }
920e13db 287 }
d270ca46 288 else if (ismn >= 12 && ismn < 24)
920e13db 289 {
d270ca46 290 if ((irow >= 0 && irow < 48) && (icol >= 0 && icol < 96))
291 {
292 celldataTr[ihit] = celltrack[irow][icol];
293 celldataPid[ihit] = cellpid[irow][icol];
294 celldataAdc[ihit] = (Float_t) celladc[irow][icol];
e6ba3040 295
d270ca46 296 }
297 else
298 {
299 celldataTr[ihit] = -1;
300 celldataPid[ihit] = -1;
301 celldataAdc[ihit] = -1;
302 }
920e13db 303 }
d270ca46 304
920e13db 305 }
e6ba3040 306
920e13db 307 pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY,
b0e4d1e1 308 celldataTr, celldataPid, celldataAdc);
3edbbba2 309 pmdcont->Add(pmdcl);
310 }
fd30f88e 311
91e6e2a0 312 fPMDclucont->Delete();
3edbbba2 313}
314// ------------------------------------------------------------------------ //
562718f9 315Int_t AliPMDClusteringV1::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
316 Int_t iord1[], Double_t edepcell[])
3edbbba2 317{
fd30f88e 318 // Does crude clustering
3edbbba2 319 // Finds out only the big patch by just searching the
320 // connected cells
321 //
2c1131dd 322 const Int_t kndim = 4609;
323 Int_t i,j,k,id1,id2,icl, numcell, clust[2][kndim];
fd30f88e 324 Int_t jd1,jd2, icell, cellcount;
325 static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
326
3edbbba2 327 AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
328
fd30f88e 329 for (j = 0; j < kNDIMX; j++)
330 {
331 for(k = 0; k < kNDIMY; k++)
332 {
333 fInfocl[0][j][k] = 0;
334 fInfocl[1][j][k] = 0;
335 }
336 }
337 for(i=0; i < kNMX; i++)
338 {
339 fInfcl[0][i] = -1;
562718f9 340
341 j = iord1[i];
342 id2 = j/kNDIMX;
343 id1 = j-id2*kNDIMX;
344
345 if(edepcell[j] <= cutoff)
fd30f88e 346 {
347 fInfocl[0][id1][id2] = -1;
348 }
3edbbba2 349 }
562718f9 350
3edbbba2 351 // ---------------------------------------------------------------
352 // crude clustering begins. Start with cell having largest adc
353 // count and loop over the cells in descending order of adc count
354 // ---------------------------------------------------------------
562718f9 355
fd30f88e 356 icl = -1;
357 cellcount = -1;
3edbbba2 358
fd30f88e 359 for(icell = 0; icell <= nmx1; icell++)
360 {
562718f9 361 j = iord1[icell];
362 id2 = j/kNDIMX;
363 id1 = j-id2*kNDIMX;
364
fd30f88e 365 if(fInfocl[0][id1][id2] == 0 )
366 {
367 icl++;
368 numcell = 0;
369 cellcount++;
370 fInfocl[0][id1][id2] = 1;
371 fInfocl[1][id1][id2] = icl;
372 fInfcl[0][cellcount] = icl;
373 fInfcl[1][cellcount] = id1;
374 fInfcl[2][cellcount] = id2;
375
376 clust[0][numcell] = id1;
377 clust[1][numcell] = id2;
378
2c1131dd 379 for(i = 1; i < kndim; i++)
fd30f88e 380 {
381 clust[0][i]=0;
382 }
383 // ---------------------------------------------------------------
384 // check for adc count in neib. cells. If ne 0 put it in this clust
385 // ---------------------------------------------------------------
386 for(i = 0; i < 6; i++)
387 {
388 jd1 = id1 + neibx[i];
389 jd2 = id2 + neiby[i];
390 if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
391 fInfocl[0][jd1][jd2] == 0)
392 {
393 numcell++;
394 fInfocl[0][jd1][jd2] = 2;
395 fInfocl[1][jd1][jd2] = icl;
396 clust[0][numcell] = jd1;
397 clust[1][numcell] = jd2;
398 cellcount++;
399 fInfcl[0][cellcount] = icl;
400 fInfcl[1][cellcount] = jd1;
401 fInfcl[2][cellcount] = jd2;
402 }
403 }
404 // ---------------------------------------------------------------
405 // check adc count for neighbour's neighbours recursively and
406 // if nonzero, add these to the cluster.
407 // ---------------------------------------------------------------
2c1131dd 408 for(i = 1; i < kndim;i++)
fd30f88e 409 {
410 if(clust[0][i] != 0)
411 {
412 id1 = clust[0][i];
413 id2 = clust[1][i];
414 for(j = 0; j < 6 ; j++)
415 {
416 jd1 = id1 + neibx[j];
417 jd2 = id2 + neiby[j];
418 if( (jd1 >= 0 && jd1 < kNDIMX) &&
419 (jd2 >= 0 && jd2 < kNDIMY) &&
420 fInfocl[0][jd1][jd2] == 0 )
421 {
422 fInfocl[0][jd1][jd2] = 2;
423 fInfocl[1][jd1][jd2] = icl;
424 numcell++;
425 clust[0][numcell] = jd1;
426 clust[1][numcell] = jd2;
427 cellcount++;
428 fInfcl[0][cellcount] = icl;
429 fInfcl[1][cellcount] = jd1;
430 fInfcl[2][cellcount] = jd2;
431 }
432 }
433 }
3edbbba2 434 }
3edbbba2 435 }
3edbbba2 436 }
3edbbba2 437 return cellcount;
438}
439// ------------------------------------------------------------------------ //
562718f9 440void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
3edbbba2 441{
442 // Does the refining of clusters
443 // Takes the big patch and does gaussian fitting and
444 // finds out the more refined clusters
445 //
22d01768 446
22d01768 447 AliPMDcludata *pmdcludata = 0;
2c1131dd 448
449 const Int_t kNmaxCell = 19; // # of cells surrounding a cluster center
450
451 Int_t ndim = incr + 1;
452
453 Int_t *ncl = 0x0;
454 Int_t *clxy = 0x0;
455 Int_t i12, i22;
456 Int_t i, j, k, i1, i2, id, icl, itest,ihld, ig, nsupcl,clno, t1, t2;
562718f9 457 Float_t clusdata[6];
2c1131dd 458 Double_t x1, y1, z1, x2, y2, z2, rr;
459
460 ncl = new Int_t [ndim];
461 clxy = new Int_t [kNmaxCell];
462
562718f9 463 // Initialisation
2c1131dd 464 for(i = 0; i<ndim; i++)
465 {
466 ncl[i] = -1;
562718f9 467 if (i < 6) clusdata[i] = 0.;
2c1131dd 468 if (i < kNmaxCell) clxy[i] = 0;
fd30f88e 469 }
562718f9 470
fd30f88e 471 // clno counts the final clusters
3edbbba2 472 // nsupcl = # of superclusters; ncl[i]= # of cells in supercluster i
473 // x, y and z store (x,y) coordinates of and energy deposited in a cell
474 // xc, yc store (x,y) coordinates of the cluster center
475 // zc stores the energy deposited in a cluster
476 // rc is cluster radius
fd30f88e 477
478 clno = -1;
3edbbba2 479 nsupcl = -1;
562718f9 480
fd30f88e 481 for(i = 0; i <= incr; i++)
482 {
483 if(fInfcl[0][i] != nsupcl)
01c4d84a 484 {
fd30f88e 485 nsupcl++;
01c4d84a 486 }
2c1131dd 487 if (nsupcl > ndim)
01c4d84a 488 {
fd30f88e 489 AliWarning("RefClust: Too many superclusters!");
2c1131dd 490 nsupcl = ndim;
fd30f88e 491 break;
01c4d84a 492 }
fd30f88e 493 ncl[nsupcl]++;
3edbbba2 494 }
fd30f88e 495
496 AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
497 incr+1,nsupcl+1));
498 id = -1;
499 icl = -1;
562718f9 500
fd30f88e 501 for(i = 0; i <= nsupcl; i++)
502 {
503 if(ncl[i] == 0)
01c4d84a 504 {
fd30f88e 505 id++;
506 icl++;
2c1131dd 507 if (clno >= 4608)
fd30f88e 508 {
2c1131dd 509 AliWarning("RefClust: Too many clusters! more than 4608");
fd30f88e 510 return;
511 }
512 clno++;
513 i1 = fInfcl[1][id];
514 i2 = fInfcl[2][id];
515
78fc1b96 516 i12 = i1 + i2*kNDIMX;
562718f9 517
fd30f88e 518 clusdata[0] = fCoord[0][i1][i2];
519 clusdata[1] = fCoord[1][i1][i2];
562718f9 520 clusdata[2] = edepcell[i12];
fd30f88e 521 clusdata[3] = 1.;
522 clusdata[4] = 0.5;
523 clusdata[5] = 0.0;
562718f9 524
2c1131dd 525 clxy[0] = i1*10000 + i2;
526
527 for(Int_t icltr = 1; icltr < kNmaxCell; icltr++)
fd30f88e 528 {
562718f9 529 clxy[icltr] = -1;
3edbbba2 530 }
e6ba3040 531
562718f9 532 pmdcludata = new AliPMDcludata(clusdata,clxy);
533 fPMDclucont->Add(pmdcludata);
3edbbba2 534 }
fd30f88e 535 else if(ncl[i] == 1)
536 {
537 id++;
538 icl++;
2c1131dd 539 if (clno >= 4608)
fd30f88e 540 {
2c1131dd 541 AliWarning("RefClust: Too many clusters! more than 4608");
fd30f88e 542 return;
543 }
544 clno++;
545 i1 = fInfcl[1][id];
546 i2 = fInfcl[2][id];
78fc1b96 547 i12 = i1 + i2*kNDIMX;
562718f9 548
fd30f88e 549 x1 = fCoord[0][i1][i2];
550 y1 = fCoord[1][i1][i2];
562718f9 551 z1 = edepcell[i12];
2c1131dd 552
553 clxy[0] = i1*10000 + i2;
554
562718f9 555 id++;
fd30f88e 556 i1 = fInfcl[1][id];
557 i2 = fInfcl[2][id];
562718f9 558
2c1131dd 559 i22 = i1 + i2*kNDIMX;
fd30f88e 560 x2 = fCoord[0][i1][i2];
561 y2 = fCoord[1][i1][i2];
562718f9 562 z2 = edepcell[i22];
2c1131dd 563
564 clxy[1] = i1*10000 + i2;
565
566
567 for(Int_t icltr = 2; icltr < kNmaxCell; icltr++)
fd30f88e 568 {
562718f9 569 clxy[icltr] = -1;
fd30f88e 570 }
fd30f88e 571
572 clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
573 clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
574 clusdata[2] = z1+z2;
575 clusdata[3] = 2.;
576 clusdata[4] = 0.5;
577 clusdata[5] = 0.0;
562718f9 578 pmdcludata = new AliPMDcludata(clusdata,clxy);
579 fPMDclucont->Add(pmdcludata);
3edbbba2 580 }
fd30f88e 581 else
582 {
2c1131dd 583
584 Int_t *iord, *tc, *t;
585 Double_t *x, *y, *z, *xc, *yc, *zc;
586
587 iord = new Int_t [ncl[i]+1];
588 tc = new Int_t [ncl[i]+1];
589 t = new Int_t [ncl[i]+1];
590
591 x = new Double_t [ncl[i]+1];
592 y = new Double_t [ncl[i]+1];
593 z = new Double_t [ncl[i]+1];
594 xc = new Double_t [ncl[i]+1];
595 yc = new Double_t [ncl[i]+1];
596 zc = new Double_t [ncl[i]+1];
597
598 for( k = 0; k < ncl[i]+1; k++)
599 {
600 iord[k] = -1;
601 t[k] = -1;
602 tc[k] = -1;
603 x[k] = -1;
604 y[k] = -1;
605 z[k] = -1;
606 xc[k] = -1;
607 yc[k] = -1;
608 zc[k] = -1;
609 }
fd30f88e 610 id++;
fd30f88e 611 // super-cluster of more than two cells - broken up into smaller
612 // clusters gaussian centers computed. (peaks separated by > 1 cell)
613 // Begin from cell having largest energy deposited This is first
614 // cluster center
615 i1 = fInfcl[1][id];
616 i2 = fInfcl[2][id];
78fc1b96 617 i12 = i1 + i2*kNDIMX;
562718f9 618
fd30f88e 619 x[0] = fCoord[0][i1][i2];
620 y[0] = fCoord[1][i1][i2];
562718f9 621 z[0] = edepcell[i12];
2c1131dd 622 t[0] = i1*10000 + i2;
fd30f88e 623
2c1131dd 624
fd30f88e 625 iord[0] = 0;
626 for(j = 1; j <= ncl[i]; j++)
627 {
628 id++;
629 i1 = fInfcl[1][id];
630 i2 = fInfcl[2][id];
78fc1b96 631 i12 = i1 + i2*kNDIMX;
562718f9 632
fd30f88e 633 iord[j] = j;
634 x[j] = fCoord[0][i1][i2];
635 y[j] = fCoord[1][i1][i2];
562718f9 636 z[j] = edepcell[i12];
2c1131dd 637 t[j] = i1*10000 + i2;
562718f9 638
3edbbba2 639 }
fd30f88e 640
fd30f88e 641 // arranging cells within supercluster in decreasing order
2c1131dd 642
fd30f88e 643 for(j = 1;j <= ncl[i]; j++)
644 {
645 itest = 0;
646 ihld = iord[j];
647 for(i1 = 0; i1 < j; i1++)
648 {
649 if(itest == 0 && z[iord[i1]] < z[ihld])
650 {
651 itest = 1;
652 for(i2 = j-1; i2 >= i1; i2--)
653 {
654 iord[i2+1] = iord[i2];
655 }
656 iord[i1] = ihld;
657 }
658 }
3edbbba2 659 }
2c1131dd 660 /* MODIFICATION PART STARTS (Tapan July 2008)
661 iord[0] is the cell with highest ADC in the crude-cluster
662 ig is the number of local maxima in the crude-cluster
663 For the higest peak we make ig=0 which means first local maximum.
664 Next we go down in terms of the ADC sequence and find out if any
665 more of the cells form local maxima. The definition of local
666 maxima is that all its neighbours are of less ADC compared to it.
667 */
562718f9 668 ig = 0;
fd30f88e 669 xc[ig] = x[iord[0]];
670 yc[ig] = y[iord[0]];
671 zc[ig] = z[iord[0]];
2c1131dd 672 tc[ig] = t[iord[0]];
673 Int_t ivalid = 0, icount = 0;
674
675 for(j=1;j<=ncl[i];j++)
fd30f88e 676 {
2c1131dd 677 x1 = x[iord[j]];
678 y1 = y[iord[j]];
679 z1 = z[iord[j]];
680 t1 = t[iord[j]];
681 rr=Distance(x1,y1,xc[ig],yc[ig]);
682
683 // Check the cells which are outside the neighbours (rr>1.2)
684 if(rr>1.2 )
fd30f88e 685 {
2c1131dd 686 ivalid=0;
687 icount=0;
688 for(Int_t j1=1;j1<j;j1++)
689 {
690 icount++;
691 Float_t rr1=Distance(x1,y1,x[iord[j1]],y[iord[j1]]);
692 if(rr1>1.2) ivalid++;
693 }
694 if(ivalid == icount && z1>0.5*zc[ig])
695 {
696 ig++;
697 xc[ig]=x1;
698 yc[ig]=y1;
699 zc[ig]=z1;
700 tc[ig]=t1;
701 }
702 }
3edbbba2 703 }
2c1131dd 704
705 icl=icl+ig+1;
706
707 // We use simple Gaussian weighting. (Tapan Jan 2005)
fd30f88e 708 // compute the number of cells belonging to each cluster.
2c1131dd 709 // cell can be shared between several clusters
710 // in the ratio of cluster energy deposition
711 // To calculate:
712 // (1) number of cells belonging to a cluster (ig) and
713 // (2) total ADC of the cluster (ig)
714 // (3) x and y positions of the cluster
715
716
717 Int_t *cellCount;
718 Int_t **cellXY;
719
720 Int_t *status;
721 Double_t *totaladc, *totaladc2, *ncell,*weight;
722 Double_t *xclust, *yclust, *sigxclust, *sigyclust;
723 Double_t *ax, *ay, *ax2, *ay2;
724
725
726 status = new Int_t [ncl[i]+1];
727 cellXY = new Int_t *[ncl[i]+1];
728
562718f9 729 cellCount = new Int_t [ig+1];
2c1131dd 730 totaladc = new Double_t [ig+1];
731 totaladc2 = new Double_t [ig+1];
732 ncell = new Double_t [ig+1];
733 weight = new Double_t [ig+1];
734 xclust = new Double_t [ig+1];
735 yclust = new Double_t [ig+1];
736 sigxclust = new Double_t [ig+1];
737 sigyclust = new Double_t [ig+1];
738 ax = new Double_t [ig+1];
739 ay = new Double_t [ig+1];
740 ax2 = new Double_t [ig+1];
741 ay2 = new Double_t [ig+1];
742
743 for(j = 0; j < ncl[i]+1; j++)
01c4d84a 744 {
2c1131dd 745 status[j] = 0;
746 cellXY[j] = new Int_t[ig+1];
fd30f88e 747 }
2c1131dd 748 //initialization
749 for(Int_t kcl = 0; kcl < ig+1; kcl++)
750 {
751 cellCount[kcl] = 0;
752 totaladc[kcl] = 0.;
753 totaladc2[kcl] = 0.;
754 ncell[kcl] = 0.;
755 weight[kcl] = 0.;
756 xclust[kcl] = 0.;
757 yclust[kcl] = 0.;
758 sigxclust[kcl] = 0.;
759 sigyclust[kcl] = 0.;
760 ax[kcl] = 0.;
761 ay[kcl] = 0.;
762 ax2[kcl] = 0.;
763 ay2[kcl] = 0.;
764 for(j = 0; j < ncl[i]+1; j++)
765 {
766 cellXY[j][kcl] = 0;
767 }
768 }
769 Double_t sumweight, gweight;
562718f9 770
2c1131dd 771 for(j = 0;j <= ncl[i]; j++)
fd30f88e 772 {
2c1131dd 773 x1 = x[iord[j]];
774 y1 = y[iord[j]];
775 z1 = z[iord[j]];
776 t1 = t[iord[j]];
777
778 for(Int_t kcl=0; kcl<=ig; kcl++)
fd30f88e 779 {
2c1131dd 780 x2 = xc[kcl];
781 y2 = yc[kcl];
782 rr = Distance(x1,y1,x2,y2);
783 t2 = tc[kcl];
784
785 if(rr==0)
fd30f88e 786 {
2c1131dd 787 ncell[kcl] = 1.;
788 totaladc[kcl] = z1;
789 totaladc2[kcl] = pow(z1,2);
790 ax[kcl] = x1 * z1;
791 ay[kcl] = y1 * z1;
792 ax2[kcl]= 0.;
793 ay2[kcl]= 0.;
794 status[j] = 1;
fd30f88e 795 }
2c1131dd 796 }
797 }
798
799 for(j = 0; j <= ncl[i]; j++)
800 {
801 Int_t maxweight = 0;
802 Double_t max = 0.;
803
804 if(status[j] == 0)
805 {
806 x1 = x[iord[j]];
807 y1 = y[iord[j]];
808 z1 = z[iord[j]];
809 t1 = t[iord[j]];
810 sumweight = 0.;
811
812 for(Int_t kcl = 0; kcl <= ig; kcl++)
fd30f88e 813 {
2c1131dd 814 x2 = xc[kcl];
815 y2 = yc[kcl];
816 rr = Distance(x1,y1,x2,y2);
817 gweight = exp(-(rr*rr)/(2*(1.2*1.2)));
818 weight[kcl] = zc[kcl] * gweight;
819 sumweight = sumweight + weight[kcl];
820
821 if(weight[kcl] > max)
fd30f88e 822 {
2c1131dd 823 max = weight[kcl];
824 maxweight = kcl;
fd30f88e 825 }
826 }
2c1131dd 827
828 cellXY[cellCount[maxweight]][maxweight] = iord[j];
829
830 cellCount[maxweight]++;
831
d270ca46 832 x2 = xc[maxweight];
833 y2 = yc[maxweight];
834 totaladc[maxweight] += z1;
835 ax[maxweight] += x1 * z1;
836 ay[maxweight] += y1 * z1;
837 totaladc2[maxweight] += pow(z1,2);
838 ax2[maxweight] += z1 * pow((x1-x2),2);
839 ay2[maxweight] += z1 * pow((y1-y2),2);
840 ncell[maxweight]++;
841
fd30f88e 842 }
843 }
844
2c1131dd 845 for(Int_t kcl = 0; kcl <= ig; kcl++)
fd30f88e 846 {
2c1131dd 847
848 if(totaladc[kcl]>0){
849 if(totaladc[kcl]>0.)xclust[kcl] = (ax[kcl])/ totaladc[kcl];
850 if(totaladc[kcl]>0.)yclust[kcl] = (ay[kcl])/ totaladc[kcl];
851
e6ba3040 852 //natasha
853 if(totaladc2[kcl] >= pow(totaladc[kcl],2))
854 {
855 sigxclust[kcl] = 0.25;
856 sigyclust[kcl] = 0.25;
857 }
858 else
859 {
860 sigxclust[kcl] = (totaladc[kcl]/(pow(totaladc[kcl],2)-totaladc2[kcl]))*ax2[kcl];
861 sigyclust[kcl] = (totaladc[kcl]/(pow(totaladc[kcl],2)-totaladc2[kcl]))*ay2[kcl];
862 }
863 }
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];
e6ba3040 876
877
2c1131dd 878 if(sigxclust[kcl] > sigyclust[kcl])
fd30f88e 879 {
2c1131dd 880 clusdata[4] = pow(sigxclust[kcl],0.5);
881 clusdata[5] = pow(sigyclust[kcl],0.5);
fd30f88e 882 }
883 else
884 {
2c1131dd 885 clusdata[4] = pow(sigyclust[kcl],0.5);
886 clusdata[5] = pow(sigxclust[kcl],0.5);
fd30f88e 887 }
e6ba3040 888
2c1131dd 889 clxy[0] = tc[kcl];
562718f9 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 }
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];
970 Float_t adc = (Float_t) celladc[id1][jd1];
971 if(adc == 0.)
972 {
973 isocount++;
974 if(isocount == kCellNeighbour)
975 {
976 Float_t cadc = (Float_t) celladc[irow][icol];
977
978 isocell = new AliPMDisocell(idet,ismn,irow,icol,cadc);
979 pmdisocell->Add(isocell);
980
981 }
982 }
983 } // neigh cell cond.
984 }
985 }
986 }
987
988
989}
990// ------------------------------------------------------------------------ //
3edbbba2 991void AliPMDClusteringV1::SetEdepCut(Float_t decut)
992{
993 fCutoff = decut;
994}
995// ------------------------------------------------------------------------ //