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