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