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