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