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