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