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