]>
Commit | Line | Data |
---|---|---|
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 | 40 | AliMillepede::AliMillepede() |
41 | : TObject(), | |
42 | fIndexLocEq(fgkMaxGlobalPar+fgkMaxLocalPar), | |
43 | fDerivLocEq(fgkMaxGlobalPar+fgkMaxLocalPar), | |
44 | fIndexAllEqs(1000*fgkMaxGlobalPar+fgkMaxLocalPar), | |
45 | fDerivAllEqs(1000*fgkMaxGlobalPar+fgkMaxLocalPar), | |
46 | fLocEqPlace(1000), | |
47 | fNIndexLocEq(0), | |
48 | fNDerivLocEq(0), | |
49 | fNIndexAllEqs(0), | |
50 | fNDerivAllEqs(0), | |
51 | fNLocEqPlace(0), | |
52 | fNLocalEquations(0), | |
53 | fResCutInit(0.), | |
54 | fResCut(0.), | |
55 | fChi2CutFactor(1.0), | |
56 | fChi2CutRef(1.0), | |
57 | fIter(0), | |
58 | fMaxIter(10), | |
59 | fNStdDev(3), | |
60 | fNGlobalConstraints(0), | |
61 | fNLocalFits(0), | |
62 | fNLocalFitsRejected(0), | |
63 | fNGlobalPar(0), | |
64 | fNLocalPar(0) | |
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 | //============================================================================= | |
81 | AliMillepede::~AliMillepede() { | |
82 | /// Destructor | |
83 | ||
227a3364 | 84 | } |
010eb601 | 85 | |
86 | //============================================================================= | |
87 | Int_t AliMillepede::InitMille(int nGlo, int nLoc, int lNStdDev, | |
88 | double lResCut, double lResCutInit) | |
89 | { | |
90 | /// Initialization of millepede | |
91 | AliDebug(1,""); | |
92 | AliDebug(1,"----------------------------------------------------"); | |
93 | AliDebug(1,""); | |
94 | AliDebug(1," Entering InitMille"); | |
95 | AliDebug(1,""); | |
96 | AliDebug(1,"-----------------------------------------------------"); | |
97 | AliDebug(1,""); | |
98 | ||
99 | fNGlobalConstraints = 0; | |
100 | fNLocalFits = 0; // Total number of local fits | |
101 | fNLocalFitsRejected = 0; // Total number of local fits rejected | |
102 | fChi2CutRef = 1.0; // Reference value for Chi^2/ndof cut | |
103 | ||
104 | AliMillepede::SetNLocalEquations(0); // Number of local fits (starts at 0) | |
105 | ||
106 | fIter = 0; // By default iterations are turned off, turned on by use SetIterations | |
107 | fMaxIter = 10; | |
108 | ||
109 | fResCutInit = lResCutInit; | |
110 | fResCut = lResCut; | |
111 | ||
112 | fNGlobalPar = nGlo; // Number of global derivatives | |
113 | fNLocalPar = nLoc; // Number of local derivatives | |
114 | fNStdDev = lNStdDev; // Number of StDev for local fit chisquare cut | |
115 | ||
116 | AliDebug(1,Form("Number of global parameters : %d", fNGlobalPar)); | |
117 | AliDebug(1,Form("Number of local parameters : %d", fNLocalPar)); | |
118 | AliDebug(1,Form("Number of standard deviations : %d", fNStdDev)); | |
119 | ||
120 | if (fNGlobalPar>fgkMaxGlobalPar || fNLocalPar>fgkMaxLocalPar) { | |
121 | AliDebug(1,"Two many parameters !!!!!"); | |
122 | return 0; | |
123 | } | |
124 | ||
125 | // Global parameters initializations | |
126 | for (int i=0; i<fNGlobalPar; i++) { | |
127 | fVecBGlo[i]=0.; | |
128 | fInitPar[i]=0.; | |
129 | fDeltaPar[i]=0.; | |
130 | fSigmaPar[i]=-1.; | |
131 | fIsNonLinear[i]=0; | |
132 | fGlo2CGLRow[i]=-1; | |
133 | fCGLRow2Glo[i]=-1; | |
134 | ||
135 | for (int j=0; j<fNGlobalPar;j++) { | |
136 | fMatCGlo[i][j]=0.; | |
137 | } | |
138 | } | |
139 | ||
140 | // Local parameters initializations | |
141 | ||
142 | for (int i=0; i<fNLocalPar; i++) { | |
143 | fVecBLoc[i]=0.; | |
144 | ||
145 | for (int j=0; j<fNLocalPar;j++) { | |
146 | fMatCLoc[i][j]=0.; | |
147 | } | |
148 | } | |
149 | ||
150 | // Then we fix all parameters... | |
151 | ||
152 | for (int j=0; j<fNGlobalPar; j++) { | |
153 | AliMillepede::SetParSigma(j,0.0); | |
154 | } | |
155 | ||
156 | fDerivLocEq.Reset(); fNDerivLocEq=0; | |
157 | fIndexLocEq.Reset(); fNIndexLocEq=0; | |
158 | ||
159 | fIndexAllEqs.Reset(); fNIndexAllEqs=0; | |
160 | fDerivAllEqs.Reset(); fNDerivAllEqs=0; | |
161 | fLocEqPlace.Reset(); fNLocEqPlace=0; | |
162 | ||
163 | AliDebug(1,""); | |
164 | AliDebug(1,"----------------------------------------------------"); | |
165 | AliDebug(1,""); | |
166 | AliDebug(1," InitMille has been successfully called!"); | |
167 | AliDebug(1,""); | |
168 | AliDebug(1,"-----------------------------------------------------"); | |
169 | AliDebug(1,""); | |
170 | ||
171 | return 1; | |
172 | } | |
173 | ||
174 | /* | |
175 | ----------------------------------------------------------- | |
176 | PARGLO: initialization of global parameters | |
177 | ----------------------------------------------------------- | |
178 | ||
179 | param = array of starting values | |
180 | ||
181 | ----------------------------------------------------------- | |
182 | */ | |
183 | Int_t AliMillepede::SetGlobalParameters(double *param) | |
184 | { | |
185 | /// initialization of global parameters | |
186 | for(Int_t iPar=0; iPar<fNGlobalPar; iPar++){ | |
187 | fInitPar[iPar] = param[iPar]; | |
188 | } | |
189 | ||
190 | return 1; | |
191 | } | |
192 | ||
193 | /* | |
194 | ----------------------------------------------------------- | |
195 | PARGLO: initialization of global parameters | |
196 | ----------------------------------------------------------- | |
197 | ||
198 | iPar = the index of the global parameter in the | |
199 | result array (equivalent to fDeltaPar[]). | |
200 | ||
201 | param = the starting value | |
202 | ||
203 | ----------------------------------------------------------- | |
204 | */ | |
205 | Int_t AliMillepede::SetGlobalParameter(int iPar, double param) | |
206 | { | |
207 | /// initialization of global parameter iPar | |
208 | if (iPar<0 || iPar>=fNGlobalPar) { | |
209 | return 0; | |
210 | } | |
211 | else { | |
212 | fInitPar[iPar] = param; | |
213 | } | |
214 | return 1; | |
215 | } | |
216 | ||
217 | ||
218 | /* | |
219 | ----------------------------------------------------------- | |
220 | PARSIG: define a constraint for a single global param | |
221 | param is 'encouraged' to vary within [-sigma;sigma] | |
222 | range | |
223 | ----------------------------------------------------------- | |
224 | ||
225 | iPar = the index of the global parameter in the | |
226 | result array (equivalent to fDeltaPar[]). | |
227 | ||
228 | sigma = value of the constraint (sigma <= 0. will | |
229 | mean that parameter is FIXED !!!) | |
230 | ||
231 | ----------------------------------------------------------- | |
232 | */ | |
233 | Int_t AliMillepede::SetParSigma(int iPar, double sigma) | |
234 | { | |
235 | /// Define a range [-sigma;sigma] where iPar is encourage to vary | |
236 | if (iPar>=fNGlobalPar) { | |
237 | return 0; | |
238 | } | |
239 | else { | |
240 | fSigmaPar[iPar] = sigma; | |
241 | } | |
242 | ||
243 | return 1; | |
244 | } | |
245 | ||
246 | ||
247 | /* | |
248 | ----------------------------------------------------------- | |
249 | NONLIN: set nonlinear flag for a single global param | |
250 | update of param durin iterations will not | |
251 | consider initial starting value | |
252 | ----------------------------------------------------------- | |
253 | ||
254 | iPar = the index of the global parameter in the | |
255 | result array (equivalent to fDeltaPar[]). | |
256 | ||
257 | ----------------------------------------------------------- | |
258 | */ | |
259 | Int_t AliMillepede::SetNonLinear(int iPar) | |
260 | { | |
261 | /// Set nonlinear flag for iPar | |
262 | if (iPar<0 || iPar>=fNGlobalPar) { | |
263 | return 0; | |
264 | } | |
265 | else { | |
266 | fIsNonLinear[iPar] = 1; | |
267 | } | |
268 | ||
269 | return 1; | |
270 | } | |
271 | ||
272 | ||
273 | /* | |
274 | ----------------------------------------------------------- | |
275 | INITUN: unit for iteration | |
276 | ----------------------------------------------------------- | |
277 | ||
278 | lChi2CutFac is used by Fitloc to define the Chi^2/ndof cut value | |
279 | ||
280 | A large cutfac value enables to take a wider range of tracks | |
281 | for first iterations, which might be useful if misalignments | |
282 | are large. | |
283 | ||
284 | As soon as cutfac differs from 0 iteration are requested. | |
285 | cutfac is then reduced, from one iteration to the other, | |
286 | and iterations are stopped when it reaches the value 1. | |
287 | ||
288 | At least one more iteration is often needed in order to remove | |
289 | tracks containing outliers. | |
290 | ||
291 | ----------------------------------------------------------- | |
292 | */ | |
293 | Int_t AliMillepede::SetIterations(double lChi2CutFac) | |
294 | { | |
295 | /// Number of iterations is calculated from lChi2CutFac | |
296 | fChi2CutFactor = TMath::Max(1.0, lChi2CutFac); | |
297 | ||
298 | AliInfo(Form("Initial cut factor is %f",fChi2CutFactor)); | |
299 | fIter = 1; // Initializes the iteration process | |
300 | return 1; | |
301 | } | |
302 | ||
303 | /* | |
304 | ----------------------------------------------------------- | |
305 | CONSTF: define a constraint equation in AliMillepede | |
306 | ----------------------------------------------------------- | |
307 | ||
308 | dercs = the row containing constraint equation | |
309 | derivatives (put into the final matrix) | |
310 | ||
311 | lLagMult = the lagrange multiplier value (sum of equation) | |
312 | ||
313 | ----------------------------------------------------------- | |
314 | */ | |
315 | Int_t AliMillepede::SetGlobalConstraint(double dercs[], double lLagMult) | |
316 | { | |
317 | /// Define a constraint equation | |
318 | if (fNGlobalConstraints>=fgkMaxGloCsts) { | |
319 | AliInfo("Too many constraints !!!"); | |
320 | return 0; | |
321 | } | |
322 | ||
323 | for (int i=0; i<fNGlobalPar; i++) { | |
324 | fMatDerConstr[fNGlobalConstraints][i] = dercs[i]; | |
325 | } | |
326 | ||
327 | fLagMult[fNGlobalConstraints] = lLagMult; | |
328 | fNGlobalConstraints++ ; | |
329 | AliInfo(Form("Number of constraints increased to %d",fNGlobalConstraints)); | |
330 | return 1; | |
331 | } | |
332 | ||
333 | ||
334 | /* | |
335 | ----------------------------------------------------------- | |
336 | EQULOC: write ONE equation in the matrices | |
337 | ----------------------------------------------------------- | |
338 | ||
339 | dergb[1..fNGlobalPar] = global parameters derivatives | |
340 | derlc[1..fNLocalPar] = local parameters derivatives | |
341 | rmeas = measured value | |
342 | sigma = error on measured value (nothing to do with SetParSigma!!!) | |
343 | ||
344 | ----------------------------------------------------------- | |
345 | */ | |
346 | Int_t AliMillepede::SetLocalEquation(double dergb[], double derlc[], double lMeas, double lSigma) | |
347 | { | |
348 | /// Write one local equation | |
349 | if (lSigma<=0.0) { // If parameter is fixed, then no equation | |
350 | for (int i=0; i<fNLocalPar; i++) { | |
351 | derlc[i] = 0.0; | |
352 | } | |
353 | for (int i=0; i<fNGlobalPar; i++) { | |
354 | dergb[i] = 0.0; | |
355 | } | |
356 | return 1; | |
357 | } | |
358 | ||
359 | // Serious equation, initialize parameters | |
360 | double lWeight = 1.0/(lSigma*lSigma); | |
361 | int iLocFirst = -1; | |
362 | int iLocLast = -1; | |
363 | int iGlobFirst = -1; | |
364 | int iGlobLast = -1; | |
365 | ||
366 | for (int i=0; i<fNLocalPar; i++) { // Retrieve local param interesting indices | |
367 | if (derlc[i]!=0.0) { | |
368 | if (iLocFirst == -1) { | |
369 | iLocFirst = i; // first index | |
370 | } | |
371 | iLocLast = i; // last index | |
372 | } | |
373 | } | |
374 | AliDebug(2,Form("%d / %d",iLocFirst, iLocLast)); | |
375 | ||
376 | for (int i=0; i<fNGlobalPar; i++) { // Idem for global parameters | |
377 | if (dergb[i]!=0.0) { | |
378 | if (iGlobFirst == -1) { | |
379 | iGlobFirst = i; // first index | |
380 | } | |
381 | iGlobLast = i; // last index | |
382 | } | |
383 | } | |
384 | AliDebug(2,Form("%d / %d",iGlobFirst,iGlobLast)); | |
385 | ||
386 | if (fNIndexLocEq==fIndexLocEq.GetSize()) fIndexLocEq.Set(2*fNIndexLocEq); | |
387 | fIndexLocEq.AddAt(-1,fNIndexLocEq++); | |
388 | if (fNDerivLocEq==fDerivLocEq.GetSize()) fDerivLocEq.Set(2*fNDerivLocEq); | |
389 | fDerivLocEq.AddAt(lMeas,fNDerivLocEq++); | |
390 | ||
391 | ||
392 | for (int i=iLocFirst; i<=iLocLast; i++) { // Store interesting local parameters | |
393 | if (derlc[i]!=0.0) { | |
394 | if (fNIndexLocEq==fIndexLocEq.GetSize()) fIndexLocEq.Set(2*fNIndexLocEq); | |
395 | fIndexLocEq.AddAt(i,fNIndexLocEq++); | |
396 | if (fNDerivLocEq==fDerivLocEq.GetSize()) fDerivLocEq.Set(2*fNDerivLocEq); | |
397 | fDerivLocEq.AddAt(derlc[i],fNDerivLocEq++); | |
398 | derlc[i] = 0.0; | |
399 | } | |
400 | } | |
401 | ||
402 | if (fNIndexLocEq==fIndexLocEq.GetSize()) fIndexLocEq.Set(2*fNIndexLocEq); | |
403 | fIndexLocEq.AddAt(-1,fNIndexLocEq++); | |
404 | if (fNDerivLocEq==fDerivLocEq.GetSize()) fDerivLocEq.Set(2*fNDerivLocEq); | |
405 | fDerivLocEq.AddAt(lWeight,fNDerivLocEq++); | |
406 | ||
407 | for (int i=iGlobFirst; i<=iGlobLast; i++) { // Store interesting global parameters | |
408 | if (dergb[i]!=0.0) { | |
409 | if (fNIndexLocEq==fIndexLocEq.GetSize()) fIndexLocEq.Set(2*fNIndexLocEq); | |
410 | fIndexLocEq.AddAt(i,fNIndexLocEq++); | |
411 | if (fNDerivLocEq==fDerivLocEq.GetSize()) fDerivLocEq.Set(2*fNDerivLocEq); | |
412 | fDerivLocEq.AddAt(dergb[i],fNDerivLocEq++); | |
413 | dergb[i] = 0.0; | |
414 | } | |
415 | } | |
416 | ||
417 | AliDebug(2,Form("Out Equloc -- NST = %d",fNDerivLocEq)); | |
418 | return 1; | |
419 | } | |
420 | ||
421 | /* | |
422 | ----------------------------------------------------------- | |
423 | FITLOC: perform local params fit, once all the equations | |
424 | have been written by EquLoc | |
425 | ----------------------------------------------------------- | |
426 | ||
427 | iFit = number of the fit, it is used to store | |
428 | fit parameters and then retrieve them | |
429 | for iterations (via FINDEXALLEQS and FDERIVALLEQS) | |
430 | ||
431 | localParams = contains the fitted track parameters and | |
432 | related errors | |
433 | ||
434 | bSingleFit = is an option, if it is set to 1, we don't | |
435 | perform the last loop. It is used to update | |
436 | the track parameters without modifying global | |
437 | matrices | |
438 | ||
439 | ----------------------------------------------------------- | |
440 | */ | |
1dd0ce54 | 441 | Int_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 | */ | |
848 | Int_t AliMillepede::GlobalFit(double par[], double error[], double pull[]) | |
849 | { | |
850 | /// perform global parameters fit once all the local equations have been fitted | |
851 | int i, j; | |
852 | int nVar = 0; | |
853 | int nGloFix = 0; | |
854 | double lConstraint; | |
855 | ||
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 | 1058 | Double_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 | */ | |
1080 | int 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 | // | |
1209 | int 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 | */ | |
1320 | Int_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 | */ | |
1361 | Int_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 | */ | |
1384 | Int_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 | */ | |
1422 | double 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 | } |