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