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