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