From fd30f88ef1c7f7c03cb419ebaac8a9cbc01bfbc5 Mon Sep 17 00:00:00 2001 From: bnandi Date: Sat, 10 Feb 2007 07:24:52 +0000 Subject: [PATCH] method Ranmar is removed --- PMD/AliPMDClusteringV1.cxx | 1132 ++++++++++++++++++------------------ PMD/AliPMDClusteringV1.h | 41 +- 2 files changed, 588 insertions(+), 585 deletions(-) diff --git a/PMD/AliPMDClusteringV1.cxx b/PMD/AliPMDClusteringV1.cxx index b80f78fd7ec..c20acd257de 100644 --- a/PMD/AliPMDClusteringV1.cxx +++ b/PMD/AliPMDClusteringV1.cxx @@ -30,11 +30,11 @@ Bhubaneswar 751 005 ( phatak@iopb.res.in ) Given the energy deposited ( or ADC value ) in each cell of supermodule ( pmd or cpv ), the code builds up superclusters and breaks them into clusters. The input is - in array fEdepCell[kNDIMX][kNDIMY] and cluster information is in array - fClusters[5][5000]. integer fClno gives total number of clusters in the + in array fEdepCell[kNDIMX][kNDIMY] and cluster information is in a + TObjarray. Integer clno gives total number of clusters in the supermodule. - fEdepCell, fClno and fClusters are the only global ( public ) variables. + fEdepCell and fClusters are the only global ( public ) variables. Others are local ( private ) to the code. At the moment, the data is read for whole detector ( all supermodules and pmd as well as cpv. This will have to be modify later ) @@ -47,22 +47,24 @@ #include #include +#include "AliPMDcludata.h" #include "AliPMDcluster.h" #include "AliPMDClustering.h" #include "AliPMDClusteringV1.h" #include "AliLog.h" +#include "TRandom.h" ClassImp(AliPMDClusteringV1) const Double_t AliPMDClusteringV1::fgkSqroot3by2=0.8660254; // sqrt(3.)/2. AliPMDClusteringV1::AliPMDClusteringV1(): - fClno(0), + pmdclucont(new TObjArray()), fCutoff(0.0) { - for(int i = 0; i < kNDIMX; i++) + for(Int_t i = 0; i < kNDIMX; i++) { - for(int j = 0; j < kNDIMY; j++) + for(Int_t j = 0; j < kNDIMY; j++) { fCoord[0][i][j] = i+j/2.; fCoord[1][i][j] = fgkSqroot3by2*j; @@ -73,10 +75,11 @@ AliPMDClusteringV1::AliPMDClusteringV1(): // ------------------------------------------------------------------------ // AliPMDClusteringV1::~AliPMDClusteringV1() { - + delete pmdclucont; } // ------------------------------------------------------------------------ // -void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn, Double_t celladc[48][96], TObjArray *pmdcont) +void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn, + Double_t celladc[48][96], TObjArray *pmdcont) { // main function to call other necessary functions to do clustering // @@ -87,7 +90,7 @@ void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn, Double_t celladc[48][96 and 0 <= jd <=96 */ - int i, i1, i2, j, nmx1, incr, id, jd; + Int_t i, i1, i2, j, nmx1, incr, id, jd; Int_t celldataX[15], celldataY[15]; Float_t clusdata[6]; @@ -97,8 +100,8 @@ void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn, Double_t celladc[48][96 // ndimXr and ndimYr are different because of different module size - Int_t ndimXr =0; - Int_t ndimYr =0; + Int_t ndimXr = 0; + Int_t ndimYr = 0; if (ismn < 12) { @@ -111,9 +114,9 @@ void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn, Double_t celladc[48][96 ndimYr = 96; } - for (Int_t i =0; i < kNDIMX; i++) + for (Int_t i = 0; i < kNDIMX; i++) { - for (Int_t j =0; j < kNDIMY; j++) + for (Int_t j = 0; j < kNDIMY; j++) { fEdepCell[i][j] = 0; fCellTrNo[i][j] = -1; @@ -124,8 +127,8 @@ void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn, Double_t celladc[48][96 { for (jd = 0; jd < ndimYr; jd++) { - j=jd; - i=id+(ndimYr/2-1)-(jd/2); + j = jd; + i = id+(ndimYr/2-1)-(jd/2); if (ismn < 12) { @@ -140,44 +143,52 @@ void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn, Double_t celladc[48][96 } } - Order(); // order the data + Order(); // order the data cutoff = fCutoff; // cutoff used to discard cells having ener. dep. - ave=0.; - nmx1=-1; + ave = 0.; + nmx1 = -1; - for(j=0;j 0.) {ave = ave + fEdepCell[i1][i2];} - if (fEdepCell[i1][i2] > cutoff ) nmx1 = nmx1 + 1; + if(fEdepCell[i1][i2] > 0.) + { + ave += fEdepCell[i1][i2]; + } + if(fEdepCell[i1][i2] > cutoff ) + { + nmx1++; + } } AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1)); if (nmx1 == 0) nmx1 = 1; - ave=ave/nmx1; + ave = ave/nmx1; AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f", kNMX,ave)); incr = CrClust(ave, cutoff, nmx1); RefClust(incr); - - AliDebug(1,Form("Detector Plane = %d Serial Module No = %d Number of clusters = %d",idet, ismn, fClno)); - - - for(i1=0; i1<=fClno; i1++) + Int_t nentries1 = pmdclucont->GetEntries(); + AliDebug(1,Form("Detector Plane = %d Serial Module No = %d Number of clusters = %d",idet, ismn, nentries1)); + AliDebug(1,Form("Total number of clusters/module = %d",nentries1)); + for (Int_t ient1 = 0; ient1 < nentries1; ient1++) { - Float_t cluXC = (Float_t) fClusters[0][i1]; - Float_t cluYC = (Float_t) fClusters[1][i1]; - Float_t cluADC = (Float_t) fClusters[2][i1]; - Float_t cluCELLS = (Float_t) fClusters[3][i1]; - Float_t cluRAD = (Float_t) fClusters[4][i1]; + AliPMDcludata *pmdcludata = + (AliPMDcludata*)pmdclucont->UncheckedAt(ient1); + Float_t cluXC = pmdcludata->GetClusX(); + Float_t cluYC = pmdcludata->GetClusY(); + Float_t cluADC = pmdcludata->GetClusADC(); + Float_t cluCELLS = pmdcludata->GetClusCells(); + Float_t cluSIGX = pmdcludata->GetClusSigmaX(); + Float_t cluSIGY = pmdcludata->GetClusSigmaY(); + Float_t cluY0 = ktwobysqrt3*cluYC; Float_t cluX0 = cluXC - cluY0/2.; - - + // // Cluster X centroid is back transformed // @@ -185,23 +196,23 @@ void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn, Double_t celladc[48][96 { clusdata[0] = cluX0 - (24-1) + cluY0/2.; } - else if (ismn >= 12 && ismn <= 23) + else if ( ismn >= 12 && ismn <= 23) { clusdata[0] = cluX0 - (48-1) + cluY0/2.; } - + clusdata[1] = cluY0; clusdata[2] = cluADC; clusdata[3] = cluCELLS; - clusdata[4] = cluRAD; - clusdata[5] = 0.; - + clusdata[4] = cluSIGX; + clusdata[5] = cluSIGY; + // // Cells associated with a cluster // for (Int_t ihit = 0; ihit < 15; ihit++) { - + if (ismn < 12) { celldataX[ihit] = fClTr[ihit][i1]%10000; @@ -213,10 +224,12 @@ void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn, Double_t celladc[48][96 celldataY[ihit] = fClTr[ihit][i1]%10000; } } - //printf("%d %f %f\n",idet,cluXC,cluYC ); pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY); pmdcont->Add(pmdcl); } + + pmdclucont->Clear(); + } // ------------------------------------------------------------------------ // void AliPMDClusteringV1::Order() @@ -224,14 +237,10 @@ void AliPMDClusteringV1::Order() // Sorting algorithm // sorts the ADC values from higher to lower // - double dd[kNMX]; - // matrix fEdepCell converted into - // one dimensional array dd. adum a place holder for double - int i, j, i1, i2, iord1[kNMX]; - // information of - // ordering is stored in iord1, original array not ordered - // - // define arrays dd and iord1 + Int_t i, j, i1, i2; + Int_t iord1[kNMX]; + Double_t dd[kNMX]; + for(i1=0; i1 < kNDIMX; i1++) { for(i2=0; i2 < kNDIMY; i2++) @@ -243,7 +252,6 @@ void AliPMDClusteringV1::Order() } TMath::Sort(kNMX,dd,iord1); //PH Using much better algorithm... - // store the sorted information in fIord for later use for(i=0; i= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) && - fInfocl[0][jd1][jd2] == 0){ - numcell=numcell+1; - fInfocl[0][jd1][jd2]=2; - fInfocl[1][jd1][jd2]=icl; - clust[0][numcell]=jd1; - clust[1][numcell]=jd2; - cellcount=cellcount+1; - fInfcl[0][cellcount]=icl; - fInfcl[1][cellcount]=jd1; - fInfcl[2][cellcount]=jd2; - } - } - // --------------------------------------------------------------- - // check adc count for neighbour's neighbours recursively and - // if nonzero, add these to the cluster. - // --------------------------------------------------------------- - for(i=1;i < 5000;i++){ - if(clust[0][i] != 0){ - id1=clust[0][i]; - id2=clust[1][i]; - for(j=0; j<6 ; j++){ - jd1=id1+neibx[j]; - jd2=id2+neiby[j]; - if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) && - fInfocl[0][jd1][jd2] == 0 ){ - fInfocl[0][jd1][jd2] = 2; - fInfocl[1][jd1][jd2] = icl; - numcell = numcell + 1; - clust[0][numcell] = jd1; - clust[1][numcell] = jd2; - cellcount = cellcount+1; - fInfcl[0][cellcount] = icl; - fInfcl[1][cellcount] = jd1; - fInfcl[2][cellcount] = jd2; + for(icell = 0; icell <= nmx1; icell++) + { + id1 = fIord[0][icell]; + id2 = fIord[1][icell]; + if(fInfocl[0][id1][id2] == 0 ) + { + icl++; + numcell = 0; + cellcount++; + fInfocl[0][id1][id2] = 1; + fInfocl[1][id1][id2] = icl; + fInfcl[0][cellcount] = icl; + fInfcl[1][cellcount] = id1; + fInfcl[2][cellcount] = id2; + + clust[0][numcell] = id1; + clust[1][numcell] = id2; + + for(i = 1; i < 5000; i++) + { + clust[0][i]=0; + } + // --------------------------------------------------------------- + // check for adc count in neib. cells. If ne 0 put it in this clust + // --------------------------------------------------------------- + for(i = 0; i < 6; i++) + { + jd1 = id1 + neibx[i]; + jd2 = id2 + neiby[i]; + if( (jd1 >= 0 && jd1 < kNDIMX) && (jd2 >= 0 && jd2 < kNDIMY) && + fInfocl[0][jd1][jd2] == 0) + { + numcell++; + fInfocl[0][jd1][jd2] = 2; + fInfocl[1][jd1][jd2] = icl; + clust[0][numcell] = jd1; + clust[1][numcell] = jd2; + cellcount++; + fInfcl[0][cellcount] = icl; + fInfcl[1][cellcount] = jd1; + fInfcl[2][cellcount] = jd2; + } + } + // --------------------------------------------------------------- + // check adc count for neighbour's neighbours recursively and + // if nonzero, add these to the cluster. + // --------------------------------------------------------------- + for(i = 1; i < 5000;i++) + { + if(clust[0][i] != 0) + { + id1 = clust[0][i]; + id2 = clust[1][i]; + for(j = 0; j < 6 ; j++) + { + jd1 = id1 + neibx[j]; + jd2 = id2 + neiby[j]; + if( (jd1 >= 0 && jd1 < kNDIMX) && + (jd2 >= 0 && jd2 < kNDIMY) && + fInfocl[0][jd1][jd2] == 0 ) + { + fInfocl[0][jd1][jd2] = 2; + fInfocl[1][jd1][jd2] = icl; + numcell++; + clust[0][numcell] = jd1; + clust[1][numcell] = jd2; + cellcount++; + fInfcl[0][cellcount] = icl; + fInfcl[1][cellcount] = jd1; + fInfcl[2][cellcount] = jd2; + } + } + } } - } } - } } - } - - // for(icell=0; icell<=cellcount; icell++){ - // ofl0 << fInfcl[0][icell] << " " << fInfcl[1][icell] << " " << - // fInfcl[2][icell] << endl; - //} - return cellcount; } // ------------------------------------------------------------------------ // -void AliPMDClusteringV1::RefClust(int incr) +void AliPMDClusteringV1::RefClust(Int_t incr) { // Does the refining of clusters // Takes the big patch and does gaussian fitting and // finds out the more refined clusters // - int i, j, k, i1, i2, id, icl, ncl[4500], iord[4500], itest; - int ihld; - int ig, nsupcl, lev1[20], lev2[20]; - double x[4500], y[4500], z[4500], x1, y1, z1, x2, y2, z2, dist; - double xc[4500], yc[4500], zc[4500], cells[4500], sum, rc[4500], rr; + const Int_t kdim = 4500; + Int_t i, j, k, i1, i2, id, icl, itest,ihld, ig, nsupcl,clno; + Double_t x1, y1, z1, x2, y2, z2, dist,rr,sum; - //asso - Int_t t[4500],cellCount[4500]; - for(i=0; i<4500; i++) + Int_t t[kdim],cellCount[kdim]; + Int_t ncl[kdim], iord[kdim], lev1[20], lev2[20]; + Double_t x[kdim], y[kdim], z[kdim]; + Double_t xc[kdim], yc[kdim], zc[kdim], cells[kdim], rc[kdim]; + + Float_t clusdata[6]; + for(Int_t kk = 0; kk < 6; kk++) { - t[i]=-1; - cellCount[i]=0; + clusdata[kk] = 0.; } - - // fClno counts the final clusters + //asso + for(i = 0; i 4500) { - AliWarning("RefClust: Too many superclusters!"); - nsupcl = 4500; - break; + for(i = 0; i < kdim; i++) + { + ncl[i] = -1; } - - ncl[nsupcl]=ncl[nsupcl]+1; - } - - AliDebug(1,Form("Number of cells = %d Number of Superclusters = %d", - incr+1,nsupcl+1)); - id=-1; - icl=-1; - for(i=0; i<=nsupcl; i++) { - if(ncl[i] == 0){ - id=id+1; - icl=icl+1; - // one cell super-clusters --> single cluster - // cluster center at the centyer of the cell - // cluster radius = half cell dimension - if (fClno >= 5000) { - AliWarning("RefClust: Too many clusters! more than 5000"); - return; - } - fClno = fClno + 1; - i1 = fInfcl[1][id]; - i2 = fInfcl[2][id]; - fClusters[0][fClno] = fCoord[0][i1][i2]; - fClusters[1][fClno] = fCoord[1][i1][i2]; - fClusters[2][fClno] = fEdepCell[i1][i2]; - fClusters[3][fClno] = 1.; - fClusters[4][fClno] = 0.5; - - //association - - fClTr[0][fClno]=fCellTrNo[i1][i2]; - for(Int_t icltr=1;icltr<14;icltr++) + for(i = 0; i <= incr; i++) + { + if(fInfcl[0][i] != nsupcl) { - fClTr[icltr][fClno]=-1; + nsupcl++; } - - //ofl1 << icl << " " << fCoord[0][i1][i2] << " " << fCoord[1][i1][i2] << - //" " << fEdepCell[i1][i2] << " " << fClusters[3][fClno] < single cluster - // cluster center is at ener. dep.-weighted mean of two cells - // cluster radius == half cell dimension - id = id + 1; - icl = icl+1; - if (fClno >= 5000) { - AliWarning("RefClust: Too many clusters! more than 5000"); - return; - } - fClno = fClno+1; - i1 = fInfcl[1][id]; - i2 = fInfcl[2][id]; - x1 = fCoord[0][i1][i2]; - y1 = fCoord[1][i1][i2]; - z1 = fEdepCell[i1][i2]; - - //asso - fClTr[0][fClno]=fCellTrNo[i1][i2]; - // - - id = id+1; - i1 = fInfcl[1][id]; - i2 = fInfcl[2][id]; - x2 = fCoord[0][i1][i2]; - y2 = fCoord[1][i1][i2]; - z2 = fEdepCell[i1][i2]; - - //asso - - fClTr[1][fClno]=fCellTrNo[i1][i2]; - for(Int_t icltr=2;icltr<14;icltr++) + if (nsupcl > kdim) { - fClTr[icltr][fClno] = -1; + AliWarning("RefClust: Too many superclusters!"); + nsupcl = kdim; + break; } - // - - fClusters[0][fClno] = (x1*z1+x2*z2)/(z1+z2); - fClusters[1][fClno] = (y1*z1+y2*z2)/(z1+z2); - fClusters[2][fClno] = z1+z2; - fClusters[3][fClno] = 2.; - fClusters[4][fClno] = 0.5; - - - //ofl1 << icl << " " << fClusters[0][fClno] << " " << fClusters[1][fClno] - // << " " << fClusters[2][fClno] << " " < 1 cell) - // Begin from cell having largest energy deposited This is first - // cluster center - i1 = fInfcl[1][id]; - i2 = fInfcl[2][id]; - x[0] = fCoord[0][i1][i2]; - y[0] = fCoord[1][i1][i2]; - z[0] = fEdepCell[i1][i2]; - - //asso - t[0]=fCellTrNo[i1][i2]; - // - - iord[0] = 0; - for(j=1;j<=ncl[i];j++){ - - id = id + 1; - i1 = fInfcl[1][id]; - i2 = fInfcl[2][id]; - iord[j] = j; - x[j] = fCoord[0][i1][i2]; - y[j] = fCoord[1][i1][i2]; - z[j] = fEdepCell[i1][i2]; - - //asso - t[j]=fCellTrNo[i1][i2]; - // - - - } - // arranging cells within supercluster in decreasing order - for(j=1;j<=ncl[i];j++){ - itest=0; - ihld=iord[j]; - for(i1=0;i1=i1;i2--){ - iord[i2+1]=iord[i2]; + id++; + icl++; + if (clno >= 5000) + { + AliWarning("RefClust: Too many clusters! more than 5000"); + return; + } + clno++; + i1 = fInfcl[1][id]; + i2 = fInfcl[2][id]; + + clusdata[0] = fCoord[0][i1][i2]; + clusdata[1] = fCoord[1][i1][i2]; + clusdata[2] = fEdepCell[i1][i2]; + clusdata[3] = 1.; + clusdata[4] = 0.5; + clusdata[5] = 0.0; + pmdcludata = new AliPMDcludata(clusdata); + pmdclucont->Add(pmdcludata); + + //association + + fClTr[0][clno] = fCellTrNo[i1][i2]; + for(Int_t icltr = 1; icltr < 14; icltr++) + { + fClTr[icltr][clno] = -1; } - iord[i1]=ihld; - } - } - } - - // compute the number of Gaussians and their centers ( first - // guess ) - // centers must be separated by cells having smaller ener. dep. - // neighbouring centers should be either strong or well-separated - ig=0; - xc[ig]=x[iord[0]]; - yc[ig]=y[iord[0]]; - zc[ig]=z[iord[0]]; - for(j=1;j<=ncl[i];j++){ - itest=-1; - x1=x[iord[j]]; - y1=y[iord[j]]; - for(k=0;k<=ig;k++){ - x2=xc[k]; y2=yc[k]; - rr=Distance(x1,y1,x2,y2); - if( rr >= 1.1 && rr < 1.8 && z[iord[j]] > zc[k]/4.) - itest=itest+1; - if( rr >= 1.8 && rr < 2.1 && z[iord[j]] > zc[k]/10.) - itest=itest+1; - if( rr >= 2.1)itest=itest+1; } - if(itest == ig){ - ig=ig+1; - xc[ig]=x1; - yc[ig]=y1; - zc[ig]=z[iord[j]]; + else if(ncl[i] == 1) + { + id++; + icl++; + if (clno >= 5000) + { + AliWarning("RefClust: Too many clusters! more than 5000"); + return; + } + clno++; + i1 = fInfcl[1][id]; + i2 = fInfcl[2][id]; + x1 = fCoord[0][i1][i2]; + y1 = fCoord[1][i1][i2]; + z1 = fEdepCell[i1][i2]; + + //asso + fClTr[0][clno] = fCellTrNo[i1][i2]; + // + + id = id+1; + i1 = fInfcl[1][id]; + i2 = fInfcl[2][id]; + x2 = fCoord[0][i1][i2]; + y2 = fCoord[1][i1][i2]; + z2 = fEdepCell[i1][i2]; + + //asso + + fClTr[1][clno] = fCellTrNo[i1][i2]; + for(Int_t icltr = 2; icltr < 14; icltr++) + { + fClTr[icltr][clno] = -1; + } + // + + clusdata[0] = (x1*z1+x2*z2)/(z1+z2); + clusdata[1] = (y1*z1+y2*z2)/(z1+z2); + clusdata[2] = z1+z2; + clusdata[3] = 2.; + clusdata[4] = 0.5; + clusdata[5] = 0.0; + pmdcludata = new AliPMDcludata(clusdata); + pmdclucont->Add(pmdcludata); } - } - // for(j=0; j<=ig; j++){ - //ofl1 << icl+j+1 << " " << xc[j] << " " < 0){ - for(j=0; j<=ncl[i]; j++){ - lev1[j]=0; - lev2[j]=0; - for(k=0; k<=ig; k++){ - dist=Distance(x[j], y[j], xc[k], yc[k]); - if(dist < sqrt(3.) ){ - + else + { + //asso + for(Int_t icg = 0; icg < kdim; icg++) + { + cellCount[icg]=0; + } + // + + id++; + iord[0] = 0; + // super-cluster of more than two cells - broken up into smaller + // clusters gaussian centers computed. (peaks separated by > 1 cell) + // Begin from cell having largest energy deposited This is first + // cluster center + i1 = fInfcl[1][id]; + i2 = fInfcl[2][id]; + x[0] = fCoord[0][i1][i2]; + y[0] = fCoord[1][i1][i2]; + z[0] = fEdepCell[i1][i2]; + + //asso + t[0] = fCellTrNo[i1][i2]; + // + + iord[0] = 0; + for(j = 1; j <= ncl[i]; j++) + { + id++; + i1 = fInfcl[1][id]; + i2 = fInfcl[2][id]; + iord[j] = j; + x[j] = fCoord[0][i1][i2]; + y[j] = fCoord[1][i1][i2]; + z[j] = fEdepCell[i1][i2]; //asso - fClTr[cellCount[k]][fClno+k+1]=t[j]; - cellCount[k]++; + t[j] = fCellTrNo[i1][i2]; // - - lev1[0]++; - i1=lev1[0]; - lev1[i1]=k; - }else{ - if(dist < 2.1){ - lev2[0]++; - i1=lev2[0]; - lev2[i1]=k; - } } - } - if(lev1[0] != 0){ - if(lev1[0] == 1){cells[lev1[1]]=cells[lev1[1]]+1.;} - else{ - sum=0.; - for(k=1; k<=lev1[0]; k++){ - sum=sum+zc[lev1[k]]; - } - for(k=1; k<=lev1[0]; k++){ - cells[lev1[k]]=cells[lev1[k]]+zc[lev1[k]]/sum; - } + + + // arranging cells within supercluster in decreasing order + + for(j = 1;j <= ncl[i]; j++) + { + itest = 0; + ihld = iord[j]; + for(i1 = 0; i1 < j; i1++) + { + if(itest == 0 && z[iord[i1]] < z[ihld]) + { + itest = 1; + for(i2 = j-1; i2 >= i1; i2--) + { + iord[i2+1] = iord[i2]; + } + iord[i1] = ihld; + } + } } - }else{ - if(lev2[0] == 0){cells[lev2[1]]=cells[lev2[1]]+1.;} - else{ - sum=0.; - for(k=1; k<=lev2[0]; k++){ - sum=sum+zc[lev2[k]]; - } - for(k=1; k<=lev2[0]; k++){ - cells[lev2[k]]=cells[lev2[k]]+zc[lev2[k]]/sum; - } + // compute the number of Gaussians and their centers ( first + // guess ) + // centers must be separated by cells having smaller ener. dep. + // neighbouring centers should be either strong or well-separated + ig=0; + xc[ig] = x[iord[0]]; + yc[ig] = y[iord[0]]; + zc[ig] = z[iord[0]]; + for(j = 1; j <= ncl[i]; j++) + { + itest = -1; + x1 = x[iord[j]]; + y1 = y[iord[j]]; + for(k = 0; k <= ig; k++) + { + x2 = xc[k]; + y2 = yc[k]; + rr = Distance(x1,y1,x2,y2); + if( rr >= 1.1 && rr < 1.8 && z[iord[j]] > zc[k]/4.) + { + itest++; + } + if( rr >= 1.8 && rr < 2.1 && z[iord[j]] > zc[k]/10.) + { + itest++; + } + if( rr >= 2.1) + { + itest++; + } + } + if(itest == ig) + { + ig++; + xc[ig] = x1; + yc[ig] = y1; + zc[ig] = z[iord[j]]; + } } - } - } - } - - // zero rest of the cell array - //asso - for(k=0; k<=ig; k++) - { - for(Int_t icltr=cellCount[k];icltr<14;icltr++) + + GaussFit(ncl[i], ig, x[0], y[0] ,z[0], xc[0], yc[0], zc[0], rc[0]); + icl += ig+1; + // compute the number of cells belonging to each cluster. + // cell is shared between several clusters ( if they are equidistant + // from it ) in the ratio of cluster energy deposition + for(j = 0; j <= ig; j++) { - fClTr[icltr][fClno]=-1; + cells[j]=0.; + } + if(ig > 0) + { + for(j = 0; j <= ncl[i]; j++) + { + lev1[j] = 0; + lev2[j] = 0; + for(k = 0; k <= ig; k++) + { + dist = Distance(x[j], y[j], xc[k], yc[k]); + if(dist < TMath::Sqrt(3.) ) + { + //asso + fClTr[cellCount[k]][clno+k+1] = t[j]; + cellCount[k]++; + // + lev1[0]++; + i1 = lev1[0]; + lev1[i1] = k; + } + else + { + if(dist < 2.1) + { + lev2[0]++; + i1 = lev2[0]; + lev2[i1] = k; + } + } + } + if(lev1[0] != 0) + { + if(lev1[0] == 1) + { + cells[lev1[1]]++; + } + else + { + sum=0.; + for(k = 1; k <= lev1[0]; k++) + { + sum += zc[lev1[k]]; + } + for(k = 1; k <= lev1[0]; k++) + { + cells[lev1[k]] += zc[lev1[k]]/sum; + } + } + } + else + { + if(lev2[0] == 0) + { + cells[lev2[1]]++; + } + else + { + sum=0.; + for( k = 1; k <= lev2[0]; k++) + { + sum += zc[lev2[k]]; + } + for(k = 1; k <= lev2[0]; k++) + { + cells[lev2[k]] += zc[lev2[k]]/sum; + } + } + } + } + } + + // zero rest of the cell array + //asso + for( k = 0; k <= ig; k++) + { + for(Int_t icltr = cellCount[k]; icltr < 14; icltr++) + { + fClTr[icltr][clno] = -1; + } + } + // + + for(j = 0; j <= ig; j++) + { + clno++; + if (clno >= 5000) + { + AliWarning("RefClust: Too many clusters! more than 5000"); + return; + } + clusdata[0] = xc[j]; + clusdata[1] = yc[j]; + clusdata[2] = zc[j]; + clusdata[4] = rc[j]; + clusdata[5] = 0.0; + if(ig == 0) + { + clusdata[3] = ncl[i]; + } + else + { + clusdata[3] = cells[j]; + } + pmdcludata = new AliPMDcludata(clusdata); + pmdclucont->Add(pmdcludata); } } - // - - - - for(j=0; j<=ig; j++){ - if (fClno >= 5000) { - AliWarning("RefClust: Too many clusters! more than 5000"); - return; - } - fClno = fClno + 1; - fClusters[0][fClno] = xc[j]; - fClusters[1][fClno] = yc[j]; - fClusters[2][fClno] = zc[j]; - fClusters[4][fClno] = rc[j]; - if(ig == 0){ - fClusters[3][fClno] = ncl[i]; - }else{ - fClusters[3][fClno] = cells[j]; - } - } } - } } // ------------------------------------------------------------------------ // -void AliPMDClusteringV1::GaussFit(Int_t ncell, Int_t nclust, Double_t &x, Double_t &y ,Double_t &z, Double_t &xc, Double_t &yc, Double_t &zc, Double_t &rc) +void AliPMDClusteringV1::GaussFit(Int_t ncell, Int_t nclust, Double_t &x, + Double_t &y ,Double_t &z, Double_t &xc, + Double_t &yc, Double_t &zc, Double_t &rc) { // Does gaussian fitting // - int i, j, i1, i2, novar, idd, jj; - double xx[4500], yy[4500], zz[4500], xxc[4500], yyc[4500]; - double a[4500], b[4500], c[4500], d[4500], ha[4500], hb[4500]; - double hc[4500], hd[4500], zzc[4500], rrc[4500]; - int neib[4500][50]; - double sum, dx, dy, str, str1, aint, sum1, rr, dum; - double x1, x2, y1, y2; + const Int_t kdim = 4500; + Int_t i, j, i1, i2, novar, idd, jj; + Double_t sum, dx, dy, str, str1, aint, sum1, rr, dum; + Double_t x1, x2, y1, y2; + Double_t xx[kdim], yy[kdim], zz[kdim], xxc[kdim], yyc[kdim]; + Double_t a[kdim], b[kdim], c[kdim], d[kdim], ha[kdim], hb[kdim]; + Double_t hc[kdim], hd[kdim], zzc[kdim], rrc[kdim]; + Int_t neib[kdim][50]; + + TRandom rnd; + str = 0.; str1 = 0.; rr = 0.3; novar = 0; - j = 0; // Just put not to see the compiler warning, BKN + j = 0; - for(i=0; i<=ncell; i++) + for(i = 0; i <= ncell; i++) { xx[i] = *(&x+i); yy[i] = *(&y+i); zz[i] = *(&z+i); - str = str + zz[i]; + str += zz[i]; } for(i=0; i<=nclust; i++) { xxc[i] = *(&xc+i); yyc[i] = *(&yc+i); zzc[i] = *(&zc+i); - str1 = str1 + zzc[i]; + str1 += zzc[i]; rrc[i] = 0.5; } - for(i=0; i<=nclust; i++) + for(i = 0; i <= nclust; i++) { zzc[i] = str/str1*zzc[i]; ha[i] = xxc[i]; @@ -716,148 +786,96 @@ void AliPMDClusteringV1::GaussFit(Int_t ncell, Int_t nclust, Double_t &x, Double x1 = xxc[i]; y1 = yyc[i]; } - for(i=0; i<=ncell; i++){ - idd=0; - x1=xx[i]; - y1=yy[i]; - for(j=0; j<=nclust; j++){ - x2=xxc[j]; - y2=yyc[j]; - if(Distance(x1,y1,x2,y2) <= 3.){ idd=idd+1; neib[i][idd]=j; } + for(i = 0; i <= ncell; i++) + { + idd = 0; + x1 = xx[i]; + y1 = yy[i]; + for(j = 0; j <= nclust; j++) + { + x2 = xxc[j]; + y2 = yyc[j]; + if(Distance(x1,y1,x2,y2) <= 3.) + { + idd++; + neib[i][idd] = j; + } + } + neib[i][0] = idd; } - neib[i][0]=idd; - } - sum=0.; - for(i1=0; i1<=ncell; i1++){ - aint=0.; - idd=neib[i1][0]; - for(i2=1; i2<=idd; i2++){ - jj=neib[i1][i2]; - dx=xx[i1]-xxc[jj]; - dy=yy[i1]-yyc[jj]; - dum=rrc[j]*rrc[jj]+rr*rr; - aint=aint+exp(-(dx*dx+dy*dy)/dum)*zzc[idd]*rr*rr/dum; + sum = 0.; + for(i1 = 0; i1 <= ncell; i1++) + { + aint = 0.; + idd = neib[i1][0]; + for(i2 = 1; i2 <= idd; i2++) + { + jj = neib[i1][i2]; + dx = xx[i1] - xxc[jj]; + dy = yy[i1] - yyc[jj]; + dum = rrc[j]*rrc[jj] + rr*rr; + aint += exp(-(dx*dx+dy*dy)/dum)*zzc[idd]*rr*rr/dum; + } + sum += (aint - zz[i1])*(aint - zz[i1])/str; + } + str1 = 0.; + for(i = 0; i <= nclust; i++) + { + a[i] = xxc[i] + 0.6*(rnd.Uniform() - 0.5); + b[i] = yyc[i] + 0.6*(rnd.Uniform() - 0.5); + c[i] = zzc[i]*(1.+ ( rnd.Uniform() - 0.5)*0.2); + str1 += zzc[i]; + d[i] = rrc[i]*(1.+ ( rnd.Uniform() - 0.5)*0.1); + + if(d[i] < 0.25) + { + d[i]=0.25; + } } - sum=sum+(aint-zz[i1])*(aint-zz[i1])/str; - } -// jmax=nclust*1000; -// if(nclust > 20)jmax=20000; -// for(j=0; j 31328 ) ii = ii - ( ii / 31328 ) * 31328; - if(jj > 30081 ) jj = jj - ( jj / 30081 ) * 30081; - itest=itest+1; - if((( ii > 0 ) && ( ii <= 31328 )) && (( jj > 0 ) && - ( jj <= 30081 ))){ - i1=ii/177+2; i2=ii-(i1-2)*177+2; i3=jj/169+1; i4=jj-(i3-1)*169; - i4 = jj - (i3-1)*169; - count1=0; - while ( count1 < 97 ){ - s=0.; - t=0.5; - count2=0; - while( count2 < 24 ){ - idum=i1*i2/179; - idum=( i1*i2 - (i1*i2/179)*179 ) * i3; - i5=idum-(idum/179)*179; - i1=i2; i2=i3; i3=i5; idum=53*i4+1; i4=idum-(idum/169)*169; - if( i4*i5-((i4*i5)/64)*64 >= 32 ) s=s+t; - t=0.5*t; - count2=count2+1; - } - u[count1] = s; - count1 = count1 +1; - } - c = 362436./16777216.; cd = 7654321./16777216.; - cm = 16777213./16777216.; - } - else{ - AliWarning("Wrong initialization"); - } - } - else{ - uni = u[i] - u[j]; - if( uni < 0.) uni = uni + 1; - u[i] = uni; - i = i -1; - if( i < 0 ) i = 96; - j = j - 1; - if ( j < 0 ) j = 96; - c = c - cd; - if( c < 0. ) c = c+cm; - uni = uni-c ; - if( uni < 0. )uni = uni+1.; - } - return uni; + return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); } // ------------------------------------------------------------------------ // void AliPMDClusteringV1::SetEdepCut(Float_t decut) diff --git a/PMD/AliPMDClusteringV1.h b/PMD/AliPMDClusteringV1.h index c20f1f4bf5a..d48ec65f578 100644 --- a/PMD/AliPMDClusteringV1.h +++ b/PMD/AliPMDClusteringV1.h @@ -31,10 +31,9 @@ class TNtuple; class TObjArray; class AliPMDcluster; - +class AliPMDcludata; class AliPMDClusteringV1: public AliPMDClustering { - public: AliPMDClusteringV1(); virtual ~AliPMDClusteringV1(); @@ -50,45 +49,31 @@ class AliPMDClusteringV1: public AliPMDClustering Double_t &yc, Double_t &zc, Double_t &rc); Double_t Distance(Double_t x1, Double_t y1, Double_t x2, Double_t y2); - Double_t Ranmar() const; void SetEdepCut(Float_t decut); protected: - + + TObjArray *pmdclucont; + AliPMDcludata *pmdcludata; + static const Double_t fgkSqroot3by2; // fgkSqroot3by2 = sqrt(3.)/2. enum { - kNMX = 11424, - kNDIMX = 119, - kNDIMY = 96 + kNMX = 11424, // no. of cells in a module + kNDIMX = 119, // max no. of cells along x direction + kNDIMY = 96 // max no. of cells along axis at 60 deg with x axis }; - - /* - kNMX : # of cells in a supermodule - kNDIMX : maximum number of cells along x direction (origin at one corner) - kNDIMY : maximum number of cells along axis at 60 degrees with x axis - */ Double_t fEdepCell[kNDIMX][kNDIMY]; //energy(ADC) in each cell - Double_t fClusters[5][5000]; // Cluster informations - Int_t fClno; // number of clusters in a supermodule - - /* - clusters[0][i] --- x position of the cluster center - clusters[1][i] --- y position of the cluster center - clusters[2][i] --- total energy in the cluster - clusters[3][i] --- number of cells forming the cluster - ( possibly fractional ) - clusters[4][i] --- cluster radius - */ - + //Variables for association - Int_t fCellTrNo[kNDIMX][kNDIMY]; //id x-y value of cells + Int_t fCellTrNo[kNDIMX][kNDIMY]; // id x-y value of cells Int_t fClTr[15][5000]; // 1d x-y cell info of attached cells Int_t fIord[2][kNMX]; // ordered list of i and j according // to decreasing energy dep. - Int_t fInfocl[2][kNDIMX][kNDIMY]; // cellwise information on the cluster to which the cell + Int_t fInfocl[2][kNDIMX][kNDIMY]; // cellwise information on the + // cluster to which the cell Int_t fInfcl[3][kNMX]; // cluster information [0][i] // -- cluster number Double_t fCoord[2][kNDIMX][kNDIMY]; @@ -106,6 +91,6 @@ class AliPMDClusteringV1: public AliPMDClustering Float_t fCutoff; // Energy(ADC) cutoff per cell before clustering - ClassDef(AliPMDClusteringV1,2) // Does clustering for PMD + ClassDef(AliPMDClusteringV1,3) // Does clustering for PMD }; #endif -- 2.31.1