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