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