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 "AliMUONClusterFinderVS.h"
  19 #include "AliMUONDigit.h"
  20 #include "AliMUONRawCluster.h"
  21 #include "AliMUONGeometrySegmentation.h"
  22 #include "AliMUONMathieson.h"
  23 #include "AliMUONClusterInput.h"
  24 #include "AliMUONDigitMapA1.h"
  25
  26 #include "AliLog.h"
  27
  28 #include <TMinuit.h>
  29 #include <TF1.h>
  30 #include <TMinuit.h>
  31 #include <Riostream.h>
  32
  33
  34 //_____________________________________________________________________
  35 // This function is minimized in the double-Mathieson fit
  36 void fcnS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
  37 void fcnS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
  38 void fcnCombiS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
  39 void fcnCombiS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
  40
  41 ClassImp(AliMUONClusterFinderVS)
  42
  43 AliMUONClusterFinderVS::AliMUONClusterFinderVS()
  44   : TObject()
  45 {
  46 // Default constructor
  47     fInput=AliMUONClusterInput::Instance();
  48     fDigitMap[0] = 0;
  49     fDigitMap[1] = 0;
  50     fTrack[0]=fTrack[1]=-1;
  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 (fDigitMap[cath]->TestHit(x[j], y[j])==kEmpty) continue;
 990           digt=(AliMUONDigit*) fDigitMap[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 (fDigitMap[cath]->TestHit(ix, iy)!=kEmpty) {
1065                     iNN++;
1066                     digt=(AliMUONDigit*) fDigitMap[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 (fDigitMap[cath]->TestHit(ix, iy)!=kEmpty) {
1130                     iNN++;
1131                     digt=(AliMUONDigit*) fDigitMap[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 = fDigitMap[cath]->GetHitIndex(i,j);
1303     AliMUONDigit* dig = (AliMUONDigit*) fDigitMap[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     fDigitMap[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 (fDigitMap[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 (fDigitMap[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 (fDigitMap[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     fSeg2[0] = fInput->Segmentation2(0);
1441     fSeg2[1] = fInput->Segmentation2(1);
1442
1443     Int_t detElemId = fInput->DetElemId();
1444
1445     Int_t npx0  = fSeg2[0]->Npx(detElemId)+1;
1446     Int_t npy0  = fSeg2[0]->Npy(detElemId)+1;
1447     fDigitMap[0]  = new AliMUONDigitMapA1(detElemId, npx0, npy0);
1448
1449     Int_t npx1  = fSeg2[0]->Npx(detElemId)+1;
1450     Int_t npy1  = fSeg2[0]->Npy(detElemId)+1;
1451     fDigitMap[1]  = new AliMUONDigitMapA1(detElemId, npx1, npy1);
1452
1453     AliMUONDigit *dig;
1454
1455     Int_t ndig, cath;
1456     Int_t nskip=0;
1457     Int_t ncls=0;
1458
1459     fDigitMap[0]->FillHits(fInput->Digits(0));
1460     fDigitMap[1]->FillHits(fInput->Digits(1));
1461 //
1462 //  Outer Loop over Cathodes
1463     for (cath = 0; cath < 2; cath++) {
1464
1465         for (ndig=0; ndig<fInput->NDigits(cath); ndig++) {
1466           dig = fInput->Digit(cath, ndig);
1467           Int_t padx = dig->PadX();
1468           Int_t pady = dig->PadY();
1469           if (fDigitMap[cath]->TestHit(padx,pady)==kUsed ||fDigitMap[0]->TestHit(padx,pady)==kEmpty) {
1470             nskip++;
1471             continue;
1472           }
1473           AliDebug(1,Form("\n CATHODE %d CLUSTER %d\n",cath,ncls));
1474           AliMUONRawCluster clus;
1475           clus.SetMultiplicity(0, 0);
1476           clus.SetMultiplicity(1, 0);
1477           clus.SetPeakSignal(cath,dig->Signal());
1478           clus.SetTrack(0, dig->Hit());
1479           clus.SetTrack(1, dig->Track(0));
1480           clus.SetTrack(2, dig->Track(1));
1481
1482           AliDebug(1,Form("idDE %d Padx %d Pady %d", fInput->DetElemId(), padx, pady));
1483
1484           // tag the beginning of cluster list in a raw cluster
1485           clus.SetNcluster(0,-1);
1486           Float_t xcu, ycu;
1487           fSeg2[cath]->GetPadC(fInput->DetElemId(), padx, pady, xcu, ycu, fZPlane);
1488           fSector= fSeg2[cath]->Sector(fInput->DetElemId(), padx, pady)/100;
1489
1490
1491
1492
1493           FindCluster(padx,pady,cath,clus);
1494           //^^^^^^^^^^^^^^^^^^^^^^^^
1495           // center of gravity
1496           if (clus.GetX(0)!=0.) clus.SetX(0, clus.GetX(0)/clus.GetCharge(0)); // clus.fX[0] /= clus.fQ[0];
1497
1498           // Force on anod
1499           clus.SetX(0,fSeg2[0]->GetAnod(fInput->DetElemId(), clus.GetX(0)));
1500           if (clus.GetY(0)!=0.) clus.SetY(0, clus.GetY(0)/clus.GetCharge(0)); // clus.fY[0] /= clus.fQ[0];
1501
1502           if(clus.GetCharge(1)!=0.) clus.SetX(1, clus.GetX(1)/clus.GetCharge(1));  // clus.fX[1] /= clus.fQ[1];
1503
1504           // Force on anod
1505           clus.SetX(1, fSeg2[0]->GetAnod(fInput->DetElemId(),clus.GetX(1)));
1506           if(clus.GetCharge(1)!=0.) clus.SetY(1, clus.GetY(1)/clus.GetCharge(1));// clus.fY[1] /= clus.fQ[1];
1507
1508           clus.SetZ(0, fZPlane);
1509           clus.SetZ(1, fZPlane);
1510
1511           AliDebug(1,Form("\n Cathode 1 multiplicite %d X(CG) %f Y(CG) %f\n",
1512                           clus.GetMultiplicity(0),clus.GetX(0),clus.GetY(0)));
1513           AliDebug(1,Form(" Cathode 2 multiplicite %d X(CG) %f Y(CG) %f\n",
1514                           clus.GetMultiplicity(1),clus.GetX(1),clus.GetY(1)));
1515           //      Analyse cluster and decluster if necessary
1516           //
1517           ncls++;
1518           clus.SetNcluster(1,fNRawClusters);
1519           clus.SetClusterType(clus.PhysicsContribution());
1520
1521           fNPeaks=0;
1522           //
1523           //
1524           Decluster(&clus);
1525           //
1526           //      reset Cluster object
1527           { // begin local scope
1528             for (int k=0;k<clus.GetMultiplicity(0);k++) clus.SetIndex(k, 0, 0);
1529           } // end local scope
1530
1531           { // begin local scope
1532             for (int k=0;k<clus.GetMultiplicity(1);k++) clus.SetIndex(k, 1, 0);
1533           } // end local scope
1534
1535           clus.SetMultiplicity(0,0);
1536           clus.SetMultiplicity(1,0);
1537
1538
1539         } // end loop ndig
1540     } // end loop cathodes
1541     delete fDigitMap[0];
1542     delete fDigitMap[1];
1543 }
1544
1545 Float_t AliMUONClusterFinderVS::SingleMathiesonFit(AliMUONRawCluster *c, Int_t cath)
1546 {
1547 // Performs a single Mathieson fit on one cathode
1548 //
1549     Double_t arglist[20];
1550     Int_t ierflag=0;
1551     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1552
1553     clusterInput.Fitter()->SetFCN(fcnS1);
1554     clusterInput.Fitter()->mninit(2,10,7);
1555     clusterInput.Fitter()->SetPrintLevel(-1 + AliLog::GetGlobalDebugLevel());
1556     arglist[0]=-1;
1557     clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1558 // Set starting values
1559     static Double_t vstart[2];
1560     vstart[0]=c->GetX(1);
1561     vstart[1]=c->GetY(0);
1562
1563
1564 // lower and upper limits
1565     static Double_t lower[2], upper[2];
1566     Int_t ix,iy, isec;
1567     fSeg2[cath]->GetPadI(fInput->DetElemId(), c->GetX(cath), c->GetY(cath), fZPlane, ix, iy);
1568     isec=fSeg2[cath]->Sector(fInput->DetElemId(), ix, iy);
1569
1570     lower[0]=vstart[0]-fSeg2[cath]->Dpx(fInput->DetElemId(), isec)/2;
1571     lower[1]=vstart[1]-fSeg2[cath]->Dpy(fInput->DetElemId(), isec)/2;
1572
1573     upper[0]=lower[0]+fSeg2[cath]->Dpx(fInput->DetElemId(), isec);
1574     upper[1]=lower[1]+fSeg2[cath]->Dpy(fInput->DetElemId(), isec);
1575
1576
1577 // step sizes
1578     static Double_t step[2]={0.0005, 0.0005};
1579
1580     clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1581     clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1582 // ready for minimisation
1583     arglist[0]= -1;
1584     arglist[1]= 0;
1585
1586     clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1587     clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1588     //    clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1589     Double_t fmin, fedm, errdef;
1590     Int_t   npari, nparx, istat;
1591
1592     clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1593     fFitStat=istat;
1594
1595 // Print results
1596 // Get fitted parameters
1597     Double_t xrec, yrec;
1598     TString chname;
1599     Double_t epxz, b1, b2;
1600     Int_t ierflg;
1601     clusterInput.Fitter()->mnpout(0, chname, xrec, epxz, b1, b2, ierflg);
1602     clusterInput.Fitter()->mnpout(1, chname, yrec, epxz, b1, b2, ierflg);
1603     fXFit[cath]=xrec;
1604     fYFit[cath]=yrec;
1605     return fmin;
1606 }
1607
1608 Float_t AliMUONClusterFinderVS::CombiSingleMathiesonFit(AliMUONRawCluster * /*c*/)
1609 {
1610 // Perform combined Mathieson fit on both cathode planes
1611 //
1612     Double_t arglist[20];
1613     Int_t ierflag=0;
1614     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1615     clusterInput.Fitter()->SetFCN(fcnCombiS1);
1616     clusterInput.Fitter()->mninit(2,10,7);
1617     clusterInput.Fitter()->SetPrintLevel(-1 + AliLog::GetGlobalDebugLevel());
1618     arglist[0]=-1;
1619     clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1620     static Double_t vstart[2];
1621     vstart[0]=fXInit[0];
1622     vstart[1]=fYInit[0];
1623
1624
1625 // lower and upper limits
1626     static Float_t lower[2], upper[2];
1627     Int_t ix,iy,isec;
1628     Float_t dpy, dpx;
1629
1630     fSeg2[0]->GetPadI(fInput->DetElemId(), fXInit[0], fYInit[0], fZPlane, ix, iy);
1631     isec=fSeg2[0]->Sector(fInput->DetElemId(), ix, iy);
1632     dpy=fSeg2[0]->Dpy(fInput->DetElemId(), isec);
1633     fSeg2[1]->GetPadI(fInput->DetElemId(), fXInit[0], fYInit[0], fZPlane, ix, iy);
1634     isec=fSeg2[1]->Sector(fInput->DetElemId(), ix, iy);
1635     dpx=fSeg2[1]->Dpx(fInput->DetElemId(), isec);
1636
1637     Int_t icount;
1638     Float_t xdum, ydum, zdum;
1639
1640 //  Find save upper and lower limits
1641
1642     icount = 0;
1643     for (fSeg2[1]->FirstPad(fInput->DetElemId(),fXInit[0], fYInit[0], fZPlane, dpx, 0.);
1644          fSeg2[1]->MorePads(fInput->DetElemId());
1645          fSeg2[1]->NextPad(fInput->DetElemId()))
1646         {
1647           ix=fSeg2[1]->Ix(); iy=fSeg2[1]->Iy();
1648           fSeg2[1]->GetPadC(fInput->DetElemId(), ix,iy, upper[0], ydum, zdum);
1649           if (icount ==0) lower[0]=upper[0];
1650           icount++;
1651         }
1652
1653     if (lower[0]>upper[0]) {xdum=lower[0]; lower[0]=upper[0]; upper[0]=xdum;}
1654
1655     icount=0;
1656     AliDebug(1,Form("\n single y %f %f", fXInit[0], fYInit[0]));
1657
1658     for (fSeg2[0]->FirstPad(fInput->DetElemId(), fXInit[0], fYInit[0], fZPlane, 0., dpy);
1659          fSeg2[0]->MorePads(fInput->DetElemId());
1660          fSeg2[0]->NextPad(fInput->DetElemId()))
1661         {
1662           ix=fSeg2[0]->Ix(); iy=fSeg2[0]->Iy();
1663           fSeg2[0]->GetPadC(fInput->DetElemId(), ix,iy,xdum,upper[1],zdum);
1664           if (icount ==0) lower[1]=upper[1];
1665           icount++;
1666           AliDebug(1,Form("\n upper lower %d %f %f", icount, upper[1], lower[1]));
1667         }
1668
1669     if (lower[1]>upper[1]) {xdum=lower[1]; lower[1]=upper[1]; upper[1]=xdum;}
1670
1671 // step sizes
1672     static Double_t step[2]={0.00001, 0.0001};
1673
1674     clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1675     clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1676 // ready for minimisation
1677     arglist[0]= -1;
1678     arglist[1]= 0;
1679
1680     clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1681     clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1682     //    clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1683     Double_t fmin, fedm, errdef;
1684     Int_t   npari, nparx, istat;
1685
1686     clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1687     fFitStat=istat;
1688
1689 // Print results
1690 // Get fitted parameters
1691     Double_t xrec, yrec;
1692     TString chname;
1693     Double_t epxz, b1, b2;
1694     Int_t ierflg;
1695     clusterInput.Fitter()->mnpout(0, chname, xrec, epxz, b1, b2, ierflg);
1696     clusterInput.Fitter()->mnpout(1, chname, yrec, epxz, b1, b2, ierflg);
1697     fXFit[0]=xrec;
1698     fYFit[0]=yrec;
1699     return fmin;
1700 }
1701
1702 Bool_t AliMUONClusterFinderVS::DoubleMathiesonFit(AliMUONRawCluster * /*c*/, Int_t cath)
1703 {
1704 // Performs a double Mathieson fit on one cathode
1705 //
1706
1707 //
1708 //  Initialise global variables for fit
1709     Double_t arglist[20];
1710     Int_t ierflag=0;
1711     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1712     clusterInput.Fitter()->SetFCN(fcnS2);
1713     clusterInput.Fitter()->mninit(5,10,7);
1714     clusterInput.Fitter()->SetPrintLevel(-1 + AliLog::GetGlobalDebugLevel());
1715     arglist[0]=-1;
1716     clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1717 // Set starting values
1718     static Double_t vstart[5];
1719     vstart[0]=fX[fIndLocal[0][cath]][cath];
1720     vstart[1]=fY[fIndLocal[0][cath]][cath];
1721     vstart[2]=fX[fIndLocal[1][cath]][cath];
1722     vstart[3]=fY[fIndLocal[1][cath]][cath];
1723     vstart[4]=Float_t(fQ[fIndLocal[0][cath]][cath])/
1724         Float_t(fQ[fIndLocal[0][cath]][cath]+fQ[fIndLocal[1][cath]][cath]);
1725 // lower and upper limits
1726     static Float_t lower[5], upper[5];
1727     Int_t isec;
1728
1729     isec=fSeg2[cath]->Sector(fInput->DetElemId(),fIx[fIndLocal[0][cath]][cath],
1730                              fIy[fIndLocal[0][cath]][cath]);
1731     lower[0]=vstart[0]-fSeg2[cath]->Dpx(fInput->DetElemId(),isec);
1732     lower[1]=vstart[1]-fSeg2[cath]->Dpy(fInput->DetElemId(),isec);
1733
1734     upper[0]=lower[0]+2.*fSeg2[cath]->Dpx(fInput->DetElemId(),isec);
1735     upper[1]=lower[1]+2.*fSeg2[cath]->Dpy(fInput->DetElemId(),isec);
1736
1737     isec=fSeg2[cath]->Sector(fInput->DetElemId(),fIx[fIndLocal[1][cath]][cath],
1738                              fIy[fIndLocal[1][cath]][cath]);
1739     lower[2]=vstart[2]-fSeg2[cath]->Dpx(fInput->DetElemId(),isec)/2;
1740     lower[3]=vstart[3]-fSeg2[cath]->Dpy(fInput->DetElemId(),isec)/2;
1741
1742     upper[2]=lower[2]+fSeg2[cath]->Dpx(fInput->DetElemId(),isec);
1743     upper[1]=lower[1]+2.*fSeg2[cath]->Dpy(fInput->DetElemId(),isec);
1744
1745
1746
1747     lower[4]=0.;
1748     upper[4]=1.;
1749 // step sizes
1750     static Double_t step[5]={0.0005, 0.0005, 0.0005, 0.0005, 0.0001};
1751
1752     clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1753     clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1754     clusterInput.Fitter()->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
1755     clusterInput.Fitter()->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
1756     clusterInput.Fitter()->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
1757 // ready for minimisation
1758     arglist[0]= -1;
1759     arglist[1]= 0;
1760
1761     clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1762     clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1763     //    clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1764 // Get fitted parameters
1765     Double_t xrec[2], yrec[2], qfrac;
1766     TString chname;
1767     Double_t epxz, b1, b2;
1768     Int_t ierflg;
1769     clusterInput.Fitter()->mnpout(0, chname, xrec[0], epxz, b1, b2, ierflg);
1770     clusterInput.Fitter()->mnpout(1, chname, yrec[0], epxz, b1, b2, ierflg);
1771     clusterInput.Fitter()->mnpout(2, chname, xrec[1], epxz, b1, b2, ierflg);
1772     clusterInput.Fitter()->mnpout(3, chname, yrec[1], epxz, b1, b2, ierflg);
1773     clusterInput.Fitter()->mnpout(4, chname, qfrac,   epxz, b1, b2, ierflg);
1774
1775     Double_t fmin, fedm, errdef;
1776     Int_t   npari, nparx, istat;
1777
1778     clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1779     fFitStat=istat;
1780     return kTRUE;
1781 }
1782
1783 Float_t AliMUONClusterFinderVS::CombiDoubleMathiesonFit(AliMUONRawCluster * /*c*/)
1784 {
1785 //
1786 // Perform combined double Mathieson fit on both cathode planes
1787 //
1788     Double_t arglist[20];
1789     Int_t ierflag=0;
1790     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1791     clusterInput.Fitter()->SetFCN(fcnCombiS2);
1792     clusterInput.Fitter()->mninit(6,10,7);
1793     clusterInput.Fitter()->SetPrintLevel(-1 + AliLog::GetGlobalDebugLevel());
1794     arglist[0]=-1;
1795     clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1796 // Set starting values
1797     static Double_t vstart[6];
1798     vstart[0]=fXInit[0];
1799     vstart[1]=fYInit[0];
1800     vstart[2]=fXInit[1];
1801     vstart[3]=fYInit[1];
1802     vstart[4]=fQrInit[0];
1803     vstart[5]=fQrInit[1];
1804 // lower and upper limits
1805     static Float_t lower[6], upper[6];
1806     Int_t ix,iy,isec;
1807     Float_t dpx, dpy;
1808
1809     fSeg2[1]->GetPadI(fInput->DetElemId(),fXInit[0], fYInit[0], fZPlane, ix, iy);
1810     isec=fSeg2[1]->Sector(fInput->DetElemId(),ix, iy);
1811     dpx=fSeg2[1]->Dpx(fInput->DetElemId(), isec);
1812
1813     fSeg2[0]->GetPadI(fInput->DetElemId(), fXInit[0], fYInit[0], fZPlane, ix, iy);
1814     isec=fSeg2[0]->Sector(fInput->DetElemId(), ix, iy);
1815     dpy=fSeg2[0]->Dpy(fInput->DetElemId(), isec);
1816
1817
1818
1819     Int_t icount;
1820     Float_t xdum, ydum, zdum;
1821     AliDebug(1,Form("\n Cluster Finder: %f %f %f %f  ", fXInit[0], fXInit[1],fYInit[0], fYInit[1] ));
1822
1823     //  Find save upper and lower limits
1824     icount = 0;
1825
1826     for (fSeg2[1]->FirstPad(fInput->DetElemId(),fXInit[0], fYInit[0], fZPlane, dpx, 0.);
1827          fSeg2[1]->MorePads(fInput->DetElemId());
1828          fSeg2[1]->NextPad(fInput->DetElemId()))
1829       {
1830         ix=fSeg2[1]->Ix(); iy=fSeg2[1]->Iy();
1831         //      if (fDigitMap[1]->TestHit(ix, iy) == kEmpty) continue;
1832         fSeg2[1]->GetPadC(fInput->DetElemId(),ix,iy,upper[0],ydum,zdum);
1833         if (icount ==0) lower[0]=upper[0];
1834         icount++;
1835       }
1836     if (lower[0]>upper[0]) {xdum=lower[0]; lower[0]=upper[0]; upper[0]=xdum;}
1837     //    vstart[0] = 0.5*(lower[0]+upper[0]);
1838
1839
1840     icount=0;
1841
1842     for (fSeg2[0]->FirstPad(fInput->DetElemId(),fXInit[0], fYInit[0], fZPlane, 0., dpy);
1843          fSeg2[0]->MorePads(fInput->DetElemId());
1844          fSeg2[0]->NextPad(fInput->DetElemId()))
1845       {
1846         ix=fSeg2[0]->Ix(); iy=fSeg2[0]->Iy();
1847         //      if (fDigitMap[0]->TestHit(ix, iy) == kEmpty) continue;
1848         fSeg2[0]->GetPadC(fInput->DetElemId(),ix,iy,xdum,upper[1],zdum);
1849         if (icount ==0) lower[1]=upper[1];
1850         icount++;
1851       }
1852
1853     if (lower[1]>upper[1]) {xdum=lower[1]; lower[1]=upper[1]; upper[1]=xdum;}
1854     //     vstart[1] = 0.5*(lower[1]+upper[1]);
1855
1856
1857     fSeg2[1]->GetPadI(fInput->DetElemId(),fXInit[1], fYInit[1], fZPlane, ix, iy);
1858     isec=fSeg2[1]->Sector(fInput->DetElemId(),ix, iy);
1859     dpx=fSeg2[1]->Dpx(fInput->DetElemId(),isec);
1860     fSeg2[0]->GetPadI(fInput->DetElemId(),fXInit[1], fYInit[1], fZPlane, ix, iy);
1861     isec=fSeg2[0]->Sector(fInput->DetElemId(),ix, iy);
1862     dpy=fSeg2[0]->Dpy(fInput->DetElemId(),isec);
1863
1864
1865     //  Find save upper and lower limits
1866
1867     icount=0;
1868
1869     for (fSeg2[1]->FirstPad(fInput->DetElemId(),fXInit[1], fYInit[1], fZPlane, dpx, 0);
1870          fSeg2[1]->MorePads(fInput->DetElemId());
1871          fSeg2[1]->NextPad(fInput->DetElemId()))
1872       {
1873         ix=fSeg2[1]->Ix(); iy=fSeg2[1]->Iy();
1874         //      if (fDigitMap[1]->TestHit(ix, iy) == kEmpty) continue;
1875         fSeg2[1]->GetPadC(fInput->DetElemId(),ix,iy,upper[2],ydum,zdum);
1876         if (icount ==0) lower[2]=upper[2];
1877         icount++;
1878       }
1879     if (lower[2]>upper[2]) {xdum=lower[2]; lower[2]=upper[2]; upper[2]=xdum;}
1880     //    vstart[2] = 0.5*(lower[2]+upper[2]);
1881
1882     icount=0;
1883
1884     for (fSeg2[0]->FirstPad(fInput->DetElemId(),fXInit[1], fYInit[1], fZPlane, 0, dpy);
1885          fSeg2[0]-> MorePads(fInput->DetElemId());
1886          fSeg2[0]->NextPad(fInput->DetElemId()))
1887       {
1888         ix=fSeg2[0]->Ix(); iy=fSeg2[0]->Iy();
1889         //      if (fDigitMap[0]->TestHit(ix, iy) != kEmpty) continue;
1890
1891         fSeg2[0]->GetPadC(fInput->DetElemId(),ix,iy,xdum,upper[3],zdum);
1892         if (icount ==0) lower[3]=upper[3];
1893         icount++;
1894
1895       }
1896     if (lower[3]>upper[3]) {xdum=lower[3]; lower[3]=upper[3]; upper[3]=xdum;}
1897
1898     lower[4]=0.;
1899     upper[4]=1.;
1900     lower[5]=0.;
1901     upper[5]=1.;
1902
1903 // step sizes
1904     static Double_t step[6]={0.0005, 0.0005, 0.0005, 0.0005, 0.001, 0.001};
1905     clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1906     clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1907     clusterInput.Fitter()->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
1908     clusterInput.Fitter()->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
1909     clusterInput.Fitter()->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
1910     clusterInput.Fitter()->mnparm(5,"a1",vstart[5],step[5],lower[5],upper[5],ierflag);
1911 // ready for minimisation
1912     arglist[0]= -1;
1913     arglist[1]= 0;
1914
1915     clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1916     clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1917     //    clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1918 // Get fitted parameters
1919     TString chname;
1920     Double_t epxz, b1, b2;
1921     Int_t ierflg;
1922     clusterInput.Fitter()->mnpout(0, chname, fXFit[0],  epxz, b1, b2, ierflg);
1923     clusterInput.Fitter()->mnpout(1, chname, fYFit[0],  epxz, b1, b2, ierflg);
1924     clusterInput.Fitter()->mnpout(2, chname, fXFit[1],  epxz, b1, b2, ierflg);
1925     clusterInput.Fitter()->mnpout(3, chname, fYFit[1],  epxz, b1, b2, ierflg);
1926     clusterInput.Fitter()->mnpout(4, chname, fQrFit[0], epxz, b1, b2, ierflg);
1927     clusterInput.Fitter()->mnpout(5, chname, fQrFit[1], epxz, b1, b2, ierflg);
1928
1929     Double_t fmin, fedm, errdef;
1930     Int_t   npari, nparx, istat;
1931
1932     clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1933     fFitStat=istat;
1934
1935     fChi2[0]=fmin;
1936     fChi2[1]=fmin;
1937     return fmin;
1938 }
1939
1940 void AliMUONClusterFinderVS::Split(AliMUONRawCluster* c)
1941 {
1942 //
1943 // One cluster for each maximum
1944 //
1945     Int_t i, j, cath;
1946     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1947     for (j=0; j<2; j++) {
1948         AliMUONRawCluster cnew;
1949         cnew.SetGhost(c->GetGhost());
1950         for (cath=0; cath<2; cath++) {
1951             cnew.SetChi2(cath,fChi2[0]);
1952             // ?? why not cnew.fChi2[cath]=fChi2[cath];
1953
1954             if (fNPeaks == 0) {
1955                 cnew.SetNcluster(0,-1);
1956                 cnew.SetNcluster(1,fNRawClusters);
1957             } else {
1958                 cnew.SetNcluster(0,fNPeaks);
1959                 cnew.SetNcluster(1,0);
1960             }
1961             cnew.SetMultiplicity(cath,0);
1962             cnew.SetX(cath, Float_t(fXFit[j]));
1963             cnew.SetY(cath, Float_t(fYFit[j]));
1964             cnew.SetZ(cath, fZPlane);
1965             if (j==0) {
1966                 cnew.SetCharge(cath, Int_t(clusterInput.TotalCharge(cath)*fQrFit[cath]));
1967             } else {
1968                 cnew.SetCharge(cath, Int_t(clusterInput.TotalCharge(cath)*(1-fQrFit[cath])));
1969             }
1970             fSeg2[cath]->SetHit(fInput->DetElemId(), fXFit[j],fYFit[j],fZPlane);
1971
1972             for (i=0; i<fMul[cath]; i++) {
1973               Float_t q1;
1974                 cnew.SetIndex(cnew.GetMultiplicity(cath), cath, c->GetIndex(i,cath));
1975
1976                 fSeg2[cath]->SetPad(fInput->DetElemId(),fIx[i][cath], fIy[i][cath]);
1977                 q1 = fInput->Mathieson()->IntXY(fInput->DetElemId(),fSeg2[cath]);
1978
1979                 cnew.SetContrib(i, cath, q1*Float_t(cnew.GetCharge(cath))/Float_t(fQ[i][cath]));
1980                 cnew.SetMultiplicity(cath, cnew.GetMultiplicity(cath)+1 );
1981             }
1982             FillCluster(&cnew,0,cath);
1983         } // cathode loop
1984         cnew.SetClusterType(cnew.PhysicsContribution());
1985         if (cnew.GetCharge(0)>0 && cnew.GetCharge(1)>0) AddRawCluster(cnew);
1986         fNPeaks++;
1987     }
1988 }
1989 void AliMUONClusterFinderVS::AddRawCluster(AliMUONRawCluster& c)
1990 {
1991   //
1992   // Add a raw cluster copy to the list
1993   //
1994   //     AliMUON *pMUON=(AliMUON*)gAlice->GetModule("MUON");
1995   //     pMUON->GetMUONData()->AddRawCluster(fInput->Chamber(),c);
1996   //     fNRawClusters++;
1997
1998   // Setting detection element in raw cluster for alignment
1999   // BB 19/05/05
2000   c.SetDetElemId(fInput->DetElemId());
2001
2002   TClonesArray &lrawcl = *fRawClusters;
2003   new(lrawcl[fNRawClusters++]) AliMUONRawCluster(c);
2004   AliDebug(1,Form("\nfNRawClusters %d\n",fNRawClusters));
2005 }
2006
2007 AliMUONClusterFinderVS& AliMUONClusterFinderVS
2008 ::operator = (const AliMUONClusterFinderVS& rhs)
2009 {
2010 // Protected assignement operator
2011
2012   if (this == &rhs) return *this;
2013
2014   AliFatal("Not implemented.");
2015
2016   return *this;
2017 }
2018
2019 //
2020 // Minimisation functions
2021 // Single Mathieson
2022 void fcnS1(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2023 {
2024     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2025     Int_t i;
2026     Float_t delta;
2027     Float_t chisq=0;
2028     Float_t qcont=0;
2029     Float_t qtot=0;
2030
2031     for (i=0; i<clusterInput.Nmul(0); i++) {
2032         Float_t q0=clusterInput.Charge(i,0);
2033         Float_t q1=clusterInput.DiscrChargeS1(i,par);
2034         delta=(q0-q1)/q0;
2035         chisq+=delta*delta;
2036         qcont+=q1;
2037         qtot+=q0;
2038     }
2039     f=chisq;
2040 }
2041
2042 void fcnCombiS1(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2043 {
2044     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2045     Int_t i, cath;
2046     Float_t delta;
2047     Float_t chisq=0;
2048     Float_t qcont=0;
2049     Float_t qtot=0;
2050
2051     for (cath=0; cath<2; cath++) {
2052         for (i=0; i<clusterInput.Nmul(cath); i++) {
2053             Float_t q0=clusterInput.Charge(i,cath);
2054             Float_t q1=clusterInput.DiscrChargeCombiS1(i,par,cath);
2055             delta=(q0-q1)/q0;
2056             chisq+=delta*delta;
2057             qcont+=q1;
2058             qtot+=q0;
2059         }
2060     }
2061     f=chisq;
2062 }
2063
2064 // Double Mathieson
2065 void fcnS2(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2066 {
2067     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2068     Int_t i;
2069     Float_t delta;
2070     Float_t chisq=0;
2071     Float_t qcont=0;
2072     Float_t qtot=0;
2073
2074     for (i=0; i<clusterInput.Nmul(0); i++) {
2075
2076         Float_t q0=clusterInput.Charge(i,0);
2077         Float_t q1=clusterInput.DiscrChargeS2(i,par);
2078         delta=(q0-q1)/q0;
2079         chisq+=delta*delta;
2080         qcont+=q1;
2081         qtot+=q0;
2082     }
2083     f=chisq;
2084 }
2085
2086 // Double Mathieson
2087 void fcnCombiS2(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2088 {
2089     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2090     Int_t i, cath;
2091     Float_t delta;
2092     Float_t chisq=0;
2093     Float_t qcont=0;
2094     Float_t qtot=0;
2095     for (cath=0; cath<2; cath++) {
2096         for (i=0; i<clusterInput.Nmul(cath); i++) {
2097             Float_t q0=clusterInput.Charge(i,cath);
2098             Float_t q1=clusterInput.DiscrChargeCombiS2(i,par,cath);
2099             delta=(q0-q1)/q0;
2100             chisq+=delta*delta;
2101             qcont+=q1;
2102             qtot+=q0;
2103         }
2104     }
2105     f=chisq;
2106 }