Corrected protection.
[u/mrichter/AliRoot.git] / 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
27
28 class AliMillePede2: public TObject
29 {
30  public:
31   //
32   enum {gkFailed,gkInvert,gkNoInversion};    // used global matrix solution methods
33   //
34   AliMillePede2();
35   AliMillePede2(const AliMillePede2& src);
36   virtual ~AliMillePede2();
37   AliMillePede2& operator=(const AliMillePede2& )             {printf("Dummy\n"); return *this;}
38   //
39   Int_t                InitMille(int nGlo, int nLoc, Int_t lNStdDev=-1,double lResCut=-1., double lResCutInit=-1.);
40   //
41   Int_t                GetNGloPar()                     const {return fNGloPar;}
42   Int_t                GetNLocPar()                     const {return fNLocPar;}
43   Long_t               GetNLocalEquations()             const {return fNLocEquations;}
44   Int_t                GetCurrentIteration()            const {return fIter;}
45   Int_t                GetNMaxIterations()              const {return fMaxIter;}
46   Int_t                GetNStdDev()                     const {return fNStdDev;} 
47   Int_t                GetNGlobalConstraints()          const {return fNGloConstraints;}
48   Int_t                GetNLagrangeConstraints()        const {return fNLagrangeConstraints;}
49   Long_t               GetNLocalFits()                  const {return fNLocFits;}
50   Long_t               GetNLocalFitsRejected()          const {return fNLocFitsRejected;}
51   Int_t                GetNGlobalsFixed()               const {return fNGloFix;}
52   Int_t                GetGlobalSolveStatus()           const {return fGloSolveStatus;}
53   Float_t              GetChi2CutFactor()               const {return fChi2CutFactor;}
54   Float_t              GetChi2CutRef()                  const {return fChi2CutRef;}
55   Float_t              GetResCurInit()                  const {return fResCutInit;}
56   Float_t              GetResCut()                      const {return fResCut;}
57   Int_t                GetMinPntValid()                 const {return fMinPntValid;}
58   Int_t                GetProcessedPoints(Int_t i)      const {return fProcPnt[i];}
59   Int_t*               GetProcessedPoints()             const {return fProcPnt;}
60   Int_t                GetParamGrID(Int_t i)            const {return fParamGrID[i];}
61   //
62   AliMatrixSq*         GetGlobalMatrix()                const {return fMatCGlo;}
63   AliSymMatrix*        GetLocalMatrix()                 const {return fMatCLoc;}
64   Double_t*            GetGlobals()                     const {return fVecBGlo;}
65   Double_t*            GetDeltaPars()                   const {return fDeltaPar;}
66   Double_t*            GetInitPars()                    const {return fInitPar;}
67   Double_t*            GetSigmaPars()                   const {return fSigmaPar;}
68   Bool_t*              GetIsLinear()                    const {return fIsLinear;}
69   Double_t             GetFinalParam(int i)             const {return fDeltaPar[i]+fInitPar[i];}
70   Double_t             GetFinalError(int i)             const {return GetParError(i);}
71   //
72   Double_t             GetGlobal(Int_t i)               const {return fVecBGlo[i];}
73   Double_t             GetInitPar(Int_t i)              const {return fInitPar[i];}
74   Double_t             GetSigmaPar(Int_t i)             const {return fSigmaPar[i];}
75   Bool_t               GetIsLinear(Int_t i)             const {return fIsLinear[i];}
76   static Bool_t        IsGlobalMatSparse()                    {return fgIsMatGloSparse;}
77   static Bool_t        IsWeightSigma()                        {return fgWeightSigma;}
78   //
79   void                 SetParamGrID(Int_t grID,Int_t i)       {fParamGrID[i] = grID; if(fNGroupsSet<grID)fNGroupsSet=grID;}
80   void                 SetNGloPar(Int_t n)                    {fNGloPar = n;}
81   void                 SetNLocPar(Int_t n)                    {fNLocPar = n;}
82   void                 SetNMaxIterations(Int_t n=10)          {fMaxIter = n;}
83   void                 SetNStdDev(Int_t n)                    {fNStdDev = n;}
84   void                 SetChi2CutFactor(Float_t v)            {fChi2CutFactor = v;}
85   void                 SetChi2CutRef(Float_t v)               {fChi2CutRef = v;}
86   void                 SetResCurInit(Float_t v)               {fResCutInit = v;}
87   void                 SetResCut(Float_t v)                   {fResCut = v;}
88   void                 SetMinPntValid(Int_t n)                {fMinPntValid = n>0 ? n:1;}
89   static void          SetGlobalMatSparse(Bool_t v=kTRUE)     {fgIsMatGloSparse = v;}
90   static void          SetWeightSigma(Bool_t v=kTRUE)         {fgWeightSigma = v;}
91   //
92   void                 SetInitPars(const Double_t* par)       {memcpy(fInitPar,par,fNGloPar*sizeof(Double_t));}
93   void                 SetSigmaPars(const Double_t* par)      {memcpy(fSigmaPar,par,fNGloPar*sizeof(Double_t));}
94   void                 SetInitPar(Int_t i,Double_t par)       {fInitPar[i] = par;;}
95   void                 SetSigmaPar(Int_t i,Double_t par)      {fSigmaPar[i] = par;}
96   //
97   Int_t                GlobalFit(Double_t *par=0, Double_t *error=0, Double_t *pull=0);
98   Int_t                GlobalFitIteration();
99   Int_t                SolveGlobalMatEq();
100   static void          SetInvChol(Bool_t v=kTRUE)             {fgInvChol = v;}
101   static void          SetMinResPrecondType(Int_t tp=0)       {fgMinResCondType = tp;}
102   static void          SetMinResTol(Double_t val=1e-12)       {fgMinResTol = val;}
103   static void          SetMinResMaxIter(Int_t val=2000)       {fgMinResMaxIter = val;}
104   static void          SetIterSolverType(Int_t val=AliMinResSolve::kSolMinRes) {fgIterSol = val;}
105   static void          SetNKrylovV(Int_t val=60)              {fgNKrylovV = val;}
106   //
107   static Bool_t        GetInvChol()                           {return fgInvChol;}
108   static Int_t         GetMinResPrecondType()                 {return fgMinResCondType;}
109   static Double_t      GetMinResTol()                         {return fgMinResTol;}
110   static Int_t         GetMinResMaxIter()                     {return fgMinResMaxIter;}
111   static Int_t         GetIterSolverType()                    {return fgIterSol;}
112   static Int_t         GetNKrylovV()                          {return fgNKrylovV;}
113   //
114   Double_t             GetParError(int iPar)           const;
115   Int_t                PrintGlobalParameters()         const;
116   void                 SetRejRunList(const UInt_t *runs, Int_t nruns);
117   void                 SetAccRunList(const UInt_t *runs, Int_t nruns);
118   Bool_t               IsRecordAcceptable() const;
119   //
120   //
121   Int_t                SetIterations(double lChi2CutFac);
122
123   //
124   // constraints
125   void                 SetGlobalConstraint(double *dergb, double val, double sigma=0);
126   void                 SetGlobalConstraint(const int *indgb,double *dergb, int ngb, double val, double sigma=0);
127   //
128   // processing of the local measurement
129   void                 SetRecordRun(Int_t run);
130   void                 SetRecordWeight(double wgh);
131   void                 SetLocalEquation(double *dergb, double *derlc, double lMeas, double lSigma);
132   void                 SetLocalEquation(int *indgb, double *dergb, int ngb, int *indlc, 
133                                         double *derlc,int nlc,double lMeas,double lSigma);
134   //
135   // manipilation with processed data and costraints records and its buffer
136   void                 SetDataRecFName(const char* flname)   {fDataRecFName = flname;}
137   const Char_t*        GetDataRecFName()               const {return fDataRecFName.Data();}
138   void                 SetConsRecFName(const char* flname)   {fConstrRecFName = flname;}
139   const Char_t*        GetConsRecFName()               const {return fConstrRecFName.Data();}
140
141   Bool_t               InitDataRecStorage(Bool_t read=kFALSE);
142   Bool_t               InitConsRecStorage(Bool_t read=kFALSE);
143   Bool_t               ImposeDataRecFile(const char* fname);
144   Bool_t               ImposeConsRecFile(const char* fname);
145   void                 CloseDataRecStorage();
146   void                 CloseConsRecStorage();
147   void                 ReadRecordData(Long_t recID);
148   void                 ReadRecordConstraint(Long_t recID);
149   Bool_t               ReadNextRecordData();
150   Bool_t               ReadNextRecordConstraint();
151   void                 SaveRecordData();
152   void                 SaveRecordConstraint();
153   AliMillePedeRecord*  GetRecord()                      const {return fRecord;}
154   Long_t               GetSelFirst()                    const {return fSelFirst;}
155   Long_t               GetSelLast()                     const {return fSelLast;}
156   void                 SetSelFirst(Long_t v)                  {fSelFirst = v;}
157   void                 SetSelLast(Long_t v)                   {fSelLast = v;}
158   //
159   Float_t              Chi2DoFLim(int nSig, int nDoF)   const;
160   //
161   // aliases for compatibility with millipede1
162   void                 SetParSigma(Int_t i,Double_t par)      {SetSigmaPar(i,par);}
163   void                 SetGlobalParameters(Double_t *par)     {SetInitPars(par);}
164   //
165  protected:
166   //
167   Int_t                LocalFit(double *localParams=0);
168   Bool_t               IsZero(Double_t v,Double_t eps=1e-16)   const {return TMath::Abs(v)<eps;}                  
169   //
170  protected:
171   //
172   Int_t                 fNLocPar;                        // number of local parameters
173   Int_t                 fNGloPar;                        // number of global parameters
174   Int_t                 fNGloSize;                       // final size of the global matrix (NGloPar+NConstraints)
175   //
176   Long_t                fNLocEquations;                  // Number of local equations 
177   Int_t                 fIter;                           // Current iteration
178   Int_t                 fMaxIter;                        // Maximum number of iterations
179   Int_t                 fNStdDev;                        // Number of standard deviations for chi2 cut 
180   Int_t                 fNGloConstraints;                // Number of constraint equations
181   Int_t                 fNLagrangeConstraints;           // Number of constraint equations requiring Lagrange multiplier
182   Long_t                fNLocFits;                       // Number of local fits
183   Long_t                fNLocFitsRejected;               // Number of local fits rejected
184   Int_t                 fNGloFix;                        // Number of globals fixed by user
185   Int_t                 fGloSolveStatus;                 // Status of global solver at current step
186   //
187   Float_t               fChi2CutFactor;                  // Cut factor for chi2 cut to accept local fit 
188   Float_t               fChi2CutRef;                     // Reference cut for chi2 cut to accept local fit 
189   Float_t               fResCutInit;                     // Cut in residual for first iterartion
190   Float_t               fResCut;                         // Cut in residual for other iterartiona
191   Int_t                 fMinPntValid;                    // min number of points for global to vary
192   //
193   Int_t                 fNGroupsSet;                     // number of groups set
194   Int_t                *fParamGrID;                      // group id for the every parameter
195   Int_t                *fProcPnt;                        // N of processed points per global variable
196   Double_t             *fVecBLoc;                        // Vector B local (parameters) 
197   Double_t             *fDiagCGlo;                       // Initial diagonal elements of C global matrix
198   Double_t             *fVecBGlo;                        // Vector B global (parameters)
199   //
200   Double_t             *fInitPar;                        // Initial global parameters
201   Double_t             *fDeltaPar;                       // Variation of global parameters
202   Double_t             *fSigmaPar;                       // Sigma of allowed variation of global parameter
203   //
204   Bool_t               *fIsLinear;                       // Flag for linear parameters
205   Bool_t               *fConstrUsed;                     // Flag for used constraints
206   //
207   Int_t                *fGlo2CGlo;                       // global ID to compressed ID buffer
208   Int_t                *fCGlo2Glo;                       // compressed ID to global ID buffer
209   //
210   // Matrices
211   AliSymMatrix         *fMatCLoc;                        // Matrix C local
212   AliMatrixSq          *fMatCGlo;                        // Matrix C global
213   AliRectMatrix        *fMatCGloLoc;                     // Rectangular matrix C g*l 
214   Int_t                *fFillIndex;                      // auxilary index array for fast matrix fill
215   Double_t             *fFillValue;                      // auxilary value array for fast matrix fill
216   //
217   // processed data record bufferization   
218   TString               fDataRecFName;                   // Name of File for data records               
219   AliMillePedeRecord   *fRecord;                         // Buffer of measurements records
220   TFile                *fDataRecFile;                    // File of processed measurements records
221   TTree                *fTreeData;                       // Tree of processed measurements records
222   Int_t                 fRecFileStatus;                  // state of the record file (0-no, 1-read, 2-rw)
223   //
224   TString               fConstrRecFName;                 // Name of File for constraints records               
225   TTree                *fTreeConstr;                     // Tree of constraint records
226   TFile                *fConsRecFile;                    // File of processed constraints records
227   Long_t                fCurrRecDataID;                  // ID of the current data record
228   Long_t                fCurrRecConstrID;                // ID of the current constraint record
229   Bool_t                fLocFitAdd;                      // Add contribution of carrent track (and not eliminate it)
230   Int_t                 fSelFirst;                       // event selection start
231   Int_t                 fSelLast;                        // event selection end
232   TArrayL*              fRejRunList;                     // list of runs to reject (if any)
233   TArrayL*              fAccRunList;                     // list of runs to select (if any)
234   //
235   static Bool_t         fgInvChol;                       // Invert global matrix in Cholesky solver
236   static Bool_t         fgWeightSigma;                   // weight parameter constraint by statistics
237   static Bool_t         fgIsMatGloSparse;                // Type of the global matrix (sparse ...)
238   static Int_t          fgMinResCondType;                // Type of the preconditioner for MinRes method 
239   static Double_t       fgMinResTol;                     // Tolerance for MinRes solution
240   static Int_t          fgMinResMaxIter;                 // Max number of iterations for the MinRes method
241   static Int_t          fgIterSol;                       // type of iterative solution: MinRes or FGMRES
242   static Int_t          fgNKrylovV;                      // size of Krylov vectors buffer in FGMRES
243   //
244   ClassDef(AliMillePede2,0)
245 };
246
247 //_____________________________________________________________________________________________
248 inline void AliMillePede2::ReadRecordData(Long_t recID)       {fTreeData->GetEntry(recID); fCurrRecDataID=recID;}
249
250 //_____________________________________________________________________________________________
251 inline void AliMillePede2::ReadRecordConstraint(Long_t recID) {fTreeConstr->GetEntry(recID); fCurrRecConstrID=recID;}
252
253 //_____________________________________________________________________________________________
254 inline void AliMillePede2::SaveRecordData()                   {fTreeData->Fill(); fRecord->Reset(); fCurrRecDataID++;}
255
256 //_____________________________________________________________________________________________
257 inline void AliMillePede2::SaveRecordConstraint()             {fTreeConstr->Fill(); fRecord->Reset();fCurrRecConstrID++;}
258
259 #endif