1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 //-----------------------------------------------------------------------------
19 /// \class AliMillepede
20 /// Detector independent alignment class
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
28 /// \author Javier Castillo
29 //-----------------------------------------------------------------------------
36 #include "AliMillepede.h"
38 //=============================================================================
39 AliMillepede::AliMillepede()
41 fIndexLocEq(fgkMaxGlobalPar+fgkMaxLocalPar),
42 fDerivLocEq(fgkMaxGlobalPar+fgkMaxLocalPar),
43 fIndexAllEqs(1000*fgkMaxGlobalPar+fgkMaxLocalPar),
44 fDerivAllEqs(1000*fgkMaxGlobalPar+fgkMaxLocalPar),
59 fNGlobalConstraints(0),
61 fNLocalFitsRejected(0),
65 /// Standard constructor
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");
79 //=============================================================================
80 AliMillepede::~AliMillepede() {
85 //=============================================================================
86 Int_t AliMillepede::InitMille(int nGlo, int nLoc, int lNStdDev,
87 double lResCut, double lResCutInit)
89 /// Initialization of millepede
91 AliDebug(1,"----------------------------------------------------");
93 AliDebug(1," Entering InitMille");
95 AliDebug(1,"-----------------------------------------------------");
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
103 AliMillepede::SetNLocalEquations(0); // Number of local fits (starts at 0)
105 fIter = 0; // By default iterations are turned off, turned on by use SetIterations
108 fResCutInit = lResCutInit;
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
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));
119 if (fNGlobalPar>fgkMaxGlobalPar || fNLocalPar>fgkMaxLocalPar) {
120 AliDebug(1,"Two many parameters !!!!!");
124 // Global parameters initializations
125 for (int i=0; i<fNGlobalPar; i++) {
134 for (int j=0; j<fNGlobalPar;j++) {
139 // Local parameters initializations
141 for (int i=0; i<fNLocalPar; i++) {
144 for (int j=0; j<fNLocalPar;j++) {
149 // Then we fix all parameters...
151 for (int j=0; j<fNGlobalPar; j++) {
152 AliMillepede::SetParSigma(j,0.0);
155 fDerivLocEq.Reset(); fNDerivLocEq=0;
156 fIndexLocEq.Reset(); fNIndexLocEq=0;
158 fIndexAllEqs.Reset(); fNIndexAllEqs=0;
159 fDerivAllEqs.Reset(); fNDerivAllEqs=0;
160 fLocEqPlace.Reset(); fNLocEqPlace=0;
163 AliDebug(1,"----------------------------------------------------");
165 AliDebug(1," InitMille has been successfully called!");
167 AliDebug(1,"-----------------------------------------------------");
174 -----------------------------------------------------------
175 PARGLO: initialization of global parameters
176 -----------------------------------------------------------
178 param = array of starting values
180 -----------------------------------------------------------
182 Int_t AliMillepede::SetGlobalParameters(double *param)
184 /// initialization of global parameters
185 for(Int_t iPar=0; iPar<fNGlobalPar; iPar++){
186 fInitPar[iPar] = param[iPar];
193 -----------------------------------------------------------
194 PARGLO: initialization of global parameters
195 -----------------------------------------------------------
197 iPar = the index of the global parameter in the
198 result array (equivalent to fDeltaPar[]).
200 param = the starting value
202 -----------------------------------------------------------
204 Int_t AliMillepede::SetGlobalParameter(int iPar, double param)
206 /// initialization of global parameter iPar
207 if (iPar<0 || iPar>=fNGlobalPar) {
211 fInitPar[iPar] = param;
218 -----------------------------------------------------------
219 PARSIG: define a constraint for a single global param
220 param is 'encouraged' to vary within [-sigma;sigma]
222 -----------------------------------------------------------
224 iPar = the index of the global parameter in the
225 result array (equivalent to fDeltaPar[]).
227 sigma = value of the constraint (sigma <= 0. will
228 mean that parameter is FIXED !!!)
230 -----------------------------------------------------------
232 Int_t AliMillepede::SetParSigma(int iPar, double sigma)
234 /// Define a range [-sigma;sigma] where iPar is encourage to vary
235 if (iPar>=fNGlobalPar) {
239 fSigmaPar[iPar] = sigma;
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 -----------------------------------------------------------
253 iPar = the index of the global parameter in the
254 result array (equivalent to fDeltaPar[]).
256 -----------------------------------------------------------
258 Int_t AliMillepede::SetNonLinear(int iPar)
260 /// Set nonlinear flag for iPar
261 if (iPar<0 || iPar>=fNGlobalPar) {
265 fIsNonLinear[iPar] = 1;
273 -----------------------------------------------------------
274 INITUN: unit for iteration
275 -----------------------------------------------------------
277 lChi2CutFac is used by Fitloc to define the Chi^2/ndof cut value
279 A large cutfac value enables to take a wider range of tracks
280 for first iterations, which might be useful if misalignments
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.
287 At least one more iteration is often needed in order to remove
288 tracks containing outliers.
290 -----------------------------------------------------------
292 Int_t AliMillepede::SetIterations(double lChi2CutFac)
294 /// Number of iterations is calculated from lChi2CutFac
295 fChi2CutFactor = TMath::Max(1.0, lChi2CutFac);
297 AliInfo(Form("Initial cut factor is %f",fChi2CutFactor));
298 fIter = 1; // Initializes the iteration process
303 -----------------------------------------------------------
304 CONSTF: define a constraint equation in AliMillepede
305 -----------------------------------------------------------
307 dercs = the row containing constraint equation
308 derivatives (put into the final matrix)
310 lLagMult = the lagrange multiplier value (sum of equation)
312 -----------------------------------------------------------
314 Int_t AliMillepede::SetGlobalConstraint(double dercs[], double lLagMult)
316 /// Define a constraint equation
317 if (fNGlobalConstraints>=fgkMaxGloCsts) {
318 AliInfo("Too many constraints !!!");
322 for (int i=0; i<fNGlobalPar; i++) {
323 fMatDerConstr[fNGlobalConstraints][i] = dercs[i];
326 fLagMult[fNGlobalConstraints] = lLagMult;
327 fNGlobalConstraints++ ;
328 AliInfo(Form("Number of constraints increased to %d",fNGlobalConstraints));
334 -----------------------------------------------------------
335 EQULOC: write ONE equation in the matrices
336 -----------------------------------------------------------
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!!!)
343 -----------------------------------------------------------
345 Int_t AliMillepede::SetLocalEquation(double dergb[], double derlc[], double lMeas, double lSigma)
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++) {
352 for (int i=0; i<fNGlobalPar; i++) {
358 // Serious equation, initialize parameters
359 double lWeight = 1.0/(lSigma*lSigma);
365 for (int i=0; i<fNLocalPar; i++) { // Retrieve local param interesting indices
367 if (iLocFirst == -1) {
368 iLocFirst = i; // first index
370 iLocLast = i; // last index
373 AliDebug(2,Form("%d / %d",iLocFirst, iLocLast));
375 for (int i=0; i<fNGlobalPar; i++) { // Idem for global parameters
377 if (iGlobFirst == -1) {
378 iGlobFirst = i; // first index
380 iGlobLast = i; // last index
383 AliDebug(2,Form("%d / %d",iGlobFirst,iGlobLast));
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++);
391 for (int i=iLocFirst; i<=iLocLast; i++) { // Store interesting local parameters
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++);
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++);
406 for (int i=iGlobFirst; i<=iGlobLast; i++) { // Store interesting global parameters
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++);
416 AliDebug(2,Form("Out Equloc -- NST = %d",fNDerivLocEq));
421 -----------------------------------------------------------
422 FITLOC: perform local params fit, once all the equations
423 have been written by EquLoc
424 -----------------------------------------------------------
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)
430 localParams = contains the fitted track parameters and
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
438 -----------------------------------------------------------
440 Int_t AliMillepede::LocalFit(int iFit, double localParams[], bool bSingleFit)
442 /// Perform local parameters fit once all the local equations have been set
444 // Few initializations
447 int iIdx, jIdx, kIdx;
448 int iIdxIdx, jIdxIdx;
458 double lWeight = 0.0;
461 double lRedChi2 = 0.0;
462 double lChi2Cut = 0.0;
465 int nEqTerms = fNDerivLocEq; // Number of terms (local + global derivatives,
466 // measurement and weight) involved in this local equation
469 // Fill the track store at first pass
471 if (fIter < 2 && !bSingleFit) { // Do it only once
472 AliDebug(1,Form("Store equation no: %d", iFit));
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++);
480 if (fNLocEqPlace==fLocEqPlace.GetSize()) fLocEqPlace.Set(2*fNLocEqPlace);
481 fLocEqPlace.AddAt(fNIndexAllEqs,fNLocEqPlace++);
484 AliDebug(2,Form("FLocEqPlace size = %d",fLocEqPlace[iFit]));
485 AliDebug(2,Form("FIndexAllEqs size = %d",fNIndexAllEqs));
489 for (i=0; i<fNLocalPar; i++) { // reset local params
492 for (j=0; j<fNLocalPar; j++) {
493 fMatCLoc[i][j] = 0.0;
497 for (i=0; i<fNGlobalPar; i++) {
498 fGlo2CGLRow[i] = -1; // reset mixed params
504 LOOPS : HOW DOES IT WORKS ?
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:
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.
516 We took them and store them into matrices.
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)
527 // FIRST LOOP : local track fit
531 // fIndexLocEq.push_back(-1);
532 if (fNIndexLocEq==fIndexLocEq.GetSize()) fIndexLocEq.Set(2*fNIndexLocEq);
533 fIndexLocEq.AddAt(-1,fNIndexLocEq++);
535 while (iEqTerm <= nEqTerms) {
536 if (fIndexLocEq[iEqTerm] == -1) {
537 if (iMeas == -1) { // First -1 : lMeas
539 iLocFirst = iEqTerm+1;
541 else if (iWeight == 0) { // Second -1 : weight
543 iLocLast = iEqTerm-1;
544 iGloFirst = iEqTerm+1;
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));
553 // Now suppress the global part (only relevant with iterations)
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
562 lMeas -= fDerivLocEq[i]*(fDeltaPar[iIdx]); // nonlinear parameter
564 // AliDebug(2,Form("lMeas after global stuff removal = %f", lMeas));
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]));
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]));
579 iEqTerm--; // end of one equation is the beginning of next
580 } // End of "end of equation" operations
581 } // End of loop on equation
583 } // End of loop on all equations used in the fit
587 // Local params matrix is completed, now invert to solve...
590 Int_t nRank = AliMillepede::SpmInv(fMatCLoc, fVecBLoc, fNLocalPar);
591 // nRank is the number of nonzero diagonal elements
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)");
598 for (i=0; i<fNLocalPar; i++) {
599 AliDebug(1,Form("%d / %.6f / %.6f", i, fVecBLoc[i], TMath::Sqrt(fMatCLoc[i][i])));
603 // Store the track params and errors
605 for (i=0; i<fNLocalPar; i++) {
606 localParams[2*i] = fVecBLoc[i];
607 localParams[2*i+1] = TMath::Sqrt(TMath::Abs(fMatCLoc[i][i]));
612 // SECOND LOOP : residual calculation
619 while (iEqTerm <= nEqTerms) {
620 if (fIndexLocEq[iEqTerm] == -1) {
621 if (iMeas == -1) { // First -1 : lMeas
623 iLocFirst = iEqTerm+1;
625 else if (iWeight == 0) { // Second -1 : weight
627 iLocLast = iEqTerm-1;
628 iGloFirst = iEqTerm+1;
630 else { // Third -1 : end of equation; start of next
631 iGloLast = iEqTerm-1;
632 lMeas = fDerivLocEq[iMeas];
633 lWeight = fDerivLocEq[iWeight];
635 // Print all (for debugging purposes)
637 // int nDerLoc = iLocLast-iLocFirst+1; // Number of local derivatives involved
638 // int nDerGlo = iGloLast-iGloFirst+1; // Number of global derivatives involved
641 // AliDebug(2,Form(". equation: measured value %.6f +/- %.6f", lMeas, 1.0/TMath::Sqrt(lWeight)));
642 // AliDebug(2,Form("Number of derivatives (global, local): %d, %d",nDerGlo,nDerLoc));
643 // AliDebug(2,"Global derivatives are: (index/derivative/parvalue) ");
645 // for (i=iGloFirst; i<=iGloLast; i++) {
646 // AliDebug(2,Form("%d / %.6f / %.6f",fIndexLocEq[i],fDerivLocEq[i],fInitPar[fIndexLocEq[i]]));
649 // AliDebug(2,"Local derivatives are: (index/derivative) ");
651 // for (i=(ja+1); i<jb; i++) {AliDebug(2,Form("%d / %.6f",fIndexLocEq[i], fDerivLocEq[i]));}
653 // Now suppress local and global parts to LMEAS;
655 // First the local part
656 for (i=iLocFirst; i<=iLocLast; i++) {
657 iIdx = fIndexLocEq[i];
658 lMeas -= fDerivLocEq[i]*fVecBLoc[iIdx];
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
666 lMeas -= fDerivLocEq[i]*(fDeltaPar[iIdx]); // nonlinear parameter
669 // lMeas contains now the residual value
670 AliDebug(2,Form("Residual value : %.6f", lMeas));
672 // reject the track if lMeas is too important (outlier)
673 if (TMath::Abs(lMeas) >= fResCutInit && fIter <= 1) {
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;
681 if (TMath::Abs(lMeas) >= fResCut && fIter > 1) {
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;
689 lChi2 += lWeight*lMeas*lMeas ; // total chi^2
690 nEq++; // number of equations
694 } // End of "end of equation" operations
695 } // End of loop on equation
697 } // End of loop on all equations used in the fit
702 AliDebug(1,Form("Final chi square / degrees of freedom %.2f / %d",lChi2, nDoF));
704 if (nDoF > 0) lRedChi2 = lChi2/float(nDoF); // Chi^2/dof
708 if (fNStdDev != 0 && nDoF > 0 && !bSingleFit) // Chisquare cut
710 lChi2Cut = AliMillepede::Chi2DoFLim(fNStdDev, nDoF)*fChi2CutFactor;
712 AliDebug(1,Form("Reject if Chisq/Ndf = %.4f > %.4f",lRedChi2,lChi2Cut));
714 if (lRedChi2 > lChi2Cut) // Reject the track if too much...
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;
724 if (bSingleFit) // Stop here if just updating the track parameters
726 fIndexLocEq.Reset(); fNIndexLocEq=0; // Reset store for the next track
727 fDerivLocEq.Reset(); fNDerivLocEq=0;
732 // THIRD LOOP: local operations are finished, track is accepted
733 // We now update the global parameters (other matrices)
740 while (iEqTerm <= nEqTerms)
742 if (fIndexLocEq[iEqTerm] == -1)
744 if (iMeas == -1) { // First -1 : lMeas
746 iLocFirst = iEqTerm+1;
748 else if (iWeight == 0) { // Second -1 : weight
750 iLocLast = iEqTerm-1;
751 iGloFirst = iEqTerm+1;
753 else { // Third -1 : end of equation; start of next
754 iGloLast = iEqTerm-1;
755 lMeas = fDerivLocEq[iMeas];
756 lWeight = fDerivLocEq[iWeight];
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
764 lMeas -= fDerivLocEq[i]*(fDeltaPar[iIdx]); // nonlinear parameter
767 for (i=iGloFirst; i<=iGloLast; i++) {
768 iIdx = fIndexLocEq[i]; // Global param indice (the matrix line)
770 fVecBGlo[iIdx] += lWeight*lMeas*fDerivLocEq[i];
771 // AliDebug(2,Form("fVecBGlo[%d] = %.6f", j, fVecBGlo[j] ));
773 // First of all, the global/global terms (exactly like local matrix)
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]));
781 // Now we have also rectangular matrices containing global/local terms.
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
788 fGlo2CGLRow[iIdx] = nGloInFit;
789 fCGLRow2Glo[nGloInFit] = iIdx;
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]));
804 } // End of "end of equation" operations
805 } // End of loop on equation
807 } // End of loop on all equations used in the fit
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);
813 for (iIdxIdx=0; iIdxIdx<nGloInFit; iIdxIdx++) {
814 iIdx = fCGLRow2Glo[iIdxIdx];
815 fVecBGlo[iIdx] -= fVecBGloCorr[iIdxIdx];
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];
824 fIndexLocEq.Reset(); fNIndexLocEq=0; // Reset store for the next track
825 fDerivLocEq.Reset(); fNDerivLocEq=0;
832 -----------------------------------------------------------
833 MAKEGLOBALFIT: perform global params fit, once all the 'tracks'
834 have been fitted by FitLoc
835 -----------------------------------------------------------
837 par[] = array containing the computed global
838 parameters (the misalignment constants)
840 error[] = array containing the error on global
841 parameters (estimated by AliMillepede)
843 pull[] = array containing the corresponding pulls
845 -----------------------------------------------------------
847 Int_t AliMillepede::GlobalFit(double par[], double error[], double pull[])
849 /// perform global parameters fit once all the local equations have been fitted
857 double localPars[2*fgkMaxLocalPar];
859 int nLocFitsGood = 0;
863 AliInfo("..... Making global fit .....");
865 nLocFitsTot = AliMillepede::GetNLocalEquations();
867 while (fIter < fMaxIter) // Iteration for the final loop
870 nLocFits = AliMillepede::GetNLocalEquations();
871 AliInfo(Form("...using %d local fits...",nLocFits));
873 // Start by saving the diagonal elements
875 for (i=0; i<fNGlobalPar; i++) {
876 fDiagCGlo[i] = fMatCGlo[i][i];
879 // Then we retrieve the different constraints: fixed parameter or global equation
881 nGloFix = 0; // First look at the fixed global params
883 for (i=0; i<fNGlobalPar; i++) {
884 if (fSigmaPar[i] <= 0.0) { // fixed global param
886 for (j=0; j<fNGlobalPar; j++) {
887 fMatCGlo[i][j] = 0.0; // Reset row and column
888 fMatCGlo[j][i] = 0.0;
892 fMatCGlo[i][i] += 1.0/(fSigmaPar[i]*fSigmaPar[i]);
896 nVar = fNGlobalPar; // Current number of equations
897 AliDebug(1,Form("Number of constraint equations : %d", fNGlobalConstraints));
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]);
907 fMatCGlo[nVar][nVar] = 0.0;
908 fVecBGlo[nVar] = float(nLocFits)*lConstraint;
913 // Intended to compute the final global chisquare
915 double lFinalCor = 0.0;
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]);
929 AliInfo(Form(" Final coeff is %.6f",lFinalCor));
930 AliInfo(Form(" Final NDOFs = %d", fNGlobalPar));
932 // The final matrix inversion
934 Int_t nRank = AliMillepede::SpmInv(fMatCGlo, fVecBGlo, nVar);
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]));
940 AliDebug(1,Form("err = %.6f", TMath::Sqrt(TMath::Abs(fMatCGlo[i][i]))));
942 step[i] = fVecBGlo[i];
944 if (fIter == 1) error[i] = fMatCGlo[i][i]; // Unfitted error
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));
951 AliInfo(Form("Total : %d local fits, %d rejected.", fNLocalFits, fNLocalFitsRejected));
952 if (fIter == 0) break; // No iterations set
955 if (fIter == fMaxIter) break; // End of story
957 // Reinitialize parameters for iteration
960 fNLocalFitsRejected = 0;
962 if (fChi2CutFactor != fChi2CutRef) {
963 fChi2CutFactor = TMath::Sqrt(fChi2CutFactor);
964 if (fChi2CutFactor < 1.2*fChi2CutRef) {
965 fChi2CutFactor = fChi2CutRef;
966 fIter = fMaxIter - 1; // Last iteration
969 AliInfo(Form("Iteration %d with cut factor %.2f", fIter, fChi2CutFactor));
971 // Reset global variables
973 for (i=0; i<nVar; i++) {
975 for (j=0; j<nVar; j++) {
976 fMatCGlo[i][j] = 0.0;
981 // We start a new iteration
984 // First we read the stores for retrieving the local params
988 for (i=0; i<nLocFitsTot; i++) {
992 (i>0) ? iEqFirst = fLocEqPlace[i-1] : iEqFirst = 0;
993 iEqLast = fLocEqPlace[i];
995 AliDebug(2,Form("Track %d : ",i));
996 AliDebug(2,Form("Starts at %d", iEqFirst));
997 AliDebug(2,Form("Ends at %d",iEqLast));
999 if (fIndexAllEqs[iEqFirst] != -999) { // Fit is still OK
1000 fIndexLocEq.Reset(); fNIndexLocEq=0;
1001 fDerivLocEq.Reset(); fNDerivLocEq=0;
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++);
1010 for (j=0; j<2*fNLocalPar; j++) {
1014 // Int_t sc = AliMillepede::LocalFit(i,localPars,0);
1015 // (sc) ? nLocFitsGood++ : fIndexAllEqs[iEqFirst] = -999;
1016 AliMillepede::LocalFit(i,localPars,0);
1019 } // End of loop on fits
1021 AliMillepede::SetNLocalEquations(nLocFitsGood);
1023 } // End of iteration loop
1025 AliMillepede::PrintGlobalParameters(); // Print the final results
1027 for (j=0; j<fNGlobalPar; j++) {
1028 par[j] = fInitPar[j]+fDeltaPar[j];
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]));
1034 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.");
1047 -----------------------------------------------------------
1048 ERRPAR: return error for parameter iPar
1049 -----------------------------------------------------------
1051 iPar = the index of the global parameter in the
1052 result array (equivalent to fDeltaPar[]).
1054 -----------------------------------------------------------
1056 Double_t AliMillepede::GetParError(Int_t iPar) const
1058 /// return error for parameter iPar
1059 Double_t lErr = -1.;
1060 if (iPar>=0 && iPar<fNGlobalPar) {
1061 lErr = TMath::Sqrt(TMath::Abs(fMatCGlo[iPar][iPar]));
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 -----------------------------------------------------------
1073 Solve the equation : V * X = B
1075 V is replaced by inverse matrix and B by X, the solution vector
1076 -----------------------------------------------------------
1078 int AliMillepede::SpmInv(double matV[][fgkMaxGloPC], double vecB[], int nGlo)
1080 /// Obtain solution of a system of linear equations with symmetric matrix
1081 /// and the inverse (using 'singular-value friendly' GAUSS pivot)
1086 double eps = 0.00000000000001;
1088 bool *bUnUsed = new bool[nGlo];
1089 double *diagV = new double[nGlo];
1090 double *rowMax = new double[nGlo];
1091 double *colMax = new double[nGlo];
1093 double *temp = new double[nGlo];
1095 for (Int_t i=0; i<nGlo; i++) {
1100 for (Int_t j=0; j<i; j++) {
1101 if (matV[j][i] == 0) {
1102 matV[j][i] = matV[i][j];
1107 // Small loop for matrix equilibration (gives a better conditioning)
1109 for (Int_t i=0; i<nGlo; i++) {
1110 for (Int_t j=0; j<nGlo; j++) {
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
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
1121 for (Int_t i=0; i<nGlo; i++) {
1122 for (Int_t j=0; j<nGlo; j++) {
1123 matV[i][j] = TMath::Sqrt(rowMax[i])*matV[i][j]*TMath::Sqrt(colMax[j]); // Equilibrate the V matrix
1125 diagV[i] = TMath::Abs(matV[i][i]); // save diagonal elem absolute values
1129 for (Int_t i=0; i<nGlo; i++) {
1133 for (Int_t j=0; j<nGlo; j++) { // First look for the pivot, ie max unused diagonal element
1134 if (bUnUsed[j] && (TMath::Abs(matV[j][j])>std::max(TMath::Abs(vPivot),eps*diagV[j]))) {
1135 vPivot = matV[j][j];
1140 if (iPivot >= 0) { // pivot found
1142 bUnUsed[iPivot] = false; // This value is used
1143 vPivot = 1.0/vPivot;
1144 matV[iPivot][iPivot] = -vPivot; // Replace pivot by its inverse
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];
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
1161 else { // No more pivot value (clear those elements)
1162 for (Int_t j=0; j<nGlo; j++) {
1166 for (Int_t k=0; k<nGlo; k++) {
1172 break; // No more pivots anyway, stop here
1176 for (Int_t i=0; i<nGlo; i++) {
1177 for (Int_t j=0; j<nGlo; j++) {
1178 matV[i][j] = TMath::Sqrt(colMax[i])*matV[i][j]*TMath::Sqrt(rowMax[j]); // Correct matrix V
1182 for (Int_t j=0; j<nGlo; j++) {
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];
1191 for (Int_t j=0; j<nGlo; j++) {
1192 vecB[j] = temp[j]; // The final result
1205 // Same method but for local fit, so heavily simplified
1207 int AliMillepede::SpmInv(double matV[][fgkMaxLocalPar], double vecB[], int nLoc)
1209 /// Obtain solution of a system of linear equations with symmetric matrix
1210 /// and the inverse (using 'singular-value friendly' GAUSS pivot)
1215 double eps = 0.0000000000001;
1217 bool *bUnUsed = new bool[nLoc];
1218 double *diagV = new double[nLoc];
1219 double *temp = new double[nLoc];
1221 for (Int_t i=0; i<nLoc; i++) {
1223 diagV[i] = TMath::Abs(matV[i][i]); // save diagonal elem absolute values
1224 for (Int_t j=0; j<i; j++) {
1225 matV[j][i] = matV[i][j] ;
1229 for (Int_t i=0; i<nLoc; i++) {
1233 for (Int_t j=0; j<nLoc; j++) { // First look for the pivot, ie max unused diagonal element
1234 if (bUnUsed[j] && (TMath::Abs(matV[j][j])>TMath::Max(TMath::Abs(vPivot),eps*diagV[j]))) {
1235 vPivot = matV[j][j];
1240 if (iPivot >= 0) { // pivot found
1242 bUnUsed[iPivot] = false;
1243 vPivot = 1.0/vPivot;
1244 matV[iPivot][iPivot] = -vPivot; // Replace pivot by its inverse
1246 for (Int_t j=0; j<nLoc; j++) {
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];
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];
1264 else { // No more pivot value (clear those elements)
1265 for (Int_t j=0; j<nLoc; j++) {
1269 for (Int_t k=0; k<j; k++) {
1276 break; // No more pivots anyway, stop here
1280 for (Int_t j=0; j<nLoc; j++) {
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];
1288 for (Int_t j=0; j<nLoc; j++) {
1301 -----------------------------------------------------------
1303 -----------------------------------------------------------
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.
1310 CALL SPAVAT(V,A,W,N,M) W = A * V * A
1313 where V = symmetric N-by-N matrix
1314 A = general N-by-M matrix
1315 W = symmetric M-by-M matrix
1316 -----------------------------------------------------------
1318 Int_t AliMillepede::SpAVAt(double matV[][fgkMaxLocalPar], double matA[][fgkMaxLocalPar], double matW[][fgkMaxGlobalPar], int nLoc, int nGlo)
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.
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
1335 matW[j][i] = matW[i][j]; // Matrix w is symmetric
1345 -----------------------------------------------------------
1347 -----------------------------------------------------------
1349 multiply general M-by-N matrix A and N-vector X
1351 CALL SPAX(A,X,Y,M,N) Y = A * X
1354 where A = general M-by-N matrix (A11 A12 ... A1N A21 A22 ...)
1357 -----------------------------------------------------------
1359 Int_t AliMillepede::SpAX(double matA[][fgkMaxLocalPar], double vecX[], double vecY[], int nCol, int nRow)
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
1374 -----------------------------------------------------------
1376 -----------------------------------------------------------
1378 Print the final results into the logfile
1380 -----------------------------------------------------------
1382 Int_t AliMillepede::PrintGlobalParameters() const
1384 /// Print the final results into the logfile
1386 double lGlobalCor =0.;
1389 AliInfo(" Result of fit for global parameters");
1390 AliInfo(" ===================================");
1391 AliInfo(" I initial final differ lastcor error gcor");
1392 AliInfo("-----------------------------------------------------------------------------------");
1394 for (int i=0; i<fNGlobalPar; i++) {
1395 lError = TMath::Sqrt(TMath::Abs(fMatCGlo[i][i]));
1396 if (fMatCGlo[i][i] < 0.0) lError = -lError;
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])));
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));
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]));
1413 ----------------------------------------------------------------
1414 CHI2DOFLIM: return the limit in chi^2/nd for n sigmas stdev authorized
1415 ----------------------------------------------------------------
1417 Only n=1, 2, and 3 are expected in input
1418 ----------------------------------------------------------------
1420 double AliMillepede::Chi2DoFLim(int nSig, int nDoF)
1422 /// return the limit in chi^2/nd for n sigmas stdev authorized
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}};
1442 lNSig = TMath::Max(1,TMath::Min(nSig,3));
1445 return table[lNSig-1][nDoF-1];
1447 else { // approximation
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);