1 /***************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 //-----------------------------------------------------//
20 // Source File : PMDClusteringV2.cxx //
22 // clustering code for alice pmd //
24 //-----------------------------------------------------//
26 /* --------------------------------------------------------------------
27 Code developed by S. C. Phatak, Institute of Physics,
28 Bhubaneswar 751 005 ( phatak@iopb.res.in ) Given the energy deposited
29 ( or ADC value ) in each cell of supermodule ( pmd or cpv ), the code
30 builds up superclusters and breaks them into clusters. The input is
31 in TObjarray and cluster information is in TObjArray.
32 integer clno gives total number of clusters in the supermodule.
33 fClusters is the global ( public ) variables.
34 Others are local ( private ) to the code.
35 At the moment, the data is read for whole detector ( all supermodules
36 and pmd as well as cpv. This will have to be modify later )
37 LAST UPDATE : October 23, 2002
38 -----------------------------------------------------------------------*/
40 #include <Riostream.h>
42 #include <TObjArray.h>
45 #include "AliPMDcludata.h"
46 #include "AliPMDcluster.h"
47 #include "AliPMDClustering.h"
48 #include "AliPMDClusteringV2.h"
51 ClassImp(AliPMDClusteringV2)
53 const Double_t AliPMDClusteringV2::fgkSqroot3by2=0.8660254; // sqrt(3.)/2.
55 AliPMDClusteringV2::AliPMDClusteringV2():
56 fPMDclucont(new TObjArray()),
59 for(int i = 0; i < kNDIMX; i++)
61 for(int j = 0; j < kNDIMY; j++)
63 fCoord[0][i][j] = i+j/2.;
64 fCoord[1][i][j] = fgkSqroot3by2*j;
68 // ------------------------------------------------------------------------ //
71 AliPMDClusteringV2::AliPMDClusteringV2(const AliPMDClusteringV2& pmdclv2):
72 AliPMDClustering(pmdclv2),
77 AliError("Copy constructor not allowed ");
80 // ------------------------------------------------------------------------ //
81 AliPMDClusteringV2 &AliPMDClusteringV2::operator=(const AliPMDClusteringV2& /*pmdclv2*/)
84 AliError("Assignment operator not allowed ");
87 // ------------------------------------------------------------------------ //
88 AliPMDClusteringV2::~AliPMDClusteringV2()
92 // ------------------------------------------------------------------------ //
94 void AliPMDClusteringV2::DoClust(Int_t idet, Int_t ismn,
95 Double_t celladc[48][96], TObjArray *pmdcont)
97 // main function to call other necessary functions to do clustering
99 AliPMDcluster *pmdcl = 0;
101 const Float_t ktwobysqrt3 = 1.1547; // 2./sqrt(3.)
102 const Int_t kNmaxCell = 19; // # of cells surrounding a cluster center
103 Int_t i, j, nmx1, incr, id, jd;
106 Int_t celldataX[kNmaxCell], celldataY[kNmaxCell];
108 Double_t cutoff, ave;
109 Double_t edepcell[kNMX];
117 else if (ismn >= 12 && ismn <= 23)
123 for (i =0; i < kNMX; i++)
128 for (id = 0; id < ndimXr; id++)
130 for (jd = 0; jd < ndimYr; jd++)
133 i = id + (ndimYr/2-1) - (jd/2);
134 Int_t ij = i + j*kNDIMX;
137 edepcell[ij] = celladc[jd][id];
139 else if (ismn >= 12 && ismn <= 23)
141 edepcell[ij] = celladc[id][jd];
148 TMath::Sort((Int_t)kNMX,edepcell,iord1);// order the data
149 cutoff = fCutoff; // cutoff used to discard cells having ener. dep.
153 for(i = 0;i < kNMX; i++)
159 if(edepcell[i] > cutoff )
165 AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1));
173 AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f",
176 incr = CrClust(ave, cutoff, nmx1,iord1, edepcell);
177 RefClust(incr,edepcell );
179 Int_t nentries1 = fPMDclucont->GetEntries();
180 AliDebug(1,Form("Detector Plane = %d Serial Module No = %d Number of clusters = %d",idet, ismn, nentries1));
181 AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
182 for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
184 AliPMDcludata *pmdcludata =
185 (AliPMDcludata*)fPMDclucont->UncheckedAt(ient1);
186 Float_t cluXC = pmdcludata->GetClusX();
187 Float_t cluYC = pmdcludata->GetClusY();
188 Float_t cluADC = pmdcludata->GetClusADC();
189 Float_t cluCELLS = pmdcludata->GetClusCells();
190 Float_t cluSIGX = pmdcludata->GetClusSigmaX();
191 Float_t cluSIGY = pmdcludata->GetClusSigmaY();
193 Float_t cluY0 = ktwobysqrt3*cluYC;
194 Float_t cluX0 = cluXC - cluY0/2.;
197 // Cluster X centroid is back transformed
201 clusdata[0] = cluX0 - (24-1) + cluY0/2.;
203 else if (ismn >= 12 && ismn <= 23)
205 clusdata[0] = cluX0 - (48-1) + cluY0/2.;
209 clusdata[2] = cluADC;
210 clusdata[3] = cluCELLS;
211 clusdata[4] = cluSIGX;
212 clusdata[5] = cluSIGY;
214 // Cells associated with a cluster
216 for (Int_t ihit = 0; ihit < kNmaxCell; ihit++)
218 Int_t dummyXY = pmdcludata->GetCellXY(ihit);
220 Int_t celldumY = dummyXY%10000;
221 Int_t celldumX = dummyXY/10000;
222 Float_t cellY = (Float_t) celldumY/10;
223 Float_t cellX = (Float_t) celldumX/10;
226 // Cell X centroid is back transformed
230 celldataX[ihit] = (Int_t) ((cellX - (24-1) + cellY/2.) + 0.5);
232 else if (ismn >= 12 && ismn <= 23)
234 celldataX[ihit] = (Int_t) ((cellX - (48-1) + cellY/2.) + 0.5 );
236 celldataY[ihit] = (Int_t) (cellY + 0.5);
239 pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY);
242 fPMDclucont->Delete();
244 // ------------------------------------------------------------------------ //
245 Int_t AliPMDClusteringV2::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
246 Int_t iord1[], Double_t edepcell[])
248 // Does crude clustering
249 // Finds out only the big patch by just searching the
253 Int_t i,j,k,id1,id2,icl, numcell;
254 Int_t jd1,jd2, icell, cellcount;
255 Int_t clust[2][5000];
256 static Int_t neibx[6] = {1,0,-1,-1,0,1}, neiby[6] = {0,1,1,0,-1,-1};
258 // neibx and neiby define ( incremental ) (i,j) for the neighbours of a
259 // cell. There are six neighbours.
260 // cellcount --- total number of cells having nonzero ener dep
261 // numcell --- number of cells in a given supercluster
263 AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
265 for (j=0; j < kNDIMX; j++)
267 for(k=0; k < kNDIMY; k++)
269 fInfocl[0][j][k] = 0;
270 fInfocl[1][j][k] = 0;
274 for(i=0; i < kNMX; i++)
282 if(edepcell[j] <= cutoff)
284 fInfocl[0][id1][id2] = -1;
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 // ---------------------------------------------------------------
293 for(icell=0; icell <= nmx1; icell++)
298 if(fInfocl[0][id1][id2] == 0 )
300 // ---------------------------------------------------------------
301 // icl -- cluster #, numcell -- # of cells in it, clust -- stores
302 // coordinates of the cells in a cluster, fInfocl[0][i1][i2] is 1 for
303 // primary and 2 for secondary cells,
304 // fInfocl[1][i1][i2] stores cluster #
305 // ---------------------------------------------------------------
309 fInfocl[0][id1][id2] = 1;
310 fInfocl[1][id1][id2] = icl;
311 fInfcl[0][cellcount] = icl;
312 fInfcl[1][cellcount] = id1;
313 fInfcl[2][cellcount] = id2;
315 clust[0][numcell] = id1;
316 clust[1][numcell] = id2;
317 for(i = 1; i < 5000; i++)
321 // ---------------------------------------------------------------
322 // check for adc count in neib. cells. If ne 0 put it in this clust
323 // ---------------------------------------------------------------
324 for(i = 0; i < 6; i++)
326 jd1 = id1 + neibx[i];
327 jd2 = id2 + neiby[i];
328 if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) &&
329 fInfocl[0][jd1][jd2] == 0)
332 fInfocl[0][jd1][jd2] = 2;
333 fInfocl[1][jd1][jd2] = icl;
334 clust[0][numcell] = jd1;
335 clust[1][numcell] = jd2;
337 fInfcl[0][cellcount] = icl;
338 fInfcl[1][cellcount] = jd1;
339 fInfcl[2][cellcount] = jd2;
342 // ---------------------------------------------------------------
343 // check adc count for neighbour's neighbours recursively and
344 // if nonzero, add these to the cluster.
345 // ---------------------------------------------------------------
346 for(i = 1;i < 5000; i++)
348 if(clust[0][i] != -1)
352 for(j = 0; j < 6 ; j++)
354 jd1 = id1 + neibx[j];
355 jd2 = id2 + neiby[j];
356 if( (jd1 >= 0 && jd1 < kNDIMX) &&
357 (jd2 >= 0 && jd2 < kNDIMY)
358 && fInfocl[0][jd1][jd2] == 0 )
360 fInfocl[0][jd1][jd2] = 2;
361 fInfocl[1][jd1][jd2] = icl;
363 clust[0][numcell] = jd1;
364 clust[1][numcell] = jd2;
366 fInfcl[0][cellcount] = icl;
367 fInfcl[1][cellcount] = jd1;
368 fInfcl[2][cellcount] = jd2;
377 // ------------------------------------------------------------------------ //
378 void AliPMDClusteringV2::RefClust(Int_t incr, Double_t edepcell[])
380 // Does the refining of clusters
381 // Takes the big patch and does gaussian fitting and
382 // finds out the more refined clusters
384 const Float_t ktwobysqrt3 = 1.1547;
385 const Int_t kNmaxCell = 19;
387 AliPMDcludata *pmdcludata = 0;
390 Int_t i, j, k, i1, i2, id, icl, itest, ihld;
391 Int_t ig, nsupcl, clno, clX,clY;
392 Int_t clxy[kNmaxCell];
395 Double_t x1, y1, z1, x2, y2, z2, rr;
397 Int_t kndim = incr + 1;
404 Double_t *x, *y, *z, *xc, *yc, *zc, *cells, *rcl, *rcs;
406 ncl = new Int_t [kndim];
407 iord = new Int_t [kndim];
408 x = new Double_t [kndim];
409 y = new Double_t [kndim];
410 z = new Double_t [kndim];
411 xc = new Double_t [kndim];
412 yc = new Double_t [kndim];
413 zc = new Double_t [kndim];
414 cells = new Double_t [kndim];
415 rcl = new Double_t [kndim];
416 rcs = new Double_t [kndim];
418 for(Int_t kk = 0; kk < 15; kk++)
420 if( kk < 6 )clusdata[kk] = 0.;
423 // nsupcl = # of superclusters; ncl[i]= # of cells in supercluster i
424 // x, y and z store (x,y) coordinates of and energy deposited in a cell
425 // xc, yc store (x,y) coordinates of the cluster center
426 // zc stores the energy deposited in a cluster, rc is cluster radius
431 for(i = 0; i < kndim; i++)
435 for(i = 0; i <= incr; i++)
437 if(fInfcl[0][i] != nsupcl)
443 AliWarning("RefClust: Too many superclusters!");
450 AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d",
455 for(i = 0; i <= nsupcl; i++)
461 // one cell super-clusters --> single cluster
462 // cluster center at the centyer of the cell
463 // cluster radius = half cell dimension
466 AliWarning("RefClust: Too many clusters! more than 5000");
472 i12 = i1 + i2*kNDIMX;
473 clusdata[0] = fCoord[0][i1][i2];
474 clusdata[1] = fCoord[1][i1][i2];
475 clusdata[2] = edepcell[i12];
482 clY = (Int_t)((ktwobysqrt3*fCoord[1][i1][i2])*10);
483 clX = (Int_t)((fCoord[0][i1][i2] - clY/20.)*10);
484 clxy[0] = clX*10000 + clY ;
486 for(Int_t icltr = 1; icltr < kNmaxCell; icltr++)
490 pmdcludata = new AliPMDcludata(clusdata,clxy);
491 fPMDclucont->Add(pmdcludata);
497 // two cell super-cluster --> single cluster
498 // cluster center is at ener. dep.-weighted mean of two cells
499 // cluster radius == half cell dimension
504 AliWarning("RefClust: Too many clusters! more than 5000");
510 i12 = i1 + i2*kNDIMX;
512 x1 = fCoord[0][i1][i2];
513 y1 = fCoord[1][i1][i2];
519 i12 = i1 + i2*kNDIMX;
521 x2 = fCoord[0][i1][i2];
522 y2 = fCoord[1][i1][i2];
525 clusdata[0] = (x1*z1+x2*z2)/(z1+z2);
526 clusdata[1] = (y1*z1+y2*z2)/(z1+z2);
529 clusdata[4] = (TMath::Sqrt(z1*z2))/(z1+z2);
532 clY = (Int_t)((ktwobysqrt3*y1)*10);
533 clX = (Int_t)((x1 - clY/20.)*10);
534 clxy[0] = clX*10000 + clY ;
536 clY = (Int_t)((ktwobysqrt3*y2)*10);
537 clX = (Int_t)((x2 - clY/20.)*10);
538 clxy[1] = clX*10000 + clY ;
540 for(Int_t icltr = 2; icltr < kNmaxCell; icltr++)
544 pmdcludata = new AliPMDcludata(clusdata, clxy);
545 fPMDclucont->Add(pmdcludata);
550 // super-cluster of more than two cells - broken up into smaller
551 // clusters gaussian centers computed. (peaks separated by > 1 cell)
552 // Begin from cell having largest energy deposited This is first
554 // *****************************************************************
555 // NOTE --- POSSIBLE MODIFICATION: ONE MAY NOT BREAKING SUPERCLUSTERS
556 // IF NO. OF CELLS IS NOT TOO LARGE ( SAY 5 OR 6 )
557 // SINCE WE EXPECT THE SUPERCLUSTER
558 // TO BE A SINGLE CLUSTER
559 //*******************************************************************
563 i12 = i1 + i2*kNDIMX;
565 x[0] = fCoord[0][i1][i2];
566 y[0] = fCoord[1][i1][i2];
567 z[0] = edepcell[i12];
570 for(j = 1; j <= ncl[i]; j++)
576 i12 = i1 + i2*kNDIMX;
578 x[j] = fCoord[0][i1][i2];
579 y[j] = fCoord[1][i1][i2];
580 z[j] = edepcell[i12];
583 // arranging cells within supercluster in decreasing order
584 for(j = 1; j <= ncl[i];j++)
588 for(i1 = 0; i1 < j; i1++)
590 if(itest == 0 && z[iord[i1]] < z[ihld])
593 for(i2 = j-1;i2 >= i1;i2--)
595 iord[i2+1] = iord[i2];
603 // compute the number of clusters and their centers ( first
605 // centers must be separated by cells having smaller ener. dep.
606 // neighbouring centers should be either strong or well-separated
611 for(j = 1; j <= ncl[i]; j++)
616 for(k = 0; k <= ig; k++)
620 rr = Distance(x1,y1,x2,y2);
621 //************************************************************
622 // finetuning cluster splitting
623 // the numbers zc/4 and zc/10 may need to be changed.
624 // Also one may need to add one more layer because our
625 // cells are smaller in absolute scale
626 //************************************************************
629 if( rr >= 1.1 && rr < 1.8 && z[iord[j]] > zc[k]/4.) itest++;
630 if( rr >= 1.8 && rr < 2.1 && z[iord[j]] > zc[k]/10.) itest++;
631 if( rr >= 2.1)itest++;
642 ClustDetails(ncl[i], ig, x, y ,z, xc, yc, zc, rcl, rcs, cells,
646 for(j = 0; j <= ig; j++)
651 AliWarning("RefClust: Too many clusters! more than 5000");
657 clusdata[4] = rcl[j];
658 clusdata[5] = rcs[j];
661 clusdata[3] = ncl[i] + 1;
665 clusdata[3] = cells[j];
668 Int_t ncellcls = testncl[j];
669 if( ncellcls < kNmaxCell )
671 for(Int_t kk = 1; kk <= ncellcls; kk++)
673 Int_t ll = testindex[pp];
674 clY = (Int_t)((ktwobysqrt3*y[ll])*10);
675 clX = (Int_t)((x[ll] - clY/20.)*10);
676 clxy[kk-1] = clX*10000 + clY ;
680 for(Int_t icltr = ncellcls ; icltr < kNmaxCell; icltr++)
685 pmdcludata = new AliPMDcludata(clusdata, clxy);
686 fPMDclucont->Add(pmdcludata);
704 // ------------------------------------------------------------------------ //
705 void AliPMDClusteringV2::ClustDetails(Int_t ncell, Int_t nclust, Double_t x[],
706 Double_t y[], Double_t z[],Double_t xc[],
707 Double_t yc[], Double_t zc[],
708 Double_t rcl[], Double_t rcs[],
709 Double_t cells[], TArrayI &testncl,
715 Int_t kndim1 = ncell + 1;//ncell
717 Int_t kndim3 = nclust + 1;//nclust
719 Int_t i, j, k, i1, i2;
720 Double_t x1, y1, x2, y2, rr, b, c, r1, r2;
721 Double_t sumx, sumy, sumxy, sumxx, sum, sum1, sumyy;
723 Double_t *str, *str1, *xcl, *ycl, *cln;
726 Double_t **clustcell;
727 str = new Double_t [kndim3];
728 str1 = new Double_t [kndim3];
729 xcl = new Double_t [kndim3];
730 ycl = new Double_t [kndim3];
731 cln = new Double_t [kndim3];
733 clustcell = new Double_t *[kndim3];
734 cell = new Int_t *[kndim3];
735 cluster = new Int_t *[kndim1];
736 for(i = 0; i < kndim1; i++)
738 cluster[i] = new Int_t [kndim2];
741 for(i = 0; i < kndim3; i++)
749 cell[i] = new Int_t [kndim2];
750 clustcell[i] = new Double_t [kndim1];
751 for(j = 0; j < kndim1; j++)
755 for(j = 0; j < kndim2; j++)
764 // more than one cluster
765 // checking cells shared between several clusters.
766 // First check if the cell is within
767 // one cell unit ( nearest neighbour). Else,
768 // if it is within 1.74 cell units ( next nearest )
769 // Else if it is upto 2 cell units etc.
771 for (i = 0; i <= ncell; i++)
777 // distance <= 1 cell unit
779 for(j = 0; j <= nclust; j++)
783 rr = Distance(x1, y1, x2, y2);
791 // next nearest neighbour
792 if(cluster[i][0] == 0)
794 for(j=0; j<=nclust; j++)
798 rr = Distance(x1, y1, x2, y2);
799 if(rr <= TMath::Sqrt(3.))
807 // next-to-next nearest neighbour
808 if(cluster[i][0] == 0)
810 for(j=0; j<=nclust; j++)
814 rr = Distance(x1, y1, x2, y2);
824 if(cluster[i][0] == 0)
826 for(j = 0; j <= nclust; j++)
830 rr = Distance(x1, y1, x2, y2);
841 // computing cluster strength. Some cells are shared.
842 for(i = 0; i <= ncell; i++)
844 if(cluster[i][0] != 0)
847 for(j = 1; j <= i1; j++)
855 for(k = 0; k < 5; k++)
857 for(i = 0; i <= ncell; i++)
859 if(cluster[i][0] != 0)
863 for(j = 1; j <= i1; j++)
865 sum += str[cluster[i][j]];
868 for(j = 1; j <= i1; j++)
871 str1[i2] += z[i]*str[i2]/sum;
872 clustcell[i2][i] = z[i]*str[i2]/sum;
878 for(j = 0; j <= nclust; j++)
885 for(i = 0; i <= nclust; i++)
891 for(j = 0; j <= ncell; j++)
893 if(clustcell[i][j] != 0)
895 sumx += clustcell[i][j]*x[j];
896 sumy += clustcell[i][j]*y[j];
897 sum += clustcell[i][j];
898 sum1 += clustcell[i][j]/z[j];
901 //** xcl and ycl are cluster centroid positions ( center of gravity )
909 for(j = 0; j <= ncell; j++)
911 sumxx += clustcell[i][j]*(x[j]-xcl[i])*(x[j]-xcl[i])/sum;
912 sumyy += clustcell[i][j]*(y[j]-ycl[i])*(y[j]-ycl[i])/sum;
913 sumxy += clustcell[i][j]*(x[j]-xcl[i])*(y[j]-ycl[i])/sum;
916 c = sumxx*sumyy-sumxy*sumxy;
917 // ******************r1 and r2 are major and minor axes ( r1 > r2 ).
918 r1 = b/2.+TMath::Sqrt(b*b/4.-c);
919 r2 = b/2.-TMath::Sqrt(b*b/4.-c);
920 // final assignments to proper external variables
930 //To get the cell position in a cluster
932 for(Int_t ii=0; ii<= ncell; ii++)
934 Int_t jj = cluster[ii][0];
935 for(Int_t kk=1; kk<= jj; kk++)
937 Int_t ll = cluster[ii][kk];
939 cell[ll][cell[ll][0]] = ii;
943 testncl.Set(nclust+1);
946 for(Int_t ii=0; ii <= nclust; ii++)
948 testncl[ii] = cell[ii][0];
949 counter += testncl[ii];
951 testindex.Set(counter);
953 for(Int_t ii=0; ii<= nclust; ii++)
955 for(Int_t jj = 1; jj<= testncl[ii]; jj++)
957 Int_t kk = cell[ii][jj];
971 for(j = 0; j <= ncell; j++)
984 for(j = 0; j <= ncell; j++)
986 sumxx += clustcell[i][j]*(x[j]-xcl[i])*(x[j]-xcl[i])/sum;
987 sumyy += clustcell[i][j]*(y[j]-ycl[i])*(y[j]-ycl[i])/sum;
988 sumxy += clustcell[i][j]*(x[j]-xcl[i])*(y[j]-ycl[i])/sum;
991 c = sumxx*sumyy-sumxy*sumxy;
992 r1 = b/2.+ TMath::Sqrt(b*b/4.-c);
993 r2 = b/2.- TMath::Sqrt(b*b/4.-c);
995 // To get the cell position in a cluster
996 testncl.Set(nclust+1);
997 testindex.Set(ncell+1);
998 cell[0][0] = ncell + 1;
999 testncl[0] = cell[0][0];
1001 for(Int_t ii = 1; ii <= ncell; ii++)
1004 Int_t kk = cell[0][ii];
1008 // final assignments
1016 for(i = 0; i < kndim3; i++)
1018 delete [] clustcell[i];
1021 delete [] clustcell;
1023 for(i = 0; i <kndim1 ; i++)
1025 delete [] cluster[i];
1035 // ------------------------------------------------------------------------ //
1036 Double_t AliPMDClusteringV2::Distance(Double_t x1, Double_t y1,
1037 Double_t x2, Double_t y2)
1039 return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
1041 // ------------------------------------------------------------------------ //
1042 void AliPMDClusteringV2::SetEdepCut(Float_t decut)
1046 // ------------------------------------------------------------------------ //