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