#include <TArrayI.h>
#include <TArrayD.h>
#include <TMath.h>
#include "AliLog.h"
#include "AliMillepede.h"
AliMillepede::AliMillepede()
: TObject(),
fIndexLocEq(fgkMaxGlobalPar+fgkMaxLocalPar),
fDerivLocEq(fgkMaxGlobalPar+fgkMaxLocalPar),
fIndexAllEqs(1000*fgkMaxGlobalPar+fgkMaxLocalPar),
fDerivAllEqs(1000*fgkMaxGlobalPar+fgkMaxLocalPar),
fLocEqPlace(1000),
fNIndexLocEq(0),
fNDerivLocEq(0),
fNIndexAllEqs(0),
fNDerivAllEqs(0),
fNLocEqPlace(0),
fNLocalEquations(0),
fResCutInit(0.),
fResCut(0.),
fChi2CutFactor(1.0),
fChi2CutRef(1.0),
fIter(0),
fMaxIter(10),
fNStdDev(3),
fNGlobalConstraints(0),
fNLocalFits(0),
fNLocalFitsRejected(0),
fNGlobalPar(0),
fNLocalPar(0)
{
AliInfo(" ");
AliInfo(" * o o o ");
AliInfo(" o o o ");
AliInfo(" o ooooo o o o oo ooo oo ooo oo ");
AliInfo(" o o o o o o o o o o o o o o o o ");
AliInfo(" o o o o o o oooo o o oooo o o oooo ");
AliInfo(" o o o o o o o ooo o o o o ");
AliInfo(" o o o o o o oo o oo ooo oo ++ starts");
AliInfo(" ");
}
AliMillepede::~AliMillepede() {
}
Int_t AliMillepede::InitMille(int nGlo, int nLoc, int lNStdDev,
double lResCut, double lResCutInit)
{
AliDebug(1,"");
AliDebug(1,"----------------------------------------------------");
AliDebug(1,"");
AliDebug(1," Entering InitMille");
AliDebug(1,"");
AliDebug(1,"-----------------------------------------------------");
AliDebug(1,"");
fNGlobalConstraints = 0;
fNLocalFits = 0;
fNLocalFitsRejected = 0;
fChi2CutRef = 1.0;
AliMillepede::SetNLocalEquations(0);
fIter = 0;
fMaxIter = 10;
fResCutInit = lResCutInit;
fResCut = lResCut;
fNGlobalPar = nGlo;
fNLocalPar = nLoc;
fNStdDev = lNStdDev;
AliDebug(1,Form("Number of global parameters : %d", fNGlobalPar));
AliDebug(1,Form("Number of local parameters : %d", fNLocalPar));
AliDebug(1,Form("Number of standard deviations : %d", fNStdDev));
if (fNGlobalPar>fgkMaxGlobalPar || fNLocalPar>fgkMaxLocalPar) {
AliDebug(1,"Two many parameters !!!!!");
return 0;
}
for (int i=0; i<fNGlobalPar; i++) {
fVecBGlo[i]=0.;
fInitPar[i]=0.;
fDeltaPar[i]=0.;
fSigmaPar[i]=-1.;
fIsNonLinear[i]=0;
fGlo2CGLRow[i]=-1;
fCGLRow2Glo[i]=-1;
for (int j=0; j<fNGlobalPar;j++) {
fMatCGlo[i][j]=0.;
}
}
for (int i=0; i<fNLocalPar; i++) {
fVecBLoc[i]=0.;
for (int j=0; j<fNLocalPar;j++) {
fMatCLoc[i][j]=0.;
}
}
for (int j=0; j<fNGlobalPar; j++) {
AliMillepede::SetParSigma(j,0.0);
}
fDerivLocEq.Reset(); fNDerivLocEq=0;
fIndexLocEq.Reset(); fNIndexLocEq=0;
fIndexAllEqs.Reset(); fNIndexAllEqs=0;
fDerivAllEqs.Reset(); fNDerivAllEqs=0;
fLocEqPlace.Reset(); fNLocEqPlace=0;
AliDebug(1,"");
AliDebug(1,"----------------------------------------------------");
AliDebug(1,"");
AliDebug(1," InitMille has been successfully called!");
AliDebug(1,"");
AliDebug(1,"-----------------------------------------------------");
AliDebug(1,"");
return 1;
}
Int_t AliMillepede::SetGlobalParameters(double *param)
{
for(Int_t iPar=0; iPar<fNGlobalPar; iPar++){
fInitPar[iPar] = param[iPar];
}
return 1;
}
Int_t AliMillepede::SetGlobalParameter(int iPar, double param)
{
if (iPar<0 || iPar>=fNGlobalPar) {
return 0;
}
else {
fInitPar[iPar] = param;
}
return 1;
}
Int_t AliMillepede::SetParSigma(int iPar, double sigma)
{
if (iPar>=fNGlobalPar) {
return 0;
}
else {
fSigmaPar[iPar] = sigma;
}
return 1;
}
Int_t AliMillepede::SetNonLinear(int iPar)
{
if (iPar<0 || iPar>=fNGlobalPar) {
return 0;
}
else {
fIsNonLinear[iPar] = 1;
}
return 1;
}
Int_t AliMillepede::SetIterations(double lChi2CutFac)
{
fChi2CutFactor = TMath::Max(1.0, lChi2CutFac);
AliInfo(Form("Initial cut factor is %f",fChi2CutFactor));
fIter = 1;
return 1;
}
Int_t AliMillepede::SetGlobalConstraint(double dercs[], double lLagMult)
{
if (fNGlobalConstraints>=fgkMaxGloCsts) {
AliInfo("Too many constraints !!!");
return 0;
}
for (int i=0; i<fNGlobalPar; i++) {
fMatDerConstr[fNGlobalConstraints][i] = dercs[i];
}
fLagMult[fNGlobalConstraints] = lLagMult;
fNGlobalConstraints++ ;
AliInfo(Form("Number of constraints increased to %d",fNGlobalConstraints));
return 1;
}
Int_t AliMillepede::SetLocalEquation(double dergb[], double derlc[], double lMeas, double lSigma)
{
if (lSigma<=0.0) {
for (int i=0; i<fNLocalPar; i++) {
derlc[i] = 0.0;
}
for (int i=0; i<fNGlobalPar; i++) {
dergb[i] = 0.0;
}
return 1;
}
double lWeight = 1.0/(lSigma*lSigma);
int iLocFirst = -1;
int iLocLast = -1;
int iGlobFirst = -1;
int iGlobLast = -1;
for (int i=0; i<fNLocalPar; i++) {
if (derlc[i]!=0.0) {
if (iLocFirst == -1) {
iLocFirst = i;
}
iLocLast = i;
}
}
AliDebug(2,Form("%d / %d",iLocFirst, iLocLast));
for (int i=0; i<fNGlobalPar; i++) {
if (dergb[i]!=0.0) {
if (iGlobFirst == -1) {
iGlobFirst = i;
}
iGlobLast = i;
}
}
AliDebug(2,Form("%d / %d",iGlobFirst,iGlobLast));
if (fNIndexLocEq==fIndexLocEq.GetSize()) fIndexLocEq.Set(2*fNIndexLocEq);
fIndexLocEq.AddAt(-1,fNIndexLocEq++);
if (fNDerivLocEq==fDerivLocEq.GetSize()) fDerivLocEq.Set(2*fNDerivLocEq);
fDerivLocEq.AddAt(lMeas,fNDerivLocEq++);
for (int i=iLocFirst; i<=iLocLast; i++) {
if (derlc[i]!=0.0) {
if (fNIndexLocEq==fIndexLocEq.GetSize()) fIndexLocEq.Set(2*fNIndexLocEq);
fIndexLocEq.AddAt(i,fNIndexLocEq++);
if (fNDerivLocEq==fDerivLocEq.GetSize()) fDerivLocEq.Set(2*fNDerivLocEq);
fDerivLocEq.AddAt(derlc[i],fNDerivLocEq++);
derlc[i] = 0.0;
}
}
if (fNIndexLocEq==fIndexLocEq.GetSize()) fIndexLocEq.Set(2*fNIndexLocEq);
fIndexLocEq.AddAt(-1,fNIndexLocEq++);
if (fNDerivLocEq==fDerivLocEq.GetSize()) fDerivLocEq.Set(2*fNDerivLocEq);
fDerivLocEq.AddAt(lWeight,fNDerivLocEq++);
for (int i=iGlobFirst; i<=iGlobLast; i++) {
if (dergb[i]!=0.0) {
if (fNIndexLocEq==fIndexLocEq.GetSize()) fIndexLocEq.Set(2*fNIndexLocEq);
fIndexLocEq.AddAt(i,fNIndexLocEq++);
if (fNDerivLocEq==fDerivLocEq.GetSize()) fDerivLocEq.Set(2*fNDerivLocEq);
fDerivLocEq.AddAt(dergb[i],fNDerivLocEq++);
dergb[i] = 0.0;
}
}
AliDebug(2,Form("Out Equloc -- NST = %d",fNDerivLocEq));
return 1;
}
Int_t AliMillepede::LocalFit(int iFit, double localParams[], Bool_t bSingleFit)
{
int iEqTerm = 0;
int i, j, k;
int iIdx, jIdx, kIdx;
int iIdxIdx, jIdxIdx;
int iMeas = -1;
int iWeight = 0;
int iLocFirst = 0;
int iLocLast = 0;
int iGloFirst = 0;
int iGloLast = 0;
int nGloInFit = 0;
double lMeas = 0.0;
double lWeight = 0.0;
double lChi2 = 0.0;
double lRedChi2 = 0.0;
double lChi2Cut = 0.0;
int nEq = 0;
int nDoF = 0;
int nEqTerms = fNDerivLocEq;
if (fIter < 2 && !bSingleFit) {
AliDebug(1,Form("Store equation no: %d", iFit));
for (i=0; i<nEqTerms; i++) {
if (fNIndexAllEqs==fIndexAllEqs.GetSize()) fIndexAllEqs.Set(2*fNIndexAllEqs);
fIndexAllEqs.AddAt(fIndexLocEq[i],fNIndexAllEqs++);
if (fNDerivAllEqs==fDerivAllEqs.GetSize()) fDerivAllEqs.Set(2*fNDerivAllEqs);
fDerivAllEqs.AddAt(fDerivLocEq[i],fNDerivAllEqs++);
}
if (fNLocEqPlace==fLocEqPlace.GetSize()) fLocEqPlace.Set(2*fNLocEqPlace);
fLocEqPlace.AddAt(fNIndexAllEqs,fNLocEqPlace++);
AliDebug(2,Form("FLocEqPlace size = %d",fLocEqPlace[iFit]));
AliDebug(2,Form("FIndexAllEqs size = %d",fNIndexAllEqs));
}
for (i=0; i<fNLocalPar; i++) {
fVecBLoc[i] = 0.0;
for (j=0; j<fNLocalPar; j++) {
fMatCLoc[i][j] = 0.0;
}
}
for (i=0; i<fNGlobalPar; i++) {
fGlo2CGLRow[i] = -1;
}
iEqTerm = 0;
if (fNIndexLocEq==fIndexLocEq.GetSize()) fIndexLocEq.Set(2*fNIndexLocEq);
fIndexLocEq.AddAt(-1,fNIndexLocEq++);
while (iEqTerm <= nEqTerms) {
if (fIndexLocEq[iEqTerm] == -1) {
if (iMeas == -1) {
iMeas = iEqTerm;
iLocFirst = iEqTerm+1;
}
else if (iWeight == 0) {
iWeight = iEqTerm;
iLocLast = iEqTerm-1;
iGloFirst = iEqTerm+1;
}
else {
iGloLast = iEqTerm-1;
lMeas = fDerivLocEq[iMeas];
lWeight = fDerivLocEq[iWeight];
for (i=iGloFirst; i<=iGloLast; i++) {
iIdx = fIndexLocEq[i];
if (fIsNonLinear[iIdx] == 0)
lMeas -= fDerivLocEq[i]*(fInitPar[iIdx]+fDeltaPar[iIdx]);
else
lMeas -= fDerivLocEq[i]*(fDeltaPar[iIdx]);
}
for (i=iLocFirst; i<=iLocLast; i++) {
iIdx = fIndexLocEq[i];
fVecBLoc[iIdx] += lWeight*lMeas*fDerivLocEq[i];
for (j=iLocFirst; j<=i ; j++) {
jIdx = fIndexLocEq[j];
fMatCLoc[iIdx][jIdx] += lWeight*fDerivLocEq[i]*fDerivLocEq[j];
}
}
iMeas = -1;
iWeight = 0;
iEqTerm--;
}
}
iEqTerm++;
}
Int_t nRank = AliMillepede::SpmInv(fMatCLoc, fVecBLoc, fNLocalPar);
AliDebug(1,"");
AliDebug(1," __________________________________________________");
AliDebug(1,Form(" Printout of local fit (FITLOC) with rank= %d", nRank));
AliDebug(1," Result of local fit : (index/parameter/error)");
for (i=0; i<fNLocalPar; i++) {
AliDebug(1,Form("%d / %.6f / %.6f", i, fVecBLoc[i], TMath::Sqrt(fMatCLoc[i][i])));
}
for (i=0; i<fNLocalPar; i++) {
localParams[2*i] = fVecBLoc[i];
localParams[2*i+1] = TMath::Sqrt(TMath::Abs(fMatCLoc[i][i]));
}
iEqTerm = 0;
iMeas = -1;
iWeight = 0;
while (iEqTerm <= nEqTerms) {
if (fIndexLocEq[iEqTerm] == -1) {
if (iMeas == -1) {
iMeas = iEqTerm;
iLocFirst = iEqTerm+1;
}
else if (iWeight == 0) {
iWeight = iEqTerm;
iLocLast = iEqTerm-1;
iGloFirst = iEqTerm+1;
}
else {
iGloLast = iEqTerm-1;
lMeas = fDerivLocEq[iMeas];
lWeight = fDerivLocEq[iWeight];
for (i=iLocFirst; i<=iLocLast; i++) {
iIdx = fIndexLocEq[i];
lMeas -= fDerivLocEq[i]*fVecBLoc[iIdx];
}
for (i=iGloFirst; i<=iGloLast; i++) {
iIdx = fIndexLocEq[i];
if (fIsNonLinear[iIdx] == 0)
lMeas -= fDerivLocEq[i]*(fInitPar[iIdx]+fDeltaPar[iIdx]);
else
lMeas -= fDerivLocEq[i]*(fDeltaPar[iIdx]);
}
AliDebug(2,Form("Residual value : %.6f", lMeas));
if (TMath::Abs(lMeas) >= fResCutInit && fIter <= 1) {
AliDebug(2,"Rejected track !!!!!");
fNLocalFitsRejected++;
fIndexLocEq.Reset(); fNIndexLocEq=0;
fDerivLocEq.Reset(); fNDerivLocEq=0;
return 0;
}
if (TMath::Abs(lMeas) >= fResCut && fIter > 1) {
AliDebug(2,"Rejected track !!!!!");
fNLocalFitsRejected++;
fIndexLocEq.Reset(); fNIndexLocEq=0;
fDerivLocEq.Reset(); fNDerivLocEq=0;
return 0;
}
lChi2 += lWeight*lMeas*lMeas ;
nEq++;
iMeas = -1;
iWeight = 0;
iEqTerm--;
}
}
iEqTerm++;
}
nDoF = nEq-nRank;
lRedChi2 = 0.0;
AliDebug(1,Form("Final chi square / degrees of freedom %.2f / %d",lChi2, nDoF));
if (nDoF > 0) lRedChi2 = lChi2/float(nDoF);
fNLocalFits++;
if (fNStdDev != 0 && nDoF > 0 && !bSingleFit)
{
lChi2Cut = AliMillepede::Chi2DoFLim(fNStdDev, nDoF)*fChi2CutFactor;
AliDebug(1,Form("Reject if Chisq/Ndf = %.4f > %.4f",lRedChi2,lChi2Cut));
if (lRedChi2 > lChi2Cut)
{
AliDebug(2,"Rejected track !!!!!");
fNLocalFitsRejected++;
fIndexLocEq.Reset(); fNIndexLocEq=0;
fDerivLocEq.Reset(); fNDerivLocEq=0;
return 0;
}
}
if (bSingleFit)
{
fIndexLocEq.Reset(); fNIndexLocEq=0;
fDerivLocEq.Reset(); fNDerivLocEq=0;
return 1;
}
iEqTerm = 0;
iMeas = -1;
iWeight = 0;
while (iEqTerm <= nEqTerms)
{
if (fIndexLocEq[iEqTerm] == -1)
{
if (iMeas == -1) {
iMeas = iEqTerm;
iLocFirst = iEqTerm+1;
}
else if (iWeight == 0) {
iWeight = iEqTerm;
iLocLast = iEqTerm-1;
iGloFirst = iEqTerm+1;
}
else {
iGloLast = iEqTerm-1;
lMeas = fDerivLocEq[iMeas];
lWeight = fDerivLocEq[iWeight];
for (i=iGloFirst; i<=iGloLast; i++) {
iIdx = fIndexLocEq[i];
if (fIsNonLinear[iIdx] == 0)
lMeas -= fDerivLocEq[i]*(fInitPar[iIdx]+fDeltaPar[iIdx]);
else
lMeas -= fDerivLocEq[i]*(fDeltaPar[iIdx]);
}
for (i=iGloFirst; i<=iGloLast; i++) {
iIdx = fIndexLocEq[i];
fVecBGlo[iIdx] += lWeight*lMeas*fDerivLocEq[i];
for (j=iGloFirst; j<=iGloLast; j++) {
jIdx = fIndexLocEq[j];
fMatCGlo[iIdx][jIdx] += lWeight*fDerivLocEq[i]*fDerivLocEq[j];
}
iIdxIdx = fGlo2CGLRow[iIdx];
if (iIdxIdx == -1) {
for (k=0; k<fNLocalPar; k++) {
fMatCGloLoc[nGloInFit][k] = 0.0;
}
fGlo2CGLRow[iIdx] = nGloInFit;
fCGLRow2Glo[nGloInFit] = iIdx;
iIdxIdx = nGloInFit;
nGloInFit++;
}
for (k=iLocFirst; k<=iLocLast ; k++) {
kIdx = fIndexLocEq[k];
fMatCGloLoc[iIdxIdx][kIdx] += lWeight*fDerivLocEq[i]*fDerivLocEq[k];
}
}
iMeas = -1;
iWeight = 0;
iEqTerm--;
}
}
iEqTerm++;
}
AliMillepede::SpAVAt(fMatCLoc, fMatCGloLoc, fMatCGloCorr, fNLocalPar, nGloInFit);
AliMillepede::SpAX(fMatCGloLoc, fVecBLoc, fVecBGloCorr, fNLocalPar, nGloInFit);
for (iIdxIdx=0; iIdxIdx<nGloInFit; iIdxIdx++) {
iIdx = fCGLRow2Glo[iIdxIdx];
fVecBGlo[iIdx] -= fVecBGloCorr[iIdxIdx];
for (jIdxIdx=0; jIdxIdx<=iIdxIdx; jIdxIdx++) {
jIdx = fCGLRow2Glo[jIdxIdx];
fMatCGlo[iIdx][jIdx] -= fMatCGloCorr[iIdxIdx][jIdxIdx];
fMatCGlo[jIdx][iIdx] = fMatCGlo[iIdx][jIdx];
}
}
fIndexLocEq.Reset(); fNIndexLocEq=0;
fDerivLocEq.Reset(); fNDerivLocEq=0;
return 1;
}
Int_t AliMillepede::GlobalFit(double par[], double error[], double pull[])
{
int i, j;
int nVar = 0;
int nGloFix = 0;
double lConstraint;
double step[fgkMaxGlobalPar]={0};
double localPars[2*fgkMaxLocalPar];
int nLocFitsGood = 0;
int nLocFitsTot = 0;
int nLocFits = 0;
AliInfo("..... Making global fit .....");
AliInfo(Form("First Iteration with cut factor %.2f", fChi2CutFactor));
nLocFitsTot = AliMillepede::GetNLocalEquations();
while (fIter < fMaxIter)
{
nLocFits = AliMillepede::GetNLocalEquations();
AliInfo(Form("...using %d local fits...",nLocFits));
for (i=0; i<fNGlobalPar; i++) {
fDiagCGlo[i] = fMatCGlo[i][i];
}
nGloFix = 0;
for (i=0; i<fNGlobalPar; i++) {
if (fSigmaPar[i] <= 0.0) {
nGloFix++;
for (j=0; j<fNGlobalPar; j++) {
fMatCGlo[i][j] = 0.0;
fMatCGlo[j][i] = 0.0;
}
}
else {
fMatCGlo[i][i] += 1.0/(fSigmaPar[i]*fSigmaPar[i]);
}
}
nVar = fNGlobalPar;
AliDebug(1,Form("Number of constraint equations : %d", fNGlobalConstraints));
for (i=0; i<fNGlobalConstraints; i++) {
lConstraint = fLagMult[i];
for (j=0; j<fNGlobalPar; j++) {
fMatCGlo[nVar][j] = float(nLocFits)*fMatDerConstr[i][j];
fMatCGlo[j][nVar] = float(nLocFits)*fMatDerConstr[i][j];
lConstraint -= fMatDerConstr[i][j]*(fInitPar[j]+fDeltaPar[j]);
}
fMatCGlo[nVar][nVar] = 0.0;
fVecBGlo[nVar] = float(nLocFits)*lConstraint;
nVar++;
}
double lFinalCor = 0.0;
if (fIter > 1) {
for (i=0; i<fNGlobalPar; i++) {
for (j=0; j<fNGlobalPar; j++) {
lFinalCor += step[i]*fMatCGlo[i][j]*step[j];
if (i == j && fSigmaPar[i] != 0) {
lFinalCor -= step[i]*step[i]/(fSigmaPar[i]*fSigmaPar[i]);
}
}
}
}
AliInfo(Form(" Final coeff is %.6f",lFinalCor));
AliInfo(Form(" Final NDOFs = %d", fNGlobalPar));
Int_t nRank = AliMillepede::SpmInv(fMatCGlo, fVecBGlo, nVar);
for (i=0; i<fNGlobalPar; i++) {
fDeltaPar[i] += fVecBGlo[i];
AliDebug(1,Form("fDeltaPar[%d] = %.6f", i, fDeltaPar[i]));
AliDebug(1,Form("fMatCGlo[%d][%d] = %.6f", i, i, fMatCGlo[i][i]));
AliDebug(1,Form("err = %.6f", TMath::Sqrt(TMath::Abs(fMatCGlo[i][i]))));
step[i] = fVecBGlo[i];
if (fIter == 1) error[i] = fMatCGlo[i][i];
}
AliInfo(Form("The rank defect of the symmetric %d by %d matrix is %d (bad if non 0)",
nVar, nVar, nVar-nGloFix-nRank));
AliInfo(Form("Total : %d local fits, %d rejected.", fNLocalFits, fNLocalFitsRejected));
if (fIter == 0) break;
fIter++;
if (fIter == fMaxIter) break;
fNLocalFits = 0;
fNLocalFitsRejected = 0;
if (fChi2CutFactor != fChi2CutRef) {
fChi2CutFactor = TMath::Sqrt(fChi2CutFactor);
if (fChi2CutFactor < 1.2*fChi2CutRef) {
fChi2CutFactor = fChi2CutRef;
fIter = fMaxIter - 1;
}
}
AliInfo("");
AliInfo("----------------------------------------------");
AliInfo(Form("Iteration %d with cut factor %.2f", fIter, fChi2CutFactor));
for (i=0; i<nVar; i++) {
fVecBGlo[i] = 0.0;
for (j=0; j<nVar; j++) {
fMatCGlo[i][j] = 0.0;
}
}
nLocFitsGood = 0;
for (i=0; i<nLocFitsTot; i++) {
int iEqFirst = 0;
int iEqLast = 0;
(i>0) ? iEqFirst = fLocEqPlace[i-1] : iEqFirst = 0;
iEqLast = fLocEqPlace[i];
AliDebug(2,Form("Track %d : ",i));
AliDebug(2,Form("Starts at %d", iEqFirst));
AliDebug(2,Form("Ends at %d",iEqLast));
if (fIndexAllEqs[iEqFirst] != -999) {
fIndexLocEq.Reset(); fNIndexLocEq=0;
fDerivLocEq.Reset(); fNDerivLocEq=0;
for (j=iEqFirst; j<iEqLast; j++) {
if (fNIndexLocEq==fIndexLocEq.GetSize()) fIndexLocEq.Set(2*fNIndexLocEq);
fIndexLocEq.AddAt(fIndexAllEqs[j],fNIndexLocEq++);
if (fNDerivLocEq==fDerivLocEq.GetSize()) fDerivLocEq.Set(2*fNDerivLocEq);
fDerivLocEq.AddAt(fDerivAllEqs[j],fNDerivLocEq++);
}
for (j=0; j<2*fNLocalPar; j++) {
localPars[j] = 0.;
}
AliMillepede::LocalFit(i,localPars,0);
nLocFitsGood++;
}
}
AliMillepede::SetNLocalEquations(nLocFitsGood);
}
AliMillepede::PrintGlobalParameters();
for (j=0; j<fNGlobalPar; j++) {
par[j] = fInitPar[j]+fDeltaPar[j];
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]);
error[j] = TMath::Sqrt(TMath::Abs(fMatCGlo[j][j]));
}
AliInfo(" ");
AliInfo(" * o o o ");
AliInfo(" o o o ");
AliInfo(" o ooooo o o o oo ooo oo ooo oo ");
AliInfo(" o o o o o o o o o o o o o o o o ");
AliInfo(" o o o o o o oooo o o oooo o o oooo ");
AliInfo(" o o o o o o o ooo o o o o ");
AliInfo(" o o o o o o oo o oo ooo oo ++ ends.");
AliInfo(" o ");
return 1;
}
Double_t AliMillepede::GetParError(int iPar) const
{
Double_t lErr = -1.;
if (iPar>=0 && iPar<fNGlobalPar) {
lErr = TMath::Sqrt(TMath::Abs(fMatCGlo[iPar][iPar]));
}
return lErr;
}
int AliMillepede::SpmInv(double matV[][fgkMaxGloPC], double vecB[], int nGlo)
{
Int_t nRank = 0;
int iPivot;
double vPivot = 0.;
double eps = 0.00000000000001;
bool *bUnUsed = new bool[nGlo];
double *diagV = new double[nGlo];
double *rowMax = new double[nGlo];
double *colMax = new double[nGlo];
double *temp = new double[nGlo];
for (Int_t i=0; i<nGlo; i++) {
rowMax[i] = 0.0;
colMax[i] = 0.0;
bUnUsed[i] = true;
for (Int_t j=0; j<i; j++) {
if (matV[j][i] == 0) {
matV[j][i] = matV[i][j];
}
}
}
for (Int_t i=0; i<nGlo; i++) {
for (Int_t j=0; j<nGlo; j++) {
if (TMath::Abs(matV[i][j]) >= rowMax[i]) rowMax[i] = TMath::Abs(matV[i][j]);
if (TMath::Abs(matV[j][i]) >= colMax[i]) colMax[i] = TMath::Abs(matV[j][i]);
}
}
for (Int_t i=0; i<nGlo; i++) {
if (0.0 != rowMax[i]) rowMax[i] = 1./rowMax[i];
if (0.0 != colMax[i]) colMax[i] = 1./colMax[i];
}
for (Int_t i=0; i<nGlo; i++) {
for (Int_t j=0; j<nGlo; j++) {
matV[i][j] = TMath::Sqrt(rowMax[i])*matV[i][j]*TMath::Sqrt(colMax[j]);
}
diagV[i] = TMath::Abs(matV[i][i]);
}
for (Int_t i=0; i<nGlo; i++) {
vPivot = 0.0;
iPivot = -1;
for (Int_t j=0; j<nGlo; j++) {
if (bUnUsed[j] && (TMath::Abs(matV[j][j])>TMath::Max(TMath::Abs(vPivot),eps*diagV[j]))) {
vPivot = matV[j][j];
iPivot = j;
}
}
if (iPivot >= 0) {
nRank++;
bUnUsed[iPivot] = false;
vPivot = 1.0/vPivot;
matV[iPivot][iPivot] = -vPivot;
for (Int_t j=0; j<nGlo; j++) {
for (Int_t jj=0; jj<nGlo; jj++) {
if (j != iPivot && jj != iPivot) {
matV[j][jj] = matV[j][jj] - vPivot*matV[j][iPivot]*matV[iPivot][jj];
}
}
}
for (Int_t j=0; j<nGlo; j++) {
if (j != iPivot) {
matV[j][iPivot] = matV[j][iPivot]*vPivot;
matV[iPivot][j] = matV[iPivot][j]*vPivot;
}
}
}
else {
for (Int_t j=0; j<nGlo; j++) {
if (bUnUsed[j]) {
vecB[j] = 0.0;
for (Int_t k=0; k<nGlo; k++) {
matV[j][k] = 0.0;
matV[k][j] = 0.0;
}
}
}
break;
}
}
for (Int_t i=0; i<nGlo; i++) {
for (Int_t j=0; j<nGlo; j++) {
matV[i][j] = TMath::Sqrt(colMax[i])*matV[i][j]*TMath::Sqrt(rowMax[j]);
}
}
for (Int_t j=0; j<nGlo; j++) {
temp[j] = 0.0;
for (Int_t jj=0; jj<nGlo; jj++) {
matV[j][jj] = -matV[j][jj];
temp[j] += matV[j][jj]*vecB[jj];
}
}
for (Int_t j=0; j<nGlo; j++) {
vecB[j] = temp[j];
}
delete [] temp;
delete [] bUnUsed;
delete [] diagV;
delete [] rowMax;
delete [] colMax;
return nRank;
}
int AliMillepede::SpmInv(double matV[][fgkMaxLocalPar], double vecB[], int nLoc)
{
Int_t nRank = 0;
Int_t iPivot = -1;
double vPivot = 0.;
double eps = 0.0000000000001;
bool *bUnUsed = new bool[nLoc];
double *diagV = new double[nLoc];
double *temp = new double[nLoc];
for (Int_t i=0; i<nLoc; i++) {
bUnUsed[i] = true;
diagV[i] = TMath::Abs(matV[i][i]);
for (Int_t j=0; j<i; j++) {
matV[j][i] = matV[i][j] ;
}
}
for (Int_t i=0; i<nLoc; i++) {
vPivot = 0.0;
iPivot = -1;
for (Int_t j=0; j<nLoc; j++) {
if (bUnUsed[j] && (TMath::Abs(matV[j][j])>TMath::Max(TMath::Abs(vPivot),eps*diagV[j]))) {
vPivot = matV[j][j];
iPivot = j;
}
}
if (iPivot >= 0) {
nRank++;
bUnUsed[iPivot] = false;
vPivot = 1.0/vPivot;
matV[iPivot][iPivot] = -vPivot;
for (Int_t j=0; j<nLoc; j++) {
if (j != iPivot) {
for (Int_t jj=0; jj<=j; jj++) {
if (jj != iPivot) {
matV[j][jj] = matV[j][jj] - vPivot*matV[j][iPivot]*matV[iPivot][jj];
matV[jj][j] = matV[j][jj];
}
}
}
}
for (Int_t j=0; j<nLoc; j++) {
if (j != iPivot) {
matV[j][iPivot] = matV[j][iPivot]*vPivot;
matV[iPivot][j] = matV[j][iPivot];
}
}
}
else {
for (Int_t j=0; j<nLoc; j++) {
if (bUnUsed[j]) {
vecB[j] = 0.0;
matV[j][j] = 0.0;
for (Int_t k=0; k<j; k++) {
matV[j][k] = 0.0;
matV[k][j] = 0.0;
}
}
}
break;
}
}
for (Int_t j=0; j<nLoc; j++) {
temp[j] = 0.0;
for (Int_t jj=0; jj<nLoc; jj++) {
matV[j][jj] = -matV[j][jj];
temp[j] += matV[j][jj]*vecB[jj];
}
}
for (Int_t j=0; j<nLoc; j++) {
vecB[j] = temp[j];
}
delete [] bUnUsed;
delete [] diagV;
delete [] temp;
return nRank;
}
Int_t AliMillepede::SpAVAt(double matV[][fgkMaxLocalPar], double matA[][fgkMaxLocalPar], double matW[][fgkMaxGlobalPar], int nLoc, int nGlo)
{
for (Int_t i=0; i<nGlo; i++) {
for (Int_t j=0; j<=i; j++) {
matW[i][j] = 0.0;
for (Int_t k=0; k<nLoc; k++) {
matW[i][j] += matA[i][k]*matV[k][k]*matA[j][k];
for (Int_t l=0; l<k; l++) {
matW[i][j] += matA[i][k]*matV[k][l]*matA[j][l];
matW[i][j] += matA[i][l]*matV[k][l]*matA[j][k];
}
}
if (i!=j){
matW[j][i] = matW[i][j];
}
}
}
return 1;
}
Int_t AliMillepede::SpAX(double matA[][fgkMaxLocalPar], double vecX[], double vecY[], int nCol, int nRow)
{
for (Int_t i=0; i<nRow; i++) {
vecY[i] = 0.0;
for (Int_t j=0; j<nCol; j++) {
vecY[i] += matA[i][j]*vecX[j];
}
}
return 1;
}
Int_t AliMillepede::PrintGlobalParameters() const
{
double lError = 0.;
double lGlobalCor =0.;
AliInfo("");
AliInfo(" Result of fit for global parameters");
AliInfo(" ===================================");
AliInfo(" I initial final differ lastcor error gcor");
AliInfo("-----------------------------------------------------------------------------------");
for (int i=0; i<fNGlobalPar; i++) {
lError = TMath::Sqrt(TMath::Abs(fMatCGlo[i][i]));
if (fMatCGlo[i][i] < 0.0) lError = -lError;
lGlobalCor = 0.0;
if (TMath::Abs(fMatCGlo[i][i]*fDiagCGlo[i]) > 0) {
lGlobalCor = TMath::Sqrt(TMath::Abs(1.0-1.0/(fMatCGlo[i][i]*fDiagCGlo[i])));
AliInfo(Form("%d\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f",
i,fInitPar[i],fInitPar[i]+fDeltaPar[i],fDeltaPar[i],fVecBGlo[i],lError,lGlobalCor));
}
else {
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]));
}
}
return 1;
}
double AliMillepede::Chi2DoFLim(int nSig, int nDoF)
{
int lNSig;
double sn[3] = {0.47523, 1.690140, 2.782170};
double table[3][30] = {{1.0000, 1.1479, 1.1753, 1.1798, 1.1775, 1.1730, 1.1680, 1.1630,
1.1581, 1.1536, 1.1493, 1.1454, 1.1417, 1.1383, 1.1351, 1.1321,
1.1293, 1.1266, 1.1242, 1.1218, 1.1196, 1.1175, 1.1155, 1.1136,
1.1119, 1.1101, 1.1085, 1.1070, 1.1055, 1.1040},
{4.0000, 3.0900, 2.6750, 2.4290, 2.2628, 2.1415, 2.0481, 1.9736,
1.9124, 1.8610, 1.8171, 1.7791, 1.7457, 1.7161, 1.6897, 1.6658,
1.6442, 1.6246, 1.6065, 1.5899, 1.5745, 1.5603, 1.5470, 1.5346,
1.5230, 1.5120, 1.5017, 1.4920, 1.4829, 1.4742},
{9.0000, 5.9146, 4.7184, 4.0628, 3.6410, 3.3436, 3.1209, 2.9468,
2.8063, 2.6902, 2.5922, 2.5082, 2.4352, 2.3711, 2.3143, 2.2635,
2.2178, 2.1764, 2.1386, 2.1040, 2.0722, 2.0428, 2.0155, 1.9901,
1.9665, 1.9443, 1.9235, 1.9040, 1.8855, 1.8681}};
if (nDoF < 1) {
return 0.0;
}
else {
lNSig = TMath::Max(1,TMath::Min(nSig,3));
if (nDoF <= 30) {
return table[lNSig-1][nDoF-1];
}
else {
return ((sn[lNSig-1]+TMath::Sqrt(float(2*nDoF-3)))*
(sn[lNSig-1]+TMath::Sqrt(float(2*nDoF-3))))/float(2*nDoF-2);
}
}
}