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