Adding a fast parametric fitter
authorshahoian <shahoian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 28 Oct 2009 17:17:04 +0000 (17:17 +0000)
committershahoian <shahoian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 28 Oct 2009 17:17:04 +0000 (17:17 +0000)
STEER/AliParamSolver.cxx [new file with mode: 0644]
STEER/AliParamSolver.h [new file with mode: 0644]
STEER/STEERLinkDef.h
STEER/libSTEER.pkg

diff --git a/STEER/AliParamSolver.cxx b/STEER/AliParamSolver.cxx
new file mode 100644 (file)
index 0000000..c209f56
--- /dev/null
@@ -0,0 +1,316 @@
+#include "AliParamSolver.h"
+#include "AliSymMatrix.h"
+#include "AliLog.h"
+#include <TMath.h>
+
+ClassImp(AliParamSolver)
+
+//______________________________________________________________________________________
+AliParamSolver::AliParamSolver()
+: fMatrix(0),fSolGlo(0),fSolLoc(0),fMaxGlobal(0),fNGlobal(0),fNPoints(0),fMaxPoints(0),
+  fRHSGlo(0),fRHSLoc(0),fMatGamma(0),fMatG(0),fCovDGl(0)
+{ 
+  // default constructor
+}
+
+//______________________________________________________________________________________
+AliParamSolver::AliParamSolver(Int_t maxglo,Int_t locbuff)
+: fMatrix(0),fSolGlo(0),fSolLoc(0),fMaxGlobal(maxglo),fNGlobal(maxglo),fNPoints(0),fMaxPoints(0),
+  fRHSGlo(0),fRHSLoc(0),fMatGamma(0),fMatG(0),fCovDGl(0)
+{ 
+  // constructor for nglo globals
+  Init(locbuff);
+}
+
+//______________________________________________________________________________________
+AliParamSolver::AliParamSolver(AliParamSolver& src)
+  : TObject(src),fMatrix(0),fSolGlo(0),fSolLoc(0),fMaxGlobal(src.fMaxGlobal),fNGlobal(src.fNGlobal),
+    fNPoints(0),fMaxPoints(0),fRHSGlo(0),fRHSLoc(0),fMatGamma(0),fMatG(0),fCovDGl(0)
+{ 
+  // copy constructor 
+  if (src.fMatrix) {
+    Init(src.fMaxPoints);
+    (*this) = src;
+  }
+}
+
+//______________________________________________________________________________________
+AliParamSolver& AliParamSolver::operator=(const AliParamSolver& src)
+{
+  // assignment operator
+  if (this==&src) return *this;
+  TObject::operator=(src);
+  if (src.fMatrix && (fNGlobal!=src.fNGlobal || fMaxPoints<src.fNPoints)) {
+    fNGlobal   = src.fNGlobal;
+    fMaxGlobal = src.fMaxGlobal;
+    if (fMatrix)   delete   fMatrix;
+    if (fSolGlo)   delete[] fSolGlo;
+    if (fSolLoc)   delete[] fSolLoc;
+    if (fRHSGlo)   delete[] fRHSGlo;
+    if (fRHSLoc)   delete[] fRHSLoc;
+    if (fMatGamma) delete[] fMatGamma;
+    if (fMatG)     delete[] fMatG;
+    if (fCovDGl)   delete[] fCovDGl;
+    Init(src.fMaxPoints);
+  }
+  if (src.fMatrix) {
+    (*fMatrix) = *(src.fMatrix);
+    memcpy(fSolGlo,src.fSolGlo,fNGlobal*sizeof(double));
+    memcpy(fSolLoc,src.fSolLoc,fNPoints*sizeof(double));
+    memcpy(fRHSGlo,src.fRHSGlo,fNGlobal*sizeof(double));
+    memcpy(fRHSLoc,src.fRHSLoc,fNPoints*sizeof(double));
+    memcpy(fMatGamma,src.fMatGamma,fNPoints*sizeof(double));
+    memcpy(fMatG,src.fMatG,fNPoints*fNGlobal*sizeof(double));
+  }
+  //
+  return *this;
+}
+
+//______________________________________________________________________________________
+AliParamSolver::~AliParamSolver()
+{ 
+  // destructor
+  delete fMatrix;
+  delete[] fSolGlo;
+  delete[] fSolLoc;
+  delete[] fRHSGlo;
+  delete[] fRHSLoc;
+  delete[] fMatGamma;
+  delete[] fMatG;
+  delete[] fCovDGl;
+}
+
+//______________________________________________________________________________________
+Bool_t AliParamSolver::SolveGlobals(Bool_t obtainCov)
+{
+  if (fNPoints<fNGlobal/3) {
+    AliError(Form("Number of points: %d is not enough for %d globals",fNPoints,fNGlobal));
+    return kFALSE;
+  }
+  //
+  if (!TestBit(kBitGloSol)) {
+    if (!fMatrix->SolveChol(fRHSGlo, fSolGlo, obtainCov)) {
+      AliError("Solution Failed");
+      return kFALSE;
+    }
+    SetBit(kBitGloSol);
+    if (obtainCov) SetBit(kBitCInv);
+  }
+  return kTRUE;
+}
+
+//______________________________________________________________________________________
+Bool_t AliParamSolver::SolveLocals()
+{
+  if (TestBit(kBitLocSol)) return kTRUE;
+  if (!TestBit(kBitGloSol)) {
+    AliError("Cannot solve for Locals before SolveGlobals is called");
+    return kFALSE;
+  }
+  for (int i=fNPoints;i--;) {
+    double beta = fRHSLoc[i];
+    double *mtG = fMatG + i*fNGlobal;  // G_i
+    for (int j=fNGlobal;j--;) beta -= mtG[j]*fSolGlo[j];
+    fSolLoc[i] = beta/fMatGamma[i];   // Gamma^-1 * (beta - G_i * a)
+  }
+  SetBit(kBitLocSol);
+  return kTRUE;
+}
+
+//______________________________________________________________________________________
+AliSymMatrix* AliParamSolver::GetCovMatrix()
+{
+  if (!TestBit(kBitGloSol)) {
+    AliError("Cannot obtain Cov.Matrix before SolveGlobals is called");
+    return 0;
+  }
+  if (TestBit(kBitCInv)) return fMatrix;
+  //
+  if (fMatrix->InvertChol()) {
+    SetBit(kBitCInv);
+    return fMatrix;
+  }
+  return 0;
+}
+
+//______________________________________________________________________________________
+Bool_t AliParamSolver::Solve(Bool_t obtainCov)
+{
+  return (SolveGlobals(obtainCov) && SolveLocals()) ? kTRUE : kFALSE; 
+}
+
+//______________________________________________________________________________________
+void AliParamSolver::AddEquation(const Double_t* dGl,const Double_t *dLc,const Double_t *res, const Double_t *covI)
+{
+  // add the measured point to chi2 normal equations
+  // Input: 
+  // dGl : NGlo x 3 matrix of derivative for each of the 3 coordinates vs global params
+  // dLc : 3-vector of derivative for each coordinate vs local param
+  // res : residual of the point (extrapolated - measured)
+  // covI: 3 x (3+1)/2 matrix of the (inverted) errors (symmetric)
+  //
+  // The contribution of the point to chi2 is  V * covI * V 
+  // with component of the vector V(i) = dGl[i]*glob + dLc[i]*loc - res[i] 
+  //
+  // Instead of the NGlo + NMeasPoints size matrix we create only NGlo size matrix using the 
+  // reduction a la millepede : http://www.desy.de/~blobel/millepede1.ps
+  const double kTiny = 1e-16;
+  //
+  if (fNPoints+1 == fMaxPoints) ExpandStorage((fNPoints+1)*2);
+  ResetBit(kBitGloSol|kBitLocSol);
+  //
+  if (TestBit(kBitCInv)) { // solution was obtained and the matrix was inverted for previous points
+    fMatrix->InvertChol();
+    ResetBit(kBitCInv);
+  }
+  //
+  double *mtG   = fMatG + fNPoints*fNGlobal;  // G_i
+  double &beta  = fRHSLoc[fNPoints];
+  double &gamma = fMatGamma[fNPoints];
+  double cDl[3];
+  //
+  // cov * dR/dloc
+  cDl[kX] = covI[kXX]*dLc[kX] + covI[kXY]*dLc[kY] + covI[kXZ]*dLc[kZ];
+  cDl[kY] = covI[kXY]*dLc[kX] + covI[kYY]*dLc[kY] + covI[kYZ]*dLc[kZ];
+  cDl[kZ] = covI[kXZ]*dLc[kX] + covI[kYZ]*dLc[kY] + covI[kZZ]*dLc[kZ];
+  //
+  for (int i=fNGlobal;i--;) {
+    const double *dGli = dGl+i*3;  // derivatives of XYZ vs i-th global param
+    double *cDgi = fCovDGl+i*3;    // cov * dR/dGl_i
+    cDgi[kX] = covI[kXX]*dGli[kX] + covI[kXY]*dGli[kY] + covI[kXZ]*dGli[kZ];
+    cDgi[kY] = covI[kXY]*dGli[kX] + covI[kYY]*dGli[kY] + covI[kYZ]*dGli[kZ];
+    cDgi[kZ] = covI[kXZ]*dGli[kX] + covI[kYZ]*dGli[kY] + covI[kZZ]*dGli[kZ];
+    //
+    mtG[i]   = cDl[kX]*dGli[kX] + cDl[kY]*dGli[kY] + cDl[kZ]*dGli[kZ];  // dR/dGl_i * cov * dR/dLoc
+  }
+  beta    = res[kX]*cDl[kX] + res[kY]*cDl[kY] + res[kZ]*cDl[kZ];  //RHS: res*cov*dR/dloc
+  gamma   = dLc[kX]*cDl[kX] + dLc[kY]*cDl[kY] + dLc[kZ]*cDl[kZ];  //RHS: dR/dloc*cov*dR/dloc
+  double locSol  = TMath::Abs(gamma)<kTiny ? 0. : beta/gamma;     //local solution: gamma^-1 beta
+  //
+  AliSymMatrix &matC = *fMatrix;
+  for (int i=fNGlobal;i--;) {
+    const double *cDgi = fCovDGl+i*3;    // cov * dR/dGl_i
+    //
+    fRHSGlo[i] += cDgi[kX]*res[kX] + cDgi[kY]*res[kY] + cDgi[kZ]*res[kZ]      // b_i = dR/dGi * cov * res
+      -           mtG[i]*locSol;                                              //     - [G gamma^-1 beta ]_i
+    //
+    for (int j=i+1;j--;) {
+      //      const double *cDgj = fCovDGl+j*3;  // cov * dR/dGl_j
+      const double *dGlj = dGl+j*3;        // derivatives of XYZ vs i-th global param
+      matC(i,j) += dGlj[kX]*cDgi[kX] + dGlj[kY]*cDgi[kY] + dGlj[kZ]*cDgi[kZ]  // C_ij = dR/dGi * cov * dR/dGj
+       - mtG[i]*mtG[j]/gamma;                                                //        - [G gamma^-1 T(G) ]_ij      
+    }
+  }
+  //
+  fNPoints++;
+  //
+}
+
+//______________________________________________________________________________________
+void AliParamSolver::AddConstraint(Int_t parID, Double_t val, Double_t err2inv)
+{
+  // add gassian constriant to parameter parID
+  if (parID>=fNGlobal) {
+    AliError(Form("Attempt to constraint non-existing parameter %d",parID));
+    return;
+  }
+  //
+  (*fMatrix)(parID,parID) += err2inv;
+  fRHSGlo[parID] += val*err2inv;
+  //
+}
+
+//______________________________________________________________________________________
+void AliParamSolver::Init(Int_t npini)
+{
+  // create storage assuming maximum npini measured points
+  fMatrix = new AliSymMatrix(fMaxGlobal);
+  fSolGlo = new Double_t[fMaxGlobal];
+  fRHSGlo = new Double_t[fMaxGlobal];
+  //
+  fCovDGl = new Double_t[3*fMaxGlobal];
+  ExpandStorage(npini);
+  fMatrix->SetSizeUsed(fNGlobal);
+  //
+}
+
+//______________________________________________________________________________________
+void AliParamSolver::ExpandStorage(Int_t newSize)
+{
+  // increase space to newSize measured points
+  newSize = newSize>fMaxPoints ? newSize : fMaxPoints+1;
+  double* tmp;
+  tmp = new Double_t[newSize];
+  if (fSolLoc) delete[] fSolLoc; 
+  fSolLoc = tmp;
+  //
+  tmp = new Double_t[newSize];
+  if (fMatGamma) {
+    memcpy(tmp,fMatGamma,fNPoints*sizeof(Double_t));
+    delete[] fMatGamma;
+  }
+  fMatGamma = tmp;
+  //
+  tmp = new Double_t[newSize];
+  if (fRHSLoc) {
+    memcpy(tmp,fRHSLoc,fNPoints*sizeof(Double_t));
+    delete[] fRHSLoc;
+  }
+  fRHSLoc = tmp;
+  //
+  tmp = new Double_t[newSize*fMaxGlobal];
+  if (fMatG) {
+    memcpy(tmp,fMatG,fNPoints*fMaxGlobal*sizeof(Double_t));
+    delete[] fMatG;
+  }
+  fMatG = tmp;
+  //
+  fMaxPoints = newSize;
+  //
+}
+
+//______________________________________________________________________________________
+void AliParamSolver::Clear(Option_t*)
+{
+  fNPoints = 0;
+  fMatrix->Reset();
+  for (int i=fNGlobal;i--;) fRHSGlo[i] = 0;
+  ResetBit(kBitGloSol | kBitLocSol | kBitCInv);
+}
+
+//______________________________________________________________________________________
+void AliParamSolver::Print(Option_t*) const
+{
+  AliInfo(Form("Solver with %d globals for %d points",fNGlobal,fNPoints));
+}
+
+//______________________________________________________________________________________
+void AliParamSolver::SetNGlobal(Int_t n) 
+{
+  if (n>fMaxGlobal) {
+    AliError(Form("Maximum number of globals was set to %n",fMaxGlobal));
+    return;
+  }
+  fNGlobal = n;
+  fMatrix->SetSizeUsed(fNGlobal);
+}
+
+//______________________________________________________________________________________
+void AliParamSolver::SetMaxGlobal(Int_t n) 
+{
+  if (n>0 && n==fMaxGlobal) return;
+  fMaxGlobal = n;
+  fNGlobal = n;
+  if (fMatrix)   delete   fMatrix;
+  if (fSolGlo)   delete[] fSolGlo;
+  if (fSolLoc)   delete[] fSolLoc;
+  if (fRHSGlo)   delete[] fRHSGlo;
+  if (fRHSLoc)   delete[] fRHSLoc;
+  if (fMatGamma) delete[] fMatGamma;
+  if (fMatG)     delete[] fMatG;
+  if (fCovDGl)   delete[] fCovDGl;
+  n = TMath::Max(16,fMaxPoints);
+  Init(n);
+  //
+}
+
diff --git a/STEER/AliParamSolver.h b/STEER/AliParamSolver.h
new file mode 100644 (file)
index 0000000..9c082ef
--- /dev/null
@@ -0,0 +1,73 @@
+#ifndef ALIPARAMSOLVER_H
+#define ALIPARAMSOLVER_H
+
+/* ----------------------------------------------------------------------------------------
+   Class to solve a set of N linearized parametric equations of the type
+   Eq(k): sum_i=0^n { g_i G_ik }  + t_k T_k = res_k
+   whith n "global" parameters gi and one "local" parameter (per equation) t_k.
+   Each measured points provides 3 measured coordinates, with proper covariance matrix.
+
+   Used for Newton-Raphson iteration step in solution of non-linear parametric equations
+   F(g,t_k) - res_k = 0, with G_ik = dF(g,t_k)/dg_i and T_k = dF(g,t_k)/dt_k
+   Solution is obtained by elimination of local parameters via large (n+N) matrix partitioning 
+
+   Author: ruben.shahoyan@cern.ch
+-------------------------------------------------------------------------------------------*/ 
+
+#include <TObject.h>
+class AliSymMatrix;
+
+class AliParamSolver: public TObject
+{
+ public:
+  enum {kBitGloSol=BIT(14),kBitLocSol=BIT(15),kBitCInv=BIT(16)};
+  enum {kXX=0,kXY=1,kXZ=2,kYX=kXY,kYY=3,kYZ=4,kZX=kXZ,kZY=kYZ,kZZ=5};
+  enum {kX,kY,kZ};
+
+  AliParamSolver();
+  AliParamSolver(Int_t maxglo,Int_t locsize=16);
+  AliParamSolver(AliParamSolver& src);
+  AliParamSolver& operator=(const AliParamSolver& src);
+  ~AliParamSolver();
+  //
+  void    AddEquation(const Double_t* dGl,const Double_t *dLc,const Double_t *res,const Double_t *covI);
+  void    AddConstraint(Int_t parID, Double_t val, Double_t err2inv);
+  Bool_t  Solve(Bool_t obtainCov=kFALSE);
+  Bool_t  SolveGlobals(Bool_t obtainCov=kFALSE);
+  Bool_t  SolveLocals();
+  void    SetMaxGlobal(Int_t n);
+  void    SetNGlobal(Int_t n);
+  void   Clear(Option_t* = "");
+  void   Print(Option_t* = "") const;
+  //
+  Int_t   GetNGlobal()          const {return fNGlobal;}
+  Int_t   GetMaxGlobal()        const {return fMaxGlobal;}
+  AliSymMatrix* GetCovMatrix();
+  Double_t*     GetLocals()     const {return (Double_t*)fSolLoc;}
+  Double_t*     GetGlobals()    const {return (Double_t*)fSolGlo;}
+
+ protected:
+  void    Init(Int_t npini=16);
+  void    ExpandStorage(Int_t newSize);
+
+ protected:
+  AliSymMatrix*    fMatrix;            // final matrix for global parameters (C in MP)
+  Double_t*        fSolGlo;            // solution for globals ( vector a in MP)
+  Double_t*        fSolLoc;            // solution for locals  ( vector alpha in MP)
+  Int_t            fMaxGlobal;         // max number of globals can process
+  Int_t            fNGlobal;           // number of globals
+  Int_t            fNPoints;           // number of added points (=number of local parameters)
+  //
+  // temp storage
+  Int_t            fMaxPoints;         // buffer size for storage
+  Double_t*        fRHSGlo;            // RHS of globals (vector b in MP)
+  Double_t*        fRHSLoc;            // RHS of locals (vector beta in MP)
+  Double_t*        fMatGamma;          // diagonals of local partition (Gamma_i in MP)
+  Double_t*        fMatG;              // off-diagonals of local partition (G_i in MP)
+  Double_t*        fCovDGl;            // temporary matrix of cov*dR/dGlo
+  //
+  ClassDef(AliParamSolver,0)
+};
+
+
+#endif
index cd58085..00f3208 100644 (file)
 #pragma link C++ class AliSymMatrix+;
 #pragma link C++ class AliSymBDMatrix+;
 #pragma link C++ class AliRectMatrix+;
+#pragma link C++ class AliParamSolver+;
 
 #pragma link C++ class AliGRPManager+;
 #pragma link C++ class AliDCSArray+;
index 33ed5bf..a7b9e0a 100644 (file)
@@ -71,6 +71,7 @@ AliRecoParam.cxx AliDetectorRecoParam.cxx \
 AliMillePedeRecord.cxx AliMillePede2.cxx AliMatrixSq.cxx \
 AliVectorSparse.cxx AliMatrixSparse.cxx \
 AliSymMatrix.cxx AliSymBDMatrix.cxx AliRectMatrix.cxx AliMinResSolve.cxx \
+AliParamSolver.cxx \
 AliGRPManager.cxx \
 AliDCSArray.cxx AliLHCReader.cxx \
 AliCTPTimeParams.cxx AliCTPInputTimeParams.cxx