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