RICH/AliRICHClusterFinder.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /*
  17   $Log$
  18   Revision 1.11.4.1  2003/07/10 17:09:50  hristov
  19   Merging to v3-10-00
  20
  21   Revision 1.12  2003/06/27 10:03:00  rdm
  22   avoid assignment of TF1 objects. In SinoidalFit() pass TF1 by pointer.
  23
  24   Revision 1.11  2001/10/23 13:03:35  hristov
  25   The access to several data members was changed from public to protected. The digitisation was adapted to the multi-event case (J.Chudoba)
  26
  27   Revision 1.10  2001/02/27 15:21:58  jbarbosa
  28   Transition to SDigits.
  29
  30   Revision 1.9  2001/01/26 20:00:27  hristov
  31   Major upgrade of AliRoot code
  32
  33   Revision 1.8  2000/11/02 09:11:12  jbarbosa
  34   Removed AliRICHRecHit.h from include.
  35
  36   Revision 1.7  2000/10/03 21:44:09  morsch
  37   Use AliSegmentation and AliHit abstract base classes.
  38
  39   Revision 1.6  2000/10/02 21:28:12  fca
  40   Removal of useless dependecies via forward declarations
  41
  42   Revision 1.5  2000/10/02 15:45:58  jbarbosa
  43   Fixed forward declarations.
  44
  45   Revision 1.4  2000/06/12 19:01:29  morsch
  46   Clean-up bug in Centered() corrected.
  47
  48   Revision 1.3  2000/06/12 15:49:44  jbarbosa
  49   Removed verbose output.
  50
  51   Revision 1.2  2000/06/12 15:18:19  jbarbosa
  52   Cleaned up version.
  53
  54   Revision 1.1  2000/04/19 13:01:48  morsch
  55   A cluster finder and hit reconstruction class for RICH (adapted from MUON).
  56   Cluster Finders for MUON and RICH should derive from the same class in the
  57   future (JB, AM).
  58
  59 */
  60
  61
  62 #include "AliRICHClusterFinder.h"
  63 #include "AliRun.h"
  64 #include "AliRICH.h"
  65 #include "AliRICHHit.h"
  66 #include "AliRICHHitMapA1.h"
  67 #include "AliRICHCerenkov.h"
  68 #include "AliRICHSDigit.h"
  69 #include "AliRICHDigit.h"
  70 #include "AliRICHRawCluster.h"
  71
  72 #include <TTree.h>
  73 #include <TCanvas.h>
  74 #include <TH1.h>
  75 #include <TF1.h>
  76 #include <TPad.h>
  77 #include <TGraph.h>
  78 #include <TPostScript.h>
  79 #include <TMinuit.h>
  80
  81 //----------------------------------------------------------
  82 static AliSegmentation* gSegmentation;
  83 static AliRICHResponse*     gResponse;
  84 static Int_t                gix[500];
  85 static Int_t                giy[500];
  86 static Float_t              gCharge[500];
  87 static Int_t                gNbins;
  88 static Int_t                gFirst=kTRUE;
  89 static TMinuit *gMyMinuit ;
  90 void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
  91 static Int_t                gChargeTot;
  92
  93 ClassImp(AliRICHClusterFinder)
  94
  95 AliRICHClusterFinder::AliRICHClusterFinder
  96 (AliSegmentation *segmentation, AliRICHResponse *response,
  97  TClonesArray *digits, Int_t chamber)
  98 {
  99
 100 // Constructor for Cluster Finder object
 101
 102     fSegmentation=segmentation;
 103     fResponse=response;
 104
 105     fDigits=digits;
 106     fNdigits = fDigits->GetEntriesFast();
 107     fChamber=chamber;
 108     fRawClusters=new TClonesArray("AliRICHRawCluster",10000);
 109     fNRawClusters=0;
 110     fCogCorr = 0;
 111     SetNperMax();
 112     SetClusterSize();
 113     SetDeclusterFlag();
 114     fNPeaks=-1;
 115 }
 116
 117 AliRICHClusterFinder::AliRICHClusterFinder()
 118 {
 119
 120 // Default constructor
 121
 122     fSegmentation=0;
 123     fResponse=0;
 124
 125     fDigits=0;
 126     fNdigits = 0;
 127     fChamber=-1;
 128     fRawClusters=new TClonesArray("AliRICHRawCluster",10000);
 129     fNRawClusters=0;
 130     fHitMap = 0;
 131     fCogCorr = 0;
 132     SetNperMax();
 133     SetClusterSize();
 134     SetDeclusterFlag();
 135     fNPeaks=-1;
 136 }
 137
 138 AliRICHClusterFinder::AliRICHClusterFinder(const AliRICHClusterFinder& ClusterFinder)
 139 {
 140 // Copy Constructor
 141 }
 142
 143 AliRICHClusterFinder::~AliRICHClusterFinder()
 144 {
 145
 146 // Destructor
 147
 148 delete fRawClusters;
 149 }
 150
 151 void AliRICHClusterFinder::AddRawCluster(const AliRICHRawCluster c)
 152 {
 153   //
 154   // Add a raw cluster copy to the list
 155   //
 156     AliRICH *pRICH=(AliRICH*)gAlice->GetModule("RICH");
 157     pRICH->AddRawCluster(fChamber,c);
 158     fNRawClusters++;
 159 }
 160
 161
 162
 163 void AliRICHClusterFinder::Decluster(AliRICHRawCluster *cluster)
 164 {
 165
 166 //
 167 // Decluster algorithm
 168
 169     Int_t mul = cluster->fMultiplicity;
 170 //    printf("Decluster - multiplicity   %d \n",mul);
 171
 172     if (mul == 1 || mul ==2) {
 173 //
 174 // Nothing special for 1- and 2-clusters
 175         if (fNPeaks != 0) {
 176             cluster->fNcluster[0]=fNPeaks;
 177             cluster->fNcluster[1]=0;
 178         }
 179         AddRawCluster(*cluster);
 180         fNPeaks++;
 181     } else if (mul ==3) {
 182 //
 183 // 3-cluster, check topology
 184 //      printf("\n 3-cluster, check topology \n");
 185         if (fDeclusterFlag) {
 186             if (Centered(cluster)) {
 187                 // ok, cluster is centered
 188             } else {
 189                 // cluster is not centered, split into 2+1
 190             }
 191         } else {
 192             if (fNPeaks != 0) {
 193                 cluster->fNcluster[0]=fNPeaks;
 194                 cluster->fNcluster[1]=0;
 195             }
 196             AddRawCluster(*cluster);
 197             fNPeaks++;
 198         }
 199     } else {
 200 //
 201 // 4-and more-pad clusters
 202 //
 203         if (mul <= fClusterSize) {
 204             if (fDeclusterFlag) {
 205                 SplitByLocalMaxima(cluster);
 206             } else {
 207                 if (fNPeaks != 0) {
 208                     cluster->fNcluster[0]=fNPeaks;
 209                     cluster->fNcluster[1]=0;
 210                 }
 211                 AddRawCluster(*cluster);
 212                 fNPeaks++;
 213             }
 214         }
 215     } // multiplicity
 216 }
 217
 218
 219 Bool_t AliRICHClusterFinder::Centered(AliRICHRawCluster *cluster)
 220 {
 221
 222 // Is the cluster centered?
 223
 224     AliRICHDigit* dig;
 225     dig= (AliRICHDigit*)fDigits->UncheckedAt(cluster->fIndexMap[0]);
 226     Int_t ix=dig->PadX();
 227     Int_t iy=dig->PadY();
 228     Int_t nn;
 229     Int_t x[kMaxNeighbours], y[kMaxNeighbours], xN[kMaxNeighbours], yN[kMaxNeighbours];
 230
 231     fSegmentation->Neighbours(ix,iy,&nn,x,y);
 232     Int_t nd=0;
 233     for (Int_t i=0; i<nn; i++) {
 234         if (fHitMap->TestHit(x[i],y[i]) == kUsed) {
 235             xN[nd]=x[i];
 236             yN[nd]=y[i];
 237             nd++;
 238
 239             //printf("Getting: %d %d %d\n",i,x[i],y[i]);
 240         }
 241     }
 242     if (nd==2) {
 243 //
 244 // cluster is centered !
 245         if (fNPeaks != 0) {
 246             cluster->fNcluster[0]=fNPeaks;
 247             cluster->fNcluster[1]=0;
 248         }
 249         cluster->fCtype=0;
 250         AddRawCluster(*cluster);
 251         fNPeaks++;
 252         return kTRUE;
 253     } else if (nd ==1) {
 254 //
 255 // Highest signal on an edge, split cluster into 2+1
 256 //
 257 // who is the neighbour ?
 258
 259       //printf("Calling GetIndex with x:%d y:%d\n",xN[0], yN[0]);
 260
 261         Int_t nind=fHitMap->GetHitIndex(xN[0], yN[0]);
 262         Int_t i1= (nind==cluster->fIndexMap[1]) ? 1:2;
 263         Int_t i2= (nind==cluster->fIndexMap[1]) ? 2:1;
 264 //
 265 // 2-cluster
 266         AliRICHRawCluster cnew;
 267         if (fNPeaks == 0) {
 268             cnew.fNcluster[0]=-1;
 269             cnew.fNcluster[1]=fNRawClusters;
 270         } else {
 271             cnew.fNcluster[0]=fNPeaks;
 272             cnew.fNcluster[1]=0;
 273         }
 274         cnew.fMultiplicity=2;
 275         cnew.fIndexMap[0]=cluster->fIndexMap[0];
 276         cnew.fIndexMap[1]=cluster->fIndexMap[i1];
 277         FillCluster(&cnew);
 278         cnew.fClusterType=cnew.PhysicsContribution();
 279         AddRawCluster(cnew);
 280         fNPeaks++;
 281 //
 282 // 1-cluster
 283         cluster->fMultiplicity=1;
 284         cluster->fIndexMap[0]=cluster->fIndexMap[i2];
 285         cluster->fIndexMap[1]=0;
 286         cluster->fIndexMap[2]=0;
 287         FillCluster(cluster);
 288         if (fNPeaks != 0) {
 289             cluster->fNcluster[0]=fNPeaks;
 290             cluster->fNcluster[1]=0;
 291         }
 292         cluster->fClusterType=cluster->PhysicsContribution();
 293         AddRawCluster(*cluster);
 294         fNPeaks++;
 295         return kFALSE;
 296     } else {
 297         printf("\n Completely screwed up %d !! \n",nd);
 298
 299     }
 300
 301         return kFALSE;
 302 }
 303 void AliRICHClusterFinder::SplitByLocalMaxima(AliRICHRawCluster *c)
 304 {
 305
 306 //
 307 // Split the cluster according to the number of maxima inside
 308
 309
 310     AliRICHDigit* dig[100], *digt;
 311     Int_t ix[100], iy[100], q[100];
 312     Float_t x[100], y[100], zdum;
 313     Int_t i; // loops over digits
 314     Int_t j; // loops over local maxima
 315     //    Float_t xPeak[2];
 316     //    Float_t yPeak[2];
 317     //    Int_t threshold=500;
 318     Int_t mul=c->fMultiplicity;
 319 //
 320 //  dump digit information into arrays
 321 //
 322     for (i=0; i<mul; i++)
 323     {
 324         dig[i]= (AliRICHDigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
 325         ix[i]= dig[i]->PadX();
 326         iy[i]= dig[i]->PadY();
 327         q[i] = dig[i]->Signal();
 328         fSegmentation->GetPadC(ix[i], iy[i], x[i], y[i], zdum);
 329     }
 330 //
 331 //  Find local maxima
 332 //
 333     Bool_t isLocal[100];
 334     Int_t nLocal=0;
 335     Int_t associatePeak[100];
 336     Int_t indLocal[100];
 337     Int_t nn;
 338     Int_t xNei[kMaxNeighbours], yNei[kMaxNeighbours];
 339     for (i=0; i<mul; i++) {
 340         fSegmentation->Neighbours(ix[i], iy[i], &nn, xNei, yNei);
 341         isLocal[i]=kTRUE;
 342         for (j=0; j<nn; j++) {
 343             if (fHitMap->TestHit(xNei[j], yNei[j])==kEmpty) continue;
 344             digt=(AliRICHDigit*) fHitMap->GetHit(xNei[j], yNei[j]);
 345             if (digt->Signal() > q[i]) {
 346                 isLocal[i]=kFALSE;
 347                 break;
 348 //
 349 // handle special case of neighbouring pads with equal signal
 350             } else if (digt->Signal() == q[i]) {
 351                 if (nLocal >0) {
 352                     for (Int_t k=0; k<nLocal; k++) {
 353                         if (xNei[j]==ix[indLocal[k]] && yNei[j]==iy[indLocal[k]]){
 354                             isLocal[i]=kFALSE;
 355                         }
 356                     }
 357                 }
 358             }
 359         } // loop over next neighbours
 360         // Maxima should not be on the edge
 361         if (isLocal[i]) {
 362             indLocal[nLocal]=i;
 363             nLocal++;
 364         }
 365     } // loop over all digits
 366 //    printf("Found %d local Maxima",nLocal);
 367 //
 368 // If only one local maximum found but multiplicity is high
 369 // take global maximum from the list of digits.
 370     if (nLocal==1 && mul>5) {
 371         Int_t nnew=0;
 372         for (i=0; i<mul; i++) {
 373             if (!isLocal[i]) {
 374                 indLocal[nLocal]=i;
 375                 isLocal[i]=kTRUE;
 376                 nLocal++;
 377                 nnew++;
 378             }
 379             if (nnew==1) break;
 380         }
 381     }
 382
 383 // If number of local maxima is 2 try to fit a double gaussian
 384     if (nLocal==-100) {
 385 //
 386 //  Initialise global variables for fit
 387         gFirst=1;
 388         gSegmentation=fSegmentation;
 389         gResponse    =fResponse;
 390         gNbins=mul;
 391
 392         for (i=0; i<mul; i++) {
 393             gix[i]=ix[i];
 394             giy[i]=iy[i];
 395             gCharge[i]=Float_t(q[i]);
 396         }
 397 //
 398         if (gFirst) {
 399             gFirst=kFALSE;
 400             gMyMinuit = new TMinuit(5);
 401         }
 402         gMyMinuit->SetFCN(fcn);
 403         gMyMinuit->mninit(5,10,7);
 404         Double_t arglist[20];
 405         Int_t ierflag=0;
 406         arglist[0]=1;
 407 //      gMyMinuit->mnexcm("SET ERR",arglist,1,ierflag);
 408 // Set starting values
 409         static Double_t vstart[5];
 410         vstart[0]=x[indLocal[0]];
 411         vstart[1]=y[indLocal[0]];
 412         vstart[2]=x[indLocal[1]];
 413         vstart[3]=y[indLocal[1]];
 414         vstart[4]=Float_t(q[indLocal[0]])/
 415             Float_t(q[indLocal[0]]+q[indLocal[1]]);
 416 // lower and upper limits
 417         static Double_t lower[5], upper[5];
 418         Int_t isec=fSegmentation->Sector(ix[indLocal[0]], iy[indLocal[0]]);
 419         lower[0]=vstart[0]-fSegmentation->Dpx(isec)/2;
 420         lower[1]=vstart[1]-fSegmentation->Dpy(isec)/2;
 421 //      lower[1]=vstart[1];
 422
 423         upper[0]=lower[0]+fSegmentation->Dpx(isec);
 424         upper[1]=lower[1]+fSegmentation->Dpy(isec);
 425 //      upper[1]=vstart[1];
 426
 427         isec=fSegmentation->Sector(ix[indLocal[1]], iy[indLocal[1]]);
 428         lower[2]=vstart[2]-fSegmentation->Dpx(isec)/2;
 429         lower[3]=vstart[3]-fSegmentation->Dpy(isec)/2;
 430 //      lower[3]=vstart[3];
 431
 432         upper[2]=lower[2]+fSegmentation->Dpx(isec);
 433         upper[3]=lower[3]+fSegmentation->Dpy(isec);
 434 //      upper[3]=vstart[3];
 435
 436         lower[4]=0.;
 437         upper[4]=1.;
 438 // step sizes
 439         static Double_t step[5]={0.005, 0.03, 0.005, 0.03, 0.01};
 440
 441         gMyMinuit->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
 442         gMyMinuit->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
 443         gMyMinuit->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
 444         gMyMinuit->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
 445         gMyMinuit->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
 446 // ready for minimisation
 447         gMyMinuit->SetPrintLevel(-1);
 448         gMyMinuit->mnexcm("SET OUT", arglist, 0, ierflag);
 449         arglist[0]= -1;
 450         arglist[1]= 0;
 451
 452         gMyMinuit->mnexcm("SET NOGR", arglist, 0, ierflag);
 453         gMyMinuit->mnexcm("SCAN", arglist, 0, ierflag);
 454         gMyMinuit->mnexcm("EXIT" , arglist, 0, ierflag);
 455 // Print results
 456 //      Double_t amin,edm,errdef;
 457 //      Int_t nvpar,nparx,icstat;
 458 //      gMyMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
 459 //      gMyMinuit->mnprin(3,amin);
 460 // Get fitted parameters
 461
 462         Double_t xrec[2], yrec[2], qfrac;
 463         TString chname;
 464         Double_t epxz, b1, b2;
 465         Int_t ierflg;
 466         gMyMinuit->mnpout(0, chname, xrec[0], epxz, b1, b2, ierflg);
 467         gMyMinuit->mnpout(1, chname, yrec[0], epxz, b1, b2, ierflg);
 468         gMyMinuit->mnpout(2, chname, xrec[1], epxz, b1, b2, ierflg);
 469         gMyMinuit->mnpout(3, chname, yrec[1], epxz, b1, b2, ierflg);
 470         gMyMinuit->mnpout(4, chname, qfrac,   epxz, b1, b2, ierflg);
 471         //printf("\n %f %f %f %f %f\n", xrec[0], yrec[0], xrec[1], yrec[1],qfrac);
 472 //      delete gMyMinuit;
 473
 474
 475  //
 476  // One cluster for each maximum
 477  //
 478         for (j=0; j<2; j++) {
 479             AliRICHRawCluster cnew;
 480             if (fNPeaks == 0) {
 481                 cnew.fNcluster[0]=-1;
 482                 cnew.fNcluster[1]=fNRawClusters;
 483             } else {
 484                 cnew.fNcluster[0]=fNPeaks;
 485                 cnew.fNcluster[1]=0;
 486             }
 487             cnew.fMultiplicity=0;
 488             cnew.fX=Float_t(xrec[j]);
 489             cnew.fY=Float_t(yrec[j]);
 490             if (j==0) {
 491                 cnew.fQ=Int_t(gChargeTot*qfrac);
 492             } else {
 493                 cnew.fQ=Int_t(gChargeTot*(1-qfrac));
 494             }
 495             gSegmentation->SetHit(xrec[j],yrec[j],0);
 496             for (i=0; i<mul; i++) {
 497                 cnew.fIndexMap[cnew.fMultiplicity]=c->fIndexMap[i];
 498                 gSegmentation->SetPad(gix[i], giy[i]);
 499                 Float_t q1=gResponse->IntXY(gSegmentation);
 500                 cnew.fContMap[cnew.fMultiplicity]=Float_t(q[i])/(q1*cnew.fQ);
 501                 cnew.fMultiplicity++;
 502             }
 503             FillCluster(&cnew,0);
 504             //printf("\n x,y %f %f ", cnew.fX, cnew.fY);
 505             cnew.fClusterType=cnew.PhysicsContribution();
 506             AddRawCluster(cnew);
 507             fNPeaks++;
 508         }
 509     }
 510
 511     Bool_t fitted=kTRUE;
 512
 513     if (nLocal !=-100 || !fitted) {
 514         // Check if enough local clusters have been found,
 515         // if not add global maxima to the list
 516         //
 517         Int_t nPerMax;
 518         if (nLocal!=0) {
 519             nPerMax=mul/nLocal;
 520         } else {
 521             printf("\n Warning, no local maximum found \n");
 522             nPerMax=fNperMax+1;
 523         }
 524
 525         if (nPerMax > fNperMax) {
 526             Int_t nGlob=mul/fNperMax-nLocal+1;
 527             if (nGlob > 0) {
 528                 Int_t nnew=0;
 529                 for (i=0; i<mul; i++) {
 530                     if (!isLocal[i]) {
 531                         indLocal[nLocal]=i;
 532                         isLocal[i]=kTRUE;
 533                         nLocal++;
 534                         nnew++;
 535                     }
 536                     if (nnew==nGlob) break;
 537                 }
 538             }
 539         }
 540         //
 541         // Associate hits to peaks
 542         //
 543         for (i=0; i<mul; i++) {
 544             Float_t dmin=1.E10;
 545             Float_t qmax=0;
 546             if (isLocal[i]) continue;
 547             for (j=0; j<nLocal; j++) {
 548                 Int_t il=indLocal[j];
 549                 Float_t d=TMath::Sqrt((x[i]-x[il])*(x[i]-x[il])
 550                                       +(y[i]-y[il])*(y[i]-y[il]));
 551                 Float_t ql=q[il];
 552                 //
 553                 // Select nearest peak
 554                 //
 555                 if (d<dmin) {
 556                     dmin=d;
 557                     qmax=ql;
 558                     associatePeak[i]=j;
 559                 } else if (d==dmin) {
 560                     //
 561                     // If more than one take highest peak
 562                     //
 563                     if (ql>qmax) {
 564                         dmin=d;
 565                         qmax=ql;
 566                         associatePeak[i]=j;
 567                     }
 568                 }
 569             }
 570         }
 571
 572
 573  //
 574  // One cluster for each maximum
 575  //
 576         for (j=0; j<nLocal; j++) {
 577             AliRICHRawCluster cnew;
 578             if (fNPeaks == 0) {
 579                 cnew.fNcluster[0]=-1;
 580                 cnew.fNcluster[1]=fNRawClusters;
 581             } else {
 582                 cnew.fNcluster[0]=fNPeaks;
 583                 cnew.fNcluster[1]=0;
 584             }
 585             cnew.fIndexMap[0]=c->fIndexMap[indLocal[j]];
 586             cnew.fMultiplicity=1;
 587             for (i=0; i<mul; i++) {
 588                 if (isLocal[i]) continue;
 589                 if (associatePeak[i]==j) {
 590                     cnew.fIndexMap[cnew.fMultiplicity]=c->fIndexMap[i];
 591                     cnew.fMultiplicity++;
 592                 }
 593             }
 594             FillCluster(&cnew);
 595             cnew.fClusterType=cnew.PhysicsContribution();
 596             AddRawCluster(cnew);
 597             fNPeaks++;
 598         }
 599     }
 600 }
 601
 602
 603 void  AliRICHClusterFinder::FillCluster(AliRICHRawCluster* c, Int_t flag)
 604 {
 605 //
 606 //  Completes cluster information starting from list of digits
 607 //
 608     AliRICHDigit* dig;
 609     Float_t x, y, z;
 610     Int_t  ix, iy;
 611     Float_t frac=0;
 612
 613     c->fPeakSignal=0;
 614     if (flag) {
 615         c->fX=0;
 616         c->fY=0;
 617         c->fQ=0;
 618     }
 619     //c->fQ=0;
 620
 621
 622     for (Int_t i=0; i<c->fMultiplicity; i++)
 623     {
 624         dig= (AliRICHDigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
 625         ix=dig->PadX()+c->fOffsetMap[i];
 626         iy=dig->PadY();
 627         Int_t q=dig->Signal();
 628         if (dig->Physics() >= dig->Signal()) {
 629           c->fPhysicsMap[i]=2;
 630         } else if (dig->Physics() == 0) {
 631           c->fPhysicsMap[i]=0;
 632         } else  c->fPhysicsMap[i]=1;
 633 //
 634 //
 635 // peak signal and track list
 636         if (flag) {
 637            if (q>c->fPeakSignal) {
 638               c->fPeakSignal=q;
 639 /*
 640             c->fTracks[0]=dig->Track(0);
 641             c->fTracks[1]=dig->Track(1);
 642             c->fTracks[2]=dig->Track(2);
 643 */
 644               //c->fTracks[0]=dig->fTrack;
 645             c->fTracks[0]=dig->Hit();
 646             c->fTracks[1]=dig->Track(0);
 647             c->fTracks[2]=dig->Track(1);
 648            }
 649         } else {
 650            if (c->fContMap[i] > frac) {
 651               frac=c->fContMap[i];
 652               c->fPeakSignal=q;
 653 /*
 654             c->fTracks[0]=dig->Track(0);
 655             c->fTracks[1]=dig->Track(1);
 656             c->fTracks[2]=dig->Track(2);
 657 */
 658               //c->fTracks[0]=dig->fTrack;
 659             c->fTracks[0]=dig->Hit();
 660             c->fTracks[1]=dig->Track(0);
 661             c->fTracks[2]=dig->Track(1);
 662            }
 663         }
 664 //
 665         if (flag) {
 666             fSegmentation->GetPadC(ix, iy, x, y, z);
 667             c->fX += q*x;
 668             c->fY += q*y;
 669             c->fQ += q;
 670         }
 671
 672     } // loop over digits
 673
 674  if (flag) {
 675
 676      c->fX/=c->fQ;
 677      c->fX=fSegmentation->GetAnod(c->fX);
 678      c->fY/=c->fQ;
 679 //
 680 //  apply correction to the coordinate along the anode wire
 681 //
 682      x=c->fX;
 683      y=c->fY;
 684      fSegmentation->GetPadI(x, y, 0, ix, iy);
 685      fSegmentation->GetPadC(ix, iy, x, y, z);
 686      Int_t isec=fSegmentation->Sector(ix,iy);
 687      TF1* cogCorr = fSegmentation->CorrFunc(isec-1);
 688
 689      if (cogCorr) {
 690          Float_t yOnPad=(c->fY-y)/fSegmentation->Dpy(isec);
 691          c->fY=c->fY-cogCorr->Eval(yOnPad, 0, 0);
 692      }
 693  }
 694 }
 695
 696
 697 void  AliRICHClusterFinder::FindCluster(Int_t i, Int_t j, AliRICHRawCluster &c){
 698 //
 699 //  Find clusters
 700 //
 701 //
 702 //  Add i,j as element of the cluster
 703 //
 704
 705     Int_t idx = fHitMap->GetHitIndex(i,j);
 706     AliRICHDigit* dig = (AliRICHDigit*) fHitMap->GetHit(i,j);
 707     Int_t q=dig->Signal();
 708     if (q > TMath::Abs(c.fPeakSignal)) {
 709         c.fPeakSignal=q;
 710 /*
 711         c.fTracks[0]=dig->fTracks[0];
 712         c.fTracks[1]=dig->fTracks[1];
 713         c.fTracks[2]=dig->fTracks[2];
 714 */
 715         //c.fTracks[0]=dig->fTrack;
 716         c.fTracks[0]=dig->Hit();
 717         c.fTracks[1]=dig->Track(0);
 718         c.fTracks[2]=dig->Track(1);
 719     }
 720 //
 721 //  Make sure that list of digits is ordered
 722 //
 723     Int_t mu=c.fMultiplicity;
 724     c.fIndexMap[mu]=idx;
 725
 726     if (dig->Physics() >= dig->Signal()) {
 727         c.fPhysicsMap[mu]=2;
 728     } else if (dig->Physics() == 0) {
 729         c.fPhysicsMap[mu]=0;
 730     } else  c.fPhysicsMap[mu]=1;
 731
 732     if (mu > 0) {
 733         for (Int_t ind=mu-1; ind>=0; ind--) {
 734             Int_t ist=(c.fIndexMap)[ind];
 735             Int_t ql=((AliRICHDigit*)fDigits
 736                       ->UncheckedAt(ist))->Signal();
 737             if (q>ql) {
 738                 c.fIndexMap[ind]=idx;
 739                 c.fIndexMap[ind+1]=ist;
 740             } else {
 741                 break;
 742             }
 743         }
 744     }
 745
 746     c.fMultiplicity++;
 747
 748     if (c.fMultiplicity >= 50 ) {
 749         printf("FindCluster - multiplicity >50  %d \n",c.fMultiplicity);
 750         c.fMultiplicity=49;
 751     }
 752
 753 // Prepare center of gravity calculation
 754     Float_t x, y, z;
 755     fSegmentation->GetPadC(i, j, x, y, z);
 756     c.fX += q*x;
 757     c.fY += q*y;
 758     c.fQ += q;
 759 // Flag hit as taken
 760     fHitMap->FlagHit(i,j);
 761 //
 762 //  Now look recursively for all neighbours
 763 //
 764     Int_t nn;
 765     Int_t xList[kMaxNeighbours], yList[kMaxNeighbours];
 766     fSegmentation->Neighbours(i,j,&nn,xList,yList);
 767     for (Int_t in=0; in<nn; in++) {
 768         Int_t ix=xList[in];
 769         Int_t iy=yList[in];
 770         if (fHitMap->TestHit(ix,iy)==kUnused) FindCluster(ix, iy, c);
 771     }
 772 }
 773
 774 //_____________________________________________________________________________
 775
 776 void AliRICHClusterFinder::FindRawClusters()
 777 {
 778   //
 779   // simple RICH cluster finder from digits -- finds neighbours and
 780   // fill the tree with raw clusters
 781   //
 782     if (!fNdigits) return;
 783
 784     fHitMap = new AliRICHHitMapA1(fSegmentation, fDigits);
 785
 786     AliRICHDigit *dig;
 787
 788     //printf ("Now I'm here");
 789
 790     Int_t ndig;
 791     Int_t nskip=0;
 792     Int_t ncls=0;
 793     fHitMap->FillHits();
 794     for (ndig=0; ndig<fNdigits; ndig++) {
 795         dig = (AliRICHDigit*)fDigits->UncheckedAt(ndig);
 796         Int_t i=dig->PadX();
 797         Int_t j=dig->PadY();
 798         if (fHitMap->TestHit(i,j)==kUsed ||fHitMap->TestHit(i,j)==kEmpty) {
 799             nskip++;
 800             continue;
 801         }
 802         AliRICHRawCluster c;
 803         c.fMultiplicity=0;
 804         c.fPeakSignal=dig->Signal();
 805 /*
 806         c.fTracks[0]=dig->fTracks[0];
 807         c.fTracks[1]=dig->fTracks[1];
 808         c.fTracks[2]=dig->fTracks[2];
 809 */
 810         //c.fTracks[0]=dig->fTrack;
 811         c.fTracks[0]=dig->Hit();
 812         c.fTracks[1]=dig->Track(0);
 813         c.fTracks[2]=dig->Track(1);
 814         // tag the beginning of cluster list in a raw cluster
 815         c.fNcluster[0]=-1;
 816         FindCluster(i,j, c);
 817         // center of gravity
 818         c.fX /= c.fQ;
 819         c.fX=fSegmentation->GetAnod(c.fX);
 820         c.fY /= c.fQ;
 821 //
 822 //  apply correction to the coordinate along the anode wire
 823 //
 824         Int_t ix,iy;
 825         Float_t x=c.fX;
 826         Float_t y=c.fY;
 827         Float_t z;
 828
 829         fSegmentation->GetPadI(x, y, 0, ix, iy);
 830         fSegmentation->GetPadC(ix, iy, x, y, z);
 831         Int_t isec=fSegmentation->Sector(ix,iy);
 832         TF1* cogCorr=fSegmentation->CorrFunc(isec-1);
 833         if (cogCorr) {
 834             Float_t yOnPad=(c.fY-y)/fSegmentation->Dpy(isec);
 835             c.fY=c.fY-cogCorr->Eval(yOnPad,0,0);
 836         }
 837
 838 //
 839 //      Analyse cluster and decluster if necessary
 840 //
 841     ncls++;
 842     c.fNcluster[1]=fNRawClusters;
 843     c.fClusterType=c.PhysicsContribution();
 844     Decluster(&c);
 845     fNPeaks=0;
 846 //
 847 //
 848 //
 849 //      reset Cluster object
 850         for (int k=0;k<c.fMultiplicity;k++) {
 851             c.fIndexMap[k]=0;
 852         }
 853         c.fMultiplicity=0;
 854     } // end loop ndig
 855     delete fHitMap;
 856 }
 857
 858 void AliRICHClusterFinder::
 859 CalibrateCOG()
 860 {
 861
 862 // Calibration
 863
 864     Float_t x[5];
 865     Float_t y[5];
 866     Int_t n, i;
 867     if (fSegmentation) {
 868         TF1 *func;
 869         fSegmentation->GiveTestPoints(n, x, y);
 870         for (i=0; i<n; i++) {
 871             func = 0;
 872             Float_t xtest=x[i];
 873             Float_t ytest=y[i];
 874             SinoidalFit(xtest, ytest, func);
 875             if (func) fSegmentation->SetCorrFunc(i, new TF1(*func));
 876         }
 877     }
 878 }
 879
 880
 881 void AliRICHClusterFinder::
 882 SinoidalFit(Float_t x, Float_t y, TF1 *func)
 883 {
 884 // Sinoidal fit
 885
 886
 887     static Int_t count=0;
 888     char canvasname[3];
 889     Float_t z;
 890
 891     count++;
 892     sprintf(canvasname,"c%d",count);
 893
 894     const Int_t kNs=101;
 895     Float_t xg[kNs], yg[kNs], xrg[kNs], yrg[kNs];
 896     Float_t xsig[kNs], ysig[kNs];
 897
 898     AliSegmentation *segmentation=fSegmentation;
 899
 900     Int_t ix,iy;
 901     segmentation->GetPadI(x,y,0,ix,iy);
 902     segmentation->GetPadC(ix,iy,x,y,z);
 903     Int_t isec=segmentation->Sector(ix,iy);
 904 // Pad Limits
 905     Float_t xmin = x-segmentation->Dpx(isec)/2;
 906     Float_t ymin = y-segmentation->Dpy(isec)/2;
 907 //
 908 //      Integration Limits
 909     Float_t dxI=fResponse->SigmaIntegration()*fResponse->ChargeSpreadX();
 910     Float_t dyI=fResponse->SigmaIntegration()*fResponse->ChargeSpreadY();
 911
 912 //
 913 //  Scanning
 914 //
 915     Int_t i;
 916     Float_t qp;
 917 //
 918 //  y-position
 919     Float_t yscan=ymin;
 920     Float_t dy=segmentation->Dpy(isec)/(kNs-1);
 921
 922     for (i=0; i<kNs; i++) {
 923 //
 924 //      Pad Loop
 925 //
 926         Float_t sum=0;
 927         Float_t qcheck=0;
 928         segmentation->SigGenInit(x, yscan, 0);
 929
 930         for (segmentation->FirstPad(x, yscan,0, dxI, dyI);
 931              segmentation->MorePads();
 932              segmentation->NextPad())
 933         {
 934             qp=fResponse->IntXY(segmentation);
 935             qp=TMath::Abs(qp);
 936 //
 937 //
 938             if (qp > 1.e-4) {
 939                 qcheck+=qp;
 940                 Int_t ixs=segmentation->Ix();
 941                 Int_t iys=segmentation->Iy();
 942                 Float_t xs,ys,zs;
 943                 segmentation->GetPadC(ixs,iys,xs,ys,zs);
 944                 sum+=qp*ys;
 945             }
 946         } // Pad loop
 947         Float_t ycog=sum/qcheck;
 948         yg[i]=(yscan-y)/segmentation->Dpy(isec);
 949         yrg[i]=(ycog-y)/segmentation->Dpy(isec);
 950         ysig[i]=ycog-yscan;
 951         yscan+=dy;
 952     } // scan loop
 953 //
 954 //  x-position
 955     Float_t xscan=xmin;
 956     Float_t dx=segmentation->Dpx(isec)/(kNs-1);
 957
 958     for (i=0; i<kNs; i++) {
 959 //
 960 //      Pad Loop
 961 //
 962         Float_t sum=0;
 963         Float_t qcheck=0;
 964         segmentation->SigGenInit(xscan, y, 0);
 965
 966         for (segmentation->FirstPad(xscan, y, 0, dxI, dyI);
 967              segmentation->MorePads();
 968              segmentation->NextPad())
 969         {
 970             qp=fResponse->IntXY(segmentation);
 971             qp=TMath::Abs(qp);
 972 //
 973 //
 974             if (qp > 1.e-2) {
 975                 qcheck+=qp;
 976                 Int_t ixs=segmentation->Ix();
 977                 Int_t iys=segmentation->Iy();
 978                 Float_t xs,ys,zs;
 979                 segmentation->GetPadC(ixs,iys,xs,ys,zs);
 980                 sum+=qp*xs;
 981             }
 982         } // Pad loop
 983         Float_t xcog=sum/qcheck;
 984         xcog=segmentation->GetAnod(xcog);
 985
 986         xg[i]=(xscan-x)/segmentation->Dpx(isec);
 987         xrg[i]=(xcog-x)/segmentation->Dpx(isec);
 988         xsig[i]=xcog-xscan;
 989         xscan+=dx;
 990     }
 991 //
 992 // Creates a Root function based on function sinoid above
 993 // and perform the fit
 994 //
 995     //    TGraph *graphx = new TGraph(kNs,xg ,xsig);
 996     //    TGraph *graphxr= new TGraph(kNs,xrg,xsig);
 997     //    TGraph *graphy = new TGraph(kNs,yg ,ysig);
 998     TGraph *graphyr= new TGraph(kNs,yrg,ysig);
 999
1000     Double_t sinoid(Double_t *x, Double_t *par);
1001     new TF1("sinoidf",sinoid,0.5,0.5,5);
1002     graphyr->Fit("sinoidf","Q");
1003     func = (TF1*)graphyr->GetListOfFunctions()->At(0);
1004 /*
1005
1006     TCanvas *c1=new TCanvas(canvasname,canvasname,400,10,600,700);
1007     TPad* pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
1008     TPad* pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
1009     TPad* pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
1010     TPad* pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
1011     pad11->SetFillColor(11);
1012     pad12->SetFillColor(11);
1013     pad13->SetFillColor(11);
1014     pad14->SetFillColor(11);
1015     pad11->Draw();
1016     pad12->Draw();
1017     pad13->Draw();
1018     pad14->Draw();
1019
1020 //
1021     pad11->cd();
1022     graphx->SetFillColor(42);
1023     graphx->SetMarkerColor(4);
1024     graphx->SetMarkerStyle(21);
1025     graphx->Draw("AC");
1026     graphx->GetHistogram()->SetXTitle("x on pad");
1027     graphx->GetHistogram()->SetYTitle("xcog-x");
1028
1029
1030     pad12->cd();
1031     graphxr->SetFillColor(42);
1032     graphxr->SetMarkerColor(4);
1033     graphxr->SetMarkerStyle(21);
1034     graphxr->Draw("AP");
1035     graphxr->GetHistogram()->SetXTitle("xcog on pad");
1036     graphxr->GetHistogram()->SetYTitle("xcog-x");
1037
1038
1039     pad13->cd();
1040     graphy->SetFillColor(42);
1041     graphy->SetMarkerColor(4);
1042     graphy->SetMarkerStyle(21);
1043     graphy->Draw("AF");
1044     graphy->GetHistogram()->SetXTitle("y on pad");
1045     graphy->GetHistogram()->SetYTitle("ycog-y");
1046
1047
1048
1049     pad14->cd();
1050     graphyr->SetFillColor(42);
1051     graphyr->SetMarkerColor(4);
1052     graphyr->SetMarkerStyle(21);
1053     graphyr->Draw("AF");
1054     graphyr->GetHistogram()->SetXTitle("ycog on pad");
1055     graphyr->GetHistogram()->SetYTitle("ycog-y");
1056
1057     c1->Update();
1058 */
1059 }
1060
1061 Double_t sinoid(Double_t *x, Double_t *par)
1062 {
1063
1064 // Sinoid function
1065
1066     Double_t arg = -2*TMath::Pi()*x[0];
1067     Double_t fitval= par[0]*TMath::Sin(arg)+
1068         par[1]*TMath::Sin(2*arg)+
1069         par[2]*TMath::Sin(3*arg)+
1070         par[3]*TMath::Sin(4*arg)+
1071         par[4]*TMath::Sin(5*arg);
1072     return fitval;
1073  }
1074
1075
1076 Double_t DoubleGauss(Double_t *x, Double_t *par)
1077 {
1078
1079 // Doublr gaussian function
1080
1081     Double_t arg1 = (x[0]-par[1])/0.18;
1082     Double_t arg2 = (x[0]-par[3])/0.18;
1083     Double_t fitval= par[0]*TMath::Exp(-arg1*arg1/2)
1084         +par[2]*TMath::Exp(-arg2*arg2/2);
1085     return fitval;
1086  }
1087
1088 Float_t DiscrCharge(Int_t i,Double_t *par)
1089 {
1090 // par[0]    x-position of first  cluster
1091 // par[1]    y-position of first  cluster
1092 // par[2]    x-position of second cluster
1093 // par[3]    y-position of second cluster
1094 // par[4]    charge fraction of first  cluster
1095 // 1-par[4]  charge fraction of second cluster
1096
1097     static Float_t qtot;
1098     if (gFirst) {
1099         qtot=0;
1100         for (Int_t jbin=0; jbin<gNbins; jbin++) {
1101             qtot+=gCharge[jbin];
1102         }
1103         gFirst=0;
1104         //printf("\n sum of charge from DiscrCharge %f\n", qtot);
1105         gChargeTot=Int_t(qtot);
1106
1107     }
1108     gSegmentation->SetPad(gix[i], giy[i]);
1109 //  First Cluster
1110     gSegmentation->SetHit(par[0],par[1],0);
1111     Float_t q1=gResponse->IntXY(gSegmentation);
1112
1113 //  Second Cluster
1114     gSegmentation->SetHit(par[2],par[3],0);
1115     Float_t q2=gResponse->IntXY(gSegmentation);
1116
1117     Float_t value = qtot*(par[4]*q1+(1.-par[4])*q2);
1118     return value;
1119 }
1120
1121 //
1122 // Minimisation function
1123 void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
1124 {
1125     Int_t i;
1126     Float_t delta;
1127     Float_t chisq=0;
1128     Float_t qcont=0;
1129     Float_t qtot=0;
1130
1131     for (i=0; i<gNbins; i++) {
1132         Float_t q0=gCharge[i];
1133         Float_t q1=DiscrCharge(i,par);
1134         delta=(q0-q1)/TMath::Sqrt(q0);
1135         chisq+=delta*delta;
1136         qcont+=q1;
1137         qtot+=q0;
1138     }
1139     chisq=chisq+=(qtot-qcont)*(qtot-qcont)*0.5;
1140     f=chisq;
1141 }
1142
1143
1144 void AliRICHClusterFinder::SetDigits(TClonesArray *RICHdigits)
1145 {
1146
1147 // Get all the digits
1148
1149     fDigits=RICHdigits;
1150     fNdigits = fDigits->GetEntriesFast();
1151 }
1152
1153 AliRICHClusterFinder& AliRICHClusterFinder::operator=(const AliRICHClusterFinder& rhs)
1154 {
1155 // Assignment operator
1156     return *this;
1157
1158 }