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