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