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