]>
Commit | Line | Data |
---|---|---|
de34b538 | 1 | /**********************************************************************************************/ |
2 | /* General class for alignment with large number of degrees of freedom */ | |
3 | /* Based on the original milliped2 by Volker Blobel */ | |
4 | /* http://www.desy.de/~blobel/mptalks.html */ | |
5 | /* */ | |
6 | /* Author: ruben.shahoyan@cern.ch */ | |
7 | /* */ | |
8 | /**********************************************************************************************/ | |
9 | ||
8a9ab0eb | 10 | #include "AliMillePede2.h" |
11 | #include "AliLog.h" | |
12 | #include <TStopwatch.h> | |
de34b538 | 13 | #include <TFile.h> |
14 | #include <TMath.h> | |
7c3070ec | 15 | #include <TVectorD.h> |
de34b538 | 16 | #include "AliMatrixSq.h" |
17 | #include "AliSymMatrix.h" | |
18 | #include "AliRectMatrix.h" | |
19 | #include "AliMatrixSparse.h" | |
8a9ab0eb | 20 | |
8a9ab0eb | 21 | |
de34b538 | 22 | using namespace std; |
8a9ab0eb | 23 | |
8a9ab0eb | 24 | |
8a9ab0eb | 25 | ClassImp(AliMillePede2) |
26 | ||
de34b538 | 27 | Bool_t AliMillePede2::fgInvChol = kTRUE; // Invert global matrix with Cholesky solver |
28 | Bool_t AliMillePede2::fgWeightSigma = kTRUE; // weight local constraint by module statistics | |
29 | Bool_t AliMillePede2::fgIsMatGloSparse = kFALSE; // use faster dense matrix by default | |
30 | Int_t AliMillePede2::fgMinResCondType = 1; // Jacoby preconditioner by default | |
31 | Double_t AliMillePede2::fgMinResTol = 1.e-11; // default tolerance | |
32 | Int_t AliMillePede2::fgMinResMaxIter = 10000; // default max number of iterations | |
8a9ab0eb | 33 | Int_t AliMillePede2::fgIterSol = AliMinResSolve::kSolMinRes; // default iterative solver |
de34b538 | 34 | Int_t AliMillePede2::fgNKrylovV = 240; // default number of Krylov vectors to keep |
8a9ab0eb | 35 | |
36 | //_____________________________________________________________________________________________ | |
37 | AliMillePede2::AliMillePede2() | |
38 | : fNLocPar(0), | |
39 | fNGloPar(0), | |
40 | fNGloSize(0), | |
41 | // | |
42 | fNLocEquations(0), | |
43 | fIter(0), | |
44 | fMaxIter(10), | |
45 | fNStdDev(3), | |
46 | fNGloConstraints(0), | |
de34b538 | 47 | fNLagrangeConstraints(0), |
8a9ab0eb | 48 | fNLocFits(0), |
49 | fNLocFitsRejected(0), | |
50 | fNGloFix(0), | |
51 | fGloSolveStatus(gkFailed), | |
52 | // | |
53 | fChi2CutFactor(1.), | |
54 | fChi2CutRef(1.), | |
55 | fResCutInit(100.), | |
56 | fResCut(100.), | |
de34b538 | 57 | fMinPntValid(1), |
8a9ab0eb | 58 | // |
de34b538 | 59 | fNGroupsSet(0), |
60 | fParamGrID(0), | |
8a9ab0eb | 61 | fProcPnt(0), |
62 | fVecBLoc(0), | |
63 | fDiagCGlo(0), | |
64 | fVecBGlo(0), | |
65 | fInitPar(0), | |
66 | fDeltaPar(0), | |
67 | fSigmaPar(0), | |
68 | fIsLinear(0), | |
69 | fConstrUsed(0), | |
70 | // | |
71 | fGlo2CGlo(0), | |
72 | fCGlo2Glo(0), | |
73 | // | |
74 | fMatCLoc(0), | |
75 | fMatCGlo(0), | |
76 | fMatCGloLoc(0), | |
de34b538 | 77 | // |
78 | fFillIndex(0), | |
79 | fFillValue(0), | |
80 | // | |
8a9ab0eb | 81 | fDataRecFName("/tmp/mp2_data_records.root"), |
82 | fRecord(0), | |
83 | fDataRecFile(0), | |
84 | fTreeData(0), | |
a8c5b94c | 85 | fRecFileStatus(0), |
8a9ab0eb | 86 | // |
87 | fConstrRecFName("/tmp/mp2_constraints_records.root"), | |
88 | fTreeConstr(0), | |
89 | fConsRecFile(0), | |
90 | fCurrRecDataID(0), | |
de34b538 | 91 | fCurrRecConstrID(0), |
92 | fLocFitAdd(kTRUE) | |
8a9ab0eb | 93 | {} |
94 | ||
95 | //_____________________________________________________________________________________________ | |
96 | AliMillePede2::AliMillePede2(const AliMillePede2& src) : | |
97 | TObject(src),fNLocPar(0),fNGloPar(0),fNGloSize(0),fNLocEquations(0),fIter(0), | |
de34b538 | 98 | fMaxIter(10),fNStdDev(3),fNGloConstraints(0),fNLagrangeConstraints(0), |
99 | fNLocFits(0),fNLocFitsRejected(0), | |
8a9ab0eb | 100 | fNGloFix(0),fGloSolveStatus(0),fChi2CutFactor(0),fChi2CutRef(0),fResCutInit(0), |
de34b538 | 101 | fResCut(0),fMinPntValid(1),fNGroupsSet(0),fParamGrID(0),fProcPnt(0),fVecBLoc(0),fDiagCGlo(0),fVecBGlo(0), |
8a9ab0eb | 102 | fInitPar(0),fDeltaPar(0),fSigmaPar(0),fIsLinear(0),fConstrUsed(0),fGlo2CGlo(0),fCGlo2Glo(0), |
de34b538 | 103 | fMatCLoc(0),fMatCGlo(0),fMatCGloLoc(0),fFillIndex(0),fFillValue(0), |
104 | fDataRecFName(0),fRecord(0),fDataRecFile(0), | |
a8c5b94c | 105 | fTreeData(0),fRecFileStatus(0),fConstrRecFName(0),fTreeConstr(0),fConsRecFile(0),fCurrRecDataID(0), |
de34b538 | 106 | fCurrRecConstrID(0),fLocFitAdd(kTRUE) |
8a9ab0eb | 107 | {printf("Dummy\n");} |
108 | ||
109 | //_____________________________________________________________________________________________ | |
110 | AliMillePede2::~AliMillePede2() | |
111 | { | |
112 | CloseDataRecStorage(); | |
113 | CloseConsRecStorage(); | |
114 | // | |
de34b538 | 115 | delete[] fParamGrID; |
8a9ab0eb | 116 | delete[] fProcPnt; |
117 | delete[] fVecBLoc; | |
118 | delete[] fDiagCGlo; | |
119 | delete[] fVecBGlo; | |
120 | delete[] fInitPar; | |
121 | delete[] fDeltaPar; | |
122 | delete[] fSigmaPar; | |
123 | delete[] fGlo2CGlo; | |
124 | delete[] fCGlo2Glo; | |
125 | delete[] fIsLinear; | |
126 | delete[] fConstrUsed; | |
de34b538 | 127 | delete[] fFillIndex; |
128 | delete[] fFillValue; | |
8a9ab0eb | 129 | // |
130 | delete fRecord; | |
131 | delete fMatCLoc; | |
132 | delete fMatCGlo; | |
133 | delete fMatCGloLoc; | |
134 | } | |
135 | ||
136 | //_____________________________________________________________________________________________ | |
137 | Int_t AliMillePede2::InitMille(int nGlo, int nLoc, int lNStdDev,double lResCut, double lResCutInit) | |
138 | { | |
139 | // | |
140 | if (nLoc>0) fNLocPar = nLoc; | |
141 | if (nGlo>0) fNGloPar = nGlo; | |
142 | if (lResCutInit>0) fResCutInit = lResCutInit; | |
143 | if (lResCut>0) fResCut = lResCut; | |
144 | if (lNStdDev>0) fNStdDev = lNStdDev; | |
145 | // | |
146 | fNGloSize = fNGloPar; | |
147 | // | |
148 | try { | |
149 | // | |
150 | if (fgIsMatGloSparse) {fMatCGlo = new AliMatrixSparse(fNGloPar); fMatCGlo->SetSymmetric(kTRUE);} | |
151 | else fMatCGlo = new AliSymMatrix(fNGloPar); | |
152 | // | |
de34b538 | 153 | fFillIndex = new Int_t[fNGloPar]; |
154 | fFillValue = new Double_t[fNGloPar]; | |
155 | // | |
8a9ab0eb | 156 | fMatCLoc = new AliSymMatrix(fNLocPar); |
de34b538 | 157 | fMatCGloLoc = new AliRectMatrix(fNGloPar,fNLocPar); |
8a9ab0eb | 158 | // |
de34b538 | 159 | fParamGrID = new Int_t[fNGloPar]; |
8a9ab0eb | 160 | fProcPnt = new Int_t[fNGloPar]; |
161 | fVecBLoc = new Double_t[fNLocPar]; | |
162 | fDiagCGlo = new Double_t[fNGloPar]; | |
163 | // | |
164 | fInitPar = new Double_t[fNGloPar]; | |
165 | fDeltaPar = new Double_t[fNGloPar]; | |
166 | fSigmaPar = new Double_t[fNGloPar]; | |
167 | fIsLinear = new Bool_t[fNGloPar]; | |
168 | // | |
169 | fGlo2CGlo = new Int_t[fNGloPar]; | |
170 | fCGlo2Glo = new Int_t[fNGloPar]; | |
171 | } | |
172 | catch(bad_alloc&) { | |
173 | AliInfo(Form("Failed to allocate the memory for %d global and %d local parameters",fNGloPar,fNLocPar)); | |
174 | return 0; | |
175 | } | |
176 | // | |
177 | memset(fVecBLoc ,0,fNLocPar*sizeof(Double_t)); | |
178 | memset(fDiagCGlo ,0,fNGloPar*sizeof(Double_t)); | |
179 | memset(fInitPar ,0,fNGloPar*sizeof(Double_t)); | |
180 | memset(fDeltaPar ,0,fNGloPar*sizeof(Double_t)); | |
181 | memset(fSigmaPar ,0,fNGloPar*sizeof(Double_t)); | |
182 | memset(fProcPnt ,0,fNGloPar*sizeof(Int_t)); | |
183 | // | |
184 | for (int i=fNGloPar;i--;) { | |
de34b538 | 185 | fGlo2CGlo[i] = fCGlo2Glo[i] = -1; |
8a9ab0eb | 186 | fIsLinear[i] = kTRUE; |
de34b538 | 187 | fParamGrID[i] = -1; |
8a9ab0eb | 188 | } |
189 | // | |
190 | return 1; | |
191 | } | |
192 | ||
8a9ab0eb | 193 | //_____________________________________________________________________________________________ |
194 | Bool_t AliMillePede2::ImposeDataRecFile(const char* fname) | |
195 | { | |
196 | CloseDataRecStorage(); | |
197 | SetDataRecFName(fname); | |
198 | return InitDataRecStorage(kTRUE); // open in read mode | |
199 | } | |
200 | ||
201 | //_____________________________________________________________________________________________ | |
202 | Bool_t AliMillePede2::ImposeConsRecFile(const char* fname) | |
203 | { | |
204 | CloseConsRecStorage(); | |
205 | SetConsRecFName(fname); | |
206 | return InitConsRecStorage(kTRUE); // open in read mode | |
207 | } | |
208 | ||
209 | //_____________________________________________________________________________________________ | |
210 | Bool_t AliMillePede2::InitDataRecStorage(Bool_t read) | |
211 | { | |
212 | // initialize the buffer for processed measurements records | |
213 | // | |
214 | if (fDataRecFile) {AliInfo("Data Records File is already initialized"); return kFALSE;} | |
215 | // | |
216 | if (!fRecord) fRecord = new AliMillePedeRecord(); | |
217 | // | |
218 | fDataRecFile = TFile::Open(GetDataRecFName(),read ? "":"recreate"); | |
a8c5b94c | 219 | if (!fDataRecFile) {AliFatal(Form("Failed to initialize data records file %s",GetDataRecFName())); return kFALSE;} |
8a9ab0eb | 220 | // |
221 | AliInfo(Form("File %s used for derivatives records",GetDataRecFName())); | |
222 | if (read) { | |
223 | fTreeData = (TTree*)fDataRecFile->Get("AliMillePedeRecords_Data"); | |
a8c5b94c | 224 | if (!fTreeData) {AliFatal(Form("Did not find data records tree in %s",GetDataRecFName())); return kFALSE;} |
8a9ab0eb | 225 | fTreeData->SetBranchAddress("Record_Data",&fRecord); |
226 | AliInfo(Form("Found %d derivatives records",fTreeData->GetEntries())); | |
227 | } | |
228 | else { | |
229 | fTreeData = new TTree("AliMillePedeRecords_Data","Data Records for AliMillePede2"); | |
230 | fTreeData->Branch("Record_Data","AliMillePedeRecord",&fRecord,32000,99); | |
231 | } | |
232 | fCurrRecDataID = -1; | |
a8c5b94c | 233 | fRecFileStatus = read ? 1:2; |
8a9ab0eb | 234 | // |
235 | return kTRUE; | |
236 | } | |
237 | ||
238 | //_____________________________________________________________________________________________ | |
239 | Bool_t AliMillePede2::InitConsRecStorage(Bool_t read) | |
240 | { | |
241 | // initialize the buffer for processed measurements records | |
242 | // | |
243 | if (fConsRecFile) {AliInfo("Constraints Records File is already initialized"); return kFALSE;} | |
244 | // | |
245 | if (!fRecord) fRecord = new AliMillePedeRecord(); | |
246 | // | |
247 | fConsRecFile = TFile::Open(GetConsRecFName(),read ? "":"recreate"); | |
248 | if (!fConsRecFile) {AliInfo(Form("Failed to initialize constraints records file %s",GetConsRecFName())); return kFALSE;} | |
249 | // | |
250 | AliInfo(Form("File %s used for constraints records",GetConsRecFName())); | |
251 | if (read) { | |
252 | fTreeConstr = (TTree*)fConsRecFile->Get("AliMillePedeRecords_Constraints"); | |
253 | if (!fTreeConstr) {AliInfo(Form("Did not find constraints records tree in %s",GetConsRecFName())); return kFALSE;} | |
254 | fTreeConstr->SetBranchAddress("Record_Constraints",&fRecord); | |
255 | AliInfo(Form("Found %d constraints records",fTreeConstr->GetEntries())); | |
256 | // | |
257 | } | |
258 | else { | |
259 | // | |
260 | fTreeConstr = new TTree("AliMillePedeRecords_Constraints","Constraints Records for AliMillePede2"); | |
261 | fTreeConstr->Branch("Record_Constraints","AliMillePedeRecord",&fRecord,32000,99); | |
262 | } | |
263 | fCurrRecConstrID = -1; | |
264 | // | |
265 | return kTRUE; | |
266 | } | |
267 | ||
268 | //_____________________________________________________________________________________________ | |
269 | void AliMillePede2::CloseDataRecStorage() | |
270 | { | |
271 | if (fTreeData) { | |
272 | if (fDataRecFile->IsWritable()) { | |
273 | fDataRecFile->cd(); | |
274 | fTreeData->Write(); | |
275 | } | |
276 | delete fTreeData; | |
277 | fTreeData = 0; | |
278 | fDataRecFile->Close(); | |
279 | delete fDataRecFile; | |
280 | fDataRecFile = 0; | |
281 | } | |
a8c5b94c | 282 | fRecFileStatus = 0; |
8a9ab0eb | 283 | // |
284 | } | |
285 | ||
286 | //_____________________________________________________________________________________________ | |
287 | void AliMillePede2::CloseConsRecStorage() | |
288 | { | |
289 | if (fTreeConstr) { | |
290 | if (fConsRecFile->IsWritable()) { | |
291 | fConsRecFile->cd(); | |
292 | fTreeConstr->Write(); | |
293 | } | |
294 | delete fTreeConstr; | |
295 | fTreeConstr = 0; | |
296 | fConsRecFile->Close(); | |
297 | delete fConsRecFile; | |
298 | fConsRecFile = 0; | |
299 | } | |
300 | // | |
301 | } | |
302 | ||
303 | //_____________________________________________________________________________________________ | |
304 | Bool_t AliMillePede2::ReadNextRecordData() | |
305 | { | |
306 | // read next data record (if any) | |
307 | if (!fTreeData || ++fCurrRecDataID >= fTreeData->GetEntries()) { fCurrRecDataID--; return kFALSE;} | |
308 | fTreeData->GetEntry(fCurrRecDataID); | |
309 | return kTRUE; | |
310 | } | |
311 | ||
312 | //_____________________________________________________________________________________________ | |
313 | Bool_t AliMillePede2::ReadNextRecordConstraint() | |
314 | { | |
315 | // read next constraint record (if any) | |
316 | if (!fTreeConstr || ++fCurrRecConstrID >= fTreeConstr->GetEntries()) { fCurrRecConstrID--; return kFALSE;} | |
317 | fTreeConstr->GetEntry(fCurrRecConstrID); | |
318 | return kTRUE; | |
319 | } | |
320 | ||
a8c5b94c | 321 | //_____________________________________________________________________________________________ |
322 | void AliMillePede2::SetRecordWeight(double wgh) | |
323 | { | |
324 | if (fRecFileStatus<2) InitDataRecStorage(); // create a buffer to store the data | |
325 | fRecord->SetWeight(wgh); | |
326 | } | |
327 | ||
8a9ab0eb | 328 | //_____________________________________________________________________________________________ |
329 | void AliMillePede2::SetLocalEquation(double *dergb, double *derlc, double lMeas, double lSigma) | |
330 | { | |
a8c5b94c | 331 | if (fRecFileStatus<2) InitDataRecStorage(); // create a buffer to store the data |
8a9ab0eb | 332 | // |
333 | // write data of single measurement | |
334 | if (lSigma<=0.0) { // If parameter is fixed, then no equation | |
335 | for (int i=fNLocPar; i--;) derlc[i] = 0.0; | |
336 | for (int i=fNGloPar; i--;) dergb[i] = 0.0; | |
337 | return; | |
338 | } | |
339 | // | |
340 | fRecord->AddResidual(lMeas); | |
341 | // | |
342 | // Retrieve local param interesting indices | |
cc9bec47 | 343 | for (int i=0;i<fNLocPar;i++) if (!IsZero(derlc[i])) {fRecord->AddIndexValue(i,derlc[i]); derlc[i] = 0.0;} |
8a9ab0eb | 344 | // |
345 | fRecord->AddWeight( 1.0/lSigma/lSigma ); | |
346 | // | |
347 | // Idem for global parameters | |
cc9bec47 | 348 | for (int i=0;i<fNGloPar;i++) if (!IsZero(dergb[i])) { |
de34b538 | 349 | fRecord->AddIndexValue(i,dergb[i]); dergb[i] = 0.0; |
350 | fRecord->MarkGroup(fParamGrID[i]); | |
351 | } | |
8a9ab0eb | 352 | // |
353 | } | |
354 | ||
355 | //_____________________________________________________________________________________________ | |
356 | void AliMillePede2::SetLocalEquation(int *indgb, double *dergb, int ngb, int *indlc, | |
357 | double *derlc,int nlc,double lMeas,double lSigma) | |
358 | { | |
359 | // write data of single measurement | |
360 | if (lSigma<=0.0) { // If parameter is fixed, then no equation | |
361 | for (int i=nlc;i--;) derlc[i] = 0.0; | |
362 | for (int i=ngb;i--;) dergb[i] = 0.0; | |
363 | return; | |
364 | } | |
365 | // | |
a8c5b94c | 366 | if (fRecFileStatus<2) InitDataRecStorage(); // create a buffer to store the data |
8a9ab0eb | 367 | // |
368 | fRecord->AddResidual(lMeas); | |
369 | // | |
370 | // Retrieve local param interesting indices | |
cc9bec47 | 371 | for (int i=0;i<nlc;i++) if (!IsZero(derlc[i])) {fRecord->AddIndexValue(indlc[i],derlc[i]); derlc[i]=0.; indlc[i]=0;} |
8a9ab0eb | 372 | // |
373 | fRecord->AddWeight( 1./lSigma/lSigma ); | |
374 | // | |
375 | // Idem for global parameters | |
cc9bec47 | 376 | for (int i=0;i<ngb;i++) if (!IsZero(dergb[i])) {fRecord->AddIndexValue(indgb[i],dergb[i]); dergb[i]=0.; indgb[i]=0;} |
8a9ab0eb | 377 | // |
378 | } | |
379 | ||
de34b538 | 380 | |
8a9ab0eb | 381 | //_____________________________________________________________________________________________ |
de34b538 | 382 | void AliMillePede2::SetGlobalConstraint(double *dergb, double val, double sigma) |
8a9ab0eb | 383 | { |
384 | // Define a constraint equation. | |
385 | if (!fConsRecFile || !fConsRecFile->IsWritable()) InitConsRecStorage(); // create a buffer to store the data | |
386 | // | |
387 | fRecord->Reset(); | |
388 | fRecord->AddResidual(val); | |
de34b538 | 389 | fRecord->AddWeight(sigma); |
cc9bec47 | 390 | for (int i=0; i<fNGloPar; i++) if (!IsZero(dergb[i])) fRecord->AddIndexValue(i,dergb[i]); |
8a9ab0eb | 391 | fNGloConstraints++; |
cc9bec47 | 392 | if (IsZero(sigma)) fNLagrangeConstraints++; |
de34b538 | 393 | // printf("NewConstraint:\n"); fRecord->Print(); //RRR |
8a9ab0eb | 394 | SaveRecordConstraint(); |
395 | // | |
396 | } | |
397 | ||
398 | //_____________________________________________________________________________________________ | |
de34b538 | 399 | void AliMillePede2::SetGlobalConstraint(const int *indgb, double *dergb, int ngb, double val,double sigma) |
8a9ab0eb | 400 | { |
401 | // Define a constraint equation. | |
402 | if (!fConsRecFile || !fConsRecFile->IsWritable()) InitConsRecStorage(); // create a buffer to store the data | |
403 | fRecord->Reset(); | |
404 | fRecord->AddResidual(val); | |
de34b538 | 405 | fRecord->AddWeight(sigma); // dummy |
cc9bec47 | 406 | for (int i=0; i<ngb; i++) if (!IsZero(dergb[i])) fRecord->AddIndexValue(indgb[i],dergb[i]); |
8a9ab0eb | 407 | fNGloConstraints++; |
cc9bec47 | 408 | if (IsZero(sigma)) fNLagrangeConstraints++; |
8a9ab0eb | 409 | SaveRecordConstraint(); |
410 | // | |
411 | } | |
412 | ||
413 | //_____________________________________________________________________________________________ | |
414 | Int_t AliMillePede2::LocalFit(double *localParams) | |
415 | { | |
416 | /* | |
417 | Perform local parameters fit once all the local equations have been set | |
418 | ----------------------------------------------------------- | |
419 | localParams = (if !=0) will contain the fitted track parameters and | |
420 | related errors | |
421 | */ | |
de34b538 | 422 | static int nrefSize = 0; |
423 | // static TArrayI refLoc,refGlo,nrefLoc,nrefGlo; | |
424 | static Int_t *refLoc=0,*refGlo=0,*nrefLoc=0,*nrefGlo=0; | |
8a9ab0eb | 425 | int nPoints = 0; |
426 | // | |
de34b538 | 427 | AliSymMatrix &matCLoc = *fMatCLoc; |
428 | AliMatrixSq &matCGlo = *fMatCGlo; | |
429 | AliRectMatrix &matCGloLoc = *fMatCGloLoc; | |
8a9ab0eb | 430 | // |
431 | memset(fVecBLoc,0,fNLocPar*sizeof(double)); | |
de34b538 | 432 | matCLoc.Reset(); |
8a9ab0eb | 433 | // |
434 | int cnt=0; | |
435 | int recSz = fRecord->GetSize(); | |
436 | // | |
437 | while(cnt<recSz) { // Transfer the measurement records to matrices | |
438 | // | |
439 | // extract addresses of residual, weight and pointers on local and global derivatives for each point | |
de34b538 | 440 | if (nrefSize<=nPoints) { |
441 | int *tmpA = 0; | |
442 | nrefSize = 2*(nPoints+1); | |
443 | tmpA = refLoc; refLoc = new Int_t[nrefSize]; if (tmpA) memcpy(refLoc,tmpA,nPoints*sizeof(int)); | |
444 | tmpA = refGlo; refGlo = new Int_t[nrefSize]; if (tmpA) memcpy(refGlo,tmpA,nPoints*sizeof(int)); | |
445 | tmpA = nrefLoc; nrefLoc = new Int_t[nrefSize]; if (tmpA) memcpy(nrefLoc,tmpA,nPoints*sizeof(int)); | |
446 | tmpA = nrefGlo; nrefGlo = new Int_t[nrefSize]; if (tmpA) memcpy(nrefGlo,tmpA,nPoints*sizeof(int)); | |
8a9ab0eb | 447 | } |
448 | // | |
de34b538 | 449 | refLoc[nPoints] = ++cnt; |
8a9ab0eb | 450 | int nLoc = 0; |
451 | while(!fRecord->IsWeight(cnt)) {nLoc++; cnt++;} | |
de34b538 | 452 | nrefLoc[nPoints] = nLoc; |
8a9ab0eb | 453 | // |
de34b538 | 454 | refGlo[nPoints] = ++cnt; |
8a9ab0eb | 455 | int nGlo = 0; |
456 | while(!fRecord->IsResidual(cnt) && cnt<recSz) {nGlo++; cnt++;} | |
de34b538 | 457 | nrefGlo[nPoints] = nGlo; |
8a9ab0eb | 458 | // |
459 | nPoints++; | |
460 | } | |
461 | double vl; | |
462 | // | |
a8c5b94c | 463 | double gloWgh = fRecord->GetWeight(); // global weight for this set |
7c3070ec | 464 | Int_t maxLocUsed = 0; |
465 | // | |
8a9ab0eb | 466 | for (int ip=nPoints;ip--;) { // Transfer the measurement records to matrices |
de34b538 | 467 | double resid = fRecord->GetValue( refLoc[ip]-1 ); |
a8c5b94c | 468 | double weight = fRecord->GetValue( refGlo[ip]-1 )*gloWgh; |
de34b538 | 469 | double *derLoc = fRecord->GetValue()+refLoc[ip]; |
470 | double *derGlo = fRecord->GetValue()+refGlo[ip]; | |
471 | int *indLoc = fRecord->GetIndex()+refLoc[ip]; | |
472 | int *indGlo = fRecord->GetIndex()+refGlo[ip]; | |
473 | // | |
a8c5b94c | 474 | for (int i=nrefGlo[ip];i--;) { // suppress the global part (only relevant with iterations) |
8a9ab0eb | 475 | int iID = indGlo[i]; // Global param indice |
de34b538 | 476 | if (fSigmaPar[iID]<=0.) continue; // fixed parameter RRRCheck |
8a9ab0eb | 477 | if (fIsLinear[iID]) resid -= derGlo[i]*(fInitPar[iID]+fDeltaPar[iID]); // linear parameter |
478 | else resid -= derGlo[i]*fDeltaPar[iID]; // nonlinear parameter | |
479 | } | |
480 | // | |
de34b538 | 481 | // Symmetric matrix, don't bother j>i coeffs |
482 | for (int i=nrefLoc[ip];i--;) { // Fill local matrix and vector | |
483 | fVecBLoc[ indLoc[i] ] += weight*resid*derLoc[i]; | |
7c3070ec | 484 | if (indLoc[i]>maxLocUsed) maxLocUsed = indLoc[i]; |
de34b538 | 485 | for (int j=i+1;j--;) matCLoc(indLoc[i] ,indLoc[j]) += weight*derLoc[i]*derLoc[j]; |
7c3070ec | 486 | } |
8a9ab0eb | 487 | // |
488 | } // end of the transfer of the measurement record to matrices | |
489 | // | |
7c3070ec | 490 | matCLoc.SetSizeUsed(++maxLocUsed); // data with B=0 may use less than declared nLocals |
491 | // | |
8a9ab0eb | 492 | // first try to solve by faster Cholesky decomposition, then by Gaussian elimination |
de34b538 | 493 | if (!matCLoc.SolveChol(fVecBLoc,kTRUE)) { |
a8c5b94c | 494 | AliInfo("Failed to solve locals by Cholesky, trying Gaussian Elimination"); |
de34b538 | 495 | if (!matCLoc.SolveSpmInv(fVecBLoc,kTRUE)) { |
8a9ab0eb | 496 | AliInfo("Failed to solve locals by Gaussian Elimination, skip..."); |
de34b538 | 497 | matCLoc.Print("d"); |
8a9ab0eb | 498 | return 0; // failed to solve |
499 | } | |
500 | } | |
501 | // | |
502 | // If requested, store the track params and errors | |
de34b538 | 503 | // printf("locfit: "); for (int i=0;i<fNLocPar;i++) printf("%+e |",fVecBLoc[i]); printf("\n"); |
7c3070ec | 504 | if (localParams) for (int i=maxLocUsed; i--;) { |
8a9ab0eb | 505 | localParams[2*i] = fVecBLoc[i]; |
de34b538 | 506 | localParams[2*i+1] = TMath::Sqrt(TMath::Abs(matCLoc.QueryDiag(i))); |
8a9ab0eb | 507 | } |
508 | // | |
509 | float lChi2 = 0; | |
510 | int nEq = 0; | |
511 | // | |
512 | for (int ip=nPoints;ip--;) { // Calculate residuals | |
de34b538 | 513 | double resid = fRecord->GetValue( refLoc[ip]-1 ); |
a8c5b94c | 514 | double weight = fRecord->GetValue( refGlo[ip]-1 )*gloWgh; |
de34b538 | 515 | double *derLoc = fRecord->GetValue()+refLoc[ip]; |
516 | double *derGlo = fRecord->GetValue()+refGlo[ip]; | |
517 | int *indLoc = fRecord->GetIndex()+refLoc[ip]; | |
518 | int *indGlo = fRecord->GetIndex()+refGlo[ip]; | |
8a9ab0eb | 519 | // |
520 | // Suppress local and global contribution in residuals; | |
de34b538 | 521 | for (int i=nrefLoc[ip];i--;) resid -= derLoc[i]*fVecBLoc[ indLoc[i] ]; // local part |
8a9ab0eb | 522 | // |
de34b538 | 523 | for (int i=nrefGlo[ip];i--;) { // global part |
8a9ab0eb | 524 | int iID = indGlo[i]; |
525 | if ( fSigmaPar[iID] <= 0.) continue; // fixed parameter RRRCheck | |
526 | if (fIsLinear[iID]) resid -= derGlo[i]*(fInitPar[iID]+fDeltaPar[iID]); // linear parameter | |
527 | else resid -= derGlo[i]*fDeltaPar[iID]; // nonlinear parameter | |
528 | } | |
529 | // | |
530 | // reject the track if the residual is too large (outlier) | |
de34b538 | 531 | double absres = TMath::Abs(resid); |
532 | if ( (absres >= fResCutInit && fIter ==1 ) || | |
533 | (absres >= fResCut && fIter > 1)) { | |
534 | if (fLocFitAdd) fNLocFitsRejected++; | |
cc9bec47 | 535 | // printf("reject res %5ld %+e\n",fCurrRecDataID,resid); |
8a9ab0eb | 536 | return 0; |
537 | } | |
538 | // | |
539 | lChi2 += weight*resid*resid ; // total chi^2 | |
540 | nEq++; // number of equations | |
541 | } // end of Calculate residuals | |
542 | // | |
a8c5b94c | 543 | lChi2 /= gloWgh; |
7c3070ec | 544 | int nDoF = nEq-maxLocUsed; |
8a9ab0eb | 545 | lChi2 = (nDoF>0) ? lChi2/nDoF : 0; // Chi^2/dof |
546 | // | |
547 | if (fNStdDev != 0 && nDoF>0 && lChi2 > Chi2DoFLim(fNStdDev,nDoF)*fChi2CutFactor) { // check final chi2 | |
de34b538 | 548 | if (fLocFitAdd) fNLocFitsRejected++; |
cc9bec47 | 549 | // printf("reject chi2 %5ld: %+e\n",fCurrRecDataID, lChi2); |
8a9ab0eb | 550 | return 0; |
551 | } | |
552 | // | |
de34b538 | 553 | if (fLocFitAdd) { |
554 | fNLocFits++; | |
555 | fNLocEquations += nEq; | |
556 | } | |
557 | else { | |
558 | fNLocFits--; | |
559 | fNLocEquations -= nEq; | |
560 | } | |
8a9ab0eb | 561 | // |
562 | // local operations are finished, track is accepted | |
563 | // We now update the global parameters (other matrices) | |
564 | // | |
565 | int nGloInFit = 0; | |
566 | // | |
567 | for (int ip=nPoints;ip--;) { // Update matrices | |
de34b538 | 568 | double resid = fRecord->GetValue( refLoc[ip]-1 ); |
a8c5b94c | 569 | double weight = fRecord->GetValue( refGlo[ip]-1 )*gloWgh; |
de34b538 | 570 | double *derLoc = fRecord->GetValue()+refLoc[ip]; |
571 | double *derGlo = fRecord->GetValue()+refGlo[ip]; | |
572 | int *indLoc = fRecord->GetIndex()+refLoc[ip]; | |
573 | int *indGlo = fRecord->GetIndex()+refGlo[ip]; | |
574 | // | |
575 | for (int i=nrefGlo[ip];i--;) { // suppress the global part | |
8a9ab0eb | 576 | int iID = indGlo[i]; // Global param indice |
577 | if ( fSigmaPar[iID] <= 0.) continue; // fixed parameter RRRCheck | |
de34b538 | 578 | if (fIsLinear[iID]) resid -= derGlo[i]*(fInitPar[iID]+fDeltaPar[iID]); // linear parameter |
579 | else resid -= derGlo[i]*fDeltaPar[iID]; // nonlinear parameter | |
8a9ab0eb | 580 | } |
581 | // | |
de34b538 | 582 | for (int ig=nrefGlo[ip];ig--;) { |
8a9ab0eb | 583 | int iIDg = indGlo[ig]; // Global param indice (the matrix line) |
584 | if ( fSigmaPar[iIDg] <= 0.) continue; // fixed parameter RRRCheck | |
de34b538 | 585 | if (fLocFitAdd) fVecBGlo[ iIDg ] += weight*resid*derGlo[ig]; //!!! |
586 | else fVecBGlo[ iIDg ] -= weight*resid*derGlo[ig]; //!!! | |
8a9ab0eb | 587 | // |
588 | // First of all, the global/global terms (exactly like local matrix) | |
de34b538 | 589 | int nfill = 0; |
590 | for (int jg=ig+1;jg--;) { // matCGlo is symmetric by construction | |
8a9ab0eb | 591 | int jIDg = indGlo[jg]; |
592 | if ( fSigmaPar[jIDg] <= 0.) continue; // fixed parameter RRRCheck | |
cc9bec47 | 593 | if ( !IsZero(vl = weight*derGlo[ig]*derGlo[jg]) ) { |
de34b538 | 594 | fFillIndex[nfill] = jIDg; |
595 | fFillValue[nfill++] = fLocFitAdd ? vl:-vl; | |
596 | } | |
597 | } | |
598 | if (nfill) matCGlo.AddToRow(iIDg,fFillValue,fFillIndex,nfill); | |
8a9ab0eb | 599 | // |
600 | // Now we have also rectangular matrices containing global/local terms. | |
601 | int iCIDg = fGlo2CGlo[iIDg]; // compressed Index of index | |
602 | if (iCIDg == -1) { | |
de34b538 | 603 | Double_t *rowGL = matCGloLoc(nGloInFit); |
7c3070ec | 604 | for (int k=maxLocUsed;k--;) rowGL[k] = 0.0; // reset the row |
8a9ab0eb | 605 | iCIDg = fGlo2CGlo[iIDg] = nGloInFit; |
606 | fCGlo2Glo[nGloInFit++] = iIDg; | |
607 | } | |
608 | // | |
de34b538 | 609 | Double_t *rowGLIDg = matCGloLoc(iCIDg); |
610 | for (int il=nrefLoc[ip];il--;) rowGLIDg[ indLoc[il] ] += weight*derGlo[ig]*derLoc[il]; | |
611 | fProcPnt[iIDg] += fLocFitAdd ? 1:-1; // update counter | |
8a9ab0eb | 612 | // |
613 | } | |
614 | } // end of Update matrices | |
615 | // | |
616 | // calculate fMatCGlo -= fMatCGloLoc * fMatCLoc * fMatCGloLoc^T | |
617 | // and fVecBGlo -= fMatCGloLoc * fVecBLoc | |
618 | // | |
de34b538 | 619 | //-------------------------------------------------------------- >>> |
8a9ab0eb | 620 | double vll; |
621 | for (int iCIDg=0; iCIDg<nGloInFit; iCIDg++) { | |
622 | int iIDg = fCGlo2Glo[iCIDg]; | |
623 | // | |
624 | vl = 0; | |
de34b538 | 625 | Double_t *rowGLIDg = matCGloLoc(iCIDg); |
7c3070ec | 626 | for (int kl=0;kl<maxLocUsed;kl++) if (rowGLIDg[kl]) vl += rowGLIDg[kl]*fVecBLoc[kl]; |
cc9bec47 | 627 | if (!IsZero(vl)) fVecBGlo[iIDg] -= fLocFitAdd ? vl : -vl; |
8a9ab0eb | 628 | // |
de34b538 | 629 | int nfill = 0; |
8a9ab0eb | 630 | for (int jCIDg=0;jCIDg<=iCIDg; jCIDg++) { |
631 | int jIDg = fCGlo2Glo[jCIDg]; | |
632 | // | |
633 | vl = 0; | |
de34b538 | 634 | Double_t *rowGLJDg = matCGloLoc(jCIDg); |
7c3070ec | 635 | for (int kl=0;kl<maxLocUsed;kl++) { |
8a9ab0eb | 636 | // diag terms |
cc9bec47 | 637 | if ( (!IsZero(vll=rowGLIDg[kl]*rowGLJDg[kl])) ) vl += matCLoc.QueryDiag(kl)*vll; |
8a9ab0eb | 638 | // |
639 | // off-diag terms | |
640 | for (int ll=0;ll<kl;ll++) { | |
cc9bec47 | 641 | if ( !IsZero(vll=rowGLIDg[kl]*rowGLJDg[ll]) ) vl += matCLoc(kl,ll)*vll; |
642 | if ( !IsZero(vll=rowGLIDg[ll]*rowGLJDg[kl]) ) vl += matCLoc(kl,ll)*vll; | |
8a9ab0eb | 643 | } |
644 | } | |
cc9bec47 | 645 | if (!IsZero(vl)) { |
de34b538 | 646 | fFillIndex[nfill] = jIDg; |
647 | fFillValue[nfill++] = fLocFitAdd ? -vl : vl; | |
648 | } | |
8a9ab0eb | 649 | } |
de34b538 | 650 | if (nfill) matCGlo.AddToRow(iIDg,fFillValue,fFillIndex,nfill); |
8a9ab0eb | 651 | } |
652 | // | |
653 | // reset compressed index array | |
654 | // | |
de34b538 | 655 | for (int i=nGloInFit;i--;) { |
656 | fGlo2CGlo[ fCGlo2Glo[i] ] = -1; | |
657 | fCGlo2Glo[i] = -1; | |
658 | } | |
8a9ab0eb | 659 | // |
de34b538 | 660 | //---------------------------------------------------- <<< |
8a9ab0eb | 661 | return 1; |
662 | } | |
663 | ||
664 | //_____________________________________________________________________________________________ | |
665 | Int_t AliMillePede2::GlobalFit(Double_t *par, Double_t *error, Double_t *pull) | |
666 | { | |
667 | // performs a requested number of global iterations | |
668 | fIter = 1; | |
669 | // | |
670 | TStopwatch sw; sw.Start(); | |
671 | // | |
672 | int res = 0; | |
673 | AliInfo("Starting Global fit."); | |
674 | while (fIter<=fMaxIter) { | |
675 | // | |
676 | res = GlobalFitIteration(); | |
677 | if (!res) break; | |
678 | // | |
cc9bec47 | 679 | if (!IsZero(fChi2CutFactor-fChi2CutRef)) { |
8a9ab0eb | 680 | fChi2CutFactor = TMath::Sqrt(fChi2CutFactor); |
681 | if (fChi2CutFactor < 1.2*fChi2CutRef) { | |
682 | fChi2CutFactor = fChi2CutRef; | |
683 | fIter = fMaxIter - 1; // Last iteration | |
684 | } | |
685 | } | |
686 | fIter++; | |
687 | } | |
688 | // | |
689 | sw.Stop(); | |
690 | AliInfo(Form("Global fit %s, CPU time: %.1f",res ? "Converged":"Failed",sw.CpuTime())); | |
691 | if (!res) return 0; | |
692 | // | |
693 | if (par) for (int i=fNGloPar;i--;) par[i] = fInitPar[i]+fDeltaPar[i]; | |
694 | // | |
695 | if (fGloSolveStatus==gkInvert) { // errors on params are available | |
de34b538 | 696 | if (error) for (int i=fNGloPar;i--;) error[i] = fProcPnt[i]>0 ? TMath::Sqrt(TMath::Abs(fMatCGlo->QueryDiag(i))) : 0.; |
697 | if (pull) for (int i=fNGloPar;i--;) pull[i] = fProcPnt[i]>0 && (fSigmaPar[i]*fSigmaPar[i]-fMatCGlo->QueryDiag(i))>0. && fSigmaPar[i]>0 | |
698 | ? fDeltaPar[i]/TMath::Sqrt(fSigmaPar[i]*fSigmaPar[i]-fMatCGlo->QueryDiag(i)) : 0; | |
8a9ab0eb | 699 | } |
700 | // | |
701 | return 1; | |
702 | } | |
703 | ||
704 | //_____________________________________________________________________________________________ | |
705 | Int_t AliMillePede2::GlobalFitIteration() | |
706 | { | |
707 | // perform global parameters fit once all the local equations have been fitted | |
708 | // | |
709 | AliInfo(Form("Global Fit Iteration#%2d (Local Fit Chi^2 cut factor: %.2f)",fIter,fChi2CutFactor)); | |
710 | // | |
711 | if (!fNGloPar || !fTreeData) { | |
712 | AliInfo("No data was stored, aborting iteration"); | |
713 | return 0; | |
714 | } | |
de34b538 | 715 | TStopwatch sw,sws; |
716 | sw.Start(); | |
717 | sws.Stop(); | |
8a9ab0eb | 718 | // |
719 | if (!fConstrUsed) { | |
720 | fConstrUsed = new Bool_t[fNGloConstraints]; | |
721 | memset(fConstrUsed,0,fNGloConstraints*sizeof(Bool_t)); | |
722 | } | |
723 | // Reset all info specific for this step | |
de34b538 | 724 | AliMatrixSq& matCGlo = *fMatCGlo; |
725 | matCGlo.Reset(); | |
8a9ab0eb | 726 | memset(fProcPnt,0,fNGloPar*sizeof(Int_t)); |
727 | // | |
728 | fNGloConstraints = fTreeConstr ? fTreeConstr->GetEntries() : 0; | |
729 | // | |
de34b538 | 730 | // count number of Lagrange constraints: they need new row/cols to be added |
731 | fNLagrangeConstraints = 0; | |
732 | for (int i=0; i<fNGloConstraints; i++) { | |
733 | ReadRecordConstraint(i); | |
cc9bec47 | 734 | if ( IsZero(fRecord->GetValue(1)) ) fNLagrangeConstraints++; // exact constraint (no error) -> Lagrange multiplier |
de34b538 | 735 | } |
736 | // | |
8a9ab0eb | 737 | // if needed, readjust the size of the global vector (for matrices this is done automatically) |
de34b538 | 738 | if (!fVecBGlo || fNGloSize!=fNGloPar+fNLagrangeConstraints) { |
8a9ab0eb | 739 | delete[] fVecBGlo; // in case some constraint was added between the two manual iterations |
de34b538 | 740 | fNGloSize = fNGloPar+fNLagrangeConstraints; |
8a9ab0eb | 741 | fVecBGlo = new Double_t[fNGloSize]; |
742 | } | |
743 | memset(fVecBGlo,0,fNGloSize*sizeof(double)); | |
744 | // | |
745 | fNLocFits = 0; | |
746 | fNLocFitsRejected = 0; | |
de34b538 | 747 | fNLocEquations = 0; |
8a9ab0eb | 748 | // |
749 | // Process data records and build the matrices | |
750 | Long_t ndr = fTreeData->GetEntries(); | |
751 | AliInfo(Form("Building the Global matrix from %d data records",ndr)); | |
de34b538 | 752 | if (!ndr) return 0; |
753 | // | |
754 | TStopwatch swt; swt.Start(); | |
755 | fLocFitAdd = kTRUE; // add contributions of matching tracks | |
8a9ab0eb | 756 | for (Long_t i=0;i<ndr;i++) { |
757 | ReadRecordData(i); | |
758 | LocalFit(); | |
7c3070ec | 759 | if ( (i%int(0.2*ndr)) == 0) printf("%.1f%% of local fits done\n", double(100.*i)/ndr); |
8a9ab0eb | 760 | } |
de34b538 | 761 | swt.Stop(); |
762 | printf("%ld local fits done: ", ndr); | |
763 | swt.Print(); | |
764 | sw.Start(kFALSE); | |
765 | // | |
8a9ab0eb | 766 | // |
de34b538 | 767 | // ---------------------- Reject parameters with low statistics ------------>> |
8a9ab0eb | 768 | fNGloFix = 0; |
de34b538 | 769 | if (fMinPntValid>1 && fNGroupsSet) { |
770 | // | |
771 | printf("Checking parameters with statistics < %d\n",fMinPntValid); | |
772 | TStopwatch swsup; | |
773 | swsup.Start(); | |
774 | // 1) build the list of parameters to fix | |
775 | Int_t fixArrSize = 10; | |
776 | Int_t nFixedGroups = 0; | |
777 | TArrayI fixGroups(fixArrSize); | |
778 | // | |
91e1aa2c | 779 | int grIDold = -2; |
780 | int oldStart = -1; | |
781 | double oldMin = 1.e20; | |
782 | double oldMax =-1.e20; | |
783 | // | |
de34b538 | 784 | for (int i=fNGloPar;i--;) { // // Reset row and column of fixed params and add 1/sig^2 to free ones |
785 | int grID = fParamGrID[i]; | |
786 | if (grID<0) continue; // not in the group | |
91e1aa2c | 787 | // |
788 | if (grID!=grIDold) { // starting new group | |
789 | if (grIDold>=0) { // decide if the group has enough statistics | |
790 | if (oldMin<fMinPntValid && oldMax<2*fMinPntValid) { // suppress group | |
791 | for (int iold=oldStart;iold>i;iold--) fProcPnt[iold] = 0; | |
792 | Bool_t fnd = kFALSE; // check if the group is already accounted | |
793 | for (int j=nFixedGroups;j--;) if (fixGroups[j]==grIDold) {fnd=kTRUE; break;} | |
794 | if (!fnd) { | |
795 | if (nFixedGroups>=fixArrSize) {fixArrSize*=2; fixGroups.Set(fixArrSize);} | |
796 | fixGroups[nFixedGroups++] = grIDold; // add group to fix | |
797 | } | |
798 | } | |
799 | } | |
800 | grIDold = grID; // mark the start of the new group | |
801 | oldStart = i; | |
802 | oldMin = 1.e20; | |
803 | oldMax = -1.e20; | |
804 | } | |
805 | if (oldMin>fProcPnt[i]) oldMin = fProcPnt[i]; | |
806 | if (oldMax<fProcPnt[i]) oldMax = fProcPnt[i]; | |
de34b538 | 807 | // |
808 | } | |
91e1aa2c | 809 | // extra check for the last group |
810 | if (grIDold>=0 && oldMin<fMinPntValid && oldMax<2*fMinPntValid) { // suppress group | |
811 | for (int iold=oldStart;iold--;) fProcPnt[iold] = 0; | |
812 | Bool_t fnd = kFALSE; // check if the group is already accounted | |
813 | for (int j=nFixedGroups;j--;) if (fixGroups[j]==grIDold) {fnd=kTRUE; break;} | |
814 | if (!fnd) { | |
815 | if (nFixedGroups>=fixArrSize) {fixArrSize*=2; fixGroups.Set(fixArrSize);} | |
816 | fixGroups[nFixedGroups++] = grIDold; // add group to fix | |
817 | } | |
818 | } | |
819 | // | |
de34b538 | 820 | // 2) loop over records and add contributions of fixed groups with negative sign |
821 | fLocFitAdd = kFALSE; | |
822 | // | |
823 | for (Long_t i=0;i<ndr;i++) { | |
824 | ReadRecordData(i); | |
825 | Bool_t suppr = kFALSE; | |
826 | for (int ifx=nFixedGroups;ifx--;)if (fRecord->IsGroupPresent(fixGroups[ifx])) suppr = kTRUE; | |
827 | if (suppr) LocalFit(); | |
828 | } | |
829 | fLocFitAdd = kTRUE; | |
830 | // | |
831 | if (nFixedGroups) { | |
832 | printf("Suppressed contributions of groups with NPoints<%d :\n",fMinPntValid); | |
833 | for (int i=0;i<nFixedGroups;i++) printf("%d ",fixGroups[i]); printf("\n"); | |
834 | } | |
835 | swsup.Stop(); | |
836 | swsup.Print(); | |
837 | } | |
838 | // ---------------------- Reject parameters with low statistics ------------<< | |
839 | // | |
840 | // add large number to diagonal of fixed params | |
8a9ab0eb | 841 | // |
8a9ab0eb | 842 | for (int i=fNGloPar;i--;) { // // Reset row and column of fixed params and add 1/sig^2 to free ones |
7c3070ec | 843 | // printf("#%3d : Nproc : %5d grp: %d\n",i,fProcPnt[i],fParamGrID[i]); |
de34b538 | 844 | if (fProcPnt[i]<1) { |
8a9ab0eb | 845 | fNGloFix++; |
8a9ab0eb | 846 | fVecBGlo[i] = 0.; |
a8c5b94c | 847 | matCGlo.DiagElem(i) = 1.;//float(fNLocEquations*fNLocEquations); |
8a9ab0eb | 848 | } |
de34b538 | 849 | else matCGlo.DiagElem(i) += (fgWeightSigma ? fProcPnt[i] : 1.)/(fSigmaPar[i]*fSigmaPar[i]); |
8a9ab0eb | 850 | } |
851 | // | |
de34b538 | 852 | for (int i=fNGloPar;i--;) fDiagCGlo[i] = matCGlo.QueryDiag(i); // save the diagonal elements |
853 | // | |
8a9ab0eb | 854 | // add constraint equations |
855 | int nVar = fNGloPar; // Current size of global matrix | |
856 | for (int i=0; i<fNGloConstraints; i++) { | |
857 | ReadRecordConstraint(i); | |
858 | double val = fRecord->GetValue(0); | |
de34b538 | 859 | double sig = fRecord->GetValue(1); |
8a9ab0eb | 860 | int *indV = fRecord->GetIndex()+2; |
861 | double *der = fRecord->GetValue()+2; | |
862 | int csize = fRecord->GetSize()-2; | |
863 | // | |
de34b538 | 864 | // check if after suppression of fixed variables there are non-0 derivatives |
865 | // and determine the max statistics of involved params | |
866 | int nSuppressed = 0; | |
867 | int maxStat = 1; | |
868 | for (int j=csize;j--;) { | |
869 | if (fProcPnt[indV[j]]<1) nSuppressed++; | |
870 | else { | |
871 | maxStat = TMath::Max(maxStat,fProcPnt[indV[j]]); | |
872 | } | |
873 | } | |
8a9ab0eb | 874 | // |
de34b538 | 875 | if (nSuppressed==csize) { |
876 | // AliInfo(Form("Neglecting constraint %d of %d derivatives since no free parameters left",i,csize)); | |
8a9ab0eb | 877 | // |
878 | // was this constraint ever created ? | |
de34b538 | 879 | if ( sig==0 && fConstrUsed[i] ) { // this is needed only for constraints with Lagrange multiplier |
8a9ab0eb | 880 | // to avoid empty row impose dummy constraint on "Lagrange multiplier" |
de34b538 | 881 | matCGlo.DiagElem(nVar) = 1.; |
8a9ab0eb | 882 | fVecBGlo[nVar++] = 0; |
883 | } | |
884 | continue; | |
885 | } | |
886 | // | |
de34b538 | 887 | // account for already accumulated corrections |
888 | for (int j=csize; j--;) val -= der[j]*(fInitPar[ indV[j] ]+fDeltaPar[ indV[j] ]); | |
8a9ab0eb | 889 | // |
de34b538 | 890 | if (sig > 0) { // this is a gaussian constriant: no Lagrange multipliers are added |
891 | // | |
892 | double sig2i = (fgWeightSigma ? TMath::Sqrt(maxStat) : 1.)/sig/sig; | |
893 | for (int ir=0;ir<csize;ir++) { | |
894 | int iID = indV[ir]; | |
895 | for (int ic=0;ic<=ir;ic++) { // matrix is symmetric | |
896 | int jID = indV[ic]; | |
897 | double vl = der[ir]*der[ic]*sig2i; | |
cc9bec47 | 898 | if (!IsZero(vl)) matCGlo(iID,jID) += vl; |
de34b538 | 899 | } |
900 | fVecBGlo[iID] += val*der[ir]*sig2i; | |
901 | } | |
902 | } | |
903 | else { // this is exact constriant: Lagrange multipliers must be added | |
904 | for (int j=csize; j--;) { | |
905 | int jID = indV[j]; | |
906 | if (fProcPnt[jID]<1) continue; // this parameter was fixed, don't put it into constraint | |
907 | matCGlo(nVar,jID) = float(fNLocEquations)*der[j]; // fMatCGlo is symmetric, only lower triangle is filled | |
908 | } | |
909 | // | |
910 | if (matCGlo.QueryDiag(nVar)) matCGlo.DiagElem(nVar) = 0.0; | |
911 | fVecBGlo[nVar++] = float(fNLocEquations)*val; //RS ? should we use here fNLocFits ? | |
912 | fConstrUsed[i] = kTRUE; | |
913 | } | |
8a9ab0eb | 914 | } |
915 | // | |
916 | AliInfo(Form("Obtained %-7ld equations from %-7ld records (%-7ld rejected). Fixed %-4d globals", | |
917 | fNLocEquations,fNLocFits,fNLocFitsRejected,fNGloFix)); | |
918 | ||
919 | // | |
de34b538 | 920 | sws.Start(); |
8a9ab0eb | 921 | fGloSolveStatus = SolveGlobalMatEq(); // obtain solution for this step |
de34b538 | 922 | sws.Stop(); |
923 | printf("Solve %d |",fIter); sws.Print(); | |
8a9ab0eb | 924 | // |
925 | sw.Stop(); | |
926 | AliInfo(Form("Iteration#%2d %s. CPU time: %.1f",fIter,fGloSolveStatus==gkFailed ? "Failed":"Converged",sw.CpuTime())); | |
927 | if (fGloSolveStatus==gkFailed) return 0; | |
928 | // | |
929 | for (int i=fNGloPar;i--;) fDeltaPar[i] += fVecBGlo[i]; // Update global parameters values (for iterations) | |
930 | // | |
de34b538 | 931 | // PrintGlobalParameters(); |
8a9ab0eb | 932 | return 1; |
933 | } | |
934 | ||
935 | //_____________________________________________________________________________________________ | |
936 | Int_t AliMillePede2::SolveGlobalMatEq() | |
937 | { | |
938 | // | |
939 | // solve global matrix equation MatCGlob*X=VecBGlo and store the result in the VecBGlo | |
940 | // | |
de34b538 | 941 | /* |
942 | printf("GlobalMatrix\n"); | |
943 | fMatCGlo->Print(); | |
944 | printf("RHS\n"); | |
945 | for (int i=0;i<fNGloPar;i++) printf("%d %+e\n",i,fVecBGlo[i]); | |
946 | */ | |
947 | // | |
8a9ab0eb | 948 | if (!fgIsMatGloSparse) { |
949 | // | |
de34b538 | 950 | if (fNLagrangeConstraints==0) { // pos-def systems are faster to solve by Cholesky |
951 | if ( ((AliSymMatrix*)fMatCGlo)->SolveChol(fVecBGlo, fgInvChol) ) return fgInvChol ? gkInvert:gkNoInversion; | |
8a9ab0eb | 952 | else AliInfo("Solution of Global Dense System by Cholesky failed, trying Gaussian Elimiation"); |
953 | } | |
954 | // | |
955 | if (((AliSymMatrix*)fMatCGlo)->SolveSpmInv(fVecBGlo, kTRUE)) return gkInvert; | |
956 | else AliInfo("Solution of Global Dense System by Gaussian Elimiation failed, trying iterative methods"); | |
957 | } | |
958 | // try to solve by minres | |
de34b538 | 959 | TVectorD sol(fNGloSize); |
8a9ab0eb | 960 | // |
961 | AliMinResSolve *slv = new AliMinResSolve(fMatCGlo,fVecBGlo); | |
962 | if (!slv) return gkFailed; | |
963 | // | |
964 | Bool_t res = kFALSE; | |
965 | if (fgIterSol == AliMinResSolve::kSolMinRes) | |
de34b538 | 966 | res = slv->SolveMinRes(sol,fgMinResCondType,fgMinResMaxIter,fgMinResTol); |
8a9ab0eb | 967 | else if (fgIterSol == AliMinResSolve::kSolFGMRes) |
de34b538 | 968 | res = slv->SolveFGMRES(sol,fgMinResCondType,fgMinResMaxIter,fgMinResTol,fgNKrylovV); |
8a9ab0eb | 969 | else |
970 | AliInfo(Form("Undefined Iteritive Solver ID=%d, only %d are defined",fgIterSol,AliMinResSolve::kNSolvers)); | |
971 | // | |
972 | if (!res) return gkFailed; | |
de34b538 | 973 | for (int i=fNGloSize;i--;) fVecBGlo[i] = sol[i]; |
974 | return gkNoInversion; | |
8a9ab0eb | 975 | // |
976 | } | |
977 | ||
978 | //_____________________________________________________________________________________________ | |
979 | Float_t AliMillePede2::Chi2DoFLim(int nSig, int nDoF) const | |
980 | { | |
981 | /// return the limit in chi^2/nd for n sigmas stdev authorized | |
982 | // Only n=1, 2, and 3 are expected in input | |
983 | int lNSig; | |
984 | float sn[3] = {0.47523, 1.690140, 2.782170}; | |
985 | float table[3][30] = {{1.0000, 1.1479, 1.1753, 1.1798, 1.1775, 1.1730, 1.1680, 1.1630, | |
986 | 1.1581, 1.1536, 1.1493, 1.1454, 1.1417, 1.1383, 1.1351, 1.1321, | |
987 | 1.1293, 1.1266, 1.1242, 1.1218, 1.1196, 1.1175, 1.1155, 1.1136, | |
988 | 1.1119, 1.1101, 1.1085, 1.1070, 1.1055, 1.1040}, | |
989 | {4.0000, 3.0900, 2.6750, 2.4290, 2.2628, 2.1415, 2.0481, 1.9736, | |
990 | 1.9124, 1.8610, 1.8171, 1.7791, 1.7457, 1.7161, 1.6897, 1.6658, | |
991 | 1.6442, 1.6246, 1.6065, 1.5899, 1.5745, 1.5603, 1.5470, 1.5346, | |
992 | 1.5230, 1.5120, 1.5017, 1.4920, 1.4829, 1.4742}, | |
993 | {9.0000, 5.9146, 4.7184, 4.0628, 3.6410, 3.3436, 3.1209, 2.9468, | |
994 | 2.8063, 2.6902, 2.5922, 2.5082, 2.4352, 2.3711, 2.3143, 2.2635, | |
995 | 2.2178, 2.1764, 2.1386, 2.1040, 2.0722, 2.0428, 2.0155, 1.9901, | |
996 | 1.9665, 1.9443, 1.9235, 1.9040, 1.8855, 1.8681}}; | |
997 | ||
998 | if (nDoF < 1) { | |
999 | return 0.0; | |
1000 | } | |
1001 | else { | |
1002 | lNSig = TMath::Max(1,TMath::Min(nSig,3)); | |
1003 | ||
1004 | if (nDoF <= 30) { | |
1005 | return table[lNSig-1][nDoF-1]; | |
1006 | } | |
1007 | else { // approximation | |
1008 | return ((sn[lNSig-1]+TMath::Sqrt(float(2*nDoF-3)))* | |
1009 | (sn[lNSig-1]+TMath::Sqrt(float(2*nDoF-3))))/float(2*nDoF-2); | |
1010 | } | |
1011 | } | |
1012 | } | |
1013 | ||
1014 | //_____________________________________________________________________________________________ | |
1015 | Int_t AliMillePede2::SetIterations(double lChi2CutFac) | |
1016 | { | |
1017 | // Number of iterations is calculated from lChi2CutFac | |
1018 | fChi2CutFactor = TMath::Max(1.0, lChi2CutFac); | |
1019 | // | |
1020 | AliInfo(Form("Initial cut factor is %f",fChi2CutFactor)); | |
1021 | fIter = 1; // Initializes the iteration process | |
1022 | return 1; | |
1023 | } | |
1024 | ||
1025 | //_____________________________________________________________________________________________ | |
1026 | Double_t AliMillePede2::GetParError(int iPar) const | |
1027 | { | |
1028 | // return error for parameter iPar | |
1029 | if (fGloSolveStatus==gkInvert) { | |
de34b538 | 1030 | double res = fMatCGlo->QueryDiag(iPar); |
8a9ab0eb | 1031 | if (res>=0) return TMath::Sqrt(res); |
1032 | } | |
1033 | return -1.; | |
1034 | } | |
1035 | ||
1036 | ||
1037 | //_____________________________________________________________________________________________ | |
1038 | Int_t AliMillePede2::PrintGlobalParameters() const | |
1039 | { | |
1040 | /// Print the final results into the logfile | |
1041 | double lError = 0.; | |
1042 | double lGlobalCor =0.; | |
1043 | ||
1044 | AliInfo(""); | |
1045 | AliInfo(" Result of fit for global parameters"); | |
1046 | AliInfo(" ==================================="); | |
1047 | AliInfo(" I initial final differ lastcor error gcor Npnt"); | |
1048 | AliInfo("----------------------------------------------------------------------------------------------"); | |
1049 | // | |
1050 | for (int i=0; i<fNGloPar; i++) { | |
1051 | lError = GetParError(i); | |
1052 | lGlobalCor = 0.0; | |
1053 | // | |
1054 | double dg; | |
de34b538 | 1055 | if (fGloSolveStatus==gkInvert && TMath::Abs( (dg=fMatCGlo->QueryDiag(i)) *fDiagCGlo[i]) > 0) { |
8a9ab0eb | 1056 | lGlobalCor = TMath::Sqrt(TMath::Abs(1.0-1.0/(dg*fDiagCGlo[i]))); |
1057 | AliInfo(Form("%d\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f\t %6d", | |
1058 | i,fInitPar[i],fInitPar[i]+fDeltaPar[i],fDeltaPar[i],fVecBGlo[i],lError,lGlobalCor,fProcPnt[i])); | |
1059 | } | |
1060 | else { | |
1061 | AliInfo(Form("%d\t %.6f\t %.6f\t %.6f\t %.6f\t OFF\t OFF\t %6d",i,fInitPar[i],fInitPar[i]+fDeltaPar[i], | |
1062 | fDeltaPar[i],fVecBGlo[i],fProcPnt[i])); | |
1063 | } | |
1064 | } | |
1065 | return 1; | |
1066 | } |