MUON/AliMUONClusterFinderVS.cxx

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