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