bug fixed for alignment, removed alignment database access from AliPMDUtility class
[u/mrichter/AliRoot.git] / PMD / AliPMDClusteringV2.cxx
CommitLineData
8c7250c5 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
8c7250c5 18//-----------------------------------------------------//
19// //
20// Source File : PMDClusteringV2.cxx //
21// //
22// clustering code for alice pmd //
23// //
24//-----------------------------------------------------//
25
26/* --------------------------------------------------------------------
27 Code developed by S. C. Phatak, Institute of Physics,
28 Bhubaneswar 751 005 ( phatak@iopb.res.in ) Given the energy deposited
29 ( or ADC value ) in each cell of supermodule ( pmd or cpv ), the code
30 builds up superclusters and breaks them into clusters. The input is
562718f9 31 in TObjarray and cluster information is in TObjArray.
32 integer clno gives total number of clusters in the supermodule.
33 fClusters is the global ( public ) variables.
8c7250c5 34 Others are local ( private ) to the code.
35 At the moment, the data is read for whole detector ( all supermodules
36 and pmd as well as cpv. This will have to be modify later )
37 LAST UPDATE : October 23, 2002
38-----------------------------------------------------------------------*/
39
090026bf 40#include <Riostream.h>
41#include <TMath.h>
8c7250c5 42#include <TObjArray.h>
562718f9 43#include <TArrayI.h>
8c7250c5 44
562718f9 45#include "AliPMDcludata.h"
8c7250c5 46#include "AliPMDcluster.h"
47#include "AliPMDClustering.h"
48#include "AliPMDClusteringV2.h"
49#include "AliLog.h"
50
51ClassImp(AliPMDClusteringV2)
52
53const Double_t AliPMDClusteringV2::fgkSqroot3by2=0.8660254; // sqrt(3.)/2.
54
55AliPMDClusteringV2::AliPMDClusteringV2():
562718f9 56 fPMDclucont(new TObjArray()),
4755c3f1 57 fCutoff(0.0),
58 fClusParam(0)
8c7250c5 59{
60 for(int i = 0; i < kNDIMX; i++)
61 {
62 for(int j = 0; j < kNDIMY; j++)
63 {
64 fCoord[0][i][j] = i+j/2.;
65 fCoord[1][i][j] = fgkSqroot3by2*j;
8c7250c5 66 }
67 }
68}
69// ------------------------------------------------------------------------ //
562718f9 70
71
72AliPMDClusteringV2::AliPMDClusteringV2(const AliPMDClusteringV2& pmdclv2):
73 AliPMDClustering(pmdclv2),
74 fPMDclucont(0),
4755c3f1 75 fCutoff(0),
76 fClusParam(0)
562718f9 77{
78 // copy constructor
79 AliError("Copy constructor not allowed ");
80
81}
82// ------------------------------------------------------------------------ //
83AliPMDClusteringV2 &AliPMDClusteringV2::operator=(const AliPMDClusteringV2& /*pmdclv2*/)
84{
85 // copy constructor
86 AliError("Assignment operator not allowed ");
87 return *this;
88}
89// ------------------------------------------------------------------------ //
8c7250c5 90AliPMDClusteringV2::~AliPMDClusteringV2()
91{
562718f9 92 delete fPMDclucont;
8c7250c5 93}
94// ------------------------------------------------------------------------ //
562718f9 95
96void AliPMDClusteringV2::DoClust(Int_t idet, Int_t ismn,
920e13db 97 Int_t celltrack[48][96],
98 Int_t cellpid[48][96],
99 Double_t celladc[48][96],
22bd512d 100 TObjArray *pmdcont)
8c7250c5 101{
102 // main function to call other necessary functions to do clustering
103 //
104 AliPMDcluster *pmdcl = 0;
105
2c1131dd 106 const Float_t ktwobysqrt3 = 1.1547; // 2./sqrt(3.)
107 const Int_t kNmaxCell = 19; // # of cells surrounding a cluster center
be8b7039 108 Int_t i = 0, j = 0, nmx1 = 0;
109 Int_t incr = 0, id = 0, jd = 0;
562718f9 110 Int_t ndimXr = 0;
111 Int_t ndimYr = 0;
2c1131dd 112 Int_t celldataX[kNmaxCell], celldataY[kNmaxCell];
920e13db 113 Int_t celldataTr[kNmaxCell], celldataPid[kNmaxCell];
b0e4d1e1 114 Float_t celldataAdc[kNmaxCell];
be8b7039 115 Float_t clusdata[6] = {0.,0.,0.,0.,0.,0.};
116 Double_t cutoff = 0., ave = 0.;
562718f9 117 Double_t edepcell[kNMX];
8c7250c5 118
8c7250c5 119
120 if (ismn < 12)
121 {
122 ndimXr = 96;
123 ndimYr = 48;
124 }
125 else if (ismn >= 12 && ismn <= 23)
126 {
127 ndimXr = 48;
128 ndimYr = 96;
129 }
562718f9 130
78fc1b96 131 for (i =0; i < kNMX; i++)
8c7250c5 132 {
562718f9 133 edepcell[i] = 0.;
8c7250c5 134 }
562718f9 135
8c7250c5 136 for (id = 0; id < ndimXr; id++)
137 {
138 for (jd = 0; jd < ndimYr; jd++)
139 {
562718f9 140 j = jd;
141 i = id + (ndimYr/2-1) - (jd/2);
142 Int_t ij = i + j*kNDIMX;
8c7250c5 143 if (ismn < 12)
144 {
562718f9 145 edepcell[ij] = celladc[jd][id];
8c7250c5 146 }
147 else if (ismn >= 12 && ismn <= 23)
148 {
562718f9 149 edepcell[ij] = celladc[id][jd];
8c7250c5 150 }
151
152 }
153 }
154
562718f9 155 Int_t iord1[kNMX];
df4e6759 156 TMath::Sort((Int_t)kNMX,edepcell,iord1);// order the data
8c7250c5 157 cutoff = fCutoff; // cutoff used to discard cells having ener. dep.
562718f9 158 ave = 0.;
159 nmx1 = -1;
8c7250c5 160
562718f9 161 for(i = 0;i < kNMX; i++)
8c7250c5 162 {
562718f9 163 if(edepcell[i] > 0.)
164 {
165 ave += edepcell[i];
166 }
167 if(edepcell[i] > cutoff )
168 {
169 nmx1++;
170 }
8c7250c5 171 }
562718f9 172
8c7250c5 173 AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
562718f9 174
175 if (nmx1 == 0)
176 {
177 nmx1 = 1;
178 }
179 ave = ave/nmx1;
180
8c7250c5 181 AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
182 kNMX,ave));
8c7250c5 183
562718f9 184 incr = CrClust(ave, cutoff, nmx1,iord1, edepcell);
185 RefClust(incr,edepcell );
186
187 Int_t nentries1 = fPMDclucont->GetEntries();
188 AliDebug(1,Form("Detector Plane = %d Serial Module No = %d Number of clusters = %d",idet, ismn, nentries1));
189 AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
190 for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
8c7250c5 191 {
562718f9 192 AliPMDcludata *pmdcludata =
193 (AliPMDcludata*)fPMDclucont->UncheckedAt(ient1);
194 Float_t cluXC = pmdcludata->GetClusX();
195 Float_t cluYC = pmdcludata->GetClusY();
196 Float_t cluADC = pmdcludata->GetClusADC();
197 Float_t cluCELLS = pmdcludata->GetClusCells();
198 Float_t cluSIGX = pmdcludata->GetClusSigmaX();
199 Float_t cluSIGY = pmdcludata->GetClusSigmaY();
200
8c7250c5 201 Float_t cluY0 = ktwobysqrt3*cluYC;
202 Float_t cluX0 = cluXC - cluY0/2.;
562718f9 203
8c7250c5 204 //
205 // Cluster X centroid is back transformed
206 //
207 if (ismn < 12)
208 {
209 clusdata[0] = cluX0 - (24-1) + cluY0/2.;
210 }
562718f9 211 else if (ismn >= 12 && ismn <= 23)
8c7250c5 212 {
213 clusdata[0] = cluX0 - (48-1) + cluY0/2.;
214 }
215
562718f9 216 clusdata[1] = cluY0;
217 clusdata[2] = cluADC;
218 clusdata[3] = cluCELLS;
219 clusdata[4] = cluSIGX;
220 clusdata[5] = cluSIGY;
8c7250c5 221 //
222 // Cells associated with a cluster
223 //
2c1131dd 224 for (Int_t ihit = 0; ihit < kNmaxCell; ihit++)
8c7250c5 225 {
562718f9 226 Int_t dummyXY = pmdcludata->GetCellXY(ihit);
227
228 Int_t celldumY = dummyXY%10000;
229 Int_t celldumX = dummyXY/10000;
c1339151 230 Float_t cellY = (Float_t) celldumY/10;
231 Float_t cellX = (Float_t) celldumX/10;
232
562718f9 233 //
234 // Cell X centroid is back transformed
235 //
236 if (ismn < 12)
237 {
c1339151 238 celldataX[ihit] = (Int_t) ((cellX - (24-1) + cellY/2.) + 0.5);
562718f9 239 }
240 else if (ismn >= 12 && ismn <= 23)
241 {
c1339151 242 celldataX[ihit] = (Int_t) ((cellX - (48-1) + cellY/2.) + 0.5 );
562718f9 243 }
b0e4d1e1 244 celldataY[ihit] = (Int_t) (cellY + 0.5);
b6d6a9b5 245
246 Int_t irow = celldataX[ihit];
247 Int_t icol = celldataY[ihit];
248
249 if ((irow >= 0 && irow < 48) && (icol >= 0 && icol < 96))
250 {
251 celldataTr[ihit] = celltrack[irow][icol];
252 celldataPid[ihit] = cellpid[irow][icol];
253 celldataAdc[ihit] = (Float_t) celladc[irow][icol];
254 }
255 else
256 {
257 celldataTr[ihit] = -1;
258 celldataPid[ihit] = -1;
259 celldataAdc[ihit] = -1;
260 }
261
8c7250c5 262 }
263
920e13db 264 pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY,
b0e4d1e1 265 celldataTr, celldataPid, celldataAdc);
8c7250c5 266 pmdcont->Add(pmdcl);
267 }
2c1131dd 268 fPMDclucont->Delete();
8c7250c5 269}
270// ------------------------------------------------------------------------ //
562718f9 271Int_t AliPMDClusteringV2::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
272 Int_t iord1[], Double_t edepcell[])
8c7250c5 273{
274 // Does crude clustering
275 // Finds out only the big patch by just searching the
276 // connected cells
277 //
278
be8b7039 279 Int_t i = 0, j = 0, k = 0, id1 =0, id2 = 0, icl = 0, numcell = 0;
280 Int_t jd1 = 0, jd2 = 0, icell = 0, cellcount = 0;
562718f9 281 Int_t clust[2][5000];
282 static Int_t neibx[6] = {1,0,-1,-1,0,1}, neiby[6] = {0,1,1,0,-1,-1};
8c7250c5 283
284 // neibx and neiby define ( incremental ) (i,j) for the neighbours of a
285 // cell. There are six neighbours.
286 // cellcount --- total number of cells having nonzero ener dep
287 // numcell --- number of cells in a given supercluster
562718f9 288
8c7250c5 289 AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
290
562718f9 291 for (j=0; j < kNDIMX; j++)
292 {
293 for(k=0; k < kNDIMY; k++)
294 {
295 fInfocl[0][j][k] = 0;
296 fInfocl[1][j][k] = 0;
297 }
298 }
299
300 for(i=0; i < kNMX; i++)
301 {
302 fInfcl[0][i] = -1;
303
304 j = iord1[i];
305 id2 = j/kNDIMX;
306 id1 = j-id2*kNDIMX;
307
308 if(edepcell[j] <= cutoff)
309 {
310 fInfocl[0][id1][id2] = -1;
311 }
8c7250c5 312 }
8c7250c5 313 // ---------------------------------------------------------------
314 // crude clustering begins. Start with cell having largest adc
315 // count and loop over the cells in descending order of adc count
316 // ---------------------------------------------------------------
562718f9 317 icl = -1;
318 cellcount = -1;
319 for(icell=0; icell <= nmx1; icell++)
320 {
321 j = iord1[icell];
322 id2 = j/kNDIMX;
323 id1 = j-id2*kNDIMX;
324 if(fInfocl[0][id1][id2] == 0 )
325 {
326 // ---------------------------------------------------------------
327 // icl -- cluster #, numcell -- # of cells in it, clust -- stores
328 // coordinates of the cells in a cluster, fInfocl[0][i1][i2] is 1 for
329 // primary and 2 for secondary cells,
330 // fInfocl[1][i1][i2] stores cluster #
331 // ---------------------------------------------------------------
332 icl++;
333 numcell = 0;
334 cellcount++;
335 fInfocl[0][id1][id2] = 1;
336 fInfocl[1][id1][id2] = icl;
337 fInfcl[0][cellcount] = icl;
338 fInfcl[1][cellcount] = id1;
339 fInfcl[2][cellcount] = id2;
340
341 clust[0][numcell] = id1;
342 clust[1][numcell] = id2;
343 for(i = 1; i < 5000; i++)
344 {
345 clust[0][i] = -1;
346 }
347 // ---------------------------------------------------------------
348 // check for adc count in neib. cells. If ne 0 put it in this clust
349 // ---------------------------------------------------------------
350 for(i = 0; i < 6; i++)
351 {
352 jd1 = id1 + neibx[i];
353 jd2 = id2 + neiby[i];
8c7250c5 354 if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
562718f9 355 fInfocl[0][jd1][jd2] == 0)
356 {
357 numcell++;
358 fInfocl[0][jd1][jd2] = 2;
359 fInfocl[1][jd1][jd2] = icl;
360 clust[0][numcell] = jd1;
361 clust[1][numcell] = jd2;
362 cellcount++;
363 fInfcl[0][cellcount] = icl;
364 fInfcl[1][cellcount] = jd1;
365 fInfcl[2][cellcount] = jd2;
366 }
367 }
368 // ---------------------------------------------------------------
369 // check adc count for neighbour's neighbours recursively and
370 // if nonzero, add these to the cluster.
371 // ---------------------------------------------------------------
372 for(i = 1;i < 5000; i++)
373 {
374 if(clust[0][i] != -1)
375 {
376 id1 = clust[0][i];
377 id2 = clust[1][i];
378 for(j = 0; j < 6 ; j++)
379 {
380 jd1 = id1 + neibx[j];
381 jd2 = id2 + neiby[j];
382 if( (jd1 >= 0 && jd1 < kNDIMX) &&
383 (jd2 >= 0 && jd2 < kNDIMY)
384 && fInfocl[0][jd1][jd2] == 0 )
385 {
386 fInfocl[0][jd1][jd2] = 2;
387 fInfocl[1][jd1][jd2] = icl;
388 numcell++;
389 clust[0][numcell] = jd1;
390 clust[1][numcell] = jd2;
391 cellcount++;
392 fInfcl[0][cellcount] = icl;
393 fInfcl[1][cellcount] = jd1;
394 fInfcl[2][cellcount] = jd2;
395 }
396 }
397 }
8c7250c5 398 }
8c7250c5 399 }
8c7250c5 400 }
8c7250c5 401 return cellcount;
402}
403// ------------------------------------------------------------------------ //
562718f9 404 void AliPMDClusteringV2::RefClust(Int_t incr, Double_t edepcell[])
8c7250c5 405{
406 // Does the refining of clusters
407 // Takes the big patch and does gaussian fitting and
408 // finds out the more refined clusters
8c7250c5 409
c1339151 410 const Float_t ktwobysqrt3 = 1.1547;
2c1131dd 411 const Int_t kNmaxCell = 19;
412
562718f9 413 AliPMDcludata *pmdcludata = 0;
c1339151 414
be8b7039 415 Int_t i12 = 0;
416 Int_t i = 0, j = 0, k = 0;
417 Int_t i1 = 0, i2 = 0, id = 0, icl = 0, itest = 0, ihld = 0;
418 Int_t ig = 0, nsupcl = 0, clno = 0, clX = 0, clY = 0;
419 Int_t clxy[kNmaxCell];
c1339151 420
be8b7039 421 Float_t clusdata[6] = {0.,0.,0.,0.,0.,0.};
422 Double_t x1 = 0., y1 = 0., z1 = 0., x2 = 0., y2 = 0., z2 = 0., rr = 0.;
c1339151 423
424 Int_t kndim = incr + 1;
425
426 TArrayI testncl;
427 TArrayI testindex;
428
429 Int_t *ncl, *iord;
430
431 Double_t *x, *y, *z, *xc, *yc, *zc, *cells, *rcl, *rcs;
432
433 ncl = new Int_t [kndim];
434 iord = new Int_t [kndim];
435 x = new Double_t [kndim];
436 y = new Double_t [kndim];
437 z = new Double_t [kndim];
438 xc = new Double_t [kndim];
439 yc = new Double_t [kndim];
440 zc = new Double_t [kndim];
441 cells = new Double_t [kndim];
442 rcl = new Double_t [kndim];
443 rcs = new Double_t [kndim];
562718f9 444
445 for(Int_t kk = 0; kk < 15; kk++)
446 {
447 if( kk < 6 )clusdata[kk] = 0.;
448 }
449
8c7250c5 450 // nsupcl = # of superclusters; ncl[i]= # of cells in supercluster i
451 // x, y and z store (x,y) coordinates of and energy deposited in a cell
452 // xc, yc store (x,y) coordinates of the cluster center
562718f9 453 // zc stores the energy deposited in a cluster, rc is cluster radius
8c7250c5 454
562718f9 455 clno = -1;
8c7250c5 456 nsupcl = -1;
c1339151 457
458 for(i = 0; i < kndim; i++)
562718f9 459 {
460 ncl[i] = -1;
8c7250c5 461 }
c1339151 462 for(i = 0; i <= incr; i++)
562718f9 463 {
464 if(fInfcl[0][i] != nsupcl)
465 {
466 nsupcl++;
467 }
468 if (nsupcl > 4500)
469 {
470 AliWarning("RefClust: Too many superclusters!");
471 nsupcl = 4500;
472 break;
473 }
474 ncl[nsupcl]++;
475 }
476
8c7250c5 477 AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
478 incr+1,nsupcl+1));
562718f9 479
480 id = -1;
481 icl = -1;
c1339151 482 for(i = 0; i <= nsupcl; i++)
562718f9 483 {
484 if(ncl[i] == 0)
8c7250c5 485 {
562718f9 486 id++;
487 icl++;
488 // one cell super-clusters --> single cluster
489 // cluster center at the centyer of the cell
490 // cluster radius = half cell dimension
491 if (clno >= 5000)
8c7250c5 492 {
562718f9 493 AliWarning("RefClust: Too many clusters! more than 5000");
494 return;
8c7250c5 495 }
562718f9 496 clno++;
497 i1 = fInfcl[1][id];
498 i2 = fInfcl[2][id];
78fc1b96 499 i12 = i1 + i2*kNDIMX;
562718f9 500 clusdata[0] = fCoord[0][i1][i2];
501 clusdata[1] = fCoord[1][i1][i2];
502 clusdata[2] = edepcell[i12];
503 clusdata[3] = 1.;
504 clusdata[4] = 0.0;
505 clusdata[5] = 0.0;
506
507 //cell information
c1339151 508
509 clY = (Int_t)((ktwobysqrt3*fCoord[1][i1][i2])*10);
510 clX = (Int_t)((fCoord[0][i1][i2] - clY/20.)*10);
511 clxy[0] = clX*10000 + clY ;
512
2c1131dd 513 for(Int_t icltr = 1; icltr < kNmaxCell; icltr++)
562718f9 514 {
515 clxy[icltr] = -1;
516 }
517 pmdcludata = new AliPMDcludata(clusdata,clxy);
518 fPMDclucont->Add(pmdcludata);
519
520
8c7250c5 521 }
562718f9 522 else if(ncl[i] == 1)
8c7250c5 523 {
562718f9 524 // two cell super-cluster --> single cluster
525 // cluster center is at ener. dep.-weighted mean of two cells
526 // cluster radius == half cell dimension
527 id++;
528 icl++;
529 if (clno >= 5000)
8c7250c5 530 {
531 AliWarning("RefClust: Too many clusters! more than 5000");
532 return;
533 }
562718f9 534 clno++;
535 i1 = fInfcl[1][id];
536 i2 = fInfcl[2][id];
78fc1b96 537 i12 = i1 + i2*kNDIMX;
562718f9 538
539 x1 = fCoord[0][i1][i2];
540 y1 = fCoord[1][i1][i2];
541 z1 = edepcell[i12];
542
543 id++;
544 i1 = fInfcl[1][id];
545 i2 = fInfcl[2][id];
546 i12 = i1 + i2*kNDIMX;
547
548 x2 = fCoord[0][i1][i2];
549 y2 = fCoord[1][i1][i2];
550 z2 = edepcell[i12];
551
552 clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
553 clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
554 clusdata[2] = z1+z2;
555 clusdata[3] = 2.;
556 clusdata[4] = (TMath::Sqrt(z1*z2))/(z1+z2);
557 clusdata[5] = 0.0;
558
c1339151 559 clY = (Int_t)((ktwobysqrt3*y1)*10);
560 clX = (Int_t)((x1 - clY/20.)*10);
561 clxy[0] = clX*10000 + clY ;
562
c1339151 563 clY = (Int_t)((ktwobysqrt3*y2)*10);
564 clX = (Int_t)((x2 - clY/20.)*10);
565 clxy[1] = clX*10000 + clY ;
566
2c1131dd 567 for(Int_t icltr = 2; icltr < kNmaxCell; icltr++)
8c7250c5 568 {
562718f9 569 clxy[icltr] = -1;
8c7250c5 570 }
562718f9 571 pmdcludata = new AliPMDcludata(clusdata, clxy);
572 fPMDclucont->Add(pmdcludata);
8c7250c5 573 }
562718f9 574 else{
575 id++;
576 iord[0] = 0;
577 // super-cluster of more than two cells - broken up into smaller
578 // clusters gaussian centers computed. (peaks separated by > 1 cell)
579 // Begin from cell having largest energy deposited This is first
580 // cluster center
581 // *****************************************************************
582 // NOTE --- POSSIBLE MODIFICATION: ONE MAY NOT BREAKING SUPERCLUSTERS
583 // IF NO. OF CELLS IS NOT TOO LARGE ( SAY 5 OR 6 )
584 // SINCE WE EXPECT THE SUPERCLUSTER
585 // TO BE A SINGLE CLUSTER
586 //*******************************************************************
587
588 i1 = fInfcl[1][id];
589 i2 = fInfcl[2][id];
78fc1b96 590 i12 = i1 + i2*kNDIMX;
562718f9 591
592 x[0] = fCoord[0][i1][i2];
593 y[0] = fCoord[1][i1][i2];
594 z[0] = edepcell[i12];
595
596 iord[0] = 0;
597 for(j = 1; j <= ncl[i]; j++)
598 {
599
600 id++;
601 i1 = fInfcl[1][id];
602 i2 = fInfcl[2][id];
78fc1b96 603 i12 = i1 + i2*kNDIMX;
562718f9 604 iord[j] = j;
605 x[j] = fCoord[0][i1][i2];
606 y[j] = fCoord[1][i1][i2];
607 z[j] = edepcell[i12];
608 }
609
610 // arranging cells within supercluster in decreasing order
611 for(j = 1; j <= ncl[i];j++)
612 {
613 itest = 0;
614 ihld = iord[j];
615 for(i1 = 0; i1 < j; i1++)
616 {
617 if(itest == 0 && z[iord[i1]] < z[ihld])
618 {
619 itest = 1;
620 for(i2 = j-1;i2 >= i1;i2--)
621 {
622 iord[i2+1] = iord[i2];
623 }
624 iord[i1] = ihld;
625 }
626 }
627 }
628
629
630 // compute the number of clusters and their centers ( first
631 // guess )
632 // centers must be separated by cells having smaller ener. dep.
633 // neighbouring centers should be either strong or well-separated
634 ig = 0;
635 xc[ig] = x[iord[0]];
636 yc[ig] = y[iord[0]];
637 zc[ig] = z[iord[0]];
638 for(j = 1; j <= ncl[i]; j++)
639 {
640 itest = -1;
641 x1 = x[iord[j]];
642 y1 = y[iord[j]];
643 for(k = 0; k <= ig; k++)
644 {
645 x2 = xc[k];
646 y2 = yc[k];
647 rr = Distance(x1,y1,x2,y2);
648 //************************************************************
649 // finetuning cluster splitting
650 // the numbers zc/4 and zc/10 may need to be changed.
651 // Also one may need to add one more layer because our
652 // cells are smaller in absolute scale
653 //************************************************************
654
655
656 if( rr >= 1.1 && rr < 1.8 && z[iord[j]] > zc[k]/4.) itest++;
657 if( rr >= 1.8 && rr < 2.1 && z[iord[j]] > zc[k]/10.) itest++;
658 if( rr >= 2.1)itest++;
659 }
660
661 if(itest == ig)
662 {
663 ig++;
664 xc[ig] = x1;
665 yc[ig] = y1;
666 zc[ig] = z[iord[j]];
667 }
668 }
c1339151 669 ClustDetails(ncl[i], ig, x, y ,z, xc, yc, zc, rcl, rcs, cells,
670 testncl, testindex);
562718f9 671
672 Int_t pp = 0;
673 for(j = 0; j <= ig; j++)
674 {
675 clno++;
676 if (clno >= 5000)
677 {
678 AliWarning("RefClust: Too many clusters! more than 5000");
679 return;
680 }
681 clusdata[0] = xc[j];
682 clusdata[1] = yc[j];
683 clusdata[2] = zc[j];
684 clusdata[4] = rcl[j];
685 clusdata[5] = rcs[j];
686 if(ig == 0)
687 {
2c1131dd 688 clusdata[3] = ncl[i] + 1;
562718f9 689 }
690 else
691 {
692 clusdata[3] = cells[j];
693 }
694 // cell information
695 Int_t ncellcls = testncl[j];
2c1131dd 696 if( ncellcls < kNmaxCell )
562718f9 697 {
698 for(Int_t kk = 1; kk <= ncellcls; kk++)
699 {
700 Int_t ll = testindex[pp];
c1339151 701 clY = (Int_t)((ktwobysqrt3*y[ll])*10);
702 clX = (Int_t)((x[ll] - clY/20.)*10);
703 clxy[kk-1] = clX*10000 + clY ;
c1339151 704
562718f9 705 pp++;
706 }
2c1131dd 707 for(Int_t icltr = ncellcls ; icltr < kNmaxCell; icltr++)
562718f9 708 {
709 clxy[icltr] = -1;
710 }
711 }
712 pmdcludata = new AliPMDcludata(clusdata, clxy);
713 fPMDclucont->Add(pmdcludata);
714 }
715 testncl.Set(0);
716 testindex.Set(0);
717 }
8c7250c5 718 }
c1339151 719 delete [] ncl;
720 delete [] iord;
721 delete [] x;
722 delete [] y;
723 delete [] z;
724 delete [] xc;
725 delete [] yc;
726 delete [] zc;
727 delete [] cells;
728 delete [] rcl;
729 delete [] rcs;
8c7250c5 730}
8c7250c5 731// ------------------------------------------------------------------------ //
c1339151 732void AliPMDClusteringV2::ClustDetails(Int_t ncell, Int_t nclust, Double_t x[],
733 Double_t y[], Double_t z[],Double_t xc[],
734 Double_t yc[], Double_t zc[],
735 Double_t rcl[], Double_t rcs[],
736 Double_t cells[], TArrayI &testncl,
562718f9 737 TArrayI &testindex)
8c7250c5 738{
739 // function begins
740 //
8c7250c5 741
c1339151 742 Int_t kndim1 = ncell + 1;//ncell
743 Int_t kndim2 = 20;
744 Int_t kndim3 = nclust + 1;//nclust
78fc1b96 745
be8b7039 746 Int_t i = 0, j = 0, k = 0, i1 = 0, i2 = 0;
747 Double_t x1 = 0., y1 = 0., x2 = 0., y2 = 0.;
748 Double_t rr = 0., b = 0., c = 0., r1 = 0., r2 = 0.;
749 Double_t sumx = 0., sumy = 0., sumxy = 0.;
750 Double_t sumxx = 0., sum = 0., sum1 = 0., sumyy = 0.;
8c7250c5 751
c1339151 752 Double_t *str, *str1, *xcl, *ycl, *cln;
753 Int_t **cell;
754 Int_t ** cluster;
755 Double_t **clustcell;
756 str = new Double_t [kndim3];
757 str1 = new Double_t [kndim3];
758 xcl = new Double_t [kndim3];
759 ycl = new Double_t [kndim3];
760 cln = new Double_t [kndim3];
562718f9 761
c1339151 762 clustcell = new Double_t *[kndim3];
763 cell = new Int_t *[kndim3];
764 cluster = new Int_t *[kndim1];
562718f9 765 for(i = 0; i < kndim1; i++)
766 {
c1339151 767 cluster[i] = new Int_t [kndim2];
768 }
769
770 for(i = 0; i < kndim3; i++)
771 {
772 str[i] = 0;
773 str1[i] = 0;
774 xcl[i] = 0;
775 ycl[i] = 0;
776 cln[i] = 0;
777
778 cell[i] = new Int_t [kndim2];
779 clustcell[i] = new Double_t [kndim1];
780 for(j = 0; j < kndim1; j++)
8c7250c5 781 {
c1339151 782 clustcell[i][j] = 0;
562718f9 783 }
c1339151 784 for(j = 0; j < kndim2; j++)
562718f9 785 {
c1339151 786 cluster[i][j] = 0;
787 cell[i][j] = 0;
8c7250c5 788 }
562718f9 789 }
c1339151 790
562718f9 791 if(nclust > 0)
792 {
793 // more than one cluster
794 // checking cells shared between several clusters.
795 // First check if the cell is within
796 // one cell unit ( nearest neighbour). Else,
797 // if it is within 1.74 cell units ( next nearest )
798 // Else if it is upto 2 cell units etc.
799
800 for (i = 0; i <= ncell; i++)
8c7250c5 801 {
c1339151 802 x1 = x[i];
803 y1 = y[i];
562718f9 804 cluster[i][0] = 0;
2c1131dd 805
562718f9 806 // distance <= 1 cell unit
2c1131dd 807
562718f9 808 for(j = 0; j <= nclust; j++)
8c7250c5 809 {
c1339151 810 x2 = xc[j];
811 y2 = yc[j];
8c7250c5 812 rr = Distance(x1, y1, x2, y2);
562718f9 813 if(rr <= 1.)
8c7250c5 814 {
815 cluster[i][0]++;
816 i1 = cluster[i][0];
817 cluster[i][i1] = j;
818 }
819 }
562718f9 820 // next nearest neighbour
821 if(cluster[i][0] == 0)
822 {
823 for(j=0; j<=nclust; j++)
824 {
c1339151 825 x2 = xc[j];
826 y2 = yc[j];
562718f9 827 rr = Distance(x1, y1, x2, y2);
828 if(rr <= TMath::Sqrt(3.))
829 {
830 cluster[i][0]++;
831 i1 = cluster[i][0];
832 cluster[i][i1] = j;
833 }
834 }
835 }
836 // next-to-next nearest neighbour
837 if(cluster[i][0] == 0)
838 {
839 for(j=0; j<=nclust; j++)
840 {
c1339151 841 x2 = xc[j];
842 y2 = yc[j];
562718f9 843 rr = Distance(x1, y1, x2, y2);
844 if(rr <= 2.)
845 {
846 cluster[i][0]++;
847 i1 = cluster[i][0];
848 cluster[i][i1] = j;
849 }
850 }
851 }
852 // one more
853 if(cluster[i][0] == 0)
854 {
855 for(j = 0; j <= nclust; j++)
856 {
c1339151 857 x2 = xc[j];
858 y2 = yc[j];
562718f9 859 rr = Distance(x1, y1, x2, y2);
860 if(rr <= 2.7)
861 {
862 cluster[i][0]++;
863 i1 = cluster[i][0];
864 cluster[i][i1] = j;
865 }
866 }
867 }
8c7250c5 868 }
562718f9 869
870 // computing cluster strength. Some cells are shared.
871 for(i = 0; i <= ncell; i++)
8c7250c5 872 {
562718f9 873 if(cluster[i][0] != 0)
8c7250c5 874 {
562718f9 875 i1 = cluster[i][0];
876 for(j = 1; j <= i1; j++)
8c7250c5 877 {
562718f9 878 i2 = cluster[i][j];
c1339151 879 str[i2] += z[i]/i1;
8c7250c5 880 }
881 }
882 }
562718f9 883
884 for(k = 0; k < 5; k++)
8c7250c5 885 {
562718f9 886 for(i = 0; i <= ncell; i++)
8c7250c5 887 {
562718f9 888 if(cluster[i][0] != 0)
8c7250c5 889 {
562718f9 890 i1=cluster[i][0];
891 sum=0.;
892 for(j = 1; j <= i1; j++)
893 {
894 sum += str[cluster[i][j]];
895 }
896
897 for(j = 1; j <= i1; j++)
898 {
899 i2 = cluster[i][j];
c1339151 900 str1[i2] += z[i]*str[i2]/sum;
901 clustcell[i2][i] = z[i]*str[i2]/sum;
562718f9 902 }
8c7250c5 903 }
904 }
562718f9 905
906
907 for(j = 0; j <= nclust; j++)
908 {
909 str[j] = str1[j];
910 str1[j] = 0.;
911 }
8c7250c5 912 }
562718f9 913
914 for(i = 0; i <= nclust; i++)
915 {
916 sumx = 0.;
917 sumy = 0.;
918 sum = 0.;
919 sum1 = 0.;
920 for(j = 0; j <= ncell; j++)
921 {
922 if(clustcell[i][j] != 0)
923 {
c1339151 924 sumx += clustcell[i][j]*x[j];
925 sumy += clustcell[i][j]*y[j];
562718f9 926 sum += clustcell[i][j];
c1339151 927 sum1 += clustcell[i][j]/z[j];
562718f9 928 }
929 }
930 //** xcl and ycl are cluster centroid positions ( center of gravity )
931
932 xcl[i] = sumx/sum;
933 ycl[i] = sumy/sum;
934 cln[i] = sum1;
935 sumxx = 0.;
936 sumyy = 0.;
937 sumxy = 0.;
938 for(j = 0; j <= ncell; j++)
939 {
c1339151 940 sumxx += clustcell[i][j]*(x[j]-xcl[i])*(x[j]-xcl[i])/sum;
941 sumyy += clustcell[i][j]*(y[j]-ycl[i])*(y[j]-ycl[i])/sum;
942 sumxy += clustcell[i][j]*(x[j]-xcl[i])*(y[j]-ycl[i])/sum;
562718f9 943 }
944 b = sumxx+sumyy;
945 c = sumxx*sumyy-sumxy*sumxy;
946 // ******************r1 and r2 are major and minor axes ( r1 > r2 ).
947 r1 = b/2.+TMath::Sqrt(b*b/4.-c);
948 r2 = b/2.-TMath::Sqrt(b*b/4.-c);
949 // final assignments to proper external variables
c1339151 950 xc[i] = xcl[i];
951 yc[i] = ycl[i];
952 zc[i] = str[i];
953 cells[i] = cln[i];
954 rcl[i] = r1;
955 rcs[i] = r2;
956
8c7250c5 957 }
562718f9 958
959 //To get the cell position in a cluster
960
961 for(Int_t ii=0; ii<= ncell; ii++)
962 {
963 Int_t jj = cluster[ii][0];
964 for(Int_t kk=1; kk<= jj; kk++)
8c7250c5 965 {
562718f9 966 Int_t ll = cluster[ii][kk];
967 cell[ll][0]++;
968 cell[ll][cell[ll][0]] = ii;
8c7250c5 969 }
562718f9 970 }
971
972 testncl.Set(nclust+1);
973 Int_t counter = 0;
974
975 for(Int_t ii=0; ii <= nclust; ii++)
976 {
977 testncl[ii] = cell[ii][0];
978 counter += testncl[ii];
979 }
980 testindex.Set(counter);
981 Int_t ll = 0;
982 for(Int_t ii=0; ii<= nclust; ii++)
983 {
984 for(Int_t jj = 1; jj<= testncl[ii]; jj++)
985 {
986 Int_t kk = cell[ii][jj];
987 testindex[ll] = kk;
988 ll++;
989 }
990 }
991
992 }
c1339151 993 else if(nclust == 0)
562718f9 994 {
8c7250c5 995 sumx = 0.;
996 sumy = 0.;
997 sum = 0.;
998 sum1 = 0.;
562718f9 999 i = 0;
1000 for(j = 0; j <= ncell; j++)
1001 {
c1339151 1002 sumx += z[j]*x[j];
1003 sumy += z[j]*y[j];
1004 sum += z[j];
562718f9 1005 sum1++;
8c7250c5 1006 }
8c7250c5 1007 xcl[i] = sumx/sum;
1008 ycl[i] = sumy/sum;
1009 cln[i] = sum1;
562718f9 1010 sumxx = 0.;
1011 sumyy = 0.;
1012 sumxy = 0.;
1013 for(j = 0; j <= ncell; j++)
1014 {
c1339151 1015 sumxx += clustcell[i][j]*(x[j]-xcl[i])*(x[j]-xcl[i])/sum;
1016 sumyy += clustcell[i][j]*(y[j]-ycl[i])*(y[j]-ycl[i])/sum;
1017 sumxy += clustcell[i][j]*(x[j]-xcl[i])*(y[j]-ycl[i])/sum;
562718f9 1018 }
1019 b = sumxx+sumyy;
1020 c = sumxx*sumyy-sumxy*sumxy;
1021 r1 = b/2.+ TMath::Sqrt(b*b/4.-c);
1022 r2 = b/2.- TMath::Sqrt(b*b/4.-c);
1023
1024 // To get the cell position in a cluster
1025 testncl.Set(nclust+1);
c1339151 1026 testindex.Set(ncell+1);
2c1131dd 1027 cell[0][0] = ncell + 1;
562718f9 1028 testncl[0] = cell[0][0];
1029 Int_t ll = 0;
c1339151 1030 for(Int_t ii = 1; ii <= ncell; ii++)
562718f9 1031 {
1032 cell[0][ii]=ii;
562718f9 1033 Int_t kk = cell[0][ii];
1034 testindex[ll] = kk;
1035 ll++;
1036 }
1037 // final assignments
c1339151 1038 xc[i] = xcl[i];
1039 yc[i] = ycl[i];
c1339151 1040 zc[i] = sum;
1041 cells[i] = cln[i];
1042 rcl[i] = r1;
1043 rcs[i] = r2;
1044 }
1045 for(i = 0; i < kndim3; i++)
1046 {
1047 delete [] clustcell[i];
1048 delete [] cell[i];
1049 }
1050 delete [] clustcell;
1051 delete [] cell;
1052 for(i = 0; i <kndim1 ; i++)
1053 {
1054 delete [] cluster[i];
8c7250c5 1055 }
c1339151 1056 delete [] cluster;
1057 delete [] str;
1058 delete [] str1;
1059 delete [] xcl;
1060 delete [] ycl;
1061 delete [] cln;
8c7250c5 1062}
1063
1064// ------------------------------------------------------------------------ //
1065Double_t AliPMDClusteringV2::Distance(Double_t x1, Double_t y1,
1066 Double_t x2, Double_t y2)
1067{
562718f9 1068 return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
8c7250c5 1069}
1070// ------------------------------------------------------------------------ //
1071void AliPMDClusteringV2::SetEdepCut(Float_t decut)
1072{
1073 fCutoff = decut;
1074}
1075// ------------------------------------------------------------------------ //
4755c3f1 1076void AliPMDClusteringV2::SetClusteringParam(Int_t cluspar)
1077{
1078 fClusParam = cluspar;
1079}
1080// ------------------------------------------------------------------------ //