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 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 int 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
858 double localPars[2*fgkMaxLocalPar];
860 int nLocFitsGood = 0;
864 AliInfo("..... Making global fit .....");
866 nLocFitsTot = AliMillepede::GetNLocalEquations();
868 while (fIter < fMaxIter) // Iteration for the final loop
871 nLocFits = AliMillepede::GetNLocalEquations();
872 AliInfo(Form("...using %d local fits...",nLocFits));
874 // Start by saving the diagonal elements
876 for (i=0; i<fNGlobalPar; i++) {
877 fDiagCGlo[i] = fMatCGlo[i][i];
880 // Then we retrieve the different constraints: fixed parameter or global equation
882 nGloFix = 0; // First look at the fixed global params
884 for (i=0; i<fNGlobalPar; i++) {
885 if (fSigmaPar[i] <= 0.0) { // fixed global param
887 for (j=0; j<fNGlobalPar; j++) {
888 fMatCGlo[i][j] = 0.0; // Reset row and column
889 fMatCGlo[j][i] = 0.0;
893 fMatCGlo[i][i] += 1.0/(fSigmaPar[i]*fSigmaPar[i]);
897 nVar = fNGlobalPar; // Current number of equations
898 AliDebug(1,Form("Number of constraint equations : %d", fNGlobalConstraints));
900 for (i=0; i<fNGlobalConstraints; i++) { // Then the constraint equation
901 lConstraint = fLagMult[i];
902 for (j=0; j<fNGlobalPar; j++) {
903 fMatCGlo[nVar][j] = float(nLocFits)*fMatDerConstr[i][j];
904 fMatCGlo[j][nVar] = float(nLocFits)*fMatDerConstr[i][j];
905 lConstraint -= fMatDerConstr[i][j]*(fInitPar[j]+fDeltaPar[j]);
908 fMatCGlo[nVar][nVar] = 0.0;
909 fVecBGlo[nVar] = float(nLocFits)*lConstraint;
914 // Intended to compute the final global chisquare
916 double lFinalCor = 0.0;
919 for (i=0; i<fNGlobalPar; i++) {
920 for (j=0; j<fNGlobalPar; j++) {
921 // printf("%d, %d, %.6f %.6f %.6f\n",i,j,step[i],fMatCGlo[i][j],step[j]);
922 lFinalCor += step[i]*fMatCGlo[i][j]*step[j];
923 if (i == j && fSigmaPar[i] != 0) {
924 lFinalCor -= step[i]*step[i]/(fSigmaPar[i]*fSigmaPar[i]);
930 AliInfo(Form(" Final coeff is %.6f",lFinalCor));
931 AliInfo(Form(" Final NDOFs = %d", fNGlobalPar));
933 // The final matrix inversion
935 Int_t nRank = AliMillepede::SpmInv(fMatCGlo, fVecBGlo, nVar);
937 for (i=0; i<fNGlobalPar; i++) {
938 fDeltaPar[i] += fVecBGlo[i]; // Update global parameters values (for iterations)
939 AliDebug(1,Form("fDeltaPar[%d] = %.6f", i, fDeltaPar[i]));
940 AliDebug(1,Form("fMatCGlo[%d][%d] = %.6f", i, i, fMatCGlo[i][i]));
941 AliDebug(1,Form("err = %.6f", TMath::Sqrt(TMath::Abs(fMatCGlo[i][i]))));
943 step[i] = fVecBGlo[i];
945 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));
952 AliInfo(Form("Total : %d local fits, %d rejected.", fNLocalFits, fNLocalFitsRejected));
953 if (fIter == 0) break; // No iterations set
956 if (fIter == fMaxIter) break; // End of story
958 // Reinitialize parameters for iteration
961 fNLocalFitsRejected = 0;
963 if (fChi2CutFactor != fChi2CutRef) {
964 fChi2CutFactor = TMath::Sqrt(fChi2CutFactor);
965 if (fChi2CutFactor < 1.2*fChi2CutRef) {
966 fChi2CutFactor = fChi2CutRef;
967 fIter = fMaxIter - 1; // Last iteration
970 AliInfo(Form("Iteration %d with cut factor %.2f", fIter, fChi2CutFactor));
972 // Reset global variables
974 for (i=0; i<nVar; i++) {
976 for (j=0; j<nVar; j++) {
977 fMatCGlo[i][j] = 0.0;
982 // We start a new iteration
985 // First we read the stores for retrieving the local params
989 for (i=0; i<nLocFitsTot; i++) {
993 (i>0) ? iEqFirst = fLocEqPlace[i-1] : iEqFirst = 0;
994 iEqLast = fLocEqPlace[i];
996 AliDebug(2,Form("Track %d : ",i));
997 AliDebug(2,Form("Starts at %d", iEqFirst));
998 AliDebug(2,Form("Ends at %d",iEqLast));
1000 if (fIndexAllEqs[iEqFirst] != -999) { // Fit is still OK
1001 fIndexLocEq.Reset(); fNIndexLocEq=0;
1002 fDerivLocEq.Reset(); fNDerivLocEq=0;
1004 for (j=iEqFirst; j<iEqLast; j++) {
1005 if (fNIndexLocEq==fIndexLocEq.GetSize()) fIndexLocEq.Set(2*fNIndexLocEq);
1006 fIndexLocEq.AddAt(fIndexAllEqs[j],fNIndexLocEq++);
1007 if (fNDerivLocEq==fDerivLocEq.GetSize()) fDerivLocEq.Set(2*fNDerivLocEq);
1008 fDerivLocEq.AddAt(fDerivAllEqs[j],fNDerivLocEq++);
1011 for (j=0; j<2*fNLocalPar; j++) {
1015 // Int_t sc = AliMillepede::LocalFit(i,localPars,0);
1016 // (sc) ? nLocFitsGood++ : fIndexAllEqs[iEqFirst] = -999;
1017 AliMillepede::LocalFit(i,localPars,0);
1020 } // End of loop on fits
1022 AliMillepede::SetNLocalEquations(nLocFitsGood);
1024 } // End of iteration loop
1026 AliMillepede::PrintGlobalParameters(); // Print the final results
1028 for (j=0; j<fNGlobalPar; j++) {
1029 par[j] = fInitPar[j]+fDeltaPar[j];
1030 pull[j] = (fSigmaPar[j] <= 0.) ? 0. : fDeltaPar[j]/TMath::Sqrt(fSigmaPar[j]*fSigmaPar[j]-fMatCGlo[j][j]);
1031 error[j] = TMath::Sqrt(TMath::Abs(fMatCGlo[j][j]));
1035 AliInfo(" * o o o ");
1037 AliInfo(" o ooooo o o o oo ooo oo ooo oo ");
1038 AliInfo(" o o o o o o o o o o o o o o o o ");
1039 AliInfo(" o o o o o o oooo o o oooo o o oooo ");
1040 AliInfo(" o o o o o o o ooo o o o o ");
1041 AliInfo(" o o o o o o oo o oo ooo oo ++ ends.");
1048 -----------------------------------------------------------
1049 ERRPAR: return error for parameter iPar
1050 -----------------------------------------------------------
1052 iPar = the index of the global parameter in the
1053 result array (equivalent to fDeltaPar[]).
1055 -----------------------------------------------------------
1057 Double_t AliMillepede::GetParError(Int_t iPar) const
1059 /// return error for parameter iPar
1060 Double_t lErr = -1.;
1061 if (iPar>=0 && iPar<fNGlobalPar) {
1062 lErr = TMath::Sqrt(TMath::Abs(fMatCGlo[iPar][iPar]));
1069 -----------------------------------------------------------
1070 SPMINV: obtain solution of a system of linear equations with symmetric matrix
1071 and the inverse (using 'singular-value friendly' GAUSS pivot)
1072 -----------------------------------------------------------
1074 Solve the equation : V * X = B
1076 V is replaced by inverse matrix and B by X, the solution vector
1077 -----------------------------------------------------------
1079 int AliMillepede::SpmInv(double matV[][fgkMaxGloPC], double vecB[], int nGlo)
1081 /// Obtain solution of a system of linear equations with symmetric matrix
1082 /// and the inverse (using 'singular-value friendly' GAUSS pivot)
1087 double eps = 0.00000000000001;
1089 bool *bUnUsed = new bool[nGlo];
1090 double *diagV = new double[nGlo];
1091 double *rowMax = new double[nGlo];
1092 double *colMax = new double[nGlo];
1094 double *temp = new double[nGlo];
1096 for (Int_t i=0; i<nGlo; i++) {
1101 for (Int_t j=0; j<i; j++) {
1102 if (matV[j][i] == 0) {
1103 matV[j][i] = matV[i][j];
1108 // Small loop for matrix equilibration (gives a better conditioning)
1110 for (Int_t i=0; i<nGlo; i++) {
1111 for (Int_t j=0; j<nGlo; j++) {
1112 if (TMath::Abs(matV[i][j]) >= rowMax[i]) rowMax[i] = TMath::Abs(matV[i][j]); // Max elemt of row i
1113 if (TMath::Abs(matV[j][i]) >= colMax[i]) colMax[i] = TMath::Abs(matV[j][i]); // Max elemt of column i
1117 for (Int_t i=0; i<nGlo; i++) {
1118 if (0.0 != rowMax[i]) rowMax[i] = 1./rowMax[i]; // Max elemt of row i
1119 if (0.0 != colMax[i]) colMax[i] = 1./colMax[i]; // Max elemt of column i
1122 for (Int_t i=0; i<nGlo; i++) {
1123 for (Int_t j=0; j<nGlo; j++) {
1124 matV[i][j] = TMath::Sqrt(rowMax[i])*matV[i][j]*TMath::Sqrt(colMax[j]); // Equilibrate the V matrix
1126 diagV[i] = TMath::Abs(matV[i][i]); // save diagonal elem absolute values
1130 for (Int_t i=0; i<nGlo; i++) {
1134 for (Int_t j=0; j<nGlo; j++) { // First look for the pivot, ie max unused diagonal element
1135 if (bUnUsed[j] && (TMath::Abs(matV[j][j])>TMath::Max(TMath::Abs(vPivot),eps*diagV[j]))) {
1136 vPivot = matV[j][j];
1141 if (iPivot >= 0) { // pivot found
1143 bUnUsed[iPivot] = false; // This value is used
1144 vPivot = 1.0/vPivot;
1145 matV[iPivot][iPivot] = -vPivot; // Replace pivot by its inverse
1147 for (Int_t j=0; j<nGlo; j++) {
1148 for (Int_t jj=0; jj<nGlo; jj++) {
1149 if (j != iPivot && jj != iPivot) {// Other elements (!!! do them first as you use old matV[k][j]'s !!!)
1150 matV[j][jj] = matV[j][jj] - vPivot*matV[j][iPivot]*matV[iPivot][jj];
1155 for (Int_t j=0; j<nGlo; j++) {
1156 if (j != iPivot) { // Pivot row or column elements
1157 matV[j][iPivot] = matV[j][iPivot]*vPivot; // Column
1158 matV[iPivot][j] = matV[iPivot][j]*vPivot; // Line
1162 else { // No more pivot value (clear those elements)
1163 for (Int_t j=0; j<nGlo; j++) {
1167 for (Int_t k=0; k<nGlo; k++) {
1173 break; // No more pivots anyway, stop here
1177 for (Int_t i=0; i<nGlo; i++) {
1178 for (Int_t j=0; j<nGlo; j++) {
1179 matV[i][j] = TMath::Sqrt(colMax[i])*matV[i][j]*TMath::Sqrt(rowMax[j]); // Correct matrix V
1183 for (Int_t j=0; j<nGlo; j++) {
1186 for (Int_t jj=0; jj<nGlo; jj++) { // Reverse matrix elements
1187 matV[j][jj] = -matV[j][jj];
1188 temp[j] += matV[j][jj]*vecB[jj];
1192 for (Int_t j=0; j<nGlo; j++) {
1193 vecB[j] = temp[j]; // The final result
1206 // Same method but for local fit, so heavily simplified
1208 int AliMillepede::SpmInv(double matV[][fgkMaxLocalPar], double vecB[], int nLoc)
1210 /// Obtain solution of a system of linear equations with symmetric matrix
1211 /// and the inverse (using 'singular-value friendly' GAUSS pivot)
1216 double eps = 0.0000000000001;
1218 bool *bUnUsed = new bool[nLoc];
1219 double *diagV = new double[nLoc];
1220 double *temp = new double[nLoc];
1222 for (Int_t i=0; i<nLoc; i++) {
1224 diagV[i] = TMath::Abs(matV[i][i]); // save diagonal elem absolute values
1225 for (Int_t j=0; j<i; j++) {
1226 matV[j][i] = matV[i][j] ;
1230 for (Int_t i=0; i<nLoc; i++) {
1234 for (Int_t j=0; j<nLoc; j++) { // First look for the pivot, ie max unused diagonal element
1235 if (bUnUsed[j] && (TMath::Abs(matV[j][j])>TMath::Max(TMath::Abs(vPivot),eps*diagV[j]))) {
1236 vPivot = matV[j][j];
1241 if (iPivot >= 0) { // pivot found
1243 bUnUsed[iPivot] = false;
1244 vPivot = 1.0/vPivot;
1245 matV[iPivot][iPivot] = -vPivot; // Replace pivot by its inverse
1247 for (Int_t j=0; j<nLoc; j++) {
1249 for (Int_t jj=0; jj<=j; jj++) {
1250 if (jj != iPivot) {// Other elements (!!! do them first as you use old matV[k][j]'s !!!)
1251 matV[j][jj] = matV[j][jj] - vPivot*matV[j][iPivot]*matV[iPivot][jj];
1252 matV[jj][j] = matV[j][jj];
1258 for (Int_t j=0; j<nLoc; j++) {
1259 if (j != iPivot) { // Pivot row or column elements
1260 matV[j][iPivot] = matV[j][iPivot]*vPivot; // Column
1261 matV[iPivot][j] = matV[j][iPivot];
1265 else { // No more pivot value (clear those elements)
1266 for (Int_t j=0; j<nLoc; j++) {
1270 for (Int_t k=0; k<j; k++) {
1277 break; // No more pivots anyway, stop here
1281 for (Int_t j=0; j<nLoc; j++) {
1283 for (Int_t jj=0; jj<nLoc; jj++) { // Reverse matrix elements
1284 matV[j][jj] = -matV[j][jj];
1285 temp[j] += matV[j][jj]*vecB[jj];
1289 for (Int_t j=0; j<nLoc; j++) {
1302 -----------------------------------------------------------
1304 -----------------------------------------------------------
1306 multiply symmetric N-by-N matrix from the left with general M-by-N
1307 matrix and from the right with the transposed of the same general
1308 matrix to form symmetric M-by-M matrix.
1311 CALL SPAVAT(V,A,W,N,M) W = A * V * A
1314 where V = symmetric N-by-N matrix
1315 A = general N-by-M matrix
1316 W = symmetric M-by-M matrix
1317 -----------------------------------------------------------
1319 Int_t AliMillepede::SpAVAt(double matV[][fgkMaxLocalPar], double matA[][fgkMaxLocalPar], double matW[][fgkMaxGlobalPar], int nLoc, int nGlo)
1321 /// multiply symmetric N-by-N matrix from the left with general M-by-N
1322 /// matrix and from the right with the transposed of the same general
1323 /// matrix to form symmetric M-by-M matrix.
1325 for (Int_t i=0; i<nGlo; i++) {
1326 for (Int_t j=0; j<=i; j++) { // Matrix w is symmetric
1327 matW[i][j] = 0.0; // Reset final matrix
1328 for (Int_t k=0; k<nLoc; k++) {
1329 matW[i][j] += matA[i][k]*matV[k][k]*matA[j][k]; // diagonal terms of v
1330 for (Int_t l=0; l<k; l++) {
1331 matW[i][j] += matA[i][k]*matV[k][l]*matA[j][l]; // Use symmetric properties
1332 matW[i][j] += matA[i][l]*matV[k][l]*matA[j][k]; // of matrix v
1336 matW[j][i] = matW[i][j]; // Matrix w is symmetric
1346 -----------------------------------------------------------
1348 -----------------------------------------------------------
1350 multiply general M-by-N matrix A and N-vector X
1352 CALL SPAX(A,X,Y,M,N) Y = A * X
1355 where A = general M-by-N matrix (A11 A12 ... A1N A21 A22 ...)
1358 -----------------------------------------------------------
1360 Int_t AliMillepede::SpAX(double matA[][fgkMaxLocalPar], double vecX[], double vecY[], int nCol, int nRow)
1362 /// multiply general M-by-N matrix A and N-vector X
1363 for (Int_t i=0; i<nRow; i++) {
1364 vecY[i] = 0.0; // Reset final vector
1365 for (Int_t j=0; j<nCol; j++) {
1366 vecY[i] += matA[i][j]*vecX[j]; // fill the vector
1375 -----------------------------------------------------------
1377 -----------------------------------------------------------
1379 Print the final results into the logfile
1381 -----------------------------------------------------------
1383 Int_t AliMillepede::PrintGlobalParameters() const
1385 /// Print the final results into the logfile
1387 double lGlobalCor =0.;
1390 AliInfo(" Result of fit for global parameters");
1391 AliInfo(" ===================================");
1392 AliInfo(" I initial final differ lastcor error gcor");
1393 AliInfo("-----------------------------------------------------------------------------------");
1395 for (int i=0; i<fNGlobalPar; i++) {
1396 lError = TMath::Sqrt(TMath::Abs(fMatCGlo[i][i]));
1397 if (fMatCGlo[i][i] < 0.0) lError = -lError;
1400 if (TMath::Abs(fMatCGlo[i][i]*fDiagCGlo[i]) > 0) {
1401 lGlobalCor = TMath::Sqrt(TMath::Abs(1.0-1.0/(fMatCGlo[i][i]*fDiagCGlo[i])));
1402 AliInfo(Form("%d\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f",
1403 i,fInitPar[i],fInitPar[i]+fDeltaPar[i],fDeltaPar[i],fVecBGlo[i],lError,lGlobalCor));
1406 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]));
1414 ----------------------------------------------------------------
1415 CHI2DOFLIM: return the limit in chi^2/nd for n sigmas stdev authorized
1416 ----------------------------------------------------------------
1418 Only n=1, 2, and 3 are expected in input
1419 ----------------------------------------------------------------
1421 double AliMillepede::Chi2DoFLim(int nSig, int nDoF)
1423 /// return the limit in chi^2/nd for n sigmas stdev authorized
1425 double sn[3] = {0.47523, 1.690140, 2.782170};
1426 double table[3][30] = {{1.0000, 1.1479, 1.1753, 1.1798, 1.1775, 1.1730, 1.1680, 1.1630,
1427 1.1581, 1.1536, 1.1493, 1.1454, 1.1417, 1.1383, 1.1351, 1.1321,
1428 1.1293, 1.1266, 1.1242, 1.1218, 1.1196, 1.1175, 1.1155, 1.1136,
1429 1.1119, 1.1101, 1.1085, 1.1070, 1.1055, 1.1040},
1430 {4.0000, 3.0900, 2.6750, 2.4290, 2.2628, 2.1415, 2.0481, 1.9736,
1431 1.9124, 1.8610, 1.8171, 1.7791, 1.7457, 1.7161, 1.6897, 1.6658,
1432 1.6442, 1.6246, 1.6065, 1.5899, 1.5745, 1.5603, 1.5470, 1.5346,
1433 1.5230, 1.5120, 1.5017, 1.4920, 1.4829, 1.4742},
1434 {9.0000, 5.9146, 4.7184, 4.0628, 3.6410, 3.3436, 3.1209, 2.9468,
1435 2.8063, 2.6902, 2.5922, 2.5082, 2.4352, 2.3711, 2.3143, 2.2635,
1436 2.2178, 2.1764, 2.1386, 2.1040, 2.0722, 2.0428, 2.0155, 1.9901,
1437 1.9665, 1.9443, 1.9235, 1.9040, 1.8855, 1.8681}};
1443 lNSig = TMath::Max(1,TMath::Min(nSig,3));
1446 return table[lNSig-1][nDoF-1];
1448 else { // approximation
1449 return ((sn[lNSig-1]+TMath::Sqrt(float(2*nDoF-3)))*
1450 (sn[lNSig-1]+TMath::Sqrt(float(2*nDoF-3))))/float(2*nDoF-2);