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