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