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