]>
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 */ | |
339fbe23 | 4 | /* and AliMillepede class by Javier */ |
5 | /* Allows operations with large sparse matrices */ | |
de34b538 | 6 | /* http://www.desy.de/~blobel/mptalks.html */ |
7 | /* */ | |
8 | /* Author: ruben.shahoyan@cern.ch */ | |
9 | /* */ | |
10 | /**********************************************************************************************/ | |
11 | ||
8a9ab0eb | 12 | #include "AliMillePede2.h" |
13 | #include "AliLog.h" | |
14 | #include <TStopwatch.h> | |
de34b538 | 15 | #include <TFile.h> |
f9efbbf6 | 16 | #include <TChain.h> |
de34b538 | 17 | #include <TMath.h> |
7c3070ec | 18 | #include <TVectorD.h> |
e30a812f | 19 | #include <TArrayL.h> |
1f6eade8 | 20 | #include <TArrayF.h> |
f9efbbf6 | 21 | #include <TSystem.h> |
de34b538 | 22 | #include "AliMatrixSq.h" |
23 | #include "AliSymMatrix.h" | |
24 | #include "AliRectMatrix.h" | |
25 | #include "AliMatrixSparse.h" | |
e30a812f | 26 | #include <stdio.h> |
27 | #include <stdlib.h> | |
28 | #include <unistd.h> | |
29 | #include <sys/types.h> | |
30 | #include <sys/stat.h> | |
31 | #include <fcntl.h> | |
f9efbbf6 | 32 | #include <fstream> |
8a9ab0eb | 33 | |
7fb585bc | 34 | //#define _DUMP_EQ_BEFORE_ |
35 | //#define _DUMP_EQ_AFTER_ | |
e041af50 | 36 | |
c9f619c7 | 37 | //#define _DUMPEQ_BEFORE_ |
38 | //#define _DUMPEQ_AFTER_ | |
8a9ab0eb | 39 | |
a0ea48e0 | 40 | using std::ifstream; |
8a9ab0eb | 41 | ClassImp(AliMillePede2) |
42 | ||
de34b538 | 43 | Bool_t AliMillePede2::fgInvChol = kTRUE; // Invert global matrix with Cholesky solver |
44 | Bool_t AliMillePede2::fgWeightSigma = kTRUE; // weight local constraint by module statistics | |
45 | Bool_t AliMillePede2::fgIsMatGloSparse = kFALSE; // use faster dense matrix by default | |
46 | Int_t AliMillePede2::fgMinResCondType = 1; // Jacoby preconditioner by default | |
47 | Double_t AliMillePede2::fgMinResTol = 1.e-11; // default tolerance | |
48 | Int_t AliMillePede2::fgMinResMaxIter = 10000; // default max number of iterations | |
8a9ab0eb | 49 | Int_t AliMillePede2::fgIterSol = AliMinResSolve::kSolMinRes; // default iterative solver |
de34b538 | 50 | Int_t AliMillePede2::fgNKrylovV = 240; // default number of Krylov vectors to keep |
8a9ab0eb | 51 | |
52 | //_____________________________________________________________________________________________ | |
53 | AliMillePede2::AliMillePede2() | |
54 | : fNLocPar(0), | |
55 | fNGloPar(0), | |
9a919f12 | 56 | fNGloParIni(0), |
8a9ab0eb | 57 | fNGloSize(0), |
58 | // | |
59 | fNLocEquations(0), | |
60 | fIter(0), | |
61 | fMaxIter(10), | |
62 | fNStdDev(3), | |
63 | fNGloConstraints(0), | |
de34b538 | 64 | fNLagrangeConstraints(0), |
8a9ab0eb | 65 | fNLocFits(0), |
66 | fNLocFitsRejected(0), | |
67 | fNGloFix(0), | |
339fbe23 | 68 | fGloSolveStatus(kFailed), |
8a9ab0eb | 69 | // |
70 | fChi2CutFactor(1.), | |
71 | fChi2CutRef(1.), | |
72 | fResCutInit(100.), | |
73 | fResCut(100.), | |
de34b538 | 74 | fMinPntValid(1), |
8a9ab0eb | 75 | // |
de34b538 | 76 | fNGroupsSet(0), |
77 | fParamGrID(0), | |
8a9ab0eb | 78 | fProcPnt(0), |
79 | fVecBLoc(0), | |
80 | fDiagCGlo(0), | |
81 | fVecBGlo(0), | |
82 | fInitPar(0), | |
83 | fDeltaPar(0), | |
84 | fSigmaPar(0), | |
85 | fIsLinear(0), | |
86 | fConstrUsed(0), | |
87 | // | |
88 | fGlo2CGlo(0), | |
89 | fCGlo2Glo(0), | |
90 | // | |
91 | fMatCLoc(0), | |
92 | fMatCGlo(0), | |
93 | fMatCGloLoc(0), | |
de34b538 | 94 | // |
95 | fFillIndex(0), | |
96 | fFillValue(0), | |
97 | // | |
58e8ba3a | 98 | fRecDataTreeName("AliMillePedeRecords_Data"), |
99 | fRecConsTreeName("AliMillePedeRecords_Consaints"), | |
100 | fRecDataBranchName("Record_Data"), | |
101 | fRecConsBranchName("Record_Consaints"), | |
102 | // | |
8a9ab0eb | 103 | fDataRecFName("/tmp/mp2_data_records.root"), |
104 | fRecord(0), | |
105 | fDataRecFile(0), | |
106 | fTreeData(0), | |
a8c5b94c | 107 | fRecFileStatus(0), |
8a9ab0eb | 108 | // |
109 | fConstrRecFName("/tmp/mp2_constraints_records.root"), | |
110 | fTreeConstr(0), | |
111 | fConsRecFile(0), | |
112 | fCurrRecDataID(0), | |
de34b538 | 113 | fCurrRecConstrID(0), |
e30a812f | 114 | fLocFitAdd(kTRUE), |
c9f619c7 | 115 | fUseRecordWeight(kTRUE), |
116 | fMinRecordLength(1), | |
e30a812f | 117 | fSelFirst(1), |
118 | fSelLast(-1), | |
119 | fRejRunList(0), | |
9a919f12 | 120 | fAccRunList(0), |
1f6eade8 | 121 | fAccRunListWgh(0), |
122 | fRunWgh(1), | |
9a919f12 | 123 | fkReGroup(0) |
e041af50 | 124 | { |
125 | fWghScl[0] = fWghScl[1] = -1; | |
126 | } | |
8a9ab0eb | 127 | |
128 | //_____________________________________________________________________________________________ | |
129 | AliMillePede2::AliMillePede2(const AliMillePede2& src) : | |
9a919f12 | 130 | TObject(src),fNLocPar(0),fNGloPar(0),fNGloParIni(0),fNGloSize(0),fNLocEquations(0),fIter(0), |
de34b538 | 131 | fMaxIter(10),fNStdDev(3),fNGloConstraints(0),fNLagrangeConstraints(0), |
132 | fNLocFits(0),fNLocFitsRejected(0), | |
8a9ab0eb | 133 | fNGloFix(0),fGloSolveStatus(0),fChi2CutFactor(0),fChi2CutRef(0),fResCutInit(0), |
de34b538 | 134 | fResCut(0),fMinPntValid(1),fNGroupsSet(0),fParamGrID(0),fProcPnt(0),fVecBLoc(0),fDiagCGlo(0),fVecBGlo(0), |
8a9ab0eb | 135 | fInitPar(0),fDeltaPar(0),fSigmaPar(0),fIsLinear(0),fConstrUsed(0),fGlo2CGlo(0),fCGlo2Glo(0), |
de34b538 | 136 | fMatCLoc(0),fMatCGlo(0),fMatCGloLoc(0),fFillIndex(0),fFillValue(0), |
58e8ba3a | 137 | fRecDataTreeName(0),fRecConsTreeName(0),fRecDataBranchName(0),fRecConsBranchName(0), |
de34b538 | 138 | fDataRecFName(0),fRecord(0),fDataRecFile(0), |
a8c5b94c | 139 | fTreeData(0),fRecFileStatus(0),fConstrRecFName(0),fTreeConstr(0),fConsRecFile(0),fCurrRecDataID(0), |
e30a812f | 140 | fCurrRecConstrID(0),fLocFitAdd(kTRUE), |
c9f619c7 | 141 | fUseRecordWeight(kTRUE), |
142 | fMinRecordLength(1), | |
e30a812f | 143 | fSelFirst(1), |
144 | fSelLast(-1), | |
145 | fRejRunList(0), | |
9a919f12 | 146 | fAccRunList(0), |
1f6eade8 | 147 | fAccRunListWgh(0), |
148 | fRunWgh(1), | |
9a919f12 | 149 | fkReGroup(0) |
e041af50 | 150 | { |
151 | fWghScl[0] = src.fWghScl[0]; | |
152 | fWghScl[1] = src.fWghScl[1]; | |
153 | printf("Dummy\n"); | |
154 | } | |
8a9ab0eb | 155 | |
156 | //_____________________________________________________________________________________________ | |
157 | AliMillePede2::~AliMillePede2() | |
158 | { | |
339fbe23 | 159 | // destructor |
8a9ab0eb | 160 | CloseDataRecStorage(); |
161 | CloseConsRecStorage(); | |
162 | // | |
de34b538 | 163 | delete[] fParamGrID; |
8a9ab0eb | 164 | delete[] fProcPnt; |
165 | delete[] fVecBLoc; | |
166 | delete[] fDiagCGlo; | |
167 | delete[] fVecBGlo; | |
168 | delete[] fInitPar; | |
169 | delete[] fDeltaPar; | |
170 | delete[] fSigmaPar; | |
171 | delete[] fGlo2CGlo; | |
172 | delete[] fCGlo2Glo; | |
173 | delete[] fIsLinear; | |
174 | delete[] fConstrUsed; | |
de34b538 | 175 | delete[] fFillIndex; |
176 | delete[] fFillValue; | |
8a9ab0eb | 177 | // |
178 | delete fRecord; | |
179 | delete fMatCLoc; | |
180 | delete fMatCGlo; | |
181 | delete fMatCGloLoc; | |
e30a812f | 182 | delete fRejRunList; |
183 | delete fAccRunList; | |
1f6eade8 | 184 | delete fAccRunListWgh; |
8a9ab0eb | 185 | } |
186 | ||
187 | //_____________________________________________________________________________________________ | |
9a919f12 | 188 | Int_t AliMillePede2::InitMille(int nGlo, int nLoc, int lNStdDev,double lResCut, double lResCutInit, const Int_t* regroup) |
8a9ab0eb | 189 | { |
339fbe23 | 190 | // init all |
9a919f12 | 191 | // |
192 | fNGloParIni = nGlo; | |
193 | if (regroup) { // regrouping is requested | |
194 | fkReGroup = regroup; | |
195 | int ng = 0; // recalculate N globals | |
196 | int maxPID = -1; | |
197 | for (int i=0;i<nGlo;i++) if (regroup[i]>=0) {ng++; if (regroup[i]>maxPID) maxPID = regroup[i];} | |
f91a16b5 | 198 | maxPID++; |
199 | AliInfo(Form("Regrouping is requested: from %d raw to %d formal globals grouped to %d real globals",nGlo,ng,maxPID)); | |
200 | nGlo = maxPID; | |
9a919f12 | 201 | } |
8a9ab0eb | 202 | if (nLoc>0) fNLocPar = nLoc; |
203 | if (nGlo>0) fNGloPar = nGlo; | |
204 | if (lResCutInit>0) fResCutInit = lResCutInit; | |
205 | if (lResCut>0) fResCut = lResCut; | |
206 | if (lNStdDev>0) fNStdDev = lNStdDev; | |
207 | // | |
d56a0575 | 208 | AliInfo(Form("NLoc: %d NGlo: %d",fNLocPar,fNGloPar)); |
209 | ||
8a9ab0eb | 210 | fNGloSize = fNGloPar; |
211 | // | |
3c5b4cc8 | 212 | if (fgIsMatGloSparse) {fMatCGlo = new AliMatrixSparse(fNGloPar); fMatCGlo->SetSymmetric(kTRUE);} |
213 | else fMatCGlo = new AliSymMatrix(fNGloPar); | |
214 | // | |
215 | fFillIndex = new Int_t[fNGloPar]; | |
216 | fFillValue = new Double_t[fNGloPar]; | |
217 | // | |
218 | fMatCLoc = new AliSymMatrix(fNLocPar); | |
219 | fMatCGloLoc = new AliRectMatrix(fNGloPar,fNLocPar); | |
220 | // | |
221 | fParamGrID = new Int_t[fNGloPar]; | |
222 | fProcPnt = new Int_t[fNGloPar]; | |
223 | fVecBLoc = new Double_t[fNLocPar]; | |
224 | fDiagCGlo = new Double_t[fNGloPar]; | |
225 | // | |
226 | fInitPar = new Double_t[fNGloPar]; | |
227 | fDeltaPar = new Double_t[fNGloPar]; | |
228 | fSigmaPar = new Double_t[fNGloPar]; | |
229 | fIsLinear = new Bool_t[fNGloPar]; | |
230 | // | |
231 | fGlo2CGlo = new Int_t[fNGloPar]; | |
232 | fCGlo2Glo = new Int_t[fNGloPar]; | |
8a9ab0eb | 233 | // |
234 | memset(fVecBLoc ,0,fNLocPar*sizeof(Double_t)); | |
235 | memset(fDiagCGlo ,0,fNGloPar*sizeof(Double_t)); | |
236 | memset(fInitPar ,0,fNGloPar*sizeof(Double_t)); | |
237 | memset(fDeltaPar ,0,fNGloPar*sizeof(Double_t)); | |
238 | memset(fSigmaPar ,0,fNGloPar*sizeof(Double_t)); | |
239 | memset(fProcPnt ,0,fNGloPar*sizeof(Int_t)); | |
240 | // | |
241 | for (int i=fNGloPar;i--;) { | |
de34b538 | 242 | fGlo2CGlo[i] = fCGlo2Glo[i] = -1; |
8a9ab0eb | 243 | fIsLinear[i] = kTRUE; |
de34b538 | 244 | fParamGrID[i] = -1; |
8a9ab0eb | 245 | } |
246 | // | |
c9f619c7 | 247 | fWghScl[0] = -1; |
248 | fWghScl[1] = -1; | |
8a9ab0eb | 249 | return 1; |
250 | } | |
251 | ||
8a9ab0eb | 252 | //_____________________________________________________________________________________________ |
253 | Bool_t AliMillePede2::ImposeDataRecFile(const char* fname) | |
254 | { | |
339fbe23 | 255 | // set filename for records |
8a9ab0eb | 256 | CloseDataRecStorage(); |
257 | SetDataRecFName(fname); | |
258 | return InitDataRecStorage(kTRUE); // open in read mode | |
259 | } | |
260 | ||
261 | //_____________________________________________________________________________________________ | |
262 | Bool_t AliMillePede2::ImposeConsRecFile(const char* fname) | |
263 | { | |
339fbe23 | 264 | // set filename for constraints |
8a9ab0eb | 265 | CloseConsRecStorage(); |
266 | SetConsRecFName(fname); | |
267 | return InitConsRecStorage(kTRUE); // open in read mode | |
268 | } | |
269 | ||
270 | //_____________________________________________________________________________________________ | |
271 | Bool_t AliMillePede2::InitDataRecStorage(Bool_t read) | |
272 | { | |
273 | // initialize the buffer for processed measurements records | |
274 | // | |
f9efbbf6 | 275 | if (fTreeData) {AliInfo("Data Records File is already initialized"); return kFALSE;} |
8a9ab0eb | 276 | // |
277 | if (!fRecord) fRecord = new AliMillePedeRecord(); | |
278 | // | |
f9efbbf6 | 279 | if (!read) { // write mode: cannot use chain |
280 | fDataRecFile = TFile::Open(GetDataRecFName(),"recreate"); | |
281 | if (!fDataRecFile) {AliFatal(Form("Failed to initialize data records file %s",GetDataRecFName())); return kFALSE;} | |
282 | AliInfo(Form("File %s used for derivatives records",GetDataRecFName())); | |
58e8ba3a | 283 | fTreeData = new TTree(GetRecDataTreeName(),"Data Records for AliMillePede2"); |
284 | fTreeData->Branch(GetRecDataBranchName(),"AliMillePedeRecord",&fRecord,32000,99); | |
8a9ab0eb | 285 | } |
f9efbbf6 | 286 | else { // use chain |
58e8ba3a | 287 | TChain* ch = new TChain(GetRecDataTreeName()); |
f9efbbf6 | 288 | // |
289 | if (fDataRecFName.EndsWith(".root")) ch->AddFile(fDataRecFName); | |
290 | else { // assume text file with list of filenames | |
291 | // | |
292 | ifstream inpf(fDataRecFName.Data()); | |
293 | if (!inpf.good()) {AliInfo(Form("Failed on input records list %s\n",fDataRecFName.Data())); return kFALSE;} | |
294 | // | |
295 | TString recfName; | |
f9efbbf6 | 296 | while ( !(recfName.ReadLine(inpf)).eof() ) { |
297 | recfName = recfName.Strip(TString::kBoth,' '); | |
298 | if (recfName.BeginsWith("//") || recfName.BeginsWith("#") || !recfName.EndsWith(".root")) { // comment | |
299 | AliInfo(Form("Skip %s\n",recfName.Data())); | |
300 | continue; | |
301 | } | |
302 | // | |
303 | recfName = recfName.Strip(TString::kBoth,','); | |
304 | recfName = recfName.Strip(TString::kBoth,'"'); | |
305 | gSystem->ExpandPathName(recfName); | |
306 | printf("Adding %s\n",recfName.Data()); | |
307 | ch->AddFile(recfName.Data()); | |
308 | } | |
309 | } | |
310 | // | |
311 | Long64_t nent = ch->GetEntries(); | |
312 | if (nent<1) { AliInfo("Obtained chain is empty"); return kFALSE;} | |
313 | fTreeData = ch; | |
58e8ba3a | 314 | fTreeData->SetBranchAddress(GetRecDataBranchName(),&fRecord); |
f9efbbf6 | 315 | AliInfo(Form("Found %lld derivatives records",nent)); |
316 | } | |
8a9ab0eb | 317 | fCurrRecDataID = -1; |
a8c5b94c | 318 | fRecFileStatus = read ? 1:2; |
8a9ab0eb | 319 | // |
320 | return kTRUE; | |
321 | } | |
322 | ||
323 | //_____________________________________________________________________________________________ | |
324 | Bool_t AliMillePede2::InitConsRecStorage(Bool_t read) | |
325 | { | |
326 | // initialize the buffer for processed measurements records | |
327 | // | |
328 | if (fConsRecFile) {AliInfo("Constraints Records File is already initialized"); return kFALSE;} | |
329 | // | |
330 | if (!fRecord) fRecord = new AliMillePedeRecord(); | |
331 | // | |
332 | fConsRecFile = TFile::Open(GetConsRecFName(),read ? "":"recreate"); | |
333 | if (!fConsRecFile) {AliInfo(Form("Failed to initialize constraints records file %s",GetConsRecFName())); return kFALSE;} | |
334 | // | |
335 | AliInfo(Form("File %s used for constraints records",GetConsRecFName())); | |
336 | if (read) { | |
58e8ba3a | 337 | fTreeConstr = (TTree*)fConsRecFile->Get(GetRecConsTreeName()); |
8a9ab0eb | 338 | if (!fTreeConstr) {AliInfo(Form("Did not find constraints records tree in %s",GetConsRecFName())); return kFALSE;} |
58e8ba3a | 339 | fTreeConstr->SetBranchAddress(GetRecConsBranchName(),&fRecord); |
73bbf779 | 340 | AliInfo(Form("Found %lld constraints records",fTreeConstr->GetEntries())); |
8a9ab0eb | 341 | // |
342 | } | |
343 | else { | |
344 | // | |
58e8ba3a | 345 | fTreeConstr = new TTree(GetRecConsTreeName(),"Constraints Records for AliMillePede2"); |
346 | fTreeConstr->Branch(GetRecConsBranchName(),"AliMillePedeRecord",&fRecord,32000,99); | |
8a9ab0eb | 347 | } |
348 | fCurrRecConstrID = -1; | |
349 | // | |
350 | return kTRUE; | |
351 | } | |
352 | ||
353 | //_____________________________________________________________________________________________ | |
354 | void AliMillePede2::CloseDataRecStorage() | |
355 | { | |
339fbe23 | 356 | // close records file |
8a9ab0eb | 357 | if (fTreeData) { |
f9efbbf6 | 358 | if (fDataRecFile && fDataRecFile->IsWritable()) { |
8a9ab0eb | 359 | fDataRecFile->cd(); |
360 | fTreeData->Write(); | |
361 | } | |
362 | delete fTreeData; | |
363 | fTreeData = 0; | |
f9efbbf6 | 364 | if (fDataRecFile) { |
365 | fDataRecFile->Close(); | |
366 | delete fDataRecFile; | |
367 | fDataRecFile = 0; | |
368 | } | |
8a9ab0eb | 369 | } |
a8c5b94c | 370 | fRecFileStatus = 0; |
8a9ab0eb | 371 | // |
372 | } | |
373 | ||
374 | //_____________________________________________________________________________________________ | |
375 | void AliMillePede2::CloseConsRecStorage() | |
376 | { | |
339fbe23 | 377 | // close constraints file |
8a9ab0eb | 378 | if (fTreeConstr) { |
379 | if (fConsRecFile->IsWritable()) { | |
380 | fConsRecFile->cd(); | |
381 | fTreeConstr->Write(); | |
382 | } | |
383 | delete fTreeConstr; | |
384 | fTreeConstr = 0; | |
385 | fConsRecFile->Close(); | |
386 | delete fConsRecFile; | |
387 | fConsRecFile = 0; | |
388 | } | |
389 | // | |
390 | } | |
391 | ||
392 | //_____________________________________________________________________________________________ | |
393 | Bool_t AliMillePede2::ReadNextRecordData() | |
394 | { | |
395 | // read next data record (if any) | |
396 | if (!fTreeData || ++fCurrRecDataID >= fTreeData->GetEntries()) { fCurrRecDataID--; return kFALSE;} | |
397 | fTreeData->GetEntry(fCurrRecDataID); | |
398 | return kTRUE; | |
399 | } | |
400 | ||
401 | //_____________________________________________________________________________________________ | |
402 | Bool_t AliMillePede2::ReadNextRecordConstraint() | |
403 | { | |
404 | // read next constraint record (if any) | |
405 | if (!fTreeConstr || ++fCurrRecConstrID >= fTreeConstr->GetEntries()) { fCurrRecConstrID--; return kFALSE;} | |
406 | fTreeConstr->GetEntry(fCurrRecConstrID); | |
407 | return kTRUE; | |
408 | } | |
409 | ||
a8c5b94c | 410 | //_____________________________________________________________________________________________ |
411 | void AliMillePede2::SetRecordWeight(double wgh) | |
412 | { | |
339fbe23 | 413 | // assign weight |
a8c5b94c | 414 | if (fRecFileStatus<2) InitDataRecStorage(); // create a buffer to store the data |
415 | fRecord->SetWeight(wgh); | |
416 | } | |
417 | ||
e30a812f | 418 | //_____________________________________________________________________________________________ |
419 | void AliMillePede2::SetRecordRun(Int_t run) | |
420 | { | |
339fbe23 | 421 | // assign run |
e30a812f | 422 | if (fRecFileStatus<2) InitDataRecStorage(); // create a buffer to store the data |
423 | fRecord->SetRunID(run); | |
424 | } | |
425 | ||
8a9ab0eb | 426 | //_____________________________________________________________________________________________ |
427 | void AliMillePede2::SetLocalEquation(double *dergb, double *derlc, double lMeas, double lSigma) | |
428 | { | |
339fbe23 | 429 | // assing derivs of loc.eq. |
a8c5b94c | 430 | if (fRecFileStatus<2) InitDataRecStorage(); // create a buffer to store the data |
8a9ab0eb | 431 | // |
432 | // write data of single measurement | |
433 | if (lSigma<=0.0) { // If parameter is fixed, then no equation | |
434 | for (int i=fNLocPar; i--;) derlc[i] = 0.0; | |
9a919f12 | 435 | for (int i=fNGloParIni; i--;) dergb[i] = 0.0; |
8a9ab0eb | 436 | return; |
437 | } | |
438 | // | |
439 | fRecord->AddResidual(lMeas); | |
440 | // | |
441 | // Retrieve local param interesting indices | |
cc9bec47 | 442 | for (int i=0;i<fNLocPar;i++) if (!IsZero(derlc[i])) {fRecord->AddIndexValue(i,derlc[i]); derlc[i] = 0.0;} |
8a9ab0eb | 443 | // |
444 | fRecord->AddWeight( 1.0/lSigma/lSigma ); | |
445 | // | |
446 | // Idem for global parameters | |
9a919f12 | 447 | for (int i=0;i<fNGloParIni;i++) if (!IsZero(dergb[i])) { |
de34b538 | 448 | fRecord->AddIndexValue(i,dergb[i]); dergb[i] = 0.0; |
9a919f12 | 449 | int idrg = GetRGId(i); |
450 | fRecord->MarkGroup(idrg<0 ? -1 : fParamGrID[i]); | |
de34b538 | 451 | } |
d56a0575 | 452 | // fRecord->Print(); |
8a9ab0eb | 453 | // |
454 | } | |
455 | ||
456 | //_____________________________________________________________________________________________ | |
457 | void AliMillePede2::SetLocalEquation(int *indgb, double *dergb, int ngb, int *indlc, | |
458 | double *derlc,int nlc,double lMeas,double lSigma) | |
459 | { | |
9a919f12 | 460 | // write data of single measurement. Note: the records ignore regrouping, store direct parameters |
8a9ab0eb | 461 | if (lSigma<=0.0) { // If parameter is fixed, then no equation |
462 | for (int i=nlc;i--;) derlc[i] = 0.0; | |
463 | for (int i=ngb;i--;) dergb[i] = 0.0; | |
464 | return; | |
465 | } | |
466 | // | |
a8c5b94c | 467 | if (fRecFileStatus<2) InitDataRecStorage(); // create a buffer to store the data |
8a9ab0eb | 468 | // |
469 | fRecord->AddResidual(lMeas); | |
470 | // | |
471 | // Retrieve local param interesting indices | |
cc9bec47 | 472 | for (int i=0;i<nlc;i++) if (!IsZero(derlc[i])) {fRecord->AddIndexValue(indlc[i],derlc[i]); derlc[i]=0.; indlc[i]=0;} |
8a9ab0eb | 473 | // |
474 | fRecord->AddWeight( 1./lSigma/lSigma ); | |
475 | // | |
476 | // Idem for global parameters | |
cc9bec47 | 477 | for (int i=0;i<ngb;i++) if (!IsZero(dergb[i])) {fRecord->AddIndexValue(indgb[i],dergb[i]); dergb[i]=0.; indgb[i]=0;} |
8a9ab0eb | 478 | // |
479 | } | |
480 | ||
de34b538 | 481 | |
8a9ab0eb | 482 | //_____________________________________________________________________________________________ |
339fbe23 | 483 | void AliMillePede2::SetGlobalConstraint(const double *dergb, double val, double sigma) |
8a9ab0eb | 484 | { |
485 | // Define a constraint equation. | |
486 | if (!fConsRecFile || !fConsRecFile->IsWritable()) InitConsRecStorage(); // create a buffer to store the data | |
487 | // | |
488 | fRecord->Reset(); | |
489 | fRecord->AddResidual(val); | |
de34b538 | 490 | fRecord->AddWeight(sigma); |
9a919f12 | 491 | for (int i=0; i<fNGloParIni; i++) if (!IsZero(dergb[i])) fRecord->AddIndexValue(i,dergb[i]); |
8a9ab0eb | 492 | fNGloConstraints++; |
cc9bec47 | 493 | if (IsZero(sigma)) fNLagrangeConstraints++; |
de34b538 | 494 | // printf("NewConstraint:\n"); fRecord->Print(); //RRR |
8a9ab0eb | 495 | SaveRecordConstraint(); |
496 | // | |
497 | } | |
498 | ||
499 | //_____________________________________________________________________________________________ | |
339fbe23 | 500 | void AliMillePede2::SetGlobalConstraint(const int *indgb, const double *dergb, int ngb, double val,double sigma) |
8a9ab0eb | 501 | { |
502 | // Define a constraint equation. | |
503 | if (!fConsRecFile || !fConsRecFile->IsWritable()) InitConsRecStorage(); // create a buffer to store the data | |
504 | fRecord->Reset(); | |
505 | fRecord->AddResidual(val); | |
de34b538 | 506 | fRecord->AddWeight(sigma); // dummy |
cc9bec47 | 507 | for (int i=0; i<ngb; i++) if (!IsZero(dergb[i])) fRecord->AddIndexValue(indgb[i],dergb[i]); |
8a9ab0eb | 508 | fNGloConstraints++; |
cc9bec47 | 509 | if (IsZero(sigma)) fNLagrangeConstraints++; |
8a9ab0eb | 510 | SaveRecordConstraint(); |
511 | // | |
512 | } | |
513 | ||
514 | //_____________________________________________________________________________________________ | |
515 | Int_t AliMillePede2::LocalFit(double *localParams) | |
516 | { | |
517 | /* | |
518 | Perform local parameters fit once all the local equations have been set | |
519 | ----------------------------------------------------------- | |
520 | localParams = (if !=0) will contain the fitted track parameters and | |
521 | related errors | |
522 | */ | |
de34b538 | 523 | static int nrefSize = 0; |
524 | // static TArrayI refLoc,refGlo,nrefLoc,nrefGlo; | |
525 | static Int_t *refLoc=0,*refGlo=0,*nrefLoc=0,*nrefGlo=0; | |
8a9ab0eb | 526 | int nPoints = 0; |
527 | // | |
de34b538 | 528 | AliSymMatrix &matCLoc = *fMatCLoc; |
529 | AliMatrixSq &matCGlo = *fMatCGlo; | |
530 | AliRectMatrix &matCGloLoc = *fMatCGloLoc; | |
8a9ab0eb | 531 | // |
532 | memset(fVecBLoc,0,fNLocPar*sizeof(double)); | |
de34b538 | 533 | matCLoc.Reset(); |
8a9ab0eb | 534 | // |
535 | int cnt=0; | |
536 | int recSz = fRecord->GetSize(); | |
537 | // | |
538 | while(cnt<recSz) { // Transfer the measurement records to matrices | |
539 | // | |
540 | // extract addresses of residual, weight and pointers on local and global derivatives for each point | |
de34b538 | 541 | if (nrefSize<=nPoints) { |
542 | int *tmpA = 0; | |
543 | nrefSize = 2*(nPoints+1); | |
544 | tmpA = refLoc; refLoc = new Int_t[nrefSize]; if (tmpA) memcpy(refLoc,tmpA,nPoints*sizeof(int)); | |
545 | tmpA = refGlo; refGlo = new Int_t[nrefSize]; if (tmpA) memcpy(refGlo,tmpA,nPoints*sizeof(int)); | |
546 | tmpA = nrefLoc; nrefLoc = new Int_t[nrefSize]; if (tmpA) memcpy(nrefLoc,tmpA,nPoints*sizeof(int)); | |
547 | tmpA = nrefGlo; nrefGlo = new Int_t[nrefSize]; if (tmpA) memcpy(nrefGlo,tmpA,nPoints*sizeof(int)); | |
8a9ab0eb | 548 | } |
549 | // | |
de34b538 | 550 | refLoc[nPoints] = ++cnt; |
8a9ab0eb | 551 | int nLoc = 0; |
552 | while(!fRecord->IsWeight(cnt)) {nLoc++; cnt++;} | |
de34b538 | 553 | nrefLoc[nPoints] = nLoc; |
8a9ab0eb | 554 | // |
de34b538 | 555 | refGlo[nPoints] = ++cnt; |
8a9ab0eb | 556 | int nGlo = 0; |
557 | while(!fRecord->IsResidual(cnt) && cnt<recSz) {nGlo++; cnt++;} | |
de34b538 | 558 | nrefGlo[nPoints] = nGlo; |
8a9ab0eb | 559 | // |
560 | nPoints++; | |
561 | } | |
c9f619c7 | 562 | if (fMinRecordLength>0 && nPoints < fMinRecordLength) return 0; // ignore |
563 | // | |
8a9ab0eb | 564 | double vl; |
565 | // | |
c9f619c7 | 566 | double gloWgh = fRunWgh; |
567 | if (fUseRecordWeight) gloWgh *= fRecord->GetWeight(); // global weight for this set | |
7c3070ec | 568 | Int_t maxLocUsed = 0; |
569 | // | |
8a9ab0eb | 570 | for (int ip=nPoints;ip--;) { // Transfer the measurement records to matrices |
de34b538 | 571 | double resid = fRecord->GetValue( refLoc[ip]-1 ); |
c9f619c7 | 572 | double weight = fRecord->GetValue( refGlo[ip]-1 )*gloWgh; |
573 | int odd = (ip&0x1); | |
574 | if (fWghScl[odd]>0) weight *= fWghScl[odd]; | |
de34b538 | 575 | double *derLoc = fRecord->GetValue()+refLoc[ip]; |
576 | double *derGlo = fRecord->GetValue()+refGlo[ip]; | |
577 | int *indLoc = fRecord->GetIndex()+refLoc[ip]; | |
578 | int *indGlo = fRecord->GetIndex()+refGlo[ip]; | |
579 | // | |
a8c5b94c | 580 | for (int i=nrefGlo[ip];i--;) { // suppress the global part (only relevant with iterations) |
9a919f12 | 581 | // |
582 | // if regrouping was requested, do it here | |
583 | if (fkReGroup) { | |
584 | int idtmp = fkReGroup[ indGlo[i] ]; | |
585 | if (idtmp == kFixParID) indGlo[i] = kFixParID; // fixed param in regrouping | |
ebf239ae | 586 | else indGlo[i] = idtmp; |
9a919f12 | 587 | } |
588 | // | |
8a9ab0eb | 589 | int iID = indGlo[i]; // Global param indice |
9a919f12 | 590 | if (iID<0 || fSigmaPar[iID]<=0.) continue; // fixed parameter RRRCheck |
8a9ab0eb | 591 | if (fIsLinear[iID]) resid -= derGlo[i]*(fInitPar[iID]+fDeltaPar[iID]); // linear parameter |
592 | else resid -= derGlo[i]*fDeltaPar[iID]; // nonlinear parameter | |
593 | } | |
594 | // | |
de34b538 | 595 | // Symmetric matrix, don't bother j>i coeffs |
596 | for (int i=nrefLoc[ip];i--;) { // Fill local matrix and vector | |
597 | fVecBLoc[ indLoc[i] ] += weight*resid*derLoc[i]; | |
7c3070ec | 598 | if (indLoc[i]>maxLocUsed) maxLocUsed = indLoc[i]; |
de34b538 | 599 | for (int j=i+1;j--;) matCLoc(indLoc[i] ,indLoc[j]) += weight*derLoc[i]*derLoc[j]; |
7c3070ec | 600 | } |
8a9ab0eb | 601 | // |
602 | } // end of the transfer of the measurement record to matrices | |
603 | // | |
7c3070ec | 604 | matCLoc.SetSizeUsed(++maxLocUsed); // data with B=0 may use less than declared nLocals |
605 | // | |
f9efbbf6 | 606 | /* //RRR |
607 | fRecord->Print("l"); | |
608 | printf("\nBefore\nLocalMatrix: "); matCLoc.Print("l"); | |
609 | printf("RHSLoc: "); for (int i=0;i<fNLocPar;i++) printf("%+e |",fVecBLoc[i]); printf("\n"); | |
610 | */ | |
8a9ab0eb | 611 | // first try to solve by faster Cholesky decomposition, then by Gaussian elimination |
de34b538 | 612 | if (!matCLoc.SolveChol(fVecBLoc,kTRUE)) { |
a8c5b94c | 613 | AliInfo("Failed to solve locals by Cholesky, trying Gaussian Elimination"); |
de34b538 | 614 | if (!matCLoc.SolveSpmInv(fVecBLoc,kTRUE)) { |
8a9ab0eb | 615 | AliInfo("Failed to solve locals by Gaussian Elimination, skip..."); |
de34b538 | 616 | matCLoc.Print("d"); |
8a9ab0eb | 617 | return 0; // failed to solve |
618 | } | |
619 | } | |
620 | // | |
621 | // If requested, store the track params and errors | |
f9efbbf6 | 622 | //RRR printf("locfit: "); for (int i=0;i<fNLocPar;i++) printf("%+e |",fVecBLoc[i]); printf("\n"); |
623 | ||
7c3070ec | 624 | if (localParams) for (int i=maxLocUsed; i--;) { |
8a9ab0eb | 625 | localParams[2*i] = fVecBLoc[i]; |
de34b538 | 626 | localParams[2*i+1] = TMath::Sqrt(TMath::Abs(matCLoc.QueryDiag(i))); |
8a9ab0eb | 627 | } |
628 | // | |
629 | float lChi2 = 0; | |
630 | int nEq = 0; | |
631 | // | |
632 | for (int ip=nPoints;ip--;) { // Calculate residuals | |
de34b538 | 633 | double resid = fRecord->GetValue( refLoc[ip]-1 ); |
c9f619c7 | 634 | double weight = fRecord->GetValue( refGlo[ip]-1 )*gloWgh; |
635 | int odd = (ip&0x1); | |
636 | if (fWghScl[odd]>0) weight *= fWghScl[odd]; | |
de34b538 | 637 | double *derLoc = fRecord->GetValue()+refLoc[ip]; |
638 | double *derGlo = fRecord->GetValue()+refGlo[ip]; | |
639 | int *indLoc = fRecord->GetIndex()+refLoc[ip]; | |
640 | int *indGlo = fRecord->GetIndex()+refGlo[ip]; | |
8a9ab0eb | 641 | // |
642 | // Suppress local and global contribution in residuals; | |
de34b538 | 643 | for (int i=nrefLoc[ip];i--;) resid -= derLoc[i]*fVecBLoc[ indLoc[i] ]; // local part |
8a9ab0eb | 644 | // |
de34b538 | 645 | for (int i=nrefGlo[ip];i--;) { // global part |
8a9ab0eb | 646 | int iID = indGlo[i]; |
9a919f12 | 647 | if ( iID<0 || fSigmaPar[iID] <= 0.) continue; // fixed parameter RRRCheck |
8a9ab0eb | 648 | if (fIsLinear[iID]) resid -= derGlo[i]*(fInitPar[iID]+fDeltaPar[iID]); // linear parameter |
649 | else resid -= derGlo[i]*fDeltaPar[iID]; // nonlinear parameter | |
650 | } | |
651 | // | |
652 | // reject the track if the residual is too large (outlier) | |
de34b538 | 653 | double absres = TMath::Abs(resid); |
654 | if ( (absres >= fResCutInit && fIter ==1 ) || | |
655 | (absres >= fResCut && fIter > 1)) { | |
656 | if (fLocFitAdd) fNLocFitsRejected++; | |
cc9bec47 | 657 | // printf("reject res %5ld %+e\n",fCurrRecDataID,resid); |
8a9ab0eb | 658 | return 0; |
659 | } | |
660 | // | |
661 | lChi2 += weight*resid*resid ; // total chi^2 | |
662 | nEq++; // number of equations | |
663 | } // end of Calculate residuals | |
664 | // | |
a8c5b94c | 665 | lChi2 /= gloWgh; |
7c3070ec | 666 | int nDoF = nEq-maxLocUsed; |
8a9ab0eb | 667 | lChi2 = (nDoF>0) ? lChi2/nDoF : 0; // Chi^2/dof |
668 | // | |
669 | if (fNStdDev != 0 && nDoF>0 && lChi2 > Chi2DoFLim(fNStdDev,nDoF)*fChi2CutFactor) { // check final chi2 | |
de34b538 | 670 | if (fLocFitAdd) fNLocFitsRejected++; |
cc9bec47 | 671 | // printf("reject chi2 %5ld: %+e\n",fCurrRecDataID, lChi2); |
8a9ab0eb | 672 | return 0; |
673 | } | |
674 | // | |
de34b538 | 675 | if (fLocFitAdd) { |
676 | fNLocFits++; | |
677 | fNLocEquations += nEq; | |
678 | } | |
679 | else { | |
680 | fNLocFits--; | |
681 | fNLocEquations -= nEq; | |
682 | } | |
8a9ab0eb | 683 | // |
684 | // local operations are finished, track is accepted | |
685 | // We now update the global parameters (other matrices) | |
686 | // | |
687 | int nGloInFit = 0; | |
688 | // | |
689 | for (int ip=nPoints;ip--;) { // Update matrices | |
de34b538 | 690 | double resid = fRecord->GetValue( refLoc[ip]-1 ); |
c9f619c7 | 691 | double weight = fRecord->GetValue( refGlo[ip]-1 )*gloWgh; |
692 | int odd = (ip&0x1); | |
693 | if (fWghScl[odd]>0) weight *= fWghScl[odd]; | |
de34b538 | 694 | double *derLoc = fRecord->GetValue()+refLoc[ip]; |
695 | double *derGlo = fRecord->GetValue()+refGlo[ip]; | |
696 | int *indLoc = fRecord->GetIndex()+refLoc[ip]; | |
697 | int *indGlo = fRecord->GetIndex()+refGlo[ip]; | |
698 | // | |
9a919f12 | 699 | for (int i=nrefGlo[ip];i--;) { // suppress the global part |
8a9ab0eb | 700 | int iID = indGlo[i]; // Global param indice |
9a919f12 | 701 | if ( iID<0 || fSigmaPar[iID] <= 0.) continue; // fixed parameter RRRCheck |
de34b538 | 702 | if (fIsLinear[iID]) resid -= derGlo[i]*(fInitPar[iID]+fDeltaPar[iID]); // linear parameter |
703 | else resid -= derGlo[i]*fDeltaPar[iID]; // nonlinear parameter | |
8a9ab0eb | 704 | } |
705 | // | |
de34b538 | 706 | for (int ig=nrefGlo[ip];ig--;) { |
8a9ab0eb | 707 | int iIDg = indGlo[ig]; // Global param indice (the matrix line) |
9a919f12 | 708 | if ( iIDg<0 || fSigmaPar[iIDg] <= 0.) continue; // fixed parameter RRRCheck |
de34b538 | 709 | if (fLocFitAdd) fVecBGlo[ iIDg ] += weight*resid*derGlo[ig]; //!!! |
710 | else fVecBGlo[ iIDg ] -= weight*resid*derGlo[ig]; //!!! | |
8a9ab0eb | 711 | // |
712 | // First of all, the global/global terms (exactly like local matrix) | |
de34b538 | 713 | int nfill = 0; |
714 | for (int jg=ig+1;jg--;) { // matCGlo is symmetric by construction | |
8a9ab0eb | 715 | int jIDg = indGlo[jg]; |
9a919f12 | 716 | if ( jIDg<0 || fSigmaPar[jIDg] <= 0.) continue; // fixed parameter RRRCheck |
cc9bec47 | 717 | if ( !IsZero(vl = weight*derGlo[ig]*derGlo[jg]) ) { |
de34b538 | 718 | fFillIndex[nfill] = jIDg; |
719 | fFillValue[nfill++] = fLocFitAdd ? vl:-vl; | |
720 | } | |
721 | } | |
722 | if (nfill) matCGlo.AddToRow(iIDg,fFillValue,fFillIndex,nfill); | |
8a9ab0eb | 723 | // |
724 | // Now we have also rectangular matrices containing global/local terms. | |
725 | int iCIDg = fGlo2CGlo[iIDg]; // compressed Index of index | |
726 | if (iCIDg == -1) { | |
de34b538 | 727 | Double_t *rowGL = matCGloLoc(nGloInFit); |
7c3070ec | 728 | for (int k=maxLocUsed;k--;) rowGL[k] = 0.0; // reset the row |
8a9ab0eb | 729 | iCIDg = fGlo2CGlo[iIDg] = nGloInFit; |
730 | fCGlo2Glo[nGloInFit++] = iIDg; | |
731 | } | |
732 | // | |
de34b538 | 733 | Double_t *rowGLIDg = matCGloLoc(iCIDg); |
734 | for (int il=nrefLoc[ip];il--;) rowGLIDg[ indLoc[il] ] += weight*derGlo[ig]*derLoc[il]; | |
735 | fProcPnt[iIDg] += fLocFitAdd ? 1:-1; // update counter | |
8a9ab0eb | 736 | // |
737 | } | |
738 | } // end of Update matrices | |
739 | // | |
f9efbbf6 | 740 | /*//RRR |
741 | printf("After GLO\n"); | |
742 | printf("MatCLoc: "); fMatCLoc->Print("l"); | |
743 | printf("MatCGlo: "); fMatCGlo->Print("l"); | |
744 | printf("MatCGlLc:"); fMatCGloLoc->Print("l"); | |
745 | printf("BGlo: "); for (int i=0; i<fNGloPar; i++) printf("%+e |",fVecBGlo[i]); printf("\n"); | |
746 | */ | |
8a9ab0eb | 747 | // calculate fMatCGlo -= fMatCGloLoc * fMatCLoc * fMatCGloLoc^T |
748 | // and fVecBGlo -= fMatCGloLoc * fVecBLoc | |
749 | // | |
de34b538 | 750 | //-------------------------------------------------------------- >>> |
8a9ab0eb | 751 | double vll; |
752 | for (int iCIDg=0; iCIDg<nGloInFit; iCIDg++) { | |
753 | int iIDg = fCGlo2Glo[iCIDg]; | |
754 | // | |
755 | vl = 0; | |
de34b538 | 756 | Double_t *rowGLIDg = matCGloLoc(iCIDg); |
7c3070ec | 757 | for (int kl=0;kl<maxLocUsed;kl++) if (rowGLIDg[kl]) vl += rowGLIDg[kl]*fVecBLoc[kl]; |
cc9bec47 | 758 | if (!IsZero(vl)) fVecBGlo[iIDg] -= fLocFitAdd ? vl : -vl; |
8a9ab0eb | 759 | // |
de34b538 | 760 | int nfill = 0; |
8a9ab0eb | 761 | for (int jCIDg=0;jCIDg<=iCIDg; jCIDg++) { |
762 | int jIDg = fCGlo2Glo[jCIDg]; | |
763 | // | |
764 | vl = 0; | |
de34b538 | 765 | Double_t *rowGLJDg = matCGloLoc(jCIDg); |
7c3070ec | 766 | for (int kl=0;kl<maxLocUsed;kl++) { |
8a9ab0eb | 767 | // diag terms |
cc9bec47 | 768 | if ( (!IsZero(vll=rowGLIDg[kl]*rowGLJDg[kl])) ) vl += matCLoc.QueryDiag(kl)*vll; |
8a9ab0eb | 769 | // |
770 | // off-diag terms | |
771 | for (int ll=0;ll<kl;ll++) { | |
cc9bec47 | 772 | if ( !IsZero(vll=rowGLIDg[kl]*rowGLJDg[ll]) ) vl += matCLoc(kl,ll)*vll; |
773 | if ( !IsZero(vll=rowGLIDg[ll]*rowGLJDg[kl]) ) vl += matCLoc(kl,ll)*vll; | |
8a9ab0eb | 774 | } |
775 | } | |
cc9bec47 | 776 | if (!IsZero(vl)) { |
de34b538 | 777 | fFillIndex[nfill] = jIDg; |
778 | fFillValue[nfill++] = fLocFitAdd ? -vl : vl; | |
779 | } | |
8a9ab0eb | 780 | } |
de34b538 | 781 | if (nfill) matCGlo.AddToRow(iIDg,fFillValue,fFillIndex,nfill); |
8a9ab0eb | 782 | } |
783 | // | |
784 | // reset compressed index array | |
785 | // | |
f9efbbf6 | 786 | /*//RRR |
787 | printf("After GLOLoc\n"); | |
788 | printf("MatCGlo: "); fMatCGlo->Print(""); | |
789 | printf("BGlo: "); for (int i=0; i<fNGloPar; i++) printf("%+e |",fVecBGlo[i]); printf("\n"); | |
790 | */ | |
de34b538 | 791 | for (int i=nGloInFit;i--;) { |
792 | fGlo2CGlo[ fCGlo2Glo[i] ] = -1; | |
793 | fCGlo2Glo[i] = -1; | |
794 | } | |
8a9ab0eb | 795 | // |
de34b538 | 796 | //---------------------------------------------------- <<< |
8a9ab0eb | 797 | return 1; |
798 | } | |
799 | ||
800 | //_____________________________________________________________________________________________ | |
801 | Int_t AliMillePede2::GlobalFit(Double_t *par, Double_t *error, Double_t *pull) | |
802 | { | |
803 | // performs a requested number of global iterations | |
804 | fIter = 1; | |
805 | // | |
806 | TStopwatch sw; sw.Start(); | |
807 | // | |
808 | int res = 0; | |
809 | AliInfo("Starting Global fit."); | |
810 | while (fIter<=fMaxIter) { | |
811 | // | |
812 | res = GlobalFitIteration(); | |
813 | if (!res) break; | |
814 | // | |
cc9bec47 | 815 | if (!IsZero(fChi2CutFactor-fChi2CutRef)) { |
8a9ab0eb | 816 | fChi2CutFactor = TMath::Sqrt(fChi2CutFactor); |
817 | if (fChi2CutFactor < 1.2*fChi2CutRef) { | |
818 | fChi2CutFactor = fChi2CutRef; | |
e30a812f | 819 | //RRR fIter = fMaxIter - 1; // Last iteration |
8a9ab0eb | 820 | } |
821 | } | |
822 | fIter++; | |
823 | } | |
824 | // | |
825 | sw.Stop(); | |
826 | AliInfo(Form("Global fit %s, CPU time: %.1f",res ? "Converged":"Failed",sw.CpuTime())); | |
827 | if (!res) return 0; | |
828 | // | |
9a919f12 | 829 | if (par) for (int i=fNGloParIni;i--;) par[i] = GetFinalParam(i); |
8a9ab0eb | 830 | // |
339fbe23 | 831 | if (fGloSolveStatus==kInvert) { // errors on params are available |
9a919f12 | 832 | if (error) for (int i=fNGloParIni;i--;) error[i] = GetFinalError(i); |
833 | if (pull) for (int i=fNGloParIni;i--;) pull[i] = GetPull(i); | |
8a9ab0eb | 834 | } |
835 | // | |
836 | return 1; | |
837 | } | |
838 | ||
839 | //_____________________________________________________________________________________________ | |
840 | Int_t AliMillePede2::GlobalFitIteration() | |
841 | { | |
842 | // perform global parameters fit once all the local equations have been fitted | |
843 | // | |
844 | AliInfo(Form("Global Fit Iteration#%2d (Local Fit Chi^2 cut factor: %.2f)",fIter,fChi2CutFactor)); | |
845 | // | |
846 | if (!fNGloPar || !fTreeData) { | |
339fbe23 | 847 | AliInfo("No data was stored, stopping iteration"); |
8a9ab0eb | 848 | return 0; |
849 | } | |
de34b538 | 850 | TStopwatch sw,sws; |
851 | sw.Start(); | |
852 | sws.Stop(); | |
8a9ab0eb | 853 | // |
854 | if (!fConstrUsed) { | |
855 | fConstrUsed = new Bool_t[fNGloConstraints]; | |
856 | memset(fConstrUsed,0,fNGloConstraints*sizeof(Bool_t)); | |
857 | } | |
858 | // Reset all info specific for this step | |
de34b538 | 859 | AliMatrixSq& matCGlo = *fMatCGlo; |
860 | matCGlo.Reset(); | |
8a9ab0eb | 861 | memset(fProcPnt,0,fNGloPar*sizeof(Int_t)); |
862 | // | |
863 | fNGloConstraints = fTreeConstr ? fTreeConstr->GetEntries() : 0; | |
864 | // | |
de34b538 | 865 | // count number of Lagrange constraints: they need new row/cols to be added |
866 | fNLagrangeConstraints = 0; | |
867 | for (int i=0; i<fNGloConstraints; i++) { | |
868 | ReadRecordConstraint(i); | |
cc9bec47 | 869 | if ( IsZero(fRecord->GetValue(1)) ) fNLagrangeConstraints++; // exact constraint (no error) -> Lagrange multiplier |
de34b538 | 870 | } |
871 | // | |
8a9ab0eb | 872 | // if needed, readjust the size of the global vector (for matrices this is done automatically) |
de34b538 | 873 | if (!fVecBGlo || fNGloSize!=fNGloPar+fNLagrangeConstraints) { |
8a9ab0eb | 874 | delete[] fVecBGlo; // in case some constraint was added between the two manual iterations |
de34b538 | 875 | fNGloSize = fNGloPar+fNLagrangeConstraints; |
8a9ab0eb | 876 | fVecBGlo = new Double_t[fNGloSize]; |
877 | } | |
878 | memset(fVecBGlo,0,fNGloSize*sizeof(double)); | |
879 | // | |
880 | fNLocFits = 0; | |
881 | fNLocFitsRejected = 0; | |
de34b538 | 882 | fNLocEquations = 0; |
8a9ab0eb | 883 | // |
884 | // Process data records and build the matrices | |
885 | Long_t ndr = fTreeData->GetEntries(); | |
e30a812f | 886 | Long_t first = fSelFirst>0 ? fSelFirst : 0; |
887 | Long_t last = fSelLast<1 ? ndr : (fSelLast>=ndr ? ndr : fSelLast+Long_t(1)); | |
888 | ndr = last - first; | |
889 | // | |
73bbf779 | 890 | AliInfo(Form("Building the Global matrix from data records %ld : %ld",first,last)); |
e30a812f | 891 | if (ndr<1) return 0; |
de34b538 | 892 | // |
893 | TStopwatch swt; swt.Start(); | |
894 | fLocFitAdd = kTRUE; // add contributions of matching tracks | |
8a9ab0eb | 895 | for (Long_t i=0;i<ndr;i++) { |
e30a812f | 896 | Long_t iev = i+first; |
897 | ReadRecordData(iev); | |
898 | if (!IsRecordAcceptable()) continue; | |
8a9ab0eb | 899 | LocalFit(); |
7c3070ec | 900 | if ( (i%int(0.2*ndr)) == 0) printf("%.1f%% of local fits done\n", double(100.*i)/ndr); |
8a9ab0eb | 901 | } |
de34b538 | 902 | swt.Stop(); |
903 | printf("%ld local fits done: ", ndr); | |
f9efbbf6 | 904 | /* |
905 | printf("MatCGlo: "); fMatCGlo->Print("l"); | |
906 | printf("BGlo: "); for (int i=0; i<fNGloPar; i++) printf("%+e |",fVecBGlo[i]); printf("\n"); | |
de34b538 | 907 | swt.Print(); |
f9efbbf6 | 908 | */ |
de34b538 | 909 | sw.Start(kFALSE); |
910 | // | |
8a9ab0eb | 911 | // |
de34b538 | 912 | // ---------------------- Reject parameters with low statistics ------------>> |
8a9ab0eb | 913 | fNGloFix = 0; |
de34b538 | 914 | if (fMinPntValid>1 && fNGroupsSet) { |
915 | // | |
916 | printf("Checking parameters with statistics < %d\n",fMinPntValid); | |
917 | TStopwatch swsup; | |
918 | swsup.Start(); | |
919 | // 1) build the list of parameters to fix | |
920 | Int_t fixArrSize = 10; | |
921 | Int_t nFixedGroups = 0; | |
922 | TArrayI fixGroups(fixArrSize); | |
923 | // | |
91e1aa2c | 924 | int grIDold = -2; |
925 | int oldStart = -1; | |
926 | double oldMin = 1.e20; | |
927 | double oldMax =-1.e20; | |
928 | // | |
de34b538 | 929 | for (int i=fNGloPar;i--;) { // // Reset row and column of fixed params and add 1/sig^2 to free ones |
930 | int grID = fParamGrID[i]; | |
931 | if (grID<0) continue; // not in the group | |
91e1aa2c | 932 | // |
933 | if (grID!=grIDold) { // starting new group | |
934 | if (grIDold>=0) { // decide if the group has enough statistics | |
935 | if (oldMin<fMinPntValid && oldMax<2*fMinPntValid) { // suppress group | |
936 | for (int iold=oldStart;iold>i;iold--) fProcPnt[iold] = 0; | |
937 | Bool_t fnd = kFALSE; // check if the group is already accounted | |
938 | for (int j=nFixedGroups;j--;) if (fixGroups[j]==grIDold) {fnd=kTRUE; break;} | |
939 | if (!fnd) { | |
940 | if (nFixedGroups>=fixArrSize) {fixArrSize*=2; fixGroups.Set(fixArrSize);} | |
941 | fixGroups[nFixedGroups++] = grIDold; // add group to fix | |
942 | } | |
943 | } | |
944 | } | |
945 | grIDold = grID; // mark the start of the new group | |
946 | oldStart = i; | |
947 | oldMin = 1.e20; | |
948 | oldMax = -1.e20; | |
949 | } | |
950 | if (oldMin>fProcPnt[i]) oldMin = fProcPnt[i]; | |
951 | if (oldMax<fProcPnt[i]) oldMax = fProcPnt[i]; | |
de34b538 | 952 | // |
953 | } | |
91e1aa2c | 954 | // extra check for the last group |
955 | if (grIDold>=0 && oldMin<fMinPntValid && oldMax<2*fMinPntValid) { // suppress group | |
956 | for (int iold=oldStart;iold--;) fProcPnt[iold] = 0; | |
957 | Bool_t fnd = kFALSE; // check if the group is already accounted | |
958 | for (int j=nFixedGroups;j--;) if (fixGroups[j]==grIDold) {fnd=kTRUE; break;} | |
959 | if (!fnd) { | |
960 | if (nFixedGroups>=fixArrSize) {fixArrSize*=2; fixGroups.Set(fixArrSize);} | |
961 | fixGroups[nFixedGroups++] = grIDold; // add group to fix | |
962 | } | |
963 | } | |
964 | // | |
de34b538 | 965 | // 2) loop over records and add contributions of fixed groups with negative sign |
966 | fLocFitAdd = kFALSE; | |
967 | // | |
968 | for (Long_t i=0;i<ndr;i++) { | |
e30a812f | 969 | Long_t iev = i+first; |
970 | ReadRecordData(iev); | |
971 | if (!IsRecordAcceptable()) continue; | |
de34b538 | 972 | Bool_t suppr = kFALSE; |
973 | for (int ifx=nFixedGroups;ifx--;)if (fRecord->IsGroupPresent(fixGroups[ifx])) suppr = kTRUE; | |
974 | if (suppr) LocalFit(); | |
975 | } | |
976 | fLocFitAdd = kTRUE; | |
977 | // | |
978 | if (nFixedGroups) { | |
979 | printf("Suppressed contributions of groups with NPoints<%d :\n",fMinPntValid); | |
980 | for (int i=0;i<nFixedGroups;i++) printf("%d ",fixGroups[i]); printf("\n"); | |
981 | } | |
982 | swsup.Stop(); | |
983 | swsup.Print(); | |
984 | } | |
985 | // ---------------------- Reject parameters with low statistics ------------<< | |
986 | // | |
987 | // add large number to diagonal of fixed params | |
8a9ab0eb | 988 | // |
8a9ab0eb | 989 | for (int i=fNGloPar;i--;) { // // Reset row and column of fixed params and add 1/sig^2 to free ones |
7c3070ec | 990 | // printf("#%3d : Nproc : %5d grp: %d\n",i,fProcPnt[i],fParamGrID[i]); |
de34b538 | 991 | if (fProcPnt[i]<1) { |
8a9ab0eb | 992 | fNGloFix++; |
8a9ab0eb | 993 | fVecBGlo[i] = 0.; |
a8c5b94c | 994 | matCGlo.DiagElem(i) = 1.;//float(fNLocEquations*fNLocEquations); |
e30a812f | 995 | // matCGlo.DiagElem(i) = float(fNLocEquations*fNLocEquations); |
8a9ab0eb | 996 | } |
de34b538 | 997 | else matCGlo.DiagElem(i) += (fgWeightSigma ? fProcPnt[i] : 1.)/(fSigmaPar[i]*fSigmaPar[i]); |
8a9ab0eb | 998 | } |
999 | // | |
de34b538 | 1000 | for (int i=fNGloPar;i--;) fDiagCGlo[i] = matCGlo.QueryDiag(i); // save the diagonal elements |
1001 | // | |
8a9ab0eb | 1002 | // add constraint equations |
1003 | int nVar = fNGloPar; // Current size of global matrix | |
1004 | for (int i=0; i<fNGloConstraints; i++) { | |
1005 | ReadRecordConstraint(i); | |
1006 | double val = fRecord->GetValue(0); | |
de34b538 | 1007 | double sig = fRecord->GetValue(1); |
8a9ab0eb | 1008 | int *indV = fRecord->GetIndex()+2; |
1009 | double *der = fRecord->GetValue()+2; | |
1010 | int csize = fRecord->GetSize()-2; | |
1011 | // | |
9a919f12 | 1012 | if (fkReGroup) { |
1013 | for (int jp=csize;jp--;) { | |
1014 | int idp = indV[jp]; | |
1015 | if (fkReGroup[idp]<0) AliFatal(Form("Constain is requested for suppressed parameter #%d",indV[jp])); | |
1016 | indV[jp] = idp; | |
1017 | } | |
1018 | } | |
de34b538 | 1019 | // check if after suppression of fixed variables there are non-0 derivatives |
1020 | // and determine the max statistics of involved params | |
1021 | int nSuppressed = 0; | |
1022 | int maxStat = 1; | |
1023 | for (int j=csize;j--;) { | |
1024 | if (fProcPnt[indV[j]]<1) nSuppressed++; | |
1025 | else { | |
1026 | maxStat = TMath::Max(maxStat,fProcPnt[indV[j]]); | |
1027 | } | |
1028 | } | |
8a9ab0eb | 1029 | // |
de34b538 | 1030 | if (nSuppressed==csize) { |
1031 | // AliInfo(Form("Neglecting constraint %d of %d derivatives since no free parameters left",i,csize)); | |
8a9ab0eb | 1032 | // |
1033 | // was this constraint ever created ? | |
de34b538 | 1034 | if ( sig==0 && fConstrUsed[i] ) { // this is needed only for constraints with Lagrange multiplier |
8a9ab0eb | 1035 | // to avoid empty row impose dummy constraint on "Lagrange multiplier" |
de34b538 | 1036 | matCGlo.DiagElem(nVar) = 1.; |
8a9ab0eb | 1037 | fVecBGlo[nVar++] = 0; |
1038 | } | |
1039 | continue; | |
1040 | } | |
1041 | // | |
de34b538 | 1042 | // account for already accumulated corrections |
1043 | for (int j=csize; j--;) val -= der[j]*(fInitPar[ indV[j] ]+fDeltaPar[ indV[j] ]); | |
8a9ab0eb | 1044 | // |
de34b538 | 1045 | if (sig > 0) { // this is a gaussian constriant: no Lagrange multipliers are added |
1046 | // | |
1047 | double sig2i = (fgWeightSigma ? TMath::Sqrt(maxStat) : 1.)/sig/sig; | |
1048 | for (int ir=0;ir<csize;ir++) { | |
1049 | int iID = indV[ir]; | |
1050 | for (int ic=0;ic<=ir;ic++) { // matrix is symmetric | |
1051 | int jID = indV[ic]; | |
1052 | double vl = der[ir]*der[ic]*sig2i; | |
cc9bec47 | 1053 | if (!IsZero(vl)) matCGlo(iID,jID) += vl; |
de34b538 | 1054 | } |
1055 | fVecBGlo[iID] += val*der[ir]*sig2i; | |
1056 | } | |
1057 | } | |
1058 | else { // this is exact constriant: Lagrange multipliers must be added | |
1059 | for (int j=csize; j--;) { | |
1060 | int jID = indV[j]; | |
1061 | if (fProcPnt[jID]<1) continue; // this parameter was fixed, don't put it into constraint | |
1062 | matCGlo(nVar,jID) = float(fNLocEquations)*der[j]; // fMatCGlo is symmetric, only lower triangle is filled | |
1063 | } | |
1064 | // | |
1065 | if (matCGlo.QueryDiag(nVar)) matCGlo.DiagElem(nVar) = 0.0; | |
1066 | fVecBGlo[nVar++] = float(fNLocEquations)*val; //RS ? should we use here fNLocFits ? | |
1067 | fConstrUsed[i] = kTRUE; | |
1068 | } | |
8a9ab0eb | 1069 | } |
1070 | // | |
1071 | AliInfo(Form("Obtained %-7ld equations from %-7ld records (%-7ld rejected). Fixed %-4d globals", | |
1072 | fNLocEquations,fNLocFits,fNLocFitsRejected,fNGloFix)); | |
1073 | ||
1074 | // | |
de34b538 | 1075 | sws.Start(); |
f9efbbf6 | 1076 | |
e041af50 | 1077 | #ifdef _DUMP_EQ_BEFORE_ |
1078 | const char* faildumpB = Form("mp2eq_before%d.dat",fIter); | |
1079 | int defoutB = dup(1); | |
c8396ea8 | 1080 | if (defoutB<0) {AliFatal("Failed on dup"); exit(1);} |
e041af50 | 1081 | int slvDumpB = open(faildumpB, O_RDWR|O_CREAT, 0666); |
1082 | if (slvDumpB>=0) { | |
1083 | dup2(slvDumpB,1); | |
1084 | printf("Solving%d for %d params\n",fIter,fNGloSize); | |
1085 | matCGlo.Print("10"); | |
1086 | for (int i=0;i<fNGloSize;i++) printf("b%2d : %+.10f\n",i,fVecBGlo[i]); | |
1087 | } | |
1088 | dup2(defoutB,1); | |
1089 | close(slvDumpB); | |
c8396ea8 | 1090 | close(defoutB); |
1091 | ||
e041af50 | 1092 | #endif |
f9efbbf6 | 1093 | /* |
1094 | printf("Solving:\n"); | |
ebf239ae | 1095 | matCGlo.Print("l"); |
f9efbbf6 | 1096 | for (int i=0;i<fNGloSize;i++) printf("b%2d : %+e\n",i,fVecBGlo[i]); |
1097 | */ | |
c9f619c7 | 1098 | #ifdef _DUMPEQ_BEFORE_ |
1099 | const char* faildumpB = Form("mp2eq_before%d.dat",fIter); | |
1100 | int defoutB = dup(1); | |
1101 | int slvDumpB = open(faildumpB, O_RDWR|O_CREAT, 0666); | |
1102 | dup2(slvDumpB,1); | |
1103 | // | |
1104 | printf("#Equation before step %d\n",fIter); | |
1105 | fMatCGlo->Print("10"); | |
1106 | printf("#RHS/STAT : NGlo:%d NGloSize:%d\n",fNGloPar,fNGloSize); | |
1107 | for (int i=0;i<fNGloSize;i++) printf("%d %+.10f %d\n",i,fVecBGlo[i],fProcPnt[i]); | |
1108 | // | |
1109 | dup2(defoutB,1); | |
1110 | close(slvDumpB); | |
c8396ea8 | 1111 | close(defoutB); |
c9f619c7 | 1112 | #endif |
1113 | // | |
8a9ab0eb | 1114 | fGloSolveStatus = SolveGlobalMatEq(); // obtain solution for this step |
c9f619c7 | 1115 | #ifdef _DUMPEQ_AFTER_ |
1116 | const char* faildumpA = Form("mp2eq_after%d.dat",fIter); | |
1117 | int defoutA = dup(1); | |
1118 | int slvDumpA = open(faildumpA, O_RDWR|O_CREAT, 0666); | |
1119 | dup2(slvDumpA,1); | |
1120 | // | |
1121 | printf("#Matrix after step %d\n",fIter); | |
1122 | fMatCGlo->Print("10"); | |
1123 | printf("#RHS/STAT : NGlo:%d NGloSize:%d\n",fNGloPar,fNGloSize); | |
1124 | for (int i=0;i<fNGloSize;i++) printf("%d %+.10f %d\n",i,fVecBGlo[i],fProcPnt[i]); | |
1125 | // | |
1126 | dup2(defoutA,1); | |
1127 | close(slvDumpA); | |
c8396ea8 | 1128 | close(defoutA); |
c9f619c7 | 1129 | #endif |
1130 | // | |
de34b538 | 1131 | sws.Stop(); |
1132 | printf("Solve %d |",fIter); sws.Print(); | |
8a9ab0eb | 1133 | // |
1134 | sw.Stop(); | |
339fbe23 | 1135 | AliInfo(Form("Iteration#%2d %s. CPU time: %.1f",fIter,fGloSolveStatus==kFailed ? "Failed":"Converged",sw.CpuTime())); |
1136 | if (fGloSolveStatus==kFailed) return 0; | |
8a9ab0eb | 1137 | // |
1138 | for (int i=fNGloPar;i--;) fDeltaPar[i] += fVecBGlo[i]; // Update global parameters values (for iterations) | |
e041af50 | 1139 | |
1140 | #ifdef _DUMP_EQ_AFTER_ | |
1141 | const char* faildumpA = Form("mp2eq_after%d.dat",fIter); | |
1142 | int defoutA = dup(1); | |
c8396ea8 | 1143 | if (defoutA<0) {AliFatal("Failed on dup"); exit(1);} |
e041af50 | 1144 | int slvDumpA = open(faildumpA, O_RDWR|O_CREAT, 0666); |
1145 | if (slvDumpA>=0) { | |
1146 | dup2(slvDumpA,1); | |
1147 | printf("Solving%d for %d params\n",fIter,fNGloSize); | |
1148 | matCGlo.Print("10"); | |
1149 | for (int i=0;i<fNGloSize;i++) printf("b%2d : %+.10f\n",i,fVecBGlo[i]); | |
1150 | } | |
1151 | dup2(defoutA,1); | |
1152 | close(slvDumpA); | |
c8396ea8 | 1153 | close(defoutA); |
e041af50 | 1154 | #endif |
8a9ab0eb | 1155 | // |
ebf239ae | 1156 | /* |
1157 | printf("Solved:\n"); | |
1158 | matCGlo.Print("l"); | |
1159 | for (int i=0;i<fNGloSize;i++) printf("b%2d : %+e (->%+e)\n",i,fVecBGlo[i], fDeltaPar[i]); | |
1160 | */ | |
1161 | ||
c9f619c7 | 1162 | PrintGlobalParameters(); |
8a9ab0eb | 1163 | return 1; |
1164 | } | |
1165 | ||
1166 | //_____________________________________________________________________________________________ | |
1167 | Int_t AliMillePede2::SolveGlobalMatEq() | |
1168 | { | |
1169 | // | |
1170 | // solve global matrix equation MatCGlob*X=VecBGlo and store the result in the VecBGlo | |
1171 | // | |
de34b538 | 1172 | /* |
1173 | printf("GlobalMatrix\n"); | |
ebf239ae | 1174 | fMatCGlo->Print("l"); |
de34b538 | 1175 | printf("RHS\n"); |
1176 | for (int i=0;i<fNGloPar;i++) printf("%d %+e\n",i,fVecBGlo[i]); | |
1177 | */ | |
1178 | // | |
8a9ab0eb | 1179 | if (!fgIsMatGloSparse) { |
1180 | // | |
de34b538 | 1181 | if (fNLagrangeConstraints==0) { // pos-def systems are faster to solve by Cholesky |
339fbe23 | 1182 | if ( ((AliSymMatrix*)fMatCGlo)->SolveChol(fVecBGlo, fgInvChol) ) return fgInvChol ? kInvert:kNoInversion; |
8a9ab0eb | 1183 | else AliInfo("Solution of Global Dense System by Cholesky failed, trying Gaussian Elimiation"); |
1184 | } | |
1185 | // | |
339fbe23 | 1186 | if (((AliSymMatrix*)fMatCGlo)->SolveSpmInv(fVecBGlo, kTRUE)) return kInvert; |
e30a812f | 1187 | else AliInfo("Solution of Global Dense System by Gaussian Elimination failed, trying iterative methods"); |
8a9ab0eb | 1188 | } |
1189 | // try to solve by minres | |
de34b538 | 1190 | TVectorD sol(fNGloSize); |
8a9ab0eb | 1191 | // |
1192 | AliMinResSolve *slv = new AliMinResSolve(fMatCGlo,fVecBGlo); | |
339fbe23 | 1193 | if (!slv) return kFailed; |
8a9ab0eb | 1194 | // |
1195 | Bool_t res = kFALSE; | |
1196 | if (fgIterSol == AliMinResSolve::kSolMinRes) | |
de34b538 | 1197 | res = slv->SolveMinRes(sol,fgMinResCondType,fgMinResMaxIter,fgMinResTol); |
8a9ab0eb | 1198 | else if (fgIterSol == AliMinResSolve::kSolFGMRes) |
de34b538 | 1199 | res = slv->SolveFGMRES(sol,fgMinResCondType,fgMinResMaxIter,fgMinResTol,fgNKrylovV); |
8a9ab0eb | 1200 | else |
1201 | AliInfo(Form("Undefined Iteritive Solver ID=%d, only %d are defined",fgIterSol,AliMinResSolve::kNSolvers)); | |
1202 | // | |
e30a812f | 1203 | if (!res) { |
1204 | const char* faildump = "fgmr_failed.dat"; | |
1205 | int defout = dup(1); | |
c8f37c50 | 1206 | if (defout<0) { |
1207 | AliInfo("Failed on dup"); | |
339fbe23 | 1208 | return kFailed; |
c8f37c50 | 1209 | } |
642724dc | 1210 | int slvDump = open(faildump, O_RDWR|O_CREAT, 0666); |
f9efbbf6 | 1211 | if (slvDump>=0) { |
1212 | dup2(slvDump,1); | |
1213 | // | |
1214 | printf("#Failed to solve using solver %d with PreCond: %d MaxIter: %d Tol: %e NKrylov: %d\n", | |
1215 | fgIterSol,fgMinResCondType,fgMinResMaxIter,fgMinResTol,fgNKrylovV); | |
1216 | printf("#Dump of matrix:\n"); | |
1217 | fMatCGlo->Print("10"); | |
1218 | printf("#Dump of RHS:\n"); | |
1219 | for (int i=0;i<fNGloSize;i++) printf("%d %+.10f\n",i,fVecBGlo[i]); | |
1220 | // | |
1221 | dup2(defout,1); | |
1222 | close(slvDump); | |
c8396ea8 | 1223 | close(defout); |
f9efbbf6 | 1224 | printf("#Dumped failed matrix and RHS to %s\n",faildump); |
1225 | } | |
1226 | else AliInfo("Failed on file open for matrix dumping"); | |
4c234df8 | 1227 | close(defout); |
339fbe23 | 1228 | return kFailed; |
e30a812f | 1229 | } |
de34b538 | 1230 | for (int i=fNGloSize;i--;) fVecBGlo[i] = sol[i]; |
c9f619c7 | 1231 | // |
339fbe23 | 1232 | return kNoInversion; |
8a9ab0eb | 1233 | // |
1234 | } | |
1235 | ||
1236 | //_____________________________________________________________________________________________ | |
1237 | Float_t AliMillePede2::Chi2DoFLim(int nSig, int nDoF) const | |
1238 | { | |
1239 | /// return the limit in chi^2/nd for n sigmas stdev authorized | |
1240 | // Only n=1, 2, and 3 are expected in input | |
1241 | int lNSig; | |
1242 | float sn[3] = {0.47523, 1.690140, 2.782170}; | |
1243 | float table[3][30] = {{1.0000, 1.1479, 1.1753, 1.1798, 1.1775, 1.1730, 1.1680, 1.1630, | |
1244 | 1.1581, 1.1536, 1.1493, 1.1454, 1.1417, 1.1383, 1.1351, 1.1321, | |
1245 | 1.1293, 1.1266, 1.1242, 1.1218, 1.1196, 1.1175, 1.1155, 1.1136, | |
1246 | 1.1119, 1.1101, 1.1085, 1.1070, 1.1055, 1.1040}, | |
1247 | {4.0000, 3.0900, 2.6750, 2.4290, 2.2628, 2.1415, 2.0481, 1.9736, | |
1248 | 1.9124, 1.8610, 1.8171, 1.7791, 1.7457, 1.7161, 1.6897, 1.6658, | |
1249 | 1.6442, 1.6246, 1.6065, 1.5899, 1.5745, 1.5603, 1.5470, 1.5346, | |
1250 | 1.5230, 1.5120, 1.5017, 1.4920, 1.4829, 1.4742}, | |
1251 | {9.0000, 5.9146, 4.7184, 4.0628, 3.6410, 3.3436, 3.1209, 2.9468, | |
1252 | 2.8063, 2.6902, 2.5922, 2.5082, 2.4352, 2.3711, 2.3143, 2.2635, | |
1253 | 2.2178, 2.1764, 2.1386, 2.1040, 2.0722, 2.0428, 2.0155, 1.9901, | |
1254 | 1.9665, 1.9443, 1.9235, 1.9040, 1.8855, 1.8681}}; | |
1255 | ||
1256 | if (nDoF < 1) { | |
1257 | return 0.0; | |
1258 | } | |
1259 | else { | |
1260 | lNSig = TMath::Max(1,TMath::Min(nSig,3)); | |
1261 | ||
1262 | if (nDoF <= 30) { | |
1263 | return table[lNSig-1][nDoF-1]; | |
1264 | } | |
1265 | else { // approximation | |
1266 | return ((sn[lNSig-1]+TMath::Sqrt(float(2*nDoF-3)))* | |
1267 | (sn[lNSig-1]+TMath::Sqrt(float(2*nDoF-3))))/float(2*nDoF-2); | |
1268 | } | |
1269 | } | |
1270 | } | |
1271 | ||
1272 | //_____________________________________________________________________________________________ | |
1273 | Int_t AliMillePede2::SetIterations(double lChi2CutFac) | |
1274 | { | |
1275 | // Number of iterations is calculated from lChi2CutFac | |
1276 | fChi2CutFactor = TMath::Max(1.0, lChi2CutFac); | |
1277 | // | |
1278 | AliInfo(Form("Initial cut factor is %f",fChi2CutFactor)); | |
1279 | fIter = 1; // Initializes the iteration process | |
1280 | return 1; | |
1281 | } | |
1282 | ||
1283 | //_____________________________________________________________________________________________ | |
1284 | Double_t AliMillePede2::GetParError(int iPar) const | |
1285 | { | |
1286 | // return error for parameter iPar | |
339fbe23 | 1287 | if (fGloSolveStatus==kInvert) { |
9a919f12 | 1288 | if (fkReGroup) iPar = fkReGroup[iPar]; |
f91a16b5 | 1289 | if (iPar<0) { |
1290 | // AliDebug(2,Form("Parameter %d was suppressed in the regrouping",iPar)); | |
1291 | return 0; | |
1292 | } | |
de34b538 | 1293 | double res = fMatCGlo->QueryDiag(iPar); |
8a9ab0eb | 1294 | if (res>=0) return TMath::Sqrt(res); |
1295 | } | |
9a919f12 | 1296 | return 0.; |
1297 | } | |
1298 | ||
1299 | //_____________________________________________________________________________________________ | |
1300 | Double_t AliMillePede2::GetPull(int iPar) const | |
1301 | { | |
1302 | // return pull for parameter iPar | |
1303 | if (fGloSolveStatus==kInvert) { | |
1304 | if (fkReGroup) iPar = fkReGroup[iPar]; | |
f91a16b5 | 1305 | if (iPar<0) { |
1306 | // AliDebug(2,Form("Parameter %d was suppressed in the regrouping",iPar)); | |
1307 | return 0; | |
1308 | } | |
9a919f12 | 1309 | // |
1310 | return fProcPnt[iPar]>0 && (fSigmaPar[iPar]*fSigmaPar[iPar]-fMatCGlo->QueryDiag(iPar))>0. && fSigmaPar[iPar]>0 | |
1311 | ? fDeltaPar[iPar]/TMath::Sqrt(fSigmaPar[iPar]*fSigmaPar[iPar]-fMatCGlo->QueryDiag(iPar)) : 0; | |
1312 | } | |
1313 | return 0.; | |
8a9ab0eb | 1314 | } |
1315 | ||
1316 | ||
1317 | //_____________________________________________________________________________________________ | |
1318 | Int_t AliMillePede2::PrintGlobalParameters() const | |
1319 | { | |
1320 | /// Print the final results into the logfile | |
1321 | double lError = 0.; | |
1322 | double lGlobalCor =0.; | |
1323 | ||
f91a16b5 | 1324 | printf("\nMillePede2 output\n"); |
1325 | printf(" Result of fit for global parameters\n"); | |
1326 | printf(" ===================================\n"); | |
1327 | printf(" I initial final differ lastcor error gcor Npnt\n"); | |
1328 | printf("----------------------------------------------------------------------------------------------\n"); | |
8a9ab0eb | 1329 | // |
f91a16b5 | 1330 | int lastPrintedId = -1; |
9a919f12 | 1331 | for (int i0=0; i0<fNGloParIni; i0++) { |
1332 | int i = GetRGId(i0); if (i<0) continue; | |
f91a16b5 | 1333 | if (i!=i0 && lastPrintedId>=0 && i<=lastPrintedId) continue; // grouped param |
1334 | lastPrintedId = i; | |
1335 | lError = GetParError(i0); | |
8a9ab0eb | 1336 | lGlobalCor = 0.0; |
1337 | // | |
1338 | double dg; | |
339fbe23 | 1339 | if (fGloSolveStatus==kInvert && TMath::Abs( (dg=fMatCGlo->QueryDiag(i)) *fDiagCGlo[i]) > 0) { |
8a9ab0eb | 1340 | lGlobalCor = TMath::Sqrt(TMath::Abs(1.0-1.0/(dg*fDiagCGlo[i]))); |
f91a16b5 | 1341 | printf("%4d(%4d)\t %+.6f\t %+.6f\t %+.6f\t %.6f\t %.6f\t %.6f\t %6d\n", |
1342 | i,i0,fInitPar[i],fInitPar[i]+fDeltaPar[i],fDeltaPar[i],fVecBGlo[i],lError,lGlobalCor,fProcPnt[i]); | |
8a9ab0eb | 1343 | } |
1344 | else { | |
f91a16b5 | 1345 | printf("%4d (%4d)\t %+.6f\t %+.6f\t %+.6f\t %.6f\t OFF\t OFF\t %6d\n",i,i0,fInitPar[i],fInitPar[i]+fDeltaPar[i], |
1346 | fDeltaPar[i],fVecBGlo[i],fProcPnt[i]); | |
8a9ab0eb | 1347 | } |
1348 | } | |
1349 | return 1; | |
1350 | } | |
e30a812f | 1351 | |
1352 | //_____________________________________________________________________________________________ | |
1f6eade8 | 1353 | Bool_t AliMillePede2::IsRecordAcceptable() |
e30a812f | 1354 | { |
1355 | // validate record according run lists set by the user | |
1a39bef4 | 1356 | static Long_t prevRunID = kMaxInt; |
e30a812f | 1357 | static Bool_t prevAns = kTRUE; |
1358 | Long_t runID = fRecord->GetRunID(); | |
1359 | if (runID!=prevRunID) { | |
1360 | int n = 0; | |
1f6eade8 | 1361 | fRunWgh = 1.; |
e30a812f | 1362 | prevRunID = runID; |
1363 | // is run to be rejected? | |
1364 | if (fRejRunList && (n=fRejRunList->GetSize())) { | |
1365 | prevAns = kTRUE; | |
1366 | for (int i=n;i--;) if (runID == (*fRejRunList)[i]) { | |
1367 | prevAns = kFALSE; | |
7fb585bc | 1368 | AliInfo(Form("New Run to reject: %ld",runID)); |
e30a812f | 1369 | break; |
1370 | } | |
1371 | } | |
1372 | else if (fAccRunList && (n=fAccRunList->GetSize())) { // is run specifically selected | |
1373 | prevAns = kFALSE; | |
7fb585bc | 1374 | for (int i=n;i--;) if (runID == (*fAccRunList)[i]) { |
1375 | prevAns = kTRUE; | |
1376 | if (fAccRunListWgh) fRunWgh = (*fAccRunListWgh)[i]; | |
1377 | AliInfo(Form("New Run to accept explicitly: %ld, weight=%f",runID,fRunWgh)); | |
1378 | break; | |
1379 | } | |
1380 | if (!prevAns) AliInfo(Form("New Run is not in the list to accept: %ld",runID)); | |
e30a812f | 1381 | } |
1382 | } | |
1383 | // | |
1384 | return prevAns; | |
1385 | // | |
1386 | } | |
1387 | ||
1388 | //_____________________________________________________________________________________________ | |
1389 | void AliMillePede2::SetRejRunList(const UInt_t *runs, Int_t nruns) | |
1390 | { | |
1391 | // set the list of runs to be rejected | |
1392 | if (fRejRunList) delete fRejRunList; | |
1393 | fRejRunList = 0; | |
1394 | if (nruns<1 || !runs) return; | |
1395 | fRejRunList = new TArrayL(nruns); | |
1396 | for (int i=0;i<nruns;i++) (*fRejRunList)[i] = runs[i]; | |
1397 | } | |
1398 | ||
1399 | //_____________________________________________________________________________________________ | |
1f6eade8 | 1400 | void AliMillePede2::SetAccRunList(const UInt_t *runs, Int_t nruns, const Float_t* wghList) |
e30a812f | 1401 | { |
1402 | // set the list of runs to be selected | |
1403 | if (fAccRunList) delete fAccRunList; | |
1f6eade8 | 1404 | if (fAccRunListWgh) delete fAccRunListWgh; |
e30a812f | 1405 | fAccRunList = 0; |
1406 | if (nruns<1 || !runs) return; | |
1407 | fAccRunList = new TArrayL(nruns); | |
1f6eade8 | 1408 | fAccRunListWgh = new TArrayF(nruns); |
1409 | for (int i=0;i<nruns;i++) { | |
1410 | (*fAccRunList)[i] = runs[i]; | |
1411 | (*fAccRunListWgh)[i] =wghList ? wghList[i] : 1.0; | |
1412 | } | |
e30a812f | 1413 | } |
9a919f12 | 1414 | |
1415 | //_____________________________________________________________________________________________ | |
1416 | void AliMillePede2::SetInitPars(const Double_t* par) | |
1417 | { | |
1418 | // initialize parameters, account for eventual grouping | |
1419 | for (int i=0;i<fNGloParIni;i++) { | |
1420 | int id = GetRGId(i); if (id<0) continue; | |
1421 | fInitPar[id] = par[i]; | |
1422 | } | |
1423 | } | |
1424 | ||
1425 | //_____________________________________________________________________________________________ | |
1426 | void AliMillePede2::SetSigmaPars(const Double_t* par) | |
1427 | { | |
1428 | // initialize sigmas, account for eventual grouping | |
1429 | for (int i=0;i<fNGloParIni;i++) { | |
1430 | int id = GetRGId(i); if (id<0) continue; | |
1431 | fSigmaPar[id] = par[i]; | |
1432 | } | |
1433 | } | |
1434 | ||
1435 | //_____________________________________________________________________________________________ | |
1436 | void AliMillePede2::SetInitPar(Int_t i,Double_t par) | |
1437 | { | |
1438 | // initialize param, account for eventual grouping | |
1439 | int id = GetRGId(i); if (id<0) return; | |
1440 | fInitPar[id] = par; | |
1441 | } | |
1442 | ||
1443 | //_____________________________________________________________________________________________ | |
1444 | void AliMillePede2::SetSigmaPar(Int_t i,Double_t par) | |
1445 | { | |
1446 | // initialize sigma, account for eventual grouping | |
1447 | int id = GetRGId(i); if (id<0) return; | |
1448 | fSigmaPar[id] = par; | |
1449 | } | |
1450 | ||
1451 |