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