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