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