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