Avoid global variables in AliMUONClusterFinderVS by seperating the input data for...
[u/mrichter/AliRoot.git] / MUON / AliMUONClusterFinderVS.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 $Log$
17 Revision 1.4  2000/06/27 16:18:47  gosset
18 Finally correct implementation of xm, ym, ixm, iym sizes
19 when at least three local maxima on cathode 1 or on cathode 2
20
21 Revision 1.3  2000/06/22 14:02:45  morsch
22 Parameterised size of xm[], ym[], ixm[], iym[] correctly implemented (PH)
23 Some HP scope problems corrected (PH)
24
25 Revision 1.2  2000/06/15 07:58:48  morsch
26 Code from MUON-dev joined
27
28 Revision 1.1.2.3  2000/06/09 21:58:33  morsch
29 Most coding rule violations corrected.
30
31 Revision 1.1.2.2  2000/02/15 08:33:52  morsch
32 Error in calculation of contribution map for double clusters (Split method) corrected   (A.M.)
33 Error in determination of track list for double cluster (FillCluster method) corrected  (A.M.)
34 Revised and extended SplitByLocalMaxima method (Isabelle Chevrot):
35         - For clusters with more than 2 maxima on one of the cathode planes all valid
36         combinations of maxima on the two cathodes are preserved. The position of the maxima is
37         taken as the hit position.
38         - New FillCluster method with 2 arguments to find tracks associated to the clusters
39         defined above added. (Method destinction by argument list not very elegant in this case,
40         should be revides (A.M.)
41         - Bug in if-statement to handle maximum 1 maximum per plane corrected
42         - Two cluster per cathode but only 1 combination valid is handled.
43         - More rigerous treatment of 1-2 and 2-1 combinations of maxima.
44
45 */
46
47 #include "AliMUONClusterFinderVS.h"
48 #include "AliMUONDigit.h"
49 #include "AliMUONRawCluster.h"
50 #include "AliMUONSegmentation.h"
51 #include "AliMUONResponse.h"
52 #include "AliMUONHitMap.h"
53 #include "AliMUONHitMapA1.h"
54 #include "AliRun.h"
55 #include "AliMUON.h"
56 #include "AliMUONClusterInput.h"
57
58 #include <TTree.h>
59 #include <TCanvas.h>
60 #include <TH1.h>
61 #include <TPad.h>
62 #include <TGraph.h> 
63 #include <TPostScript.h> 
64 #include <TMinuit.h> 
65 #include <stdio.h>
66 #include <iostream.h>
67
68 //_____________________________________________________________________
69 /*
70 static AliMUONSegmentation* fgSegmentation[2];
71 static AliMUONResponse*     fgResponse;
72 static Int_t                fgix[500][2];
73 static Int_t                fgiy[500][2];
74 static Float_t              fgCharge[500][2];
75 static Int_t                fgNbins[2];
76 static Int_t                fgFirst=kTRUE;
77 static Int_t                fgChargeTot[2];
78 static Float_t              fgQtot[2];
79 static TMinuit*             fgMyMinuit ;
80 */
81 // This function is minimized in the double-Mathieson fit
82 void fcnS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
83 void fcnS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
84 void fcnCombiS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
85 void fcnCombiS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
86
87 ClassImp(AliMUONClusterFinderVS)
88
89     AliMUONClusterFinderVS::AliMUONClusterFinderVS
90 (AliMUONSegmentation *seg1, AliMUONSegmentation *seg2,
91  AliMUONResponse *response, 
92  TClonesArray *digits1, TClonesArray *digits2, 
93  Int_t chamber)   
94     :AliMUONClusterFinder(seg1, response, digits1, chamber)
95 {
96 // Constructor
97     fSegmentation[1]=seg2;
98     fDigits2=digits2;
99     fNdigits2 = fDigits2->GetEntriesFast();
100     fHitMap2=0;
101     fTrack[0]=fTrack[1]=-1;
102     
103 }
104
105     AliMUONClusterFinderVS::AliMUONClusterFinderVS()
106         :AliMUONClusterFinder()
107 {
108 // Default constructor
109     fSegmentation[1]=0;
110     fDigits2=0;
111     fNdigits2 = 0;
112     fHitMap2 = 0;
113     fTrack[0]=fTrack[1]=-1;
114 }
115
116 AliMUONClusterFinderVS::AliMUONClusterFinderVS(
117     const AliMUONClusterFinderVS & clusterFinder)
118 {
119 // Dummy copy Constructor
120     ;
121 }
122
123 void AliMUONClusterFinderVS::SetDigits(TClonesArray *MUONdigits1, TClonesArray *MUONdigits2) {
124 // Set pointers to digit lists 
125     fDigits=MUONdigits1;
126     fNdigits = fDigits->GetEntriesFast();
127     fDigits2=MUONdigits2;
128     fNdigits2 = fDigits2->GetEntriesFast();
129 }
130
131 // Get Segmentation
132 AliMUONSegmentation*  AliMUONClusterFinderVS::Segmentation(Int_t i)
133 {
134 // Return pointer to segmentation of cathode plane number 1 (i=0) or 2 (i=1)
135     return ((i==0)? fSegmentation[0] : fSegmentation[1]);
136 }
137
138 // Get Number of Digits
139 Int_t   AliMUONClusterFinderVS::NDigits(Int_t i)
140 {
141 // Return number of digits for cathode plane i+1
142     return ((i==0)? fNdigits : fNdigits2);
143 }
144
145 // Get Digits
146 TClonesArray*  AliMUONClusterFinderVS::Digits(Int_t i)
147 {
148 // Return pointer to digits for cathode plane i+1
149     return ((i==0)? fDigits : fDigits2);
150 }
151     
152
153 AliMUONHitMap*   AliMUONClusterFinderVS::HitMap(Int_t i)
154 {
155 // Return pointer to  HitMap
156     return ((i==0)? fHitMap : fHitMap2);
157 }
158
159 void AliMUONClusterFinderVS::Decluster(AliMUONRawCluster *cluster)
160 {
161 // Decluster by local maxima
162     SplitByLocalMaxima(cluster);
163 }
164
165 void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
166 {
167 // Split complex cluster by local maxima 
168
169     Int_t cath, i;
170
171     AliMUONClusterInput::Instance()->SetCluster(c);
172
173     fMul[0]=c->fMultiplicity[0];
174     fMul[1]=c->fMultiplicity[1];
175
176 //
177 //  dump digit information into arrays
178 //
179
180     Float_t qtot;
181     
182     for (cath=0; cath<2; cath++) {
183         qtot=0;
184         for (i=0; i<fMul[cath]; i++)
185         {
186             // pointer to digit
187             fDig[i][cath]=(AliMUONDigit*)
188                 (Digits(cath)->UncheckedAt(c->fIndexMap[i][cath]));
189             // pad coordinates
190             fIx[i][cath]= fDig[i][cath]->fPadX;
191             fIy[i][cath]= fDig[i][cath]->fPadY;
192             // pad charge
193             fQ[i][cath] = fDig[i][cath]->fSignal;
194             // pad centre coordinates
195             Segmentation(cath)->
196                 GetPadCxy(fIx[i][cath], fIy[i][cath], fX[i][cath], fY[i][cath]);
197         } // loop over cluster digits
198     }  // loop over cathodes
199
200
201     FindLocalMaxima(c);
202
203 //
204 //  Initialise and perform mathieson fits
205     Float_t chi2, oldchi2;
206 //  ++++++++++++++++++*************+++++++++++++++++++++
207 //  (1) No more than one local maximum per cathode plane 
208 //  +++++++++++++++++++++++++++++++*************++++++++
209     if ((fNLocal[0]==1 && (fNLocal[1]==0 ||  fNLocal[1]==1)) || 
210         (fNLocal[0]==0 && fNLocal[1]==1)) {
211
212 // Perform combined single Mathieson fit
213 // Initial values for coordinates (x,y) 
214
215         // One local maximum on cathodes 1 and 2 (X->cathode 2, Y->cathode 1)
216         if (fNLocal[0]==1 &&  fNLocal[1]==1) {
217             fXInit[0]=c->fX[1];
218             fYInit[0]=c->fY[0];
219             // One local maximum on cathode 1 (X,Y->cathode 1)
220         } else if (fNLocal[0]==1) {
221             fXInit[0]=c->fX[0];
222             fYInit[0]=c->fY[0];
223             // One local maximum on cathode 2  (X,Y->cathode 2)
224         } else {
225             fXInit[0]=c->fX[1];
226             fYInit[0]=c->fY[1];
227         }
228         fprintf(stderr,"\n cas (1) CombiSingleMathiesonFit(c)\n");
229         chi2=CombiSingleMathiesonFit(c);
230 //      Int_t ndf = fgNbins[0]+fgNbins[1]-2;
231 //      Float_t prob = TMath::Prob(Double_t(chi2),ndf);
232 //      prob1->Fill(prob);
233 //      chi2_1->Fill(chi2);
234         oldchi2=chi2;
235         fprintf(stderr," chi2 %f ",chi2);
236
237         c->fX[0]=fXFit[0];
238         c->fY[0]=fYFit[0];
239
240         c->fX[1]=fXFit[0];
241         c->fY[1]=fYFit[0];
242         c->fChi2[0]=chi2;
243         c->fChi2[1]=chi2;
244         c->fX[0]=Segmentation(0)->GetAnod(c->fX[0]);
245         c->fX[1]=Segmentation(1)->GetAnod(c->fX[1]);
246         
247 // If reasonable chi^2 add result to the list of rawclusters
248         //      if (chi2 < 50) {
249         if (chi2 < 0.3) {
250             AddRawCluster(*c);
251 // If not try combined double Mathieson Fit
252         } else {
253             fprintf(stderr," MAUVAIS CHI2 !!!\n");
254             if (fNLocal[0]==1 &&  fNLocal[1]==1) {
255                 fXInit[0]=fX[fIndLocal[0][1]][1];
256                 fYInit[0]=fY[fIndLocal[0][0]][0];
257                 fXInit[1]=fX[fIndLocal[0][1]][1];
258                 fYInit[1]=fY[fIndLocal[0][0]][0];
259             } else if (fNLocal[0]==1) {
260                 fXInit[0]=fX[fIndLocal[0][0]][0];
261                 fYInit[0]=fY[fIndLocal[0][0]][0];
262                 fXInit[1]=fX[fIndLocal[0][0]][0];
263                 fYInit[1]=fY[fIndLocal[0][0]][0];
264             } else {
265                 fXInit[0]=fX[fIndLocal[0][1]][1];
266                 fYInit[0]=fY[fIndLocal[0][1]][1];
267                 fXInit[1]=fX[fIndLocal[0][1]][1];
268                 fYInit[1]=fY[fIndLocal[0][1]][1];
269             }
270             
271 //  Initial value for charge ratios
272             fQrInit[0]=0.5;
273             fQrInit[1]=0.5;
274             fprintf(stderr,"\n cas (1) CombiDoubleMathiesonFit(c)\n");
275             chi2=CombiDoubleMathiesonFit(c);
276 //          Int_t ndf = fgNbins[0]+fgNbins[1]-6;
277 //          Float_t prob = TMath::Prob(chi2,ndf);
278 //          prob2->Fill(prob);
279 //          chi2_2->Fill(chi2);
280             
281 // Was this any better ??
282             fprintf(stderr," Old and new chi2 %f %f ", oldchi2, chi2);
283             if (fFitStat!=0 && chi2>0 && (2.*chi2 < oldchi2)) {
284                 fprintf(stderr," Split\n");
285                 // Split cluster into two according to fit result
286                 Split(c);
287             } else {
288                 fprintf(stderr," Don't Split\n");
289                 // Don't split
290                 AddRawCluster(*c);
291             }
292         }
293
294 //  +++++++++++++++++++++++++++++++++++++++
295 //  (2) Two local maxima per cathode plane 
296 //  +++++++++++++++++++++++++++++++++++++++
297     } else if (fNLocal[0]==2 &&  fNLocal[1]==2) {
298 //
299 //  Let's look for ghosts first 
300 //
301         Float_t xm[4][2], ym[4][2];
302         Float_t dpx, dpy, dx, dy;
303         Int_t ixm[4][2], iym[4][2];
304         Int_t isec, im1, im2, ico;
305 //
306 //  Form the 2x2 combinations
307 //  0-0, 0-1, 1-0, 1-1  
308         ico=0;
309         for (im1=0; im1<2; im1++) {
310             for (im2=0; im2<2; im2++) {     
311                 xm[ico][0]=fX[fIndLocal[im1][0]][0];
312                 ym[ico][0]=fY[fIndLocal[im1][0]][0];
313                 xm[ico][1]=fX[fIndLocal[im2][1]][1];
314                 ym[ico][1]=fY[fIndLocal[im2][1]][1];
315
316                 ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
317                 iym[ico][0]=fIy[fIndLocal[im1][0]][0];
318                 ixm[ico][1]=fIx[fIndLocal[im2][1]][1];
319                 iym[ico][1]=fIy[fIndLocal[im2][1]][1];
320                 ico++;
321             }
322         }
323 // ico = 0 : first local maximum on cathodes 1 and 2
324 // ico = 1 : fisrt local maximum on cathode 1 and second on cathode 2
325 // ico = 2 : second local maximum on cathode 1 and first on cathode 1
326 // ico = 3 : second local maximum on cathodes 1 and 2
327
328 // Analyse the combinations and keep those that are possible !
329 // For each combination check consistency in x and y    
330         Int_t iacc;
331         Bool_t accepted[4];
332         iacc=0;
333         
334         for (ico=0; ico<4; ico++) {
335             accepted[ico]=kFALSE;
336 // cathode one: x-coordinate
337             isec=Segmentation(0)->Sector(ixm[ico][0], iym[ico][0]);
338             dpx=Segmentation(0)->Dpx(isec)/2.;
339             dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
340 // cathode two: y-coordinate
341             isec=Segmentation(1)->Sector(ixm[ico][1], iym[ico][1]);
342             dpy=Segmentation(1)->Dpy(isec)/2.;
343             dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
344 //          printf("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy );
345             if ((dx <= dpx) && (dy <= dpy)) {
346                 // consistent
347                 accepted[ico]=kTRUE;
348                 iacc++;
349             } else {
350                 // reject
351                 accepted[ico]=kFALSE;
352             }
353         }
354
355         if (iacc==2) {
356             fprintf(stderr,"\n iacc=2: No problem ! \n");
357         } else if (iacc==4) {
358             fprintf(stderr,"\n iacc=4: Ok, but ghost problem !!! \n");
359         } else if (iacc==0) {
360             fprintf(stderr,"\n iacc=0: I don't know what to do with this !!!!!!!!! \n");
361         }
362
363 //  Initial value for charge ratios
364         fQrInit[0]=Float_t(fQ[fIndLocal[0][0]][0])/
365             Float_t(fQ[fIndLocal[0][0]][0]+fQ[fIndLocal[1][0]][0]);
366         fQrInit[1]=Float_t(fQ[fIndLocal[0][1]][1])/
367             Float_t(fQ[fIndLocal[0][1]][1]+fQ[fIndLocal[1][1]][1]);
368         
369 // ******* iacc = 0 *******
370 // No combinations found between the 2 cathodes
371 // We keep the center of gravity of the cluster
372         if (iacc==0) {
373             AddRawCluster(*c);
374         }
375
376 // ******* iacc = 1 *******
377 // Only one combination found between the 2 cathodes
378         if (iacc==1) {
379
380 // Initial values for the 2 maxima (x,y)
381
382 // 1 maximum is initialised with the maximum of the combination found (X->cathode 2, Y->cathode 1)
383 // 1 maximum is initialised with the other maximum of the first cathode  
384             if (accepted[0]){
385                 fprintf(stderr,"ico=0\n");
386                 fXInit[0]=xm[0][1];
387                 fYInit[0]=ym[0][0];
388                 fXInit[1]=xm[3][0];
389                 fYInit[1]=ym[3][0];
390             } else if (accepted[1]){
391                 fprintf(stderr,"ico=1\n");
392                 fXInit[0]=xm[1][1];
393                 fYInit[0]=ym[1][0];
394                 fXInit[1]=xm[2][0];
395                 fYInit[1]=ym[2][0];
396             } else if (accepted[2]){
397                 fprintf(stderr,"ico=2\n");
398                 fXInit[0]=xm[2][1];
399                 fYInit[0]=ym[2][0];
400                 fXInit[1]=xm[1][0];
401                 fYInit[1]=ym[1][0];
402             } else if (accepted[3]){
403                 fprintf(stderr,"ico=3\n");
404                 fXInit[0]=xm[3][1];
405                 fYInit[0]=ym[3][0];
406                 fXInit[1]=xm[0][0];
407                 fYInit[1]=ym[0][0];
408             }
409             fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
410             chi2=CombiDoubleMathiesonFit(c);
411 //          Int_t ndf = fgNbins[0]+fgNbins[1]-6;
412 //          Float_t prob = TMath::Prob(chi2,ndf);
413 //          prob2->Fill(prob);
414 //          chi2_2->Fill(chi2);
415             fprintf(stderr," chi2 %f\n",chi2);
416
417 // If reasonable chi^2 add result to the list of rawclusters
418             if (chi2<10) {
419                 Split(c);
420
421             } else {
422 // 1 maximum is initialised with the maximum of the combination found (X->cathode 2, Y->cathode 1)
423 // 1 maximum is initialised with the other maximum of the second cathode  
424                 if (accepted[0]){
425                     fprintf(stderr,"ico=0\n");
426                     fXInit[0]=xm[0][1];
427                     fYInit[0]=ym[0][0];
428                     fXInit[1]=xm[3][1];
429                     fYInit[1]=ym[3][1];
430                 } else if (accepted[1]){
431                     fprintf(stderr,"ico=1\n");
432                     fXInit[0]=xm[1][1];
433                     fYInit[0]=ym[1][0];
434                     fXInit[1]=xm[2][1];
435                     fYInit[1]=ym[2][1];
436                 } else if (accepted[2]){
437                     fprintf(stderr,"ico=2\n");
438                     fXInit[0]=xm[2][1];
439                     fYInit[0]=ym[2][0];
440                     fXInit[1]=xm[1][1];
441                     fYInit[1]=ym[1][1];
442                 } else if (accepted[3]){
443                     fprintf(stderr,"ico=3\n");
444                     fXInit[0]=xm[3][1];
445                     fYInit[0]=ym[3][0];
446                     fXInit[1]=xm[0][1];
447                     fYInit[1]=ym[0][1];
448                 }
449                 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
450                 chi2=CombiDoubleMathiesonFit(c);
451 //              Int_t ndf = fgNbins[0]+fgNbins[1]-6;
452 //              Float_t prob = TMath::Prob(chi2,ndf);
453 //              prob2->Fill(prob);
454 //              chi2_2->Fill(chi2);
455                 fprintf(stderr," chi2 %f\n",chi2);
456
457 // If reasonable chi^2 add result to the list of rawclusters
458                 if (chi2<10) {
459                     Split(c);
460                 } else {
461 //We keep only the combination found (X->cathode 2, Y->cathode 1)
462                     for (Int_t ico=0; ico<2; ico++) {
463                         if (accepted[ico]) {
464                             AliMUONRawCluster cnew;
465                             Int_t cath;    
466                             for (cath=0; cath<2; cath++) {
467                                 cnew.fX[cath]=Float_t(xm[ico][1]);
468                                 cnew.fY[cath]=Float_t(ym[ico][0]);
469                                 cnew.fMultiplicity[cath]=c->fMultiplicity[cath];
470                                 for (i=0; i<fMul[cath]; i++) {
471                                     cnew.fIndexMap[i][cath]=c->fIndexMap[i][cath];
472                                     fSegmentation[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
473                                 }
474                                 fprintf(stderr,"\nRawCluster %d cath %d\n",ico,cath);
475                                 fprintf(stderr,"mult_av %d\n",c->fMultiplicity[cath]);
476                                 FillCluster(&cnew,cath);
477                             } 
478                             cnew.fClusterType=cnew.PhysicsContribution();
479                             AddRawCluster(cnew);
480                             fNPeaks++;
481                         }
482                     }
483                 }
484             }
485         }
486         
487 // ******* iacc = 2 *******
488 // Two combinations found between the 2 cathodes
489         if (iacc==2) {
490
491 // Was the same maximum taken twice
492             if ((accepted[0]&&accepted[1]) || (accepted[2]&&accepted[3])) {
493                 fprintf(stderr,"\n Maximum taken twice !!!\n");
494
495 // Have a try !! with that 
496                 if (accepted[0]&&accepted[3]) {
497                     fXInit[0]=xm[0][1];
498                     fYInit[0]=ym[0][0];
499                     fXInit[1]=xm[1][1];
500                     fYInit[1]=ym[1][0];
501                 } else {
502                     fXInit[0]=xm[2][1];
503                     fYInit[0]=ym[2][0];
504                     fXInit[1]=xm[3][1];
505                     fYInit[1]=ym[3][0];
506                 }
507                 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
508                 chi2=CombiDoubleMathiesonFit(c);
509 //                  Int_t ndf = fgNbins[0]+fgNbins[1]-6;
510 //                  Float_t prob = TMath::Prob(chi2,ndf);
511 //                  prob2->Fill(prob);
512 //                  chi2_2->Fill(chi2);
513                 Split(c);
514                 
515             } else {
516 // No ghosts ! No Problems ! -  Perform one fit only !
517                 if (accepted[0]&&accepted[3]) {
518                     fXInit[0]=xm[0][1];
519                     fYInit[0]=ym[0][0];
520                     fXInit[1]=xm[3][1];
521                     fYInit[1]=ym[3][0];
522                 } else {
523                     fXInit[0]=xm[1][1];
524                     fYInit[0]=ym[1][0];
525                     fXInit[1]=xm[2][1];
526                     fYInit[1]=ym[2][0];
527                 }
528                 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
529                 chi2=CombiDoubleMathiesonFit(c);
530 //                  Int_t ndf = fgNbins[0]+fgNbins[1]-6;
531 //                  Float_t prob = TMath::Prob(chi2,ndf);
532 //                  prob2->Fill(prob);
533 //                  chi2_2->Fill(chi2);
534                 fprintf(stderr," chi2 %f\n",chi2);
535                 Split(c);
536             }
537             
538 // ******* iacc = 4 *******
539 // Four combinations found between the 2 cathodes
540 // Ghost !!
541         } else if (iacc==4) {
542 // Perform fits for the two possibilities !!    
543             fXInit[0]=xm[0][1];
544             fYInit[0]=ym[0][0];
545             fXInit[1]=xm[3][1];
546             fYInit[1]=ym[3][0];
547             fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
548             chi2=CombiDoubleMathiesonFit(c);
549 //              Int_t ndf = fgNbins[0]+fgNbins[1]-6;
550 //              Float_t prob = TMath::Prob(chi2,ndf);
551 //              prob2->Fill(prob);
552 //              chi2_2->Fill(chi2);
553             fprintf(stderr," chi2 %f\n",chi2);
554             Split(c);
555             fXInit[0]=xm[1][1];
556             fYInit[0]=ym[1][0];
557             fXInit[1]=xm[2][1];
558             fYInit[1]=ym[2][0];
559             fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
560             chi2=CombiDoubleMathiesonFit(c);
561 //              ndf = fgNbins[0]+fgNbins[1]-6;
562 //              prob = TMath::Prob(chi2,ndf);
563 //              prob2->Fill(prob);
564 //              chi2_2->Fill(chi2);
565             fprintf(stderr," chi2 %f\n",chi2);
566             Split(c);
567         }
568
569     } else if (fNLocal[0]==2 &&  fNLocal[1]==1) {
570 //  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
571 //  (3) Two local maxima on cathode 1 and one maximum on cathode 2 
572 //  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
573 //
574         Float_t xm[4][2], ym[4][2];
575         Float_t dpx, dpy, dx, dy;
576         Int_t ixm[4][2], iym[4][2];
577         Int_t isec, im1, ico;
578 //
579 //  Form the 2x2 combinations
580 //  0-0, 0-1, 1-0, 1-1  
581         ico=0;
582         for (im1=0; im1<2; im1++) {
583             xm[ico][0]=fX[fIndLocal[im1][0]][0];
584             ym[ico][0]=fY[fIndLocal[im1][0]][0];
585             xm[ico][1]=fX[fIndLocal[0][1]][1];
586             ym[ico][1]=fY[fIndLocal[0][1]][1];
587             
588             ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
589             iym[ico][0]=fIy[fIndLocal[im1][0]][0];
590             ixm[ico][1]=fIx[fIndLocal[0][1]][1];
591             iym[ico][1]=fIy[fIndLocal[0][1]][1];
592             ico++;
593         }
594 // ico = 0 : first local maximum on cathodes 1 and 2
595 // ico = 1 : second local maximum on cathode 1 and first on cathode 2
596
597 // Analyse the combinations and keep those that are possible !
598 // For each combination check consistency in x and y    
599         Int_t iacc;
600         Bool_t accepted[4];
601         iacc=0;
602         
603         for (ico=0; ico<2; ico++) {
604             accepted[ico]=kFALSE;
605             isec=Segmentation(0)->Sector(ixm[ico][0], iym[ico][0]);
606             dpx=Segmentation(0)->Dpx(isec)/2.;
607             dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
608             isec=Segmentation(1)->Sector(ixm[ico][1], iym[ico][1]);
609             dpy=Segmentation(1)->Dpy(isec)/2.;
610             dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
611 //          printf("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy );
612             if ((dx <= dpx) && (dy <= dpy)) {
613                 // consistent
614                 accepted[ico]=kTRUE;
615                 iacc++;
616             } else {
617                 // reject
618                 accepted[ico]=kFALSE;
619             }
620         }
621         
622         Float_t chi21 = 100;
623         Float_t chi22 = 100;
624         
625         if (accepted[0]) {
626             fXInit[0]=xm[0][1];
627             fYInit[0]=ym[0][0];
628             fXInit[1]=xm[1][0];
629             fYInit[1]=ym[1][0];
630             chi21=CombiDoubleMathiesonFit(c);
631 //          Int_t ndf = fgNbins[0]+fgNbins[1]-6;
632 //          Float_t prob = TMath::Prob(chi2,ndf);
633 //          prob2->Fill(prob);
634 //          chi2_2->Fill(chi21);
635             fprintf(stderr," chi2 %f\n",chi21);
636             if (chi21<10) Split(c);
637         } else if (accepted[1]) {
638             fXInit[0]=xm[1][1];
639             fYInit[0]=ym[1][0];
640             fXInit[1]=xm[0][0];
641             fYInit[1]=ym[0][0];
642             chi22=CombiDoubleMathiesonFit(c);
643 //          Int_t ndf = fgNbins[0]+fgNbins[1]-6;
644 //          Float_t prob = TMath::Prob(chi2,ndf);
645 //          prob2->Fill(prob);
646 //          chi2_2->Fill(chi22);
647             fprintf(stderr," chi2 %f\n",chi22);
648             if (chi22<10) Split(c);
649         }
650
651         if (chi21 > 10 && chi22 > 10) {
652 // We keep only the combination found (X->cathode 2, Y->cathode 1)
653             for (Int_t ico=0; ico<2; ico++) {
654                 if (accepted[ico]) {
655                     AliMUONRawCluster cnew;
656                     Int_t cath;    
657                     for (cath=0; cath<2; cath++) {
658                         cnew.fX[cath]=Float_t(xm[ico][1]);
659                         cnew.fY[cath]=Float_t(ym[ico][0]);
660                         cnew.fMultiplicity[cath]=c->fMultiplicity[cath];
661                         for (i=0; i<fMul[cath]; i++) {
662                             cnew.fIndexMap[i][cath]=c->fIndexMap[i][cath];
663                             fSegmentation[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
664                         }
665                         fprintf(stderr,"\nRawCluster %d cath %d\n",ico,cath);
666                         fprintf(stderr,"mult_av %d\n",c->fMultiplicity[cath]);
667                         FillCluster(&cnew,cath);
668                     } 
669                     cnew.fClusterType=cnew.PhysicsContribution();
670                     AddRawCluster(cnew);
671                     fNPeaks++;
672                 }
673             }
674         }
675         
676 //  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
677 //  (3') One local maximum on cathode 1 and two maxima on cathode 2 
678 //  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
679     } else if (fNLocal[0]==1 && fNLocal[1]==2) {
680         
681         Float_t xm[4][2], ym[4][2];
682         Float_t dpx, dpy, dx, dy;
683         Int_t ixm[4][2], iym[4][2];
684         Int_t isec, im1, ico;
685 //
686 //  Form the 2x2 combinations
687 //  0-0, 0-1, 1-0, 1-1  
688         ico=0;
689         for (im1=0; im1<2; im1++) {
690             xm[ico][0]=fX[fIndLocal[0][0]][0];
691             ym[ico][0]=fY[fIndLocal[0][0]][0];
692             xm[ico][1]=fX[fIndLocal[im1][1]][1];
693             ym[ico][1]=fY[fIndLocal[im1][1]][1];
694             
695             ixm[ico][0]=fIx[fIndLocal[0][0]][0];
696             iym[ico][0]=fIy[fIndLocal[0][0]][0];
697             ixm[ico][1]=fIx[fIndLocal[im1][1]][1];
698             iym[ico][1]=fIy[fIndLocal[im1][1]][1];
699             ico++;
700         }
701 // ico = 0 : first local maximum on cathodes 1 and 2
702 // ico = 1 : first local maximum on cathode 1 and second on cathode 2
703
704 // Analyse the combinations and keep those that are possible !
705 // For each combination check consistency in x and y    
706         Int_t iacc;
707         Bool_t accepted[4];
708         iacc=0;
709         
710         for (ico=0; ico<2; ico++) {
711             accepted[ico]=kFALSE;
712             isec=Segmentation(0)->Sector(ixm[ico][0], iym[ico][0]);
713             dpx=Segmentation(0)->Dpx(isec)/2.;
714             dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
715             isec=Segmentation(1)->Sector(ixm[ico][1], iym[ico][1]);
716             dpy=Segmentation(1)->Dpy(isec)/2.;
717             dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
718 //          printf("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy );
719             if ((dx <= dpx) && (dy <= dpy)) {
720                 // consistent
721                 accepted[ico]=kTRUE;
722                 fprintf(stderr,"ico %d\n",ico);
723                 iacc++;
724             } else {
725                 // reject
726                 accepted[ico]=kFALSE;
727             }
728         }
729
730         Float_t chi21 = 100;
731         Float_t chi22 = 100;
732
733         if (accepted[0]) {
734             fXInit[0]=xm[0][0];
735             fYInit[0]=ym[0][1];
736             fXInit[1]=xm[1][1];
737             fYInit[1]=ym[1][1];
738             chi21=CombiDoubleMathiesonFit(c);
739 //          Int_t ndf = fgNbins[0]+fgNbins[1]-6;
740 //          Float_t prob = TMath::Prob(chi2,ndf);
741 //          prob2->Fill(prob);
742 //          chi2_2->Fill(chi21);
743             fprintf(stderr," chi2 %f\n",chi21);
744             if (chi21<10) Split(c);
745         } else if (accepted[1]) {
746             fXInit[0]=xm[1][0];
747             fYInit[0]=ym[1][1];
748             fXInit[1]=xm[0][1];
749             fYInit[1]=ym[0][1];
750             chi22=CombiDoubleMathiesonFit(c);
751 //          Int_t ndf = fgNbins[0]+fgNbins[1]-6;
752 //          Float_t prob = TMath::Prob(chi2,ndf);
753 //          prob2->Fill(prob);
754 //          chi2_2->Fill(chi22);
755             fprintf(stderr," chi2 %f\n",chi22);
756             if (chi22<10) Split(c);
757         }
758
759         if (chi21 > 10 && chi22 > 10) {
760 //We keep only the combination found (X->cathode 2, Y->cathode 1)
761             for (Int_t ico=0; ico<2; ico++) {
762                 if (accepted[ico]) {
763                     AliMUONRawCluster cnew;
764                     Int_t cath;    
765                     for (cath=0; cath<2; cath++) {
766                         cnew.fX[cath]=Float_t(xm[ico][1]);
767                         cnew.fY[cath]=Float_t(ym[ico][0]);
768                         cnew.fMultiplicity[cath]=c->fMultiplicity[cath];
769                         for (i=0; i<fMul[cath]; i++) {
770                             cnew.fIndexMap[i][cath]=c->fIndexMap[i][cath];
771                             fSegmentation[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
772                         }
773                         fprintf(stderr,"\nRawCluster %d cath %d\n",ico,cath);
774                         fprintf(stderr,"mult_av %d\n",c->fMultiplicity[cath]);
775                         FillCluster(&cnew,cath);
776                     } 
777                     cnew.fClusterType=cnew.PhysicsContribution();
778                     AddRawCluster(cnew);
779                     fNPeaks++;
780                 }
781             }
782         }
783
784 //  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
785 //  (4) At least three local maxima on cathode 1 or on cathode 2 
786 //  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
787     } else if (fNLocal[0]>2 || fNLocal[1]>2) {
788         
789         Int_t param = fNLocal[0]*fNLocal[1];
790         Int_t ii;
791
792         Float_t ** xm = new Float_t * [param];
793         for (ii=0; ii<param; ii++) xm[ii]=new Float_t [2];
794         Float_t ** ym = new Float_t * [param];
795         for (ii=0; ii<param; ii++) ym[ii]=new Float_t [2];
796         Int_t ** ixm = new Int_t * [param];
797         for (ii=0; ii<param; ii++) ixm[ii]=new Int_t [2];
798         Int_t ** iym = new Int_t * [param];
799         for (ii=0; ii<param; ii++) iym[ii]=new Int_t [2];
800         
801         Int_t isec, ico;
802         Float_t dpx, dpy, dx, dy;
803
804         ico=0;
805         for (Int_t im1=0; im1<fNLocal[0]; im1++) {
806             for (Int_t im2=0; im2<fNLocal[1]; im2++) {
807                 xm[ico][0]=fX[fIndLocal[im1][0]][0];
808                 ym[ico][0]=fY[fIndLocal[im1][0]][0];
809                 xm[ico][1]=fX[fIndLocal[im2][1]][1];
810                 ym[ico][1]=fY[fIndLocal[im2][1]][1];
811
812                 ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
813                 iym[ico][0]=fIy[fIndLocal[im1][0]][0];
814                 ixm[ico][1]=fIx[fIndLocal[im2][1]][1];
815                 iym[ico][1]=fIy[fIndLocal[im2][1]][1];
816                 ico++;
817             }
818         }
819         
820         Int_t nIco = ico;
821         
822         fprintf(stderr,"nIco %d\n",nIco);
823         for (ico=0; ico<nIco; ico++) {
824             fprintf(stderr,"ico = %d\n",ico);
825             isec=Segmentation(0)->Sector(ixm[ico][0], iym[ico][0]);
826             dpx=Segmentation(0)->Dpx(isec)/2.;
827             dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
828             isec=Segmentation(1)->Sector(ixm[ico][1], iym[ico][1]);
829             dpy=Segmentation(1)->Dpy(isec)/2.;
830             dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
831
832             fprintf(stderr,"dx %f dpx %f dy %f dpy %f\n",dx,dpx,dy,dpy);
833             fprintf(stderr,"  X %f Y %f\n",xm[ico][1],ym[ico][0]);
834             if ((dx <= dpx) && (dy <= dpy)) {
835                 fprintf(stderr,"ok\n");
836                 Int_t cath;    
837                 AliMUONRawCluster cnew;
838                 for (cath=0; cath<2; cath++) {
839                     cnew.fX[cath]=Float_t(xm[ico][1]);
840                     cnew.fY[cath]=Float_t(ym[ico][0]);
841                     cnew.fMultiplicity[cath]=c->fMultiplicity[cath];
842                     for (i=0; i<fMul[cath]; i++) {
843                         cnew.fIndexMap[i][cath]=c->fIndexMap[i][cath];
844                         fSegmentation[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
845                     }
846                     FillCluster(&cnew,cath);
847                 } 
848                 cnew.fClusterType=cnew.PhysicsContribution();
849                 AddRawCluster(cnew);
850                 fNPeaks++;
851             }
852         }
853         delete [] xm;
854         delete [] ym;
855         delete [] ixm;
856         delete [] iym;
857     }
858 }
859
860 void AliMUONClusterFinderVS::FindLocalMaxima(AliMUONRawCluster* c)
861 {
862 // Find all local maxima of a cluster
863    
864     AliMUONDigit* digt;
865     
866     Int_t cath, cath1; // loops over cathodes
867     Int_t i;           // loops over digits
868     Int_t j;           // loops over cathodes
869 //
870 //  Find local maxima
871 //
872 //  counters for number of local maxima
873     fNLocal[0]=fNLocal[1]=0;
874 //  flags digits as local maximum
875     Bool_t isLocal[100][2];
876     for (i=0; i<100;i++) {
877         isLocal[i][0]=isLocal[i][1]=kFALSE;
878     }
879 //  number of next neighbours and arrays to store them 
880     Int_t nn;
881     Int_t x[kMaxNeighbours], y[kMaxNeighbours];
882 // loop over cathodes
883     for (cath=0; cath<2; cath++) {
884 // loop over cluster digits
885         for (i=0; i<fMul[cath]; i++) {
886 // get neighbours for that digit and assume that it is local maximum        
887             Segmentation(cath)->Neighbours(fIx[i][cath], fIy[i][cath], &nn, x, y);
888             isLocal[i][cath]=kTRUE;
889             Int_t isec= Segmentation(cath)->Sector(fIx[i][cath], fIy[i][cath]);
890             Float_t a0 = Segmentation(cath)->Dpx(isec)*Segmentation(cath)->Dpy(isec);
891 // loop over next neighbours, if at least one neighbour has higher charger assumption
892 // digit is not local maximum 
893             for (j=0; j<nn; j++) {
894                 if (HitMap(cath)->TestHit(x[j], y[j])==kEmpty) continue;
895                 digt=(AliMUONDigit*) HitMap(cath)->GetHit(x[j], y[j]);
896                 isec=Segmentation(cath)->Sector(x[j], y[j]);
897                 Float_t a1 = Segmentation(cath)->Dpx(isec)*Segmentation(cath)->Dpy(isec);
898                 if (digt->fSignal/a1 > fQ[i][cath]/a0) {
899                     isLocal[i][cath]=kFALSE;
900                     break;
901 //
902 // handle special case of neighbouring pads with equal signal
903                 } else if (digt->fSignal == fQ[i][cath]) {
904                     if (fNLocal[cath]>0) {
905                         for (Int_t k=0; k<fNLocal[cath]; k++) {
906                             if (x[j]==fIx[fIndLocal[k][cath]][cath] 
907                                 && y[j]==fIy[fIndLocal[k][cath]][cath])
908                             {
909                                 isLocal[i][cath]=kFALSE;
910                             } 
911                         } // loop over local maxima
912                     } // are there already local maxima
913                 } // same charge ? 
914             } // loop over next neighbours
915             if (isLocal[i][cath]) {
916                 fIndLocal[fNLocal[cath]][cath]=i;
917                 fNLocal[cath]++;
918             } 
919         } // loop over all digits
920     } // loop over cathodes
921     
922     printf("\n Found %d %d %d %d local Maxima\n",
923            fNLocal[0], fNLocal[1], fMul[0], fMul[1]);
924     fprintf(stderr,"\n Cathode 1 local Maxima %d Multiplicite %d\n",fNLocal[0], fMul[0]);
925     fprintf(stderr," Cathode 2 local Maxima %d Multiplicite %d\n",fNLocal[1], fMul[1]);
926     Int_t ix, iy, isec;
927     Float_t dpx, dpy;
928     
929     
930     if (fNLocal[1]==2 &&  (fNLocal[0]==1 || fNLocal[0]==0)) {
931         Int_t iback=fNLocal[0];
932         
933 //  Two local maxima on cathode 2 and one maximum on cathode 1 
934 //  Look for local maxima considering up and down neighbours on the 1st cathode only
935 //
936 //  Loop over cluster digits
937         cath=0;
938         cath1=1;
939         
940         for (i=0; i<fMul[cath]; i++) {
941             isec=Segmentation(cath)->Sector(fIx[i][cath],fIy[i][cath]);
942             dpy=Segmentation(cath)->Dpy(isec);
943             dpx=Segmentation(cath)->Dpx(isec);
944             if (isLocal[i][cath]) continue;
945 // Pad position should be consistent with position of local maxima on the opposite cathode
946             if ((TMath::Abs(fX[i][cath]-fX[fIndLocal[0][cath1]][cath1]) > dpx/2.) && 
947                 (TMath::Abs(fX[i][cath]-fX[fIndLocal[1][cath1]][cath1]) > dpx/2.))
948                 continue;
949
950 // get neighbours for that digit and assume that it is local maximum        
951             isLocal[i][cath]=kTRUE;
952 // compare signal to that on the two neighbours on the left and on the right
953             Segmentation(cath)->GetPadIxy(fX[i][cath],fY[i][cath]+dpy,ix,iy);
954 // iNN counts the number of neighbours with signal, it should be 1 or 2
955             Int_t iNN=0;
956             if (HitMap(cath)->TestHit(ix, iy)!=kEmpty) {
957                 iNN++;
958                 digt=(AliMUONDigit*) HitMap(cath)->GetHit(ix,iy);
959                 if (digt->fSignal > fQ[i][cath]) isLocal[i][cath]=kFALSE;
960             }
961             Segmentation(cath)->GetPadIxy(fX[i][cath],fY[i][cath]-dpy,ix,iy);
962             if (HitMap(cath)->TestHit(ix, iy)!=kEmpty) {
963                 iNN++;
964                 digt=(AliMUONDigit*) HitMap(cath)->GetHit(ix,iy);
965                 if (digt->fSignal > fQ[i][cath]) isLocal[i][cath]=kFALSE;
966             }
967             if (isLocal[i][cath] && iNN>0) {
968                 fIndLocal[fNLocal[cath]][cath]=i;
969                 fNLocal[cath]++;
970             } 
971         } // loop over all digits
972 // if one additional maximum has been found we are happy 
973 // if more maxima have been found restore the previous situation
974         fprintf(stderr,"\n New search gives %d local maxima for cathode 1 \n",fNLocal[0]);
975         fprintf(stderr,"                  %d local maxima for cathode 2 \n",fNLocal[1]);
976         if (fNLocal[cath]>2) {
977             fNLocal[cath]=iback;
978         }
979         
980     } // 1,2 local maxima
981     
982     if (fNLocal[0]==2 &&  (fNLocal[1]==1 || fNLocal[1]==0)) {
983         Int_t iback=fNLocal[1];
984         
985 //  Two local maxima on cathode 1 and one maximum on cathode 2 
986 //  Look for local maxima considering left and right neighbours on the 2nd cathode only
987         cath=1;
988         Int_t cath1=0;
989         
990
991 //
992 //  Loop over cluster digits
993         for (i=0; i<fMul[cath]; i++) {
994             isec=Segmentation(cath)->Sector(fIx[i][cath],fIy[i][cath]);
995             dpx=Segmentation(cath)->Dpx(isec);
996             dpy=Segmentation(cath)->Dpy(isec);
997             if (isLocal[i][cath]) continue;
998 // Pad position should be consistent with position of local maxima on the opposite cathode
999             if ((TMath::Abs(fY[i][cath]-fY[fIndLocal[0][cath1]][cath1]) > dpy/2.) && 
1000                 (TMath::Abs(fY[i][cath]-fY[fIndLocal[1][cath1]][cath1]) > dpy/2.))
1001                 continue;
1002 //
1003 // get neighbours for that digit and assume that it is local maximum        
1004             isLocal[i][cath]=kTRUE;
1005 // compare signal to that on the two neighbours on the left and on the right
1006             Segmentation(cath)->GetPadIxy(fX[i][cath]+dpx,fY[i][cath],ix,iy);
1007 // iNN counts the number of neighbours with signal, it should be 1 or 2
1008             Int_t iNN=0;
1009             if (HitMap(cath)->TestHit(ix, iy)!=kEmpty) {
1010                 iNN++;
1011                 digt=(AliMUONDigit*) HitMap(cath)->GetHit(ix,iy);
1012                 if (digt->fSignal > fQ[i][cath]) isLocal[i][cath]=kFALSE;
1013             }
1014             Segmentation(cath)->GetPadIxy(fX[i][cath]-dpx,fY[i][cath],ix,iy);
1015             if (HitMap(cath)->TestHit(ix, iy)!=kEmpty) {
1016                 iNN++;
1017                 digt=(AliMUONDigit*) HitMap(cath)->GetHit(ix,iy);
1018                 if (digt->fSignal > fQ[i][cath]) isLocal[i][cath]=kFALSE;
1019             }
1020             if (isLocal[i][cath] && iNN>0) {
1021                 fIndLocal[fNLocal[cath]][cath]=i;
1022                 fNLocal[cath]++;
1023             } 
1024         } // loop over all digits
1025 // if one additional maximum has been found we are happy 
1026 // if more maxima have been found restore the previous situation
1027         fprintf(stderr,"\n New search gives %d local maxima for cathode 1 \n",fNLocal[0]);
1028         fprintf(stderr,"\n                  %d local maxima for cathode 2 \n",fNLocal[1]);
1029 //      printf("\n New search gives %d %d \n",fNLocal[0],fNLocal[1]);
1030         if (fNLocal[cath]>2) {
1031             fNLocal[cath]=iback;
1032         }
1033
1034
1035
1036     } // 2,1 local maxima
1037 }
1038
1039
1040 void  AliMUONClusterFinderVS::FillCluster(AliMUONRawCluster* c, Int_t flag, Int_t cath) 
1041 {
1042 //
1043 //  Completes cluster information starting from list of digits
1044 //
1045     AliMUONDigit* dig;
1046     Float_t x, y;
1047     Int_t  ix, iy;
1048     
1049     if (cath==1) {
1050         c->fPeakSignal[cath]=c->fPeakSignal[0]; 
1051     } else {
1052         c->fPeakSignal[cath]=0;
1053     }
1054     
1055     
1056     if (flag) {
1057         c->fX[cath]=0;
1058         c->fY[cath]=0;
1059         c->fQ[cath]=0;
1060     }
1061
1062 //    fprintf(stderr,"\n fPeakSignal %d\n",c->fPeakSignal[cath]);
1063     for (Int_t i=0; i<c->fMultiplicity[cath]; i++)
1064     {
1065         dig= (AliMUONDigit*)Digits(cath)->UncheckedAt(c->fIndexMap[i][cath]);
1066         ix=dig->fPadX+c->fOffsetMap[i][cath];
1067         iy=dig->fPadY;
1068         Int_t q=dig->fSignal;
1069         if (!flag) q=Int_t(q*c->fContMap[i][cath]);
1070 //      fprintf(stderr,"q %d c->fPeakSignal[ %d ] %d\n",q,cath,c->fPeakSignal[cath]);
1071         if (dig->fPhysics >= dig->fSignal) {
1072             c->fPhysicsMap[i]=2;
1073         } else if (dig->fPhysics == 0) {
1074             c->fPhysicsMap[i]=0;
1075         } else  c->fPhysicsMap[i]=1;
1076 //
1077 // 
1078 //      fprintf(stderr,"q %d c->fPeakSignal[cath] %d\n",q,c->fPeakSignal[cath]);
1079 // peak signal and track list
1080         if (q>c->fPeakSignal[cath]) {
1081             c->fPeakSignal[cath]=q;
1082             c->fTracks[0]=dig->fHit;
1083             c->fTracks[1]=dig->fTracks[0];
1084             c->fTracks[2]=dig->fTracks[1];
1085 //          fprintf(stderr," c->fTracks[0] %d c->fTracks[1] %d\n",dig->fHit,dig->fTracks[0]);
1086         }
1087 //
1088         if (flag) {
1089             Segmentation(cath)->GetPadCxy(ix, iy, x, y);
1090             c->fX[cath] += q*x;
1091             c->fY[cath] += q*y;
1092             c->fQ[cath] += q;
1093         }
1094     } // loop over digits
1095 //    fprintf(stderr," fin du cluster c\n");
1096
1097
1098     if (flag) {
1099         c->fX[cath]/=c->fQ[cath];
1100         c->fX[cath]=Segmentation(cath)->GetAnod(c->fX[cath]);
1101         c->fY[cath]/=c->fQ[cath]; 
1102 //
1103 //  apply correction to the coordinate along the anode wire
1104 //
1105         x=c->fX[cath];   
1106         y=c->fY[cath];
1107         Segmentation(cath)->GetPadIxy(x, y, ix, iy);
1108         Segmentation(cath)->GetPadCxy(ix, iy, x, y);
1109         Int_t isec=Segmentation(cath)->Sector(ix,iy);
1110         TF1* cogCorr = Segmentation(cath)->CorrFunc(isec-1);
1111         
1112         if (cogCorr) {
1113             Float_t yOnPad=(c->fY[cath]-y)/Segmentation(cath)->Dpy(isec);
1114             c->fY[cath]=c->fY[cath]-cogCorr->Eval(yOnPad, 0, 0);
1115         }
1116     }
1117 }
1118
1119 void  AliMUONClusterFinderVS::FillCluster(AliMUONRawCluster* c, Int_t cath) 
1120 {
1121 //
1122 //  Completes cluster information starting from list of digits
1123 //
1124     static Float_t dr0;
1125
1126     AliMUONDigit* dig;
1127
1128     if (cath==0) {
1129         dr0 = 10000;
1130     }
1131     
1132     Float_t xpad, ypad;
1133     Float_t dx, dy, dr;
1134
1135     for (Int_t i=0; i<c->fMultiplicity[cath]; i++)
1136     {
1137       dig= (AliMUONDigit*)Digits(cath)->UncheckedAt(c->fIndexMap[i][cath]);
1138       Segmentation(cath)->
1139         GetPadCxy(dig->fPadX,dig->fPadY,xpad,ypad);
1140         fprintf(stderr,"x %f y %f cx %f cy %f\n",xpad,ypad,c->fX[0],c->fY[0]);
1141         dx = xpad - c->fX[0];
1142         dy = ypad - c->fY[0];
1143         dr = TMath::Sqrt(dx*dx+dy*dy);
1144
1145         if (dr < dr0) {
1146             dr0 = dr;
1147             fprintf(stderr," dr %f\n",dr);
1148             Int_t q=dig->fSignal;
1149             if (dig->fPhysics >= dig->fSignal) {
1150                 c->fPhysicsMap[i]=2;
1151             } else if (dig->fPhysics == 0) {
1152                 c->fPhysicsMap[i]=0;
1153             } else  c->fPhysicsMap[i]=1;
1154             c->fPeakSignal[cath]=q;
1155             c->fTracks[0]=dig->fHit;
1156             c->fTracks[1]=dig->fTracks[0];
1157             c->fTracks[2]=dig->fTracks[1];
1158             fprintf(stderr," c->fTracks[0] %d c->fTracks[1] %d\n",dig->fHit,dig->fTracks[0]);
1159         }
1160 //
1161     } // loop over digits
1162
1163 //  apply correction to the coordinate along the anode wire
1164     c->fX[cath]=Segmentation(cath)->GetAnod(c->fX[cath]);
1165 }
1166
1167 void  AliMUONClusterFinderVS::FindCluster(Int_t i, Int_t j, Int_t cath, AliMUONRawCluster &c){
1168 //
1169 //  Find clusterset
1170 //
1171 //
1172 //  Add i,j as element of the cluster
1173 //
1174
1175     Int_t idx = HitMap(cath)->GetHitIndex(i,j);
1176     AliMUONDigit* dig = (AliMUONDigit*) HitMap(cath)->GetHit(i,j);
1177     Int_t q=dig->fSignal;
1178     Int_t theX=dig->fPadX;
1179     Int_t theY=dig->fPadY;    
1180     if (q > TMath::Abs(c.fPeakSignal[0]) && q > TMath::Abs(c.fPeakSignal[1])) {
1181         c.fPeakSignal[cath]=q;
1182         c.fTracks[0]=dig->fHit;
1183         c.fTracks[1]=dig->fTracks[0];
1184         c.fTracks[2]=dig->fTracks[1];
1185     }
1186
1187 //
1188 //  Make sure that list of digits is ordered 
1189 // 
1190     Int_t mu=c.fMultiplicity[cath];
1191     c.fIndexMap[mu][cath]=idx;
1192     
1193     if (dig->fPhysics >= dig->fSignal) {
1194         c.fPhysicsMap[mu]=2;
1195     } else if (dig->fPhysics == 0) {
1196         c.fPhysicsMap[mu]=0;
1197     } else  c.fPhysicsMap[mu]=1;
1198     if (mu > 0) {
1199         for (Int_t ind=mu-1; ind>=0; ind--) {
1200             Int_t ist=(c.fIndexMap)[ind][cath];
1201             Int_t ql=((AliMUONDigit*)Digits(cath)
1202                       ->UncheckedAt(ist))->fSignal;
1203             Int_t ix=((AliMUONDigit*)Digits(cath)
1204                       ->UncheckedAt(ist))->fPadX;
1205             Int_t iy=((AliMUONDigit*)Digits(cath)
1206                       ->UncheckedAt(ist))->fPadY;
1207
1208             if (q>ql || (q==ql && theX > ix && theY < iy)) {
1209                 c.fIndexMap[ind][cath]=idx;
1210                 c.fIndexMap[ind+1][cath]=ist;
1211             } else {
1212                 break;
1213             }
1214         }
1215     }
1216     
1217     c.fMultiplicity[cath]++;
1218     if (c.fMultiplicity[cath] >= 50 ) {
1219         printf("FindCluster - multiplicity >50  %d \n",c.fMultiplicity[0]);
1220         c.fMultiplicity[cath]=49;
1221     }
1222
1223 // Prepare center of gravity calculation
1224     Float_t x, y;
1225     Segmentation(cath)->GetPadCxy(i, j, x, y);
1226             
1227     c.fX[cath] += q*x;
1228     c.fY[cath] += q*y;
1229     c.fQ[cath] += q;
1230 // Flag hit as taken  
1231     HitMap(cath)->FlagHit(i,j);
1232 //
1233 //  Now look recursively for all neighbours and pad hit on opposite cathode
1234 //
1235 //  Loop over neighbours
1236     Int_t ix,iy;
1237     Int_t nn;
1238     Int_t xList[kMaxNeighbours], yList[kMaxNeighbours];
1239     Segmentation(cath)->Neighbours(i,j,&nn,xList,yList);
1240     for (Int_t in=0; in<nn; in++) {
1241         ix=xList[in];
1242         iy=yList[in];
1243         if (HitMap(cath)->TestHit(ix,iy)==kUnused) FindCluster(ix, iy, cath, c);
1244     }
1245 //  Neighbours on opposite cathode 
1246 //  Take into account that several pads can overlap with the present pad
1247     Float_t xmin, xmax, ymin, ymax, xc, yc;
1248     Int_t iop;
1249     Int_t isec=Segmentation(cath)->Sector(i,j);    
1250     if (cath==0) {
1251         iop=1;
1252         xmin=x-Segmentation(cath)->Dpx(isec);
1253         xmax=x+Segmentation(cath)->Dpx(isec);
1254         xc=xmin+.001;
1255         while (xc < xmax) {
1256             xc+=Segmentation(iop)->Dpx(isec);
1257             Segmentation(iop)->GetPadIxy(xc,y,ix,iy);
1258             if (ix>=(Segmentation(iop)->Npx()) || (iy>=Segmentation(iop)->Npy())) continue;
1259             if (HitMap(iop)->TestHit(ix,iy)==kUnused) FindCluster(ix, iy, iop, c);
1260         }
1261     } else {
1262         iop=0;
1263         ymin=y-Segmentation(cath)->Dpy(isec);
1264         ymax=y+Segmentation(cath)->Dpy(isec);
1265         yc=ymin+.001;
1266         while (yc < ymax) {
1267             yc+=Segmentation(iop)->Dpy(isec);
1268             Segmentation(iop)->GetPadIxy(x,yc,ix,iy);
1269             if (ix>=(Segmentation(iop)->Npx()) || (iy>=Segmentation(iop)->Npy())) continue;
1270             if (HitMap(iop)->TestHit(ix,iy)==kUnused) FindCluster(ix, iy, iop, c);
1271         }
1272     }
1273 }
1274
1275 //_____________________________________________________________________________
1276
1277 void AliMUONClusterFinderVS::FindRawClusters()
1278 {
1279   //
1280   // MUON cluster finder from digits -- finds neighbours on both cathodes and 
1281   // fills the tree with raw clusters
1282   //
1283
1284     if (!NDigits(0) && !NDigits(1)) return;
1285
1286     fHitMap  = new AliMUONHitMapA1(fSegmentation[0] , fDigits);
1287     fHitMap2 = new AliMUONHitMapA1(fSegmentation[1], fDigits2);
1288
1289     AliMUONDigit *dig;
1290
1291     Int_t ndig, cath;
1292     Int_t nskip=0;
1293     Int_t ncls=0;
1294     HitMap(0)->FillHits();
1295     HitMap(1)->FillHits();
1296 //
1297 //  Outer Loop over Cathodes
1298     for (cath=0; cath<2; cath++) {
1299         for (ndig=0; ndig<NDigits(cath); ndig++) {
1300             dig = (AliMUONDigit*)Digits(cath)->UncheckedAt(ndig);
1301             Int_t i=dig->fPadX;
1302             Int_t j=dig->fPadY;
1303             if (HitMap(cath)->TestHit(i,j)==kUsed ||fHitMap->TestHit(i,j)==kEmpty) {
1304                 nskip++;
1305                 continue;
1306             }
1307             fprintf(stderr,"\n CATHODE %d CLUSTER %d\n",cath,ncls);
1308             AliMUONRawCluster c;
1309             c.fMultiplicity[0]=0;
1310             c.fMultiplicity[1]=0;
1311             c.fPeakSignal[cath]=dig->fSignal;
1312             c.fTracks[0]=dig->fHit;
1313             c.fTracks[1]=dig->fTracks[0];
1314             c.fTracks[2]=dig->fTracks[1];
1315             // tag the beginning of cluster list in a raw cluster
1316             c.fNcluster[0]=-1;
1317
1318             FindCluster(i,j,cath,c);
1319
1320             // center of gravity
1321             c.fX[0] /= c.fQ[0];
1322             c.fX[0]=Segmentation(0)->GetAnod(c.fX[0]);
1323             c.fY[0] /= c.fQ[0];
1324             c.fX[1] /= c.fQ[1];
1325             c.fX[1]=Segmentation(0)->GetAnod(c.fX[1]);
1326             c.fY[1] /= c.fQ[1];
1327             fprintf(stderr,"\n Cathode 1 multiplicite %d X(CG) %f Y(CG) %f\n",c.fMultiplicity[0],c.fX[0],c.fY[0]);
1328             fprintf(stderr," Cathode 2 multiplicite %d X(CG) %f Y(CG) %f\n",c.fMultiplicity[1],c.fX[1],c.fY[1]);
1329
1330 //     Mathieson Fit
1331 /*
1332             Bool_t fitted;
1333             
1334             fitted=SingleMathiesonFit(&c, 0);
1335             c.fX[0]=Segmentation(0)->GetAnod(c.fX[0]);
1336             fitted=SingleMathiesonFit(&c, 1);       
1337             c.fX[1]=Segmentation(1)->GetAnod(c.fX[1]);
1338 */ 
1339 //
1340 //      Analyse cluster and decluster if necessary
1341 //      
1342         ncls++;
1343         c.fNcluster[1]=fNRawClusters;
1344         c.fClusterType=c.PhysicsContribution();
1345
1346         fNPeaks=0;
1347 //
1348 //
1349         Decluster(&c);
1350 //      AddRawCluster(c);
1351
1352 //
1353 //      reset Cluster object
1354         { // begin local scope
1355             for (int k=0;k<c.fMultiplicity[0];k++) c.fIndexMap[k][0]=0;
1356         } // end local scope
1357
1358         { // begin local scope
1359             for (int k=0;k<c.fMultiplicity[1];k++) c.fIndexMap[k][1]=0;
1360         } // end local scope
1361         
1362         c.fMultiplicity[0]=c.fMultiplicity[0]=0;
1363
1364         
1365         } // end loop ndig
1366     } // end loop cathodes
1367     delete fHitMap;
1368     delete fHitMap2;
1369 }
1370
1371 Float_t AliMUONClusterFinderVS::SingleMathiesonFit(AliMUONRawCluster *c, Int_t cath)
1372 {
1373 //
1374     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1375     
1376     clusterInput.Fitter()->SetFCN(fcnS1);
1377     clusterInput.Fitter()->mninit(2,10,7);
1378     Double_t arglist[20];
1379     Int_t ierflag=0;
1380     arglist[0]=1;
1381 //      clusterInput.Fitter()->mnexcm("SET ERR",arglist,1,ierflag);
1382 // Set starting values 
1383     static Double_t vstart[2];
1384     vstart[0]=c->fX[1];
1385     vstart[1]=c->fY[0];
1386     
1387     
1388 // lower and upper limits
1389     static Double_t lower[2], upper[2];
1390     Int_t ix,iy;
1391     Segmentation(cath)->GetPadIxy(c->fX[cath], c->fY[cath], ix, iy);
1392     Int_t isec=Segmentation(cath)->Sector(ix, iy);
1393     lower[0]=vstart[0]-Segmentation(cath)->Dpx(isec)/2;
1394     lower[1]=vstart[1]-Segmentation(cath)->Dpy(isec)/2;
1395     
1396     upper[0]=lower[0]+Segmentation(cath)->Dpx(isec);
1397     upper[1]=lower[1]+Segmentation(cath)->Dpy(isec);
1398     
1399 // step sizes
1400     static Double_t step[2]={0.0005, 0.0005};
1401     
1402     clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1403     clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1404 // ready for minimisation       
1405     clusterInput.Fitter()->SetPrintLevel(1);
1406     clusterInput.Fitter()->mnexcm("SET OUT", arglist, 0, ierflag);
1407     arglist[0]= -1;
1408     arglist[1]= 0;
1409     
1410     clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1411     clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1412     clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1413     Double_t fmin, fedm, errdef;
1414     Int_t   npari, nparx, istat;
1415       
1416     clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);  
1417     fFitStat=istat;
1418     
1419 // Print results
1420 // Get fitted parameters
1421     Double_t xrec, yrec;
1422     TString chname;
1423     Double_t epxz, b1, b2;
1424     Int_t ierflg;
1425     clusterInput.Fitter()->mnpout(0, chname, xrec, epxz, b1, b2, ierflg);       
1426     clusterInput.Fitter()->mnpout(1, chname, yrec, epxz, b1, b2, ierflg);       
1427     fXFit[cath]=xrec;
1428     fYFit[cath]=yrec;
1429     return fmin;
1430 }
1431
1432 Float_t AliMUONClusterFinderVS::CombiSingleMathiesonFit(AliMUONRawCluster *c)
1433 {
1434 // Perform combined Mathieson fit on both cathode planes
1435 //
1436     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1437     clusterInput.Fitter()->SetFCN(fcnCombiS1);
1438     clusterInput.Fitter()->mninit(2,10,7);
1439     Double_t arglist[20];
1440     Int_t ierflag=0;
1441     arglist[0]=1;
1442     static Double_t vstart[2];
1443     vstart[0]=fXInit[0];
1444     vstart[1]=fYInit[0];
1445     
1446     
1447 // lower and upper limits
1448     static Double_t lower[2], upper[2];
1449     Int_t ix,iy,isec;
1450     Segmentation(0)->GetPadIxy(fXInit[0], fYInit[0], ix, iy);
1451     isec=Segmentation(0)->Sector(ix, iy);
1452     Float_t dpy=Segmentation(0)->Dpy(isec)/2;
1453     Segmentation(1)->GetPadIxy(fXInit[0], fYInit[0], ix, iy);
1454     isec=Segmentation(1)->Sector(ix, iy);
1455     Float_t dpx=Segmentation(1)->Dpx(isec)/2;
1456
1457
1458     lower[0]=vstart[0]-dpx;
1459     lower[1]=vstart[1]-dpy;
1460     
1461     upper[0]=vstart[0]+dpx;
1462     upper[1]=vstart[1]+dpy;
1463     
1464 // step sizes
1465     static Double_t step[2]={0.00001, 0.0001};
1466     
1467     clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1468     clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1469 // ready for minimisation       
1470     clusterInput.Fitter()->SetPrintLevel(1);
1471     clusterInput.Fitter()->mnexcm("SET OUT", arglist, 0, ierflag);
1472     arglist[0]= -1;
1473     arglist[1]= 0;
1474     
1475     clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1476     clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1477     clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1478     Double_t fmin, fedm, errdef;
1479     Int_t   npari, nparx, istat;
1480       
1481     clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);  
1482     fFitStat=istat;
1483     
1484 // Print results
1485 // Get fitted parameters
1486     Double_t xrec, yrec;
1487     TString chname;
1488     Double_t epxz, b1, b2;
1489     Int_t ierflg;
1490     clusterInput.Fitter()->mnpout(0, chname, xrec, epxz, b1, b2, ierflg);       
1491     clusterInput.Fitter()->mnpout(1, chname, yrec, epxz, b1, b2, ierflg);       
1492     fXFit[0]=xrec;
1493     fYFit[0]=yrec;
1494     return fmin;
1495 }
1496
1497 Bool_t AliMUONClusterFinderVS::DoubleMathiesonFit(AliMUONRawCluster *c, Int_t cath)
1498 {
1499 //
1500 //  Initialise global variables for fit
1501     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1502     clusterInput.Fitter()->SetFCN(fcnS2);
1503     clusterInput.Fitter()->mninit(5,10,7);
1504     Double_t arglist[20];
1505     Int_t ierflag=0;
1506     arglist[0]=1;
1507 // Set starting values 
1508     static Double_t vstart[5];
1509     vstart[0]=fX[fIndLocal[0][cath]][cath];
1510     vstart[1]=fY[fIndLocal[0][cath]][cath];     
1511     vstart[2]=fX[fIndLocal[1][cath]][cath];
1512     vstart[3]=fY[fIndLocal[1][cath]][cath];     
1513     vstart[4]=Float_t(fQ[fIndLocal[0][cath]][cath])/
1514         Float_t(fQ[fIndLocal[0][cath]][cath]+fQ[fIndLocal[1][cath]][cath]);
1515 // lower and upper limits
1516     static Double_t lower[5], upper[5];
1517     Int_t isec=Segmentation(cath)->Sector(fIx[fIndLocal[0][cath]][cath], fIy[fIndLocal[0][cath]][cath]);
1518     lower[0]=vstart[0]-Segmentation(cath)->Dpx(isec);
1519     lower[1]=vstart[1]-Segmentation(cath)->Dpy(isec);
1520     
1521     upper[0]=lower[0]+2.*Segmentation(cath)->Dpx(isec);
1522     upper[1]=lower[1]+2.*Segmentation(cath)->Dpy(isec);
1523     
1524     isec=Segmentation(cath)->Sector(fIx[fIndLocal[1][cath]][cath], fIy[fIndLocal[1][cath]][cath]);
1525     lower[2]=vstart[2]-Segmentation(cath)->Dpx(isec)/2;
1526     lower[3]=vstart[3]-Segmentation(cath)->Dpy(isec)/2;
1527     
1528     upper[2]=lower[2]+Segmentation(cath)->Dpx(isec);
1529     upper[3]=lower[3]+Segmentation(cath)->Dpy(isec);
1530     
1531     lower[4]=0.;
1532     upper[4]=1.;
1533 // step sizes
1534     static Double_t step[5]={0.0005, 0.0005, 0.0005, 0.0005, 0.0001};
1535     
1536     clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1537     clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1538     clusterInput.Fitter()->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
1539     clusterInput.Fitter()->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
1540     clusterInput.Fitter()->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
1541 // ready for minimisation       
1542     clusterInput.Fitter()->SetPrintLevel(-1);
1543     clusterInput.Fitter()->mnexcm("SET OUT", arglist, 0, ierflag);
1544     arglist[0]= -1;
1545     arglist[1]= 0;
1546     
1547     clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1548     clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1549     clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1550 // Get fitted parameters
1551     Double_t xrec[2], yrec[2], qfrac;
1552     TString chname;
1553     Double_t epxz, b1, b2;
1554     Int_t ierflg;
1555     clusterInput.Fitter()->mnpout(0, chname, xrec[0], epxz, b1, b2, ierflg);    
1556     clusterInput.Fitter()->mnpout(1, chname, yrec[0], epxz, b1, b2, ierflg);    
1557     clusterInput.Fitter()->mnpout(2, chname, xrec[1], epxz, b1, b2, ierflg);    
1558     clusterInput.Fitter()->mnpout(3, chname, yrec[1], epxz, b1, b2, ierflg);    
1559     clusterInput.Fitter()->mnpout(4, chname, qfrac,   epxz, b1, b2, ierflg);    
1560
1561     Double_t fmin, fedm, errdef;
1562     Int_t   npari, nparx, istat;
1563       
1564     clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);  
1565     fFitStat=istat;
1566     return kTRUE;
1567 }
1568
1569 Float_t AliMUONClusterFinderVS::CombiDoubleMathiesonFit(AliMUONRawCluster *c)
1570 {
1571 //
1572 // Perform combined double Mathieson fit on both cathode planes
1573 //
1574     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1575     clusterInput.Fitter()->SetFCN(fcnCombiS2);
1576     clusterInput.Fitter()->mninit(6,10,7);
1577     Double_t arglist[20];
1578     Int_t ierflag=0;
1579     arglist[0]=1;
1580 // Set starting values 
1581     static Double_t vstart[6];
1582     vstart[0]=fXInit[0];
1583     vstart[1]=fYInit[0];
1584     vstart[2]=fXInit[1];
1585     vstart[3]=fYInit[1];
1586     vstart[4]=fQrInit[0];
1587     vstart[5]=fQrInit[1];
1588 // lower and upper limits
1589     static Double_t lower[6], upper[6];
1590     Int_t ix,iy,isec;
1591     Float_t dpx, dpy;
1592     
1593     Segmentation(1)->GetPadIxy(fXInit[0], fYInit[0], ix, iy);
1594     isec=Segmentation(1)->Sector(ix, iy);
1595     dpx=Segmentation(1)->Dpx(isec);
1596
1597     Segmentation(0)->GetPadIxy(fXInit[0], fYInit[0], ix, iy);
1598     isec=Segmentation(0)->Sector(ix, iy);
1599     dpy=Segmentation(0)->Dpy(isec);
1600
1601     lower[0]=vstart[0]-dpx;
1602     lower[1]=vstart[1]-dpy;
1603     upper[0]=vstart[0]+dpx;
1604     upper[1]=vstart[1]+dpy;
1605
1606
1607     Segmentation(1)->GetPadIxy(fXInit[1], fYInit[1], ix, iy);
1608     isec=Segmentation(1)->Sector(ix, iy);
1609     dpx=Segmentation(1)->Dpx(isec);
1610     Segmentation(0)->GetPadIxy(fXInit[1], fYInit[1], ix, iy);
1611     isec=Segmentation(0)->Sector(ix, iy);
1612     dpy=Segmentation(0)->Dpy(isec);
1613
1614     lower[2]=vstart[2]-dpx;
1615     lower[3]=vstart[3]-dpy;
1616     upper[2]=vstart[2]+dpx;
1617     upper[3]=vstart[3]+dpy;
1618
1619
1620     lower[4]=0.;
1621     upper[4]=1.;
1622     lower[5]=0.;
1623     upper[5]=1.;
1624
1625 // step sizes
1626     static Double_t step[6]={0.0005, 0.0005, 0.0005, 0.0005, 0.001, 0.001};
1627     
1628     clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1629     clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1630     clusterInput.Fitter()->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
1631     clusterInput.Fitter()->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
1632     clusterInput.Fitter()->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
1633     clusterInput.Fitter()->mnparm(5,"a1",vstart[5],step[5],lower[5],upper[5],ierflag);
1634 // ready for minimisation       
1635     clusterInput.Fitter()->SetPrintLevel(-1);
1636     clusterInput.Fitter()->mnexcm("SET OUT", arglist, 0, ierflag);
1637     arglist[0]= -1;
1638     arglist[1]= 0;
1639     
1640     clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1641     clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1642     clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1643 // Get fitted parameters
1644     TString chname;
1645     Double_t epxz, b1, b2;
1646     Int_t ierflg;
1647     clusterInput.Fitter()->mnpout(0, chname, fXFit[0],  epxz, b1, b2, ierflg);  
1648     clusterInput.Fitter()->mnpout(1, chname, fYFit[0],  epxz, b1, b2, ierflg);  
1649     clusterInput.Fitter()->mnpout(2, chname, fXFit[1],  epxz, b1, b2, ierflg);  
1650     clusterInput.Fitter()->mnpout(3, chname, fYFit[1],  epxz, b1, b2, ierflg);  
1651     clusterInput.Fitter()->mnpout(4, chname, fQrFit[0], epxz, b1, b2, ierflg);  
1652     clusterInput.Fitter()->mnpout(5, chname, fQrFit[1], epxz, b1, b2, ierflg);  
1653
1654     Double_t fmin, fedm, errdef;
1655     Int_t   npari, nparx, istat;
1656       
1657     clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);  
1658     fFitStat=istat;
1659     
1660     fChi2[0]=fmin;
1661     fChi2[1]=fmin;
1662     return fmin;
1663 }
1664
1665 void AliMUONClusterFinderVS::Split(AliMUONRawCluster* c)
1666 {
1667 //
1668 // One cluster for each maximum
1669 //
1670     Int_t i, j, cath;
1671     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1672     for (j=0; j<2; j++) {
1673         AliMUONRawCluster cnew;
1674         for (cath=0; cath<2; cath++) {
1675             cnew.fChi2[cath]=fChi2[0];
1676             
1677             if (fNPeaks == 0) {
1678                 cnew.fNcluster[0]=-1;
1679                 cnew.fNcluster[1]=fNRawClusters;
1680             } else {
1681                 cnew.fNcluster[0]=fNPeaks;
1682                 cnew.fNcluster[1]=0;
1683             }
1684             cnew.fMultiplicity[cath]=0;
1685             cnew.fX[cath]=Float_t(fXFit[j]);
1686             cnew.fY[cath]=Float_t(fYFit[j]);
1687             if (j==0) {
1688                 cnew.fQ[cath]=Int_t(clusterInput.TotalCharge(cath)*fQrFit[cath]);
1689             } else {
1690                 cnew.fQ[cath]=Int_t(clusterInput.TotalCharge(cath)*(1-fQrFit[cath]));
1691             }
1692             fSegmentation[cath]->SetHit(fXFit[j],fYFit[j]);
1693             for (i=0; i<fMul[cath]; i++) {
1694                 cnew.fIndexMap[cnew.fMultiplicity[cath]][cath]=
1695                     c->fIndexMap[i][cath];
1696                 fSegmentation[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
1697                 Float_t q1=fResponse->IntXY(fSegmentation[cath]);
1698                 cnew.fContMap[i][cath]
1699                     =(q1*Float_t(cnew.fQ[cath]))/Float_t(fQ[i][cath]);
1700                 cnew.fMultiplicity[cath]++;
1701 //              printf(" fXFIT %f fYFIT %f Multiplicite %d\n",cnew.fX[cath],cnew.fY[cath],cnew.fMultiplicity[cath]);
1702             }
1703             FillCluster(&cnew,0,cath);
1704         } // cathode loop
1705         
1706         cnew.fClusterType=cnew.PhysicsContribution();
1707         if (cnew.fQ[0]>0 && cnew.fQ[1]>0) AddRawCluster(cnew);
1708         fNPeaks++;
1709     }
1710 }
1711
1712
1713 //
1714 // Minimisation functions
1715 // Single Mathieson
1716 void fcnS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
1717 {
1718     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());    
1719     Int_t i;
1720     Float_t delta;
1721     Float_t chisq=0;
1722     Float_t qcont=0;
1723     Float_t qtot=0;
1724
1725     for (i=0; i<clusterInput.Nmul(0); i++) {
1726         Float_t q0=clusterInput.Charge(i,0);
1727         Float_t q1=clusterInput.DiscrChargeS1(i,par);
1728         delta=(q0-q1)/q0;
1729         chisq+=delta*delta;
1730         qcont+=q1;
1731         qtot+=q0;
1732     }
1733     f=chisq;
1734 }
1735
1736 void fcnCombiS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
1737 {
1738     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());    
1739     Int_t i, cath;
1740     Float_t delta;
1741     Float_t chisq=0;
1742     Float_t qcont=0;
1743     Float_t qtot=0;
1744     //    Float_t chi2temp=0;
1745
1746     for (cath=0; cath<2; cath++) {
1747 //      chisq=0;
1748         for (i=0; i<clusterInput.Nmul(cath); i++) {
1749             Float_t q0=clusterInput.Charge(i,cath);
1750             Float_t q1=clusterInput.DiscrChargeCombiS1(i,par,cath);
1751             //      delta=(q0-q1);
1752             delta=(q0-q1)/q0;
1753             chisq+=delta*delta;
1754             qcont+=q1;
1755             qtot+=q0;
1756         }
1757 //      if (cath == 0) chi2temp=chisq/clusterInput.Nbins[cath];
1758     }
1759 //    chisq = chisq/clusterInput.Nbins[1]+chi2temp; 
1760     f=chisq;
1761 }
1762
1763 // Double Mathieson
1764 void fcnS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
1765 {
1766     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());    
1767     Int_t i;
1768     Float_t delta;
1769     Float_t chisq=0;
1770     Float_t qcont=0;
1771     Float_t qtot=0;
1772     
1773     for (i=0; i<clusterInput.Nmul(0); i++) {
1774
1775         Float_t q0=clusterInput.Charge(i,0);
1776         Float_t q1=clusterInput.DiscrChargeS2(i,par);
1777         delta=(q0-q1)/q0;
1778         chisq+=delta*delta;
1779         qcont+=q1;
1780         qtot+=q0;
1781     }
1782 //    chisq=chisq+=(qtot-qcont)*(qtot-qcont)*0.5;
1783     f=chisq;
1784 }
1785
1786 // Double Mathieson
1787 void fcnCombiS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
1788 {
1789     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());    
1790     Int_t i, cath;
1791     Float_t delta;
1792     Float_t chisq=0;
1793     Float_t qcont=0;
1794     Float_t qtot=0;
1795     //    Float_t chi2temp=0;
1796
1797     for (cath=0; cath<2; cath++) {
1798 //      chisq=0;
1799         for (i=0; i<clusterInput.Nmul(cath); i++) {
1800             Float_t q0=clusterInput.Charge(i,cath);
1801             Float_t q1=clusterInput.DiscrChargeCombiS2(i,par,cath);
1802             //      delta=(q0-q1);
1803             delta=(q0-q1)/q0;
1804             chisq+=delta*delta;
1805             qcont+=q1;
1806             qtot+=q0;
1807         }
1808 //      if (cath == 0) chi2temp=chisq/clusterInput.Nbins[cath];
1809     }
1810 //    chisq = chisq/clusterInput.Nbins[1]+chi2temp;     
1811     f=chisq;
1812 }
1813
1814 void AliMUONClusterFinderVS::AddRawCluster(const AliMUONRawCluster c)
1815 {
1816   //
1817   // Add a raw cluster copy to the list
1818   //
1819     AliMUON *pMUON=(AliMUON*)gAlice->GetModule("MUON");
1820     pMUON->AddRawCluster(fChamber,c); 
1821     fNRawClusters++;
1822     fprintf(stderr,"\nfNRawClusters %d\n",fNRawClusters);
1823 }
1824
1825
1826 AliMUONClusterFinderVS& AliMUONClusterFinderVS
1827 ::operator = (const AliMUONClusterFinderVS& rhs)
1828 {
1829 // Dummy assignment operator
1830     return *this;
1831 }
1832
1833
1834
1835
1836
1837
1838
1839
1840