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),
 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 Chi2
1280   SetFitFMin(fChi2);
1281
1282   // Set track parameters at vertex
1283   AliMUONTrackParam trackParam;
1284   SetTrackParam(&trackParam, fTrackPar, fPosition);
1285   SetTrackParamAtVertex(&trackParam);
1286
1287   // Set track parameters at hits
1288   for (Int_t i = fNmbTrackHits-1; i>=0; i--) {
1289     if ((*fChi2Smooth)[i] < 0) {
1290       // Propagate through last chambers
1291       trackParam.ExtrapToZ(((AliMUONHitForRec*)((*fTrackHits)[i]))->GetZ());
1292     } else {
1293       // Take saved info
1294       SetTrackParam(&trackParam, (TMatrixD*)fParSmooth->UncheckedAt(i), ((AliMUONHitForRec*)((*fTrackHits)[i]))->GetZ());
1295     }
1296     AddTrackParamAtHit(&trackParam,(AliMUONHitForRec*)fTrackHits->UncheckedAt(i));
1297     // Fill array of HitForRec's
1298     AddHitForRecAtHit((AliMUONHitForRec*)fTrackHits->UncheckedAt(i));
1299   }
1300 }
1301
1302   //__________________________________________________________________________
1303 void AliMUONTrackK::SetTrackParam(AliMUONTrackParam *trackParam, TMatrixD *par, Double_t z)
1304 {
1305 /// Fill AliMUONTrackParam object
1306
1307   trackParam->SetBendingCoor((*par)(0,0));
1308   trackParam->SetNonBendingCoor((*par)(1,0));
1309   trackParam->SetBendingSlope(TMath::Tan((*par)(2,0)));
1310   trackParam->SetNonBendingSlope(TMath::Tan((*par)(3,0))/TMath::Cos((*par)(2,0)));
1311   trackParam->SetInverseBendingMomentum((*par)(4,0)/TMath::Cos((*par)(3,0)));
1312   trackParam->SetZ(z);
1313 }
1314
1315   //__________________________________________________________________________
1316 void AliMUONTrackK::Branson(void)
1317 {
1318 /// Propagates track to the vertex thru absorber using Branson correction
1319 /// (makes use of the AliMUONTrackParam class)
1320
1321   //AliMUONTrackParam *trackParam = new AliMUONTrackParam();
1322   AliMUONTrackParam trackParam = AliMUONTrackParam();
1323   /*
1324   trackParam->SetBendingCoor((*fTrackPar)(0,0));
1325   trackParam->SetNonBendingCoor((*fTrackPar)(1,0));
1326   trackParam->SetBendingSlope(TMath::Tan((*fTrackPar)(2,0)));
1327   trackParam->SetNonBendingSlope(TMath::Tan((*fTrackPar)(3,0))/TMath::Cos((*fTrackPar)(2,0)));
1328   trackParam->SetInverseBendingMomentum((*fTrackPar)(4,0)/TMath::Cos((*fTrackPar)(3,0)));
1329   trackParam->SetZ(fPosition);
1330   */
1331   SetTrackParam(&trackParam, fTrackPar, fPosition);
1332
1333   trackParam.ExtrapToVertex(Double_t(0.), Double_t(0.), Double_t(0.));
1334   //trackParam.ExtrapToVertex();
1335
1336   (*fTrackPar)(0,0) = trackParam.GetBendingCoor();
1337   (*fTrackPar)(1,0) = trackParam.GetNonBendingCoor();
1338   (*fTrackPar)(2,0) = TMath::ATan(trackParam.GetBendingSlope());
1339   (*fTrackPar)(3,0) = TMath::ATan(TMath::Cos((*fTrackPar)(2,0))*trackParam.GetNonBendingSlope());
1340   (*fTrackPar)(4,0) = TMath::Cos((*fTrackPar)(3,0))*trackParam.GetInverseBendingMomentum();
1341   fPosition = trackParam.GetZ();
1342   //delete trackParam;
1343   if (fgDebug > 0) cout << 1/(*fTrackPar)(4,0) << " " << fPosition << " " << (*fTrackPar)(0,0) << endl;
1344
1345   // Get covariance matrix
1346   *fCovariance = *fWeight;
1347   if (fCovariance->Determinant() != 0) {
1348     Int_t ifail;
1349     mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1350   } else {
1351     AliWarning(" Determinant fCovariance=0:");
1352   }
1353 }
1354
1355   //__________________________________________________________________________
1356 void AliMUONTrackK::GoToZ(Double_t zEnd)
1357 {
1358 /// Propagates track to given Z
1359
1360   ParPropagation(zEnd);
1361   MSThin(1); // multiple scattering in the chamber
1362   WeightPropagation(zEnd, kFALSE);
1363   fPosition = fPositionNew;
1364   *fTrackPar = *fTrackParNew;
1365 }
1366
1367   //__________________________________________________________________________
1368 void AliMUONTrackK::GoToVertex(Int_t iflag)
1369 {
1370 /// Version 3.08
1371 /// Propagates track to the vertex
1372 /// All material constants are taken from AliRoot
1373
1374     static Double_t x01[5] = { 24.282,  // C
1375                                24.282,  // C
1376                                11.274,  // Concrete
1377                                 1.758,  // Fe
1378                                 1.758}; // Fe (cm)
1379   // inner part theta < 3 degrees
1380     static Double_t x02[5] = { 30413,  // Air
1381                                24.282, // C
1382                                11.274, // Concrete
1383                                1.758,  // Fe
1384                                0.369}; // W (cm)
1385   // z positions of the materials inside the absober outer part theta > 3 degres
1386   static Double_t zPos[10] = {-90, -105, -315, -443, -468};
1387
1388   Double_t dZ, r0Norm, x0, deltaP, dChi2, pTotal, pOld;
1389   AliMUONHitForRec *hit;
1390   AliMUONRawCluster *clus;
1391   TClonesArray *rawclusters;
1392
1393   // First step to the rear end of the absorber
1394   Double_t zRear = -503;
1395   GoToZ(zRear);
1396   Double_t tan3 = TMath::Tan(3./180*TMath::Pi());
1397
1398   // Go through absorber
1399   pOld = 1/(*fTrackPar)(4,0);
1400   Double_t r0Rear = (*fTrackPar)(0,0)*(*fTrackPar)(0,0) +
1401                     (*fTrackPar)(1,0)*(*fTrackPar)(1,0);
1402   r0Rear = TMath::Sqrt(r0Rear)/TMath::Abs(fPosition)/tan3;
1403   r0Norm = r0Rear;
1404   Double_t p0, cos25, cos60;
1405   if (!iflag) goto vertex;
1406
1407   for (Int_t i=4; i>=0; i--) {
1408     ParPropagation(zPos[i]);
1409     WeightPropagation(zPos[i], kFALSE);
1410     dZ = TMath::Abs (fPositionNew-fPosition);
1411     if (r0Norm > 1) x0 = x01[i];
1412     else x0 = x02[i];
1413     MSLine(dZ,x0); // multiple scattering in the medium (linear approximation)
1414     fPosition = fPositionNew;
1415     *fTrackPar = *fTrackParNew;
1416     r0Norm = (*fTrackPar)(0,0)*(*fTrackPar)(0,0) +
1417              (*fTrackPar)(1,0)*(*fTrackPar)(1,0);
1418     r0Norm = TMath::Sqrt(r0Norm)/TMath::Abs(fPosition)/tan3;
1419   }
1420   // Correct momentum for energy losses
1421   pTotal = 1/TMath::Abs((*fTrackPar)(4,0));
1422   p0 = pTotal;
1423   cos25 = TMath::Cos(2.5/180*TMath::Pi());
1424   cos60 = TMath::Cos(6.0/180*TMath::Pi());
1425   for (Int_t j=0; j<1; j++) {
1426     /*
1427     if (r0Rear > 1) {
1428       if (p0 < 20) {
1429         deltaP = 2.164 + 0.145e-1*p0 - 0.417e-3*p0*p0;
1430       } else {
1431         deltaP = 2.275 + 0.102e-2*p0 - 0.674e-6*p0*p0;
1432       }
1433     } else {
1434       if (p0 < 20) {
1435         deltaP = 2.581 + 0.188e-1*p0 - 0.398e-3*p0*p0;
1436       } else {
1437         deltaP = 2.727 + 0.356e-2*p0 + 0.242e-5*p0*p0;
1438       }
1439     }
1440     */
1441     /*
1442     if (r0Rear < 1) {
1443       //W
1444       if (p0<15) {
1445         deltaP = 2.737 + 0.0494*p0 - 0.001123*p0*p0;
1446       } else {
1447         deltaP = 3.0643 + 0.01346*p0;
1448       }
1449       deltaP *= 0.95;
1450     } else {
1451       //Pb
1452       if (p0<15) {
1453         deltaP  = 2.1380 + 0.0351*p0 - 0.000853*p0*p0;
1454       } else {
1455         deltaP = 2.407 + 0.00702*p0;
1456       }
1457       deltaP *= 0.95;
1458     }
1459     */
1460     /*
1461     if (r0Rear < 1) {
1462       //W
1463       if (p0<18) {
1464         deltaP = 2.439 + 0.806e-1*p0 - 0.500e-2*p0*p0 + 0.106e-3*p0*p0*p0;
1465       } else {
1466         deltaP = 2.767 + 0.742e-2*p0 - 0.196e-4*p0*p0 + 0.403e-7*p0*p0*p0;
1467       }
1468       //deltaP += 0.2;
1469       deltaP *= cos25;
1470     } else {
1471       //Pb
1472       if (p0<18) {
1473         deltaP  = 2.209 + 0.800e-2*p0;
1474       } else {
1475         deltaP = 2.285 + 0.141e-2*p0 - 0.446e-6*p0*p0;
1476       }
1477       //deltaP += 0.2;
1478       deltaP *= cos60;
1479     }
1480     deltaP *= 1.1;
1481     */
1482     //*
1483     if (r0Rear  < 1) {
1484       if (p0 < 20) {
1485         deltaP = 2.5938 + 0.0570 * p0 - 0.001151 * p0 * p0;
1486       } else {
1487         deltaP = 3.0714 + 0.011767 * p0;
1488       }
1489       deltaP *= 0.75;
1490     } else {
1491       if (p0 < 20) {
1492         deltaP  = 2.1207 + 0.05478 * p0 - 0.00145079 * p0 * p0;
1493       } else {
1494         deltaP = 2.6069 + 0.0051705 * p0;
1495       }
1496       deltaP *= 0.9;
1497     }
1498     //*/
1499
1500     p0 = pTotal + deltaP/TMath::Abs(TMath::Cos((*fTrackPar)(2,0))/TMath::Cos((*fTrackPar)(3,0)));
1501   }
1502   (*fTrackPar)(4,0) = 1/p0*TMath::Sign((Double_t)1.,(*fTrackPar)(4,0));
1503
1504   // Go to the vertex
1505 vertex:
1506   ParPropagation((Double_t)0.);
1507   WeightPropagation((Double_t)0., kFALSE);
1508   fPosition = fPositionNew;
1509   //*fTrackPar = *fTrackParNew;
1510   // Add vertex as a hit
1511   TMatrixD pointWeight(fgkSize,fgkSize);
1512   TMatrixD point(fgkSize,1);
1513   TMatrixD trackParTmp = point;
1514   point(0,0) = 0; // vertex coordinate - should be taken somewhere
1515   point(1,0) = 0; // vertex coordinate - should be taken somewhere
1516   pointWeight(0,0) = 1/1.e-3/1.e-3; // 10 um error
1517   pointWeight(1,1) = 1/1.e-3/1.e-3; // 10 um error
1518   TryPoint(point,pointWeight,trackParTmp,dChi2);
1519   *fTrackPar = trackParTmp;
1520   *fWeight += pointWeight;
1521   fChi2 += dChi2; // Chi2
1522   if (fgDebug < 0) return; // no output
1523
1524   cout << pOld << " " << 1/(*fTrackPar)(4,0) << " " << dChi2 << " " << fChi2 << " " << fNmbTrackHits << endl;
1525   for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
1526     hit =  (AliMUONHitForRec*) ((*fTrackHits)[i1]);
1527     printf ("%5d", hit->GetChamberNumber());
1528   }
1529   cout << endl;
1530   if (fgDebug > 0) {
1531     for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
1532       hit =  (AliMUONHitForRec*) ((*fTrackHits)[i1]);
1533       //cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetHitNumber() << " ";
1534       //cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetZ() << " ";
1535       printf ("%5d", fgHitForRec->IndexOf(hit));
1536     }
1537     cout << endl;
1538   }
1539
1540   // from raw clusters
1541   for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
1542     hit =  (AliMUONHitForRec*) ((*fTrackHits)[i1]);
1543     if (hit->GetHitNumber() < 0) { // combined cluster / track finder
1544       Int_t index = -hit->GetHitNumber() / 100000;
1545       Int_t iPos = -hit->GetHitNumber() - index * 100000;
1546       clus = (AliMUONRawCluster*) fgCombi->DetElem(index-1)->RawClusters()->UncheckedAt(iPos);
1547     } else {
1548       rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
1549       clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
1550     }
1551     printf ("%5d", clus->GetTrack(1)%10000000);
1552
1553     cout << endl;
1554     for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
1555       hit =  (AliMUONHitForRec*) ((*fTrackHits)[i1]);
1556       if (hit->GetHitNumber() < 0) { // combined cluster / track finder
1557         Int_t index = -hit->GetHitNumber() / 100000;
1558         Int_t iPos = -hit->GetHitNumber() - index * 100000;
1559         clus = (AliMUONRawCluster*) fgCombi->DetElem(index-1)->RawClusters()->UncheckedAt(iPos);
1560       } else {
1561         rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
1562         clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
1563       }
1564       if (clus->GetTrack(2) != -1) printf ("%5d", clus->GetTrack(2)%10000000);
1565       else printf ("%5s", "    ");
1566     }
1567   }
1568   cout << endl;
1569   for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
1570     //cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetHitNumber() << " ";
1571     cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetZ() << " ";
1572     //cout << fgHitForRec->IndexOf(((AliMUONHitForRec*)((*fTrackHits)[i1]))) << " ";
1573   }
1574   cout << endl;
1575   for (Int_t i1=0; i1<fNmbTrackHits; i1++) printf("%8.4f", (*fChi2Smooth)[i1]);
1576   cout << endl;
1577   cout << "---------------------------------------------------" << endl;
1578
1579   // Get covariance matrix
1580   /* Not needed - covariance matrix is not interesting to anybody
1581   *fCovariance = *fWeight;
1582   if (fCovariance->Determinant() != 0) {
1583     Int_t ifail;
1584     mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1585   } else {
1586     AliWarning(" Determinant fCovariance=0:" );
1587   }
1588   */
1589 }
1590
1591   //__________________________________________________________________________
1592 void AliMUONTrackK::MSLine(Double_t dZ, Double_t x0)
1593 {
1594 /// Adds multiple scattering in a thick layer for linear propagation
1595
1596   Double_t cosAlph = TMath::Cos((*fTrackPar)(2,0));
1597   Double_t tanAlph = TMath::Tan((*fTrackPar)(2,0));
1598   Double_t cosBeta = TMath::Cos((*fTrackPar)(3,0));
1599   Double_t sinBeta;
1600   sinBeta = TMath::Sin((*fTrackPar)(3,0));
1601   Double_t tanBeta = TMath::Tan((*fTrackPar)(3,0));
1602   Double_t momentum = 1/(*fTrackPar)(4,0);
1603   Double_t velo = 1; // relativistic velocity
1604   Double_t step = TMath::Abs(dZ/cosAlph/cosBeta); // step length
1605
1606   // Projected scattering angle
1607   Double_t theta0 = 0.0136/velo/momentum/TMath::Sqrt(x0)*(1+0.038*TMath::Log(step/x0));
1608   Double_t theta02 = theta0*theta0;
1609   Double_t dl2 = step*step/2*theta02;
1610   Double_t dl3 = dl2*step*2/3;
1611
1612   //Derivatives
1613   Double_t dYdT = 1/cosAlph;
1614   Double_t dYdB = 0; //(*fTrackPar)(2,0)*sinBeta/cosAlph;
1615   Double_t dXdT = tanAlph*tanBeta;
1616   //Double_t dXdB = (1+(*fTrackPar)(2,0)*tanAlph*sinBeta*sinBeta)/cosBeta;
1617   Double_t dXdB = 1/cosBeta;
1618   Double_t dAdT = 1/cosBeta;
1619   Double_t dAdB = 0; //(*fTrackPar)(2,0)*tanBeta;
1620
1621   // Get covariance matrix
1622   *fCovariance = *fWeight;
1623   if (fCovariance->Determinant() != 0) {
1624     //   fCovariance->Invert();
1625     Int_t ifail;
1626     mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1627   } else {
1628     AliWarning(" Determinant fCovariance=0:" );
1629   }
1630
1631   (*fCovariance)(0,0) += dl3*(dYdT*dYdT+dYdB*dYdB); // <yy>
1632   (*fCovariance)(1,1) += dl3*(dXdT*dXdT+dXdB*dXdB); // <xx>
1633   (*fCovariance)(2,2) += theta02*step*(dAdT*dAdT+dAdB*dAdB); // <aa>
1634   (*fCovariance)(3,3) += theta02*step; // <bb>
1635
1636   (*fCovariance)(0,1) += dl3*(dYdT*dXdT+dYdB*dXdB); // <yx>
1637   (*fCovariance)(1,0) = (*fCovariance)(0,1);
1638
1639   (*fCovariance)(0,2) += dl2*(dYdT*dAdT+dYdB*dAdB); // <ya>
1640   (*fCovariance)(2,0) = (*fCovariance)(0,2);
1641
1642   (*fCovariance)(0,3) += dl2*dYdB; // <yb>
1643   (*fCovariance)(3,0) = (*fCovariance)(0,3);
1644
1645   (*fCovariance)(1,2) += dl2*(-dXdT*dAdT+dXdB*dAdB); // <xa>
1646   (*fCovariance)(2,1) = (*fCovariance)(1,2);
1647
1648   (*fCovariance)(1,3) += dl2*dXdB; // <xb>
1649   (*fCovariance)(3,1) = (*fCovariance)(1,3);
1650
1651   (*fCovariance)(2,3) += theta02*step*dAdB; // <ab>
1652   (*fCovariance)(3,2) = (*fCovariance)(2,3);
1653
1654   // Get weight matrix
1655   *fWeight = *fCovariance;
1656   if (fWeight->Determinant() != 0) {
1657     Int_t ifail;
1658     mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1659   } else {
1660     AliWarning(" Determinant fWeight=0:");
1661   }
1662 }
1663
1664   //__________________________________________________________________________
1665 Bool_t AliMUONTrackK::Recover(void)
1666 {
1667 /// Adds new failed track(s) which can be tried to be recovered
1668   Int_t nRecTracks;
1669   TClonesArray *trackPtr;
1670   AliMUONTrackK *trackK;
1671
1672   if (fgDebug > 0) cout << " ******** Enter Recover " << endl;
1673   trackPtr = fgTrackReconstructor->GetRecTracksPtr();
1674
1675   // Remove hit with the highest chi2
1676   Double_t chi2 = 0;
1677   if (fgDebug > 0) {
1678     for (Int_t i=0; i<fNmbTrackHits; i++) {
1679       chi2 = fChi2Smooth ? (*fChi2Smooth)[i] : (*fChi2Array)[i];
1680       printf("%10.4f", chi2);
1681     }
1682     printf("\n");
1683     for (Int_t i=0; i<fNmbTrackHits; i++) {
1684       printf("%10d", ((AliMUONHitForRec*)fTrackHits->UncheckedAt(i))->GetChamberNumber());
1685     }
1686     printf("\n");
1687   }
1688   Double_t chi2max = 0;
1689   Int_t imax = 0;
1690   for (Int_t i=0; i<fNmbTrackHits; i++) {
1691     chi2 = fChi2Smooth ? (*fChi2Smooth)[i] : (*fChi2Array)[i];
1692     if (chi2 < chi2max) continue;
1693     chi2max = chi2;
1694     imax = i;
1695   }
1696   //if (chi2max < 10) return kFALSE; // !!!
1697   //if (chi2max < 25) imax = fNmbTrackHits - 1;
1698   if (chi2max < 15) imax = fNmbTrackHits - 1; // discard the last point
1699   // Check if the outlier is not from the seed segment
1700   AliMUONHitForRec *skipHit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(imax);
1701   if (skipHit == fStartSegment->GetHitForRec1() || skipHit == fStartSegment->GetHitForRec2()) {
1702     //DropBranches(fStartSegment); // drop all tracks with the same seed segment
1703     return kFALSE; // to be changed probably
1704   }
1705
1706   // Make a copy of track hit collection
1707   TObjArray *hits = new TObjArray(*fTrackHits);
1708   Int_t imax0;
1709   imax0 = imax;
1710
1711   // Hits after the found one will be removed
1712   if (GetStation0() == 3 && skipHit->GetChamberNumber() >= 7) {
1713     SortHits(1, fTrackHits); // unsort hits
1714     imax = fTrackHits->IndexOf(skipHit);
1715   }
1716   Int_t nTrackHits = fNmbTrackHits;
1717   for (Int_t i=imax+1; i<nTrackHits; i++) {
1718     AliMUONHitForRec *hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(i);
1719     fTrackHits->Remove(hit);
1720     hit->SetNTrackHits(hit->GetNTrackHits()-1); // unmark hit
1721     fNmbTrackHits--;
1722   }
1723
1724   // Check if the track candidate doesn't exist yet
1725   if (ExistDouble()) { delete hits; return kFALSE; }
1726
1727   //DropBranches(imax0, hits); // drop branches downstream the discarded hit
1728   delete hits;
1729
1730   nRecTracks = fgTrackReconstructor->GetNRecTracks();
1731   skipHit = (AliMUONHitForRec*) ((*fTrackHits)[fNmbTrackHits-1]);
1732   // Remove all saved steps and smoother matrices after the skipped hit
1733   RemoveMatrices(skipHit->GetZ());
1734
1735   //AZ(z->-z) if (skipHit->GetZ() > fStartSegment->GetHitForRec2()->GetZ() || !fNSteps) {
1736   if (TMath::Abs(skipHit->GetZ()) > TMath::Abs(fStartSegment->GetHitForRec2()->GetZ()) || !fNSteps) {
1737     // Propagation toward high Z or skipped hit next to segment -
1738     // start track from segment
1739     trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(fStartSegment);
1740     fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
1741     trackK->fRecover = 1;
1742     trackK->fSkipHit = skipHit;
1743     trackK->fNmbTrackHits = fNmbTrackHits;
1744     delete trackK->fTrackHits; // not efficient ?
1745     trackK->fTrackHits = new TObjArray(*fTrackHits);
1746     if (fgDebug > 0) cout << nRecTracks << " " << trackK->fRecover << endl;
1747     return kTRUE;
1748   }
1749
1750   trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(NULL, NULL);
1751   *trackK = *this;
1752   fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
1753   //AZ(z->-z) trackK->fTrackDir = -1;
1754   trackK->fTrackDir = 1;
1755   trackK->fRecover = 1;
1756   trackK->fSkipHit = skipHit;
1757   Int_t iD = trackK->fParFilter->Last()->GetUniqueID();
1758   if (iD > 1) {
1759     trackK->fParFilter->Last()->SetUniqueID(iD-1);
1760     CreateMatrix(trackK->fParFilter);
1761   }
1762   *((TMatrixD*)trackK->fParFilter->Last()) = *((TMatrixD*)fParExtrap->Last());
1763   trackK->fParFilter->Last()->SetUniqueID(1);
1764   *(trackK->fTrackPar) = *((TMatrixD*)trackK->fParFilter->Last());
1765   iD = trackK->fCovFilter->Last()->GetUniqueID();
1766   if (iD > 1) {
1767     trackK->fCovFilter->Last()->SetUniqueID(iD-1);
1768     CreateMatrix(trackK->fCovFilter);
1769   }
1770   *((TMatrixD*)trackK->fCovFilter->Last()) = *((TMatrixD*)fCovExtrap->Last());
1771   trackK->fCovFilter->Last()->SetUniqueID(1);
1772   *(trackK->fWeight) = *((TMatrixD*)trackK->fCovFilter->Last());
1773   if (trackK->fWeight->Determinant() != 0) {
1774     Int_t ifail;
1775     mnvertLocal(&((*(trackK->fWeight))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1776   } else {
1777     AliWarning(" Determinant fWeight=0:");
1778   }
1779   trackK->fPosition = trackK->fPositionNew = (*fSteps)[fNSteps-1];
1780   trackK->fChi2 = 0;
1781   for (Int_t i=0; i<fNmbTrackHits-1; i++) trackK->fChi2 += (*fChi2Array)[i];
1782   if (fgDebug > 0) cout << nRecTracks << " " << trackK->fRecover << endl;
1783   return kTRUE;
1784 }
1785
1786   //__________________________________________________________________________
1787 void AliMUONTrackK::AddMatrices(AliMUONTrackK *trackK, Double_t dChi2, AliMUONHitForRec *hitAdd)
1788 {
1789 /// Adds matrices for the smoother and keep Chi2 for the point
1790 /// Track parameters
1791   //trackK->fParFilter->Last()->Print();
1792   Int_t iD = trackK->fParFilter->Last()->GetUniqueID();
1793   if (iD > 1) {
1794     trackK->fParFilter->Last()->SetUniqueID(iD-1);
1795     CreateMatrix(trackK->fParFilter);
1796     iD = 1;
1797   }
1798   *((TMatrixD*)(trackK->fParFilter->Last())) = *(trackK->fTrackPar);
1799   trackK->fParFilter->Last()->SetUniqueID(iD);
1800   if (fgDebug > 1) {
1801     cout << " Add matrices" << " " << fPosition << " " << dChi2 << " " << fgHitForRec->IndexOf(hitAdd) << endl;
1802     //trackK->fTrackPar->Print();
1803     //trackK->fTrackParNew->Print();
1804     trackK->fParFilter->Last()->Print();
1805     cout << " Add matrices" << endl;
1806   }
1807   // Covariance
1808   *(trackK->fCovariance) = *(trackK->fWeight);
1809   if (trackK->fCovariance->Determinant() != 0) {
1810     Int_t ifail;
1811     mnvertLocal(&((*(trackK->fCovariance))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1812   } else {
1813     AliWarning(" Determinant fCovariance=0:");
1814   }
1815   iD = trackK->fCovFilter->Last()->GetUniqueID();
1816   if (iD > 1) {
1817     trackK->fCovFilter->Last()->SetUniqueID(iD-1);
1818     CreateMatrix(trackK->fCovFilter);
1819     iD = 1;
1820   }
1821   *((TMatrixD*)(trackK->fCovFilter->Last())) = *(trackK->fCovariance);
1822   trackK->fCovFilter->Last()->SetUniqueID(iD);
1823
1824   // Save Chi2-increment for point
1825   Int_t indx = trackK->fTrackHits->IndexOf(hitAdd);
1826   if (indx < 0) indx = fNmbTrackHits;
1827   if (trackK->fChi2Array->fN <= indx) trackK->fChi2Array->Set(indx+10);
1828   trackK->fChi2Array->AddAt(dChi2,indx);
1829 }
1830
1831   //__________________________________________________________________________
1832 void AliMUONTrackK::CreateMatrix(TObjArray *objArray) const
1833 {
1834 /// Create new matrix and add it to TObjArray
1835
1836   TMatrixD *matrix = (TMatrixD*) objArray->First();
1837   TMatrixD *tmp = new TMatrixD(*matrix);
1838   objArray->AddAtAndExpand(tmp,objArray->GetLast());
1839   //cout << " !!! " << tmp << " " << tmp->GetUniqueID() << " " << tmp->GetNoElements() << endl;
1840 }
1841
1842   //__________________________________________________________________________
1843 void AliMUONTrackK::RemoveMatrices(Double_t zEnd)
1844 {
1845 /// Remove matrices (and saved steps) in the smoother part with abs(z) < abs(zEnd)
1846
1847   for (Int_t i=fNSteps-1; i>=0; i--) {
1848     if (fgDebug > 1) printf("%10.4f %10.4f \n", (*fSteps)[i], zEnd);
1849     if (TMath::Abs((*fSteps)[i]) > TMath::Abs(zEnd) - 0.01) break;
1850     RemoveMatrices(this);
1851   } // for (Int_t i=fNSteps-1;
1852 }
1853
1854   //__________________________________________________________________________
1855 void AliMUONTrackK::RemoveMatrices(AliMUONTrackK* trackK)
1856 {
1857 /// Remove last saved matrices and steps in the smoother part
1858
1859   trackK->fNSteps--;
1860   Int_t i = trackK->fNSteps;
1861
1862   Int_t id = 0;
1863   // Delete only matrices not shared by several tracks
1864   id = trackK->fParExtrap->Last()->GetUniqueID();
1865   if (id > 1) {
1866     trackK->fParExtrap->Last()->SetUniqueID(id-1);
1867     trackK->fParExtrap->RemoveAt(i);
1868   }
1869   else ((TMatrixD*)(trackK->fParExtrap->RemoveAt(i)))->Delete();
1870   id = fParFilter->Last()->GetUniqueID();
1871   if (id > 1) {
1872     trackK->fParFilter->Last()->SetUniqueID(id-1);
1873     trackK->fParFilter->RemoveAt(i);
1874   }
1875   else ((TMatrixD*)(trackK->fParFilter->RemoveAt(i)))->Delete();
1876   id = trackK->fCovExtrap->Last()->GetUniqueID();
1877   if (id > 1) {
1878     trackK->fCovExtrap->Last()->SetUniqueID(id-1);
1879     trackK->fCovExtrap->RemoveAt(i);
1880   }
1881   else ((TMatrixD*)(trackK->fCovExtrap->RemoveAt(i)))->Delete();
1882   id = trackK->fCovFilter->Last()->GetUniqueID();
1883   if (id > 1) {
1884     trackK->fCovFilter->Last()->SetUniqueID(id-1);
1885     trackK->fCovFilter->RemoveAt(i);
1886   }
1887   else ((TMatrixD*)(trackK->fCovFilter->RemoveAt(i)))->Delete();
1888   id = trackK->fJacob->Last()->GetUniqueID();
1889   if (id > 1) {
1890     trackK->fJacob->Last()->SetUniqueID(id-1);
1891     trackK->fJacob->RemoveAt(i);
1892   }
1893   else ((TMatrixD*)(trackK->fJacob->RemoveAt(i)))->Delete();
1894 }
1895
1896   //__________________________________________________________________________
1897 Bool_t AliMUONTrackK::Smooth(void)
1898 {
1899 /// Apply smoother
1900   Int_t ihit = fNmbTrackHits - 1;
1901   AliMUONHitForRec *hit = (AliMUONHitForRec*) (*fTrackHits)[ihit];
1902   fChi2Smooth = new TArrayD(fNmbTrackHits);
1903   fChi2Smooth->Reset(-1);
1904   fChi2 = 0;
1905   fParSmooth = new TObjArray(15);
1906   fParSmooth->Clear();
1907
1908   if (fgDebug > 0) {
1909     cout << " ******** Enter Smooth " << endl;
1910     cout << (*fSteps)[fNSteps-1] << " " << fPosition << " " << fPositionNew << endl;
1911     /*
1912     for (Int_t i=fNSteps-1; i>=0; i--) cout << (*fSteps)[i] << " ";
1913     cout << endl;
1914     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();}
1915     */
1916     for (Int_t i=ihit; i>=0; i--) cout << ((AliMUONHitForRec*)(*fTrackHits)[i])->GetZ() << " ";
1917     cout << endl;
1918   }
1919
1920   // Find last point corresponding to the last hit
1921   Int_t iLast = fNSteps - 1;
1922   for ( ; iLast>=0; iLast--) {
1923     //AZ(z->-z) if ((*fSteps)[iLast] + 0.01 > ((AliMUONHitForRec*)(*fTrackHits)[fNmbTrackHits-1])->GetZ()) break;
1924     if (TMath::Abs((*fSteps)[iLast]) + 0.01 > TMath::Abs(((AliMUONHitForRec*)(*fTrackHits)[fNmbTrackHits-1])->GetZ())) break;
1925   }
1926
1927   TMatrixD parSmooth = *((TMatrixD*)(fParFilter->UncheckedAt(iLast))); // last filtered = first smoothed
1928   //parSmooth.Dump();
1929   TMatrixD covSmooth = *((TMatrixD*)(fCovFilter->UncheckedAt(iLast))); // last filtered = first smoothed
1930   TMatrixD tmp=covSmooth, weight=covSmooth, gain = covSmooth;
1931   TMatrixD tmpPar = *fTrackPar;
1932   //parSmooth.Print(); ((TMatrixD*)(fParExtrap->Last()))->Print();
1933
1934   Bool_t found;
1935   Double_t chi2max = 0;
1936   for (Int_t i=iLast+1; i>0; i--) {
1937     if (i == iLast + 1) goto L33; // temporary fix
1938
1939     // Smoother gain matrix
1940     weight = *((TMatrixD*)(fCovExtrap->UncheckedAt(i)));
1941     if (weight.Determinant() != 0) {
1942       Int_t ifail;
1943       mnvertLocal(&(weight(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1944     } else {
1945       AliWarning(" Determinant weight=0:");
1946     }
1947     // Fk'Wkk+1
1948     tmp = TMatrixD(*((TMatrixD*)(fJacob->UncheckedAt(i))),TMatrixD::kTransposeMult,weight);
1949     gain = TMatrixD(*((TMatrixD*)(fCovFilter->UncheckedAt(i-1))),TMatrixD::kMult,tmp);
1950     //cout << (*fSteps)[i-1] << " " << (*fSteps)[i] << endl; gain.Print();
1951
1952     // Smoothed parameter vector
1953     //((TMatrixD*)(fParExtrap->UncheckedAt(i)))->Dump();
1954     parSmooth -= *((TMatrixD*)(fParExtrap->UncheckedAt(i)));
1955     tmpPar = *((TMatrixD*)(fParFilter->UncheckedAt(i-1)));
1956     tmpPar += TMatrixD(gain,TMatrixD::kMult,parSmooth);
1957     parSmooth = tmpPar;
1958
1959     // Smoothed covariance
1960     covSmooth -= *((TMatrixD*)(fCovExtrap->UncheckedAt(i)));
1961     weight = TMatrixD(TMatrixD::kTransposed,gain);
1962     tmp = TMatrixD(covSmooth,TMatrixD::kMult,weight);
1963     covSmooth = *((TMatrixD*)(fCovFilter->UncheckedAt(i-1)));
1964     covSmooth += TMatrixD(gain,TMatrixD::kMult,tmp);
1965
1966     // Check if there was a measurement at given z
1967     found = kFALSE;
1968     for ( ; ihit>=0; ihit--) {
1969       hit = (AliMUONHitForRec*) (*fTrackHits)[ihit];
1970       if (TMath::Abs(hit->GetZ()-(*fSteps)[i-1]) < 0.1) { found = kTRUE; break; }
1971       //AZ(z->-z) else if (ihit < fNmbTrackHits-1 && hit->GetZ() > (*fSteps)[i-1]) { ihit++; break; }
1972       else if (ihit < fNmbTrackHits-1 && TMath::Abs(hit->GetZ()) > TMath::Abs((*fSteps)[i-1])) { ihit++; break; }
1973     }
1974     if (!found) continue; // no measurement - skip the rest
1975     else if (fgDebug > 1) cout << i << " " << ihit << " " << hit->GetZ() << endl;
1976     if (ihit == 0) continue; // the first hit - skip the rest
1977
1978 L33:
1979     // Smoothed residuals
1980     tmpPar = 0;
1981     tmpPar(0,0) = hit->GetBendingCoor() - parSmooth(0,0); // bending coordinate
1982     tmpPar(1,0) = hit->GetNonBendingCoor() - parSmooth(1,0); // non-bending coordinate
1983     if (fgDebug > 1) {
1984       cout << hit->GetBendingCoor() << " " << parSmooth(0,0) << " " << tmpPar(0,0) << endl;
1985       cout << hit->GetNonBendingCoor() << " " << parSmooth(1,0) << " " << tmpPar(1,0) << endl;
1986     }
1987     // Cov. matrix of smoothed residuals
1988     tmp = 0;
1989     tmp(0,0) = hit->GetBendingReso2() - covSmooth(0,0);
1990     tmp(1,1) = hit->GetNonBendingReso2() - covSmooth(1,1);
1991     tmp(0,1) = tmp(1,0) = -covSmooth(0,1);
1992     tmp(2,2) = tmp(3,3) = tmp(4,4) = 1;
1993
1994     // Calculate Chi2 of smoothed residuals
1995     if (tmp.Determinant() != 0) {
1996       Int_t ifail;
1997       mnvertLocal(&(tmp(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1998     } else {
1999       AliWarning(" Determinant tmp=0:");
2000     }
2001     TMatrixD vector(tmp,TMatrixD::kMult,tmpPar);
2002     TMatrixD chi2(tmpPar,TMatrixD::kTransposeMult,vector);
2003     if (fgDebug > 1) chi2.Print();
2004     (*fChi2Smooth)[ihit] = chi2(0,0);
2005     if (chi2(0,0) > chi2max) chi2max = chi2(0,0);
2006     fChi2 += chi2(0,0);
2007     if (chi2(0,0) < 0) {
2008       //chi2.Print();
2009       AliError(Form(" *** chi2 < 0: %d %d ", i, iLast));
2010     }
2011     // Save smoothed parameters
2012     TMatrixD *par = new TMatrixD(parSmooth);
2013     fParSmooth->AddAtAndExpand(par, ihit);
2014
2015   } // for (Int_t i=iLast+1;
2016
2017   //if (chi2max > 16) {
2018   //if (chi2max > 25) {
2019   //if (chi2max > 50) {
2020   //if (chi2max > 100) {
2021   if (chi2max > fgkChi2max) {
2022     //if (Recover()) DropBranches();
2023     //Recover();
2024     Outlier();
2025     return kFALSE;
2026   }
2027   return kTRUE;
2028 }
2029
2030   //__________________________________________________________________________
2031 void AliMUONTrackK::Outlier()
2032 {
2033 /// Adds new track with removed hit having the highest chi2
2034
2035   if (fgDebug > 0) {
2036     cout << " ******** Enter Outlier " << endl;
2037     for (Int_t i=0; i<fNmbTrackHits; i++) printf("%10.4f", (*fChi2Smooth)[i]);
2038     printf("\n");
2039     for (Int_t i=0; i<fNmbTrackHits; i++) {
2040       printf("%10d", ((AliMUONHitForRec*)fTrackHits->UncheckedAt(i))->GetChamberNumber());
2041     }
2042     printf("\n");
2043   }
2044
2045   Double_t chi2max = 0;
2046   Int_t imax = 0;
2047   for (Int_t i=0; i<fNmbTrackHits; i++) {
2048     if ((*fChi2Smooth)[i] < chi2max) continue;
2049     chi2max = (*fChi2Smooth)[i];
2050     imax = i;
2051   }
2052   // Check if the outlier is not from the seed segment
2053   AliMUONHitForRec *hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(imax);
2054   if (hit == fStartSegment->GetHitForRec1() || hit == fStartSegment->GetHitForRec2()) return; // to be changed probably
2055
2056   // Check for missing station
2057   Int_t ok = 1;
2058   if (imax == 0) {
2059     if (((AliMUONHitForRec*)fTrackHits->UncheckedAt(1))->GetChamberNumber() < 8) ok--;
2060   } else if (imax == fNmbTrackHits-1) {
2061     if (((AliMUONHitForRec*)fTrackHits->UncheckedAt(fNmbTrackHits-2))->GetChamberNumber() > 1) ok--;
2062   }
2063   else if (((AliMUONHitForRec*)fTrackHits->UncheckedAt(imax-1))->GetChamberNumber()/2 - ((AliMUONHitForRec*)fTrackHits->UncheckedAt(imax+1))->GetChamberNumber()/2 > 1) ok--;
2064   if (!ok) { Recover(); return; } // try to recover track
2065   //AZ if (!ok) { if (fgDebug >= 0) cout << imax << endl; DropBranches(imax, 0); return; }
2066
2067   // Remove saved steps and smoother matrices after the outlier
2068   RemoveMatrices(hit->GetZ());
2069
2070   // Check for possible double track candidates
2071   //if (ExistDouble(hit)) return;
2072
2073   TClonesArray *trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2074   Int_t nRecTracks = fgTrackReconstructor->GetNRecTracks();
2075
2076   AliMUONTrackK *trackK = 0;
2077   if (!fNSteps || GetStation0() == 3 && hit->GetChamberNumber() > 7) {
2078     // start track from segment
2079     trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(fStartSegment);
2080     fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
2081     trackK->fRecover = 2;
2082     trackK->fSkipHit = hit;
2083     trackK->fNmbTrackHits = fNmbTrackHits;
2084
2085     hit = (AliMUONHitForRec*) trackK->fTrackHits->UncheckedAt(0);
2086     hit->SetNTrackHits(hit->GetNTrackHits()-1);
2087     hit = (AliMUONHitForRec*) trackK->fTrackHits->UncheckedAt(1);
2088     hit->SetNTrackHits(hit->GetNTrackHits()-1);
2089     delete trackK->fTrackHits; // not efficient ?
2090     trackK->fTrackHits = new TObjArray(*fTrackHits);
2091     for (Int_t i = 0; i < fNmbTrackHits; i++) {
2092       hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(i);
2093       hit->SetNTrackHits(hit->GetNTrackHits()+1);
2094     }
2095
2096     if (GetStation0() == 3) trackK->SortHits(1, trackK->fTrackHits);
2097     if (fgDebug > 0) cout << nRecTracks << " Outlier" << endl;
2098     return;
2099   }
2100   trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(NULL, NULL);
2101   *trackK = *this;
2102   fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
2103   trackK->fTrackDir = 1;
2104   trackK->fRecover = 2;
2105   trackK->fSkipHit = hit;
2106   Int_t iD = trackK->fParFilter->Last()->GetUniqueID();
2107   if (iD > 1) {
2108     trackK->fParFilter->Last()->SetUniqueID(iD-1);
2109     CreateMatrix(trackK->fParFilter);
2110   }
2111   *((TMatrixD*)trackK->fParFilter->Last()) = *((TMatrixD*)fParExtrap->Last());
2112   trackK->fParFilter->Last()->SetUniqueID(1);
2113   *(trackK->fTrackPar) = *((TMatrixD*)trackK->fParFilter->Last());
2114   iD = trackK->fCovFilter->Last()->GetUniqueID();
2115   if (iD > 1) {
2116     trackK->fCovFilter->Last()->SetUniqueID(iD-1);
2117     CreateMatrix(trackK->fCovFilter);
2118   }
2119   *((TMatrixD*)trackK->fCovFilter->Last()) = *((TMatrixD*)fCovExtrap->Last());
2120   trackK->fCovFilter->Last()->SetUniqueID(1);
2121   *(trackK->fWeight) = *((TMatrixD*)trackK->fCovFilter->Last());
2122   if (trackK->fWeight->Determinant() != 0) {
2123     Int_t ifail;
2124     mnvertLocal(&((*(trackK->fWeight))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
2125   } else {
2126     AliWarning(" Determinant fWeight=0:");
2127   }
2128   trackK->fPosition = trackK->fPositionNew = (*fSteps)[fNSteps-1];
2129   trackK->fChi2 = 0;
2130   for (Int_t i=0; i<fNmbTrackHits-1; i++) trackK->fChi2 += (*fChi2Array)[i];
2131   if (fgDebug > 0) cout << nRecTracks << " Outlier " << trackK->fChi2 << endl;
2132 }
2133
2134   //__________________________________________________________________________
2135 Double_t AliMUONTrackK::GetChi2PerPoint(Int_t iPoint) const
2136 {
2137 /// Return Chi2 at point
2138   return fChi2Smooth ? (*fChi2Smooth)[iPoint] : (*fChi2Array)[iPoint];
2139   //return 0.;
2140 }
2141
2142   //__________________________________________________________________________
2143 void AliMUONTrackK::Print(FILE *lun) const
2144 {
2145 /// Print out track information
2146
2147   Int_t flag = 1;
2148   AliMUONHitForRec *hit = 0;
2149     // from raw clusters
2150   AliMUONRawCluster *clus = 0;
2151   TClonesArray *rawclusters = 0;
2152   for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
2153     hit =  (AliMUONHitForRec*) ((*fTrackHits)[i1]);
2154     rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
2155     clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
2156     if (TMath::Abs(clus->GetTrack(1)-1) < 2) {
2157       if (clus->GetTrack(2)) flag = 2;
2158       continue;
2159     }
2160     if (clus->GetTrack(2) && TMath::Abs(clus->GetTrack(2)-1) < 2) {
2161       flag = 3;
2162       continue;
2163     }
2164     flag = 0;
2165     break;
2166
2167     Int_t sig[2]={1,1}, tid[2]={0};
2168     for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
2169       if (GetChi2PerPoint(i1) < -0.1) continue;
2170       hit =  (AliMUONHitForRec*) ((*fTrackHits)[i1]);
2171       rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
2172       clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
2173       for (Int_t j=0; j<2; j++) {
2174         tid[j] = clus->GetTrack(j+1) - 1;
2175         if (clus->GetTrack(j+1) < 0) { sig[j] = 0; tid[j] = 999; }
2176       }
2177       fprintf(lun,"%3d %3d %10.4f", gAlice->GetEvNumber(), hit->GetChamberNumber(), GetChi2PerPoint(i1));
2178       if (!(clus->GetTrack(2))) fprintf(lun, "%3d %3d", sig[0], tid[0]); // simple cluster
2179       else { // track overlap
2180         fprintf(lun, "%3d %3d", TMath::Max(sig[0],sig[1]), TMath::Min(tid[0],tid[1]));
2181         //if (tid[0] < 2) flag *= 2;
2182         //else if (tid[1] < 2) flag *= 3;
2183       }
2184       fprintf (lun, "%3d \n", flag);
2185     }
2186   }
2187 }
2188
2189   //__________________________________________________________________________
2190 void AliMUONTrackK::DropBranches(Int_t imax, TObjArray *hits)
2191 {
2192 /// Drop branches downstream of the skipped hit
2193   Int_t nRecTracks;
2194   TClonesArray *trackPtr;
2195   AliMUONTrackK *trackK;
2196
2197   trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2198   nRecTracks = fgTrackReconstructor->GetNRecTracks();
2199   Int_t icand = trackPtr->IndexOf(this);
2200   if (!hits) hits = fTrackHits;
2201
2202   // Check if the track candidate doesn't exist yet
2203   for (Int_t i=icand+1; i<nRecTracks; i++) {
2204     trackK = (AliMUONTrackK*) ((*trackPtr)[i]);
2205     if (trackK->fNmbTrackHits == 2 && trackK->GetRecover() == 0) continue;
2206     if (trackK->GetRecover() < 0) continue;
2207
2208     if (trackK->fNmbTrackHits >= imax + 1) {
2209       for (Int_t j=0; j<=imax; j++) {
2210         //if (j != fNmbTrackHits-1 && (*trackK->fTrackHits)[j] != (*fTrackHits)[j]) break;
2211         if ((*trackK->fTrackHits)[j] != (*hits)[j]) break;
2212         if (j == imax) {
2213           if (hits != fTrackHits) {
2214             // Drop all branches downstream the hit (for Recover)
2215             trackK->SetRecover(-1);
2216             if (fgDebug >= 0 )cout << " Recover " << i << endl;
2217             continue;
2218           }
2219           // Check if the removal of the hit breaks the track
2220           Int_t ok = 1;
2221           if (imax == 0) {
2222             if (((AliMUONHitForRec*)trackK->fTrackHits->UncheckedAt(1))->GetChamberNumber() < 8) ok--; }
2223           else if (imax == trackK->fNmbTrackHits-1) continue;
2224             // else if (imax == trackK->fNmbTrackHits-1) {
2225             //if (((AliMUONHitForRec*)trackK->fTrackHits->UncheckedAt(trackK->fNmbTrackHits-2))->GetChamberNumber() > 1) ok--;
2226             //}
2227           else if (((AliMUONHitForRec*)trackK->fTrackHits->UncheckedAt(imax-1))->GetChamberNumber()/2 - ((AliMUONHitForRec*)trackK->fTrackHits->UncheckedAt(imax+1))->GetChamberNumber()/2 > 1) ok--;
2228           if (!ok) trackK->SetRecover(-1);
2229         }
2230       } // for (Int_t j=0;
2231     }
2232   } // for (Int_t i=0;
2233 }
2234
2235   //__________________________________________________________________________
2236 void AliMUONTrackK::DropBranches(AliMUONSegment *segment)
2237 {
2238 /// Drop all candidates with the same seed segment
2239   Int_t nRecTracks;
2240   TClonesArray *trackPtr;
2241   AliMUONTrackK *trackK;
2242
2243   trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2244   nRecTracks = fgTrackReconstructor->GetNRecTracks();
2245   Int_t icand = trackPtr->IndexOf(this);
2246
2247   for (Int_t i=icand+1; i<nRecTracks; i++) {
2248     trackK = (AliMUONTrackK*) ((*trackPtr)[i]);
2249     if (trackK->fNmbTrackHits == 2 && trackK->GetRecover() == 0) continue;
2250     if (trackK->GetRecover() < 0) continue;
2251     if (trackK->fStartSegment == segment) trackK->SetRecover(-1);
2252   }
2253   if (fgDebug >= 0) cout << " Drop segment " << endl;
2254 }
2255
2256   //__________________________________________________________________________
2257 AliMUONHitForRec* AliMUONTrackK::GetHitLastOk(void)
2258 {
2259 /// Return the hit where track stopped
2260
2261   if (!fNSteps) return (AliMUONHitForRec*)((*fTrackHits)[1]);
2262   return fSkipHit;
2263 }
2264
2265   //__________________________________________________________________________
2266 Int_t AliMUONTrackK::GetStation0(void)
2267 {
2268 /// Return seed station number
2269   return fStartSegment->GetHitForRec1()->GetChamberNumber() / 2;
2270 }
2271
2272   //__________________________________________________________________________
2273 Bool_t AliMUONTrackK::ExistDouble(AliMUONHitForRec *hit)
2274 {
2275 /// Check if the track will make a double after outlier removal
2276
2277   TClonesArray *trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2278   Int_t nRecTracks = fgTrackReconstructor->GetNRecTracks();
2279   TObjArray *hitArray = new TObjArray(*fTrackHits);
2280   TObjArray *hitArray1 = new TObjArray(*hitArray);
2281   hitArray1->Remove(hit);
2282   hitArray1->Compress();
2283
2284   Bool_t same = kFALSE;
2285   for (Int_t i=0; i<nRecTracks; i++) {
2286     AliMUONTrackK *trackK = (AliMUONTrackK*) ((*trackPtr)[i]);
2287     if (trackK->fNmbTrackHits == 2 && trackK->GetRecover() == 0) continue;
2288     if (trackK == this) continue;
2289     if (trackK->fNmbTrackHits == fNmbTrackHits || trackK->fNmbTrackHits == fNmbTrackHits-1) {
2290       TObjArray *hits = new TObjArray(*trackK->fTrackHits);
2291       same = kTRUE;
2292       if (trackK->fNmbTrackHits == fNmbTrackHits) {
2293         for (Int_t j=0; j<fNmbTrackHits; j++) {
2294           if (hits->UncheckedAt(j) != hitArray->UncheckedAt(j)) { same = kFALSE; break; }
2295         }
2296         if (same) { delete hits; break; }
2297         if (trackK->fSkipHit) {
2298           TObjArray *hits1 = new TObjArray(*hits);
2299           if (hits1->Remove(trackK->fSkipHit) > 0) {
2300             hits1->Compress();
2301             same = kTRUE;
2302             for (Int_t j=0; j<fNmbTrackHits-1; j++) {
2303               if (hits1->UncheckedAt(j) != hitArray1->UncheckedAt(j)) { same = kFALSE; break; }
2304             }
2305             if (same) { delete hits1; break; }
2306           }
2307           delete hits1;
2308         }
2309       } else {
2310         // Check with removed outlier
2311         same = kTRUE;
2312         for (Int_t j=0; j<fNmbTrackHits-1; j++) {
2313           if (hits->UncheckedAt(j) != hitArray1->UncheckedAt(j)) { same = kFALSE; break; }
2314         }
2315         if (same) { delete hits; break; }
2316       }
2317       delete hits;
2318     }
2319   } // for (Int_t i=0; i<nRecTracks;
2320   delete hitArray; delete hitArray1;
2321   if (same && fgDebug >= 0) cout << " Same" << endl;
2322   return same;
2323 }
2324
2325   //__________________________________________________________________________
2326 Bool_t AliMUONTrackK::ExistDouble(void)
2327 {
2328 /// Check if the track will make a double after recovery
2329
2330   TClonesArray *trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2331   Int_t nRecTracks = fgTrackReconstructor->GetNRecTracks();
2332
2333   TObjArray *hitArray = new TObjArray(*fTrackHits);
2334   if (GetStation0() == 3) SortHits(0, hitArray); // sort
2335   //if (GetStation0() == 3) SortHits(1, hitArray); // unsort
2336
2337   Bool_t same = kFALSE;
2338   for (Int_t i=0; i<nRecTracks; i++) {
2339     AliMUONTrackK *trackK = (AliMUONTrackK*) ((*trackPtr)[i]);
2340     if (trackK->fNmbTrackHits == 2 && trackK->GetRecover() == 0) continue;
2341     if (trackK == this) continue;
2342     //AZ if (trackK->GetRecover() < 0) continue; //
2343     if (trackK->fNmbTrackHits >= fNmbTrackHits) {
2344       TObjArray *hits = new TObjArray(*trackK->fTrackHits);
2345       if (trackK->GetStation0() == 3) SortHits(0, hits); // sort
2346       //if (trackK->GetStation0() == 3) SortHits(1, hits); // unsort
2347       for (Int_t j=0; j<fNmbTrackHits; j++) {
2348         //cout << fNmbTrackHits << " " << i << " " << j << " " << (*hitArray)[j] << " " << (*hits)[j] << " " << trackK->fSkipHit << endl;
2349         if (j != fNmbTrackHits-1 && (*hitArray)[j] != (*hits)[j]) break;
2350         if (j == fNmbTrackHits-1) {
2351           if (trackK->fSkipHit && TMath::Abs(((AliMUONHitForRec*)((*hitArray)[j]))->GetZ()-trackK->fSkipHit->GetZ()) < 0.5) same = kTRUE;
2352           //if (trackK->fNmbTrackHits > fNmbTrackHits &&
2353           //if (trackK->fSkipHit) cout << ((AliMUONHitForRec*)((*hitArray)[j]))->GetZ() << " " << trackK->fSkipHit->GetZ() << endl;
2354         }
2355       } // for (Int_t j=0;
2356       delete hits;
2357       if (same) break;
2358     }
2359   } // for (Int_t i=0;
2360   delete hitArray;
2361   return same;
2362 }