STEER/AliMillePede2.h

   1 #ifndef ALIMILLEPEDE2_H
   2 #define ALIMILLEPEDE2_H
   3
   4 #include <TObject.h>
   5 #include <TTree.h>
   6 #include <TFile.h>
   7 #include <TMatrixD.h>
   8 #include "AliSymMatrix.h"
   9 #include "AliMatrixSparse.h"
  10 #include "AliMinResSolve.h"
  11 class AliLog;
  12 class TStopwatch;
  13
  14 /**********************************************************************************************/
  15 /* AliMillePedeRecords: class to store the data of single track processing                    */
  16 /* Format: for each measured point the data is stored consequtively                           */
  17 /* INDEX                                                      VALUE                           */
  18 /* -1                                                         residual                        */
  19 /* Local_param_id                                             dResidual/dLocal_param          */
  20 /* ...                                                        ...                             */
  21 /* -2                                                         weight of the measurement       */
  22 /* Global_param_od                                            dResidual/dGlobal_param         */
  23 /* ...                                                        ...                             */
  24 /*                                                                                            */
  25 /* The records for all processed tracks are stored in the temporary tree in orgder to be      */
  26 /* reused for multiple iterations of MillePede                                                */
  27 /*                                                                                            */
  28 /**********************************************************************************************/
  29
  30 class AliMillePedeRecord : public TObject
  31 {
  32  public:
  33   AliMillePedeRecord();
  34   AliMillePedeRecord(const AliMillePedeRecord& src);
  35   AliMillePedeRecord& operator=(const AliMillePedeRecord& rhs);
  36   //
  37   virtual    ~AliMillePedeRecord();
  38   void       Reset()                                             {fSize = 0;}
  39   void       Print(const Option_t *opt="")                 const;
  40   //
  41   Int_t      GetSize()                                     const {return fSize;}
  42   Int_t     *GetIndex()                                    const {return fIndex;}
  43   Int_t      GetIndex(int i)                               const {return fIndex[i];}
  44   //
  45   void       GetIndexValue(int i,int &ind,double &val)     const {ind=fIndex[i]; val=fValue[i];}
  46   void       AddIndexValue(int ind, double val);
  47   void       AddResidual(double val)                             {AddIndexValue(-1,val);}
  48   void       AddWeight(double val)                               {AddIndexValue(-2,val);}
  49   Bool_t     IsResidual(int i)                             const {return fIndex[i]==-1;}
  50   Bool_t     IsWeight(int i)                               const {return fIndex[i]==-2;}
  51   //
  52   Double_t  *GetValue()                                    const {return fValue;}
  53   Double_t   GetValue(int i)                               const {return fValue[i];}
  54   //
  55  protected:
  56   Int_t      GetBufferSize()                               const {return GetUniqueID();}
  57   void       SetBufferSize(int sz)                               {SetUniqueID(sz);}
  58   void       Expand(int bfsize);
  59   //
  60  protected:
  61   Int_t      fSize;                             // size of the record
  62   Int_t   *  fIndex;                            //[fSize] index of variables
  63   Double_t*  fValue;                            //[fSize] array of values: derivs,residuals
  64   //
  65   ClassDef(AliMillePedeRecord,1)                // Record of track residuals and local/global deriavtives
  66 };
  67
  68 inline void  AliMillePedeRecord::AddIndexValue(int ind, double val)
  69 {
  70   if (fSize>=GetBufferSize()) Expand(2*(fSize+1));
  71   fIndex[fSize]=ind;
  72   fValue[fSize++]=val;
  73 }
  74
  75 //----------------------------------------------------------------------------------------------
  76 //----------------------------------------------------------------------------------------------
  77 //----------------------------------------------------------------------------------------------
  78 //----------------------------------------------------------------------------------------------
  79
  80 class AliMillePede2: public TObject
  81 {
  82  public:
  83   //
  84   enum {gkFailed,gkInvert,gkMinRes};    // used global matrix solution methods
  85   //
  86   AliMillePede2();
  87   AliMillePede2(const AliMillePede2& src);
  88   virtual ~AliMillePede2();
  89   AliMillePede2& operator=(const AliMillePede2& )             {printf("Dummy\n"); return *this;}
  90   //
  91   Int_t                InitMille(int nGlo, int nLoc, Int_t lNStdDev=-1,double lResCut=-1., double lResCutInit=-1.);
  92   //
  93   Int_t                GetNGloPar()                     const {return fNGloPar;}
  94   Int_t                GetNLocPar()                     const {return fNLocPar;}
  95   Long_t               GetNLocalEquations()             const {return fNLocEquations;}
  96   Int_t                GetCurrentIteration()            const {return fIter;}
  97   Int_t                GetNMaxIterations()              const {return fMaxIter;}
  98   Int_t                GetNStdDev()                     const {return fNStdDev;}
  99   Int_t                GetNGlobalConstraints()          const {return fNGloConstraints;}
 100   Long_t               GetNLocalFits()                  const {return fNLocFits;}
 101   Long_t               GetNLocalFitsRejected()          const {return fNLocFitsRejected;}
 102   Int_t                GetNGlobalsFixed()               const {return fNGloFix;}
 103   Int_t                GetGlobalSolveStatus()           const {return fGloSolveStatus;}
 104   Float_t              GetChi2CutFactor()               const {return fChi2CutFactor;}
 105   Float_t              GetChi2CutRef()                  const {return fChi2CutRef;}
 106   Float_t              GetResCurInit()                  const {return fResCutInit;}
 107   Float_t              GetResCut()                      const {return fResCut;}
 108   Int_t                GetMinPntValid()                 const {return fMinPntValid;}
 109   Int_t                GetProcessedPoints(Int_t i)      const {return fProcPnt[i];}
 110   Int_t*               GetProcessedPoints()             const {return fProcPnt;}
 111   //
 112   AliMatrixSq*         GetGlobalMatrix()                const {return fMatCGlo;}
 113   Double_t*            GetGlobals()                     const {return fVecBGlo;}
 114   Double_t*            GetDeltaPars()                   const {return fDeltaPar;}
 115   Double_t*            GetInitPars()                    const {return fInitPar;}
 116   Double_t*            GetSigmaPars()                   const {return fSigmaPar;}
 117   Bool_t*              GetIsLinear()                    const {return fIsLinear;}
 118   Double_t             GetGlobal(Int_t i)               const {return fVecBGlo[i];}
 119   Double_t             GetInitPar(Int_t i)              const {return fInitPar[i];}
 120   Double_t             GetSigmaPar(Int_t i)             const {return fSigmaPar[i];}
 121   Bool_t               GetIsLinear(Int_t i)             const {return fIsLinear[i];}
 122   static Bool_t        IsGlobalMatSparse()                    {return fgIsMatGloSparse;}
 123   //
 124   void                 SetNGloPar(Int_t n)                    {fNGloPar = n;}
 125   void                 SetNLocPar(Int_t n)                    {fNLocPar = n;}
 126   void                 SetNMaxIterations(Int_t n=10)          {fMaxIter = n;}
 127   void                 SetNStdDev(Int_t n)                    {fNStdDev = n;}
 128   void                 SetChi2CutFactor(Float_t v)            {fChi2CutFactor = v;}
 129   void                 SetChi2CutRef(Float_t v)               {fChi2CutRef = v;}
 130   void                 SetResCurInit(Float_t v)               {fResCutInit = v;}
 131   void                 SetResCut(Float_t v)                   {fResCut = v;}
 132   void                 SetMinPntValid(Int_t n)                {fMinPntValid = n;}
 133   static void          SetGlobalMatSparse(Bool_t v=kTRUE)     {fgIsMatGloSparse = v;}
 134   //
 135   void                 SetInitPars(const Double_t* par)       {memcpy(fInitPar,par,fNGloPar*sizeof(Double_t));}
 136   void                 SetSigmaPars(const Double_t* par)      {memcpy(fSigmaPar,par,fNGloPar*sizeof(Double_t));}
 137   void                 SetInitPar(Int_t i,Double_t par)       {fInitPar[i] = par;;}
 138   void                 SetSigmaPar(Int_t i,Double_t par)      {fSigmaPar[i] = par;}
 139   //
 140   Int_t                GlobalFit(Double_t *par=0, Double_t *error=0, Double_t *pull=0);
 141   Int_t                GlobalFitIteration();
 142   Int_t                SolveGlobalMatEq();
 143   static void          SetMinResPrecondType(Int_t tp=0)       {fgMinResCondType = tp;}
 144   static void          SetMinResTol(Double_t val=1e-12)       {fgMinResTol = val;}
 145   static void          SetMinResMaxIter(Int_t val=2000)       {fgMinResMaxIter = val;}
 146   static void          SetIterSolverType(Int_t val=AliMinResSolve::kSolMinRes) {fgIterSol = val;}
 147   static void          SetNKrylovV(Int_t val=60)              {fgNKrylovV = val;}
 148   //
 149   static Int_t         GetMinResPrecondType()                 {return fgMinResCondType;}
 150   static Double_t      GetMinResTol()                         {return fgMinResTol;}
 151   static Int_t         GetMinResMaxIter()                     {return fgMinResMaxIter;}
 152   static Int_t         GetIterSolverType()                    {return fgIterSol;}
 153   static Int_t         GetNKrylovV()                          {return fgNKrylovV;}
 154   //
 155   Double_t             GetParError(int iPar)           const;
 156   Int_t                PrintGlobalParameters()         const;
 157   //
 158   //
 159   Int_t                SetIterations(double lChi2CutFac);
 160
 161   //
 162   // constraints
 163   void                 SetGlobalConstraint(double *dergb, double val);
 164   void                 SetGlobalConstraint(const int *indgb,double *dergb, int ngb, double val);
 165   //
 166   // processing of the local measurement
 167   void                 SetLocalEquation(double *dergb, double *derlc, double lMeas, double lSigma);
 168   void                 SetLocalEquation(int *indgb, double *dergb, int ngb, int *indlc,
 169                                         double *derlc,int nlc,double lMeas,double lSigma);
 170   //
 171   // manipilation with processed data and costraints records and its buffer
 172   void                 SetDataRecFName(const char* flname)   {fDataRecFName = flname;}
 173   const Char_t*        GetDataRecFName()               const {return fDataRecFName.Data();}
 174   void                 SetConsRecFName(const char* flname)   {fConstrRecFName = flname;}
 175   const Char_t*        GetConsRecFName()               const {return fConstrRecFName.Data();}
 176
 177   Bool_t               InitDataRecStorage(Bool_t read=kFALSE);
 178   Bool_t               InitConsRecStorage(Bool_t read=kFALSE);
 179   Bool_t               ImposeDataRecFile(const char* fname);
 180   Bool_t               ImposeConsRecFile(const char* fname);
 181   void                 CloseDataRecStorage();
 182   void                 CloseConsRecStorage();
 183   void                 ReadRecordData(Long_t recID);
 184   void                 ReadRecordConstraint(Long_t recID);
 185   Bool_t               ReadNextRecordData();
 186   Bool_t               ReadNextRecordConstraint();
 187   void                 SaveRecordData();
 188   void                 SaveRecordConstraint();
 189   void                 ResetConstraints();
 190   AliMillePedeRecord*  GetRecord()                      const {return fRecord;}
 191   //
 192   Float_t              Chi2DoFLim(int nSig, int nDoF)   const;
 193   //
 194   // aliases for compatibility with millipede1
 195   void                 SetParSigma(Int_t i,Double_t par)      {SetSigmaPar(i,par);}
 196   void                 SetGlobalParameters(Double_t *par)     {SetInitPars(par);}
 197   //
 198  protected:
 199   //
 200   Int_t                LocalFit(double *localParams=0);
 201   //
 202  protected:
 203   //
 204   Int_t                 fNLocPar;                        // number of local parameters
 205   Int_t                 fNGloPar;                        // number of global parameters
 206   Int_t                 fNGloSize;                       // final size of the global matrix (NGloPar+NConstraints)
 207   //
 208   Long_t                fNLocEquations;                  // Number of local equations
 209   Int_t                 fIter;                           // Current iteration
 210   Int_t                 fMaxIter;                        // Maximum number of iterations
 211   Int_t                 fNStdDev;                        // Number of standard deviations for chi2 cut
 212   Int_t                 fNGloConstraints;                // Number of constraint equations
 213   Long_t                fNLocFits;                       // Number of local fits
 214   Long_t                fNLocFitsRejected;               // Number of local fits rejected
 215   Int_t                 fNGloFix;                        // Number of globals fixed by user
 216   Int_t                 fGloSolveStatus;                 // Status of global solver at current step
 217   //
 218   Float_t               fChi2CutFactor;                  // Cut factor for chi2 cut to accept local fit
 219   Float_t               fChi2CutRef;                     // Reference cut for chi2 cut to accept local fit
 220   Float_t               fResCutInit;                     // Cut in residual for first iterartion
 221   Float_t               fResCut;                         // Cut in residual for other iterartiona
 222   Int_t                 fMinPntValid;                    // min number of points for global to vary
 223   //
 224   Int_t                *fProcPnt;                        // N of processed points per global variable
 225   Double_t             *fVecBLoc;                        // Vector B local (parameters)
 226   Double_t             *fDiagCGlo;                       // Initial diagonal elements of C global matrix
 227   Double_t             *fVecBGlo;                        // Vector B global (parameters)
 228   //
 229   Double_t             *fInitPar;                        // Initial global parameters
 230   Double_t             *fDeltaPar;                       // Variation of global parameters
 231   Double_t             *fSigmaPar;                       // Sigma of allowed variation of global parameter
 232   //
 233   Bool_t               *fIsLinear;                       // Flag for linear parameters
 234   Bool_t               *fConstrUsed;                     // Flag for used constraints
 235   //
 236   Int_t                *fGlo2CGlo;                       // global ID to compressed ID buffer
 237   Int_t                *fCGlo2Glo;                       // compressed ID to global ID buffer
 238   //
 239   // Matrices
 240   AliSymMatrix         *fMatCLoc;                        // Matrix C local
 241   AliMatrixSq          *fMatCGlo;                        // Matrix C global
 242   TMatrixD             *fMatCGloLoc;                     // Rectangular matrix C g*l
 243   //
 244   // processed data record bufferization
 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   //
 250   TString               fConstrRecFName;                 // Name of File for constraints records
 251   TTree                *fTreeConstr;                     // Tree of constraint records
 252   TFile                *fConsRecFile;                    // File of processed constraints records
 253   Long_t                fCurrRecDataID;                  // ID of the current data record
 254   Long_t                fCurrRecConstrID;                // ID of the current constraint record
 255   //
 256   static Bool_t         fgIsMatGloSparse;                // Type of the global matrix (sparse ...)
 257   static Int_t          fgMinResCondType;                // Type of the preconditioner for MinRes method
 258   static Double_t       fgMinResTol;                     // Tolerance for MinRes solution
 259   static Int_t          fgMinResMaxIter;                 // Max number of iterations for the MinRes method
 260   static Int_t          fgIterSol;                       // type of iterative solution: MinRes or FGMRES
 261   static Int_t          fgNKrylovV;                      // size of Krylov vectors buffer in FGMRES
 262   //
 263   ClassDef(AliMillePede2,0)
 264 };
 265
 266 //_____________________________________________________________________________________________
 267 inline void AliMillePede2::ReadRecordData(Long_t recID)       {fTreeData->GetEntry(recID); fCurrRecDataID=recID;}
 268
 269 //_____________________________________________________________________________________________
 270 inline void AliMillePede2::ReadRecordConstraint(Long_t recID) {fTreeConstr->GetEntry(recID); fCurrRecConstrID=recID;}
 271
 272 //_____________________________________________________________________________________________
 273 inline void AliMillePede2::SaveRecordData()                   {fTreeData->Fill(); fRecord->Reset(); fCurrRecDataID++;}
 274
 275 //_____________________________________________________________________________________________
 276 inline void AliMillePede2::SaveRecordConstraint()             {fTreeConstr->Fill(); fRecord->Reset();fCurrRecConstrID++;}
 277
 278 #endif