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