]> git.uio.no Git - u/mrichter/AliRoot.git/blame - MUON/AliMillepede.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Updated version (Gustavo)
[u/mrichter/AliRoot.git] / MUON / AliMillepede.cxx
CommitLineData
010eb601 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>
090026bf 33#include <TMath.h>
010eb601 34
35#include "AliLog.h"
36
37#include "AliMillepede.h"
38
39//=============================================================================
8395bff1 40AliMillepede::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)
010eb601 65{
66 /// Standard constructor
67
010eb601 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//=============================================================================
81AliMillepede::~AliMillepede() {
82 /// Destructor
83
227a3364 84}
010eb601 85
86//=============================================================================
87Int_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*/
183Int_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*/
205Int_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*/
233Int_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*/
259Int_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*/
293Int_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*/
315Int_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*/
346Int_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*/
441Int_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++) {
1702b054 600 AliDebug(1,Form("%d / %.6f / %.6f", i, fVecBLoc[i], TMath::Sqrt(fMatCLoc[i][i])));
010eb601 601 }
602
603
604// Store the track params and errors
605
606 for (i=0; i<fNLocalPar; i++) {
607 localParams[2*i] = fVecBLoc[i];
1702b054 608 localParams[2*i+1] = TMath::Sqrt(TMath::Abs(fMatCLoc[i][i]));
010eb601 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,"");
1702b054 642// AliDebug(2,Form(". equation: measured value %.6f +/- %.6f", lMeas, 1.0/TMath::Sqrt(lWeight)));
010eb601 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)
1702b054 674 if (TMath::Abs(lMeas) >= fResCutInit && fIter <= 1) {
010eb601 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
1702b054 682 if (TMath::Abs(lMeas) >= fResCut && fIter > 1) {
010eb601 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*/
848Int_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]));
1702b054 941 AliDebug(1,Form("err = %.6f", TMath::Sqrt(TMath::Abs(fMatCGlo[i][i]))));
010eb601 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) {
1702b054 964 fChi2CutFactor = TMath::Sqrt(fChi2CutFactor);
010eb601 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];
1702b054 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]));
010eb601 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*/
1057Double_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) {
1702b054 1062 lErr = TMath::Sqrt(TMath::Abs(fMatCGlo[iPar][iPar]));
010eb601 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*/
1079int 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++) {
1702b054 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
010eb601 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++) {
1702b054 1124 matV[i][j] = TMath::Sqrt(rowMax[i])*matV[i][j]*TMath::Sqrt(colMax[j]); // Equilibrate the V matrix
010eb601 1125 }
1702b054 1126 diagV[i] = TMath::Abs(matV[i][i]); // save diagonal elem absolute values
010eb601 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
1702b054 1135 if (bUnUsed[j] && (TMath::Abs(matV[j][j])>std::max(TMath::Abs(vPivot),eps*diagV[j]))) {
010eb601 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++) {
1702b054 1179 matV[i][j] = TMath::Sqrt(colMax[i])*matV[i][j]*TMath::Sqrt(rowMax[j]); // Correct matrix V
010eb601 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//
1208int 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;
1702b054 1224 diagV[i] = TMath::Abs(matV[i][i]); // save diagonal elem absolute values
010eb601 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
1702b054 1235 if (bUnUsed[j] && (TMath::Abs(matV[j][j])>TMath::Max(TMath::Abs(vPivot),eps*diagV[j]))) {
010eb601 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*/
1319Int_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*/
1360Int_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*/
1383Int_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++) {
1702b054 1396 lError = TMath::Sqrt(TMath::Abs(fMatCGlo[i][i]));
010eb601 1397 if (fMatCGlo[i][i] < 0.0) lError = -lError;
1398 lGlobalCor = 0.0;
1399
1702b054 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])));
010eb601 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*/
1421double 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
1702b054 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);
010eb601 1451 }
1452 }
1453}