]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PMD/AliPMDClusteringV1.cxx
Updating DP list for simulation.
[u/mrichter/AliRoot.git] / PMD / AliPMDClusteringV1.cxx
index b96eb60d7d184f13bef2b48eaf6cebec543520aa..3d5582f5767691b2dbf679f40d3328331d243c82 100644 (file)
@@ -113,17 +113,18 @@ void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn,
   Int_t    i,  j, nmx1, incr, id, jd;
   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 = new Double_t [11424];
+  Double_t cellenergy[11424];
   
 
   // call the isolated cell search method
 
-  CalculateIsoCell(idet, ismn, celladc, pmdisocell);
+  FindIsoCell(idet, ismn, celladc, pmdisocell);
 
   // ndimXr and ndimYr are different because of different module size
 
@@ -177,11 +178,9 @@ 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
@@ -245,7 +244,7 @@ void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn,
       clusdata[3]     = cluCELLS;
       clusdata[4]     = cluSIGX;
       clusdata[5]     = cluSIGY;
-      
+
       //
       // Cells associated with a cluster
       //
@@ -269,20 +268,42 @@ void AliPMDClusteringV1::DoClust(Int_t idet, Int_t ismn,
          Int_t irow = celldataX[ihit];
          Int_t icol = celldataY[ihit];
 
-         if ((irow >= 0 && irow < 48) && (icol >= 0 && icol < 96))
+         if (ismn < 12)
            {
-             celldataTr[ihit]  = celltrack[irow][icol];
-             celldataPid[ihit] = cellpid[irow][icol];
+             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
+         else if (ismn >= 12 && ismn < 24)
            {
-             celldataTr[ihit] = -1;
-             celldataPid[ihit] = -1;
+             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,
-                               celldataTr, celldataPid);
+                               celldataTr, celldataPid, celldataAdc);
       pmdcont->Add(pmdcl);
     }
   
@@ -505,6 +526,7 @@ void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
            {
              clxy[icltr] = -1;
            }
+
          pmdcludata  = new AliPMDcludata(clusdata,clxy);
          fPMDclucont->Add(pmdcludata);
        }
@@ -760,14 +782,14 @@ void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
                  
                  if(rr==0)
                    {
-                     ncell[kcl] = 1.;
+                     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;
+                     totaladc2[kcl] = z1*z1;
+                     ax[kcl]        = x1 * z1;
+                     ay[kcl]        = y1 * z1;
+                     ax2[kcl]       = 0.;
+                     ay2[kcl]       = 0.;
+                     status[j]      = 1;
                    }
                }
            }
@@ -805,44 +827,41 @@ void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
                                  
                  cellCount[maxweight]++;
                  
-                 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)
-                     
-                     if(sumweight>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);
-                       }
-                   }
+                 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]++;
+
                }
            }
          
          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];
-             } 
-           
+             
+             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) 
@@ -854,19 +873,21 @@ void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
              clusdata[1] = yclust[kcl];
              clusdata[2] = totaladc[kcl];
              clusdata[3] = ncell[kcl];
+             
+             
              if(sigxclust[kcl] > sigyclust[kcl]) 
                {
-                 clusdata[4] = pow(sigxclust[kcl],0.5);
-                 clusdata[5] = pow(sigyclust[kcl],0.5);
+                 clusdata[4] = TMath::Sqrt(sigxclust[kcl]);
+                 clusdata[5] = TMath::Sqrt(sigyclust[kcl]);
                }
              else
                {
-                 clusdata[4] = pow(sigyclust[kcl],0.5);
-                 clusdata[5] = pow(sigxclust[kcl],0.5);
+                 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++)
                {
@@ -876,7 +897,7 @@ void AliPMDClusteringV1::RefClust(Int_t incr, Double_t edepcell[])
                      Ncell++;
                    }
                } 
-
+             
              pmdcludata = new AliPMDcludata(clusdata,clxy);
              fPMDclucont->Add(pmdcludata);
            }
@@ -919,7 +940,7 @@ Double_t AliPMDClusteringV1::Distance(Double_t x1, Double_t y1,
   return TMath::Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
 }
 // ------------------------------------------------------------------------ //
-void AliPMDClusteringV1::CalculateIsoCell(Int_t idet, Int_t ismn, Double_t celladc[][96], TObjArray *pmdisocell)
+void AliPMDClusteringV1::FindIsoCell(Int_t idet, Int_t ismn, Double_t celladc[][96], TObjArray *pmdisocell)
 {
   // Does isolated cell search for offline calibration
 
@@ -946,8 +967,12 @@ void AliPMDClusteringV1::CalculateIsoCell(Int_t idet, Int_t ismn, Double_t cella
                {
                  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 == 0.)
+                 if(adc < 1.)
                    {
                      isocount++;
                      if(isocount == kCellNeighbour)