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