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