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