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