]>
Commit | Line | Data |
---|---|---|
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 | ||
57 | ClassImp(AliPMDClusteringV1) | |
58 | ||
59 | const Double_t AliPMDClusteringV1::fgkSqroot3by2=0.8660254; // sqrt(3.)/2. | |
60 | ||
61 | AliPMDClusteringV1::AliPMDClusteringV1(): | |
562718f9 | 62 | fPMDclucont(new TObjArray()), |
3edbbba2 | 63 | fCutoff(0.0) |
64 | { | |
fd30f88e | 65 | for(Int_t i = 0; i < kNDIMX; i++) |
3edbbba2 | 66 | { |
fd30f88e | 67 | for(Int_t j = 0; j < kNDIMY; j++) |
3edbbba2 | 68 | { |
69 | fCoord[0][i][j] = i+j/2.; | |
70 | fCoord[1][i][j] = fgkSqroot3by2*j; | |
3edbbba2 | 71 | } |
72 | } | |
73 | } | |
74 | // ------------------------------------------------------------------------ // | |
562718f9 | 75 | AliPMDClusteringV1::AliPMDClusteringV1(const AliPMDClusteringV1& pmdclv1): |
76 | AliPMDClustering(pmdclv1), | |
77 | fPMDclucont(0), | |
78 | fCutoff(0) | |
79 | { | |
80 | // copy constructor | |
81 | AliError("Copy constructor not allowed "); | |
82 | ||
83 | } | |
84 | // ------------------------------------------------------------------------ // | |
85 | AliPMDClusteringV1 &AliPMDClusteringV1::operator=(const AliPMDClusteringV1& /*pmdclv1*/) | |
86 | { | |
87 | // copy constructor | |
88 | AliError("Assignment operator not allowed "); | |
89 | return *this; | |
90 | } | |
91 | // ------------------------------------------------------------------------ // | |
3edbbba2 | 92 | AliPMDClusteringV1::~AliPMDClusteringV1() |
93 | { | |
562718f9 | 94 | delete fPMDclucont; |
3edbbba2 | 95 | } |
96 | // ------------------------------------------------------------------------ // | |
fd30f88e | 97 | void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn, |
98 | Double_t celladc[48][96], TObjArray *pmdcont) | |
3edbbba2 | 99 | { |
100 | // main function to call other necessary functions to do clustering | |
101 | // | |
01c4d84a | 102 | |
562718f9 | 103 | AliPMDcluster *pmdcl = 0; |
3edbbba2 | 104 | |
562718f9 | 105 | Int_t i, j, nmx1, incr, id, jd; |
106 | Int_t celldataX[15], celldataY[15]; | |
107 | Float_t clusdata[6]; | |
108 | Double_t cutoff, ave; | |
109 | Double_t edepcell[kNMX]; | |
3edbbba2 | 110 | |
111 | const float ktwobysqrt3 = 1.1547; // 2./sqrt(3.) | |
112 | ||
01c4d84a | 113 | // ndimXr and ndimYr are different because of different module size |
114 | ||
fd30f88e | 115 | Int_t ndimXr = 0; |
116 | Int_t ndimYr = 0; | |
01c4d84a | 117 | |
118 | if (ismn < 12) | |
119 | { | |
120 | ndimXr = 96; | |
121 | ndimYr = 48; | |
122 | } | |
123 | else if (ismn >= 12 && ismn <= 23) | |
124 | { | |
125 | ndimXr = 48; | |
126 | ndimYr = 96; | |
127 | } | |
128 | ||
562718f9 | 129 | Int_t kk = 0; |
fd30f88e | 130 | for (Int_t i = 0; i < kNDIMX; i++) |
3edbbba2 | 131 | { |
fd30f88e | 132 | for (Int_t j = 0; j < kNDIMY; j++) |
01c4d84a | 133 | { |
01c4d84a | 134 | fCellTrNo[i][j] = -1; |
562718f9 | 135 | edepcell[kk] = 0.; |
136 | kk++; | |
01c4d84a | 137 | } |
138 | } | |
139 | ||
140 | for (id = 0; id < ndimXr; id++) | |
141 | { | |
142 | for (jd = 0; jd < ndimYr; jd++) | |
3edbbba2 | 143 | { |
fd30f88e | 144 | j = jd; |
145 | i = id+(ndimYr/2-1)-(jd/2); | |
01c4d84a | 146 | |
75de4f00 | 147 | //PH Int_t ij = i + j*kNDIMX; |
148 | Int_t ij = i + j*ndimXr; | |
562718f9 | 149 | |
01c4d84a | 150 | if (ismn < 12) |
151 | { | |
562718f9 | 152 | edepcell[ij] = celladc[jd][id]; |
153 | fCellTrNo[i][j] = jd*10000+id; // for association | |
01c4d84a | 154 | } |
155 | else if (ismn >= 12 && ismn <= 23) | |
156 | { | |
562718f9 | 157 | edepcell[ij] = celladc[id][jd]; |
158 | fCellTrNo[i][j] = id*10000+jd; // for association | |
01c4d84a | 159 | } |
160 | ||
3edbbba2 | 161 | } |
162 | } | |
562718f9 | 163 | |
164 | Int_t iord1[kNMX]; | |
165 | TMath::Sort(kNMX,edepcell,iord1);// order the data | |
166 | cutoff = fCutoff; // cutoff to discard cells | |
167 | ave = 0.; | |
fd30f88e | 168 | nmx1 = -1; |
562718f9 | 169 | for(i = 0;i < kNMX; i++) |
3edbbba2 | 170 | { |
562718f9 | 171 | if(edepcell[i] > 0.) |
fd30f88e | 172 | { |
562718f9 | 173 | ave += edepcell[i]; |
fd30f88e | 174 | } |
562718f9 | 175 | if(edepcell[i] > cutoff ) |
fd30f88e | 176 | { |
177 | nmx1++; | |
178 | } | |
3edbbba2 | 179 | } |
3edbbba2 | 180 | |
181 | AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1)); | |
182 | ||
3edbbba2 | 183 | if (nmx1 == 0) nmx1 = 1; |
fd30f88e | 184 | ave = ave/nmx1; |
3edbbba2 | 185 | |
186 | AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f", | |
187 | kNMX,ave)); | |
562718f9 | 188 | |
189 | incr = CrClust(ave, cutoff, nmx1,iord1, edepcell ); | |
562718f9 | 190 | RefClust(incr,edepcell); |
191 | Int_t nentries1 = fPMDclucont->GetEntries(); | |
fd30f88e | 192 | AliDebug(1,Form("Detector Plane = %d Serial Module No = %d Number of clusters = %d",idet, ismn, nentries1)); |
193 | AliDebug(1,Form("Total number of clusters/module = %d",nentries1)); | |
562718f9 | 194 | |
fd30f88e | 195 | for (Int_t ient1 = 0; ient1 < nentries1; ient1++) |
3edbbba2 | 196 | { |
fd30f88e | 197 | AliPMDcludata *pmdcludata = |
562718f9 | 198 | (AliPMDcludata*)fPMDclucont->UncheckedAt(ient1); |
fd30f88e | 199 | Float_t cluXC = pmdcludata->GetClusX(); |
200 | Float_t cluYC = pmdcludata->GetClusY(); | |
201 | Float_t cluADC = pmdcludata->GetClusADC(); | |
202 | Float_t cluCELLS = pmdcludata->GetClusCells(); | |
203 | Float_t cluSIGX = pmdcludata->GetClusSigmaX(); | |
204 | Float_t cluSIGY = pmdcludata->GetClusSigmaY(); | |
205 | ||
3edbbba2 | 206 | Float_t cluY0 = ktwobysqrt3*cluYC; |
207 | Float_t cluX0 = cluXC - cluY0/2.; | |
fd30f88e | 208 | |
3edbbba2 | 209 | // |
210 | // Cluster X centroid is back transformed | |
211 | // | |
01c4d84a | 212 | if (ismn < 12) |
213 | { | |
214 | clusdata[0] = cluX0 - (24-1) + cluY0/2.; | |
215 | } | |
fd30f88e | 216 | else if ( ismn >= 12 && ismn <= 23) |
01c4d84a | 217 | { |
218 | clusdata[0] = cluX0 - (48-1) + cluY0/2.; | |
219 | } | |
fd30f88e | 220 | |
01c4d84a | 221 | clusdata[1] = cluY0; |
222 | clusdata[2] = cluADC; | |
223 | clusdata[3] = cluCELLS; | |
fd30f88e | 224 | clusdata[4] = cluSIGX; |
225 | clusdata[5] = cluSIGY; | |
226 | ||
3edbbba2 | 227 | // |
228 | // Cells associated with a cluster | |
229 | // | |
562718f9 | 230 | |
3edbbba2 | 231 | for (Int_t ihit = 0; ihit < 15; ihit++) |
232 | { | |
01c4d84a | 233 | if (ismn < 12) |
234 | { | |
562718f9 | 235 | celldataX[ihit] = pmdcludata->GetCellXY(ihit)%10000; |
236 | celldataY[ihit] = pmdcludata->GetCellXY(ihit)/10000; | |
01c4d84a | 237 | } |
238 | else if (ismn >= 12 && ismn <= 23) | |
239 | { | |
562718f9 | 240 | celldataX[ihit] = pmdcludata->GetCellXY(ihit)/10000; |
241 | celldataY[ihit] = pmdcludata->GetCellXY(ihit)%10000; | |
01c4d84a | 242 | } |
3edbbba2 | 243 | } |
3edbbba2 | 244 | pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY); |
245 | pmdcont->Add(pmdcl); | |
246 | } | |
fd30f88e | 247 | |
562718f9 | 248 | fPMDclucont->Clear(); |
3edbbba2 | 249 | |
3edbbba2 | 250 | } |
251 | // ------------------------------------------------------------------------ // | |
562718f9 | 252 | Int_t AliPMDClusteringV1::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1, |
253 | Int_t iord1[], Double_t edepcell[]) | |
3edbbba2 | 254 | { |
fd30f88e | 255 | // Does crude clustering |
3edbbba2 | 256 | // Finds out only the big patch by just searching the |
257 | // connected cells | |
258 | // | |
fd30f88e | 259 | Int_t i,j,k,id1,id2,icl, numcell, clust[2][5000]; |
260 | Int_t jd1,jd2, icell, cellcount; | |
261 | static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1}; | |
262 | ||
3edbbba2 | 263 | AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff)); |
264 | ||
fd30f88e | 265 | for (j = 0; j < kNDIMX; j++) |
266 | { | |
267 | for(k = 0; k < kNDIMY; k++) | |
268 | { | |
269 | fInfocl[0][j][k] = 0; | |
270 | fInfocl[1][j][k] = 0; | |
271 | } | |
272 | } | |
273 | for(i=0; i < kNMX; i++) | |
274 | { | |
275 | fInfcl[0][i] = -1; | |
562718f9 | 276 | |
277 | j = iord1[i]; | |
278 | id2 = j/kNDIMX; | |
279 | id1 = j-id2*kNDIMX; | |
280 | ||
281 | if(edepcell[j] <= cutoff) | |
fd30f88e | 282 | { |
283 | fInfocl[0][id1][id2] = -1; | |
284 | } | |
3edbbba2 | 285 | } |
562718f9 | 286 | |
3edbbba2 | 287 | // --------------------------------------------------------------- |
288 | // crude clustering begins. Start with cell having largest adc | |
289 | // count and loop over the cells in descending order of adc count | |
290 | // --------------------------------------------------------------- | |
562718f9 | 291 | |
fd30f88e | 292 | icl = -1; |
293 | cellcount = -1; | |
3edbbba2 | 294 | |
fd30f88e | 295 | for(icell = 0; icell <= nmx1; icell++) |
296 | { | |
562718f9 | 297 | j = iord1[icell]; |
298 | id2 = j/kNDIMX; | |
299 | id1 = j-id2*kNDIMX; | |
300 | ||
fd30f88e | 301 | if(fInfocl[0][id1][id2] == 0 ) |
302 | { | |
303 | icl++; | |
304 | numcell = 0; | |
305 | cellcount++; | |
306 | fInfocl[0][id1][id2] = 1; | |
307 | fInfocl[1][id1][id2] = icl; | |
308 | fInfcl[0][cellcount] = icl; | |
309 | fInfcl[1][cellcount] = id1; | |
310 | fInfcl[2][cellcount] = id2; | |
311 | ||
312 | clust[0][numcell] = id1; | |
313 | clust[1][numcell] = id2; | |
314 | ||
315 | for(i = 1; i < 5000; i++) | |
316 | { | |
317 | clust[0][i]=0; | |
318 | } | |
319 | // --------------------------------------------------------------- | |
320 | // check for adc count in neib. cells. If ne 0 put it in this clust | |
321 | // --------------------------------------------------------------- | |
322 | for(i = 0; i < 6; i++) | |
323 | { | |
324 | jd1 = id1 + neibx[i]; | |
325 | jd2 = id2 + neiby[i]; | |
326 | if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) && | |
327 | fInfocl[0][jd1][jd2] == 0) | |
328 | { | |
329 | numcell++; | |
330 | fInfocl[0][jd1][jd2] = 2; | |
331 | fInfocl[1][jd1][jd2] = icl; | |
332 | clust[0][numcell] = jd1; | |
333 | clust[1][numcell] = jd2; | |
334 | cellcount++; | |
335 | fInfcl[0][cellcount] = icl; | |
336 | fInfcl[1][cellcount] = jd1; | |
337 | fInfcl[2][cellcount] = jd2; | |
338 | } | |
339 | } | |
340 | // --------------------------------------------------------------- | |
341 | // check adc count for neighbour's neighbours recursively and | |
342 | // if nonzero, add these to the cluster. | |
343 | // --------------------------------------------------------------- | |
344 | for(i = 1; i < 5000;i++) | |
345 | { | |
346 | if(clust[0][i] != 0) | |
347 | { | |
348 | id1 = clust[0][i]; | |
349 | id2 = clust[1][i]; | |
350 | for(j = 0; j < 6 ; j++) | |
351 | { | |
352 | jd1 = id1 + neibx[j]; | |
353 | jd2 = id2 + neiby[j]; | |
354 | if( (jd1 >= 0 && jd1 < kNDIMX) && | |
355 | (jd2 >= 0 && jd2 < kNDIMY) && | |
356 | fInfocl[0][jd1][jd2] == 0 ) | |
357 | { | |
358 | fInfocl[0][jd1][jd2] = 2; | |
359 | fInfocl[1][jd1][jd2] = icl; | |
360 | numcell++; | |
361 | clust[0][numcell] = jd1; | |
362 | clust[1][numcell] = jd2; | |
363 | cellcount++; | |
364 | fInfcl[0][cellcount] = icl; | |
365 | fInfcl[1][cellcount] = jd1; | |
366 | fInfcl[2][cellcount] = jd2; | |
367 | } | |
368 | } | |
369 | } | |
3edbbba2 | 370 | } |
3edbbba2 | 371 | } |
3edbbba2 | 372 | } |
3edbbba2 | 373 | return cellcount; |
374 | } | |
375 | // ------------------------------------------------------------------------ // | |
562718f9 | 376 | void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[]) |
3edbbba2 | 377 | { |
378 | // Does the refining of clusters | |
379 | // Takes the big patch and does gaussian fitting and | |
380 | // finds out the more refined clusters | |
381 | // | |
22d01768 | 382 | |
562718f9 | 383 | |
562718f9 | 384 | |
22d01768 | 385 | AliPMDcludata *pmdcludata = 0; |
562718f9 | 386 | |
973b7456 | 387 | Int_t *cellCount = 0x0; |
388 | Int_t **cellXY = 0x0; | |
22d01768 | 389 | const Int_t kdim = 4500; |
562718f9 | 390 | |
22d01768 | 391 | Int_t i, j, k, i1, i2, id, icl, itest,ihld, ig, nsupcl,clno; |
392 | Int_t t[kdim]; | |
973b7456 | 393 | Int_t ncl[kdim], iord[kdim], lev1[kdim], lev2[kdim]; |
22d01768 | 394 | Int_t clxy[15]; |
562718f9 | 395 | Float_t clusdata[6]; |
396 | Double_t x1, y1, z1, x2, y2, z2, dist,rr,sum; | |
fd30f88e | 397 | Double_t x[kdim], y[kdim], z[kdim]; |
398 | Double_t xc[kdim], yc[kdim], zc[kdim], cells[kdim], rc[kdim]; | |
562718f9 | 399 | |
400 | // Initialisation | |
fd30f88e | 401 | for(i = 0; i<kdim; i++) |
402 | { | |
403 | t[i] = -1; | |
562718f9 | 404 | ncl[i] = -1; |
405 | if (i < 6) clusdata[i] = 0.; | |
406 | if (i < 15) clxy[i] = 0; | |
fd30f88e | 407 | } |
562718f9 | 408 | |
fd30f88e | 409 | // clno counts the final clusters |
3edbbba2 | 410 | // nsupcl = # of superclusters; ncl[i]= # of cells in supercluster i |
411 | // x, y and z store (x,y) coordinates of and energy deposited in a cell | |
412 | // xc, yc store (x,y) coordinates of the cluster center | |
413 | // zc stores the energy deposited in a cluster | |
414 | // rc is cluster radius | |
fd30f88e | 415 | |
416 | clno = -1; | |
3edbbba2 | 417 | nsupcl = -1; |
562718f9 | 418 | |
fd30f88e | 419 | for(i = 0; i <= incr; i++) |
420 | { | |
421 | if(fInfcl[0][i] != nsupcl) | |
01c4d84a | 422 | { |
fd30f88e | 423 | nsupcl++; |
01c4d84a | 424 | } |
fd30f88e | 425 | if (nsupcl > kdim) |
01c4d84a | 426 | { |
fd30f88e | 427 | AliWarning("RefClust: Too many superclusters!"); |
428 | nsupcl = kdim; | |
429 | break; | |
01c4d84a | 430 | } |
fd30f88e | 431 | ncl[nsupcl]++; |
3edbbba2 | 432 | } |
fd30f88e | 433 | |
434 | AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d", | |
435 | incr+1,nsupcl+1)); | |
436 | id = -1; | |
437 | icl = -1; | |
562718f9 | 438 | |
fd30f88e | 439 | for(i = 0; i <= nsupcl; i++) |
440 | { | |
441 | if(ncl[i] == 0) | |
01c4d84a | 442 | { |
fd30f88e | 443 | id++; |
444 | icl++; | |
445 | if (clno >= 5000) | |
446 | { | |
447 | AliWarning("RefClust: Too many clusters! more than 5000"); | |
448 | return; | |
449 | } | |
450 | clno++; | |
451 | i1 = fInfcl[1][id]; | |
452 | i2 = fInfcl[2][id]; | |
453 | ||
562718f9 | 454 | Int_t i12 = i1 + i2*kNDIMX; |
455 | ||
fd30f88e | 456 | clusdata[0] = fCoord[0][i1][i2]; |
457 | clusdata[1] = fCoord[1][i1][i2]; | |
562718f9 | 458 | clusdata[2] = edepcell[i12]; |
fd30f88e | 459 | clusdata[3] = 1.; |
460 | clusdata[4] = 0.5; | |
461 | clusdata[5] = 0.0; | |
562718f9 | 462 | |
463 | clxy[0] = fCellTrNo[i1][i2]; //association | |
464 | for(Int_t icltr = 1; icltr < 15; icltr++) | |
fd30f88e | 465 | { |
562718f9 | 466 | clxy[icltr] = -1; |
3edbbba2 | 467 | } |
562718f9 | 468 | pmdcludata = new AliPMDcludata(clusdata,clxy); |
469 | fPMDclucont->Add(pmdcludata); | |
3edbbba2 | 470 | } |
fd30f88e | 471 | else if(ncl[i] == 1) |
472 | { | |
473 | id++; | |
474 | icl++; | |
475 | if (clno >= 5000) | |
476 | { | |
477 | AliWarning("RefClust: Too many clusters! more than 5000"); | |
478 | return; | |
479 | } | |
480 | clno++; | |
481 | i1 = fInfcl[1][id]; | |
482 | i2 = fInfcl[2][id]; | |
562718f9 | 483 | Int_t i12 = i1 + i2*kNDIMX; |
484 | ||
fd30f88e | 485 | x1 = fCoord[0][i1][i2]; |
486 | y1 = fCoord[1][i1][i2]; | |
562718f9 | 487 | z1 = edepcell[i12]; |
488 | clxy[0] = fCellTrNo[i1][i2]; //asso | |
489 | id++; | |
fd30f88e | 490 | i1 = fInfcl[1][id]; |
491 | i2 = fInfcl[2][id]; | |
562718f9 | 492 | |
493 | Int_t i22 = i1 + i2*kNDIMX; | |
fd30f88e | 494 | x2 = fCoord[0][i1][i2]; |
495 | y2 = fCoord[1][i1][i2]; | |
562718f9 | 496 | z2 = edepcell[i22]; |
497 | clxy[1] = fCellTrNo[i1][i2]; //asso | |
498 | for(Int_t icltr = 2; icltr < 15; icltr++) | |
fd30f88e | 499 | { |
562718f9 | 500 | clxy[icltr] = -1; |
fd30f88e | 501 | } |
fd30f88e | 502 | |
503 | clusdata[0] = (x1*z1+x2*z2)/(z1+z2); | |
504 | clusdata[1] = (y1*z1+y2*z2)/(z1+z2); | |
505 | clusdata[2] = z1+z2; | |
506 | clusdata[3] = 2.; | |
507 | clusdata[4] = 0.5; | |
508 | clusdata[5] = 0.0; | |
562718f9 | 509 | pmdcludata = new AliPMDcludata(clusdata,clxy); |
510 | fPMDclucont->Add(pmdcludata); | |
3edbbba2 | 511 | } |
fd30f88e | 512 | else |
513 | { | |
fd30f88e | 514 | id++; |
515 | iord[0] = 0; | |
516 | // super-cluster of more than two cells - broken up into smaller | |
517 | // clusters gaussian centers computed. (peaks separated by > 1 cell) | |
518 | // Begin from cell having largest energy deposited This is first | |
519 | // cluster center | |
520 | i1 = fInfcl[1][id]; | |
521 | i2 = fInfcl[2][id]; | |
562718f9 | 522 | Int_t i12 = i1 + i2*kNDIMX; |
523 | ||
fd30f88e | 524 | x[0] = fCoord[0][i1][i2]; |
525 | y[0] = fCoord[1][i1][i2]; | |
562718f9 | 526 | z[0] = edepcell[i12]; |
527 | ||
528 | t[0] = fCellTrNo[i1][i2]; //asso | |
fd30f88e | 529 | |
530 | iord[0] = 0; | |
531 | for(j = 1; j <= ncl[i]; j++) | |
532 | { | |
533 | id++; | |
534 | i1 = fInfcl[1][id]; | |
535 | i2 = fInfcl[2][id]; | |
562718f9 | 536 | Int_t i12 = i1 + i2*kNDIMX; |
537 | ||
fd30f88e | 538 | iord[j] = j; |
539 | x[j] = fCoord[0][i1][i2]; | |
540 | y[j] = fCoord[1][i1][i2]; | |
562718f9 | 541 | z[j] = edepcell[i12]; |
542 | ||
543 | t[j] = fCellTrNo[i1][i2]; //asso | |
3edbbba2 | 544 | } |
fd30f88e | 545 | |
fd30f88e | 546 | // arranging cells within supercluster in decreasing order |
547 | ||
548 | for(j = 1;j <= ncl[i]; j++) | |
549 | { | |
550 | itest = 0; | |
551 | ihld = iord[j]; | |
552 | for(i1 = 0; i1 < j; i1++) | |
553 | { | |
554 | if(itest == 0 && z[iord[i1]] < z[ihld]) | |
555 | { | |
556 | itest = 1; | |
557 | for(i2 = j-1; i2 >= i1; i2--) | |
558 | { | |
559 | iord[i2+1] = iord[i2]; | |
560 | } | |
561 | iord[i1] = ihld; | |
562 | } | |
563 | } | |
3edbbba2 | 564 | } |
fd30f88e | 565 | // compute the number of Gaussians and their centers ( first |
566 | // guess ) | |
567 | // centers must be separated by cells having smaller ener. dep. | |
568 | // neighbouring centers should be either strong or well-separated | |
562718f9 | 569 | ig = 0; |
fd30f88e | 570 | xc[ig] = x[iord[0]]; |
571 | yc[ig] = y[iord[0]]; | |
572 | zc[ig] = z[iord[0]]; | |
573 | for(j = 1; j <= ncl[i]; j++) | |
574 | { | |
575 | itest = -1; | |
576 | x1 = x[iord[j]]; | |
577 | y1 = y[iord[j]]; | |
578 | for(k = 0; k <= ig; k++) | |
579 | { | |
580 | x2 = xc[k]; | |
581 | y2 = yc[k]; | |
582 | rr = Distance(x1,y1,x2,y2); | |
562718f9 | 583 | if(rr >= 1.1 && rr < 1.8 && z[iord[j]] > zc[k]/4.)itest++; |
584 | if(rr >= 1.8 && rr < 2.1 && z[iord[j]] > zc[k]/10.)itest++; | |
585 | if( rr >= 2.1)itest++; | |
fd30f88e | 586 | } |
587 | if(itest == ig) | |
588 | { | |
589 | ig++; | |
590 | xc[ig] = x1; | |
591 | yc[ig] = y1; | |
592 | zc[ig] = z[iord[j]]; | |
593 | } | |
3edbbba2 | 594 | } |
fd30f88e | 595 | GaussFit(ncl[i], ig, x[0], y[0] ,z[0], xc[0], yc[0], zc[0], rc[0]); |
596 | icl += ig+1; | |
597 | // compute the number of cells belonging to each cluster. | |
598 | // cell is shared between several clusters ( if they are equidistant | |
599 | // from it ) in the ratio of cluster energy deposition | |
562718f9 | 600 | |
601 | Int_t jj = 15; | |
602 | cellCount = new Int_t [ig+1]; | |
603 | cellXY = new Int_t *[jj]; | |
604 | for(Int_t ij = 0; ij < 15; ij++) cellXY[ij] = new Int_t [ig+1]; | |
605 | ||
fd30f88e | 606 | for(j = 0; j <= ig; j++) |
01c4d84a | 607 | { |
562718f9 | 608 | cellCount[j] = 0; |
609 | cells[j] = 0.; | |
fd30f88e | 610 | } |
562718f9 | 611 | |
fd30f88e | 612 | if(ig > 0) |
613 | { | |
614 | for(j = 0; j <= ncl[i]; j++) | |
615 | { | |
616 | lev1[j] = 0; | |
617 | lev2[j] = 0; | |
618 | for(k = 0; k <= ig; k++) | |
619 | { | |
620 | dist = Distance(x[j], y[j], xc[k], yc[k]); | |
621 | if(dist < TMath::Sqrt(3.) ) | |
622 | { | |
623 | //asso | |
562718f9 | 624 | if (cellCount[k] < 15) |
625 | { | |
626 | cellXY[cellCount[k]][k] = t[j]; | |
627 | } | |
fd30f88e | 628 | cellCount[k]++; |
629 | // | |
630 | lev1[0]++; | |
631 | i1 = lev1[0]; | |
632 | lev1[i1] = k; | |
633 | } | |
634 | else | |
635 | { | |
636 | if(dist < 2.1) | |
637 | { | |
638 | lev2[0]++; | |
639 | i1 = lev2[0]; | |
640 | lev2[i1] = k; | |
641 | } | |
642 | } | |
643 | } | |
644 | if(lev1[0] != 0) | |
645 | { | |
646 | if(lev1[0] == 1) | |
647 | { | |
648 | cells[lev1[1]]++; | |
649 | } | |
650 | else | |
651 | { | |
652 | sum=0.; | |
653 | for(k = 1; k <= lev1[0]; k++) | |
654 | { | |
655 | sum += zc[lev1[k]]; | |
656 | } | |
657 | for(k = 1; k <= lev1[0]; k++) | |
658 | { | |
659 | cells[lev1[k]] += zc[lev1[k]]/sum; | |
660 | } | |
661 | } | |
662 | } | |
663 | else | |
664 | { | |
665 | if(lev2[0] == 0) | |
666 | { | |
667 | cells[lev2[1]]++; | |
668 | } | |
669 | else | |
670 | { | |
671 | sum=0.; | |
672 | for( k = 1; k <= lev2[0]; k++) | |
673 | { | |
674 | sum += zc[lev2[k]]; | |
675 | } | |
676 | for(k = 1; k <= lev2[0]; k++) | |
677 | { | |
678 | cells[lev2[k]] += zc[lev2[k]]/sum; | |
679 | } | |
680 | } | |
681 | } | |
682 | } | |
683 | } | |
684 | ||
685 | // zero rest of the cell array | |
686 | //asso | |
687 | for( k = 0; k <= ig; k++) | |
688 | { | |
562718f9 | 689 | for(Int_t icltr = cellCount[k]; icltr < 15; icltr++) |
fd30f88e | 690 | { |
562718f9 | 691 | cellXY[icltr][k] = -1; |
fd30f88e | 692 | } |
693 | } | |
694 | // | |
695 | ||
696 | for(j = 0; j <= ig; j++) | |
697 | { | |
698 | clno++; | |
699 | if (clno >= 5000) | |
700 | { | |
701 | AliWarning("RefClust: Too many clusters! more than 5000"); | |
702 | return; | |
703 | } | |
704 | clusdata[0] = xc[j]; | |
705 | clusdata[1] = yc[j]; | |
706 | clusdata[2] = zc[j]; | |
707 | clusdata[4] = rc[j]; | |
708 | clusdata[5] = 0.0; | |
709 | if(ig == 0) | |
710 | { | |
711 | clusdata[3] = ncl[i]; | |
712 | } | |
713 | else | |
714 | { | |
715 | clusdata[3] = cells[j]; | |
716 | } | |
562718f9 | 717 | |
718 | ||
719 | for (Int_t ii=0; ii < 15; ii++) | |
720 | { | |
721 | clxy[ii] = cellXY[ii][j]; | |
722 | } | |
723 | pmdcludata = new AliPMDcludata(clusdata,clxy); | |
724 | fPMDclucont->Add(pmdcludata); | |
01c4d84a | 725 | } |
562718f9 | 726 | delete [] cellCount; |
727 | for(Int_t jj = 0; jj < 15; jj++)delete [] cellXY[jj]; | |
728 | delete [] cellXY; | |
01c4d84a | 729 | } |
3edbbba2 | 730 | } |
3edbbba2 | 731 | } |
732 | // ------------------------------------------------------------------------ // | |
fd30f88e | 733 | void AliPMDClusteringV1::GaussFit(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, Double_t &rc) | |
3edbbba2 | 736 | { |
737 | // Does gaussian fitting | |
738 | // | |
562718f9 | 739 | |
fd30f88e | 740 | const Int_t kdim = 4500; |
741 | Int_t i, j, i1, i2, novar, idd, jj; | |
562718f9 | 742 | Int_t neib[kdim][50]; |
743 | ||
fd30f88e | 744 | Double_t sum, dx, dy, str, str1, aint, sum1, rr, dum; |
745 | Double_t x1, x2, y1, y2; | |
746 | Double_t xx[kdim], yy[kdim], zz[kdim], xxc[kdim], yyc[kdim]; | |
747 | Double_t a[kdim], b[kdim], c[kdim], d[kdim], ha[kdim], hb[kdim]; | |
748 | Double_t hc[kdim], hd[kdim], zzc[kdim], rrc[kdim]; | |
fd30f88e | 749 | |
750 | TRandom rnd; | |
751 | ||
3edbbba2 | 752 | str = 0.; |
753 | str1 = 0.; | |
754 | rr = 0.3; | |
755 | novar = 0; | |
562718f9 | 756 | j = 0; |
3edbbba2 | 757 | |
fd30f88e | 758 | for(i = 0; i <= ncell; i++) |
3edbbba2 | 759 | { |
760 | xx[i] = *(&x+i); | |
761 | yy[i] = *(&y+i); | |
762 | zz[i] = *(&z+i); | |
fd30f88e | 763 | str += zz[i]; |
3edbbba2 | 764 | } |
765 | for(i=0; i<=nclust; i++) | |
766 | { | |
767 | xxc[i] = *(&xc+i); | |
768 | yyc[i] = *(&yc+i); | |
769 | zzc[i] = *(&zc+i); | |
fd30f88e | 770 | str1 += zzc[i]; |
3edbbba2 | 771 | rrc[i] = 0.5; |
772 | } | |
fd30f88e | 773 | for(i = 0; i <= nclust; i++) |
3edbbba2 | 774 | { |
775 | zzc[i] = str/str1*zzc[i]; | |
776 | ha[i] = xxc[i]; | |
777 | hb[i] = yyc[i]; | |
778 | hc[i] = zzc[i]; | |
779 | hd[i] = rrc[i]; | |
780 | x1 = xxc[i]; | |
781 | y1 = yyc[i]; | |
782 | } | |
562718f9 | 783 | |
fd30f88e | 784 | for(i = 0; i <= ncell; i++) |
785 | { | |
786 | idd = 0; | |
787 | x1 = xx[i]; | |
788 | y1 = yy[i]; | |
789 | for(j = 0; j <= nclust; j++) | |
790 | { | |
791 | x2 = xxc[j]; | |
792 | y2 = yyc[j]; | |
793 | if(Distance(x1,y1,x2,y2) <= 3.) | |
794 | { | |
795 | idd++; | |
796 | neib[i][idd] = j; | |
797 | } | |
798 | } | |
799 | neib[i][0] = idd; | |
3edbbba2 | 800 | } |
fd30f88e | 801 | sum = 0.; |
802 | for(i1 = 0; i1 <= ncell; i1++) | |
803 | { | |
804 | aint = 0.; | |
805 | idd = neib[i1][0]; | |
806 | for(i2 = 1; i2 <= idd; i2++) | |
807 | { | |
808 | jj = neib[i1][i2]; | |
809 | dx = xx[i1] - xxc[jj]; | |
810 | dy = yy[i1] - yyc[jj]; | |
811 | dum = rrc[j]*rrc[jj] + rr*rr; | |
812 | aint += exp(-(dx*dx+dy*dy)/dum)*zzc[idd]*rr*rr/dum; | |
813 | } | |
814 | sum += (aint - zz[i1])*(aint - zz[i1])/str; | |
815 | } | |
816 | str1 = 0.; | |
562718f9 | 817 | |
fd30f88e | 818 | for(i = 0; i <= nclust; i++) |
819 | { | |
820 | a[i] = xxc[i] + 0.6*(rnd.Uniform() - 0.5); | |
821 | b[i] = yyc[i] + 0.6*(rnd.Uniform() - 0.5); | |
822 | c[i] = zzc[i]*(1.+ ( rnd.Uniform() - 0.5)*0.2); | |
823 | str1 += zzc[i]; | |
824 | d[i] = rrc[i]*(1.+ ( rnd.Uniform() - 0.5)*0.1); | |
825 | ||
826 | if(d[i] < 0.25) | |
827 | { | |
828 | d[i]=0.25; | |
829 | } | |
3edbbba2 | 830 | } |
fd30f88e | 831 | for(i = 0; i <= nclust; i++) |
832 | { | |
833 | c[i] = c[i]*str/str1; | |
3edbbba2 | 834 | } |
fd30f88e | 835 | sum1=0.; |
562718f9 | 836 | |
fd30f88e | 837 | for(i1 = 0; i1 <= ncell; i1++) |
838 | { | |
839 | aint = 0.; | |
840 | idd = neib[i1][0]; | |
841 | for(i2 = 1; i2 <= idd; i2++) | |
842 | { | |
843 | jj = neib[i1][i2]; | |
844 | dx = xx[i1] - a[jj]; | |
845 | dy = yy[i1] - b[jj]; | |
846 | dum = d[jj]*d[jj]+rr*rr; | |
847 | aint += exp(-(dx*dx+dy*dy)/dum)*c[i2]*rr*rr/dum; | |
848 | } | |
849 | sum1 += (aint - zz[i1])*(aint - zz[i1])/str; | |
3edbbba2 | 850 | } |
851 | ||
fd30f88e | 852 | if(sum1 < sum) |
853 | { | |
854 | for(i2 = 0; i2 <= nclust; i2++) | |
855 | { | |
856 | xxc[i2] = a[i2]; | |
857 | yyc[i2] = b[i2]; | |
858 | zzc[i2] = c[i2]; | |
859 | rrc[i2] = d[i2]; | |
860 | sum = sum1; | |
861 | } | |
862 | } | |
863 | for(j = 0; j <= nclust; j++) | |
864 | { | |
865 | *(&xc+j) = xxc[j]; | |
866 | *(&yc+j) = yyc[j]; | |
867 | *(&zc+j) = zzc[j]; | |
868 | *(&rc+j) = rrc[j]; | |
3edbbba2 | 869 | } |
3edbbba2 | 870 | } |
871 | // ------------------------------------------------------------------------ // | |
fd30f88e | 872 | Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1, |
873 | Double_t x2, Double_t y2) | |
3edbbba2 | 874 | { |
fd30f88e | 875 | return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); |
3edbbba2 | 876 | } |
877 | // ------------------------------------------------------------------------ // | |
878 | void AliPMDClusteringV1::SetEdepCut(Float_t decut) | |
879 | { | |
880 | fCutoff = decut; | |
881 | } | |
882 | // ------------------------------------------------------------------------ // |