Update for MillePede2 and related classes
authorshahoian <shahoian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 8 May 2009 20:24:15 +0000 (20:24 +0000)
committershahoian <shahoian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 8 May 2009 20:24:15 +0000 (20:24 +0000)
17 files changed:
STEER/AliMatrixSparse.cxx
STEER/AliMatrixSparse.h
STEER/AliMatrixSq.cxx
STEER/AliMatrixSq.h
STEER/AliMillePede2.cxx
STEER/AliMillePede2.h
STEER/AliMillePedeRecord.cxx [new file with mode: 0644]
STEER/AliMillePedeRecord.h [new file with mode: 0644]
STEER/AliMinResSolve.cxx
STEER/AliRectMatrix.cxx [new file with mode: 0644]
STEER/AliRectMatrix.h [new file with mode: 0644]
STEER/AliSymMatrix.cxx
STEER/AliSymMatrix.h
STEER/AliVectorSparse.cxx [new file with mode: 0644]
STEER/AliVectorSparse.h [new file with mode: 0644]
STEER/STEERLinkDef.h
STEER/libSTEER.pkg

index 4cf5252..26f4c42 100644 (file)
 #include "AliMatrixSparse.h"
+#include "TStopwatch.h"
 
-//___________________________________________________________
-AliVectorSparse::AliVectorSparse()
-  : fNElems(0),fIndex(0),fElems(0) {}
-
-//___________________________________________________________
-AliVectorSparse::AliVectorSparse(const AliVectorSparse& src)
-  : fNElems(src.fNElems),fIndex(0),fElems(0)
-{
-  fIndex = new UShort_t[fNElems];
-  fElems = new Double_t[fNElems];
-  memcpy(fIndex,src.fIndex,fNElems*sizeof(UShort_t));
-  memcpy(fElems,src.fElems,fNElems*sizeof(Double_t));
-}
-
-//___________________________________________________________
-void AliVectorSparse::Clear()
-{
-  delete[] fIndex; fIndex = 0;
-  delete[] fElems; fElems = 0;
-  fNElems = 0;
-}
-
-//___________________________________________________________
-Float_t AliMatrixSparse::GetDensity() const
-{
-  // get fraction of non-zero elements
-  Int_t nel = 0;
-  for (int i=GetSize();i--;) nel += GetRow(i)->GetNElems();
-  int den = IsSymmetric() ? (GetSize()+1)*GetSize()/2 : GetSize()*GetSize();
-  return float(nel)/den;
-}
-
-//___________________________________________________________
-AliVectorSparse& AliVectorSparse::operator=(const AliVectorSparse& src)
-{
-  if (&src==this) return *this;
-  Clear();
-  fNElems = src.fNElems;
-  fIndex = new UShort_t[fNElems];
-  fElems = new Double_t[fNElems];
-  memcpy(fIndex,src.fIndex,fNElems*sizeof(UShort_t));
-  memcpy(fElems,src.fElems,fNElems*sizeof(Double_t));
-  //
-  return *this;
-}
-
-//___________________________________________________________
-Double_t AliVectorSparse::FindIndex(Int_t ind) const
-{
-  // return an element with given index
-  //printf("V: findindex\n");
-  int first = 0;
-  int last = fNElems-1;
-  while (first<=last) {
-    int mid = (first+last)>>1;
-    if (ind>fIndex[mid]) first = mid+1;
-    else if (ind<fIndex[mid]) last = mid-1;
-    else return fElems[mid];
-  }
-  return 0.0;
-}
-
-//___________________________________________________________
-void AliVectorSparse::SetToZero(Int_t ind)
-{
-  // set element to 0 if it was already defined
-  int first = 0;
-  int last = fNElems-1;
-  while (first<=last) {
-    int mid = (first+last)>>1;
-    if (ind>fIndex[mid]) first = mid+1;
-    else if (ind<fIndex[mid]) last = mid-1;
-    else {fElems[mid] = 0.; return;}
-  }
-}
-
-//___________________________________________________________
-Double_t& AliVectorSparse::FindIndexAdd(Int_t ind)
-{
-  // increment an element with given index
-  //printf("V: findindexAdd\n");
-  int first = 0;
-  int last = fNElems-1;
-  while (first<=last) {
-    int mid = (first+last)>>1;
-    if (ind>fIndex[mid]) first = mid+1;
-    else if (ind<fIndex[mid]) last = mid-1;
-    else return fElems[mid];
-  }
-  // need to insert a new element
-  UShort_t *arrI = new UShort_t[fNElems+1];
-  memcpy(arrI,fIndex,first*sizeof(UShort_t));
-  arrI[first] = ind;
-  memcpy(arrI+first+1,fIndex+first,(fNElems-first)*sizeof(UShort_t));
-  delete[] fIndex;
-  fIndex = arrI;
-  //
-  Double_t   *arrE = new Double_t[fNElems+1];
-  memcpy(arrE,fElems,first*sizeof(Double_t));
-  arrE[first] = 0;
-  memcpy(arrE+first+1,fElems+first,(fNElems-first)*sizeof(Double_t));
-  delete[] fElems;
-  fElems = arrE;
-  //
-  fNElems++;
-  return fElems[first];
-  //
-}
-
-//__________________________________________________________
-void AliVectorSparse::ReSize(Int_t sz,Bool_t copy)
-{
-  if (sz<1) {Clear(); return;}
-    // need to insert a new element
-  UShort_t *arrI = new UShort_t[sz];
-  Double_t *arrE = new Double_t[sz];
-  memset(arrI,0,sz*sizeof(UShort_t));
-  memset(arrE,0,sz*sizeof(Double_t));
-  //
-  if (copy && fIndex) {
-    int cpsz = TMath::Min(fNElems,sz);
-    memcpy(arrI,fIndex,cpsz*sizeof(UShort_t));
-    memcpy(arrE,fElems,cpsz*sizeof(Double_t));
-  }
-  delete[] fIndex;
-  delete[] fElems;
-  fIndex = arrI;
-  fElems = arrE;
-  fNElems = sz;
-  //
-}
+/**********************************************************************************************/
+/* Sparse matrix class, used as a global matrix for AliMillePede2                             */
+/*                                                                                            */ 
+/* Author: ruben.shahoyan@cern.ch                                                             */
+/*                                                                                            */ 
+/**********************************************************************************************/
 
-//__________________________________________________________
-void AliVectorSparse::SortIndices(Bool_t valuesToo)
-{
-  // sort indices in increasing order. Used to fix the row after ILUk decomposition
-  for (int i=fNElems;i--;) for (int j=i;j--;) if (fIndex[i]<fIndex[j]) { //swap
-       UShort_t tmpI = fIndex[i]; fIndex[i] = fIndex[j]; fIndex[j]=tmpI;
-       if (valuesToo) {Double_t tmpV = fElems[i];fElems[i]=fElems[j];fElems[j]=tmpV;}
-      }
-}
-
-//__________________________________________________________
-void AliVectorSparse::Print(Option_t* )  const
-{
-  printf("|");
-  for (int i=0;i<fNElems;i++) printf("%2d:%+.2e|",fIndex[i],fElems[i]);
-  printf("|\n");
-}
-
-
-///////////////////////////////////////////////////////////////////////////////////////////
 //___________________________________________________________
 ClassImp(AliMatrixSparse)
 
@@ -256,6 +112,52 @@ void AliMatrixSparse::MultiplyByVec(Double_t* vecIn, Double_t* vecOut) const
 //___________________________________________________________
 void AliMatrixSparse::SortIndices(Bool_t valuesToo)
 {
+  TStopwatch sw; 
+  sw.Start();
+  printf("AliMatrixSparse:sort>>\n");
   // sort columns in increasing order. Used to fix the matrix after ILUk decompostion
   for (int i=GetSize();i--;) GetRow(i)->SortIndices(valuesToo);
+  sw.Stop();
+  sw.Print();
+  printf("AliMatrixSparse:sort<<\n");
+}
+
+//___________________________________________________________
+void AliMatrixSparse::AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n)
+{
+  // for sym. matrix count how many elems to add have row>=col and assign excplicitly 
+  // those which have row<col
+  //   
+  // range in increasing order of indices
+  for (int i=n;i--;) {
+    for (int j=i;j>=0;j--) {
+      if (indc[j]>indc[i]) { // swap
+       int    ti = indc[i]; indc[i] = indc[j]; indc[j] = ti;
+       double tv = valc[i]; valc[i] = valc[j]; valc[j] = tv;
+      }
+    }
+  }
+  //
+  int ni=n;
+  if (IsSymmetric()) {
+    while(ni--) {
+      if (indc[ni]>r) (*this)(indc[ni],r) += valc[ni]; 
+      else break;  // use the fact that the indices are ranged in increasing order
+    }
+  }
+  //
+  if (ni<0) return;
+  AliVectorSparse* row = GetRowAdd(r);
+  row->Add(valc,indc,ni+1);
+}
+
+//___________________________________________________________
+Float_t AliMatrixSparse::GetDensity() const
+{
+  // get fraction of non-zero elements
+  Int_t nel = 0;
+  for (int i=GetSize();i--;) nel += GetRow(i)->GetNElems();
+  int den = IsSymmetric() ? (GetSize()+1)*GetSize()/2 : GetSize()*GetSize();
+  return float(nel)/den;
 }
+
index 20874a5..9dcf519 100644 (file)
@@ -2,61 +2,17 @@
 #define ALIMATRIXSPARSE_H
 
 #include "AliMatrixSq.h"
+#include "AliVectorSparse.h"
 
 
-///////////////////////////////////////////////////////////////////////////////////////
-class AliVectorSparse {
- public:
-  AliVectorSparse();
-  AliVectorSparse(const AliVectorSparse& src);
-  virtual ~AliVectorSparse() {Clear();}
-  virtual void Print(Option_t* option="")                 const;
-  //
-  Int_t     GetNElems()                                   const {return fNElems;}
-  UShort_t *GetIndices()                                  const {return fIndex;}
-  Double_t *GetElems()                                    const {return fElems;}
-  UShort_t& GetIndex(Int_t i)                                   {return fIndex[i];}
-  Double_t& GetElem(Int_t i)                              const {return fElems[i];}
-  void      Clear();
-  void      Reset()                                             {memset(fElems,0,fNElems*sizeof(Double_t));}
-  void      ReSize(Int_t sz,Bool_t copy=kFALSE);
-  void      SortIndices(Bool_t valuesToo=kFALSE);
-  //
-  AliVectorSparse& operator=(const AliVectorSparse& src);
-  //
-  virtual Double_t         operator()(Int_t ind)         const;
-  virtual Double_t&        operator()(Int_t ind);
-  virtual void             SetToZero(Int_t ind);
-  Double_t                 FindIndex(Int_t ind)          const;
-  Double_t&                FindIndexAdd(Int_t ind);
-  //
-  Int_t     GetLastIndex()                               const {return fIndex[fNElems-1];}
-  Double_t  GetLastElem()                                const {return fElems[fNElems-1];}
-  Double_t &GetLastElem()                                      {return fElems[fNElems-1];}
-  //
- protected:
-  Int_t            fNElems;   // 
-  UShort_t*        fIndex;    // Index of stored elems
-  Double_t*        fElems;    // pointer on elements
-};
-
+/**********************************************************************************************/
+/* Sparse matrix class, used as a global matrix for AliMillePede2                             */
+/*                                                                                            */ 
+/* Author: ruben.shahoyan@cern.ch                                                             */
+/*                                                                                            */ 
+/**********************************************************************************************/
 
-//___________________________________________________
-inline Double_t AliVectorSparse::operator()(Int_t ind) const
-{
-  return FindIndex(ind);
-}
 
-//___________________________________________________
-inline Double_t& AliVectorSparse::operator()(Int_t ind)
-{
-  return FindIndexAdd(ind);
-}
-
-//////////////////////////////////////////////////////////////////////////////////////
-
-/////////////////////////////////////////////////////////////////////////////////////////////
-// Sparse matrix class
 //
 class AliMatrixSparse : public AliMatrixSq 
 {
@@ -82,8 +38,10 @@ class AliMatrixSparse : public AliMatrixSq
   Double_t&        DiagElem(Int_t r);
   void             SortIndices(Bool_t valuesToo=kFALSE);
   //
+  void MultiplyByVec(TVectorD &vecIn, TVectorD &vecOut) const; 
   void MultiplyByVec(Double_t* vecIn, Double_t* vecOut) const;
-  void MultiplyByVec(TVectorD &vecIn, TVectorD &vecOut) const {MultiplyByVec((Double_t*)vecIn.GetMatrixArray(),(Double_t*)vecOut.GetMatrixArray());}
+  //
+  void AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n);
   //
  protected:
   //
@@ -92,6 +50,12 @@ class AliMatrixSparse : public AliMatrixSq
   ClassDef(AliMatrixSparse,0)
 };
 
+//___________________________________________________
+inline void AliMatrixSparse::MultiplyByVec(TVectorD &vecIn, TVectorD &vecOut) const 
+{
+  MultiplyByVec((Double_t*)vecIn.GetMatrixArray(),(Double_t*)vecOut.GetMatrixArray());
+}
+
 //___________________________________________________
 inline void AliMatrixSparse::SetToZero(Int_t row,Int_t col)
 {
@@ -140,5 +104,6 @@ inline Double_t &AliMatrixSparse::DiagElem(Int_t row)
   //
 }
 
+
 #endif
 
index d315e54..f5e7232 100644 (file)
@@ -34,10 +34,10 @@ void AliMatrixSq::PrintCOO() const
   // get number of non-zero elements
   int nnz = 0;
   int sz = GetSize();
-  for (int ir=0;ir<sz;ir++) for (int ic=0;ic<sz;ic++) if (Querry(ir,ic)!=0) nnz++;
+  for (int ir=0;ir<sz;ir++) for (int ic=0;ic<sz;ic++) if (Query(ir,ic)!=0) nnz++;
   //
   printf("%d %d %d\n",sz,sz,nnz);
   double vl;
-  for (int ir=0;ir<sz;ir++) for (int ic=0;ic<sz;ic++) if ((vl=Querry(ir,ic))!=0) printf("%d %d %f\n",ir,ic,vl);
+  for (int ir=0;ir<sz;ir++) for (int ic=0;ic<sz;ic++) if ((vl=Query(ir,ic))!=0) printf("%d %d %f\n",ir,ic,vl);
   //
 }
index 09f0828..0cc7701 100644 (file)
@@ -15,13 +15,14 @@ class AliMatrixSq : public TMatrixDBase {
   //
   virtual void  Clear(Option_t* option="")                       = 0;//{Error("Clear","Dummy");}
   //
-  virtual       Double_t      Querry(Int_t rown, Int_t coln)     const {return operator()(rown,coln);}
+  virtual       Double_t      Query(Int_t rown, Int_t coln)     const {return operator()(rown,coln);}
   virtual       Double_t      operator()(Int_t rown, Int_t coln) const = 0;//{Error("(i,j)","Dummy");return 0;}
   virtual       Double_t&     operator()(Int_t rown, Int_t coln) = 0;//{Error("(i,j)","Dummy");return 0;}
   //
-  virtual       Double_t      QuerryDiag(Int_t rc)               const {return DiagElem(rc);}
+  virtual       Double_t      QueryDiag(Int_t rc)               const {return DiagElem(rc);}
   virtual       Double_t      DiagElem(Int_t r)                  const = 0;
   virtual       Double_t&     DiagElem(Int_t r)                  = 0;
+  virtual       void          AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n) = 0;
   //
   virtual void  Print(Option_t* option="")           const       = 0;//{Error("Print","Dummy");}
   virtual void  Reset()                                          = 0;
index aca79e2..62a3f94 100644 (file)
+/**********************************************************************************************/
+/* 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 <TStopwatch.h>
+#include <TFile.h>
+#include <TMath.h>
+#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;i<rhs.GetSize();i++) {
-      double val;
-      int    ind;
-      rhs.GetIndexValue(i,ind,val);
-      AddIndexValue(ind,val);
-    }
-  }
-  return *this;
-}
 
-//_____________________________________________________________________________________________
-AliMillePedeRecord::~AliMillePedeRecord() {delete[] fIndex; delete[] fValue;}
-
-//_____________________________________________________________________________________________
-void AliMillePedeRecord::Print(const Option_t *) const
-{
-  if (!fSize) {AliInfo("No data"); return;}
-  int cnt=0,point=0;
-  //  
-  while(cnt<fSize) {
-    //
-    double resid = fValue[cnt++];
-    double *derLoc = GetValue()+cnt;
-    int    *indLoc = GetIndex()+cnt;
-    int     nLoc = 0;
-    while(!IsWeight(cnt)) {nLoc++;cnt++;}
-    double weight = GetValue(cnt++);
-    double *derGlo = GetValue()+cnt;
-    int    *indGlo = GetIndex()+cnt;
-    int     nGlo = 0;
-    while(!IsResidual(cnt) && cnt<fSize) {nGlo++; cnt++;} 
-    //
-    printf("\n*** Point#%2d | Residual = %+.4e | Weight = %+.4e\n",point++,resid,weight);
-    printf("Locals : "); 
-    for (int i=0;i<nLoc;i++) printf("[%5d] %+.4e|",indLoc[i],derLoc[i]); printf("\n");
-    printf("Globals: "); 
-    for (int i=0;i<nGlo;i++) printf("[%5d] %+.4e|",indGlo[i],derGlo[i]); printf("\n");
-    //
-  }
-  //
-}
+using namespace std;
 
-//_____________________________________________________________________________________________
-void AliMillePedeRecord::Expand(int bfsize)
-{
-  // add extra space 
-  bfsize = TMath::Max(bfsize,GetBufferSize());
-  int *tmpI = new int[bfsize];
-  memcpy(tmpI,fIndex, fSize*sizeof(int));
-  delete fIndex;
-  fIndex = tmpI;
-  //
-  double *tmpD = new double[bfsize];
-  memcpy(tmpD,fValue, fSize*sizeof(double));
-  delete fValue;
-  fValue = tmpD;
-  //
-  SetBufferSize(bfsize);
-}
 
-//----------------------------------------------------------------------------------------------
-//----------------------------------------------------------------------------------------------
-//----------------------------------------------------------------------------------------------
-//----------------------------------------------------------------------------------------------
 ClassImp(AliMillePede2)
 
-Bool_t   AliMillePede2::fgIsMatGloSparse = kFALSE;   // use faster dense matrix by default
-Int_t    AliMillePede2::fgMinResCondType = kTRUE;        // No preconditioner by default
-Double_t AliMillePede2::fgMinResTol      = 1.e-8;   // default tolerance
-Int_t    AliMillePede2::fgMinResMaxIter  = 3000;     // default max number of iterations 
+Bool_t   AliMillePede2::fgInvChol        = kTRUE;     // Invert global matrix with Cholesky solver
+Bool_t   AliMillePede2::fgWeightSigma    = kTRUE;     // weight local constraint by module statistics
+Bool_t   AliMillePede2::fgIsMatGloSparse = kFALSE;    // use faster dense matrix by default
+Int_t    AliMillePede2::fgMinResCondType = 1;         // Jacoby preconditioner by default
+Double_t AliMillePede2::fgMinResTol      = 1.e-11;    // default tolerance
+Int_t    AliMillePede2::fgMinResMaxIter  = 10000;     // default max number of iterations 
 Int_t    AliMillePede2::fgIterSol        = AliMinResSolve::kSolMinRes; // default iterative solver
-Int_t    AliMillePede2::fgNKrylovV       = 30;       // default number of Krylov vectors to keep
+Int_t    AliMillePede2::fgNKrylovV       = 240;        // default number of Krylov vectors to keep
 
 //_____________________________________________________________________________________________
 AliMillePede2::AliMillePede2() 
@@ -109,6 +43,7 @@ AliMillePede2::AliMillePede2()
   fMaxIter(10),
   fNStdDev(3),
   fNGloConstraints(0),
+  fNLagrangeConstraints(0),
   fNLocFits(0),
   fNLocFitsRejected(0),
   fNGloFix(0),
@@ -118,8 +53,10 @@ AliMillePede2::AliMillePede2()
   fChi2CutRef(1.),
   fResCutInit(100.),
   fResCut(100.),
-  fMinPntValid(0),
+  fMinPntValid(1),
 //
+  fNGroupsSet(0),
+  fParamGrID(0),
   fProcPnt(0),
   fVecBLoc(0),
   fDiagCGlo(0),
@@ -136,7 +73,10 @@ AliMillePede2::AliMillePede2()
   fMatCLoc(0),
   fMatCGlo(0),
   fMatCGloLoc(0),
-//
+  //
+  fFillIndex(0),
+  fFillValue(0),
+  //
   fDataRecFName("/tmp/mp2_data_records.root"),
   fRecord(0),
   fDataRecFile(0),  
@@ -146,18 +86,22 @@ AliMillePede2::AliMillePede2()
   fTreeConstr(0),
   fConsRecFile(0),
   fCurrRecDataID(0),
-  fCurrRecConstrID(0)
+  fCurrRecConstrID(0),
+  fLocFitAdd(kTRUE)
 {}
 
 //_____________________________________________________________________________________________
 AliMillePede2::AliMillePede2(const AliMillePede2& src) : 
   TObject(src),fNLocPar(0),fNGloPar(0),fNGloSize(0),fNLocEquations(0),fIter(0),
-  fMaxIter(10),fNStdDev(3),fNGloConstraints(0),fNLocFits(0),fNLocFitsRejected(0),
+  fMaxIter(10),fNStdDev(3),fNGloConstraints(0),fNLagrangeConstraints(0),
+  fNLocFits(0),fNLocFitsRejected(0),
   fNGloFix(0),fGloSolveStatus(0),fChi2CutFactor(0),fChi2CutRef(0),fResCutInit(0),
-  fResCut(0),fMinPntValid(0),fProcPnt(0),fVecBLoc(0),fDiagCGlo(0),fVecBGlo(0),
+  fResCut(0),fMinPntValid(1),fNGroupsSet(0),fParamGrID(0),fProcPnt(0),fVecBLoc(0),fDiagCGlo(0),fVecBGlo(0),
   fInitPar(0),fDeltaPar(0),fSigmaPar(0),fIsLinear(0),fConstrUsed(0),fGlo2CGlo(0),fCGlo2Glo(0),
-  fMatCLoc(0),fMatCGlo(0),fMatCGloLoc(0),fDataRecFName(0),fRecord(0),fDataRecFile(0),
-  fTreeData(0),fConstrRecFName(0),fTreeConstr(0),fConsRecFile(0),fCurrRecDataID(0),fCurrRecConstrID(0)
+  fMatCLoc(0),fMatCGlo(0),fMatCGloLoc(0),fFillIndex(0),fFillValue(0),
+  fDataRecFName(0),fRecord(0),fDataRecFile(0),
+  fTreeData(0),fConstrRecFName(0),fTreeConstr(0),fConsRecFile(0),fCurrRecDataID(0),
+  fCurrRecConstrID(0),fLocFitAdd(kTRUE)
 {printf("Dummy\n");}
 
 //_____________________________________________________________________________________________
@@ -166,6 +110,7 @@ AliMillePede2::~AliMillePede2()
   CloseDataRecStorage();
   CloseConsRecStorage();
   //
+  delete[] fParamGrID;
   delete[] fProcPnt;
   delete[] fVecBLoc;
   delete[] fDiagCGlo;
@@ -177,6 +122,8 @@ AliMillePede2::~AliMillePede2()
   delete[] fCGlo2Glo;
   delete[] fIsLinear;
   delete[] fConstrUsed;
+  delete[] fFillIndex;
+  delete[] fFillValue;
   //
   delete fRecord;
   delete fMatCLoc;
@@ -201,9 +148,13 @@ Int_t AliMillePede2::InitMille(int nGlo, int nLoc, int lNStdDev,double lResCut,
     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 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 +180,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)
 {
@@ -392,7 +334,10 @@ void AliMillePede2::SetLocalEquation(double *dergb, double *derlc, double lMeas,
   fRecord->AddWeight( 1.0/lSigma/lSigma );
   //
   // Idem for global parameters
-  for (int i=0;i<fNGloPar;i++) if (dergb[i]!=0.0) {fRecord->AddIndexValue(i,dergb[i]); dergb[i] = 0.0;}
+  for (int i=0;i<fNGloPar;i++) if (dergb[i]!=0.0) {
+    fRecord->AddIndexValue(i,dergb[i]); dergb[i] = 0.0;
+    fRecord->MarkGroup(fParamGrID[i]);
+  }
   //
 }
 
@@ -421,31 +366,35 @@ void AliMillePede2::SetLocalEquation(int *indgb, double *dergb, int ngb, int *in
   //
 }
 
+
 //_____________________________________________________________________________________________
-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
+  fRecord->AddWeight(sigma);
   for (int i=0; i<fNGloPar; i++) if (dergb[i]!=0)  fRecord->AddIndexValue(i,dergb[i]);
   fNGloConstraints++;
+  if (sigma==0) 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
+  fRecord->AddWeight(sigma);   // dummy
   for (int i=0; i<ngb; i++) if (dergb[i]!=0)  fRecord->AddIndexValue(indgb[i],dergb[i]);
   fNGloConstraints++;
+  if (sigma==0) fNLagrangeConstraints++;
   SaveRecordConstraint();
   //
 }
@@ -459,16 +408,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 +426,84 @@ Int_t AliMillePede2::LocalFit(double *localParams)
   while(cnt<recSz) {  // Transfer the measurement records to matrices
     //
     // extract addresses of residual, weight and pointers on local and global derivatives for each point
-    if (nRefSize<=nPoints) {
-      nRefSize = 2*(nPoints+1); 
-      RefLoc.Set(nRefSize); 
-      RefGlo.Set(nRefSize);
-      nRefLoc.Set(nRefSize); 
-      nRefGlo.Set(nRefSize);
+    if (nrefSize<=nPoints) {
+      int *tmpA = 0;
+      nrefSize = 2*(nPoints+1); 
+      tmpA = refLoc;  refLoc  = new Int_t[nrefSize]; if (tmpA) memcpy(refLoc,tmpA,nPoints*sizeof(int));
+      tmpA = refGlo;  refGlo  = new Int_t[nrefSize]; if (tmpA) memcpy(refGlo,tmpA,nPoints*sizeof(int));
+      tmpA = nrefLoc; nrefLoc = new Int_t[nrefSize]; if (tmpA) memcpy(nrefLoc,tmpA,nPoints*sizeof(int));
+      tmpA = nrefGlo; nrefGlo = new Int_t[nrefSize]; if (tmpA) memcpy(nrefGlo,tmpA,nPoints*sizeof(int));
     }
     //
-    RefLoc[nPoints]  = ++cnt;
+    refLoc[nPoints]  = ++cnt;
     int nLoc = 0;
     while(!fRecord->IsWeight(cnt)) {nLoc++; cnt++;}
-    nRefLoc[nPoints] = nLoc;
+    nrefLoc[nPoints] = nLoc;
     //
-    RefGlo[nPoints]  = ++cnt;
+    refGlo[nPoints]  = ++cnt;
     int nGlo = 0;
     while(!fRecord->IsResidual(cnt) && cnt<recSz) {nGlo++; cnt++;} 
-    nRefGlo[nPoints] = nGlo;
+    nrefGlo[nPoints] = nGlo;
     //
     nPoints++;
   }
   double vl;
   //
   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 );    
+    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];
+      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
   //
   // first try to solve by faster Cholesky decomposition, then by Gaussian elimination
-  if (!MatCLoc.SolveChol(fVecBLoc,kTRUE)) {
+  if (!matCLoc.SolveChol(fVecBLoc,kTRUE)) {
     AliInfo("Failed to solve locals by Cholesky, traying Gaussian Elimination");
-    if (!MatCLoc.SolveSpmInv(fVecBLoc,kTRUE)) { 
+    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
+  //  printf("locfit: "); for (int i=0;i<fNLocPar;i++) printf("%+e |",fVecBLoc[i]); printf("\n");
   if (localParams) for (int i=fNLocPar; i--;) {
       localParams[2*i]   = fVecBLoc[i];
-      localParams[2*i+1] = TMath::Sqrt(TMath::Abs(MatCLoc.QuerryDiag(i)));
+      localParams[2*i+1] = TMath::Sqrt(TMath::Abs(matCLoc.QueryDiag(i)));
     }
   //
   float lChi2 = 0;
   int   nEq   = 0;
   //
   for (int ip=nPoints;ip--;) {  // Calculate residuals
-    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];
+    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];
     //
     // 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 +511,10 @@ 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++;      
       return 0;
     }
     // 
@@ -573,17 +522,22 @@ Int_t AliMillePede2::LocalFit(double *localParams)
     nEq++;                        // number of equations                       
   } // end of Calculate residuals
   //
-  //
   int nDoF = nEq-fNLocPar;
   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++;      
     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 +545,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 );    
+    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 ( (vl = weight*derGlo[ig]*derGlo[jg])!=0 ) {
+         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=fNLocPar;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 +596,48 @@ Int_t AliMillePede2::LocalFit(double *localParams)
   // calculate fMatCGlo -= fMatCGloLoc * fMatCLoc * fMatCGloLoc^T
   // and       fVecBGlo -= fMatCGloLoc * fVecBLoc
   //
+  //-------------------------------------------------------------- >>>
   double vll;
   for (int iCIDg=0; iCIDg<nGloInFit; iCIDg++) {
     int iIDg = fCGlo2Glo[iCIDg];
     //
     vl = 0;
-    for (int kl=0;kl<fNLocPar;kl++) if ( (vll = MatCGloLoc(iCIDg,kl))!=0) vl += vll*fVecBLoc[kl];
-    if  (vl!=0) fVecBGlo[iIDg] -= vl;
+    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;
     //
+    int nfill = 0;
     for (int jCIDg=0;jCIDg<=iCIDg; jCIDg++) {  
       int jIDg = fCGlo2Glo[jCIDg];
       //
       vl = 0;
+      Double_t *rowGLJDg =  matCGloLoc(jCIDg);
       for (int kl=0;kl<fNLocPar;kl++) {
        // diag terms
-       if ( (vll=MatCGloLoc(iCIDg,kl)*MatCGloLoc(jCIDg,kl))!=0 ) vl += MatCLoc.QuerryDiag(kl)*vll;
+       if ( (vll=rowGLIDg[kl]*rowGLJDg[kl])!=0 ) vl += matCLoc.QueryDiag(kl)*vll;
        //
        // off-diag terms
        for (int ll=0;ll<kl;ll++) {
-         if ( (vll=MatCGloLoc(iCIDg,kl)*MatCGloLoc(jCIDg,ll))!=0 ) vl += MatCLoc(kl,ll)*vll;
-         if ( (vll=MatCGloLoc(iCIDg,ll)*MatCGloLoc(jCIDg,kl))!=0 ) vl += MatCLoc(kl,ll)*vll;
+         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 (vl!=0) MatCGlo(iIDg,jIDg) -= vl; 
-      //
+      if (vl!=0) {
+       fFillIndex[nfill]   = jIDg;
+       fFillValue[nfill++] = fLocFitAdd ? -vl : vl;
+      }
     }
+    if (nfill)         matCGlo.AddToRow(iIDg,fFillValue,fFillIndex,nfill);
   }
   //
   // reset compressed index array
   //
-  for (int i=nGloInFit;i--;) { fGlo2CGlo[ fCGlo2Glo[i] ] = -1; fCGlo2Glo[i] = -1;}
+  for (int i=nGloInFit;i--;) { 
+    fGlo2CGlo[ fCGlo2Glo[i] ] = -1;
+    fCGlo2Glo[i] = -1;
+  }
   //
+  //---------------------------------------------------- <<<
   return 1;
 }
 
@@ -704,9 +673,9 @@ Int_t AliMillePede2::GlobalFit(Double_t *par, Double_t *error, Double_t *pull)
   if (par) for (int i=fNGloPar;i--;) par[i] = fInitPar[i]+fDeltaPar[i]; 
   //
   if (fGloSolveStatus==gkInvert) { // errors on params are available
-    if (error) for (int i=fNGloPar;i--;) error[i] = fProcPnt[i]>0 ? 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 +692,186 @@ 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; i<fNGloConstraints; i++) {
+    ReadRecordConstraint(i);
+    if ( fRecord->GetValue(1)==0 ) 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<ndr;i++) {
     ReadRecordData(i);
     LocalFit();
+    if ( (i%int(0.1*ndr)) == 0) printf("%.1f%% of local fits done\n", double(100.*i)/ndr);
   }
+  swt.Stop();
+  printf("%ld local fits done: ", ndr);
+  swt.Print();
+  sw.Start(kFALSE);
+  //
   //
+  // ---------------------- Reject parameters with low statistics ------------>>
   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);
+    //
+    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 (fProcPnt[i]>=fMinPntValid) continue;
+      fProcPnt[i] = 0;
+      // check if the group is already accounted
+      Bool_t fnd = kFALSE;
+      for (int j=nFixedGroups;j--;) if (fixGroups[j]==grID) {fnd=kTRUE; break;}
+      if (fnd) continue;
+      if (nFixedGroups>=fixArrSize) {fixArrSize*=2; fixGroups.Set(fixArrSize);}
+      fixGroups[nFixedGroups++] = grID; // add group to fix
+      //
+    }
+    // 2) loop over records and add contributions of fixed groups with negative sign
+    fLocFitAdd = kFALSE;
+    //
+    for (Long_t i=0;i<ndr;i++) {
+      ReadRecordData(i);
+      Bool_t suppr = kFALSE;
+      for (int ifx=nFixedGroups;ifx--;)if (fRecord->IsGroupPresent(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;i<nFixedGroups;i++) printf("%d ",fixGroups[i]); printf("\n");
+    }
+    swsup.Stop();
+    swsup.Print();
+  }
+  // ---------------------- Reject parameters with low statistics ------------<<
+  //
+  // add large number to diagonal of fixed params  
   //
-  //  Account for fixed parameters
   for (int i=fNGloPar;i--;) { // // Reset row and column of fixed params and add 1/sig^2 to free ones
-    if ( fSigmaPar[i] <= 0. || fProcPnt[i]<fMinPntValid) {
+    if (fProcPnt[i]<1) {
       fNGloFix++; 
-      fProcPnt[i] *= -1;
       fVecBGlo[i] = 0.;
-      if (IsGlobalMatSparse()) {
-       AliVectorSparse& row = *((AliMatrixSparse*)fMatCGlo)->GetRow(i);
-       row.Clear();
-       row(i) = float(fNLocEquations);
-       for (int j=i+1;j<fNGloPar;j++) ((AliMatrixSparse*)fMatCGlo)->SetToZero(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) = 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; i<fNGloConstraints; i++) {
     ReadRecordConstraint(i);
     double val   = fRecord->GetValue(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;ir<csize;ir++) {
+       int iID = indV[ir];
+       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;
+       }
+       fVecBGlo[iID] += val*der[ir]*sig2i;
+      }
+    }
+    else {   // this is exact constriant:  Lagrange multipliers must be added
+      for (int j=csize; j--;) {
+       int jID = indV[j];
+       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  
+      }
+      //
+      if (matCGlo.QueryDiag(nVar)) matCGlo.DiagElem(nVar) = 0.0;
+      fVecBGlo[nVar++] = float(fNLocEquations)*val; //RS ? should we use here fNLocFits ? 
+      fConstrUsed[i] = kTRUE;
+    }
   }
   //
   AliInfo(Form("Obtained %-7ld equations from %-7ld records (%-7ld rejected). Fixed %-4d globals",
               fNLocEquations,fNLocFits,fNLocFitsRejected,fNGloFix));
 
   //
+  sws.Start();
   fGloSolveStatus = SolveGlobalMatEq();                     // obtain solution for this step
+  sws.Stop();
+  printf("Solve %d |",fIter); sws.Print();
   //
   sw.Stop();
   AliInfo(Form("Iteration#%2d %s. CPU time: %.1f",fIter,fGloSolveStatus==gkFailed ? "Failed":"Converged",sw.CpuTime()));
@@ -827,6 +879,7 @@ Int_t AliMillePede2::GlobalFitIteration()
   //
   for (int i=fNGloPar;i--;) fDeltaPar[i] += fVecBGlo[i];    // Update global parameters values (for iterations)
   //
+  //  PrintGlobalParameters();
   return 1;
 }
 
@@ -836,10 +889,17 @@ Int_t AliMillePede2::SolveGlobalMatEq()
   //
   // solve global matrix equation MatCGlob*X=VecBGlo and store the result in the VecBGlo
   //
+  /*
+  printf("GlobalMatrix\n");
+  fMatCGlo->Print();
+  printf("RHS\n");
+  for (int i=0;i<fNGloPar;i++) printf("%d %+e\n",i,fVecBGlo[i]);
+  */
+  //
   if (!fgIsMatGloSparse) {
     //
-    if (fNGloConstraints==0) { // pos-def systems are faster to solve by Cholesky
-      if ( ((AliSymMatrix*)fMatCGlo)->SolveChol(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 +907,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 +978,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 +1003,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]));
index 9c91225..9d81293 100644 (file)
@@ -1,87 +1,34 @@
 #ifndef ALIMILLEPEDE2_H
 #define ALIMILLEPEDE2_H
 
+/**********************************************************************************************/
+/* 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 <TObject.h>
 #include <TTree.h>
-#include <TFile.h>
-#include <TMatrixD.h>
-#include "AliSymMatrix.h"
-#include "AliMatrixSparse.h"
 #include "AliMinResSolve.h"
+#include "AliMillePedeRecord.h"
+class TFile;
+class AliMatrixSq;
+class AliSymMatrix;
+class AliRectMatrix;
+class AliMatrixSparse;
 class AliLog;
 class TStopwatch;
 
-/**********************************************************************************************/
-/* AliMillePedeRecords: class to store the data of single track processing                    */
-/* Format: for each measured point the data is stored consequtively                           */
-/* INDEX                                                      VALUE                           */
-/* -1                                                         residual                        */
-/* Local_param_id                                             dResidual/dLocal_param          */
-/* ...                                                        ...                             */
-/* -2                                                         weight of the measurement       */
-/* Global_param_od                                            dResidual/dGlobal_param         */
-/* ...                                                        ...                             */
-/*                                                                                            */
-/* The records for all processed tracks are stored in the temporary tree in orgder to be      */
-/* reused for multiple iterations of MillePede                                                */
-/*                                                                                            */
-/**********************************************************************************************/
 
-class AliMillePedeRecord : public TObject
-{
- public:
-  AliMillePedeRecord();
-  AliMillePedeRecord(const AliMillePedeRecord& src);
-  AliMillePedeRecord& operator=(const AliMillePedeRecord& rhs);
-  //
-  virtual    ~AliMillePedeRecord();
-  void       Reset()                                             {fSize = 0;}
-  void       Print(const Option_t *opt="")                 const;
-  //
-  Int_t      GetSize()                                     const {return fSize;}
-  Int_t     *GetIndex()                                    const {return fIndex;}
-  Int_t      GetIndex(int i)                               const {return fIndex[i];}
-  //
-  void       GetIndexValue(int i,int &ind,double &val)     const {ind=fIndex[i]; val=fValue[i];}
-  void       AddIndexValue(int ind, double val);
-  void       AddResidual(double val)                             {AddIndexValue(-1,val);}
-  void       AddWeight(double val)                               {AddIndexValue(-2,val);}
-  Bool_t     IsResidual(int i)                             const {return fIndex[i]==-1;}
-  Bool_t     IsWeight(int i)                               const {return fIndex[i]==-2;}
-  //
-  Double_t  *GetValue()                                    const {return fValue;}
-  Double_t   GetValue(int i)                               const {return fValue[i];}
-  //
- protected:
-  Int_t      GetBufferSize()                               const {return GetUniqueID();}
-  void       SetBufferSize(int sz)                               {SetUniqueID(sz);}
-  void       Expand(int bfsize);
-  //
- protected:
-  Int_t      fSize;                             // size of the record
-  Int_t   *  fIndex;                            //[fSize] index of variables
-  Double_t*  fValue;                            //[fSize] array of values: derivs,residuals
-  //
-  ClassDef(AliMillePedeRecord,1)                // Record of track residuals and local/global deriavtives
-};
-
-inline void  AliMillePedeRecord::AddIndexValue(int ind, double val) 
-{
-  if (fSize>=GetBufferSize()) Expand(2*(fSize+1));
-  fIndex[fSize]=ind; 
-  fValue[fSize++]=val;
-}
-
-//----------------------------------------------------------------------------------------------
-//----------------------------------------------------------------------------------------------
-//----------------------------------------------------------------------------------------------
-//----------------------------------------------------------------------------------------------
 
 class AliMillePede2: public TObject
 {
  public:
   //
-  enum {gkFailed,gkInvert,gkMinRes};    // used global matrix solution methods
+  enum {gkFailed,gkInvert,gkNoInversion};    // used global matrix solution methods
   //
   AliMillePede2();
   AliMillePede2(const AliMillePede2& src);
@@ -97,6 +44,7 @@ class AliMillePede2: public TObject
   Int_t                GetNMaxIterations()              const {return fMaxIter;}
   Int_t                GetNStdDev()                     const {return fNStdDev;} 
   Int_t                GetNGlobalConstraints()          const {return fNGloConstraints;}
+  Int_t                GetNLagrangeConstraints()        const {return fNLagrangeConstraints;}
   Long_t               GetNLocalFits()                  const {return fNLocFits;}
   Long_t               GetNLocalFitsRejected()          const {return fNLocFitsRejected;}
   Int_t                GetNGlobalsFixed()               const {return fNGloFix;}
@@ -108,19 +56,26 @@ class AliMillePede2: public TObject
   Int_t                GetMinPntValid()                 const {return fMinPntValid;}
   Int_t                GetProcessedPoints(Int_t i)      const {return fProcPnt[i];}
   Int_t*               GetProcessedPoints()             const {return fProcPnt;}
+  Int_t                GetParamGrID(Int_t i)            const {return fParamGrID[i];}
   //
   AliMatrixSq*         GetGlobalMatrix()                const {return fMatCGlo;}
+  AliSymMatrix*        GetLocalMatrix()                 const {return fMatCLoc;}
   Double_t*            GetGlobals()                     const {return fVecBGlo;}
   Double_t*            GetDeltaPars()                   const {return fDeltaPar;}
   Double_t*            GetInitPars()                    const {return fInitPar;}
   Double_t*            GetSigmaPars()                   const {return fSigmaPar;}
-  Bool_t*              GetIsLinear()                    const {return fIsLinear;}        
+  Bool_t*              GetIsLinear()                    const {return fIsLinear;}
+  Double_t             GetFinalParam(int i)             const {return fDeltaPar[i]+fInitPar[i];}
+  Double_t             GetFinalError(int i)             const {return GetParError(i);}
+  //
   Double_t             GetGlobal(Int_t i)               const {return fVecBGlo[i];}
   Double_t             GetInitPar(Int_t i)              const {return fInitPar[i];}
   Double_t             GetSigmaPar(Int_t i)             const {return fSigmaPar[i];}
   Bool_t               GetIsLinear(Int_t i)             const {return fIsLinear[i];}
   static Bool_t        IsGlobalMatSparse()                    {return fgIsMatGloSparse;}
+  static Bool_t        IsWeightSigma()                        {return fgWeightSigma;}
   //
+  void                 SetParamGrID(Int_t grID,Int_t i)       {fParamGrID[i] = grID; if(fNGroupsSet<grID)fNGroupsSet=grID;}
   void                 SetNGloPar(Int_t n)                    {fNGloPar = n;}
   void                 SetNLocPar(Int_t n)                    {fNLocPar = n;}
   void                 SetNMaxIterations(Int_t n=10)          {fMaxIter = n;}
@@ -129,8 +84,9 @@ class AliMillePede2: public TObject
   void                 SetChi2CutRef(Float_t v)               {fChi2CutRef = v;}
   void                 SetResCurInit(Float_t v)               {fResCutInit = v;}
   void                 SetResCut(Float_t v)                   {fResCut = v;}
-  void                 SetMinPntValid(Int_t n)                {fMinPntValid = n;}
+  void                 SetMinPntValid(Int_t n)                {fMinPntValid = n>0 ? n:1;}
   static void          SetGlobalMatSparse(Bool_t v=kTRUE)     {fgIsMatGloSparse = v;}
+  static void          SetWeightSigma(Bool_t v=kTRUE)         {fgWeightSigma = v;}
   //
   void                 SetInitPars(const Double_t* par)       {memcpy(fInitPar,par,fNGloPar*sizeof(Double_t));}
   void                 SetSigmaPars(const Double_t* par)      {memcpy(fSigmaPar,par,fNGloPar*sizeof(Double_t));}
@@ -140,12 +96,14 @@ class AliMillePede2: public TObject
   Int_t                GlobalFit(Double_t *par=0, Double_t *error=0, Double_t *pull=0);
   Int_t                GlobalFitIteration();
   Int_t                SolveGlobalMatEq();
+  static void          SetInvChol(Bool_t v=kTRUE)             {fgInvChol = v;}
   static void          SetMinResPrecondType(Int_t tp=0)       {fgMinResCondType = tp;}
   static void          SetMinResTol(Double_t val=1e-12)       {fgMinResTol = val;}
   static void          SetMinResMaxIter(Int_t val=2000)       {fgMinResMaxIter = val;}
   static void          SetIterSolverType(Int_t val=AliMinResSolve::kSolMinRes) {fgIterSol = val;}
   static void          SetNKrylovV(Int_t val=60)              {fgNKrylovV = val;}
   //
+  static Bool_t        GetInvChol()                           {return fgInvChol;}
   static Int_t         GetMinResPrecondType()                 {return fgMinResCondType;}
   static Double_t      GetMinResTol()                         {return fgMinResTol;}
   static Int_t         GetMinResMaxIter()                     {return fgMinResMaxIter;}
@@ -160,8 +118,8 @@ class AliMillePede2: public TObject
 
   //
   // constraints
-  void                 SetGlobalConstraint(double *dergb, double val);
-  void                 SetGlobalConstraint(const int *indgb,double *dergb, int ngb, double val);
+  void                 SetGlobalConstraint(double *dergb, double val, double sigma=0);
+  void                 SetGlobalConstraint(const int *indgb,double *dergb, int ngb, double val, double sigma=0);
   //
   // processing of the local measurement
   void                 SetLocalEquation(double *dergb, double *derlc, double lMeas, double lSigma);
@@ -186,7 +144,6 @@ class AliMillePede2: public TObject
   Bool_t               ReadNextRecordConstraint();
   void                 SaveRecordData();
   void                 SaveRecordConstraint();
-  void                 ResetConstraints();
   AliMillePedeRecord*  GetRecord()                      const {return fRecord;}
   //
   Float_t              Chi2DoFLim(int nSig, int nDoF)   const;
@@ -210,6 +167,7 @@ class AliMillePede2: public TObject
   Int_t                 fMaxIter;                        // Maximum number of iterations
   Int_t                 fNStdDev;                        // Number of standard deviations for chi2 cut 
   Int_t                 fNGloConstraints;                // Number of constraint equations
+  Int_t                 fNLagrangeConstraints;           // Number of constraint equations requiring Lagrange multiplier
   Long_t                fNLocFits;                       // Number of local fits
   Long_t                fNLocFitsRejected;               // Number of local fits rejected
   Int_t                 fNGloFix;                        // Number of globals fixed by user
@@ -221,6 +179,8 @@ class AliMillePede2: public TObject
   Float_t               fResCut;                         // Cut in residual for other iterartiona
   Int_t                 fMinPntValid;                    // min number of points for global to vary
   //
+  Int_t                 fNGroupsSet;                     // number of groups set
+  Int_t                *fParamGrID;                      // group id for the every parameter
   Int_t                *fProcPnt;                        // N of processed points per global variable
   Double_t             *fVecBLoc;                        // Vector B local (parameters) 
   Double_t             *fDiagCGlo;                       // Initial diagonal elements of C global matrix
@@ -239,7 +199,9 @@ class AliMillePede2: public TObject
   // Matrices
   AliSymMatrix         *fMatCLoc;                        // Matrix C local
   AliMatrixSq          *fMatCGlo;                        // Matrix C global
-  TMatrixD             *fMatCGloLoc;                     // Rectangular matrix C g*l 
+  AliRectMatrix        *fMatCGloLoc;                     // Rectangular matrix C g*l 
+  Int_t                *fFillIndex;                      // auxilary index array for fast matrix fill
+  Double_t             *fFillValue;                      // auxilary value array for fast matrix fill
   //
   // processed data record bufferization   
   TString               fDataRecFName;                   // Name of File for data records               
@@ -252,7 +214,10 @@ class AliMillePede2: public TObject
   TFile                *fConsRecFile;                    // File of processed constraints records
   Long_t                fCurrRecDataID;                  // ID of the current data record
   Long_t                fCurrRecConstrID;                // ID of the current constraint record
+  Bool_t                fLocFitAdd;                      // Add contribution of carrent track (and not eliminate it)
   //
+  static Bool_t         fgInvChol;                       // Invert global matrix in Cholesky solver
+  static Bool_t         fgWeightSigma;                   // weight parameter constraint by statistics
   static Bool_t         fgIsMatGloSparse;                // Type of the global matrix (sparse ...)
   static Int_t          fgMinResCondType;                // Type of the preconditioner for MinRes method 
   static Double_t       fgMinResTol;                     // Tolerance for MinRes solution
diff --git a/STEER/AliMillePedeRecord.cxx b/STEER/AliMillePedeRecord.cxx
new file mode 100644 (file)
index 0000000..4cba214
--- /dev/null
@@ -0,0 +1,145 @@
+#include "AliMillePedeRecord.h"
+#include <TMath.h>
+#include "AliLog.h"
+
+/**********************************************************************************************/
+/* AliMillePedeRecords: class to store the data of single track processing                    */
+/* Format: for each measured point the data is stored consequtively                           */
+/* INDEX                                                      VALUE                           */
+/* -1                                                         residual                        */
+/* Local_param_id                                             dResidual/dLocal_param          */
+/* ...                                                        ...                             */
+/* -2                                                         weight of the measurement       */
+/* Global_param_od                                            dResidual/dGlobal_param         */
+/* ...                                                        ...                             */
+/*                                                                                            */
+/* The records for all processed tracks are stored in the temporary tree in orgder to be      */
+/* reused for multiple iterations of MillePede                                                */
+/*                                                                                            */ 
+/* Author: ruben.shahoyan@cern.ch                                                             */
+/*                                                                                            */ 
+/**********************************************************************************************/
+
+ClassImp(AliMillePedeRecord)
+
+//_____________________________________________________________________________________________
+AliMillePedeRecord::AliMillePedeRecord() : 
+fSize(0),fNGroups(0),fGroupID(0),fIndex(0),fValue(0) {SetUniqueID(0);}
+
+//_____________________________________________________________________________________________
+AliMillePedeRecord::AliMillePedeRecord(const AliMillePedeRecord& src) : 
+  TObject(src),fSize(src.fSize),fNGroups(src.fNGroups),fGroupID(0),fIndex(0),fValue(0)
+{
+  fIndex = new Int_t[GetDtBufferSize()];
+  memcpy(fIndex,src.fIndex,fSize*sizeof(Int_t));
+  fValue = new Double_t[GetDtBufferSize()];
+  memcpy(fValue,src.fValue,fSize*sizeof(Double_t));
+  fGroupID = new Int_t[GetGrBufferSize()];
+  memcpy(fGroupID,src.fGroupID,GetGrBufferSize()*sizeof(Int_t));
+}
+
+//_____________________________________________________________________________________________
+AliMillePedeRecord& AliMillePedeRecord::operator=(const AliMillePedeRecord& rhs)
+{
+  if (this!=&rhs) {
+    Reset();
+    for (int i=0;i<rhs.GetSize();i++) {
+      Double_t val;
+      Int_t    ind;
+      rhs.GetIndexValue(i,ind,val);
+      AddIndexValue(ind,val);
+    }
+    for (int i=0;i<rhs.GetNGroups();i++) MarkGroup(rhs.GetGroupID(i));
+  }
+  return *this;
+}
+
+//_____________________________________________________________________________________________
+AliMillePedeRecord::~AliMillePedeRecord() {delete[] fIndex; delete[] fValue; delete[] fGroupID;}
+
+//_____________________________________________________________________________________________
+void AliMillePedeRecord::Reset()
+{
+  fSize = 0;
+  for (int i=fNGroups;i--;) fGroupID[i] = -1;
+  fNGroups = 0;
+}
+
+//_____________________________________________________________________________________________
+void AliMillePedeRecord::Print(const Option_t *) const
+{
+  if (!fSize) {AliInfo("No data"); return;}
+  int cnt=0,point=0;
+  //  
+  if (fNGroups) printf("Groups: ");
+  for (int i=0;i<fNGroups;i++) printf("%4d |",GetGroupID(i)); printf("\n");
+  while(cnt<fSize) {
+    //
+    Double_t resid = fValue[cnt++];
+    Double_t *derLoc = GetValue()+cnt;
+    int    *indLoc = GetIndex()+cnt;
+    int     nLoc = 0;
+    while(!IsWeight(cnt)) {nLoc++;cnt++;}
+    Double_t weight = GetValue(cnt++);
+    Double_t *derGlo = GetValue()+cnt;
+    int    *indGlo = GetIndex()+cnt;
+    int     nGlo = 0;
+    while(!IsResidual(cnt) && cnt<fSize) {nGlo++; cnt++;} 
+    //
+    printf("\n*** Point#%2d | Residual = %+.4e | Weight = %+.4e\n",point++,resid,weight);
+    printf("Locals : "); 
+    for (int i=0;i<nLoc;i++) printf("[%5d] %+.4e|",indLoc[i],derLoc[i]); printf("\n");
+    printf("Globals: "); 
+    for (int i=0;i<nGlo;i++) printf("[%5d] %+.4e|",indGlo[i],derGlo[i]); printf("\n");
+    //
+  }
+  //
+}
+
+//_____________________________________________________________________________________________
+void AliMillePedeRecord::ExpandDtBuffer(Int_t bfsize)
+{
+  // add extra space for derivatives data
+  bfsize = TMath::Max(bfsize,GetDtBufferSize());
+  Int_t *tmpI = new Int_t[bfsize];
+  memcpy(tmpI,fIndex, fSize*sizeof(Int_t));
+  delete fIndex;
+  fIndex = tmpI;
+  //
+  Double_t *tmpD = new Double_t[bfsize];
+  memcpy(tmpD,fValue, fSize*sizeof(Double_t));
+  delete fValue;
+  fValue = tmpD;
+  //
+  SetDtBufferSize(bfsize);
+}
+
+//_____________________________________________________________________________________________
+void AliMillePedeRecord::ExpandGrBuffer(Int_t bfsize)
+{
+  // add extra space for groupID data 
+  bfsize = TMath::Max(bfsize,GetGrBufferSize());
+  Int_t *tmpI = new Int_t[bfsize];
+  memcpy(tmpI,fGroupID, fNGroups*sizeof(Int_t));
+  delete fGroupID;
+  fGroupID = tmpI;
+  for (int i=fNGroups;i<bfsize;i++) fGroupID[i] = -1;
+  //
+  SetGrBufferSize(bfsize);
+}
+
+//_____________________________________________________________________________________________
+void AliMillePedeRecord::MarkGroup(Int_t id)
+{
+  // mark the presence of the detector group
+  if (fNGroups>0 && fGroupID[fNGroups-1]==id) return; // already there
+  if (fNGroups>=GetGrBufferSize()) ExpandGrBuffer(2*(fNGroups+1));
+  fGroupID[fNGroups++] = id;  
+}
+
+//_____________________________________________________________________________________________
+Bool_t AliMillePedeRecord::IsGroupPresent(Int_t id) const
+{
+  for (int i=fNGroups;i--;) if (GetGroupID(i)==id) return kTRUE;
+  return kFALSE;
+}
diff --git a/STEER/AliMillePedeRecord.h b/STEER/AliMillePedeRecord.h
new file mode 100644 (file)
index 0000000..2dee67f
--- /dev/null
@@ -0,0 +1,78 @@
+#ifndef ALIMILLEPEDERECORD_H
+#define ALIMILLEPEDERECORD_H
+
+/**********************************************************************************************/
+/* AliMillePedeRecords: class to store the data of single track processing                    */
+/* Format: for each measured point the data is stored consequtively                           */
+/* INDEX                                                      VALUE                           */
+/* -1                                                         residual                        */
+/* Local_param_id                                             dResidual/dLocal_param          */
+/* ...                                                        ...                             */
+/* -2                                                         weight of the measurement       */
+/* Global_param_od                                            dResidual/dGlobal_param         */
+/* ...                                                        ...                             */
+/*                                                                                            */
+/* The records for all processed tracks are stored in the temporary tree in orgder to be      */
+/* reused for multiple iterations of MillePede                                                */
+/*                                                                                            */
+/* Author: ruben.shahoyan@cern.ch                                                             */
+/*                                                                                            */
+/**********************************************************************************************/
+#include <TObject.h>
+
+class AliMillePedeRecord : public TObject
+{
+ public:
+  AliMillePedeRecord();
+  AliMillePedeRecord(const AliMillePedeRecord& src);
+  AliMillePedeRecord& operator=(const AliMillePedeRecord& rhs);
+  //
+  virtual    ~AliMillePedeRecord();
+  void       Reset();
+  void       Print(const Option_t *opt="")                   const;
+  //
+  Int_t      GetSize()                                       const {return fSize;}
+  Int_t     *GetIndex()                                      const {return fIndex;}
+  Int_t      GetIndex(int i)                                 const {return fIndex[i];}
+  //
+  void       GetIndexValue(Int_t i,Int_t &ind,Double_t &val) const {ind=fIndex[i]; val=fValue[i];}
+  void       AddIndexValue(Int_t ind, Double_t val);
+  void       AddResidual(Double_t val)                             {AddIndexValue(-1,val);}
+  void       AddWeight(Double_t val)                               {AddIndexValue(-2,val);}
+  Bool_t     IsResidual(Int_t i)                             const {return fIndex[i]==-1;}
+  Bool_t     IsWeight(Int_t i)                               const {return fIndex[i]==-2;}
+  //
+  Double_t  *GetValue()                                      const {return fValue;}
+  Double_t   GetValue(Int_t i)                               const {return fValue[i];}
+  //
+  void       MarkGroup(Int_t id);
+  Int_t      GetNGroups()                                    const {return fNGroups;}
+  Int_t      GetGroupID(Int_t i)                             const {return fGroupID[i];}
+  Bool_t     IsGroupPresent(Int_t id)                        const;
+  //
+ protected:
+  Int_t      GetDtBufferSize()                               const {return GetUniqueID()&0x0000ffff;}
+  Int_t      GetGrBufferSize()                               const {return GetUniqueID()>>16;}
+  void       SetDtBufferSize(Int_t sz)                             {SetUniqueID((GetGrBufferSize()<<16)+sz);}
+  void       SetGrBufferSize(Int_t sz)                             {SetUniqueID(GetDtBufferSize()+(sz<<16));}
+  void       ExpandDtBuffer(Int_t bfsize);
+  void       ExpandGrBuffer(Int_t bfsize);
+  //
+ protected:
+  Int_t      fSize;                             // size of the record
+  Int_t      fNGroups;                          // number of groups (e.g. detectors) contributing
+  Int_t   *  fGroupID;                          //[fNGroups] groups id's (in increasing order)
+  Int_t   *  fIndex;                            //[fSize] index of variables
+  Double_t*  fValue;                            //[fSize] array of values: derivs,residuals
+  //
+  ClassDef(AliMillePedeRecord,1)                // Record of track residuals and local/global deriavtives
+};
+
+inline void  AliMillePedeRecord::AddIndexValue(Int_t ind, Double_t val) 
+{
+  if (fSize>=GetDtBufferSize()) ExpandDtBuffer(2*(fSize+1));
+  fIndex[fSize]=ind; 
+  fValue[fSize++]=val;
+}
+
+#endif
index bf821be..54d47f0 100644 (file)
@@ -70,12 +70,12 @@ Int_t AliMinResSolve::BuildPrecon(Int_t prec)
     fPrecAux  = new Double_t[ fSize ];
     //
     // copy inverse diagonal
-    for (int i=0;i<fSize;i++) fPrecDiag[i] = fMatrix->QuerryDiag(i);
+    for (int i=0;i<fSize;i++) fPrecDiag[i] = fMatrix->QueryDiag(i);
     //
     for (int i=0;i<fSize;i++) {
       fPrecDiag[i] = TMath::Abs(fPrecDiag[i])>kTiny  ? 1./fPrecDiag[i] : 1./TMath::Sign(kTiny,fPrecDiag[i]);
       for (int j=i+1;j<fSize;j++) {
-       double vl = fMatrix->Querry(j,i);
+       double vl = fMatrix->Query(j,i);
        if (vl!=0) fPrecDiag[j] -= fPrecDiag[i]*vl*vl; 
       }
     }
@@ -311,7 +311,7 @@ Bool_t AliMinResSolve::SolveMinRes(double *VecSol,Int_t precon,int itnlim,double
   //   Set up y and v for the first Lanczos vector v1.
   //   y  =  beta1 P' v1,  where  P = C**(-1). v is really P' v1.
   //
-  for (int i=fSize;i--;) fPVecY[i] = fPVecR1[i] = fRHS[i];
+  for (int i=fSize;i--;) fPVecY[i]  = fPVecR1[i] = fRHS[i];
   //
   if ( precon>0 ) ApplyPrecon( fRHS, fPVecY);
   beta1 = 0; for (int i=fSize;i--;) beta1 += fRHS[i]*fPVecY[i]; //
@@ -508,7 +508,7 @@ void AliMinResSolve::ApplyPrecon(const double* vecRHS, double* vecOut) const
   const Double_t kTiny = 1E-12;
   if (fPrecon==kPreconDiag) {
     for (int i=fSize;i--;) {
-      double d = fMatrix->QuerryDiag(i);
+      double d = fMatrix->QueryDiag(i);
       vecOut[i] =  vecRHS[i]/(TMath::Abs(d)>kTiny ? d : kTiny);
     }
     //    return;
@@ -517,12 +517,12 @@ void AliMinResSolve::ApplyPrecon(const double* vecRHS, double* vecOut) const
   else if (fPrecon==kPreconDILU) {
     for (int i=0;i<fSize;i++) {
       double el = 0;
-      for (int j=i;j--;) {double vl = fMatrix->Querry(i,j);if (vl!=0) el += fPrecAux[j]*vl;}
+      for (int j=i;j--;) {double vl = fMatrix->Query(i,j);if (vl!=0) el += fPrecAux[j]*vl;}
       fPrecAux[i] = fPrecDiag[i]*(vecRHS[i]-el); 
     }
     for (int i=fSize;i--;) {
       double el = 0;
-      for (int j=i+1;j<fSize;j++) {double vl = fMatrix->Querry(i,j);if (vl!=0) el += vl*vecOut[j];}
+      for (int j=i+1;j<fSize;j++) {double vl = fMatrix->Query(i,j);if (vl!=0) el += vl*vecOut[j];}
       vecOut[i] = fPrecAux[i] - fPrecDiag[i]*el;
     }
     //    return;
@@ -646,6 +646,12 @@ Int_t  AliMinResSolve::BuildPreconILUK(Int_t lofM)
   for(int j=fSize;j--;) jw[j] = -1;     // set indicator array jw to -1 
   //
   for(int i=0; i<fSize; i++ ) {            // beginning of main loop
+    if ( (i%int(0.1*fSize)) == 0) {
+      printf("BuildPrecon: row %d of %d\n",i,fSize);
+      sw.Stop();
+      sw.Print();
+      sw.Start(kFALSE);
+    }
     /* setup array jw[], and initial i-th row */
     AliVectorSparse& rowLi = *fMatL->GetRow(i);
     AliVectorSparse& rowUi = *fMatU->GetRow(i);
@@ -673,7 +679,7 @@ Int_t  AliMinResSolve::BuildPreconILUK(Int_t lofM)
       else rowUi.GetElem(jw[col]) = rowM.GetElem(j);
     }
     if (fMatrix->IsSymmetric()) for (int col=i+1;col<fSize;col++) {      // part of the row I on the right of diagonal is stored as 
-       double vl = fMatrix->Querry(col,i);    // the lower part of the column I
+       double vl = fMatrix->Query(col,i);    // the lower part of the column I
        if (vl==0) continue;                   
        rowUi.GetElem(jw[col]) = vl;
       }
@@ -764,7 +770,7 @@ Int_t  AliMinResSolve::BuildPreconILUKDense(Int_t lofM)
     }
     // copy row from csmat into L,U D
     for(int j=fSize; j--;) {  // L and D part 
-      double vl = fMatrix->Querry(i,j);
+      double vl = fMatrix->Query(i,j);
       if (vl==0) continue;
       if( j < i )   rowLi.GetElem(jw[j]) = vl;
       else if(j==i) fPrecDiag[i] = vl;
@@ -818,6 +824,9 @@ Int_t  AliMinResSolve::PreconILUKsymb(Int_t lofM)
    * Adapted from iluk.c: lofC of ITSOL_1 package by Y.Saad: http://www-users.cs.umn.edu/~saad/software/
    *----------------------------------------------------------------------------*/
   //
+  TStopwatch sw;
+  printf("PreconILUKsymb>>\n");
+  sw.Start();
   AliMatrixSparse* Matrix = (AliMatrixSparse*)fMatrix;
   //
   UChar_t **ulvl=0,*levls=0;
@@ -862,7 +871,7 @@ Int_t  AliMinResSolve::PreconILUKsymb(Int_t lofM)
       }
     }
     if (Matrix->IsSymmetric()) for (int col=i+1;col<fSize;col++) { // U-part of symmetric matrix
-       if (fMatrix->Querry(col,i)==0) continue;    // Due to the symmetry  == matrix(i,col)
+       if (fMatrix->Query(col,i)==0) continue;    // Due to the symmetry  == matrix(i,col)
        jbuf[incu] = col;
        levls[incu] = 0;
        iw[col] = incu++;
@@ -938,6 +947,9 @@ Int_t  AliMinResSolve::PreconILUKsymb(Int_t lofM)
   //
   fMatL->SortIndices();
   fMatU->SortIndices();
+  sw.Stop();
+  sw.Print();
+  printf("PreconILUKsymb<<\n");
   return 0;
 }
 
@@ -978,7 +990,7 @@ Int_t  AliMinResSolve::PreconILUKsymbDense(Int_t lofM)
     //
     // assign lof = 0 for matrix elements
     for(int j=0;j<fSize; j++) {
-      if (fMatrix->Querry(i,j)==0) continue;
+      if (fMatrix->Query(i,j)==0) continue;
       if (j<i) {                      // L-part
        jbuf[incl] = j;
        levls[incl] = 0;
diff --git a/STEER/AliRectMatrix.cxx b/STEER/AliRectMatrix.cxx
new file mode 100644 (file)
index 0000000..3490912
--- /dev/null
@@ -0,0 +1,81 @@
+#include "AliRectMatrix.h"
+#include <TString.h>
+//
+
+ClassImp(AliRectMatrix)
+
+
+//___________________________________________________________
+AliRectMatrix::AliRectMatrix() 
+: fNRows(0),fNCols(0),fRows(0)
+{}
+
+//___________________________________________________________
+AliRectMatrix::AliRectMatrix(Int_t nrow,Int_t ncol)
+  : fNRows(nrow),fNCols(ncol),fRows(0)
+{
+  //
+  fRows = new Double_t*[fNRows];
+  for (int i=fNRows;i--;) {
+    fRows[i] = new Double_t[fNCols];
+    memset(fRows[i],0,fNCols*sizeof(Double_t));
+  }
+  //
+}
+
+//___________________________________________________________
+AliRectMatrix::AliRectMatrix(const AliRectMatrix &src)
+  : TObject(src),fNRows(src.fNRows), fNCols(src.fNCols), fRows(0)
+{
+  for (int i=fNRows;i--;) {
+    fRows[i] = new Double_t[fNCols];
+    memcpy(fRows[i], src.fRows[i], fNCols*sizeof(Double_t));
+  }
+}
+
+//___________________________________________________________
+AliRectMatrix::~AliRectMatrix()
+{
+  if (fNRows) for (int i=fNRows;i--;) delete[] fRows[i];
+  delete[] fRows;
+}
+
+//___________________________________________________________
+AliRectMatrix& AliRectMatrix::operator=(const AliRectMatrix& src)
+{
+  //
+  if (&src == this) return *this;
+  if (fNRows) for (int i=fNRows;i--;) delete[] fRows[i];
+  delete[] fRows;
+  fNRows = src.fNRows;
+  fNCols = src.fNCols;
+  fRows = new Double_t*[fNRows];
+  for (int i=fNRows;i--;) {
+    fRows[i] = new Double_t[fNCols];
+    memcpy(fRows[i], src.fRows[i], fNCols*sizeof(Double_t));
+  }
+  //
+  return *this;
+}
+
+//___________________________________________________________
+void AliRectMatrix::Print(Option_t* option) const
+{
+  printf("Rectangular Matrix:  %d rows %d columns\n",fNRows,fNCols);
+  TString opt = option; opt.ToLower();
+  if (opt.IsNull()) return;
+  for (int i=0;i<fNRows;i++) {
+    for (Int_t j=0;j<=fNCols;j++) printf("%+.3e|",Query(i,j));
+    printf("\n");
+  }
+}
+
+
+//___________________________________________________________
+void AliRectMatrix::Reset()
+{
+  for (int i=fNRows;i--;) {
+    double *row = GetRow(i); 
+    for (int j=fNCols;j--;)  row[j] = 0.;
+  }
+}
diff --git a/STEER/AliRectMatrix.h b/STEER/AliRectMatrix.h
new file mode 100644 (file)
index 0000000..0885b01
--- /dev/null
@@ -0,0 +1,52 @@
+#ifndef ALIRECTMATRIX_H
+#define ALIRECTMATRIX_H
+
+#include "TObject.h"
+class TString;
+
+class AliRectMatrix : public TObject {
+  //
+ public:
+  AliRectMatrix();
+  AliRectMatrix(Int_t nrow,Int_t ncol);
+  AliRectMatrix(const AliRectMatrix &src);
+  virtual ~AliRectMatrix();
+  //
+  Int_t         GetNRows()                            const {return fNRows;}
+  Int_t         GetNCols()                            const {return fNCols;}
+  //
+  Double_t      Query(Int_t rown, Int_t coln)         const {return operator()(rown,coln);}
+
+  AliRectMatrix& operator=(const AliRectMatrix& src);
+  Double_t      operator()(Int_t rown, Int_t coln)    const;
+  Double_t&     operator()(Int_t rown, Int_t coln);
+  Double_t*     operator()(Int_t row)                 const {return GetRow(row);}
+  Double_t*     GetRow(Int_t row)                     const {return fRows[row];}
+  //
+  void          Reset();
+  //
+  virtual void  Print(Option_t* option="")           const;
+  //
+ protected:
+  //
+  Int_t   fNRows;       // Number of rows
+  Int_t   fNCols;       // Number of columns
+  Double_t **fRows;     // pointers on rows
+  //
+  ClassDef(AliRectMatrix,0) //Rectangular Matrix Class
+};
+
+//___________________________________________________________
+inline Double_t AliRectMatrix::operator()(Int_t row, Int_t col) const
+{
+  return (const Double_t&) GetRow(row)[col];
+}
+
+//___________________________________________________________
+inline Double_t& AliRectMatrix::operator()(Int_t row, Int_t col)
+{
+  return (Double_t&) fRows[row][col];
+}
+
+
+#endif
index 9152a88..50bed50 100644 (file)
@@ -82,7 +82,8 @@ AliSymMatrix&  AliSymMatrix::operator=(const AliSymMatrix& src)
       for (int i=0;i<fNrows;i++) delete[] fElemsAdd[i]; 
       delete[] fElemsAdd;
       //
-      fNrowIndex = fNcols = src.GetSize();
+      fNrowIndex = src.GetSize(); 
+      fNcols = src.GetSize();
       fNrows = 0;
       fElems     = new Double_t[fNcols*(fNcols+1)/2];
       int nmainel = src.fNcols*(src.fNcols+1);
@@ -121,7 +122,9 @@ void AliSymMatrix::Clear(Option_t*)
     delete[] fElemsAdd;
     fElemsAdd = 0;
   }
-  fNrowIndex = fNcols = fNrows = 0;
+  fNrowIndex = 0;
+  fNcols = 0;
+  fNrows = 0;
   //
 }
 
@@ -164,12 +167,12 @@ void AliSymMatrix::MultiplyByVec(Double_t *vecIn,Double_t *vecOut) const
 AliSymMatrix* AliSymMatrix::DecomposeChol() 
 {
   // Return a matrix with Choleski decomposition
-  /*Adopted from Numerical Recipes in C, ch.2-9, http://www.nr.com
-    consturcts Cholesky decomposition of SYMMETRIC and
-    POSITIVELY-DEFINED matrix a (a=L*Lt)
-    Only upper triangle of the matrix has to be filled.
-    In opposite to function from the book, the matrix is modified:
-    lower triangle and diagonal are refilled. */
+  /Adopted from Numerical Recipes in C, ch.2-9, http://www.nr.com
+  // consturcts Cholesky decomposition of SYMMETRIC and
+  // POSITIVELY-DEFINED matrix a (a=L*Lt)
+  // Only upper triangle of the matrix has to be filled.
+  // In opposite to function from the book, the matrix is modified:
+  // lower triangle and diagonal are refilled.
   //
   if (!fgBuffer || fgBuffer->GetSize()!=GetSize()) {
     delete fgBuffer; 
@@ -186,18 +189,20 @@ AliSymMatrix* AliSymMatrix::DecomposeChol()
   AliSymMatrix& mchol = *fgBuffer;
   //
   for (int i=0;i<fNrowIndex;i++) {
+    Double_t *rowi = mchol.GetRow(i);
     for (int j=i;j<fNrowIndex;j++) {
-      double sum=(*this)(i,j);
-      for (int k=i-1;k>=0;k--) sum -= mchol(i,k)*mchol(j,k);
+      Double_t *rowj = mchol.GetRow(j);
+      double sum = rowj[i];
+      for (int k=i-1;k>=0;k--) if (rowi[k]&&rowj[k]) sum -= rowi[k]*rowj[k];
       if (i == j) {
        if (sum <= 0.0) { // not positive-definite
          printf("The matrix is not positive definite [%e]\n"
                 "Choleski decomposition is not possible\n",sum);
          return 0;
        }
-       mchol(i,i) = TMath::Sqrt(sum);
+       rowi[i] = TMath::Sqrt(sum);
        //
-      } else mchol(j,i) = sum/mchol(i,i);
+      } else rowj[i] = sum/rowi[i];
     }
   }
   return fgBuffer;
@@ -218,6 +223,7 @@ Bool_t AliSymMatrix::InvertChol()
   return kTRUE;
   //
 }
+
 //___________________________________________________________
 void AliSymMatrix::InvertChol(AliSymMatrix* pmchol) 
 {
@@ -230,9 +236,12 @@ void AliSymMatrix::InvertChol(AliSymMatrix* pmchol)
   for (int i=0;i<fNrowIndex;i++) { 
     mchol(i,i) =  1.0/mchol(i,i);
     for (int j=i+1;j<fNrowIndex;j++) { 
+      Double_t *rowj = mchol.GetRow(j);
       sum = 0.0; 
-      for (int k=i;k<j;k++) sum -= mchol(j,k)*mchol(k,i);
-      mchol(j,i) = sum/mchol(j,j);
+      for (int k=i;k<j;k++) if (rowj[k]) { 
+       double &mki = mchol(k,i); if (mki) sum -= rowj[k]*mki;
+      }
+      rowj[i] = sum/rowj[j];
     } 
   }
   //
@@ -240,7 +249,13 @@ void AliSymMatrix::InvertChol(AliSymMatrix* pmchol)
   for (int i=fNrowIndex;i--;) {
     for (int j=i+1;j--;) {
       sum = 0;
-      for (int k=i;k<fNrowIndex;k++) sum += mchol(i,k)*mchol(j,k);
+      for (int k=i;k<fNrowIndex;k++) {
+       double &mik = mchol(i,k); 
+       if (mik) {
+         double &mjk = mchol(j,k);
+         if (mjk) sum += mik*mjk;
+       }
+      }
       (*this)(j,i) = sum;
     }
   }
@@ -251,12 +266,12 @@ void AliSymMatrix::InvertChol(AliSymMatrix* pmchol)
 //___________________________________________________________
 Bool_t AliSymMatrix::SolveChol(Double_t *b, Bool_t invert) 
 {
-  /* Adopted from Numerical Recipes in C, ch.2-9, http://www.nr.com
-     Solves the set of n linear equations A x = b, 
-     where a is a positive-definite symmetric matrix. 
-     a[1..n][1..n] is the output of the routine CholDecomposw. 
-     Only the lower triangle of a is accessed. b[1..n] is input as the 
-     right-hand side vector. The solution vector is returned in b[1..n].*/
+  // Adopted from Numerical Recipes in C, ch.2-9, http://www.nr.com
+  // Solves the set of n linear equations A x = b, 
+  // where a is a positive-definite symmetric matrix. 
+  // a[1..n][1..n] is the output of the routine CholDecomposw. 
+  // Only the lower triangle of a is accessed. b[1..n] is input as the 
+  // right-hand side vector. The solution vector is returned in b[1..n].
   //
   Int_t i,k;
   Double_t sum;
@@ -270,11 +285,15 @@ Bool_t AliSymMatrix::SolveChol(Double_t *b, Bool_t invert)
   AliSymMatrix& mchol = *pmchol;
   //
   for (i=0;i<fNrowIndex;i++) {
-    for (sum=b[i],k=i-1;k>=0;k--) sum -= mchol(i,k)*b[k];
-    b[i]=sum/mchol(i,i);
+    Double_t *rowi = mchol.GetRow(i);
+    for (sum=b[i],k=i-1;k>=0;k--) if (rowi[k]&&b[k]) sum -= rowi[k]*b[k];
+    b[i]=sum/rowi[i];
   }
+  //
   for (i=fNrowIndex-1;i>=0;i--) {
-    for (sum=b[i],k=i+1;k<fNrowIndex;k++) sum -= mchol(k,i)*b[k];
+    for (sum=b[i],k=i+1;k<fNrowIndex;k++) if (b[k]) {
+      double &mki=mchol(k,i); if (mki) sum -= mki*b[k];
+    }
     b[i]=sum/mchol(i,i);
   }
   //
@@ -338,6 +357,51 @@ void AliSymMatrix::Reset()
   //
 }
 
+//___________________________________________________________
+/*
+void AliSymMatrix::AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n)
+{
+  //   for (int i=n;i--;) {
+  //     (*this)(indc[i],r) += valc[i];
+  //   }
+  //   return;
+
+  double *row;
+  if (r>=fNrowIndex) {
+    AddRows(r-fNrowIndex+1); 
+    row = &((fElemsAdd[r-fNcols])[0]);
+  }
+  else row = &fElems[GetIndex(r,0)];
+  //
+  int nadd = 0;
+  for (int i=n;i--;) {
+    if (indc[i]>r) continue;
+    row[indc[i]] += valc[i];
+    nadd++;
+  }
+  if (nadd == n) return;
+  //
+  // add to col>row
+  for (int i=n;i--;) {
+    if (indc[i]>r) (*this)(indc[i],r) += valc[i];
+  }
+  //
+}
+*/
+
+//___________________________________________________________
+Double_t* AliSymMatrix::GetRow(Int_t r)
+{
+  if (r>=fNrowIndex) {
+    int nn = fNrowIndex;
+    AddRows(r-fNrowIndex+1); 
+    printf("create %d of %d\n",r, nn);
+    return &((fElemsAdd[r-fNcols])[0]);
+  }
+  else return &fElems[GetIndex(r,0)];
+}
+
+
 //___________________________________________________________
 int AliSymMatrix::SolveSpmInv(double *vecB, Bool_t stabilize)
 {
@@ -359,7 +423,7 @@ int AliSymMatrix::SolveSpmInv(double *vecB, Bool_t stabilize)
     colMax   = new double[nGlo];
     for (Int_t i=nGlo; i--;) rowMax[i] = colMax[i] = 0.0;
     for (Int_t i=nGlo; i--;) for (Int_t j=i+1;j--;) { 
-       double vl = TMath::Abs(Querry(i,j));
+       double vl = TMath::Abs(Query(i,j));
        if (vl==0) continue;
        if (vl > rowMax[i]) rowMax[i] = vl; // Max elemt of row i
        if (vl > colMax[j]) colMax[j] = vl; // Max elemt of column j
@@ -391,16 +455,16 @@ int AliSymMatrix::SolveSpmInv(double *vecB, Bool_t stabilize)
   //
   if (stabilize) for (int i=0;i<nGlo; i++) { // Small loop for matrix equilibration (gives a better conditioning) 
       for (int j=0;j<=i; j++) {
-       double vl = Querry(i,j);
+       double vl = Query(i,j);
        if (vl!=0) SetEl(i,j, TMath::Sqrt(rowMax[i])*vl*TMath::Sqrt(colMax[j]) ); // Equilibrate the V matrix
       }
       for (int j=i+1;j<nGlo;j++) {
-       double vl = Querry(j,i);
+       double vl = Query(j,i);
        if (vl!=0) fgBuffer->SetEl(j,i,TMath::Sqrt(rowMax[i])*vl*TMath::Sqrt(colMax[j]) ); // Equilibrate the V matrix
       }
     }
   //
-  for (Int_t j=nGlo; j--;) fgBuffer->DiagElem(j) = TMath::Abs(QuerryDiag(j)); // save diagonal elem absolute values 
+  for (Int_t j=nGlo; j--;) fgBuffer->DiagElem(j) = TMath::Abs(QueryDiag(j)); // save diagonal elem absolute values 
   //
   for (Int_t i=0; i<nGlo; i++) {
     vPivot = 0.0;
@@ -408,7 +472,7 @@ int AliSymMatrix::SolveSpmInv(double *vecB, Bool_t stabilize)
     //
     for (Int_t j=0; j<nGlo; j++) { // First look for the pivot, ie max unused diagonal element       
       double vl;
-      if (bUnUsed[j] && (TMath::Abs(vl=QuerryDiag(j))>TMath::Max(TMath::Abs(vPivot),eps*fgBuffer->QuerryDiag(j)))) {    
+      if (bUnUsed[j] && (TMath::Abs(vl=QueryDiag(j))>TMath::Max(TMath::Abs(vPivot),eps*fgBuffer->QueryDiag(j)))) {    
        vPivot = vl;
        iPivot = j;
       }
@@ -424,8 +488,8 @@ int AliSymMatrix::SolveSpmInv(double *vecB, Bool_t stabilize)
        for (Int_t jj=0; jj<nGlo; jj++) {  
          if (j != iPivot && jj != iPivot) {// Other elements (!!! do them first as you use old matV[k][j]'s !!!)         
            double &r = j>=jj ? (*this)(j,jj) : (*fgBuffer)(jj,j);
-           r -= vPivot* ( j>iPivot  ? Querry(j,iPivot)  : fgBuffer->Querry(iPivot,j) )
-             *          ( iPivot>jj ? Querry(iPivot,jj) : fgBuffer->Querry(jj,iPivot));
+           r -= vPivot* ( j>iPivot  ? Query(j,iPivot)  : fgBuffer->Query(iPivot,j) )
+             *          ( iPivot>jj ? Query(iPivot,jj) : fgBuffer->Query(jj,iPivot));
          }
        }
       }
@@ -460,8 +524,8 @@ int AliSymMatrix::SolveSpmInv(double *vecB, Bool_t stabilize)
     rowMax[j] = 0.0;
     for (Int_t jj=0; jj<nGlo; jj++) { // Reverse matrix elements
       double vl;
-      if (j>=jj) vl = (*this)(j,jj)     = -Querry(j,jj);
-      else       vl = (*fgBuffer)(j,jj) = -fgBuffer->Querry(j,jj);
+      if (j>=jj) vl = (*this)(j,jj)     = -Query(j,jj);
+      else       vl = (*fgBuffer)(j,jj) = -fgBuffer->Query(j,jj);
       rowMax[j] += vl*vecB[jj];
     }          
   }
index b154e32..f22d5d5 100644 (file)
@@ -28,12 +28,15 @@ class AliSymMatrix : public AliMatrixSq {
   Double_t      DiagElem(Int_t r)                                const {return (*(const AliSymMatrix*)this)(r,r);}
   Double_t&     DiagElem(Int_t r)                                      {return (*this)(r,r);}
   //
+  Double_t*     GetRow(Int_t r);
+  //
   void          Print(Option_t* option="")                       const;
   void          AddRows(int nrows=1);
   //
   void          Scale(Double_t coeff);
   void          MultiplyByVec(Double_t* vecIn, Double_t* vecOut) const;
   void          MultiplyByVec(TVectorD &vecIn, TVectorD &vecOut) const;
+  void          AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n);
   //
   // ---------------------------------- Dummy methods of MatrixBase
   virtual       const Double_t   *GetMatrixArray  () const {return fElems;};
@@ -111,5 +114,10 @@ inline void AliSymMatrix::Scale(Double_t coeff)
   for (int i=fNrowIndex;i--;) for (int j=i;j--;) { double& el = operator()(i,j); if (el) el *= coeff;}
 }
 
+//___________________________________________________________
+inline void AliSymMatrix::AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n)
+{
+  for (int i=n;i--;) (*this)(indc[i],r) += valc[i];
+}
 
 #endif
diff --git a/STEER/AliVectorSparse.cxx b/STEER/AliVectorSparse.cxx
new file mode 100644 (file)
index 0000000..2fc86f1
--- /dev/null
@@ -0,0 +1,212 @@
+#include "AliVectorSparse.h"
+
+/**********************************************************************************************/
+/* Sparse vector class, used as row of the AliMatrixSparse class                              */
+/*                                                                                            */ 
+/* Author: ruben.shahoyan@cern.ch                                                             */
+/*                                                                                            */ 
+/**********************************************************************************************/
+
+ClassImp(AliVectorSparse)
+
+//___________________________________________________________
+AliVectorSparse::AliVectorSparse()
+  : fNElems(0),fIndex(0),fElems(0) {}
+
+//___________________________________________________________
+AliVectorSparse::AliVectorSparse(const AliVectorSparse& src)
+  : TObject(src),fNElems(src.fNElems),fIndex(0),fElems(0)
+{
+  fIndex = new UShort_t[fNElems];
+  fElems = new Double_t[fNElems];
+  memcpy(fIndex,src.fIndex,fNElems*sizeof(UShort_t));
+  memcpy(fElems,src.fElems,fNElems*sizeof(Double_t));
+}
+
+//___________________________________________________________
+void AliVectorSparse::Clear(Option_t*)
+{
+  delete[] fIndex; fIndex = 0;
+  delete[] fElems; fElems = 0;
+  fNElems = 0;
+}
+
+//___________________________________________________________
+AliVectorSparse& AliVectorSparse::operator=(const AliVectorSparse& src)
+{
+  if (&src==this) return *this;
+  Clear();
+  TObject::operator=(src);
+  fNElems = src.fNElems;
+  fIndex = new UShort_t[fNElems];
+  fElems = new Double_t[fNElems];
+  memcpy(fIndex,src.fIndex,fNElems*sizeof(UShort_t));
+  memcpy(fElems,src.fElems,fNElems*sizeof(Double_t));
+  //
+  return *this;
+}
+
+//___________________________________________________________
+Double_t AliVectorSparse::FindIndex(Int_t ind) const
+{
+  // return an element with given index
+  //printf("V: findindex\n");
+  int first = 0;
+  int last = fNElems-1;
+  while (first<=last) {
+    int mid = (first+last)>>1;
+    if (ind>fIndex[mid]) first = mid+1;
+    else if (ind<fIndex[mid]) last = mid-1;
+    else return fElems[mid];
+  }
+  return 0.0;
+}
+
+//___________________________________________________________
+void AliVectorSparse::SetToZero(Int_t ind)
+{
+  // set element to 0 if it was already defined
+  int first = 0;
+  int last = fNElems-1;
+  while (first<=last) {
+    int mid = (first+last)>>1;
+    if (ind>fIndex[mid]) first = mid+1;
+    else if (ind<fIndex[mid]) last = mid-1;
+    else {fElems[mid] = 0.; return;}
+  }
+}
+
+//___________________________________________________________
+Double_t& AliVectorSparse::FindIndexAdd(Int_t ind)
+{
+  // increment an element with given index
+  //printf("V: findindexAdd\n");
+  int first = 0;
+  int last = fNElems-1;
+  while (first<=last) {
+    int mid = (first+last)>>1;
+    if (ind>fIndex[mid]) first = mid+1;
+    else if (ind<fIndex[mid]) last = mid-1;
+    else return fElems[mid];
+  }
+  // need to insert a new element
+  UShort_t *arrI = new UShort_t[fNElems+1];
+  memcpy(arrI,fIndex,first*sizeof(UShort_t));
+  arrI[first] = ind;
+  memcpy(arrI+first+1,fIndex+first,(fNElems-first)*sizeof(UShort_t));
+  delete[] fIndex;
+  fIndex = arrI;
+  //
+  Double_t   *arrE = new Double_t[fNElems+1];
+  memcpy(arrE,fElems,first*sizeof(Double_t));
+  arrE[first] = 0;
+  memcpy(arrE+first+1,fElems+first,(fNElems-first)*sizeof(Double_t));
+  delete[] fElems;
+  fElems = arrE;
+  //
+  fNElems++;
+  return fElems[first];
+  //
+}
+
+//__________________________________________________________
+void AliVectorSparse::ReSize(Int_t sz,Bool_t copy)
+{
+  if (sz<1) {Clear(); return;}
+    // need to insert a new element
+  UShort_t *arrI = new UShort_t[sz];
+  Double_t *arrE = new Double_t[sz];
+  memset(arrI,0,sz*sizeof(UShort_t));
+  memset(arrE,0,sz*sizeof(Double_t));
+  //
+  if (copy && fIndex) {
+    int cpsz = TMath::Min(fNElems,sz);
+    memcpy(arrI,fIndex,cpsz*sizeof(UShort_t));
+    memcpy(arrE,fElems,cpsz*sizeof(Double_t));
+  }
+  delete[] fIndex;
+  delete[] fElems;
+  fIndex = arrI;
+  fElems = arrE;
+  fNElems = sz;
+  //
+}
+
+//__________________________________________________________
+void AliVectorSparse::SortIndices(Bool_t valuesToo)
+{
+  // sort indices in increasing order. Used to fix the row after ILUk decomposition
+  for (int i=fNElems;i--;) for (int j=i;j--;) if (fIndex[i]<fIndex[j]) { //swap
+       UShort_t tmpI = fIndex[i]; fIndex[i] = fIndex[j]; fIndex[j]=tmpI;
+       if (valuesToo) {Double_t tmpV = fElems[i];fElems[i]=fElems[j];fElems[j]=tmpV;}
+      }
+}
+
+//__________________________________________________________
+void AliVectorSparse::Print(Option_t* )  const
+{
+  printf("|");
+  for (int i=0;i<fNElems;i++) printf("%2d:%+.2e|",fIndex[i],fElems[i]);
+  printf("|\n");
+}
+
+//___________________________________________________________
+void AliVectorSparse::Add(Double_t *valc,Int_t *indc,Int_t n)
+{
+  // add indiced array to row. Indices must be in increasing order
+  int indx;
+  int nadd = 0;
+  //
+  int last = fNElems-1;
+  int mid = 0;
+  for (int i=n;i--;) {
+    // if the element with this index is already defined, just add the value
+    int first = 0;
+    Bool_t toAdd = kTRUE;
+    indx = indc[i];
+    while (first<=last) {
+      mid = (first+last)>>1;
+      if (indx>fIndex[mid]) first = mid+1;
+      else if (indx<fIndex[mid]) last = mid-1;
+      else {
+       fElems[mid] += valc[i];
+       indc[i] = -1;
+       toAdd = kFALSE;
+       last = mid-1;   // profit from the indices being ordered
+       break;
+      }
+    }
+    if (toAdd) nadd++;
+  }
+  //
+  if (nadd<1) return; // nothing to do anymore
+  // 
+  // need to expand the row
+  UShort_t *arrI = new UShort_t[fNElems+nadd];
+  Double_t *arrE = new Double_t[fNElems+nadd];
+  // copy old elems embedding the new ones
+  int inew=0,iold=0;
+  for (int i=0;i<n;i++) {  
+    if ( (indx=indc[i])<0) continue;
+    while (iold<fNElems && fIndex[iold]<indx) {
+      arrI[inew]   = fIndex[iold];
+      arrE[inew++] = fElems[iold++];
+    }
+    arrI[inew] = indx;
+    arrE[inew++] = valc[i];
+  }
+  // copy the rest
+  while (iold<fNElems) {
+    arrI[inew]   = fIndex[iold];
+    arrE[inew++] = fElems[iold++];
+  }
+  //
+  delete[] fIndex;
+  delete[] fElems;
+  fIndex = arrI;
+  fElems = arrE;
+  //
+  fNElems += nadd;
+  //
+}
+
diff --git a/STEER/AliVectorSparse.h b/STEER/AliVectorSparse.h
new file mode 100644 (file)
index 0000000..edc6a74
--- /dev/null
@@ -0,0 +1,68 @@
+#ifndef ALIVECTORSPARSE_H
+#define ALIVECTORSPARSE_H
+
+#include <TObject.h>
+#include <TMath.h>
+
+/**********************************************************************************************/
+/* Sparse vector class, used as row of the AliMatrixSparse class                              */
+/*                                                                                            */ 
+/* Author: ruben.shahoyan@cern.ch                                                             */
+/*                                                                                            */ 
+/**********************************************************************************************/
+
+///////////////////////////////////////////////////////////////////////////////////////
+class AliVectorSparse : public TObject {
+ public:
+  AliVectorSparse();
+  AliVectorSparse(const AliVectorSparse& src);
+  virtual ~AliVectorSparse() {Clear();}
+  virtual void Print(Option_t* option="")                 const;
+  //
+  Int_t     GetNElems()                                   const {return fNElems;}
+  UShort_t *GetIndices()                                  const {return fIndex;}
+  Double_t *GetElems()                                    const {return fElems;}
+  UShort_t& GetIndex(Int_t i)                                   {return fIndex[i];}
+  Double_t& GetElem(Int_t i)                              const {return fElems[i];}
+  void      Clear(Option_t* option="");
+  void      Reset()                                             {memset(fElems,0,fNElems*sizeof(Double_t));}
+  void      ReSize(Int_t sz,Bool_t copy=kFALSE);
+  void      SortIndices(Bool_t valuesToo=kFALSE);
+  void      Add(Double_t *valc,Int_t *indc,Int_t n);
+  //
+  AliVectorSparse& operator=(const AliVectorSparse& src);
+  //
+  virtual Double_t         operator()(Int_t ind)         const;
+  virtual Double_t&        operator()(Int_t ind);
+  virtual void             SetToZero(Int_t ind);
+  Double_t                 FindIndex(Int_t ind)          const;
+  Double_t&                FindIndexAdd(Int_t ind);
+  //
+  Int_t     GetLastIndex()                               const {return fIndex[fNElems-1];}
+  Double_t  GetLastElem()                                const {return fElems[fNElems-1];}
+  Double_t &GetLastElem()                                      {return fElems[fNElems-1];}
+  //
+ protected:
+  Int_t            fNElems;   // 
+  UShort_t*        fIndex;    // Index of stored elems
+  Double_t*        fElems;    // pointer on elements
+  //
+  ClassDef(AliVectorSparse,0)
+};
+
+
+//___________________________________________________
+inline Double_t AliVectorSparse::operator()(Int_t ind) const
+{
+  return FindIndex(ind);
+}
+
+//___________________________________________________
+inline Double_t& AliVectorSparse::operator()(Int_t ind)
+{
+  return FindIndexAdd(ind);
+}
+
+//////////////////////////////////////////////////////////////////////////////////////
+
+#endif
index 1c3f288..3c7ec34 100644 (file)
 #pragma link C++ class AliMillePedeRecord+;
 #pragma link C++ class AliMinResSolve+;
 #pragma link C++ class AliMatrixSparse+;
+#pragma link C++ class AliVectorSparse+;
 #pragma link C++ class AliMatrixSq+;
 #pragma link C++ class AliSymMatrix+;
-
+#pragma link C++ class AliRectMatrix+;
 
 #endif
index 9be9129..13f0dfc 100644 (file)
@@ -68,7 +68,9 @@ AliGRPObject.cxx \
 AliQAv1.cxx \
 AliRunInfo.cxx AliEventInfo.cxx \
 AliRecoParam.cxx AliDetectorRecoParam.cxx \
-AliMillePede2.cxx AliMatrixSq.cxx AliMatrixSparse.cxx AliSymMatrix.cxx AliMinResSolve.cxx
+AliMillePedeRecord.cxx AliMillePede2.cxx AliMatrixSq.cxx \
+AliVectorSparse.cxx AliMatrixSparse.cxx \
+AliSymMatrix.cxx AliRectMatrix.cxx AliMinResSolve.cxx
 
 HDRS:= $(SRCS:.cxx=.h)