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