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