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 : PMDClusteringV1.cxx, Version 00 //
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
33 : in array edepcell[kNMX] and cluster information is in a
34 : TObjarray. Integer clno gives total number of clusters in the
35 : supermodule.
36 :
37 : fClusters is the only global ( public ) variables.
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 :
44 : #include <Riostream.h>
45 : #include <TMath.h>
46 : #include <TNtuple.h>
47 : #include <TObjArray.h>
48 : #include "TRandom.h"
49 : #include <stdio.h>
50 :
51 : #include "AliPMDcludata.h"
52 : #include "AliPMDcluster.h"
53 : #include "AliPMDClustering.h"
54 : #include "AliPMDClusteringV1.h"
55 : #include "AliLog.h"
56 :
57 :
58 12 : ClassImp(AliPMDClusteringV1)
59 :
60 : const Double_t AliPMDClusteringV1::fgkSqroot3by2=0.8660254; // sqrt(3.)/2.
61 :
62 8 : AliPMDClusteringV1::AliPMDClusteringV1():
63 24 : fPMDclucont(new TObjArray()),
64 8 : fCutoff(0.0),
65 8 : fClusParam(0)
66 40 : {
67 1920 : for(Int_t i = 0; i < kNDIMX; i++)
68 : {
69 184688 : for(Int_t j = 0; j < kNDIMY; j++)
70 : {
71 91392 : fCoord[0][i][j] = i+j/2.;
72 91392 : fCoord[1][i][j] = fgkSqroot3by2*j;
73 : }
74 : }
75 16 : }
76 : // ------------------------------------------------------------------------ //
77 : AliPMDClusteringV1::AliPMDClusteringV1(const AliPMDClusteringV1& pmdclv1):
78 0 : AliPMDClustering(pmdclv1),
79 0 : fPMDclucont(0),
80 0 : fCutoff(0),
81 0 : fClusParam(0)
82 0 : {
83 : // copy constructor
84 0 : AliError("Copy constructor not allowed ");
85 :
86 0 : }
87 : // ------------------------------------------------------------------------ //
88 : AliPMDClusteringV1 &AliPMDClusteringV1::operator=(const AliPMDClusteringV1& /*pmdclv1*/)
89 : {
90 : // copy constructor
91 0 : AliError("Assignment operator not allowed ");
92 0 : return *this;
93 : }
94 : // ------------------------------------------------------------------------ //
95 : AliPMDClusteringV1::~AliPMDClusteringV1()
96 48 : {
97 16 : delete fPMDclucont;
98 24 : }
99 : // ------------------------------------------------------------------------ //
100 : void 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],
104 : TObjArray *pmdcont)
105 : {
106 : // main function to call other necessary functions to do clustering
107 : //
108 :
109 : AliPMDcluster *pmdcl = 0;
110 :
111 : const float ktwobysqrt3 = 1.1547; // 2./sqrt(3.)
112 : const Int_t kNmaxCell = 19; // # of cells surrounding a cluster center
113 :
114 : Int_t i = 0, j = 0, nmx1 = 0;
115 : Int_t incr = 0, id = 0, jd = 0;
116 196 : Int_t celldataX[kNmaxCell], celldataY[kNmaxCell];
117 98 : Int_t celldataTr[kNmaxCell], celldataPid[kNmaxCell];
118 98 : Float_t celldataAdc[kNmaxCell];
119 98 : Float_t clusdata[6] = {0.,0.,0.,0.,0.,0.};
120 : Double_t cutoff, ave;
121 98 : Double_t edepcell[kNMX];
122 98 : Double_t cellenergy[kNMX];
123 :
124 : // ndimXr and ndimYr are different because of different module size
125 :
126 : Int_t ndimXr = 0;
127 : Int_t ndimYr = 0;
128 :
129 98 : if (ismn < 12)
130 : {
131 : ndimXr = 96;
132 : ndimYr = 48;
133 62 : }
134 36 : else if (ismn >= 12 && ismn <= 23)
135 : {
136 : ndimXr = 48;
137 : ndimYr = 96;
138 36 : }
139 :
140 2239300 : for (i = 0; i < kNMX; i++)
141 : {
142 1119552 : cellenergy[i] = 0.;
143 : }
144 :
145 : Int_t kk = 0;
146 23520 : for (i = 0; i < kNDIMX; i++)
147 : {
148 2262428 : for (j = 0; j < kNDIMY; j++)
149 : {
150 1119552 : edepcell[kk] = 0.;
151 1119552 : kk++;
152 : }
153 : }
154 :
155 15556 : for (id = 0; id < ndimXr; id++)
156 : {
157 918528 : for (jd = 0; jd < ndimYr; jd++)
158 : {
159 : j = jd;
160 451584 : i = id+(ndimYr/2-1)-(jd/2);
161 :
162 451584 : Int_t ij = i + j*kNDIMX;
163 :
164 451584 : if (ismn < 12)
165 : {
166 285696 : cellenergy[ij] = celladc[jd][id];
167 285696 : }
168 165888 : else if (ismn >= 12 && ismn <= 23)
169 : {
170 165888 : cellenergy[ij] = celladc[id][jd];
171 165888 : }
172 : }
173 : }
174 :
175 2239300 : for (i = 0; i < kNMX; i++)
176 : {
177 1119552 : edepcell[i] = cellenergy[i];
178 : }
179 :
180 : Bool_t jsort = true;
181 : // the dimension of iord1 is increased twice
182 98 : Int_t iord1[2*kNMX];
183 98 : TMath::Sort((Int_t)kNMX,edepcell,iord1,jsort);// order the data
184 98 : cutoff = fCutoff; // cutoff to discard cells
185 : ave = 0.;
186 : nmx1 = -1;
187 2239300 : for(i = 0;i < kNMX; i++)
188 : {
189 1119552 : if(edepcell[i] > 0.)
190 : {
191 360 : ave += edepcell[i];
192 360 : }
193 1119552 : if(edepcell[i] > cutoff )
194 : {
195 352 : nmx1++;
196 352 : }
197 : }
198 :
199 294 : AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
200 :
201 146 : if (nmx1 == 0) nmx1 = 1;
202 98 : ave = ave/nmx1;
203 294 : AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
204 : kNMX,ave));
205 :
206 98 : incr = CrClust(ave, cutoff, nmx1,iord1, edepcell );
207 98 : RefClust(incr,edepcell);
208 98 : Int_t nentries1 = fPMDclucont->GetEntries();
209 294 : AliDebug(1,Form("Detector Plane = %d Serial Module No = %d Number of clusters = %d",idet, ismn, nentries1));
210 294 : AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
211 :
212 540 : for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
213 : {
214 : AliPMDcludata *pmdcludata =
215 172 : (AliPMDcludata*)fPMDclucont->UncheckedAt(ient1);
216 172 : Float_t cluXC = pmdcludata->GetClusX();
217 172 : Float_t cluYC = pmdcludata->GetClusY();
218 172 : Float_t cluADC = pmdcludata->GetClusADC();
219 172 : Float_t cluCELLS = pmdcludata->GetClusCells();
220 172 : Float_t cluSIGX = pmdcludata->GetClusSigmaX();
221 172 : Float_t cluSIGY = pmdcludata->GetClusSigmaY();
222 :
223 172 : Float_t cluY0 = ktwobysqrt3*cluYC;
224 172 : Float_t cluX0 = cluXC - cluY0/2.;
225 :
226 : //
227 : // Cluster X centroid is back transformed
228 : //
229 :
230 172 : if (ismn < 12)
231 : {
232 128 : clusdata[0] = cluX0 - (24-1) + cluY0/2.;
233 128 : }
234 44 : else if ( ismn >= 12 && ismn <= 23)
235 : {
236 44 : clusdata[0] = cluX0 - (48-1) + cluY0/2.;
237 44 : }
238 :
239 172 : clusdata[1] = cluY0;
240 172 : clusdata[2] = cluADC;
241 172 : clusdata[3] = cluCELLS;
242 172 : clusdata[4] = cluSIGX;
243 172 : clusdata[5] = cluSIGY;
244 :
245 : //
246 : // Cells associated with a cluster
247 : //
248 :
249 6880 : for (Int_t ihit = 0; ihit < kNmaxCell; ihit++)
250 : {
251 3268 : Int_t cellrow = pmdcludata->GetCellXY(ihit)/10000;
252 3268 : Int_t cellcol = pmdcludata->GetCellXY(ihit)%10000;
253 :
254 3268 : if (ismn < 12)
255 : {
256 2432 : celldataX[ihit] = cellrow - (24-1) + int(cellcol/2.);
257 2432 : }
258 836 : else if (ismn >= 12 && ismn <= 23)
259 : {
260 836 : celldataX[ihit] = cellrow - (48-1) + int(cellcol/2.);
261 836 : }
262 :
263 3268 : celldataY[ihit] = cellcol;
264 :
265 3268 : Int_t irow = celldataX[ihit];
266 : Int_t icol = celldataY[ihit];
267 :
268 3268 : if (ismn < 12)
269 : {
270 2432 : if ((irow >= 0 && irow < 96) && (icol >= 0 && icol < 48))
271 : {
272 274 : celldataTr[ihit] = celltrack[icol][irow];
273 274 : celldataPid[ihit] = cellpid[icol][irow];
274 274 : celldataAdc[ihit] = (Float_t) celladc[icol][irow];
275 274 : }
276 : else
277 : {
278 2158 : celldataTr[ihit] = -1;
279 2158 : celldataPid[ihit] = -1;
280 2158 : celldataAdc[ihit] = -1;
281 : }
282 : }
283 836 : else if (ismn >= 12 && ismn < 24)
284 : {
285 836 : if ((irow >= 0 && irow < 48) && (icol >= 0 && icol < 96))
286 : {
287 78 : celldataTr[ihit] = celltrack[irow][icol];
288 78 : celldataPid[ihit] = cellpid[irow][icol];
289 78 : celldataAdc[ihit] = (Float_t) celladc[irow][icol];
290 :
291 78 : }
292 : else
293 : {
294 758 : celldataTr[ihit] = -1;
295 758 : celldataPid[ihit] = -1;
296 758 : celldataAdc[ihit] = -1;
297 : }
298 : }
299 :
300 : }
301 :
302 344 : pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY,
303 172 : celldataTr, celldataPid, celldataAdc);
304 172 : pmdcont->Add(pmdcl);
305 : }
306 :
307 98 : fPMDclucont->Delete();
308 98 : }
309 : // ------------------------------------------------------------------------ //
310 : Int_t AliPMDClusteringV1::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
311 : Int_t iord1[], Double_t edepcell[])
312 : {
313 : // Does crude clustering
314 : // Finds out only the big patch by just searching the
315 : // connected cells
316 : //
317 : const Int_t kndim = 4609;
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 196 : Int_t clust[2][kndim];
321 : static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
322 :
323 294 : AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
324 :
325 23520 : for (j = 0; j < kNDIMX; j++)
326 : {
327 2262428 : for(k = 0; k < kNDIMY; k++)
328 : {
329 1119552 : fInfocl[0][j][k] = 0;
330 1119552 : fInfocl[1][j][k] = 0;
331 : }
332 : }
333 2239300 : for(i=0; i < kNMX; i++)
334 : {
335 1119552 : fInfcl[0][i] = -1;
336 :
337 1119552 : j = iord1[i];
338 1119552 : id2 = j/kNDIMX;
339 1119552 : id1 = j-id2*kNDIMX;
340 :
341 1119552 : if(edepcell[j] <= cutoff)
342 : {
343 1119200 : fInfocl[0][id1][id2] = -1;
344 1119200 : }
345 : }
346 :
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 : // ---------------------------------------------------------------
351 :
352 : icl = -1;
353 : cellcount = -1;
354 :
355 996 : for(icell = 0; icell <= nmx1; icell++)
356 : {
357 400 : j = iord1[icell];
358 400 : id2 = j/kNDIMX;
359 400 : id1 = j-id2*kNDIMX;
360 :
361 400 : if(fInfocl[0][id1][id2] == 0 )
362 : {
363 172 : icl++;
364 : numcell = 0;
365 172 : cellcount++;
366 172 : fInfocl[0][id1][id2] = 1;
367 172 : fInfocl[1][id1][id2] = icl;
368 172 : fInfcl[0][cellcount] = icl;
369 172 : fInfcl[1][cellcount] = id1;
370 172 : fInfcl[2][cellcount] = id2;
371 :
372 172 : clust[0][numcell] = id1;
373 172 : clust[1][numcell] = id2;
374 :
375 1585496 : for(i = 1; i < kndim; i++)
376 : {
377 792576 : clust[0][i]=0;
378 : }
379 : // ---------------------------------------------------------------
380 : // check for adc count in neib. cells. If ne 0 put it in this clust
381 : // ---------------------------------------------------------------
382 2408 : for(i = 0; i < 6; i++)
383 : {
384 1032 : jd1 = id1 + neibx[i];
385 1032 : jd2 = id2 + neiby[i];
386 2064 : if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
387 1032 : fInfocl[0][jd1][jd2] == 0)
388 : {
389 92 : numcell++;
390 92 : fInfocl[0][jd1][jd2] = 2;
391 92 : fInfocl[1][jd1][jd2] = icl;
392 92 : clust[0][numcell] = jd1;
393 92 : clust[1][numcell] = jd2;
394 92 : cellcount++;
395 92 : fInfcl[0][cellcount] = icl;
396 92 : fInfcl[1][cellcount] = jd1;
397 92 : fInfcl[2][cellcount] = jd2;
398 92 : }
399 : }
400 : // ---------------------------------------------------------------
401 : // check adc count for neighbour's neighbours recursively and
402 : // if nonzero, add these to the cluster.
403 : // ---------------------------------------------------------------
404 1585496 : for(i = 1; i < kndim;i++)
405 : {
406 792576 : if(clust[0][i] != 0)
407 : {
408 : id1 = clust[0][i];
409 180 : id2 = clust[1][i];
410 2520 : for(j = 0; j < 6 ; j++)
411 : {
412 1080 : jd1 = id1 + neibx[j];
413 1080 : jd2 = id2 + neiby[j];
414 4320 : if( (jd1 >= 0 && jd1 < kNDIMX) &&
415 2160 : (jd2 >= 0 && jd2 < kNDIMY) &&
416 1080 : fInfocl[0][jd1][jd2] == 0 )
417 : {
418 88 : fInfocl[0][jd1][jd2] = 2;
419 88 : fInfocl[1][jd1][jd2] = icl;
420 88 : numcell++;
421 88 : clust[0][numcell] = jd1;
422 88 : clust[1][numcell] = jd2;
423 88 : cellcount++;
424 88 : fInfcl[0][cellcount] = icl;
425 88 : fInfcl[1][cellcount] = jd1;
426 88 : fInfcl[2][cellcount] = jd2;
427 88 : }
428 : }
429 : }
430 : }
431 : }
432 : }
433 98 : return cellcount;
434 98 : }
435 : // ------------------------------------------------------------------------ //
436 : void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
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 : //
442 :
443 : AliPMDcludata *pmdcludata = 0;
444 :
445 : const Int_t kNmaxCell = 19; // # of cells surrounding a cluster center
446 :
447 196 : Int_t ndim = incr + 1;
448 :
449 : Int_t *ncl = 0x0;
450 : Int_t *clxy = 0x0;
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;
456 98 : Float_t clusdata[6];
457 : Double_t x1 = 0, y1 = 0, z1 = 0, x2 = 0, y2 = 0, z2 = 0, rr = 0;
458 :
459 98 : ncl = new Int_t [ndim];
460 98 : clxy = new Int_t [kNmaxCell];
461 :
462 : // Initialisation
463 900 : for(i = 0; i<ndim; i++)
464 : {
465 352 : ncl[i] = -1;
466 : }
467 1372 : for(i = 0; i<6; i++)
468 : {
469 588 : clusdata[i] = 0.;
470 : }
471 3920 : for(i = 0; i<19; i++)
472 : {
473 1862 : clxy[i] = 0;
474 : }
475 :
476 : // clno counts the final clusters
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
482 :
483 : clno = -1;
484 : nsupcl = -1;
485 :
486 900 : for(i = 0; i <= incr; i++)
487 : {
488 352 : if(fInfcl[0][i] != nsupcl)
489 : {
490 172 : nsupcl++;
491 172 : }
492 352 : if (nsupcl >= ndim)
493 : {
494 0 : AliWarning("RefClust: Too many superclusters!");
495 : nsupcl = ndim - 1;
496 0 : break;
497 : }
498 352 : ncl[nsupcl]++;
499 : }
500 :
501 294 : AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
502 : incr+1,nsupcl+1));
503 : id = -1;
504 : icl = -1;
505 :
506 540 : for(i = 0; i <= nsupcl; i++)
507 : {
508 172 : if(ncl[i] == 0)
509 : {
510 122 : id++;
511 122 : icl++;
512 122 : if (clno >= 4608)
513 : {
514 0 : AliWarning("RefClust: Too many clusters! more than 4608");
515 0 : delete [] ncl;
516 0 : delete [] clxy;
517 0 : return;
518 : }
519 122 : clno++;
520 122 : i1 = fInfcl[1][id];
521 122 : i2 = fInfcl[2][id];
522 :
523 122 : i12 = i1 + i2*kNDIMX;
524 :
525 122 : clusdata[0] = fCoord[0][i1][i2];
526 122 : clusdata[1] = fCoord[1][i1][i2];
527 122 : clusdata[2] = edepcell[i12];
528 122 : clusdata[3] = 1.;
529 122 : clusdata[4] = 999.5;
530 122 : clusdata[5] = 999.5;
531 :
532 122 : clxy[0] = i1*10000 + i2;
533 :
534 4636 : for(Int_t icltr = 1; icltr < kNmaxCell; icltr++)
535 : {
536 2196 : clxy[icltr] = -1;
537 : }
538 :
539 122 : pmdcludata = new AliPMDcludata(clusdata,clxy);
540 122 : fPMDclucont->Add(pmdcludata);
541 122 : }
542 50 : else if(ncl[i] == 1)
543 : {
544 24 : id++;
545 24 : icl++;
546 24 : if (clno >= 4608)
547 : {
548 0 : AliWarning("RefClust: Too many clusters! more than 4608");
549 0 : delete [] ncl;
550 0 : delete [] clxy;
551 :
552 0 : return;
553 : }
554 24 : clno++;
555 24 : i1 = fInfcl[1][id];
556 24 : i2 = fInfcl[2][id];
557 24 : i12 = i1 + i2*kNDIMX;
558 :
559 24 : x1 = fCoord[0][i1][i2];
560 24 : y1 = fCoord[1][i1][i2];
561 24 : z1 = edepcell[i12];
562 :
563 24 : clxy[0] = i1*10000 + i2;
564 :
565 24 : id++;
566 24 : i1 = fInfcl[1][id];
567 24 : i2 = fInfcl[2][id];
568 :
569 24 : i22 = i1 + i2*kNDIMX;
570 24 : x2 = fCoord[0][i1][i2];
571 24 : y2 = fCoord[1][i1][i2];
572 24 : z2 = edepcell[i22];
573 :
574 24 : clxy[1] = i1*10000 + i2;
575 :
576 :
577 864 : for(Int_t icltr = 2; icltr < kNmaxCell; icltr++)
578 : {
579 408 : clxy[icltr] = -1;
580 : }
581 :
582 24 : clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
583 24 : clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
584 24 : clusdata[2] = z1+z2;
585 24 : clusdata[3] = 2.;
586 24 : clusdata[4] = 0.5;
587 24 : clusdata[5] = 0.0;
588 24 : pmdcludata = new AliPMDcludata(clusdata,clxy);
589 24 : fPMDclucont->Add(pmdcludata);
590 24 : }
591 : else
592 : {
593 : Int_t *iord, *tc, *t;
594 : Double_t *x, *y, *z, *xc, *yc, *zc;
595 :
596 26 : iord = new Int_t [ncl[i]+1];
597 26 : tc = new Int_t [ncl[i]+1];
598 26 : t = new Int_t [ncl[i]+1];
599 26 : x = new Double_t [ncl[i]+1];
600 26 : y = new Double_t [ncl[i]+1];
601 26 : z = new Double_t [ncl[i]+1];
602 26 : xc = new Double_t [ncl[i]+1];
603 26 : yc = new Double_t [ncl[i]+1];
604 26 : zc = new Double_t [ncl[i]+1];
605 :
606 416 : for( k = 0; k < ncl[i]+1; k++)
607 : {
608 182 : iord[k] = -1;
609 182 : t[k] = -1;
610 182 : tc[k] = -1;
611 182 : x[k] = -1;
612 182 : y[k] = -1;
613 182 : z[k] = -1;
614 182 : xc[k] = -1;
615 182 : yc[k] = -1;
616 182 : zc[k] = -1;
617 : }
618 26 : id++;
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 26 : i1 = fInfcl[1][id];
624 26 : i2 = fInfcl[2][id];
625 26 : i12 = i1 + i2*kNDIMX;
626 :
627 26 : x[0] = fCoord[0][i1][i2];
628 26 : y[0] = fCoord[1][i1][i2];
629 26 : z[0] = edepcell[i12];
630 26 : t[0] = i1*10000 + i2;
631 :
632 :
633 26 : iord[0] = 0;
634 364 : for(j = 1; j <= ncl[i]; j++)
635 : {
636 156 : id++;
637 156 : i1 = fInfcl[1][id];
638 156 : i2 = fInfcl[2][id];
639 156 : i12 = i1 + i2*kNDIMX;
640 156 : iord[j] = j;
641 156 : x[j] = fCoord[0][i1][i2];
642 156 : y[j] = fCoord[1][i1][i2];
643 156 : z[j] = edepcell[i12];
644 156 : t[j] = i1*10000 + i2;
645 : }
646 :
647 : // arranging cells within supercluster in decreasing order
648 :
649 364 : for(j = 1;j <= ncl[i]; j++)
650 : {
651 : itest = 0;
652 156 : ihld = iord[j];
653 2008 : for(i1 = 0; i1 < j; i1++)
654 : {
655 1542 : if(itest == 0 && z[iord[i1]] < z[ihld])
656 : {
657 : itest = 1;
658 692 : for(i2 = j-1; i2 >= i1; i2--)
659 : {
660 250 : iord[i2+1] = iord[i2];
661 : }
662 96 : iord[i1] = ihld;
663 96 : }
664 : }
665 : }
666 :
667 26 : 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 416 : for(j = 0;j <= ncl[i]; j++)
676 : {
677 182 : supx += x[iord[j]]*z[iord[j]];
678 182 : supy += y[iord[j]]*z[iord[j]];
679 182 : supz += z[iord[j]];
680 182 : if(j < 19)
681 : {
682 182 : clxy[j] = t[iord[j]];
683 182 : }
684 : }
685 :
686 26 : if( ncl[i] + 1 < 19)
687 : {
688 672 : for(Int_t ncel = ncl[i] + 1; ncel < kNmaxCell; ncel ++ )
689 : {
690 312 : clxy[ncel] = -1;
691 : }
692 24 : }
693 26 : clusdata[0] = supx/supz;
694 26 : clusdata[1] = supy/supz;
695 26 : clusdata[2] = supz;
696 26 : clusdata[3] = ncl[i]+1;
697 26 : clusdata[4] = 0.5;
698 26 : clusdata[5] = 0.0;
699 26 : pmdcludata = new AliPMDcludata(clusdata,clxy);
700 26 : fPMDclucont->Add(pmdcludata);
701 26 : }
702 :
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 : */
711 :
712 26 : if (fClusParam == 2)
713 : {
714 : // This part is to split the supercluster
715 : //
716 : ig = 0;
717 0 : xc[ig] = x[iord[0]];
718 0 : yc[ig] = y[iord[0]];
719 0 : zc[ig] = z[iord[0]];
720 0 : tc[ig] = t[iord[0]];
721 : Int_t ivalid = 0, icount = 0;
722 :
723 0 : for(j=1;j<=ncl[i];j++)
724 : {
725 0 : x1 = x[iord[j]];
726 0 : y1 = y[iord[j]];
727 0 : z1 = z[iord[j]];
728 0 : t1 = t[iord[j]];
729 0 : rr=Distance(x1,y1,xc[ig],yc[ig]);
730 :
731 : // Check the cells which are outside the neighbours (rr>1.2)
732 0 : if(rr>1.2 )
733 : {
734 : ivalid=0;
735 : icount=0;
736 0 : for(Int_t j1=1;j1<j;j1++)
737 : {
738 0 : icount++;
739 0 : Float_t rr1=Distance(x1,y1,x[iord[j1]],y[iord[j1]]);
740 0 : if(rr1>1.2) ivalid++;
741 : }
742 0 : if(ivalid == icount && z1>0.5*zc[ig])
743 : {
744 0 : ig++;
745 0 : xc[ig]=x1;
746 0 : yc[ig]=y1;
747 0 : zc[ig]=z1;
748 0 : tc[ig]=t1;
749 0 : }
750 : }
751 : }
752 :
753 0 : icl=icl+ig+1;
754 :
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
763 :
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 0 : status = new Int_t [ncl[i]+1];
775 0 : cellXY = new Int_t *[ncl[i]+1];
776 :
777 0 : cellCount = new Int_t [ig+1];
778 0 : totaladc = new Double_t [ig+1];
779 0 : totaladc2 = new Double_t [ig+1];
780 0 : ncell = new Double_t [ig+1];
781 0 : weight = new Double_t [ig+1];
782 0 : xclust = new Double_t [ig+1];
783 0 : yclust = new Double_t [ig+1];
784 0 : sigxclust = new Double_t [ig+1];
785 0 : sigyclust = new Double_t [ig+1];
786 0 : ax = new Double_t [ig+1];
787 0 : ay = new Double_t [ig+1];
788 0 : ax2 = new Double_t [ig+1];
789 0 : ay2 = new Double_t [ig+1];
790 :
791 0 : for(j = 0; j < ncl[i]+1; j++)
792 : {
793 0 : status[j] = 0;
794 0 : cellXY[j] = new Int_t[ig+1];
795 : }
796 : //initialization
797 0 : for(Int_t kcl = 0; kcl < ig+1; kcl++)
798 : {
799 0 : cellCount[kcl] = 0;
800 0 : totaladc[kcl] = 0.;
801 0 : totaladc2[kcl] = 0.;
802 0 : ncell[kcl] = 0.;
803 0 : weight[kcl] = 0.;
804 0 : xclust[kcl] = 0.;
805 0 : yclust[kcl] = 0.;
806 0 : sigxclust[kcl] = 0.;
807 0 : sigyclust[kcl] = 0.;
808 0 : ax[kcl] = 0.;
809 0 : ay[kcl] = 0.;
810 0 : ax2[kcl] = 0.;
811 0 : ay2[kcl] = 0.;
812 0 : for(j = 0; j < ncl[i]+1; j++)
813 : {
814 0 : cellXY[j][kcl] = 0;
815 : }
816 : }
817 : Double_t sumweight, gweight;
818 :
819 0 : for(j = 0;j <= ncl[i]; j++)
820 : {
821 0 : x1 = x[iord[j]];
822 0 : y1 = y[iord[j]];
823 0 : z1 = z[iord[j]];
824 0 : t1 = t[iord[j]];
825 :
826 0 : for(Int_t kcl=0; kcl<=ig; kcl++)
827 : {
828 0 : x2 = xc[kcl];
829 0 : y2 = yc[kcl];
830 0 : rr = Distance(x1,y1,x2,y2);
831 0 : t2 = tc[kcl];
832 :
833 0 : if(rr==0)
834 : {
835 0 : ncell[kcl] = 1.;
836 0 : totaladc[kcl] = z1;
837 0 : totaladc2[kcl] = z1*z1;
838 0 : ax[kcl] = x1 * z1;
839 0 : ay[kcl] = y1 * z1;
840 0 : ax2[kcl] = 0.;
841 0 : ay2[kcl] = 0.;
842 0 : status[j] = 1;
843 0 : }
844 : }
845 : }
846 :
847 0 : for(j = 0; j <= ncl[i]; j++)
848 : {
849 : Int_t maxweight = 0;
850 : Double_t max = 0.;
851 :
852 0 : if(status[j] == 0)
853 : {
854 0 : x1 = x[iord[j]];
855 0 : y1 = y[iord[j]];
856 0 : z1 = z[iord[j]];
857 0 : t1 = t[iord[j]];
858 : sumweight = 0.;
859 :
860 0 : for(Int_t kcl = 0; kcl <= ig; kcl++)
861 : {
862 0 : x2 = xc[kcl];
863 0 : y2 = yc[kcl];
864 0 : rr = Distance(x1,y1,x2,y2);
865 0 : gweight = exp(-(rr*rr)/(2*(1.2*1.2)));
866 0 : weight[kcl] = zc[kcl] * gweight;
867 0 : sumweight = sumweight + weight[kcl];
868 :
869 0 : if(weight[kcl] > max)
870 : {
871 : max = weight[kcl];
872 : maxweight = kcl;
873 0 : }
874 : }
875 :
876 0 : cellXY[cellCount[maxweight]][maxweight] = iord[j];
877 :
878 0 : cellCount[maxweight]++;
879 :
880 0 : x2 = xc[maxweight];
881 0 : y2 = yc[maxweight];
882 0 : totaladc[maxweight] += z1;
883 0 : ax[maxweight] += x1*z1;
884 0 : ay[maxweight] += y1*z1;
885 0 : totaladc2[maxweight] += z1*z1;
886 0 : ax2[maxweight] += z1*(x1-x2)*(x1-x2);
887 0 : ay2[maxweight] += z1*(y1-y2)*(y1-y2);
888 0 : ncell[maxweight]++;
889 :
890 0 : }
891 : }
892 :
893 0 : for(Int_t kcl = 0; kcl <= ig; kcl++)
894 : {
895 0 : if(totaladc[kcl] > 0.)
896 : {
897 0 : xclust[kcl] = (ax[kcl])/ totaladc[kcl];
898 0 : yclust[kcl] = (ay[kcl])/ totaladc[kcl];
899 :
900 : //natasha
901 0 : Float_t sqtotadc = totaladc[kcl]*totaladc[kcl];
902 0 : if(totaladc2[kcl] >= sqtotadc)
903 : {
904 0 : sigxclust[kcl] = 0.25;
905 0 : sigyclust[kcl] = 0.25;
906 0 : }
907 : else
908 : {
909 0 : sigxclust[kcl] = (totaladc[kcl]/(sqtotadc-totaladc2[kcl]))*ax2[kcl];
910 0 : sigyclust[kcl] = (totaladc[kcl]/(sqtotadc-totaladc2[kcl]))*ay2[kcl];
911 : }
912 0 : }
913 :
914 0 : for(j = 0; j < cellCount[kcl]; j++) clno++;
915 :
916 0 : if (clno >= 4608)
917 : {
918 0 : AliWarning("RefClust: Too many clusters! more than 4608");
919 :
920 0 : delete [] cellCount;
921 0 : for(Int_t jj = 0; jj < ncl[i]+1; jj++) delete [] cellXY[jj];
922 0 : delete [] cellXY;
923 :
924 0 : delete [] status;
925 0 : delete [] totaladc;
926 0 : delete [] totaladc2;
927 0 : delete [] ncell;
928 0 : delete [] xclust;
929 0 : delete [] yclust;
930 0 : delete [] sigxclust;
931 0 : delete [] sigyclust;
932 0 : delete [] ax;
933 0 : delete [] ay;
934 0 : delete [] ax2;
935 0 : delete [] ay2;
936 0 : delete [] weight;
937 :
938 0 : delete [] iord;
939 0 : delete [] tc;
940 0 : delete [] t;
941 0 : delete [] x;
942 0 : delete [] y;
943 0 : delete [] z;
944 0 : delete [] xc;
945 0 : delete [] yc;
946 0 : delete [] zc;
947 :
948 :
949 0 : delete [] ncl;
950 0 : delete [] clxy;
951 :
952 0 : return;
953 : }
954 0 : clusdata[0] = xclust[kcl];
955 0 : clusdata[1] = yclust[kcl];
956 0 : clusdata[2] = totaladc[kcl];
957 0 : clusdata[3] = ncell[kcl];
958 :
959 0 : if(sigxclust[kcl] > sigyclust[kcl])
960 : {
961 0 : clusdata[4] = TMath::Sqrt(sigxclust[kcl]);
962 0 : clusdata[5] = TMath::Sqrt(sigyclust[kcl]);
963 0 : }
964 : else
965 : {
966 0 : clusdata[4] = TMath::Sqrt(sigyclust[kcl]);
967 0 : clusdata[5] = TMath::Sqrt(sigxclust[kcl]);
968 : }
969 :
970 0 : clxy[0] = tc[kcl];
971 :
972 : Int_t Ncell=1;
973 0 : for (Int_t ii = 0; ii < cellCount[kcl]; ii++)
974 : {
975 0 : if(ii<18)
976 : {
977 0 : clxy[Ncell] = t[cellXY[ii][kcl]];
978 0 : Ncell++;
979 0 : }
980 : }
981 :
982 0 : pmdcludata = new AliPMDcludata(clusdata,clxy);
983 0 : fPMDclucont->Add(pmdcludata);
984 : }
985 0 : delete [] cellCount;
986 0 : for(Int_t jj = 0; jj < ncl[i]+1; jj++) delete [] cellXY[jj];
987 0 : delete [] cellXY;
988 :
989 0 : delete [] status;
990 0 : delete [] totaladc;
991 0 : delete [] totaladc2;
992 0 : delete [] ncell;
993 0 : delete [] xclust;
994 0 : delete [] yclust;
995 0 : delete [] sigxclust;
996 0 : delete [] sigyclust;
997 0 : delete [] ax;
998 0 : delete [] ay;
999 0 : delete [] ax2;
1000 0 : delete [] ay2;
1001 0 : delete [] weight;
1002 :
1003 0 : }
1004 :
1005 52 : delete [] iord;
1006 52 : delete [] tc;
1007 52 : delete [] t;
1008 52 : delete [] x;
1009 52 : delete [] y;
1010 52 : delete [] z;
1011 52 : delete [] xc;
1012 52 : delete [] yc;
1013 52 : delete [] zc;
1014 :
1015 :
1016 26 : }
1017 : }
1018 196 : delete [] ncl;
1019 196 : delete [] clxy;
1020 196 : }
1021 : // ------------------------------------------------------------------------ //
1022 : Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1,
1023 : Double_t x2, Double_t y2)
1024 : {
1025 0 : return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
1026 : }
1027 : // ------------------------------------------------------------------------ //
1028 : void AliPMDClusteringV1::SetEdepCut(Float_t decut)
1029 : {
1030 196 : fCutoff = decut;
1031 98 : }
1032 : // ------------------------------------------------------------------------ //
1033 : void AliPMDClusteringV1::SetClusteringParam(Int_t cluspar)
1034 : {
1035 196 : fClusParam = cluspar;
1036 98 : }
1037 : // ------------------------------------------------------------------------ //
|
