MUON/AliMUONTrackK.cxx

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