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