]>
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; |
1f6eade8 | 26 | class TArrayF; |
8a9ab0eb | 27 | |
8a9ab0eb | 28 | |
29 | class 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 | //_____________________________________________________________________________________________ | |
281 | inline void AliMillePede2::ReadRecordData(Long_t recID) {fTreeData->GetEntry(recID); fCurrRecDataID=recID;} | |
282 | ||
283 | //_____________________________________________________________________________________________ | |
284 | inline void AliMillePede2::ReadRecordConstraint(Long_t recID) {fTreeConstr->GetEntry(recID); fCurrRecConstrID=recID;} | |
285 | ||
286 | //_____________________________________________________________________________________________ | |
287 | inline void AliMillePede2::SaveRecordData() {fTreeData->Fill(); fRecord->Reset(); fCurrRecDataID++;} | |
288 | ||
289 | //_____________________________________________________________________________________________ | |
290 | inline void AliMillePede2::SaveRecordConstraint() {fTreeConstr->Fill(); fRecord->Reset();fCurrRecConstrID++;} | |
291 | ||
292 | #endif |