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