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