STEER/STEER/AliMillePede2.h

   1 #ifndef ALIMILLEPEDE2_H
   2 #define ALIMILLEPEDE2_H
   3
   4 /**********************************************************************************************/
   5 /* General class for alignment with large number of degrees of freedom                        */
   6 /* Based on the original milliped2 by Volker Blobel                                           */
   7 /* http://www.desy.de/~blobel/mptalks.html                                                    */
   8 /*                                                                                            */
   9 /* Author: ruben.shahoyan@cern.ch                                                             */
  10 /*                                                                                            */
  11 /**********************************************************************************************/
  12
  13 #include <TObject.h>
  14 #include <TString.h>
  15 #include <TTree.h>
  16 #include "AliMinResSolve.h"
  17 #include "AliMillePedeRecord.h"
  18 class TFile;
  19 class AliMatrixSq;
  20 class AliSymMatrix;
  21 class AliRectMatrix;
  22 class AliMatrixSparse;
  23 class AliLog;
  24 class TStopwatch;
  25 class TArrayL;
  26 class TArrayF;
  27
  28
  29 class AliMillePede2: public TObject
  30 {
  31  public:
  32   //
  33   enum {kFailed,kInvert,kNoInversion};    // used global matrix solution methods
  34   enum {kFixParID=-1};                    // dummy id for fixed param
  35   //
  36   AliMillePede2();
  37   AliMillePede2(const AliMillePede2& src);
  38   virtual ~AliMillePede2();
  39   AliMillePede2& operator=(const AliMillePede2& )             {printf("Dummy\n"); return *this;}
  40   //
  41   Int_t                InitMille(int nGlo, int nLoc, Int_t lNStdDev=-1,double lResCut=-1., double lResCutInit=-1., const Int_t* regroup=0);
  42   //
  43   Int_t                GetNGloPar()                     const {return fNGloPar;}
  44   Int_t                GetNGloParIni()                  const {return fNGloParIni;}
  45   const Int_t*         GetRegrouping()                  const {return fkReGroup;}
  46   Int_t                GetNLocPar()                     const {return fNLocPar;}
  47   Long_t               GetNLocalEquations()             const {return fNLocEquations;}
  48   Int_t                GetCurrentIteration()            const {return fIter;}
  49   Int_t                GetNMaxIterations()              const {return fMaxIter;}
  50   Int_t                GetNStdDev()                     const {return fNStdDev;}
  51   Int_t                GetNGlobalConstraints()          const {return fNGloConstraints;}
  52   Int_t                GetNLagrangeConstraints()        const {return fNLagrangeConstraints;}
  53   Long_t               GetNLocalFits()                  const {return fNLocFits;}
  54   Long_t               GetNLocalFitsRejected()          const {return fNLocFitsRejected;}
  55   Int_t                GetNGlobalsFixed()               const {return fNGloFix;}
  56   Int_t                GetGlobalSolveStatus()           const {return fGloSolveStatus;}
  57   Float_t              GetChi2CutFactor()               const {return fChi2CutFactor;}
  58   Float_t              GetChi2CutRef()                  const {return fChi2CutRef;}
  59   Float_t              GetResCurInit()                  const {return fResCutInit;}
  60   Float_t              GetResCut()                      const {return fResCut;}
  61   Int_t                GetMinPntValid()                 const {return fMinPntValid;}
  62   Int_t                GetRGId(Int_t i)                 const {return fkReGroup ? (fkReGroup[i]<0 ? -1:fkReGroup[i]) : i;}
  63   Int_t                GetProcessedPoints(Int_t i)      const {int ir=GetRGId(i); return ir<=0 ? 0:fProcPnt[ir];}
  64   Int_t*               GetProcessedPoints()             const {return fProcPnt;}
  65   Int_t                GetParamGrID(Int_t i)            const {int ir=GetRGId(i); return ir<=0 ? 0:fParamGrID[ir];}
  66   //
  67   AliMatrixSq*         GetGlobalMatrix()                const {return fMatCGlo;}
  68   AliSymMatrix*        GetLocalMatrix()                 const {return fMatCLoc;}
  69   Double_t*            GetGlobals()                     const {return fVecBGlo;}
  70   Double_t*            GetDeltaPars()                   const {return fDeltaPar;}
  71   Double_t*            GetInitPars()                    const {return fInitPar;}
  72   Double_t*            GetSigmaPars()                   const {return fSigmaPar;}
  73   Bool_t*              GetIsLinear()                    const {return fIsLinear;}
  74   Double_t             GetFinalParam(int i)             const {int ir=GetRGId(i); return ir<0 ? 0:fDeltaPar[ir]+fInitPar[ir];}
  75   Double_t             GetFinalError(int i)             const {return GetParError(i);}
  76   Double_t             GetPull(int i)                   const;
  77   //
  78   Double_t             GetGlobal(Int_t i)               const {int ir=GetRGId(i); return ir<0 ? 0:fVecBGlo[ir];}
  79   Double_t             GetInitPar(Int_t i)              const {int ir=GetRGId(i); return ir<0 ? 0:fInitPar[ir];}
  80   Double_t             GetSigmaPar(Int_t i)             const {int ir=GetRGId(i); return ir<0 ? 0:fSigmaPar[ir];}
  81   Bool_t               GetIsLinear(Int_t i)             const {int ir=GetRGId(i); return ir<0 ? 0:fIsLinear[ir];}
  82   static Bool_t        IsGlobalMatSparse()                    {return fgIsMatGloSparse;}
  83   static Bool_t        IsWeightSigma()                        {return fgWeightSigma;}
  84   void                 SetWghScale(Double_t wOdd=1,Double_t wEven=1)    {fWghScl[0] = wOdd; fWghScl[1] = wEven;}
  85   void                 SetUseRecordWeight(Bool_t v=kTRUE)     {fUseRecordWeight=v;}
  86   Bool_t               GetUseRecordWeight()             const {return fUseRecordWeight;}
  87   void                 SetMinRecordLength(Int_t v=1)          {fMinRecordLength = v;}
  88   Int_t                GetMinRecordLength()             const {return fMinRecordLength;}
  89   //
  90   void                 SetParamGrID(Int_t grID,Int_t i)       {int ir=GetRGId(i); if(ir<0) return; fParamGrID[ir] = grID; if(fNGroupsSet<grID)fNGroupsSet=grID;}
  91   void                 SetNGloPar(Int_t n)                    {fNGloPar = n;}
  92   void                 SetNLocPar(Int_t n)                    {fNLocPar = n;}
  93   void                 SetNMaxIterations(Int_t n=10)          {fMaxIter = n;}
  94   void                 SetNStdDev(Int_t n)                    {fNStdDev = n;}
  95   void                 SetChi2CutFactor(Float_t v)            {fChi2CutFactor = v;}
  96   void                 SetChi2CutRef(Float_t v)               {fChi2CutRef = v;}
  97   void                 SetResCurInit(Float_t v)               {fResCutInit = v;}
  98   void                 SetResCut(Float_t v)                   {fResCut = v;}
  99   void                 SetMinPntValid(Int_t n)                {fMinPntValid = n>0 ? n:1;}
 100   static void          SetGlobalMatSparse(Bool_t v=kTRUE)     {fgIsMatGloSparse = v;}
 101   static void          SetWeightSigma(Bool_t v=kTRUE)         {fgWeightSigma = v;}
 102   //
 103   void                 SetInitPars(const Double_t* par);
 104   void                 SetSigmaPars(const Double_t* par);
 105   void                 SetInitPar(Int_t i,Double_t par);
 106   void                 SetSigmaPar(Int_t i,Double_t par);
 107   //
 108   Int_t                GlobalFit(Double_t *par=0, Double_t *error=0, Double_t *pull=0);
 109   Int_t                GlobalFitIteration();
 110   Int_t                SolveGlobalMatEq();
 111   static void          SetInvChol(Bool_t v=kTRUE)             {fgInvChol = v;}
 112   static void          SetMinResPrecondType(Int_t tp=0)       {fgMinResCondType = tp;}
 113   static void          SetMinResTol(Double_t val=1e-12)       {fgMinResTol = val;}
 114   static void          SetMinResMaxIter(Int_t val=2000)       {fgMinResMaxIter = val;}
 115   static void          SetIterSolverType(Int_t val=AliMinResSolve::kSolMinRes) {fgIterSol = val;}
 116   static void          SetNKrylovV(Int_t val=60)              {fgNKrylovV = val;}
 117   //
 118   static Bool_t        GetInvChol()                           {return fgInvChol;}
 119   static Int_t         GetMinResPrecondType()                 {return fgMinResCondType;}
 120   static Double_t      GetMinResTol()                         {return fgMinResTol;}
 121   static Int_t         GetMinResMaxIter()                     {return fgMinResMaxIter;}
 122   static Int_t         GetIterSolverType()                    {return fgIterSol;}
 123   static Int_t         GetNKrylovV()                          {return fgNKrylovV;}
 124   //
 125   Double_t             GetParError(int iPar)           const;
 126   Int_t                PrintGlobalParameters()         const;
 127   void                 SetRejRunList(const UInt_t *runs, Int_t nruns);
 128   void                 SetAccRunList(const UInt_t *runs, Int_t nruns, const Float_t* wghList=0);
 129   Bool_t               IsRecordAcceptable();
 130   //
 131   //
 132   Int_t                SetIterations(double lChi2CutFac);
 133
 134   //
 135   // constraints
 136   void                 SetGlobalConstraint(const double *dergb, double val, double sigma=0);
 137   void                 SetGlobalConstraint(const int *indgb, const double *dergb, int ngb, double val, double sigma=0);
 138   //
 139   // processing of the local measurement
 140   void                 SetRecordRun(Int_t run);
 141   void                 SetRecordWeight(double wgh);
 142   void                 SetLocalEquation(double *dergb, double *derlc, double lMeas, double lSigma);
 143   void                 SetLocalEquation(int *indgb, double *dergb, int ngb, int *indlc,
 144                                         double *derlc,int nlc,double lMeas,double lSigma);
 145   //
 146   // manipilation with processed data and costraints records and its buffer
 147   void                 SetDataRecFName(const char* flname)   {fDataRecFName = flname;}
 148   const Char_t*        GetDataRecFName()               const {return fDataRecFName.Data();}
 149   void                 SetConsRecFName(const char* flname)   {fConstrRecFName = flname;}
 150   const Char_t*        GetConsRecFName()               const {return fConstrRecFName.Data();}
 151   //
 152   void   SetRecDataTreeName(const char* name=0)     {fRecDataTreeName = name;   if (fRecDataTreeName.IsNull()) fRecDataTreeName = "AliMillePedeRecords_Data";}
 153   void   SetRecConsTreeName(const char* name=0)     {fRecConsTreeName = name;   if (fRecConsTreeName.IsNull()) fRecConsTreeName = "AliMillePedeRecords_Consaints";}
 154   void   SetRecDataBranchName(const char* name=0)   {fRecDataBranchName = name; if (fRecDataBranchName.IsNull()) fRecDataBranchName = "Record_Data";}
 155   void   SetRecConsBranchName(const char* name=0)   {fRecConsBranchName = name; if (fRecConsBranchName.IsNull()) fRecConsBranchName = "Record_Consaints";}
 156   const char* GetRecDataTreeName()     const {return fRecDataTreeName.Data();}
 157   const char* GetRecConsTreeName()     const {return fRecConsTreeName.Data();}
 158   const char* GetRecDataBranchName()   const {return fRecDataBranchName.Data();}
 159   const char* GetRecConsBranchName()   const {return fRecConsBranchName.Data();}
 160   //
 161   Bool_t               InitDataRecStorage(Bool_t read=kFALSE);
 162   Bool_t               InitConsRecStorage(Bool_t read=kFALSE);
 163   Bool_t               ImposeDataRecFile(const char* fname);
 164   Bool_t               ImposeConsRecFile(const char* fname);
 165   void                 CloseDataRecStorage();
 166   void                 CloseConsRecStorage();
 167   void                 ReadRecordData(Long_t recID);
 168   void                 ReadRecordConstraint(Long_t recID);
 169   Bool_t               ReadNextRecordData();
 170   Bool_t               ReadNextRecordConstraint();
 171   void                 SaveRecordData();
 172   void                 SaveRecordConstraint();
 173   AliMillePedeRecord*  GetRecord()                      const {return fRecord;}
 174   Long_t               GetSelFirst()                    const {return fSelFirst;}
 175   Long_t               GetSelLast()                     const {return fSelLast;}
 176   void                 SetSelFirst(Long_t v)                  {fSelFirst = v;}
 177   void                 SetSelLast(Long_t v)                   {fSelLast = v;}
 178   //
 179   Float_t              Chi2DoFLim(int nSig, int nDoF)   const;
 180   //
 181   // aliases for compatibility with millipede1
 182   void                 SetParSigma(Int_t i,Double_t par)      {SetSigmaPar(i,par);}
 183   void                 SetGlobalParameters(Double_t *par)     {SetInitPars(par);}
 184   void                 SetNonLinear(int index, Bool_t v=kTRUE) {int id = GetRGId(index); if (id<0) return; fIsLinear[id] = !v;}
 185   //
 186  protected:
 187   //
 188   Int_t                LocalFit(double *localParams=0);
 189   Bool_t               IsZero(Double_t v,Double_t eps=1e-16)   const {return TMath::Abs(v)<eps;}
 190   //
 191  protected:
 192   //
 193   Int_t                 fNLocPar;                        // number of local parameters
 194   Int_t                 fNGloPar;                        // number of global parameters
 195   Int_t                 fNGloParIni;                     // number of global parameters before grouping
 196   Int_t                 fNGloSize;                       // final size of the global matrix (NGloPar+NConstraints)
 197   //
 198   Long_t                fNLocEquations;                  // Number of local equations
 199   Int_t                 fIter;                           // Current iteration
 200   Int_t                 fMaxIter;                        // Maximum number of iterations
 201   Int_t                 fNStdDev;                        // Number of standard deviations for chi2 cut
 202   Int_t                 fNGloConstraints;                // Number of constraint equations
 203   Int_t                 fNLagrangeConstraints;           // Number of constraint equations requiring Lagrange multiplier
 204   Long_t                fNLocFits;                       // Number of local fits
 205   Long_t                fNLocFitsRejected;               // Number of local fits rejected
 206   Int_t                 fNGloFix;                        // Number of globals fixed by user
 207   Int_t                 fGloSolveStatus;                 // Status of global solver at current step
 208   //
 209   Float_t               fChi2CutFactor;                  // Cut factor for chi2 cut to accept local fit
 210   Float_t               fChi2CutRef;                     // Reference cut for chi2 cut to accept local fit
 211   Float_t               fResCutInit;                     // Cut in residual for first iterartion
 212   Float_t               fResCut;                         // Cut in residual for other iterartiona
 213   Int_t                 fMinPntValid;                    // min number of points for global to vary
 214   //
 215   Int_t                 fNGroupsSet;                     // number of groups set
 216   Int_t                *fParamGrID;                      //[fNGloPar] group id for the every parameter
 217   Int_t                *fProcPnt;                        //[fNGloPar] N of processed points per global variable
 218   Double_t             *fVecBLoc;                        //[fNLocPar] Vector B local (parameters)
 219   Double_t             *fDiagCGlo;                       //[fNGloPar] Initial diagonal elements of C global matrix
 220   Double_t             *fVecBGlo;                        //! Vector B global (parameters)
 221   //
 222   Double_t             *fInitPar;                        //[fNGloPar] Initial global parameters
 223   Double_t             *fDeltaPar;                       //[fNGloPar] Variation of global parameters
 224   Double_t             *fSigmaPar;                       //[fNGloPar] Sigma of allowed variation of global parameter
 225   //
 226   Bool_t               *fIsLinear;                       //[fNGloPar] Flag for linear parameters
 227   Bool_t               *fConstrUsed;                     //! Flag for used constraints
 228   //
 229   Int_t                *fGlo2CGlo;                       //[fNGloPar] global ID to compressed ID buffer
 230   Int_t                *fCGlo2Glo;                       //[fNGloPar] compressed ID to global ID buffer
 231   //
 232   // Matrices
 233   AliSymMatrix         *fMatCLoc;                        // Matrix C local
 234   AliMatrixSq          *fMatCGlo;                        // Matrix C global
 235   AliRectMatrix        *fMatCGloLoc;                     // Rectangular matrix C g*l
 236   Int_t                *fFillIndex;                      //[fNGloPar] auxilary index array for fast matrix fill
 237   Double_t             *fFillValue;                      //[fNGloPar] auxilary value array for fast matrix fill
 238   //
 239   // processed data record bufferization
 240   TString               fRecDataTreeName;                // Name of data records tree
 241   TString               fRecConsTreeName;                // Name of constraints records tree
 242   TString               fRecDataBranchName;              // Name of data records branch name
 243   TString               fRecConsBranchName;              // Name of constraints records branch name
 244
 245   TString               fDataRecFName;                   // Name of File for data records
 246   AliMillePedeRecord   *fRecord;                         // Buffer of measurements records
 247   TFile                *fDataRecFile;                    // File of processed measurements records
 248   TTree                *fTreeData;                       // Tree of processed measurements records
 249   Int_t                 fRecFileStatus;                  // state of the record file (0-no, 1-read, 2-rw)
 250   //
 251   TString               fConstrRecFName;                 // Name of File for constraints records
 252   TTree                *fTreeConstr;                     //! Tree of constraint records
 253   TFile                *fConsRecFile;                    //! File of processed constraints records
 254   Long_t                fCurrRecDataID;                  // ID of the current data record
 255   Long_t                fCurrRecConstrID;                // ID of the current constraint record
 256   Bool_t                fLocFitAdd;                      // Add contribution of carrent track (and not eliminate it)
 257   Bool_t                fUseRecordWeight;                // force or ignore the record weight
 258   Int_t                 fMinRecordLength;                // ignore shorter records
 259   Int_t                 fSelFirst;                       // event selection start
 260   Int_t                 fSelLast;                        // event selection end
 261   TArrayL*              fRejRunList;                     // list of runs to reject (if any)
 262   TArrayL*              fAccRunList;                     // list of runs to select (if any)
 263   TArrayF*              fAccRunListWgh;                  // optional weights for data of accepted runs (if any)
 264   Double_t              fRunWgh;                         // run weight
 265   Double_t              fWghScl[2];                      // optional rescaling for odd/even residual weights (see its usage in LocalFit)
 266   const Int_t*          fkReGroup;                       // optional regrouping of parameters wrt ID's from the records
 267   //
 268   static Bool_t         fgInvChol;                       // Invert global matrix in Cholesky solver
 269   static Bool_t         fgWeightSigma;                   // weight parameter constraint by statistics
 270   static Bool_t         fgIsMatGloSparse;                // Type of the global matrix (sparse ...)
 271   static Int_t          fgMinResCondType;                // Type of the preconditioner for MinRes method
 272   static Double_t       fgMinResTol;                     // Tolerance for MinRes solution
 273   static Int_t          fgMinResMaxIter;                 // Max number of iterations for the MinRes method
 274   static Int_t          fgIterSol;                       // type of iterative solution: MinRes or FGMRES
 275   static Int_t          fgNKrylovV;                      // size of Krylov vectors buffer in FGMRES
 276   //
 277   ClassDef(AliMillePede2,1)
 278 };
 279
 280 //_____________________________________________________________________________________________
 281 inline void AliMillePede2::ReadRecordData(Long_t recID)       {fTreeData->GetEntry(recID); fCurrRecDataID=recID;}
 282
 283 //_____________________________________________________________________________________________
 284 inline void AliMillePede2::ReadRecordConstraint(Long_t recID) {fTreeConstr->GetEntry(recID); fCurrRecConstrID=recID;}
 285
 286 //_____________________________________________________________________________________________
 287 inline void AliMillePede2::SaveRecordData()                   {fTreeData->Fill(); fRecord->Reset(); fCurrRecDataID++;}
 288
 289 //_____________________________________________________________________________________________
 290 inline void AliMillePede2::SaveRecordConstraint()             {fTreeConstr->Fill(); fRecord->Reset();fCurrRecConstrID++;}
 291
 292 #endif