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