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