MUON/AliMUONTrackK.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 /* $Id$ */
  17
  18 // ---------------------
  19 // Class AliMUONTrackK
  20 // ---------------------
  21 // Reconstructed track in the muons system based on the extended
  22 // Kalman filter approach
  23 // Author: Alexander Zinchenko, JINR Dubna
  24
  25 #include <Riostream.h>
  26 #include <TClonesArray.h>
  27 #include <TArrayD.h>
  28 #include <TMatrixD.h>
  29
  30 #include "AliMUONTrackK.h"
  31 #include "AliMUON.h"
  32
  33 #include "AliMUONTrackReconstructor.h"
  34 #include "AliMagF.h"
  35 #include "AliMUONSegment.h"
  36 #include "AliMUONHitForRec.h"
  37 #include "AliMUONRawCluster.h"
  38 #include "AliMUONTrackParam.h"
  39 #include "AliMUONEventRecoCombi.h"
  40 #include "AliMUONDetElement.h"
  41 #include "AliRun.h"
  42 #include "AliLog.h"
  43
  44 const Int_t AliMUONTrackK::fgkSize = 5;
  45 const Int_t AliMUONTrackK::fgkNSigma = 12;
  46 const Double_t AliMUONTrackK::fgkChi2max = 100;
  47 const Int_t AliMUONTrackK::fgkTriesMax = 10000;
  48 const Double_t AliMUONTrackK::fgkEpsilon = 0.002;
  49
  50 void mnvertLocal(Double_t* a, Int_t l, Int_t m, Int_t n, Int_t& ifail); // from AliMUONTrack
  51
  52   Int_t AliMUONTrackK::fgDebug = -1; //-1;
  53 Int_t AliMUONTrackK::fgNOfPoints = 0;
  54 AliMUONTrackReconstructor* AliMUONTrackK::fgTrackReconstructor = NULL;
  55 TClonesArray* AliMUONTrackK::fgHitForRec = NULL;
  56 AliMUONEventRecoCombi *AliMUONTrackK::fgCombi = NULL;
  57
  58 ClassImp(AliMUONTrackK) // Class implementation in ROOT context
  59
  60   //__________________________________________________________________________
  61 AliMUONTrackK::AliMUONTrackK()
  62   : AliMUONTrack(),
  63     fStartSegment(0x0),
  64     fPosition(0.),
  65     fPositionNew(0.),
  66     fChi2(0.),
  67     fTrackHits(0x0),
  68     fNmbTrackHits(0),
  69     fTrackDir(1),
  70     fBPFlag(kFALSE),
  71     fRecover(0),
  72     fSkipHit(0x0),
  73     fTrackPar(0x0),
  74     fTrackParNew(0x0),
  75     fCovariance(0x0),
  76     fWeight(0x0),
  77     fParExtrap(0x0),
  78     fParFilter(0x0),
  79     fParSmooth(0x0),
  80     fCovExtrap(0x0),
  81     fCovFilter(0x0),
  82     fJacob(0x0),
  83     fNSteps(0),
  84     fSteps(0x0),
  85     fChi2Array(0x0),
  86     fChi2Smooth(0x0)
  87 {
  88 /// Default constructor
  89
  90   fgTrackReconstructor = NULL; // pointer to event reconstructor
  91   fgHitForRec = NULL; // pointer to points
  92   fgNOfPoints = 0; // number of points
  93
  94   return;
  95 }
  96
  97   //__________________________________________________________________________
  98 AliMUONTrackK::AliMUONTrackK(AliMUONTrackReconstructor *TrackReconstructor, TClonesArray *hitForRec)
  99   : AliMUONTrack(),
 100     fStartSegment(0x0),
 101     fPosition(0.),
 102     fPositionNew(0.),
 103     fChi2(0.),
 104     fTrackHits(0x0),
 105     fNmbTrackHits(0),
 106     fTrackDir(1),
 107     fBPFlag(kFALSE),
 108     fRecover(0),
 109     fSkipHit(0x0),
 110     fTrackPar(0x0),
 111     fTrackParNew(0x0),
 112     fCovariance(0x0),
 113     fWeight(0x0),
 114     fParExtrap(0x0),
 115     fParFilter(0x0),
 116     fParSmooth(0x0),
 117     fCovExtrap(0x0),
 118     fCovFilter(0x0),
 119     fJacob(0x0),
 120     fNSteps(0),
 121     fSteps(0x0),
 122     fChi2Array(0x0),
 123     fChi2Smooth(0x0)
 124 {
 125 /// Constructor
 126
 127   if (!TrackReconstructor) return;
 128   fgTrackReconstructor = TrackReconstructor; // pointer to event reconstructor
 129   fgHitForRec = hitForRec; // pointer to points
 130   fgNOfPoints = fgHitForRec->GetEntriesFast(); // number of points
 131   fgCombi = NULL;
 132   if (fgTrackReconstructor->GetTrackMethod() == 3) fgCombi = AliMUONEventRecoCombi::Instance();
 133 }
 134
 135   //__________________________________________________________________________
 136 AliMUONTrackK::AliMUONTrackK(AliMUONSegment *segment)
 137   //: AliMUONTrack(segment, segment, fgTrackReconstructor)
 138   : AliMUONTrack(NULL, segment, fgTrackReconstructor),
 139     fStartSegment(segment),
 140     fPosition(0.),
 141     fPositionNew(0.),
 142     fChi2(0.),
 143     fTrackHits(new TObjArray(13)),
 144     fNmbTrackHits(2),
 145     fTrackDir(1),
 146     fBPFlag(kFALSE),
 147     fRecover(0),
 148     fSkipHit(0x0),
 149     fTrackPar(new TMatrixD(fgkSize,1)),
 150     fTrackParNew(new TMatrixD(fgkSize,1)),
 151     fCovariance(new TMatrixD(fgkSize,fgkSize)),
 152     fWeight(new TMatrixD(fgkSize,fgkSize)),
 153     fParExtrap(new TObjArray(15)),
 154     fParFilter(new TObjArray(15)),
 155     fParSmooth(0x0),
 156     fCovExtrap(new TObjArray(15)),
 157     fCovFilter(new TObjArray(15)),
 158     fJacob(new TObjArray(15)),
 159     fNSteps(0),
 160     fSteps(new TArrayD(15)),
 161     fChi2Array(new TArrayD(13)),
 162     fChi2Smooth(0x0)
 163 {
 164 /// Constructor from a segment
 165   Double_t dX, dY, dZ;
 166   AliMUONHitForRec *hit1, *hit2;
 167   AliMUONRawCluster *clus;
 168   TClonesArray *rawclusters;
 169
 170   // Pointers to hits from the segment
 171   hit1 = segment->GetHitForRec1();
 172   hit2 = segment->GetHitForRec2();
 173   hit1->SetNTrackHits(hit1->GetNTrackHits()+1); // mark hit as being on track
 174   hit2->SetNTrackHits(hit2->GetNTrackHits()+1); // mark hit as being on track
 175   // check sorting in Z
 176   if (TMath::Abs(hit1->GetZ()) > TMath::Abs(hit2->GetZ())) {
 177     hit1 = hit2;
 178     hit2 = segment->GetHitForRec1();
 179   }
 180
 181   // Fill array of track parameters
 182   if (hit1->GetChamberNumber() > 7) {
 183     // last tracking station
 184     (*fTrackPar)(0,0) = hit1->GetBendingCoor(); // y
 185     (*fTrackPar)(1,0) = hit1->GetNonBendingCoor(); // x
 186     fPosition = hit1->GetZ(); // z
 187     fTrackHits->Add(hit2); // add hit 2
 188     fTrackHits->Add(hit1); // add hit 1
 189     //AZ(Z->-Z) fTrackDir = -1;
 190     fTrackDir = 1;
 191   } else {
 192     // last but one tracking station
 193     (*fTrackPar)(0,0) = hit2->GetBendingCoor(); // y
 194     (*fTrackPar)(1,0) = hit2->GetNonBendingCoor(); // x
 195     fPosition = hit2->GetZ(); // z
 196     fTrackHits->Add(hit1); // add hit 1
 197     fTrackHits->Add(hit2); // add hit 2
 198     //AZ(Z->-Z) fTrackDir = 1;
 199     fTrackDir = -1;
 200   }
 201   dZ = hit2->GetZ() - hit1->GetZ();
 202   dY = hit2->GetBendingCoor() - hit1->GetBendingCoor();
 203   dX = hit2->GetNonBendingCoor() - hit1->GetNonBendingCoor();
 204   (*fTrackPar)(2,0) = TMath::ATan2(dY,dZ); // alpha
 205   if ((*fTrackPar)(2,0) < 0.) (*fTrackPar)(2,0) += 2*TMath::Pi(); // from 0 to 2*pi
 206   (*fTrackPar)(3,0) = TMath::ATan2(-dX,dZ/TMath::Cos((*fTrackPar)(2,0))); // beta
 207   (*fTrackPar)(2,0) -= TMath::Pi();
 208   (*fTrackPar)(4,0) = 1/fgTrackReconstructor->GetBendingMomentumFromImpactParam(segment->GetBendingImpact()); // 1/Pt
 209   (*fTrackPar)(4,0) *= TMath::Cos((*fTrackPar)(3,0)); // 1/p
 210   // Evaluate covariance (and weight) matrix
 211   EvalCovariance(dZ);
 212
 213   if (fgDebug < 0 ) return;
 214   cout << fgTrackReconstructor->GetNRecTracks()-1 << " " << fgTrackReconstructor->GetBendingMomentumFromImpactParam(segment->GetBendingImpact()) << " " << 1/(*fTrackPar)(4,0) << " ";
 215   if (fgTrackReconstructor->GetRecTrackRefHits()) {
 216     // from track ref. hits
 217     cout << ((AliMUONHitForRec*)((*fTrackHits)[0]))->GetTTRTrack() << "<-->" << ((AliMUONHitForRec*)((*fTrackHits)[1]))->GetTTRTrack() << " @ " << fStartSegment->GetHitForRec1()->GetChamberNumber() << endl;
 218   } else {
 219     // from raw clusters
 220     for (Int_t i=0; i<2; i++) {
 221       hit1 = (AliMUONHitForRec*) ((*fTrackHits)[i]);
 222       rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit1->GetChamberNumber());
 223       clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit1->GetHitNumber());
 224       cout << clus->GetTrack(1);
 225       if (clus->GetTrack(2) != -1) cout << " " << clus->GetTrack(2);
 226       if (i == 0) cout << " <--> ";
 227     }
 228     cout << " @ " << fStartSegment->GetHitForRec1()->GetChamberNumber() << endl;
 229   }
 230 }
 231
 232   //__________________________________________________________________________
 233 AliMUONTrackK::~AliMUONTrackK()
 234 {
 235 /// Destructor
 236
 237   if (fTrackHits) {
 238     //cout << fNmbTrackHits << endl;
 239     for (Int_t i = 0; i < fNmbTrackHits; i++) {
 240       AliMUONHitForRec *hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(i);
 241       hit->SetNTrackHits(hit->GetNTrackHits()-1);
 242     }
 243     delete fTrackHits; // delete the TObjArray of pointers to TrackHit's
 244     fTrackHits = NULL;
 245   }
 246   if (fTrackPar) {
 247     delete fTrackPar; delete fTrackParNew; delete fCovariance;
 248     delete fWeight;
 249   }
 250
 251   if (fSteps) delete fSteps;
 252   if (fChi2Array) delete fChi2Array;
 253   if (fChi2Smooth) delete fChi2Smooth;
 254   if (fParSmooth) {fParSmooth->Delete(); delete fParSmooth; }
 255   // Delete only matrices not shared by several tracks
 256   if (!fParExtrap) return;
 257
 258   Int_t id = 0;
 259   for (Int_t i=0; i<fNSteps; i++) {
 260     //if (fParExtrap->UncheckedAt(i)->GetUniqueID() > 1)
 261     //  fParExtrap->UncheckedAt(i)->SetUniqueID(fParExtrap->RemoveAt(i)->GetUniqueID()-1);
 262     id = fParExtrap->UncheckedAt(i)->GetUniqueID();
 263     if (id > 1) {
 264       fParExtrap->UncheckedAt(i)->SetUniqueID(id-1);
 265       fParExtrap->RemoveAt(i);
 266     }
 267     //if (fParFilter->UncheckedAt(i)->GetUniqueID() > 1)
 268     //  fParFilter->UncheckedAt(i)->SetUniqueID(fParFilter->RemoveAt(i)->GetUniqueID()-1);
 269     id = fParFilter->UncheckedAt(i)->GetUniqueID();
 270     if (id > 1) {
 271       fParFilter->UncheckedAt(i)->SetUniqueID(id-1);
 272       fParFilter->RemoveAt(i);
 273     }
 274     //if (fCovExtrap->UncheckedAt(i)->GetUniqueID() > 1)
 275     //  fCovExtrap->UncheckedAt(i)->SetUniqueID(fCovExtrap->RemoveAt(i)->GetUniqueID()-1);
 276     id = fCovExtrap->UncheckedAt(i)->GetUniqueID();
 277     if (id > 1) {
 278       fCovExtrap->UncheckedAt(i)->SetUniqueID(id-1);
 279       fCovExtrap->RemoveAt(i);
 280     }
 281
 282     //if (fCovFilter->UncheckedAt(i)->GetUniqueID() > 1)
 283     //  fCovFilter->UncheckedAt(i)->SetUniqueID(fCovFilter->RemoveAt(i)->GetUniqueID()-1);
 284     id = fCovFilter->UncheckedAt(i)->GetUniqueID();
 285     if (id > 1) {
 286       fCovFilter->UncheckedAt(i)->SetUniqueID(id-1);
 287       fCovFilter->RemoveAt(i);
 288     }
 289     //if (fJacob->UncheckedAt(i)->GetUniqueID() > 1)
 290     //  fJacob->UncheckedAt(i)->SetUniqueID(fJacob->RemoveAt(i)->GetUniqueID()-1);
 291     id = fJacob->UncheckedAt(i)->GetUniqueID();
 292     if (id > 1) {
 293       fJacob->UncheckedAt(i)->SetUniqueID(id-1);
 294       fJacob->RemoveAt(i);
 295     }
 296   }
 297   /*
 298   for (Int_t i=0; i<fNSteps; i++) {
 299     if (fParExtrap->UncheckedAt(i)) ((TMatrixD*)fParExtrap->UncheckedAt(i))->Delete();
 300     if (fParFilter->UncheckedAt(i)) ((TMatrixD*)fParFilter->UncheckedAt(i))->Delete();
 301     if (fCovExtrap->UncheckedAt(i)) ((TMatrixD*)fCovExtrap->UncheckedAt(i))->Delete();
 302     cout << fCovFilter->UncheckedAt(i) << " " << (*fSteps)[i] << endl;
 303     if (fCovFilter->UncheckedAt(i)) ((TMatrixD*)fCovFilter->UncheckedAt(i))->Delete();
 304     if (fJacob->UncheckedAt(i)) ((TMatrixD*)fJacob->UncheckedAt(i))->Delete();
 305   }
 306   */
 307   fParExtrap->Delete(); fParFilter->Delete();
 308   fCovExtrap->Delete(); fCovFilter->Delete();
 309   fJacob->Delete();
 310   delete fParExtrap; delete fParFilter;
 311   delete fCovExtrap; delete fCovFilter;
 312   delete fJacob;
 313 }
 314
 315   //__________________________________________________________________________
 316 AliMUONTrackK & AliMUONTrackK::operator=(const AliMUONTrackK& source)
 317 {
 318 /// Assignment operator
 319
 320   // Members
 321   if(&source == this) return *this;
 322
 323   // base class assignement
 324   //AZ TObject::operator=(source);
 325   AliMUONTrack::operator=(source);
 326
 327   fStartSegment = source.fStartSegment;
 328   fNmbTrackHits = source.fNmbTrackHits;
 329   fChi2 = source.fChi2;
 330   fPosition = source.fPosition;
 331   fPositionNew = source.fPositionNew;
 332   fTrackDir = source.fTrackDir;
 333   fBPFlag = source.fBPFlag;
 334   fRecover = source.fRecover;
 335   fSkipHit = source.fSkipHit;
 336
 337   // Pointers
 338   fTrackHits = new TObjArray(*source.fTrackHits);
 339   //source.fTrackHits->Dump();
 340   //fTrackHits->Dump();
 341   for (Int_t i = 0; i < fNmbTrackHits; i++) {
 342     AliMUONHitForRec *hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(i);
 343     hit->SetNTrackHits(hit->GetNTrackHits()+1);
 344   }
 345
 346   fTrackPar = new TMatrixD(*source.fTrackPar); // track parameters
 347   fTrackParNew = new TMatrixD(*source.fTrackParNew); // track parameters
 348   fCovariance = new TMatrixD(*source.fCovariance); // covariance matrix
 349   fWeight = new TMatrixD(*source.fWeight); // weight matrix (inverse of covariance)
 350
 351   // For smoother
 352   fParExtrap = new TObjArray(*source.fParExtrap);
 353   fParFilter = new TObjArray(*source.fParFilter);
 354   fCovExtrap = new TObjArray(*source.fCovExtrap);
 355   fCovFilter = new TObjArray(*source.fCovFilter);
 356   fJacob = new TObjArray(*source.fJacob);
 357   fSteps = new TArrayD(*source.fSteps);
 358   fNSteps = source.fNSteps;
 359   fChi2Array = NULL;
 360   if (source.fChi2Array) fChi2Array = new TArrayD(*source.fChi2Array);
 361   fChi2Smooth = NULL;
 362   fParSmooth = NULL;
 363
 364   for (Int_t i=0; i<fParExtrap->GetEntriesFast(); i++) {
 365     fParExtrap->UncheckedAt(i)->SetUniqueID(fParExtrap->UncheckedAt(i)->GetUniqueID()+1);
 366     fParFilter->UncheckedAt(i)->SetUniqueID(fParFilter->UncheckedAt(i)->GetUniqueID()+1);
 367     fCovExtrap->UncheckedAt(i)->SetUniqueID(fCovExtrap->UncheckedAt(i)->GetUniqueID()+1);
 368     fCovFilter->UncheckedAt(i)->SetUniqueID(fCovFilter->UncheckedAt(i)->GetUniqueID()+1);
 369     fJacob->UncheckedAt(i)->SetUniqueID(fJacob->UncheckedAt(i)->GetUniqueID()+1);
 370   }
 371   return *this;
 372 }
 373
 374   //__________________________________________________________________________
 375 void AliMUONTrackK::EvalCovariance(Double_t dZ)
 376 {
 377 /// Evaluate covariance (and weight) matrix for track candidate
 378   Double_t sigmaB, sigmaNonB, tanA, tanB, dAdY, rad, dBdX, dBdY;
 379
 380   dZ = -dZ;
 381   sigmaB = fgTrackReconstructor->GetBendingResolution(); // bending resolution
 382   sigmaNonB = fgTrackReconstructor->GetNonBendingResolution(); // non-bending resolution
 383
 384   (*fWeight)(0,0) = sigmaB*sigmaB; // <yy>
 385
 386   (*fWeight)(1,1) = sigmaNonB*sigmaNonB; // <xx>
 387
 388   tanA = TMath::Tan((*fTrackPar)(2,0));
 389   dAdY = 1/(1+tanA*tanA)/dZ;
 390   (*fWeight)(2,2) = dAdY*dAdY*(*fWeight)(0,0)*2; // <aa>
 391   (*fWeight)(0,2) = dAdY*(*fWeight)(0,0); // <ya>
 392   (*fWeight)(2,0) = (*fWeight)(0,2);
 393
 394   rad = dZ/TMath::Cos((*fTrackPar)(2,0));
 395   tanB = TMath::Tan((*fTrackPar)(3,0));
 396   dBdX = 1/(1+tanB*tanB)/rad;
 397   dBdY = 0; // neglect
 398   (*fWeight)(3,3) = dBdX*dBdX*(*fWeight)(1,1)*2; // <bb>
 399   (*fWeight)(1,3) = dBdX*(*fWeight)(1,1); // <xb>
 400   (*fWeight)(3,1) = (*fWeight)(1,3);
 401
 402   (*fWeight)(4,4) = ((*fTrackPar)(4,0)*0.5)*((*fTrackPar)(4,0)*0.5); // error 50%
 403
 404   // check whether the Invert method returns flag if matrix cannot be inverted,
 405   // and do not calculate the Determinant in that case !!!!
 406   if (fWeight->Determinant() != 0) {
 407     // fWeight->Invert();
 408     Int_t ifail;
 409     mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
 410   } else {
 411     AliWarning(" Determinant fWeight=0:");
 412   }
 413   return;
 414 }
 415
 416   //__________________________________________________________________________
 417 Bool_t AliMUONTrackK::KalmanFilter(Int_t ichamBeg, Int_t ichamEnd, Bool_t Back, Double_t zDipole1, Double_t zDipole2)
 418 {
 419 /// Follows track through detector stations
 420   Bool_t miss, success;
 421   Int_t ichamb, iFB, iMin, iMax, dChamb, ichambOK;
 422   Int_t ihit, firstIndx = -1, lastIndx = -1, currIndx = -1, iz0 = -1, iz = -1;
 423   Double_t zEnd, dChi2;
 424   AliMUONHitForRec *hitAdd, *firstHit = NULL, *lastHit = NULL, *hit;
 425   AliMUONRawCluster *clus;
 426   TClonesArray *rawclusters;
 427   hit = 0; clus = 0; rawclusters = 0;
 428
 429   miss = success = kTRUE;
 430   Int_t endOfProp = 0;
 431   //iFB = (ichamEnd == ichamBeg) ? fTrackDir : TMath::Sign(1,ichamEnd-ichamBeg);
 432   iFB = (ichamEnd == ichamBeg) ? -fTrackDir : TMath::Sign(1,ichamEnd-ichamBeg);
 433   iMin = TMath::Min(ichamEnd,ichamBeg);
 434   iMax = TMath::Max(ichamEnd,ichamBeg);
 435   ichamb = ichamBeg;
 436   ichambOK = ichamb;
 437
 438   if (Back) {
 439     // backpropagation
 440     currIndx = 1;
 441     if (((AliMUONHitForRec*)fTrackHits->First())->GetChamberNumber() != ichamb) currIndx = 0;
 442   } else if (fRecover) {
 443     hit = GetHitLastOk();
 444     currIndx = fTrackHits->IndexOf(hit);
 445     if (currIndx < 0) hit = fStartSegment->GetHitForRec1(); // for station 3
 446     Back = kTRUE;
 447     ichamb = hit->GetChamberNumber();
 448     if (hit == fSkipHit || fRecover == 2 && currIndx >= 0) {
 449       // start from the last point or outlier
 450       // remove the skipped hit
 451       fTrackHits->Remove(fSkipHit); // remove hit
 452       fNmbTrackHits --;
 453       fSkipHit->SetNTrackHits(fSkipHit->GetNTrackHits()-1); // unmark hit
 454       if (fRecover == 1) {
 455         // recovery
 456         Back = kFALSE;
 457         fRecover = 0;
 458         ichambOK = ((AliMUONHitForRec*)((*fTrackHits)[fNmbTrackHits-1]))->GetChamberNumber();
 459         //if (ichambOK >= 8 && ((AliMUONHitForRec*)((*fTrackHits)[0]))->GetChamberNumber() == 6) ichambOK = 6;
 460         if (fgTrackReconstructor->GetTrackMethod() == 3 &&
 461             fSkipHit->GetHitNumber() < 0) {
 462           iz0 = fgCombi->IZfromHit(fSkipHit);
 463           currIndx = -1;
 464         }
 465         else currIndx = fgHitForRec->IndexOf(fSkipHit);
 466       } else {
 467         // outlier
 468         fTrackHits->Compress();
 469       }
 470     } // if (hit == fSkipHit)
 471     else if (currIndx < 0) currIndx = fTrackHits->IndexOf(hit);
 472   } // else if (fRecover)
 473   else {
 474     // Get indices of the 1'st and last hits on the track candidate
 475     firstHit = (AliMUONHitForRec*) fTrackHits->First();
 476     lastHit = (AliMUONHitForRec*) fTrackHits->Last();
 477     if (fgTrackReconstructor->GetTrackMethod() == 3 &&
 478         lastHit->GetHitNumber() < 0) iz0 = fgCombi->IZfromHit(lastHit);
 479     else {
 480       firstIndx = fgHitForRec->IndexOf(firstHit);
 481       lastIndx = fgHitForRec->IndexOf(lastHit);
 482       currIndx = TMath::Abs (TMath::Max(firstIndx*iFB,lastIndx*iFB));
 483     }
 484   }
 485
 486   if (iz0 < 0) iz0 = iFB;
 487   while (ichamb >= iMin && ichamb <= iMax) {
 488   // Find the closest hit in Z, not belonging to the current plane
 489     if (Back) {
 490       // backpropagation
 491       if (currIndx < fNmbTrackHits) {
 492         hitAdd = (AliMUONHitForRec*) fTrackHits->UncheckedAt(currIndx);
 493         zEnd = hitAdd->GetZ();
 494         //AZ(z->-z) } else zEnd = -9999;
 495       } else zEnd = 9999;
 496     } else {
 497       //AZ(Z->-Z) zEnd = -9999;
 498       zEnd = 9999;
 499       for (ihit=currIndx+iFB; ihit>=0 && ihit<fgNOfPoints; ihit+=iFB) {
 500         hitAdd = (AliMUONHitForRec*) ((*fgHitForRec)[ihit]);
 501         //if (TMath::Abs(hitAdd->GetZ()-fPosition) > 0.1) {
 502         if (TMath::Abs(hitAdd->GetZ()-fPosition) > 0.5) {
 503           zEnd = hitAdd->GetZ();
 504           currIndx = ihit;
 505           break;
 506         }
 507       }
 508
 509       // Combined cluster / track finder
 510       if (zEnd > 999 && iFB < 0 && fgTrackReconstructor->GetTrackMethod() == 3) {
 511         currIndx = -2;
 512         AliMUONHitForRec hitTmp;
 513         for (iz = iz0 - iFB; iz < fgCombi->Nz(); iz++) {
 514           if (TMath::Abs(fgCombi->Z(iz)-fPosition) < 0.5) continue;
 515           Int_t *pDEatZ = fgCombi->DEatZ(iz);
 516           //cout << iz << " " << fgCombi->Z(iz) << endl;
 517           zEnd = fgCombi->Z(iz);
 518           iz0 = iz;
 519           AliMUONDetElement *detElem = fgCombi->DetElem(pDEatZ[0]);
 520           hitAdd = &hitTmp;
 521           hitAdd->SetChamberNumber(detElem->Chamber());
 522           //hitAdd = (AliMUONHitForRec*) detElem->HitsForRec()->First();
 523           if (hitAdd) break;
 524         }
 525       }
 526     }
 527     if (zEnd>999 && ichamb==ichamEnd) endOfProp = 1; // end-of-propagation
 528     else {
 529       // Check if there is a chamber without hits
 530       if (zEnd>999 || TMath::Abs(hitAdd->GetChamberNumber()-ichamb) > 1) {
 531         if (!Back && zEnd<999) currIndx -= iFB;
 532         ichamb += iFB;
 533         zEnd = AliMUONConstants::DefaultChamberZ(ichamb);
 534         miss = kTRUE;
 535       } else {
 536         ichamb = hitAdd->GetChamberNumber();
 537         miss = kFALSE;
 538       }
 539     }
 540     if (ichamb<iMin || ichamb>iMax) break;
 541     // Check for missing station
 542     if (!Back) {
 543       dChamb = TMath::Abs(ichamb-ichambOK);
 544       //cout << dChamb << " " << ichambOK << " " << fgNOfPoints << endl;
 545       Int_t dStatMiss = TMath::Abs (ichamb/2 - ichambOK/2);
 546       if (zEnd > 999) dStatMiss++;
 547       if (dStatMiss > 1) {
 548       //if (dStatMiss == 2 && ichambOK/2 != 3 || dStatMiss > 2) { // AZ - missing st. 3
 549         // missing station - abandon track
 550         //cout << dChamb << " " << ichambOK << " " << fgNOfPoints << " " << 1/(*fTrackPar)(4,0) << endl;
 551         if (fgDebug >= 10) {
 552           for (Int_t i1=0; i1<fgNOfPoints; i1++) {
 553             cout << " Hit #" << ((AliMUONHitForRec*)((*fgHitForRec)[i1]))->GetChamberNumber() << " ";
 554             cout << ((AliMUONHitForRec*)((*fgHitForRec)[i1]))->GetBendingCoor() << " ";
 555             cout << ((AliMUONHitForRec*)((*fgHitForRec)[i1]))->GetNonBendingCoor() << " ";
 556             cout << ((AliMUONHitForRec*)((*fgHitForRec)[i1]))->GetZ() << " " << " ";
 557             cout << ((AliMUONHitForRec*)((*fgHitForRec)[i1]))->GetTTRTrack() << endl;
 558           }
 559           cout << endl;
 560           cout << fNmbTrackHits << endl;
 561           for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
 562             hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
 563             printf(" * %d %10.4f %10.4f %10.4f",
 564                    hit->GetChamberNumber(), hit->GetBendingCoor(),
 565                    hit->GetNonBendingCoor(), hit->GetZ());
 566             if (fgTrackReconstructor->GetRecTrackRefHits()) {
 567               // from track ref. hits
 568               printf(" %3d %3d \n", hit->GetTrackRefSignal(), hit->GetTTRTrack());
 569             } else {
 570               // from raw clusters
 571               rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
 572               clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
 573               printf("%3d", clus->GetTrack(1)-1);
 574               if (clus->GetTrack(2) != 0) printf("%3d \n", clus->GetTrack(2)-1);
 575               else printf("\n");
 576             }
 577           }
 578         } // if (fgDebug >= 10)
 579         if (fNmbTrackHits>2 && fRecover==0) Recover(); // try to recover track later
 580         return kFALSE;
 581       } // if (dStatMiss > 1)
 582     } // if (!Back)
 583     if (endOfProp != 0) break;
 584
 585     // propagate to the found Z
 586
 587     // Check if track steps into dipole
 588     //AZ(z->-z) if (fPosition>zDipole2 && zEnd<zDipole2) {
 589     if (fPosition<zDipole2 && zEnd>zDipole2) {
 590       //LinearPropagation(zDipole2-zBeg);
 591       ParPropagation(zDipole2);
 592       MSThin(1); // multiple scattering in the chamber
 593       WeightPropagation(zDipole2, kTRUE); // propagate weight matrix
 594       fPosition = fPositionNew;
 595       *fTrackPar = *fTrackParNew;
 596       //MagnetPropagation(zEnd);
 597       ParPropagation(zEnd);
 598       WeightPropagation(zEnd, kTRUE);
 599       fPosition = fPositionNew;
 600     }
 601     // Check if track steps out of dipole
 602     //AZ(z->-z) else if (fPosition>zDipole1 && zEnd<zDipole1) {
 603     else if (fPosition<zDipole1 && zEnd>zDipole1) {
 604       //MagnetPropagation(zDipole1-zBeg);
 605       ParPropagation(zDipole1);
 606       MSThin(1); // multiple scattering in the chamber
 607       WeightPropagation(zDipole1, kTRUE);
 608       fPosition = fPositionNew;
 609       *fTrackPar = *fTrackParNew;
 610       //LinearPropagation(zEnd-zDipole1);
 611       ParPropagation(zEnd);
 612       WeightPropagation(zEnd, kTRUE);
 613       fPosition = fPositionNew;
 614     } else {
 615       ParPropagation(zEnd);
 616       //MSThin(1); // multiple scattering in the chamber
 617       if (TMath::Abs(zEnd-fPosition) > 5) MSThin(1); // multiple scattering in the chamber
 618       WeightPropagation(zEnd, kTRUE);
 619       fPosition = fPositionNew;
 620     }
 621
 622     // Add measurement
 623     if (fRecover != 0 && hitAdd == fSkipHit && !miss) {
 624       // recovered track - remove the hit
 625       miss = kTRUE;
 626       ichamb = fSkipHit->GetChamberNumber();
 627       // remove the skipped hit
 628       fTrackHits->Remove(fSkipHit);
 629       fNmbTrackHits --;
 630       //AZ fSkipHit->SetNTrackHits(fSkipHit->GetNTrackHits()-1); // unmark hit
 631       Back = kFALSE;
 632       fRecover = 0;
 633       ichambOK = ((AliMUONHitForRec*)((*fTrackHits)[fNmbTrackHits-1]))->GetChamberNumber();
 634       //? if (fgTrackReconstructor->GetTrackMethod() != 3) currIndx = fgHitForRec->IndexOf(fSkipHit);
 635       currIndx = fgHitForRec->IndexOf(fSkipHit);
 636     }
 637
 638     if (Back && !miss) {
 639       // backward propagator
 640       if (currIndx) {
 641         TMatrixD pointWeight(fgkSize,fgkSize);
 642         TMatrixD point(fgkSize,1);
 643         TMatrixD trackParTmp = point;
 644         point(0,0) = hitAdd->GetBendingCoor();
 645         point(1,0) = hitAdd->GetNonBendingCoor();
 646         pointWeight(0,0) = 1/hitAdd->GetBendingReso2();
 647         pointWeight(1,1) = 1/hitAdd->GetNonBendingReso2();
 648         TryPoint(point,pointWeight,trackParTmp,dChi2);
 649         *fTrackPar = trackParTmp;
 650         *fWeight += pointWeight;
 651         if (fTrackDir > 0) AddMatrices (this, dChi2, hitAdd);
 652         fChi2 += dChi2; // Chi2
 653       } else *fTrackPar = *fTrackParNew; // adjust end point to the last hit
 654       if (ichamb==ichamEnd) break;
 655       currIndx ++;
 656     } else {
 657       // forward propagator
 658       if (miss || !FindPoint(ichamb, zEnd, currIndx, iFB, hitAdd, iz)) {
 659         // missing point
 660         *fTrackPar = *fTrackParNew;
 661       } else {
 662         //add point
 663         fTrackHits->Add(hitAdd); // add hit
 664         fNmbTrackHits ++;
 665         hitAdd->SetNTrackHits(hitAdd->GetNTrackHits()+1); // mark hit as being on track
 666         ichambOK = ichamb;
 667         currIndx = fgHitForRec->IndexOf(hitAdd); // Check
 668       }
 669     }
 670   } // while
 671   if (fgDebug > 0) cout << fNmbTrackHits << " " << fChi2 << " " << 1/(*fTrackPar)(4,0) << " " << fPosition << endl;
 672   if (1/TMath::Abs((*fTrackPar)(4,0)) < fgTrackReconstructor->GetMinBendingMomentum()) success = kFALSE; // p < p_min
 673   return success;
 674 }
 675
 676   //__________________________________________________________________________
 677 void AliMUONTrackK::ParPropagation(Double_t zEnd)
 678 {
 679 /// Propagation of track parameters to zEnd
 680   Int_t iFB, nTries;
 681   Double_t dZ, step, distance, charge;
 682   Double_t vGeant3[7], vGeant3New[7];
 683   AliMUONTrackParam trackParam;
 684
 685   nTries = 0;
 686   // First step using linear extrapolation
 687   dZ = zEnd - fPosition;
 688   fPositionNew = fPosition;
 689   *fTrackParNew = *fTrackPar;
 690   if (dZ == 0) return; //AZ ???
 691   iFB = (Int_t)TMath::Sign(Double_t(1.0),dZ);
 692   step = dZ/TMath::Cos((*fTrackPar)(2,0))/TMath::Cos((*fTrackPar)(3,0)); // linear estimate
 693   //AZ(z->-z) charge = iFB*TMath::Sign(Double_t(1.0),(*fTrackPar)(4,0));
 694   charge = -iFB*TMath::Sign(Double_t(1.0),(*fTrackPar)(4,0));
 695   SetGeantParam(vGeant3,iFB);
 696   //fTrackParNew->Print();
 697
 698   // Check if overstep
 699   do {
 700     step = TMath::Abs(step);
 701     // Propagate parameters
 702     trackParam.ExtrapOneStepRungekutta(charge,step,vGeant3,vGeant3New);
 703     //extrap_onestep_rungekutta(charge,step,vGeant3,vGeant3New);
 704     distance = zEnd - vGeant3New[2];
 705     step *= dZ/(vGeant3New[2]-fPositionNew);
 706     nTries ++;
 707   } while (distance*iFB < 0 && TMath::Abs(distance) > fgkEpsilon);
 708
 709   GetFromGeantParam(vGeant3New,iFB);
 710   //fTrackParNew->Print();
 711
 712   // Position adjustment (until within tolerance)
 713   while (TMath::Abs(distance) > fgkEpsilon) {
 714     dZ = zEnd - fPositionNew;
 715     iFB = (Int_t)TMath::Sign(Double_t(1.0),dZ);
 716     step = dZ/TMath::Cos((*fTrackParNew)(2,0))/TMath::Cos((*fTrackParNew)(3,0));
 717     step = TMath::Abs(step);
 718     SetGeantParam(vGeant3,iFB);
 719     do {
 720       // binary search
 721       // Propagate parameters
 722       trackParam.ExtrapOneStepRungekutta(charge,step,vGeant3,vGeant3New);
 723       //extrap_onestep_rungekutta(charge,step,vGeant3,vGeant3New);
 724       distance = zEnd - vGeant3New[2];
 725       step /= 2;
 726       nTries ++;
 727       if (nTries > fgkTriesMax) AliError(Form(" Too many tries: %d", nTries));
 728     } while (distance*iFB < 0);
 729
 730     GetFromGeantParam(vGeant3New,iFB);
 731   }
 732   //cout << nTries << endl;
 733   //fTrackParNew->Print();
 734   return;
 735 }
 736
 737   //__________________________________________________________________________
 738 void AliMUONTrackK::WeightPropagation(Double_t zEnd, Bool_t smooth)
 739 {
 740 /// Propagation of the weight matrix
 741 /// W = DtWD, where D is Jacobian
 742   Int_t i, j;
 743   Double_t dPar;
 744
 745   TMatrixD jacob(fgkSize,fgkSize);
 746   jacob = 0;
 747
 748   // Save initial and propagated parameters
 749   TMatrixD trackPar0 = *fTrackPar;
 750   TMatrixD trackParNew0 = *fTrackParNew;
 751
 752   // Get covariance matrix
 753   *fCovariance = *fWeight;
 754   // check whether the Invert method returns flag if matrix cannot be inverted,
 755   // and do not calculate the Determinant in that case !!!!
 756   if (fCovariance->Determinant() != 0) {
 757     Int_t ifail;
 758     mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
 759     //fCovariance->Print();
 760   } else {
 761     AliWarning(" Determinant fCovariance=0:");
 762   }
 763
 764   // Loop over parameters to find change of the propagated vs initial ones
 765   for (i=0; i<fgkSize; i++) {
 766     dPar = TMath::Sqrt((*fCovariance)(i,i));
 767     *fTrackPar = trackPar0;
 768     if (i == 4) dPar *= TMath::Sign(1.,-trackPar0(4,0)); // 1/p
 769     (*fTrackPar)(i,0) += dPar;
 770     ParPropagation(zEnd);
 771     for (j=0; j<fgkSize; j++) {
 772       jacob(j,i) = ((*fTrackParNew)(j,0)-trackParNew0(j,0))/dPar;
 773     }
 774   }
 775
 776   //trackParNew0.Print();
 777   //TMatrixD par1(jacob,TMatrixD::kMult,trackPar0); //
 778   //par1.Print();
 779   TMatrixD jacob0 = jacob;
 780   if (jacob.Determinant() != 0) {
 781     jacob.Invert();
 782   } else {
 783     AliWarning(" Determinant jacob=0:");
 784   }
 785   TMatrixD weight1(*fWeight,TMatrixD::kMult,jacob); // WD
 786   *fWeight = TMatrixD(jacob,TMatrixD::kTransposeMult,weight1); // DtWD
 787
 788   // Restore initial and propagated parameters
 789   *fTrackPar = trackPar0;
 790   *fTrackParNew = trackParNew0;
 791
 792   // Save for smoother
 793   if (!smooth) return; // do not use smoother
 794   if (fTrackDir < 0) return; // only for propagation towards int. point
 795   TMatrixD *tmp = new TMatrixD(*fTrackParNew); // extrapolated parameters
 796   fParExtrap->Add(tmp);
 797   tmp->SetUniqueID(1);
 798   tmp = new TMatrixD(*fTrackParNew); // filtered parameters (if no measurement will be found)
 799   fParFilter->Add(tmp);
 800   tmp->SetUniqueID(1);
 801   *fCovariance = *fWeight;
 802   if (fCovariance->Determinant() != 0) {
 803     Int_t ifail;
 804     mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
 805   } else {
 806     AliWarning(" Determinant fCovariance=0:");
 807   }
 808   tmp = new TMatrixD(*fCovariance); // extrapolated covariance
 809   fCovExtrap->Add(tmp);
 810   tmp->SetUniqueID(1);
 811   tmp = new TMatrixD(*fCovariance); // filtered covariance (if no measurement will be found)
 812   fCovFilter->Add(tmp);
 813   tmp->SetUniqueID(1);
 814   tmp = new TMatrixD(jacob0); // Jacobian
 815   fJacob->Add(tmp);
 816   tmp->SetUniqueID(1);
 817   if (fSteps->fN <= fNSteps) fSteps->Set(fSteps->fN+10);
 818   fSteps->AddAt(fPositionNew,fNSteps++);
 819   if (fgDebug > 0) cout << " WeightPropagation " << fNSteps << " " << fPositionNew << endl;
 820   return;
 821 }
 822
 823   //__________________________________________________________________________
 824 Bool_t AliMUONTrackK::FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int_t iFB, AliMUONHitForRec *&hitAdd, Int_t iz)
 825 {
 826 /// Picks up point within a window for the chamber No ichamb
 827 /// Split the track if there are more than 1 hit
 828   Int_t ihit, nRecTracks;
 829   Double_t windowB, windowNonB, dChi2Tmp=0, dChi2, y, x, savePosition=0;
 830   TClonesArray *trackPtr;
 831   AliMUONHitForRec *hit, *hitLoop;
 832   AliMUONTrackK *trackK;
 833   AliMUONDetElement *detElem = NULL;
 834
 835   Bool_t ok = kFALSE;
 836
 837   if (fgTrackReconstructor->GetTrackMethod() == 3 && iz >= 0) {
 838     // Combined cluster / track finder
 839     // Check if extrapolated track passes thru det. elems. at iz
 840     Int_t *pDEatZ = fgCombi->DEatZ(iz);
 841     Int_t nDetElem = pDEatZ[-1];
 842     //cout << fgCombi->Z(iz) << " " << nDetElem << endl;
 843     for (Int_t i = 0; i < nDetElem; i++) {
 844       detElem = fgCombi->DetElem(pDEatZ[i]);
 845       if (detElem->Inside((*fTrackParNew)(1,0), (*fTrackParNew)(0,0), fPosition)) {
 846         detElem->ClusterReco((*fTrackParNew)(1,0), (*fTrackParNew)(0,0));
 847         hitAdd = (AliMUONHitForRec*) detElem->HitsForRec()->First();
 848         ok = kTRUE;
 849         break;
 850       }
 851     }
 852     if (!ok) return ok; // outside det. elem.
 853     ok = kFALSE;
 854   }
 855
 856   //sigmaB = fgTrackReconstructor->GetBendingResolution(); // bending resolution
 857   //sigmaNonB = fgTrackReconstructor->GetNonBendingResolution(); // non-bending resolution
 858   *fCovariance = *fWeight;
 859   // check whether the Invert method returns flag if matrix cannot be inverted,
 860   // and do not calculate the Determinant in that case !!!!
 861   if (fCovariance->Determinant() != 0) {
 862       Int_t ifail;
 863       mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
 864   } else {
 865     AliWarning(" Determinant fCovariance=0:");
 866   }
 867   //windowB = fgkNSigma*TMath::Sqrt((*fCovariance)(0,0)+sigmaB*sigmaB);
 868   //windowNonB = fgkNSigma*TMath::Sqrt((*fCovariance)(1,1)+sigmaNonB*sigmaNonB);
 869   // Loop over all hits and take hits from the chamber
 870   TMatrixD pointWeight(fgkSize,fgkSize);
 871   TMatrixD saveWeight = pointWeight;
 872   TMatrixD pointWeightTmp = pointWeight;
 873   TMatrixD point(fgkSize,1);
 874   TMatrixD trackPar = point;
 875   TMatrixD trackParTmp = point;
 876   Int_t nHitsOK = 0, ihitB = currIndx, ihitE = fgNOfPoints, iDhit = iFB;
 877   Double_t zLast;
 878   zLast = ((AliMUONHitForRec*)fTrackHits->Last())->GetZ();
 879   if (fgTrackReconstructor->GetTrackMethod() == 3 && detElem) {
 880     ihitB = 0;
 881     ihitE = detElem->NHitsForRec();
 882     iDhit = 1;
 883   }
 884
 885
 886   for (ihit = ihitB; ihit >= 0 && ihit < ihitE; ihit+=iDhit) {
 887     if (fgTrackReconstructor->GetTrackMethod() != 3 || !detElem)
 888       hit = (AliMUONHitForRec*) ((*fgHitForRec)[ihit]);
 889     else hit = (AliMUONHitForRec*) detElem->HitsForRec()->UncheckedAt(ihit);
 890     if (hit->GetChamberNumber() == ichamb) {
 891       //if (TMath::Abs(hit->GetZ()-zEnd) < 0.1) {
 892       if (TMath::Abs(hit->GetZ()-zEnd) < 0.5) {
 893         y = hit->GetBendingCoor();
 894         x = hit->GetNonBendingCoor();
 895         if (hit->GetBendingReso2() < 0) {
 896           // Combined cluster / track finder
 897           hit->SetBendingReso2(fgTrackReconstructor->GetBendingResolution()*
 898                                fgTrackReconstructor->GetBendingResolution());
 899           hit->SetNonBendingReso2(fgTrackReconstructor->GetNonBendingResolution()*
 900                                   fgTrackReconstructor->GetNonBendingResolution());
 901         }
 902         windowB = fgkNSigma*TMath::Sqrt((*fCovariance)(0,0)+hit->GetBendingReso2());
 903         windowNonB = fgkNSigma*TMath::Sqrt((*fCovariance)(1,1)+hit->GetNonBendingReso2());
 904
 905         // windowB = TMath::Min (windowB,5.);
 906         if (fgkNSigma > 6) windowB = TMath::Min (windowB,5.);
 907         /*
 908         if (TMath::Abs(hit->GetZ()-zLast) < 50) {
 909           windowB = TMath::Min (windowB,0.5);
 910           windowNonB = TMath::Min (windowNonB,3.);
 911         } else if (TMath::Abs(hit->GetZ()-zLast) < 200) {
 912           windowB = TMath::Min (windowB,1.5);
 913           windowNonB = TMath::Min (windowNonB,3.);
 914         } else if (TMath::Abs(hit->GetZ()-zLast) < 350) {
 915           windowB = TMath::Min (windowB,4.);
 916           windowNonB = TMath::Min (windowNonB,6.);
 917         }
 918         */
 919
 920         // Test
 921         /*
 922         if (TMath::Abs(hit->GetZ()-zLast) < 50) {
 923           windowB = 5*TMath::Sqrt((*fCovariance)(0,0)+hit->GetBendingReso2());
 924           windowNonB = 5*TMath::Sqrt((*fCovariance)(1,1)+hit->GetNonBendingReso2());
 925         } else {
 926           windowB = fgkNSigma*TMath::Sqrt((*fCovariance)(0,0)+hit->GetBendingReso2());
 927           windowNonB = fgkNSigma*TMath::Sqrt((*fCovariance)(1,1)+hit->GetNonBendingReso2());
 928         }
 929         */
 930
 931         //cout << TMath::Abs((*fTrackParNew)(0,0)-y) << " " << windowB << " " << TMath::Abs((*fTrackParNew)(1,0)-x) << " " << windowNonB << " " << zEnd << endl;
 932         if (TMath::Abs((*fTrackParNew)(0,0)-y) <= windowB &&
 933             TMath::Abs((*fTrackParNew)(1,0)-x) <= windowNonB) {
 934         //if (TMath::Abs((*fTrackParNew)(0,0)-y) <= windowB &&
 935           //    TMath::Abs((*fTrackParNew)(1,0)-x) <= windowNonB &&
 936           //  hit->GetTrackRefSignal() == 1) { // just for test
 937           // Vector of measurements and covariance matrix
 938           //fprintf(lun1,"%3d %3d %10.4f %10.4f \n", gAlice->GetEvNumber(), ichamb, x, y);
 939           if (TMath::Abs(hit->GetZ()-zEnd) > 0.05) {
 940             // Adjust position: for multiple hits in the chamber or misalignment (Z as a function of X or Y)
 941             //AliWarning(Form(" *** adjust %f %f ", zEnd, hit->GetZ()));
 942             zEnd = hit->GetZ();
 943             *fTrackPar = *fTrackParNew;
 944             ParPropagation(zEnd);
 945             WeightPropagation(zEnd, kTRUE);
 946             fPosition = fPositionNew;
 947             *fTrackPar = *fTrackParNew;
 948             // Get covariance
 949             *fCovariance = *fWeight;
 950             if (fCovariance->Determinant() != 0) {
 951               Int_t ifail;
 952               mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
 953             } else {
 954               AliWarning(" Determinant fCovariance=0:" );
 955             }
 956           }
 957           point.Zero();
 958           point(0,0) = y;
 959           point(1,0) = x;
 960           pointWeight(0,0) = 1/hit->GetBendingReso2();
 961           pointWeight(1,1) = 1/hit->GetNonBendingReso2();
 962           TryPoint(point,pointWeight,trackPar,dChi2);
 963           if (TMath::Abs(1./(trackPar)(4,0)) < fgTrackReconstructor->GetMinBendingMomentum()) continue; // p < p_min - next hit
 964           // if (TMath::Sign(1.,(trackPar)(4,0)*(*fTrackPar)(4,0)) < 0) continue; // change of sign
 965           ok = kTRUE;
 966           nHitsOK++;
 967           //if (nHitsOK > -1) {
 968           if (nHitsOK == 1) {
 969             // Save current members
 970             saveWeight = *fWeight;
 971             savePosition = fPosition;
 972             // temporary storage for the current track
 973             dChi2Tmp = dChi2;
 974             trackParTmp = trackPar;
 975             pointWeightTmp = pointWeight;
 976             hitAdd = hit;
 977             if (fgDebug > 0) cout << " Added point: " << x << " " << y << " " << dChi2 << endl;
 978           } else {
 979             // branching: create a new track
 980             trackPtr = fgTrackReconstructor->GetRecTracksPtr();
 981             nRecTracks = fgTrackReconstructor->GetNRecTracks();
 982             trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(NULL, NULL);
 983             *trackK = *this;
 984             fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
 985             if (fgDebug > 0) cout << " ******** New track: " << ichamb << " " << hit->GetTTRTrack() << " " << 1/(trackPar)(4,0) << " " << hit->GetBendingCoor() << " " << hit->GetNonBendingCoor() << " " << fNmbTrackHits << " " << nRecTracks << endl;
 986             trackK->fRecover = 0;
 987             *(trackK->fTrackPar) = trackPar;
 988             *(trackK->fWeight) += pointWeight;
 989             trackK->fChi2 += dChi2;
 990             // Check
 991             /*
 992             *(trackK->fCovariance) = *(trackK->fWeight);
 993             if (trackK->fCovariance->Determinant() != 0) {
 994               Int_t ifail;
 995               mnvertLocal(&((*(trackK->fCovariance))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
 996             }
 997             cout << (*(trackK->fCovariance))(0,0) << " " << (*(trackK->fCovariance))(1,1) << " " << (*fCovariance)(0,0) << " " << (*fCovariance)(1,1) << endl;
 998             */
 999             // Add filtered matrices for the smoother
1000             if (fTrackDir > 0) {
1001               if (nHitsOK > 2) { // check for position adjustment
1002                 for (Int_t i=trackK->fNSteps-1; i>=0; i--) {
1003                   if (TMath::Abs(hit->GetZ()-(*trackK->fSteps)[i]) > 0.1) {
1004                   //if (TMath::Abs(hit->GetZ()-(trackK->fSteps)[i]) > 0.1) {
1005                     RemoveMatrices(trackK);
1006                     AliError(Form(" *** Position adjustment 1: %f %f", hit->GetZ(), (*trackK->fSteps)[i]));
1007                   }
1008                   else break;
1009                 }
1010               }
1011               AddMatrices (trackK, dChi2, hit);
1012             }
1013             // Mark hits as being on 2 tracks
1014             for (Int_t i=0; i<fNmbTrackHits; i++) {
1015               hitLoop = (AliMUONHitForRec*) ((*fTrackHits)[i]);
1016               //AZ hitLoop->SetNTrackHits(hitLoop->GetNTrackHits()+1);
1017               if (fgDebug >=10) {
1018                 cout << " ** ";
1019                 cout << hitLoop->GetChamberNumber() << " ";
1020                 cout << hitLoop->GetBendingCoor() << " ";
1021                 cout << hitLoop->GetNonBendingCoor() << " ";
1022                 cout << hitLoop->GetZ() << " " << " ";
1023                 cout << hitLoop->GetTrackRefSignal() << " " << " ";
1024                 cout << hitLoop->GetTTRTrack() << endl;
1025                 printf(" ** %d %10.4f %10.4f %10.4f %d %d \n",
1026                        hitLoop->GetChamberNumber(), hitLoop->GetBendingCoor(),
1027                        hitLoop->GetNonBendingCoor(), hitLoop->GetZ(),
1028                        hitLoop->GetTrackRefSignal(), hitLoop->GetTTRTrack());
1029               }
1030             }
1031             //add point
1032             trackK->fTrackHits->Add(hit); // add hit
1033             trackK->fNmbTrackHits ++;
1034             hit->SetNTrackHits(hit->GetNTrackHits()+1); // mark hit as being on track
1035             if (ichamb == 9) {
1036               // the last chamber
1037               trackK->fTrackDir = 1;
1038               trackK->fBPFlag = kTRUE;
1039             }
1040           }
1041         }
1042       }
1043     } else break; // different chamber
1044   } // for (ihit=currIndx;
1045   if (ok) {
1046     // Restore members
1047     *fTrackPar = trackParTmp;
1048     *fWeight = saveWeight;
1049     *fWeight += pointWeightTmp;
1050     fChi2 += dChi2Tmp; // Chi2
1051     fPosition = savePosition;
1052     // Add filtered matrices for the smoother
1053     if (fTrackDir > 0) {
1054       for (Int_t i=fNSteps-1; i>=0; i--) {
1055         if (TMath::Abs(fPosition-(*fSteps)[i]) > 0.1) {
1056         //if (TMath::Abs(fPosition-fSteps[i]) > 0.1) {
1057           RemoveMatrices(this);
1058           if (fgDebug > 0) cout << " *** Position adjustment 2 " << fPosition << " " << (*fSteps)[i] << endl;
1059         }
1060         else break;
1061       } // for (Int_t i=fNSteps-1;
1062       AddMatrices (this, dChi2Tmp, hitAdd);
1063       /*
1064       if (nHitsOK > 1) {
1065         for (Int_t i=0; i<fNSteps; i++) cout << (*fSteps)[i] << " "; cout << endl;
1066         for (Int_t i=0; i<trackK->fNSteps; i++) cout << (*trackK->fSteps)[i] << " "; cout << endl;
1067       }
1068       */
1069     } // if (fTrackDir > 0)
1070   }
1071   return ok;
1072 }
1073
1074   //__________________________________________________________________________
1075 void AliMUONTrackK::TryPoint(TMatrixD &point, const TMatrixD &pointWeight, TMatrixD &trackParTmp, Double_t &dChi2)
1076 {
1077 /// Adds a measurement point (modifies track parameters and computes
1078 /// change of Chi2)
1079
1080   // Solving linear system (W+U)p' = U(m-p) + (W+U)p
1081   TMatrixD wu = *fWeight;
1082   wu += pointWeight; // W+U
1083   trackParTmp = point;
1084   trackParTmp -= *fTrackParNew; // m-p
1085   TMatrixD right(pointWeight,TMatrixD::kMult,trackParTmp); // U(m-p)
1086   TMatrixD right1(wu,TMatrixD::kMult,*fTrackParNew); // (W+U)p
1087   right += right1; // U(m-p) + (W+U)p
1088
1089   // check whether the Invert method returns flag if matrix cannot be inverted,
1090   // and do not calculate the Determinant in that case !!!!
1091   if (wu.Determinant() != 0) {
1092     Int_t ifail;
1093     mnvertLocal(&((wu)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1094   } else {
1095     AliWarning(" Determinant wu=0:");
1096   }
1097   trackParTmp = TMatrixD(wu,TMatrixD::kMult,right);
1098
1099   right1 = trackParTmp;
1100   right1 -= point; // p'-m
1101   point = trackParTmp;
1102   point -= *fTrackParNew; // p'-p
1103   right = TMatrixD(*fWeight,TMatrixD::kMult,point); // W(p'-p)
1104   TMatrixD value(point,TMatrixD::kTransposeMult,right); // (p'-p)'W(p'-p)
1105   dChi2 = value(0,0);
1106   right = TMatrixD(pointWeight,TMatrixD::kMult,right1); // U(p'-m)
1107   value = TMatrixD(right1,TMatrixD::kTransposeMult,right); // (p'-m)'U(p'-m)
1108   dChi2 += value(0,0);
1109   return;
1110 }
1111
1112   //__________________________________________________________________________
1113 void AliMUONTrackK::MSThin(Int_t sign)
1114 {
1115 /// Adds multiple scattering in a thin layer (only angles are affected)
1116   Double_t cosAlph, cosBeta, momentum, velo, path, theta0;
1117
1118   // check whether the Invert method returns flag if matrix cannot be inverted,
1119   // and do not calculate the Determinant in that case !!!!
1120   if (fWeight->Determinant() != 0) {
1121     Int_t ifail;
1122     mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1123   } else {
1124     AliWarning(" Determinant fWeight=0:");
1125   }
1126
1127   cosAlph = TMath::Cos((*fTrackParNew)(2,0));
1128   cosBeta = TMath::Cos((*fTrackParNew)(3,0));
1129   momentum = 1/(*fTrackParNew)(4,0); // particle momentum
1130   //velo = momentum/TMath::Sqrt(momentum*momentum+muonMass*muonMass); // velocity/c for muon hypothesis
1131   velo = 1; // relativistic
1132   path = TMath::Abs(fgTrackReconstructor->GetChamberThicknessInX0()/cosAlph/cosBeta); // path length
1133   theta0 = 0.0136/velo/momentum*TMath::Sqrt(path)*(1+0.038*TMath::Log(path)); // projected scattering angle
1134
1135   (*fWeight)(2,2) += sign*theta0/cosBeta*theta0/cosBeta; // alpha
1136   (*fWeight)(3,3) += sign*theta0*theta0; // beta
1137   Int_t ifail;
1138   mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1139   return;
1140 }
1141
1142   //__________________________________________________________________________
1143 void AliMUONTrackK::StartBack(void)
1144 {
1145 /// Starts backpropagator
1146
1147   fBPFlag = kTRUE;
1148   fChi2 = 0;
1149   for (Int_t i=0; i<fgkSize; i++) {
1150     for (Int_t j=0; j<fgkSize; j++) {
1151       if (j==i) (*fWeight)(i,i) /= 100;
1152       //if (j==i) (*fWeight)(i,i) /= fNmbTrackHits*fNmbTrackHits;
1153       else (*fWeight)(j,i) = 0;
1154     }
1155   }
1156   // Sort hits on track in descending order in abs(z)
1157   SortHits(0, fTrackHits);
1158 }
1159
1160   //__________________________________________________________________________
1161 void AliMUONTrackK::SortHits(Int_t iflag, TObjArray *array)
1162 {
1163 /// Sort hits in Z if the seed segment in the last but one station
1164 /// (if iflag==0 in descending order in abs(z), if !=0 - unsort)
1165
1166   if (iflag && ((AliMUONHitForRec*)(array->UncheckedAt(0)))->GetChamberNumber() == 6) return;
1167   Double_t z = 0, zmax = TMath::Abs(((AliMUONHitForRec*)(array->UncheckedAt(0)))->GetZ());
1168   Int_t i = 1, entries = array->GetEntriesFast();
1169   for ( ; i<entries; i++) {
1170     if (iflag) {
1171       if (((AliMUONHitForRec*)(array->UncheckedAt(i)))->GetChamberNumber() == 6) break;
1172     } else {
1173       z = TMath::Abs(((AliMUONHitForRec*)(array->UncheckedAt(i)))->GetZ());
1174       if (z < zmax) break;
1175       zmax = z;
1176       if (fgDebug >= 10) cout << " " << zmax << " " << z << " " << i << endl;
1177     }
1178   }
1179   if (!iflag) i--;
1180   for (Int_t j=0; j<=(i-1)/2; j++) {
1181     TObject *hit = array->UncheckedAt(j);
1182     array->AddAt(array->UncheckedAt(i-j),j);
1183     array->AddAt(hit,i-j);
1184   }
1185   if (fgDebug >= 10) {
1186     for (i=0; i<entries; i++)
1187       cout << ((AliMUONHitForRec*)(array->UncheckedAt(i)))->GetZ() << " ";
1188     cout << " - Sort" << endl;
1189   }
1190 }
1191
1192   //__________________________________________________________________________
1193 void AliMUONTrackK::SetGeantParam(Double_t *VGeant3, Int_t iFB)
1194 {
1195 /// Set vector of Geant3 parameters pointed to by "VGeant3"
1196 /// from track parameters
1197
1198   VGeant3[0] = (*fTrackParNew)(1,0); // X
1199   VGeant3[1] = (*fTrackParNew)(0,0); // Y
1200   VGeant3[2] = fPositionNew; // Z
1201   VGeant3[3] = iFB*TMath::Sin((*fTrackParNew)(3,0)); // Px/Ptot
1202   VGeant3[4] = iFB*TMath::Cos((*fTrackParNew)(3,0))*TMath::Sin((*fTrackParNew)(2,0)); // Py/Ptot
1203   VGeant3[5] = iFB*TMath::Sqrt(1.0-VGeant3[3]*VGeant3[3]-VGeant3[4]*VGeant3[4]); // Pz/Ptot
1204   VGeant3[6] = 1/TMath::Abs((*fTrackParNew)(4,0)); // Ptot
1205 }
1206
1207   //__________________________________________________________________________
1208 void AliMUONTrackK::GetFromGeantParam(Double_t *VGeant3, Int_t iFB)
1209 {
1210 /// Get track parameters from vector of Geant3 parameters pointed
1211 /// to by "VGeant3"
1212
1213   fPositionNew = VGeant3[2]; // Z
1214   (*fTrackParNew)(0,0) = VGeant3[1]; // Y
1215   (*fTrackParNew)(1,0) = VGeant3[0]; // X
1216   (*fTrackParNew)(3,0) = TMath::ASin(iFB*VGeant3[3]); // beta
1217   (*fTrackParNew)(2,0) = TMath::ASin(iFB*VGeant3[4]/TMath::Cos((*fTrackParNew)(3,0))); // alpha
1218   (*fTrackParNew)(4,0) = 1/VGeant3[6]*TMath::Sign(Double_t(1.0),(*fTrackPar)(4,0)); // 1/Ptot
1219 }
1220
1221   //__________________________________________________________________________
1222 void AliMUONTrackK::SetTrackQuality(Int_t iChi2)
1223 {
1224 /// Computes "track quality" from Chi2 (if iChi2==0) or vice versa
1225
1226   if (fChi2 > 500) {
1227     AliWarning(Form(" *** Too high Chi2: %f ", fChi2));
1228     fChi2 = 500;
1229   }
1230   if (iChi2 == 0) fChi2 = fNmbTrackHits + (500.-fChi2)/501;
1231   else fChi2 = 500 - (fChi2-fNmbTrackHits)*501;
1232 }
1233
1234   //__________________________________________________________________________
1235 Int_t AliMUONTrackK::Compare(const TObject* trackK) const
1236 {
1237 /// "Compare" function to sort with decreasing "track quality".
1238 /// Returns +1 (0, -1) if quality of current track
1239 /// is smaller than (equal to, larger than) quality of trackK
1240
1241   if (fChi2 < ((AliMUONTrackK*)trackK)->fChi2) return(+1);
1242   else if (fChi2 == ((AliMUONTrackK*)trackK)->fChi2) return(0);
1243   else return(-1);
1244 }
1245
1246   //__________________________________________________________________________
1247 Bool_t AliMUONTrackK::KeepTrack(AliMUONTrackK* track0) const
1248 {
1249 /// Check whether or not to keep current track
1250 /// (keep, if it has less than half of common hits with track0)
1251   Int_t hitsInCommon, nHits0, i, j, nTrackHits2;
1252   AliMUONHitForRec *hit0, *hit1;
1253
1254   hitsInCommon = 0;
1255   nHits0 = track0->fNmbTrackHits;
1256   nTrackHits2 = fNmbTrackHits/2;
1257
1258   for (i=0; i<nHits0; i++) {
1259     // Check if hit belongs to several tracks
1260     hit0 = (AliMUONHitForRec*) (*track0->fTrackHits)[i];
1261     if (hit0->GetNTrackHits() == 1) continue;
1262     for (j=0; j<fNmbTrackHits; j++) {
1263       hit1 = (AliMUONHitForRec*) (*fTrackHits)[j];
1264       if (hit1->GetNTrackHits() == 1) continue;
1265       if (hit0 == hit1) {
1266         hitsInCommon++;
1267         if (hitsInCommon >= nTrackHits2) return kFALSE;
1268         break;
1269       }
1270     } // for (j=0;
1271   } // for (i=0;
1272   return kTRUE;
1273 }
1274
1275   //__________________________________________________________________________
1276 void AliMUONTrackK::Kill(void)
1277 {
1278 /// Kill track candidate
1279   fgTrackReconstructor->GetRecTracksPtr()->Remove(this);
1280 }
1281
1282   //__________________________________________________________________________
1283 void AliMUONTrackK::FillMUONTrack(void)
1284 {
1285 /// Compute track parameters at hit positions (as for AliMUONTrack)
1286
1287   // Set number of hits per track
1288   SetNTrackHits(fNmbTrackHits);
1289   // Set Chi2
1290   SetFitFMin(fChi2);
1291
1292   // Set track parameters at vertex
1293   AliMUONTrackParam trackParam;
1294   SetTrackParam(&trackParam, fTrackPar, fPosition);
1295   SetTrackParamAtVertex(&trackParam);
1296
1297   // Set track parameters at hits
1298   for (Int_t i = fNmbTrackHits-1; i>=0; i--) {
1299     if ((*fChi2Smooth)[i] < 0) {
1300       // Propagate through last chambers
1301       trackParam.ExtrapToZ(((AliMUONHitForRec*)((*fTrackHits)[i]))->GetZ());
1302     } else {
1303       // Take saved info
1304       SetTrackParam(&trackParam, (TMatrixD*)fParSmooth->UncheckedAt(i), ((AliMUONHitForRec*)((*fTrackHits)[i]))->GetZ());
1305     }
1306     AddTrackParamAtHit(&trackParam);
1307     // Fill array of HitForRec's
1308     AddHitForRecAtHit((AliMUONHitForRec*)fTrackHits->UncheckedAt(i));
1309   }
1310 }
1311
1312   //__________________________________________________________________________
1313 void AliMUONTrackK::SetTrackParam(AliMUONTrackParam *trackParam, TMatrixD *par, Double_t z)
1314 {
1315 /// Fill AliMUONTrackParam object
1316
1317   trackParam->SetBendingCoor((*par)(0,0));
1318   trackParam->SetNonBendingCoor((*par)(1,0));
1319   trackParam->SetBendingSlope(TMath::Tan((*par)(2,0)));
1320   trackParam->SetNonBendingSlope(TMath::Tan((*par)(3,0))/TMath::Cos((*par)(2,0)));
1321   trackParam->SetInverseBendingMomentum((*par)(4,0)/TMath::Cos((*par)(3,0)));
1322   trackParam->SetZ(z);
1323 }
1324
1325   //__________________________________________________________________________
1326 void AliMUONTrackK::Branson(void)
1327 {
1328 /// Propagates track to the vertex thru absorber using Branson correction
1329 /// (makes use of the AliMUONTrackParam class)
1330
1331   //AliMUONTrackParam *trackParam = new AliMUONTrackParam();
1332   AliMUONTrackParam trackParam = AliMUONTrackParam();
1333   /*
1334   trackParam->SetBendingCoor((*fTrackPar)(0,0));
1335   trackParam->SetNonBendingCoor((*fTrackPar)(1,0));
1336   trackParam->SetBendingSlope(TMath::Tan((*fTrackPar)(2,0)));
1337   trackParam->SetNonBendingSlope(TMath::Tan((*fTrackPar)(3,0))/TMath::Cos((*fTrackPar)(2,0)));
1338   trackParam->SetInverseBendingMomentum((*fTrackPar)(4,0)/TMath::Cos((*fTrackPar)(3,0)));
1339   trackParam->SetZ(fPosition);
1340   */
1341   SetTrackParam(&trackParam, fTrackPar, fPosition);
1342
1343   trackParam.ExtrapToVertex(Double_t(0.), Double_t(0.), Double_t(0.));
1344   //trackParam.ExtrapToVertex();
1345
1346   (*fTrackPar)(0,0) = trackParam.GetBendingCoor();
1347   (*fTrackPar)(1,0) = trackParam.GetNonBendingCoor();
1348   (*fTrackPar)(2,0) = TMath::ATan(trackParam.GetBendingSlope());
1349   (*fTrackPar)(3,0) = TMath::ATan(TMath::Cos((*fTrackPar)(2,0))*trackParam.GetNonBendingSlope());
1350   (*fTrackPar)(4,0) = TMath::Cos((*fTrackPar)(3,0))*trackParam.GetInverseBendingMomentum();
1351   fPosition = trackParam.GetZ();
1352   //delete trackParam;
1353   if (fgDebug > 0) cout << 1/(*fTrackPar)(4,0) << " " << fPosition << " " << (*fTrackPar)(0,0) << endl;
1354
1355   // Get covariance matrix
1356   *fCovariance = *fWeight;
1357   if (fCovariance->Determinant() != 0) {
1358     Int_t ifail;
1359     mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1360   } else {
1361     AliWarning(" Determinant fCovariance=0:");
1362   }
1363 }
1364
1365   //__________________________________________________________________________
1366 void AliMUONTrackK::GoToZ(Double_t zEnd)
1367 {
1368 /// Propagates track to given Z
1369
1370   ParPropagation(zEnd);
1371   MSThin(1); // multiple scattering in the chamber
1372   WeightPropagation(zEnd, kFALSE);
1373   fPosition = fPositionNew;
1374   *fTrackPar = *fTrackParNew;
1375 }
1376
1377   //__________________________________________________________________________
1378 void AliMUONTrackK::GoToVertex(Int_t iflag)
1379 {
1380 /// Version 3.08
1381 /// Propagates track to the vertex
1382 /// All material constants are taken from AliRoot
1383
1384     static Double_t x01[5] = { 24.282,  // C
1385                                24.282,  // C
1386                                11.274,  // Concrete
1387                                 1.758,  // Fe
1388                                 1.758}; // Fe (cm)
1389   // inner part theta < 3 degrees
1390     static Double_t x02[5] = { 30413,  // Air
1391                                24.282, // C
1392                                11.274, // Concrete
1393                                1.758,  // Fe
1394                                0.369}; // W (cm)
1395   // z positions of the materials inside the absober outer part theta > 3 degres
1396   static Double_t zPos[10] = {-90, -105, -315, -443, -468};
1397
1398   Double_t dZ, r0Norm, x0, deltaP, dChi2, pTotal, pOld;
1399   AliMUONHitForRec *hit;
1400   AliMUONRawCluster *clus;
1401   TClonesArray *rawclusters;
1402
1403   // First step to the rear end of the absorber
1404   Double_t zRear = -503;
1405   GoToZ(zRear);
1406   Double_t tan3 = TMath::Tan(3./180*TMath::Pi());
1407
1408   // Go through absorber
1409   pOld = 1/(*fTrackPar)(4,0);
1410   Double_t r0Rear = (*fTrackPar)(0,0)*(*fTrackPar)(0,0) +
1411                     (*fTrackPar)(1,0)*(*fTrackPar)(1,0);
1412   r0Rear = TMath::Sqrt(r0Rear)/TMath::Abs(fPosition)/tan3;
1413   r0Norm = r0Rear;
1414   Double_t p0, cos25, cos60;
1415   if (!iflag) goto vertex;
1416
1417   for (Int_t i=4; i>=0; i--) {
1418     ParPropagation(zPos[i]);
1419     WeightPropagation(zPos[i], kFALSE);
1420     dZ = TMath::Abs (fPositionNew-fPosition);
1421     if (r0Norm > 1) x0 = x01[i];
1422     else x0 = x02[i];
1423     MSLine(dZ,x0); // multiple scattering in the medium (linear approximation)
1424     fPosition = fPositionNew;
1425     *fTrackPar = *fTrackParNew;
1426     r0Norm = (*fTrackPar)(0,0)*(*fTrackPar)(0,0) +
1427              (*fTrackPar)(1,0)*(*fTrackPar)(1,0);
1428     r0Norm = TMath::Sqrt(r0Norm)/TMath::Abs(fPosition)/tan3;
1429   }
1430   // Correct momentum for energy losses
1431   pTotal = 1/TMath::Abs((*fTrackPar)(4,0));
1432   p0 = pTotal;
1433   cos25 = TMath::Cos(2.5/180*TMath::Pi());
1434   cos60 = TMath::Cos(6.0/180*TMath::Pi());
1435   for (Int_t j=0; j<1; j++) {
1436     /*
1437     if (r0Rear > 1) {
1438       if (p0 < 20) {
1439         deltaP = 2.164 + 0.145e-1*p0 - 0.417e-3*p0*p0;
1440       } else {
1441         deltaP = 2.275 + 0.102e-2*p0 - 0.674e-6*p0*p0;
1442       }
1443     } else {
1444       if (p0 < 20) {
1445         deltaP = 2.581 + 0.188e-1*p0 - 0.398e-3*p0*p0;
1446       } else {
1447         deltaP = 2.727 + 0.356e-2*p0 + 0.242e-5*p0*p0;
1448       }
1449     }
1450     */
1451     /*
1452     if (r0Rear < 1) {
1453       //W
1454       if (p0<15) {
1455         deltaP = 2.737 + 0.0494*p0 - 0.001123*p0*p0;
1456       } else {
1457         deltaP = 3.0643 + 0.01346*p0;
1458       }
1459       deltaP *= 0.95;
1460     } else {
1461       //Pb
1462       if (p0<15) {
1463         deltaP  = 2.1380 + 0.0351*p0 - 0.000853*p0*p0;
1464       } else {
1465         deltaP = 2.407 + 0.00702*p0;
1466       }
1467       deltaP *= 0.95;
1468     }
1469     */
1470     /*
1471     if (r0Rear < 1) {
1472       //W
1473       if (p0<18) {
1474         deltaP = 2.439 + 0.806e-1*p0 - 0.500e-2*p0*p0 + 0.106e-3*p0*p0*p0;
1475       } else {
1476         deltaP = 2.767 + 0.742e-2*p0 - 0.196e-4*p0*p0 + 0.403e-7*p0*p0*p0;
1477       }
1478       //deltaP += 0.2;
1479       deltaP *= cos25;
1480     } else {
1481       //Pb
1482       if (p0<18) {
1483         deltaP  = 2.209 + 0.800e-2*p0;
1484       } else {
1485         deltaP = 2.285 + 0.141e-2*p0 - 0.446e-6*p0*p0;
1486       }
1487       //deltaP += 0.2;
1488       deltaP *= cos60;
1489     }
1490     deltaP *= 1.1;
1491     */
1492     //*
1493     if (r0Rear  < 1) {
1494       if (p0 < 20) {
1495         deltaP = 2.5938 + 0.0570 * p0 - 0.001151 * p0 * p0;
1496       } else {
1497         deltaP = 3.0714 + 0.011767 * p0;
1498       }
1499       deltaP *= 0.75;
1500     } else {
1501       if (p0 < 20) {
1502         deltaP  = 2.1207 + 0.05478 * p0 - 0.00145079 * p0 * p0;
1503       } else {
1504         deltaP = 2.6069 + 0.0051705 * p0;
1505       }
1506       deltaP *= 0.9;
1507     }
1508     //*/
1509
1510     p0 = pTotal + deltaP/TMath::Abs(TMath::Cos((*fTrackPar)(2,0))/TMath::Cos((*fTrackPar)(3,0)));
1511   }
1512   (*fTrackPar)(4,0) = 1/p0*TMath::Sign((Double_t)1.,(*fTrackPar)(4,0));
1513
1514   // Go to the vertex
1515 vertex:
1516   ParPropagation((Double_t)0.);
1517   WeightPropagation((Double_t)0., kFALSE);
1518   fPosition = fPositionNew;
1519   //*fTrackPar = *fTrackParNew;
1520   // Add vertex as a hit
1521   TMatrixD pointWeight(fgkSize,fgkSize);
1522   TMatrixD point(fgkSize,1);
1523   TMatrixD trackParTmp = point;
1524   point(0,0) = 0; // vertex coordinate - should be taken somewhere
1525   point(1,0) = 0; // vertex coordinate - should be taken somewhere
1526   pointWeight(0,0) = 1/1.e-3/1.e-3; // 10 um error
1527   pointWeight(1,1) = 1/1.e-3/1.e-3; // 10 um error
1528   TryPoint(point,pointWeight,trackParTmp,dChi2);
1529   *fTrackPar = trackParTmp;
1530   *fWeight += pointWeight;
1531   fChi2 += dChi2; // Chi2
1532   if (fgDebug < 0) return; // no output
1533
1534   cout << pOld << " " << 1/(*fTrackPar)(4,0) << " " << dChi2 << " " << fChi2 << " " << fNmbTrackHits << endl;
1535   for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
1536     hit =  (AliMUONHitForRec*) ((*fTrackHits)[i1]);
1537     printf ("%5d", hit->GetChamberNumber());
1538   }
1539   cout << endl;
1540   if (fgDebug > 0) {
1541     for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
1542       hit =  (AliMUONHitForRec*) ((*fTrackHits)[i1]);
1543       //cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetHitNumber() << " ";
1544       //cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetZ() << " ";
1545       printf ("%5d", fgHitForRec->IndexOf(hit));
1546     }
1547     cout << endl;
1548   }
1549   if (fgTrackReconstructor->GetRecTrackRefHits()) {
1550       // from track ref. hits
1551     for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
1552       hit =  (AliMUONHitForRec*) ((*fTrackHits)[i1]);
1553       cout << hit->GetTTRTrack() + hit->GetTrackRefSignal()*10000 << " ";
1554     }
1555   } else {
1556     // from raw clusters
1557     for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
1558       hit =  (AliMUONHitForRec*) ((*fTrackHits)[i1]);
1559       if (hit->GetHitNumber() < 0) { // combined cluster / track finder
1560         Int_t index = -hit->GetHitNumber() / 100000;
1561         Int_t iPos = -hit->GetHitNumber() - index * 100000;
1562         clus = (AliMUONRawCluster*) fgCombi->DetElem(index-1)->RawClusters()->UncheckedAt(iPos);
1563       } else {
1564         rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
1565         clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
1566       }
1567       printf ("%5d", clus->GetTrack(1)%10000000);
1568     }
1569     cout << endl;
1570     for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
1571       hit =  (AliMUONHitForRec*) ((*fTrackHits)[i1]);
1572       if (hit->GetHitNumber() < 0) { // combined cluster / track finder
1573         Int_t index = -hit->GetHitNumber() / 100000;
1574         Int_t iPos = -hit->GetHitNumber() - index * 100000;
1575         clus = (AliMUONRawCluster*) fgCombi->DetElem(index-1)->RawClusters()->UncheckedAt(iPos);
1576       } else {
1577         rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
1578         clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
1579       }
1580       if (clus->GetTrack(2) != -1) printf ("%5d", clus->GetTrack(2)%10000000);
1581       else printf ("%5s", "    ");
1582     }
1583   }
1584   cout << endl;
1585   for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
1586     //cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetHitNumber() << " ";
1587     cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetZ() << " ";
1588     //cout << fgHitForRec->IndexOf(((AliMUONHitForRec*)((*fTrackHits)[i1]))) << " ";
1589   }
1590   cout << endl;
1591   for (Int_t i1=0; i1<fNmbTrackHits; i1++) printf("%8.4f", (*fChi2Smooth)[i1]);
1592   cout << endl;
1593   cout << "---------------------------------------------------" << endl;
1594
1595   // Get covariance matrix
1596   /* Not needed - covariance matrix is not interesting to anybody
1597   *fCovariance = *fWeight;
1598   if (fCovariance->Determinant() != 0) {
1599     Int_t ifail;
1600     mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1601   } else {
1602     AliWarning(" Determinant fCovariance=0:" );
1603   }
1604   */
1605 }
1606
1607   //__________________________________________________________________________
1608 void AliMUONTrackK::MSLine(Double_t dZ, Double_t x0)
1609 {
1610 /// Adds multiple scattering in a thick layer for linear propagation
1611
1612   Double_t cosAlph = TMath::Cos((*fTrackPar)(2,0));
1613   Double_t tanAlph = TMath::Tan((*fTrackPar)(2,0));
1614   Double_t cosBeta = TMath::Cos((*fTrackPar)(3,0));
1615   Double_t sinBeta;
1616   sinBeta = TMath::Sin((*fTrackPar)(3,0));
1617   Double_t tanBeta = TMath::Tan((*fTrackPar)(3,0));
1618   Double_t momentum = 1/(*fTrackPar)(4,0);
1619   Double_t velo = 1; // relativistic velocity
1620   Double_t step = TMath::Abs(dZ/cosAlph/cosBeta); // step length
1621
1622   // Projected scattering angle
1623   Double_t theta0 = 0.0136/velo/momentum/TMath::Sqrt(x0)*(1+0.038*TMath::Log(step/x0));
1624   Double_t theta02 = theta0*theta0;
1625   Double_t dl2 = step*step/2*theta02;
1626   Double_t dl3 = dl2*step*2/3;
1627
1628   //Derivatives
1629   Double_t dYdT = 1/cosAlph;
1630   Double_t dYdB = 0; //(*fTrackPar)(2,0)*sinBeta/cosAlph;
1631   Double_t dXdT = tanAlph*tanBeta;
1632   //Double_t dXdB = (1+(*fTrackPar)(2,0)*tanAlph*sinBeta*sinBeta)/cosBeta;
1633   Double_t dXdB = 1/cosBeta;
1634   Double_t dAdT = 1/cosBeta;
1635   Double_t dAdB = 0; //(*fTrackPar)(2,0)*tanBeta;
1636
1637   // Get covariance matrix
1638   *fCovariance = *fWeight;
1639   if (fCovariance->Determinant() != 0) {
1640     //   fCovariance->Invert();
1641     Int_t ifail;
1642     mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1643   } else {
1644     AliWarning(" Determinant fCovariance=0:" );
1645   }
1646
1647   (*fCovariance)(0,0) += dl3*(dYdT*dYdT+dYdB*dYdB); // <yy>
1648   (*fCovariance)(1,1) += dl3*(dXdT*dXdT+dXdB*dXdB); // <xx>
1649   (*fCovariance)(2,2) += theta02*step*(dAdT*dAdT+dAdB*dAdB); // <aa>
1650   (*fCovariance)(3,3) += theta02*step; // <bb>
1651
1652   (*fCovariance)(0,1) += dl3*(dYdT*dXdT+dYdB*dXdB); // <yx>
1653   (*fCovariance)(1,0) = (*fCovariance)(0,1);
1654
1655   (*fCovariance)(0,2) += dl2*(dYdT*dAdT+dYdB*dAdB); // <ya>
1656   (*fCovariance)(2,0) = (*fCovariance)(0,2);
1657
1658   (*fCovariance)(0,3) += dl2*dYdB; // <yb>
1659   (*fCovariance)(3,0) = (*fCovariance)(0,3);
1660
1661   (*fCovariance)(1,2) += dl2*(-dXdT*dAdT+dXdB*dAdB); // <xa>
1662   (*fCovariance)(2,1) = (*fCovariance)(1,2);
1663
1664   (*fCovariance)(1,3) += dl2*dXdB; // <xb>
1665   (*fCovariance)(3,1) = (*fCovariance)(1,3);
1666
1667   (*fCovariance)(2,3) += theta02*step*dAdB; // <ab>
1668   (*fCovariance)(3,2) = (*fCovariance)(2,3);
1669
1670   // Get weight matrix
1671   *fWeight = *fCovariance;
1672   if (fWeight->Determinant() != 0) {
1673     Int_t ifail;
1674     mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1675   } else {
1676     AliWarning(" Determinant fWeight=0:");
1677   }
1678 }
1679
1680   //__________________________________________________________________________
1681 Bool_t AliMUONTrackK::Recover(void)
1682 {
1683 /// Adds new failed track(s) which can be tried to be recovered
1684   Int_t nRecTracks;
1685   TClonesArray *trackPtr;
1686   AliMUONTrackK *trackK;
1687
1688   if (fgDebug > 0) cout << " ******** Enter Recover " << endl;
1689   trackPtr = fgTrackReconstructor->GetRecTracksPtr();
1690
1691   // Remove hit with the highest chi2
1692   Double_t chi2 = 0;
1693   if (fgDebug > 0) {
1694     for (Int_t i=0; i<fNmbTrackHits; i++) {
1695       chi2 = fChi2Smooth ? (*fChi2Smooth)[i] : (*fChi2Array)[i];
1696       printf("%10.4f", chi2);
1697     }
1698     printf("\n");
1699     for (Int_t i=0; i<fNmbTrackHits; i++) {
1700       printf("%10d", ((AliMUONHitForRec*)fTrackHits->UncheckedAt(i))->GetChamberNumber());
1701     }
1702     printf("\n");
1703   }
1704   Double_t chi2max = 0;
1705   Int_t imax = 0;
1706   for (Int_t i=0; i<fNmbTrackHits; i++) {
1707     chi2 = fChi2Smooth ? (*fChi2Smooth)[i] : (*fChi2Array)[i];
1708     if (chi2 < chi2max) continue;
1709     chi2max = chi2;
1710     imax = i;
1711   }
1712   //if (chi2max < 10) return kFALSE; // !!!
1713   //if (chi2max < 25) imax = fNmbTrackHits - 1;
1714   if (chi2max < 15) imax = fNmbTrackHits - 1; // discard the last point
1715   // Check if the outlier is not from the seed segment
1716   AliMUONHitForRec *skipHit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(imax);
1717   if (skipHit == fStartSegment->GetHitForRec1() || skipHit == fStartSegment->GetHitForRec2()) {
1718     //DropBranches(fStartSegment); // drop all tracks with the same seed segment
1719     return kFALSE; // to be changed probably
1720   }
1721
1722   // Make a copy of track hit collection
1723   TObjArray *hits = new TObjArray(*fTrackHits);
1724   Int_t imax0;
1725   imax0 = imax;
1726
1727   // Hits after the found one will be removed
1728   if (GetStation0() == 3 && skipHit->GetChamberNumber() >= 7) {
1729     SortHits(1, fTrackHits); // unsort hits
1730     imax = fTrackHits->IndexOf(skipHit);
1731   }
1732   Int_t nTrackHits = fNmbTrackHits;
1733   for (Int_t i=imax+1; i<nTrackHits; i++) {
1734     AliMUONHitForRec *hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(i);
1735     fTrackHits->Remove(hit);
1736     hit->SetNTrackHits(hit->GetNTrackHits()-1); // unmark hit
1737     fNmbTrackHits--;
1738   }
1739
1740   // Check if the track candidate doesn't exist yet
1741   if (ExistDouble()) { delete hits; return kFALSE; }
1742
1743   //DropBranches(imax0, hits); // drop branches downstream the discarded hit
1744   delete hits;
1745
1746   nRecTracks = fgTrackReconstructor->GetNRecTracks();
1747   skipHit = (AliMUONHitForRec*) ((*fTrackHits)[fNmbTrackHits-1]);
1748   // Remove all saved steps and smoother matrices after the skipped hit
1749   RemoveMatrices(skipHit->GetZ());
1750
1751   //AZ(z->-z) if (skipHit->GetZ() > fStartSegment->GetHitForRec2()->GetZ() || !fNSteps) {
1752   if (TMath::Abs(skipHit->GetZ()) > TMath::Abs(fStartSegment->GetHitForRec2()->GetZ()) || !fNSteps) {
1753     // Propagation toward high Z or skipped hit next to segment -
1754     // start track from segment
1755     trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(fStartSegment);
1756     fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
1757     trackK->fRecover = 1;
1758     trackK->fSkipHit = skipHit;
1759     trackK->fNmbTrackHits = fNmbTrackHits;
1760     delete trackK->fTrackHits; // not efficient ?
1761     trackK->fTrackHits = new TObjArray(*fTrackHits);
1762     if (fgDebug > 0) cout << nRecTracks << " " << trackK->fRecover << endl;
1763     return kTRUE;
1764   }
1765
1766   trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(NULL, NULL);
1767   *trackK = *this;
1768   fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
1769   //AZ(z->-z) trackK->fTrackDir = -1;
1770   trackK->fTrackDir = 1;
1771   trackK->fRecover = 1;
1772   trackK->fSkipHit = skipHit;
1773   Int_t iD = trackK->fParFilter->Last()->GetUniqueID();
1774   if (iD > 1) {
1775     trackK->fParFilter->Last()->SetUniqueID(iD-1);
1776     CreateMatrix(trackK->fParFilter);
1777   }
1778   *((TMatrixD*)trackK->fParFilter->Last()) = *((TMatrixD*)fParExtrap->Last());
1779   trackK->fParFilter->Last()->SetUniqueID(1);
1780   *(trackK->fTrackPar) = *((TMatrixD*)trackK->fParFilter->Last());
1781   iD = trackK->fCovFilter->Last()->GetUniqueID();
1782   if (iD > 1) {
1783     trackK->fCovFilter->Last()->SetUniqueID(iD-1);
1784     CreateMatrix(trackK->fCovFilter);
1785   }
1786   *((TMatrixD*)trackK->fCovFilter->Last()) = *((TMatrixD*)fCovExtrap->Last());
1787   trackK->fCovFilter->Last()->SetUniqueID(1);
1788   *(trackK->fWeight) = *((TMatrixD*)trackK->fCovFilter->Last());
1789   if (trackK->fWeight->Determinant() != 0) {
1790     Int_t ifail;
1791     mnvertLocal(&((*(trackK->fWeight))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1792   } else {
1793     AliWarning(" Determinant fWeight=0:");
1794   }
1795   trackK->fPosition = trackK->fPositionNew = (*fSteps)[fNSteps-1];
1796   trackK->fChi2 = 0;
1797   for (Int_t i=0; i<fNmbTrackHits-1; i++) trackK->fChi2 += (*fChi2Array)[i];
1798   if (fgDebug > 0) cout << nRecTracks << " " << trackK->fRecover << endl;
1799   return kTRUE;
1800 }
1801
1802   //__________________________________________________________________________
1803 void AliMUONTrackK::AddMatrices(AliMUONTrackK *trackK, Double_t dChi2, AliMUONHitForRec *hitAdd)
1804 {
1805 /// Adds matrices for the smoother and keep Chi2 for the point
1806 /// Track parameters
1807   //trackK->fParFilter->Last()->Print();
1808   Int_t iD = trackK->fParFilter->Last()->GetUniqueID();
1809   if (iD > 1) {
1810     trackK->fParFilter->Last()->SetUniqueID(iD-1);
1811     CreateMatrix(trackK->fParFilter);
1812     iD = 1;
1813   }
1814   *((TMatrixD*)(trackK->fParFilter->Last())) = *(trackK->fTrackPar);
1815   trackK->fParFilter->Last()->SetUniqueID(iD);
1816   if (fgDebug > 1) {
1817     cout << " Add matrices" << " " << fPosition << " " << dChi2 << " " << fgHitForRec->IndexOf(hitAdd) << endl;
1818     //trackK->fTrackPar->Print();
1819     //trackK->fTrackParNew->Print();
1820     trackK->fParFilter->Last()->Print();
1821     cout << " Add matrices" << endl;
1822   }
1823   // Covariance
1824   *(trackK->fCovariance) = *(trackK->fWeight);
1825   if (trackK->fCovariance->Determinant() != 0) {
1826     Int_t ifail;
1827     mnvertLocal(&((*(trackK->fCovariance))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1828   } else {
1829     AliWarning(" Determinant fCovariance=0:");
1830   }
1831   iD = trackK->fCovFilter->Last()->GetUniqueID();
1832   if (iD > 1) {
1833     trackK->fCovFilter->Last()->SetUniqueID(iD-1);
1834     CreateMatrix(trackK->fCovFilter);
1835     iD = 1;
1836   }
1837   *((TMatrixD*)(trackK->fCovFilter->Last())) = *(trackK->fCovariance);
1838   trackK->fCovFilter->Last()->SetUniqueID(iD);
1839
1840   // Save Chi2-increment for point
1841   Int_t indx = trackK->fTrackHits->IndexOf(hitAdd);
1842   if (indx < 0) indx = fNmbTrackHits;
1843   if (trackK->fChi2Array->fN <= indx) trackK->fChi2Array->Set(indx+10);
1844   trackK->fChi2Array->AddAt(dChi2,indx);
1845 }
1846
1847   //__________________________________________________________________________
1848 void AliMUONTrackK::CreateMatrix(TObjArray *objArray) const
1849 {
1850 /// Create new matrix and add it to TObjArray
1851
1852   TMatrixD *matrix = (TMatrixD*) objArray->First();
1853   TMatrixD *tmp = new TMatrixD(*matrix);
1854   objArray->AddAtAndExpand(tmp,objArray->GetLast());
1855   //cout << " !!! " << tmp << " " << tmp->GetUniqueID() << " " << tmp->GetNoElements() << endl;
1856 }
1857
1858   //__________________________________________________________________________
1859 void AliMUONTrackK::RemoveMatrices(Double_t zEnd)
1860 {
1861 /// Remove matrices (and saved steps) in the smoother part with abs(z) < abs(zEnd)
1862
1863   for (Int_t i=fNSteps-1; i>=0; i--) {
1864     if (fgDebug > 1) printf("%10.4f %10.4f \n", (*fSteps)[i], zEnd);
1865     if (TMath::Abs((*fSteps)[i]) > TMath::Abs(zEnd) - 0.01) break;
1866     RemoveMatrices(this);
1867   } // for (Int_t i=fNSteps-1;
1868 }
1869
1870   //__________________________________________________________________________
1871 void AliMUONTrackK::RemoveMatrices(AliMUONTrackK* trackK)
1872 {
1873 /// Remove last saved matrices and steps in the smoother part
1874
1875   trackK->fNSteps--;
1876   Int_t i = trackK->fNSteps;
1877
1878   Int_t id = 0;
1879   // Delete only matrices not shared by several tracks
1880   id = trackK->fParExtrap->Last()->GetUniqueID();
1881   if (id > 1) {
1882     trackK->fParExtrap->Last()->SetUniqueID(id-1);
1883     trackK->fParExtrap->RemoveAt(i);
1884   }
1885   else ((TMatrixD*)(trackK->fParExtrap->RemoveAt(i)))->Delete();
1886   id = fParFilter->Last()->GetUniqueID();
1887   if (id > 1) {
1888     trackK->fParFilter->Last()->SetUniqueID(id-1);
1889     trackK->fParFilter->RemoveAt(i);
1890   }
1891   else ((TMatrixD*)(trackK->fParFilter->RemoveAt(i)))->Delete();
1892   id = trackK->fCovExtrap->Last()->GetUniqueID();
1893   if (id > 1) {
1894     trackK->fCovExtrap->Last()->SetUniqueID(id-1);
1895     trackK->fCovExtrap->RemoveAt(i);
1896   }
1897   else ((TMatrixD*)(trackK->fCovExtrap->RemoveAt(i)))->Delete();
1898   id = trackK->fCovFilter->Last()->GetUniqueID();
1899   if (id > 1) {
1900     trackK->fCovFilter->Last()->SetUniqueID(id-1);
1901     trackK->fCovFilter->RemoveAt(i);
1902   }
1903   else ((TMatrixD*)(trackK->fCovFilter->RemoveAt(i)))->Delete();
1904   id = trackK->fJacob->Last()->GetUniqueID();
1905   if (id > 1) {
1906     trackK->fJacob->Last()->SetUniqueID(id-1);
1907     trackK->fJacob->RemoveAt(i);
1908   }
1909   else ((TMatrixD*)(trackK->fJacob->RemoveAt(i)))->Delete();
1910 }
1911
1912   //__________________________________________________________________________
1913 Bool_t AliMUONTrackK::Smooth(void)
1914 {
1915 /// Apply smoother
1916   Int_t ihit = fNmbTrackHits - 1;
1917   AliMUONHitForRec *hit = (AliMUONHitForRec*) (*fTrackHits)[ihit];
1918   fChi2Smooth = new TArrayD(fNmbTrackHits);
1919   fChi2Smooth->Reset(-1);
1920   fChi2 = 0;
1921   fParSmooth = new TObjArray(15);
1922   fParSmooth->Clear();
1923
1924   if (fgDebug > 0) {
1925     cout << " ******** Enter Smooth " << endl;
1926     cout << (*fSteps)[fNSteps-1] << " " << fPosition << " " << fPositionNew << endl;
1927     /*
1928     for (Int_t i=fNSteps-1; i>=0; i--) cout << (*fSteps)[i] << " ";
1929     cout << endl;
1930     for (Int_t i=fNSteps-1; i>=0; i--) {cout << i << " " << (*fSteps)[i]; ((TMatrixD*)fParFilter->UncheckedAt(i))->Print(); ((TMatrixD*)fParExtrap->UncheckedAt(i))->Print(); ((TMatrixD*)fJacob->UncheckedAt(i))->Print();}
1931     */
1932     for (Int_t i=ihit; i>=0; i--) cout << ((AliMUONHitForRec*)(*fTrackHits)[i])->GetZ() << " ";
1933     cout << endl;
1934   }
1935
1936   // Find last point corresponding to the last hit
1937   Int_t iLast = fNSteps - 1;
1938   for ( ; iLast>=0; iLast--) {
1939     //AZ(z->-z) if ((*fSteps)[iLast] + 0.01 > ((AliMUONHitForRec*)(*fTrackHits)[fNmbTrackHits-1])->GetZ()) break;
1940     if (TMath::Abs((*fSteps)[iLast]) + 0.01 > TMath::Abs(((AliMUONHitForRec*)(*fTrackHits)[fNmbTrackHits-1])->GetZ())) break;
1941   }
1942
1943   TMatrixD parSmooth = *((TMatrixD*)(fParFilter->UncheckedAt(iLast))); // last filtered = first smoothed
1944   //parSmooth.Dump();
1945   TMatrixD covSmooth = *((TMatrixD*)(fCovFilter->UncheckedAt(iLast))); // last filtered = first smoothed
1946   TMatrixD tmp=covSmooth, weight=covSmooth, gain = covSmooth;
1947   TMatrixD tmpPar = *fTrackPar;
1948   //parSmooth.Print(); ((TMatrixD*)(fParExtrap->Last()))->Print();
1949
1950   Bool_t found;
1951   Double_t chi2max = 0;
1952   for (Int_t i=iLast+1; i>0; i--) {
1953     if (i == iLast + 1) goto L33; // temporary fix
1954
1955     // Smoother gain matrix
1956     weight = *((TMatrixD*)(fCovExtrap->UncheckedAt(i)));
1957     if (weight.Determinant() != 0) {
1958       Int_t ifail;
1959       mnvertLocal(&(weight(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1960     } else {
1961       AliWarning(" Determinant weight=0:");
1962     }
1963     // Fk'Wkk+1
1964     tmp = TMatrixD(*((TMatrixD*)(fJacob->UncheckedAt(i))),TMatrixD::kTransposeMult,weight);
1965     gain = TMatrixD(*((TMatrixD*)(fCovFilter->UncheckedAt(i-1))),TMatrixD::kMult,tmp);
1966     //cout << (*fSteps)[i-1] << " " << (*fSteps)[i] << endl; gain.Print();
1967
1968     // Smoothed parameter vector
1969     //((TMatrixD*)(fParExtrap->UncheckedAt(i)))->Dump();
1970     parSmooth -= *((TMatrixD*)(fParExtrap->UncheckedAt(i)));
1971     tmpPar = *((TMatrixD*)(fParFilter->UncheckedAt(i-1)));
1972     tmpPar += TMatrixD(gain,TMatrixD::kMult,parSmooth);
1973     parSmooth = tmpPar;
1974
1975     // Smoothed covariance
1976     covSmooth -= *((TMatrixD*)(fCovExtrap->UncheckedAt(i)));
1977     weight = TMatrixD(TMatrixD::kTransposed,gain);
1978     tmp = TMatrixD(covSmooth,TMatrixD::kMult,weight);
1979     covSmooth = *((TMatrixD*)(fCovFilter->UncheckedAt(i-1)));
1980     covSmooth += TMatrixD(gain,TMatrixD::kMult,tmp);
1981
1982     // Check if there was a measurement at given z
1983     found = kFALSE;
1984     for ( ; ihit>=0; ihit--) {
1985       hit = (AliMUONHitForRec*) (*fTrackHits)[ihit];
1986       if (TMath::Abs(hit->GetZ()-(*fSteps)[i-1]) < 0.1) { found = kTRUE; break; }
1987       //AZ(z->-z) else if (ihit < fNmbTrackHits-1 && hit->GetZ() > (*fSteps)[i-1]) { ihit++; break; }
1988       else if (ihit < fNmbTrackHits-1 && TMath::Abs(hit->GetZ()) > TMath::Abs((*fSteps)[i-1])) { ihit++; break; }
1989     }
1990     if (!found) continue; // no measurement - skip the rest
1991     else if (fgDebug > 1) cout << i << " " << ihit << " " << hit->GetZ() << endl;
1992     if (ihit == 0) continue; // the first hit - skip the rest
1993
1994 L33:
1995     // Smoothed residuals
1996     tmpPar = 0;
1997     tmpPar(0,0) = hit->GetBendingCoor() - parSmooth(0,0); // bending coordinate
1998     tmpPar(1,0) = hit->GetNonBendingCoor() - parSmooth(1,0); // non-bending coordinate
1999     if (fgDebug > 1) {
2000       cout << hit->GetBendingCoor() << " " << parSmooth(0,0) << " " << tmpPar(0,0) << endl;
2001       cout << hit->GetNonBendingCoor() << " " << parSmooth(1,0) << " " << tmpPar(1,0) << endl;
2002     }
2003     // Cov. matrix of smoothed residuals
2004     tmp = 0;
2005     tmp(0,0) = hit->GetBendingReso2() - covSmooth(0,0);
2006     tmp(1,1) = hit->GetNonBendingReso2() - covSmooth(1,1);
2007     tmp(0,1) = tmp(1,0) = -covSmooth(0,1);
2008     tmp(2,2) = tmp(3,3) = tmp(4,4) = 1;
2009
2010     // Calculate Chi2 of smoothed residuals
2011     if (tmp.Determinant() != 0) {
2012       Int_t ifail;
2013       mnvertLocal(&(tmp(0,0)), fgkSize,fgkSize,fgkSize,ifail);
2014     } else {
2015       AliWarning(" Determinant tmp=0:");
2016     }
2017     TMatrixD vector(tmp,TMatrixD::kMult,tmpPar);
2018     TMatrixD chi2(tmpPar,TMatrixD::kTransposeMult,vector);
2019     if (fgDebug > 1) chi2.Print();
2020     (*fChi2Smooth)[ihit] = chi2(0,0);
2021     if (chi2(0,0) > chi2max) chi2max = chi2(0,0);
2022     fChi2 += chi2(0,0);
2023     if (chi2(0,0) < 0) {
2024       //chi2.Print();
2025       AliError(Form(" *** chi2 < 0: %d %d ", i, iLast));
2026     }
2027     // Save smoothed parameters
2028     TMatrixD *par = new TMatrixD(parSmooth);
2029     fParSmooth->AddAtAndExpand(par, ihit);
2030
2031   } // for (Int_t i=iLast+1;
2032
2033   //if (chi2max > 16) {
2034   //if (chi2max > 25) {
2035   //if (chi2max > 50) {
2036   //if (chi2max > 100) {
2037   if (chi2max > fgkChi2max) {
2038     //if (Recover()) DropBranches();
2039     //Recover();
2040     Outlier();
2041     return kFALSE;
2042   }
2043   return kTRUE;
2044 }
2045
2046   //__________________________________________________________________________
2047 void AliMUONTrackK::Outlier()
2048 {
2049 /// Adds new track with removed hit having the highest chi2
2050
2051   if (fgDebug > 0) {
2052     cout << " ******** Enter Outlier " << endl;
2053     for (Int_t i=0; i<fNmbTrackHits; i++) printf("%10.4f", (*fChi2Smooth)[i]);
2054     printf("\n");
2055     for (Int_t i=0; i<fNmbTrackHits; i++) {
2056       printf("%10d", ((AliMUONHitForRec*)fTrackHits->UncheckedAt(i))->GetChamberNumber());
2057     }
2058     printf("\n");
2059   }
2060
2061   Double_t chi2max = 0;
2062   Int_t imax = 0;
2063   for (Int_t i=0; i<fNmbTrackHits; i++) {
2064     if ((*fChi2Smooth)[i] < chi2max) continue;
2065     chi2max = (*fChi2Smooth)[i];
2066     imax = i;
2067   }
2068   // Check if the outlier is not from the seed segment
2069   AliMUONHitForRec *hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(imax);
2070   if (hit == fStartSegment->GetHitForRec1() || hit == fStartSegment->GetHitForRec2()) return; // to be changed probably
2071
2072   // Check for missing station
2073   Int_t ok = 1;
2074   if (imax == 0) {
2075     if (((AliMUONHitForRec*)fTrackHits->UncheckedAt(1))->GetChamberNumber() < 8) ok--;
2076   } else if (imax == fNmbTrackHits-1) {
2077     if (((AliMUONHitForRec*)fTrackHits->UncheckedAt(fNmbTrackHits-2))->GetChamberNumber() > 1) ok--;
2078   }
2079   else if (((AliMUONHitForRec*)fTrackHits->UncheckedAt(imax-1))->GetChamberNumber()/2 - ((AliMUONHitForRec*)fTrackHits->UncheckedAt(imax+1))->GetChamberNumber()/2 > 1) ok--;
2080   if (!ok) { Recover(); return; } // try to recover track
2081   //AZ if (!ok) { if (fgDebug >= 0) cout << imax << endl; DropBranches(imax, 0); return; }
2082
2083   // Remove saved steps and smoother matrices after the outlier
2084   RemoveMatrices(hit->GetZ());
2085
2086   // Check for possible double track candidates
2087   //if (ExistDouble(hit)) return;
2088
2089   TClonesArray *trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2090   Int_t nRecTracks = fgTrackReconstructor->GetNRecTracks();
2091
2092   AliMUONTrackK *trackK = 0;
2093   if (!fNSteps || GetStation0() == 3 && hit->GetChamberNumber() > 7) {
2094     // start track from segment
2095     trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(fStartSegment);
2096     fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
2097     trackK->fRecover = 2;
2098     trackK->fSkipHit = hit;
2099     trackK->fNmbTrackHits = fNmbTrackHits;
2100
2101     hit = (AliMUONHitForRec*) trackK->fTrackHits->UncheckedAt(0);
2102     hit->SetNTrackHits(hit->GetNTrackHits()-1);
2103     hit = (AliMUONHitForRec*) trackK->fTrackHits->UncheckedAt(1);
2104     hit->SetNTrackHits(hit->GetNTrackHits()-1);
2105     delete trackK->fTrackHits; // not efficient ?
2106     trackK->fTrackHits = new TObjArray(*fTrackHits);
2107     for (Int_t i = 0; i < fNmbTrackHits; i++) {
2108       hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(i);
2109       hit->SetNTrackHits(hit->GetNTrackHits()+1);
2110     }
2111
2112     if (GetStation0() == 3) trackK->SortHits(1, trackK->fTrackHits);
2113     if (fgDebug > 0) cout << nRecTracks << " Outlier" << endl;
2114     return;
2115   }
2116   trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(NULL, NULL);
2117   *trackK = *this;
2118   fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
2119   trackK->fTrackDir = 1;
2120   trackK->fRecover = 2;
2121   trackK->fSkipHit = hit;
2122   Int_t iD = trackK->fParFilter->Last()->GetUniqueID();
2123   if (iD > 1) {
2124     trackK->fParFilter->Last()->SetUniqueID(iD-1);
2125     CreateMatrix(trackK->fParFilter);
2126   }
2127   *((TMatrixD*)trackK->fParFilter->Last()) = *((TMatrixD*)fParExtrap->Last());
2128   trackK->fParFilter->Last()->SetUniqueID(1);
2129   *(trackK->fTrackPar) = *((TMatrixD*)trackK->fParFilter->Last());
2130   iD = trackK->fCovFilter->Last()->GetUniqueID();
2131   if (iD > 1) {
2132     trackK->fCovFilter->Last()->SetUniqueID(iD-1);
2133     CreateMatrix(trackK->fCovFilter);
2134   }
2135   *((TMatrixD*)trackK->fCovFilter->Last()) = *((TMatrixD*)fCovExtrap->Last());
2136   trackK->fCovFilter->Last()->SetUniqueID(1);
2137   *(trackK->fWeight) = *((TMatrixD*)trackK->fCovFilter->Last());
2138   if (trackK->fWeight->Determinant() != 0) {
2139     Int_t ifail;
2140     mnvertLocal(&((*(trackK->fWeight))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
2141   } else {
2142     AliWarning(" Determinant fWeight=0:");
2143   }
2144   trackK->fPosition = trackK->fPositionNew = (*fSteps)[fNSteps-1];
2145   trackK->fChi2 = 0;
2146   for (Int_t i=0; i<fNmbTrackHits-1; i++) trackK->fChi2 += (*fChi2Array)[i];
2147   if (fgDebug > 0) cout << nRecTracks << " Outlier " << trackK->fChi2 << endl;
2148 }
2149
2150   //__________________________________________________________________________
2151 Double_t AliMUONTrackK::GetChi2PerPoint(Int_t iPoint) const
2152 {
2153 /// Return Chi2 at point
2154   return fChi2Smooth ? (*fChi2Smooth)[iPoint] : (*fChi2Array)[iPoint];
2155   //return 0.;
2156 }
2157
2158   //__________________________________________________________________________
2159 void AliMUONTrackK::Print(FILE *lun) const
2160 {
2161 /// Print out track information
2162
2163   Int_t flag = 1;
2164   AliMUONHitForRec *hit = 0;
2165   if (fgTrackReconstructor->GetRecTrackRefHits()) {
2166     // from track ref. hits
2167     for (Int_t j=0; j<fNmbTrackHits; j++) {
2168       hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(j);
2169       if (hit->GetTTRTrack() > 1) { flag = 0; break; }
2170     }
2171     for (Int_t j=0; j<fNmbTrackHits; j++) {
2172       printf("%10.4f", GetChi2PerPoint(j));
2173       if (GetChi2PerPoint(j) > -0.1) {
2174         hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(j);
2175         fprintf(lun,"%3d %3d %10.4f", gAlice->GetEvNumber(), hit->GetChamberNumber(), GetChi2PerPoint(j));
2176         fprintf(lun, "%3d %3d %3d \n", hit->GetTrackRefSignal(), hit->GetTTRTrack(), flag);
2177       }
2178     }
2179     printf("\n");
2180   } else {
2181     // from raw clusters
2182     AliMUONRawCluster *clus = 0;
2183     TClonesArray *rawclusters = 0;
2184     for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
2185       hit =  (AliMUONHitForRec*) ((*fTrackHits)[i1]);
2186       rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
2187       clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
2188       if (TMath::Abs(clus->GetTrack(1)-1) < 2) {
2189         if (clus->GetTrack(2)) flag = 2;
2190         continue;
2191       }
2192       if (clus->GetTrack(2) && TMath::Abs(clus->GetTrack(2)-1) < 2) {
2193         flag = 3;
2194         continue;
2195       }
2196       flag = 0;
2197       break;
2198     }
2199     Int_t sig[2]={1,1}, tid[2]={0};
2200     for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
2201       if (GetChi2PerPoint(i1) < -0.1) continue;
2202       hit =  (AliMUONHitForRec*) ((*fTrackHits)[i1]);
2203       rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
2204       clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
2205       for (Int_t j=0; j<2; j++) {
2206         tid[j] = clus->GetTrack(j+1) - 1;
2207         if (clus->GetTrack(j+1) < 0) { sig[j] = 0; tid[j] = 999; }
2208       }
2209       fprintf(lun,"%3d %3d %10.4f", gAlice->GetEvNumber(), hit->GetChamberNumber(), GetChi2PerPoint(i1));
2210       if (!(clus->GetTrack(2))) fprintf(lun, "%3d %3d", sig[0], tid[0]); // simple cluster
2211       else { // track overlap
2212         fprintf(lun, "%3d %3d", TMath::Max(sig[0],sig[1]), TMath::Min(tid[0],tid[1]));
2213         //if (tid[0] < 2) flag *= 2;
2214         //else if (tid[1] < 2) flag *= 3;
2215       }
2216       fprintf (lun, "%3d \n", flag);
2217     }
2218   }
2219 }
2220
2221   //__________________________________________________________________________
2222 void AliMUONTrackK::DropBranches(Int_t imax, TObjArray *hits)
2223 {
2224 /// Drop branches downstream of the skipped hit
2225   Int_t nRecTracks;
2226   TClonesArray *trackPtr;
2227   AliMUONTrackK *trackK;
2228
2229   trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2230   nRecTracks = fgTrackReconstructor->GetNRecTracks();
2231   Int_t icand = trackPtr->IndexOf(this);
2232   if (!hits) hits = fTrackHits;
2233
2234   // Check if the track candidate doesn't exist yet
2235   for (Int_t i=icand+1; i<nRecTracks; i++) {
2236     trackK = (AliMUONTrackK*) ((*trackPtr)[i]);
2237     if (trackK->fNmbTrackHits == 2 && trackK->GetRecover() == 0) continue;
2238     if (trackK->GetRecover() < 0) continue;
2239
2240     if (trackK->fNmbTrackHits >= imax + 1) {
2241       for (Int_t j=0; j<=imax; j++) {
2242         //if (j != fNmbTrackHits-1 && (*trackK->fTrackHits)[j] != (*fTrackHits)[j]) break;
2243         if ((*trackK->fTrackHits)[j] != (*hits)[j]) break;
2244         if (j == imax) {
2245           if (hits != fTrackHits) {
2246             // Drop all branches downstream the hit (for Recover)
2247             trackK->SetRecover(-1);
2248             if (fgDebug >= 0 )cout << " Recover " << i << endl;
2249             continue;
2250           }
2251           // Check if the removal of the hit breaks the track
2252           Int_t ok = 1;
2253           if (imax == 0) {
2254             if (((AliMUONHitForRec*)trackK->fTrackHits->UncheckedAt(1))->GetChamberNumber() < 8) ok--; }
2255           else if (imax == trackK->fNmbTrackHits-1) continue;
2256             // else if (imax == trackK->fNmbTrackHits-1) {
2257             //if (((AliMUONHitForRec*)trackK->fTrackHits->UncheckedAt(trackK->fNmbTrackHits-2))->GetChamberNumber() > 1) ok--;
2258             //}
2259           else if (((AliMUONHitForRec*)trackK->fTrackHits->UncheckedAt(imax-1))->GetChamberNumber()/2 - ((AliMUONHitForRec*)trackK->fTrackHits->UncheckedAt(imax+1))->GetChamberNumber()/2 > 1) ok--;
2260           if (!ok) trackK->SetRecover(-1);
2261         }
2262       } // for (Int_t j=0;
2263     }
2264   } // for (Int_t i=0;
2265 }
2266
2267   //__________________________________________________________________________
2268 void AliMUONTrackK::DropBranches(AliMUONSegment *segment)
2269 {
2270 /// Drop all candidates with the same seed segment
2271   Int_t nRecTracks;
2272   TClonesArray *trackPtr;
2273   AliMUONTrackK *trackK;
2274
2275   trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2276   nRecTracks = fgTrackReconstructor->GetNRecTracks();
2277   Int_t icand = trackPtr->IndexOf(this);
2278
2279   for (Int_t i=icand+1; i<nRecTracks; i++) {
2280     trackK = (AliMUONTrackK*) ((*trackPtr)[i]);
2281     if (trackK->fNmbTrackHits == 2 && trackK->GetRecover() == 0) continue;
2282     if (trackK->GetRecover() < 0) continue;
2283     if (trackK->fStartSegment == segment) trackK->SetRecover(-1);
2284   }
2285   if (fgDebug >= 0) cout << " Drop segment " << endl;
2286 }
2287
2288   //__________________________________________________________________________
2289 AliMUONHitForRec* AliMUONTrackK::GetHitLastOk(void)
2290 {
2291 /// Return the hit where track stopped
2292
2293   if (!fNSteps) return (AliMUONHitForRec*)((*fTrackHits)[1]);
2294   return fSkipHit;
2295 }
2296
2297   //__________________________________________________________________________
2298 Int_t AliMUONTrackK::GetStation0(void)
2299 {
2300 /// Return seed station number
2301   return fStartSegment->GetHitForRec1()->GetChamberNumber() / 2;
2302 }
2303
2304   //__________________________________________________________________________
2305 Bool_t AliMUONTrackK::ExistDouble(AliMUONHitForRec *hit)
2306 {
2307 /// Check if the track will make a double after outlier removal
2308
2309   TClonesArray *trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2310   Int_t nRecTracks = fgTrackReconstructor->GetNRecTracks();
2311   TObjArray *hitArray = new TObjArray(*fTrackHits);
2312   TObjArray *hitArray1 = new TObjArray(*hitArray);
2313   hitArray1->Remove(hit);
2314   hitArray1->Compress();
2315
2316   Bool_t same = kFALSE;
2317   for (Int_t i=0; i<nRecTracks; i++) {
2318     AliMUONTrackK *trackK = (AliMUONTrackK*) ((*trackPtr)[i]);
2319     if (trackK->fNmbTrackHits == 2 && trackK->GetRecover() == 0) continue;
2320     if (trackK == this) continue;
2321     if (trackK->fNmbTrackHits == fNmbTrackHits || trackK->fNmbTrackHits == fNmbTrackHits-1) {
2322       TObjArray *hits = new TObjArray(*trackK->fTrackHits);
2323       same = kTRUE;
2324       if (trackK->fNmbTrackHits == fNmbTrackHits) {
2325         for (Int_t j=0; j<fNmbTrackHits; j++) {
2326           if (hits->UncheckedAt(j) != hitArray->UncheckedAt(j)) { same = kFALSE; break; }
2327         }
2328         if (same) { delete hits; break; }
2329         if (trackK->fSkipHit) {
2330           TObjArray *hits1 = new TObjArray(*hits);
2331           if (hits1->Remove(trackK->fSkipHit) > 0) {
2332             hits1->Compress();
2333             same = kTRUE;
2334             for (Int_t j=0; j<fNmbTrackHits-1; j++) {
2335               if (hits1->UncheckedAt(j) != hitArray1->UncheckedAt(j)) { same = kFALSE; break; }
2336             }
2337             if (same) { delete hits1; break; }
2338           }
2339           delete hits1;
2340         }
2341       } else {
2342         // Check with removed outlier
2343         same = kTRUE;
2344         for (Int_t j=0; j<fNmbTrackHits-1; j++) {
2345           if (hits->UncheckedAt(j) != hitArray1->UncheckedAt(j)) { same = kFALSE; break; }
2346         }
2347         if (same) { delete hits; break; }
2348       }
2349       delete hits;
2350     }
2351   } // for (Int_t i=0; i<nRecTracks;
2352   delete hitArray; delete hitArray1;
2353   if (same && fgDebug >= 0) cout << " Same" << endl;
2354   return same;
2355 }
2356
2357   //__________________________________________________________________________
2358 Bool_t AliMUONTrackK::ExistDouble(void)
2359 {
2360 /// Check if the track will make a double after recovery
2361
2362   TClonesArray *trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2363   Int_t nRecTracks = fgTrackReconstructor->GetNRecTracks();
2364
2365   TObjArray *hitArray = new TObjArray(*fTrackHits);
2366   if (GetStation0() == 3) SortHits(0, hitArray); // sort
2367   //if (GetStation0() == 3) SortHits(1, hitArray); // unsort
2368
2369   Bool_t same = kFALSE;
2370   for (Int_t i=0; i<nRecTracks; i++) {
2371     AliMUONTrackK *trackK = (AliMUONTrackK*) ((*trackPtr)[i]);
2372     if (trackK->fNmbTrackHits == 2 && trackK->GetRecover() == 0) continue;
2373     if (trackK == this) continue;
2374     //AZ if (trackK->GetRecover() < 0) continue; //
2375     if (trackK->fNmbTrackHits >= fNmbTrackHits) {
2376       TObjArray *hits = new TObjArray(*trackK->fTrackHits);
2377       if (trackK->GetStation0() == 3) SortHits(0, hits); // sort
2378       //if (trackK->GetStation0() == 3) SortHits(1, hits); // unsort
2379       for (Int_t j=0; j<fNmbTrackHits; j++) {
2380         //cout << fNmbTrackHits << " " << i << " " << j << " " << (*hitArray)[j] << " " << (*hits)[j] << " " << trackK->fSkipHit << endl;
2381         if (j != fNmbTrackHits-1 && (*hitArray)[j] != (*hits)[j]) break;
2382         if (j == fNmbTrackHits-1) {
2383           if (trackK->fSkipHit && TMath::Abs(((AliMUONHitForRec*)((*hitArray)[j]))->GetZ()-trackK->fSkipHit->GetZ()) < 0.5) same = kTRUE;
2384           //if (trackK->fNmbTrackHits > fNmbTrackHits &&
2385           //if (trackK->fSkipHit) cout << ((AliMUONHitForRec*)((*hitArray)[j]))->GetZ() << " " << trackK->fSkipHit->GetZ() << endl;
2386         }
2387       } // for (Int_t j=0;
2388       delete hits;
2389       if (same) break;
2390     }
2391   } // for (Int_t i=0;
2392   delete hitArray;
2393   return same;
2394 }