PatRec adapted to new IO
[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
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 //__________________________________________________________________________________________________