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