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