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