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