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