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