Setting of aliases to rawReader done only once. Base decision on cosmic or calib...
[u/mrichter/AliRoot.git] / STEER / STEER / AliMillepede.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 //-----------------------------------------------------------------------------
19 /// \class AliMillepede
20 /// Detector independent alignment class
21 ///
22 /// This modified C++ version is based on a C++ translation of Millepede used 
23 /// for LHCb Vertex Detector alignment (lhcb-2005-101), available here 
24 /// http://isscvs.cern.ch/cgi-bin/cvsweb.cgi/Alignment/AlignmentTools/src/?cvsroot=lhcb
25 /// The original millepede fortran package is available at:  
26 /// http://www.desy.de/~blobel/wwwmille.html  
27 ///
28 /// \author Javier Castillo
29 //-----------------------------------------------------------------------------
30
31 #include <TArrayI.h>
32 #include <TArrayD.h>
33 #include <TMath.h>
34
35 #include "AliLog.h"
36
37 #include "AliMillepede.h"
38
39 //=============================================================================
40 AliMillepede::AliMillepede()
41   : TObject(),
42     fIndexLocEq(fgkMaxGlobalPar+fgkMaxLocalPar),
43     fDerivLocEq(fgkMaxGlobalPar+fgkMaxLocalPar),
44     fIndexAllEqs(1000*fgkMaxGlobalPar+fgkMaxLocalPar),
45     fDerivAllEqs(1000*fgkMaxGlobalPar+fgkMaxLocalPar),
46     fLocEqPlace(1000), 
47     fNIndexLocEq(0),
48     fNDerivLocEq(0),
49     fNIndexAllEqs(0),
50     fNDerivAllEqs(0),
51     fNLocEqPlace(0),
52     fNLocalEquations(0),
53     fResCutInit(0.),
54     fResCut(0.),
55     fChi2CutFactor(1.0),
56     fChi2CutRef(1.0),
57     fIter(0),
58     fMaxIter(10),
59     fNStdDev(3),
60     fNGlobalConstraints(0),
61     fNLocalFits(0),
62     fNLocalFitsRejected(0),
63     fNGlobalPar(0),
64     fNLocalPar(0)
65 {
66   /// Standard constructor
67
68   AliInfo("                                           ");
69   AliInfo("            * o o                   o      ");
70   AliInfo("              o o                   o      ");
71   AliInfo("   o ooooo  o o o  oo  ooo   oo   ooo  oo  ");
72   AliInfo("    o  o  o o o o o  o o  o o  o o  o o  o ");
73   AliInfo("    o  o  o o o o oooo o  o oooo o  o oooo ");
74   AliInfo("    o  o  o o o o o    ooo  o    o  o o    ");
75   AliInfo("    o  o  o o o o  oo  o     oo   ooo  oo  ++ starts");         
76   AliInfo("                                           ");
77
78 }
79
80 //=============================================================================
81 AliMillepede::~AliMillepede() {
82   /// Destructor
83
84 }
85
86 //=============================================================================
87 Int_t AliMillepede::InitMille(int nGlo, int nLoc, int lNStdDev,
88                               double lResCut, double lResCutInit)
89 {
90   /// Initialization of millepede
91   AliDebug(1,"");
92   AliDebug(1,"----------------------------------------------------");
93   AliDebug(1,"");
94   AliDebug(1,"    Entering InitMille");
95   AliDebug(1,"");
96   AliDebug(1,"-----------------------------------------------------");
97   AliDebug(1,"");
98
99   fNGlobalConstraints = 0;
100   fNLocalFits  = 0;                      // Total number of local fits
101   fNLocalFitsRejected  = 0;              // Total number of local fits rejected
102   fChi2CutRef  = 1.0;                    // Reference value for Chi^2/ndof cut
103
104   AliMillepede::SetNLocalEquations(0);     // Number of local fits (starts at 0)
105
106   fIter    = 0;  // By default iterations are turned off, turned on by use SetIterations
107   fMaxIter = 10;
108
109   fResCutInit = lResCutInit; 
110   fResCut = lResCut;
111
112   fNGlobalPar = nGlo;     // Number of global derivatives
113   fNLocalPar  = nLoc;     // Number of local derivatives
114   fNStdDev    = lNStdDev; // Number of StDev for local fit chisquare cut
115
116   AliDebug(1,Form("Number of global parameters   : %d", fNGlobalPar));
117   AliDebug(1,Form("Number of local parameters    : %d", fNLocalPar));
118   AliDebug(1,Form("Number of standard deviations : %d", fNStdDev));
119
120   if (fNGlobalPar>fgkMaxGlobalPar || fNLocalPar>fgkMaxLocalPar) {
121     AliDebug(1,"Two many parameters !!!!!");
122     return 0;
123   }
124
125   // Global parameters initializations
126   for (int i=0; i<fNGlobalPar; i++) {  
127     fVecBGlo[i]=0.;
128     fInitPar[i]=0.;
129     fDeltaPar[i]=0.;
130     fSigmaPar[i]=-1.;
131     fIsNonLinear[i]=0;
132     fGlo2CGLRow[i]=-1;
133     fCGLRow2Glo[i]=-1;
134     
135     for (int j=0; j<fNGlobalPar;j++) {    
136       fMatCGlo[i][j]=0.;
137     }
138   }
139
140   // Local parameters initializations
141   
142   for (int i=0; i<fNLocalPar; i++) {
143     fVecBLoc[i]=0.;
144     
145     for (int j=0; j<fNLocalPar;j++) {
146       fMatCLoc[i][j]=0.;
147     }
148   }
149
150   // Then we fix all parameters...
151
152   for (int j=0; j<fNGlobalPar; j++) {
153     AliMillepede::SetParSigma(j,0.0);
154   }
155
156   fDerivLocEq.Reset();  fNDerivLocEq=0;  
157   fIndexLocEq.Reset();  fNIndexLocEq=0; 
158
159   fIndexAllEqs.Reset();  fNIndexAllEqs=0;
160   fDerivAllEqs.Reset();  fNDerivAllEqs=0;
161   fLocEqPlace.Reset();  fNLocEqPlace=0;
162
163   AliDebug(1,"");
164   AliDebug(1,"----------------------------------------------------");
165   AliDebug(1,"");
166   AliDebug(1,"    InitMille has been successfully called!");
167   AliDebug(1,"");
168   AliDebug(1,"-----------------------------------------------------");
169   AliDebug(1,"");
170         
171   return 1;
172 }
173
174 /*
175 -----------------------------------------------------------
176   PARGLO: initialization of global parameters
177 -----------------------------------------------------------
178
179   param    = array of starting values
180
181 -----------------------------------------------------------
182 */
183 Int_t AliMillepede::SetGlobalParameters(double *param)
184 {
185   /// initialization of global parameters
186   for(Int_t iPar=0; iPar<fNGlobalPar; iPar++){
187     fInitPar[iPar] = param[iPar];
188   }
189
190  return 1;
191 }
192
193 /*
194 -----------------------------------------------------------
195   PARGLO: initialization of global parameters
196 -----------------------------------------------------------
197
198   iPar    = the index of the global parameter in the 
199              result array (equivalent to fDeltaPar[]).
200
201   param    = the starting value
202
203 -----------------------------------------------------------
204 */
205 Int_t AliMillepede::SetGlobalParameter(int iPar, double param)
206 {
207   /// initialization of global parameter iPar
208   if (iPar<0 || iPar>=fNGlobalPar) {
209     return 0;
210   }
211   else {
212     fInitPar[iPar] = param;
213   }
214   return 1;
215 }
216
217
218 /*
219 -----------------------------------------------------------
220   PARSIG: define a constraint for a single global param
221           param is 'encouraged' to vary within [-sigma;sigma] 
222           range
223 -----------------------------------------------------------
224
225   iPar    = the index of the global parameter in the 
226              result array (equivalent to fDeltaPar[]).
227
228   sigma    = value of the constraint (sigma <= 0. will 
229              mean that parameter is FIXED !!!) 
230  
231 -----------------------------------------------------------
232 */ 
233 Int_t AliMillepede::SetParSigma(int iPar, double sigma)
234 {
235   /// Define a range [-sigma;sigma] where iPar is encourage to vary 
236   if (iPar>=fNGlobalPar) {
237     return 0;
238   }
239   else {
240     fSigmaPar[iPar] = sigma;
241   }
242
243   return 1;
244 }
245
246
247 /*
248 -----------------------------------------------------------
249   NONLIN: set nonlinear flag for a single global param
250           update of param durin iterations will not
251           consider initial starting value
252 -----------------------------------------------------------
253
254   iPar    = the index of the global parameter in the 
255              result array (equivalent to fDeltaPar[]).
256
257 -----------------------------------------------------------
258 */
259 Int_t AliMillepede::SetNonLinear(int iPar)
260 {
261   /// Set nonlinear flag for iPar 
262   if (iPar<0 || iPar>=fNGlobalPar) {
263     return 0;
264   }
265   else {
266     fIsNonLinear[iPar] = 1;
267   }
268   
269   return 1;
270 }
271
272
273 /*
274 -----------------------------------------------------------
275   INITUN: unit for iteration
276 -----------------------------------------------------------
277   
278   lChi2CutFac is used by Fitloc to define the Chi^2/ndof cut value
279
280   A large cutfac value enables to take a wider range of tracks 
281   for first iterations, which might be useful if misalignments
282   are large.
283
284   As soon as cutfac differs from 0 iteration are requested.
285   cutfac is then reduced, from one iteration to the other,
286   and iterations are stopped when it reaches the value 1.
287
288   At least one more iteration is often needed in order to remove
289   tracks containing outliers.
290   
291 -----------------------------------------------------------
292 */
293 Int_t AliMillepede::SetIterations(double lChi2CutFac)
294 {
295   /// Number of iterations is calculated from lChi2CutFac 
296   fChi2CutFactor = TMath::Max(1.0, lChi2CutFac);
297
298   AliInfo(Form("Initial cut factor is %f",fChi2CutFactor));
299   fIter = 1; // Initializes the iteration process
300   return 1;
301 }
302
303 /*
304 -----------------------------------------------------------
305   CONSTF: define a constraint equation in AliMillepede
306 -----------------------------------------------------------
307
308   dercs    = the row containing constraint equation 
309              derivatives (put into the final matrix)
310
311   lLagMult      = the lagrange multiplier value (sum of equation)            
312
313 -----------------------------------------------------------
314 */
315 Int_t AliMillepede::SetGlobalConstraint(double dercs[], double lLagMult)
316
317   /// Define a constraint equation
318   if (fNGlobalConstraints>=fgkMaxGloCsts) {  
319     AliInfo("Too many constraints !!!");
320     return 0;
321   }
322         
323   for (int i=0; i<fNGlobalPar; i++) {
324     fMatDerConstr[fNGlobalConstraints][i] = dercs[i];
325   }
326         
327   fLagMult[fNGlobalConstraints] = lLagMult;
328   fNGlobalConstraints++ ;
329   AliInfo(Form("Number of constraints increased to %d",fNGlobalConstraints));
330   return 1;
331 }
332
333
334 /*
335 -----------------------------------------------------------
336   EQULOC: write ONE equation in the matrices
337 -----------------------------------------------------------
338
339   dergb[1..fNGlobalPar] = global parameters derivatives
340   derlc[1..fNLocalPar]  = local parameters derivatives
341   rmeas                 = measured value
342   sigma                 = error on measured value (nothing to do with SetParSigma!!!)
343
344 -----------------------------------------------------------
345 */
346 Int_t AliMillepede::SetLocalEquation(double dergb[], double derlc[], double lMeas, double lSigma)
347 {       
348   /// Write one local equation
349   if (lSigma<=0.0) { // If parameter is fixed, then no equation
350     for (int i=0; i<fNLocalPar; i++) {
351       derlc[i] = 0.0;
352     }
353     for (int i=0; i<fNGlobalPar; i++) {
354       dergb[i] = 0.0;
355     }
356     return 1;
357   }
358     
359   // Serious equation, initialize parameters    
360   double lWeight =  1.0/(lSigma*lSigma);
361   int iLocFirst  = -1;
362   int iLocLast   = -1;
363   int iGlobFirst = -1;
364   int iGlobLast  = -1;
365  
366   for (int i=0; i<fNLocalPar; i++) { // Retrieve local param interesting indices
367     if (derlc[i]!=0.0) {
368       if (iLocFirst == -1) { 
369         iLocFirst = i;  // first index
370       }
371       iLocLast = i;     // last index
372     }
373   }
374   AliDebug(2,Form("%d / %d",iLocFirst, iLocLast));
375         
376   for (int i=0; i<fNGlobalPar; i++) { // Idem for global parameters
377     if (dergb[i]!=0.0) {
378       if (iGlobFirst == -1) {
379         iGlobFirst = i;  // first index
380       }
381       iGlobLast = i;     // last index
382     }
383   }
384   AliDebug(2,Form("%d / %d",iGlobFirst,iGlobLast));
385
386   if (fNIndexLocEq==fIndexLocEq.GetSize()) fIndexLocEq.Set(2*fNIndexLocEq);
387   fIndexLocEq.AddAt(-1,fNIndexLocEq++);    
388   if (fNDerivLocEq==fDerivLocEq.GetSize()) fDerivLocEq.Set(2*fNDerivLocEq);
389   fDerivLocEq.AddAt(lMeas,fNDerivLocEq++);    
390
391
392   for (int i=iLocFirst; i<=iLocLast; i++) { // Store interesting local parameters
393     if (derlc[i]!=0.0) {
394       if (fNIndexLocEq==fIndexLocEq.GetSize()) fIndexLocEq.Set(2*fNIndexLocEq);
395       fIndexLocEq.AddAt(i,fNIndexLocEq++);    
396       if (fNDerivLocEq==fDerivLocEq.GetSize()) fDerivLocEq.Set(2*fNDerivLocEq);
397       fDerivLocEq.AddAt(derlc[i],fNDerivLocEq++);    
398       derlc[i] = 0.0;
399     }
400   }
401
402   if (fNIndexLocEq==fIndexLocEq.GetSize()) fIndexLocEq.Set(2*fNIndexLocEq);
403   fIndexLocEq.AddAt(-1,fNIndexLocEq++);    
404   if (fNDerivLocEq==fDerivLocEq.GetSize()) fDerivLocEq.Set(2*fNDerivLocEq);
405   fDerivLocEq.AddAt(lWeight,fNDerivLocEq++);    
406
407   for (int i=iGlobFirst; i<=iGlobLast; i++) { // Store interesting global parameters
408     if (dergb[i]!=0.0) {
409       if (fNIndexLocEq==fIndexLocEq.GetSize()) fIndexLocEq.Set(2*fNIndexLocEq);
410       fIndexLocEq.AddAt(i,fNIndexLocEq++);    
411       if (fNDerivLocEq==fDerivLocEq.GetSize()) fDerivLocEq.Set(2*fNDerivLocEq);
412       fDerivLocEq.AddAt(dergb[i],fNDerivLocEq++);    
413       dergb[i] = 0.0;
414     }
415   }
416   
417   AliDebug(2,Form("Out Equloc --  NST = %d",fNDerivLocEq));
418   return 1;     
419 }
420
421 /*
422 -----------------------------------------------------------
423   FITLOC:  perform local params fit, once all the equations
424            have been written by EquLoc
425 -----------------------------------------------------------
426
427   iFit        = number of the fit, it is used to store 
428                 fit parameters and then retrieve them 
429                 for iterations (via FINDEXALLEQS and FDERIVALLEQS)
430
431   localParams = contains the fitted track parameters and
432                 related errors
433
434   bSingleFit  = is an option, if it is set to 1, we don't 
435                 perform the last loop. It is used to update 
436                 the track parameters without modifying global
437                 matrices
438
439 -----------------------------------------------------------
440 */
441 Int_t AliMillepede::LocalFit(int iFit, double localParams[], Bool_t bSingleFit)
442 {
443   /// Perform local parameters fit once all the local equations have been set
444
445   // Few initializations
446   int iEqTerm = 0;
447   int i, j, k;
448   int iIdx, jIdx, kIdx;
449   int iIdxIdx, jIdxIdx;
450   int iMeas   = -1;
451   int iWeight = 0;
452   int iLocFirst = 0;
453   int iLocLast  = 0;
454   int iGloFirst = 0;
455   int iGloLast  = 0;
456   int nGloInFit = 0;
457         
458   double lMeas    = 0.0;
459   double lWeight  = 0.0;
460
461   double lChi2    = 0.0;
462   double lRedChi2 = 0.0;
463   double lChi2Cut = 0.0;
464   int nEq  = 0;
465   int nDoF = 0;
466   int nEqTerms = fNDerivLocEq; // Number of terms (local + global derivatives,           
467                                // measurement and weight) involved in this local equation
468
469
470   // Fill the track store at first pass
471
472   if (fIter < 2 && !bSingleFit) { // Do it only once 
473     AliDebug(1,Form("Store equation no: %d", iFit)); 
474
475     for (i=0; i<nEqTerms; i++) { // Store the track parameters
476       if (fNIndexAllEqs==fIndexAllEqs.GetSize()) fIndexAllEqs.Set(2*fNIndexAllEqs);
477       fIndexAllEqs.AddAt(fIndexLocEq[i],fNIndexAllEqs++);
478       if (fNDerivAllEqs==fDerivAllEqs.GetSize()) fDerivAllEqs.Set(2*fNDerivAllEqs);
479       fDerivAllEqs.AddAt(fDerivLocEq[i],fNDerivAllEqs++);
480     }
481     if (fNLocEqPlace==fLocEqPlace.GetSize()) fLocEqPlace.Set(2*fNLocEqPlace);
482     fLocEqPlace.AddAt(fNIndexAllEqs,fNLocEqPlace++);
483
484     
485     AliDebug(2,Form("FLocEqPlace size = %d",fLocEqPlace[iFit])); 
486     AliDebug(2,Form("FIndexAllEqs size   = %d",fNIndexAllEqs)); 
487   }
488
489
490   for (i=0; i<fNLocalPar; i++) { // reset local params
491     fVecBLoc[i] = 0.0;
492
493     for (j=0; j<fNLocalPar; j++) {
494       fMatCLoc[i][j] = 0.0;
495     }
496   }
497         
498   for (i=0; i<fNGlobalPar; i++) {
499     fGlo2CGLRow[i] = -1;  // reset mixed params
500   }
501
502
503 /*
504
505   LOOPS : HOW DOES IT WORKS ?   
506
507   Now start by reading the informations stored with EquLoc.
508   Those informations are in vector FINDEXSTORE and FDERIVSTORE.
509   Each -1 in FINDEXSTORE delimits the equation parameters:
510   
511   First -1  ---> lMeas in FDERIVSTORE 
512   Then we have indices of local eq in FINDEXSTORE, and derivatives in FDERIVSTORE
513   Second -1 ---> weight in FDERIVSTORE
514   Then follows indices and derivatives of global eq.
515   ....
516   
517   We took them and store them into matrices.
518   
519   As we want ONLY local params, we substract the part of the estimated value
520   due to global params. Indeed we could have already an idea of these params,
521   with previous alignment constants for example (set with PARGLO). Also if there
522   are more than one iteration (FITLOC could be called by FITGLO)
523
524 */
525
526     
527 //
528 // FIRST LOOP : local track fit
529 //
530         
531   iEqTerm = 0;
532 //   fIndexLocEq.push_back(-1);
533   if (fNIndexLocEq==fIndexLocEq.GetSize()) fIndexLocEq.Set(2*fNIndexLocEq);
534   fIndexLocEq.AddAt(-1,fNIndexLocEq++);    
535   
536   while (iEqTerm <= nEqTerms) {
537     if (fIndexLocEq[iEqTerm] == -1) {
538       if (iMeas == -1) {        // First  -1 : lMeas
539         iMeas = iEqTerm;
540         iLocFirst = iEqTerm+1;
541       }  
542       else if (iWeight == 0) {  // Second -1 : weight 
543         iWeight = iEqTerm;
544         iLocLast = iEqTerm-1;
545         iGloFirst = iEqTerm+1;
546       } 
547       else {                    // Third  -1 : end of equation; start of next  
548         iGloLast = iEqTerm-1;
549         lMeas   = fDerivLocEq[iMeas];
550         lWeight         = fDerivLocEq[iWeight];
551 //      AliDebug(1,Form("lMeas = %f", lMeas));
552 //      AliDebug(1,Form("lWeight = %f", lWeight));
553         
554         // Now suppress the global part (only relevant with iterations)
555         // 
556         for (i=iGloFirst; i<=iGloLast; i++) {
557           iIdx = fIndexLocEq[i];              // Global param indice
558 //        AliDebug(2,Form("fDeltaPar[%d] = %f", iIdx, fDeltaPar[iIdx]));        
559 //        AliDebug(2,Form("Starting misalignment = %f",fInitPar[iIdx]));        
560           if (fIsNonLinear[iIdx] == 0)
561             lMeas -= fDerivLocEq[i]*(fInitPar[iIdx]+fDeltaPar[iIdx]); // linear parameter
562           else
563             lMeas -= fDerivLocEq[i]*(fDeltaPar[iIdx]); // nonlinear parameter
564         }
565 //      AliDebug(2,Form("lMeas after global stuff removal = %f", lMeas));
566                                 
567         for (i=iLocFirst; i<=iLocLast; i++) { // Finally fill local matrix and vector
568           iIdx = fIndexLocEq[i];   // Local param indice (the matrix line) 
569           fVecBLoc[iIdx] += lWeight*lMeas*fDerivLocEq[i];  
570 //        AliDebug(2,Form("fVecBLoc[%d] = %f", iIdx, fVecBLoc[iIdx]));
571                                         
572           for (j=iLocFirst; j<=i ; j++) { // Symmetric matrix, don't bother j>i coeffs
573             jIdx = fIndexLocEq[j];                                              
574             fMatCLoc[iIdx][jIdx] += lWeight*fDerivLocEq[i]*fDerivLocEq[j];          
575 //          AliDebug(2,Form("fMatCLoc[%d][%d] = ", iIdx, jIdx, fMatCLoc[iIdx][jIdx]));
576           }
577         }
578         iMeas   = -1;
579         iWeight = 0;
580         iEqTerm--;   // end of one equation is the beginning of next
581       } // End of "end of equation" operations
582     } // End of loop on equation
583     iEqTerm++;
584   } // End of loop on all equations used in the fit
585
586
587 //
588 // Local params matrix is completed, now invert to solve...
589 //
590         
591   Int_t nRank = AliMillepede::SpmInv(fMatCLoc, fVecBLoc, fNLocalPar);
592   // nRank is the number of nonzero diagonal elements 
593         
594   AliDebug(1,"");
595   AliDebug(1," __________________________________________________");
596   AliDebug(1,Form(" Printout of local fit  (FITLOC)  with rank= %d", nRank));
597   AliDebug(1," Result of local fit :      (index/parameter/error)");
598   
599   for (i=0; i<fNLocalPar; i++) {
600     AliDebug(1,Form("%d   /   %.6f   /   %.6f", i, fVecBLoc[i], TMath::Sqrt(fMatCLoc[i][i])));  
601   }
602   
603
604 // Store the track params and errors
605
606   for (i=0; i<fNLocalPar; i++) {
607     localParams[2*i] = fVecBLoc[i];
608     localParams[2*i+1] = TMath::Sqrt(TMath::Abs(fMatCLoc[i][i]));
609   }
610
611     
612 //
613 // SECOND LOOP : residual calculation
614 //
615   
616   iEqTerm = 0;
617   iMeas = -1;
618   iWeight = 0;
619
620   while (iEqTerm <= nEqTerms) {
621     if (fIndexLocEq[iEqTerm] == -1) {
622       if (iMeas == -1) {        // First  -1 : lMeas
623         iMeas = iEqTerm;
624         iLocFirst = iEqTerm+1;
625       }  
626       else if (iWeight == 0) {  // Second -1 : weight 
627         iWeight = iEqTerm;
628         iLocLast = iEqTerm-1;
629         iGloFirst = iEqTerm+1;
630       } 
631       else {                    // Third  -1 : end of equation; start of next  
632         iGloLast = iEqTerm-1;
633         lMeas   = fDerivLocEq[iMeas];
634         lWeight         = fDerivLocEq[iWeight];
635         
636         // Print all (for debugging purposes)
637
638 //      int nDerLoc = iLocLast-iLocFirst+1;   // Number of local derivatives involved
639 //      int nDerGlo = iGloLast-iGloFirst+1;   // Number of global derivatives involved
640
641 //      AliDebug(2,"");
642 //      AliDebug(2,Form(". equation:  measured value %.6f +/- %.6f", lMeas, 1.0/TMath::Sqrt(lWeight)));
643 //      AliDebug(2,Form("Number of derivatives (global, local): %d, %d",nDerGlo,nDerLoc));
644 //      AliDebug(2,"Global derivatives are: (index/derivative/parvalue) ");
645         
646 //      for (i=iGloFirst; i<=iGloLast; i++) {
647 //        AliDebug(2,Form("%d / %.6f / %.6f",fIndexLocEq[i],fDerivLocEq[i],fInitPar[fIndexLocEq[i]]));
648 //      } 
649
650 //      AliDebug(2,"Local derivatives are: (index/derivative) ");
651         
652 //      for (i=(ja+1); i<jb; i++) {AliDebug(2,Form("%d / %.6f",fIndexLocEq[i], fDerivLocEq[i]));}         
653
654         // Now suppress local and global parts to LMEAS;
655         //
656         // First the local part 
657         for (i=iLocFirst; i<=iLocLast; i++) { 
658           iIdx = fIndexLocEq[i];
659           lMeas -= fDerivLocEq[i]*fVecBLoc[iIdx];
660         }
661         // Then the global part
662         for (i=iGloFirst; i<=iGloLast; i++) {
663           iIdx = fIndexLocEq[i];
664           if (fIsNonLinear[iIdx] == 0)
665             lMeas -= fDerivLocEq[i]*(fInitPar[iIdx]+fDeltaPar[iIdx]); // linear parameter
666           else
667             lMeas -= fDerivLocEq[i]*(fDeltaPar[iIdx]); // nonlinear parameter
668         }
669
670         // lMeas contains now the residual value
671         AliDebug(2,Form("Residual value : %.6f", lMeas));
672
673         // reject the track if lMeas is too important (outlier)
674         if (TMath::Abs(lMeas) >= fResCutInit && fIter <= 1) {
675           AliDebug(2,"Rejected track !!!!!");
676           fNLocalFitsRejected++;      
677           fIndexLocEq.Reset();  fNIndexLocEq=0; // reset stores and go to the next track 
678           fDerivLocEq.Reset();  fNDerivLocEq=0;   
679           return 0;
680         }
681
682         if (TMath::Abs(lMeas) >= fResCut && fIter > 1) {
683           AliDebug(2,"Rejected track !!!!!");
684           fNLocalFitsRejected++;      
685           fIndexLocEq.Reset();  fNIndexLocEq=0; // reset stores and go to the next track 
686           fDerivLocEq.Reset();  fNDerivLocEq=0;   
687           return 0;
688         }
689
690         lChi2 += lWeight*lMeas*lMeas ; // total chi^2
691         nEq++;                    // number of equations                        
692         iMeas   = -1;
693         iWeight = 0;
694         iEqTerm--;
695       } // End of "end of equation" operations
696     }   // End of loop on equation
697     iEqTerm++;
698   } // End of loop on all equations used in the fit
699
700   nDoF = nEq-nRank;     
701   lRedChi2 = 0.0;
702
703   AliDebug(1,Form("Final chi square / degrees of freedom %.2f / %d",lChi2, nDoF));
704   
705   if (nDoF > 0) lRedChi2 = lChi2/float(nDoF);  // Chi^2/dof
706         
707   fNLocalFits++;
708
709   if (fNStdDev != 0 && nDoF > 0 && !bSingleFit) // Chisquare cut
710   {
711     lChi2Cut = AliMillepede::Chi2DoFLim(fNStdDev, nDoF)*fChi2CutFactor;
712     
713     AliDebug(1,Form("Reject if Chisq/Ndf = %.4f > %.4f",lRedChi2,lChi2Cut));
714  
715     if (lRedChi2 > lChi2Cut) // Reject the track if too much...
716     {
717       AliDebug(2,"Rejected track !!!!!");
718       fNLocalFitsRejected++;      
719       fIndexLocEq.Reset();  fNIndexLocEq=0; // reset stores and go to the next track 
720       fDerivLocEq.Reset();  fNDerivLocEq=0;
721       return 0;
722     }
723   }
724
725   if (bSingleFit) // Stop here if just updating the track parameters
726   {
727     fIndexLocEq.Reset();  fNIndexLocEq=0; // Reset store for the next track 
728     fDerivLocEq.Reset();  fNDerivLocEq=0;
729     return 1;
730   }
731
732 //  
733 // THIRD LOOP: local operations are finished, track is accepted 
734 // We now update the global parameters (other matrices)
735 //
736
737   iEqTerm = 0;
738   iMeas = -1;
739   iWeight = 0;
740
741   while (iEqTerm <= nEqTerms)
742   {
743     if (fIndexLocEq[iEqTerm] == -1)
744     {
745       if (iMeas == -1) {        // First  -1 : lMeas
746         iMeas = iEqTerm;
747         iLocFirst = iEqTerm+1;
748       }  
749       else if (iWeight == 0) {  // Second -1 : weight 
750         iWeight = iEqTerm;
751         iLocLast = iEqTerm-1;
752         iGloFirst = iEqTerm+1;
753       } 
754       else {                    // Third  -1 : end of equation; start of next  
755         iGloLast = iEqTerm-1;
756         lMeas   = fDerivLocEq[iMeas];
757         lWeight         = fDerivLocEq[iWeight];
758
759         // Now suppress the global part
760         for (i=iGloFirst; i<=iGloLast; i++) {
761           iIdx = fIndexLocEq[i];   // Global param indice
762           if (fIsNonLinear[iIdx] == 0)
763             lMeas -= fDerivLocEq[i]*(fInitPar[iIdx]+fDeltaPar[iIdx]); // linear parameter
764           else
765             lMeas -= fDerivLocEq[i]*(fDeltaPar[iIdx]); // nonlinear parameter
766         }
767         
768         for (i=iGloFirst; i<=iGloLast; i++) {
769           iIdx = fIndexLocEq[i];   // Global param indice (the matrix line)          
770         
771           fVecBGlo[iIdx] += lWeight*lMeas*fDerivLocEq[i];  
772 //        AliDebug(2,Form("fVecBGlo[%d] = %.6f", j, fVecBGlo[j] ));
773
774           // First of all, the global/global terms (exactly like local matrix)
775           //      
776           for (j=iGloFirst; j<=iGloLast; j++) {   
777             jIdx = fIndexLocEq[j];                      
778             fMatCGlo[iIdx][jIdx] += lWeight*fDerivLocEq[i]*fDerivLocEq[j];
779 //          AliDebug(2,Form("fMatCGlo[%d][%d] = %.6f",iIdx,jIdx,fMatCGlo[iIdx][jIdx]));
780           } 
781
782           // Now we have also rectangular matrices containing global/local terms.
783           //
784           iIdxIdx = fGlo2CGLRow[iIdx];  // Index of index          
785           if (iIdxIdx == -1) {    // New global variable         
786             for (k=0; k<fNLocalPar; k++) {
787               fMatCGloLoc[nGloInFit][k] = 0.0;  // Initialize the row
788             }
789             fGlo2CGLRow[iIdx] = nGloInFit;
790             fCGLRow2Glo[nGloInFit] = iIdx;
791             iIdxIdx = nGloInFit;
792             nGloInFit++;
793           }
794
795           // Now fill the rectangular matrix
796           for (k=iLocFirst; k<=iLocLast ; k++) {
797             kIdx = fIndexLocEq[k];                                              
798             fMatCGloLoc[iIdxIdx][kIdx] += lWeight*fDerivLocEq[i]*fDerivLocEq[k];
799 //          AliDebug(2,Form("fMatCGloLoc[%d][%d] = %.6f",iIdxIdx,kIdx,fMatCGloLoc[iIdxIdx][kIdx]));
800           } 
801         }
802         iMeas   = -1;
803         iWeight =  0;
804         iEqTerm--;
805       } // End of "end of equation" operations
806     }   // End of loop on equation
807     iEqTerm++;
808   } // End of loop on all equations used in the fit
809         
810   // Third loop is finished, now we update the correction matrices
811   AliMillepede::SpAVAt(fMatCLoc, fMatCGloLoc, fMatCGloCorr, fNLocalPar, nGloInFit);
812   AliMillepede::SpAX(fMatCGloLoc, fVecBLoc, fVecBGloCorr, fNLocalPar, nGloInFit);
813
814   for (iIdxIdx=0; iIdxIdx<nGloInFit; iIdxIdx++) {
815     iIdx = fCGLRow2Glo[iIdxIdx];
816     fVecBGlo[iIdx] -= fVecBGloCorr[iIdxIdx];
817     
818     for (jIdxIdx=0; jIdxIdx<=iIdxIdx; jIdxIdx++) {    
819       jIdx = fCGLRow2Glo[jIdxIdx];
820       fMatCGlo[iIdx][jIdx] -= fMatCGloCorr[iIdxIdx][jIdxIdx];
821       fMatCGlo[jIdx][iIdx] = fMatCGlo[iIdx][jIdx];
822     }
823   }
824         
825   fIndexLocEq.Reset();  fNIndexLocEq=0; // Reset store for the next track 
826   fDerivLocEq.Reset();  fNDerivLocEq=0;
827
828   return 1;
829 }
830  
831
832 /*
833 -----------------------------------------------------------
834   MAKEGLOBALFIT:  perform global params fit, once all the 'tracks'
835                   have been fitted by FitLoc
836 -----------------------------------------------------------
837
838   par[]        = array containing the computed global 
839                  parameters (the misalignment constants)
840
841   error[]      = array containing the error on global 
842                  parameters (estimated by AliMillepede)
843
844   pull[]        = array containing the corresponding pulls 
845
846 -----------------------------------------------------------
847 */
848 Int_t AliMillepede::GlobalFit(double par[], double error[], double pull[])
849 {
850   /// perform global parameters fit once all the local equations have been fitted
851   int i, j;
852   int nVar    = 0;
853   int nGloFix = 0;
854   double lConstraint;
855
856   double step[fgkMaxGlobalPar]={0};
857
858   double localPars[2*fgkMaxLocalPar];
859
860   int nLocFitsGood = 0;
861   int nLocFitsTot  = 0;
862   int nLocFits     = 0;
863
864   AliInfo("..... Making global fit .....");
865   AliInfo(Form("First Iteration with cut factor %.2f", fChi2CutFactor));
866
867   nLocFitsTot = AliMillepede::GetNLocalEquations();
868         
869   while (fIter < fMaxIter)  // Iteration for the final loop
870   {
871
872     nLocFits = AliMillepede::GetNLocalEquations();
873     AliInfo(Form("...using %d local fits...",nLocFits));
874
875 // Start by saving the diagonal elements
876     
877     for (i=0; i<fNGlobalPar; i++) {
878       fDiagCGlo[i] = fMatCGlo[i][i];
879     }
880
881 //  Then we retrieve the different constraints: fixed parameter or global equation
882
883     nGloFix = 0; // First look at the fixed global params
884     
885     for (i=0; i<fNGlobalPar; i++) {    
886       if (fSigmaPar[i] <= 0.0) {  // fixed global param
887         nGloFix++;
888         for (j=0; j<fNGlobalPar; j++) {
889           fMatCGlo[i][j] = 0.0;  // Reset row and column
890           fMatCGlo[j][i] = 0.0;
891         }
892       }
893       else {
894         fMatCGlo[i][i] += 1.0/(fSigmaPar[i]*fSigmaPar[i]);
895       }
896     }
897         
898     nVar = fNGlobalPar;  // Current number of equations 
899     AliDebug(1,Form("Number of constraint equations : %d", fNGlobalConstraints));
900     
901     for (i=0; i<fNGlobalConstraints; i++) { // Then the constraint equation    
902       lConstraint = fLagMult[i];
903       for (j=0; j<fNGlobalPar; j++) {   
904         fMatCGlo[nVar][j] = float(nLocFits)*fMatDerConstr[i][j];
905         fMatCGlo[j][nVar] = float(nLocFits)*fMatDerConstr[i][j];          
906         lConstraint -= fMatDerConstr[i][j]*(fInitPar[j]+fDeltaPar[j]);
907       }
908         
909       fMatCGlo[nVar][nVar] = 0.0;
910       fVecBGlo[nVar] = float(nLocFits)*lConstraint;
911       nVar++;
912     }
913
914
915     // Intended to compute the final global chisquare
916
917     double lFinalCor = 0.0;
918
919     if (fIter > 1) {    
920       for (i=0; i<fNGlobalPar; i++) {   
921         for (j=0; j<fNGlobalPar; j++) {
922 //        printf("%d, %d, %.6f  %.6f  %.6f\n",i,j,step[i],fMatCGlo[i][j],step[j]);
923           lFinalCor += step[i]*fMatCGlo[i][j]*step[j]; 
924           if (i == j && fSigmaPar[i] != 0) {
925             lFinalCor -= step[i]*step[i]/(fSigmaPar[i]*fSigmaPar[i]);
926           }
927         }
928       }
929     }
930
931     AliInfo(Form(" Final coeff is %.6f",lFinalCor));            
932     AliInfo(Form(" Final NDOFs = %d", fNGlobalPar));
933
934     //  The final matrix inversion
935
936     Int_t nRank = AliMillepede::SpmInv(fMatCGlo, fVecBGlo, nVar);
937
938     for (i=0; i<fNGlobalPar; i++) {    
939       fDeltaPar[i] += fVecBGlo[i];    // Update global parameters values (for iterations)
940       AliDebug(1,Form("fDeltaPar[%d] = %.6f", i, fDeltaPar[i]));
941       AliDebug(1,Form("fMatCGlo[%d][%d] = %.6f", i, i, fMatCGlo[i][i]));
942       AliDebug(1,Form("err = %.6f", TMath::Sqrt(TMath::Abs(fMatCGlo[i][i]))));
943
944       step[i] = fVecBGlo[i];
945
946       if (fIter == 1) error[i] = fMatCGlo[i][i]; // Unfitted error
947     }
948     AliInfo(Form("The rank defect of the symmetric %d by %d matrix is %d (bad if non 0)",
949             nVar, nVar, nVar-nGloFix-nRank));
950                 
951     AliInfo(Form("Total : %d local fits, %d rejected.", fNLocalFits, fNLocalFitsRejected));
952     if (fIter == 0)  break;  // No iterations set     
953     fIter++;
954         
955     if (fIter == fMaxIter)  break;  // End of story         
956
957     // Reinitialize parameters for iteration
958     //
959     fNLocalFits = 0;
960     fNLocalFitsRejected = 0;
961
962     if (fChi2CutFactor != fChi2CutRef) {    
963       fChi2CutFactor = TMath::Sqrt(fChi2CutFactor);
964       if (fChi2CutFactor < 1.2*fChi2CutRef) {
965         fChi2CutFactor = fChi2CutRef;
966         fIter = fMaxIter - 1;     // Last iteration
967       }
968     }
969     AliInfo("");
970     AliInfo("----------------------------------------------");
971     AliInfo(Form("Iteration %d with cut factor %.2f", fIter, fChi2CutFactor));
972     
973     // Reset global variables
974     //    
975     for (i=0; i<nVar; i++) {    
976       fVecBGlo[i] = 0.0;
977       for (j=0; j<nVar; j++) {     
978         fMatCGlo[i][j] = 0.0;
979       }
980     }
981
982     //
983     // We start a new iteration
984     //
985
986     // First we read the stores for retrieving the local params
987     //
988     nLocFitsGood = 0;
989
990     for (i=0; i<nLocFitsTot; i++) {
991       int iEqFirst = 0;
992       int iEqLast = 0;
993
994       (i>0) ? iEqFirst = fLocEqPlace[i-1] : iEqFirst = 0;
995       iEqLast = fLocEqPlace[i];
996
997       AliDebug(2,Form("Track %d : ",i));
998       AliDebug(2,Form("Starts at %d", iEqFirst));
999       AliDebug(2,Form("Ends at %d",iEqLast));
1000
1001       if (fIndexAllEqs[iEqFirst] != -999) { // Fit is still OK      
1002         fIndexLocEq.Reset();  fNIndexLocEq=0;
1003         fDerivLocEq.Reset();  fNDerivLocEq=0;
1004         
1005         for (j=iEqFirst; j<iEqLast; j++) {
1006           if (fNIndexLocEq==fIndexLocEq.GetSize()) fIndexLocEq.Set(2*fNIndexLocEq);
1007           fIndexLocEq.AddAt(fIndexAllEqs[j],fNIndexLocEq++);    
1008           if (fNDerivLocEq==fDerivLocEq.GetSize()) fDerivLocEq.Set(2*fNDerivLocEq);
1009           fDerivLocEq.AddAt(fDerivAllEqs[j],fNDerivLocEq++);    
1010         }
1011
1012         for (j=0; j<2*fNLocalPar; j++) {
1013           localPars[j] = 0.;
1014         }       
1015
1016 //      Int_t sc = AliMillepede::LocalFit(i,localPars,0);
1017 //      (sc) ? nLocFitsGood++ : fIndexAllEqs[iEqFirst] = -999;   
1018         AliMillepede::LocalFit(i,localPars,0);
1019         nLocFitsGood++;
1020       }
1021     } // End of loop on fits
1022
1023     AliMillepede::SetNLocalEquations(nLocFitsGood);
1024
1025   } // End of iteration loop
1026         
1027   AliMillepede::PrintGlobalParameters(); // Print the final results
1028
1029   for (j=0; j<fNGlobalPar; j++) {  
1030     par[j]   = fInitPar[j]+fDeltaPar[j];
1031     pull[j]  = (fSigmaPar[j] <= 0. || fSigmaPar[j]*fSigmaPar[j]-fMatCGlo[j][j] <=0.) ? 0. : fDeltaPar[j]/TMath::Sqrt(fSigmaPar[j]*fSigmaPar[j]-fMatCGlo[j][j]);
1032     error[j] = TMath::Sqrt(TMath::Abs(fMatCGlo[j][j]));
1033   }
1034
1035   AliInfo(" ");
1036   AliInfo("            * o o                   o      ");
1037   AliInfo("              o o                   o      ");
1038   AliInfo("   o ooooo  o o o  oo  ooo   oo   ooo  oo  ");
1039   AliInfo("    o  o  o o o o o  o o  o o  o o  o o  o ");
1040   AliInfo("    o  o  o o o o oooo o  o oooo o  o oooo ");
1041   AliInfo("    o  o  o o o o o    ooo  o    o  o o    ");
1042   AliInfo("    o  o  o o o o  oo  o     oo   ooo  oo ++ ends.");
1043   AliInfo("                       o                   ");         
1044
1045   return 1;
1046 }
1047
1048 /*
1049 -----------------------------------------------------------
1050   ERRPAR: return error for parameter iPar
1051 -----------------------------------------------------------
1052
1053   iPar     = the index of the global parameter in the 
1054              result array (equivalent to fDeltaPar[]).
1055  
1056 -----------------------------------------------------------
1057 */ 
1058 Double_t AliMillepede::GetParError(int iPar) const
1059 {
1060   /// return error for parameter iPar
1061   Double_t lErr = -1.;
1062   if (iPar>=0 && iPar<fNGlobalPar) {
1063     lErr = TMath::Sqrt(TMath::Abs(fMatCGlo[iPar][iPar]));
1064   }
1065   return lErr;
1066 }
1067
1068
1069 /*
1070 -----------------------------------------------------------
1071   SPMINV:  obtain solution of a system of linear equations with symmetric matrix 
1072            and the inverse (using 'singular-value friendly' GAUSS pivot)
1073 -----------------------------------------------------------
1074
1075         Solve the equation :  V * X = B
1076         
1077         V is replaced by inverse matrix and B by X, the solution vector
1078 -----------------------------------------------------------
1079 */
1080 int AliMillepede::SpmInv(double matV[][fgkMaxGloPC], double vecB[], int nGlo)
1081 {
1082   ///  Obtain solution of a system of linear equations with symmetric matrix 
1083   ///  and the inverse (using 'singular-value friendly' GAUSS pivot)
1084     
1085   Int_t nRank = 0;
1086   int iPivot;
1087   double vPivot = 0.;
1088   double eps = 0.00000000000001;
1089
1090   bool *bUnUsed = new bool[nGlo];
1091   double *diagV = new double[nGlo];
1092   double *rowMax = new double[nGlo];
1093   double *colMax = new double[nGlo];
1094
1095   double *temp = new double[nGlo];
1096
1097   for (Int_t i=0; i<nGlo; i++) {
1098     rowMax[i] = 0.0;
1099     colMax[i] = 0.0;
1100     bUnUsed[i] = true;
1101
1102     for (Int_t j=0; j<i; j++) {
1103       if (matV[j][i] == 0) {
1104         matV[j][i] = matV[i][j];
1105       }
1106     }
1107   }
1108   
1109   // Small loop for matrix equilibration (gives a better conditioning) 
1110
1111   for (Int_t i=0; i<nGlo; i++) {
1112     for (Int_t j=0; j<nGlo; j++) { 
1113       if (TMath::Abs(matV[i][j]) >= rowMax[i]) rowMax[i] = TMath::Abs(matV[i][j]); // Max elemt of row i
1114       if (TMath::Abs(matV[j][i]) >= colMax[i]) colMax[i] = TMath::Abs(matV[j][i]); // Max elemt of column i
1115     }
1116   }
1117
1118   for (Int_t i=0; i<nGlo; i++) {
1119     if (0.0 != rowMax[i]) rowMax[i] = 1./rowMax[i]; // Max elemt of row i
1120     if (0.0 != colMax[i]) colMax[i] = 1./colMax[i]; // Max elemt of column i
1121   }
1122
1123   for (Int_t i=0; i<nGlo; i++) {
1124     for (Int_t j=0; j<nGlo; j++) {
1125       matV[i][j] = TMath::Sqrt(rowMax[i])*matV[i][j]*TMath::Sqrt(colMax[j]); // Equilibrate the V matrix
1126     }
1127     diagV[i] = TMath::Abs(matV[i][i]); // save diagonal elem absolute values    
1128   }
1129
1130
1131   for (Int_t i=0; i<nGlo; i++) {
1132     vPivot = 0.0;
1133     iPivot = -1;
1134     
1135     for (Int_t j=0; j<nGlo; j++) { // First look for the pivot, ie max unused diagonal element       
1136       if (bUnUsed[j] && (TMath::Abs(matV[j][j])>TMath::Max(TMath::Abs(vPivot),eps*diagV[j]))) {    
1137         vPivot = matV[j][j];
1138         iPivot = j;
1139       }
1140     }
1141     
1142     if (iPivot >= 0) {   // pivot found          
1143       nRank++;
1144       bUnUsed[iPivot] = false; // This value is used
1145       vPivot = 1.0/vPivot;
1146       matV[iPivot][iPivot] = -vPivot; // Replace pivot by its inverse
1147       
1148       for (Int_t j=0; j<nGlo; j++) {      
1149         for (Int_t jj=0; jj<nGlo; jj++) {  
1150           if (j != iPivot && jj != iPivot) {// Other elements (!!! do them first as you use old matV[k][j]'s !!!)         
1151             matV[j][jj] = matV[j][jj] - vPivot*matV[j][iPivot]*matV[iPivot][jj];
1152           }
1153         }
1154       }
1155       
1156       for (Int_t j=0; j<nGlo; j++) {      
1157         if (j != iPivot) { // Pivot row or column elements 
1158           matV[j][iPivot] = matV[j][iPivot]*vPivot;     // Column
1159           matV[iPivot][j] = matV[iPivot][j]*vPivot;     // Line
1160         }
1161       }
1162     }
1163     else {  // No more pivot value (clear those elements)
1164       for (Int_t j=0; j<nGlo; j++) {
1165         if (bUnUsed[j]) {
1166           vecB[j] = 0.0;
1167
1168           for (Int_t k=0; k<nGlo; k++) {
1169             matV[j][k] = 0.0;
1170             matV[k][j] = 0.0;
1171           }
1172         }
1173       }
1174       break;  // No more pivots anyway, stop here
1175     }
1176   }
1177   
1178   for (Int_t i=0; i<nGlo; i++) {
1179     for (Int_t j=0; j<nGlo; j++) {
1180       matV[i][j] = TMath::Sqrt(colMax[i])*matV[i][j]*TMath::Sqrt(rowMax[j]); // Correct matrix V
1181     }
1182   }
1183   
1184   for (Int_t j=0; j<nGlo; j++) {
1185     temp[j] = 0.0;
1186     
1187     for (Int_t jj=0; jj<nGlo; jj++) { // Reverse matrix elements
1188       matV[j][jj] = -matV[j][jj];
1189       temp[j] += matV[j][jj]*vecB[jj];
1190     }           
1191   }
1192
1193   for (Int_t j=0; j<nGlo; j++) {
1194     vecB[j] = temp[j]; // The final result
1195   }
1196   
1197   delete [] temp;
1198   delete [] bUnUsed;
1199   delete [] diagV;
1200   delete [] rowMax;
1201   delete [] colMax;
1202
1203   return nRank;
1204 }
1205
1206 //
1207 // Same method but for local fit, so heavily simplified
1208 //
1209 int AliMillepede::SpmInv(double matV[][fgkMaxLocalPar], double vecB[], int nLoc)
1210 {
1211   ///  Obtain solution of a system of linear equations with symmetric matrix 
1212   ///  and the inverse (using 'singular-value friendly' GAUSS pivot)
1213
1214   Int_t nRank = 0;
1215   Int_t iPivot = -1;
1216   double vPivot = 0.;
1217   double eps = 0.0000000000001;
1218
1219   bool *bUnUsed = new bool[nLoc];
1220   double *diagV  = new double[nLoc];
1221   double *temp  = new double[nLoc];
1222   
1223   for (Int_t i=0; i<nLoc; i++) {
1224     bUnUsed[i] = true;
1225     diagV[i] = TMath::Abs(matV[i][i]);     // save diagonal elem absolute values
1226     for (Int_t j=0; j<i; j++) {
1227       matV[j][i] = matV[i][j] ;
1228     }
1229   }
1230         
1231   for (Int_t i=0; i<nLoc; i++) {
1232     vPivot = 0.0;
1233     iPivot = -1;
1234                 
1235     for (Int_t j=0; j<nLoc; j++) { // First look for the pivot, ie max unused diagonal element 
1236       if (bUnUsed[j] && (TMath::Abs(matV[j][j])>TMath::Max(TMath::Abs(vPivot),eps*diagV[j]))) {
1237         vPivot = matV[j][j];
1238         iPivot = j;
1239       }
1240     }
1241                 
1242     if (iPivot >= 0) {   // pivot found
1243       nRank++;
1244       bUnUsed[iPivot] = false;
1245       vPivot = 1.0/vPivot;
1246       matV[iPivot][iPivot] = -vPivot; // Replace pivot by its inverse
1247       
1248       for (Int_t j=0; j<nLoc; j++) {
1249         if (j != iPivot) {
1250           for (Int_t jj=0; jj<=j; jj++) {       
1251             if (jj != iPivot) {// Other elements (!!! do them first as you use old matV[k][j]'s !!!)
1252               matV[j][jj] = matV[j][jj] - vPivot*matV[j][iPivot]*matV[iPivot][jj];
1253               matV[jj][j] = matV[j][jj];
1254             }
1255           }
1256         }
1257       }
1258
1259       for (Int_t j=0; j<nLoc; j++) {      
1260         if (j != iPivot) {      // Pivot row or column elements         
1261           matV[j][iPivot] = matV[j][iPivot]*vPivot; // Column
1262           matV[iPivot][j] = matV[j][iPivot];
1263         }
1264       }
1265     }
1266     else { // No more pivot value (clear those elements)
1267       for (Int_t j=0; j<nLoc; j++) {
1268         if (bUnUsed[j]) {
1269           vecB[j] = 0.0;
1270           matV[j][j] = 0.0;
1271           for (Int_t k=0; k<j; k++) {     
1272             matV[j][k] = 0.0;
1273             matV[k][j] = 0.0;
1274           }
1275         }
1276       }
1277
1278       break;  // No more pivots anyway, stop here
1279     }
1280   }
1281
1282   for (Int_t j=0; j<nLoc; j++) {  
1283     temp[j] = 0.0;    
1284     for (Int_t jj=0; jj<nLoc; jj++) { // Reverse matrix elements
1285       matV[j][jj] = -matV[j][jj];
1286       temp[j] += matV[j][jj]*vecB[jj];
1287     }                   
1288   }
1289
1290   for (Int_t j=0; j<nLoc; j++) {
1291     vecB[j] = temp[j];
1292   }
1293
1294   delete [] bUnUsed;
1295   delete [] diagV;
1296   delete [] temp;
1297   
1298   return nRank;
1299 }
1300
1301
1302 /*
1303 -----------------------------------------------------------
1304   SPAVAT
1305 -----------------------------------------------------------
1306
1307   multiply symmetric N-by-N matrix from the left with general M-by-N
1308   matrix and from the right with the transposed of the same  general
1309   matrix  to  form  symmetric  M-by-M   matrix.
1310   
1311                                                        T
1312   CALL SPAVAT(V,A,W,N,M)      W   =   A   *   V   *   A
1313                            M*M     M*N     N*N     N*M
1314   
1315   where V = symmetric N-by-N matrix
1316         A = general N-by-M matrix
1317         W = symmetric M-by-M matrix
1318 -----------------------------------------------------------
1319 */
1320 Int_t AliMillepede::SpAVAt(double matV[][fgkMaxLocalPar], double matA[][fgkMaxLocalPar], double matW[][fgkMaxGlobalPar], int nLoc, int nGlo)
1321 {
1322   ///  multiply symmetric N-by-N matrix from the left with general M-by-N
1323   ///  matrix and from the right with the transposed of the same general
1324   ///  matrix to form symmetric M-by-M matrix.
1325
1326   for (Int_t i=0; i<nGlo; i++) {
1327     for (Int_t j=0; j<=i; j++) {  // Matrix w is symmetric
1328       matW[i][j] = 0.0;      // Reset final matrix                      
1329       for (Int_t k=0; k<nLoc; k++) {    
1330         matW[i][j] += matA[i][k]*matV[k][k]*matA[j][k];    // diagonal terms of v
1331         for (Int_t l=0; l<k; l++) {       
1332           matW[i][j] += matA[i][k]*matV[k][l]*matA[j][l];  // Use symmetric properties
1333           matW[i][j] += matA[i][l]*matV[k][l]*matA[j][k];  // of matrix v
1334         }
1335       }
1336       if (i!=j){
1337         matW[j][i] = matW[i][j]; // Matrix w is symmetric
1338       }
1339     }
1340   }
1341         
1342   return 1;
1343 }
1344
1345
1346 /*
1347 -----------------------------------------------------------
1348   SPAX
1349 -----------------------------------------------------------
1350
1351   multiply general M-by-N matrix A and N-vector X
1352  
1353   CALL  SPAX(A,X,Y,M,N)          Y =  A * X
1354                                  M   M*N  N
1355  
1356   where A = general M-by-N matrix (A11 A12 ... A1N  A21 A22 ...)
1357         X = N vector
1358         Y = M vector
1359 -----------------------------------------------------------
1360 */
1361 Int_t AliMillepede::SpAX(double matA[][fgkMaxLocalPar], double vecX[], double vecY[], int nCol, int nRow)
1362 {
1363   ///   multiply general M-by-N matrix A and N-vector X
1364   for (Int_t i=0; i<nRow; i++) {
1365     vecY[i] = 0.0;          // Reset final vector                       
1366     for (Int_t j=0; j<nCol; j++) {
1367       vecY[i] += matA[i][j]*vecX[j];  // fill the vector
1368     }
1369   }
1370         
1371   return 1;
1372 }
1373
1374
1375 /*
1376 -----------------------------------------------------------
1377   PRTGLO
1378 -----------------------------------------------------------
1379
1380   Print the final results into the logfile
1381
1382 -----------------------------------------------------------
1383 */
1384 Int_t AliMillepede::PrintGlobalParameters() const
1385 {
1386   ///  Print the final results into the logfile
1387   double lError = 0.;
1388   double lGlobalCor =0.;
1389         
1390   AliInfo("");
1391   AliInfo("   Result of fit for global parameters");
1392   AliInfo("   ===================================");
1393   AliInfo("    I       initial       final       differ        lastcor        error       gcor");
1394   AliInfo("-----------------------------------------------------------------------------------");
1395         
1396   for (int i=0; i<fNGlobalPar; i++) {
1397     lError = TMath::Sqrt(TMath::Abs(fMatCGlo[i][i]));
1398     if (fMatCGlo[i][i] < 0.0) lError = -lError;
1399     lGlobalCor = 0.0;
1400                 
1401     if (TMath::Abs(fMatCGlo[i][i]*fDiagCGlo[i]) > 0) {    
1402       lGlobalCor = TMath::Sqrt(TMath::Abs(1.0-1.0/(fMatCGlo[i][i]*fDiagCGlo[i])));
1403       AliInfo(Form("%d\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f",
1404                    i,fInitPar[i],fInitPar[i]+fDeltaPar[i],fDeltaPar[i],fVecBGlo[i],lError,lGlobalCor));
1405     }
1406     else {    
1407       AliInfo(Form("%d\t %.6f\t %.6f\t %.6f\t %.6f\t OFF\t OFF",i,fInitPar[i],fInitPar[i]+fDeltaPar[i],fDeltaPar[i],fVecBGlo[i]));
1408     }
1409   }
1410   return 1;
1411 }
1412
1413
1414 /*
1415 ----------------------------------------------------------------
1416   CHI2DOFLIM:  return the limit in chi^2/nd for n sigmas stdev authorized
1417 ----------------------------------------------------------------
1418
1419   Only n=1, 2, and 3 are expected in input
1420 ----------------------------------------------------------------
1421 */
1422 double AliMillepede::Chi2DoFLim(int nSig, int nDoF)
1423 {
1424   /// return the limit in chi^2/nd for n sigmas stdev authorized
1425   int lNSig;
1426   double sn[3]        = {0.47523, 1.690140, 2.782170};
1427   double table[3][30] = {{1.0000, 1.1479, 1.1753, 1.1798, 1.1775, 1.1730, 1.1680, 1.1630,
1428                           1.1581, 1.1536, 1.1493, 1.1454, 1.1417, 1.1383, 1.1351, 1.1321,
1429                           1.1293, 1.1266, 1.1242, 1.1218, 1.1196, 1.1175, 1.1155, 1.1136,
1430                           1.1119, 1.1101, 1.1085, 1.1070, 1.1055, 1.1040},
1431                          {4.0000, 3.0900, 2.6750, 2.4290, 2.2628, 2.1415, 2.0481, 1.9736,
1432                           1.9124, 1.8610, 1.8171, 1.7791, 1.7457, 1.7161, 1.6897, 1.6658,
1433                           1.6442, 1.6246, 1.6065, 1.5899, 1.5745, 1.5603, 1.5470, 1.5346,
1434                           1.5230, 1.5120, 1.5017, 1.4920, 1.4829, 1.4742},
1435                          {9.0000, 5.9146, 4.7184, 4.0628, 3.6410, 3.3436, 3.1209, 2.9468,
1436                           2.8063, 2.6902, 2.5922, 2.5082, 2.4352, 2.3711, 2.3143, 2.2635,
1437                           2.2178, 2.1764, 2.1386, 2.1040, 2.0722, 2.0428, 2.0155, 1.9901,
1438                           1.9665, 1.9443, 1.9235, 1.9040, 1.8855, 1.8681}};
1439
1440   if (nDoF < 1) {
1441     return 0.0;
1442   }
1443   else {  
1444     lNSig = TMath::Max(1,TMath::Min(nSig,3));
1445
1446     if (nDoF <= 30) {    
1447       return table[lNSig-1][nDoF-1];
1448     }
1449     else { // approximation
1450       return ((sn[lNSig-1]+TMath::Sqrt(float(2*nDoF-3)))*
1451               (sn[lNSig-1]+TMath::Sqrt(float(2*nDoF-3))))/float(2*nDoF-2);
1452     }
1453   }
1454 }