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