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