X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PMD%2FAliPMDClusteringV1.cxx;h=3d5582f5767691b2dbf679f40d3328331d243c82;hb=cb4de549675559a942698c4053e0f64bd49ec9cf;hp=0c75b7c98fa042cb81b33ed9fbb98198d1d9b31a;hpb=562718f9132ba0846826c005b5475689b32e102d;p=u%2Fmrichter%2FAliRoot.git diff --git a/PMD/AliPMDClusteringV1.cxx b/PMD/AliPMDClusteringV1.cxx index 0c75b7c98fa..3d5582f5767 100644 --- a/PMD/AliPMDClusteringV1.cxx +++ b/PMD/AliPMDClusteringV1.cxx @@ -50,10 +50,12 @@ #include "AliPMDcludata.h" #include "AliPMDcluster.h" +#include "AliPMDisocell.h" #include "AliPMDClustering.h" #include "AliPMDClusteringV1.h" #include "AliLog.h" + ClassImp(AliPMDClusteringV1) const Double_t AliPMDClusteringV1::fgkSqroot3by2=0.8660254; // sqrt(3.)/2. @@ -94,21 +96,35 @@ AliPMDClusteringV1::~AliPMDClusteringV1() delete fPMDclucont; } // ------------------------------------------------------------------------ // -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, + Int_t celltrack[48][96], + Int_t cellpid[48][96], + Double_t celladc[48][96], + TObjArray *pmdisocell, TObjArray *pmdcont) { // main function to call other necessary functions to do clustering // AliPMDcluster *pmdcl = 0; + const float ktwobysqrt3 = 1.1547; // 2./sqrt(3.) + const Int_t kNmaxCell = 19; // # of cells surrounding a cluster center + Int_t i, j, nmx1, incr, id, jd; - Int_t celldataX[15], celldataY[15]; + Int_t celldataX[kNmaxCell], celldataY[kNmaxCell]; + Int_t celldataTr[kNmaxCell], celldataPid[kNmaxCell]; + Float_t celldataAdc[kNmaxCell]; Float_t clusdata[6]; Double_t cutoff, ave; Double_t edepcell[kNMX]; + + + Double_t cellenergy[11424]; + - const float ktwobysqrt3 = 1.1547; // 2./sqrt(3.) + // call the isolated cell search method + + FindIsoCell(idet, ismn, celladc, pmdisocell); // ndimXr and ndimYr are different because of different module size @@ -126,12 +142,16 @@ void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn, ndimYr = 96; } + for (i =0; i < 11424; i++) + { + cellenergy[i] = 0.; + } + Int_t kk = 0; - for (Int_t i = 0; i < kNDIMX; i++) + for (i = 0; i < kNDIMX; i++) { - for (Int_t j = 0; j < kNDIMY; j++) + for (j = 0; j < kNDIMY; j++) { - fCellTrNo[i][j] = -1; edepcell[kk] = 0.; kk++; } @@ -145,24 +165,26 @@ void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn, i = id+(ndimYr/2-1)-(jd/2); Int_t ij = i + j*kNDIMX; - + if (ismn < 12) { - edepcell[ij] = celladc[jd][id]; - fCellTrNo[i][j] = jd*10000+id; // for association + cellenergy[ij] = celladc[jd][id]; } else if (ismn >= 12 && ismn <= 23) { - edepcell[ij] = celladc[id][jd]; - fCellTrNo[i][j] = id*10000+jd; // for association + cellenergy[ij] = celladc[id][jd]; } - } } + for (i = 0; i < kNMX; i++) + { + edepcell[i] = cellenergy[i]; + } + Int_t iord1[kNMX]; - TMath::Sort(kNMX,edepcell,iord1);// order the data - cutoff = fCutoff; // cutoff to discard cells + TMath::Sort((Int_t)kNMX,edepcell,iord1);// order the data + cutoff = fCutoff; // cutoff to discard cells ave = 0.; nmx1 = -1; for(i = 0;i < kNMX; i++) @@ -176,17 +198,15 @@ void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn, nmx1++; } } - + AliDebug(1,Form("Number of cells having energy >= %f are %d",cutoff,nmx1)); if (nmx1 == 0) nmx1 = 1; ave = ave/nmx1; - AliDebug(1,Form("Number of cells in a SuperM = %d and Average = %f", kNMX,ave)); incr = CrClust(ave, cutoff, nmx1,iord1, edepcell ); - RefClust(incr,edepcell); Int_t nentries1 = fPMDclucont->GetEntries(); AliDebug(1,Form("Detector Plane = %d Serial Module No = %d Number of clusters = %d",idet, ismn, nentries1)); @@ -209,6 +229,7 @@ void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn, // // Cluster X centroid is back transformed // + if (ismn < 12) { clusdata[0] = cluX0 - (24-1) + cluY0/2.; @@ -217,36 +238,76 @@ void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn, { clusdata[0] = cluX0 - (48-1) + cluY0/2.; } - + clusdata[1] = cluY0; clusdata[2] = cluADC; clusdata[3] = cluCELLS; clusdata[4] = cluSIGX; clusdata[5] = cluSIGY; - + // // Cells associated with a cluster // - for (Int_t ihit = 0; ihit < 15; ihit++) + for (Int_t ihit = 0; ihit < kNmaxCell; ihit++) { + Int_t cellrow = pmdcludata->GetCellXY(ihit)/10000; + Int_t cellcol = pmdcludata->GetCellXY(ihit)%10000; + if (ismn < 12) { - celldataX[ihit] = pmdcludata->GetCellXY(ihit)%10000; - celldataY[ihit] = pmdcludata->GetCellXY(ihit)/10000; + celldataX[ihit] = cellrow - (24-1) + int(cellcol/2.); } else if (ismn >= 12 && ismn <= 23) { - celldataX[ihit] = pmdcludata->GetCellXY(ihit)/10000; - celldataY[ihit] = pmdcludata->GetCellXY(ihit)%10000; + celldataX[ihit] = cellrow - (48-1) + int(cellcol/2.); + } + + celldataY[ihit] = cellcol; + + Int_t irow = celldataX[ihit]; + Int_t icol = celldataY[ihit]; + + if (ismn < 12) + { + if ((irow >= 0 && irow < 96) && (icol >= 0 && icol < 48)) + { + celldataTr[ihit] = celltrack[icol][irow]; + celldataPid[ihit] = cellpid[icol][irow]; + celldataAdc[ihit] = (Float_t) celladc[icol][irow]; + } + else + { + celldataTr[ihit] = -1; + celldataPid[ihit] = -1; + celldataAdc[ihit] = -1; + } + } + else if (ismn >= 12 && ismn < 24) + { + if ((irow >= 0 && irow < 48) && (icol >= 0 && icol < 96)) + { + celldataTr[ihit] = celltrack[irow][icol]; + celldataPid[ihit] = cellpid[irow][icol]; + celldataAdc[ihit] = (Float_t) celladc[irow][icol]; + + } + else + { + celldataTr[ihit] = -1; + celldataPid[ihit] = -1; + celldataAdc[ihit] = -1; + } } + } - pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY); + + pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY, + celldataTr, celldataPid, celldataAdc); pmdcont->Add(pmdcl); } - fPMDclucont->Clear(); - + fPMDclucont->Delete(); } // ------------------------------------------------------------------------ // Int_t AliPMDClusteringV1::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1, @@ -256,7 +317,8 @@ Int_t AliPMDClusteringV1::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1, // Finds out only the big patch by just searching the // connected cells // - Int_t i,j,k,id1,id2,icl, numcell, clust[2][5000]; + const Int_t kndim = 4609; + Int_t i,j,k,id1,id2,icl, numcell, clust[2][kndim]; Int_t jd1,jd2, icell, cellcount; static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1}; @@ -312,7 +374,7 @@ Int_t AliPMDClusteringV1::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1, clust[0][numcell] = id1; clust[1][numcell] = id2; - for(i = 1; i < 5000; i++) + for(i = 1; i < kndim; i++) { clust[0][i]=0; } @@ -341,7 +403,7 @@ Int_t AliPMDClusteringV1::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1, // check adc count for neighbour's neighbours recursively and // if nonzero, add these to the cluster. // --------------------------------------------------------------- - for(i = 1; i < 5000;i++) + for(i = 1; i < kndim;i++) { if(clust[0][i] != 0) { @@ -379,32 +441,29 @@ void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[]) // Takes the big patch and does gaussian fitting and // finds out the more refined clusters // - + AliPMDcludata *pmdcludata = 0; - - const Int_t kdim = 4500; - - Int_t i, j, k, i1, i2, id, icl, itest,ihld, ig, nsupcl,clno; - Int_t lev1[20], lev2[20]; - Int_t ncl[kdim], iord[kdim], t[kdim]; - Int_t clxy[15]; - - Int_t *cellCount; - Int_t **cellXY; - + + const Int_t kNmaxCell = 19; // # of cells surrounding a cluster center + + Int_t ndim = incr + 1; + + Int_t *ncl = 0x0; + Int_t *clxy = 0x0; + Int_t i12, i22; + Int_t i, j, k, i1, i2, id, icl, itest,ihld, ig, nsupcl,clno, t1, t2; Float_t clusdata[6]; - Double_t x1, y1, z1, x2, y2, z2, dist,rr,sum; - Double_t x[kdim], y[kdim], z[kdim]; - Double_t xc[kdim], yc[kdim], zc[kdim], cells[kdim], rc[kdim]; + Double_t x1, y1, z1, x2, y2, z2, rr; + + ncl = new Int_t [ndim]; + clxy = new Int_t [kNmaxCell]; - // Initialisation - for(i = 0; i kdim) + if (nsupcl > ndim) { AliWarning("RefClust: Too many superclusters!"); - nsupcl = kdim; + nsupcl = ndim; break; } ncl[nsupcl]++; @@ -443,16 +502,16 @@ void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[]) { id++; icl++; - if (clno >= 5000) + if (clno >= 4608) { - AliWarning("RefClust: Too many clusters! more than 5000"); + AliWarning("RefClust: Too many clusters! more than 4608"); return; } clno++; i1 = fInfcl[1][id]; i2 = fInfcl[2][id]; - Int_t i12 = i1 + i2*kNDIMX; + i12 = i1 + i2*kNDIMX; clusdata[0] = fCoord[0][i1][i2]; clusdata[1] = fCoord[1][i1][i2]; @@ -461,11 +520,13 @@ void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[]) clusdata[4] = 0.5; clusdata[5] = 0.0; - clxy[0] = fCellTrNo[i1][i2]; //association - for(Int_t icltr = 1; icltr < 15; icltr++) + clxy[0] = i1*10000 + i2; + + for(Int_t icltr = 1; icltr < kNmaxCell; icltr++) { clxy[icltr] = -1; } + pmdcludata = new AliPMDcludata(clusdata,clxy); fPMDclucont->Add(pmdcludata); } @@ -473,30 +534,35 @@ void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[]) { id++; icl++; - if (clno >= 5000) + if (clno >= 4608) { - AliWarning("RefClust: Too many clusters! more than 5000"); + AliWarning("RefClust: Too many clusters! more than 4608"); return; } clno++; i1 = fInfcl[1][id]; i2 = fInfcl[2][id]; - Int_t i12 = i1 + i2*kNDIMX; + i12 = i1 + i2*kNDIMX; x1 = fCoord[0][i1][i2]; y1 = fCoord[1][i1][i2]; z1 = edepcell[i12]; - clxy[0] = fCellTrNo[i1][i2]; //asso + + clxy[0] = i1*10000 + i2; + id++; i1 = fInfcl[1][id]; i2 = fInfcl[2][id]; - Int_t i22 = i1 + i2*kNDIMX; + i22 = i1 + i2*kNDIMX; x2 = fCoord[0][i1][i2]; y2 = fCoord[1][i1][i2]; z2 = edepcell[i22]; - clxy[1] = fCellTrNo[i1][i2]; //asso - for(Int_t icltr = 2; icltr < 15; icltr++) + + clxy[1] = i1*10000 + i2; + + + for(Int_t icltr = 2; icltr < kNmaxCell; icltr++) { clxy[icltr] = -1; } @@ -512,36 +578,62 @@ void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[]) } else { + + Int_t *iord, *tc, *t; + Double_t *x, *y, *z, *xc, *yc, *zc; + + iord = new Int_t [ncl[i]+1]; + tc = new Int_t [ncl[i]+1]; + t = new Int_t [ncl[i]+1]; + + x = new Double_t [ncl[i]+1]; + y = new Double_t [ncl[i]+1]; + z = new Double_t [ncl[i]+1]; + xc = new Double_t [ncl[i]+1]; + yc = new Double_t [ncl[i]+1]; + zc = new Double_t [ncl[i]+1]; + + for( k = 0; k < ncl[i]+1; k++) + { + iord[k] = -1; + t[k] = -1; + tc[k] = -1; + x[k] = -1; + y[k] = -1; + z[k] = -1; + xc[k] = -1; + yc[k] = -1; + zc[k] = -1; + } 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]; - Int_t i12 = i1 + i2*kNDIMX; + i12 = i1 + i2*kNDIMX; x[0] = fCoord[0][i1][i2]; y[0] = fCoord[1][i1][i2]; z[0] = edepcell[i12]; - - t[0] = fCellTrNo[i1][i2]; //asso + t[0] = i1*10000 + i2; + iord[0] = 0; for(j = 1; j <= ncl[i]; j++) { id++; i1 = fInfcl[1][id]; i2 = fInfcl[2][id]; - Int_t i12 = i1 + i2*kNDIMX; + i12 = i1 + i2*kNDIMX; iord[j] = j; x[j] = fCoord[0][i1][i2]; y[j] = fCoord[1][i1][i2]; z[j] = edepcell[i12]; + t[j] = i1*10000 + i2; - t[j] = fCellTrNo[i1][i2]; //asso } // arranging cells within supercluster in decreasing order @@ -563,317 +655,341 @@ void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[]) } } } - // 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 + /* MODIFICATION PART STARTS (Tapan July 2008) + iord[0] is the cell with highest ADC in the crude-cluster + ig is the number of local maxima in the crude-cluster + For the higest peak we make ig=0 which means first local maximum. + Next we go down in terms of the ADC sequence and find out if any + more of the cells form local maxima. The definition of local + maxima is that all its neighbours are of less ADC compared to it. + */ ig = 0; xc[ig] = x[iord[0]]; yc[ig] = y[iord[0]]; zc[ig] = z[iord[0]]; - for(j = 1; j <= ncl[i]; j++) + tc[ig] = t[iord[0]]; + Int_t ivalid = 0, icount = 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) + x1 = x[iord[j]]; + y1 = y[iord[j]]; + z1 = z[iord[j]]; + t1 = t[iord[j]]; + rr=Distance(x1,y1,xc[ig],yc[ig]); + + // Check the cells which are outside the neighbours (rr>1.2) + if(rr>1.2 ) { - ig++; - xc[ig] = x1; - yc[ig] = y1; - zc[ig] = z[iord[j]]; - } + ivalid=0; + icount=0; + for(Int_t j1=1;j11.2) ivalid++; + } + if(ivalid == icount && z1>0.5*zc[ig]) + { + ig++; + xc[ig]=x1; + yc[ig]=y1; + zc[ig]=z1; + tc[ig]=t1; + } + } } - GaussFit(ncl[i], ig, x[0], y[0] ,z[0], xc[0], yc[0], zc[0], rc[0]); - icl += ig+1; + + icl=icl+ig+1; + + // We use simple Gaussian weighting. (Tapan Jan 2005) // 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 - - Int_t jj = 15; + // cell can be shared between several clusters + // in the ratio of cluster energy deposition + // To calculate: + // (1) number of cells belonging to a cluster (ig) and + // (2) total ADC of the cluster (ig) + // (3) x and y positions of the cluster + + + Int_t *cellCount; + Int_t **cellXY; + + Int_t *status; + Double_t *totaladc, *totaladc2, *ncell,*weight; + Double_t *xclust, *yclust, *sigxclust, *sigyclust; + Double_t *ax, *ay, *ax2, *ay2; + + + status = new Int_t [ncl[i]+1]; + cellXY = new Int_t *[ncl[i]+1]; + cellCount = new Int_t [ig+1]; - cellXY = new Int_t *[jj]; - for(Int_t ij = 0; ij < 15; ij++) cellXY[ij] = new Int_t [ig+1]; - - for(j = 0; j <= ig; j++) + totaladc = new Double_t [ig+1]; + totaladc2 = new Double_t [ig+1]; + ncell = new Double_t [ig+1]; + weight = new Double_t [ig+1]; + xclust = new Double_t [ig+1]; + yclust = new Double_t [ig+1]; + sigxclust = new Double_t [ig+1]; + sigyclust = new Double_t [ig+1]; + ax = new Double_t [ig+1]; + ay = new Double_t [ig+1]; + ax2 = new Double_t [ig+1]; + ay2 = new Double_t [ig+1]; + + for(j = 0; j < ncl[i]+1; j++) { - cellCount[j] = 0; - cells[j] = 0.; + status[j] = 0; + cellXY[j] = new Int_t[ig+1]; } + //initialization + for(Int_t kcl = 0; kcl < ig+1; kcl++) + { + cellCount[kcl] = 0; + totaladc[kcl] = 0.; + totaladc2[kcl] = 0.; + ncell[kcl] = 0.; + weight[kcl] = 0.; + xclust[kcl] = 0.; + yclust[kcl] = 0.; + sigxclust[kcl] = 0.; + sigyclust[kcl] = 0.; + ax[kcl] = 0.; + ay[kcl] = 0.; + ax2[kcl] = 0.; + ay2[kcl] = 0.; + for(j = 0; j < ncl[i]+1; j++) + { + cellXY[j][kcl] = 0; + } + } + Double_t sumweight, gweight; - if(ig > 0) + for(j = 0;j <= ncl[i]; j++) { - for(j = 0; j <= ncl[i]; j++) + x1 = x[iord[j]]; + y1 = y[iord[j]]; + z1 = z[iord[j]]; + t1 = t[iord[j]]; + + for(Int_t kcl=0; kcl<=ig; kcl++) { - 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 - if (cellCount[k] < 15) - { - cellXY[cellCount[k]][k] = 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) + x2 = xc[kcl]; + y2 = yc[kcl]; + rr = Distance(x1,y1,x2,y2); + t2 = tc[kcl]; + + if(rr==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; - } - } + ncell[kcl] = 1.; + totaladc[kcl] = z1; + totaladc2[kcl] = z1*z1; + ax[kcl] = x1 * z1; + ay[kcl] = y1 * z1; + ax2[kcl] = 0.; + ay2[kcl] = 0.; + status[j] = 1; } - else + } + } + + for(j = 0; j <= ncl[i]; j++) + { + Int_t maxweight = 0; + Double_t max = 0.; + + if(status[j] == 0) + { + x1 = x[iord[j]]; + y1 = y[iord[j]]; + z1 = z[iord[j]]; + t1 = t[iord[j]]; + sumweight = 0.; + + for(Int_t kcl = 0; kcl <= ig; kcl++) { - if(lev2[0] == 0) - { - cells[lev2[1]]++; - } - else + x2 = xc[kcl]; + y2 = yc[kcl]; + rr = Distance(x1,y1,x2,y2); + gweight = exp(-(rr*rr)/(2*(1.2*1.2))); + weight[kcl] = zc[kcl] * gweight; + sumweight = sumweight + weight[kcl]; + + if(weight[kcl] > max) { - 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; - } + max = weight[kcl]; + maxweight = kcl; } } + + cellXY[cellCount[maxweight]][maxweight] = iord[j]; + + cellCount[maxweight]++; + + x2 = xc[maxweight]; + y2 = yc[maxweight]; + totaladc[maxweight] += z1; + ax[maxweight] += x1*z1; + ay[maxweight] += y1*z1; + totaladc2[maxweight] += z1*z1; + ax2[maxweight] += z1*(x1-x2)*(x1-x2); + ay2[maxweight] += z1*(y1-y2)*(y1-y2); + ncell[maxweight]++; + } } - // zero rest of the cell array - //asso - for( k = 0; k <= ig; k++) + for(Int_t kcl = 0; kcl <= ig; kcl++) { - for(Int_t icltr = cellCount[k]; icltr < 15; icltr++) + + if(totaladc[kcl] > 0.) { - cellXY[icltr][k] = -1; + xclust[kcl] = (ax[kcl])/ totaladc[kcl]; + yclust[kcl] = (ay[kcl])/ totaladc[kcl]; + + //natasha + Float_t sqtotadc = totaladc[kcl]*totaladc[kcl]; + if(totaladc2[kcl] >= sqtotadc) + { + sigxclust[kcl] = 0.25; + sigyclust[kcl] = 0.25; + } + else + { + sigxclust[kcl] = (totaladc[kcl]/(sqtotadc-totaladc2[kcl]))*ax2[kcl]; + sigyclust[kcl] = (totaladc[kcl]/(sqtotadc-totaladc2[kcl]))*ay2[kcl]; + } } - } - // - - for(j = 0; j <= ig; j++) - { - clno++; - if (clno >= 5000) + + for(j = 0; j < cellCount[kcl]; j++) clno++; + + if (clno >= 4608) { - AliWarning("RefClust: Too many clusters! more than 5000"); + AliWarning("RefClust: Too many clusters! more than 4608"); 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[0] = xclust[kcl]; + clusdata[1] = yclust[kcl]; + clusdata[2] = totaladc[kcl]; + clusdata[3] = ncell[kcl]; + + + if(sigxclust[kcl] > sigyclust[kcl]) { - clusdata[3] = ncl[i]; + clusdata[4] = TMath::Sqrt(sigxclust[kcl]); + clusdata[5] = TMath::Sqrt(sigyclust[kcl]); } else { - clusdata[3] = cells[j]; + clusdata[4] = TMath::Sqrt(sigyclust[kcl]); + clusdata[5] = TMath::Sqrt(sigxclust[kcl]); } - - - for (Int_t ii=0; ii < 15; ii++) + + clxy[0] = tc[kcl]; + + Int_t Ncell=1; + for (Int_t ii = 0; ii < cellCount[kcl]; ii++) { - clxy[ii] = cellXY[ii][j]; - } + if(ii<18) + { + clxy[Ncell] = t[cellXY[ii][kcl]]; + Ncell++; + } + } + pmdcludata = new AliPMDcludata(clusdata,clxy); fPMDclucont->Add(pmdcludata); } + + delete [] iord; + delete [] tc; + delete [] t; + delete [] x; + delete [] y; + delete [] z; + delete [] xc; + delete [] yc; + delete [] zc; + delete [] cellCount; - for(Int_t jj = 0; jj < 15; jj++)delete [] cellXY[jj]; - delete [] cellXY; + for(Int_t jj = 0; jj < ncl[i]+1; jj++) delete [] cellXY[jj]; + + delete [] status; + delete [] totaladc; + delete [] totaladc2; + delete [] ncell; + delete [] xclust; + delete [] yclust; + delete [] sigxclust; + delete [] sigyclust; + delete [] ax; + delete [] ay; + delete [] ax2; + delete [] ay2; + delete [] weight; } } + delete [] ncl; + delete [] clxy; } // ------------------------------------------------------------------------ // -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) +Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1, + Double_t x2, Double_t y2) { - // Does gaussian fitting - // + return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); +} +// ------------------------------------------------------------------------ // +void AliPMDClusteringV1::FindIsoCell(Int_t idet, Int_t ismn, Double_t celladc[][96], TObjArray *pmdisocell) +{ + // Does isolated cell search for offline calibration - const Int_t kdim = 4500; - Int_t i, j, i1, i2, novar, idd, jj; - Int_t neib[kdim][50]; + AliPMDisocell *isocell = 0; - 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]; - - TRandom rnd; - - str = 0.; - str1 = 0.; - rr = 0.3; - novar = 0; - j = 0; + const Int_t kMaxRow = 48; + const Int_t kMaxCol = 96; + const Int_t kCellNeighbour = 6; - for(i = 0; i <= ncell; i++) - { - xx[i] = *(&x+i); - yy[i] = *(&y+i); - zz[i] = *(&z+i); - str += zz[i]; - } - for(i=0; i<=nclust; i++) - { - xxc[i] = *(&xc+i); - yyc[i] = *(&yc+i); - zzc[i] = *(&zc+i); - str1 += zzc[i]; - rrc[i] = 0.5; - } - for(i = 0; i <= nclust; i++) - { - zzc[i] = str/str1*zzc[i]; - ha[i] = xxc[i]; - hb[i] = yyc[i]; - hc[i] = zzc[i]; - hd[i] = rrc[i]; - x1 = xxc[i]; - y1 = yyc[i]; - } - - for(i = 0; i <= ncell; i++) + Int_t id1, jd1; + + Int_t neibx[6] = {1,0,-1,-1,0,1}; + Int_t neiby[6] = {0,1,1,0,-1,-1}; + + + for(Int_t irow = 0; irow < kMaxRow; irow++) { - idd = 0; - x1 = xx[i]; - y1 = yy[i]; - for(j = 0; j <= nclust; j++) + for(Int_t icol = 0; icol < kMaxCol; icol++) { - x2 = xxc[j]; - y2 = yyc[j]; - if(Distance(x1,y1,x2,y2) <= 3.) - { - idd++; - neib[i][idd] = j; + if(celladc[irow][icol] > 0) + { + Int_t isocount = 0; + for(Int_t ii = 0; ii < kCellNeighbour; ii++) + { + id1 = irow + neibx[ii]; + jd1 = icol + neiby[ii]; + if (id1 < 0) id1 = 0; + if (id1 > kMaxRow-1) id1 = kMaxRow - 1; + if (jd1 < 0) jd1 = 0; + if (jd1 > kMaxCol-1) jd1 = kMaxCol - 1; + Float_t adc = (Float_t) celladc[id1][jd1]; + if(adc < 1.) + { + isocount++; + if(isocount == kCellNeighbour) + { + Float_t cadc = (Float_t) celladc[irow][icol]; + + isocell = new AliPMDisocell(idet,ismn,irow,icol,cadc); + pmdisocell->Add(isocell); + + } + } + } // neigh cell cond. } } - 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 += 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; - } - } - for(i = 0; i <= nclust; i++) - { - c[i] = c[i]*str/str1; - } - sum1=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] - a[jj]; - dy = yy[i1] - b[jj]; - dum = d[jj]*d[jj]+rr*rr; - aint += exp(-(dx*dx+dy*dy)/dum)*c[i2]*rr*rr/dum; - } - sum1 += (aint - zz[i1])*(aint - zz[i1])/str; } - if(sum1 < sum) - { - for(i2 = 0; i2 <= nclust; i2++) - { - xxc[i2] = a[i2]; - yyc[i2] = b[i2]; - zzc[i2] = c[i2]; - rrc[i2] = d[i2]; - sum = sum1; - } - } - for(j = 0; j <= nclust; j++) - { - *(&xc+j) = xxc[j]; - *(&yc+j) = yyc[j]; - *(&zc+j) = zzc[j]; - *(&rc+j) = rrc[j]; - } -} -// ------------------------------------------------------------------------ // -Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1, - Double_t x2, Double_t y2) -{ - return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); + } // ------------------------------------------------------------------------ // void AliPMDClusteringV1::SetEdepCut(Float_t decut)