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