]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliMillePede2.h
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
New macro to keep track of timing performances of the segmentation methods (Laurent)
[u/mrichter/AliRoot.git] / STEER / AliMillePede2.h
CommitLineData
8a9ab0eb 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"
11class AliLog;
12class 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
30class 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
68inline 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
80class 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//_____________________________________________________________________________________________
267inline void AliMillePede2::ReadRecordData(Long_t recID) {fTreeData->GetEntry(recID); fCurrRecDataID=recID;}
268
269//_____________________________________________________________________________________________
270inline void AliMillePede2::ReadRecordConstraint(Long_t recID) {fTreeConstr->GetEntry(recID); fCurrRecConstrID=recID;}
271
272//_____________________________________________________________________________________________
273inline void AliMillePede2::SaveRecordData() {fTreeData->Fill(); fRecord->Reset(); fCurrRecDataID++;}
274
275//_____________________________________________________________________________________________
276inline void AliMillePede2::SaveRecordConstraint() {fTreeConstr->Fill(); fRecord->Reset();fCurrRecConstrID++;}
277
278#endif