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