X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliMillePede2.cxx;h=4963e0d17023c26a3c0c170899cbd6568e205888;hb=d9719df10e380ba1c38bd18363f951d63d3e1d2a;hp=aca79e251d6af6f0898082428bc23c3e495c7537;hpb=f748efd7abe787c78622a18a13c9a18db6f2f9ed;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliMillePede2.cxx b/STEER/AliMillePede2.cxx index aca79e251d6..4963e0d1702 100644 --- a/STEER/AliMillePede2.cxx +++ b/STEER/AliMillePede2.cxx @@ -1,102 +1,37 @@ +/**********************************************************************************************/ +/* General class for alignment with large number of degrees of freedom */ +/* Based on the original milliped2 by Volker Blobel */ +/* http://www.desy.de/~blobel/mptalks.html */ +/* */ +/* Author: ruben.shahoyan@cern.ch */ +/* */ +/**********************************************************************************************/ + #include "AliMillePede2.h" #include "AliLog.h" #include +#include +#include +#include +#include "AliMatrixSq.h" +#include "AliSymMatrix.h" +#include "AliRectMatrix.h" +#include "AliMatrixSparse.h" -using namespace std; - -ClassImp(AliMillePedeRecord) -//_____________________________________________________________________________________________ -AliMillePedeRecord::AliMillePedeRecord() : -fSize(0),fIndex(0),fValue(0) {SetBufferSize(0);} - -//_____________________________________________________________________________________________ -AliMillePedeRecord::AliMillePedeRecord(const AliMillePedeRecord& src) : - TObject(src),fSize(src.fSize),fIndex(0),fValue(0) -{ - fIndex = new int[GetBufferSize()]; - memcpy(fIndex,src.fIndex,fSize*sizeof(int)); - fValue = new double[GetBufferSize()]; - memcpy(fValue,src.fValue,fSize*sizeof(double)); -} - -//_____________________________________________________________________________________________ -AliMillePedeRecord& AliMillePedeRecord::operator=(const AliMillePedeRecord& rhs) -{ - if (this!=&rhs) { - Reset(); - for (int i=0;iSetSymmetric(kTRUE);} else fMatCGlo = new AliSymMatrix(fNGloPar); // + fFillIndex = new Int_t[fNGloPar]; + fFillValue = new Double_t[fNGloPar]; + // fMatCLoc = new AliSymMatrix(fNLocPar); - fMatCGloLoc = new TMatrixD(fNGloPar,fNLocPar); + fMatCGloLoc = new AliRectMatrix(fNGloPar,fNLocPar); // + fParamGrID = new Int_t[fNGloPar]; fProcPnt = new Int_t[fNGloPar]; fVecBLoc = new Double_t[fNLocPar]; fDiagCGlo = new Double_t[fNGloPar]; @@ -229,23 +182,14 @@ Int_t AliMillePede2::InitMille(int nGlo, int nLoc, int lNStdDev,double lResCut, memset(fProcPnt ,0,fNGloPar*sizeof(Int_t)); // for (int i=fNGloPar;i--;) { - fGlo2CGlo[i]=fCGlo2Glo[i] = -1; + fGlo2CGlo[i] = fCGlo2Glo[i] = -1; fIsLinear[i] = kTRUE; + fParamGrID[i] = -1; } // return 1; } -//_____________________________________________________________________________________________ -void AliMillePede2::ResetConstraints() -{ - // reset constraints tree. Allows to redefine the constraints for preprocessed data - CloseConsRecStorage(); - InitConsRecStorage(kFALSE); - fNGloConstraints = 0; - AliInfo("Constraints are reset"); -} - //_____________________________________________________________________________________________ Bool_t AliMillePede2::ImposeDataRecFile(const char* fname) { @@ -272,12 +216,12 @@ Bool_t AliMillePede2::InitDataRecStorage(Bool_t read) if (!fRecord) fRecord = new AliMillePedeRecord(); // fDataRecFile = TFile::Open(GetDataRecFName(),read ? "":"recreate"); - if (!fDataRecFile) {AliInfo(Form("Failed to initialize data records file %s",GetDataRecFName())); return kFALSE;} + if (!fDataRecFile) {AliFatal(Form("Failed to initialize data records file %s",GetDataRecFName())); return kFALSE;} // AliInfo(Form("File %s used for derivatives records",GetDataRecFName())); if (read) { fTreeData = (TTree*)fDataRecFile->Get("AliMillePedeRecords_Data"); - if (!fTreeData) {AliInfo(Form("Did not find data records tree in %s",GetDataRecFName())); return kFALSE;} + if (!fTreeData) {AliFatal(Form("Did not find data records tree in %s",GetDataRecFName())); return kFALSE;} fTreeData->SetBranchAddress("Record_Data",&fRecord); AliInfo(Form("Found %d derivatives records",fTreeData->GetEntries())); } @@ -286,6 +230,7 @@ Bool_t AliMillePede2::InitDataRecStorage(Bool_t read) fTreeData->Branch("Record_Data","AliMillePedeRecord",&fRecord,32000,99); } fCurrRecDataID = -1; + fRecFileStatus = read ? 1:2; // return kTRUE; } @@ -334,6 +279,7 @@ void AliMillePede2::CloseDataRecStorage() delete fDataRecFile; fDataRecFile = 0; } + fRecFileStatus = 0; // } @@ -372,10 +318,17 @@ Bool_t AliMillePede2::ReadNextRecordConstraint() return kTRUE; } +//_____________________________________________________________________________________________ +void AliMillePede2::SetRecordWeight(double wgh) +{ + if (fRecFileStatus<2) InitDataRecStorage(); // create a buffer to store the data + fRecord->SetWeight(wgh); +} + //_____________________________________________________________________________________________ void AliMillePede2::SetLocalEquation(double *dergb, double *derlc, double lMeas, double lSigma) { - if (!fDataRecFile || !fDataRecFile->IsWritable()) InitDataRecStorage(); // create a buffer to store the data + if (fRecFileStatus<2) InitDataRecStorage(); // create a buffer to store the data // // write data of single measurement if (lSigma<=0.0) { // If parameter is fixed, then no equation @@ -387,12 +340,15 @@ void AliMillePede2::SetLocalEquation(double *dergb, double *derlc, double lMeas, fRecord->AddResidual(lMeas); // // Retrieve local param interesting indices - for (int i=0;iAddIndexValue(i,derlc[i]); derlc[i] = 0.0;} + for (int i=0;iAddIndexValue(i,derlc[i]); derlc[i] = 0.0;} // fRecord->AddWeight( 1.0/lSigma/lSigma ); // // Idem for global parameters - for (int i=0;iAddIndexValue(i,dergb[i]); dergb[i] = 0.0;} + for (int i=0;iAddIndexValue(i,dergb[i]); dergb[i] = 0.0; + fRecord->MarkGroup(fParamGrID[i]); + } // } @@ -407,45 +363,49 @@ void AliMillePede2::SetLocalEquation(int *indgb, double *dergb, int ngb, int *in return; } // - if (!fDataRecFile || !fDataRecFile->IsWritable()) InitDataRecStorage(); // create a buffer to store the data + if (fRecFileStatus<2) InitDataRecStorage(); // create a buffer to store the data // fRecord->AddResidual(lMeas); // // Retrieve local param interesting indices - for (int i=0;iAddIndexValue(indlc[i],derlc[i]); derlc[i]=0.; indlc[i]=0;} + for (int i=0;iAddIndexValue(indlc[i],derlc[i]); derlc[i]=0.; indlc[i]=0;} // fRecord->AddWeight( 1./lSigma/lSigma ); // // Idem for global parameters - for (int i=0;iAddIndexValue(indgb[i],dergb[i]); dergb[i]=0.; indgb[i]=0;} + for (int i=0;iAddIndexValue(indgb[i],dergb[i]); dergb[i]=0.; indgb[i]=0;} // } + //_____________________________________________________________________________________________ -void AliMillePede2::SetGlobalConstraint(double *dergb, double val) +void AliMillePede2::SetGlobalConstraint(double *dergb, double val, double sigma) { // Define a constraint equation. if (!fConsRecFile || !fConsRecFile->IsWritable()) InitConsRecStorage(); // create a buffer to store the data // fRecord->Reset(); fRecord->AddResidual(val); - fRecord->AddWeight(val); // dummy - for (int i=0; iAddIndexValue(i,dergb[i]); + fRecord->AddWeight(sigma); + for (int i=0; iAddIndexValue(i,dergb[i]); fNGloConstraints++; + if (IsZero(sigma)) fNLagrangeConstraints++; + // printf("NewConstraint:\n"); fRecord->Print(); //RRR SaveRecordConstraint(); // } //_____________________________________________________________________________________________ -void AliMillePede2::SetGlobalConstraint(const int *indgb, double *dergb, int ngb, double val) +void AliMillePede2::SetGlobalConstraint(const int *indgb, double *dergb, int ngb, double val,double sigma) { // Define a constraint equation. if (!fConsRecFile || !fConsRecFile->IsWritable()) InitConsRecStorage(); // create a buffer to store the data fRecord->Reset(); fRecord->AddResidual(val); - fRecord->AddWeight(val); // dummy - for (int i=0; iAddIndexValue(indgb[i],dergb[i]); + fRecord->AddWeight(sigma); // dummy + for (int i=0; iAddIndexValue(indgb[i],dergb[i]); fNGloConstraints++; + if (IsZero(sigma)) fNLagrangeConstraints++; SaveRecordConstraint(); // } @@ -459,16 +419,17 @@ Int_t AliMillePede2::LocalFit(double *localParams) localParams = (if !=0) will contain the fitted track parameters and related errors */ - static int nRefSize = 0; - static TArrayI RefLoc,RefGlo,nRefLoc,nRefGlo; + static int nrefSize = 0; + // static TArrayI refLoc,refGlo,nrefLoc,nrefGlo; + static Int_t *refLoc=0,*refGlo=0,*nrefLoc=0,*nrefGlo=0; int nPoints = 0; // - AliSymMatrix &MatCLoc = *fMatCLoc; - AliMatrixSq &MatCGlo = *fMatCGlo; - TMatrixD &MatCGloLoc = *fMatCGloLoc; + AliSymMatrix &matCLoc = *fMatCLoc; + AliMatrixSq &matCGlo = *fMatCGlo; + AliRectMatrix &matCGloLoc = *fMatCGloLoc; // memset(fVecBLoc,0,fNLocPar*sizeof(double)); - MatCLoc.Reset(); + matCLoc.Reset(); // int cnt=0; int recSz = fRecord->GetSize(); @@ -476,86 +437,90 @@ Int_t AliMillePede2::LocalFit(double *localParams) while(cntIsWeight(cnt)) {nLoc++; cnt++;} - nRefLoc[nPoints] = nLoc; + nrefLoc[nPoints] = nLoc; // - RefGlo[nPoints] = ++cnt; + refGlo[nPoints] = ++cnt; int nGlo = 0; while(!fRecord->IsResidual(cnt) && cntGetWeight(); // global weight for this set + 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 ); - double *derLoc = fRecord->GetValue()+RefLoc[ip]; - double *derGlo = fRecord->GetValue()+RefGlo[ip]; - int *indLoc = fRecord->GetIndex()+RefLoc[ip]; - int *indGlo = fRecord->GetIndex()+RefGlo[ip]; - // - for (int i=nRefGlo[ip];i--;) { // suppress the global part (only relevant with iterations) + double resid = fRecord->GetValue( refLoc[ip]-1 ); + double weight = fRecord->GetValue( refGlo[ip]-1 )*gloWgh; + double *derLoc = fRecord->GetValue()+refLoc[ip]; + double *derGlo = fRecord->GetValue()+refGlo[ip]; + int *indLoc = fRecord->GetIndex()+refLoc[ip]; + int *indGlo = fRecord->GetIndex()+refGlo[ip]; + // + for (int i=nrefGlo[ip];i--;) { // suppress the global part (only relevant with iterations) int iID = indGlo[i]; // Global param indice - if ( fSigmaPar[iID] <= 0.) continue; // fixed parameter RRRCheck + if (fSigmaPar[iID]<=0.) continue; // fixed parameter RRRCheck if (fIsLinear[iID]) resid -= derGlo[i]*(fInitPar[iID]+fDeltaPar[iID]); // linear parameter else resid -= derGlo[i]*fDeltaPar[iID]; // nonlinear parameter } // - for (int i=nRefLoc[ip];i--;) { // Fill local matrix and vector - int iID = indLoc[i]; // Local param indice (the matrix line) - if ( (vl=weight*resid*derLoc[i])!=0 ) fVecBLoc[iID] += vl; - // - for (int j=i+1;j--;) { // Symmetric matrix, don't bother j>i coeffs - int jID = indLoc[j]; - if ( (vl=weight*derLoc[i]*derLoc[j])!=0 ) MatCLoc(iID,jID) += vl; - } - } + // 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 (!MatCLoc.SolveSpmInv(fVecBLoc,kTRUE)) { + if (!matCLoc.SolveChol(fVecBLoc,kTRUE)) { + AliInfo("Failed to solve locals by Cholesky, trying Gaussian Elimination"); + if (!matCLoc.SolveSpmInv(fVecBLoc,kTRUE)) { AliInfo("Failed to solve locals by Gaussian Elimination, skip..."); - MatCLoc.Print("d"); + matCLoc.Print("d"); return 0; // failed to solve } } // // If requested, store the track params and errors - if (localParams) for (int i=fNLocPar; i--;) { + // printf("locfit: "); for (int i=0;iGetValue( RefLoc[ip]-1 ); - double weight = fRecord->GetValue( RefGlo[ip]-1 ); - double *derLoc = fRecord->GetValue()+RefLoc[ip]; - double *derGlo = fRecord->GetValue()+RefGlo[ip]; - int *indLoc = fRecord->GetIndex()+RefLoc[ip]; - int *indGlo = fRecord->GetIndex()+RefGlo[ip]; + double resid = fRecord->GetValue( refLoc[ip]-1 ); + double weight = fRecord->GetValue( refGlo[ip]-1 )*gloWgh; + double *derLoc = fRecord->GetValue()+refLoc[ip]; + double *derGlo = fRecord->GetValue()+refGlo[ip]; + int *indLoc = fRecord->GetIndex()+refLoc[ip]; + int *indGlo = fRecord->GetIndex()+refGlo[ip]; // // Suppress local and global contribution in residuals; - for (int i=nRefLoc[ip];i--;) resid -= derLoc[i]*fVecBLoc[ indLoc[i] ]; // local part + for (int i=nrefLoc[ip];i--;) resid -= derLoc[i]*fVecBLoc[ indLoc[i] ]; // local part // - for (int i=nRefGlo[ip];i--;) { // global part + for (int i=nrefGlo[ip];i--;) { // global part int iID = indGlo[i]; if ( fSigmaPar[iID] <= 0.) continue; // fixed parameter RRRCheck if (fIsLinear[iID]) resid -= derGlo[i]*(fInitPar[iID]+fDeltaPar[iID]); // linear parameter @@ -563,9 +528,11 @@ Int_t AliMillePede2::LocalFit(double *localParams) } // // reject the track if the residual is too large (outlier) - if ( (TMath::Abs(resid) >= fResCutInit && fIter ==1 ) || - (TMath::Abs(resid) >= fResCut && fIter > 1)) { - fNLocFitsRejected++; + double absres = TMath::Abs(resid); + if ( (absres >= fResCutInit && fIter ==1 ) || + (absres >= fResCut && fIter > 1)) { + if (fLocFitAdd) fNLocFitsRejected++; + // printf("reject res %5ld %+e\n",fCurrRecDataID,resid); return 0; } // @@ -573,17 +540,24 @@ Int_t AliMillePede2::LocalFit(double *localParams) nEq++; // number of equations } // end of Calculate residuals // - // - int nDoF = nEq-fNLocPar; + lChi2 /= gloWgh; + 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 - fNLocFitsRejected++; + if (fLocFitAdd) fNLocFitsRejected++; + // printf("reject chi2 %5ld: %+e\n",fCurrRecDataID, lChi2); return 0; } // - fNLocFits++; - fNLocEquations += nEq; + if (fLocFitAdd) { + fNLocFits++; + fNLocEquations += nEq; + } + else { + fNLocFits--; + fNLocEquations -= nEq; + } // // local operations are finished, track is accepted // We now update the global parameters (other matrices) @@ -591,46 +565,50 @@ Int_t AliMillePede2::LocalFit(double *localParams) int nGloInFit = 0; // for (int ip=nPoints;ip--;) { // Update matrices - double resid = fRecord->GetValue( RefLoc[ip]-1 ); - double weight = fRecord->GetValue( RefGlo[ip]-1 ); - double *derLoc = fRecord->GetValue()+RefLoc[ip]; - double *derGlo = fRecord->GetValue()+RefGlo[ip]; - int *indLoc = fRecord->GetIndex()+RefLoc[ip]; - int *indGlo = fRecord->GetIndex()+RefGlo[ip]; - // - for (int i=nRefGlo[ip];i--;) { // suppress the global part + double resid = fRecord->GetValue( refLoc[ip]-1 ); + double weight = fRecord->GetValue( refGlo[ip]-1 )*gloWgh; + double *derLoc = fRecord->GetValue()+refLoc[ip]; + double *derGlo = fRecord->GetValue()+refGlo[ip]; + int *indLoc = fRecord->GetIndex()+refLoc[ip]; + int *indGlo = fRecord->GetIndex()+refGlo[ip]; + // + for (int i=nrefGlo[ip];i--;) { // suppress the global part int iID = indGlo[i]; // Global param indice if ( fSigmaPar[iID] <= 0.) continue; // fixed parameter RRRCheck - if (fIsLinear[iID] == 0) resid -= derGlo[i]*(fInitPar[iID]+fDeltaPar[iID]); // linear parameter - else resid -= derGlo[i]*fDeltaPar[iID]; // nonlinear parameter + if (fIsLinear[iID]) resid -= derGlo[i]*(fInitPar[iID]+fDeltaPar[iID]); // linear parameter + else resid -= derGlo[i]*fDeltaPar[iID]; // nonlinear parameter } // - for (int ig=nRefGlo[ip];ig--;) { + for (int ig=nrefGlo[ip];ig--;) { int iIDg = indGlo[ig]; // Global param indice (the matrix line) if ( fSigmaPar[iIDg] <= 0.) continue; // fixed parameter RRRCheck - if ( (vl = weight*resid*derGlo[ig])!=0 ) fVecBGlo[ iIDg ] += vl; + if (fLocFitAdd) fVecBGlo[ iIDg ] += weight*resid*derGlo[ig]; //!!! + else fVecBGlo[ iIDg ] -= weight*resid*derGlo[ig]; //!!! // // First of all, the global/global terms (exactly like local matrix) - for (int jg=ig+1;jg--;) { // MatCGlo is symmetric by construction + int nfill = 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 ) MatCGlo(iIDg,jIDg) += vl; - } + if ( !IsZero(vl = weight*derGlo[ig]*derGlo[jg]) ) { + fFillIndex[nfill] = jIDg; + fFillValue[nfill++] = fLocFitAdd ? vl:-vl; + } + } + if (nfill) matCGlo.AddToRow(iIDg,fFillValue,fFillIndex,nfill); // // Now we have also rectangular matrices containing global/local terms. int iCIDg = fGlo2CGlo[iIDg]; // compressed Index of index if (iCIDg == -1) { - for (int k=fNLocPar;k--;) MatCGloLoc(nGloInFit,k) = 0.0; // reset the row + Double_t *rowGL = matCGloLoc(nGloInFit); + for (int k=maxLocUsed;k--;) rowGL[k] = 0.0; // reset the row iCIDg = fGlo2CGlo[iIDg] = nGloInFit; fCGlo2Glo[nGloInFit++] = iIDg; } // - for (int il=nRefLoc[ip];il--;) { - int iIDl = indLoc[il]; - if ( (vl = weight*derGlo[ig]*derLoc[il])!=0) MatCGloLoc(iCIDg,iIDl) += vl; - } - // update counter - fProcPnt[iIDg]++; + Double_t *rowGLIDg = matCGloLoc(iCIDg); + for (int il=nrefLoc[ip];il--;) rowGLIDg[ indLoc[il] ] += weight*derGlo[ig]*derLoc[il]; + fProcPnt[iIDg] += fLocFitAdd ? 1:-1; // update counter // } } // end of Update matrices @@ -638,37 +616,48 @@ Int_t AliMillePede2::LocalFit(double *localParams) // calculate fMatCGlo -= fMatCGloLoc * fMatCLoc * fMatCGloLoc^T // and fVecBGlo -= fMatCGloLoc * fVecBLoc // + //-------------------------------------------------------------- >>> double vll; for (int iCIDg=0; iCIDg0 ? TMath::Sqrt(TMath::Abs(fMatCGlo->QuerryDiag(i))) : 0.; - if (pull) for (int i=fNGloPar;i--;) pull[i] = fProcPnt[i]>0 && (fSigmaPar[i]*fSigmaPar[i]-fMatCGlo->QuerryDiag(i))>0. && fSigmaPar[i]>0 - ? fDeltaPar[i]/TMath::Sqrt(fSigmaPar[i]*fSigmaPar[i]-fMatCGlo->QuerryDiag(i)) : 0; + if (error) for (int i=fNGloPar;i--;) error[i] = fProcPnt[i]>0 ? TMath::Sqrt(TMath::Abs(fMatCGlo->QueryDiag(i))) : 0.; + if (pull) for (int i=fNGloPar;i--;) pull[i] = fProcPnt[i]>0 && (fSigmaPar[i]*fSigmaPar[i]-fMatCGlo->QueryDiag(i))>0. && fSigmaPar[i]>0 + ? fDeltaPar[i]/TMath::Sqrt(fSigmaPar[i]*fSigmaPar[i]-fMatCGlo->QueryDiag(i)) : 0; } // return 1; @@ -723,103 +712,215 @@ Int_t AliMillePede2::GlobalFitIteration() AliInfo("No data was stored, aborting iteration"); return 0; } - TStopwatch sw; sw.Start(); + TStopwatch sw,sws; + sw.Start(); + sws.Stop(); // if (!fConstrUsed) { fConstrUsed = new Bool_t[fNGloConstraints]; memset(fConstrUsed,0,fNGloConstraints*sizeof(Bool_t)); } // Reset all info specific for this step - AliMatrixSq& MatCGlo = *fMatCGlo; - MatCGlo.Reset(); + AliMatrixSq& matCGlo = *fMatCGlo; + matCGlo.Reset(); memset(fProcPnt,0,fNGloPar*sizeof(Int_t)); // fNGloConstraints = fTreeConstr ? fTreeConstr->GetEntries() : 0; // + // count number of Lagrange constraints: they need new row/cols to be added + fNLagrangeConstraints = 0; + for (int i=0; iGetValue(1)) ) fNLagrangeConstraints++; // exact constraint (no error) -> Lagrange multiplier + } + // // if needed, readjust the size of the global vector (for matrices this is done automatically) - if (!fVecBGlo || fNGloSize!=fNGloPar+fNGloConstraints) { + if (!fVecBGlo || fNGloSize!=fNGloPar+fNLagrangeConstraints) { delete[] fVecBGlo; // in case some constraint was added between the two manual iterations - fNGloSize = fNGloPar+fNGloConstraints; + fNGloSize = fNGloPar+fNLagrangeConstraints; fVecBGlo = new Double_t[fNGloSize]; } memset(fVecBGlo,0,fNGloSize*sizeof(double)); // fNLocFits = 0; fNLocFitsRejected = 0; - fNLocEquations = 0; + fNLocEquations = 0; // // Process data records and build the matrices Long_t ndr = fTreeData->GetEntries(); AliInfo(Form("Building the Global matrix from %d data records",ndr)); + if (!ndr) return 0; + // + TStopwatch swt; swt.Start(); + fLocFitAdd = kTRUE; // add contributions of matching tracks for (Long_t i=0;i> fNGloFix = 0; - for (int i=fNGloPar;i--;) fDiagCGlo[i] = MatCGlo.QuerryDiag(i); // save the diagonal elements + if (fMinPntValid>1 && fNGroupsSet) { + // + printf("Checking parameters with statistics < %d\n",fMinPntValid); + TStopwatch swsup; + swsup.Start(); + // 1) build the list of parameters to fix + Int_t fixArrSize = 10; + Int_t nFixedGroups = 0; + TArrayI fixGroups(fixArrSize); + // + int grIDold = -2; + int oldStart = -1; + double oldMin = 1.e20; + double oldMax =-1.e20; + // + for (int i=fNGloPar;i--;) { // // Reset row and column of fixed params and add 1/sig^2 to free ones + int grID = fParamGrID[i]; + if (grID<0) continue; // not in the group + // + if (grID!=grIDold) { // starting new group + if (grIDold>=0) { // decide if the group has enough statistics + if (oldMini;iold--) fProcPnt[iold] = 0; + Bool_t fnd = kFALSE; // check if the group is already accounted + for (int j=nFixedGroups;j--;) if (fixGroups[j]==grIDold) {fnd=kTRUE; break;} + if (!fnd) { + if (nFixedGroups>=fixArrSize) {fixArrSize*=2; fixGroups.Set(fixArrSize);} + fixGroups[nFixedGroups++] = grIDold; // add group to fix + } + } + } + grIDold = grID; // mark the start of the new group + oldStart = i; + oldMin = 1.e20; + oldMax = -1.e20; + } + if (oldMin>fProcPnt[i]) oldMin = fProcPnt[i]; + if (oldMax=0 && oldMin=fixArrSize) {fixArrSize*=2; fixGroups.Set(fixArrSize);} + fixGroups[nFixedGroups++] = grIDold; // add group to fix + } + } + // + // 2) loop over records and add contributions of fixed groups with negative sign + fLocFitAdd = kFALSE; + // + for (Long_t i=0;iIsGroupPresent(fixGroups[ifx])) suppr = kTRUE; + if (suppr) LocalFit(); + } + fLocFitAdd = kTRUE; + // + if (nFixedGroups) { + printf("Suppressed contributions of groups with NPoints<%d :\n",fMinPntValid); + for (int i=0;iGetRow(i); - row.Clear(); - row(i) = float(fNLocEquations); - for (int j=i+1;jSetToZero(i,j); - } - else - for (int j=fNGloPar;j--;) if (MatCGlo.Querry(i,j)) MatCGlo(i,j)=0; - MatCGlo.DiagElem(i) = float(fNLocEquations); + matCGlo.DiagElem(i) = 1.;//float(fNLocEquations*fNLocEquations); } - else MatCGlo.DiagElem(i) += 1.0/(fSigmaPar[i]*fSigmaPar[i]); + else matCGlo.DiagElem(i) += (fgWeightSigma ? fProcPnt[i] : 1.)/(fSigmaPar[i]*fSigmaPar[i]); } // + for (int i=fNGloPar;i--;) fDiagCGlo[i] = matCGlo.QueryDiag(i); // save the diagonal elements + // // add constraint equations int nVar = fNGloPar; // Current size of global matrix for (int i=0; iGetValue(0); + double sig = fRecord->GetValue(1); int *indV = fRecord->GetIndex()+2; double *der = fRecord->GetValue()+2; int csize = fRecord->GetSize()-2; // - // check if after suppression of fixed variables this there are non-0 derivatives - int NSuppressed = 0; - for (int j=csize;j--;) if (fProcPnt[indV[j]]<1) NSuppressed++; + // check if after suppression of fixed variables there are non-0 derivatives + // and determine the max statistics of involved params + int nSuppressed = 0; + int maxStat = 1; + for (int j=csize;j--;) { + if (fProcPnt[indV[j]]<1) nSuppressed++; + else { + maxStat = TMath::Max(maxStat,fProcPnt[indV[j]]); + } + } // - if (NSuppressed==csize) { - AliInfo(Form("Neglecting constraint %d of %d derivatives since no free parameters left",i,csize)); + if (nSuppressed==csize) { + // AliInfo(Form("Neglecting constraint %d of %d derivatives since no free parameters left",i,csize)); // // was this constraint ever created ? - if ( fConstrUsed[i] ) { + if ( sig==0 && fConstrUsed[i] ) { // this is needed only for constraints with Lagrange multiplier // to avoid empty row impose dummy constraint on "Lagrange multiplier" - MatCGlo.DiagElem(nVar) = 1.; + matCGlo.DiagElem(nVar) = 1.; fVecBGlo[nVar++] = 0; } continue; } // - for (int j=csize; j--;) { - int jID = indV[j]; - val -= der[j]*(fInitPar[jID]+fDeltaPar[jID]); - if (fProcPnt[jID]<1) continue; // this parameter was fixed, don't put it into constraint - MatCGlo(nVar,jID) = float(fNLocEquations)*der[j]; // fMatCGlo is symmetric, only lower triangle is filled - } + // account for already accumulated corrections + for (int j=csize; j--;) val -= der[j]*(fInitPar[ indV[j] ]+fDeltaPar[ indV[j] ]); // - if (MatCGlo.QuerryDiag(nVar)) MatCGlo.DiagElem(nVar) = 0.0; - fVecBGlo[nVar++] = float(fNLocEquations)*val; //RS ? should we use here fNLocFits ? - fConstrUsed[i] = kTRUE; + if (sig > 0) { // this is a gaussian constriant: no Lagrange multipliers are added + // + double sig2i = (fgWeightSigma ? TMath::Sqrt(maxStat) : 1.)/sig/sig; + for (int ir=0;irPrint(); + printf("RHS\n"); + for (int i=0;iSolveChol(fVecBGlo, kTRUE) ) return gkInvert; + if (fNLagrangeConstraints==0) { // pos-def systems are faster to solve by Cholesky + if ( ((AliSymMatrix*)fMatCGlo)->SolveChol(fVecBGlo, fgInvChol) ) return fgInvChol ? gkInvert:gkNoInversion; else AliInfo("Solution of Global Dense System by Cholesky failed, trying Gaussian Elimiation"); } // @@ -847,22 +956,22 @@ Int_t AliMillePede2::SolveGlobalMatEq() else AliInfo("Solution of Global Dense System by Gaussian Elimiation failed, trying iterative methods"); } // try to solve by minres - TVectorD SOL(fNGloSize); + TVectorD sol(fNGloSize); // AliMinResSolve *slv = new AliMinResSolve(fMatCGlo,fVecBGlo); if (!slv) return gkFailed; // Bool_t res = kFALSE; if (fgIterSol == AliMinResSolve::kSolMinRes) - res = slv->SolveMinRes(SOL,fgMinResCondType,fgMinResMaxIter,fgMinResTol); + res = slv->SolveMinRes(sol,fgMinResCondType,fgMinResMaxIter,fgMinResTol); else if (fgIterSol == AliMinResSolve::kSolFGMRes) - res = slv->SolveFGMRES(SOL,fgMinResCondType,fgMinResMaxIter,fgMinResTol,fgNKrylovV); + res = slv->SolveFGMRES(sol,fgMinResCondType,fgMinResMaxIter,fgMinResTol,fgNKrylovV); else AliInfo(Form("Undefined Iteritive Solver ID=%d, only %d are defined",fgIterSol,AliMinResSolve::kNSolvers)); // if (!res) return gkFailed; - for (int i=fNGloSize;i--;) fVecBGlo[i] = SOL[i]; - return gkMinRes; + for (int i=fNGloSize;i--;) fVecBGlo[i] = sol[i]; + return gkNoInversion; // } @@ -918,7 +1027,7 @@ Double_t AliMillePede2::GetParError(int iPar) const { // return error for parameter iPar if (fGloSolveStatus==gkInvert) { - double res = fMatCGlo->QuerryDiag(iPar); + double res = fMatCGlo->QueryDiag(iPar); if (res>=0) return TMath::Sqrt(res); } return -1.; @@ -943,7 +1052,7 @@ Int_t AliMillePede2::PrintGlobalParameters() const lGlobalCor = 0.0; // double dg; - if (fGloSolveStatus==gkInvert && TMath::Abs( (dg=fMatCGlo->QuerryDiag(i)) *fDiagCGlo[i]) > 0) { + if (fGloSolveStatus==gkInvert && TMath::Abs( (dg=fMatCGlo->QueryDiag(i)) *fDiagCGlo[i]) > 0) { lGlobalCor = TMath::Sqrt(TMath::Abs(1.0-1.0/(dg*fDiagCGlo[i]))); AliInfo(Form("%d\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f\t %6d", i,fInitPar[i],fInitPar[i]+fDeltaPar[i],fDeltaPar[i],fVecBGlo[i],lError,lGlobalCor,fProcPnt[i]));