]>
Commit | Line | Data |
---|---|---|
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" | |
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 |