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