MUON/AliMUONClusterFinderVS.cxx

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