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