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 //-----------------------------------------------------------------------------
37 #include "AliMillepede.h"
39 //=============================================================================
40 AliMillepede::AliMillepede()
42 fIndexLocEq(fgkMaxGlobalPar+fgkMaxLocalPar),
43 fDerivLocEq(fgkMaxGlobalPar+fgkMaxLocalPar),
44 fIndexAllEqs(1000*fgkMaxGlobalPar+fgkMaxLocalPar),
45 fDerivAllEqs(1000*fgkMaxGlobalPar+fgkMaxLocalPar),
60 fNGlobalConstraints(0),
62 fNLocalFitsRejected(0),
66 /// Standard constructor
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");
80 //=============================================================================
81 AliMillepede::~AliMillepede() {
86 //=============================================================================
87 Int_t AliMillepede::InitMille(int nGlo, int nLoc, int lNStdDev,
88 double lResCut, double lResCutInit)
90 /// Initialization of millepede
92 AliDebug(1,"----------------------------------------------------");
94 AliDebug(1," Entering InitMille");
96 AliDebug(1,"-----------------------------------------------------");
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
104 AliMillepede::SetNLocalEquations(0); // Number of local fits (starts at 0)
106 fIter = 0; // By default iterations are turned off, turned on by use SetIterations
109 fResCutInit = lResCutInit;
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
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));
120 if (fNGlobalPar>fgkMaxGlobalPar || fNLocalPar>fgkMaxLocalPar) {
121 AliDebug(1,"Two many parameters !!!!!");
125 // Global parameters initializations
126 for (int i=0; i<fNGlobalPar; i++) {
135 for (int j=0; j<fNGlobalPar;j++) {
140 // Local parameters initializations
142 for (int i=0; i<fNLocalPar; i++) {
145 for (int j=0; j<fNLocalPar;j++) {
150 // Then we fix all parameters...
152 for (int j=0; j<fNGlobalPar; j++) {
153 AliMillepede::SetParSigma(j,0.0);
156 fDerivLocEq.Reset(); fNDerivLocEq=0;
157 fIndexLocEq.Reset(); fNIndexLocEq=0;
159 fIndexAllEqs.Reset(); fNIndexAllEqs=0;
160 fDerivAllEqs.Reset(); fNDerivAllEqs=0;
161 fLocEqPlace.Reset(); fNLocEqPlace=0;
164 AliDebug(1,"----------------------------------------------------");
166 AliDebug(1," InitMille has been successfully called!");
168 AliDebug(1,"-----------------------------------------------------");
175 -----------------------------------------------------------
176 PARGLO: initialization of global parameters
177 -----------------------------------------------------------
179 param = array of starting values
181 -----------------------------------------------------------
183 Int_t AliMillepede::SetGlobalParameters(double *param)
185 /// initialization of global parameters
186 for(Int_t iPar=0; iPar<fNGlobalPar; iPar++){
187 fInitPar[iPar] = param[iPar];
194 -----------------------------------------------------------
195 PARGLO: initialization of global parameters
196 -----------------------------------------------------------
198 iPar = the index of the global parameter in the
199 result array (equivalent to fDeltaPar[]).
201 param = the starting value
203 -----------------------------------------------------------
205 Int_t AliMillepede::SetGlobalParameter(int iPar, double param)
207 /// initialization of global parameter iPar
208 if (iPar<0 || iPar>=fNGlobalPar) {
212 fInitPar[iPar] = param;
219 -----------------------------------------------------------
220 PARSIG: define a constraint for a single global param
221 param is 'encouraged' to vary within [-sigma;sigma]
223 -----------------------------------------------------------
225 iPar = the index of the global parameter in the
226 result array (equivalent to fDeltaPar[]).
228 sigma = value of the constraint (sigma <= 0. will
229 mean that parameter is FIXED !!!)
231 -----------------------------------------------------------
233 Int_t AliMillepede::SetParSigma(int iPar, double sigma)
235 /// Define a range [-sigma;sigma] where iPar is encourage to vary
236 if (iPar>=fNGlobalPar) {
240 fSigmaPar[iPar] = sigma;
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 -----------------------------------------------------------
254 iPar = the index of the global parameter in the
255 result array (equivalent to fDeltaPar[]).
257 -----------------------------------------------------------
259 Int_t AliMillepede::SetNonLinear(int iPar)
261 /// Set nonlinear flag for iPar
262 if (iPar<0 || iPar>=fNGlobalPar) {
266 fIsNonLinear[iPar] = 1;
274 -----------------------------------------------------------
275 INITUN: unit for iteration
276 -----------------------------------------------------------
278 lChi2CutFac is used by Fitloc to define the Chi^2/ndof cut value
280 A large cutfac value enables to take a wider range of tracks
281 for first iterations, which might be useful if misalignments
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.
288 At least one more iteration is often needed in order to remove
289 tracks containing outliers.
291 -----------------------------------------------------------
293 Int_t AliMillepede::SetIterations(double lChi2CutFac)
295 /// Number of iterations is calculated from lChi2CutFac
296 fChi2CutFactor = TMath::Max(1.0, lChi2CutFac);
298 AliInfo(Form("Initial cut factor is %f",fChi2CutFactor));
299 fIter = 1; // Initializes the iteration process
304 -----------------------------------------------------------
305 CONSTF: define a constraint equation in AliMillepede
306 -----------------------------------------------------------
308 dercs = the row containing constraint equation
309 derivatives (put into the final matrix)
311 lLagMult = the lagrange multiplier value (sum of equation)
313 -----------------------------------------------------------
315 Int_t AliMillepede::SetGlobalConstraint(double dercs[], double lLagMult)
317 /// Define a constraint equation
318 if (fNGlobalConstraints>=fgkMaxGloCsts) {
319 AliInfo("Too many constraints !!!");
323 for (int i=0; i<fNGlobalPar; i++) {
324 fMatDerConstr[fNGlobalConstraints][i] = dercs[i];
327 fLagMult[fNGlobalConstraints] = lLagMult;
328 fNGlobalConstraints++ ;
329 AliInfo(Form("Number of constraints increased to %d",fNGlobalConstraints));
335 -----------------------------------------------------------
336 EQULOC: write ONE equation in the matrices
337 -----------------------------------------------------------
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!!!)
344 -----------------------------------------------------------
346 Int_t AliMillepede::SetLocalEquation(double dergb[], double derlc[], double lMeas, double lSigma)
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++) {
353 for (int i=0; i<fNGlobalPar; i++) {
359 // Serious equation, initialize parameters
360 double lWeight = 1.0/(lSigma*lSigma);
366 for (int i=0; i<fNLocalPar; i++) { // Retrieve local param interesting indices
368 if (iLocFirst == -1) {
369 iLocFirst = i; // first index
371 iLocLast = i; // last index
374 AliDebug(2,Form("%d / %d",iLocFirst, iLocLast));
376 for (int i=0; i<fNGlobalPar; i++) { // Idem for global parameters
378 if (iGlobFirst == -1) {
379 iGlobFirst = i; // first index
381 iGlobLast = i; // last index
384 AliDebug(2,Form("%d / %d",iGlobFirst,iGlobLast));
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++);
392 for (int i=iLocFirst; i<=iLocLast; i++) { // Store interesting local parameters
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++);
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++);
407 for (int i=iGlobFirst; i<=iGlobLast; i++) { // Store interesting global parameters
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++);
417 AliDebug(2,Form("Out Equloc -- NST = %d",fNDerivLocEq));
422 -----------------------------------------------------------
423 FITLOC: perform local params fit, once all the equations
424 have been written by EquLoc
425 -----------------------------------------------------------
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)
431 localParams = contains the fitted track parameters and
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
439 -----------------------------------------------------------
441 Int_t AliMillepede::LocalFit(int iFit, double localParams[], Bool_t bSingleFit)
443 /// Perform local parameters fit once all the local equations have been set
445 // Few initializations
448 int iIdx, jIdx, kIdx;
449 int iIdxIdx, jIdxIdx;
459 double lWeight = 0.0;
462 double lRedChi2 = 0.0;
463 double lChi2Cut = 0.0;
466 int nEqTerms = fNDerivLocEq; // Number of terms (local + global derivatives,
467 // measurement and weight) involved in this local equation
470 // Fill the track store at first pass
472 if (fIter < 2 && !bSingleFit) { // Do it only once
473 AliDebug(1,Form("Store equation no: %d", iFit));
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++);
481 if (fNLocEqPlace==fLocEqPlace.GetSize()) fLocEqPlace.Set(2*fNLocEqPlace);
482 fLocEqPlace.AddAt(fNIndexAllEqs,fNLocEqPlace++);
485 AliDebug(2,Form("FLocEqPlace size = %d",fLocEqPlace[iFit]));
486 AliDebug(2,Form("FIndexAllEqs size = %d",fNIndexAllEqs));
490 for (i=0; i<fNLocalPar; i++) { // reset local params
493 for (j=0; j<fNLocalPar; j++) {
494 fMatCLoc[i][j] = 0.0;
498 for (i=0; i<fNGlobalPar; i++) {
499 fGlo2CGLRow[i] = -1; // reset mixed params
505 LOOPS : HOW DOES IT WORKS ?
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:
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.
517 We took them and store them into matrices.
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)
528 // FIRST LOOP : local track fit
532 // fIndexLocEq.push_back(-1);
533 if (fNIndexLocEq==fIndexLocEq.GetSize()) fIndexLocEq.Set(2*fNIndexLocEq);
534 fIndexLocEq.AddAt(-1,fNIndexLocEq++);
536 while (iEqTerm <= nEqTerms) {
537 if (fIndexLocEq[iEqTerm] == -1) {
538 if (iMeas == -1) { // First -1 : lMeas
540 iLocFirst = iEqTerm+1;
542 else if (iWeight == 0) { // Second -1 : weight
544 iLocLast = iEqTerm-1;
545 iGloFirst = iEqTerm+1;
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));
554 // Now suppress the global part (only relevant with iterations)
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
563 lMeas -= fDerivLocEq[i]*(fDeltaPar[iIdx]); // nonlinear parameter
565 // AliDebug(2,Form("lMeas after global stuff removal = %f", lMeas));
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]));
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]));
580 iEqTerm--; // end of one equation is the beginning of next
581 } // End of "end of equation" operations
582 } // End of loop on equation
584 } // End of loop on all equations used in the fit
588 // Local params matrix is completed, now invert to solve...
591 Int_t nRank = AliMillepede::SpmInv(fMatCLoc, fVecBLoc, fNLocalPar);
592 // nRank is the number of nonzero diagonal elements
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)");
599 for (i=0; i<fNLocalPar; i++) {
600 AliDebug(1,Form("%d / %.6f / %.6f", i, fVecBLoc[i], TMath::Sqrt(fMatCLoc[i][i])));
604 // Store the track params and errors
606 for (i=0; i<fNLocalPar; i++) {
607 localParams[2*i] = fVecBLoc[i];
608 localParams[2*i+1] = TMath::Sqrt(TMath::Abs(fMatCLoc[i][i]));
613 // SECOND LOOP : residual calculation
620 while (iEqTerm <= nEqTerms) {
621 if (fIndexLocEq[iEqTerm] == -1) {
622 if (iMeas == -1) { // First -1 : lMeas
624 iLocFirst = iEqTerm+1;
626 else if (iWeight == 0) { // Second -1 : weight
628 iLocLast = iEqTerm-1;
629 iGloFirst = iEqTerm+1;
631 else { // Third -1 : end of equation; start of next
632 iGloLast = iEqTerm-1;
633 lMeas = fDerivLocEq[iMeas];
634 lWeight = fDerivLocEq[iWeight];
636 // Print all (for debugging purposes)
638 // int nDerLoc = iLocLast-iLocFirst+1; // Number of local derivatives involved
639 // int nDerGlo = iGloLast-iGloFirst+1; // Number of global derivatives involved
642 // AliDebug(2,Form(". equation: measured value %.6f +/- %.6f", lMeas, 1.0/TMath::Sqrt(lWeight)));
643 // AliDebug(2,Form("Number of derivatives (global, local): %d, %d",nDerGlo,nDerLoc));
644 // AliDebug(2,"Global derivatives are: (index/derivative/parvalue) ");
646 // for (i=iGloFirst; i<=iGloLast; i++) {
647 // AliDebug(2,Form("%d / %.6f / %.6f",fIndexLocEq[i],fDerivLocEq[i],fInitPar[fIndexLocEq[i]]));
650 // AliDebug(2,"Local derivatives are: (index/derivative) ");
652 // for (i=(ja+1); i<jb; i++) {AliDebug(2,Form("%d / %.6f",fIndexLocEq[i], fDerivLocEq[i]));}
654 // Now suppress local and global parts to LMEAS;
656 // First the local part
657 for (i=iLocFirst; i<=iLocLast; i++) {
658 iIdx = fIndexLocEq[i];
659 lMeas -= fDerivLocEq[i]*fVecBLoc[iIdx];
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
667 lMeas -= fDerivLocEq[i]*(fDeltaPar[iIdx]); // nonlinear parameter
670 // lMeas contains now the residual value
671 AliDebug(2,Form("Residual value : %.6f", lMeas));
673 // reject the track if lMeas is too important (outlier)
674 if (TMath::Abs(lMeas) >= fResCutInit && fIter <= 1) {
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;
682 if (TMath::Abs(lMeas) >= fResCut && fIter > 1) {
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;
690 lChi2 += lWeight*lMeas*lMeas ; // total chi^2
691 nEq++; // number of equations
695 } // End of "end of equation" operations
696 } // End of loop on equation
698 } // End of loop on all equations used in the fit
703 AliDebug(1,Form("Final chi square / degrees of freedom %.2f / %d",lChi2, nDoF));
705 if (nDoF > 0) lRedChi2 = lChi2/float(nDoF); // Chi^2/dof
709 if (fNStdDev != 0 && nDoF > 0 && !bSingleFit) // Chisquare cut
711 lChi2Cut = AliMillepede::Chi2DoFLim(fNStdDev, nDoF)*fChi2CutFactor;
713 AliDebug(1,Form("Reject if Chisq/Ndf = %.4f > %.4f",lRedChi2,lChi2Cut));
715 if (lRedChi2 > lChi2Cut) // Reject the track if too much...
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;
725 if (bSingleFit) // Stop here if just updating the track parameters
727 fIndexLocEq.Reset(); fNIndexLocEq=0; // Reset store for the next track
728 fDerivLocEq.Reset(); fNDerivLocEq=0;
733 // THIRD LOOP: local operations are finished, track is accepted
734 // We now update the global parameters (other matrices)
741 while (iEqTerm <= nEqTerms)
743 if (fIndexLocEq[iEqTerm] == -1)
745 if (iMeas == -1) { // First -1 : lMeas
747 iLocFirst = iEqTerm+1;
749 else if (iWeight == 0) { // Second -1 : weight
751 iLocLast = iEqTerm-1;
752 iGloFirst = iEqTerm+1;
754 else { // Third -1 : end of equation; start of next
755 iGloLast = iEqTerm-1;
756 lMeas = fDerivLocEq[iMeas];
757 lWeight = fDerivLocEq[iWeight];
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
765 lMeas -= fDerivLocEq[i]*(fDeltaPar[iIdx]); // nonlinear parameter
768 for (i=iGloFirst; i<=iGloLast; i++) {
769 iIdx = fIndexLocEq[i]; // Global param indice (the matrix line)
771 fVecBGlo[iIdx] += lWeight*lMeas*fDerivLocEq[i];
772 // AliDebug(2,Form("fVecBGlo[%d] = %.6f", j, fVecBGlo[j] ));
774 // First of all, the global/global terms (exactly like local matrix)
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]));
782 // Now we have also rectangular matrices containing global/local terms.
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
789 fGlo2CGLRow[iIdx] = nGloInFit;
790 fCGLRow2Glo[nGloInFit] = iIdx;
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]));
805 } // End of "end of equation" operations
806 } // End of loop on equation
808 } // End of loop on all equations used in the fit
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);
814 for (iIdxIdx=0; iIdxIdx<nGloInFit; iIdxIdx++) {
815 iIdx = fCGLRow2Glo[iIdxIdx];
816 fVecBGlo[iIdx] -= fVecBGloCorr[iIdxIdx];
818 for (jIdxIdx=0; jIdxIdx<=iIdxIdx; jIdxIdx++) {
819 jIdx = fCGLRow2Glo[jIdxIdx];
820 fMatCGlo[iIdx][jIdx] -= fMatCGloCorr[iIdxIdx][jIdxIdx];
821 fMatCGlo[jIdx][iIdx] = fMatCGlo[iIdx][jIdx];
825 fIndexLocEq.Reset(); fNIndexLocEq=0; // Reset store for the next track
826 fDerivLocEq.Reset(); fNDerivLocEq=0;
833 -----------------------------------------------------------
834 MAKEGLOBALFIT: perform global params fit, once all the 'tracks'
835 have been fitted by FitLoc
836 -----------------------------------------------------------
838 par[] = array containing the computed global
839 parameters (the misalignment constants)
841 error[] = array containing the error on global
842 parameters (estimated by AliMillepede)
844 pull[] = array containing the corresponding pulls
846 -----------------------------------------------------------
848 Int_t AliMillepede::GlobalFit(double par[], double error[], double pull[])
850 /// perform global parameters fit once all the local equations have been fitted
856 double step[fgkMaxGlobalPar];
858 double localPars[2*fgkMaxLocalPar];
860 int nLocFitsGood = 0;
864 AliInfo("..... Making global fit .....");
865 AliInfo(Form("First Iteration with cut factor %.2f", fChi2CutFactor));
867 nLocFitsTot = AliMillepede::GetNLocalEquations();
869 while (fIter < fMaxIter) // Iteration for the final loop
872 nLocFits = AliMillepede::GetNLocalEquations();
873 AliInfo(Form("...using %d local fits...",nLocFits));
875 // Start by saving the diagonal elements
877 for (i=0; i<fNGlobalPar; i++) {
878 fDiagCGlo[i] = fMatCGlo[i][i];
881 // Then we retrieve the different constraints: fixed parameter or global equation
883 nGloFix = 0; // First look at the fixed global params
885 for (i=0; i<fNGlobalPar; i++) {
886 if (fSigmaPar[i] <= 0.0) { // fixed global param
888 for (j=0; j<fNGlobalPar; j++) {
889 fMatCGlo[i][j] = 0.0; // Reset row and column
890 fMatCGlo[j][i] = 0.0;
894 fMatCGlo[i][i] += 1.0/(fSigmaPar[i]*fSigmaPar[i]);
898 nVar = fNGlobalPar; // Current number of equations
899 AliDebug(1,Form("Number of constraint equations : %d", fNGlobalConstraints));
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]);
909 fMatCGlo[nVar][nVar] = 0.0;
910 fVecBGlo[nVar] = float(nLocFits)*lConstraint;
915 // Intended to compute the final global chisquare
917 double lFinalCor = 0.0;
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]);
931 AliInfo(Form(" Final coeff is %.6f",lFinalCor));
932 AliInfo(Form(" Final NDOFs = %d", fNGlobalPar));
934 // The final matrix inversion
936 Int_t nRank = AliMillepede::SpmInv(fMatCGlo, fVecBGlo, nVar);
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]));
942 AliDebug(1,Form("err = %.6f", TMath::Sqrt(TMath::Abs(fMatCGlo[i][i]))));
944 step[i] = fVecBGlo[i];
946 if (fIter == 1) error[i] = fMatCGlo[i][i]; // Unfitted error
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));
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
970 AliInfo("----------------------------------------------");
971 AliInfo(Form("Iteration %d with cut factor %.2f", fIter, fChi2CutFactor));
973 // Reset global variables
975 for (i=0; i<nVar; i++) {
977 for (j=0; j<nVar; j++) {
978 fMatCGlo[i][j] = 0.0;
983 // We start a new iteration
986 // First we read the stores for retrieving the local params
990 for (i=0; i<nLocFitsTot; i++) {
994 (i>0) ? iEqFirst = fLocEqPlace[i-1] : iEqFirst = 0;
995 iEqLast = fLocEqPlace[i];
997 AliDebug(2,Form("Track %d : ",i));
998 AliDebug(2,Form("Starts at %d", iEqFirst));
999 AliDebug(2,Form("Ends at %d",iEqLast));
1001 if (fIndexAllEqs[iEqFirst] != -999) { // Fit is still OK
1002 fIndexLocEq.Reset(); fNIndexLocEq=0;
1003 fDerivLocEq.Reset(); fNDerivLocEq=0;
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++);
1012 for (j=0; j<2*fNLocalPar; j++) {
1016 // Int_t sc = AliMillepede::LocalFit(i,localPars,0);
1017 // (sc) ? nLocFitsGood++ : fIndexAllEqs[iEqFirst] = -999;
1018 AliMillepede::LocalFit(i,localPars,0);
1021 } // End of loop on fits
1023 AliMillepede::SetNLocalEquations(nLocFitsGood);
1025 } // End of iteration loop
1027 AliMillepede::PrintGlobalParameters(); // Print the final results
1029 for (j=0; j<fNGlobalPar; j++) {
1030 par[j] = fInitPar[j]+fDeltaPar[j];
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]);
1032 error[j] = TMath::Sqrt(TMath::Abs(fMatCGlo[j][j]));
1036 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.");
1049 -----------------------------------------------------------
1050 ERRPAR: return error for parameter iPar
1051 -----------------------------------------------------------
1053 iPar = the index of the global parameter in the
1054 result array (equivalent to fDeltaPar[]).
1056 -----------------------------------------------------------
1058 Double_t AliMillepede::GetParError(int iPar) const
1060 /// return error for parameter iPar
1061 Double_t lErr = -1.;
1062 if (iPar>=0 && iPar<fNGlobalPar) {
1063 lErr = TMath::Sqrt(TMath::Abs(fMatCGlo[iPar][iPar]));
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 -----------------------------------------------------------
1075 Solve the equation : V * X = B
1077 V is replaced by inverse matrix and B by X, the solution vector
1078 -----------------------------------------------------------
1080 int AliMillepede::SpmInv(double matV[][fgkMaxGloPC], double vecB[], int nGlo)
1082 /// Obtain solution of a system of linear equations with symmetric matrix
1083 /// and the inverse (using 'singular-value friendly' GAUSS pivot)
1088 double eps = 0.00000000000001;
1090 bool *bUnUsed = new bool[nGlo];
1091 double *diagV = new double[nGlo];
1092 double *rowMax = new double[nGlo];
1093 double *colMax = new double[nGlo];
1095 double *temp = new double[nGlo];
1097 for (Int_t i=0; i<nGlo; i++) {
1102 for (Int_t j=0; j<i; j++) {
1103 if (matV[j][i] == 0) {
1104 matV[j][i] = matV[i][j];
1109 // Small loop for matrix equilibration (gives a better conditioning)
1111 for (Int_t i=0; i<nGlo; i++) {
1112 for (Int_t j=0; j<nGlo; j++) {
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
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
1123 for (Int_t i=0; i<nGlo; i++) {
1124 for (Int_t j=0; j<nGlo; j++) {
1125 matV[i][j] = TMath::Sqrt(rowMax[i])*matV[i][j]*TMath::Sqrt(colMax[j]); // Equilibrate the V matrix
1127 diagV[i] = TMath::Abs(matV[i][i]); // save diagonal elem absolute values
1131 for (Int_t i=0; i<nGlo; i++) {
1135 for (Int_t j=0; j<nGlo; j++) { // First look for the pivot, ie max unused diagonal element
1136 if (bUnUsed[j] && (TMath::Abs(matV[j][j])>TMath::Max(TMath::Abs(vPivot),eps*diagV[j]))) {
1137 vPivot = matV[j][j];
1142 if (iPivot >= 0) { // pivot found
1144 bUnUsed[iPivot] = false; // This value is used
1145 vPivot = 1.0/vPivot;
1146 matV[iPivot][iPivot] = -vPivot; // Replace pivot by its inverse
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];
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
1163 else { // No more pivot value (clear those elements)
1164 for (Int_t j=0; j<nGlo; j++) {
1168 for (Int_t k=0; k<nGlo; k++) {
1174 break; // No more pivots anyway, stop here
1178 for (Int_t i=0; i<nGlo; i++) {
1179 for (Int_t j=0; j<nGlo; j++) {
1180 matV[i][j] = TMath::Sqrt(colMax[i])*matV[i][j]*TMath::Sqrt(rowMax[j]); // Correct matrix V
1184 for (Int_t j=0; j<nGlo; j++) {
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];
1193 for (Int_t j=0; j<nGlo; j++) {
1194 vecB[j] = temp[j]; // The final result
1207 // Same method but for local fit, so heavily simplified
1209 int AliMillepede::SpmInv(double matV[][fgkMaxLocalPar], double vecB[], int nLoc)
1211 /// Obtain solution of a system of linear equations with symmetric matrix
1212 /// and the inverse (using 'singular-value friendly' GAUSS pivot)
1217 double eps = 0.0000000000001;
1219 bool *bUnUsed = new bool[nLoc];
1220 double *diagV = new double[nLoc];
1221 double *temp = new double[nLoc];
1223 for (Int_t i=0; i<nLoc; i++) {
1225 diagV[i] = TMath::Abs(matV[i][i]); // save diagonal elem absolute values
1226 for (Int_t j=0; j<i; j++) {
1227 matV[j][i] = matV[i][j] ;
1231 for (Int_t i=0; i<nLoc; i++) {
1235 for (Int_t j=0; j<nLoc; j++) { // First look for the pivot, ie max unused diagonal element
1236 if (bUnUsed[j] && (TMath::Abs(matV[j][j])>TMath::Max(TMath::Abs(vPivot),eps*diagV[j]))) {
1237 vPivot = matV[j][j];
1242 if (iPivot >= 0) { // pivot found
1244 bUnUsed[iPivot] = false;
1245 vPivot = 1.0/vPivot;
1246 matV[iPivot][iPivot] = -vPivot; // Replace pivot by its inverse
1248 for (Int_t j=0; j<nLoc; j++) {
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];
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];
1266 else { // No more pivot value (clear those elements)
1267 for (Int_t j=0; j<nLoc; j++) {
1271 for (Int_t k=0; k<j; k++) {
1278 break; // No more pivots anyway, stop here
1282 for (Int_t j=0; j<nLoc; j++) {
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];
1290 for (Int_t j=0; j<nLoc; j++) {
1303 -----------------------------------------------------------
1305 -----------------------------------------------------------
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.
1312 CALL SPAVAT(V,A,W,N,M) W = A * V * A
1315 where V = symmetric N-by-N matrix
1316 A = general N-by-M matrix
1317 W = symmetric M-by-M matrix
1318 -----------------------------------------------------------
1320 Int_t AliMillepede::SpAVAt(double matV[][fgkMaxLocalPar], double matA[][fgkMaxLocalPar], double matW[][fgkMaxGlobalPar], int nLoc, int nGlo)
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.
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
1337 matW[j][i] = matW[i][j]; // Matrix w is symmetric
1347 -----------------------------------------------------------
1349 -----------------------------------------------------------
1351 multiply general M-by-N matrix A and N-vector X
1353 CALL SPAX(A,X,Y,M,N) Y = A * X
1356 where A = general M-by-N matrix (A11 A12 ... A1N A21 A22 ...)
1359 -----------------------------------------------------------
1361 Int_t AliMillepede::SpAX(double matA[][fgkMaxLocalPar], double vecX[], double vecY[], int nCol, int nRow)
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
1376 -----------------------------------------------------------
1378 -----------------------------------------------------------
1380 Print the final results into the logfile
1382 -----------------------------------------------------------
1384 Int_t AliMillepede::PrintGlobalParameters() const
1386 /// Print the final results into the logfile
1388 double lGlobalCor =0.;
1391 AliInfo(" Result of fit for global parameters");
1392 AliInfo(" ===================================");
1393 AliInfo(" I initial final differ lastcor error gcor");
1394 AliInfo("-----------------------------------------------------------------------------------");
1396 for (int i=0; i<fNGlobalPar; i++) {
1397 lError = TMath::Sqrt(TMath::Abs(fMatCGlo[i][i]));
1398 if (fMatCGlo[i][i] < 0.0) lError = -lError;
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])));
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));
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]));
1415 ----------------------------------------------------------------
1416 CHI2DOFLIM: return the limit in chi^2/nd for n sigmas stdev authorized
1417 ----------------------------------------------------------------
1419 Only n=1, 2, and 3 are expected in input
1420 ----------------------------------------------------------------
1422 double AliMillepede::Chi2DoFLim(int nSig, int nDoF)
1424 /// return the limit in chi^2/nd for n sigmas stdev authorized
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}};
1444 lNSig = TMath::Max(1,TMath::Min(nSig,3));
1447 return table[lNSig-1][nDoF-1];
1449 else { // approximation
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);