#include <TStopwatch.h>
#include <TFile.h>
#include <TMath.h>
+#include <TVectorD.h>
#include "AliMatrixSq.h"
#include "AliSymMatrix.h"
#include "AliRectMatrix.h"
fRecord->AddResidual(lMeas);
//
// Retrieve local param interesting indices
- for (int i=0;i<fNLocPar;i++) if (derlc[i]!=0.0) {fRecord->AddIndexValue(i,derlc[i]); derlc[i] = 0.0;}
+ for (int i=0;i<fNLocPar;i++) if (!IsZero(derlc[i])) {fRecord->AddIndexValue(i,derlc[i]); derlc[i] = 0.0;}
//
fRecord->AddWeight( 1.0/lSigma/lSigma );
//
// Idem for global parameters
- for (int i=0;i<fNGloPar;i++) if (dergb[i]!=0.0) {
+ for (int i=0;i<fNGloPar;i++) if (!IsZero(dergb[i])) {
fRecord->AddIndexValue(i,dergb[i]); dergb[i] = 0.0;
fRecord->MarkGroup(fParamGrID[i]);
}
fRecord->AddResidual(lMeas);
//
// Retrieve local param interesting indices
- for (int i=0;i<nlc;i++) if (derlc[i]!=0.0) {fRecord->AddIndexValue(indlc[i],derlc[i]); derlc[i]=0.; indlc[i]=0;}
+ for (int i=0;i<nlc;i++) if (!IsZero(derlc[i])) {fRecord->AddIndexValue(indlc[i],derlc[i]); derlc[i]=0.; indlc[i]=0;}
//
fRecord->AddWeight( 1./lSigma/lSigma );
//
// Idem for global parameters
- for (int i=0;i<ngb;i++) if (dergb[i]!=0.0) {fRecord->AddIndexValue(indgb[i],dergb[i]); dergb[i]=0.; indgb[i]=0;}
+ for (int i=0;i<ngb;i++) if (!IsZero(dergb[i])) {fRecord->AddIndexValue(indgb[i],dergb[i]); dergb[i]=0.; indgb[i]=0;}
//
}
fRecord->Reset();
fRecord->AddResidual(val);
fRecord->AddWeight(sigma);
- for (int i=0; i<fNGloPar; i++) if (dergb[i]!=0) fRecord->AddIndexValue(i,dergb[i]);
+ for (int i=0; i<fNGloPar; i++) if (!IsZero(dergb[i])) fRecord->AddIndexValue(i,dergb[i]);
fNGloConstraints++;
- if (sigma==0) fNLagrangeConstraints++;
+ if (IsZero(sigma)) fNLagrangeConstraints++;
// printf("NewConstraint:\n"); fRecord->Print(); //RRR
SaveRecordConstraint();
//
fRecord->Reset();
fRecord->AddResidual(val);
fRecord->AddWeight(sigma); // dummy
- for (int i=0; i<ngb; i++) if (dergb[i]!=0) fRecord->AddIndexValue(indgb[i],dergb[i]);
+ for (int i=0; i<ngb; i++) if (!IsZero(dergb[i])) fRecord->AddIndexValue(indgb[i],dergb[i]);
fNGloConstraints++;
- if (sigma==0) fNLagrangeConstraints++;
+ if (IsZero(sigma)) fNLagrangeConstraints++;
SaveRecordConstraint();
//
}
}
double vl;
//
+ Int_t maxLocUsed = 0;
+ //
for (int ip=nPoints;ip--;) { // Transfer the measurement records to matrices
double resid = fRecord->GetValue( refLoc[ip]-1 );
double weight = fRecord->GetValue( refGlo[ip]-1 );
// Symmetric matrix, don't bother j>i coeffs
for (int i=nrefLoc[ip];i--;) { // Fill local matrix and vector
fVecBLoc[ indLoc[i] ] += weight*resid*derLoc[i];
+ if (indLoc[i]>maxLocUsed) maxLocUsed = indLoc[i];
for (int j=i+1;j--;) matCLoc(indLoc[i] ,indLoc[j]) += weight*derLoc[i]*derLoc[j];
- }
+ }
//
} // end of the transfer of the measurement record to matrices
//
+ matCLoc.SetSizeUsed(++maxLocUsed); // data with B=0 may use less than declared nLocals
+ //
// first try to solve by faster Cholesky decomposition, then by Gaussian elimination
if (!matCLoc.SolveChol(fVecBLoc,kTRUE)) {
AliInfo("Failed to solve locals by Cholesky, traying Gaussian Elimination");
//
// If requested, store the track params and errors
// printf("locfit: "); for (int i=0;i<fNLocPar;i++) printf("%+e |",fVecBLoc[i]); printf("\n");
- if (localParams) for (int i=fNLocPar; i--;) {
+ if (localParams) for (int i=maxLocUsed; i--;) {
localParams[2*i] = fVecBLoc[i];
localParams[2*i+1] = TMath::Sqrt(TMath::Abs(matCLoc.QueryDiag(i)));
}
if ( (absres >= fResCutInit && fIter ==1 ) ||
(absres >= fResCut && fIter > 1)) {
if (fLocFitAdd) fNLocFitsRejected++;
+ // printf("reject res %5ld %+e\n",fCurrRecDataID,resid);
return 0;
}
//
nEq++; // number of equations
} // end of Calculate residuals
//
- int nDoF = nEq-fNLocPar;
+ int nDoF = nEq-maxLocUsed;
lChi2 = (nDoF>0) ? lChi2/nDoF : 0; // Chi^2/dof
//
if (fNStdDev != 0 && nDoF>0 && lChi2 > Chi2DoFLim(fNStdDev,nDoF)*fChi2CutFactor) { // check final chi2
if (fLocFitAdd) fNLocFitsRejected++;
+ // printf("reject chi2 %5ld: %+e\n",fCurrRecDataID, lChi2);
return 0;
}
//
for (int jg=ig+1;jg--;) { // matCGlo is symmetric by construction
int jIDg = indGlo[jg];
if ( fSigmaPar[jIDg] <= 0.) continue; // fixed parameter RRRCheck
- if ( (vl = weight*derGlo[ig]*derGlo[jg])!=0 ) {
+ if ( !IsZero(vl = weight*derGlo[ig]*derGlo[jg]) ) {
fFillIndex[nfill] = jIDg;
fFillValue[nfill++] = fLocFitAdd ? vl:-vl;
}
int iCIDg = fGlo2CGlo[iIDg]; // compressed Index of index
if (iCIDg == -1) {
Double_t *rowGL = matCGloLoc(nGloInFit);
- for (int k=fNLocPar;k--;) rowGL[k] = 0.0; // reset the row
+ for (int k=maxLocUsed;k--;) rowGL[k] = 0.0; // reset the row
iCIDg = fGlo2CGlo[iIDg] = nGloInFit;
fCGlo2Glo[nGloInFit++] = iIDg;
}
//
vl = 0;
Double_t *rowGLIDg = matCGloLoc(iCIDg);
- for (int kl=0;kl<fNLocPar;kl++) if (rowGLIDg[kl]) vl += rowGLIDg[kl]*fVecBLoc[kl];
- if (vl!=0) fVecBGlo[iIDg] -= fLocFitAdd ? vl : -vl;
+ for (int kl=0;kl<maxLocUsed;kl++) if (rowGLIDg[kl]) vl += rowGLIDg[kl]*fVecBLoc[kl];
+ if (!IsZero(vl)) fVecBGlo[iIDg] -= fLocFitAdd ? vl : -vl;
//
int nfill = 0;
for (int jCIDg=0;jCIDg<=iCIDg; jCIDg++) {
//
vl = 0;
Double_t *rowGLJDg = matCGloLoc(jCIDg);
- for (int kl=0;kl<fNLocPar;kl++) {
+ for (int kl=0;kl<maxLocUsed;kl++) {
// diag terms
- if ( (vll=rowGLIDg[kl]*rowGLJDg[kl])!=0 ) vl += matCLoc.QueryDiag(kl)*vll;
+ if ( (!IsZero(vll=rowGLIDg[kl]*rowGLJDg[kl])) ) vl += matCLoc.QueryDiag(kl)*vll;
//
// off-diag terms
for (int ll=0;ll<kl;ll++) {
- if ( (vll=rowGLIDg[kl]*rowGLJDg[ll])!=0 ) vl += matCLoc(kl,ll)*vll;
- if ( (vll=rowGLIDg[ll]*rowGLJDg[kl])!=0 ) vl += matCLoc(kl,ll)*vll;
+ if ( !IsZero(vll=rowGLIDg[kl]*rowGLJDg[ll]) ) vl += matCLoc(kl,ll)*vll;
+ if ( !IsZero(vll=rowGLIDg[ll]*rowGLJDg[kl]) ) vl += matCLoc(kl,ll)*vll;
}
}
- if (vl!=0) {
+ if (!IsZero(vl)) {
fFillIndex[nfill] = jIDg;
fFillValue[nfill++] = fLocFitAdd ? -vl : vl;
}
res = GlobalFitIteration();
if (!res) break;
//
- if (fChi2CutFactor != fChi2CutRef) {
+ if (!IsZero(fChi2CutFactor-fChi2CutRef)) {
fChi2CutFactor = TMath::Sqrt(fChi2CutFactor);
if (fChi2CutFactor < 1.2*fChi2CutRef) {
fChi2CutFactor = fChi2CutRef;
fNLagrangeConstraints = 0;
for (int i=0; i<fNGloConstraints; i++) {
ReadRecordConstraint(i);
- if ( fRecord->GetValue(1)==0 ) fNLagrangeConstraints++; // exact constraint (no error) -> Lagrange multiplier
+ if ( IsZero(fRecord->GetValue(1)) ) fNLagrangeConstraints++; // exact constraint (no error) -> Lagrange multiplier
}
//
// if needed, readjust the size of the global vector (for matrices this is done automatically)
for (Long_t i=0;i<ndr;i++) {
ReadRecordData(i);
LocalFit();
- if ( (i%int(0.1*ndr)) == 0) printf("%.1f%% of local fits done\n", double(100.*i)/ndr);
+ if ( (i%int(0.2*ndr)) == 0) printf("%.1f%% of local fits done\n", double(100.*i)/ndr);
}
swt.Stop();
printf("%ld local fits done: ", ndr);
// add large number to diagonal of fixed params
//
for (int i=fNGloPar;i--;) { // // Reset row and column of fixed params and add 1/sig^2 to free ones
+ // printf("#%3d : Nproc : %5d grp: %d\n",i,fProcPnt[i],fParamGrID[i]);
if (fProcPnt[i]<1) {
fNGloFix++;
fVecBGlo[i] = 0.;
for (int ic=0;ic<=ir;ic++) { // matrix is symmetric
int jID = indV[ic];
double vl = der[ir]*der[ic]*sig2i;
- if (vl!=0) matCGlo(iID,jID) += vl;
+ if (!IsZero(vl)) matCGlo(iID,jID) += vl;
}
fVecBGlo[iID] += val*der[ir]*sig2i;
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff