Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 2006-07, ALICE Experiment at CERN, All rights reserved. *
3 : * *
4 : * Author: The ALICE Off-line Project. *
5 : * Contributors are mentioned in the code where appropriate. *
6 : * *
7 : * Permission to use, copy, modify and distribute this software and its *
8 : * documentation strictly for non-commercial purposes is hereby granted *
9 : * without fee, provided that the above copyright notice appears in all *
10 : * copies and that both the copyright notice and this permission notice *
11 : * appear in the supporting documentation. The authors make no claims *
12 : * about the suitability of this software for any purpose. It is *
13 : * provided "as is" without express or implied warranty. *
14 : **************************************************************************/
15 :
16 : //-------------------------------------------------------------------------
17 : // Implementation of the AliNDLocalRegression class
18 : //-------------------------------------------------------------------------
19 :
20 : /*
21 : Related task: https://alice.its.cern.ch/jira/browse/ATO-193
22 :
23 : Algorithm secription - see:
24 : Kernel_smoother: Local polynomial regression
25 : http://en.wikipedia.org/w/index.php?title=Kernel_smoother&oldid=627785784
26 :
27 :
28 : Formally, the local polynomial regression is computed by solving a weighted least square problem.
29 : Weights are provided as a width of the gausian kernel.
30 : Local fit parameters are computed on the grid defined by axis set defiend by THn.
31 : For example use please check UnitTest:
32 : .L $ALICE_ROOT/../src/STAT/test/AliNDLocalRegressionTest.C+
33 : //
34 : Init:
35 : AliNDLocalRegression *pfitNDIdeal=0;
36 : pfitNDIdeal->SetHistogram((THn*)(hN->Clone()));
37 : pfitNDIdeal->SetCuts(3,0.8,1); // outlier rejection setting see /AliNDLocalRegressionTest.C:UnitTestGaussNoisePlusOutliers() for motivation
38 : pfitNDIdeal->MakeFit(treeIn, "val:err", "xyz0:xyz1","Entry$%2==1", "0.05:0.05","2:2",0.001);
39 :
40 : Usage:
41 :
42 : Double_t xyz[2]={1,2}
43 : pfitNDIdeal->Eval(xyz);
44 :
45 : In TFormulas:
46 : pfitNDGaus0->AddVisualCorrection(pfitNDGaus0,2);
47 : pfitNDGaus1->AddVisualCorrection(pfitNDGaus0,3);
48 : treeIn->Draw("(AliNDLocalRegression::GetCorrND(3,xyz0,xyz1)-AliNDLocalRegression::GetCorrND(2,xyz0,xyz1))/sqrt(AliNDLocalRegression::GetCorrNDError(3,xyz0,xyz1)**2+AliNDLocalRegression::GetCorrNDError(2,xyz0,xyz1)**2)>>pullsGaus01(200,-20,20)","","");
49 :
50 :
51 : To do:
52 : 1.) Statistical error of the local interpolation ignores Gaussian kernel weights
53 : errors are overestimated - find a proper mathematical formula to estimate statistical error of estimator
54 : 2.) Implent regularization for smoothing - requesting approximate smoothnes in values and derivative
55 :
56 :
57 : author: marian.ivanov@cern.ch
58 : */
59 :
60 : #include <TVectorD.h>
61 :
62 : #include "AliNDLocalRegression.h"
63 : #include "AliLog.h"
64 :
65 : #include "THn.h"
66 : #include "TObjString.h"
67 : #include "TTreeStream.h"
68 : #include "AliMathBase.h"
69 : #include "TMatrixD.h"
70 : #include "TRobustEstimator.h"
71 : #include "AliMathBase.h"
72 : #include "TStopwatch.h"
73 :
74 128 : ClassImp(AliNDLocalRegression)
75 :
76 : TObjArray *AliNDLocalRegression::fgVisualCorrection=0;
77 : // instance of correction for visualization
78 : Int_t AliNDLocalRegression::fgVerboseLevel=1000;
79 :
80 : AliNDLocalRegression::AliNDLocalRegression():
81 0 : TNamed(),
82 0 : fHistPoints(0), // ND histogram defining regression granularity
83 0 : fRobustFractionLTS(0), // fraction of data used for the robust mean and robust rms estimator (LTS https://en.wikipedia.org/wiki/Least_trimmed_squares)
84 0 : fRobustRMSLTSCut(0), // cut on the robust RMS |value-localmean|<fRobustRMSLTSCut*localRMS
85 0 : fCutType(0), // type of the cut 0- no cut 1-cut localmean=median, 2-cut localmen=rosbut mean
86 0 : fInputTree(0), // input tree - object is not owner
87 0 : fStreamer(0), // optional streamer
88 0 : fFormulaVal(0), // value:err definition formula
89 0 : fSelection(0), // point selector formula
90 0 : fFormulaVar(0), //: separated variable definition formula
91 0 : fKernelWidthFormula(0), //: separated - kernel width for the regression
92 0 : fPolDimensionFormula(0), //: separated - polynom for the regression
93 0 : fNParameters(0), // number of local paramters to fit
94 0 : fLocalFitParam(0), // local fit parameters
95 0 : fLocalFitQuality(0), // local fit quality
96 0 : fLocalFitCovar(0), // local fit covariance matrix
97 0 : fLocalRobustStat(0), // local robust statistic
98 0 : fBinIndex(0), //[fNParameters] working arrays current bin index
99 0 : fBinCenter(0), //[fNParameters] working current local variables - bin center
100 0 : fBinDelta(0), //[fNParameters] working current local variables - bin delta
101 0 : fBinWidth(0), //[fNParameters] working current local variables - bin delta
102 0 : fUseBinNorm(kFALSE) // switch make polynom in units of bins (kTRUE) or in natural units (kFALSE)
103 0 : {
104 0 : if (!fgVisualCorrection) fgVisualCorrection= new TObjArray;
105 0 : }
106 :
107 : AliNDLocalRegression::AliNDLocalRegression(const char* name, const char* title):
108 0 : TNamed(name,title),
109 0 : fHistPoints(0), // ND histogram defining regression granularity
110 0 : fRobustFractionLTS(0), // fraction of data used for the robust mean and robust rms estimator (LTS https://en.wikipedia.org/wiki/Least_trimmed_squares)
111 0 : fRobustRMSLTSCut(0), // cut on the robust RMS |value-localmean|<fRobustRMSLTSCut*localRMS
112 0 : fCutType(0), // type of the cut 0- no cut 1-cut localmean=median, 2-cut localmen=rosbut mean
113 0 : fInputTree(0), // input tree - object is not owner
114 0 : fStreamer(0), // optional streamer
115 0 : fFormulaVal(0), // value:err definition formula
116 0 : fSelection(0), // point selector formula
117 0 : fFormulaVar(0), //: separated variable definition formula
118 0 : fKernelWidthFormula(0), //: separated - kernel width for the regression
119 0 : fPolDimensionFormula(0), //: separated - polynom for the regression
120 0 : fNParameters(0), // number of local paramters to fit
121 0 : fLocalFitParam(0), // local fit parameters
122 0 : fLocalFitQuality(0), // local fit quality
123 0 : fLocalFitCovar(0), // local fit covariance matrix
124 0 : fLocalRobustStat(0), // local robust statistic
125 0 : fBinIndex(0), //[fNParameters] working arrays current bin index
126 0 : fBinCenter(0), //[fNParameters] working current local variables - bin center
127 0 : fBinDelta(0), //[fNParameters] working current local variables - bin delta
128 0 : fBinWidth(0), //[fNParameters] working current local variables - bin delta
129 0 : fUseBinNorm(kFALSE) // switch make polynom in units of bins (kTRUE) or in natural units (kFALSE)
130 0 : {
131 0 : }
132 :
133 0 : AliNDLocalRegression::~AliNDLocalRegression(){
134 : //
135 : // destructor
136 : //
137 0 : if (fHistPoints) delete fHistPoints;
138 0 : if (fStreamer) delete fStreamer;
139 : // fInputTree(0), //! input tree - object is not owner
140 0 : delete fFormulaVal; // value:err definition formula
141 0 : delete fSelection; // point selector formula
142 0 : delete fFormulaVar; //: separated variable definition formula
143 0 : delete fKernelWidthFormula; //: separated - kernel width for the regression
144 0 : delete fPolDimensionFormula; //: separated - polynom for the regression
145 : //
146 0 : if (fLocalFitParam) fLocalFitParam->Delete();
147 0 : if (fLocalFitQuality) fLocalFitQuality->Delete();
148 0 : if (fLocalFitCovar) fLocalFitCovar->Delete();
149 0 : delete fLocalFitParam; // local fit parameters
150 0 : delete fLocalFitQuality; // local fit quality
151 0 : delete fLocalFitCovar; // local fit covariance matrix
152 : //
153 0 : delete fLocalRobustStat;
154 : //
155 0 : delete[] fBinIndex;
156 0 : delete[] fBinCenter;
157 0 : delete[] fBinDelta;
158 0 : }
159 :
160 : void AliNDLocalRegression::SetHistogram(THn* histo ){
161 : //
162 : // Setup the local regression ayout according THn hitogram binning
163 : //
164 0 : if (fHistPoints!=0){
165 0 : AliError("Hostogram initialized\n");
166 0 : return ;
167 : }
168 0 : fHistPoints=histo;
169 0 : fLocalFitParam = new TObjArray(fHistPoints->GetNbins());
170 0 : fLocalFitParam->SetOwner(kTRUE);
171 0 : fLocalFitQuality = new TObjArray(fHistPoints->GetNbins());
172 0 : fLocalFitQuality->SetOwner(kTRUE);
173 0 : fLocalFitCovar = new TObjArray(fHistPoints->GetNbins());
174 0 : fLocalFitCovar->SetOwner(kTRUE);
175 :
176 : //
177 : // Check histogram
178 : //
179 0 : Int_t ndim = histo->GetNdimensions();
180 : Bool_t isOK=kTRUE;
181 0 : for (Int_t idim=0; idim<ndim; idim++){
182 0 : TAxis * axis = histo->GetAxis(idim);
183 0 : if (axis->GetNbins()<2) {
184 0 : AliError(TString::Format("Invalid binning nbins<2 %d", axis->GetNbins()).Data());
185 0 : }
186 0 : if (axis->GetXmin()>=axis->GetXmax()) {
187 0 : AliError(TString::Format("Invalid range <%f,%f", axis->GetXmin(),axis->GetXmax()).Data());
188 0 : }
189 : }
190 :
191 :
192 0 : }
193 : void AliNDLocalRegression::SetCuts(Double_t nSigma, Double_t robustFraction, Int_t estimator){
194 : //
195 : //
196 : //
197 0 : fRobustFractionLTS=robustFraction; // fraction of data used for the robust mean and robust rms estimator (LTS https://en.wikipedia.org/wiki/Least_trimmed_squares)
198 0 : fRobustRMSLTSCut=nSigma; // cut on the robust RMS |value-localmean|<fRobustRMSLTSCut*localRMS
199 0 : fCutType=estimator; // type of the cut 0- no cut 1-cut localmean=median, 2-cut localmen=rosbut mean
200 :
201 0 : }
202 :
203 : Bool_t AliNDLocalRegression::CleanCovariance(){
204 : //
205 : // Clean covariance matrix if not needed anymore
206 : //
207 0 : if (fLocalFitCovar) delete fLocalFitCovar;
208 0 : fLocalFitCovar=0;
209 0 : };
210 :
211 :
212 : Bool_t AliNDLocalRegression::MakeFit(TTree * tree , const char* formulaVal, const char * formulaVar, const char*selection, const char * formulaKernel, const char * dimensionFormula, Double_t weightCut, Int_t entries, Bool_t useBinNorm){
213 : //
214 : // Make a local fit in grid as specified by the input THn histogram
215 : // Histogram has to be set before invocation of method
216 : //
217 : // Output:
218 : // array of fit parameters and covariance matrices
219 : //
220 : // Input Parameters:
221 : // tree - input tree
222 : // formulaVal - : separated variable:error string
223 : // formulaVar - : separate varaible list
224 : // selection - selection (cut) for TTreeDraw
225 : // kernelWidth - : separated list of width of kernel for local fitting
226 : // dimenstionFormula - dummy for the moment
227 : //
228 : //Algorithm:
229 : // 1.) Check consistency of input data
230 : //
231 : // 2.) Cache input data from tree to the array of vector TVectorD
232 : //
233 : // 3.) Calculate robust local mean and robust local RMS in case outlier removal algorithm specified
234 : //
235 : // 4.) Make local fit
236 : //
237 : // const Double_t kEpsilon=1e-6;
238 : const Int_t kMaxDim=100;
239 0 : Int_t binRange[kMaxDim]={0};
240 0 : Double_t nchi2Cut=-2*TMath::Log(weightCut); // transform probability to nsigma cut
241 0 : if (fHistPoints==NULL){
242 0 : AliError("ND histogram not initialized\n");
243 0 : return kFALSE;
244 : }
245 0 : if (tree==NULL || tree->GetEntries()==0){
246 0 : AliError("Empty tree\n");
247 0 : return kFALSE;
248 : }
249 0 : if (formulaVar==NULL || formulaVar==0) {
250 0 : AliError("Empty variable list\n");
251 0 : return kFALSE;
252 : }
253 0 : if (formulaKernel==NULL) {
254 0 : AliError("Kernel width not specified\n");
255 0 : return kFALSE;
256 : }
257 0 : fUseBinNorm=useBinNorm;
258 : //
259 0 : fInputTree= tree; // should be better TRef?
260 0 : fFormulaVal = new TObjString(formulaVal);
261 0 : fFormulaVar = new TObjString(formulaVar);
262 0 : fSelection = new TObjString(selection);
263 0 : fKernelWidthFormula = new TObjString(formulaKernel);
264 0 : fPolDimensionFormula = new TObjString(dimensionFormula);
265 0 : TObjArray * arrayFormulaVar=fFormulaVar->String().Tokenize(":");
266 0 : Int_t nvarFormula = arrayFormulaVar->GetEntries();
267 0 : if (nvarFormula!=fHistPoints->GetNdimensions()){
268 0 : AliError("Histogram/points mismatch\n");
269 0 : return kFALSE;
270 : }
271 0 : TObjArray * arrayKernel=fKernelWidthFormula->String().Tokenize(":");
272 0 : Int_t nwidthFormula = arrayKernel->GetEntries();
273 0 : if (nvarFormula!=nwidthFormula){
274 0 : delete arrayKernel;
275 0 : delete arrayFormulaVar;
276 0 : AliError("Variable/Kernel mismath\n");
277 0 : return kFALSE;
278 : }
279 0 : fNParameters=nvarFormula;
280 : //
281 : // 2.) Load input data
282 : //
283 : //
284 0 : Int_t entriesVal = tree->Draw(formulaVal,selection,"goffpara",entries);
285 0 : if (entriesVal==0) {
286 0 : AliError(TString::Format("Empty point list\t%s\t%s\n",formulaVal,selection).Data());
287 0 : return kFALSE;
288 : }
289 0 : if (tree->GetVal(0)==NULL || (tree->GetVal(1)==NULL)){
290 0 : AliError(TString::Format("Wrong selection\t%s\t%s\n",formulaVar,selection).Data());
291 0 : return kFALSE;
292 : }
293 0 : TVectorD values(entriesVal,tree->GetVal(0));
294 0 : TVectorD errors(entriesVal,tree->GetVal(1));
295 0 : Double_t *pvalues = values.GetMatrixArray();
296 0 : Double_t *perrors = errors.GetMatrixArray();
297 :
298 :
299 : // 2.b) variables
300 0 : TObjArray pointArray(fNParameters);
301 0 : Int_t entriesVar = tree->Draw(formulaVar,selection,"goffpara",entries);
302 0 : if (entriesVal!=entriesVar) {
303 0 : AliError(TString::Format("Wrong selection\t%s\t%s\n",formulaVar,selection).Data());
304 0 : return kFALSE;
305 : }
306 0 : for (Int_t ipar=0; ipar<fNParameters; ipar++) pointArray.AddAt(new TVectorD(entriesVar,tree->GetVal(ipar)),ipar);
307 : // 2.c) kernel array
308 0 : TObjArray kernelArrayI2(fNParameters);
309 0 : Int_t entriesKernel = tree->Draw(formulaKernel,selection,"goffpara",entries);
310 0 : for (Int_t ipar=0; ipar<fNParameters; ipar++) {
311 0 : TVectorD* vdI2 = new TVectorD(entriesVar,tree->GetVal(ipar));
312 0 : for (int k=entriesVar;k--;) { // to speed up, precalculate inverse squared
313 0 : double kv = (*vdI2)[k];
314 0 : if (TMath::Abs(kv)<1e-12) AliFatalClassF("Kernel width=%f for entry %d of par:%d",kv,k,ipar);
315 0 : (*vdI2)[k] = 1./(kv*kv);
316 : }
317 0 : kernelArrayI2.AddAt(vdI2,ipar);
318 : }
319 : //
320 0 : Double_t *pvecVar[kMaxDim]={0};
321 0 : Double_t *pvecKernelI2[kMaxDim]={0};
322 0 : for (Int_t idim=0; idim<pointArray.GetEntries(); idim++){
323 0 : pvecVar[idim]=((TVectorD*)(pointArray.At(idim)))->GetMatrixArray();
324 0 : pvecKernelI2[idim]=((TVectorD*)(kernelArrayI2.At(idim)))->GetMatrixArray();
325 0 : binRange[idim]=fHistPoints->GetAxis(idim)->GetNbins();
326 : }
327 : //
328 : //
329 : //
330 0 : Int_t nbins = fHistPoints->GetNbins();
331 0 : fBinIndex = new Int_t[fHistPoints->GetNdimensions()];
332 0 : fBinCenter = new Double_t[fHistPoints->GetNdimensions()];
333 0 : fBinDelta = new Double_t[fHistPoints->GetNdimensions()];
334 0 : fBinWidth = new Double_t[fHistPoints->GetNdimensions()];
335 :
336 : //
337 : // 3.)
338 : //
339 0 : if (fCutType>0 && fRobustRMSLTSCut>0){
340 0 : MakeRobustStatistic(values, errors, pointArray, kernelArrayI2, weightCut, fRobustFractionLTS);
341 : }
342 : //
343 :
344 : // 4.) Make local fits
345 : //
346 :
347 0 : Double_t *binHypFit = new Double_t[2*fHistPoints->GetNdimensions()];
348 : //
349 0 : TLinearFitter fitter(1+2*fNParameters,TString::Format("hyp%d",2*fNParameters).Data());
350 0 : for (Int_t ibin=0; ibin<nbins; ibin++) {
351 0 : fHistPoints->GetBinContent(ibin,fBinIndex);
352 : Bool_t isUnderFlowBin=kFALSE;
353 : Bool_t isOverFlowBin=kFALSE;
354 0 : for (Int_t idim=0; idim<fNParameters; idim++) {
355 0 : if (fBinIndex[idim]==0) isUnderFlowBin=kTRUE;
356 0 : if (fBinIndex[idim]>binRange[idim]) isOverFlowBin=kTRUE;
357 0 : fBinCenter[idim]=fHistPoints->GetAxis(idim)->GetBinCenter(fBinIndex[idim]);
358 0 : fBinWidth[idim]=fHistPoints->GetAxis(idim)->GetBinWidth(fBinIndex[idim]);
359 : }
360 0 : if (isUnderFlowBin || isOverFlowBin) continue;
361 0 : fitter.ClearPoints();
362 : // add fit points
363 0 : for (Int_t ipoint=0; ipoint<entriesVal; ipoint++){
364 : Double_t sumChi2=0;
365 0 : if (fCutType>0 && fRobustRMSLTSCut>0){
366 0 : Double_t localRMS=(*fLocalRobustStat)(ibin,2);
367 0 : Double_t localMean=(*fLocalRobustStat)(ibin,1);
368 0 : Double_t localMedian=(*fLocalRobustStat)(ibin,0);
369 0 : if (fCutType==1){
370 0 : if (TMath::Abs(pvalues[ipoint]-localMedian)>fRobustRMSLTSCut*localRMS) continue;
371 : }
372 0 : if (fCutType==2){
373 0 : if (TMath::Abs(pvalues[ipoint]-localMean)>fRobustRMSLTSCut*localRMS) continue;
374 : }
375 0 : }
376 0 : for (Int_t idim=0; idim<fNParameters; idim++){
377 : //TVectorD &vecVar=*((TVectorD*)(pointArray.UncheckedAt(idim)));
378 : //TVectorD &vecKernel=*((TVectorD*)(kernelArray.UncheckedAt(idim)));
379 0 : fBinDelta[idim]=pvecVar[idim][ipoint]-fBinCenter[idim];
380 0 : sumChi2+= (fBinDelta[idim]*fBinDelta[idim]) * pvecKernelI2[idim][ipoint];
381 0 : if (sumChi2>nchi2Cut) break;//continue;
382 0 : if (fUseBinNorm){
383 0 : binHypFit[2*idim]=fBinDelta[idim]/fBinWidth[idim];
384 0 : binHypFit[2*idim+1]=binHypFit[2*idim]*binHypFit[2*idim];
385 0 : }else{
386 0 : binHypFit[2*idim]=fBinDelta[idim];
387 0 : binHypFit[2*idim+1]=fBinDelta[idim]*fBinDelta[idim];
388 : }
389 : }
390 0 : if (sumChi2>nchi2Cut) continue;
391 : // Double_t weight=TMath::Exp(-sumChi2*0.5);
392 : // fitter.AddPoint(binHypFit,pvalues[ipoint], perrors[ipoint]/weight);
393 0 : Double_t weightI=TMath::Exp(sumChi2*0.5);
394 0 : fitter.AddPoint(binHypFit,pvalues[ipoint], perrors[ipoint]*weightI);
395 0 : }
396 0 : TVectorD * fitParam=new TVectorD(fNParameters*2+1);
397 0 : TVectorD * fitQuality=new TVectorD(3);
398 0 : TMatrixD * fitCovar=new TMatrixD(fNParameters*2+1,fNParameters*2+1);
399 0 : Double_t normRMS=0;
400 0 : Int_t nBinPoints=fitter.GetNpoints();
401 0 : Bool_t fitOK=kFALSE;
402 0 : (*fitQuality)[0]=0;
403 0 : (*fitQuality)[1]=0;
404 0 : (*fitQuality)[2]=0;
405 :
406 0 : if (fitter.GetNpoints()>fNParameters*2+2){
407 0 : fitOK = (fitter.Eval()==0);
408 0 : if (fitOK){
409 0 : normRMS=fitter.GetChisquare()/(fitter.GetNpoints()-fitter.GetNumberFreeParameters());
410 0 : fitter.GetParameters(*fitParam);
411 0 : fitter.GetCovarianceMatrix(*fitCovar);
412 0 : (*fitQuality)[0]=nBinPoints;
413 0 : (*fitQuality)[1]=normRMS;
414 0 : (*fitQuality)[2]=ibin;
415 0 : fLocalFitParam->AddAt(fitParam,ibin);
416 0 : fLocalFitQuality->AddAt(fitQuality,ibin);
417 0 : fLocalFitCovar->AddAt(fitCovar,ibin);
418 : }
419 : }
420 0 : if (fStreamer){
421 0 : TVectorD pfBinCenter(fNParameters, fBinCenter);
422 0 : Double_t median=0,mean=0,rms=0;
423 0 : if (fLocalRobustStat){
424 0 : median=(*fLocalRobustStat)(ibin,0);
425 0 : mean=(*fLocalRobustStat)(ibin,1);
426 0 : rms=(*fLocalRobustStat)(ibin,2);
427 0 : }
428 0 : (*fStreamer)<<"localFit"<<
429 0 : "ibin="<<ibin<< // bin index
430 0 : "fitOK="<<fitOK<<
431 0 : "localMedian="<<median<<
432 0 : "localMean="<<mean<<
433 0 : "localRMS="<<rms<<
434 0 : "nBinPoints="<<nBinPoints<< // center of the bin
435 0 : "binCenter.="<<&pfBinCenter<< //
436 0 : "normRMS="<<normRMS<<
437 0 : "fitParam.="<<fitParam<<
438 0 : "fitCovar.="<<fitCovar<<
439 0 : "fitOK="<<fitOK<<
440 : "\n";
441 0 : }
442 0 : if (!fitOK) { // avoid memory leak for failed fits
443 0 : delete fitParam;
444 0 : delete fitQuality;
445 0 : delete fitCovar;
446 : }
447 0 : }
448 :
449 : return kTRUE;
450 0 : }
451 :
452 :
453 : Bool_t AliNDLocalRegression::MakeRobustStatistic(TVectorD &values,TVectorD &errors, TObjArray &pointArray, TObjArray &kernelArrayI2, Double_t weightCut, Double_t robustFraction){
454 : //
455 : // Calculate robust statistic information
456 : // use raw array to make faster calcualtion
457 : const Int_t kMaxDim=100;
458 0 : Double_t *pvalues=values.GetMatrixArray();
459 0 : Double_t *pvecVar[kMaxDim]={0};
460 0 : Double_t *pvecKernelI2[kMaxDim]={0};
461 0 : for (Int_t idim=0; idim<pointArray.GetEntries(); idim++){
462 0 : pvecVar[idim]=((TVectorD*)(pointArray.At(idim)))->GetMatrixArray();
463 0 : pvecKernelI2[idim]=((TVectorD*)(kernelArrayI2.At(idim)))->GetMatrixArray();
464 :
465 : }
466 :
467 :
468 0 : Double_t nchi2Cut=-2*TMath::Log(weightCut); // transform probability to nsigma cut
469 0 : if (robustFraction>1) robustFraction=1;
470 :
471 0 : Int_t nbins = fHistPoints->GetNbins(); //
472 0 : Int_t npoints= values.GetNrows(); // number of points for fit
473 0 : if (fLocalRobustStat){
474 0 : delete fLocalRobustStat;
475 : }
476 0 : fLocalRobustStat=new TMatrixD(nbins,3);
477 :
478 0 : TVectorD valueLocal(npoints);
479 0 : for (Int_t ibin=0; ibin<nbins; ibin++){
480 0 : fHistPoints->GetBinContent(ibin,fBinIndex); //
481 0 : for (Int_t idim=0; idim<fNParameters; idim++){
482 0 : fBinCenter[idim]=fHistPoints->GetAxis(idim)->GetBinCenter(fBinIndex[idim]);
483 0 : fBinWidth[idim]=fHistPoints->GetAxis(idim)->GetBinWidth(fBinIndex[idim]);
484 : }
485 : Int_t indexLocal=0;
486 0 : for (Int_t ipoint=0; ipoint<npoints; ipoint++){
487 : Double_t sumChi2=0;
488 0 : for (Int_t idim=0; idim<fNParameters; idim++){
489 : //TVectorD &vecVar=*((TVectorD*)(pointArray.UncheckedAt(idim)));
490 : //TVectorD &vecKernel=*((TVectorD*)(kernelArray.UncheckedAt(idim)));
491 0 : fBinDelta[idim]=pvecVar[idim][ipoint]-fBinCenter[idim];
492 0 : sumChi2+= (fBinDelta[idim]*fBinDelta[idim]) * pvecKernelI2[idim][ipoint];
493 0 : if (sumChi2>nchi2Cut) break; //continue;
494 : }
495 0 : if (sumChi2>nchi2Cut) continue;
496 0 : valueLocal[indexLocal]=pvalues[ipoint];
497 0 : indexLocal++;
498 0 : }
499 0 : Double_t median=0,meanX=0, rmsX=0;
500 0 : if (indexLocal*robustFraction-1>3){
501 0 : median=TMath::Median(indexLocal,valueLocal.GetMatrixArray());
502 0 : AliMathBase::EvaluateUni(indexLocal,valueLocal.GetMatrixArray(), meanX,rmsX, indexLocal*robustFraction-1);
503 : }
504 0 : (*fLocalRobustStat)(ibin,0)=median;
505 0 : (*fLocalRobustStat)(ibin,1)=meanX;
506 0 : (*fLocalRobustStat)(ibin,2)=rmsX;
507 0 : }
508 0 : }
509 :
510 :
511 :
512 : Double_t AliNDLocalRegression::Eval(Double_t *point ){
513 : //
514 : //
515 : //
516 : const Double_t almost0=0.00000001;
517 : // backward compatibility
518 0 : if(!fBinWidth){
519 0 : fBinWidth = new Double_t[fHistPoints->GetNdimensions()];
520 0 : }
521 :
522 0 : for (Int_t iDim=0; iDim<fNParameters; iDim++){
523 0 : if (point[iDim]<= fHistPoints->GetAxis(iDim)->GetXmin()) point[iDim]=fHistPoints->GetAxis(iDim)->GetXmin()+almost0*fHistPoints->GetAxis(iDim)->GetBinWidth(0);
524 0 : if (point[iDim]>= fHistPoints->GetAxis(iDim)->GetXmax()) point[iDim]=fHistPoints->GetAxis(iDim)->GetXmax()-almost0*fHistPoints->GetAxis(iDim)->GetBinWidth(0);
525 : }
526 :
527 0 : Int_t ibin = fHistPoints->GetBin(point);
528 : Bool_t rangeOK=kTRUE;
529 0 : if (ibin>=fLocalFitParam->GetEntriesFast() ){
530 : rangeOK=kFALSE;
531 0 : }else{
532 0 : if (fLocalFitParam->UncheckedAt(ibin)==NULL) {
533 : rangeOK=kFALSE;
534 0 : }
535 : }
536 0 : if (!rangeOK) return 0;
537 :
538 0 : fHistPoints->GetBinContent(ibin,fBinIndex);
539 0 : for (Int_t idim=0; idim<fNParameters; idim++){
540 0 : fBinCenter[idim]=fHistPoints->GetAxis(idim)->GetBinCenter(fBinIndex[idim]);
541 0 : fBinWidth[idim]=fHistPoints->GetAxis(idim)->GetBinWidth(fBinIndex[idim]);
542 : }
543 0 : TVectorD &vecParam = *((TVectorD*)fLocalFitParam->At(ibin));
544 0 : Double_t value=vecParam[0];
545 0 : if (!rangeOK) return value;
546 0 : if (fUseBinNorm) {
547 0 : for (Int_t ipar=0; ipar<fNParameters; ipar++){
548 0 : Double_t delta=(point[ipar]-fBinCenter[ipar])/fBinWidth[ipar];
549 0 : value+=(vecParam[1+2*ipar]+vecParam[1+2*ipar+1]*delta)*delta;
550 : }
551 0 : } else {
552 0 : for (Int_t ipar=0; ipar<fNParameters; ipar++){
553 0 : Double_t delta=(point[ipar]-fBinCenter[ipar]);
554 0 : value+=(vecParam[1+2*ipar]+vecParam[1+2*ipar+1]*delta)*delta;
555 : }
556 : }
557 0 : return value;
558 0 : }
559 :
560 : Double_t AliNDLocalRegression::EvalError(Double_t *point ){
561 : //
562 : //
563 : //
564 0 : if (fLocalFitCovar==NULL) {
565 0 : ::Error("AliNDLocalRegression::EvalError", "Covariance matrix not available");
566 0 : return 0 ;
567 : }
568 0 : for (Int_t iDim=0; iDim<fNParameters; iDim++){
569 0 : if (point[iDim]< fHistPoints->GetAxis(iDim)->GetXmin()) point[iDim]=fHistPoints->GetAxis(iDim)->GetXmin();
570 0 : if (point[iDim]> fHistPoints->GetAxis(iDim)->GetXmax()) point[iDim]=fHistPoints->GetAxis(iDim)->GetXmax();
571 : }
572 :
573 0 : Int_t ibin = fHistPoints->GetBin(point);
574 0 : if (fLocalFitParam->At(ibin)==NULL) return 0;
575 0 : fHistPoints->GetBinContent(ibin,fBinIndex);
576 0 : for (Int_t idim=0; idim<fNParameters; idim++){
577 0 : fBinCenter[idim]=fHistPoints->GetAxis(idim)->GetBinCenter(fBinIndex[idim]);
578 : }
579 0 : TMatrixD &vecCovar = *((TMatrixD*)fLocalFitCovar->At(ibin));
580 : //TVectorD &vecQuality = *((TVectorD*)fLocalFitQuality->At(ibin));
581 0 : Double_t value=TMath::Sqrt(vecCovar(0,0)); // fill covariance to be used
582 : return value;
583 0 : }
584 :
585 : Bool_t AliNDLocalRegression::Derivative(Double_t *point, Double_t *d)
586 : {
587 : // fill d by partial derivatives
588 : //
589 : const Double_t almost0=0.00000001;
590 0 : for (Int_t iDim=0; iDim<fNParameters; iDim++){
591 0 : const TAxis* ax = fHistPoints->GetAxis(iDim);
592 0 : if (point[iDim]<=ax->GetXmin()) point[iDim] = ax->GetXmin()+almost0*ax->GetBinWidth(0);
593 0 : else if (point[iDim]>=ax->GetXmax()) point[iDim] = ax->GetXmax()-almost0*ax->GetBinWidth(0);
594 : }
595 :
596 0 : Int_t ibin = fHistPoints->GetBin(point);
597 0 : if (ibin>=fLocalFitParam->GetEntriesFast() ||
598 0 : !fLocalFitParam->UncheckedAt(ibin) ) return kFALSE;
599 : //
600 0 : fHistPoints->GetBinContent(ibin,fBinIndex);
601 0 : for (Int_t idim=0; idim<fNParameters; idim++){
602 0 : const TAxis* ax = fHistPoints->GetAxis(idim);
603 0 : fBinCenter[idim] = ax->GetBinCenter(fBinIndex[idim]);
604 0 : fBinWidth[idim] = ax->GetBinWidth(fBinIndex[idim]);
605 : }
606 0 : TVectorD &vecParam = *((TVectorD*)fLocalFitParam->At(ibin));
607 :
608 0 : if (fUseBinNorm){
609 0 : for (Int_t ipar=0; ipar<fNParameters; ipar++){
610 0 : Double_t delta=(point[ipar]-fBinCenter[ipar])/fBinWidth[ipar];
611 0 : d[ipar] = vecParam[1+2*ipar] + 2.0*vecParam[1+2*ipar+1]*delta;
612 : }
613 0 : }else{
614 0 : for (Int_t ipar=0; ipar<fNParameters; ipar++){
615 0 : Double_t delta=(point[ipar]-fBinCenter[ipar]);
616 0 : d[ipar] = vecParam[1+2*ipar] + 2.0*vecParam[1+2*ipar+1]*delta;
617 : }
618 : }
619 : return kTRUE;
620 0 : }
621 :
622 : Bool_t AliNDLocalRegression::EvalAndDerivative(Double_t *point, Double_t &val, Double_t *d)
623 : {
624 : // fill d by partial derivatives and calculate value in one go
625 : //
626 : const Double_t almost0=0.00000001;
627 0 : for (Int_t iDim=0; iDim<fNParameters; iDim++){
628 0 : const TAxis* ax = fHistPoints->GetAxis(iDim);
629 0 : if (point[iDim]<=ax->GetXmin()) point[iDim] = ax->GetXmin()+almost0*ax->GetBinWidth(0);
630 0 : else if (point[iDim]>=ax->GetXmax()) point[iDim] = ax->GetXmax()-almost0*ax->GetBinWidth(0);
631 : }
632 0 : val = 0;
633 0 : Int_t ibin = fHistPoints->GetBin(point);
634 0 : if (ibin>=fLocalFitParam->GetEntriesFast() ||
635 0 : !fLocalFitParam->UncheckedAt(ibin) ) return kFALSE;
636 : //
637 0 : fHistPoints->GetBinContent(ibin,fBinIndex);
638 0 : for (Int_t idim=0; idim<fNParameters; idim++){
639 0 : const TAxis* ax = fHistPoints->GetAxis(idim);
640 0 : fBinCenter[idim] = ax->GetBinCenter(fBinIndex[idim]);
641 0 : fBinWidth[idim] = ax->GetBinWidth(fBinIndex[idim]);
642 : }
643 0 : TVectorD &vecParam = *((TVectorD*)fLocalFitParam->At(ibin));
644 :
645 0 : val = vecParam[0];
646 0 : if (fUseBinNorm){
647 0 : for (Int_t ipar=0; ipar<fNParameters; ipar++){
648 0 : Double_t delta=(point[ipar]-fBinCenter[ipar])/fBinWidth[ipar];
649 0 : d[ipar] = vecParam[1+2*ipar] + 2.0*vecParam[1+2*ipar+1]*delta;
650 0 : val += (vecParam[1+2*ipar]+vecParam[1+2*ipar+1]*delta)*delta;
651 : }
652 0 : } else {
653 0 : for (Int_t ipar=0; ipar<fNParameters; ipar++){
654 0 : Double_t delta=(point[ipar]-fBinCenter[ipar]);
655 0 : d[ipar] = vecParam[1+2*ipar] + 2.0*vecParam[1+2*ipar+1]*delta;
656 0 : val += (vecParam[1+2*ipar]+vecParam[1+2*ipar+1]*delta)*delta;
657 : }
658 : }
659 : return kTRUE;
660 0 : }
661 :
662 :
663 : Int_t AliNDLocalRegression::GetVisualCorrectionIndex(const char *corName){
664 : //
665 0 : return TMath::Hash(corName)%1000000;
666 : }
667 :
668 :
669 : void AliNDLocalRegression::AddVisualCorrection(AliNDLocalRegression* corr, Int_t position){
670 : /// make correction available for visualization using
671 : /// TFormula, TFX and TTree::Draw
672 : /// important in order to check corrections and also compute dervied variables
673 : /// e.g correction partial derivatives
674 : ///
675 : /// NOTE - class is not owner of correction
676 0 : if (position==0) {
677 0 : position=GetVisualCorrectionIndex(corr->GetName());
678 0 : }
679 :
680 0 : if (!fgVisualCorrection) fgVisualCorrection=new TObjArray(1000000);
681 0 : if (position>=fgVisualCorrection->GetEntriesFast())
682 0 : fgVisualCorrection->Expand((position+10)*2);
683 0 : if (fgVisualCorrection->At(position)!=NULL){
684 0 : ::Error("AliNDLocalRegression::AddVisualCorrection","Correction %d already defined Old: %s New: %s",position,fgVisualCorrection->At(position)->GetName(), corr->GetName());
685 0 : }
686 0 : fgVisualCorrection->AddAt(corr, position);
687 0 : }
688 :
689 : AliNDLocalRegression* AliNDLocalRegression::GetVisualCorrection(Int_t position) {
690 : /// Get visula correction registered at index=position
691 0 : return fgVisualCorrection? (AliNDLocalRegression*)fgVisualCorrection->At(position):0;
692 : }
693 :
694 : Double_t AliNDLocalRegression::GetCorrND(Double_t index, Double_t par0){
695 : //
696 : //
697 0 : AliNDLocalRegression *corr = (AliNDLocalRegression*)fgVisualCorrection->At(index);
698 0 : if (!corr) return 0;
699 0 : return corr->Eval(&par0);
700 0 : }
701 :
702 : Double_t AliNDLocalRegression::GetCorrNDError(Double_t index, Double_t par0){
703 : //
704 : //
705 0 : AliNDLocalRegression *corr = (AliNDLocalRegression*)fgVisualCorrection->At(index);
706 0 : if (!corr) return 0;
707 0 : return corr->EvalError(&par0);
708 0 : }
709 :
710 : Double_t AliNDLocalRegression::GetCorrND(Double_t index, Double_t par0, Double_t par1){
711 : //
712 : //
713 0 : AliNDLocalRegression *corr = (AliNDLocalRegression*)fgVisualCorrection->At(index);
714 0 : if (!corr) return 0;
715 0 : Double_t par[2]={par0,par1};
716 0 : return corr->Eval(par);
717 0 : }
718 : Double_t AliNDLocalRegression::GetCorrNDError(Double_t index, Double_t par0, Double_t par1){
719 : //
720 : //
721 0 : AliNDLocalRegression *corr = (AliNDLocalRegression*)fgVisualCorrection->At(index);
722 0 : if (!corr) return 0;
723 0 : Double_t par[2]={par0,par1};
724 0 : return corr->EvalError(par);
725 0 : }
726 :
727 : Double_t AliNDLocalRegression::GetCorrND(Double_t index, Double_t par0, Double_t par1, Double_t par2){
728 : //
729 : //
730 0 : AliNDLocalRegression *corr = (AliNDLocalRegression*)fgVisualCorrection->At(index);
731 0 : if (!corr) return 0;
732 0 : Double_t par[3]={par0,par1,par2};
733 0 : return corr->Eval(par);
734 0 : }
735 :
736 : Double_t AliNDLocalRegression::GetCorrNDError(Double_t index, Double_t par0, Double_t par1, Double_t par2){
737 : //
738 : //
739 0 : AliNDLocalRegression *corr = (AliNDLocalRegression*)fgVisualCorrection->At(index);
740 0 : if (!corr) return 0;
741 0 : Double_t par[3]={par0,par1,par2};
742 0 : return corr->EvalError(par);
743 0 : }
744 :
745 :
746 :
747 : Double_t AliNDLocalRegression::GetCorrND(Double_t index, Double_t par0, Double_t par1, Double_t par2, Double_t par3){
748 : //
749 : //
750 0 : AliNDLocalRegression *corr = (AliNDLocalRegression*)fgVisualCorrection->At(index);
751 0 : if (!corr) return 0;
752 0 : Double_t par[4]={par0,par1,par2,par3};
753 0 : return corr->Eval(par);
754 0 : }
755 :
756 : Double_t AliNDLocalRegression::GetCorrNDError(Double_t index, Double_t par0, Double_t par1, Double_t par2, Double_t par3){
757 : //
758 : //
759 0 : AliNDLocalRegression *corr = (AliNDLocalRegression*)fgVisualCorrection->At(index);
760 0 : if (!corr) return 0;
761 0 : Double_t par[4]={par0,par1,par2,par3};
762 0 : return corr->EvalError(par);
763 0 : }
764 :
765 :
766 :
767 :
768 : Bool_t AliNDLocalRegression::AddWeekConstrainsAtBoundaries(Int_t nDims, Int_t *indexes, Double_t *relWeight, TTreeSRedirector* pcstream, Bool_t useCommon){
769 : //
770 : // Adding week constrain AtBoundaries
771 : //
772 : // Technique similar to "Kalman update" of measurement used at boundaries - https://en.wikipedia.org/wiki/Kalman_filter
773 : //
774 : // 1.) Make backup of original parameters
775 : // 2.) Book Kalman matrices
776 : // 3.) Loop over all measurements bins and update mesurements -adding boundary measurements as additional measurement
777 : // relWeight vector specify relative weight of such measurement (err_i=sigma_i*refWeight_i) - not yet implemented
778 : // 4.) replace original parameters with constrained parameters
779 : // procedure can be repeated
780 : /*
781 : Input parameters example:
782 : nDims=2;
783 : Int_t indexes[2]={0,1};
784 : Double_t relWeight0[6]={1,1,1,1,1,1};
785 : Double_t relWeight1[6]={1,1,10,1,1,10};
786 : pcstream=new TTreeSRedirector("constrainStream.root","recreate");
787 :
788 : AliNDLocalRegression * regression0 = ( AliNDLocalRegression *)AliNDLocalRegression::GetVisualCorrections()->FindObject("pfitNDGaus0");
789 : AliNDLocalRegression * regression1 = ( AliNDLocalRegression *)AliNDLocalRegression::GetVisualCorrections()->FindObject("pfitNDGaus1");
790 :
791 : regressionUpdate0 = (AliNDLocalRegression *)regression0->Clone();
792 : regressionUpdate1 = (AliNDLocalRegression *)regression1->Clone();
793 : AddWeekConstrainsAtBoundaries( regressionUpdate0, nDims, indexes,relWeight0, pcstream);
794 : AddWeekConstrainsAtBoundaries( regressionUpdate0, nDims, indexes,relWeight0, pcstream);
795 : AddWeekConstrainsAtBoundaries( regressionUpdate0, nDims, indexes,relWeight0, pcstream);
796 : AddWeekConstrainsAtBoundaries( regressionUpdate0, nDims, indexes,relWeight0, pcstream);
797 : AddWeekConstrainsAtBoundaries( regressionUpdate1, nDims, indexes,relWeight1, pcstream);
798 : AddWeekConstrainsAtBoundaries( regressionUpdate1, nDims, indexes,relWeight1, pcstream);
799 : AddWeekConstrainsAtBoundaries( regressionUpdate1, nDims, indexes,relWeight1, pcstream);
800 : AddWeekConstrainsAtBoundaries( regressionUpdate1, nDims, indexes,relWeight1, pcstream);
801 :
802 : regressionUpdate0->SetName("pfitNDGaus0_Updated");
803 : regressionUpdate1->SetName("pfitNDGaus1_Updated");
804 : AliNDLocalRegression::AddVisualCorrection(regressionUpdate0);
805 : AliNDLocalRegression::AddVisualCorrection(regressionUpdate1);
806 : treeIn->SetAlias( regressionUpdate0->GetName(), TString::Format("AliNDLocalRegression::GetCorrND(%d,xyz0,xyz1+0)", regressionUpdate0->GetVisualCorrectionIndex()).Data());
807 : treeIn->SetAlias( regressionUpdate1->GetName(), TString::Format("AliNDLocalRegression::GetCorrND(%d,xyz0,xyz1+0)", regressionUpdate1->GetVisualCorrectionIndex()).Data());
808 : delete pcstream;
809 :
810 :
811 : TFile *f = TFile::Open("constrainStream.root")
812 : */
813 : const Double_t kScale=0.5;
814 : const Double_t singularity_tolerance = 1e-200;
815 0 : if (fLocalFitCovar==NULL) {
816 0 : ::Error("AliNDLocalRegression::EvalError", "Covariance matrix not available");
817 0 : return 0 ;
818 : }
819 : //
820 : // 1.) Make backup of original parameters
821 : //
822 0 : const TObjArray *vecParamOrig = fLocalFitParam;
823 : const TObjArray *vecCovarOrig = fLocalFitCovar;
824 0 : TObjArray *vecParamUpdated = new TObjArray(fHistPoints->GetNbins());
825 0 : TObjArray *vecCovarUpdated = new TObjArray(fHistPoints->GetNbins());
826 : //
827 : // 2.) Book local varaibles and Kalman matrices
828 : //
829 0 : Int_t nParams= 1+2*fNParameters;
830 0 : Int_t nMeas= nDims*6; // update each dimension specified 2 ends 2 measurements (value and first derivative)
831 :
832 0 : TMatrixD matWeight(nParams,nParams); // weight matrix for side param
833 0 : TMatrixD matCovarSide(nParams,nParams); // reweighted covariance matrix for side parameters
834 :
835 0 : TMatrixD vecXk(nParams,1); // X vector - parameter of the local fit at bin
836 0 : TMatrixD covXk(nParams,nParams); // X covariance
837 0 : TMatrixD matHk(nMeas,nParams); // vector to mesurement (values at boundary of bin)
838 0 : TMatrixD measR(nMeas,nMeas); // measurement error at boundary as provided by bin in local neigberhood
839 0 : TMatrixD vecZk(nMeas,1); // measurement at boundary
840 : //
841 0 : TMatrixD measRBin(nMeas,nMeas); // measurement error bin
842 0 : TMatrixD vecZkBin(nMeas,1); // measurement bin
843 0 : TMatrixD matrixTransformBin(nMeas, nParams); // vector to measurement to calculate error matrix current bin
844 : //
845 0 : TMatrixD vecZkSide(3,1); // measurement side
846 0 : TMatrixD matrixTransformSide(3,nParams);// vector to measurement to calculate error matrix side bin
847 :
848 : //
849 0 : TMatrixD vecYk(nMeas,1); // Innovation or measurement residual
850 0 : TMatrixD matHkT(nParams,nMeas);
851 0 : TMatrixD matSk(nMeas,nMeas); // Innovation (or residual) covariance
852 0 : TMatrixD matKk(nParams,nMeas); // Optimal Kalman gain
853 0 : TMatrixD mat1(nParams,nParams); // update covariance matrix
854 0 : TMatrixD covXk2(nParams,nParams); //
855 0 : TMatrixD covOut(nParams,nParams); //
856 0 : mat1.UnitMatrix();
857 : //
858 : // 3.) Loop over all measurements bins and update mesurements -adding boundary measurements as additional measurement
859 : // relWeight vector specify relative weight of such measurement (err_i=sigma_i*refWeight_i
860 0 : const THn* his = GetHistogram();
861 0 : Int_t binIndex[999]={0};
862 0 : Int_t binIndexSide[999]={0};
863 0 : Int_t nbinsAxis[999]={0};
864 0 : Double_t binCenter[999]={0};
865 0 : Double_t binWidth[999]={0};
866 :
867 0 : if (relWeight!=NULL) for (Int_t iParam=0; iParam<nParams; iParam++){
868 : Int_t index=0;
869 0 : if (iParam<3) index=iParam;
870 0 : if (iParam>=3) {
871 0 : Int_t dim=(iParam-3)/2;
872 0 : Int_t deriv=1+(iParam-3)%2;
873 0 : index=3*dim+deriv;
874 0 : }
875 0 : matWeight(iParam,iParam)=relWeight[index];
876 0 : }
877 :
878 :
879 0 : for (Int_t iDim=0; iDim<nDims; iDim++){nbinsAxis[iDim]=his->GetAxis(iDim)->GetNbins();}
880 : // Int_t nBins=fHistPoints->GetNbins();
881 0 : Int_t nBins=fLocalFitParam->GetSize();
882 0 : for (Int_t iBin=0; iBin<nBins; iBin++){ // loop over bins
883 0 : if (iBin%fgVerboseLevel==0) printf("%d\n",iBin);
884 : //
885 0 : his->GetBinContent(iBin,binIndex);
886 0 : for (Int_t iDim=0; iDim<nDims; iDim++) { // fill common info for bin of interest
887 0 : binCenter[iDim]= his->GetAxis(iDim)->GetBinCenter(binIndex[iDim]);
888 0 : binWidth[iDim] = his->GetAxis(iDim)->GetBinWidth(binIndex[iDim]);
889 : }
890 0 : if (fLocalFitParam->UncheckedAt(iBin)==NULL) continue;
891 0 : Double_t *vecParam0 = ((TVectorD*)(fLocalFitParam->UncheckedAt(iBin)))->GetMatrixArray();
892 0 : TMatrixD matParam0(nParams,1, vecParam0);
893 0 : TMatrixD & matCovar0=*(((TMatrixD*)(fLocalFitCovar->UncheckedAt(iBin))));
894 0 : measR.Zero();
895 0 : vecZk.Zero();
896 0 : measRBin.Zero();
897 0 : vecZkBin.Zero();
898 0 : matrixTransformBin.Zero();
899 0 : covXk=matCovar0;
900 0 : vecXk=matParam0;
901 : //
902 : // neiborhood loop
903 0 : Int_t constCounter=0;
904 0 : Int_t constCounterDim[100]={0};
905 0 : for (Int_t iDim=0; iDim<nDims; iDim++){ // loop in n dim
906 0 : constCounterDim[iDim]=0; // number of constraints per dimension
907 0 : for (Int_t iSide=-1; iSide<=1; iSide+=2){ // left right loop
908 0 : for (Int_t jDim=0; jDim<nDims; jDim++) binIndexSide[jDim]= binIndex[jDim];
909 0 : vecZkSide.Zero();
910 0 : matrixTransformSide.Zero();
911 : //
912 0 : binIndexSide[iDim]+=iSide;
913 0 : if (binIndexSide[iDim]<0) binIndexSide[iDim]=0;
914 0 : if (binIndexSide[iDim]>his->GetAxis(iDim)->GetNbins()) binIndexSide[iDim]=his->GetAxis(iDim)->GetNbins();
915 0 : Bool_t isConst=binIndexSide[iDim]>0 &&binIndexSide[iDim]<=his->GetAxis(iDim)->GetNbins() && (fLocalFitParam)->UncheckedAt(his->GetBin(binIndexSide))!=NULL;
916 0 : Int_t binSide=his->GetBin(binIndexSide);
917 0 : if (binSide>=nBins ) binIndexSide[iDim]=binIndex[iDim];
918 0 : if ((fLocalFitParam)->UncheckedAt(his->GetBin(binIndexSide))==NULL) binIndexSide[iDim]=binIndex[iDim];
919 0 : if (isConst) {
920 0 : constCounter++;
921 0 : constCounterDim[iDim]++;
922 0 : }
923 0 : Double_t localCenter=his->GetAxis(iDim)->GetBinCenter(binIndex[iDim]);
924 0 : Double_t sideCenter= his->GetAxis(iDim)->GetBinCenter(binIndexSide[iDim]);
925 0 : Double_t position= (iSide<0) ? his->GetAxis(iDim)->GetBinLowEdge(binIndex[iDim]) : his->GetAxis(iDim)->GetBinUpEdge(binIndex[iDim]);
926 0 : Double_t* vecParamSide = ((TVectorD*)(fLocalFitParam)->UncheckedAt(his->GetBin(binIndexSide)))->GetMatrixArray();
927 0 : TMatrixD matParamSide(nParams,1, vecParamSide);
928 0 : if (relWeight==NULL){
929 0 : matCovarSide=*((TMatrixD*)(fLocalFitCovar->UncheckedAt(his->GetBin(binIndexSide))));
930 : }
931 0 : if (relWeight!=NULL){
932 0 : matCovarSide=TMatrixD( matWeight,TMatrixD::kMult,*((TMatrixD*)(fLocalFitCovar->UncheckedAt(his->GetBin(binIndexSide)))));
933 0 : matCovarSide*=matWeight;
934 : }
935 0 : if (!isConst) matCovarSide*=1000;
936 : //
937 0 : Double_t deltaLocal=(position-localCenter);
938 0 : Double_t deltaSide=(position-sideCenter);
939 0 : if (fUseBinNorm){
940 0 : deltaLocal/=fBinWidth[iDim];
941 0 : deltaSide/=fBinWidth[iDim];
942 0 : }
943 : //
944 0 : matrixTransformSide(0,0)=1; matrixTransformSide(0,1+2*iDim)=deltaSide; matrixTransformSide(0,1+2*iDim+1)=deltaSide*deltaSide;
945 0 : matrixTransformSide(1,1+2*iDim)=1; matrixTransformSide(1,1+2*iDim+1)=2*deltaSide;
946 0 : matrixTransformSide(2,1+2*iDim+1)=2;
947 : //
948 0 : Int_t iMeas0=6*iDim+3*(iSide+1)/2;
949 0 : matrixTransformBin(iMeas0+0,0)=1; matrixTransformBin(iMeas0+0,1+2*iDim)=deltaLocal; matrixTransformBin(iMeas0+0,1+2*iDim+1)=deltaSide*deltaLocal;
950 0 : matrixTransformBin(iMeas0+1,1+2*iDim)=1; matrixTransformBin(iMeas0+1,1+2*iDim+1)=2*deltaLocal;
951 0 : matrixTransformBin(iMeas0+2,1+2*iDim+1)=2;
952 : //
953 0 : for (Int_t iconst=0; iconst<3; iconst++){
954 0 : Int_t iMeas=iMeas0+iconst;
955 : Double_t localMeasurement=0;
956 : Double_t sideMeasurement=0;
957 0 : if (iconst==0){ // measurement - derivative 0
958 0 : localMeasurement=vecParam0[0]+deltaLocal*(vecParam0[1+2*iDim]+vecParam0[2+2*iDim]*deltaLocal);
959 0 : sideMeasurement=vecParamSide[0]+deltaSide*(vecParamSide[1+2*iDim]+vecParamSide[2+2*iDim]*deltaSide);
960 0 : }
961 0 : if (iconst==1){ // measurement -derivative 1
962 0 : localMeasurement=(vecParam0[1+2*iDim]+2*vecParam0[2+2*iDim]*deltaLocal);
963 0 : sideMeasurement=(vecParamSide[1+2*iDim]+2*vecParamSide[2+2*iDim]*deltaSide);
964 0 : }
965 0 : if (iconst==2){
966 0 : localMeasurement=2*vecParam0[2+2*iDim];
967 0 : sideMeasurement=2*vecParamSide[2+2*iDim];
968 0 : }
969 0 : vecZkSide(iconst,0)=sideMeasurement;
970 0 : vecZk(iMeas,0)=sideMeasurement;
971 0 : if (!isConst) vecZk(iMeas,0)=localMeasurement;
972 0 : vecZkBin(iMeas,0)=localMeasurement;
973 : }
974 0 : TMatrixD measRSide0(matrixTransformSide,TMatrixD::kMult,matCovarSide); // (iconst,iconst) = (iconst,nParam)*(nParams,nParams)*(nParams,iconst
975 0 : TMatrixD matrixTransformSideT(TMatrixD::kTransposed ,matrixTransformSide);
976 0 : TMatrixD measRSide(measRSide0,TMatrixD::kMult,matrixTransformSideT);
977 : // update measutement Covariance matrix for given side
978 0 : for (Int_t iconst=0; iconst<3; iconst++)
979 0 : for (Int_t jconst=0; jconst<3; jconst++){
980 0 : measR(iMeas0+iconst,iMeas0+jconst)=measRSide(iconst,jconst);
981 : }
982 0 : if (pcstream){
983 0 : TMatrixD vecZkSideCheck(matrixTransformSide,TMatrixD::kMult,matParamSide); // (iconst,1) = (iConst,nParam)*(nParams,1)
984 : //
985 0 : (*pcstream)<<"checkSide"<< // check agreement in 1D
986 0 : "iBin="<<iBin<<
987 0 : "iDim="<<iDim<<
988 0 : "iSide="<<iSide<<
989 0 : "vecZkSide.="<<&vecZkSide<<
990 0 : "vecZkSideCheck.="<<&vecZkSideCheck<<
991 0 : "measRSide.="<<&measRSide<<
992 0 : "vecZk.="<<&vecZk<<
993 0 : "vecZkBin.="<<&vecZkBin<<
994 : "\n";
995 0 : }
996 0 : }
997 : }
998 : //
999 : //
1000 0 : TMatrixD measRBin0(matrixTransformBin,TMatrixD::kMult,matCovar0); // (iconst,iconst) = (iconst,nParam)*(nParams,nParams)*(nParams,iconst
1001 0 : TMatrixD matrixTransformBinT(TMatrixD::kTransposed ,matrixTransformBin);
1002 0 : TMatrixD measRBin(measRBin0,TMatrixD::kMult,matrixTransformBinT);
1003 : //
1004 : // make Kalman Update of state vector with side mesurement
1005 : //
1006 0 : matHk=matrixTransformBin;
1007 0 : matHkT= matrixTransformBinT;
1008 : //
1009 0 : vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
1010 0 : if (useCommon) vecYk*=0.5; // in case we are using middle point use only half of delta
1011 0 : matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
1012 0 : Double_t determinant=0;
1013 : Int_t constCounter2=0;
1014 0 : for (Int_t kDim=0; kDim<nDims; kDim++) if (constCounterDim[kDim]>0) constCounter2++;
1015 0 : if (constCounter2==nDims){
1016 0 : determinant= matSk.Determinant();
1017 0 : }
1018 0 : if (TMath::Abs(determinant)<singularity_tolerance ) {
1019 0 : vecParamUpdated->AddAt(new TVectorD(*((TVectorD*)(fLocalFitParam->UncheckedAt(iBin)))),iBin);
1020 0 : vecCovarUpdated->AddAt(new TMatrixD(*((TMatrixD*)(fLocalFitCovar->UncheckedAt(iBin)))),iBin);
1021 0 : AliDebug(1,TString::Format("Update matrix not possible matSk.Determinant() too small, skipping bin %d",iBin).Data());
1022 : }else{
1023 0 : matSk.Invert();
1024 0 : matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
1025 0 : vecXk += matKk*vecYk; // updated vector
1026 0 : covXk2 = (mat1-(matKk*matHk));
1027 0 : covOut = covXk2*covXk;
1028 : //
1029 0 : vecParamUpdated->AddAt(new TVectorD(nParams,vecXk.GetMatrixArray()), iBin);
1030 0 : vecCovarUpdated->AddAt(new TMatrixD(covOut), iBin);
1031 : }
1032 0 : if (pcstream){
1033 0 : TMatrixD vecZkBinCheck(matrixTransformBin,TMatrixD::kMult,matParam0);
1034 0 : TVectorD vecPos(nDims,binCenter);
1035 0 : TVectorD *vecXk0= (TVectorD*)(fLocalFitParam->UncheckedAt(iBin));
1036 0 : TMatrixD vecYkUpdated=(vecZk-matHk*vecXk);
1037 : //
1038 0 : (*pcstream)<<"checkBin"<< // check agreement in all sides
1039 0 : "iBin="<<iBin<< // bin index
1040 0 : "vecPos.="<<&vecPos<< // bin position
1041 0 : "determinant="<<determinant<< // determinant
1042 0 : "constCounter="<<constCounter<<
1043 0 : "constCounter2="<<constCounter<<
1044 : //
1045 0 : "vecXk0.="<<vecXk0<< // original parameter vector
1046 0 : "vecXk.="<<&vecXk<< // parameter vector at bin after update
1047 0 : "covXk.="<<&covXk<< // covaraince matrix before update
1048 0 : "covOut.="<<&covOut<< // covaraince matrix after update
1049 0 : "vecZk.="<<&vecZk<< // measurement vector - values according side measurement
1050 0 : "vecZkBin.="<<&vecZkBin<< // expected vector according parameters for bin
1051 0 : "vecZkBinCheck.="<<&vecZkBinCheck<< // expected vector according parameters at bin centers - crosscheck tracsrormation matrix
1052 0 : "measRBin.="<<&measRBin<< // expected error of extrapolation
1053 0 : "measR.="<<&measR<< // error of the side measurement
1054 : // tmporary data
1055 0 : "vecYk.="<<&vecYk<< // delta vector (nparams)
1056 0 : "matSk.="<<&matSk<< // inovation covariance (nMeas,nMeas)
1057 0 : "matKk.="<<&matKk<< // optimal Kalman gain (nParams,nMeas)
1058 0 : "covXk2.="<<&covXk2<<
1059 : //
1060 0 : "vecYkUpdated.="<<&vecYkUpdated<< // diff after kalman update
1061 : "\n";
1062 0 : }
1063 0 : }
1064 : //
1065 : // 4.) replace original parameters with constrained parameters
1066 : //
1067 0 : fLocalFitParam= vecParamUpdated;
1068 0 : fLocalFitCovar= vecCovarUpdated;
1069 : return 0;
1070 0 : }
1071 :
1072 : void AliNDLocalRegression::DumpToTree(Int_t nDiv, TTreeStream & stream){
1073 : //
1074 : //
1075 : //
1076 : const Int_t kMaxDim=100;
1077 0 : Int_t nBins=fHistPoints->GetNbins();
1078 :
1079 0 : TVectorD binLowEdge(fNParameters);
1080 0 : TVectorD binLocal(fNParameters);
1081 0 : TVectorF binIndexF(fNParameters);
1082 0 : Double_t *pbinLowEdge= binLowEdge.GetMatrixArray();
1083 : //
1084 0 : for (Int_t iBin=0; iBin<nBins; iBin++){ // loop over bins
1085 0 : if (iBin%fgVerboseLevel==0) printf("%d\n",iBin);
1086 0 : fHistPoints->GetBinContent(iBin,fBinIndex);
1087 0 : for (Int_t iDim=0; iDim<fNParameters; iDim++) { // fill common info for bin of interest
1088 0 : binIndexF[iDim]=fBinIndex[iDim];
1089 0 : fBinCenter[iDim]= fHistPoints->GetAxis(iDim)->GetBinCenter(fBinIndex[iDim]);
1090 0 : binLowEdge[iDim]= fHistPoints->GetAxis(iDim)->GetBinLowEdge(fBinIndex[iDim]);
1091 0 : fBinDelta[iDim] = fHistPoints->GetAxis(iDim)->GetBinWidth(fBinIndex[iDim]);
1092 : }
1093 :
1094 0 : for (Int_t iDim=0; iDim<fNParameters; iDim++){
1095 0 : for (Int_t jDim=0; jDim<fNParameters; jDim++) binLocal[jDim]=fBinCenter[jDim];
1096 0 : for (Int_t iDiv=0; iDiv<nDiv; iDiv++){
1097 0 : binLocal[iDim]=binLowEdge[iDim]+(fBinDelta[iDim]*(iDiv-1.))/Double_t(nDiv);
1098 0 : Double_t value2= Eval(binLocal.GetMatrixArray());
1099 0 : binLocal[iDim]=binLowEdge[iDim]+(fBinDelta[iDim]*Double_t(iDiv))/Double_t(nDiv);
1100 0 : Double_t value=Eval(binLocal.GetMatrixArray());
1101 0 : Double_t derivative=nDiv*(value-value2)/fBinDelta[iDim];
1102 0 : stream<<
1103 0 : "pos.="<<&binLocal<< // position vector
1104 0 : "iBin="<<iBin<< // bin index
1105 0 : "iDim="<<iDim<< // scan bin index
1106 0 : "iDiv="<<iDiv<< // division index
1107 0 : "binIndexF.="<<&binIndexF<< // bin index
1108 0 : "value="<<value<< // value at
1109 0 : "derivative="<<derivative<< // numerical derivative
1110 : "\n";
1111 0 : }
1112 : }
1113 : }
1114 0 : }
1115 :
1116 :
1117 : Double_t AliNDLocalRegression::EvalGraphKernel(TGraph * gr, Double_t evalTime, Double_t kernelWidth, Int_t sigmaCut, Bool_t evalLog, Int_t pol, TVectorD *param, TMatrixD *covar){
1118 : ///
1119 : ///
1120 : /// Interpolate graph using local regression with kernel width
1121 : ///
1122 : Int_t kMinEntries=4;
1123 0 : Int_t npoints=gr->GetN();
1124 0 : Int_t index0= TMath::BinarySearch(npoints, gr->GetX(), evalTime-sigmaCut*kernelWidth);
1125 0 : Int_t index1= TMath::BinarySearch(npoints, gr->GetX(), evalTime+sigmaCut*kernelWidth);
1126 0 : if (index1-index0 < kMinEntries) {
1127 0 : Int_t index= TMath::BinarySearch(npoints, gr->GetX(), evalTime);
1128 0 : index0=index-kMinEntries/2;
1129 0 : index1=index+kMinEntries/2;
1130 0 : }
1131 0 : if (index0<0) index0=0;
1132 0 : if (index1>=npoints) index1=npoints-1;
1133 0 : TLinearFitter fitter(pol+1, TString::Format("pol%d",pol));
1134 0 : Double_t mkernel2=1./(kernelWidth*kernelWidth);
1135 0 : for (Int_t ipoint=index0; ipoint<=index1; ipoint++){
1136 0 : Double_t x=gr->GetX()[ipoint]-evalTime;
1137 0 : Double_t y=gr->GetY()[ipoint];
1138 0 : if (evalLog==kFALSE){
1139 0 : fitter.AddPoint(&x, gr->GetY()[ipoint], TMath::Exp(x*x*mkernel2));
1140 : }else{
1141 0 : if (y>0) fitter.AddPoint(&x, TMath::Log(y) , TMath::Exp(x*x*mkernel2));
1142 : }
1143 0 : }
1144 0 : Int_t hasFailed=(fitter.GetNpoints()>kMinEntries)? fitter.Eval():1;
1145 0 : if (hasFailed) return 0;
1146 :
1147 0 : if (param){
1148 0 : fitter.GetParameters(*param);
1149 : }
1150 0 : if (covar) fitter.GetCovarianceMatrix(*covar);
1151 0 : return (evalLog==0) ? fitter.GetParameter(0) : TMath::Exp(fitter.GetParameter(0));
1152 0 : }
1153 :
1154 :
|
