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