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