f096d8ea234367b80c926b6f3971c97afba9688e
[u/mrichter/AliRoot.git] / MUON / AliMUONTrack.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 ///////////////////////////////////////////////////
19 //
20 // Reconstructed track
21 // in
22 // ALICE
23 // dimuon
24 // spectrometer
25 //
26 ///////////////////////////////////////////////////
27
28 #include "AliMUONTrack.h"
29
30 #include "AliMUONTrackReconstructor.h" 
31 #include "AliMUONHitForRec.h" 
32 #include "AliMUONSegment.h" 
33 #include "AliMUONTrackHit.h"
34 #include "AliMUONTriggerTrack.h"
35 #include "AliMUONConstants.h"
36
37 #include "AliLog.h"
38
39 #include <Riostream.h> // for cout
40 #include <TMath.h>
41 #include <TMatrixD.h>
42 #include <TObjArray.h>
43 #include <TVirtualFitter.h>
44
45 #include <stdlib.h> // for exit()
46
47 // Functions to be minimized with Minuit
48 void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag);
49 void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag);
50
51 void mnvertLocal(Double_t* a, Int_t l, Int_t m, Int_t n, Int_t& ifail);
52
53 Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit);
54
55 TVirtualFitter* AliMUONTrack::fgFitter = NULL; 
56
57 ClassImp(AliMUONTrack) // Class implementation in ROOT context
58
59 //__________________________________________________________________________
60 AliMUONTrack::AliMUONTrack()
61   : TObject() 
62 {
63   // Default constructor
64   fgFitter = 0;
65   fTrackReconstructor = 0;
66   fTrackHitsPtr = NULL;
67   fTrackParamAtHit = NULL;
68   fHitForRecAtHit = NULL;
69   fTrackID = 0;
70 }
71
72   //__________________________________________________________________________
73 AliMUONTrack::AliMUONTrack(AliMUONSegment* BegSegment, AliMUONSegment* EndSegment, AliMUONTrackReconstructor* TrackReconstructor)
74   : TObject()
75 {
76   // Constructor from two Segment's
77   fTrackReconstructor = TrackReconstructor; // link back to TrackReconstructor
78   // memory allocation for the TObjArray of pointers to reconstructed TrackHit's
79   fTrackHitsPtr = new TObjArray(10);
80   fNTrackHits = 0;
81   if (BegSegment) { //AZ
82     AddSegment(BegSegment); // add hits from BegSegment
83     AddSegment(EndSegment); // add hits from EndSegment
84     fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z
85     SetTrackParamAtVertex(); // set track parameters at vertex
86   }
87   fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",10);
88   fHitForRecAtHit = new TClonesArray("AliMUONHitForRec",10);
89   // set fit conditions...
90   fFitMCS = 0;
91   fFitNParam = 3;
92   fFitStart = 1;
93   fFitFMin = -1.0;
94   fMatchTrigger = kFALSE;
95   fChi2MatchTrigger = 0;
96   fTrackID = 0;
97   return;
98 }
99
100   //__________________________________________________________________________
101 AliMUONTrack::AliMUONTrack(AliMUONSegment* Segment, AliMUONHitForRec* HitForRec, AliMUONTrackReconstructor* TrackReconstructor)
102   : TObject()
103 {
104   // Constructor from one Segment and one HitForRec
105   fTrackReconstructor = TrackReconstructor; // link back to TrackReconstructor
106   // memory allocation for the TObjArray of pointers to reconstructed TrackHit's
107   fTrackHitsPtr = new TObjArray(10);
108   fNTrackHits = 0;
109   AddSegment(Segment); // add hits from Segment
110   AddHitForRec(HitForRec); // add HitForRec
111   fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z
112   SetTrackParamAtVertex(); // set track parameters at vertex
113   fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",10);
114   fHitForRecAtHit = new TClonesArray("AliMUONHitForRec",10);
115   // set fit conditions...
116   fFitMCS = 0;
117   fFitNParam = 3;
118   fFitStart = 1;
119   fFitFMin = -1.0;
120   fMatchTrigger = kFALSE;
121   fChi2MatchTrigger = 0;
122   fTrackID = 0;
123   return;
124 }
125
126   //__________________________________________________________________________
127 AliMUONTrack::~AliMUONTrack()
128 {
129   // Destructor
130   if (fTrackHitsPtr) {
131     // delete the TObjArray of pointers to TrackHit's   
132     delete fTrackHitsPtr; 
133     fTrackHitsPtr = NULL;
134   }
135   
136   if (fTrackParamAtHit) {
137     // delete the TClonesArray of pointers to TrackParam
138     delete fTrackParamAtHit;
139     fTrackParamAtHit = NULL;
140   }
141
142   if (fHitForRecAtHit) {
143     // delete the TClonesArray of pointers to HitForRec
144     delete fHitForRecAtHit;
145     fHitForRecAtHit = NULL;
146   }
147 }
148
149   //__________________________________________________________________________
150 AliMUONTrack::AliMUONTrack (const AliMUONTrack& theMUONTrack)
151   :  TObject(theMUONTrack)
152 {
153   //fTrackReconstructor = new AliMUONTrackReconstructor(*MUONTrack.fTrackReconstructor);
154                                // is it right ?
155                                // NO, because it would use dummy copy constructor
156                                // and AliMUONTrack is not the owner of its TrackReconstructor 
157   fTrackReconstructor = theMUONTrack.fTrackReconstructor;
158   fTrackParamAtVertex = theMUONTrack.fTrackParamAtVertex;
159
160  // necessary to make a copy of the objects and not only the pointers in TObjArray.
161   fTrackHitsPtr = 0;
162   if (theMUONTrack.fTrackHitsPtr) {
163     fTrackHitsPtr  =  new TObjArray(10);
164     for (Int_t index = 0; index < (theMUONTrack.fTrackHitsPtr)->GetEntriesFast(); index++) {
165       AliMUONTrackHit *trackHit = new AliMUONTrackHit(*(AliMUONTrackHit*)(theMUONTrack.fTrackHitsPtr)->At(index));
166       fTrackHitsPtr->Add(trackHit);
167     }
168     fTrackHitsPtr->SetOwner(); // nedeed for deleting TClonesArray
169   }
170
171   // necessary to make a copy of the objects and not only the pointers in TClonesArray.
172   fTrackParamAtHit = 0;
173   if (theMUONTrack.fTrackParamAtHit) {
174     fTrackParamAtHit  =  new TClonesArray("AliMUONTrackParam",10);
175     for (Int_t index = 0; index < (theMUONTrack.fTrackParamAtHit)->GetEntriesFast(); index++) {
176       {new ((*fTrackParamAtHit)[fTrackParamAtHit->GetEntriesFast()]) 
177           AliMUONTrackParam(*(AliMUONTrackParam*)(theMUONTrack.fTrackParamAtHit)->At(index));}
178     }
179   }  
180
181   // necessary to make a copy of the objects and not only the pointers in TClonesArray.
182   fHitForRecAtHit = 0;
183   if (theMUONTrack.fHitForRecAtHit) {
184     fHitForRecAtHit  =  new TClonesArray("AliMUONHitForRec",10);
185     for (Int_t index = 0; index < (theMUONTrack.fHitForRecAtHit)->GetEntriesFast(); index++) {
186       {new ((*fHitForRecAtHit)[fHitForRecAtHit->GetEntriesFast()]) 
187           AliMUONHitForRec(*(AliMUONHitForRec*)(theMUONTrack.fHitForRecAtHit)->At(index));}
188     }
189   }  
190
191   fNTrackHits       =  theMUONTrack.fNTrackHits;
192   fFitMCS           =  theMUONTrack.fFitMCS;
193   fFitNParam        =  theMUONTrack.fFitNParam;
194   fFitFMin          =  theMUONTrack.fFitFMin;
195   fFitStart         =  theMUONTrack.fFitStart;
196   fMatchTrigger     =  theMUONTrack.fMatchTrigger;
197   fChi2MatchTrigger =  theMUONTrack.fChi2MatchTrigger;
198   fTrackID          =  theMUONTrack.fTrackID;
199 }
200
201   //__________________________________________________________________________
202 AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& theMUONTrack)
203 {
204
205   // check assignement to self
206   if (this == &theMUONTrack)
207     return *this;
208
209   // base class assignement
210   TObject::operator=(theMUONTrack);
211
212   // fTrackReconstructor =  new AliMUONTrackReconstructor(*MUONTrack.fTrackReconstructor); 
213   // is it right ?
214   // is it right ? NO because it would use dummy copy constructor
215   fTrackReconstructor =  theMUONTrack.fTrackReconstructor;
216   fTrackParamAtVertex =  theMUONTrack.fTrackParamAtVertex;
217
218  // necessary to make a copy of the objects and not only the pointers in TObjArray.
219   fTrackHitsPtr = 0;
220   if (theMUONTrack.fTrackHitsPtr) {
221     fTrackHitsPtr  =  new TObjArray(10);
222     for (Int_t index = 0; index < (theMUONTrack.fTrackHitsPtr)->GetEntriesFast(); index++) {
223       AliMUONTrackHit *trackHit = new AliMUONTrackHit(*(AliMUONTrackHit*)(theMUONTrack.fTrackHitsPtr)->At(index));
224       fTrackHitsPtr->Add(trackHit);
225     }
226     fTrackHitsPtr->SetOwner();  // nedeed for deleting TClonesArray
227   }  
228
229   // necessary to make a copy of the objects and not only the pointers in TClonesArray.
230   fTrackParamAtHit = 0;
231   if (theMUONTrack.fTrackParamAtHit) {
232     fTrackParamAtHit  =  new TClonesArray("AliMUONTrackParam",10);
233     for (Int_t index = 0; index < (theMUONTrack.fTrackParamAtHit)->GetEntriesFast(); index++) {
234       {new ((*fTrackParamAtHit)[fTrackParamAtHit->GetEntriesFast()]) 
235           AliMUONTrackParam(*(AliMUONTrackParam*)(theMUONTrack.fTrackParamAtHit)->At(index));}
236     }
237   }  
238
239   // necessary to make a copy of the objects and not only the pointers in TClonesArray.
240   fHitForRecAtHit = 0;
241   if (theMUONTrack.fHitForRecAtHit) {  
242     fHitForRecAtHit  =  new TClonesArray("AliMUONHitForRec",10);
243     for (Int_t index = 0; index < (theMUONTrack.fHitForRecAtHit)->GetEntriesFast(); index++) {
244       {new ((*fHitForRecAtHit)[fHitForRecAtHit->GetEntriesFast()]) 
245           AliMUONHitForRec(*(AliMUONHitForRec*)(theMUONTrack.fHitForRecAtHit)->At(index));}
246     }
247   }  
248
249   fNTrackHits         =  theMUONTrack.fNTrackHits;
250   fFitMCS             =  theMUONTrack.fFitMCS;
251   fFitNParam          =  theMUONTrack.fFitNParam;
252   fFitFMin            =  theMUONTrack.fFitFMin;
253   fFitStart           =  theMUONTrack.fFitStart;
254   fMatchTrigger       =  theMUONTrack.fMatchTrigger;
255   fChi2MatchTrigger   =  theMUONTrack.fChi2MatchTrigger;
256   fTrackID            =  theMUONTrack.fTrackID;
257
258   return *this;
259 }
260
261   //__________________________________________________________________________
262 void AliMUONTrack::Remove()
263 {
264   // Remove current track from array of tracks,
265   // and corresponding track hits from array of track hits.
266   // Compress the TClonesArray it belongs to.
267   AliMUONTrackHit *nextTrackHit;
268   AliMUONTrackReconstructor *eventRec = this->fTrackReconstructor;
269   TClonesArray *trackHitsPtr = eventRec->GetRecTrackHitsPtr();
270   // Loop over all track hits of track
271   AliMUONTrackHit *trackHit = (AliMUONTrackHit*) fTrackHitsPtr->First();
272   while (trackHit) {
273     nextTrackHit = (AliMUONTrackHit*) fTrackHitsPtr->After(trackHit);
274     // Remove TrackHit from event TClonesArray.
275     // Destructor is called,
276     // hence links between HitForRec's and TrackHit's are updated
277     trackHitsPtr->Remove(trackHit);
278     trackHit = nextTrackHit;
279   }
280   // Remove the track from event TClonesArray
281   // Destructor is called,
282   // hence space for TObjArray of pointers to TrackHit's is freed
283   eventRec->GetRecTracksPtr()->Remove(this);
284   // Number of tracks decreased by 1
285   eventRec->SetNRecTracks(eventRec->GetNRecTracks() - 1);
286   // Compress event TClonesArray of Track's:
287   // this is essential to retrieve the TClonesArray afterwards
288   eventRec->GetRecTracksPtr()->Compress();
289   // Compress event TClonesArray of TrackHit's:
290   // this is probably also essential to retrieve the TClonesArray afterwards
291   trackHitsPtr->Compress();
292 }
293
294   //__________________________________________________________________________
295 void AliMUONTrack::SetFitMCS(Int_t FitMCS)
296 {
297   // Set multiple Coulomb scattering option for track fit "fFitMCS"
298   // from "FitMCS" argument: 0 without, 1 with
299   if ((FitMCS == 0) || (FitMCS == 1)) fFitMCS = FitMCS;
300   // better implementation with enum(with, without) ????
301   else {
302     cout << "ERROR in AliMUONTrack::SetFitMCS(FitMCS)" << endl;
303     cout << "FitMCS = " << FitMCS << " is neither 0 nor 1" << endl;
304     exit(0);
305   }
306   return;
307 }
308
309   //__________________________________________________________________________
310 void AliMUONTrack::SetFitNParam(Int_t FitNParam)
311 {
312   // Set number of parameters for track fit "fFitNParam" from "FitNParam":
313   // 3 for momentum, 5 for momentum and position
314   if ((FitNParam == 3) || (FitNParam == 5)) fFitNParam = FitNParam;
315   else {
316     cout << "ERROR in AliMUONTrack::SetFitNParam(FitNParam)" << endl;
317     cout << "FitNParam = " << FitNParam << " is neither 3 nor 5" << endl;
318     exit(0);
319   }
320   return;
321 }
322
323   //__________________________________________________________________________
324 void AliMUONTrack::SetFitStart(Int_t FitStart)
325 {
326   // Set multiple Coulomb scattering option for track fit "fFitStart"
327   // from "FitStart" argument: 0 without, 1 with
328   if ((FitStart == 0) || (FitStart == 1)) fFitStart = FitStart;
329   // better implementation with enum(vertex, firstHit) ????
330   else {
331     cout << "ERROR in AliMUONTrack::SetFitStart(FitStart)" << endl;
332     cout << "FitStart = " << FitStart << " is neither 0 nor 1" << endl;
333     exit(0);
334   }
335   return;
336 }
337
338   //__________________________________________________________________________
339 AliMUONTrackParam* AliMUONTrack::GetTrackParamAtFirstHit(void) const {
340   // Get pointer to TrackParamAtFirstHit
341   return ((AliMUONTrackHit*) (fTrackHitsPtr->First()))->GetTrackParam();}
342
343   //__________________________________________________________________________
344 void AliMUONTrack::RecursiveDump(void) const
345 {
346   // Recursive dump of AliMUONTrack, i.e. with dump of TrackHit's and HitForRec's
347   AliMUONTrackHit *trackHit;
348   AliMUONHitForRec *hitForRec;
349   cout << "Recursive dump of Track: " << this << endl;
350   // Track
351   this->Dump();
352   for (Int_t trackHitIndex = 0; trackHitIndex < fNTrackHits; trackHitIndex++) {
353     trackHit = (AliMUONTrackHit*) ((*fTrackHitsPtr)[trackHitIndex]);
354     // TrackHit
355     cout << "TrackHit: " << trackHit << " (index: " << trackHitIndex << ")" << endl;
356     trackHit->Dump();
357     hitForRec = trackHit->GetHitForRecPtr();
358     // HitForRec
359     cout << "HitForRec: " << hitForRec << endl;
360     hitForRec->Dump();
361   }
362   return;
363 }
364   
365   //__________________________________________________________________________
366 Bool_t* AliMUONTrack::CompatibleTrack(AliMUONTrack * Track, Double_t Sigma2Cut) const
367 {
368   // Return kTRUE/kFALSE for each chamber if hit is compatible or not 
369   TClonesArray *hitArray, *thisHitArray;
370   AliMUONHitForRec *hit, *thisHit;
371   Int_t chamberNumber;
372   Float_t deltaZ;
373   Float_t deltaZMax = 1.; // 1 cm
374   Float_t chi2 = 0;
375   Bool_t *nCompHit = new Bool_t[AliMUONConstants::NTrackingCh()]; 
376
377   for ( Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
378     nCompHit[ch] = kFALSE;
379   }
380
381   thisHitArray = this->GetHitForRecAtHit();
382
383   hitArray =  Track->GetHitForRecAtHit();
384
385   for (Int_t iHthis = 0; iHthis < thisHitArray->GetEntriesFast(); iHthis++) {
386     thisHit = (AliMUONHitForRec*) thisHitArray->At(iHthis);
387     chamberNumber = thisHit->GetChamberNumber();
388     if (chamberNumber < 0 || chamberNumber > AliMUONConstants::NTrackingCh()) continue; 
389     nCompHit[chamberNumber] = kFALSE;
390     for (Int_t iH = 0; iH < hitArray->GetEntriesFast(); iH++) {
391       hit = (AliMUONHitForRec*) hitArray->At(iH);
392       deltaZ = TMath::Abs(thisHit->GetZ() - hit->GetZ());
393       chi2 = thisHit->NormalizedChi2WithHitForRec(hit,Sigma2Cut); // set cut to 4 sigmas
394       if (chi2 < 3. && deltaZ < deltaZMax) {
395         nCompHit[chamberNumber] = kTRUE;
396         break;
397       }
398     }  
399   }
400   
401   return nCompHit;
402 }
403
404   //__________________________________________________________________________
405 Int_t AliMUONTrack::HitsInCommon(AliMUONTrack* Track) const
406 {
407   // Returns the number of hits in common
408   // between the current track ("this")
409   // and the track pointed to by "Track".
410   Int_t hitsInCommon = 0;
411   AliMUONTrackHit *trackHit1, *trackHit2;
412   // Loop over hits of first track
413   trackHit1 = (AliMUONTrackHit*) this->GetTrackHitsPtr()->First();
414   while (trackHit1) {
415     // Loop over hits of second track
416     trackHit2 = (AliMUONTrackHit*) Track->GetTrackHitsPtr()->First();
417     while (trackHit2) {
418       // Increment "hitsInCommon" if both TrackHits point to the same HitForRec
419       if ( (trackHit1->GetHitForRecPtr()) ==
420            (trackHit2->GetHitForRecPtr())    ) hitsInCommon++;
421       trackHit2 = (AliMUONTrackHit*) Track->GetTrackHitsPtr()->After(trackHit2);
422     } // trackHit2
423     trackHit1 = (AliMUONTrackHit*) this->GetTrackHitsPtr()->After(trackHit1);
424   } // trackHit1
425   return hitsInCommon;
426 }
427
428   //__________________________________________________________________________
429 void AliMUONTrack::MatchTriggerTrack(TClonesArray *triggerTrackArray)
430 {
431   // Match this track with one trigger track if possible
432   AliMUONTrackParam trackParam; 
433   AliMUONTriggerTrack *triggerTrack;
434   Double_t xTrack, yTrack, ySlopeTrack, dTrigTrackMin2, dTrigTrack2;
435   Double_t nSigmaCut2;
436
437   Double_t distSigma[3]={1,1,0.02}; // sigma of distributions (trigger-track) X,Y,slopeY
438   Double_t distTriggerTrack[3] = {0,0,0};
439
440   fMatchTrigger = kFALSE;
441   fChi2MatchTrigger = 0;
442
443   trackParam = *((AliMUONTrackParam*) fTrackParamAtHit->Last()); 
444   trackParam.ExtrapToZ(AliMUONConstants::DefaultChamberZ(10)); // extrap to 1st trigger chamber
445
446   nSigmaCut2 =  fTrackReconstructor->GetMaxSigma2Distance(); // nb of sigma**2 for cut
447   xTrack = trackParam.GetNonBendingCoor();
448   yTrack = trackParam.GetBendingCoor();
449   ySlopeTrack = trackParam.GetBendingSlope();
450   dTrigTrackMin2 = 999;
451   
452   triggerTrack = (AliMUONTriggerTrack*) triggerTrackArray->First();
453   while(triggerTrack){
454     distTriggerTrack[0] = (triggerTrack->GetX11()-xTrack)/distSigma[0];
455     distTriggerTrack[1] = (triggerTrack->GetY11()-yTrack)/distSigma[1];
456     distTriggerTrack[2] = (TMath::Tan(triggerTrack->GetThetay())-ySlopeTrack)/distSigma[2];
457     dTrigTrack2 = 0;
458     for (Int_t iVar = 0; iVar < 3; iVar++)
459       dTrigTrack2 += distTriggerTrack[iVar]*distTriggerTrack[iVar];
460     if (dTrigTrack2 < dTrigTrackMin2 && dTrigTrack2 < nSigmaCut2) {
461       dTrigTrackMin2 = dTrigTrack2;
462       fMatchTrigger = kTRUE;
463       fChi2MatchTrigger =  dTrigTrack2/3.; // Normalized Chi2, 3 variables (X,Y,slopeY)
464     }
465     triggerTrack = (AliMUONTriggerTrack*) triggerTrackArray->After(triggerTrack);
466   }
467
468 }
469   //__________________________________________________________________________
470 void AliMUONTrack::Fit()
471 {
472   // Fit the current track ("this"),
473   // with or without multiple Coulomb scattering according to "fFitMCS",
474   // with the number of parameters given by "fFitNParam"
475   // (3 if one keeps X and Y fixed in "TrackParam", 5 if one lets them vary),
476   // starting, according to "fFitStart",
477   // with track parameters at vertex or at the first TrackHit.
478   // "fFitMCS", "fFitNParam" and "fFitStart" have to be set before
479   // by calling the corresponding Set methods.
480   Double_t arg[1], benC, errorParam, invBenP, lower, nonBenC, upper, x, y;
481   char parName[50];
482   AliMUONTrackParam *trackParam;
483   // Check if Minuit is initialized...
484   fgFitter = TVirtualFitter::Fitter(this); // add 3 or 5 for the maximum number of parameters ???
485   fgFitter->Clear(); // necessary ???? probably yes
486   // how to reset the printout number at every fit ????
487   // is there any risk to leave it like that ????
488   // how to go faster ???? choice of Minuit parameters like EDM ????
489   // choice of function to be minimized according to fFitMCS
490   if (fFitMCS == 0) fgFitter->SetFCN(TrackChi2);
491   else fgFitter->SetFCN(TrackChi2MCS);
492   // Switch off printout
493   arg[0] = -1;
494   fgFitter->ExecuteCommand("SET PRINT", arg, 1); // More printing !!!!
495   // No warnings
496   fgFitter->ExecuteCommand("SET NOW", arg, 0);
497   // Parameters according to "fFitStart"
498   // (should be a function to be used at every place where needed ????)
499   if (fFitStart == 0) trackParam = &fTrackParamAtVertex;
500   else trackParam = this->GetTrackParamAtFirstHit();
501   // set first 3 Minuit parameters
502   // could be tried with no limits for the search (min=max=0) ????
503   fgFitter->SetParameter(0, "InvBenP",
504                          trackParam->GetInverseBendingMomentum(),
505                          0.003, -0.4, 0.4);
506   fgFitter->SetParameter(1, "BenS",
507                          trackParam->GetBendingSlope(),
508                          0.001, -0.5, 0.5);
509   fgFitter->SetParameter(2, "NonBenS",
510                          trackParam->GetNonBendingSlope(),
511                          0.001, -0.5, 0.5);
512   if (fFitNParam == 5) {
513     // set last 2 Minuit parameters
514     // mandatory limits in Bending to avoid NaN values of parameters
515     fgFitter->SetParameter(3, "X",
516                            trackParam->GetNonBendingCoor(),
517                            0.03, -500.0, 500.0);
518     // mandatory limits in non Bending to avoid NaN values of parameters
519     fgFitter->SetParameter(4, "Y",
520                            trackParam->GetBendingCoor(),
521                            0.10, -500.0, 500.0);
522   }
523   // search without gradient calculation in the function
524   fgFitter->ExecuteCommand("SET NOGRADIENT", arg, 0);
525   // minimization
526   fgFitter->ExecuteCommand("MINIMIZE", arg, 0);
527   // exit from Minuit
528   //  fgFitter->ExecuteCommand("EXIT", arg, 0); // necessary ????
529   // get results into "invBenP", "benC", "nonBenC" ("x", "y")
530   fgFitter->GetParameter(0, parName, invBenP, errorParam, lower, upper);
531   fgFitter->GetParameter(1, parName, benC, errorParam, lower, upper);
532   fgFitter->GetParameter(2, parName, nonBenC, errorParam, lower, upper);
533   if (fFitNParam == 5) {
534     fgFitter->GetParameter(3, parName, x, errorParam, lower, upper);
535     fgFitter->GetParameter(4, parName, y, errorParam, lower, upper);
536   }
537   // result of the fit into track parameters
538   trackParam->SetInverseBendingMomentum(invBenP);
539   trackParam->SetBendingSlope(benC);
540   trackParam->SetNonBendingSlope(nonBenC);
541   if (fFitNParam == 5) {
542     trackParam->SetNonBendingCoor(x);
543     trackParam->SetBendingCoor(y);
544   }
545   // global result of the fit
546   Double_t fedm, errdef;
547   Int_t npari, nparx;
548   fgFitter->GetStats(fFitFMin, fedm, errdef, npari, nparx);
549 }
550
551   //__________________________________________________________________________
552 void AliMUONTrack::AddSegment(AliMUONSegment* Segment)
553 {
554   // Add Segment to the track
555   AddHitForRec(Segment->GetHitForRec1()); // 1st hit
556   AddHitForRec(Segment->GetHitForRec2()); // 2nd hit
557 }
558
559   //__________________________________________________________________________
560 void AliMUONTrack::AddHitForRec(AliMUONHitForRec* HitForRec)
561 {
562   // Add HitForRec to the track:
563   // actual TrackHit into TClonesArray of TrackHit's for the event;
564   // pointer to actual TrackHit in TObjArray of pointers to TrackHit's for the track
565   TClonesArray *recTrackHitsPtr = this->fTrackReconstructor->GetRecTrackHitsPtr();
566   Int_t eventTrackHits = this->fTrackReconstructor->GetNRecTrackHits();
567   // event
568   AliMUONTrackHit* trackHit =
569     new ((*recTrackHitsPtr)[eventTrackHits]) AliMUONTrackHit(HitForRec);
570   this->fTrackReconstructor->SetNRecTrackHits(eventTrackHits + 1);
571   // track
572   if (fTrackHitsPtr->IsOwner()) AliFatal("fTrackHitsPtr is owner");
573   fTrackHitsPtr->Add(trackHit);
574   fNTrackHits++;
575 }
576
577   //__________________________________________________________________________
578 void AliMUONTrack::SetTrackParamAtHit(Int_t indexHit, AliMUONTrackParam *TrackParam) const
579 {
580   // Set track parameters at TrackHit with index "indexHit"
581   // from the track parameters pointed to by "TrackParam".
582   //PH  AliMUONTrackHit* trackHit = (AliMUONTrackHit*) ((*fTrackHitsPtr)[indexHit]);
583   AliMUONTrackHit* trackHit = (AliMUONTrackHit*) (fTrackHitsPtr->At(indexHit));
584   trackHit->SetTrackParam(TrackParam);
585 }
586
587   //__________________________________________________________________________
588 void AliMUONTrack::SetTrackParamAtVertex()
589 {
590   // Set track parameters at vertex.
591   // TrackHit's are assumed to be only in stations(1..) 4 and 5,
592   // and sorted according to increasing Z..
593   // Parameters are calculated from information in HitForRec's
594   // of first and last TrackHit's.
595   AliMUONTrackParam *trackParam =
596     &fTrackParamAtVertex; // pointer to track parameters
597   // Pointer to HitForRec of first TrackHit
598   AliMUONHitForRec *firstHit =
599     ((AliMUONTrackHit*) (fTrackHitsPtr->First()))->GetHitForRecPtr();
600   // Pointer to HitForRec of last TrackHit
601   AliMUONHitForRec *lastHit =
602     ((AliMUONTrackHit*) (fTrackHitsPtr->Last()))->GetHitForRecPtr();
603   // Z difference between first and last hits
604   Double_t deltaZ = firstHit->GetZ() - lastHit->GetZ();
605   // bending slope in stations(1..) 4 and 5
606   Double_t bendingSlope =
607     (firstHit->GetBendingCoor() - lastHit->GetBendingCoor()) / deltaZ;
608   trackParam->SetBendingSlope(bendingSlope);
609   // impact parameter
610   Double_t impactParam =
611     firstHit->GetBendingCoor() - bendingSlope * firstHit->GetZ(); // same if from firstHit and  lastHit ????
612   // signed bending momentum
613   Double_t signedBendingMomentum =
614     fTrackReconstructor->GetBendingMomentumFromImpactParam(impactParam);
615   trackParam->SetInverseBendingMomentum(1.0 / signedBendingMomentum);
616   // bending slope at vertex
617   trackParam->
618     SetBendingSlope(bendingSlope +
619                     impactParam / fTrackReconstructor->GetSimpleBPosition());
620   // non bending slope
621   Double_t nonBendingSlope =
622     (firstHit->GetNonBendingCoor() - lastHit->GetNonBendingCoor()) / deltaZ;
623   trackParam->SetNonBendingSlope(nonBendingSlope);
624   // vertex coordinates at (0,0,0)
625   trackParam->SetZ(0.0);
626   trackParam->SetBendingCoor(0.0);
627   trackParam->SetNonBendingCoor(0.0);
628 }
629
630   //__________________________________________________________________________
631 void AliMUONTrack::AddTrackParamAtHit(const AliMUONTrackParam *trackParam) 
632 {
633   // Add track paremeters
634
635   if (!fTrackParamAtHit)
636     fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",10);  
637
638   new ((*fTrackParamAtHit)[fTrackParamAtHit->GetEntriesFast()]) AliMUONTrackParam(*trackParam);
639 }
640
641   //__________________________________________________________________________
642 void AliMUONTrack::AddHitForRecAtHit(const AliMUONHitForRec *hitForRec) 
643 {
644   if (!fHitForRecAtHit)
645     fHitForRecAtHit = new TClonesArray("AliMUONHitForRec",10); 
646
647   new ((*fHitForRecAtHit)[fHitForRecAtHit->GetEntriesFast()]) AliMUONHitForRec(*hitForRec);
648 }
649
650   //__________________________________________________________________________
651 void TrackChi2(Int_t &NParam, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/)
652 {
653   // Return the "Chi2" to be minimized with Minuit for track fitting,
654   // with "NParam" parameters
655   // and their current values in array pointed to by "Param".
656   // Assumes that the track hits are sorted according to increasing Z.
657   // Track parameters at each TrackHit are updated accordingly.
658   // Multiple Coulomb scattering is not taken into account
659   AliMUONTrack *trackBeingFitted;
660   AliMUONTrackHit* hit;
661   AliMUONTrackParam param1;
662   Int_t hitNumber;
663   Double_t zHit;
664   Chi2 = 0.0; // initialize Chi2
665   // copy of track parameters to be fitted
666   trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit();
667   if (trackBeingFitted->GetFitStart() == 0)
668     param1 = *(trackBeingFitted->GetTrackParamAtVertex());
669   else param1 = *(trackBeingFitted->GetTrackParamAtFirstHit());
670   // Minuit parameters to be fitted into this copy
671   param1.SetInverseBendingMomentum(Param[0]);
672   param1.SetBendingSlope(Param[1]);
673   param1.SetNonBendingSlope(Param[2]);
674   if (NParam == 5) {
675     param1.SetNonBendingCoor(Param[3]);
676     param1.SetBendingCoor(Param[4]);
677   }
678   // Follow track through all planes of track hits
679   for (hitNumber = 0; hitNumber < trackBeingFitted->GetNTrackHits(); hitNumber++) {
680     hit = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber];
681     zHit = hit->GetHitForRecPtr()->GetZ();
682     // do something special if 2 hits with same Z ????
683     // security against infinite loop ????
684     (&param1)->ExtrapToZ(zHit); // extrapolation
685     hit->SetTrackParam(&param1);
686     // Increment Chi2
687     // done hit per hit, with hit resolution,
688     // and not with point and angle like in "reco_muon.F" !!!!
689     // Needs to add multiple scattering contribution ????
690     Double_t dX =
691       hit->GetHitForRecPtr()->GetNonBendingCoor() - (&param1)->GetNonBendingCoor();
692     Double_t dY =
693       hit->GetHitForRecPtr()->GetBendingCoor() - (&param1)->GetBendingCoor();
694     Chi2 =
695       Chi2 +
696       dX * dX / hit->GetHitForRecPtr()->GetNonBendingReso2() +
697       dY * dY / hit->GetHitForRecPtr()->GetBendingReso2();
698   }
699 }
700
701   //__________________________________________________________________________
702 void TrackChi2MCS(Int_t &NParam, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/)
703 {
704   // Return the "Chi2" to be minimized with Minuit for track fitting,
705   // with "NParam" parameters
706   // and their current values in array pointed to by "Param".
707   // Assumes that the track hits are sorted according to increasing Z.
708   // Track parameters at each TrackHit are updated accordingly.
709   // Multiple Coulomb scattering is taken into account with covariance matrix.
710   AliMUONTrack *trackBeingFitted;
711   AliMUONTrackParam param1;
712   Chi2 = 0.0; // initialize Chi2
713   // copy of track parameters to be fitted
714   trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit();
715   if (trackBeingFitted->GetFitStart() == 0)
716     param1 = *(trackBeingFitted->GetTrackParamAtVertex());
717   else param1 = *(trackBeingFitted->GetTrackParamAtFirstHit());
718   // Minuit parameters to be fitted into this copy
719   param1.SetInverseBendingMomentum(Param[0]);
720   param1.SetBendingSlope(Param[1]);
721   param1.SetNonBendingSlope(Param[2]);
722   if (NParam == 5) {
723     param1.SetNonBendingCoor(Param[3]);
724     param1.SetBendingCoor(Param[4]);
725   }
726
727   AliMUONTrackHit *hit;
728   Int_t chCurrent, chPrev = 0, hitNumber, hitNumber1, hitNumber2, hitNumber3;
729   Double_t z, z1, z2, z3;
730   AliMUONTrackHit *hit1, *hit2, *hit3;
731   Double_t hbc1, hbc2, pbc1, pbc2;
732   Double_t hnbc1, hnbc2, pnbc1, pnbc2;
733   Int_t numberOfHit = trackBeingFitted->GetNTrackHits();
734   TMatrixD *covBending = new TMatrixD(numberOfHit, numberOfHit);
735   TMatrixD *covNonBending = new TMatrixD(numberOfHit, numberOfHit);
736   Double_t *msa2 = new Double_t[numberOfHit];
737
738   // Predicted coordinates and  multiple scattering angles are first calculated
739   for (hitNumber = 0; hitNumber < numberOfHit; hitNumber++) {
740     hit = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber];
741     z = hit->GetHitForRecPtr()->GetZ();
742     // do something special if 2 hits with same Z ????
743     // security against infinite loop ????
744     (&param1)->ExtrapToZ(z); // extrapolation
745     hit->SetTrackParam(&param1);
746     // square of multiple scattering angle at current hit, with one chamber
747     msa2[hitNumber] = MultipleScatteringAngle2(hit);
748     // correction for eventual missing hits or multiple hits in a chamber,
749     // according to the number of chambers
750     // between the current hit and the previous one
751     chCurrent = hit->GetHitForRecPtr()->GetChamberNumber();
752     if (hitNumber > 0) msa2[hitNumber] = msa2[hitNumber] * (chCurrent - chPrev);
753     chPrev = chCurrent;
754   }
755
756   // Calculates the covariance matrix
757   for (hitNumber1 = 0; hitNumber1 < numberOfHit; hitNumber1++) { 
758     hit1 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber1];
759     z1 = hit1->GetHitForRecPtr()->GetZ();
760     for (hitNumber2 = hitNumber1; hitNumber2 < numberOfHit; hitNumber2++) {
761       hit2 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber2];
762       z2 = hit2->GetHitForRecPtr()->GetZ();
763       // initialization to 0 (diagonal plus upper triangular part)
764       (*covBending)(hitNumber2, hitNumber1) = 0.0;
765       // contribution from multiple scattering in bending plane:
766       // loop over upstream hits
767       for (hitNumber3 = 0; hitNumber3 < hitNumber1; hitNumber3++) {     
768         hit3 = (AliMUONTrackHit*)
769           (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber3];
770         z3 = hit3->GetHitForRecPtr()->GetZ();
771         (*covBending)(hitNumber2, hitNumber1) =
772           (*covBending)(hitNumber2, hitNumber1) +
773           ((z1 - z3) * (z2 - z3) * msa2[hitNumber3]); 
774       }
775       // equal contribution from multiple scattering in non bending plane
776       (*covNonBending)(hitNumber2, hitNumber1) =
777         (*covBending)(hitNumber2, hitNumber1);
778       if (hitNumber1 == hitNumber2) {
779         // Diagonal elements: add contribution from position measurements
780         // in bending plane
781         (*covBending)(hitNumber2, hitNumber1) =
782           (*covBending)(hitNumber2, hitNumber1) +
783           hit1->GetHitForRecPtr()->GetBendingReso2();
784         // and in non bending plane
785         (*covNonBending)(hitNumber2, hitNumber1) =
786           (*covNonBending)(hitNumber2, hitNumber1) +
787           hit1->GetHitForRecPtr()->GetNonBendingReso2();
788       }
789       else {
790         // Non diagonal elements: symmetrization
791         // for bending plane
792         (*covBending)(hitNumber1, hitNumber2) =
793           (*covBending)(hitNumber2, hitNumber1);
794         // and non bending plane
795         (*covNonBending)(hitNumber1, hitNumber2) =
796           (*covNonBending)(hitNumber2, hitNumber1);
797       }
798     } // for (hitNumber2 = hitNumber1;...
799   } // for (hitNumber1 = 0;...
800     
801   // Inversion of covariance matrices
802   // with "mnvertLocal", local "mnvert" function of Minuit.
803   // One cannot use directly "mnvert" since "TVirtualFitter" does not know it.
804   // One will have to replace this local function by the right inversion function
805   // from a specialized Root package for symmetric positive definite matrices,
806   // when available!!!!
807   Int_t ifailBending;
808   mnvertLocal(&((*covBending)(0,0)), numberOfHit, numberOfHit, numberOfHit,
809               ifailBending);
810   Int_t ifailNonBending;
811   mnvertLocal(&((*covNonBending)(0,0)), numberOfHit, numberOfHit, numberOfHit,
812               ifailNonBending);
813
814   // It would be worth trying to calculate the inverse of the covariance matrix
815   // only once per fit, since it cannot change much in principle,
816   // and it would save a lot of computing time !!!!
817   
818   // Calculates Chi2
819   if ((ifailBending == 0) && (ifailNonBending == 0)) {
820     // with Multiple Scattering if inversion correct
821     for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) { 
822       hit1 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber1];
823       hbc1 = hit1->GetHitForRecPtr()->GetBendingCoor();
824       pbc1 = hit1->GetTrackParam()->GetBendingCoor();
825       hnbc1 = hit1->GetHitForRecPtr()->GetNonBendingCoor();
826       pnbc1 = hit1->GetTrackParam()->GetNonBendingCoor();
827       for (hitNumber2 = 0; hitNumber2 < numberOfHit; hitNumber2++) {
828         hit2 = (AliMUONTrackHit*)
829           (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber2];
830         hbc2 = hit2->GetHitForRecPtr()->GetBendingCoor();
831         pbc2 = hit2->GetTrackParam()->GetBendingCoor();
832         hnbc2 = hit2->GetHitForRecPtr()->GetNonBendingCoor();
833         pnbc2 = hit2->GetTrackParam()->GetNonBendingCoor();
834         Chi2 = Chi2 +
835           ((*covBending)(hitNumber2, hitNumber1) *
836            (hbc1 - pbc1) * (hbc2 - pbc2)) +
837           ((*covNonBending)(hitNumber2, hitNumber1) *
838            (hnbc1 - pnbc1) * (hnbc2 - pnbc2));
839       }
840     }
841   } else {
842     // without Multiple Scattering if inversion impossible
843     for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) { 
844       hit1 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber1];
845       hbc1 = hit1->GetHitForRecPtr()->GetBendingCoor();
846       pbc1 = hit1->GetTrackParam()->GetBendingCoor();
847       hnbc1 = hit1->GetHitForRecPtr()->GetNonBendingCoor();
848       pnbc1 = hit1->GetTrackParam()->GetNonBendingCoor();
849       Chi2 = Chi2 + 
850         ((hbc1 - pbc1) * (hbc1 - pbc1) /
851          hit1->GetHitForRecPtr()->GetBendingReso2()) +
852         ((hnbc1 - pnbc1) * (hnbc1 - pnbc1) /
853          hit1->GetHitForRecPtr()->GetNonBendingReso2());
854     }
855   }
856   
857   delete covBending;
858   delete covNonBending;
859   delete [] msa2;
860 }
861
862 Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit)
863 {
864   // Returns square of multiple Coulomb scattering angle
865   // at TrackHit pointed to by "TrackHit"
866   Double_t slopeBending, slopeNonBending, radiationLength, inverseBendingMomentum2, inverseTotalMomentum2;
867   Double_t varMultipleScatteringAngle;
868   AliMUONTrack *trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit();
869   AliMUONTrackParam *param = TrackHit->GetTrackParam();
870   // Better implementation in AliMUONTrack class ????
871   slopeBending = param->GetBendingSlope();
872   slopeNonBending = param->GetNonBendingSlope();
873   // thickness in radiation length for the current track,
874   // taking local angle into account
875   radiationLength =
876     trackBeingFitted->GetTrackReconstructor()->GetChamberThicknessInX0() *
877     TMath::Sqrt(1.0 +
878                 slopeBending * slopeBending + slopeNonBending * slopeNonBending);
879   inverseBendingMomentum2 = 
880     param->GetInverseBendingMomentum() * param->GetInverseBendingMomentum();
881   inverseTotalMomentum2 =
882     inverseBendingMomentum2 * (1.0 + slopeBending * slopeBending) /
883     (1.0 + slopeBending *slopeBending + slopeNonBending * slopeNonBending); 
884   varMultipleScatteringAngle = 0.0136 * (1.0 + 0.038 * TMath::Log(radiationLength));
885   // The velocity is assumed to be 1 !!!!
886   varMultipleScatteringAngle = inverseTotalMomentum2 * radiationLength *
887     varMultipleScatteringAngle * varMultipleScatteringAngle;
888   return varMultipleScatteringAngle;
889 }
890
891 //______________________________________________________________________________
892  void mnvertLocal(Double_t *a, Int_t l, Int_t, Int_t n, Int_t &ifail)
893 {
894 //*-*-*-*-*-*-*-*-*-*-*-*Inverts a symmetric matrix*-*-*-*-*-*-*-*-*-*-*-*-*
895 //*-*                    ==========================
896 //*-*        inverts a symmetric matrix.   matrix is first scaled to
897 //*-*        have all ones on the diagonal (equivalent to change of units)
898 //*-*        but no pivoting is done since matrix is positive-definite.
899 //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
900
901   // taken from TMinuit package of Root (l>=n)
902   // fVERTs, fVERTq and fVERTpp changed to localVERTs, localVERTq and localVERTpp
903   //  Double_t localVERTs[n], localVERTq[n], localVERTpp[n];
904   Double_t * localVERTs = new Double_t[n];
905   Double_t * localVERTq = new Double_t[n];
906   Double_t * localVERTpp = new Double_t[n];
907   // fMaxint changed to localMaxint
908   Int_t localMaxint = n;
909
910     /* System generated locals */
911     Int_t aOffset;
912
913     /* Local variables */
914     Double_t si;
915     Int_t i, j, k, kp1, km1;
916
917     /* Parameter adjustments */
918     aOffset = l + 1;
919     a -= aOffset;
920
921     /* Function Body */
922     ifail = 0;
923     if (n < 1) goto L100;
924     if (n > localMaxint) goto L100;
925 //*-*-                  scale matrix by sqrt of diag elements
926     for (i = 1; i <= n; ++i) {
927         si = a[i + i*l];
928         if (si <= 0) goto L100;
929         localVERTs[i-1] = 1 / TMath::Sqrt(si);
930     }
931     for (i = 1; i <= n; ++i) {
932         for (j = 1; j <= n; ++j) {
933             a[i + j*l] = a[i + j*l]*localVERTs[i-1]*localVERTs[j-1];
934         }
935     }
936 //*-*-                                       . . . start main loop . . . .
937     for (i = 1; i <= n; ++i) {
938         k = i;
939 //*-*-                  preparation for elimination step1
940         if (a[k + k*l] != 0) localVERTq[k-1] = 1 / a[k + k*l];
941         else goto L100;
942         localVERTpp[k-1] = 1;
943         a[k + k*l] = 0;
944         kp1 = k + 1;
945         km1 = k - 1;
946         if (km1 < 0) goto L100;
947         else if (km1 == 0) goto L50;
948         else               goto L40;
949 L40:
950         for (j = 1; j <= km1; ++j) {
951             localVERTpp[j-1] = a[j + k*l];
952             localVERTq[j-1]  = a[j + k*l]*localVERTq[k-1];
953             a[j + k*l]   = 0;
954         }
955 L50:
956         if (k - n < 0) goto L51;
957         else if (k - n == 0) goto L60;
958         else                goto L100;
959 L51:
960         for (j = kp1; j <= n; ++j) {
961             localVERTpp[j-1] = a[k + j*l];
962             localVERTq[j-1]  = -a[k + j*l]*localVERTq[k-1];
963             a[k + j*l]   = 0;
964         }
965 //*-*-                  elimination proper
966 L60:
967         for (j = 1; j <= n; ++j) {
968             for (k = j; k <= n; ++k) { a[j + k*l] += localVERTpp[j-1]*localVERTq[k-1]; }
969         }
970     }
971 //*-*-                  elements of left diagonal and unscaling
972     for (j = 1; j <= n; ++j) {
973         for (k = 1; k <= j; ++k) {
974             a[k + j*l] = a[k + j*l]*localVERTs[k-1]*localVERTs[j-1];
975             a[j + k*l] = a[k + j*l];
976         }
977     }
978     delete [] localVERTs;
979     delete [] localVERTq;
980     delete [] localVERTpp;
981     return;
982 //*-*-                  failure return
983 L100:
984     delete [] localVERTs;
985     delete [] localVERTq;
986     delete [] localVERTpp;
987     ifail = 1;
988 } /* mnvertLocal */
989
990 //_____________________________________________-
991 void AliMUONTrack::Print(Option_t* opt) const
992 {
993 //
994   // Printing Track information 
995   // "full" option for printing all the information about the track
996   //
997   TString sopt(opt);
998   sopt.ToUpper();
999  
1000   if ( sopt.Contains("FULL") ) { 
1001     cout << "<AliMUONTrack> No.Clusters=" << setw(2)   << GetNTrackHits() << 
1002       //      ", Bending P="<< setw(8) << setprecision(5)      << 1./GetInverseBendingMomentum() << 
1003       //", NonBendSlope=" << setw(8) << setprecision(5)  << GetNonBendingSlope()*180./TMath::Pi() <<
1004       //", BendSlope=" << setw(8) << setprecision(5)     << GetBendingSlope()*180./TMath::Pi() <<
1005       ", Match2Trig=" << setw(1) << GetMatchTrigger()  << 
1006       ", Chi2-tracking-trigger=" << setw(8) << setprecision(5) <<  GetChi2MatchTrigger() << endl ;
1007     GetTrackParamAtHit()->First()->Print("full");
1008   }
1009   else {
1010     cout << "<AliMUONTrack>";
1011     GetTrackParamAtHit()->First()->Print("");
1012
1013   }
1014     
1015 }