Filling residuals (GlobalQA) only for tracks with kTPCin
[u/mrichter/AliRoot.git] / CORRFW / AliCFUnfolding.cxx
CommitLineData
c0b10ad4 1/**************************************************************************
2 * Copyright(c) 1998-1999, 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// //
18// AliCFUnfolding Class //
19// Class to handle general unfolding procedure //
20// For the moment only bayesian unfolding is supported //
21// The next steps are to add chi2 minimisation and weighting methods //
22// //
23// Author : renaud.vernet@cern.ch //
24//--------------------------------------------------------------------//
25
26
27#include "AliCFUnfolding.h"
28#include "TMath.h"
29#include "TAxis.h"
30#include "AliLog.h"
85b6bda9 31#include "TF1.h"
32#include "TH1D.h"
33#include "TH2D.h"
34#include "TH3D.h"
c0b10ad4 35
36ClassImp(AliCFUnfolding)
37
38//______________________________________________________________
39
40AliCFUnfolding::AliCFUnfolding() :
41 TNamed(),
42 fResponse(0x0),
43 fPrior(0x0),
c0b10ad4 44 fEfficiency(0x0),
45 fMeasured(0x0),
46 fMaxNumIterations(0),
47 fNVariables(0),
48 fMaxChi2(0),
49 fUseSmoothing(kFALSE),
85b6bda9 50 fSmoothFunction(0x0),
51 fSmoothOption(""),
52 fOriginalPrior(0x0),
c0b10ad4 53 fInverseResponse(0x0),
54 fMeasuredEstimate(0x0),
55 fConditional(0x0),
56 fProjResponseInT(0x0),
57 fUnfolded(0x0),
58 fCoordinates2N(0x0),
59 fCoordinatesN_M(0x0),
60 fCoordinatesN_T(0x0)
61{
62 //
63 // default constructor
64 //
65}
66
67//______________________________________________________________
68
69AliCFUnfolding::AliCFUnfolding(const Char_t* name, const Char_t* title, const Int_t nVar,
70 const THnSparse* response, const THnSparse* efficiency, const THnSparse* measured, const THnSparse* prior) :
71 TNamed(name,title),
72 fResponse((THnSparse*)response->Clone()),
73 fPrior(0x0),
c0b10ad4 74 fEfficiency((THnSparse*)efficiency->Clone()),
75 fMeasured((THnSparse*)measured->Clone()),
76 fMaxNumIterations(0),
77 fNVariables(nVar),
78 fMaxChi2(0),
79 fUseSmoothing(kFALSE),
85b6bda9 80 fSmoothFunction(0x0),
81 fSmoothOption(""),
82 fOriginalPrior(0x0),
c0b10ad4 83 fInverseResponse(0x0),
84 fMeasuredEstimate(0x0),
85 fConditional(0x0),
86 fProjResponseInT(0x0),
87 fUnfolded(0x0),
88 fCoordinates2N(0x0),
89 fCoordinatesN_M(0x0),
90 fCoordinatesN_T(0x0)
91{
92 //
93 // named constructor
94 //
95
96 AliInfo(Form("\n\n--------------------------\nCreating an unfolder :\n--------------------------\nresponse matrix has %d dimension(s)",fResponse->GetNdimensions()));
97
98 if (!prior) CreateFlatPrior(); // if no prior distribution declared, simply use a flat distribution
99 else {
100 fPrior = (THnSparse*) prior->Clone();
101 fOriginalPrior = (THnSparse*)fPrior->Clone();
85b6bda9 102 if (fPrior->GetNdimensions() != fNVariables)
103 AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions()));
c0b10ad4 104 }
85b6bda9 105
106 if (fEfficiency->GetNdimensions() != fNVariables)
107 AliFatal(Form("The efficiency matrix should have %d dimensions, and it has actually %d",fNVariables,fEfficiency->GetNdimensions()));
108 if (fMeasured->GetNdimensions() != fNVariables)
109 AliFatal(Form("The measured matrix should have %d dimensions, and it has actually %d",fNVariables,fMeasured->GetNdimensions()));
110 if (fResponse->GetNdimensions() != 2*fNVariables)
111 AliFatal(Form("The response matrix should have %d dimensions, and it has actually %d",2*fNVariables,fResponse->GetNdimensions()));
c0b10ad4 112
85b6bda9 113
c0b10ad4 114 for (Int_t iVar=0; iVar<fNVariables; iVar++) {
115 AliInfo(Form("prior matrix has %d bins in dimension %d",fPrior ->GetAxis(iVar)->GetNbins(),iVar));
116 AliInfo(Form("efficiency matrix has %d bins in dimension %d",fEfficiency->GetAxis(iVar)->GetNbins(),iVar));
117 AliInfo(Form("measured matrix has %d bins in dimension %d",fMeasured ->GetAxis(iVar)->GetNbins(),iVar));
118 }
119 Init();
120}
121
122
123//______________________________________________________________
124
125AliCFUnfolding::AliCFUnfolding(const AliCFUnfolding& c) :
126 TNamed(c),
127 fResponse((THnSparse*)c.fResponse->Clone()),
128 fPrior((THnSparse*)c.fPrior->Clone()),
c0b10ad4 129 fEfficiency((THnSparse*)c.fEfficiency->Clone()),
130 fMeasured((THnSparse*)c.fMeasured->Clone()),
131 fMaxNumIterations(c.fMaxNumIterations),
132 fNVariables(c.fNVariables),
133 fMaxChi2(c.fMaxChi2),
134 fUseSmoothing(c.fUseSmoothing),
85b6bda9 135 fSmoothFunction((TF1*)c.fSmoothFunction->Clone()),
136 fSmoothOption(fSmoothOption),
137 fOriginalPrior((THnSparse*)c.fOriginalPrior->Clone()),
c0b10ad4 138 fInverseResponse((THnSparse*)c.fInverseResponse->Clone()),
139 fMeasuredEstimate((THnSparse*)fMeasuredEstimate->Clone()),
140 fConditional((THnSparse*)c.fConditional->Clone()),
141 fProjResponseInT((THnSparse*)c.fProjResponseInT->Clone()),
142 fUnfolded((THnSparse*)c.fUnfolded->Clone()),
143 fCoordinates2N(new Int_t(*c.fCoordinates2N)),
144 fCoordinatesN_M(new Int_t(*c.fCoordinatesN_M)),
145 fCoordinatesN_T(new Int_t(*c.fCoordinatesN_T))
146{
147 //
148 // copy constructor
149 //
150}
151
152//______________________________________________________________
153
154AliCFUnfolding& AliCFUnfolding::operator=(const AliCFUnfolding& c) {
155 //
156 // assignment operator
157 //
158
159 if (this!=&c) {
160 TNamed::operator=(c);
161 fResponse = (THnSparse*)c.fResponse->Clone() ;
162 fPrior = (THnSparse*)c.fPrior->Clone() ;
c0b10ad4 163 fEfficiency = (THnSparse*)c.fEfficiency->Clone() ;
164 fMeasured = (THnSparse*)c.fMeasured->Clone() ;
165 fMaxNumIterations = c.fMaxNumIterations ;
166 fNVariables = c.fNVariables ;
167 fMaxChi2 = c.fMaxChi2 ;
168 fUseSmoothing = c.fUseSmoothing ;
85b6bda9 169 fSmoothFunction = (TF1*)c.fSmoothFunction->Clone();
170 fSmoothOption = c.fSmoothOption ;
171 fOriginalPrior = (THnSparse*)c.fOriginalPrior->Clone() ;
c0b10ad4 172 fInverseResponse = (THnSparse*)c.fInverseResponse->Clone() ;
173 fMeasuredEstimate = (THnSparse*)fMeasuredEstimate->Clone() ;
174 fConditional = (THnSparse*)c.fConditional->Clone() ;
175 fProjResponseInT = (THnSparse*)c.fProjResponseInT->Clone() ;
176 fUnfolded = (THnSparse*)c.fUnfolded->Clone() ;
177 fCoordinates2N = new Int_t(*c.fCoordinates2N) ;
178 fCoordinatesN_M = new Int_t(*c.fCoordinatesN_M) ;
179 fCoordinatesN_T = new Int_t(*c.fCoordinatesN_T) ;
180 }
181 return *this;
182}
183
184//______________________________________________________________
185
186AliCFUnfolding::~AliCFUnfolding() {
187 //
188 // destructor
189 //
190 if (fResponse) delete fResponse;
191 if (fPrior) delete fPrior;
c0b10ad4 192 if (fEfficiency) delete fEfficiency;
193 if (fMeasured) delete fMeasured;
85b6bda9 194 if (fSmoothFunction) delete fSmoothFunction;
195 if (fOriginalPrior) delete fOriginalPrior;
c0b10ad4 196 if (fInverseResponse) delete fInverseResponse;
197 if (fMeasuredEstimate) delete fMeasuredEstimate;
198 if (fConditional) delete fConditional;
199 if (fProjResponseInT) delete fProjResponseInT;
200 if (fCoordinates2N) delete [] fCoordinates2N;
201 if (fCoordinatesN_M) delete [] fCoordinatesN_M;
202 if (fCoordinatesN_T) delete [] fCoordinatesN_T;
203}
204
205//______________________________________________________________
206
207void AliCFUnfolding::Init() {
208 //
209 // initialisation function : creates internal settings
210 //
211
212 fCoordinates2N = new Int_t[2*fNVariables];
213 fCoordinatesN_M = new Int_t[fNVariables];
214 fCoordinatesN_T = new Int_t[fNVariables];
215
216 // create the matrix of conditional probabilities P(M|T)
217 CreateConditional();
218
219 // create the frame of the inverse response matrix
220 fInverseResponse = (THnSparse*) fResponse->Clone();
221 // create the frame of the unfolded spectrum
222 fUnfolded = (THnSparse*) fPrior->Clone();
223 // create the frame of the measurement estimate spectrum
224 fMeasuredEstimate = (THnSparse*) fMeasured->Clone();
225}
226
227//______________________________________________________________
228
229void AliCFUnfolding::CreateEstMeasured() {
230 //
231 // This function creates a estimate (M) of the reconstructed spectrum
232 // given the a priori distribution (T) and the conditional matrix (COND)
233 //
234 // --> P(M) = SUM { P(M|T) * P(T) }
235 // --> M(i) = SUM_k { COND(i,k) * T(k) }
236 //
237 // This is needed to calculate the inverse response matrix
238 //
239
240
241 // clean the measured estimate spectrum
242 for (Long64_t i=0; i<fMeasuredEstimate->GetNbins(); i++) {
243 fMeasuredEstimate->GetBinContent(i,fCoordinatesN_M);
244 fMeasuredEstimate->SetBinContent(fCoordinatesN_M,0.);
85b6bda9 245 fMeasuredEstimate->SetBinError (fCoordinatesN_M,0.);
c0b10ad4 246 }
247
85b6bda9 248 THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone();
249 priorTimesEff->Multiply(fEfficiency);
250
c0b10ad4 251 // fill it
252 for (Int_t iBin=0; iBin<fConditional->GetNbins(); iBin++) {
253 Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N);
85b6bda9 254 Double_t conditionalError = fConditional->GetBinError (iBin);
c0b10ad4 255 GetCoordinates();
85b6bda9 256 Double_t priorTimesEffValue = priorTimesEff->GetBinContent(fCoordinatesN_T);
257 Double_t priorTimesEffError = priorTimesEff->GetBinError (fCoordinatesN_T);
258 Double_t fill = conditionalValue * priorTimesEffValue ;
259
260 if (fill>0.) {
261 fMeasuredEstimate->AddBinContent(fCoordinatesN_M,fill);
262
263 // error calculation : gaussian error propagation (may be overestimated...)
264 Double_t err2 = TMath::Power(fMeasuredEstimate->GetBinError(fCoordinatesN_M),2) ;
265 err2 += TMath::Power(conditionalValue*priorTimesEffError,2) + TMath::Power(conditionalError*priorTimesEffValue,2) ;
266 Double_t err = TMath::Sqrt(err2);
267 fMeasuredEstimate->SetBinError(fCoordinatesN_M,err);
268 }
c0b10ad4 269 }
85b6bda9 270 delete priorTimesEff ;
c0b10ad4 271}
272
273//______________________________________________________________
274
275void AliCFUnfolding::CreateInvResponse() {
276 //
277 // Creates the inverse response matrix (INV) with Bayesian method
278 // : uses the conditional matrix (COND), the prior probabilities (T) and the efficiency map (E)
279 //
280 // --> P(T|M) = P(M|T) * P(T) * eff(T) / SUM { P(M|T) * P(T) }
281 // --> INV(i,j) = COND(i,j) * T(j) * E(j) / SUM_k { COND(i,k) * T(k) }
282 //
283
85b6bda9 284 THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone();
285 priorTimesEff->Multiply(fEfficiency);
286
c0b10ad4 287 for (Int_t iBin=0; iBin<fConditional->GetNbins(); iBin++) {
288 Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N);
85b6bda9 289 Double_t conditionalError = fConditional->GetBinError (iBin);
c0b10ad4 290 GetCoordinates();
85b6bda9 291 Double_t estMeasuredValue = fMeasuredEstimate->GetBinContent(fCoordinatesN_M);
292 Double_t estMeasuredError = fMeasuredEstimate->GetBinError (fCoordinatesN_M);
293 Double_t priorTimesEffValue = priorTimesEff ->GetBinContent(fCoordinatesN_T);
294 Double_t priorTimesEffError = priorTimesEff ->GetBinError (fCoordinatesN_T);
295 Double_t fill = (estMeasuredValue>0. ? conditionalValue * priorTimesEffValue / estMeasuredValue : 0. ) ;
296 // error calculation : gaussian error propagation (may be overestimated...)
297 Double_t err = 0. ;
298 if (estMeasuredValue>0.) {
299 err = TMath::Sqrt( TMath::Power(conditionalError * priorTimesEffValue * estMeasuredValue ,2) +
300 TMath::Power(conditionalValue * priorTimesEffError * estMeasuredValue ,2) +
301 TMath::Power(conditionalValue * priorTimesEffValue * estMeasuredError ,2) )
302 / TMath::Power(estMeasuredValue,2) ;
303 }
304 if (fill>0. || fInverseResponse->GetBinContent(fCoordinates2N)>0.) {
305 fInverseResponse->SetBinContent(fCoordinates2N,fill);
306 fInverseResponse->SetBinError (fCoordinates2N,err );
307 }
308 }
309 delete priorTimesEff ;
c0b10ad4 310}
311
312//______________________________________________________________
313
314void AliCFUnfolding::Unfold() {
315 //
316 // Main routine called by the user :
317 // it calculates the unfolded spectrum from the response matrix and the measured spectrum
318 // several iterations are performed until a reasonable chi2 is reached
319 //
320
321 Int_t iIterBayes=0 ;
322 Double_t chi2=0 ;
323
324 for (iIterBayes=0; iIterBayes<fMaxNumIterations; iIterBayes++) { // bayes iterations
325 CreateEstMeasured();
326 CreateInvResponse();
327 CreateUnfolded();
328 chi2 = GetChi2();
c0b10ad4 329 if (fMaxChi2>0. && chi2<fMaxChi2) {
330 break;
331 }
332 // update the prior distribution
85b6bda9 333 if (fUseSmoothing) {
334 if (Smooth()) {
335 AliError("Couldn't smooth the unfolded spectrum!!");
336 return;
337 }
338 }
c0b10ad4 339 fPrior = (THnSparse*)fUnfolded->Clone() ; // this should be changed (memory)
340 }
85b6bda9 341 AliInfo(Form("\n\n=======================\nFinished at iteration %d : Chi2 is %e and you required it to be < %e\n=======================\n\n",iIterBayes,chi2,fMaxChi2));
c0b10ad4 342}
343
344//______________________________________________________________
345
346void AliCFUnfolding::CreateUnfolded() {
347 //
348 // Creates the unfolded (T) spectrum from the measured spectrum (M) and the inverse response matrix (INV)
349 // We have P(T) = SUM { P(T|M) * P(M) }
350 // --> T(i) = SUM_k { INV(i,k) * M(k) }
351 //
352
353
354 // clear the unfolded spectrum
355 for (Long64_t i=0; i<fUnfolded->GetNbins(); i++) {
356 fUnfolded->GetBinContent(i,fCoordinatesN_T);
357 fUnfolded->SetBinContent(fCoordinatesN_T,0.);
85b6bda9 358 fUnfolded->SetBinError (fCoordinatesN_T,0.);
c0b10ad4 359 }
360
361 for (Int_t iBin=0; iBin<fInverseResponse->GetNbins(); iBin++) {
362 Double_t invResponseValue = fInverseResponse->GetBinContent(iBin,fCoordinates2N);
85b6bda9 363 Double_t invResponseError = fInverseResponse->GetBinError (iBin);
c0b10ad4 364 GetCoordinates();
85b6bda9 365 Double_t effValue = fEfficiency->GetBinContent(fCoordinatesN_T);
366 Double_t effError = fEfficiency->GetBinError (fCoordinatesN_T);
367 Double_t measuredValue = fMeasured ->GetBinContent(fCoordinatesN_M);
368 Double_t measuredError = fMeasured ->GetBinError (fCoordinatesN_M);
369 Double_t fill = (effValue>0. ? invResponseValue * measuredValue / effValue : 0.) ;
370
371 if (fill>0.) {
372 fUnfolded->AddBinContent(fCoordinatesN_T,fill);
373
374 // error calculation : gaussian error propagation (may be overestimated...)
375 Double_t err2 = TMath::Power(fUnfolded->GetBinError(fCoordinatesN_T),2) ;
376 err2 += TMath::Power(invResponseError * measuredValue * effValue,2) / TMath::Power(effValue,4) ;
377 err2 += TMath::Power(invResponseValue * measuredError * effValue,2) / TMath::Power(effValue,4) ;
378 err2 += TMath::Power(invResponseValue * measuredValue * effError,2) / TMath::Power(effValue,4) ;
379 Double_t err = TMath::Sqrt(err2);
380 fUnfolded->SetBinError(fCoordinatesN_T,err);
381 }
c0b10ad4 382 }
383}
384
85b6bda9 385//______________________________________________________________
386
c0b10ad4 387void AliCFUnfolding::GetCoordinates() {
388 //
389 // assign coordinates in Measured and True spaces (dim=N) from coordinates in global space (dim=2N)
390 //
391 for (Int_t i = 0; i<fNVariables ; i++) {
392 fCoordinatesN_M[i] = fCoordinates2N[i];
393 fCoordinatesN_T[i] = fCoordinates2N[i+fNVariables];
394 }
395}
396
397//______________________________________________________________
398
399void AliCFUnfolding::CreateConditional() {
400 //
401 // creates the conditional probability matrix (R*) holding the P(M|T), given the reponse matrix R
402 //
403 // --> R*(i,j) = R(i,j) / SUM_k{ R(k,j) }
404 //
405
406 fConditional = (THnSparse*) fResponse->Clone(); // output of this function
407 fProjResponseInT = (THnSparse*) fPrior->Clone(); // output denominator :
408 // projection of the response matrix on the TRUE axis
85b6bda9 409 Int_t* dim = new Int_t [fNVariables];
410 for (Int_t iDim=0; iDim<fNVariables; iDim++) dim[iDim] = fNVariables+iDim ; //dimensions corresponding to TRUE values (i.e. from N to 2N-1)
411 fProjResponseInT = fConditional->Projection(fNVariables,dim,"E"); //project
412 delete [] dim;
c0b10ad4 413
414 // fill the conditional probability matrix
415 for (Int_t iBin=0; iBin<fResponse->GetNbins(); iBin++) {
416 Double_t responseValue = fResponse->GetBinContent(iBin,fCoordinates2N);
85b6bda9 417 Double_t responseError = fResponse->GetBinError (iBin);
c0b10ad4 418 GetCoordinates();
85b6bda9 419 Double_t projValue = fProjResponseInT->GetBinContent(fCoordinatesN_T);
420 Double_t projError = fProjResponseInT->GetBinError (fCoordinatesN_T);
421
422 Double_t fill = responseValue / projValue ;
423 if (fill>0. || fConditional->GetBinContent(fCoordinates2N)>0.) {
424 fConditional->SetBinContent(fCoordinates2N,fill);
425 // gaussian error for the moment
426 Double_t err2 = TMath::Power(responseError*projValue,2) + TMath::Power(responseValue*projError,2) ;
427 Double_t err = TMath::Sqrt(err2);
428 err /= TMath::Power(projValue,2) ;
429 fConditional->SetBinError (fCoordinates2N,err);
430 }
c0b10ad4 431 }
432}
433
434//______________________________________________________________
435
436Double_t AliCFUnfolding::GetChi2() {
437 //
438 // Returns the chi2 between unfolded and a priori spectrum
439 //
440
441 Double_t chi2 = 0. ;
442 for (Int_t iBin=0; iBin<fPrior->GetNbins(); iBin++) {
443 Double_t priorValue = fPrior->GetBinContent(iBin);
444 chi2 += (priorValue>0. ? TMath::Power(fUnfolded->GetBinContent(iBin) - priorValue,2) / priorValue : 0.) ;
445 }
446 return chi2;
447}
448
449//______________________________________________________________
450
451void AliCFUnfolding::SetMaxChi2PerDOF(Double_t val) {
452 //
453 // Max. chi2 per degree of freedom : user setting
454 //
455
456 Int_t nDOF = 1 ;
457 for (Int_t iDim=0; iDim<fNVariables; iDim++) {
458 nDOF *= fPrior->GetAxis(iDim)->GetNbins();
459 }
460 AliInfo(Form("Number of degrees of freedom = %d",nDOF));
461 fMaxChi2 = val * nDOF ;
462}
463
464//______________________________________________________________
465
85b6bda9 466Short_t AliCFUnfolding::Smooth() {
c0b10ad4 467 //
468 // Smoothes the unfolded spectrum
85b6bda9 469 //
470 // By default each cell content is replaced by the average with the neighbouring bins (but not diagonally-neighbouring bins)
471 // However, if a specific function fcn has been defined in UseSmoothing(fcn), the unfolded will be fit and updated using fcn
c0b10ad4 472 //
473
85b6bda9 474 if (fSmoothFunction) {
475 AliInfo(Form("Smoothing spectrum with fit function %p",fSmoothFunction));
476 return SmoothUsingFunction();
477 }
478 else return SmoothUsingNeighbours();
479}
480
481//______________________________________________________________
482
483Short_t AliCFUnfolding::SmoothUsingNeighbours() {
484 //
485 // Smoothes the unfolded spectrum using neighouring bins
486 //
487
c0b10ad4 488 Int_t* numBins = new Int_t[fNVariables];
489 for (Int_t iVar=0; iVar<fNVariables; iVar++) numBins[iVar]=fUnfolded->GetAxis(iVar)->GetNbins();
490
491 //need a copy because fUnfolded will be updated during the loop, and this creates problems
492 THnSparse* copy = (THnSparse*)fUnfolded->Clone();
493
494 for (Int_t iBin=0; iBin<copy->GetNbins(); iBin++) { //loop on non-empty bins
495 Double_t content = copy->GetBinContent(iBin,fCoordinatesN_T);
85b6bda9 496 Double_t error2 = TMath::Power(copy->GetBinContent(iBin),2);
c0b10ad4 497
498 // skip the under/overflow bins...
499 Bool_t isOutside = kFALSE ;
500 for (Int_t iVar=0; iVar<fNVariables; iVar++) {
501 if (fCoordinatesN_T[iVar]<1 || fCoordinatesN_T[iVar]>numBins[iVar]) {
502 isOutside=kTRUE;
503 break;
504 }
505 }
506 if (isOutside) continue;
507
508 Int_t neighbours = 0; // number of neighbours to average with
509
510 for (Int_t iVar=0; iVar<fNVariables; iVar++) {
511 if (fCoordinatesN_T[iVar] > 1) { // must not be on low edge border
512 fCoordinatesN_T[iVar]-- ; //get lower neighbouring bin
85b6bda9 513 content += copy->GetBinContent(fCoordinatesN_T);
514 error2 += TMath::Power(copy->GetBinError(fCoordinatesN_T),2);
c0b10ad4 515 neighbours++;
516 fCoordinatesN_T[iVar]++ ; //back to initial coordinate
517 }
518 if (fCoordinatesN_T[iVar] < numBins[iVar]) { // must not be on up edge border
519 fCoordinatesN_T[iVar]++ ; //get upper neighbouring bin
85b6bda9 520 content += copy->GetBinContent(fCoordinatesN_T);
521 error2 += TMath::Power(copy->GetBinError(fCoordinatesN_T),2);
c0b10ad4 522 neighbours++;
523 fCoordinatesN_T[iVar]-- ; //back to initial coordinate
524 }
525 }
85b6bda9 526 // make an average
527 fUnfolded->SetBinContent(fCoordinatesN_T,content/(1.+neighbours));
528 fUnfolded->SetBinError (fCoordinatesN_T,TMath::Sqrt(error2)/(1.+neighbours));
c0b10ad4 529 }
530 delete [] numBins;
531 delete copy;
85b6bda9 532 return 0;
c0b10ad4 533}
534
85b6bda9 535//______________________________________________________________
536
537Short_t AliCFUnfolding::SmoothUsingFunction() {
538 //
539 // Fits the unfolded spectrum using the function fSmoothFunction
540 //
541
542 AliInfo(Form("Smooth function is a %s with option \"%s\" and has %d parameters : ",fSmoothFunction->ClassName(),fSmoothOption,fSmoothFunction->GetNpar()));
543
544 for (Int_t iPar=0; iPar<fSmoothFunction->GetNpar(); iPar++) AliInfo(Form("par[%d]=%e",iPar,fSmoothFunction->GetParameter(iPar)));
545
546 switch (fNVariables) {
547 case 1 : fUnfolded->Projection(0) ->Fit(fSmoothFunction,fSmoothOption); break;
548 case 2 : fUnfolded->Projection(1,0) ->Fit(fSmoothFunction,fSmoothOption); break; // (1,0) instead of (0,1) -> TAxis issue
549 case 3 : fUnfolded->Projection(0,1,2)->Fit(fSmoothFunction,fSmoothOption); break;
550 default: AliError(Form("Cannot handle such fit in %d dimensions",fNVariables)) ; return 1;
551 }
552
553 Int_t nDim = fNVariables;
554 Int_t* bins = new Int_t[nDim]; // number of bins for each variable
555 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
556
557 for (Int_t iVar=0; iVar<nDim; iVar++) {
558 bins[iVar] = fUnfolded->GetAxis(iVar)->GetNbins();
559 nBins *= bins[iVar];
560 }
561
562 Int_t *bin = new Int_t[nDim]; // bin to fill the THnSparse (holding the bin coordinates)
563 Double_t x[3] = {0,0,0} ; // value in bin center (max dimension is 3 (TF3))
564
565 // loop on the bins and update of fUnfolded
566 // THnSparse::Multiply(TF1*) doesn't exist, so let's do it bin by bin
567 for (Long_t iBin=0; iBin<nBins; iBin++) {
568 Long_t bin_tmp = iBin ;
569 for (Int_t iVar=0; iVar<nDim; iVar++) {
570 bin[iVar] = 1 + bin_tmp % bins[iVar] ;
571 bin_tmp /= bins[iVar] ;
572 x[iVar] = fUnfolded->GetAxis(iVar)->GetBinCenter(bin[iVar]);
573 }
574 Double_t functionValue = fSmoothFunction->Eval(x[0],x[1],x[2]) ;
575 fUnfolded->SetBinContent(bin,functionValue);
576 fUnfolded->SetBinError (bin,functionValue*fUnfolded->GetBinError(bin));
577 }
578 return 0;
579}
c0b10ad4 580
581//______________________________________________________________
582
583void AliCFUnfolding::CreateFlatPrior() {
584 //
585 // Creates a flat prior distribution
586 //
587
588 AliInfo("Creating a flat a priori distribution");
589
590 // create the frame of the THnSparse given (for example) the one from the efficiency map
591 fPrior = (THnSparse*) fEfficiency->Clone();
592
85b6bda9 593 if (fNVariables != fPrior->GetNdimensions())
594 AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions()));
595
c0b10ad4 596 Int_t nDim = fNVariables;
597 Int_t* bins = new Int_t[nDim]; // number of bins for each variable
598 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
599
600 for (Int_t iVar=0; iVar<nDim; iVar++) {
601 bins[iVar] = fPrior->GetAxis(iVar)->GetNbins();
602 nBins *= bins[iVar];
603 }
604
605 Int_t *bin = new Int_t[nDim]; // bin to fill the THnSparse (holding the bin coordinates)
606
607 // loop that sets 1 in each bin
608 for (Long_t iBin=0; iBin<nBins; iBin++) {
609 Long_t bin_tmp = iBin ;
610 for (Int_t iVar=0; iVar<nDim; iVar++) {
611 bin[iVar] = 1 + bin_tmp % bins[iVar] ;
612 bin_tmp /= bins[iVar] ;
613 }
614 fPrior->SetBinContent(bin,1.); // put 1 everywhere
85b6bda9 615 fPrior->SetBinError (bin,0.); // put 0 everywhere
c0b10ad4 616 }
617
618 fOriginalPrior = (THnSparse*)fPrior->Clone();
619
620 delete [] bin;
621 delete [] bins;
622}