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