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