1e5701e3e0cf338420f50ebf8888e7c2f9f947c2
[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 /* $Id$ */
17
18 #include <TMinuit.h> 
19 #include <TF1.h>
20 #include <TMinuit.h> 
21 #include <Riostream.h>
22
23 #include "AliMUONClusterFinderVS.h"
24 #include "AliMUONDigit.h"
25 #include "AliMUONRawCluster.h"
26 #include "AliSegmentation.h"
27 #include "AliMUONGeometrySegmentation.h"
28 #include "AliMUONMathieson.h"
29 #include "AliMUONClusterInput.h"
30 #include "AliMUONHitMapA1.h"
31 #include "AliLog.h"
32
33 //_____________________________________________________________________
34 // This function is minimized in the double-Mathieson fit
35 void fcnS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
36 void fcnS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
37 void fcnCombiS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
38 void fcnCombiS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
39
40 ClassImp(AliMUONClusterFinderVS)
41
42 AliMUONClusterFinderVS::AliMUONClusterFinderVS()
43   : TObject()
44 {
45 // Default constructor
46     fInput=AliMUONClusterInput::Instance();
47 //     cout <<  " TYPE" << fSegmentationType << endl;
48     fHitMap[0] = 0;
49     fHitMap[1] = 0;
50     fTrack[0]=fTrack[1]=-1;
51     fDebugLevel = 0; // make silent default
52     fGhostChi2Cut = 1e6; // nothing done by default
53     fSeg[0]    = 0;
54     fSeg[1]    = 0;
55     fSeg2[0]    = 0;
56     fSeg2[1]    = 0;
57
58     for(Int_t i=0; i<100; i++) {
59       for (Int_t j=0; j<2; j++) {
60         fDig[i][j] = 0;
61       }
62     } 
63     fRawClusters = new TClonesArray("AliMUONRawCluster",1000);
64     fNRawClusters = 0;
65 }
66  //____________________________________________________________________________
67 AliMUONClusterFinderVS::~AliMUONClusterFinderVS()
68 {
69   // Reset tracks information
70    fNRawClusters = 0;
71    if (fRawClusters) {
72      fRawClusters->Delete();
73      delete fRawClusters;
74    }
75 }
76
77 AliMUONClusterFinderVS::AliMUONClusterFinderVS(const AliMUONClusterFinderVS & clusterFinder):TObject(clusterFinder)
78 {
79 // Protected copy constructor
80
81   AliFatal("Not implemented.");
82 }
83 //____________________________________________________________________________
84 void AliMUONClusterFinderVS::ResetRawClusters()
85 {
86   // Reset tracks information
87   fNRawClusters = 0;
88   if (fRawClusters) fRawClusters->Clear();
89 }
90 //____________________________________________________________________________
91 void AliMUONClusterFinderVS::Decluster(AliMUONRawCluster *cluster)
92 {
93 // Decluster by local maxima
94     SplitByLocalMaxima(cluster);
95 }
96 //____________________________________________________________________________
97 void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
98 {
99 // Split complex cluster by local maxima 
100     Int_t cath, i;
101
102     fInput->SetCluster(c);
103
104     fMul[0]=c->GetMultiplicity(0);
105     fMul[1]=c->GetMultiplicity(1);
106
107 //
108 //  dump digit information into arrays
109 //
110
111     Float_t qtot;
112     
113     for (cath=0; cath<2; cath++) {
114         qtot=0;
115         for (i=0; i<fMul[cath]; i++)
116         {
117             // pointer to digit
118             fDig[i][cath]=fInput->Digit(cath, c->GetIndex(i, cath));
119             // pad coordinates
120             fIx[i][cath]= fDig[i][cath]->PadX();
121             fIy[i][cath]= fDig[i][cath]->PadY();
122             // pad charge
123             fQ[i][cath] = fDig[i][cath]->Signal();
124             // pad centre coordinates
125             if (fSegmentationType == 1)
126               fSeg[cath]->
127                 GetPadC(fIx[i][cath], fIy[i][cath], fX[i][cath], fY[i][cath], fZ[i][cath]);
128             else
129               fSeg2[cath]->
130                 GetPadC(fInput->DetElemId(), fIx[i][cath], fIy[i][cath], fX[i][cath], fY[i][cath], fZ[i][cath]);
131         } // loop over cluster digits
132     }  // loop over cathodes
133
134
135     FindLocalMaxima(c);
136
137 //
138 //  Initialise and perform mathieson fits
139     Float_t chi2, oldchi2;
140 //  ++++++++++++++++++*************+++++++++++++++++++++
141 //  (1) No more than one local maximum per cathode plane 
142 //  +++++++++++++++++++++++++++++++*************++++++++
143     if ((fNLocal[0]==1 && (fNLocal[1]==0 ||  fNLocal[1]==1)) || 
144         (fNLocal[0]==0 && fNLocal[1]==1)) {
145 // Perform combined single Mathieson fit
146 // Initial values for coordinates (x,y) 
147
148         // One local maximum on cathodes 1 and 2 (X->cathode 2, Y->cathode 1)
149         if (fNLocal[0]==1 &&  fNLocal[1]==1) {
150             fXInit[0]=c->GetX(1);
151             fYInit[0]=c->GetY(0);
152             // One local maximum on cathode 1 (X,Y->cathode 1)
153         } else if (fNLocal[0]==1) {
154             fXInit[0]=c->GetX(0);
155             fYInit[0]=c->GetY(0);
156             // One local maximum on cathode 2  (X,Y->cathode 2)
157         } else {
158             fXInit[0]=c->GetX(1);
159             fYInit[0]=c->GetY(1);
160         }
161         AliDebug(1,"cas (1) CombiSingleMathiesonFit(c)");
162         chi2=CombiSingleMathiesonFit(c);
163 //      Int_t ndf = fgNbins[0]+fgNbins[1]-2;
164 //      Float_t prob = TMath::Prob(Double_t(chi2),ndf);
165 //      prob1->Fill(prob);
166 //      chi2_1->Fill(chi2);
167         oldchi2=chi2;
168         AliDebug(1,Form(" chi2 %f ",chi2));        
169
170         c->SetX(0, fXFit[0]);
171         c->SetY(0, fYFit[0]);
172
173         c->SetX(1,fXFit[0]);
174         c->SetY(1,fYFit[0]);
175         c->SetChi2(0,chi2);
176         c->SetChi2(1,chi2);
177         // Force on anod
178         if (fSegmentationType == 1) {
179           c->SetX(0, fSeg[0]->GetAnod(c->GetX(0)));
180           c->SetX(1, fSeg[1]->GetAnod(c->GetX(1)));
181         } else {
182           c->SetX(0, fSeg2[0]->GetAnod(fInput->DetElemId(), c->GetX(0)));
183           c->SetX(1, fSeg2[1]->GetAnod(fInput->DetElemId(), c->GetX(1)));
184         }
185         
186 // If reasonable chi^2 add result to the list of rawclusters
187         if (chi2 < 0.3) {
188             AddRawCluster(*c);
189 // If not try combined double Mathieson Fit
190         } else {
191                 AliDebug(1," MAUVAIS CHI2 !!!\n");
192             if (fNLocal[0]==1 &&  fNLocal[1]==1) {
193                 fXInit[0]=fX[fIndLocal[0][1]][1];
194                 fYInit[0]=fY[fIndLocal[0][0]][0];
195                 fXInit[1]=fX[fIndLocal[0][1]][1];
196                 fYInit[1]=fY[fIndLocal[0][0]][0];
197             } else if (fNLocal[0]==1) {
198                 fXInit[0]=fX[fIndLocal[0][0]][0];
199                 fYInit[0]=fY[fIndLocal[0][0]][0];
200                 fXInit[1]=fX[fIndLocal[0][0]][0];
201                 fYInit[1]=fY[fIndLocal[0][0]][0];
202             } else {
203                 fXInit[0]=fX[fIndLocal[0][1]][1];
204                 fYInit[0]=fY[fIndLocal[0][1]][1];
205                 fXInit[1]=fX[fIndLocal[0][1]][1];
206                 fYInit[1]=fY[fIndLocal[0][1]][1];
207             }
208             
209 //  Initial value for charge ratios
210             fQrInit[0]=0.5;
211             fQrInit[1]=0.5;
212             AliDebug(1,"\n cas (1) CombiDoubleMathiesonFit(c)\n");
213             chi2=CombiDoubleMathiesonFit(c);
214 //          Int_t ndf = fgNbins[0]+fgNbins[1]-6;
215 //          Float_t prob = TMath::Prob(chi2,ndf);
216 //          prob2->Fill(prob);
217 //          chi2_2->Fill(chi2);
218             
219 // Was this any better ??
220             AliDebug(1,Form(" Old and new chi2 %f %f ", oldchi2, chi2));
221             if (fFitStat!=0 && chi2>0 && (2.*chi2 < oldchi2)) {
222               AliDebug(1,"Split");
223                 // Split cluster into two according to fit result
224                 Split(c);
225             } else {
226               AliDebug(1,"Do not Split");
227                 // Don't split
228                 AddRawCluster(*c);
229             }
230         }
231
232 //  +++++++++++++++++++++++++++++++++++++++
233 //  (2) Two local maxima per cathode plane 
234 //  +++++++++++++++++++++++++++++++++++++++
235     } else if (fNLocal[0]==2 &&  fNLocal[1]==2) {
236 //
237 //  Let's look for ghosts first 
238
239         Float_t xm[4][2], ym[4][2];
240         Float_t dpx, dpy, dx, dy;
241         Int_t ixm[4][2], iym[4][2];
242         Int_t isec, im1, im2, ico;
243 //
244 //  Form the 2x2 combinations
245 //  0-0, 0-1, 1-0, 1-1  
246         ico=0;
247         for (im1=0; im1<2; im1++) {
248             for (im2=0; im2<2; im2++) {     
249                 xm[ico][0]=fX[fIndLocal[im1][0]][0];
250                 ym[ico][0]=fY[fIndLocal[im1][0]][0];
251                 xm[ico][1]=fX[fIndLocal[im2][1]][1];
252                 ym[ico][1]=fY[fIndLocal[im2][1]][1];
253
254                 ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
255                 iym[ico][0]=fIy[fIndLocal[im1][0]][0];
256                 ixm[ico][1]=fIx[fIndLocal[im2][1]][1];
257                 iym[ico][1]=fIy[fIndLocal[im2][1]][1];
258                 ico++;
259             }
260         }
261 // ico = 0 : first local maximum on cathodes 1 and 2
262 // ico = 1 : fisrt local maximum on cathode 1 and second on cathode 2
263 // ico = 2 : second local maximum on cathode 1 and first on cathode 1
264 // ico = 3 : second local maximum on cathodes 1 and 2
265
266 // Analyse the combinations and keep those that are possible !
267 // For each combination check consistency in x and y    
268         Int_t   iacc;
269         Bool_t  accepted[4];
270         Float_t dr[4] = {1.e4, 1.e4, 1.e4, 1.e4};
271         iacc=0;
272
273 // In case of staggering maxima are displaced by exactly half the pad-size in y. 
274 // We have to take into account the numerical precision in the consistency check;       
275         Float_t eps = 1.e-5;
276 //
277         for (ico=0; ico<4; ico++) {
278             accepted[ico]=kFALSE;
279 // cathode one: x-coordinate
280             if (fSegmentationType == 1) {
281               isec=fSeg[0]->Sector(ixm[ico][0], iym[ico][0]);
282               dpx=fSeg[0]->Dpx(isec)/2.;
283             } else {
284               isec=fSeg2[0]->Sector(fInput->DetElemId(), ixm[ico][0], iym[ico][0]);
285               dpx=fSeg2[0]->Dpx(fInput->DetElemId(), isec)/2.;
286             }
287             dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
288 // cathode two: y-coordinate
289             if (fSegmentationType == 1) {
290               isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
291               dpy=fSeg[1]->Dpy(isec)/2.;
292             } else {
293               isec=fSeg2[1]->Sector(fInput->DetElemId(), ixm[ico][1], iym[ico][1]);
294               dpy=fSeg2[1]->Dpy(fInput->DetElemId(), isec)/2.;
295             }
296             dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
297             AliDebug(2,Form("\n %i %f %f %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy, dx, dpx ));
298             if ((dx <= dpx) && (dy <= dpy+eps)) {
299                 // consistent
300                 accepted[ico]=kTRUE;
301                 dr[ico] = TMath::Sqrt(dx*dx+dy*dy);
302                 iacc++;
303             } else {
304                 // reject
305                 accepted[ico]=kFALSE;
306             }
307         }
308         AliDebug(1,Form("\n iacc= %d:\n", iacc));
309         if (iacc == 3) {
310             if (accepted[0] && accepted[1]) {
311                 if (dr[0] >= dr[1]) {
312                     accepted[0]=kFALSE;
313                 } else {
314                     accepted[1]=kFALSE;
315                 }
316             }
317
318             if (accepted[2] && accepted[3]) {
319                 if (dr[2] >= dr[3]) {
320                     accepted[2]=kFALSE;
321                 } else {
322                     accepted[3]=kFALSE;
323                 }
324             }
325 /*          
326 // eliminate one candidate
327             Float_t drmax = 0;
328             Int_t icobad = -1;
329
330             for (ico=0; ico<4; ico++) {
331                 if (accepted[ico] && dr[ico] > drmax) {
332                     icobad = ico;
333                     drmax  = dr[ico];
334                 }
335             }
336             
337             accepted[icobad] = kFALSE;
338 */
339             iacc = 2;
340         }
341         
342         
343         AliDebug(1,Form("\n iacc= %d:\n", iacc));
344         if (iacc==2) {
345                 AliDebug(1,"\n iacc=2: No problem ! \n");
346         } else if (iacc==4) {
347                 AliDebug(1,"\n iacc=4: Ok, but ghost problem !!! \n");
348         } else if (iacc==0) {
349                 AliDebug(1,"\n iacc=0: I don't know what to do with this !!!!!!!!! \n");
350         }
351
352 //  Initial value for charge ratios
353         fQrInit[0]=Float_t(fQ[fIndLocal[0][0]][0])/
354             Float_t(fQ[fIndLocal[0][0]][0]+fQ[fIndLocal[1][0]][0]);
355         fQrInit[1]=Float_t(fQ[fIndLocal[0][1]][1])/
356             Float_t(fQ[fIndLocal[0][1]][1]+fQ[fIndLocal[1][1]][1]);
357         
358 // ******* iacc = 0 *******
359 // No combinations found between the 2 cathodes
360 // We keep the center of gravity of the cluster
361         if (iacc==0) {
362             AddRawCluster(*c);
363         }
364
365 // ******* iacc = 1 *******
366 // Only one combination found between the 2 cathodes
367         if (iacc==1) {
368 // Initial values for the 2 maxima (x,y)
369
370 // 1 maximum is initialised with the maximum of the combination found (X->cathode 2, Y->cathode 1)
371 // 1 maximum is initialised with the other maximum of the first cathode  
372             if (accepted[0]){
373                 AliDebug(1,"ico=0");
374                 fXInit[0]=xm[0][1];
375                 fYInit[0]=ym[0][0];
376                 fXInit[1]=xm[3][0];
377                 fYInit[1]=ym[3][0];
378             } else if (accepted[1]){
379                 AliDebug(1,"ico=1");
380                 fXInit[0]=xm[1][1];
381                 fYInit[0]=ym[1][0];
382                 fXInit[1]=xm[2][0];
383                 fYInit[1]=ym[2][0];
384             } else if (accepted[2]){
385                 AliDebug(1,"ico=2");
386                 fXInit[0]=xm[2][1];
387                 fYInit[0]=ym[2][0];
388                 fXInit[1]=xm[1][0];
389                 fYInit[1]=ym[1][0];
390             } else if (accepted[3]){
391                 AliDebug(1,"ico=3");
392                 fXInit[0]=xm[3][1];
393                 fYInit[0]=ym[3][0];
394                 fXInit[1]=xm[0][0];
395                 fYInit[1]=ym[0][0];
396             }
397                 AliDebug(1,"cas (2) CombiDoubleMathiesonFit(c)");
398             chi2=CombiDoubleMathiesonFit(c);
399 //          Int_t ndf = fgNbins[0]+fgNbins[1]-6;
400 //          Float_t prob = TMath::Prob(chi2,ndf);
401 //          prob2->Fill(prob);
402 //          chi2_2->Fill(chi2);
403             AliDebug(1,Form(" chi2 %f\n",chi2));
404
405 // If reasonable chi^2 add result to the list of rawclusters
406             if (chi2<10) {
407                 Split(c);
408
409             } else {
410 // 1 maximum is initialised with the maximum of the combination found (X->cathode 2, Y->cathode 1)
411 // 1 maximum is initialised with the other maximum of the second cathode  
412                 if (accepted[0]){
413                         AliDebug(1,"ico=0");
414                     fXInit[0]=xm[0][1];
415                     fYInit[0]=ym[0][0];
416                     fXInit[1]=xm[3][1];
417                     fYInit[1]=ym[3][1];
418                 } else if (accepted[1]){
419                         AliDebug(1,"ico=1");
420                     fXInit[0]=xm[1][1];
421                     fYInit[0]=ym[1][0];
422                     fXInit[1]=xm[2][1];
423                     fYInit[1]=ym[2][1];
424                 } else if (accepted[2]){
425                         AliDebug(1,"ico=2");
426                     fXInit[0]=xm[2][1];
427                     fYInit[0]=ym[2][0];
428                     fXInit[1]=xm[1][1];
429                     fYInit[1]=ym[1][1];
430                 } else if (accepted[3]){
431                         AliDebug(1,"ico=3");
432                     fXInit[0]=xm[3][1];
433                     fYInit[0]=ym[3][0];
434                     fXInit[1]=xm[0][1];
435                     fYInit[1]=ym[0][1];
436                 }
437                 AliDebug(1,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
438                 chi2=CombiDoubleMathiesonFit(c);
439 //              Int_t ndf = fgNbins[0]+fgNbins[1]-6;
440 //              Float_t prob = TMath::Prob(chi2,ndf);
441 //              prob2->Fill(prob);
442 //              chi2_2->Fill(chi2);
443                 AliDebug(1,Form(" chi2 %f\n",chi2));
444
445 // If reasonable chi^2 add result to the list of rawclusters
446                 if (chi2<10) {
447                     Split(c);
448                 } else {
449 //We keep only the combination found (X->cathode 2, Y->cathode 1)
450                     for (Int_t ico=0; ico<2; ico++) {
451                         if (accepted[ico]) {
452                             AliMUONRawCluster cnew;
453                             Int_t cath;    
454                             for (cath=0; cath<2; cath++) {
455                                 cnew.SetX(cath, Float_t(xm[ico][1]));
456                                 cnew.SetY(cath, Float_t(ym[ico][0]));
457                                 cnew.SetZ(cath, fZPlane);
458                                 
459                                 cnew.SetMultiplicity(cath,c->GetMultiplicity(cath));
460                                 for (i=0; i<fMul[cath]; i++) {
461                                   cnew.SetIndex(i, cath, c->GetIndex(i,cath));
462                                   if (fSegmentationType == 1) 
463                                     fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
464                                   else 
465                                     fSeg2[cath]->SetPad(fInput->DetElemId(), fIx[i][cath], fIy[i][cath]);
466                                 }
467                                 AliDebug(1,Form("\nRawCluster %d cath %d\n",ico,cath));
468                                 AliDebug(1,Form("mult_av %d\n",c->GetMultiplicity(cath)));
469                                 FillCluster(&cnew,cath);
470                             } 
471                             cnew.SetClusterType(cnew.PhysicsContribution());
472                             AddRawCluster(cnew);
473                             fNPeaks++;
474                         }
475                     }
476                 }
477             }
478         }
479         
480 // ******* iacc = 2 *******
481 // Two combinations found between the 2 cathodes
482         if (iacc==2) {
483 // Was the same maximum taken twice
484             if ((accepted[0]&&accepted[1]) || (accepted[2]&&accepted[3])) {
485                 AliDebug(1,"\n Maximum taken twice !!!\n");
486
487 // Have a try !! with that
488                 if (accepted[0]&&accepted[3]) {
489                     fXInit[0]=xm[0][1];
490                     fYInit[0]=ym[0][0];
491                     fXInit[1]=xm[1][1];
492                     fYInit[1]=ym[1][0];
493                 } else {
494                     fXInit[0]=xm[2][1];
495                     fYInit[0]=ym[2][0];
496                     fXInit[1]=xm[3][1];
497                     fYInit[1]=ym[3][0];
498                 }
499                 AliDebug(1,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
500                 chi2=CombiDoubleMathiesonFit(c);
501 //                  Int_t ndf = fgNbins[0]+fgNbins[1]-6;
502 //                  Float_t prob = TMath::Prob(chi2,ndf);
503 //                  prob2->Fill(prob);
504 //                  chi2_2->Fill(chi2);
505                 Split(c);
506                 
507             } else {
508 // No ghosts ! No Problems ! -  Perform one fit only !
509                 if (accepted[0]&&accepted[3]) {
510                     fXInit[0]=xm[0][1];
511                     fYInit[0]=ym[0][0];
512                     fXInit[1]=xm[3][1];
513                     fYInit[1]=ym[3][0];
514                 } else {
515                     fXInit[0]=xm[1][1];
516                     fYInit[0]=ym[1][0];
517                     fXInit[1]=xm[2][1];
518                     fYInit[1]=ym[2][0];
519                 }
520                 AliDebug(1,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
521                 chi2=CombiDoubleMathiesonFit(c);
522 //                  Int_t ndf = fgNbins[0]+fgNbins[1]-6;
523 //                  Float_t prob = TMath::Prob(chi2,ndf);
524 //                  prob2->Fill(prob);
525 //                  chi2_2->Fill(chi2);
526                 AliDebug(1,Form(" chi2 %f\n",chi2));
527                 Split(c);
528             }
529             
530 // ******* iacc = 4 *******
531 // Four combinations found between the 2 cathodes
532 // Ghost !!
533         } else if (iacc==4) {
534 // Perform fits for the two possibilities !!    
535 // Accept if charges are compatible on both cathodes
536 // If none are compatible, keep everything
537             fXInit[0]=xm[0][1];
538             fYInit[0]=ym[0][0];
539             fXInit[1]=xm[3][1];
540             fYInit[1]=ym[3][0];
541             AliDebug(1,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
542             chi2=CombiDoubleMathiesonFit(c);
543 //              Int_t ndf = fgNbins[0]+fgNbins[1]-6;
544 //              Float_t prob = TMath::Prob(chi2,ndf);
545 //              prob2->Fill(prob);
546 //              chi2_2->Fill(chi2);
547             AliDebug(1,Form(" chi2 %f\n",chi2));
548             // store results of fit and postpone decision
549             Double_t sXFit[2],sYFit[2],sQrFit[2];
550             Float_t sChi2[2];
551             for (Int_t i=0;i<2;i++) {
552                 sXFit[i]=fXFit[i];
553                 sYFit[i]=fYFit[i];
554                 sQrFit[i]=fQrFit[i];
555                 sChi2[i]=fChi2[i];
556             }
557             fXInit[0]=xm[1][1];
558             fYInit[0]=ym[1][0];
559             fXInit[1]=xm[2][1];
560             fYInit[1]=ym[2][0];
561             AliDebug(1,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
562             chi2=CombiDoubleMathiesonFit(c);
563 //              ndf = fgNbins[0]+fgNbins[1]-6;
564 //              prob = TMath::Prob(chi2,ndf);
565 //              prob2->Fill(prob);
566 //              chi2_2->Fill(chi2);
567             AliDebug(1,Form(" chi2 %f\n",chi2));
568             // We have all informations to perform the decision
569             // Compute the chi2 for the 2 possibilities
570             Float_t chi2fi,chi2si,chi2f,chi2s;
571
572             chi2f = (TMath::Log(fInput->TotalCharge(0)*fQrFit[0]
573                   /  (fInput->TotalCharge(1)*fQrFit[1]) )
574                   / fInput->ChargeCorrel() );
575             chi2f *=chi2f;
576             chi2fi = (TMath::Log(fInput->TotalCharge(0)*(1-fQrFit[0])
577                   /  (fInput->TotalCharge(1)*(1-fQrFit[1])) )
578                   / fInput->ChargeCorrel() );
579             chi2f += chi2fi*chi2fi;
580
581             chi2s = (TMath::Log(fInput->TotalCharge(0)*sQrFit[0]
582                   /  (fInput->TotalCharge(1)*sQrFit[1]) )
583                   / fInput->ChargeCorrel() );
584             chi2s *=chi2s;
585             chi2si = (TMath::Log(fInput->TotalCharge(0)*(1-sQrFit[0])
586                   /  (fInput->TotalCharge(1)*(1-sQrFit[1])) )
587                   / fInput->ChargeCorrel() );
588             chi2s += chi2si*chi2si;
589
590             // usefull to store the charge matching chi2 in the cluster
591             // fChi2[0]=sChi2[1]=chi2f;
592             // fChi2[1]=sChi2[0]=chi2s;
593
594             if (chi2f<=fGhostChi2Cut && chi2s<=fGhostChi2Cut)
595                 c->SetGhost(1);
596             if   (chi2f>fGhostChi2Cut && chi2s>fGhostChi2Cut) {
597                 // we keep the ghost
598                 c->SetGhost(2);
599                 chi2s=-1;
600                 chi2f=-1;
601             }
602             if (chi2f<=fGhostChi2Cut)
603                 Split(c);
604             if (chi2s<=fGhostChi2Cut) {
605                 // retreive saved values
606                 for (Int_t i=0;i<2;i++) {
607                     fXFit[i]=sXFit[i];
608                     fYFit[i]=sYFit[i];
609                     fQrFit[i]=sQrFit[i];
610                     fChi2[i]=sChi2[i];
611                 }
612                 Split(c);
613             }
614             c->SetGhost(0);
615         }
616
617     } else if (fNLocal[0]==2 &&  fNLocal[1]==1) {
618 //  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
619 //  (3) Two local maxima on cathode 1 and one maximum on cathode 2 
620 //  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
621 //
622         Float_t xm[4][2], ym[4][2];
623         Float_t dpx, dpy, dx, dy;
624         Int_t ixm[4][2], iym[4][2];
625         Int_t isec, im1, ico;
626 //
627 //  Form the 2x2 combinations
628 //  0-0, 0-1, 1-0, 1-1  
629         ico=0;
630         for (im1=0; im1<2; im1++) {
631             xm[ico][0]=fX[fIndLocal[im1][0]][0];
632             ym[ico][0]=fY[fIndLocal[im1][0]][0];
633             xm[ico][1]=fX[fIndLocal[0][1]][1];
634             ym[ico][1]=fY[fIndLocal[0][1]][1];
635             
636             ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
637             iym[ico][0]=fIy[fIndLocal[im1][0]][0];
638             ixm[ico][1]=fIx[fIndLocal[0][1]][1];
639             iym[ico][1]=fIy[fIndLocal[0][1]][1];
640             ico++;
641         }
642 // ico = 0 : first local maximum on cathodes 1 and 2
643 // ico = 1 : second local maximum on cathode 1 and first on cathode 2
644
645 // Analyse the combinations and keep those that are possible !
646 // For each combination check consistency in x and y    
647         Int_t iacc;
648         Bool_t accepted[4];
649         iacc=0;
650         // In case of staggering maxima are displaced by exactly half the pad-size in y. 
651         // We have to take into account the numerical precision in the consistency check;
652         
653         Float_t eps = 1.e-5;
654
655         for (ico=0; ico<2; ico++) {
656             accepted[ico]=kFALSE;
657             if (fSegmentationType == 1) {
658               isec=fSeg[0]->Sector(ixm[ico][0], iym[ico][0]);
659               dpx=fSeg[0]->Dpx(isec)/2.;
660             } else {
661               isec=fSeg2[0]->Sector(fInput->DetElemId(), ixm[ico][0], iym[ico][0]);
662               dpx=fSeg2[0]->Dpx(fInput->DetElemId(), isec)/2.;
663             }
664             dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
665             if (fSegmentationType == 1) {
666               isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
667               dpy=fSeg[1]->Dpy(isec)/2.;
668             } else {
669               isec=fSeg2[1]->Sector(fInput->DetElemId(), ixm[ico][1], iym[ico][1]);
670               dpy=fSeg2[1]->Dpy(fInput->DetElemId(), isec)/2.;
671             }
672             dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
673             AliDebug(2,Form("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy ));
674             if ((dx <= dpx) && (dy <= dpy+eps)) {
675                 // consistent
676                 accepted[ico]=kTRUE;
677                 iacc++;
678             } else {
679                 // reject
680                 accepted[ico]=kFALSE;
681             }
682         }
683         
684         Float_t chi21 = 100;
685         Float_t chi22 = 100;
686         Float_t chi23 = 100;
687
688         //  Initial value for charge ratios
689         fQrInit[0]=Float_t(fQ[fIndLocal[0][0]][0])/
690             Float_t(fQ[fIndLocal[0][0]][0]+fQ[fIndLocal[1][0]][0]);
691         fQrInit[1]=fQrInit[0];
692         
693         if (accepted[0] && accepted[1]) {
694             
695             fXInit[0]=0.5*(xm[0][1]+xm[0][0]);
696             fYInit[0]=ym[0][0];
697             fXInit[1]=0.5*(xm[0][1]+xm[1][0]);
698             fYInit[1]=ym[1][0];
699             fQrInit[0]=0.5;
700             fQrInit[1]=0.5;
701             chi23=CombiDoubleMathiesonFit(c);
702             if (chi23<10) {
703                 Split(c);
704                 Float_t yst;
705                 yst = fYFit[0];
706                 fYFit[0] = fYFit[1];
707                 fYFit[1] = yst;
708                 Split(c);
709             }
710         } else if (accepted[0]) {
711             fXInit[0]=xm[0][1];
712             fYInit[0]=ym[0][0];
713             fXInit[1]=xm[1][0];
714             fYInit[1]=ym[1][0];
715             chi21=CombiDoubleMathiesonFit(c);
716 //          Int_t ndf = fgNbins[0]+fgNbins[1]-6;
717 //          Float_t prob = TMath::Prob(chi2,ndf);
718 //          prob2->Fill(prob);
719 //          chi2_2->Fill(chi21);
720             AliDebug(1,Form(" chi2 %f\n",chi21));
721             if (chi21<10) Split(c);
722         } else if (accepted[1]) {
723             fXInit[0]=xm[1][1];
724             fYInit[0]=ym[1][0];
725             fXInit[1]=xm[0][0];
726             fYInit[1]=ym[0][0];
727             chi22=CombiDoubleMathiesonFit(c);
728 //          Int_t ndf = fgNbins[0]+fgNbins[1]-6;
729 //          Float_t prob = TMath::Prob(chi2,ndf);
730 //          prob2->Fill(prob);
731 //          chi2_2->Fill(chi22);
732             AliDebug(1,Form(" chi2 %f\n",chi22));
733             if (chi22<10) Split(c);
734         }
735
736         if (chi21 > 10 && chi22 > 10 && chi23 > 10) {
737 // We keep only the combination found (X->cathode 2, Y->cathode 1)
738             for (Int_t ico=0; ico<2; ico++) {
739                 if (accepted[ico]) {
740                     AliMUONRawCluster cnew;
741                     Int_t cath;    
742                     for (cath=0; cath<2; cath++) {
743                         cnew.SetX(cath, Float_t(xm[ico][1]));
744                         cnew.SetY(cath, Float_t(ym[ico][0]));
745                         cnew.SetZ(cath, fZPlane);
746                         cnew.SetMultiplicity(cath, c->GetMultiplicity(cath));
747                         for (i=0; i<fMul[cath]; i++) {
748                             cnew.SetIndex(i, cath, c->GetIndex(i, cath));
749                             if (fSegmentationType == 1) 
750                               fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
751                             else
752                               fSeg2[cath]->SetPad(fInput->DetElemId(), fIx[i][cath], fIy[i][cath]);
753
754                         }
755                         AliDebug(1,Form("\nRawCluster %d cath %d\n",ico,cath));
756                         AliDebug(1,Form("mult_av %d\n",c->GetMultiplicity(cath)));
757                         
758                         FillCluster(&cnew,cath);
759                     } 
760                     cnew.SetClusterType(cnew.PhysicsContribution());
761                     AddRawCluster(cnew);
762                     fNPeaks++;
763                 }
764             }
765         }
766         
767 //  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
768 //  (3') One local maximum on cathode 1 and two maxima on cathode 2 
769 //  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
770     } else if (fNLocal[0]==1 && fNLocal[1]==2) {
771         Float_t xm[4][2], ym[4][2];
772         Float_t dpx, dpy, dx, dy;
773         Int_t ixm[4][2], iym[4][2];
774         Int_t isec, im1, ico;
775 //
776 //  Form the 2x2 combinations
777 //  0-0, 0-1, 1-0, 1-1  
778         ico=0;
779         for (im1=0; im1<2; im1++) {
780             xm[ico][0]=fX[fIndLocal[0][0]][0];
781             ym[ico][0]=fY[fIndLocal[0][0]][0];
782             xm[ico][1]=fX[fIndLocal[im1][1]][1];
783             ym[ico][1]=fY[fIndLocal[im1][1]][1];
784             
785             ixm[ico][0]=fIx[fIndLocal[0][0]][0];
786             iym[ico][0]=fIy[fIndLocal[0][0]][0];
787             ixm[ico][1]=fIx[fIndLocal[im1][1]][1];
788             iym[ico][1]=fIy[fIndLocal[im1][1]][1];
789             ico++;
790         }
791 // ico = 0 : first local maximum on cathodes 1 and 2
792 // ico = 1 : first local maximum on cathode 1 and second on cathode 2
793
794 // Analyse the combinations and keep those that are possible !
795 // For each combination check consistency in x and y    
796         Int_t iacc;
797         Bool_t accepted[4];
798         iacc=0;
799         // In case of staggering maxima are displaced by exactly half the pad-size in y. 
800         // We have to take into account the numerical precision in the consistency check;       
801         Float_t eps = 1.e-5;
802
803         
804         for (ico=0; ico<2; ico++) {
805             accepted[ico]=kFALSE;
806             if (fSegmentationType == 1) {
807               isec=fSeg[0]->Sector(ixm[ico][0], iym[ico][0]);
808               dpx=fSeg[0]->Dpx(isec)/2.;
809             } else {
810               isec=fSeg2[0]->Sector(fInput->DetElemId(), ixm[ico][0], iym[ico][0]);
811               dpx=fSeg2[0]->Dpx(fInput->DetElemId(), isec)/2.;
812             }
813             dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
814             if (fSegmentationType == 1) {
815               isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
816               dpy=fSeg[1]->Dpy(isec)/2.;
817             } else {
818               isec=fSeg2[1]->Sector(fInput->DetElemId(), ixm[ico][1], iym[ico][1]);
819               dpy=fSeg2[1]->Dpy(fInput->DetElemId(), isec)/2.;
820             }
821             dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
822             AliDebug(1,Form("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy ));
823             if ((dx <= dpx) && (dy <= dpy+eps)) {
824                 // consistent
825                 accepted[ico]=kTRUE;
826                 AliDebug(1,Form("ico %d\n",ico));
827                 iacc++;
828             } else {
829                 // reject
830                 accepted[ico]=kFALSE;
831             }
832         }
833
834         Float_t chi21 = 100;
835         Float_t chi22 = 100;
836         Float_t chi23 = 100;
837
838         fQrInit[1]=Float_t(fQ[fIndLocal[0][1]][1])/
839             Float_t(fQ[fIndLocal[0][1]][1]+fQ[fIndLocal[1][1]][1]);
840         
841         fQrInit[0]=fQrInit[1];
842
843         
844         if (accepted[0] && accepted[1]) {
845             fXInit[0]=xm[0][1];
846             fYInit[0]=0.5*(ym[0][0]+ym[0][1]);
847             fXInit[1]=xm[1][1];
848             fYInit[1]=0.5*(ym[0][0]+ym[1][1]);
849             fQrInit[0]=0.5;
850             fQrInit[1]=0.5;
851             chi23=CombiDoubleMathiesonFit(c);
852             if (chi23<10) {
853                 Split(c);
854                 Float_t yst;
855                 yst = fYFit[0];
856                 fYFit[0] = fYFit[1];
857                 fYFit[1] = yst;
858                 Split(c);
859             }
860         } else if (accepted[0]) {
861             fXInit[0]=xm[0][0];
862             fYInit[0]=ym[0][1];
863             fXInit[1]=xm[1][1];
864             fYInit[1]=ym[1][1];
865             chi21=CombiDoubleMathiesonFit(c);
866 //          Int_t ndf = fgNbins[0]+fgNbins[1]-6;
867 //          Float_t prob = TMath::Prob(chi2,ndf);
868 //          prob2->Fill(prob);
869 //          chi2_2->Fill(chi21);
870             AliDebug(1,Form(" chi2 %f\n",chi21));
871             if (chi21<10) Split(c);
872         } else if (accepted[1]) {
873             fXInit[0]=xm[1][0];
874             fYInit[0]=ym[1][1];
875             fXInit[1]=xm[0][1];
876             fYInit[1]=ym[0][1];
877             chi22=CombiDoubleMathiesonFit(c);
878 //          Int_t ndf = fgNbins[0]+fgNbins[1]-6;
879 //          Float_t prob = TMath::Prob(chi2,ndf);
880 //          prob2->Fill(prob);
881 //          chi2_2->Fill(chi22);
882             AliDebug(1,Form(" chi2 %f\n",chi22));
883             if (chi22<10) Split(c);
884         }
885
886         if (chi21 > 10 && chi22 > 10 && chi23 > 10) {
887 //We keep only the combination found (X->cathode 2, Y->cathode 1)
888             for (Int_t ico=0; ico<2; ico++) {
889                 if (accepted[ico]) {
890                     AliMUONRawCluster cnew;
891                     Int_t cath;    
892                     for (cath=0; cath<2; cath++) {
893                         cnew.SetX(cath, Float_t(xm[ico][1]));
894                         cnew.SetY(cath, Float_t(ym[ico][0]));
895                         cnew.SetZ(cath, fZPlane);
896                         cnew.SetMultiplicity(cath, c->GetMultiplicity(cath));
897                         for (i=0; i<fMul[cath]; i++) {
898                             cnew.SetIndex(i, cath, c->GetIndex(i, cath));
899                             if (fSegmentationType == 1) 
900                               fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
901                             else 
902                               fSeg2[cath]->SetPad(fInput->DetElemId(), fIx[i][cath], fIy[i][cath]);
903                         }
904                         AliDebug(1,Form("\nRawCluster %d cath %d\n",ico,cath));
905                         AliDebug(1,Form("mult_av %d\n",c->GetMultiplicity(cath)));
906                         FillCluster(&cnew,cath);
907                     } 
908                     cnew.SetClusterType(cnew.PhysicsContribution());
909                     AddRawCluster(cnew);
910                     fNPeaks++;
911                 }
912             }
913         }
914
915 //  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
916 //  (4) At least three local maxima on cathode 1 or on cathode 2 
917 //  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
918     } else if (fNLocal[0]>2 || fNLocal[1]>2) {
919         Int_t param = fNLocal[0]*fNLocal[1];
920         Int_t ii;
921
922         Float_t ** xm = new Float_t * [param];
923         for (ii=0; ii<param; ii++) xm[ii]=new Float_t [2];
924         Float_t ** ym = new Float_t * [param];
925         for (ii=0; ii<param; ii++) ym[ii]=new Float_t [2];
926         Int_t ** ixm = new Int_t * [param];
927         for (ii=0; ii<param; ii++) ixm[ii]=new Int_t [2];
928         Int_t ** iym = new Int_t * [param];
929         for (ii=0; ii<param; ii++) iym[ii]=new Int_t [2];
930         
931         Int_t isec, ico;
932         Float_t dpx, dpy, dx, dy;
933
934         ico=0;
935         for (Int_t im1=0; im1<fNLocal[0]; im1++) {
936             for (Int_t im2=0; im2<fNLocal[1]; im2++) {
937                 xm[ico][0]=fX[fIndLocal[im1][0]][0];
938                 ym[ico][0]=fY[fIndLocal[im1][0]][0];
939                 xm[ico][1]=fX[fIndLocal[im2][1]][1];
940                 ym[ico][1]=fY[fIndLocal[im2][1]][1];
941
942                 ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
943                 iym[ico][0]=fIy[fIndLocal[im1][0]][0];
944                 ixm[ico][1]=fIx[fIndLocal[im2][1]][1];
945                 iym[ico][1]=fIy[fIndLocal[im2][1]][1];
946                 ico++;
947             }
948         }
949         
950         Int_t nIco = ico;
951         AliDebug(1,Form("nIco %d\n",nIco));
952         for (ico=0; ico<nIco; ico++) {
953             AliDebug(1,Form("ico = %d\n",ico));
954             if (fSegmentationType == 1) {
955               isec=fSeg[0]->Sector(ixm[ico][0], iym[ico][0]);
956               dpx=fSeg[0]->Dpx(isec)/2.;
957             } else {
958               isec=fSeg2[0]->Sector(fInput->DetElemId(), ixm[ico][0], iym[ico][0]);
959               dpx=fSeg2[0]->Dpx(fInput->DetElemId(), isec)/2.;
960             } 
961             dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
962             if (fSegmentationType == 1) {
963               isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
964               dpy=fSeg[1]->Dpy(isec)/2.;
965             } else {
966               isec=fSeg2[1]->Sector(fInput->DetElemId(), ixm[ico][1], iym[ico][1]);
967               dpy=fSeg2[1]->Dpy(fInput->DetElemId(), isec)/2.;
968             }
969             dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
970                 AliDebug(1,Form("dx %f dpx %f dy %f dpy %f\n",dx,dpx,dy,dpy));
971                 AliDebug(1,Form("  X %f Y %f\n",xm[ico][1],ym[ico][0]));
972             if ((dx <= dpx) && (dy <= dpy)) {
973                         AliDebug(1,"ok\n");
974                 Int_t cath;    
975                 AliMUONRawCluster cnew;
976                 for (cath=0; cath<2; cath++) {
977                     cnew.SetX(cath, Float_t(xm[ico][1]));
978                     cnew.SetY(cath, Float_t(ym[ico][0]));
979                     cnew.SetZ(cath, fZPlane);
980                     cnew.SetMultiplicity(cath, c->GetMultiplicity(cath));
981                     for (i=0; i<fMul[cath]; i++) {
982                         cnew.SetIndex(i, cath, c->GetIndex(i, cath));
983                         if (fSegmentationType == 1) 
984                           fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
985                         else 
986                           fSeg2[cath]->SetPad(fInput->DetElemId(), fIx[i][cath], fIy[i][cath]);
987                     }
988                     FillCluster(&cnew,cath);
989                 } 
990                 cnew.SetClusterType(cnew.PhysicsContribution());
991                 AddRawCluster(cnew);
992                 fNPeaks++;
993             }
994         }
995         delete [] xm;
996         delete [] ym;
997         delete [] ixm;
998         delete [] iym;
999     }
1000 }
1001
1002 void AliMUONClusterFinderVS::FindLocalMaxima(AliMUONRawCluster* /*c*/)
1003 {
1004 // Find all local maxima of a cluster
1005     AliDebug(1,"\n Find Local maxima  !");
1006     
1007     AliMUONDigit* digt;
1008     
1009     Int_t cath, cath1; // loops over cathodes
1010     Int_t i;           // loops over digits
1011     Int_t j;           // loops over cathodes
1012 //
1013 //  Find local maxima
1014 //
1015 //  counters for number of local maxima
1016     fNLocal[0]=fNLocal[1]=0;
1017 //  flags digits as local maximum
1018     Bool_t isLocal[100][2];
1019     for (i=0; i<100;i++) {
1020         isLocal[i][0]=isLocal[i][1]=kFALSE;
1021     }
1022 //  number of next neighbours and arrays to store them 
1023     Int_t nn;
1024     Int_t x[10], y[10];
1025     // loop over cathodes
1026     for (cath=0; cath<2; cath++) {
1027       // loop over cluster digits
1028       for (i=0; i<fMul[cath]; i++) {
1029         // get neighbours for that digit and assume that it is local maximum
1030         Int_t isec;
1031         Float_t a0;
1032
1033         if (fSegmentationType == 1)         
1034           fSeg[cath]->Neighbours(fIx[i][cath], fIy[i][cath], &nn, x, y);
1035         else 
1036           fSeg2[cath]->Neighbours(fInput->DetElemId(), fIx[i][cath], fIy[i][cath], &nn, x, y);
1037           
1038         isLocal[i][cath]=kTRUE;
1039
1040         if (fSegmentationType == 1) {
1041           isec = fSeg[cath]->Sector(fIx[i][cath], fIy[i][cath]);
1042           a0   = fSeg[cath]->Dpx(isec)*fSeg[cath]->Dpy(isec);
1043         } else {
1044           isec = fSeg2[cath]->Sector(fInput->DetElemId(), fIx[i][cath], fIy[i][cath]);
1045           a0   = fSeg2[cath]->Dpx(fInput->DetElemId(), isec)*fSeg2[cath]->Dpy(fInput->DetElemId(), isec);
1046         }
1047         // loop over next neighbours, if at least one neighbour has higher charger assumption
1048         // digit is not local maximum 
1049         for (j=0; j<nn; j++) {
1050           if (fHitMap[cath]->TestHit(x[j], y[j])==kEmpty) continue;
1051           digt=(AliMUONDigit*) fHitMap[cath]->GetHit(x[j], y[j]);
1052           Float_t a1;
1053           if (fSegmentationType == 1) {
1054             isec=fSeg[cath]->Sector(x[j], y[j]);
1055             a1 = fSeg[cath]->Dpx(isec)*fSeg[cath]->Dpy(isec);
1056           } else {
1057             isec=fSeg2[cath]->Sector(fInput->DetElemId(), x[j], y[j]);
1058             a1 = fSeg2[cath]->Dpx(fInput->DetElemId(),isec)*fSeg2[cath]->Dpy(fInput->DetElemId(), isec);
1059           }
1060           if (digt->Signal()/a1 > fQ[i][cath]/a0) {
1061             isLocal[i][cath]=kFALSE;
1062             break;
1063             //
1064             // handle special case of neighbouring pads with equal signal
1065           } else if (digt->Signal() == fQ[i][cath]) {
1066             if (fNLocal[cath]>0) {
1067               for (Int_t k=0; k<fNLocal[cath]; k++) {
1068                 if (x[j]==fIx[fIndLocal[k][cath]][cath] 
1069                     && y[j]==fIy[fIndLocal[k][cath]][cath])
1070                   {
1071                     isLocal[i][cath]=kFALSE;
1072                   } 
1073               } // loop over local maxima
1074             } // are there already local maxima
1075           } // same charge ? 
1076         } // loop over next neighbours
1077         if (isLocal[i][cath]) {
1078           fIndLocal[fNLocal[cath]][cath]=i;
1079           fNLocal[cath]++;
1080         } 
1081       } // loop over all digits
1082     } // loop over cathodes
1083
1084     AliDebug(1,Form("\n Found %d %d %d %d local Maxima\n",
1085                fNLocal[0], fNLocal[1], fMul[0], fMul[1]));
1086         AliDebug(1,Form("\n Cathode 1 local Maxima %d Multiplicite %d\n",fNLocal[0], fMul[0]));
1087         AliDebug(1,Form(" Cathode 2 local Maxima %d Multiplicite %d\n",fNLocal[1], fMul[1]));
1088     Int_t ix, iy, isec;
1089     Float_t dpx, dpy;
1090     
1091     
1092     if (fNLocal[1]==2 &&  (fNLocal[0]==1 || fNLocal[0]==0)) {
1093         Int_t iback=fNLocal[0];
1094         
1095 //  Two local maxima on cathode 2 and one maximum on cathode 1 
1096 //  Look for local maxima considering up and down neighbours on the 1st cathode only
1097 //
1098 //  Loop over cluster digits
1099         cath=0;
1100         cath1=1;
1101         
1102         for (i=0; i<fMul[cath]; i++) {
1103           if (fSegmentationType == 1) {
1104             isec=fSeg[cath]->Sector(fIx[i][cath],fIy[i][cath]);
1105             dpy=fSeg[cath]->Dpy(isec);
1106             dpx=fSeg[cath]->Dpx(isec);
1107           } else {
1108             isec=fSeg2[cath]->Sector(fInput->DetElemId(), fIx[i][cath],fIy[i][cath]);
1109             dpy=fSeg2[cath]->Dpy(fInput->DetElemId(), isec);
1110             dpx=fSeg2[cath]->Dpx(fInput->DetElemId(), isec);
1111           }
1112             if (isLocal[i][cath]) continue;
1113 // Pad position should be consistent with position of local maxima on the opposite cathode
1114             if ((TMath::Abs(fX[i][cath]-fX[fIndLocal[0][cath1]][cath1]) > dpx/2.) && 
1115                 (TMath::Abs(fX[i][cath]-fX[fIndLocal[1][cath1]][cath1]) > dpx/2.))
1116                 continue;
1117
1118 // get neighbours for that digit and assume that it is local maximum        
1119             isLocal[i][cath]=kTRUE;
1120 // compare signal to that on the two neighbours on the left and on the right
1121 // iNN counts the number of neighbours with signal, it should be 1 or 2
1122             Int_t iNN=0;
1123             if (fSegmentationType == 1) {
1124
1125               for (fSeg[cath]->FirstPad(fX[i][cath], fY[i][cath], fZPlane, 0., dpy);
1126                    fSeg[cath]->MorePads();
1127                    fSeg[cath]->NextPad())
1128                 {
1129                   ix = fSeg[cath]->Ix();
1130                   iy = fSeg[cath]->Iy();
1131                   // skip the current pad
1132                   if (iy == fIy[i][cath]) continue;
1133                 
1134                   if (fHitMap[cath]->TestHit(ix, iy)!=kEmpty) {
1135                     iNN++;
1136                     digt=(AliMUONDigit*) fHitMap[cath]->GetHit(ix,iy);
1137                     if (digt->Signal() > fQ[i][cath]) isLocal[i][cath]=kFALSE;
1138                   }
1139                 } // Loop over pad neighbours in y
1140             } else {
1141
1142               for (fSeg2[cath]->FirstPad(fInput->DetElemId(), fX[i][cath], fY[i][cath], fZPlane, 0., dpy);
1143                    fSeg2[cath]->MorePads(fInput->DetElemId());
1144                    fSeg2[cath]->NextPad(fInput->DetElemId()))
1145                 {
1146                   ix = fSeg2[cath]->Ix();
1147                   iy = fSeg2[cath]->Iy();
1148                   // skip the current pad
1149                   if (iy == fIy[i][cath]) continue;
1150                 
1151                   if (fHitMap[cath]->TestHit(ix, iy)!=kEmpty) {
1152                     iNN++;
1153                     digt=(AliMUONDigit*) fHitMap[cath]->GetHit(ix,iy);
1154                     if (digt->Signal() > fQ[i][cath]) isLocal[i][cath]=kFALSE;
1155                   }
1156                 } // Loop over pad neighbours in y
1157             }
1158             if (isLocal[i][cath] && iNN>0) {
1159                 fIndLocal[fNLocal[cath]][cath]=i;
1160                 fNLocal[cath]++;
1161             } 
1162         } // loop over all digits
1163 // if one additional maximum has been found we are happy 
1164 // if more maxima have been found restore the previous situation
1165         AliDebug(1,Form("\n New search gives %d local maxima for cathode 1 \n",
1166                     fNLocal[0]));
1167         AliDebug(1,Form("                  %d local maxima for cathode 2 \n",
1168                     fNLocal[1]));
1169         if (fNLocal[cath]>2) {
1170             fNLocal[cath]=iback;
1171         }
1172         
1173     } // 1,2 local maxima
1174     
1175     if (fNLocal[0]==2 &&  (fNLocal[1]==1 || fNLocal[1]==0)) {
1176         Int_t iback=fNLocal[1];
1177         
1178 //  Two local maxima on cathode 1 and one maximum on cathode 2 
1179 //  Look for local maxima considering left and right neighbours on the 2nd cathode only
1180         cath=1;
1181         Int_t cath1 = 0;
1182         Float_t eps = 1.e-5;
1183         
1184 //
1185 //  Loop over cluster digits
1186         for (i=0; i<fMul[cath]; i++) {
1187           if (fSegmentationType == 1) {
1188             isec=fSeg[cath]->Sector(fIx[i][cath],fIy[i][cath]);
1189             dpx=fSeg[cath]->Dpx(isec);
1190             dpy=fSeg[cath]->Dpy(isec);
1191           } else {
1192             isec=fSeg2[cath]->Sector(fInput->DetElemId(), fIx[i][cath],fIy[i][cath]);
1193             dpx=fSeg2[cath]->Dpx(fInput->DetElemId(), isec);
1194             dpy=fSeg2[cath]->Dpy(fInput->DetElemId(), isec);
1195           }
1196         
1197             if (isLocal[i][cath]) continue;
1198 // Pad position should be consistent with position of local maxima on the opposite cathode
1199             if ((TMath::Abs(fY[i][cath]-fY[fIndLocal[0][cath1]][cath1]) > dpy/2.+eps) && 
1200                 (TMath::Abs(fY[i][cath]-fY[fIndLocal[1][cath1]][cath1]) > dpy/2.+eps))
1201                 continue;
1202             
1203 //
1204 // get neighbours for that digit and assume that it is local maximum        
1205             isLocal[i][cath]=kTRUE;
1206 // compare signal to that on the two neighbours on the left and on the right
1207
1208 // iNN counts the number of neighbours with signal, it should be 1 or 2
1209             Int_t iNN=0;
1210             if (fSegmentationType == 1) {
1211               for (fSeg[cath]->FirstPad(fX[i][cath], fY[i][cath], fZPlane, dpx, 0.);
1212                    fSeg[cath]->MorePads();
1213                    fSeg[cath]->NextPad())
1214                 {
1215
1216                   ix = fSeg[cath]->Ix();
1217                   iy = fSeg[cath]->Iy();
1218
1219                   // skip the current pad
1220                   if (ix == fIx[i][cath]) continue;
1221                 
1222                   if (fHitMap[cath]->TestHit(ix, iy)!=kEmpty) {
1223                     iNN++;
1224                     digt=(AliMUONDigit*) fHitMap[cath]->GetHit(ix,iy);
1225                     if (digt->Signal() > fQ[i][cath]) isLocal[i][cath]=kFALSE;
1226                   }
1227                 } // Loop over pad neighbours in x
1228             } else {
1229               for (fSeg2[cath]->FirstPad(fInput->DetElemId(), fX[i][cath], fY[i][cath], fZPlane, dpx, 0.);
1230                    fSeg2[cath]->MorePads(fInput->DetElemId());
1231                    fSeg2[cath]->NextPad(fInput->DetElemId()))
1232                 {
1233
1234                   ix = fSeg2[cath]->Ix();
1235                   iy = fSeg2[cath]->Iy();
1236
1237                   // skip the current pad
1238                   if (ix == fIx[i][cath]) continue;
1239                 
1240                   if (fHitMap[cath]->TestHit(ix, iy)!=kEmpty) {
1241                     iNN++;
1242                     digt=(AliMUONDigit*) fHitMap[cath]->GetHit(ix,iy);
1243                     if (digt->Signal() > fQ[i][cath]) isLocal[i][cath]=kFALSE;
1244                   }
1245                 } // Loop over pad neighbours in x
1246             }
1247             if (isLocal[i][cath] && iNN>0) {
1248                 fIndLocal[fNLocal[cath]][cath]=i;
1249                 fNLocal[cath]++;
1250             } 
1251         } // loop over all digits
1252 // if one additional maximum has been found we are happy 
1253 // if more maxima have been found restore the previous situation
1254         AliDebug(1,Form("\n New search gives %d local maxima for cathode 1 \n",fNLocal[0]));
1255         AliDebug(1,Form("\n                  %d local maxima for cathode 2 \n",fNLocal[1]));
1256         AliDebug(1,Form("\n New search gives %d %d \n",fNLocal[0],fNLocal[1]));
1257         if (fNLocal[cath]>2) {
1258             fNLocal[cath]=iback;
1259         }
1260     } // 2,1 local maxima
1261 }
1262
1263
1264 void  AliMUONClusterFinderVS::FillCluster(AliMUONRawCluster* c, Int_t flag, Int_t cath) 
1265 {
1266 //
1267 //  Completes cluster information starting from list of digits
1268 //
1269     AliMUONDigit* dig;
1270     Float_t x, y, z;
1271     Int_t  ix, iy;
1272     
1273     if (cath==1) {
1274         c->SetPeakSignal(cath,c->GetPeakSignal(0));     
1275     } else {
1276         c->SetPeakSignal(cath,0);
1277     }
1278     
1279     
1280     if (flag) {
1281         c->SetX(cath,0.);
1282         c->SetY(cath,0.);
1283         c->SetCharge(cath,0);
1284     }
1285
1286     AliDebug(1,Form("\n fPeakSignal %d\n",c->GetPeakSignal(cath)));
1287     for (Int_t i=0; i<c->GetMultiplicity(cath); i++)
1288     {
1289         dig= fInput->Digit(cath,c->GetIndex(i,cath));
1290         ix=dig->PadX()+c->GetOffset(i,cath);
1291         iy=dig->PadY();
1292         Int_t q=dig->Signal();
1293         if (!flag) q=Int_t(q*c->GetContrib(i,cath));
1294 //      fprintf(stderr,"q %d c->fPeakSignal[ %d ] %d\n",q,cath,c->fPeakSignal[cath]);
1295         if (dig->Physics() >= dig->Signal()) {
1296             c->SetPhysics(i,2);
1297         } else if (dig->Physics() == 0) {
1298             c->SetPhysics(i,0);
1299         } else  c->SetPhysics(i,1);
1300 //
1301 // 
1302         AliDebug(2,Form("q %d c->fPeakSignal[cath] %d\n",q,c->GetPeakSignal(cath)));
1303 // peak signal and track list
1304         if (q>c->GetPeakSignal(cath)) {
1305             c->SetPeakSignal(cath, q);
1306             c->SetTrack(0,dig->Hit());
1307             c->SetTrack(1,dig->Track(0));
1308             c->SetTrack(2,dig->Track(1));
1309 //          fprintf(stderr," c->fTracks[0] %d c->fTracks[1] %d\n",dig->fHit,dig->fTracks[0]);
1310         }
1311 //
1312         if (flag) {
1313           if (fSegmentationType == 1) 
1314             fSeg[cath]->GetPadC(ix, iy, x, y, z);
1315           else
1316             fSeg2[cath]->GetPadC(fInput->DetElemId(), ix, iy, x, y, z);
1317         
1318             c->AddX(cath, q*x);
1319             c->AddY(cath, q*y);
1320             c->AddCharge(cath, q);
1321         }
1322     } // loop over digits
1323     AliDebug(1," fin du cluster c\n");
1324
1325
1326     if (flag) {
1327         c->SetX(cath, c->GetX(cath)/c->GetCharge(cath));
1328 // Force on anod
1329         if (fSegmentationType == 1) 
1330           c->SetX(cath, fSeg[cath]->GetAnod(c->GetX(cath)));
1331         else
1332           c->SetX(cath, fSeg2[cath]->GetAnod(fInput->DetElemId(), c->GetX(cath)));
1333         c->SetY(cath, c->GetY(cath)/c->GetCharge(cath)); 
1334 //
1335 //  apply correction to the coordinate along the anode wire
1336 //
1337         x=c->GetX(cath);   
1338         y=c->GetY(cath);
1339         TF1* cogCorr;
1340         Int_t isec;
1341         if (fSegmentationType == 1) {
1342           fSeg[cath]->GetPadI(x, y, fZPlane, ix, iy);
1343           fSeg[cath]->GetPadC(ix, iy, x, y, z);
1344           isec=fSeg[cath]->Sector(ix,iy);
1345           cogCorr = fSeg[cath]->CorrFunc(isec-1);
1346         } else {
1347           fSeg2[cath]->GetPadI(fInput->DetElemId(), x, y, fZPlane, ix, iy);
1348           fSeg2[cath]->GetPadC(fInput->DetElemId(), ix, iy, x, y, z);
1349           isec=fSeg2[cath]->Sector(fInput->DetElemId(), ix,iy);
1350           cogCorr = fSeg2[cath]->CorrFunc(fInput->DetElemId(), isec-1);
1351         }
1352         
1353         if (cogCorr) {
1354           Float_t yOnPad;
1355           if (fSegmentationType == 1) 
1356             yOnPad=(c->GetY(cath)-y)/fSeg[cath]->Dpy(isec);
1357           else 
1358             yOnPad=(c->GetY(cath)-y)/fSeg2[cath]->Dpy(fInput->DetElemId(), isec);
1359
1360           c->SetY(cath, c->GetY(cath)-cogCorr->Eval(yOnPad, 0, 0));
1361         }
1362     }
1363 }
1364
1365 void  AliMUONClusterFinderVS::FillCluster(AliMUONRawCluster* c, Int_t cath) 
1366 {
1367 //
1368 //  Completes cluster information starting from list of digits
1369 //
1370     static Float_t dr0;
1371
1372     AliMUONDigit* dig;
1373
1374     if (cath==0) {
1375         dr0 = 10000;
1376     }
1377     
1378     Float_t xpad, ypad, zpad;
1379     Float_t dx, dy, dr;
1380
1381     for (Int_t i=0; i<c->GetMultiplicity(cath); i++)
1382     {
1383         dig = fInput->Digit(cath,c->GetIndex(i,cath));
1384         if (fSegmentationType == 1) 
1385           fSeg[cath]->
1386             GetPadC(dig->PadX(),dig->PadY(),xpad,ypad, zpad);
1387         else
1388           fSeg2[cath]->
1389             GetPadC(fInput->DetElemId(),dig->PadX(),dig->PadY(),xpad,ypad, zpad);
1390         AliDebug(1,Form("x %f y %f cx %f cy %f\n",xpad,ypad,c->GetX(0),c->GetY(0)));
1391         dx = xpad - c->GetX(0);
1392         dy = ypad - c->GetY(0);
1393         dr = TMath::Sqrt(dx*dx+dy*dy);
1394
1395         if (dr < dr0) {
1396             dr0 = dr;
1397             AliDebug(1,Form(" dr %f\n",dr));
1398             Int_t q=dig->Signal();
1399             if (dig->Physics() >= dig->Signal()) {
1400                 c->SetPhysics(i,2);
1401             } else if (dig->Physics() == 0) {
1402                 c->SetPhysics(i,0);
1403             } else  c->SetPhysics(i,1);
1404             c->SetPeakSignal(cath,q);
1405             c->SetTrack(0,dig->Hit());
1406             c->SetTrack(1,dig->Track(0));
1407             c->SetTrack(2,dig->Track(1));
1408             AliDebug(1,Form(" c->fTracks[0] %d c->fTracks[1] %d\n",dig->Hit(),
1409                     dig->Track(0)));
1410         }
1411 //
1412     } // loop over digits
1413
1414 //  apply correction to the coordinate along the anode wire
1415 // Force on anod
1416     if (fSegmentationType == 1) 
1417       c->SetX(cath,fSeg[cath]->GetAnod(c->GetX(cath)));
1418     else
1419       c->SetX(cath,fSeg2[cath]->GetAnod(fInput->DetElemId(), c->GetX(cath)));
1420 }
1421
1422 void  AliMUONClusterFinderVS::FindCluster(Int_t i, Int_t j, Int_t cath, AliMUONRawCluster &c){
1423
1424
1425 //
1426 //  Find a super cluster on both cathodes
1427 //
1428 //
1429 //  Add i,j as element of the cluster
1430 //
1431     
1432     Int_t idx = fHitMap[cath]->GetHitIndex(i,j);
1433     AliMUONDigit* dig = (AliMUONDigit*) fHitMap[cath]->GetHit(i,j);
1434     Int_t q=dig->Signal();
1435     Int_t theX=dig->PadX();
1436     Int_t theY=dig->PadY(); 
1437    
1438     if (q > TMath::Abs(c.GetPeakSignal(0)) && q > TMath::Abs(c.GetPeakSignal(1))) {
1439         c.SetPeakSignal(cath,q);
1440         c.SetTrack(0,dig->Hit());
1441         c.SetTrack(1,dig->Track(0));
1442         c.SetTrack(2,dig->Track(1));
1443     }
1444
1445 //
1446 //  Make sure that list of digits is ordered 
1447 // 
1448     Int_t mu=c.GetMultiplicity(cath);
1449     c.SetIndex(mu, cath, idx);
1450     
1451     if (dig->Physics() >= dig->Signal()) {
1452         c.SetPhysics(mu,2);
1453     } else if (dig->Physics() == 0) {
1454         c.SetPhysics(mu,0);
1455     } else  c.SetPhysics(mu,1);
1456
1457     
1458     if (mu > 0) {
1459         for (Int_t ind = mu-1; ind >= 0; ind--) {
1460             Int_t ist=c.GetIndex(ind,cath);
1461             Int_t ql=fInput->Digit(cath, ist)->Signal();
1462             Int_t ix=fInput->Digit(cath, ist)->PadX();
1463             Int_t iy=fInput->Digit(cath, ist)->PadY();
1464             
1465             if (q>ql || (q==ql && theX > ix && theY < iy)) {
1466                 c.SetIndex(ind, cath, idx);
1467                 c.SetIndex(ind+1, cath, ist);
1468             } else {
1469                 
1470                 break;
1471             }
1472         }
1473     }
1474
1475     c.SetMultiplicity(cath, c.GetMultiplicity(cath)+1);
1476     if (c.GetMultiplicity(cath) >= 50 ) {
1477       AliDebug(1,Form("FindCluster - multiplicity >50  %d \n",c.GetMultiplicity(0)));
1478         c.SetMultiplicity(cath, 49);
1479     }
1480
1481 // Prepare center of gravity calculation
1482     Float_t x, y, z;
1483     if (fSegmentationType == 1) 
1484       fSeg[cath]->GetPadC(i, j, x, y, z);
1485     else
1486       fSeg2[cath]->GetPadC(fInput->DetElemId(), i, j, x, y, z);
1487     c.AddX(cath,q*x);
1488     c.AddY(cath,q*y);
1489     c.AddCharge(cath,q);
1490 //
1491 // Flag hit as "taken"  
1492     fHitMap[cath]->FlagHit(i,j);
1493 //
1494 //  Now look recursively for all neighbours and pad hit on opposite cathode
1495 //
1496 //  Loop over neighbours
1497     Int_t ix,iy;
1498     ix=iy=0;
1499     Int_t nn;
1500     Int_t xList[10], yList[10];
1501     if (fSegmentationType == 1) 
1502       fSeg[cath]->Neighbours(i,j,&nn,xList,yList);
1503     else 
1504       fSeg2[cath]->Neighbours(fInput->DetElemId(), i,j,&nn,xList,yList);
1505     for (Int_t in=0; in<nn; in++) {
1506         ix=xList[in];
1507         iy=yList[in];
1508         
1509         if (fHitMap[cath]->TestHit(ix,iy)==kUnused) {
1510             AliDebug(2,Form("\n Neighbours %d %d %d", cath, ix, iy));
1511             FindCluster(ix, iy, cath, c);
1512         }
1513         
1514    }
1515     Int_t nOpp=0;
1516     Int_t iXopp[50], iYopp[50];
1517     
1518 //  Neighbours on opposite cathode 
1519 //  Take into account that several pads can overlap with the present pad
1520     Int_t isec;
1521     if (fSegmentationType == 1) 
1522       isec=fSeg[cath]->Sector(i,j);    
1523     else
1524       isec=fSeg2[cath]->Sector(fInput->DetElemId(), i,j);    
1525
1526     Int_t iop;
1527     Float_t dx, dy;
1528
1529     if (fSegmentationType == 1) {
1530       if (cath==0) {
1531         iop = 1;
1532         dx  = (fSeg[cath]->Dpx(isec))/2.;
1533         dy  = 0.;
1534       } else {
1535         iop = 0;
1536         dx  = 0.;
1537         dy  = (fSeg[cath]->Dpy(isec))/2;
1538       }
1539    
1540
1541     
1542       // loop over pad neighbours on opposite cathode
1543       for (fSeg[iop]->FirstPad(x, y, fZPlane, dx, dy);
1544            fSeg[iop]->MorePads();
1545            fSeg[iop]->NextPad())
1546         {
1547         
1548           ix = fSeg[iop]->Ix(); iy = fSeg[iop]->Iy();
1549           AliDebug(2,Form("\n ix, iy: %f %f %f %d %d %d", x,y,z,ix, iy, fSector));
1550           if (fHitMap[iop]->TestHit(ix,iy)==kUnused){
1551             iXopp[nOpp]=ix;
1552             iYopp[nOpp++]=iy;
1553             AliDebug(2,Form("\n Opposite %d %d %d", iop, ix, iy));
1554           }
1555         
1556         } // Loop over pad neighbours
1557       //  This had to go outside the loop since recursive calls inside the iterator are not possible
1558       //
1559       Int_t jopp;
1560       for (jopp=0; jopp<nOpp; jopp++) {
1561         if (fHitMap[iop]->TestHit(iXopp[jopp],iYopp[jopp]) == kUnused) 
1562           FindCluster(iXopp[jopp], iYopp[jopp], iop, c);
1563       }
1564     } else {
1565
1566       if (cath==0) {
1567         iop = 1;
1568         dx  = (fSeg2[cath]->Dpx(fInput->DetElemId(), isec))/2.;
1569         dy  = 0.;
1570       } else {
1571         iop = 0;
1572         dx  = 0.;
1573         dy  = (fSeg2[cath]->Dpy(fInput->DetElemId(), isec))/2;
1574       }
1575    
1576
1577     
1578       // loop over pad neighbours on opposite cathode
1579       for (fSeg2[iop]->FirstPad(fInput->DetElemId(), x, y, fZPlane, dx, dy);
1580            fSeg2[iop]->MorePads(fInput->DetElemId());
1581            fSeg2[iop]->NextPad(fInput->DetElemId()))
1582         {
1583         
1584           ix = fSeg2[iop]->Ix(); iy = fSeg2[iop]->Iy();
1585           AliDebug(2,Form("\n ix, iy: %f %f %f %d %d %d", x,y,z,ix, iy, fSector));
1586           if (fHitMap[iop]->TestHit(ix,iy)==kUnused){
1587             iXopp[nOpp]=ix;
1588             iYopp[nOpp++]=iy;
1589             AliDebug(2,Form("\n Opposite %d %d %d", iop, ix, iy));
1590           }
1591         
1592         } // Loop over pad neighbours
1593       //  This had to go outside the loop since recursive calls inside the iterator are not possible
1594       //
1595       Int_t jopp;
1596       for (jopp=0; jopp<nOpp; jopp++) {
1597         if (fHitMap[iop]->TestHit(iXopp[jopp],iYopp[jopp]) == kUnused) 
1598           FindCluster(iXopp[jopp], iYopp[jopp], iop, c);
1599       }
1600     }
1601 }
1602
1603 //_____________________________________________________________________________
1604
1605 void AliMUONClusterFinderVS::FindRawClusters()
1606 {
1607   //
1608   // MUON cluster finder from digits -- finds neighbours on both cathodes and 
1609   // fills the tree with raw clusters
1610   //
1611
1612     ResetRawClusters();
1613 //  Return if no input datad available
1614     if (!fInput->NDigits(0) && !fInput->NDigits(1)) return;
1615
1616     fSegmentationType = fInput->GetSegmentationType();
1617
1618     if (fSegmentationType == 1) {
1619       fSeg[0] = fInput->Segmentation(0);
1620       fSeg[1] = fInput->Segmentation(1);
1621
1622       fHitMap[0]  = new AliMUONHitMapA1(fSeg[0], fInput->Digits(0));
1623       fHitMap[1]  = new AliMUONHitMapA1(fSeg[1], fInput->Digits(1));
1624
1625     } else {
1626       fSeg2[0] = fInput->Segmentation2(0);
1627       fSeg2[1] = fInput->Segmentation2(1);
1628
1629       fHitMap[0]  = new AliMUONHitMapA1(fInput->DetElemId(), fSeg2[0], fInput->Digits(0));
1630       fHitMap[1]  = new AliMUONHitMapA1(fInput->DetElemId(), fSeg2[1], fInput->Digits(1));
1631     }
1632  
1633     AliMUONDigit *dig;
1634
1635     Int_t ndig, cath;
1636     Int_t nskip=0;
1637     Int_t ncls=0;
1638     fHitMap[0]->FillHits();
1639     fHitMap[1]->FillHits();
1640 //
1641 //  Outer Loop over Cathodes
1642     for (cath=0; cath<2; cath++) {
1643         for (ndig=0; ndig<fInput->NDigits(cath); ndig++) {
1644             dig = fInput->Digit(cath, ndig);
1645             Int_t i=dig->PadX();
1646             Int_t j=dig->PadY();
1647             if (fHitMap[cath]->TestHit(i,j)==kUsed ||fHitMap[0]->TestHit(i,j)==kEmpty) {
1648                 nskip++;
1649                 continue;
1650             }
1651             AliDebug(1,Form("\n CATHODE %d CLUSTER %d\n",cath,ncls));
1652             AliMUONRawCluster c;
1653             c.SetMultiplicity(0, 0);
1654             c.SetMultiplicity(1, 0);
1655             c.SetPeakSignal(cath,dig->Signal());
1656             c.SetTrack(0, dig->Hit());
1657             c.SetTrack(1, dig->Track(0));
1658             c.SetTrack(2, dig->Track(1));
1659             // tag the beginning of cluster list in a raw cluster
1660             c.SetNcluster(0,-1);
1661             Float_t xcu, ycu;
1662             if (fSegmentationType == 1) {
1663               fSeg[cath]->GetPadC(i,j,xcu, ycu, fZPlane);
1664               fSector= fSeg[cath]->Sector(i,j)/100;
1665             } else {
1666               fSeg2[cath]->GetPadC(fInput->DetElemId(), i,j,xcu, ycu, fZPlane);
1667               fSector= fSeg2[cath]->Sector(fInput->DetElemId(), i,j)/100;
1668             }
1669             AliDebug(1,Form("\n New Seed %d %d ", i,j));
1670         
1671             
1672             FindCluster(i,j,cath,c);
1673 //          ^^^^^^^^^^^^^^^^^^^^^^^^
1674             // center of gravity
1675             if (c.GetX(0)!=0.) c.SetX(0, c.GetX(0)/c.GetCharge(0)); // c.fX[0] /= c.fQ[0];
1676 // Force on anod
1677             if (fSegmentationType == 1) 
1678               c.SetX(0,fSeg[0]->GetAnod(c.GetX(0)));
1679             else 
1680               c.SetX(0,fSeg2[0]->GetAnod(fInput->DetElemId(), c.GetX(0)));
1681             if (c.GetY(0)!=0.) c.SetY(0, c.GetY(0)/c.GetCharge(0)); // c.fY[0] /= c.fQ[0];
1682             
1683             if(c.GetCharge(1)!=0.) c.SetX(1, c.GetX(1)/c.GetCharge(1));  // c.fX[1] /= c.fQ[1];
1684                                         
1685            // Force on anod
1686             if (fSegmentationType == 1) 
1687               c.SetX(1, fSeg[0]->GetAnod(c.GetX(1)));
1688             else 
1689               c.SetX(1, fSeg2[0]->GetAnod(fInput->DetElemId(),c.GetX(1)));
1690             if(c.GetCharge(1)!=0.) c.SetY(1, c.GetY(1)/c.GetCharge(1));// c.fY[1] /= c.fQ[1];
1691             
1692             c.SetZ(0, fZPlane);
1693             c.SetZ(1, fZPlane);     
1694
1695                 AliDebug(1,Form("\n Cathode 1 multiplicite %d X(CG) %f Y(CG) %f\n",
1696                         c.GetMultiplicity(0),c.GetX(0),c.GetY(0)));
1697                 AliDebug(1,Form(" Cathode 2 multiplicite %d X(CG) %f Y(CG) %f\n",
1698                         c.GetMultiplicity(1),c.GetX(1),c.GetY(1)));
1699 //      Analyse cluster and decluster if necessary
1700 //      
1701         ncls++;
1702         c.SetNcluster(1,fNRawClusters);
1703         c.SetClusterType(c.PhysicsContribution());
1704
1705         fNPeaks=0;
1706 //
1707 //
1708         Decluster(&c);
1709 //
1710 //      reset Cluster object
1711         { // begin local scope
1712             for (int k=0;k<c.GetMultiplicity(0);k++) c.SetIndex(k, 0, 0);
1713         } // end local scope
1714
1715         { // begin local scope
1716             for (int k=0;k<c.GetMultiplicity(1);k++) c.SetIndex(k, 1, 0);
1717         } // end local scope
1718         
1719         c.SetMultiplicity(0,0);
1720         c.SetMultiplicity(1,0);
1721
1722         
1723         } // end loop ndig
1724     } // end loop cathodes
1725     delete fHitMap[0];
1726     delete fHitMap[1];
1727 }
1728
1729 Float_t AliMUONClusterFinderVS::SingleMathiesonFit(AliMUONRawCluster *c, Int_t cath)
1730 {
1731 // Performs a single Mathieson fit on one cathode
1732 // 
1733     Double_t arglist[20];
1734     Int_t ierflag=0;
1735     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1736     
1737     clusterInput.Fitter()->SetFCN(fcnS1);
1738     clusterInput.Fitter()->mninit(2,10,7);
1739     clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
1740     arglist[0]=-1;
1741     clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1742 // Set starting values 
1743     static Double_t vstart[2];
1744     vstart[0]=c->GetX(1);
1745     vstart[1]=c->GetY(0);
1746     
1747     
1748 // lower and upper limits
1749     static Double_t lower[2], upper[2];
1750     Int_t ix,iy, isec;
1751     if (fSegmentationType == 1) {
1752       fSeg[cath]->GetPadI(c->GetX(cath), c->GetY(cath), fZPlane, ix, iy);
1753       isec=fSeg[cath]->Sector(ix, iy);
1754
1755       lower[0]=vstart[0]-fSeg[cath]->Dpx(isec)/2;
1756       lower[1]=vstart[1]-fSeg[cath]->Dpy(isec)/2;
1757     
1758       upper[0]=lower[0]+fSeg[cath]->Dpx(isec);
1759       upper[1]=lower[1]+fSeg[cath]->Dpy(isec);
1760
1761     } else {
1762       fSeg2[cath]->GetPadI(fInput->DetElemId(), c->GetX(cath), c->GetY(cath), fZPlane, ix, iy);
1763       isec=fSeg2[cath]->Sector(fInput->DetElemId(), ix, iy);
1764
1765       lower[0]=vstart[0]-fSeg2[cath]->Dpx(fInput->DetElemId(), isec)/2;
1766       lower[1]=vstart[1]-fSeg2[cath]->Dpy(fInput->DetElemId(), isec)/2;
1767     
1768       upper[0]=lower[0]+fSeg2[cath]->Dpx(fInput->DetElemId(), isec);
1769       upper[1]=lower[1]+fSeg2[cath]->Dpy(fInput->DetElemId(), isec);
1770     }
1771
1772 // step sizes
1773     static Double_t step[2]={0.0005, 0.0005};
1774     
1775     clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1776     clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1777 // ready for minimisation       
1778     arglist[0]= -1;
1779     arglist[1]= 0;
1780     
1781     clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1782     clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1783     //    clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1784     Double_t fmin, fedm, errdef;
1785     Int_t   npari, nparx, istat;
1786       
1787     clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);  
1788     fFitStat=istat;
1789     
1790 // Print results
1791 // Get fitted parameters
1792     Double_t xrec, yrec;
1793     TString chname;
1794     Double_t epxz, b1, b2;
1795     Int_t ierflg;
1796     clusterInput.Fitter()->mnpout(0, chname, xrec, epxz, b1, b2, ierflg);       
1797     clusterInput.Fitter()->mnpout(1, chname, yrec, epxz, b1, b2, ierflg);       
1798     fXFit[cath]=xrec;
1799     fYFit[cath]=yrec;
1800     return fmin;
1801 }
1802
1803 Float_t AliMUONClusterFinderVS::CombiSingleMathiesonFit(AliMUONRawCluster * /*c*/)
1804 {
1805 // Perform combined Mathieson fit on both cathode planes
1806 //
1807     Double_t arglist[20];
1808     Int_t ierflag=0;
1809     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1810     clusterInput.Fitter()->SetFCN(fcnCombiS1);
1811     clusterInput.Fitter()->mninit(2,10,7);
1812     clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
1813     arglist[0]=-1;
1814     clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1815     static Double_t vstart[2];
1816     vstart[0]=fXInit[0];
1817     vstart[1]=fYInit[0];
1818     
1819     
1820 // lower and upper limits
1821     static Float_t lower[2], upper[2];
1822     Int_t ix,iy,isec;
1823     Float_t dpy, dpx;
1824
1825     if (fSegmentationType == 1) {
1826       fSeg[0]->GetPadI(fXInit[0], fYInit[0], fZPlane, ix, iy);
1827       isec=fSeg[0]->Sector(ix, iy);
1828       dpy=fSeg[0]->Dpy(isec);
1829       fSeg[1]->GetPadI(fXInit[0], fYInit[0], fZPlane, ix, iy);
1830       isec=fSeg[1]->Sector(ix, iy);
1831       dpx=fSeg[1]->Dpx(isec);
1832
1833     } else {
1834       fSeg2[0]->GetPadI(fInput->DetElemId(), fXInit[0], fYInit[0], fZPlane, ix, iy);
1835       isec=fSeg2[0]->Sector(fInput->DetElemId(), ix, iy);
1836       dpy=fSeg2[0]->Dpy(fInput->DetElemId(), isec);
1837       fSeg2[1]->GetPadI(fInput->DetElemId(), fXInit[0], fYInit[0], fZPlane, ix, iy);
1838       isec=fSeg2[1]->Sector(fInput->DetElemId(), ix, iy);
1839       dpx=fSeg2[1]->Dpx(fInput->DetElemId(), isec);
1840
1841     }
1842     Int_t icount;
1843     Float_t xdum, ydum, zdum;
1844
1845 //  Find save upper and lower limits    
1846     
1847     icount = 0;
1848     if (fSegmentationType == 1) {
1849       for (fSeg[1]->FirstPad(fXInit[0], fYInit[0], fZPlane, dpx, 0.); 
1850            fSeg[1]->MorePads(); 
1851            fSeg[1]->NextPad())
1852         {
1853           ix=fSeg[1]->Ix(); iy=fSeg[1]->Iy();
1854           fSeg[1]->GetPadC(ix,iy, upper[0], ydum, zdum);        
1855           if (icount ==0) lower[0]=upper[0];
1856           icount++;
1857         }
1858     } else {
1859       for (fSeg2[1]->FirstPad(fInput->DetElemId(),fXInit[0], fYInit[0], fZPlane, dpx, 0.); 
1860            fSeg2[1]->MorePads(fInput->DetElemId()); 
1861            fSeg2[1]->NextPad(fInput->DetElemId()))
1862         {
1863           ix=fSeg2[1]->Ix(); iy=fSeg2[1]->Iy();
1864           fSeg2[1]->GetPadC(fInput->DetElemId(), ix,iy, upper[0], ydum, zdum);  
1865           if (icount ==0) lower[0]=upper[0];
1866           icount++;
1867         }
1868     }
1869     if (lower[0]>upper[0]) {xdum=lower[0]; lower[0]=upper[0]; upper[0]=xdum;}
1870         
1871     icount=0;
1872     AliDebug(1,Form("\n single y %f %f", fXInit[0], fYInit[0]));
1873     
1874     if (fSegmentationType == 1) {
1875       for (fSeg[0]->FirstPad(fXInit[0], fYInit[0], fZPlane, 0., dpy); 
1876            fSeg[0]->MorePads(); 
1877            fSeg[0]->NextPad())
1878         {
1879           ix=fSeg[0]->Ix(); iy=fSeg[0]->Iy();
1880           fSeg[0]->GetPadC(ix,iy,xdum,upper[1],zdum);   
1881           if (icount ==0) lower[1]=upper[1];
1882           icount++;
1883           AliDebug(1,Form("\n upper lower %d %f %f", icount, upper[1], lower[1]));
1884         }
1885     } else {
1886       for (fSeg2[0]->FirstPad(fInput->DetElemId(), fXInit[0], fYInit[0], fZPlane, 0., dpy); 
1887            fSeg2[0]->MorePads(fInput->DetElemId()); 
1888            fSeg2[0]->NextPad(fInput->DetElemId()))
1889         {
1890           ix=fSeg2[0]->Ix(); iy=fSeg2[0]->Iy();
1891           fSeg2[0]->GetPadC(fInput->DetElemId(), ix,iy,xdum,upper[1],zdum);     
1892           if (icount ==0) lower[1]=upper[1];
1893           icount++;
1894           AliDebug(1,Form("\n upper lower %d %f %f", icount, upper[1], lower[1]));
1895         }
1896     }
1897     if (lower[1]>upper[1]) {xdum=lower[1]; lower[1]=upper[1]; upper[1]=xdum;}
1898
1899 // step sizes
1900     static Double_t step[2]={0.00001, 0.0001};
1901     
1902     clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1903     clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1904 // ready for minimisation       
1905     arglist[0]= -1;
1906     arglist[1]= 0;
1907     
1908     clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1909     clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1910     //    clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1911     Double_t fmin, fedm, errdef;
1912     Int_t   npari, nparx, istat;
1913       
1914     clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);  
1915     fFitStat=istat;
1916     
1917 // Print results
1918 // Get fitted parameters
1919     Double_t xrec, yrec;
1920     TString chname;
1921     Double_t epxz, b1, b2;
1922     Int_t ierflg;
1923     clusterInput.Fitter()->mnpout(0, chname, xrec, epxz, b1, b2, ierflg);       
1924     clusterInput.Fitter()->mnpout(1, chname, yrec, epxz, b1, b2, ierflg);       
1925     fXFit[0]=xrec;
1926     fYFit[0]=yrec;
1927     return fmin;
1928 }
1929
1930 Bool_t AliMUONClusterFinderVS::DoubleMathiesonFit(AliMUONRawCluster * /*c*/, Int_t cath)
1931 {
1932 // Performs a double Mathieson fit on one cathode
1933 // 
1934
1935 //
1936 //  Initialise global variables for fit
1937     Double_t arglist[20];
1938     Int_t ierflag=0;
1939     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1940     clusterInput.Fitter()->SetFCN(fcnS2);
1941     clusterInput.Fitter()->mninit(5,10,7);
1942     clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
1943     arglist[0]=-1;
1944     clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1945 // Set starting values 
1946     static Double_t vstart[5];
1947     vstart[0]=fX[fIndLocal[0][cath]][cath];
1948     vstart[1]=fY[fIndLocal[0][cath]][cath];     
1949     vstart[2]=fX[fIndLocal[1][cath]][cath];
1950     vstart[3]=fY[fIndLocal[1][cath]][cath];     
1951     vstart[4]=Float_t(fQ[fIndLocal[0][cath]][cath])/
1952         Float_t(fQ[fIndLocal[0][cath]][cath]+fQ[fIndLocal[1][cath]][cath]);
1953 // lower and upper limits
1954     static Float_t lower[5], upper[5];
1955     Int_t isec;
1956
1957     if (fSegmentationType == 1) {
1958       isec=fSeg[cath]->Sector(fIx[fIndLocal[0][cath]][cath], fIy[fIndLocal[0][cath]][cath]);
1959       lower[0]=vstart[0]-fSeg[cath]->Dpx(isec);
1960       lower[1]=vstart[1]-fSeg[cath]->Dpy(isec);
1961     
1962       upper[0]=lower[0]+2.*fSeg[cath]->Dpx(isec);
1963       upper[1]=lower[1]+2.*fSeg[cath]->Dpy(isec);
1964     
1965       isec=fSeg[cath]->Sector(fIx[fIndLocal[1][cath]][cath], fIy[fIndLocal[1][cath]][cath]);
1966       lower[2]=vstart[2]-fSeg[cath]->Dpx(isec)/2;
1967       lower[3]=vstart[3]-fSeg[cath]->Dpy(isec)/2;
1968     
1969       upper[2]=lower[2]+fSeg[cath]->Dpx(isec);
1970       upper[3]=lower[3]+fSeg[cath]->Dpy(isec);
1971
1972     } else {
1973       isec=fSeg2[cath]->Sector(fInput->DetElemId(),fIx[fIndLocal[0][cath]][cath], 
1974                                fIy[fIndLocal[0][cath]][cath]);
1975       lower[0]=vstart[0]-fSeg2[cath]->Dpx(fInput->DetElemId(),isec);
1976       lower[1]=vstart[1]-fSeg2[cath]->Dpy(fInput->DetElemId(),isec);
1977     
1978       upper[0]=lower[0]+2.*fSeg2[cath]->Dpx(fInput->DetElemId(),isec);
1979       upper[1]=lower[1]+2.*fSeg2[cath]->Dpy(fInput->DetElemId(),isec);
1980     
1981       isec=fSeg2[cath]->Sector(fInput->DetElemId(),fIx[fIndLocal[1][cath]][cath], 
1982                               fIy[fIndLocal[1][cath]][cath]);
1983       lower[2]=vstart[2]-fSeg2[cath]->Dpx(fInput->DetElemId(),isec)/2;
1984       lower[3]=vstart[3]-fSeg2[cath]->Dpy(fInput->DetElemId(),isec)/2;
1985     
1986       upper[2]=lower[2]+fSeg2[cath]->Dpx(fInput->DetElemId(),isec);
1987       upper[1]=lower[1]+2.*fSeg2[cath]->Dpy(fInput->DetElemId(),isec);
1988
1989     }
1990
1991     lower[4]=0.;
1992     upper[4]=1.;
1993 // step sizes
1994     static Double_t step[5]={0.0005, 0.0005, 0.0005, 0.0005, 0.0001};
1995     
1996     clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1997     clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1998     clusterInput.Fitter()->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
1999     clusterInput.Fitter()->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
2000     clusterInput.Fitter()->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
2001 // ready for minimisation       
2002     arglist[0]= -1;
2003     arglist[1]= 0;
2004     
2005     clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
2006     clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
2007     //    clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
2008 // Get fitted parameters
2009     Double_t xrec[2], yrec[2], qfrac;
2010     TString chname;
2011     Double_t epxz, b1, b2;
2012     Int_t ierflg;
2013     clusterInput.Fitter()->mnpout(0, chname, xrec[0], epxz, b1, b2, ierflg);    
2014     clusterInput.Fitter()->mnpout(1, chname, yrec[0], epxz, b1, b2, ierflg);    
2015     clusterInput.Fitter()->mnpout(2, chname, xrec[1], epxz, b1, b2, ierflg);    
2016     clusterInput.Fitter()->mnpout(3, chname, yrec[1], epxz, b1, b2, ierflg);    
2017     clusterInput.Fitter()->mnpout(4, chname, qfrac,   epxz, b1, b2, ierflg);    
2018
2019     Double_t fmin, fedm, errdef;
2020     Int_t   npari, nparx, istat;
2021       
2022     clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);  
2023     fFitStat=istat;
2024     return kTRUE;
2025 }
2026
2027 Float_t AliMUONClusterFinderVS::CombiDoubleMathiesonFit(AliMUONRawCluster * /*c*/)
2028 {
2029 //
2030 // Perform combined double Mathieson fit on both cathode planes
2031 //
2032     Double_t arglist[20];
2033     Int_t ierflag=0;
2034     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2035     clusterInput.Fitter()->SetFCN(fcnCombiS2);
2036     clusterInput.Fitter()->mninit(6,10,7);
2037     clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
2038     arglist[0]=-1;
2039     clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
2040 // Set starting values 
2041     static Double_t vstart[6];
2042     vstart[0]=fXInit[0];
2043     vstart[1]=fYInit[0];
2044     vstart[2]=fXInit[1];
2045     vstart[3]=fYInit[1];
2046     vstart[4]=fQrInit[0];
2047     vstart[5]=fQrInit[1];
2048 // lower and upper limits
2049     static Float_t lower[6], upper[6];
2050     Int_t ix,iy,isec;
2051     Float_t dpx, dpy;
2052   if (fSegmentationType == 1) {
2053     fSeg[1]->GetPadI(fXInit[0], fYInit[0], fZPlane, ix, iy);
2054     isec=fSeg[1]->Sector(ix, iy);
2055     dpx=fSeg[1]->Dpx(isec);
2056
2057     fSeg[0]->GetPadI(fXInit[0], fYInit[0], fZPlane, ix, iy);
2058     isec=fSeg[0]->Sector(ix, iy);
2059     dpy=fSeg[0]->Dpy(isec);
2060
2061   } else {
2062     fSeg2[1]->GetPadI(fInput->DetElemId(),fXInit[0], fYInit[0], fZPlane, ix, iy);
2063     isec=fSeg2[1]->Sector(fInput->DetElemId(),ix, iy);
2064     dpx=fSeg2[1]->Dpx(fInput->DetElemId(), isec);
2065
2066     fSeg2[0]->GetPadI(fInput->DetElemId(), fXInit[0], fYInit[0], fZPlane, ix, iy);
2067     isec=fSeg2[0]->Sector(fInput->DetElemId(), ix, iy);
2068     dpy=fSeg2[0]->Dpy(fInput->DetElemId(), isec);
2069
2070   }
2071
2072     Int_t icount;
2073     Float_t xdum, ydum, zdum;
2074     AliDebug(1,Form("\n Cluster Finder: %f %f %f %f  ", fXInit[0], fXInit[1],fYInit[0], fYInit[1] ));
2075
2076     if (fSegmentationType == 1) {
2077  
2078       //  Find save upper and lower limits    
2079       icount = 0;
2080     
2081       for (fSeg[1]->FirstPad(fXInit[0], fYInit[0], fZPlane, dpx, 0.); 
2082            fSeg[1]->MorePads(); 
2083            fSeg[1]->NextPad())
2084         {
2085           ix=fSeg[1]->Ix(); iy=fSeg[1]->Iy();
2086           //    if (fHitMap[1]->TestHit(ix, iy) == kEmpty) continue;
2087           fSeg[1]->GetPadC(ix,iy,upper[0],ydum,zdum);   
2088           if (icount ==0) lower[0]=upper[0];
2089           icount++;
2090         }
2091       if (lower[0]>upper[0]) {xdum=lower[0]; lower[0]=upper[0]; upper[0]=xdum;}    
2092       //    vstart[0] = 0.5*(lower[0]+upper[0]);
2093
2094     
2095       icount=0;
2096     
2097       for (fSeg[0]->FirstPad(fXInit[0], fYInit[0], fZPlane, 0., dpy); 
2098            fSeg[0]->MorePads(); 
2099            fSeg[0]->NextPad())
2100         {
2101           ix=fSeg[0]->Ix(); iy=fSeg[0]->Iy();
2102           //    if (fHitMap[0]->TestHit(ix, iy) == kEmpty) continue;
2103           fSeg[0]->GetPadC(ix,iy,xdum,upper[1],zdum);   
2104           if (icount ==0) lower[1]=upper[1];
2105           icount++;
2106         }
2107     
2108       if (lower[1]>upper[1]) {xdum=lower[1]; lower[1]=upper[1]; upper[1]=xdum;}    
2109       //     vstart[1] = 0.5*(lower[1]+upper[1]);
2110
2111
2112       fSeg[1]->GetPadI(fXInit[1], fYInit[1], fZPlane, ix, iy);
2113       isec=fSeg[1]->Sector(ix, iy);
2114       dpx=fSeg[1]->Dpx(isec);
2115       fSeg[0]->GetPadI(fXInit[1], fYInit[1], fZPlane, ix, iy);
2116       isec=fSeg[0]->Sector(ix, iy);
2117       dpy=fSeg[0]->Dpy(isec);
2118
2119
2120       //  Find save upper and lower limits    
2121
2122       icount=0;
2123     
2124       for (fSeg[1]->FirstPad(fXInit[1], fYInit[1], fZPlane, dpx, 0); 
2125            fSeg[1]->MorePads(); fSeg[1]->NextPad())
2126         {
2127           ix=fSeg[1]->Ix(); iy=fSeg[1]->Iy();
2128           //    if (fHitMap[1]->TestHit(ix, iy) == kEmpty) continue;
2129           fSeg[1]->GetPadC(ix,iy,upper[2],ydum,zdum);   
2130           if (icount ==0) lower[2]=upper[2];
2131           icount++;
2132         }
2133       if (lower[2]>upper[2]) {xdum=lower[2]; lower[2]=upper[2]; upper[2]=xdum;}    
2134       //    vstart[2] = 0.5*(lower[2]+upper[2]);
2135
2136       icount=0;
2137     
2138       for (fSeg[0]->FirstPad(fXInit[1], fYInit[1], fZPlane, 0, dpy); 
2139            fSeg[0]-> MorePads(); fSeg[0]->NextPad())
2140         {
2141           ix=fSeg[0]->Ix(); iy=fSeg[0]->Iy();
2142           //    if (fHitMap[0]->TestHit(ix, iy) != kEmpty) continue;
2143         
2144           fSeg[0]->GetPadC(ix,iy,xdum,upper[3],zdum);   
2145           if (icount ==0) lower[3]=upper[3];
2146           icount++;
2147
2148         }
2149       if (lower[3]>upper[3]) {xdum=lower[3]; lower[3]=upper[3]; upper[3]=xdum;}    
2150     
2151       //     vstart[3] = 0.5*(lower[3]+upper[3]);
2152     } else {
2153
2154     //  Find save upper and lower limits    
2155       icount = 0;
2156     
2157       for (fSeg2[1]->FirstPad(fInput->DetElemId(),fXInit[0], fYInit[0], fZPlane, dpx, 0.); 
2158            fSeg2[1]->MorePads(fInput->DetElemId()); 
2159            fSeg2[1]->NextPad(fInput->DetElemId()))
2160         {
2161           ix=fSeg2[1]->Ix(); iy=fSeg2[1]->Iy();
2162           //    if (fHitMap[1]->TestHit(ix, iy) == kEmpty) continue;
2163           fSeg2[1]->GetPadC(fInput->DetElemId(),ix,iy,upper[0],ydum,zdum);      
2164           if (icount ==0) lower[0]=upper[0];
2165           icount++;
2166         }
2167       if (lower[0]>upper[0]) {xdum=lower[0]; lower[0]=upper[0]; upper[0]=xdum;}    
2168       //    vstart[0] = 0.5*(lower[0]+upper[0]);
2169
2170     
2171       icount=0;
2172     
2173       for (fSeg2[0]->FirstPad(fInput->DetElemId(),fXInit[0], fYInit[0], fZPlane, 0., dpy); 
2174            fSeg2[0]->MorePads(fInput->DetElemId()); 
2175            fSeg2[0]->NextPad(fInput->DetElemId()))
2176         {
2177           ix=fSeg2[0]->Ix(); iy=fSeg2[0]->Iy();
2178           //    if (fHitMap[0]->TestHit(ix, iy) == kEmpty) continue;
2179           fSeg2[0]->GetPadC(fInput->DetElemId(),ix,iy,xdum,upper[1],zdum);      
2180           if (icount ==0) lower[1]=upper[1];
2181           icount++;
2182         }
2183     
2184       if (lower[1]>upper[1]) {xdum=lower[1]; lower[1]=upper[1]; upper[1]=xdum;}    
2185       //     vstart[1] = 0.5*(lower[1]+upper[1]);
2186
2187
2188       fSeg2[1]->GetPadI(fInput->DetElemId(),fXInit[1], fYInit[1], fZPlane, ix, iy);
2189       isec=fSeg2[1]->Sector(fInput->DetElemId(),ix, iy);
2190       dpx=fSeg2[1]->Dpx(fInput->DetElemId(),isec);
2191       fSeg2[0]->GetPadI(fInput->DetElemId(),fXInit[1], fYInit[1], fZPlane, ix, iy);
2192       isec=fSeg2[0]->Sector(fInput->DetElemId(),ix, iy);
2193       dpy=fSeg2[0]->Dpy(fInput->DetElemId(),isec);
2194
2195
2196       //  Find save upper and lower limits    
2197
2198       icount=0;
2199     
2200       for (fSeg2[1]->FirstPad(fInput->DetElemId(),fXInit[1], fYInit[1], fZPlane, dpx, 0); 
2201            fSeg2[1]->MorePads(fInput->DetElemId()); 
2202            fSeg2[1]->NextPad(fInput->DetElemId()))
2203         {
2204           ix=fSeg2[1]->Ix(); iy=fSeg2[1]->Iy();
2205           //    if (fHitMap[1]->TestHit(ix, iy) == kEmpty) continue;
2206           fSeg2[1]->GetPadC(fInput->DetElemId(),ix,iy,upper[2],ydum,zdum);      
2207           if (icount ==0) lower[2]=upper[2];
2208           icount++;
2209         }
2210       if (lower[2]>upper[2]) {xdum=lower[2]; lower[2]=upper[2]; upper[2]=xdum;}    
2211       //    vstart[2] = 0.5*(lower[2]+upper[2]);
2212
2213       icount=0;
2214     
2215       for (fSeg2[0]->FirstPad(fInput->DetElemId(),fXInit[1], fYInit[1], fZPlane, 0, dpy); 
2216            fSeg2[0]-> MorePads(fInput->DetElemId()); 
2217            fSeg2[0]->NextPad(fInput->DetElemId()))
2218         {
2219           ix=fSeg2[0]->Ix(); iy=fSeg2[0]->Iy();
2220           //    if (fHitMap[0]->TestHit(ix, iy) != kEmpty) continue;
2221         
2222           fSeg2[0]->GetPadC(fInput->DetElemId(),ix,iy,xdum,upper[3],zdum);      
2223           if (icount ==0) lower[3]=upper[3];
2224           icount++;
2225
2226         }
2227       if (lower[3]>upper[3]) {xdum=lower[3]; lower[3]=upper[3]; upper[3]=xdum;}  
2228     }
2229     lower[4]=0.;
2230     upper[4]=1.;
2231     lower[5]=0.;
2232     upper[5]=1.;
2233
2234 // step sizes
2235     static Double_t step[6]={0.0005, 0.0005, 0.0005, 0.0005, 0.001, 0.001};
2236     clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
2237     clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
2238     clusterInput.Fitter()->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
2239     clusterInput.Fitter()->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
2240     clusterInput.Fitter()->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
2241     clusterInput.Fitter()->mnparm(5,"a1",vstart[5],step[5],lower[5],upper[5],ierflag);
2242 // ready for minimisation       
2243     arglist[0]= -1;
2244     arglist[1]= 0;
2245     
2246     clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
2247     clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
2248     //    clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
2249 // Get fitted parameters
2250     TString chname;
2251     Double_t epxz, b1, b2;
2252     Int_t ierflg;
2253     clusterInput.Fitter()->mnpout(0, chname, fXFit[0],  epxz, b1, b2, ierflg);  
2254     clusterInput.Fitter()->mnpout(1, chname, fYFit[0],  epxz, b1, b2, ierflg);  
2255     clusterInput.Fitter()->mnpout(2, chname, fXFit[1],  epxz, b1, b2, ierflg);  
2256     clusterInput.Fitter()->mnpout(3, chname, fYFit[1],  epxz, b1, b2, ierflg);  
2257     clusterInput.Fitter()->mnpout(4, chname, fQrFit[0], epxz, b1, b2, ierflg);  
2258     clusterInput.Fitter()->mnpout(5, chname, fQrFit[1], epxz, b1, b2, ierflg);  
2259
2260     Double_t fmin, fedm, errdef;
2261     Int_t   npari, nparx, istat;
2262       
2263     clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);  
2264     fFitStat=istat;
2265     
2266     fChi2[0]=fmin;
2267     fChi2[1]=fmin;
2268     return fmin;
2269 }
2270
2271 void AliMUONClusterFinderVS::Split(AliMUONRawCluster* c)
2272 {
2273 //
2274 // One cluster for each maximum
2275 //
2276     Int_t i, j, cath;
2277     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2278     for (j=0; j<2; j++) {
2279         AliMUONRawCluster cnew;
2280         cnew.SetGhost(c->GetGhost());
2281         for (cath=0; cath<2; cath++) {
2282             cnew.SetChi2(cath,fChi2[0]);
2283             // ?? why not cnew.fChi2[cath]=fChi2[cath];
2284             
2285             if (fNPeaks == 0) {
2286                 cnew.SetNcluster(0,-1);
2287                 cnew.SetNcluster(1,fNRawClusters);
2288             } else {
2289                 cnew.SetNcluster(0,fNPeaks);
2290                 cnew.SetNcluster(1,0);
2291             }
2292             cnew.SetMultiplicity(cath,0);
2293             cnew.SetX(cath, Float_t(fXFit[j]));
2294             cnew.SetY(cath, Float_t(fYFit[j]));
2295             cnew.SetZ(cath, fZPlane);
2296             if (j==0) {
2297                 cnew.SetCharge(cath, Int_t(clusterInput.TotalCharge(cath)*fQrFit[cath]));
2298             } else {
2299                 cnew.SetCharge(cath, Int_t(clusterInput.TotalCharge(cath)*(1-fQrFit[cath])));
2300             }
2301             if (fSegmentationType == 1) 
2302               fSeg[cath]->SetHit(fXFit[j],fYFit[j],fZPlane);
2303             else 
2304               fSeg2[cath]->SetHit(fInput->DetElemId(), fXFit[j],fYFit[j],fZPlane);
2305
2306             for (i=0; i<fMul[cath]; i++) {
2307               Float_t q1;
2308                 cnew.SetIndex(cnew.GetMultiplicity(cath), cath, c->GetIndex(i,cath));
2309                 if (fSegmentationType == 1) {
2310                   fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
2311                   q1 = fInput->Mathieson()->IntXY(fSeg[cath]);
2312                 } else {
2313                   fSeg2[cath]->SetPad(fInput->DetElemId(),fIx[i][cath], fIy[i][cath]);
2314                   q1 = fInput->Mathieson()->IntXY(fInput->DetElemId(),fSeg2[cath]);
2315                 }
2316                 cnew.SetContrib(i, cath, q1*Float_t(cnew.GetCharge(cath))/Float_t(fQ[i][cath]));
2317                 cnew.SetMultiplicity(cath, cnew.GetMultiplicity(cath)+1 );
2318             }
2319             FillCluster(&cnew,0,cath);
2320         } // cathode loop
2321         
2322         cnew.SetClusterType(cnew.PhysicsContribution());
2323         if (cnew.GetCharge(0)>0 && cnew.GetCharge(1)>0) AddRawCluster(cnew);
2324         fNPeaks++;
2325     }
2326 }
2327 void AliMUONClusterFinderVS::AddRawCluster(const AliMUONRawCluster& c)
2328 {
2329   //
2330   // Add a raw cluster copy to the list
2331   //
2332
2333 //     AliMUON *pMUON=(AliMUON*)gAlice->GetModule("MUON");
2334 //     pMUON->GetMUONData()->AddRawCluster(fInput->Chamber(),c); 
2335 //     fNRawClusters++;
2336
2337   
2338     TClonesArray &lrawcl = *fRawClusters;
2339     new(lrawcl[fNRawClusters++]) AliMUONRawCluster(c);
2340     AliDebug(1,Form("\nfNRawClusters %d\n",fNRawClusters));
2341 }
2342
2343 AliMUONClusterFinderVS& AliMUONClusterFinderVS
2344 ::operator = (const AliMUONClusterFinderVS& rhs)
2345 {
2346 // Protected assignement operator
2347
2348   if (this == &rhs) return *this;
2349
2350   AliFatal("Not implemented.");
2351     
2352   return *this;  
2353 }
2354
2355 //
2356 // Minimisation functions
2357 // Single Mathieson
2358 void fcnS1(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2359 {
2360     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());    
2361     Int_t i;
2362     Float_t delta;
2363     Float_t chisq=0;
2364     Float_t qcont=0;
2365     Float_t qtot=0;
2366
2367     for (i=0; i<clusterInput.Nmul(0); i++) {
2368         Float_t q0=clusterInput.Charge(i,0);
2369         Float_t q1=clusterInput.DiscrChargeS1(i,par);
2370         delta=(q0-q1)/q0;
2371         chisq+=delta*delta;
2372         qcont+=q1;
2373         qtot+=q0;
2374     }
2375     f=chisq;
2376 }
2377
2378 void fcnCombiS1(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2379 {
2380     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());    
2381     Int_t i, cath;
2382     Float_t delta;
2383     Float_t chisq=0;
2384     Float_t qcont=0;
2385     Float_t qtot=0;
2386
2387     for (cath=0; cath<2; cath++) {
2388         for (i=0; i<clusterInput.Nmul(cath); i++) {
2389             Float_t q0=clusterInput.Charge(i,cath);
2390             Float_t q1=clusterInput.DiscrChargeCombiS1(i,par,cath);
2391             delta=(q0-q1)/q0;
2392             chisq+=delta*delta;
2393             qcont+=q1;
2394             qtot+=q0;
2395         }
2396     }
2397     f=chisq;
2398 }
2399
2400 // Double Mathieson
2401 void fcnS2(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2402 {
2403     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());    
2404     Int_t i;
2405     Float_t delta;
2406     Float_t chisq=0;
2407     Float_t qcont=0;
2408     Float_t qtot=0;
2409     
2410     for (i=0; i<clusterInput.Nmul(0); i++) {
2411
2412         Float_t q0=clusterInput.Charge(i,0);
2413         Float_t q1=clusterInput.DiscrChargeS2(i,par);
2414         delta=(q0-q1)/q0;
2415         chisq+=delta*delta;
2416         qcont+=q1;
2417         qtot+=q0;
2418     }
2419     f=chisq;
2420 }
2421
2422 // Double Mathieson
2423 void fcnCombiS2(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2424 {
2425     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());    
2426     Int_t i, cath;
2427     Float_t delta;
2428     Float_t chisq=0;
2429     Float_t qcont=0;
2430     Float_t qtot=0;
2431     for (cath=0; cath<2; cath++) {
2432         for (i=0; i<clusterInput.Nmul(cath); i++) {
2433             Float_t q0=clusterInput.Charge(i,cath);
2434             Float_t q1=clusterInput.DiscrChargeCombiS2(i,par,cath);
2435             delta=(q0-q1)/q0;
2436             chisq+=delta*delta;
2437             qcont+=q1;
2438             qtot+=q0;
2439         }
2440     }
2441     f=chisq;
2442 }