]>
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 | |
562718f9 | 147 | Int_t ij = i + j*kNDIMX; |
148 | ||
01c4d84a | 149 | if (ismn < 12) |
150 | { | |
562718f9 | 151 | edepcell[ij] = celladc[jd][id]; |
152 | fCellTrNo[i][j] = jd*10000+id; // for association | |
01c4d84a | 153 | } |
154 | else if (ismn >= 12 && ismn <= 23) | |
155 | { | |
562718f9 | 156 | edepcell[ij] = celladc[id][jd]; |
157 | fCellTrNo[i][j] = id*10000+jd; // for association | |
01c4d84a | 158 | } |
159 | ||
3edbbba2 | 160 | } |
161 | } | |
562718f9 | 162 | |
163 | Int_t iord1[kNMX]; | |
164 | TMath::Sort(kNMX,edepcell,iord1);// order the data | |
165 | cutoff = fCutoff; // cutoff to discard cells | |
166 | ave = 0.; | |
fd30f88e | 167 | nmx1 = -1; |
562718f9 | 168 | for(i = 0;i < kNMX; i++) |
3edbbba2 | 169 | { |
562718f9 | 170 | if(edepcell[i] > 0.) |
fd30f88e | 171 | { |
562718f9 | 172 | ave += edepcell[i]; |
fd30f88e | 173 | } |
562718f9 | 174 | if(edepcell[i] > cutoff ) |
fd30f88e | 175 | { |
176 | nmx1++; | |
177 | } | |
3edbbba2 | 178 | } |
3edbbba2 | 179 | |
180 | AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1)); | |
181 | ||
3edbbba2 | 182 | if (nmx1 == 0) nmx1 = 1; |
fd30f88e | 183 | ave = ave/nmx1; |
3edbbba2 | 184 | |
185 | AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f", | |
186 | kNMX,ave)); | |
562718f9 | 187 | |
188 | incr = CrClust(ave, cutoff, nmx1,iord1, edepcell ); | |
189 | ||
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 | // | |
562718f9 | 382 | |
383 | AliPMDcludata *pmdcludata = 0; | |
384 | ||
fd30f88e | 385 | const Int_t kdim = 4500; |
562718f9 | 386 | |
387 | Int_t i, j, k, i1, i2, id, icl, itest,ihld, ig, nsupcl,clno; | |
388 | Int_t lev1[20], lev2[20]; | |
389 | Int_t ncl[kdim], iord[kdim], t[kdim]; | |
390 | Int_t clxy[15]; | |
391 | ||
392 | Int_t *cellCount; | |
393 | Int_t **cellXY; | |
394 | ||
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]; | |
399 | ||
562718f9 | 400 | |
401 | // Initialisation | |
fd30f88e | 402 | for(i = 0; i<kdim; i++) |
403 | { | |
404 | t[i] = -1; | |
562718f9 | 405 | ncl[i] = -1; |
406 | if (i < 6) clusdata[i] = 0.; | |
407 | if (i < 15) clxy[i] = 0; | |
fd30f88e | 408 | } |
562718f9 | 409 | |
fd30f88e | 410 | // clno counts the final clusters |
3edbbba2 | 411 | // nsupcl = # of superclusters; ncl[i]= # of cells in supercluster i |
412 | // x, y and z store (x,y) coordinates of and energy deposited in a cell | |
413 | // xc, yc store (x,y) coordinates of the cluster center | |
414 | // zc stores the energy deposited in a cluster | |
415 | // rc is cluster radius | |
fd30f88e | 416 | |
417 | clno = -1; | |
3edbbba2 | 418 | nsupcl = -1; |
562718f9 | 419 | |
fd30f88e | 420 | for(i = 0; i <= incr; i++) |
421 | { | |
422 | if(fInfcl[0][i] != nsupcl) | |
01c4d84a | 423 | { |
fd30f88e | 424 | nsupcl++; |
01c4d84a | 425 | } |
fd30f88e | 426 | if (nsupcl > kdim) |
01c4d84a | 427 | { |
fd30f88e | 428 | AliWarning("RefClust: Too many superclusters!"); |
429 | nsupcl = kdim; | |
430 | break; | |
01c4d84a | 431 | } |
fd30f88e | 432 | ncl[nsupcl]++; |
3edbbba2 | 433 | } |
fd30f88e | 434 | |
435 | AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d", | |
436 | incr+1,nsupcl+1)); | |
437 | id = -1; | |
438 | icl = -1; | |
562718f9 | 439 | |
fd30f88e | 440 | for(i = 0; i <= nsupcl; i++) |
441 | { | |
442 | if(ncl[i] == 0) | |
01c4d84a | 443 | { |
fd30f88e | 444 | id++; |
445 | icl++; | |
446 | if (clno >= 5000) | |
447 | { | |
448 | AliWarning("RefClust: Too many clusters! more than 5000"); | |
449 | return; | |
450 | } | |
451 | clno++; | |
452 | i1 = fInfcl[1][id]; | |
453 | i2 = fInfcl[2][id]; | |
454 | ||
562718f9 | 455 | Int_t i12 = i1 + i2*kNDIMX; |
456 | ||
fd30f88e | 457 | clusdata[0] = fCoord[0][i1][i2]; |
458 | clusdata[1] = fCoord[1][i1][i2]; | |
562718f9 | 459 | clusdata[2] = edepcell[i12]; |
fd30f88e | 460 | clusdata[3] = 1.; |
461 | clusdata[4] = 0.5; | |
462 | clusdata[5] = 0.0; | |
562718f9 | 463 | |
464 | clxy[0] = fCellTrNo[i1][i2]; //association | |
465 | for(Int_t icltr = 1; icltr < 15; icltr++) | |
fd30f88e | 466 | { |
562718f9 | 467 | clxy[icltr] = -1; |
3edbbba2 | 468 | } |
562718f9 | 469 | pmdcludata = new AliPMDcludata(clusdata,clxy); |
470 | fPMDclucont->Add(pmdcludata); | |
3edbbba2 | 471 | } |
fd30f88e | 472 | else if(ncl[i] == 1) |
473 | { | |
474 | id++; | |
475 | icl++; | |
476 | if (clno >= 5000) | |
477 | { | |
478 | AliWarning("RefClust: Too many clusters! more than 5000"); | |
479 | return; | |
480 | } | |
481 | clno++; | |
482 | i1 = fInfcl[1][id]; | |
483 | i2 = fInfcl[2][id]; | |
562718f9 | 484 | Int_t i12 = i1 + i2*kNDIMX; |
485 | ||
fd30f88e | 486 | x1 = fCoord[0][i1][i2]; |
487 | y1 = fCoord[1][i1][i2]; | |
562718f9 | 488 | z1 = edepcell[i12]; |
489 | clxy[0] = fCellTrNo[i1][i2]; //asso | |
490 | id++; | |
fd30f88e | 491 | i1 = fInfcl[1][id]; |
492 | i2 = fInfcl[2][id]; | |
562718f9 | 493 | |
494 | Int_t i22 = i1 + i2*kNDIMX; | |
fd30f88e | 495 | x2 = fCoord[0][i1][i2]; |
496 | y2 = fCoord[1][i1][i2]; | |
562718f9 | 497 | z2 = edepcell[i22]; |
498 | clxy[1] = fCellTrNo[i1][i2]; //asso | |
499 | for(Int_t icltr = 2; icltr < 15; icltr++) | |
fd30f88e | 500 | { |
562718f9 | 501 | clxy[icltr] = -1; |
fd30f88e | 502 | } |
fd30f88e | 503 | |
504 | clusdata[0] = (x1*z1+x2*z2)/(z1+z2); | |
505 | clusdata[1] = (y1*z1+y2*z2)/(z1+z2); | |
506 | clusdata[2] = z1+z2; | |
507 | clusdata[3] = 2.; | |
508 | clusdata[4] = 0.5; | |
509 | clusdata[5] = 0.0; | |
562718f9 | 510 | pmdcludata = new AliPMDcludata(clusdata,clxy); |
511 | fPMDclucont->Add(pmdcludata); | |
3edbbba2 | 512 | } |
fd30f88e | 513 | else |
514 | { | |
fd30f88e | 515 | id++; |
516 | iord[0] = 0; | |
517 | // super-cluster of more than two cells - broken up into smaller | |
518 | // clusters gaussian centers computed. (peaks separated by > 1 cell) | |
519 | // Begin from cell having largest energy deposited This is first | |
520 | // cluster center | |
521 | i1 = fInfcl[1][id]; | |
522 | i2 = fInfcl[2][id]; | |
562718f9 | 523 | Int_t i12 = i1 + i2*kNDIMX; |
524 | ||
fd30f88e | 525 | x[0] = fCoord[0][i1][i2]; |
526 | y[0] = fCoord[1][i1][i2]; | |
562718f9 | 527 | z[0] = edepcell[i12]; |
528 | ||
529 | t[0] = fCellTrNo[i1][i2]; //asso | |
fd30f88e | 530 | |
531 | iord[0] = 0; | |
532 | for(j = 1; j <= ncl[i]; j++) | |
533 | { | |
534 | id++; | |
535 | i1 = fInfcl[1][id]; | |
536 | i2 = fInfcl[2][id]; | |
562718f9 | 537 | Int_t i12 = i1 + i2*kNDIMX; |
538 | ||
fd30f88e | 539 | iord[j] = j; |
540 | x[j] = fCoord[0][i1][i2]; | |
541 | y[j] = fCoord[1][i1][i2]; | |
562718f9 | 542 | z[j] = edepcell[i12]; |
543 | ||
544 | t[j] = fCellTrNo[i1][i2]; //asso | |
3edbbba2 | 545 | } |
fd30f88e | 546 | |
fd30f88e | 547 | // arranging cells within supercluster in decreasing order |
548 | ||
549 | for(j = 1;j <= ncl[i]; j++) | |
550 | { | |
551 | itest = 0; | |
552 | ihld = iord[j]; | |
553 | for(i1 = 0; i1 < j; i1++) | |
554 | { | |
555 | if(itest == 0 && z[iord[i1]] < z[ihld]) | |
556 | { | |
557 | itest = 1; | |
558 | for(i2 = j-1; i2 >= i1; i2--) | |
559 | { | |
560 | iord[i2+1] = iord[i2]; | |
561 | } | |
562 | iord[i1] = ihld; | |
563 | } | |
564 | } | |
3edbbba2 | 565 | } |
fd30f88e | 566 | // compute the number of Gaussians and their centers ( first |
567 | // guess ) | |
568 | // centers must be separated by cells having smaller ener. dep. | |
569 | // neighbouring centers should be either strong or well-separated | |
562718f9 | 570 | ig = 0; |
fd30f88e | 571 | xc[ig] = x[iord[0]]; |
572 | yc[ig] = y[iord[0]]; | |
573 | zc[ig] = z[iord[0]]; | |
574 | for(j = 1; j <= ncl[i]; j++) | |
575 | { | |
576 | itest = -1; | |
577 | x1 = x[iord[j]]; | |
578 | y1 = y[iord[j]]; | |
579 | for(k = 0; k <= ig; k++) | |
580 | { | |
581 | x2 = xc[k]; | |
582 | y2 = yc[k]; | |
583 | rr = Distance(x1,y1,x2,y2); | |
562718f9 | 584 | if(rr >= 1.1 && rr < 1.8 && z[iord[j]] > zc[k]/4.)itest++; |
585 | if(rr >= 1.8 && rr < 2.1 && z[iord[j]] > zc[k]/10.)itest++; | |
586 | if( rr >= 2.1)itest++; | |
fd30f88e | 587 | } |
588 | if(itest == ig) | |
589 | { | |
590 | ig++; | |
591 | xc[ig] = x1; | |
592 | yc[ig] = y1; | |
593 | zc[ig] = z[iord[j]]; | |
594 | } | |
3edbbba2 | 595 | } |
fd30f88e | 596 | GaussFit(ncl[i], ig, x[0], y[0] ,z[0], xc[0], yc[0], zc[0], rc[0]); |
597 | icl += ig+1; | |
598 | // compute the number of cells belonging to each cluster. | |
599 | // cell is shared between several clusters ( if they are equidistant | |
600 | // from it ) in the ratio of cluster energy deposition | |
562718f9 | 601 | |
602 | Int_t jj = 15; | |
603 | cellCount = new Int_t [ig+1]; | |
604 | cellXY = new Int_t *[jj]; | |
605 | for(Int_t ij = 0; ij < 15; ij++) cellXY[ij] = new Int_t [ig+1]; | |
606 | ||
fd30f88e | 607 | for(j = 0; j <= ig; j++) |
01c4d84a | 608 | { |
562718f9 | 609 | cellCount[j] = 0; |
610 | cells[j] = 0.; | |
fd30f88e | 611 | } |
562718f9 | 612 | |
fd30f88e | 613 | if(ig > 0) |
614 | { | |
615 | for(j = 0; j <= ncl[i]; j++) | |
616 | { | |
617 | lev1[j] = 0; | |
618 | lev2[j] = 0; | |
619 | for(k = 0; k <= ig; k++) | |
620 | { | |
621 | dist = Distance(x[j], y[j], xc[k], yc[k]); | |
622 | if(dist < TMath::Sqrt(3.) ) | |
623 | { | |
624 | //asso | |
562718f9 | 625 | if (cellCount[k] < 15) |
626 | { | |
627 | cellXY[cellCount[k]][k] = t[j]; | |
628 | } | |
fd30f88e | 629 | cellCount[k]++; |
630 | // | |
631 | lev1[0]++; | |
632 | i1 = lev1[0]; | |
633 | lev1[i1] = k; | |
634 | } | |
635 | else | |
636 | { | |
637 | if(dist < 2.1) | |
638 | { | |
639 | lev2[0]++; | |
640 | i1 = lev2[0]; | |
641 | lev2[i1] = k; | |
642 | } | |
643 | } | |
644 | } | |
645 | if(lev1[0] != 0) | |
646 | { | |
647 | if(lev1[0] == 1) | |
648 | { | |
649 | cells[lev1[1]]++; | |
650 | } | |
651 | else | |
652 | { | |
653 | sum=0.; | |
654 | for(k = 1; k <= lev1[0]; k++) | |
655 | { | |
656 | sum += zc[lev1[k]]; | |
657 | } | |
658 | for(k = 1; k <= lev1[0]; k++) | |
659 | { | |
660 | cells[lev1[k]] += zc[lev1[k]]/sum; | |
661 | } | |
662 | } | |
663 | } | |
664 | else | |
665 | { | |
666 | if(lev2[0] == 0) | |
667 | { | |
668 | cells[lev2[1]]++; | |
669 | } | |
670 | else | |
671 | { | |
672 | sum=0.; | |
673 | for( k = 1; k <= lev2[0]; k++) | |
674 | { | |
675 | sum += zc[lev2[k]]; | |
676 | } | |
677 | for(k = 1; k <= lev2[0]; k++) | |
678 | { | |
679 | cells[lev2[k]] += zc[lev2[k]]/sum; | |
680 | } | |
681 | } | |
682 | } | |
683 | } | |
684 | } | |
685 | ||
686 | // zero rest of the cell array | |
687 | //asso | |
688 | for( k = 0; k <= ig; k++) | |
689 | { | |
562718f9 | 690 | for(Int_t icltr = cellCount[k]; icltr < 15; icltr++) |
fd30f88e | 691 | { |
562718f9 | 692 | cellXY[icltr][k] = -1; |
fd30f88e | 693 | } |
694 | } | |
695 | // | |
696 | ||
697 | for(j = 0; j <= ig; j++) | |
698 | { | |
699 | clno++; | |
700 | if (clno >= 5000) | |
701 | { | |
702 | AliWarning("RefClust: Too many clusters! more than 5000"); | |
703 | return; | |
704 | } | |
705 | clusdata[0] = xc[j]; | |
706 | clusdata[1] = yc[j]; | |
707 | clusdata[2] = zc[j]; | |
708 | clusdata[4] = rc[j]; | |
709 | clusdata[5] = 0.0; | |
710 | if(ig == 0) | |
711 | { | |
712 | clusdata[3] = ncl[i]; | |
713 | } | |
714 | else | |
715 | { | |
716 | clusdata[3] = cells[j]; | |
717 | } | |
562718f9 | 718 | |
719 | ||
720 | for (Int_t ii=0; ii < 15; ii++) | |
721 | { | |
722 | clxy[ii] = cellXY[ii][j]; | |
723 | } | |
724 | pmdcludata = new AliPMDcludata(clusdata,clxy); | |
725 | fPMDclucont->Add(pmdcludata); | |
01c4d84a | 726 | } |
562718f9 | 727 | delete [] cellCount; |
728 | for(Int_t jj = 0; jj < 15; jj++)delete [] cellXY[jj]; | |
729 | delete [] cellXY; | |
01c4d84a | 730 | } |
3edbbba2 | 731 | } |
3edbbba2 | 732 | } |
733 | // ------------------------------------------------------------------------ // | |
fd30f88e | 734 | void AliPMDClusteringV1::GaussFit(Int_t ncell, Int_t nclust, Double_t &x, |
735 | Double_t &y ,Double_t &z, Double_t &xc, | |
736 | Double_t &yc, Double_t &zc, Double_t &rc) | |
3edbbba2 | 737 | { |
738 | // Does gaussian fitting | |
739 | // | |
562718f9 | 740 | |
fd30f88e | 741 | const Int_t kdim = 4500; |
742 | Int_t i, j, i1, i2, novar, idd, jj; | |
562718f9 | 743 | Int_t neib[kdim][50]; |
744 | ||
fd30f88e | 745 | Double_t sum, dx, dy, str, str1, aint, sum1, rr, dum; |
746 | Double_t x1, x2, y1, y2; | |
747 | Double_t xx[kdim], yy[kdim], zz[kdim], xxc[kdim], yyc[kdim]; | |
748 | Double_t a[kdim], b[kdim], c[kdim], d[kdim], ha[kdim], hb[kdim]; | |
749 | Double_t hc[kdim], hd[kdim], zzc[kdim], rrc[kdim]; | |
fd30f88e | 750 | |
751 | TRandom rnd; | |
752 | ||
3edbbba2 | 753 | str = 0.; |
754 | str1 = 0.; | |
755 | rr = 0.3; | |
756 | novar = 0; | |
562718f9 | 757 | j = 0; |
3edbbba2 | 758 | |
fd30f88e | 759 | for(i = 0; i <= ncell; i++) |
3edbbba2 | 760 | { |
761 | xx[i] = *(&x+i); | |
762 | yy[i] = *(&y+i); | |
763 | zz[i] = *(&z+i); | |
fd30f88e | 764 | str += zz[i]; |
3edbbba2 | 765 | } |
766 | for(i=0; i<=nclust; i++) | |
767 | { | |
768 | xxc[i] = *(&xc+i); | |
769 | yyc[i] = *(&yc+i); | |
770 | zzc[i] = *(&zc+i); | |
fd30f88e | 771 | str1 += zzc[i]; |
3edbbba2 | 772 | rrc[i] = 0.5; |
773 | } | |
fd30f88e | 774 | for(i = 0; i <= nclust; i++) |
3edbbba2 | 775 | { |
776 | zzc[i] = str/str1*zzc[i]; | |
777 | ha[i] = xxc[i]; | |
778 | hb[i] = yyc[i]; | |
779 | hc[i] = zzc[i]; | |
780 | hd[i] = rrc[i]; | |
781 | x1 = xxc[i]; | |
782 | y1 = yyc[i]; | |
783 | } | |
562718f9 | 784 | |
fd30f88e | 785 | for(i = 0; i <= ncell; i++) |
786 | { | |
787 | idd = 0; | |
788 | x1 = xx[i]; | |
789 | y1 = yy[i]; | |
790 | for(j = 0; j <= nclust; j++) | |
791 | { | |
792 | x2 = xxc[j]; | |
793 | y2 = yyc[j]; | |
794 | if(Distance(x1,y1,x2,y2) <= 3.) | |
795 | { | |
796 | idd++; | |
797 | neib[i][idd] = j; | |
798 | } | |
799 | } | |
800 | neib[i][0] = idd; | |
3edbbba2 | 801 | } |
fd30f88e | 802 | sum = 0.; |
803 | for(i1 = 0; i1 <= ncell; i1++) | |
804 | { | |
805 | aint = 0.; | |
806 | idd = neib[i1][0]; | |
807 | for(i2 = 1; i2 <= idd; i2++) | |
808 | { | |
809 | jj = neib[i1][i2]; | |
810 | dx = xx[i1] - xxc[jj]; | |
811 | dy = yy[i1] - yyc[jj]; | |
812 | dum = rrc[j]*rrc[jj] + rr*rr; | |
813 | aint += exp(-(dx*dx+dy*dy)/dum)*zzc[idd]*rr*rr/dum; | |
814 | } | |
815 | sum += (aint - zz[i1])*(aint - zz[i1])/str; | |
816 | } | |
817 | str1 = 0.; | |
562718f9 | 818 | |
fd30f88e | 819 | for(i = 0; i <= nclust; i++) |
820 | { | |
821 | a[i] = xxc[i] + 0.6*(rnd.Uniform() - 0.5); | |
822 | b[i] = yyc[i] + 0.6*(rnd.Uniform() - 0.5); | |
823 | c[i] = zzc[i]*(1.+ ( rnd.Uniform() - 0.5)*0.2); | |
824 | str1 += zzc[i]; | |
825 | d[i] = rrc[i]*(1.+ ( rnd.Uniform() - 0.5)*0.1); | |
826 | ||
827 | if(d[i] < 0.25) | |
828 | { | |
829 | d[i]=0.25; | |
830 | } | |
3edbbba2 | 831 | } |
fd30f88e | 832 | for(i = 0; i <= nclust; i++) |
833 | { | |
834 | c[i] = c[i]*str/str1; | |
3edbbba2 | 835 | } |
fd30f88e | 836 | sum1=0.; |
562718f9 | 837 | |
fd30f88e | 838 | for(i1 = 0; i1 <= ncell; i1++) |
839 | { | |
840 | aint = 0.; | |
841 | idd = neib[i1][0]; | |
842 | for(i2 = 1; i2 <= idd; i2++) | |
843 | { | |
844 | jj = neib[i1][i2]; | |
845 | dx = xx[i1] - a[jj]; | |
846 | dy = yy[i1] - b[jj]; | |
847 | dum = d[jj]*d[jj]+rr*rr; | |
848 | aint += exp(-(dx*dx+dy*dy)/dum)*c[i2]*rr*rr/dum; | |
849 | } | |
850 | sum1 += (aint - zz[i1])*(aint - zz[i1])/str; | |
3edbbba2 | 851 | } |
852 | ||
fd30f88e | 853 | if(sum1 < sum) |
854 | { | |
855 | for(i2 = 0; i2 <= nclust; i2++) | |
856 | { | |
857 | xxc[i2] = a[i2]; | |
858 | yyc[i2] = b[i2]; | |
859 | zzc[i2] = c[i2]; | |
860 | rrc[i2] = d[i2]; | |
861 | sum = sum1; | |
862 | } | |
863 | } | |
864 | for(j = 0; j <= nclust; j++) | |
865 | { | |
866 | *(&xc+j) = xxc[j]; | |
867 | *(&yc+j) = yyc[j]; | |
868 | *(&zc+j) = zzc[j]; | |
869 | *(&rc+j) = rrc[j]; | |
3edbbba2 | 870 | } |
3edbbba2 | 871 | } |
872 | // ------------------------------------------------------------------------ // | |
fd30f88e | 873 | Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1, |
874 | Double_t x2, Double_t y2) | |
3edbbba2 | 875 | { |
fd30f88e | 876 | return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); |
3edbbba2 | 877 | } |
878 | // ------------------------------------------------------------------------ // | |
879 | void AliPMDClusteringV1::SetEdepCut(Float_t decut) | |
880 | { | |
881 | fCutoff = decut; | |
882 | } | |
883 | // ------------------------------------------------------------------------ // |