Arrays with variable size created by new and deleted at the end
[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 /*
17 $Log$
18 Revision 1.10  2001/04/09 12:25:09  gosset
19 Inversion of covariance matrices with local copy of TMinuit::mnvert,
20 for symmetric positive definite matrices, instead of TMatrixD::Invert
21
22 Revision 1.9  2001/01/17 20:59:24  hristov
23 chPrev initialised
24
25 Revision 1.8  2001/01/08 11:01:02  gosset
26 Modifications used for addendum to Dimuon TDR (JP Cussonneau):
27 *. MaxBendingMomentum to make both a segment and a track (default 500)
28 *. MaxChi2 per degree of freedom to make a track (default 100)
29 *. MinBendingMomentum used also to make a track
30    and not only a segment (default 3)
31 *. wider roads for track search in stations 1 to 3
32 *. extrapolation to actual Z instead of Z(chamber) in FollowTracks
33 *. in track fit:
34    - limits on parameters X and Y (+/-500)
35    - covariance matrices in double precision
36    - normalization of covariance matrices before inversion
37    - suppression of Minuit printouts
38 *. correction against memory leak (delete extrapHit) in FollowTracks
39 *. RMax to 10 degrees with Z(chamber) instead of fixed values;
40    RMin and Rmax cuts suppressed in NewHitForRecFromGEANT,
41    because useless with realistic geometry
42
43 Revision 1.7  2000/09/19 15:50:46  gosset
44 TrackChi2MCS function: covariance matrix better calculated,
45 taking into account missing planes...
46
47 Revision 1.6  2000/07/20 12:45:27  gosset
48 New "EventReconstructor..." structure,
49         hopefully more adapted to tree/streamer.
50 "AliMUONEventReconstructor::RemoveDoubleTracks"
51         to keep only one track among similar ones.
52
53 Revision 1.5  2000/07/18 16:04:06  gosset
54 AliMUONEventReconstructor package:
55 * a few minor modifications and more comments
56 * a few corrections
57   * right sign for Z of raw clusters
58   * right loop over chambers inside station
59   * symmetrized covariance matrix for measurements (TrackChi2MCS)
60   * right sign of charge in extrapolation (ExtrapToZ)
61   * right zEndAbsorber for Branson correction below 3 degrees
62 * use of TVirtualFitter instead of TMinuit for AliMUONTrack::Fit
63 * no parameter for AliMUONTrack::Fit() but more fit parameters in Track object
64
65 Revision 1.4  2000/06/30 10:15:48  gosset
66 Changes to EventReconstructor...:
67 precision fit with multiple Coulomb scattering;
68 extrapolation to vertex with Branson correction in absorber (JPC)
69
70 Revision 1.3  2000/06/25 13:23:28  hristov
71 stdlib.h needed for non-Linux compilation
72
73 Revision 1.2  2000/06/15 07:58:48  morsch
74 Code from MUON-dev joined
75
76 Revision 1.1.2.3  2000/06/12 10:11:34  morsch
77 Dummy copy constructor and assignment operator added
78
79 Revision 1.1.2.2  2000/06/09 12:58:05  gosset
80 Removed comment beginnings in Log sections of .cxx files
81 Suppressed most violations of coding rules
82
83 Revision 1.1.2.1  2000/06/07 14:44:53  gosset
84 Addition of files for track reconstruction in C++
85 */
86
87 //__________________________________________________________________________
88 //
89 // Reconstructed track in ALICE dimuon spectrometer
90 //__________________________________________________________________________
91
92 #include "AliMUONTrack.h"
93
94 #include <iostream.h>
95
96 #include <TClonesArray.h>
97 #include <TMath.h>
98 #include <TMatrixD.h>
99 #include <TObjArray.h>
100 #include <TVirtualFitter.h>
101
102 #include "AliMUONEventReconstructor.h" 
103 #include "AliMUONHitForRec.h" 
104 #include "AliMUONSegment.h" 
105 #include "AliMUONTrackHit.h"
106
107 #include <stdlib.h>
108
109 // Functions to be minimized with Minuit
110 void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag);
111 void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag);
112
113 void mnvertLocal(Double_t* a, Int_t l, Int_t m, Int_t n, Int_t& ifail);
114
115 Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit);
116
117 ClassImp(AliMUONTrack) // Class implementation in ROOT context
118
119 TVirtualFitter* AliMUONTrack::fgFitter = NULL; 
120
121   //__________________________________________________________________________
122 AliMUONTrack::AliMUONTrack(AliMUONSegment* BegSegment, AliMUONSegment* EndSegment, AliMUONEventReconstructor* EventReconstructor)
123 {
124   // Constructor from two Segment's
125   fEventReconstructor = EventReconstructor; // link back to EventReconstructor
126   // memory allocation for the TObjArray of pointers to reconstructed TrackHit's
127   fTrackHitsPtr = new TObjArray(10);
128   fNTrackHits = 0;
129   AddSegment(BegSegment); // add hits from BegSegment
130   AddSegment(EndSegment); // add hits from EndSegment
131   fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z
132   SetTrackParamAtVertex(); // set track parameters at vertex
133   // set fit conditions...
134   fFitMCS = 0;
135   fFitNParam = 3;
136   fFitStart = 1;
137   fFitFMin = -1.0;
138   return;
139 }
140
141   //__________________________________________________________________________
142 AliMUONTrack::AliMUONTrack(AliMUONSegment* Segment, AliMUONHitForRec* HitForRec, AliMUONEventReconstructor* EventReconstructor)
143 {
144   // Constructor from one Segment and one HitForRec
145   fEventReconstructor = EventReconstructor; // link back to EventReconstructor
146   // memory allocation for the TObjArray of pointers to reconstructed TrackHit's
147   fTrackHitsPtr = new TObjArray(10);
148   fNTrackHits = 0;
149   AddSegment(Segment); // add hits from Segment
150   AddHitForRec(HitForRec); // add HitForRec
151   fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z
152   SetTrackParamAtVertex(); // set track parameters at vertex
153   // set fit conditions...
154   fFitMCS = 0;
155   fFitNParam = 3;
156   fFitStart = 1;
157   fFitFMin = -1.0;
158   return;
159 }
160
161   //__________________________________________________________________________
162 AliMUONTrack::~AliMUONTrack()
163 {
164   // Destructor
165   if (fTrackHitsPtr) {
166     delete fTrackHitsPtr; // delete the TObjArray of pointers to TrackHit's
167     fTrackHitsPtr = NULL;
168   }
169 }
170
171   //__________________________________________________________________________
172 AliMUONTrack::AliMUONTrack (const AliMUONTrack& MUONTrack)
173 {
174 // Dummy copy constructor
175 }
176
177   //__________________________________________________________________________
178 AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& MUONTrack)
179 {
180 // Dummy assignment operator
181     return *this;
182 }
183
184   //__________________________________________________________________________
185 void AliMUONTrack::Remove()
186 {
187   // Remove current track from array of tracks,
188   // and corresponding track hits from array of track hits.
189   // Compress the TClonesArray it belongs to.
190   AliMUONTrackHit *nextTrackHit;
191   AliMUONEventReconstructor *eventRec = this->fEventReconstructor;
192   TClonesArray *trackHitsPtr = eventRec->GetRecTrackHitsPtr();
193   // Loop over all track hits of track
194   AliMUONTrackHit *trackHit = (AliMUONTrackHit*) fTrackHitsPtr->First();
195   while (trackHit) {
196     nextTrackHit = (AliMUONTrackHit*) fTrackHitsPtr->After(trackHit);
197     // Remove TrackHit from event TClonesArray.
198     // Destructor is called,
199     // hence links between HitForRec's and TrackHit's are updated
200     trackHitsPtr->Remove(trackHit);
201     trackHit = nextTrackHit;
202   }
203   // Remove the track from event TClonesArray
204   // Destructor is called,
205   // hence space for TObjArray of pointers to TrackHit's is freed
206   eventRec->GetRecTracksPtr()->Remove(this);
207   // Number of tracks decreased by 1
208   eventRec->SetNRecTracks(eventRec->GetNRecTracks() - 1);
209   // Compress event TClonesArray of Track's:
210   // this is essential to retrieve the TClonesArray afterwards
211   eventRec->GetRecTracksPtr()->Compress();
212   // Compress event TClonesArray of TrackHit's:
213   // this is probably also essential to retrieve the TClonesArray afterwards
214   trackHitsPtr->Compress();
215 }
216
217   //__________________________________________________________________________
218 void AliMUONTrack::SetFitMCS(Int_t FitMCS)
219 {
220   // Set multiple Coulomb scattering option for track fit "fFitMCS"
221   // from "FitMCS" argument: 0 without, 1 with
222   if ((FitMCS == 0) || (FitMCS == 1)) fFitMCS = FitMCS;
223   // better implementation with enum(with, without) ????
224   else {
225     cout << "ERROR in AliMUONTrack::SetFitMCS(FitMCS)" << endl;
226     cout << "FitMCS = " << FitMCS << " is neither 0 nor 1" << endl;
227     exit(0);
228   }
229   return;
230 }
231
232   //__________________________________________________________________________
233 void AliMUONTrack::SetFitNParam(Int_t FitNParam)
234 {
235   // Set number of parameters for track fit "fFitNParam" from "FitNParam":
236   // 3 for momentum, 5 for momentum and position
237   if ((FitNParam == 3) || (FitNParam == 5)) fFitNParam = FitNParam;
238   else {
239     cout << "ERROR in AliMUONTrack::SetFitNParam(FitNParam)" << endl;
240     cout << "FitNParam = " << FitNParam << " is neither 3 nor 5" << endl;
241     exit(0);
242   }
243   return;
244 }
245
246   //__________________________________________________________________________
247 void AliMUONTrack::SetFitStart(Int_t FitStart)
248 {
249   // Set multiple Coulomb scattering option for track fit "fFitStart"
250   // from "FitStart" argument: 0 without, 1 with
251   if ((FitStart == 0) || (FitStart == 1)) fFitStart = FitStart;
252   // better implementation with enum(vertex, firstHit) ????
253   else {
254     cout << "ERROR in AliMUONTrack::SetFitStart(FitStart)" << endl;
255     cout << "FitStart = " << FitStart << " is neither 0 nor 1" << endl;
256     exit(0);
257   }
258   return;
259 }
260
261   //__________________________________________________________________________
262 AliMUONTrackParam* AliMUONTrack::GetTrackParamAtFirstHit(void) {
263   // Get pointer to TrackParamAtFirstHit
264   return ((AliMUONTrackHit*) (fTrackHitsPtr->First()))->GetTrackParam();}
265
266   //__________________________________________________________________________
267 void AliMUONTrack::RecursiveDump(void)
268 {
269   // Recursive dump of AliMUONTrack, i.e. with dump of TrackHit's and HitForRec's
270   AliMUONTrackHit *trackHit;
271   AliMUONHitForRec *hitForRec;
272   cout << "Recursive dump of Track: " << this << endl;
273   // Track
274   this->Dump();
275   for (Int_t trackHitIndex = 0; trackHitIndex < fNTrackHits; trackHitIndex++) {
276     trackHit = (AliMUONTrackHit*) ((*fTrackHitsPtr)[trackHitIndex]);
277     // TrackHit
278     cout << "TrackHit: " << trackHit << " (index: " << trackHitIndex << ")" << endl;
279     trackHit->Dump();
280     hitForRec = trackHit->GetHitForRecPtr();
281     // HitForRec
282     cout << "HitForRec: " << hitForRec << endl;
283     hitForRec->Dump();
284   }
285   return;
286 }
287
288   //__________________________________________________________________________
289 Int_t AliMUONTrack::HitsInCommon(AliMUONTrack* Track)
290 {
291   // Returns the number of hits in common
292   // between the current track ("this")
293   // and the track pointed to by "Track".
294   Int_t hitsInCommon = 0;
295   AliMUONTrackHit *trackHit1, *trackHit2;
296   // Loop over hits of first track
297   trackHit1 = (AliMUONTrackHit*) this->GetTrackHitsPtr()->First();
298   while (trackHit1) {
299     // Loop over hits of second track
300     trackHit2 = (AliMUONTrackHit*) Track->GetTrackHitsPtr()->First();
301     while (trackHit2) {
302       // Increment "hitsInCommon" if both TrackHits point to the same HitForRec
303       if ( (trackHit1->GetHitForRecPtr()) ==
304            (trackHit2->GetHitForRecPtr())    ) hitsInCommon++;
305       trackHit2 = (AliMUONTrackHit*) Track->GetTrackHitsPtr()->After(trackHit2);
306     } // trackHit2
307     trackHit1 = (AliMUONTrackHit*) this->GetTrackHitsPtr()->After(trackHit1);
308   } // trackHit1
309   return hitsInCommon;
310 }
311
312   //__________________________________________________________________________
313 void AliMUONTrack::Fit()
314 {
315   // Fit the current track ("this"),
316   // with or without multiple Coulomb scattering according to "fFitMCS",
317   // with the number of parameters given by "fFitNParam"
318   // (3 if one keeps X and Y fixed in "TrackParam", 5 if one lets them vary),
319   // starting, according to "fFitStart",
320   // with track parameters at vertex or at the first TrackHit.
321   // "fFitMCS", "fFitNParam" and "fFitStart" have to be set before
322   // by calling the corresponding Set methods.
323   Double_t arg[1], benC, errorParam, invBenP, lower, nonBenC, upper, x, y;
324   char parName[50];
325   AliMUONTrackParam *trackParam;
326   // Check if Minuit is initialized...
327   fgFitter = TVirtualFitter::Fitter(this); // add 3 or 5 for the maximum number of parameters ???
328   fgFitter->Clear(); // necessary ???? probably yes
329   // how to reset the printout number at every fit ????
330   // is there any risk to leave it like that ????
331   // how to go faster ???? choice of Minuit parameters like EDM ????
332   // choice of function to be minimized according to fFitMCS
333   if (fFitMCS == 0) fgFitter->SetFCN(TrackChi2);
334   else fgFitter->SetFCN(TrackChi2MCS);
335   arg[0] = -1;
336   fgFitter->ExecuteCommand("SET PRINT", arg, 1); // More printing !!!!
337   // Parameters according to "fFitStart"
338   // (should be a function to be used at every place where needed ????)
339   if (fFitStart == 0) trackParam = &fTrackParamAtVertex;
340   else trackParam = this->GetTrackParamAtFirstHit();
341   // set first 3 Minuit parameters
342   // could be tried with no limits for the search (min=max=0) ????
343   fgFitter->SetParameter(0, "InvBenP",
344                          trackParam->GetInverseBendingMomentum(),
345                          0.003, -0.4, 0.4);
346   fgFitter->SetParameter(1, "BenS",
347                          trackParam->GetBendingSlope(),
348                          0.001, -0.5, 0.5);
349   fgFitter->SetParameter(2, "NonBenS",
350                          trackParam->GetNonBendingSlope(),
351                          0.001, -0.5, 0.5);
352   if (fFitNParam == 5) {
353     // set last 2 Minuit parameters
354     // mandatory limits in Bending to avoid NaN values of parameters
355     fgFitter->SetParameter(3, "X",
356                            trackParam->GetNonBendingCoor(),
357                            0.03, -500.0, 500.0);
358     // mandatory limits in non Bending to avoid NaN values of parameters
359     fgFitter->SetParameter(4, "Y",
360                            trackParam->GetBendingCoor(),
361                            0.10, -500.0, 500.0);
362   }
363   // search without gradient calculation in the function
364   fgFitter->ExecuteCommand("SET NOGRADIENT", arg, 0);
365   // minimization
366   fgFitter->ExecuteCommand("MINIMIZE", arg, 0);
367   // exit from Minuit
368   fgFitter->ExecuteCommand("EXIT", arg, 0); // necessary ????
369   // get results into "invBenP", "benC", "nonBenC" ("x", "y")
370   fgFitter->GetParameter(0, parName, invBenP, errorParam, lower, upper);
371   fgFitter->GetParameter(1, parName, benC, errorParam, lower, upper);
372   fgFitter->GetParameter(2, parName, nonBenC, errorParam, lower, upper);
373   if (fFitNParam == 5) {
374     fgFitter->GetParameter(3, parName, x, errorParam, lower, upper);
375     fgFitter->GetParameter(4, parName, y, errorParam, lower, upper);
376   }
377   // result of the fit into track parameters
378   trackParam->SetInverseBendingMomentum(invBenP);
379   trackParam->SetBendingSlope(benC);
380   trackParam->SetNonBendingSlope(nonBenC);
381   if (fFitNParam == 5) {
382     trackParam->SetNonBendingCoor(x);
383     trackParam->SetBendingCoor(y);
384   }
385   // global result of the fit
386   Double_t fedm, errdef;
387   Int_t npari, nparx;
388   fgFitter->GetStats(fFitFMin, fedm, errdef, npari, nparx);
389 }
390
391   //__________________________________________________________________________
392 void AliMUONTrack::AddSegment(AliMUONSegment* Segment)
393 {
394   // Add Segment to the track
395   AddHitForRec(Segment->GetHitForRec1()); // 1st hit
396   AddHitForRec(Segment->GetHitForRec2()); // 2nd hit
397 }
398
399   //__________________________________________________________________________
400 void AliMUONTrack::AddHitForRec(AliMUONHitForRec* HitForRec)
401 {
402   // Add HitForRec to the track:
403   // actual TrackHit into TClonesArray of TrackHit's for the event;
404   // pointer to actual TrackHit in TObjArray of pointers to TrackHit's for the track
405   TClonesArray *recTrackHitsPtr = this->fEventReconstructor->GetRecTrackHitsPtr();
406   Int_t eventTrackHits = this->fEventReconstructor->GetNRecTrackHits();
407   // event
408   AliMUONTrackHit* trackHit =
409     new ((*recTrackHitsPtr)[eventTrackHits]) AliMUONTrackHit(HitForRec);
410   this->fEventReconstructor->SetNRecTrackHits(eventTrackHits + 1);
411   // track
412   fTrackHitsPtr->Add(trackHit);
413   fNTrackHits++;
414 }
415
416   //__________________________________________________________________________
417 void AliMUONTrack::SetTrackParamAtHit(Int_t indexHit, AliMUONTrackParam *TrackParam)
418 {
419   // Set track parameters at TrackHit with index "indexHit"
420   // from the track parameters pointed to by "TrackParam".
421   AliMUONTrackHit* trackHit = (AliMUONTrackHit*) ((*fTrackHitsPtr)[indexHit]);
422   trackHit->SetTrackParam(TrackParam);
423 }
424
425   //__________________________________________________________________________
426 void AliMUONTrack::SetTrackParamAtVertex()
427 {
428   // Set track parameters at vertex.
429   // TrackHit's are assumed to be only in stations(1..) 4 and 5,
430   // and sorted according to increasing Z..
431   // Parameters are calculated from information in HitForRec's
432   // of first and last TrackHit's.
433   AliMUONTrackParam *trackParam =
434     &fTrackParamAtVertex; // pointer to track parameters
435   // Pointer to HitForRec of first TrackHit
436   AliMUONHitForRec *firstHit =
437     ((AliMUONTrackHit*) (fTrackHitsPtr->First()))->GetHitForRecPtr();
438   // Pointer to HitForRec of last TrackHit
439   AliMUONHitForRec *lastHit =
440     ((AliMUONTrackHit*) (fTrackHitsPtr->Last()))->GetHitForRecPtr();
441   // Z difference between first and last hits
442   Double_t deltaZ = firstHit->GetZ() - lastHit->GetZ();
443   // bending slope in stations(1..) 4 and 5
444   Double_t bendingSlope =
445     (firstHit->GetBendingCoor() - lastHit->GetBendingCoor()) / deltaZ;
446   trackParam->SetBendingSlope(bendingSlope);
447   // impact parameter
448   Double_t impactParam =
449     firstHit->GetBendingCoor() - bendingSlope * firstHit->GetZ(); // same if from firstHit and  lastHit ????
450   // signed bending momentum
451   Double_t signedBendingMomentum =
452     fEventReconstructor->GetBendingMomentumFromImpactParam(impactParam);
453   trackParam->SetInverseBendingMomentum(1.0 / signedBendingMomentum);
454   // bending slope at vertex
455   trackParam->
456     SetBendingSlope(bendingSlope +
457                     impactParam / fEventReconstructor->GetSimpleBPosition());
458   // non bending slope
459   Double_t nonBendingSlope =
460     (firstHit->GetNonBendingCoor() - lastHit->GetNonBendingCoor()) / deltaZ;
461   trackParam->SetNonBendingSlope(nonBendingSlope);
462   // vertex coordinates at (0,0,0)
463   trackParam->SetZ(0.0);
464   trackParam->SetBendingCoor(0.0);
465   trackParam->SetNonBendingCoor(0.0);
466 }
467
468   //__________________________________________________________________________
469 void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag)
470 {
471   // Return the "Chi2" to be minimized with Minuit for track fitting,
472   // with "NParam" parameters
473   // and their current values in array pointed to by "Param".
474   // Assumes that the track hits are sorted according to increasing Z.
475   // Track parameters at each TrackHit are updated accordingly.
476   // Multiple Coulomb scattering is not taken into account
477   AliMUONTrack *trackBeingFitted;
478   AliMUONTrackHit* hit;
479   AliMUONTrackParam param1;
480   Int_t hitNumber;
481   Double_t zHit;
482   Chi2 = 0.0; // initialize Chi2
483   // copy of track parameters to be fitted
484   trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit();
485   if (trackBeingFitted->GetFitStart() == 0)
486     param1 = *(trackBeingFitted->GetTrackParamAtVertex());
487   else param1 = *(trackBeingFitted->GetTrackParamAtFirstHit());
488   // Minuit parameters to be fitted into this copy
489   param1.SetInverseBendingMomentum(Param[0]);
490   param1.SetBendingSlope(Param[1]);
491   param1.SetNonBendingSlope(Param[2]);
492   if (NParam == 5) {
493     param1.SetNonBendingCoor(Param[3]);
494     param1.SetBendingCoor(Param[4]);
495   }
496   // Follow track through all planes of track hits
497   for (hitNumber = 0; hitNumber < trackBeingFitted->GetNTrackHits(); hitNumber++) {
498     hit = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber];
499     zHit = hit->GetHitForRecPtr()->GetZ();
500     // do something special if 2 hits with same Z ????
501     // security against infinite loop ????
502     (&param1)->ExtrapToZ(zHit); // extrapolation
503     hit->SetTrackParam(&param1);
504     // Increment Chi2
505     // done hit per hit, with hit resolution,
506     // and not with point and angle like in "reco_muon.F" !!!!
507     // Needs to add multiple scattering contribution ????
508     Double_t dX =
509       hit->GetHitForRecPtr()->GetNonBendingCoor() - (&param1)->GetNonBendingCoor();
510     Double_t dY =
511       hit->GetHitForRecPtr()->GetBendingCoor() - (&param1)->GetBendingCoor();
512     Chi2 =
513       Chi2 +
514       dX * dX / hit->GetHitForRecPtr()->GetNonBendingReso2() +
515       dY * dY / hit->GetHitForRecPtr()->GetBendingReso2();
516   }
517 }
518
519   //__________________________________________________________________________
520 void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag)
521 {
522   // Return the "Chi2" to be minimized with Minuit for track fitting,
523   // with "NParam" parameters
524   // and their current values in array pointed to by "Param".
525   // Assumes that the track hits are sorted according to increasing Z.
526   // Track parameters at each TrackHit are updated accordingly.
527   // Multiple Coulomb scattering is taken into account with covariance matrix.
528   AliMUONTrack *trackBeingFitted;
529   AliMUONTrackParam param1;
530   Chi2 = 0.0; // initialize Chi2
531   // copy of track parameters to be fitted
532   trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit();
533   if (trackBeingFitted->GetFitStart() == 0)
534     param1 = *(trackBeingFitted->GetTrackParamAtVertex());
535   else param1 = *(trackBeingFitted->GetTrackParamAtFirstHit());
536   // Minuit parameters to be fitted into this copy
537   param1.SetInverseBendingMomentum(Param[0]);
538   param1.SetBendingSlope(Param[1]);
539   param1.SetNonBendingSlope(Param[2]);
540   if (NParam == 5) {
541     param1.SetNonBendingCoor(Param[3]);
542     param1.SetBendingCoor(Param[4]);
543   }
544
545   AliMUONTrackHit *hit;
546   Int_t chCurrent, chPrev = 0, hitNumber, hitNumber1, hitNumber2, hitNumber3;
547   Double_t z, z1, z2, z3;
548   AliMUONTrackHit *hit1, *hit2, *hit3;
549   Double_t hbc1, hbc2, pbc1, pbc2;
550   Double_t hnbc1, hnbc2, pnbc1, pnbc2;
551   Int_t numberOfHit = trackBeingFitted->GetNTrackHits();
552   TMatrixD *covBending = new TMatrixD(numberOfHit, numberOfHit);
553   TMatrixD *covNonBending = new TMatrixD(numberOfHit, numberOfHit);
554   Double_t *msa2 = new Double_t[numberOfHit];
555
556   // Predicted coordinates and  multiple scattering angles are first calculated
557   for (hitNumber = 0; hitNumber < numberOfHit; hitNumber++) {
558     hit = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber];
559     z = hit->GetHitForRecPtr()->GetZ();
560     // do something special if 2 hits with same Z ????
561     // security against infinite loop ????
562     (&param1)->ExtrapToZ(z); // extrapolation
563     hit->SetTrackParam(&param1);
564     // square of multiple scattering angle at current hit, with one chamber
565     msa2[hitNumber] = MultipleScatteringAngle2(hit);
566     // correction for eventual missing hits or multiple hits in a chamber,
567     // according to the number of chambers
568     // between the current hit and the previous one
569     chCurrent = hit->GetHitForRecPtr()->GetChamberNumber();
570     if (hitNumber > 0) msa2[hitNumber] = msa2[hitNumber] * (chCurrent - chPrev);
571     chPrev = chCurrent;
572   }
573
574   // Calculates the covariance matrix
575   for (hitNumber1 = 0; hitNumber1 < numberOfHit; hitNumber1++) { 
576     hit1 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber1];
577     z1 = hit1->GetHitForRecPtr()->GetZ();
578     for (hitNumber2 = hitNumber1; hitNumber2 < numberOfHit; hitNumber2++) {
579       hit2 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber2];
580       z2 = hit2->GetHitForRecPtr()->GetZ();
581       // initialization to 0 (diagonal plus upper triangular part)
582       (*covBending)(hitNumber2, hitNumber1) = 0.0;
583       // contribution from multiple scattering in bending plane:
584       // loop over upstream hits
585       for (hitNumber3 = 0; hitNumber3 < hitNumber1; hitNumber3++) {     
586         hit3 = (AliMUONTrackHit*)
587           (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber3];
588         z3 = hit3->GetHitForRecPtr()->GetZ();
589         (*covBending)(hitNumber2, hitNumber1) =
590           (*covBending)(hitNumber2, hitNumber1) +
591           ((z1 - z3) * (z2 - z3) * msa2[hitNumber3]); 
592       }
593       // equal contribution from multiple scattering in non bending plane
594       (*covNonBending)(hitNumber2, hitNumber1) =
595         (*covBending)(hitNumber2, hitNumber1);
596       if (hitNumber1 == hitNumber2) {
597         // Diagonal elements: add contribution from position measurements
598         // in bending plane
599         (*covBending)(hitNumber2, hitNumber1) =
600           (*covBending)(hitNumber2, hitNumber1) +
601           hit1->GetHitForRecPtr()->GetBendingReso2();
602         // and in non bending plane
603         (*covNonBending)(hitNumber2, hitNumber1) =
604           (*covNonBending)(hitNumber2, hitNumber1) +
605           hit1->GetHitForRecPtr()->GetNonBendingReso2();
606       }
607       else {
608         // Non diagonal elements: symmetrization
609         // for bending plane
610         (*covBending)(hitNumber1, hitNumber2) =
611           (*covBending)(hitNumber2, hitNumber1);
612         // and non bending plane
613         (*covNonBending)(hitNumber1, hitNumber2) =
614           (*covNonBending)(hitNumber2, hitNumber1);
615       }
616     } // for (hitNumber2 = hitNumber1;...
617   } // for (hitNumber1 = 0;...
618     
619   // Inversion of covariance matrices
620   // with "mnvertLocal", local "mnvert" function of Minuit.
621   // One cannot use directly "mnvert" since "TVirtualFitter" does not know it.
622   // One will have to replace this local function by the right inversion function
623   // from a specialized Root package for symmetric positive definite matrices,
624   // when available!!!!
625   Int_t ifailBending;
626   mnvertLocal(&((*covBending)(0,0)), numberOfHit, numberOfHit, numberOfHit,
627               ifailBending);
628   Int_t ifailNonBending;
629   mnvertLocal(&((*covNonBending)(0,0)), numberOfHit, numberOfHit, numberOfHit,
630               ifailNonBending);
631
632   // It would be worth trying to calculate the inverse of the covariance matrix
633   // only once per fit, since it cannot change much in principle,
634   // and it would save a lot of computing time !!!!
635   
636   // Calculates Chi2
637   if ((ifailBending == 0) && (ifailNonBending == 0)) {
638     // with Multiple Scattering if inversion correct
639     for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) { 
640       hit1 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber1];
641       hbc1 = hit1->GetHitForRecPtr()->GetBendingCoor();
642       pbc1 = hit1->GetTrackParam()->GetBendingCoor();
643       hnbc1 = hit1->GetHitForRecPtr()->GetNonBendingCoor();
644       pnbc1 = hit1->GetTrackParam()->GetNonBendingCoor();
645       for (hitNumber2 = 0; hitNumber2 < numberOfHit; hitNumber2++) {
646         hit2 = (AliMUONTrackHit*)
647           (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber2];
648         hbc2 = hit2->GetHitForRecPtr()->GetBendingCoor();
649         pbc2 = hit2->GetTrackParam()->GetBendingCoor();
650         hnbc2 = hit2->GetHitForRecPtr()->GetNonBendingCoor();
651         pnbc2 = hit2->GetTrackParam()->GetNonBendingCoor();
652         Chi2 = Chi2 +
653           ((*covBending)(hitNumber2, hitNumber1) *
654            (hbc1 - pbc1) * (hbc2 - pbc2)) +
655           ((*covNonBending)(hitNumber2, hitNumber1) *
656            (hnbc1 - pnbc1) * (hnbc2 - pnbc2));
657       }
658     }
659   } else {
660     // without Multiple Scattering if inversion impossible
661     for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) { 
662       hit1 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber1];
663       hbc1 = hit1->GetHitForRecPtr()->GetBendingCoor();
664       pbc1 = hit1->GetTrackParam()->GetBendingCoor();
665       hnbc1 = hit1->GetHitForRecPtr()->GetNonBendingCoor();
666       pnbc1 = hit1->GetTrackParam()->GetNonBendingCoor();
667       Chi2 = Chi2 + 
668         ((hbc1 - pbc1) * (hbc1 - pbc1) /
669          hit1->GetHitForRecPtr()->GetBendingReso2()) +
670         ((hnbc1 - pnbc1) * (hnbc1 - pnbc1) /
671          hit1->GetHitForRecPtr()->GetNonBendingReso2());
672     }
673   }
674   
675   delete covBending;
676   delete covNonBending;
677   delete [] msa2;
678 }
679
680 Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit)
681 {
682   // Returns square of multiple Coulomb scattering angle
683   // at TrackHit pointed to by "TrackHit"
684   Double_t slopeBending, slopeNonBending, radiationLength, inverseBendingMomentum2, inverseTotalMomentum2;
685   Double_t varMultipleScatteringAngle;
686   AliMUONTrack *trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit();
687   AliMUONTrackParam *param = TrackHit->GetTrackParam();
688   // Better implementation in AliMUONTrack class ????
689   slopeBending = param->GetBendingSlope();
690   slopeNonBending = param->GetNonBendingSlope();
691   // thickness in radiation length for the current track,
692   // taking local angle into account
693   radiationLength =
694     trackBeingFitted->GetEventReconstructor()->GetChamberThicknessInX0() *
695     TMath::Sqrt(1.0 +
696                 slopeBending * slopeBending + slopeNonBending * slopeNonBending);
697   inverseBendingMomentum2 = 
698     param->GetInverseBendingMomentum() * param->GetInverseBendingMomentum();
699   inverseTotalMomentum2 =
700     inverseBendingMomentum2 * (1.0 + slopeBending * slopeBending) /
701     (1.0 + slopeBending *slopeBending + slopeNonBending * slopeNonBending); 
702   varMultipleScatteringAngle = 0.0136 * (1.0 + 0.038 * TMath::Log(radiationLength));
703   // The velocity is assumed to be 1 !!!!
704   varMultipleScatteringAngle = inverseTotalMomentum2 * radiationLength *
705     varMultipleScatteringAngle * varMultipleScatteringAngle;
706   return varMultipleScatteringAngle;
707 }
708
709 //______________________________________________________________________________
710  void mnvertLocal(Double_t *a, Int_t l, Int_t, Int_t n, Int_t &ifail)
711 {
712 //*-*-*-*-*-*-*-*-*-*-*-*Inverts a symmetric matrix*-*-*-*-*-*-*-*-*-*-*-*-*
713 //*-*                    ==========================
714 //*-*        inverts a symmetric matrix.   matrix is first scaled to
715 //*-*        have all ones on the diagonal (equivalent to change of units)
716 //*-*        but no pivoting is done since matrix is positive-definite.
717 //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
718
719   // taken from TMinuit package of Root (l>=n)
720   // fVERTs, fVERTq and fVERTpp changed to localVERTs, localVERTq and localVERTpp
721   //  Double_t localVERTs[n], localVERTq[n], localVERTpp[n];
722   Double_t * localVERTs = new Double_t[n];
723   Double_t * localVERTq = new Double_t[n];
724   Double_t * localVERTpp = new Double_t[n];
725   // fMaxint changed to localMaxint
726   Int_t localMaxint = n;
727
728     /* System generated locals */
729     Int_t a_offset;
730
731     /* Local variables */
732     Double_t si;
733     Int_t i, j, k, kp1, km1;
734
735     /* Parameter adjustments */
736     a_offset = l + 1;
737     a -= a_offset;
738
739     /* Function Body */
740     ifail = 0;
741     if (n < 1) goto L100;
742     if (n > localMaxint) goto L100;
743 //*-*-                  scale matrix by sqrt of diag elements
744     for (i = 1; i <= n; ++i) {
745         si = a[i + i*l];
746         if (si <= 0) goto L100;
747         localVERTs[i-1] = 1 / TMath::Sqrt(si);
748     }
749     for (i = 1; i <= n; ++i) {
750         for (j = 1; j <= n; ++j) {
751             a[i + j*l] = a[i + j*l]*localVERTs[i-1]*localVERTs[j-1];
752         }
753     }
754 //*-*-                                       . . . start main loop . . . .
755     for (i = 1; i <= n; ++i) {
756         k = i;
757 //*-*-                  preparation for elimination step1
758         if (a[k + k*l] != 0) localVERTq[k-1] = 1 / a[k + k*l];
759         else goto L100;
760         localVERTpp[k-1] = 1;
761         a[k + k*l] = 0;
762         kp1 = k + 1;
763         km1 = k - 1;
764         if (km1 < 0) goto L100;
765         else if (km1 == 0) goto L50;
766         else               goto L40;
767 L40:
768         for (j = 1; j <= km1; ++j) {
769             localVERTpp[j-1] = a[j + k*l];
770             localVERTq[j-1]  = a[j + k*l]*localVERTq[k-1];
771             a[j + k*l]   = 0;
772         }
773 L50:
774         if (k - n < 0) goto L51;
775         else if (k - n == 0) goto L60;
776         else                goto L100;
777 L51:
778         for (j = kp1; j <= n; ++j) {
779             localVERTpp[j-1] = a[k + j*l];
780             localVERTq[j-1]  = -a[k + j*l]*localVERTq[k-1];
781             a[k + j*l]   = 0;
782         }
783 //*-*-                  elimination proper
784 L60:
785         for (j = 1; j <= n; ++j) {
786             for (k = j; k <= n; ++k) { a[j + k*l] += localVERTpp[j-1]*localVERTq[k-1]; }
787         }
788     }
789 //*-*-                  elements of left diagonal and unscaling
790     for (j = 1; j <= n; ++j) {
791         for (k = 1; k <= j; ++k) {
792             a[k + j*l] = a[k + j*l]*localVERTs[k-1]*localVERTs[j-1];
793             a[j + k*l] = a[k + j*l];
794         }
795     }
796     delete localVERTs;
797     delete localVERTq;
798     delete localVERTpp;
799     return;
800 //*-*-                  failure return
801 L100:
802     delete localVERTs;
803     delete localVERTq;
804     delete localVERTpp;
805     ifail = 1;
806 } /* mnvertLocal */
807