Algorithm of cluster deconvolution in. It works, for the moment, up to cluster size = 2

diff --git a/RICH/AliRICH.h b/RICH/AliRICH.h
index 07b4203d895ee0198f86b3b245c02f7b5d2be83e..66b4161495fa05a41a39e7c995254d9c8296e667 100644 (file)
--- a/RICH/AliRICH.h
+++ b/RICH/AliRICH.h
@@ -227,6 +227,7 @@ public:
          void       Print(Option_t *option="")const;                                                       //virtual
   inline void       AddDigit(AliRICHdigit *pDig);                                                          //
   inline void       CoG();                                                                                 // 
+         void       Reset() {fSize=fQdc=fStatus=fChamber=fDimXY=kBad;fX=fY=kBad;delete fDigits;fDigits=0;} //                                                                               //
 protected:
   Int_t         fSize;        //how many digits belong to this cluster    
   Int_t         fDimXY;       //100*xdim+ydim box containing the cluster

diff --git a/RICH/AliRICHClusterFinder.cxx b/RICH/AliRICHClusterFinder.cxx
index 0bdfa1db3702eeb1e14a02816588be0e4f83633d..17657ea724923b22aa3ba18a2110818c06d51471 100644 (file)
--- a/RICH/AliRICHClusterFinder.cxx
+++ b/RICH/AliRICHClusterFinder.cxx
@@ -15,12 +15,13 @@
 
 
 #include "AliRICHClusterFinder.h"
-#include "AliRICH.h"
 #include "AliRICHMap.h"
-#include "AliRICHParam.h"
+#include <TMinuit.h>
+#include <TVector3.h>
 #include <AliLoader.h>
 #include <AliRun.h>
 
+void RICHMinMathieson(Int_t &npar, Double_t *gin, Double_t &chi2, Double_t *par, Int_t iflag);
 
 ClassImp(AliRICHClusterFinder)
 //__________________________________________________________________________________________________
@@ -34,19 +35,19 @@ AliRICHClusterFinder::AliRICHClusterFinder(AliRICH *pRICH)
   
 }//main ctor
 //__________________________________________________________________________________________________
-void AliRICHClusterFinder::FindLocalMaxima(AliRICHcluster &rawCluster)
+void AliRICHClusterFinder::FindLocalMaxima()
 {// Split the cluster according to the number of maxima inside
-  Info("SplitbyLocalMaxima","Start.");
+  Info("FindLocalMaxima","Start.");
   Int_t Nlocal = 0;
   Int_t localX[100],localY[100];
-  for(Int_t iDig1=0;iDig1<rawCluster.Size();iDig1++) {
+  for(Int_t iDig1=0;iDig1<fCurrentCluster.Size();iDig1++) {
     Int_t iNotMax = 0;
-    AliRICHdigit *pDig1 = (AliRICHdigit *)rawCluster.Digits()->At(iDig1);
+    AliRICHdigit *pDig1 = (AliRICHdigit *)fCurrentCluster.Digits()->At(iDig1);
     Int_t padX1 = pDig1->X();
     Int_t padY1 = pDig1->Y();
     Double_t padQ1 = pDig1->Q();
-    for(Int_t iDig2=0;iDig2<rawCluster.Size();iDig2++) {
-      AliRICHdigit *pDig2 = (AliRICHdigit *)rawCluster.Digits()->At(iDig2);
+    for(Int_t iDig2=0;iDig2<fCurrentCluster.Size();iDig2++) {
+      AliRICHdigit *pDig2 = (AliRICHdigit *)fCurrentCluster.Digits()->At(iDig2);
       Int_t padX2 = pDig2->X();
       Int_t padY2 = pDig2->Y();
       Double_t padQ2 = pDig2->Q();
@@ -105,15 +106,14 @@ void AliRICHClusterFinder::FindRawClusters(Int_t iChamber)
     Int_t i=dig->X();   Int_t j=dig->Y();
     if(fHitMap->TestHit(i,j)==kUsed) continue;
        
-    AliRICHcluster rawCluster;
-    FormRawCluster(i,j,rawCluster);
+    FormRawCluster(i,j);
        
     if(AliRICHParam::IsResolveClusters()) {
-      ResolveCluster(rawCluster); // ResolveCluster serialization will happen inside
+      ResolveCluster(); // ResolveCluster serialization will happen inside
     } else {
-      WriteRawCluster(rawCluster); // simply output of the RawCluster found without deconvolution
+      WriteRawCluster(); // simply output of the RawCluster found without deconvolution
     }    
-    
+    fCurrentCluster.Reset();
   }//digits loop
 
   delete fHitMap;
@@ -121,34 +121,116 @@ void AliRICHClusterFinder::FindRawClusters(Int_t iChamber)
   
 }//FindRawClusters()
 //__________________________________________________________________________________________________
-void  AliRICHClusterFinder::FormRawCluster(Int_t i, Int_t j, AliRICHcluster &rawCluster)
+void  AliRICHClusterFinder::FormRawCluster(Int_t i, Int_t j)
 {// Builder of the final Raw Cluster (before deconvolution)  
   Info("FormRawCluster","Start with digit(%i,%i)",i,j);
   
-  rawCluster.AddDigit((AliRICHdigit*) fHitMap->GetHit(i,j));
+  fCurrentCluster.AddDigit((AliRICHdigit*) fHitMap->GetHit(i,j));
   fHitMap->FlagHit(i,j);// Flag hit as taken  
 
   Int_t listX[4], listY[4];    //  Now look recursively for all neighbours
   for (Int_t iNeighbour=0;iNeighbour<Rich()->Param()->PadNeighbours(i,j,listX,listY);iNeighbour++)
     if(fHitMap->TestHit(listX[iNeighbour],listY[iNeighbour])==kUnused) 
-                      FormRawCluster(listX[iNeighbour],listY[iNeighbour],rawCluster);    
+                      FormRawCluster(listX[iNeighbour],listY[iNeighbour]);    
 }//AddDigit2Cluster()
 //__________________________________________________________________________________________________
-void AliRICHClusterFinder::ResolveCluster(AliRICHcluster &rawCluster)
+void AliRICHClusterFinder::ResolveCluster()
 {// Decluster algorithm
   Info("ResolveCluster","Start.");    
+
+  fCurrentCluster.CoG(); // first initial approxmation of the CoG...to start minimization.
+  fCurrentCluster.Print();
+  switch (fCurrentCluster.Size()) {
   
-  rawCluster.SetStatus(kRaw);// just dummy to compile...
-     
+  case 1:
+    WriteRawCluster(); break;
+  case 2:
+    FitCoG();
+    WriteRawCluster(); break;
+    
+  default:
+    WriteRawCluster(); break;
+  }     
 }//ResolveCluster()
 //__________________________________________________________________________________________________
-void AliRICHClusterFinder::WriteRawCluster(AliRICHcluster &rawCluster)
+void AliRICHClusterFinder::WriteRawCluster()
 {// out the current RawCluster
   Info("WriteRawCluster","Start.");
   
-  rawCluster.CoG();
-  
-  Rich()->AddCluster(rawCluster);
+  Rich()->AddCluster(fCurrentCluster);
   
 }//WriteRawCluster()
 //__________________________________________________________________________________________________
+void AliRICHClusterFinder::FitCoG()
+{// Fit cluster size 2 by single Mathieson
+  Info("FitCoG","Start.");
+  
+  TMinuit *pMinuit = new TMinuit(2);
+  
+  Double_t arglist;
+  Int_t ierflag = 0;
+
+  static Double_t vstart[2];
+  static Double_t lower[2], upper[2];
+  static Double_t step[2]={0.001,0.001};
+
+  TString chname;
+  Int_t ierflg;
+  
+  pMinuit->SetObjectFit((TObject*)this);
+  pMinuit->SetFCN(RICHMinMathieson);
+  pMinuit->mninit(5,10,7);
+
+  vstart[0] = fCurrentCluster.X();
+  vstart[1] = fCurrentCluster.Y();
+
+  lower[0] = vstart[0] - 2*AliRICHParam::PadSizeX();
+  upper[0] = vstart[0] + 2*AliRICHParam::PadSizeX();
+  lower[1] = vstart[1] - 2*AliRICHParam::PadSizeY();
+  upper[1] = vstart[1] + 2*AliRICHParam::PadSizeY();
+
+
+  pMinuit->mnparm(0," x position ",vstart[0],step[0],lower[0],upper[0],ierflag);
+  pMinuit->mnparm(1," y position ",vstart[1],step[1],lower[1],upper[1],ierflag);
+
+  arglist = -1;
+
+  pMinuit->SetPrintLevel(-1);
+  pMinuit->mnexcm("SET NOGR",&arglist, 1, ierflag);
+  pMinuit->mnexcm("SET NOW",&arglist, 1, ierflag);
+  arglist = 1;
+  pMinuit->mnexcm("SET ERR", &arglist, 1,ierflg);
+  arglist = -1;
+  pMinuit->mnexcm("SIMPLEX",&arglist, 0, ierflag);
+  pMinuit->mnexcm("MIGRAD",&arglist, 0, ierflag);
+  pMinuit->mnexcm("EXIT" ,&arglist, 0, ierflag);  
+  Double_t xCoG,yCoG;
+  Double_t eps, b1, b2;
+  pMinuit->mnpout(0,chname, xCoG, eps , b1, b2, ierflg);
+  pMinuit->mnpout(1,chname, yCoG, eps , b1, b2, ierflg);
+  delete pMinuit;
+}
+//__________________________________________________________________________________________________
+void RICHMinMathieson(Int_t &, Double_t *, Double_t &chi2, Double_t *par, Int_t iflag)
+{// Minimization function of Mathieson
+  
+  AliRICHcluster *pRawCluster = ((AliRICHClusterFinder*)gMinuit->GetObjectFit())->GetCurrentCluster();
+
+  TVector3 centroid(par[0],par[1],0);
+  
+  chi2 = 0;
+  Int_t qtot = pRawCluster->Q();
+  for(Int_t i=0;i<pRawCluster->Size();i++) {
+    Int_t    padX = ((AliRICHdigit *)pRawCluster->Digits()->At(i))->X();
+    Int_t    padY = ((AliRICHdigit *)pRawCluster->Digits()->At(i))->Y();
+    Double_t padQ = ((AliRICHdigit *)pRawCluster->Digits()->At(i))->Q();
+     chi2 += TMath::Power((qtot*AliRICHParam::Loc2PadFrac(centroid,padX,padY)-padQ),2)/padQ;
+  }   
+  if(iflag == 3)
+    {
+            cout << " --- end convergence...summary --- " << endl;
+            cout << " x position " << par[0] << endl;
+            cout << " y position " << par[1] << endl;
+            cout << " chi2       " << chi2 << endl;
+    }
+}//RICHMinMathieson()

diff --git a/RICH/AliRICHClusterFinder.h b/RICH/AliRICHClusterFinder.h
index 530174ccfa7ce9396e970e5f9ee3f6ec0acd0445..08d83ea5f62505f20d797bee4333bf06b79adf86 100644 (file)
--- a/RICH/AliRICHClusterFinder.h
+++ b/RICH/AliRICHClusterFinder.h
@@ -6,9 +6,8 @@
 
 #include "TTask.h"
 
-class AliRICH;
+#include "AliRICH.h"
 class AliHitMap;
-class AliRICHcluster;
 
 class AliRICHClusterFinder : public TTask
 {
@@ -17,16 +16,18 @@ public:
   virtual ~AliRICHClusterFinder()                                          {;}
   
   AliRICH *Rich() {return fRICH;}                                             //Pointer to RICH  
-  void     Exec();                                                              //Loop on events and chambers  
-  void     FindRawClusters(Int_t iChamber);                                     //Find raw clusters  
-  void     FormRawCluster(Int_t i, Int_t j, AliRICHcluster &cluster);           //form a raw cluster
-  void     FindLocalMaxima(AliRICHcluster &cluster);                            //Find local maxima in a cluster
-  void     ResolveCluster(AliRICHcluster  &cluster);                            //Try to resolve a cluster with maxima > 2
-  //PH  void     CoG();                                                               //Evaluate the CoG as the best  
-  void     WriteRawCluster(AliRICHcluster &cluster);                            //write in the list of the raw clusters  
-protected:      
+  void     Exec();                                                            //Loop on events and chambers  
+  void     FindRawClusters(Int_t iChamber);                                   //Find raw clusters  
+  void     FormRawCluster(Int_t i, Int_t j);                                  //form a raw cluster
+  void     FindLocalMaxima();                                                 //Find local maxima in a cluster
+  void     ResolveCluster();                                                  //Try to resolve a cluster with maxima > 2
+  void     FitCoG();                                                          //Evaluate the CoG as the best 
+  void     WriteRawCluster();                                                 //write in the list of the raw clusters  
+  AliRICHcluster *GetCurrentCluster() {return &fCurrentCluster;}              //Return pointer to the current cluster
+protected:
   AliRICH                *fRICH;                         //Pointer to RICH
-  AliHitMap              *fHitMap;                       //Hit Map with digit positions  
+  AliHitMap              *fHitMap;                       //Hit Map with digit positions
+  AliRICHcluster         fCurrentCluster;                //Current cluster to examine
   ClassDef(AliRICHClusterFinder,0) //Finds raw clusters, trasfers them to resolved clusters through declustering.
 };
 #endif

diff --git a/RICH/AliRICHParam.cxx b/RICH/AliRICHParam.cxx
index d29e61e52c3ab7c2bf815c1cd5dbddd03ae16ccf..f172b2b24c4991af22cb126fdaa80fcade08f8f6 100644 (file)
--- a/RICH/AliRICHParam.cxx
+++ b/RICH/AliRICHParam.cxx
@@ -15,7 +15,7 @@
 #include "AliRICHParam.h"
 
 Bool_t   AliRICHParam::fgIsWireSag            =kTRUE;
-Bool_t   AliRICHParam::fgIsResolveClusters    =kFALSE;
+Bool_t   AliRICHParam::fgIsResolveClusters    =kTRUE;
 Double_t AliRICHParam::fgAngleRot             =-60;
 Int_t    AliRICHParam::fgHV                   =2150;
 

diff --git a/RICH/AliRICHParam.h b/RICH/AliRICHParam.h
index 665b4436b10644f049b9aa4f434e886acb1d1f04..a976bda1ca2fce990177bd948f4e81925103f583 100644 (file)
--- a/RICH/AliRICHParam.h
+++ b/RICH/AliRICHParam.h
@@ -43,8 +43,8 @@ public:
   static Double_t InnerFreonLength()         {return 133;}   
   static Double_t InnerFreonWidth()          {return 41.3;}   
   static Double_t IonisationPotential()      {return 26.0e-9;}                            
-  static Double_t MathiensonDeltaX()         {return 5*0.18;}    
-  static Double_t MathiensonDeltaY()         {return 5*0.18;}    
+  static Double_t MathiesonDeltaX()          {return 5*0.18;}    
+  static Double_t MathiesonDeltaY()          {return 5*0.18;}    
   static Int_t    MaxQdc()                   {return 4095;}          
   static Double_t QdcSlope(Int_t sec)        {HV(sec);return 27;}
   static Double_t AlphaFeedback(Int_t sec)   {HV(sec);return 0.036;}
@@ -58,7 +58,7 @@ public:
     static void  SetHV(Int_t hv)             {fgHV       =hv;}  
     static void  SetAngleRot(Double_t rot)   {fgAngleRot =rot;}         
 
-  inline static Double_t Mathienson(Double_t lx1,Double_t lx2,Double_t ly1,Double_t ly2);   
+  inline static Double_t Mathieson(Double_t lx1,Double_t lx2,Double_t ly1,Double_t ly2);   
   inline static void    Loc2Area(TVector3 hitX3,Int_t &padxMin,Int_t &padyMin,Int_t &padxMax,Int_t &padyMax);
   inline static Int_t   PadNeighbours(Int_t iPadX,Int_t iPadY,Int_t aListX[4],Int_t aListY[4]);
   inline static Int_t   Loc2Pad(Double_t x,Double_t y,Int_t &padx,Int_t &pady); 
@@ -202,10 +202,10 @@ Double_t AliRICHParam::Loc2PadFrac(TVector3 hitX3,Int_t padx,Int_t pady)
   Double_t normYmin=(hitX3.Y()-padYcenter-PadSizeY()/2)  /AnodeCathodeGap();
   Double_t normYmax=(hitX3.Y()-padYcenter+PadSizeY()/2)  /AnodeCathodeGap();
   
-  return Mathienson(normXmin,normYmin,normXmax,normYmax);
+  return Mathieson(normXmin,normYmin,normXmax,normYmax);
 }//Loc2PadQdc()
 //__________________________________________________________________________________________________
-Double_t AliRICHParam::Mathienson(Double_t xMin,Double_t yMin,Double_t xMax,Double_t yMax)
+Double_t AliRICHParam::Mathieson(Double_t xMin,Double_t yMin,Double_t xMax,Double_t yMax)
 {//see NIM A370(1988)602-603 
   const Double_t SqrtKx3=0.77459667;const Double_t Kx2=0.962;const Double_t Kx4=0.379;
   const Double_t SqrtKy3=0.77459667;const Double_t Ky2=0.962;const Double_t Ky4=0.379;
@@ -221,7 +221,7 @@ void AliRICHParam::Loc2Area(TVector3 hitX3,Int_t &iPadXmin,Int_t &iPadYmin,Int_t
 {//calculates the area of disintegration for a given hit. Area is a rectangulare set pf pads
  //defined by its left-down and right-up coners
   //  hitX3.SetX(Shift2NearestWire(hitX3.X());
-  Loc2Pad(hitX3.X()-MathiensonDeltaX(),hitX3.Y()-MathiensonDeltaY(),iPadXmin,iPadYmin);   
-  Loc2Pad(hitX3.X()+MathiensonDeltaX(),hitX3.Y()+MathiensonDeltaY(),iPadXmax,iPadYmax);     
+  Loc2Pad(hitX3.X()-MathiesonDeltaX(),hitX3.Y()-MathiesonDeltaY(),iPadXmin,iPadYmin);   
+  Loc2Pad(hitX3.X()+MathiesonDeltaX(),hitX3.Y()+MathiesonDeltaY(),iPadXmax,iPadYmax);     
 }//
 #endif //AliRICHParam_h