Update of ExB efect simulation for SPD (G. Bruno)
[u/mrichter/AliRoot.git] / PMD / AliPMDClusteringV1.cxx
CommitLineData
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
2c1131dd 57
3edbbba2 58ClassImp(AliPMDClusteringV1)
59
60const Double_t AliPMDClusteringV1::fgkSqroot3by2=0.8660254; // sqrt(3.)/2.
61
62AliPMDClusteringV1::AliPMDClusteringV1():
562718f9 63 fPMDclucont(new TObjArray()),
4755c3f1 64 fCutoff(0.0),
65 fClusParam(0)
3edbbba2 66{
fd30f88e 67 for(Int_t i = 0; i < kNDIMX; i++)
3edbbba2 68 {
fd30f88e 69 for(Int_t j = 0; j < kNDIMY; j++)
3edbbba2 70 {
71 fCoord[0][i][j] = i+j/2.;
72 fCoord[1][i][j] = fgkSqroot3by2*j;
3edbbba2 73 }
74 }
75}
76// ------------------------------------------------------------------------ //
562718f9 77AliPMDClusteringV1::AliPMDClusteringV1(const AliPMDClusteringV1& pmdclv1):
78 AliPMDClustering(pmdclv1),
79 fPMDclucont(0),
4755c3f1 80 fCutoff(0),
81 fClusParam(0)
562718f9 82{
83 // copy constructor
84 AliError("Copy constructor not allowed ");
85
86}
87// ------------------------------------------------------------------------ //
88AliPMDClusteringV1 &AliPMDClusteringV1::operator=(const AliPMDClusteringV1& /*pmdclv1*/)
89{
90 // copy constructor
91 AliError("Assignment operator not allowed ");
92 return *this;
93}
94// ------------------------------------------------------------------------ //
3edbbba2 95AliPMDClusteringV1::~AliPMDClusteringV1()
96{
562718f9 97 delete fPMDclucont;
3edbbba2 98}
99// ------------------------------------------------------------------------ //
920e13db 100void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn,
101 Int_t celltrack[48][96],
102 Int_t cellpid[48][96],
103 Double_t celladc[48][96],
22bd512d 104 TObjArray *pmdcont)
3edbbba2 105{
106 // main function to call other necessary functions to do clustering
107 //
01c4d84a 108
562718f9 109 AliPMDcluster *pmdcl = 0;
3edbbba2 110
2c1131dd 111 const float ktwobysqrt3 = 1.1547; // 2./sqrt(3.)
112 const Int_t kNmaxCell = 19; // # of cells surrounding a cluster center
113
be8b7039 114 Int_t i = 0, j = 0, nmx1 = 0;
115 Int_t incr = 0, id = 0, jd = 0;
2c1131dd 116 Int_t celldataX[kNmaxCell], celldataY[kNmaxCell];
920e13db 117 Int_t celldataTr[kNmaxCell], celldataPid[kNmaxCell];
b0e4d1e1 118 Float_t celldataAdc[kNmaxCell];
be8b7039 119 Float_t clusdata[6] = {0.,0.,0.,0.,0.,0.};
562718f9 120 Double_t cutoff, ave;
121 Double_t edepcell[kNMX];
4e2efa25 122 Double_t cellenergy[kNMX];
9c294f00 123
01c4d84a 124 // ndimXr and ndimYr are different because of different module size
125
fd30f88e 126 Int_t ndimXr = 0;
127 Int_t ndimYr = 0;
01c4d84a 128
129 if (ismn < 12)
130 {
131 ndimXr = 96;
132 ndimYr = 48;
133 }
134 else if (ismn >= 12 && ismn <= 23)
135 {
136 ndimXr = 48;
137 ndimYr = 96;
138 }
139
4e2efa25 140 for (i = 0; i < kNMX; i++)
9c294f00 141 {
142 cellenergy[i] = 0.;
143 }
144
562718f9 145 Int_t kk = 0;
78fc1b96 146 for (i = 0; i < kNDIMX; i++)
3edbbba2 147 {
78fc1b96 148 for (j = 0; j < kNDIMY; j++)
01c4d84a 149 {
562718f9 150 edepcell[kk] = 0.;
151 kk++;
01c4d84a 152 }
153 }
154
155 for (id = 0; id < ndimXr; id++)
156 {
157 for (jd = 0; jd < ndimYr; jd++)
3edbbba2 158 {
fd30f88e 159 j = jd;
160 i = id+(ndimYr/2-1)-(jd/2);
01c4d84a 161
9c294f00 162 Int_t ij = i + j*kNDIMX;
9c294f00 163
01c4d84a 164 if (ismn < 12)
165 {
2c1131dd 166 cellenergy[ij] = celladc[jd][id];
01c4d84a 167 }
168 else if (ismn >= 12 && ismn <= 23)
169 {
2c1131dd 170 cellenergy[ij] = celladc[id][jd];
01c4d84a 171 }
3edbbba2 172 }
173 }
562718f9 174
78fc1b96 175 for (i = 0; i < kNMX; i++)
939c3b8e 176 {
177 edepcell[i] = cellenergy[i];
178 }
9c294f00 179
4e2efa25 180 Bool_t jsort = true;
562718f9 181 Int_t iord1[kNMX];
4e2efa25 182 TMath::Sort((Int_t)kNMX,edepcell,iord1,jsort);// order the data
183 cutoff = fCutoff; // cutoff to discard cells
562718f9 184 ave = 0.;
fd30f88e 185 nmx1 = -1;
562718f9 186 for(i = 0;i < kNMX; i++)
3edbbba2 187 {
562718f9 188 if(edepcell[i] > 0.)
fd30f88e 189 {
562718f9 190 ave += edepcell[i];
fd30f88e 191 }
562718f9 192 if(edepcell[i] > cutoff )
fd30f88e 193 {
194 nmx1++;
195 }
3edbbba2 196 }
9c294f00 197
3edbbba2 198 AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
199
3edbbba2 200 if (nmx1 == 0) nmx1 = 1;
fd30f88e 201 ave = ave/nmx1;
3edbbba2 202 AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
203 kNMX,ave));
562718f9 204
205 incr = CrClust(ave, cutoff, nmx1,iord1, edepcell );
562718f9 206 RefClust(incr,edepcell);
207 Int_t nentries1 = fPMDclucont->GetEntries();
fd30f88e 208 AliDebug(1,Form("Detector Plane = %d Serial Module No = %d Number of clusters = %d",idet, ismn, nentries1));
209 AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
562718f9 210
fd30f88e 211 for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
3edbbba2 212 {
fd30f88e 213 AliPMDcludata *pmdcludata =
562718f9 214 (AliPMDcludata*)fPMDclucont->UncheckedAt(ient1);
fd30f88e 215 Float_t cluXC = pmdcludata->GetClusX();
216 Float_t cluYC = pmdcludata->GetClusY();
217 Float_t cluADC = pmdcludata->GetClusADC();
218 Float_t cluCELLS = pmdcludata->GetClusCells();
219 Float_t cluSIGX = pmdcludata->GetClusSigmaX();
220 Float_t cluSIGY = pmdcludata->GetClusSigmaY();
221
3edbbba2 222 Float_t cluY0 = ktwobysqrt3*cluYC;
223 Float_t cluX0 = cluXC - cluY0/2.;
fd30f88e 224
3edbbba2 225 //
226 // Cluster X centroid is back transformed
227 //
2c1131dd 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 }
2c1131dd 237
01c4d84a 238 clusdata[1] = cluY0;
239 clusdata[2] = cluADC;
240 clusdata[3] = cluCELLS;
fd30f88e 241 clusdata[4] = cluSIGX;
242 clusdata[5] = cluSIGY;
e6ba3040 243
3edbbba2 244 //
245 // Cells associated with a cluster
246 //
562718f9 247
2c1131dd 248 for (Int_t ihit = 0; ihit < kNmaxCell; ihit++)
3edbbba2 249 {
2c1131dd 250 Int_t cellrow = pmdcludata->GetCellXY(ihit)/10000;
251 Int_t cellcol = pmdcludata->GetCellXY(ihit)%10000;
252
01c4d84a 253 if (ismn < 12)
254 {
2c1131dd 255 celldataX[ihit] = cellrow - (24-1) + int(cellcol/2.);
01c4d84a 256 }
257 else if (ismn >= 12 && ismn <= 23)
258 {
2c1131dd 259 celldataX[ihit] = cellrow - (48-1) + int(cellcol/2.);
01c4d84a 260 }
2c1131dd 261
262 celldataY[ihit] = cellcol;
920e13db 263
264 Int_t irow = celldataX[ihit];
265 Int_t icol = celldataY[ihit];
2c1131dd 266
d270ca46 267 if (ismn < 12)
920e13db 268 {
d270ca46 269 if ((irow >= 0 && irow < 96) && (icol >= 0 && icol < 48))
270 {
271 celldataTr[ihit] = celltrack[icol][irow];
272 celldataPid[ihit] = cellpid[icol][irow];
273 celldataAdc[ihit] = (Float_t) celladc[icol][irow];
274 }
275 else
276 {
277 celldataTr[ihit] = -1;
278 celldataPid[ihit] = -1;
279 celldataAdc[ihit] = -1;
280 }
920e13db 281 }
d270ca46 282 else if (ismn >= 12 && ismn < 24)
920e13db 283 {
d270ca46 284 if ((irow >= 0 && irow < 48) && (icol >= 0 && icol < 96))
285 {
286 celldataTr[ihit] = celltrack[irow][icol];
287 celldataPid[ihit] = cellpid[irow][icol];
288 celldataAdc[ihit] = (Float_t) celladc[irow][icol];
e6ba3040 289
d270ca46 290 }
291 else
292 {
293 celldataTr[ihit] = -1;
294 celldataPid[ihit] = -1;
295 celldataAdc[ihit] = -1;
296 }
920e13db 297 }
d270ca46 298
920e13db 299 }
e6ba3040 300
920e13db 301 pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY,
b0e4d1e1 302 celldataTr, celldataPid, celldataAdc);
3edbbba2 303 pmdcont->Add(pmdcl);
304 }
fd30f88e 305
91e6e2a0 306 fPMDclucont->Delete();
3edbbba2 307}
308// ------------------------------------------------------------------------ //
562718f9 309Int_t AliPMDClusteringV1::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
310 Int_t iord1[], Double_t edepcell[])
3edbbba2 311{
fd30f88e 312 // Does crude clustering
3edbbba2 313 // Finds out only the big patch by just searching the
314 // connected cells
315 //
2c1131dd 316 const Int_t kndim = 4609;
be8b7039 317 Int_t i=0,j=0,k=0,id1=0,id2=0,icl=0, numcell=0;
318 Int_t jd1=0,jd2=0, icell=0, cellcount=0;
319 Int_t clust[2][kndim];
fd30f88e 320 static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
321
3edbbba2 322 AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
323
fd30f88e 324 for (j = 0; j < kNDIMX; j++)
325 {
326 for(k = 0; k < kNDIMY; k++)
327 {
328 fInfocl[0][j][k] = 0;
329 fInfocl[1][j][k] = 0;
330 }
331 }
332 for(i=0; i < kNMX; i++)
333 {
334 fInfcl[0][i] = -1;
562718f9 335
336 j = iord1[i];
337 id2 = j/kNDIMX;
338 id1 = j-id2*kNDIMX;
339
340 if(edepcell[j] <= cutoff)
fd30f88e 341 {
342 fInfocl[0][id1][id2] = -1;
343 }
3edbbba2 344 }
562718f9 345
3edbbba2 346 // ---------------------------------------------------------------
347 // crude clustering begins. Start with cell having largest adc
348 // count and loop over the cells in descending order of adc count
349 // ---------------------------------------------------------------
562718f9 350
fd30f88e 351 icl = -1;
352 cellcount = -1;
3edbbba2 353
fd30f88e 354 for(icell = 0; icell <= nmx1; icell++)
355 {
562718f9 356 j = iord1[icell];
357 id2 = j/kNDIMX;
358 id1 = j-id2*kNDIMX;
359
fd30f88e 360 if(fInfocl[0][id1][id2] == 0 )
361 {
362 icl++;
363 numcell = 0;
364 cellcount++;
365 fInfocl[0][id1][id2] = 1;
366 fInfocl[1][id1][id2] = icl;
367 fInfcl[0][cellcount] = icl;
368 fInfcl[1][cellcount] = id1;
369 fInfcl[2][cellcount] = id2;
370
371 clust[0][numcell] = id1;
372 clust[1][numcell] = id2;
373
2c1131dd 374 for(i = 1; i < kndim; i++)
fd30f88e 375 {
376 clust[0][i]=0;
377 }
378 // ---------------------------------------------------------------
379 // check for adc count in neib. cells. If ne 0 put it in this clust
380 // ---------------------------------------------------------------
381 for(i = 0; i < 6; i++)
382 {
383 jd1 = id1 + neibx[i];
384 jd2 = id2 + neiby[i];
385 if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
386 fInfocl[0][jd1][jd2] == 0)
387 {
388 numcell++;
389 fInfocl[0][jd1][jd2] = 2;
390 fInfocl[1][jd1][jd2] = icl;
391 clust[0][numcell] = jd1;
392 clust[1][numcell] = jd2;
393 cellcount++;
394 fInfcl[0][cellcount] = icl;
395 fInfcl[1][cellcount] = jd1;
396 fInfcl[2][cellcount] = jd2;
397 }
398 }
399 // ---------------------------------------------------------------
400 // check adc count for neighbour's neighbours recursively and
401 // if nonzero, add these to the cluster.
402 // ---------------------------------------------------------------
2c1131dd 403 for(i = 1; i < kndim;i++)
fd30f88e 404 {
405 if(clust[0][i] != 0)
406 {
407 id1 = clust[0][i];
408 id2 = clust[1][i];
409 for(j = 0; j < 6 ; j++)
410 {
411 jd1 = id1 + neibx[j];
412 jd2 = id2 + neiby[j];
413 if( (jd1 >= 0 && jd1 < kNDIMX) &&
414 (jd2 >= 0 && jd2 < kNDIMY) &&
415 fInfocl[0][jd1][jd2] == 0 )
416 {
417 fInfocl[0][jd1][jd2] = 2;
418 fInfocl[1][jd1][jd2] = icl;
419 numcell++;
420 clust[0][numcell] = jd1;
421 clust[1][numcell] = jd2;
422 cellcount++;
423 fInfcl[0][cellcount] = icl;
424 fInfcl[1][cellcount] = jd1;
425 fInfcl[2][cellcount] = jd2;
426 }
427 }
428 }
3edbbba2 429 }
3edbbba2 430 }
3edbbba2 431 }
3edbbba2 432 return cellcount;
433}
434// ------------------------------------------------------------------------ //
562718f9 435void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
3edbbba2 436{
437 // Does the refining of clusters
438 // Takes the big patch and does gaussian fitting and
439 // finds out the more refined clusters
440 //
22d01768 441
22d01768 442 AliPMDcludata *pmdcludata = 0;
2c1131dd 443
444 const Int_t kNmaxCell = 19; // # of cells surrounding a cluster center
445
446 Int_t ndim = incr + 1;
447
448 Int_t *ncl = 0x0;
449 Int_t *clxy = 0x0;
be8b7039 450 Int_t i12 = 0, i22 = 0;
451 Int_t i = 0, j = 0, k = 0;
452 Int_t i1 = 0, i2 = 0, id = 0, icl = 0;
453 Int_t itest = 0, ihld = 0, ig = 0;
454 Int_t nsupcl = 0, clno = 0, t1 = 0, t2 = 0;
562718f9 455 Float_t clusdata[6];
be8b7039 456 Double_t x1 = 0, y1 = 0, z1 = 0, x2 = 0, y2 = 0, z2 = 0, rr = 0;
2c1131dd 457
458 ncl = new Int_t [ndim];
459 clxy = new Int_t [kNmaxCell];
460
562718f9 461 // Initialisation
2c1131dd 462 for(i = 0; i<ndim; i++)
463 {
464 ncl[i] = -1;
4e2efa25 465 }
466 for(i = 0; i<6; i++)
467 {
468 clusdata[i] = 0.;
469 }
470 for(i = 0; i<19; i++)
471 {
472 clxy[i] = 0;
fd30f88e 473 }
562718f9 474
fd30f88e 475 // clno counts the final clusters
3edbbba2 476 // nsupcl = # of superclusters; ncl[i]= # of cells in supercluster i
477 // x, y and z store (x,y) coordinates of and energy deposited in a cell
478 // xc, yc store (x,y) coordinates of the cluster center
479 // zc stores the energy deposited in a cluster
480 // rc is cluster radius
fd30f88e 481
482 clno = -1;
3edbbba2 483 nsupcl = -1;
562718f9 484
fd30f88e 485 for(i = 0; i <= incr; i++)
486 {
487 if(fInfcl[0][i] != nsupcl)
01c4d84a 488 {
fd30f88e 489 nsupcl++;
01c4d84a 490 }
2c1131dd 491 if (nsupcl > ndim)
01c4d84a 492 {
fd30f88e 493 AliWarning("RefClust: Too many superclusters!");
2c1131dd 494 nsupcl = ndim;
fd30f88e 495 break;
01c4d84a 496 }
fd30f88e 497 ncl[nsupcl]++;
3edbbba2 498 }
fd30f88e 499
500 AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
501 incr+1,nsupcl+1));
502 id = -1;
503 icl = -1;
562718f9 504
fd30f88e 505 for(i = 0; i <= nsupcl; i++)
506 {
507 if(ncl[i] == 0)
01c4d84a 508 {
fd30f88e 509 id++;
510 icl++;
2c1131dd 511 if (clno >= 4608)
fd30f88e 512 {
2c1131dd 513 AliWarning("RefClust: Too many clusters! more than 4608");
fd30f88e 514 return;
515 }
516 clno++;
517 i1 = fInfcl[1][id];
518 i2 = fInfcl[2][id];
519
78fc1b96 520 i12 = i1 + i2*kNDIMX;
562718f9 521
fd30f88e 522 clusdata[0] = fCoord[0][i1][i2];
523 clusdata[1] = fCoord[1][i1][i2];
562718f9 524 clusdata[2] = edepcell[i12];
fd30f88e 525 clusdata[3] = 1.;
22bd512d 526 clusdata[4] = 999.5;
527 clusdata[5] = 999.5;
562718f9 528
2c1131dd 529 clxy[0] = i1*10000 + i2;
530
531 for(Int_t icltr = 1; icltr < kNmaxCell; icltr++)
fd30f88e 532 {
562718f9 533 clxy[icltr] = -1;
3edbbba2 534 }
e6ba3040 535
562718f9 536 pmdcludata = new AliPMDcludata(clusdata,clxy);
537 fPMDclucont->Add(pmdcludata);
3edbbba2 538 }
fd30f88e 539 else if(ncl[i] == 1)
540 {
541 id++;
542 icl++;
2c1131dd 543 if (clno >= 4608)
fd30f88e 544 {
2c1131dd 545 AliWarning("RefClust: Too many clusters! more than 4608");
fd30f88e 546 return;
547 }
548 clno++;
549 i1 = fInfcl[1][id];
550 i2 = fInfcl[2][id];
78fc1b96 551 i12 = i1 + i2*kNDIMX;
562718f9 552
fd30f88e 553 x1 = fCoord[0][i1][i2];
554 y1 = fCoord[1][i1][i2];
562718f9 555 z1 = edepcell[i12];
2c1131dd 556
557 clxy[0] = i1*10000 + i2;
558
562718f9 559 id++;
fd30f88e 560 i1 = fInfcl[1][id];
561 i2 = fInfcl[2][id];
562718f9 562
2c1131dd 563 i22 = i1 + i2*kNDIMX;
fd30f88e 564 x2 = fCoord[0][i1][i2];
565 y2 = fCoord[1][i1][i2];
562718f9 566 z2 = edepcell[i22];
2c1131dd 567
568 clxy[1] = i1*10000 + i2;
569
570
571 for(Int_t icltr = 2; icltr < kNmaxCell; icltr++)
fd30f88e 572 {
562718f9 573 clxy[icltr] = -1;
fd30f88e 574 }
fd30f88e 575
576 clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
577 clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
578 clusdata[2] = z1+z2;
579 clusdata[3] = 2.;
580 clusdata[4] = 0.5;
581 clusdata[5] = 0.0;
562718f9 582 pmdcludata = new AliPMDcludata(clusdata,clxy);
583 fPMDclucont->Add(pmdcludata);
3edbbba2 584 }
fd30f88e 585 else
586 {
2c1131dd 587
588 Int_t *iord, *tc, *t;
589 Double_t *x, *y, *z, *xc, *yc, *zc;
590
591 iord = new Int_t [ncl[i]+1];
592 tc = new Int_t [ncl[i]+1];
593 t = new Int_t [ncl[i]+1];
594
595 x = new Double_t [ncl[i]+1];
596 y = new Double_t [ncl[i]+1];
597 z = new Double_t [ncl[i]+1];
598 xc = new Double_t [ncl[i]+1];
599 yc = new Double_t [ncl[i]+1];
600 zc = new Double_t [ncl[i]+1];
601
602 for( k = 0; k < ncl[i]+1; k++)
603 {
604 iord[k] = -1;
605 t[k] = -1;
606 tc[k] = -1;
607 x[k] = -1;
608 y[k] = -1;
609 z[k] = -1;
610 xc[k] = -1;
611 yc[k] = -1;
612 zc[k] = -1;
613 }
fd30f88e 614 id++;
fd30f88e 615 // super-cluster of more than two cells - broken up into smaller
616 // clusters gaussian centers computed. (peaks separated by > 1 cell)
617 // Begin from cell having largest energy deposited This is first
618 // cluster center
619 i1 = fInfcl[1][id];
620 i2 = fInfcl[2][id];
78fc1b96 621 i12 = i1 + i2*kNDIMX;
562718f9 622
fd30f88e 623 x[0] = fCoord[0][i1][i2];
624 y[0] = fCoord[1][i1][i2];
562718f9 625 z[0] = edepcell[i12];
2c1131dd 626 t[0] = i1*10000 + i2;
fd30f88e 627
2c1131dd 628
fd30f88e 629 iord[0] = 0;
630 for(j = 1; j <= ncl[i]; j++)
631 {
632 id++;
633 i1 = fInfcl[1][id];
634 i2 = fInfcl[2][id];
78fc1b96 635 i12 = i1 + i2*kNDIMX;
562718f9 636
fd30f88e 637 iord[j] = j;
638 x[j] = fCoord[0][i1][i2];
639 y[j] = fCoord[1][i1][i2];
562718f9 640 z[j] = edepcell[i12];
2c1131dd 641 t[j] = i1*10000 + i2;
562718f9 642
3edbbba2 643 }
fd30f88e 644
fd30f88e 645 // arranging cells within supercluster in decreasing order
2c1131dd 646
fd30f88e 647 for(j = 1;j <= ncl[i]; j++)
648 {
649 itest = 0;
650 ihld = iord[j];
651 for(i1 = 0; i1 < j; i1++)
652 {
653 if(itest == 0 && z[iord[i1]] < z[ihld])
654 {
655 itest = 1;
656 for(i2 = j-1; i2 >= i1; i2--)
657 {
658 iord[i2+1] = iord[i2];
659 }
660 iord[i1] = ihld;
661 }
662 }
3edbbba2 663 }
b634e7ee 664
665 if (fClusParam == 1)
666 {
667 // Clustering algorithm returns from here for PP collisions
668 // for pp, only the output of crude clusterng is taken
669 // sigx and sigy are not calculated at this moment
670
671 Double_t supx=0., supy=0., supz=0.;
672
673 for(j = 0;j <= ncl[i]; j++)
674 {
675 supx += x[iord[j]]*z[iord[j]];
676 supy += y[iord[j]]*z[iord[j]];
677 supz += z[iord[j]];
678 if(j < 19)
679 {
680 clxy[j] = t[iord[j]];
681 }
682 }
683
684 if( ncl[i] + 1 < 19)
685 {
686 for(Int_t ncel = ncl[i] + 1; ncel < kNmaxCell; ncel ++ )
687 {
688 clxy[ncel] = -1;
689 }
690 }
691 clusdata[0] = supx/supz;
692 clusdata[1] = supy/supz;
693 clusdata[2] = supz;
694 clusdata[3] = ncl[i]+1;
695 clusdata[4] = 0.5;
696 clusdata[5] = 0.0;
697 pmdcludata = new AliPMDcludata(clusdata,clxy);
698 fPMDclucont->Add(pmdcludata);
699 }
700
2c1131dd 701 /* MODIFICATION PART STARTS (Tapan July 2008)
702 iord[0] is the cell with highest ADC in the crude-cluster
703 ig is the number of local maxima in the crude-cluster
704 For the higest peak we make ig=0 which means first local maximum.
705 Next we go down in terms of the ADC sequence and find out if any
706 more of the cells form local maxima. The definition of local
707 maxima is that all its neighbours are of less ADC compared to it.
708 */
b634e7ee 709
710 if (fClusParam == 2)
fd30f88e 711 {
b634e7ee 712 // This part is to split the supercluster
713 //
714 ig = 0;
715 xc[ig] = x[iord[0]];
716 yc[ig] = y[iord[0]];
717 zc[ig] = z[iord[0]];
718 tc[ig] = t[iord[0]];
719 Int_t ivalid = 0, icount = 0;
720
721 for(j=1;j<=ncl[i];j++)
fd30f88e 722 {
b634e7ee 723 x1 = x[iord[j]];
724 y1 = y[iord[j]];
725 z1 = z[iord[j]];
726 t1 = t[iord[j]];
727 rr=Distance(x1,y1,xc[ig],yc[ig]);
728
729 // Check the cells which are outside the neighbours (rr>1.2)
730 if(rr>1.2 )
2c1131dd 731 {
b634e7ee 732 ivalid=0;
733 icount=0;
734 for(Int_t j1=1;j1<j;j1++)
735 {
736 icount++;
737 Float_t rr1=Distance(x1,y1,x[iord[j1]],y[iord[j1]]);
738 if(rr1>1.2) ivalid++;
739 }
740 if(ivalid == icount && z1>0.5*zc[ig])
741 {
742 ig++;
743 xc[ig]=x1;
744 yc[ig]=y1;
745 zc[ig]=z1;
746 tc[ig]=t1;
747 }
748 }
749 }
2c1131dd 750
b634e7ee 751 icl=icl+ig+1;
2c1131dd 752
b634e7ee 753 // We use simple Gaussian weighting. (Tapan Jan 2005)
754 // compute the number of cells belonging to each cluster.
755 // cell can be shared between several clusters
756 // in the ratio of cluster energy deposition
757 // To calculate:
758 // (1) number of cells belonging to a cluster (ig) and
759 // (2) total ADC of the cluster (ig)
760 // (3) x and y positions of the cluster
2c1131dd 761
b634e7ee 762
763 Int_t *cellCount;
764 Int_t **cellXY;
765
766 Int_t *status;
767 Double_t *totaladc, *totaladc2, *ncell,*weight;
768 Double_t *xclust, *yclust, *sigxclust, *sigyclust;
769 Double_t *ax, *ay, *ax2, *ay2;
770
771
772 status = new Int_t [ncl[i]+1];
773 cellXY = new Int_t *[ncl[i]+1];
774
775 cellCount = new Int_t [ig+1];
776 totaladc = new Double_t [ig+1];
777 totaladc2 = new Double_t [ig+1];
778 ncell = new Double_t [ig+1];
779 weight = new Double_t [ig+1];
780 xclust = new Double_t [ig+1];
781 yclust = new Double_t [ig+1];
782 sigxclust = new Double_t [ig+1];
783 sigyclust = new Double_t [ig+1];
784 ax = new Double_t [ig+1];
785 ay = new Double_t [ig+1];
786 ax2 = new Double_t [ig+1];
787 ay2 = new Double_t [ig+1];
788
789 for(j = 0; j < ncl[i]+1; j++)
2c1131dd 790 {
b634e7ee 791 status[j] = 0;
792 cellXY[j] = new Int_t[ig+1];
2c1131dd 793 }
b634e7ee 794 //initialization
795 for(Int_t kcl = 0; kcl < ig+1; kcl++)
fd30f88e 796 {
b634e7ee 797 cellCount[kcl] = 0;
798 totaladc[kcl] = 0.;
799 totaladc2[kcl] = 0.;
800 ncell[kcl] = 0.;
801 weight[kcl] = 0.;
802 xclust[kcl] = 0.;
803 yclust[kcl] = 0.;
804 sigxclust[kcl] = 0.;
805 sigyclust[kcl] = 0.;
806 ax[kcl] = 0.;
807 ay[kcl] = 0.;
808 ax2[kcl] = 0.;
809 ay2[kcl] = 0.;
810 for(j = 0; j < ncl[i]+1; j++)
fd30f88e 811 {
b634e7ee 812 cellXY[j][kcl] = 0;
fd30f88e 813 }
2c1131dd 814 }
b634e7ee 815 Double_t sumweight, gweight;
2c1131dd 816
b634e7ee 817 for(j = 0;j <= ncl[i]; j++)
818 {
819 x1 = x[iord[j]];
2c1131dd 820 y1 = y[iord[j]];
821 z1 = z[iord[j]];
822 t1 = t[iord[j]];
b634e7ee 823
824 for(Int_t kcl=0; kcl<=ig; kcl++)
fd30f88e 825 {
b634e7ee 826 x2 = xc[kcl];
827 y2 = yc[kcl];
2c1131dd 828 rr = Distance(x1,y1,x2,y2);
b634e7ee 829 t2 = tc[kcl];
2c1131dd 830
b634e7ee 831 if(rr==0)
fd30f88e 832 {
b634e7ee 833 ncell[kcl] = 1.;
834 totaladc[kcl] = z1;
835 totaladc2[kcl] = z1*z1;
836 ax[kcl] = x1 * z1;
837 ay[kcl] = y1 * z1;
838 ax2[kcl] = 0.;
839 ay2[kcl] = 0.;
840 status[j] = 1;
fd30f88e 841 }
842 }
b634e7ee 843 }
844
845 for(j = 0; j <= ncl[i]; j++)
846 {
847 Int_t maxweight = 0;
848 Double_t max = 0.;
2c1131dd 849
b634e7ee 850 if(status[j] == 0)
851 {
852 x1 = x[iord[j]];
853 y1 = y[iord[j]];
854 z1 = z[iord[j]];
855 t1 = t[iord[j]];
856 sumweight = 0.;
857
858 for(Int_t kcl = 0; kcl <= ig; kcl++)
859 {
860 x2 = xc[kcl];
861 y2 = yc[kcl];
862 rr = Distance(x1,y1,x2,y2);
863 gweight = exp(-(rr*rr)/(2*(1.2*1.2)));
864 weight[kcl] = zc[kcl] * gweight;
865 sumweight = sumweight + weight[kcl];
866
867 if(weight[kcl] > max)
868 {
869 max = weight[kcl];
870 maxweight = kcl;
871 }
872 }
873
874 cellXY[cellCount[maxweight]][maxweight] = iord[j];
875
876 cellCount[maxweight]++;
877
878 x2 = xc[maxweight];
879 y2 = yc[maxweight];
880 totaladc[maxweight] += z1;
881 ax[maxweight] += x1*z1;
882 ay[maxweight] += y1*z1;
883 totaladc2[maxweight] += z1*z1;
884 ax2[maxweight] += z1*(x1-x2)*(x1-x2);
885 ay2[maxweight] += z1*(y1-y2)*(y1-y2);
886 ncell[maxweight]++;
887
888 }
fd30f88e 889 }
fd30f88e 890
b634e7ee 891 for(Int_t kcl = 0; kcl <= ig; kcl++)
939c3b8e 892 {
b634e7ee 893 if(totaladc[kcl] > 0.)
894 {
895 xclust[kcl] = (ax[kcl])/ totaladc[kcl];
896 yclust[kcl] = (ay[kcl])/ totaladc[kcl];
897
898 //natasha
899 Float_t sqtotadc = totaladc[kcl]*totaladc[kcl];
900 if(totaladc2[kcl] >= sqtotadc)
901 {
902 sigxclust[kcl] = 0.25;
903 sigyclust[kcl] = 0.25;
904 }
905 else
906 {
907 sigxclust[kcl] = (totaladc[kcl]/(sqtotadc-totaladc2[kcl]))*ax2[kcl];
908 sigyclust[kcl] = (totaladc[kcl]/(sqtotadc-totaladc2[kcl]))*ay2[kcl];
909 }
910 }
911
912 for(j = 0; j < cellCount[kcl]; j++) clno++;
939c3b8e 913
b634e7ee 914 if (clno >= 4608)
939c3b8e 915 {
b634e7ee 916 AliWarning("RefClust: Too many clusters! more than 4608");
917 return;
918 }
919 clusdata[0] = xclust[kcl];
920 clusdata[1] = yclust[kcl];
921 clusdata[2] = totaladc[kcl];
922 clusdata[3] = ncell[kcl];
923
924 if(sigxclust[kcl] > sigyclust[kcl])
925 {
926 clusdata[4] = TMath::Sqrt(sigxclust[kcl]);
927 clusdata[5] = TMath::Sqrt(sigyclust[kcl]);
939c3b8e 928 }
929 else
930 {
b634e7ee 931 clusdata[4] = TMath::Sqrt(sigyclust[kcl]);
932 clusdata[5] = TMath::Sqrt(sigxclust[kcl]);
2c1131dd 933 }
b634e7ee 934
935 clxy[0] = tc[kcl];
936
937 Int_t Ncell=1;
938 for (Int_t ii = 0; ii < cellCount[kcl]; ii++)
939 {
940 if(ii<18)
941 {
942 clxy[Ncell] = t[cellXY[ii][kcl]];
943 Ncell++;
944 }
945 }
946
947 pmdcludata = new AliPMDcludata(clusdata,clxy);
948 fPMDclucont->Add(pmdcludata);
949 }
950 delete [] cellCount;
951 for(Int_t jj = 0; jj < ncl[i]+1; jj++) delete [] cellXY[jj];
939c3b8e 952
b634e7ee 953 delete [] status;
954 delete [] totaladc;
955 delete [] totaladc2;
956 delete [] ncell;
957 delete [] xclust;
958 delete [] yclust;
959 delete [] sigxclust;
960 delete [] sigyclust;
961 delete [] ax;
962 delete [] ay;
963 delete [] ax2;
964 delete [] ay2;
965 delete [] weight;
966
01c4d84a 967 }
b634e7ee 968
2c1131dd 969 delete [] iord;
970 delete [] tc;
971 delete [] t;
972 delete [] x;
973 delete [] y;
974 delete [] z;
975 delete [] xc;
976 delete [] yc;
977 delete [] zc;
b634e7ee 978
979
fd30f88e 980 }
3edbbba2 981 }
2c1131dd 982 delete [] ncl;
983 delete [] clxy;
3edbbba2 984}
985// ------------------------------------------------------------------------ //
fd30f88e 986Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1,
987 Double_t x2, Double_t y2)
3edbbba2 988{
fd30f88e 989 return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
3edbbba2 990}
991// ------------------------------------------------------------------------ //
992void AliPMDClusteringV1::SetEdepCut(Float_t decut)
993{
994 fCutoff = decut;
995}
996// ------------------------------------------------------------------------ //
4755c3f1 997void AliPMDClusteringV1::SetClusteringParam(Int_t cluspar)
998{
999 fClusParam = cluspar;
1000}
1001// ------------------------------------------------------------------------ //