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