]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/STEER/AliMillePede2.h
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'master' into TPCdev
[u/mrichter/AliRoot.git] / STEER / STEER / AliMillePede2.h
CommitLineData
8a9ab0eb 1#ifndef ALIMILLEPEDE2_H
2#define ALIMILLEPEDE2_H
3
de34b538 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
8a9ab0eb 13#include <TObject.h>
7c3070ec 14#include <TString.h>
8a9ab0eb 15#include <TTree.h>
8a9ab0eb 16#include "AliMinResSolve.h"
de34b538 17#include "AliMillePedeRecord.h"
18class TFile;
19class AliMatrixSq;
20class AliSymMatrix;
21class AliRectMatrix;
22class AliMatrixSparse;
8a9ab0eb 23class AliLog;
24class TStopwatch;
e30a812f 25class TArrayL;
1f6eade8 26class TArrayF;
8a9ab0eb 27
8a9ab0eb 28
29class AliMillePede2: public TObject
30{
31 public:
32 //
339fbe23 33 enum {kFailed,kInvert,kNoInversion}; // used global matrix solution methods
9a919f12 34 enum {kFixParID=-1}; // dummy id for fixed param
8a9ab0eb 35 //
36 AliMillePede2();
37 AliMillePede2(const AliMillePede2& src);
38 virtual ~AliMillePede2();
39 AliMillePede2& operator=(const AliMillePede2& ) {printf("Dummy\n"); return *this;}
40 //
9a919f12 41 Int_t InitMille(int nGlo, int nLoc, Int_t lNStdDev=-1,double lResCut=-1., double lResCutInit=-1., const Int_t* regroup=0);
8a9ab0eb 42 //
43 Int_t GetNGloPar() const {return fNGloPar;}
9a919f12 44 Int_t GetNGloParIni() const {return fNGloParIni;}
45 const Int_t* GetRegrouping() const {return fkReGroup;}
8a9ab0eb 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;}
de34b538 52 Int_t GetNLagrangeConstraints() const {return fNLagrangeConstraints;}
8a9ab0eb 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;}
9a919f12 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];}
8a9ab0eb 64 Int_t* GetProcessedPoints() const {return fProcPnt;}
9a919f12 65 Int_t GetParamGrID(Int_t i) const {int ir=GetRGId(i); return ir<=0 ? 0:fParamGrID[ir];}
8a9ab0eb 66 //
67 AliMatrixSq* GetGlobalMatrix() const {return fMatCGlo;}
de34b538 68 AliSymMatrix* GetLocalMatrix() const {return fMatCLoc;}
8a9ab0eb 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;}
de34b538 73 Bool_t* GetIsLinear() const {return fIsLinear;}
9a919f12 74 Double_t GetFinalParam(int i) const {int ir=GetRGId(i); return ir<0 ? 0:fDeltaPar[ir]+fInitPar[ir];}
de34b538 75 Double_t GetFinalError(int i) const {return GetParError(i);}
9a919f12 76 Double_t GetPull(int i) const;
de34b538 77 //
9a919f12 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];}
8a9ab0eb 82 static Bool_t IsGlobalMatSparse() {return fgIsMatGloSparse;}
de34b538 83 static Bool_t IsWeightSigma() {return fgWeightSigma;}
c9f619c7 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;}
8a9ab0eb 89 //
9a919f12 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;}
8a9ab0eb 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;}
de34b538 99 void SetMinPntValid(Int_t n) {fMinPntValid = n>0 ? n:1;}
8a9ab0eb 100 static void SetGlobalMatSparse(Bool_t v=kTRUE) {fgIsMatGloSparse = v;}
de34b538 101 static void SetWeightSigma(Bool_t v=kTRUE) {fgWeightSigma = v;}
8a9ab0eb 102 //
9a919f12 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);
8a9ab0eb 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();
de34b538 111 static void SetInvChol(Bool_t v=kTRUE) {fgInvChol = v;}
8a9ab0eb 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 //
de34b538 118 static Bool_t GetInvChol() {return fgInvChol;}
8a9ab0eb 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;
e30a812f 127 void SetRejRunList(const UInt_t *runs, Int_t nruns);
1f6eade8 128 void SetAccRunList(const UInt_t *runs, Int_t nruns, const Float_t* wghList=0);
129 Bool_t IsRecordAcceptable();
8a9ab0eb 130 //
131 //
132 Int_t SetIterations(double lChi2CutFac);
133
134 //
135 // constraints
339fbe23 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);
8a9ab0eb 138 //
139 // processing of the local measurement
e30a812f 140 void SetRecordRun(Int_t run);
a8c5b94c 141 void SetRecordWeight(double wgh);
8a9ab0eb 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();}
58e8ba3a 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 //
8a9ab0eb 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();
8a9ab0eb 173 AliMillePedeRecord* GetRecord() const {return fRecord;}
e30a812f 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;}
8a9ab0eb 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);}
9a919f12 184 void SetNonLinear(int index, Bool_t v=kTRUE) {int id = GetRGId(index); if (id<0) return; fIsLinear[id] = !v;}
8a9ab0eb 185 //
186 protected:
187 //
188 Int_t LocalFit(double *localParams=0);
cc9bec47 189 Bool_t IsZero(Double_t v,Double_t eps=1e-16) const {return TMath::Abs(v)<eps;}
8a9ab0eb 190 //
191 protected:
192 //
193 Int_t fNLocPar; // number of local parameters
194 Int_t fNGloPar; // number of global parameters
9a919f12 195 Int_t fNGloParIni; // number of global parameters before grouping
8a9ab0eb 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
de34b538 203 Int_t fNLagrangeConstraints; // Number of constraint equations requiring Lagrange multiplier
8a9ab0eb 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 //
de34b538 215 Int_t fNGroupsSet; // number of groups set
d56a0575 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)
8a9ab0eb 221 //
d56a0575 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
8a9ab0eb 225 //
d56a0575 226 Bool_t *fIsLinear; //[fNGloPar] Flag for linear parameters
227 Bool_t *fConstrUsed; //! Flag for used constraints
8a9ab0eb 228 //
d56a0575 229 Int_t *fGlo2CGlo; //[fNGloPar] global ID to compressed ID buffer
230 Int_t *fCGlo2Glo; //[fNGloPar] compressed ID to global ID buffer
8a9ab0eb 231 //
232 // Matrices
233 AliSymMatrix *fMatCLoc; // Matrix C local
234 AliMatrixSq *fMatCGlo; // Matrix C global
de34b538 235 AliRectMatrix *fMatCGloLoc; // Rectangular matrix C g*l
d56a0575 236 Int_t *fFillIndex; //[fNGloPar] auxilary index array for fast matrix fill
237 Double_t *fFillValue; //[fNGloPar] auxilary value array for fast matrix fill
8a9ab0eb 238 //
239 // processed data record bufferization
58e8ba3a 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
8a9ab0eb 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
a8c5b94c 249 Int_t fRecFileStatus; // state of the record file (0-no, 1-read, 2-rw)
8a9ab0eb 250 //
251 TString fConstrRecFName; // Name of File for constraints records
d56a0575 252 TTree *fTreeConstr; //! Tree of constraint records
253 TFile *fConsRecFile; //! File of processed constraints records
8a9ab0eb 254 Long_t fCurrRecDataID; // ID of the current data record
255 Long_t fCurrRecConstrID; // ID of the current constraint record
de34b538 256 Bool_t fLocFitAdd; // Add contribution of carrent track (and not eliminate it)
c9f619c7 257 Bool_t fUseRecordWeight; // force or ignore the record weight
258 Int_t fMinRecordLength; // ignore shorter records
e30a812f 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)
1f6eade8 263 TArrayF* fAccRunListWgh; // optional weights for data of accepted runs (if any)
264 Double_t fRunWgh; // run weight
c9f619c7 265 Double_t fWghScl[2]; // optional rescaling for odd/even residual weights (see its usage in LocalFit)
9a919f12 266 const Int_t* fkReGroup; // optional regrouping of parameters wrt ID's from the records
8a9ab0eb 267 //
de34b538 268 static Bool_t fgInvChol; // Invert global matrix in Cholesky solver
269 static Bool_t fgWeightSigma; // weight parameter constraint by statistics
8a9ab0eb 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 //
d56a0575 277 ClassDef(AliMillePede2,1)
8a9ab0eb 278};
279
280//_____________________________________________________________________________________________
281inline void AliMillePede2::ReadRecordData(Long_t recID) {fTreeData->GetEntry(recID); fCurrRecDataID=recID;}
282
283//_____________________________________________________________________________________________
284inline void AliMillePede2::ReadRecordConstraint(Long_t recID) {fTreeConstr->GetEntry(recID); fCurrRecConstrID=recID;}
285
286//_____________________________________________________________________________________________
287inline void AliMillePede2::SaveRecordData() {fTreeData->Fill(); fRecord->Reset(); fCurrRecDataID++;}
288
289//_____________________________________________________________________________________________
290inline void AliMillePede2::SaveRecordConstraint() {fTreeConstr->Fill(); fRecord->Reset();fCurrRecConstrID++;}
291
292#endif