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 #include "AliRICHClusterFinder.h"
  18 #include "AliRun.h"
  19 #include "AliRICH.h"
  20 #include "AliRICHMap.h"
  21 #include "AliRICHSDigit.h"
  22 #include "AliRICHDigit.h"
  23 #include "AliRICHRawCluster.h"
  24 #include "AliRICHParam.h"
  25
  26
  27
  28 #include <TTree.h>
  29 #include <TCanvas.h>
  30 #include <TH1.h>
  31 #include <TF1.h>
  32 #include <TPad.h>
  33 #include <TGraph.h>
  34 #include <TPostScript.h>
  35 #include <TMinuit.h>
  36
  37 //----------------------------------------------------------
  38 static AliSegmentation     *gSegmentation;
  39 static AliRICHResponse*     gResponse;
  40 static Int_t                gix[500];
  41 static Int_t                giy[500];
  42 static Float_t              gCharge[500];
  43 static Int_t                gNbins;
  44 static Int_t                gFirst=kTRUE;
  45 static TMinuit *gMyMinuit ;
  46 void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
  47 static Int_t                gChargeTot;
  48
  49 ClassImp(AliRICHClusterFinder)
  50 //__________________________________________________________________________________________________
  51 AliRICHClusterFinder::AliRICHClusterFinder(AliRICH *pRICH)
  52 {//main ctor
  53   Info("main ctor","Start.");
  54
  55   fRICH=pRICH;
  56
  57   fSegmentation=Rich()->C(1)->GetSegmentationModel();
  58   fResponse    =Rich()->C(1)->GetResponseModel();
  59
  60   fDigits=0;    fNdigits=0;
  61   fChamber=0;
  62   fHitMap=0;
  63
  64   fCogCorr = 0;
  65   SetNperMax();
  66   SetClusterSize();
  67   SetDeclusterFlag();
  68   fNPeaks=-1;
  69 }//main ctor
  70 //__________________________________________________________________________________________________
  71 void AliRICHClusterFinder::Decluster(AliRICHRawCluster *cluster)
  72 {// Decluster algorithm
  73   Info("Decluster","Start.");
  74   if(cluster->fMultiplicity==1||cluster->fMultiplicity==2){//Nothing special for 1- and 2-clusters
  75     if(fNPeaks != 0) {cluster->fNcluster[0]=fNPeaks; cluster->fNcluster[1]=0;}
  76     AddRawCluster(*cluster);
  77     fNPeaks++;
  78   }else if(cluster->fMultiplicity==3){// 3-cluster, check topology
  79     if(fDeclusterFlag)
  80       Centered(cluster);// ok, cluster is centered and added in Centered()
  81     else{//if(fDeclusterFlag)
  82       if(fNPeaks!=0){cluster->fNcluster[0]=fNPeaks;cluster->fNcluster[1]=0;}
  83       AddRawCluster(*cluster);
  84       fNPeaks++;
  85     }//if(fDeclusterFlag)
  86   }else{//4-and more-pad clusters
  87     if(cluster->fMultiplicity<= fClusterSize){
  88       if(fDeclusterFlag)
  89         SplitByLocalMaxima(cluster);
  90       else{
  91         if(fNPeaks!= 0){cluster->fNcluster[0]=fNPeaks; cluster->fNcluster[1]=0; }
  92         AddRawCluster(*cluster);
  93         fNPeaks++;
  94       }//if(fDeclusterFlag)
  95     }//if <= fClusterSize
  96   }//if multiplicity
  97 }//Decluster()
  98 //__________________________________________________________________________________________________
  99 Bool_t AliRICHClusterFinder::Centered(AliRICHRawCluster *cluster)
 100 {//Is the cluster centered?
 101
 102   AliRICHDigit* dig;
 103   dig= (AliRICHDigit*)fDigits->UncheckedAt(cluster->fIndexMap[0]);
 104   Int_t x[kMaxNeighbours], y[kMaxNeighbours], xN[kMaxNeighbours], yN[kMaxNeighbours];
 105   Int_t nn=Rich()->Param()->Neighbours(dig->PadX(),dig->PadY(),x,y);
 106
 107
 108   Int_t nd=0;
 109   for (Int_t i=0; i<nn; i++){//neighbours loop
 110     if(fHitMap->TestHit(x[i],y[i]) == kUsed){
 111       xN[nd]=x[i];
 112       yN[nd]=y[i];
 113       nd++;
 114     }
 115   }//neighbours loop
 116
 117   if(nd==2){// cluster is centered !
 118         if (fNPeaks != 0) {
 119             cluster->fNcluster[0]=fNPeaks;
 120             cluster->fNcluster[1]=0;
 121         }
 122         cluster->fCtype=0;
 123         AddRawCluster(*cluster);
 124         fNPeaks++;
 125         return kTRUE;
 126     } else if (nd ==1) {
 127 // Highest signal on an edge, split cluster into 2+1
 128 // who is the neighbour ?
 129         Int_t nind=fHitMap->GetHitIndex(xN[0], yN[0]);
 130         Int_t i1= (nind==cluster->fIndexMap[1]) ? 1:2;
 131         Int_t i2= (nind==cluster->fIndexMap[1]) ? 2:1;
 132 // 2-cluster
 133         AliRICHRawCluster cnew;
 134         if (fNPeaks == 0) {
 135             cnew.fNcluster[0]=-1;
 136             cnew.fNcluster[1]=fNRawClusters;
 137         } else {
 138             cnew.fNcluster[0]=fNPeaks;
 139             cnew.fNcluster[1]=0;
 140         }
 141         cnew.fMultiplicity=2;
 142         cnew.fIndexMap[0]=cluster->fIndexMap[0];
 143         cnew.fIndexMap[1]=cluster->fIndexMap[i1];
 144         FillCluster(&cnew);
 145         cnew.fClusterType=cnew.PhysicsContribution();
 146         AddRawCluster(cnew);
 147         fNPeaks++;
 148 // 1-cluster
 149         cluster->fMultiplicity=1;
 150         cluster->fIndexMap[0]=cluster->fIndexMap[i2];
 151         cluster->fIndexMap[1]=0;
 152         cluster->fIndexMap[2]=0;
 153         FillCluster(cluster);
 154         if (fNPeaks != 0) {
 155             cluster->fNcluster[0]=fNPeaks;
 156             cluster->fNcluster[1]=0;
 157         }
 158         cluster->fClusterType=cluster->PhysicsContribution();
 159         AddRawCluster(*cluster);
 160         fNPeaks++;
 161         return kFALSE;
 162     } else {
 163         Warning("Centered","\n Completely screwed up %d !! \n",nd);
 164
 165     }
 166   return kFALSE;
 167 }//Centered()
 168 //__________________________________________________________________________________________________
 169 void AliRICHClusterFinder::SplitByLocalMaxima(AliRICHRawCluster *c)
 170 {// Split the cluster according to the number of maxima inside
 171     AliRICHDigit* dig[100], *digt;
 172     Int_t ix[100], iy[100], q[100];
 173     Float_t x[100], y[100];
 174     Int_t i; // loops over digits
 175     Int_t j; // loops over local maxima
 176     Int_t mul=c->fMultiplicity;
 177 //  dump digit information into arrays
 178   for (i=0; i<mul; i++){
 179     dig[i]= (AliRICHDigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
 180     ix[i]= dig[i]->PadX();
 181     iy[i]= dig[i]->PadY();
 182     q[i] = dig[i]->Signal();
 183     Rich()->Param()->Pad2Local(ix[i], iy[i], x[i], y[i]);
 184   }
 185 //  Find local maxima
 186     Bool_t isLocal[100];
 187     Int_t nLocal=0;
 188     Int_t associatePeak[100];
 189     Int_t indLocal[100];
 190     Int_t nn;
 191     Int_t xNei[kMaxNeighbours], yNei[kMaxNeighbours];
 192     for (i=0; i<mul; i++) {
 193         fSegmentation->Neighbours(ix[i], iy[i], &nn, xNei, yNei);
 194         isLocal[i]=kTRUE;
 195         for (j=0; j<nn; j++) {
 196             if (fHitMap->TestHit(xNei[j], yNei[j])==kEmpty) continue;
 197             digt=(AliRICHDigit*) fHitMap->GetHit(xNei[j], yNei[j]);
 198             if (digt->Signal() > q[i]) {
 199                 isLocal[i]=kFALSE;
 200                 break;
 201 // handle special case of neighbouring pads with equal signal
 202             } else if (digt->Signal() == q[i]) {
 203                 if (nLocal >0) {
 204                     for (Int_t k=0; k<nLocal; k++) {
 205                         if (xNei[j]==ix[indLocal[k]] && yNei[j]==iy[indLocal[k]]){
 206                             isLocal[i]=kFALSE;
 207                         }
 208                     }
 209                 }
 210             }
 211         } // loop over next neighbours
 212         // Maxima should not be on the edge
 213         if (isLocal[i]) {
 214             indLocal[nLocal]=i;
 215             nLocal++;
 216         }
 217     } // loop over all digits
 218 // If only one local maximum found but multiplicity is high take global maximum from the list of digits.
 219     if (nLocal==1 && mul>5) {
 220         Int_t nnew=0;
 221         for (i=0; i<mul; i++) {
 222             if (!isLocal[i]) {
 223                 indLocal[nLocal]=i;
 224                 isLocal[i]=kTRUE;
 225                 nLocal++;
 226                 nnew++;
 227             }
 228             if (nnew==1) break;
 229         }
 230     }
 231 // If number of local maxima is 2 try to fit a double gaussian
 232     if (nLocal==-100) {
 233 //  Initialise global variables for fit
 234         gFirst=1;
 235         gSegmentation=fSegmentation;
 236         gResponse    =fResponse;
 237         gNbins=mul;
 238
 239         for (i=0; i<mul; i++) {
 240             gix[i]=ix[i];
 241             giy[i]=iy[i];
 242             gCharge[i]=Float_t(q[i]);
 243         }
 244         if (gFirst) {
 245             gFirst=kFALSE;
 246             gMyMinuit = new TMinuit(5);
 247         }
 248         gMyMinuit->SetFCN(fcn);
 249         gMyMinuit->mninit(5,10,7);
 250         Double_t arglist[20];
 251         Int_t ierflag=0;
 252         arglist[0]=1;
 253 // Set starting values
 254         static Double_t vstart[5];
 255         vstart[0]=x[indLocal[0]];
 256         vstart[1]=y[indLocal[0]];
 257         vstart[2]=x[indLocal[1]];
 258         vstart[3]=y[indLocal[1]];
 259         vstart[4]=Float_t(q[indLocal[0]])/
 260             Float_t(q[indLocal[0]]+q[indLocal[1]]);
 261 // lower and upper limits
 262         static Double_t lower[5], upper[5];
 263         Int_t isec=fSegmentation->Sector(ix[indLocal[0]], iy[indLocal[0]]);
 264         lower[0]=vstart[0]-fSegmentation->Dpx(isec)/2;
 265         lower[1]=vstart[1]-fSegmentation->Dpy(isec)/2;
 266
 267         upper[0]=lower[0]+fSegmentation->Dpx(isec);
 268         upper[1]=lower[1]+fSegmentation->Dpy(isec);
 269
 270         isec=fSegmentation->Sector(ix[indLocal[1]], iy[indLocal[1]]);
 271         lower[2]=vstart[2]-fSegmentation->Dpx(isec)/2;
 272         lower[3]=vstart[3]-fSegmentation->Dpy(isec)/2;
 273
 274         upper[2]=lower[2]+fSegmentation->Dpx(isec);
 275         upper[3]=lower[3]+fSegmentation->Dpy(isec);
 276
 277         lower[4]=0.;
 278         upper[4]=1.;
 279 // step sizes
 280         static Double_t step[5]={0.005, 0.03, 0.005, 0.03, 0.01};
 281
 282         gMyMinuit->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
 283         gMyMinuit->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
 284         gMyMinuit->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
 285         gMyMinuit->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
 286         gMyMinuit->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
 287 // ready for minimisation
 288         gMyMinuit->SetPrintLevel(-1);
 289         gMyMinuit->mnexcm("SET OUT", arglist, 0, ierflag);
 290         arglist[0]= -1;
 291         arglist[1]= 0;
 292
 293         gMyMinuit->mnexcm("SET NOGR", arglist, 0, ierflag);
 294         gMyMinuit->mnexcm("SCAN", arglist, 0, ierflag);
 295         gMyMinuit->mnexcm("EXIT" , arglist, 0, ierflag);
 296
 297         Double_t xrec[2], yrec[2], qfrac;
 298         TString chname;
 299         Double_t epxz, b1, b2;
 300         Int_t ierflg;
 301         gMyMinuit->mnpout(0, chname, xrec[0], epxz, b1, b2, ierflg);
 302         gMyMinuit->mnpout(1, chname, yrec[0], epxz, b1, b2, ierflg);
 303         gMyMinuit->mnpout(2, chname, xrec[1], epxz, b1, b2, ierflg);
 304         gMyMinuit->mnpout(3, chname, yrec[1], epxz, b1, b2, ierflg);
 305         gMyMinuit->mnpout(4, chname, qfrac,   epxz, b1, b2, ierflg);
 306  // One cluster for each maximum
 307         for (j=0; j<2; j++) {
 308             AliRICHRawCluster cnew;
 309             if (fNPeaks == 0) {
 310                 cnew.fNcluster[0]=-1;
 311                 cnew.fNcluster[1]=fNRawClusters;
 312             } else {
 313                 cnew.fNcluster[0]=fNPeaks;
 314                 cnew.fNcluster[1]=0;
 315             }
 316             cnew.fMultiplicity=0;
 317             cnew.fX=Float_t(xrec[j]);
 318             cnew.fY=Float_t(yrec[j]);
 319             if (j==0) {
 320                 cnew.fQ=Int_t(gChargeTot*qfrac);
 321             } else {
 322                 cnew.fQ=Int_t(gChargeTot*(1-qfrac));
 323             }
 324             gSegmentation->SetHit(xrec[j],yrec[j],0);
 325             for (i=0; i<mul; i++) {
 326                 cnew.fIndexMap[cnew.fMultiplicity]=c->fIndexMap[i];
 327                 gSegmentation->SetPad(gix[i], giy[i]);
 328                 Float_t q1=gResponse->IntXY(gSegmentation);
 329                 cnew.fContMap[cnew.fMultiplicity]=Float_t(q[i])/(q1*cnew.fQ);
 330                 cnew.fMultiplicity++;
 331             }
 332             FillCluster(&cnew,0);
 333             cnew.fClusterType=cnew.PhysicsContribution();
 334             AddRawCluster(cnew);
 335             fNPeaks++;
 336         }
 337     }
 338
 339     Bool_t fitted=kTRUE;
 340
 341     if (nLocal !=-100 || !fitted) {
 342         // Check if enough local clusters have been found, if not add global maxima to the list
 343         Int_t nPerMax;
 344         if (nLocal!=0) {
 345             nPerMax=mul/nLocal;
 346         } else {
 347             Warning("SplitByLocalMaxima","no local maximum found");
 348             nPerMax=fNperMax+1;
 349         }
 350
 351         if (nPerMax > fNperMax) {
 352             Int_t nGlob=mul/fNperMax-nLocal+1;
 353             if (nGlob > 0) {
 354                 Int_t nnew=0;
 355                 for (i=0; i<mul; i++) {
 356                     if (!isLocal[i]) {
 357                         indLocal[nLocal]=i;
 358                         isLocal[i]=kTRUE;
 359                         nLocal++;
 360                         nnew++;
 361                     }
 362                     if (nnew==nGlob) break;
 363                 }
 364             }
 365         }
 366         for (i=0; i<mul; i++) { // Associate hits to peaks
 367             Float_t dmin=1.E10;
 368             Float_t qmax=0;
 369             if (isLocal[i]) continue;
 370             for (j=0; j<nLocal; j++) {
 371                 Int_t il=indLocal[j];
 372                 Float_t d=TMath::Sqrt((x[i]-x[il])*(x[i]-x[il])
 373                                       +(y[i]-y[il])*(y[i]-y[il]));
 374                 Float_t ql=q[il];
 375                 if (d<dmin) {           // Select nearest peak
 376                     dmin=d;
 377                     qmax=ql;
 378                     associatePeak[i]=j;
 379                 } else if (d==dmin) {               // If more than one take highest peak
 380                     if (ql>qmax) {
 381                         dmin=d;
 382                         qmax=ql;
 383                         associatePeak[i]=j;
 384                     }
 385                 }
 386             }
 387         }
 388  // One cluster for each maximum
 389         for (j=0; j<nLocal; j++) {
 390             AliRICHRawCluster cnew;
 391             if (fNPeaks == 0) {
 392                 cnew.fNcluster[0]=-1;
 393                 cnew.fNcluster[1]=fNRawClusters;
 394             } else {
 395                 cnew.fNcluster[0]=fNPeaks;
 396                 cnew.fNcluster[1]=0;
 397             }
 398             cnew.fIndexMap[0]=c->fIndexMap[indLocal[j]];
 399             cnew.fMultiplicity=1;
 400             for (i=0; i<mul; i++) {
 401                 if (isLocal[i]) continue;
 402                 if (associatePeak[i]==j) {
 403                     cnew.fIndexMap[cnew.fMultiplicity]=c->fIndexMap[i];
 404                     cnew.fMultiplicity++;
 405                 }
 406             }
 407             FillCluster(&cnew);
 408             cnew.fClusterType=cnew.PhysicsContribution();
 409             AddRawCluster(cnew);
 410             fNPeaks++;
 411         }
 412     }
 413 }//SplitByLocalMaxima(AliRICHRawCluster *c)
 414 //__________________________________________________________________________________________________
 415 void  AliRICHClusterFinder::FillCluster(AliRICHRawCluster* c, Int_t flag)
 416 {//  Completes cluster information starting from list of digits
 417     AliRICHDigit* dig;
 418     Float_t x, y, z;
 419     Int_t  ix, iy;
 420     Float_t frac=0;
 421
 422     c->fPeakSignal=0;
 423     if (flag) {
 424         c->fX=0;
 425         c->fY=0;
 426         c->fQ=0;
 427     }
 428
 429
 430     for (Int_t i=0; i<c->fMultiplicity; i++){
 431         dig= (AliRICHDigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
 432         ix=dig->PadX()+c->fOffsetMap[i];
 433         iy=dig->PadY();
 434         Int_t q=dig->Signal();
 435         if (dig->Physics() >= dig->Signal()) {
 436           c->fPhysicsMap[i]=2;
 437         } else if (dig->Physics() == 0) {
 438           c->fPhysicsMap[i]=0;
 439         } else  c->fPhysicsMap[i]=1;
 440 // peak signal and track list
 441         if (flag) {
 442            if (q>c->fPeakSignal) {
 443               c->fPeakSignal=q;
 444             c->fTracks[0]=dig->Hit();
 445             c->fTracks[1]=dig->Track(0);
 446             c->fTracks[2]=dig->Track(1);
 447            }
 448         } else {
 449            if (c->fContMap[i] > frac) {
 450               frac=c->fContMap[i];
 451               c->fPeakSignal=q;
 452             c->fTracks[0]=dig->Hit();
 453             c->fTracks[1]=dig->Track(0);
 454             c->fTracks[2]=dig->Track(1);
 455            }
 456         }
 457         if (flag) {
 458             fSegmentation->GetPadC(ix, iy, x, y, z);
 459             c->fX += q*x;
 460             c->fY += q*y;
 461             c->fQ += q;
 462         }
 463
 464     } // loop over digits
 465
 466  if (flag) {
 467
 468      c->fX/=c->fQ;
 469      c->fX=fSegmentation->GetAnod(c->fX);
 470      c->fY/=c->fQ;
 471 //  apply correction to the coordinate along the anode wire
 472      x=c->fX;
 473      y=c->fY;
 474      Rich()->Param()->Local2Pad(x,y,ix,iy);
 475      Rich()->Param()->Pad2Local(ix,iy,x,y);
 476      Int_t isec=fSegmentation->Sector(ix,iy);
 477      TF1* cogCorr = fSegmentation->CorrFunc(isec-1);
 478
 479      if (cogCorr) {
 480          Float_t yOnPad=(c->fY-y)/fSegmentation->Dpy(isec);
 481          c->fY=c->fY-cogCorr->Eval(yOnPad, 0, 0);
 482      }
 483  }
 484 }//FillCluster(AliRICHRawCluster* c, Int_t flag)
 485 //__________________________________________________________________________________________________
 486 void  AliRICHClusterFinder::AddDigit2Cluster(Int_t i, Int_t j, AliRICHRawCluster &c)
 487 {//Find clusters Add i,j as element of the cluster
 488   Info("AddDigit2Cluster","Start with digit(%i,%i)",i,j);
 489
 490   Int_t idx = fHitMap->GetHitIndex(i,j);
 491   AliRICHDigit* dig = (AliRICHDigit*) fHitMap->GetHit(i,j);
 492   Int_t q=dig->Signal();
 493   if(q>TMath::Abs(c.fPeakSignal)){
 494         c.fPeakSignal=q;
 495         c.fTracks[0]=dig->Hit();
 496         c.fTracks[1]=dig->Track(0);
 497         c.fTracks[2]=dig->Track(1);
 498     }
 499 //  Make sure that list of digits is ordered
 500     Int_t mu=c.fMultiplicity;
 501     c.fIndexMap[mu]=idx;
 502
 503     if (dig->Physics() >= dig->Signal()) {
 504         c.fPhysicsMap[mu]=2;
 505     } else if (dig->Physics() == 0) {
 506         c.fPhysicsMap[mu]=0;
 507     } else  c.fPhysicsMap[mu]=1;
 508
 509     if (mu > 0) {
 510         for (Int_t ind=mu-1; ind>=0; ind--) {
 511             Int_t ist=(c.fIndexMap)[ind];
 512             Int_t ql=((AliRICHDigit*)fDigits->UncheckedAt(ist))->Signal();
 513             if (q>ql) {
 514                 c.fIndexMap[ind]=idx;
 515                 c.fIndexMap[ind+1]=ist;
 516             } else {
 517                 break;
 518             }
 519         }
 520     }
 521
 522   c.fMultiplicity++;
 523     if (c.fMultiplicity >= 50 ) {
 524         Info("AddDigit2CLuster","multiplicity >50  %d \n",c.fMultiplicity);
 525         c.fMultiplicity=49;
 526     }
 527   Float_t x,y;// Prepare center of gravity calculation
 528   Rich()->Param()->Pad2Local(i,j,x,y);
 529   c.fX+=q*x;    c.fY+=q*y;    c.fQ += q;
 530   fHitMap->FlagHit(i,j);// Flag hit as taken
 531
 532
 533   Int_t xList[4], yList[4];    //  Now look recursively for all neighbours
 534   for (Int_t iNei=0;iNei<Rich()->Param()->Neighbours(i,j,xList,yList);iNei++)
 535     if(fHitMap->TestHit(xList[iNei],yList[iNei])==kUnused) AddDigit2Cluster(xList[iNei],yList[iNei],c);
 536 }//AddDigit2Cluster()
 537 //__________________________________________________________________________________________________
 538 void AliRICHClusterFinder::FindRawClusters()
 539 {//finds neighbours and fill the tree with raw clusters
 540   Info("FindRawClusters","Start for Chamber %i.",fChamber);
 541
 542   if(!fNdigits)return;
 543
 544   fHitMap=new AliRICHMap(fDigits);
 545
 546   for(Int_t iDigN=0;iDigN<fNdigits;iDigN++){//digits loop
 547     AliRICHDigit *dig=(AliRICHDigit*)fDigits->UncheckedAt(iDigN);
 548     Int_t i=dig->PadX();   Int_t j=dig->PadY();
 549     if(fHitMap->TestHit(i,j)==kUsed||fHitMap->TestHit(i,j)==kEmpty) continue;
 550
 551     AliRICHRawCluster c;
 552     c.fMultiplicity=0; c.fPeakSignal=dig->Signal();
 553     c.fTracks[0]=dig->Hit();c.fTracks[1]=dig->Track(0);c.fTracks[2]=dig->Track(1);
 554     c.fNcluster[0]=-1;// tag the beginning of cluster list in a raw cluster
 555
 556     AddDigit2Cluster(i,j,c);//form initial cluster
 557
 558     c.fX /= c.fQ;       // center of gravity
 559     //c.fX=fSegmentation->GetAnod(c.fX);
 560     c.fY /= c.fQ;
 561     AddRawCluster(c);
 562
 563 //    Int_t ix,iy;//  apply correction to the coordinate along the anode wire
 564 //    Float_t x=c.fX, y=c.fY;
 565 //    Rich()->Param()->Local2Pad(x,y,ix,iy);
 566 //    Rich()->Param()->Pad2Local(ix,iy,x,y);
 567 //    Int_t isec=fSegmentation->Sector(ix,iy);
 568 //    TF1* cogCorr=fSegmentation->CorrFunc(isec-1);
 569 //    if(cogCorr){
 570 //      Float_t yOnPad=(c.fY-y)/fSegmentation->Dpy(isec);
 571 //      c.fY=c.fY-cogCorr->Eval(yOnPad,0,0);
 572 //    }
 573
 574     c.fNcluster[1]=fNRawClusters; c.fClusterType=c.PhysicsContribution();
 575
 576     Decluster(&c);
 577
 578     fNPeaks=0;
 579
 580     c.fMultiplicity=0; for(int k=0;k<c.fMultiplicity;k++) c.fIndexMap[k]=0;//reset cluster object
 581   }//digits loop
 582   delete fHitMap;
 583   Info("FindRawClusters","Stop.");
 584 }//FindRawClusters()
 585 //__________________________________________________________________________________________________
 586 void AliRICHClusterFinder::CalibrateCOG()
 587 {// Calibration
 588
 589     Float_t x[5];
 590     Float_t y[5];
 591     Int_t n, i;
 592     if (fSegmentation) {
 593         TF1 *func;
 594         fSegmentation->GiveTestPoints(n, x, y);
 595         for (i=0; i<n; i++) {
 596             func = 0;
 597             Float_t xtest=x[i];
 598             Float_t ytest=y[i];
 599             SinoidalFit(xtest, ytest, func);
 600             if (func) fSegmentation->SetCorrFunc(i, new TF1(*func));
 601         }
 602     }
 603 }//CalibrateCOG()
 604 //__________________________________________________________________________________________________
 605 void AliRICHClusterFinder::SinoidalFit(Float_t x, Float_t y, TF1 *func)
 606 {//Sinoidal fit
 607   static Int_t count=0;
 608   Float_t z;
 609
 610     count++;
 611
 612     const Int_t kNs=101;
 613     Float_t xg[kNs], yg[kNs], xrg[kNs], yrg[kNs];
 614     Float_t xsig[kNs], ysig[kNs];
 615
 616     Int_t ix,iy;
 617     fSegmentation->GetPadI(x,y,0,ix,iy);
 618     fSegmentation->GetPadC(ix,iy,x,y,z);
 619     Int_t isec=fSegmentation->Sector(ix,iy);
 620 // Pad Limits
 621     Float_t xmin = x-fSegmentation->Dpx(isec)/2;
 622     Float_t ymin = y-fSegmentation->Dpy(isec)/2;
 623 //
 624 //      Integration Limits
 625     Float_t dxI=fResponse->SigmaIntegration()*fResponse->ChargeSpreadX();
 626     Float_t dyI=fResponse->SigmaIntegration()*fResponse->ChargeSpreadY();
 627
 628 //
 629 //  Scanning
 630 //
 631     Int_t i;
 632     Float_t qp=0;
 633
 634 //  y-position
 635     Float_t yscan=ymin;
 636     Float_t dy=fSegmentation->Dpy(isec)/(kNs-1);
 637
 638     for (i=0; i<kNs; i++) {//      Pad Loop
 639         Float_t sum=0;
 640         Float_t qcheck=0;
 641         fSegmentation->SigGenInit(x, yscan, 0);
 642
 643         for (fSegmentation->FirstPad(x, yscan,0, dxI, dyI);
 644              fSegmentation->MorePads();
 645              fSegmentation->NextPad())
 646         {
 647             qp=fResponse->IntXY(fSegmentation);
 648             qp=TMath::Abs(qp);
 649             if (qp > 1.e-4) {
 650                 qcheck+=qp;
 651                 Int_t ixs=fSegmentation->Ix();
 652                 Int_t iys=fSegmentation->Iy();
 653                 Float_t xs,ys,zs;
 654                 fSegmentation->GetPadC(ixs,iys,xs,ys,zs);
 655                 sum+=qp*ys;
 656             }
 657         } // Pad loop
 658         Float_t ycog=sum/qcheck;
 659         yg[i]=(yscan-y)/fSegmentation->Dpy(isec);
 660         yrg[i]=(ycog-y)/fSegmentation->Dpy(isec);
 661         ysig[i]=ycog-yscan;
 662         yscan+=dy;
 663     } // scan loop
 664 //  x-position
 665     Float_t xscan=xmin;
 666     Float_t dx=fSegmentation->Dpx(isec)/(kNs-1);
 667
 668     for (i=0; i<kNs; i++) {//      Pad Loop
 669         Float_t sum=0;
 670         Float_t qcheck=0;
 671         fSegmentation->SigGenInit(xscan, y, 0);
 672
 673         for (fSegmentation->FirstPad(xscan, y, 0, dxI, dyI);
 674              fSegmentation->MorePads();
 675              fSegmentation->NextPad())
 676         {
 677             qp=fResponse->IntXY(fSegmentation);
 678             qp=TMath::Abs(qp);
 679             if (qp > 1.e-2) {
 680                 qcheck+=qp;
 681                 Int_t ixs=fSegmentation->Ix();
 682                 Int_t iys=fSegmentation->Iy();
 683                 Float_t xs,ys,zs;
 684                 fSegmentation->GetPadC(ixs,iys,xs,ys,zs);
 685                 sum+=qp*xs;
 686             }
 687         } // Pad loop
 688         Float_t xcog=sum/qcheck;
 689         xcog=fSegmentation->GetAnod(xcog);
 690
 691         xg[i]=(xscan-x)/fSegmentation->Dpx(isec);
 692         xrg[i]=(xcog-x)/fSegmentation->Dpx(isec);
 693         xsig[i]=xcog-xscan;
 694         xscan+=dx;
 695     }
 696 // Creates a Root function based on function sinoid above and perform the fit
 697     TGraph *graphyr= new TGraph(kNs,yrg,ysig);
 698     Double_t sinoid(Double_t *x, Double_t *par);
 699     new TF1("sinoidf",sinoid,0.5,0.5,5);
 700     graphyr->Fit("sinoidf","Q");
 701     func = (TF1*)graphyr->GetListOfFunctions()->At(0);
 702 }//SinoidalFit()
 703 //__________________________________________________________________________________________________
 704 Double_t sinoid(Double_t *x, Double_t *par)
 705 {// Sinoid function
 706
 707     Double_t arg = -2*TMath::Pi()*x[0];
 708     Double_t fitval= par[0]*TMath::Sin(arg)+
 709         par[1]*TMath::Sin(2*arg)+
 710         par[2]*TMath::Sin(3*arg)+
 711         par[3]*TMath::Sin(4*arg)+
 712         par[4]*TMath::Sin(5*arg);
 713     return fitval;
 714 }//sinoid()
 715 //__________________________________________________________________________________________________
 716 Double_t DoubleGauss(Double_t *x, Double_t *par)
 717 {//Double gaussian function
 718   Double_t arg1 = (x[0]-par[1])/0.18;
 719   Double_t arg2 = (x[0]-par[3])/0.18;
 720   return par[0]*TMath::Exp(-arg1*arg1/2)+par[2]*TMath::Exp(-arg2*arg2/2);
 721 }
 722 //__________________________________________________________________________________________________
 723 Float_t DiscrCharge(Int_t i,Double_t *par)
 724 {
 725 // par[0]    x-position of first  cluster
 726 // par[1]    y-position of first  cluster
 727 // par[2]    x-position of second cluster
 728 // par[3]    y-position of second cluster
 729 // par[4]    charge fraction of first  cluster
 730 // 1-par[4]  charge fraction of second cluster
 731
 732     static Float_t qtot;
 733     if (gFirst) {
 734         qtot=0;
 735         for (Int_t jbin=0; jbin<gNbins; jbin++) {
 736             qtot+=gCharge[jbin];
 737         }
 738         gFirst=0;
 739         //printf("\n sum of charge from DiscrCharge %f\n", qtot);
 740         gChargeTot=Int_t(qtot);
 741
 742     }
 743     gSegmentation->SetPad(gix[i], giy[i]);
 744 //  First Cluster
 745     gSegmentation->SetHit(par[0],par[1],0);
 746     Float_t q1=gResponse->IntXY(gSegmentation);
 747
 748 //  Second Cluster
 749     gSegmentation->SetHit(par[2],par[3],0);
 750     Float_t q2=gResponse->IntXY(gSegmentation);
 751
 752     Float_t value = qtot*(par[4]*q1+(1.-par[4])*q2);
 753     return value;
 754 }//DiscrCharge(Int_t i,Double_t *par)
 755 //__________________________________________________________________________________________________
 756 void fcn(Int_t &npar, Double_t */*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
 757 {// Minimisation function
 758   npar=1;
 759     Int_t i;
 760     Float_t delta;
 761     Float_t chisq=0;
 762     Float_t qcont=0;
 763     Float_t qtot=0;
 764
 765     for (i=0; i<gNbins; i++) {
 766         Float_t q0=gCharge[i];
 767         Float_t q1=DiscrCharge(i,par);
 768         delta=(q0-q1)/TMath::Sqrt(q0);
 769         chisq+=delta*delta;
 770         qcont+=q1;
 771         qtot+=q0;
 772     }
 773     chisq=chisq+=(qtot-qcont)*(qtot-qcont)*0.5;
 774     f=chisq;
 775 }//
 776 //__________________________________________________________________________________________________
 777 void AliRICHClusterFinder::Exec()
 778 {
 779   Info("Exec","Start.");
 780
 781   Rich()->GetLoader()->LoadDigits();
 782
 783   for(Int_t iEventN=0;iEventN<gAlice->GetEventsPerRun();iEventN++){//events loop
 784     gAlice->GetRunLoader()->GetEvent(iEventN);
 785
 786     Rich()->GetLoader()->MakeTree("R");  Rich()->MakeBranch("R");
 787     Rich()->ResetDigitsOld();  Rich()->ResetRawClusters();
 788
 789     Rich()->GetLoader()->TreeD()->GetEntry(0);
 790     for(fChamber=1;fChamber<=kNCH;fChamber++){//chambers loop
 791       fDigits=Rich()->DigitsOld(fChamber); fNdigits=fDigits->GetEntries();
 792
 793       FindRawClusters();
 794
 795     }//chambers loop
 796
 797     Rich()->GetLoader()->TreeR()->Fill();
 798     Rich()->GetLoader()->WriteRecPoints("OVERWRITE");
 799   }//events loop
 800   Rich()->GetLoader()->UnloadDigits(); Rich()->GetLoader()->UnloadRecPoints();
 801   Rich()->ResetDigitsOld();  Rich()->ResetRawClusters();
 802   Info("Exec","Stop.");
 803 }//Exec()
 804 //__________________________________________________________________________________________________