Another histos for lumi
[u/mrichter/AliRoot.git] / PMD / AliPMDClusteringV1.cxx
index e03a864..d9a4f9a 100644 (file)
@@ -61,7 +61,8 @@ const Double_t AliPMDClusteringV1::fgkSqroot3by2=0.8660254;  // sqrt(3.)/2.
 
 AliPMDClusteringV1::AliPMDClusteringV1():
   fPMDclucont(new TObjArray()),
-  fCutoff(0.0)
+  fCutoff(0.0),
+  fClusParam(0)
 {
   for(Int_t i = 0; i < kNDIMX; i++)
     {
@@ -76,7 +77,8 @@ AliPMDClusteringV1::AliPMDClusteringV1():
 AliPMDClusteringV1::AliPMDClusteringV1(const AliPMDClusteringV1& pmdclv1):
   AliPMDClustering(pmdclv1),
   fPMDclucont(0),
-  fCutoff(0)
+  fCutoff(0),
+  fClusParam(0)
 {
   // copy constructor
   AliError("Copy constructor not allowed ");
@@ -95,8 +97,11 @@ 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 *pmdcont)
 {
   // main function to call other necessary functions to do clustering
   //
@@ -106,16 +111,16 @@ void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn,
   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    i = 0,  j = 0, nmx1 = 0;
+  Int_t    incr = 0, id = 0, jd = 0;
   Int_t    celldataX[kNmaxCell], celldataY[kNmaxCell];
-  Float_t  clusdata[6];
+  Int_t    celldataTr[kNmaxCell], celldataPid[kNmaxCell];
+  Float_t  celldataAdc[kNmaxCell];
+  Float_t  clusdata[6] = {0.,0.,0.,0.,0.,0.};
   Double_t cutoff, ave;
   Double_t edepcell[kNMX];
+  Double_t cellenergy[kNMX];
   
-  
-  Double_t *cellenergy = new Double_t [11424];
-  
-
   // ndimXr and ndimYr are different because of different module size
 
   Int_t ndimXr = 0;
@@ -132,12 +137,11 @@ void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn,
       ndimYr = 96;
     }
 
-  for (i =0; i < 11424; i++)
+  for (i = 0; i < kNMX; i++)
   {
       cellenergy[i] = 0.;
   }
 
-
   Int_t kk = 0;
   for (i = 0; i < kNDIMX; i++)
     {
@@ -169,15 +173,15 @@ void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn,
     }
   
   for (i = 0; i < kNMX; i++)
-  {
-    edepcell[i] = cellenergy[i];
-  }
-
-  delete [] cellenergy;
+    {
+      edepcell[i] = cellenergy[i];
+    }
 
-  Int_t iord1[kNMX];
-  TMath::Sort((Int_t)kNMX,edepcell,iord1);// order the data
-  cutoff = fCutoff;                       // cutoff to discard cells
+  Bool_t jsort = true;
+  // the dimension of iord1 is increased twice
+  Int_t iord1[2*kNMX];
+  TMath::Sort((Int_t)kNMX,edepcell,iord1,jsort);// order the data
+  cutoff = fCutoff;                             // cutoff to discard cells
   ave  = 0.;
   nmx1 = -1;
   for(i = 0;i < kNMX; i++)
@@ -231,14 +235,13 @@ 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
       //
@@ -258,10 +261,46 @@ void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn,
            }
          
          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];
 
-      pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY);
+               }
+             else
+               {
+                 celldataTr[ihit]  = -1;
+                 celldataPid[ihit] = -1;
+                 celldataAdc[ihit] = -1;
+               }
+           }
+         
+       }
+      
+      pmdcl = new AliPMDcluster(idet, ismn, clusdata, celldataX, celldataY,
+                               celldataTr, celldataPid, celldataAdc);
       pmdcont->Add(pmdcl);
     }
   
@@ -276,8 +315,9 @@ Int_t AliPMDClusteringV1::CrClust(Double_t ave, Double_t cutoff, Int_t nmx1,
   // connected cells
   //
   const Int_t kndim = 4609;
-  Int_t i,j,k,id1,id2,icl, numcell, clust[2][kndim];
-  Int_t jd1,jd2, icell, cellcount;
+  Int_t i=0,j=0,k=0,id1=0,id2=0,icl=0, numcell=0;
+  Int_t jd1=0,jd2=0, icell=0, cellcount=0;
+  Int_t clust[2][kndim];
   static Int_t neibx[6]={1,0,-1,-1,0,1}, neiby[6]={0,1,1,0,-1,-1};
 
   AliDebug(1,Form("kNMX = %d nmx1 = %d kNDIMX = %d kNDIMY = %d ave = %f cutoff = %f",kNMX,nmx1,kNDIMX,kNDIMY,ave,cutoff));
@@ -408,10 +448,13 @@ void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
   
   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;
+  Int_t    i12 = 0, i22 = 0;
+  Int_t    i = 0, j = 0, k = 0;
+  Int_t    i1 = 0, i2 = 0, id = 0, icl = 0;
+  Int_t    itest = 0, ihld = 0, ig = 0;
+  Int_t    nsupcl = 0, clno = 0, t1 = 0, t2 = 0;
   Float_t  clusdata[6];
-  Double_t x1, y1, z1, x2, y2, z2, rr;
+  Double_t x1 = 0, y1 = 0, z1 = 0, x2 = 0, y2 = 0, z2 = 0, rr = 0;
   
   ncl   = new Int_t [ndim];
   clxy  = new Int_t [kNmaxCell];
@@ -420,8 +463,14 @@ void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
   for(i = 0; i<ndim; i++)
     {
       ncl[i]  = -1; 
-      if (i < 6) clusdata[i] = 0.;
-      if (i < kNmaxCell) clxy[i]    = 0;
+    }
+  for(i = 0; i<6; i++)
+    {
+      clusdata[i] = 0.;
+    }
+  for(i = 0; i<19; i++)
+    {
+      clxy[i] = 0;
     }
 
   // clno counts the final clusters
@@ -440,10 +489,10 @@ void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
        {
          nsupcl++;
        }
-      if (nsupcl > ndim) 
+      if (nsupcl >= ndim) 
        {
          AliWarning("RefClust: Too many superclusters!");
-         nsupcl = ndim;
+         nsupcl = ndim - 1;
          break;
        }
       ncl[nsupcl]++;
@@ -463,6 +512,8 @@ void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
          if (clno >= 4608) 
            {
              AliWarning("RefClust: Too many clusters! more than 4608");
+             delete [] ncl;
+             delete [] clxy;
              return;
            }
          clno++;
@@ -475,8 +526,8 @@ void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
          clusdata[1] = fCoord[1][i1][i2];
          clusdata[2] = edepcell[i12];
          clusdata[3] = 1.;
-         clusdata[4] = 0.5;
-         clusdata[5] = 0.0;
+         clusdata[4] = 999.5;
+         clusdata[5] = 999.5;
 
          clxy[0] = i1*10000 + i2;
          
@@ -484,6 +535,7 @@ void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
            {
              clxy[icltr] = -1;
            }
+
          pmdcludata  = new AliPMDcludata(clusdata,clxy);
          fPMDclucont->Add(pmdcludata);
        }
@@ -494,6 +546,9 @@ void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
          if (clno >= 4608) 
            {
              AliWarning("RefClust: Too many clusters! more than 4608");
+             delete [] ncl;
+             delete [] clxy;
+
              return;
            }
          clno++;
@@ -535,14 +590,12 @@ 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];
@@ -584,13 +637,11 @@ void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
              i1      = fInfcl[1][id];
              i2      = fInfcl[2][id];
              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;
-
            }
          
          // arranging cells within supercluster in decreasing order
@@ -612,6 +663,43 @@ void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
                    }
                }
            }
+
+         if (fClusParam == 1)
+           {
+             // Clustering algorithm returns from here for PP collisions
+             // for pp, only the output of crude clusterng is taken
+             // sigx and sigy are not calculated at this moment
+
+             Double_t supx=0., supy=0., supz=0.;
+         
+             for(j = 0;j <= ncl[i]; j++)
+               {
+                 supx += x[iord[j]]*z[iord[j]];
+                 supy += y[iord[j]]*z[iord[j]];
+                 supz += z[iord[j]];
+                 if(j < 19)
+                   {
+                     clxy[j] = t[iord[j]];
+                   }
+               }
+             
+             if( ncl[i] + 1 < 19)
+               {
+                 for(Int_t ncel = ncl[i] + 1; ncel < kNmaxCell; ncel ++ )
+                   {
+                     clxy[ncel] = -1;
+                   }
+               }
+             clusdata[0] = supx/supz;
+             clusdata[1] = supy/supz;
+             clusdata[2] = supz;
+             clusdata[3] = ncl[i]+1;
+             clusdata[4] = 0.5;
+             clusdata[5] = 0.0;
+             pmdcludata  = new AliPMDcludata(clusdata,clxy);
+             fPMDclucont->Add(pmdcludata);
+           }
+
          /* 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
@@ -620,252 +708,300 @@ void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
             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]];
-         tc[ig] = t[iord[0]];
-         Int_t ivalid = 0,  icount = 0;
-         
-         for(j=1;j<=ncl[i];j++)
+
+         if (fClusParam == 2)
            {
-             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 ) 
+             // This part is to split the supercluster
+             //
+             ig = 0;
+             xc[ig] = x[iord[0]];
+             yc[ig] = y[iord[0]];
+             zc[ig] = z[iord[0]];
+             tc[ig] = t[iord[0]];
+             Int_t ivalid = 0,  icount = 0;
+         
+             for(j=1;j<=ncl[i];j++)
                {
-                 ivalid=0;
-                 icount=0;
-                 for(Int_t j1=1;j1<j;j1++)
-                   {
-                     icount++;
-                     Float_t rr1=Distance(x1,y1,x[iord[j1]],y[iord[j1]]);              
-                     if(rr1>1.2) ivalid++;
-                   }
-                 if(ivalid == icount && z1>0.5*zc[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]=z1;
-                     tc[ig]=t1;
-                   }
-               }         
-           }
-         
-         icl=icl+ig+1;
-         
-         //  We use simple Gaussian weighting. (Tapan Jan 2005)
-         // compute the number of cells belonging to each cluster.
-         // 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
-         
+                     ivalid=0;
+                     icount=0;
+                     for(Int_t j1=1;j1<j;j1++)
+                       {
+                         icount++;
+                         Float_t rr1=Distance(x1,y1,x[iord[j1]],y[iord[j1]]);
+                         if(rr1>1.2) ivalid++;
+                       }
+                     if(ivalid == icount && z1>0.5*zc[ig])
+                       {
+                         ig++;
+                         xc[ig]=x1; 
+                         yc[ig]=y1; 
+                         zc[ig]=z1;
+                         tc[ig]=t1;
+                       }
+                   }     
+               }
          
-         Int_t *cellCount;
-         Int_t **cellXY;
+             icl=icl+ig+1;
          
-         Int_t    *status;
-         Double_t *totaladc, *totaladc2, *ncell,*weight;
-         Double_t *xclust, *yclust, *sigxclust, *sigyclust;
-         Double_t *ax, *ay, *ax2, *ay2;
+             //  We use simple Gaussian weighting. (Tapan Jan 2005)
+             // compute the number of cells belonging to each cluster.
+             // 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];
-         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];
+             status    = new Int_t [ncl[i]+1];
+             cellXY    = new Int_t *[ncl[i]+1];
          
-         for(j = 0; j < ncl[i]+1; j++) 
-           {
-             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++)
+             cellCount = new Int_t [ig+1];
+             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++) 
                {
-                 cellXY[j][kcl] = 0;
+                 status[j]     = 0;
+                 cellXY[j] = new Int_t[ig+1];
                }
-           }
-         Double_t sumweight, gweight; 
-         
-         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++)
+             //initialization
+             for(Int_t kcl = 0; kcl < ig+1; kcl++)
                {
-                 x2 = xc[kcl];
-                 y2 = yc[kcl];
-                 rr = Distance(x1,y1,x2,y2);
-                 t2 = tc[kcl];           
-                 
-                 if(rr==0)
+                 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++)
                    {
-                     ncell[kcl] = 1.;
-                     totaladc[kcl]  = z1;
-                     totaladc2[kcl]  = pow(z1,2);
-                     ax[kcl] =  x1 * z1;
-                     ay[kcl] =  y1 * z1;
-                     ax2[kcl]=  0.;
-                     ay2[kcl]=  0.;
-                     status[j] = 1;
+                     cellXY[j][kcl] = 0;
                    }
                }
-           }
-         
-         for(j = 0; j <= ncl[i]; j++)
-           {
-             Int_t   maxweight = 0;
-             Double_t     max  = 0.;
+             Double_t sumweight, gweight; 
              
-             if(status[j] == 0)
-               { 
-                 x1 = x[iord[j]]; 
+             for(j = 0;j <= ncl[i]; j++)
+               {
+                 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++)
-                   {
-                     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)
-                       {
-                         max       =  weight[kcl];
-                         maxweight =  kcl;
-                       }
-                   }
-                 
-                 cellXY[cellCount[maxweight]][maxweight] = iord[j];
-                                 
-                 cellCount[maxweight]++;
                  
-                 for(Int_t kcl = 0; kcl <= ig; kcl++)
+                 for(Int_t kcl=0; kcl<=ig; kcl++)
                    {
                      x2 = xc[kcl];
                      y2 = yc[kcl];
                      rr = Distance(x1,y1,x2,y2);
-                     t2 = tc[kcl];
-                     // For calculating weighted mean:
-                     // Weighted_mean = (\sum w_i x_i) / (\sum w_i)
+                     t2 = tc[kcl];               
                      
-                     if(sumweight>0.)
+                     if(rr==0)
                        {
-                         totaladc[kcl] = totaladc[kcl] + z1*(weight[kcl]/sumweight);
-                         ncell[kcl]    = ncell[kcl]    + (weight[kcl]/sumweight);        
-                         ax[kcl]       = ax[kcl]       + (x1 * z1 *weight[kcl]/sumweight);
-                         ay[kcl]       = ay[kcl]       + (y1 * z1 *weight[kcl]/sumweight);
-                         
-                         // For calculating weighted sigma:
-                         // Wieghted_sigma= ((\sum w_i)/((\sum w_i)^2 - \sum (w_i^2))) \sum w_i (x_i - \Hat\mu)^2
-                         // I assume here x1,y1 are cluster centers, otherwise life becomes too difficult
-                         totaladc2[kcl] = totaladc2[kcl] +  pow((z1*(weight[kcl]/sumweight)),2);
-                         ax2[kcl] = ax2[kcl] + (z1 *weight[kcl]/sumweight) * pow((x1-x2),2);
-                         ay2[kcl] = ay2[kcl] + (z1 *weight[kcl]/sumweight) * pow((y1-y2),2);
+                         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;
                        }
                    }
                }
-           }
          
-         for(Int_t kcl = 0; kcl <= ig; kcl++)
-           {
-
-             if(totaladc[kcl]>0){
-               if(totaladc[kcl]>0.)xclust[kcl] = (ax[kcl])/ totaladc[kcl];
-               if(totaladc[kcl]>0.)yclust[kcl] = (ay[kcl])/ totaladc[kcl];
-               
-               if(totaladc[kcl]>0.)sigxclust[kcl] = (totaladc[kcl]/(pow(totaladc[kcl],2)-totaladc2[kcl]))*ax2[kcl];
-               if(totaladc[kcl]>0.)sigyclust[kcl] = (totaladc[kcl]/(pow(totaladc[kcl],2)-totaladc2[kcl]))*ay2[kcl];
-             } 
-           
-             for(j = 0; j < cellCount[kcl]; j++) clno++; 
-             
-             if (clno >= 4608) 
-               {
-                 AliWarning("RefClust: Too many clusters! more than 4608");
-                 return;
-               }
-             clusdata[0] = xclust[kcl];
-             clusdata[1] = yclust[kcl];
-             clusdata[2] = totaladc[kcl];
-             clusdata[3] = ncell[kcl];
-             if(sigxclust[kcl] > sigyclust[kcl]) 
+             for(j = 0; j <= ncl[i]; j++)
                {
-                 clusdata[4] = pow(sigxclust[kcl],0.5);
-                 clusdata[5] = pow(sigyclust[kcl],0.5);
+                 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++)
+                       {
+                         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)
+                           {
+                             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]++;
+                     
+                   }
                }
-             else
+         
+             for(Int_t kcl = 0; kcl <= ig; kcl++)
                {
-                 clusdata[4] = pow(sigyclust[kcl],0.5);
-                 clusdata[5] = pow(sigxclust[kcl],0.5);
-               }
+                 if(totaladc[kcl] > 0.)
+                   {
+                     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 < cellCount[kcl]; j++) clno++; 
+                 
+                 if (clno >= 4608) 
+                   {
+                     AliWarning("RefClust: Too many clusters! more than 4608");
 
-             clxy[0] = tc[kcl];
+                     delete [] cellCount;
+                     for(Int_t jj = 0; jj < ncl[i]+1; jj++) delete [] cellXY[jj];
+                     delete [] cellXY;
 
-             Int_t Ncell=1;
-             for (Int_t ii = 0; ii < cellCount[kcl]; ii++)
-               {
-                 Float_t xx = x[cellXY[ii][kcl]];
-                 Float_t yy = y[cellXY[ii][kcl]];
-
-                 rr=Distance(xclust[kcl],yclust[kcl],xx,yy);
-                 // Natasha
-                 // We store only the nearest and nest-nearest neighbours
-                 if(rr<2.2) 
-                   {   
-                     clxy[Ncell] = t[cellXY[ii][kcl]];
-                     Ncell++;
+                     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 [] iord;
+                     delete [] tc;       
+                     delete [] t;
+                     delete [] x;
+                     delete [] y;
+                     delete [] z;
+                     delete [] xc;
+                     delete [] yc;
+                     delete [] zc;
+
+
+                     delete [] ncl;
+                     delete [] clxy;
+
+                     return;
+                   }
+                 clusdata[0] = xclust[kcl];
+                 clusdata[1] = yclust[kcl];
+                 clusdata[2] = totaladc[kcl];
+                 clusdata[3] = ncell[kcl];
+                 
+                 if(sigxclust[kcl] > sigyclust[kcl]) 
+                   {
+                     clusdata[4] = TMath::Sqrt(sigxclust[kcl]);
+                     clusdata[5] = TMath::Sqrt(sigyclust[kcl]);
                    }
-               } 
+                 else
+                   {
+                     clusdata[4] = TMath::Sqrt(sigyclust[kcl]);
+                     clusdata[5] = TMath::Sqrt(sigxclust[kcl]);
+                   }
+                 
+                 clxy[0] = tc[kcl];
+                 
+                 Int_t Ncell=1;
+                 for (Int_t ii = 0; ii < cellCount[kcl]; ii++)
+                   {
+                     if(ii<18) 
+                       {       
+                         clxy[Ncell] = t[cellXY[ii][kcl]];
+                         Ncell++;
+                       }
+                   } 
+                 
+                 pmdcludata = new AliPMDcludata(clusdata,clxy);
+                 fPMDclucont->Add(pmdcludata);
+               }
+             delete [] cellCount;
+             for(Int_t jj = 0; jj < ncl[i]+1; jj++) delete [] cellXY[jj];
+             delete [] cellXY;
+
+             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;
 
-             pmdcludata = new AliPMDcludata(clusdata,clxy);
-             fPMDclucont->Add(pmdcludata);
            }
-         
+
          delete [] iord;
          delete [] tc;   
          delete [] t;
@@ -875,23 +1011,8 @@ void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
          delete [] xc;
          delete [] yc;
          delete [] zc;
-         
-         delete [] cellCount;
-         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;
@@ -909,3 +1030,8 @@ void AliPMDClusteringV1::SetEdepCut(Float_t decut)
   fCutoff = decut;
 }
 // ------------------------------------------------------------------------ //
+void AliPMDClusteringV1::SetClusteringParam(Int_t cluspar)
+{
+  fClusParam = cluspar;
+}
+// ------------------------------------------------------------------------ //