]> git.uio.no Git - u/mrichter/AliRoot.git/blob - MUON/AliMillepede.cxx
DP:Misalignment of CPV added
[u/mrichter/AliRoot.git] / MUON / 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 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       int 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[1010];
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
866   nLocFitsTot = AliMillepede::GetNLocalEquations();
867         
868   while (fIter < fMaxIter)  // Iteration for the final loop
869   {
870
871     nLocFits = AliMillepede::GetNLocalEquations();
872     AliInfo(Form("...using %d local fits...",nLocFits));
873
874 // Start by saving the diagonal elements
875     
876     for (i=0; i<fNGlobalPar; i++) {
877       fDiagCGlo[i] = fMatCGlo[i][i];
878     }
879
880 //  Then we retrieve the different constraints: fixed parameter or global equation
881
882     nGloFix = 0; // First look at the fixed global params
883     
884     for (i=0; i<fNGlobalPar; i++) {    
885       if (fSigmaPar[i] <= 0.0) {  // fixed global param
886         nGloFix++;
887         for (j=0; j<fNGlobalPar; j++) {
888           fMatCGlo[i][j] = 0.0;  // Reset row and column
889           fMatCGlo[j][i] = 0.0;
890         }
891       }
892       else {
893         fMatCGlo[i][i] += 1.0/(fSigmaPar[i]*fSigmaPar[i]);
894       }
895     }
896         
897     nVar = fNGlobalPar;  // Current number of equations 
898     AliDebug(1,Form("Number of constraint equations : %d", fNGlobalConstraints));
899     
900     for (i=0; i<fNGlobalConstraints; i++) { // Then the constraint equation    
901       lConstraint = fLagMult[i];
902       for (j=0; j<fNGlobalPar; j++) {   
903         fMatCGlo[nVar][j] = float(nLocFits)*fMatDerConstr[i][j];
904         fMatCGlo[j][nVar] = float(nLocFits)*fMatDerConstr[i][j];          
905         lConstraint -= fMatDerConstr[i][j]*(fInitPar[j]+fDeltaPar[j]);
906       }
907         
908       fMatCGlo[nVar][nVar] = 0.0;
909       fVecBGlo[nVar] = float(nLocFits)*lConstraint;
910       nVar++;
911     }
912
913
914     // Intended to compute the final global chisquare
915
916     double lFinalCor = 0.0;
917
918     if (fIter > 1) {    
919       for (i=0; i<fNGlobalPar; i++) {   
920         for (j=0; j<fNGlobalPar; j++) {
921 //        printf("%d, %d, %.6f  %.6f  %.6f\n",i,j,step[i],fMatCGlo[i][j],step[j]);
922           lFinalCor += step[i]*fMatCGlo[i][j]*step[j]; 
923           if (i == j && fSigmaPar[i] != 0) {
924             lFinalCor -= step[i]*step[i]/(fSigmaPar[i]*fSigmaPar[i]);
925           }
926         }
927       }
928     }
929
930     AliInfo(Form(" Final coeff is %.6f",lFinalCor));            
931     AliInfo(Form(" Final NDOFs = %d", fNGlobalPar));
932
933     //  The final matrix inversion
934
935     Int_t nRank = AliMillepede::SpmInv(fMatCGlo, fVecBGlo, nVar);
936
937     for (i=0; i<fNGlobalPar; i++) {    
938       fDeltaPar[i] += fVecBGlo[i];    // Update global parameters values (for iterations)
939       AliDebug(1,Form("fDeltaPar[%d] = %.6f", i, fDeltaPar[i]));
940       AliDebug(1,Form("fMatCGlo[%d][%d] = %.6f", i, i, fMatCGlo[i][i]));
941       AliDebug(1,Form("err = %.6f", TMath::Sqrt(TMath::Abs(fMatCGlo[i][i]))));
942
943       step[i] = fVecBGlo[i];
944
945       if (fIter == 1) error[i] = fMatCGlo[i][i]; // Unfitted error
946     }
947     AliInfo("");
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("");
952     AliInfo(Form("Total : %d local fits, %d rejected.", fNLocalFits, fNLocalFitsRejected));
953     if (fIter == 0)  break;  // No iterations set     
954     fIter++;
955         
956     if (fIter == fMaxIter)  break;  // End of story         
957
958     // Reinitialize parameters for iteration
959     //
960     fNLocalFits = 0;
961     fNLocalFitsRejected = 0;
962
963     if (fChi2CutFactor != fChi2CutRef) {    
964       fChi2CutFactor = TMath::Sqrt(fChi2CutFactor);
965       if (fChi2CutFactor < 1.2*fChi2CutRef) {
966         fChi2CutFactor = fChi2CutRef;
967         fIter = fMaxIter - 1;     // Last iteration
968       }
969     }
970     AliInfo(Form("Iteration %d with cut factor %.2f", fIter, fChi2CutFactor));
971     
972     // Reset global variables
973     //    
974     for (i=0; i<nVar; i++) {    
975       fVecBGlo[i] = 0.0;
976       for (j=0; j<nVar; j++) {     
977         fMatCGlo[i][j] = 0.0;
978       }
979     }
980
981     //
982     // We start a new iteration
983     //
984
985     // First we read the stores for retrieving the local params
986     //
987     nLocFitsGood = 0;
988
989     for (i=0; i<nLocFitsTot; i++) {
990       int iEqFirst = 0;
991       int iEqLast = 0;
992
993       (i>0) ? iEqFirst = fLocEqPlace[i-1] : iEqFirst = 0;
994       iEqLast = fLocEqPlace[i];
995
996       AliDebug(2,Form("Track %d : ",i));
997       AliDebug(2,Form("Starts at %d", iEqFirst));
998       AliDebug(2,Form("Ends at %d",iEqLast));
999
1000       if (fIndexAllEqs[iEqFirst] != -999) { // Fit is still OK      
1001         fIndexLocEq.Reset();  fNIndexLocEq=0;
1002         fDerivLocEq.Reset();  fNDerivLocEq=0;
1003         
1004         for (j=iEqFirst; j<iEqLast; j++) {
1005           if (fNIndexLocEq==fIndexLocEq.GetSize()) fIndexLocEq.Set(2*fNIndexLocEq);
1006           fIndexLocEq.AddAt(fIndexAllEqs[j],fNIndexLocEq++);    
1007           if (fNDerivLocEq==fDerivLocEq.GetSize()) fDerivLocEq.Set(2*fNDerivLocEq);
1008           fDerivLocEq.AddAt(fDerivAllEqs[j],fNDerivLocEq++);    
1009         }
1010
1011         for (j=0; j<2*fNLocalPar; j++) {
1012           localPars[j] = 0.;
1013         }       
1014
1015 //      Int_t sc = AliMillepede::LocalFit(i,localPars,0);
1016 //      (sc) ? nLocFitsGood++ : fIndexAllEqs[iEqFirst] = -999;   
1017         AliMillepede::LocalFit(i,localPars,0);
1018         nLocFitsGood++;
1019       }
1020     } // End of loop on fits
1021
1022     AliMillepede::SetNLocalEquations(nLocFitsGood);
1023
1024   } // End of iteration loop
1025         
1026   AliMillepede::PrintGlobalParameters(); // Print the final results
1027
1028   for (j=0; j<fNGlobalPar; j++) {  
1029     par[j]   = fInitPar[j]+fDeltaPar[j];
1030     pull[j]  = (fSigmaPar[j] <= 0.) ? 0. : fDeltaPar[j]/TMath::Sqrt(fSigmaPar[j]*fSigmaPar[j]-fMatCGlo[j][j]);
1031     error[j] = TMath::Sqrt(TMath::Abs(fMatCGlo[j][j]));
1032   }
1033
1034   AliInfo(" ");
1035   AliInfo("            * o o                   o      ");
1036   AliInfo("              o o                   o      ");
1037   AliInfo("   o ooooo  o o o  oo  ooo   oo   ooo  oo  ");
1038   AliInfo("    o  o  o o o o o  o o  o o  o o  o o  o ");
1039   AliInfo("    o  o  o o o o oooo o  o oooo o  o oooo ");
1040   AliInfo("    o  o  o o o o o    ooo  o    o  o o    ");
1041   AliInfo("    o  o  o o o o  oo  o     oo   ooo  oo ++ ends.");
1042   AliInfo("                       o                   ");         
1043
1044   return 1;
1045 }
1046
1047 /*
1048 -----------------------------------------------------------
1049   ERRPAR: return error for parameter iPar
1050 -----------------------------------------------------------
1051
1052   iPar     = the index of the global parameter in the 
1053              result array (equivalent to fDeltaPar[]).
1054  
1055 -----------------------------------------------------------
1056 */ 
1057 Double_t AliMillepede::GetParError(Int_t iPar) const
1058 {
1059   /// return error for parameter iPar
1060   Double_t lErr = -1.;
1061   if (iPar>=0 && iPar<fNGlobalPar) {
1062     lErr = TMath::Sqrt(TMath::Abs(fMatCGlo[iPar][iPar]));
1063   }
1064   return lErr;
1065 }
1066
1067
1068 /*
1069 -----------------------------------------------------------
1070   SPMINV:  obtain solution of a system of linear equations with symmetric matrix 
1071            and the inverse (using 'singular-value friendly' GAUSS pivot)
1072 -----------------------------------------------------------
1073
1074         Solve the equation :  V * X = B
1075         
1076         V is replaced by inverse matrix and B by X, the solution vector
1077 -----------------------------------------------------------
1078 */
1079 int AliMillepede::SpmInv(double matV[][fgkMaxGloPC], double vecB[], int nGlo)
1080 {
1081   ///  Obtain solution of a system of linear equations with symmetric matrix 
1082   ///  and the inverse (using 'singular-value friendly' GAUSS pivot)
1083     
1084   Int_t nRank = 0;
1085   int iPivot;
1086   double vPivot = 0.;
1087   double eps = 0.00000000000001;
1088
1089   bool *bUnUsed = new bool[nGlo];
1090   double *diagV = new double[nGlo];
1091   double *rowMax = new double[nGlo];
1092   double *colMax = new double[nGlo];
1093
1094   double *temp = new double[nGlo];
1095
1096   for (Int_t i=0; i<nGlo; i++) {
1097     rowMax[i] = 0.0;
1098     colMax[i] = 0.0;
1099     bUnUsed[i] = true;
1100
1101     for (Int_t j=0; j<i; j++) {
1102       if (matV[j][i] == 0) {
1103         matV[j][i] = matV[i][j];
1104       }
1105     }
1106   }
1107   
1108   // Small loop for matrix equilibration (gives a better conditioning) 
1109
1110   for (Int_t i=0; i<nGlo; i++) {
1111     for (Int_t j=0; j<nGlo; j++) { 
1112       if (TMath::Abs(matV[i][j]) >= rowMax[i]) rowMax[i] = TMath::Abs(matV[i][j]); // Max elemt of row i
1113       if (TMath::Abs(matV[j][i]) >= colMax[i]) colMax[i] = TMath::Abs(matV[j][i]); // Max elemt of column i
1114     }
1115   }
1116
1117   for (Int_t i=0; i<nGlo; i++) {
1118     if (0.0 != rowMax[i]) rowMax[i] = 1./rowMax[i]; // Max elemt of row i
1119     if (0.0 != colMax[i]) colMax[i] = 1./colMax[i]; // Max elemt of column i
1120   }
1121
1122   for (Int_t i=0; i<nGlo; i++) {
1123     for (Int_t j=0; j<nGlo; j++) {
1124       matV[i][j] = TMath::Sqrt(rowMax[i])*matV[i][j]*TMath::Sqrt(colMax[j]); // Equilibrate the V matrix
1125     }
1126     diagV[i] = TMath::Abs(matV[i][i]); // save diagonal elem absolute values    
1127   }
1128
1129
1130   for (Int_t i=0; i<nGlo; i++) {
1131     vPivot = 0.0;
1132     iPivot = -1;
1133     
1134     for (Int_t j=0; j<nGlo; j++) { // First look for the pivot, ie max unused diagonal element       
1135       if (bUnUsed[j] && (TMath::Abs(matV[j][j])>TMath::Max(TMath::Abs(vPivot),eps*diagV[j]))) {    
1136         vPivot = matV[j][j];
1137         iPivot = j;
1138       }
1139     }
1140     
1141     if (iPivot >= 0) {   // pivot found          
1142       nRank++;
1143       bUnUsed[iPivot] = false; // This value is used
1144       vPivot = 1.0/vPivot;
1145       matV[iPivot][iPivot] = -vPivot; // Replace pivot by its inverse
1146       
1147       for (Int_t j=0; j<nGlo; j++) {      
1148         for (Int_t jj=0; jj<nGlo; jj++) {  
1149           if (j != iPivot && jj != iPivot) {// Other elements (!!! do them first as you use old matV[k][j]'s !!!)         
1150             matV[j][jj] = matV[j][jj] - vPivot*matV[j][iPivot]*matV[iPivot][jj];
1151           }
1152         }
1153       }
1154       
1155       for (Int_t j=0; j<nGlo; j++) {      
1156         if (j != iPivot) { // Pivot row or column elements 
1157           matV[j][iPivot] = matV[j][iPivot]*vPivot;     // Column
1158           matV[iPivot][j] = matV[iPivot][j]*vPivot;     // Line
1159         }
1160       }
1161     }
1162     else {  // No more pivot value (clear those elements)
1163       for (Int_t j=0; j<nGlo; j++) {
1164         if (bUnUsed[j]) {
1165           vecB[j] = 0.0;
1166
1167           for (Int_t k=0; k<nGlo; k++) {
1168             matV[j][k] = 0.0;
1169             matV[k][j] = 0.0;
1170           }
1171         }
1172       }
1173       break;  // No more pivots anyway, stop here
1174     }
1175   }
1176   
1177   for (Int_t i=0; i<nGlo; i++) {
1178     for (Int_t j=0; j<nGlo; j++) {
1179       matV[i][j] = TMath::Sqrt(colMax[i])*matV[i][j]*TMath::Sqrt(rowMax[j]); // Correct matrix V
1180     }
1181   }
1182   
1183   for (Int_t j=0; j<nGlo; j++) {
1184     temp[j] = 0.0;
1185     
1186     for (Int_t jj=0; jj<nGlo; jj++) { // Reverse matrix elements
1187       matV[j][jj] = -matV[j][jj];
1188       temp[j] += matV[j][jj]*vecB[jj];
1189     }           
1190   }
1191
1192   for (Int_t j=0; j<nGlo; j++) {
1193     vecB[j] = temp[j]; // The final result
1194   }
1195   
1196   delete [] temp;
1197   delete [] bUnUsed;
1198   delete [] diagV;
1199   delete [] rowMax;
1200   delete [] colMax;
1201
1202   return nRank;
1203 }
1204
1205 //
1206 // Same method but for local fit, so heavily simplified
1207 //
1208 int AliMillepede::SpmInv(double matV[][fgkMaxLocalPar], double vecB[], int nLoc)
1209 {
1210   ///  Obtain solution of a system of linear equations with symmetric matrix 
1211   ///  and the inverse (using 'singular-value friendly' GAUSS pivot)
1212
1213   Int_t nRank = 0;
1214   Int_t iPivot = -1;
1215   double vPivot = 0.;
1216   double eps = 0.0000000000001;
1217
1218   bool *bUnUsed = new bool[nLoc];
1219   double *diagV  = new double[nLoc];
1220   double *temp  = new double[nLoc];
1221   
1222   for (Int_t i=0; i<nLoc; i++) {
1223     bUnUsed[i] = true;
1224     diagV[i] = TMath::Abs(matV[i][i]);     // save diagonal elem absolute values
1225     for (Int_t j=0; j<i; j++) {
1226       matV[j][i] = matV[i][j] ;
1227     }
1228   }
1229         
1230   for (Int_t i=0; i<nLoc; i++) {
1231     vPivot = 0.0;
1232     iPivot = -1;
1233                 
1234     for (Int_t j=0; j<nLoc; j++) { // First look for the pivot, ie max unused diagonal element 
1235       if (bUnUsed[j] && (TMath::Abs(matV[j][j])>TMath::Max(TMath::Abs(vPivot),eps*diagV[j]))) {
1236         vPivot = matV[j][j];
1237         iPivot = j;
1238       }
1239     }
1240                 
1241     if (iPivot >= 0) {   // pivot found
1242       nRank++;
1243       bUnUsed[iPivot] = false;
1244       vPivot = 1.0/vPivot;
1245       matV[iPivot][iPivot] = -vPivot; // Replace pivot by its inverse
1246       
1247       for (Int_t j=0; j<nLoc; j++) {
1248         if (j != iPivot) {
1249           for (Int_t jj=0; jj<=j; jj++) {       
1250             if (jj != iPivot) {// Other elements (!!! do them first as you use old matV[k][j]'s !!!)
1251               matV[j][jj] = matV[j][jj] - vPivot*matV[j][iPivot]*matV[iPivot][jj];
1252               matV[jj][j] = matV[j][jj];
1253             }
1254           }
1255         }
1256       }
1257
1258       for (Int_t j=0; j<nLoc; j++) {      
1259         if (j != iPivot) {      // Pivot row or column elements         
1260           matV[j][iPivot] = matV[j][iPivot]*vPivot; // Column
1261           matV[iPivot][j] = matV[j][iPivot];
1262         }
1263       }
1264     }
1265     else { // No more pivot value (clear those elements)
1266       for (Int_t j=0; j<nLoc; j++) {
1267         if (bUnUsed[j]) {
1268           vecB[j] = 0.0;
1269           matV[j][j] = 0.0;
1270           for (Int_t k=0; k<j; k++) {     
1271             matV[j][k] = 0.0;
1272             matV[k][j] = 0.0;
1273           }
1274         }
1275       }
1276
1277       break;  // No more pivots anyway, stop here
1278     }
1279   }
1280
1281   for (Int_t j=0; j<nLoc; j++) {  
1282     temp[j] = 0.0;    
1283     for (Int_t jj=0; jj<nLoc; jj++) { // Reverse matrix elements
1284       matV[j][jj] = -matV[j][jj];
1285       temp[j] += matV[j][jj]*vecB[jj];
1286     }                   
1287   }
1288
1289   for (Int_t j=0; j<nLoc; j++) {
1290     vecB[j] = temp[j];
1291   }
1292
1293   delete [] bUnUsed;
1294   delete [] diagV;
1295   delete [] temp;
1296   
1297   return nRank;
1298 }
1299
1300
1301 /*
1302 -----------------------------------------------------------
1303   SPAVAT
1304 -----------------------------------------------------------
1305
1306   multiply symmetric N-by-N matrix from the left with general M-by-N
1307   matrix and from the right with the transposed of the same  general
1308   matrix  to  form  symmetric  M-by-M   matrix.
1309   
1310                                                        T
1311   CALL SPAVAT(V,A,W,N,M)      W   =   A   *   V   *   A
1312                            M*M     M*N     N*N     N*M
1313   
1314   where V = symmetric N-by-N matrix
1315         A = general N-by-M matrix
1316         W = symmetric M-by-M matrix
1317 -----------------------------------------------------------
1318 */
1319 Int_t AliMillepede::SpAVAt(double matV[][fgkMaxLocalPar], double matA[][fgkMaxLocalPar], double matW[][fgkMaxGlobalPar], int nLoc, int nGlo)
1320 {
1321   ///  multiply symmetric N-by-N matrix from the left with general M-by-N
1322   ///  matrix and from the right with the transposed of the same general
1323   ///  matrix to form symmetric M-by-M matrix.
1324
1325   for (Int_t i=0; i<nGlo; i++) {
1326     for (Int_t j=0; j<=i; j++) {  // Matrix w is symmetric
1327       matW[i][j] = 0.0;      // Reset final matrix                      
1328       for (Int_t k=0; k<nLoc; k++) {    
1329         matW[i][j] += matA[i][k]*matV[k][k]*matA[j][k];    // diagonal terms of v
1330         for (Int_t l=0; l<k; l++) {       
1331           matW[i][j] += matA[i][k]*matV[k][l]*matA[j][l];  // Use symmetric properties
1332           matW[i][j] += matA[i][l]*matV[k][l]*matA[j][k];  // of matrix v
1333         }
1334       }
1335       if (i!=j){
1336         matW[j][i] = matW[i][j]; // Matrix w is symmetric
1337       }
1338     }
1339   }
1340         
1341   return 1;
1342 }
1343
1344
1345 /*
1346 -----------------------------------------------------------
1347   SPAX
1348 -----------------------------------------------------------
1349
1350   multiply general M-by-N matrix A and N-vector X
1351  
1352   CALL  SPAX(A,X,Y,M,N)          Y =  A * X
1353                                  M   M*N  N
1354  
1355   where A = general M-by-N matrix (A11 A12 ... A1N  A21 A22 ...)
1356         X = N vector
1357         Y = M vector
1358 -----------------------------------------------------------
1359 */
1360 Int_t AliMillepede::SpAX(double matA[][fgkMaxLocalPar], double vecX[], double vecY[], int nCol, int nRow)
1361 {
1362   ///   multiply general M-by-N matrix A and N-vector X
1363   for (Int_t i=0; i<nRow; i++) {
1364     vecY[i] = 0.0;          // Reset final vector                       
1365     for (Int_t j=0; j<nCol; j++) {
1366       vecY[i] += matA[i][j]*vecX[j];  // fill the vector
1367     }
1368   }
1369         
1370   return 1;
1371 }
1372
1373
1374 /*
1375 -----------------------------------------------------------
1376   PRTGLO
1377 -----------------------------------------------------------
1378
1379   Print the final results into the logfile
1380
1381 -----------------------------------------------------------
1382 */
1383 Int_t AliMillepede::PrintGlobalParameters() const
1384 {
1385   ///  Print the final results into the logfile
1386   double lError = 0.;
1387   double lGlobalCor =0.;
1388         
1389   AliInfo("");
1390   AliInfo("   Result of fit for global parameters");
1391   AliInfo("   ===================================");
1392   AliInfo("    I       initial       final       differ        lastcor        error       gcor");
1393   AliInfo("-----------------------------------------------------------------------------------");
1394         
1395   for (int i=0; i<fNGlobalPar; i++) {
1396     lError = TMath::Sqrt(TMath::Abs(fMatCGlo[i][i]));
1397     if (fMatCGlo[i][i] < 0.0) lError = -lError;
1398     lGlobalCor = 0.0;
1399                 
1400     if (TMath::Abs(fMatCGlo[i][i]*fDiagCGlo[i]) > 0) {    
1401       lGlobalCor = TMath::Sqrt(TMath::Abs(1.0-1.0/(fMatCGlo[i][i]*fDiagCGlo[i])));
1402       AliInfo(Form("%d\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f",
1403                    i,fInitPar[i],fInitPar[i]+fDeltaPar[i],fDeltaPar[i],fVecBGlo[i],lError,lGlobalCor));
1404     }
1405     else {    
1406       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]));
1407     }
1408   }
1409   return 1;
1410 }
1411
1412
1413 /*
1414 ----------------------------------------------------------------
1415   CHI2DOFLIM:  return the limit in chi^2/nd for n sigmas stdev authorized
1416 ----------------------------------------------------------------
1417
1418   Only n=1, 2, and 3 are expected in input
1419 ----------------------------------------------------------------
1420 */
1421 double AliMillepede::Chi2DoFLim(int nSig, int nDoF)
1422 {
1423   /// return the limit in chi^2/nd for n sigmas stdev authorized
1424   int lNSig;
1425   double sn[3]        = {0.47523, 1.690140, 2.782170};
1426   double table[3][30] = {{1.0000, 1.1479, 1.1753, 1.1798, 1.1775, 1.1730, 1.1680, 1.1630,
1427                           1.1581, 1.1536, 1.1493, 1.1454, 1.1417, 1.1383, 1.1351, 1.1321,
1428                           1.1293, 1.1266, 1.1242, 1.1218, 1.1196, 1.1175, 1.1155, 1.1136,
1429                           1.1119, 1.1101, 1.1085, 1.1070, 1.1055, 1.1040},
1430                          {4.0000, 3.0900, 2.6750, 2.4290, 2.2628, 2.1415, 2.0481, 1.9736,
1431                           1.9124, 1.8610, 1.8171, 1.7791, 1.7457, 1.7161, 1.6897, 1.6658,
1432                           1.6442, 1.6246, 1.6065, 1.5899, 1.5745, 1.5603, 1.5470, 1.5346,
1433                           1.5230, 1.5120, 1.5017, 1.4920, 1.4829, 1.4742},
1434                          {9.0000, 5.9146, 4.7184, 4.0628, 3.6410, 3.3436, 3.1209, 2.9468,
1435                           2.8063, 2.6902, 2.5922, 2.5082, 2.4352, 2.3711, 2.3143, 2.2635,
1436                           2.2178, 2.1764, 2.1386, 2.1040, 2.0722, 2.0428, 2.0155, 1.9901,
1437                           1.9665, 1.9443, 1.9235, 1.9040, 1.8855, 1.8681}};
1438
1439   if (nDoF < 1) {
1440     return 0.0;
1441   }
1442   else {  
1443     lNSig = TMath::Max(1,TMath::Min(nSig,3));
1444
1445     if (nDoF <= 30) {    
1446       return table[lNSig-1][nDoF-1];
1447     }
1448     else { // approximation
1449       return ((sn[lNSig-1]+TMath::Sqrt(float(2*nDoF-3)))*
1450               (sn[lNSig-1]+TMath::Sqrt(float(2*nDoF-3))))/float(2*nDoF-2);
1451     }
1452   }
1453 }