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