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