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