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