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