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