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