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