A cluster finder and hit reconstruction class for RICH (adapted from MUON).
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 19 Apr 2000 13:02:04 +0000 (13:02 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 19 Apr 2000 13:02:04 +0000 (13:02 +0000)
Cluster Finders for MUON and RICH should derive from the same class in the
future (JB, AM).

RICH/AliRICHClusterFinder.cxx [new file with mode: 0644]
RICH/AliRICHClusterFinder.h [new file with mode: 0644]

diff --git a/RICH/AliRICHClusterFinder.cxx b/RICH/AliRICHClusterFinder.cxx
new file mode 100644 (file)
index 0000000..f03b98f
--- /dev/null
@@ -0,0 +1,1054 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/*
+  $Log$
+*/
+
+
+#include "AliRICHClusterFinder.h"
+#include "TTree.h"
+#include "AliRun.h"
+#include <TCanvas.h>
+#include <TH1.h>
+#include <TPad.h>
+#include <TGraph.h> 
+#include <TPostScript.h> 
+#include <TMinuit.h> 
+
+//----------------------------------------------------------
+static AliRICHSegmentation* gSegmentation;
+static AliRICHResponse*     gResponse;
+static Int_t                gix[500];
+static Int_t                giy[500];
+static Float_t              gCharge[500];
+static Int_t                gNbins;
+static Int_t                gFirst=kTRUE;
+static TMinuit *gMyMinuit ;
+void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
+static Int_t                gChargeTot;
+
+ClassImp(AliRICHClusterFinder)
+
+    AliRICHClusterFinder::AliRICHClusterFinder
+(AliRICHSegmentation *segmentation, AliRICHResponse *response, 
+ TClonesArray *digits, Int_t chamber)   
+{
+    fSegmentation=segmentation;
+    fResponse=response;
+    
+    fDigits=digits;
+    fNdigits = fDigits->GetEntriesFast();
+    fChamber=chamber;
+    fRawClusters=new TClonesArray("AliRICHRawCluster",10000);
+    fNRawClusters=0;
+    fCogCorr = 0;
+    SetNperMax();
+    SetClusterSize();
+    SetDeclusterFlag();
+    fNPeaks=-1;
+}
+
+    AliRICHClusterFinder::AliRICHClusterFinder()
+{
+    fSegmentation=0;
+    fResponse=0;
+    
+    fDigits=0;
+    fNdigits = 0;
+    fChamber=-1;
+    fRawClusters=new TClonesArray("AliRICHRawCluster",10000);
+    fNRawClusters=0;
+    fHitMap = 0;
+    fCogCorr = 0;
+    SetNperMax();
+    SetClusterSize();
+    SetDeclusterFlag();
+    fNPeaks=-1;
+}
+
+void AliRICHClusterFinder::AddRawCluster(const AliRICHRawCluster c)
+{
+  //
+  // Add a raw cluster copy to the list
+  //
+    AliRICH *RICH=(AliRICH*)gAlice->GetModule("RICH");
+    RICH->AddRawCluster(fChamber,c); 
+    fNRawClusters++;
+}
+
+
+
+void AliRICHClusterFinder::Decluster(AliRICHRawCluster *cluster)
+{
+//    AliRICHDigit *dig;
+//    Int_t q;
+
+    
+    Int_t mul = cluster->fMultiplicity;
+//    printf("Decluster - multiplicity   %d \n",mul);
+
+    if (mul == 1 || mul ==2) {
+//
+// Nothing special for 1- and 2-clusters
+       if (fNPeaks != 0) {
+            cluster->fNcluster[0]=fNPeaks;
+            cluster->fNcluster[1]=0;
+        } 
+       AddRawCluster(*cluster); 
+        fNPeaks++;
+    } else if (mul ==3) {
+//
+// 3-cluster, check topology
+//     printf("\n 3-cluster, check topology \n");
+       if (fDeclusterFlag) {
+           if (Centered(cluster)) {
+               // ok, cluster is centered 
+           } else {
+               // cluster is not centered, split into 2+1
+           }
+       } else {
+           if (fNPeaks != 0) {
+               cluster->fNcluster[0]=fNPeaks;
+               cluster->fNcluster[1]=0;
+           } 
+           AddRawCluster(*cluster); 
+           fNPeaks++;
+       }           
+    } else {
+// 
+// 4-and more-pad clusters
+//
+       if (mul <= fClusterSize) {
+           if (fDeclusterFlag) {
+               SplitByLocalMaxima(cluster);
+           } else {
+               if (fNPeaks != 0) {
+                   cluster->fNcluster[0]=fNPeaks;
+                   cluster->fNcluster[1]=0;
+               } 
+               AddRawCluster(*cluster);
+               fNPeaks++;
+           }   
+       }
+    } // multiplicity 
+}
+
+
+Bool_t AliRICHClusterFinder::Centered(AliRICHRawCluster *cluster)
+{
+    AliRICHDigit* dig;
+    dig= (AliRICHDigit*)fDigits->UncheckedAt(cluster->fIndexMap[0]);
+    Int_t ix=dig->fPadX;
+    Int_t iy=dig->fPadY;
+    Int_t nn;
+    Int_t X[kMaxNeighbours], Y[kMaxNeighbours], XN[kMaxNeighbours], YN[kMaxNeighbours];
+    
+    fSegmentation->Neighbours(ix,iy,&nn,X,Y);
+    Int_t nd=0;
+    for (Int_t i=0; i<nn; i++) {
+       if (fHitMap->TestHit(X[i],Y[i]) == used) {
+           XN[nd]=X[i];
+           YN[nd]=Y[i];
+           nd++;
+       }
+    }
+    if (nd==2) {
+//
+// cluster is centered !
+       if (fNPeaks != 0) {
+            cluster->fNcluster[0]=fNPeaks;
+            cluster->fNcluster[1]=0;
+        }  
+       cluster->fCtype=0;
+       AddRawCluster(*cluster);
+       fNPeaks++;
+       return kTRUE;
+    } else if (nd ==1) {
+//
+// Highest signal on an edge, split cluster into 2+1
+//
+// who is the neighbour ?
+       Int_t nind=fHitMap->GetHitIndex(XN[0], YN[0]);
+       Int_t i1= (nind==cluster->fIndexMap[1]) ? 1:2;
+       Int_t i2= (nind==cluster->fIndexMap[1]) ? 2:1;    
+//
+// 2-cluster
+       AliRICHRawCluster cnew;
+       if (fNPeaks == 0) {
+            cnew.fNcluster[0]=-1;
+            cnew.fNcluster[1]=fNRawClusters;
+        } else {
+            cnew.fNcluster[0]=fNPeaks;
+            cnew.fNcluster[1]=0;
+        }
+       cnew.fMultiplicity=2;
+       cnew.fIndexMap[0]=cluster->fIndexMap[0];
+       cnew.fIndexMap[1]=cluster->fIndexMap[i1];
+       FillCluster(&cnew);
+       cnew.fClusterType=cnew.PhysicsContribution();
+       AddRawCluster(cnew);
+        fNPeaks++;
+//
+// 1-cluster
+       cluster->fMultiplicity=1;
+       cluster->fIndexMap[0]=cluster->fIndexMap[i2];
+       cluster->fIndexMap[1]=0;
+       cluster->fIndexMap[2]=0;        
+       FillCluster(cluster);
+        if (fNPeaks != 0) {
+            cluster->fNcluster[0]=fNPeaks;
+            cluster->fNcluster[1]=0;
+        }  
+       cluster->fClusterType=cluster->PhysicsContribution();
+       AddRawCluster(*cluster);
+       fNPeaks++;
+       return kFALSE;
+    } else {
+       printf("\n Completely screwed up %d !! \n",nd);
+       
+    }
+    
+       return kFALSE;
+}
+void AliRICHClusterFinder::SplitByLocalMaxima(AliRICHRawCluster *c)
+{
+    AliRICHDigit* dig[100], *digt;
+    Int_t ix[100], iy[100], q[100];
+    Float_t x[100], y[100];
+    Int_t i; // loops over digits
+    Int_t j; // loops over local maxima
+    //    Float_t xPeak[2];
+    //    Float_t yPeak[2];
+    //    Int_t threshold=500;
+    Int_t mul=c->fMultiplicity;
+//
+//  dump digit information into arrays
+//
+    for (i=0; i<mul; i++)
+    {
+       dig[i]= (AliRICHDigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
+       ix[i]= dig[i]->fPadX;
+       iy[i]= dig[i]->fPadY;
+       q[i] = dig[i]->fSignal;
+       fSegmentation->GetPadCxy(ix[i], iy[i], x[i], y[i]);
+    }
+//
+//  Find local maxima
+//
+    Bool_t IsLocal[100];
+    Int_t NLocal=0;
+    Int_t AssocPeak[100];
+    Int_t IndLocal[100];
+    Int_t nn;
+    Int_t X[kMaxNeighbours], Y[kMaxNeighbours];
+    for (i=0; i<mul; i++) {
+       fSegmentation->Neighbours(ix[i], iy[i], &nn, X, Y);
+       IsLocal[i]=kTRUE;
+       for (j=0; j<nn; j++) {
+           if (fHitMap->TestHit(X[j], Y[j])==empty) continue;
+           digt=(AliRICHDigit*) fHitMap->GetHit(X[j], Y[j]);
+           if (digt->fSignal > q[i]) {
+               IsLocal[i]=kFALSE;
+               break;
+//
+// handle special case of neighbouring pads with equal signal
+           } else if (digt->fSignal == q[i]) {
+               if (NLocal >0) {
+                   for (Int_t k=0; k<NLocal; k++) {
+                       if (X[j]==ix[IndLocal[k]] && Y[j]==iy[IndLocal[k]]){
+                           IsLocal[i]=kFALSE;
+                       }
+                   }
+               }
+           } 
+       } // loop over next neighbours
+       // Maxima should not be on the edge
+       if (IsLocal[i]) {
+           IndLocal[NLocal]=i;
+           NLocal++;
+       } 
+    } // loop over all digits
+//    printf("Found %d local Maxima",NLocal);
+//
+// If only one local maximum found but multiplicity is high 
+// take global maximum from the list of digits.    
+    if (NLocal==1 && mul>5) {
+       Int_t nnew=0;
+       for (i=0; i<mul; i++) {
+           if (!IsLocal[i]) {
+               IndLocal[NLocal]=i;
+               IsLocal[i]=kTRUE;
+               NLocal++;
+               nnew++;
+           }
+           if (nnew==1) break;
+       }
+    }
+    
+// If number of local maxima is 2 try to fit a double gaussian
+    if (NLocal==-100) {
+//
+//  Initialise global variables for fit
+       gFirst=1;
+       gSegmentation=fSegmentation;
+       gResponse    =fResponse;
+       gNbins=mul;
+       
+       for (i=0; i<mul; i++) {
+           gix[i]=ix[i];
+           giy[i]=iy[i];
+           gCharge[i]=Float_t(q[i]);
+       }
+//
+       if (gFirst) {
+           gFirst=kFALSE;
+           gMyMinuit = new TMinuit(5);
+       }
+       gMyMinuit->SetFCN(fcn);
+       gMyMinuit->mninit(5,10,7);
+       Double_t arglist[20];
+       Int_t ierflag=0;
+       arglist[0]=1;
+//     gMyMinuit->mnexcm("SET ERR",arglist,1,ierflag);
+// Set starting values 
+       static Double_t vstart[5];
+       vstart[0]=x[IndLocal[0]];
+       vstart[1]=y[IndLocal[0]];       
+       vstart[2]=x[IndLocal[1]];
+       vstart[3]=y[IndLocal[1]];       
+       vstart[4]=Float_t(q[IndLocal[0]])/
+           Float_t(q[IndLocal[0]]+q[IndLocal[1]]);
+// lower and upper limits
+       static Double_t lower[5], upper[5];
+       Int_t isec=fSegmentation->Sector(ix[IndLocal[0]], iy[IndLocal[0]]);
+       lower[0]=vstart[0]-fSegmentation->Dpx(isec)/2;
+       lower[1]=vstart[1]-fSegmentation->Dpy(isec)/2;
+//     lower[1]=vstart[1];
+       
+       upper[0]=lower[0]+fSegmentation->Dpx(isec);
+       upper[1]=lower[1]+fSegmentation->Dpy(isec);
+//     upper[1]=vstart[1];
+       
+       isec=fSegmentation->Sector(ix[IndLocal[1]], iy[IndLocal[1]]);
+       lower[2]=vstart[2]-fSegmentation->Dpx(isec)/2;
+       lower[3]=vstart[3]-fSegmentation->Dpy(isec)/2;
+//     lower[3]=vstart[3];
+       
+       upper[2]=lower[2]+fSegmentation->Dpx(isec);
+       upper[3]=lower[3]+fSegmentation->Dpy(isec);
+//     upper[3]=vstart[3];
+       
+       lower[4]=0.;
+       upper[4]=1.;
+// step sizes
+       static Double_t step[5]={0.005, 0.03, 0.005, 0.03, 0.01};
+       
+       gMyMinuit->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
+       gMyMinuit->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
+       gMyMinuit->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
+       gMyMinuit->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
+       gMyMinuit->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
+// ready for minimisation      
+       gMyMinuit->SetPrintLevel(-1);
+       gMyMinuit->mnexcm("SET OUT", arglist, 0, ierflag);
+       arglist[0]= -1;
+       arglist[1]= 0;
+       
+       gMyMinuit->mnexcm("SET NOGR", arglist, 0, ierflag);
+       gMyMinuit->mnexcm("SCAN", arglist, 0, ierflag);
+       gMyMinuit->mnexcm("EXIT" , arglist, 0, ierflag);
+// Print results
+//     Double_t amin,edm,errdef;
+//     Int_t nvpar,nparx,icstat;
+//     gMyMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
+//     gMyMinuit->mnprin(3,amin);
+// Get fitted parameters
+
+       Double_t xrec[2], yrec[2], qfrac;
+       TString chname;
+       Double_t epxz, b1, b2;
+       Int_t ierflg;
+       gMyMinuit->mnpout(0, chname, xrec[0], epxz, b1, b2, ierflg);    
+       gMyMinuit->mnpout(1, chname, yrec[0], epxz, b1, b2, ierflg);    
+       gMyMinuit->mnpout(2, chname, xrec[1], epxz, b1, b2, ierflg);    
+       gMyMinuit->mnpout(3, chname, yrec[1], epxz, b1, b2, ierflg);    
+       gMyMinuit->mnpout(4, chname, qfrac,   epxz, b1, b2, ierflg);    
+       //printf("\n %f %f %f %f %f\n", xrec[0], yrec[0], xrec[1], yrec[1],qfrac);
+//     delete gMyMinuit;
+       
+       
+ //
+ // One cluster for each maximum
+ //
+       for (j=0; j<2; j++) {
+           AliRICHRawCluster cnew;
+           if (fNPeaks == 0) {
+               cnew.fNcluster[0]=-1;
+               cnew.fNcluster[1]=fNRawClusters;
+           } else {
+               cnew.fNcluster[0]=fNPeaks;
+               cnew.fNcluster[1]=0;
+           }
+           cnew.fMultiplicity=0;
+           cnew.fX=Float_t(xrec[j]);
+           cnew.fY=Float_t(yrec[j]);
+           if (j==0) {
+               cnew.fQ=Int_t(gChargeTot*qfrac);
+           } else {
+               cnew.fQ=Int_t(gChargeTot*(1-qfrac));
+           }
+           gSegmentation->SetHit(xrec[j],yrec[j]);
+           for (i=0; i<mul; i++) {
+               cnew.fIndexMap[cnew.fMultiplicity]=c->fIndexMap[i];
+               gSegmentation->SetPad(gix[i], giy[i]);
+               Float_t q1=gResponse->IntXY(gSegmentation);
+               cnew.fContMap[cnew.fMultiplicity]=Float_t(q[i])/(q1*cnew.fQ);
+               cnew.fMultiplicity++;
+           }
+           FillCluster(&cnew,0);
+           //printf("\n x,y %f %f ", cnew.fX, cnew.fY);
+           cnew.fClusterType=cnew.PhysicsContribution();
+           AddRawCluster(cnew);
+           fNPeaks++;
+       }
+    }
+
+    Bool_t fitted=kTRUE;
+
+    if (NLocal !=-100 || !fitted) {
+       // Check if enough local clusters have been found,
+       // if not add global maxima to the list 
+       //
+       Int_t nPerMax;
+       if (NLocal!=0) {
+           nPerMax=mul/NLocal;
+       } else {
+           printf("\n Warning, no local maximum found \n");
+           nPerMax=fNperMax+1;
+       }
+       
+       if (nPerMax > fNperMax) {
+           Int_t nGlob=mul/fNperMax-NLocal+1;
+           if (nGlob > 0) {
+               Int_t nnew=0;
+               for (i=0; i<mul; i++) {
+                   if (!IsLocal[i]) {
+                       IndLocal[NLocal]=i;
+                       IsLocal[i]=kTRUE;
+                       NLocal++;
+                       nnew++;
+                   }
+                   if (nnew==nGlob) break;
+               }
+           }
+       }
+       //
+       // Associate hits to peaks
+       //
+       for (i=0; i<mul; i++) {
+           Float_t dmin=1.E10;
+           Float_t qmax=0;
+           if (IsLocal[i]) continue;
+           for (j=0; j<NLocal; j++) {
+               Int_t il=IndLocal[j];
+               Float_t d=TMath::Sqrt((x[i]-x[il])*(x[i]-x[il])
+                                     +(y[i]-y[il])*(y[i]-y[il]));
+               Float_t ql=q[il];
+               //
+               // Select nearest peak
+               //
+               if (d<dmin) {
+                   dmin=d;
+                   qmax=ql;
+                   AssocPeak[i]=j;
+               } else if (d==dmin) {
+                   //
+                   // If more than one take highest peak
+                   //
+                   if (ql>qmax) {
+                       dmin=d;
+                       qmax=ql;
+                       AssocPeak[i]=j;
+                   }
+               }
+           }
+       }
+       
+
+ //
+ // One cluster for each maximum
+ //
+       for (j=0; j<NLocal; j++) {
+           AliRICHRawCluster cnew;
+           if (fNPeaks == 0) {
+               cnew.fNcluster[0]=-1;
+               cnew.fNcluster[1]=fNRawClusters;
+           } else {
+               cnew.fNcluster[0]=fNPeaks;
+               cnew.fNcluster[1]=0;
+           }
+           cnew.fIndexMap[0]=c->fIndexMap[IndLocal[j]];
+           cnew.fMultiplicity=1;
+           for (i=0; i<mul; i++) {
+               if (IsLocal[i]) continue;
+               if (AssocPeak[i]==j) {
+                   cnew.fIndexMap[cnew.fMultiplicity]=c->fIndexMap[i];
+                   cnew.fMultiplicity++;
+               }
+           }
+           FillCluster(&cnew);
+           cnew.fClusterType=cnew.PhysicsContribution();
+           AddRawCluster(cnew);
+           fNPeaks++;
+       }
+    }
+}
+
+
+void  AliRICHClusterFinder::FillCluster(AliRICHRawCluster* c, Int_t flag) 
+{
+//
+//  Completes cluster information starting from list of digits
+//
+    AliRICHDigit* dig;
+    Float_t x, y;
+    Int_t  ix, iy;
+    Float_t frac=0;
+    
+    c->fPeakSignal=0;
+    if (flag) {
+       c->fX=0;
+       c->fY=0;
+       c->fQ=0;
+    }
+    //c->fQ=0;
+
+    for (Int_t i=0; i<c->fMultiplicity; i++)
+    {
+       dig= (AliRICHDigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
+       ix=dig->fPadX+c->fOffsetMap[i];
+       iy=dig->fPadY;
+       Int_t q=dig->fSignal;
+       if (dig->fPhysics >= dig->fSignal) {
+         c->fPhysicsMap[i]=2;
+       } else if (dig->fPhysics == 0) {
+         c->fPhysicsMap[i]=0;
+       } else  c->fPhysicsMap[i]=1;
+//
+//
+// peak signal and track list
+       if (flag) {
+          if (q>c->fPeakSignal) {
+             c->fPeakSignal=q;
+/*
+           c->fTracks[0]=dig->fTracks[0];
+           c->fTracks[1]=dig->fTracks[1];
+           c->fTracks[2]=dig->fTracks[2];
+*/
+             //c->fTracks[0]=dig->fTrack;
+           c->fTracks[0]=dig->fHit;
+           c->fTracks[1]=dig->fTracks[0];
+           c->fTracks[2]=dig->fTracks[1];
+          }
+       } else {
+          if (c->fContMap[i] > frac) {
+              frac=c->fContMap[i];
+             c->fPeakSignal=q;
+/*
+           c->fTracks[0]=dig->fTracks[0];
+           c->fTracks[1]=dig->fTracks[1];
+           c->fTracks[2]=dig->fTracks[2];
+*/
+             //c->fTracks[0]=dig->fTrack;
+           c->fTracks[0]=dig->fHit;
+           c->fTracks[1]=dig->fTracks[0];
+           c->fTracks[2]=dig->fTracks[1];
+          }
+       }
+//
+       if (flag) {
+           fSegmentation->GetPadCxy(ix, iy, x, y);
+           c->fX += q*x;
+           c->fY += q*y;
+           c->fQ += q;
+       }
+
+    } // loop over digits
+
+ if (flag) {
+     
+     c->fX/=c->fQ;
+     c->fX=fSegmentation->GetAnod(c->fX);
+     c->fY/=c->fQ; 
+//
+//  apply correction to the coordinate along the anode wire
+//
+     x=c->fX;   
+     y=c->fY;
+     fSegmentation->GetPadIxy(x, y, ix, iy);
+     fSegmentation->GetPadCxy(ix, iy, x, y);
+     Int_t isec=fSegmentation->Sector(ix,iy);
+     TF1* CogCorr = fSegmentation->CorrFunc(isec-1);
+     
+     if (CogCorr) {
+        Float_t YonPad=(c->fY-y)/fSegmentation->Dpy(isec);
+        c->fY=c->fY-CogCorr->Eval(YonPad, 0, 0);
+     }
+ }
+}
+
+
+void  AliRICHClusterFinder::FindCluster(Int_t i, Int_t j, AliRICHRawCluster &c){
+//
+//  Find clusters
+//
+//
+//  Add i,j as element of the cluster
+//    
+    
+    Int_t idx = fHitMap->GetHitIndex(i,j);
+    AliRICHDigit* dig = (AliRICHDigit*) fHitMap->GetHit(i,j);
+    Int_t q=dig->fSignal;
+    if (q > TMath::Abs(c.fPeakSignal)) {
+       c.fPeakSignal=q;
+/*
+       c.fTracks[0]=dig->fTracks[0];
+       c.fTracks[1]=dig->fTracks[1];
+       c.fTracks[2]=dig->fTracks[2];
+*/
+       //c.fTracks[0]=dig->fTrack;
+       c.fTracks[0]=dig->fHit;
+       c.fTracks[1]=dig->fTracks[0];
+       c.fTracks[2]=dig->fTracks[1];
+    }
+//
+//  Make sure that list of digits is ordered 
+// 
+    Int_t mu=c.fMultiplicity;
+    c.fIndexMap[mu]=idx;
+
+    if (dig->fPhysics >= dig->fSignal) {
+        c.fPhysicsMap[mu]=2;
+    } else if (dig->fPhysics == 0) {
+        c.fPhysicsMap[mu]=0;
+    } else  c.fPhysicsMap[mu]=1;
+
+    if (mu > 0) {
+       for (Int_t ind=mu-1; ind>=0; ind--) {
+           Int_t ist=(c.fIndexMap)[ind];
+           Int_t ql=((AliRICHDigit*)fDigits
+                     ->UncheckedAt(ist))->fSignal;
+           if (q>ql) {
+               c.fIndexMap[ind]=idx;
+               c.fIndexMap[ind+1]=ist;
+           } else {
+               break;
+           }
+       }
+    }
+    
+    c.fMultiplicity++;
+    
+    if (c.fMultiplicity >= 50 ) {
+       printf("FindCluster - multiplicity >50  %d \n",c.fMultiplicity);
+       c.fMultiplicity=49;
+    }
+
+// Prepare center of gravity calculation
+    Float_t x, y;
+    fSegmentation->GetPadCxy(i, j, x, y);
+    c.fX += q*x;
+    c.fY += q*y;
+    c.fQ += q;
+// Flag hit as taken  
+    fHitMap->FlagHit(i,j);
+//
+//  Now look recursively for all neighbours
+//  
+    Int_t nn;
+    Int_t Xlist[kMaxNeighbours], Ylist[kMaxNeighbours];
+    fSegmentation->Neighbours(i,j,&nn,Xlist,Ylist);
+    for (Int_t in=0; in<nn; in++) {
+       Int_t ix=Xlist[in];
+       Int_t iy=Ylist[in];
+       if (fHitMap->TestHit(ix,iy)==unused) FindCluster(ix, iy, c);
+    }
+}
+
+//_____________________________________________________________________________
+
+void AliRICHClusterFinder::FindRawClusters()
+{
+  //
+  // simple RICH cluster finder from digits -- finds neighbours and 
+  // fill the tree with raw clusters
+  //
+    if (!fNdigits) return;
+
+    fHitMap = new AliRICHHitMapA1(fSegmentation, fDigits);
+
+    AliRICHDigit *dig;
+
+    //printf ("Now I'm here");
+
+    Int_t ndig;
+    Int_t nskip=0;
+    Int_t ncls=0;
+    fHitMap->FillHits();
+    for (ndig=0; ndig<fNdigits; ndig++) {
+       dig = (AliRICHDigit*)fDigits->UncheckedAt(ndig);
+       Int_t i=dig->fPadX;
+       Int_t j=dig->fPadY;
+       if (fHitMap->TestHit(i,j)==used ||fHitMap->TestHit(i,j)==empty) {
+           nskip++;
+           continue;
+       }
+       AliRICHRawCluster c;
+       c.fMultiplicity=0;
+       c.fPeakSignal=dig->fSignal;
+/*
+       c.fTracks[0]=dig->fTracks[0];
+       c.fTracks[1]=dig->fTracks[1];
+       c.fTracks[2]=dig->fTracks[2];
+*/
+       //c.fTracks[0]=dig->fTrack;
+       c.fTracks[0]=dig->fHit;
+       c.fTracks[1]=dig->fTracks[0];
+       c.fTracks[2]=dig->fTracks[1];
+        // tag the beginning of cluster list in a raw cluster
+        c.fNcluster[0]=-1;
+       FindCluster(i,j, c);
+       // center of gravity
+       c.fX /= c.fQ;
+       c.fX=fSegmentation->GetAnod(c.fX);
+       c.fY /= c.fQ;
+//
+//  apply correction to the coordinate along the anode wire
+//
+       Int_t ix,iy;
+       Float_t x=c.fX;   
+       Float_t y=c.fY;
+       fSegmentation->GetPadIxy(x, y, ix, iy);
+       fSegmentation->GetPadCxy(ix, iy, x, y);
+       Int_t isec=fSegmentation->Sector(ix,iy);
+       TF1* CogCorr=fSegmentation->CorrFunc(isec-1);
+       if (CogCorr) {
+           Float_t YonPad=(c.fY-y)/fSegmentation->Dpy(isec);
+           c.fY=c.fY-CogCorr->Eval(YonPad,0,0);
+       }
+
+//
+//      Analyse cluster and decluster if necessary
+//     
+    ncls++;
+    c.fNcluster[1]=fNRawClusters;
+    c.fClusterType=c.PhysicsContribution();
+    Decluster(&c);
+    fNPeaks=0;
+//
+//
+//
+//      reset Cluster object
+       for (int k=0;k<c.fMultiplicity;k++) {
+           c.fIndexMap[k]=0;
+       }
+       c.fMultiplicity=0;
+    } // end loop ndig    
+    delete fHitMap;
+}
+
+void AliRICHClusterFinder::
+CalibrateCOG()
+{
+    Float_t x[5];
+    Float_t y[5];
+    Int_t n, i;
+    TF1 func;
+    if (fSegmentation) {
+       fSegmentation->GiveTestPoints(n, x, y);
+       for (i=0; i<n; i++) {
+           Float_t xtest=x[i];
+           Float_t ytest=y[i];     
+           SinoidalFit(xtest, ytest, func);
+           fSegmentation->SetCorrFunc(i, new TF1(func));
+       }
+    }
+}
+
+
+void AliRICHClusterFinder::
+SinoidalFit(Float_t x, Float_t y, TF1 &func)
+{
+//
+    static Int_t count=0;
+    char canvasname[3];
+    count++;
+    sprintf(canvasname,"c%d",count);
+
+    const Int_t ns=101;
+    Float_t xg[ns], yg[ns], xrg[ns], yrg[ns];
+    Float_t xsig[ns], ysig[ns];
+   
+    AliRICHSegmentation *segmentation=fSegmentation;
+
+    Int_t ix,iy;
+    segmentation->GetPadIxy(x,y,ix,iy);   
+    segmentation->GetPadCxy(ix,iy,x,y);   
+    Int_t isec=segmentation->Sector(ix,iy);
+// Pad Limits    
+    Float_t xmin = x-segmentation->Dpx(isec)/2;
+    Float_t ymin = y-segmentation->Dpy(isec)/2;
+//             
+//      Integration Limits
+    Float_t dxI=fResponse->SigmaIntegration()*fResponse->ChargeSpreadX();
+    Float_t dyI=fResponse->SigmaIntegration()*fResponse->ChargeSpreadY();
+
+//
+//  Scanning
+//
+    Int_t i;
+    Float_t qp;
+//
+//  y-position
+    Float_t yscan=ymin;
+    Float_t dy=segmentation->Dpy(isec)/(ns-1);
+
+    for (i=0; i<ns; i++) {
+//
+//      Pad Loop
+//      
+       Float_t sum=0;
+       Float_t qcheck=0;
+       segmentation->SigGenInit(x, yscan, 0);
+       
+       for (segmentation->FirstPad(x, yscan, dxI, dyI); 
+            segmentation->MorePads(); 
+            segmentation->NextPad()) 
+       {
+           qp=fResponse->IntXY(segmentation);
+           qp=TMath::Abs(qp);
+//
+//
+           if (qp > 1.e-4) {
+               qcheck+=qp;
+               Int_t ixs=segmentation->Ix();
+               Int_t iys=segmentation->Iy();
+               Float_t xs,ys;
+               segmentation->GetPadCxy(ixs,iys,xs,ys);
+               sum+=qp*ys;
+           }
+       } // Pad loop
+       Float_t ycog=sum/qcheck;
+       yg[i]=(yscan-y)/segmentation->Dpy(isec);
+       yrg[i]=(ycog-y)/segmentation->Dpy(isec);
+       ysig[i]=ycog-yscan;
+       yscan+=dy;
+    } // scan loop
+//
+//  x-position
+    Float_t xscan=xmin;
+    Float_t dx=segmentation->Dpx(isec)/(ns-1);
+
+    for (i=0; i<ns; i++) {
+//
+//      Pad Loop
+//      
+       Float_t sum=0;
+       Float_t qcheck=0;
+       segmentation->SigGenInit(xscan, y, 0);
+       
+       for (segmentation->FirstPad(xscan, y, dxI, dyI); 
+            segmentation->MorePads(); 
+            segmentation->NextPad()) 
+       {
+           qp=fResponse->IntXY(segmentation);
+           qp=TMath::Abs(qp);
+//
+//
+           if (qp > 1.e-2) {
+               qcheck+=qp;
+               Int_t ixs=segmentation->Ix();
+               Int_t iys=segmentation->Iy();
+               Float_t xs,ys;
+               segmentation->GetPadCxy(ixs,iys,xs,ys);
+               sum+=qp*xs;
+           }
+       } // Pad loop
+       Float_t xcog=sum/qcheck;
+       xcog=segmentation->GetAnod(xcog);
+       
+       xg[i]=(xscan-x)/segmentation->Dpx(isec);
+       xrg[i]=(xcog-x)/segmentation->Dpx(isec);
+       xsig[i]=xcog-xscan;
+       xscan+=dx;
+    }
+//
+// Creates a Root function based on function sinoid above
+// and perform the fit
+//
+    //    TGraph *graphx = new TGraph(ns,xg ,xsig);
+    //    TGraph *graphxr= new TGraph(ns,xrg,xsig);   
+    //    TGraph *graphy = new TGraph(ns,yg ,ysig);
+    TGraph *graphyr= new TGraph(ns,yrg,ysig);
+
+    Double_t sinoid(Double_t *x, Double_t *par);
+    new TF1("sinoidf",sinoid,0.5,0.5,5);
+    graphyr->Fit("sinoidf","Q");
+    func = *((TF1*)((graphyr->GetListOfFunctions())->At(0)));
+/*
+    
+    TCanvas *c1=new TCanvas(canvasname,canvasname,400,10,600,700);
+    TPad* pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
+    TPad* pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
+    TPad* pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
+    TPad* pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
+    pad11->SetFillColor(11);
+    pad12->SetFillColor(11);
+    pad13->SetFillColor(11);
+    pad14->SetFillColor(11);
+    pad11->Draw();
+    pad12->Draw();
+    pad13->Draw();
+    pad14->Draw();
+    
+//
+    pad11->cd();
+    graphx->SetFillColor(42);
+    graphx->SetMarkerColor(4);
+    graphx->SetMarkerStyle(21);
+    graphx->Draw("AC");
+    graphx->GetHistogram()->SetXTitle("x on pad");
+    graphx->GetHistogram()->SetYTitle("xcog-x");   
+
+
+    pad12->cd();
+    graphxr->SetFillColor(42);
+    graphxr->SetMarkerColor(4);
+    graphxr->SetMarkerStyle(21);
+    graphxr->Draw("AP");
+    graphxr->GetHistogram()->SetXTitle("xcog on pad");
+    graphxr->GetHistogram()->SetYTitle("xcog-x");   
+    
+
+    pad13->cd();
+    graphy->SetFillColor(42);
+    graphy->SetMarkerColor(4);
+    graphy->SetMarkerStyle(21);
+    graphy->Draw("AF");
+    graphy->GetHistogram()->SetXTitle("y on pad");
+    graphy->GetHistogram()->SetYTitle("ycog-y");   
+    
+
+
+    pad14->cd();
+    graphyr->SetFillColor(42);
+    graphyr->SetMarkerColor(4);
+    graphyr->SetMarkerStyle(21);
+    graphyr->Draw("AF");
+    graphyr->GetHistogram()->SetXTitle("ycog on pad");
+    graphyr->GetHistogram()->SetYTitle("ycog-y");   
+    
+    c1->Update();
+*/
+}
+
+Double_t sinoid(Double_t *x, Double_t *par)
+{
+    Double_t arg = -2*TMath::Pi()*x[0];
+    Double_t fitval= par[0]*TMath::Sin(arg)+
+       par[1]*TMath::Sin(2*arg)+
+       par[2]*TMath::Sin(3*arg)+
+       par[3]*TMath::Sin(4*arg)+
+       par[4]*TMath::Sin(5*arg);
+    return fitval;
+ }
+
+
+Double_t DoubleGauss(Double_t *x, Double_t *par)
+{
+    Double_t arg1 = (x[0]-par[1])/0.18;
+    Double_t arg2 = (x[0]-par[3])/0.18;
+    Double_t fitval= par[0]*TMath::Exp(-arg1*arg1/2)
+       +par[2]*TMath::Exp(-arg2*arg2/2);
+    return fitval;
+ }
+
+Float_t DiscrCharge(Int_t i,Double_t *par) 
+{
+// par[0]    x-position of first  cluster
+// par[1]    y-position of first  cluster
+// par[2]    x-position of second cluster
+// par[3]    y-position of second cluster
+// par[4]    charge fraction of first  cluster
+// 1-par[4]  charge fraction of second cluster
+
+    static Float_t qtot;
+    if (gFirst) {
+       qtot=0;
+       for (Int_t jbin=0; jbin<gNbins; jbin++) {
+           qtot+=gCharge[jbin];
+       }
+       gFirst=0;
+       //printf("\n sum of charge from DiscrCharge %f\n", qtot);
+       gChargeTot=Int_t(qtot);
+       
+    }
+    gSegmentation->SetPad(gix[i], giy[i]);
+//  First Cluster
+    gSegmentation->SetHit(par[0],par[1]);
+    Float_t q1=gResponse->IntXY(gSegmentation);
+    
+//  Second Cluster
+    gSegmentation->SetHit(par[2],par[3]);
+    Float_t q2=gResponse->IntXY(gSegmentation);
+    
+    Float_t value = qtot*(par[4]*q1+(1.-par[4])*q2);
+    return value;
+}
+
+//
+// Minimisation function
+void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
+{
+    Int_t i;
+    Float_t delta;
+    Float_t chisq=0;
+    Float_t qcont=0;
+    Float_t qtot=0;
+    
+    for (i=0; i<gNbins; i++) {
+       Float_t q0=gCharge[i];
+       Float_t q1=DiscrCharge(i,par);
+       delta=(q0-q1)/TMath::Sqrt(q0);
+       chisq+=delta*delta;
+       qcont+=q1;
+       qtot+=q0;
+    }
+    chisq=chisq+=(qtot-qcont)*(qtot-qcont)*0.5;
+    f=chisq;
+}
+
+
+
+
+
+
+
+
+
+
+
diff --git a/RICH/AliRICHClusterFinder.h b/RICH/AliRICHClusterFinder.h
new file mode 100644 (file)
index 0000000..1e9afff
--- /dev/null
@@ -0,0 +1,90 @@
+#ifndef AliRICHClusterFinder_H
+#define AliRICHClusterFinder_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+
+////////////////////////////////////////////////
+//  RICH Cluster Finder Class                 //
+////////////////////////////////////////////////
+#include "AliRICHHitMap.h"
+#include "TF1.h"
+class AliRICHClusterFinder :
+ public TObject
+{
+public:
+    TClonesArray*           fDigits;
+    Int_t                   fNdigits;
+protected:
+    AliRICHSegmentation*    fSegmentation;
+    AliRICHResponse*        fResponse;
+    TClonesArray*           fRawClusters;
+    Int_t                   fChamber;
+    Int_t                   fNRawClusters;
+    AliRICHHitMapA1*        fHitMap;
+    TF1*                    fCogCorr;
+    Int_t                   fNperMax;
+    Int_t                   fDeclusterFlag;
+    Int_t                   fClusterSize;
+    Int_t                   fNPeaks; 
+ public:
+    AliRICHClusterFinder
+       (AliRICHSegmentation *segmentation,
+        AliRICHResponse *response, TClonesArray *digits, Int_t chamber);
+    AliRICHClusterFinder();
+    ~AliRICHClusterFinder(){delete fRawClusters;}
+    virtual void SetSegmentation(
+       AliRICHSegmentation *segmentation){
+       fSegmentation=segmentation;
+    }
+    virtual void SetResponse(AliRICHResponse *response) {
+       fResponse=response;
+    }
+
+    virtual void SetDigits(TClonesArray *RICHdigits) {
+       fDigits=RICHdigits;
+       fNdigits = fDigits->GetEntriesFast();
+    }
+    
+    virtual void SetChamber(Int_t ich){
+       fChamber=ich;
+    }
+    
+    virtual void AddRawCluster(const AliRICHRawCluster);
+    // Search for raw clusters
+    virtual void FindRawClusters();
+    virtual void  FindCluster(Int_t i, Int_t j, AliRICHRawCluster &c);
+    // Decluster
+    virtual void Decluster(AliRICHRawCluster *cluster);
+    // Set max. Number of pads per local cluster
+    virtual void SetNperMax(Int_t npermax=5) {fNperMax = npermax;}
+    // Decluster ?
+    virtual void SetDeclusterFlag(Int_t flag=1) {fDeclusterFlag =flag;}
+    // Set max. cluster size ; bigger clusters will be rejected
+    virtual void SetClusterSize(Int_t clsize=5) {fClusterSize = clsize;}
+    // Self Calibration of COG 
+    virtual void CalibrateCOG();
+    virtual void SinoidalFit(Float_t x, Float_t y, TF1 &func);
+    //
+    virtual void CorrectCOG(){;}
+    
+    //
+    virtual Bool_t Centered(AliRICHRawCluster *cluster);
+    virtual void   SplitByLocalMaxima(AliRICHRawCluster *cluster);
+    virtual void   FillCluster(AliRICHRawCluster *cluster, Int_t);
+    virtual void   FillCluster(AliRICHRawCluster *cluster) {
+       FillCluster(cluster,1);}
+    TClonesArray* RawClusters(){return fRawClusters;}
+    ClassDef(AliRICHClusterFinder,1) //Class for clustering and reconstruction of space points
+};
+#endif
+
+
+
+
+
+
+