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