/**********************************************************************************************/
/* General class for alignment with large number of degrees of freedom */
/* Based on the original milliped2 by Volker Blobel */
+/* and AliMillepede class by Javier */
+/* Allows operations with large sparse matrices */
/* http://www.desy.de/~blobel/mptalks.html */
/* */
/* Author: ruben.shahoyan@cern.ch */
fNLocFits(0),
fNLocFitsRejected(0),
fNGloFix(0),
- fGloSolveStatus(gkFailed),
+ fGloSolveStatus(kFailed),
//
fChi2CutFactor(1.),
fChi2CutRef(1.),
//_____________________________________________________________________________________________
AliMillePede2::~AliMillePede2()
{
+ // destructor
CloseDataRecStorage();
CloseConsRecStorage();
//
//_____________________________________________________________________________________________
Int_t AliMillePede2::InitMille(int nGlo, int nLoc, int lNStdDev,double lResCut, double lResCutInit)
{
- //
+ // init all
if (nLoc>0) fNLocPar = nLoc;
if (nGlo>0) fNGloPar = nGlo;
if (lResCutInit>0) fResCutInit = lResCutInit;
//
fNGloSize = fNGloPar;
//
- try {
- //
- if (fgIsMatGloSparse) {fMatCGlo = new AliMatrixSparse(fNGloPar); fMatCGlo->SetSymmetric(kTRUE);}
- else fMatCGlo = new AliSymMatrix(fNGloPar);
- //
- fFillIndex = new Int_t[fNGloPar];
- fFillValue = new Double_t[fNGloPar];
- //
- fMatCLoc = new AliSymMatrix(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];
- //
- fInitPar = new Double_t[fNGloPar];
- fDeltaPar = new Double_t[fNGloPar];
- fSigmaPar = new Double_t[fNGloPar];
- fIsLinear = new Bool_t[fNGloPar];
- //
- fGlo2CGlo = new Int_t[fNGloPar];
- fCGlo2Glo = new Int_t[fNGloPar];
- }
- catch(bad_alloc&) {
- AliInfo(Form("Failed to allocate the memory for %d global and %d local parameters",fNGloPar,fNLocPar));
- return 0;
- }
+ if (fgIsMatGloSparse) {fMatCGlo = new AliMatrixSparse(fNGloPar); fMatCGlo->SetSymmetric(kTRUE);}
+ else fMatCGlo = new AliSymMatrix(fNGloPar);
+ //
+ fFillIndex = new Int_t[fNGloPar];
+ fFillValue = new Double_t[fNGloPar];
+ //
+ fMatCLoc = new AliSymMatrix(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];
+ //
+ fInitPar = new Double_t[fNGloPar];
+ fDeltaPar = new Double_t[fNGloPar];
+ fSigmaPar = new Double_t[fNGloPar];
+ fIsLinear = new Bool_t[fNGloPar];
+ //
+ fGlo2CGlo = new Int_t[fNGloPar];
+ fCGlo2Glo = new Int_t[fNGloPar];
//
memset(fVecBLoc ,0,fNLocPar*sizeof(Double_t));
memset(fDiagCGlo ,0,fNGloPar*sizeof(Double_t));
//_____________________________________________________________________________________________
Bool_t AliMillePede2::ImposeDataRecFile(const char* fname)
{
+ // set filename for records
CloseDataRecStorage();
SetDataRecFName(fname);
return InitDataRecStorage(kTRUE); // open in read mode
//_____________________________________________________________________________________________
Bool_t AliMillePede2::ImposeConsRecFile(const char* fname)
{
+ // set filename for constraints
CloseConsRecStorage();
SetConsRecFName(fname);
return InitConsRecStorage(kTRUE); // open in read mode
//_____________________________________________________________________________________________
void AliMillePede2::CloseDataRecStorage()
{
+ // close records file
if (fTreeData) {
if (fDataRecFile && fDataRecFile->IsWritable()) {
fDataRecFile->cd();
//_____________________________________________________________________________________________
void AliMillePede2::CloseConsRecStorage()
{
+ // close constraints file
if (fTreeConstr) {
if (fConsRecFile->IsWritable()) {
fConsRecFile->cd();
//_____________________________________________________________________________________________
void AliMillePede2::SetRecordWeight(double wgh)
{
+ // assign weight
if (fRecFileStatus<2) InitDataRecStorage(); // create a buffer to store the data
fRecord->SetWeight(wgh);
}
//_____________________________________________________________________________________________
void AliMillePede2::SetRecordRun(Int_t run)
{
+ // assign run
if (fRecFileStatus<2) InitDataRecStorage(); // create a buffer to store the data
fRecord->SetRunID(run);
}
//_____________________________________________________________________________________________
void AliMillePede2::SetLocalEquation(double *dergb, double *derlc, double lMeas, double lSigma)
{
+ // assing derivs of loc.eq.
if (fRecFileStatus<2) InitDataRecStorage(); // create a buffer to store the data
//
// write data of single measurement
//_____________________________________________________________________________________________
-void AliMillePede2::SetGlobalConstraint(double *dergb, double val, double sigma)
+void AliMillePede2::SetGlobalConstraint(const double *dergb, double val, double sigma)
{
// Define a constraint equation.
if (!fConsRecFile || !fConsRecFile->IsWritable()) InitConsRecStorage(); // create a buffer to store the data
}
//_____________________________________________________________________________________________
-void AliMillePede2::SetGlobalConstraint(const int *indgb, double *dergb, int ngb, double val,double sigma)
+void AliMillePede2::SetGlobalConstraint(const int *indgb, const double *dergb, int ngb, double val,double sigma)
{
// Define a constraint equation.
if (!fConsRecFile || !fConsRecFile->IsWritable()) InitConsRecStorage(); // create a buffer to store the data
//
if (par) for (int i=fNGloPar;i--;) par[i] = fInitPar[i]+fDeltaPar[i];
//
- if (fGloSolveStatus==gkInvert) { // errors on params are available
+ if (fGloSolveStatus==kInvert) { // errors on params are available
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;
AliInfo(Form("Global Fit Iteration#%2d (Local Fit Chi^2 cut factor: %.2f)",fIter,fChi2CutFactor));
//
if (!fNGloPar || !fTreeData) {
- AliInfo("No data was stored, aborting iteration");
+ AliInfo("No data was stored, stopping iteration");
return 0;
}
TStopwatch sw,sws;
printf("Solve %d |",fIter); sws.Print();
//
sw.Stop();
- AliInfo(Form("Iteration#%2d %s. CPU time: %.1f",fIter,fGloSolveStatus==gkFailed ? "Failed":"Converged",sw.CpuTime()));
- if (fGloSolveStatus==gkFailed) return 0;
+ AliInfo(Form("Iteration#%2d %s. CPU time: %.1f",fIter,fGloSolveStatus==kFailed ? "Failed":"Converged",sw.CpuTime()));
+ if (fGloSolveStatus==kFailed) return 0;
//
for (int i=fNGloPar;i--;) fDeltaPar[i] += fVecBGlo[i]; // Update global parameters values (for iterations)
//
if (!fgIsMatGloSparse) {
//
if (fNLagrangeConstraints==0) { // pos-def systems are faster to solve by Cholesky
- if ( ((AliSymMatrix*)fMatCGlo)->SolveChol(fVecBGlo, fgInvChol) ) return fgInvChol ? gkInvert:gkNoInversion;
+ if ( ((AliSymMatrix*)fMatCGlo)->SolveChol(fVecBGlo, fgInvChol) ) return fgInvChol ? kInvert:kNoInversion;
else AliInfo("Solution of Global Dense System by Cholesky failed, trying Gaussian Elimiation");
}
//
- if (((AliSymMatrix*)fMatCGlo)->SolveSpmInv(fVecBGlo, kTRUE)) return gkInvert;
+ if (((AliSymMatrix*)fMatCGlo)->SolveSpmInv(fVecBGlo, kTRUE)) return kInvert;
else AliInfo("Solution of Global Dense System by Gaussian Elimination failed, trying iterative methods");
}
// try to solve by minres
TVectorD sol(fNGloSize);
//
AliMinResSolve *slv = new AliMinResSolve(fMatCGlo,fVecBGlo);
- if (!slv) return gkFailed;
+ if (!slv) return kFailed;
//
Bool_t res = kFALSE;
if (fgIterSol == AliMinResSolve::kSolMinRes)
int defout = dup(1);
if (defout<0) {
AliInfo("Failed on dup");
- return gkFailed;
+ return kFailed;
}
int slvDump = open(faildump, O_RDWR|O_CREAT, 0666);
if (slvDump>=0) {
//
dup2(defout,1);
close(slvDump);
- close(defout);
printf("#Dumped failed matrix and RHS to %s\n",faildump);
}
else AliInfo("Failed on file open for matrix dumping");
- return gkFailed;
+ close(defout);
+ return kFailed;
}
for (int i=fNGloSize;i--;) fVecBGlo[i] = sol[i];
- return gkNoInversion;
+ return kNoInversion;
//
}
Double_t AliMillePede2::GetParError(int iPar) const
{
// return error for parameter iPar
- if (fGloSolveStatus==gkInvert) {
+ if (fGloSolveStatus==kInvert) {
double res = fMatCGlo->QueryDiag(iPar);
if (res>=0) return TMath::Sqrt(res);
}
lGlobalCor = 0.0;
//
double dg;
- if (fGloSolveStatus==gkInvert && TMath::Abs( (dg=fMatCGlo->QueryDiag(i)) *fDiagCGlo[i]) > 0) {
+ if (fGloSolveStatus==kInvert && 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]));
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff