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