MUON/AliMUONClusterFinderVS.cxx

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