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