New class for dE/dx analysis (comparison with Bethe-Bloch, check of response function...
[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
7bc6a013 275 for (Long_t i=0; i<fMeasuredEstimate->GetNbins(); i++) {
c0b10ad4 276 fMeasuredEstimate->GetBinContent(i,fCoordinatesN_M);
277 fMeasuredEstimate->SetBinContent(fCoordinatesN_M,0.);
85b6bda9 278 fMeasuredEstimate->SetBinError (fCoordinatesN_M,0.);
c0b10ad4 279 }
280
85b6bda9 281 THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone();
282 priorTimesEff->Multiply(fEfficiency);
283
c0b10ad4 284 // fill it
7bc6a013 285 for (Long_t iBin=0; iBin<fConditional->GetNbins(); iBin++) {
c0b10ad4 286 Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N);
85b6bda9 287 Double_t conditionalError = fConditional->GetBinError (iBin);
c0b10ad4 288 GetCoordinates();
85b6bda9 289 Double_t priorTimesEffValue = priorTimesEff->GetBinContent(fCoordinatesN_T);
290 Double_t priorTimesEffError = priorTimesEff->GetBinError (fCoordinatesN_T);
291 Double_t fill = conditionalValue * priorTimesEffValue ;
292
293 if (fill>0.) {
294 fMeasuredEstimate->AddBinContent(fCoordinatesN_M,fill);
295
296 // error calculation : gaussian error propagation (may be overestimated...)
297 Double_t err2 = TMath::Power(fMeasuredEstimate->GetBinError(fCoordinatesN_M),2) ;
298 err2 += TMath::Power(conditionalValue*priorTimesEffError,2) + TMath::Power(conditionalError*priorTimesEffValue,2) ;
299 Double_t err = TMath::Sqrt(err2);
300 fMeasuredEstimate->SetBinError(fCoordinatesN_M,err);
301 }
c0b10ad4 302 }
85b6bda9 303 delete priorTimesEff ;
c0b10ad4 304}
305
306//______________________________________________________________
307
308void AliCFUnfolding::CreateInvResponse() {
309 //
310 // Creates the inverse response matrix (INV) with Bayesian method
311 // : uses the conditional matrix (COND), the prior probabilities (T) and the efficiency map (E)
312 //
313 // --> P(T|M) = P(M|T) * P(T) * eff(T) / SUM { P(M|T) * P(T) }
314 // --> INV(i,j) = COND(i,j) * T(j) * E(j) / SUM_k { COND(i,k) * T(k) }
315 //
316
85b6bda9 317 THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone();
318 priorTimesEff->Multiply(fEfficiency);
319
7bc6a013 320 for (Long_t iBin=0; iBin<fConditional->GetNbins(); iBin++) {
c0b10ad4 321 Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N);
85b6bda9 322 Double_t conditionalError = fConditional->GetBinError (iBin);
c0b10ad4 323 GetCoordinates();
85b6bda9 324 Double_t estMeasuredValue = fMeasuredEstimate->GetBinContent(fCoordinatesN_M);
325 Double_t estMeasuredError = fMeasuredEstimate->GetBinError (fCoordinatesN_M);
326 Double_t priorTimesEffValue = priorTimesEff ->GetBinContent(fCoordinatesN_T);
327 Double_t priorTimesEffError = priorTimesEff ->GetBinError (fCoordinatesN_T);
328 Double_t fill = (estMeasuredValue>0. ? conditionalValue * priorTimesEffValue / estMeasuredValue : 0. ) ;
329 // error calculation : gaussian error propagation (may be overestimated...)
330 Double_t err = 0. ;
331 if (estMeasuredValue>0.) {
332 err = TMath::Sqrt( TMath::Power(conditionalError * priorTimesEffValue * estMeasuredValue ,2) +
333 TMath::Power(conditionalValue * priorTimesEffError * estMeasuredValue ,2) +
334 TMath::Power(conditionalValue * priorTimesEffValue * estMeasuredError ,2) )
335 / TMath::Power(estMeasuredValue,2) ;
336 }
337 if (fill>0. || fInverseResponse->GetBinContent(fCoordinates2N)>0.) {
338 fInverseResponse->SetBinContent(fCoordinates2N,fill);
339 fInverseResponse->SetBinError (fCoordinates2N,err );
340 }
341 }
342 delete priorTimesEff ;
c0b10ad4 343}
344
345//______________________________________________________________
346
347void AliCFUnfolding::Unfold() {
348 //
349 // Main routine called by the user :
350 // it calculates the unfolded spectrum from the response matrix and the measured spectrum
351 // several iterations are performed until a reasonable chi2 is reached
352 //
353
354 Int_t iIterBayes=0 ;
355 Double_t chi2=0 ;
356
357 for (iIterBayes=0; iIterBayes<fMaxNumIterations; iIterBayes++) { // bayes iterations
358 CreateEstMeasured();
359 CreateInvResponse();
360 CreateUnfolded();
361 chi2 = GetChi2();
7bc6a013 362 AliDebug(1,Form("Chi2 at iteration %d is %e",iIterBayes,chi2));
c0b10ad4 363 if (fMaxChi2>0. && chi2<fMaxChi2) {
364 break;
365 }
366 // update the prior distribution
85b6bda9 367 if (fUseSmoothing) {
368 if (Smooth()) {
369 AliError("Couldn't smooth the unfolded spectrum!!");
7bc6a013 370 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 371 return;
372 }
373 }
c0b10ad4 374 fPrior = (THnSparse*)fUnfolded->Clone() ; // this should be changed (memory)
375 }
85b6bda9 376 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 377}
378
379//______________________________________________________________
380
381void AliCFUnfolding::CreateUnfolded() {
382 //
383 // Creates the unfolded (T) spectrum from the measured spectrum (M) and the inverse response matrix (INV)
384 // We have P(T) = SUM { P(T|M) * P(M) }
385 // --> T(i) = SUM_k { INV(i,k) * M(k) }
386 //
387
388
389 // clear the unfolded spectrum
7bc6a013 390 for (Long_t i=0; i<fUnfolded->GetNbins(); i++) {
c0b10ad4 391 fUnfolded->GetBinContent(i,fCoordinatesN_T);
392 fUnfolded->SetBinContent(fCoordinatesN_T,0.);
85b6bda9 393 fUnfolded->SetBinError (fCoordinatesN_T,0.);
c0b10ad4 394 }
395
7bc6a013 396 for (Long_t iBin=0; iBin<fInverseResponse->GetNbins(); iBin++) {
c0b10ad4 397 Double_t invResponseValue = fInverseResponse->GetBinContent(iBin,fCoordinates2N);
85b6bda9 398 Double_t invResponseError = fInverseResponse->GetBinError (iBin);
c0b10ad4 399 GetCoordinates();
85b6bda9 400 Double_t effValue = fEfficiency->GetBinContent(fCoordinatesN_T);
401 Double_t effError = fEfficiency->GetBinError (fCoordinatesN_T);
402 Double_t measuredValue = fMeasured ->GetBinContent(fCoordinatesN_M);
403 Double_t measuredError = fMeasured ->GetBinError (fCoordinatesN_M);
404 Double_t fill = (effValue>0. ? invResponseValue * measuredValue / effValue : 0.) ;
405
406 if (fill>0.) {
407 fUnfolded->AddBinContent(fCoordinatesN_T,fill);
408
409 // error calculation : gaussian error propagation (may be overestimated...)
410 Double_t err2 = TMath::Power(fUnfolded->GetBinError(fCoordinatesN_T),2) ;
411 err2 += TMath::Power(invResponseError * measuredValue * effValue,2) / TMath::Power(effValue,4) ;
412 err2 += TMath::Power(invResponseValue * measuredError * effValue,2) / TMath::Power(effValue,4) ;
413 err2 += TMath::Power(invResponseValue * measuredValue * effError,2) / TMath::Power(effValue,4) ;
414 Double_t err = TMath::Sqrt(err2);
415 fUnfolded->SetBinError(fCoordinatesN_T,err);
416 }
c0b10ad4 417 }
418}
419
85b6bda9 420//______________________________________________________________
421
c0b10ad4 422void AliCFUnfolding::GetCoordinates() {
423 //
424 // assign coordinates in Measured and True spaces (dim=N) from coordinates in global space (dim=2N)
425 //
426 for (Int_t i = 0; i<fNVariables ; i++) {
427 fCoordinatesN_M[i] = fCoordinates2N[i];
428 fCoordinatesN_T[i] = fCoordinates2N[i+fNVariables];
429 }
430}
431
432//______________________________________________________________
433
434void AliCFUnfolding::CreateConditional() {
435 //
436 // creates the conditional probability matrix (R*) holding the P(M|T), given the reponse matrix R
437 //
438 // --> R*(i,j) = R(i,j) / SUM_k{ R(k,j) }
439 //
440
441 fConditional = (THnSparse*) fResponse->Clone(); // output of this function
442 fProjResponseInT = (THnSparse*) fPrior->Clone(); // output denominator :
443 // projection of the response matrix on the TRUE axis
85b6bda9 444 Int_t* dim = new Int_t [fNVariables];
445 for (Int_t iDim=0; iDim<fNVariables; iDim++) dim[iDim] = fNVariables+iDim ; //dimensions corresponding to TRUE values (i.e. from N to 2N-1)
446 fProjResponseInT = fConditional->Projection(fNVariables,dim,"E"); //project
447 delete [] dim;
c0b10ad4 448
449 // fill the conditional probability matrix
7bc6a013 450 for (Long_t iBin=0; iBin<fResponse->GetNbins(); iBin++) {
c0b10ad4 451 Double_t responseValue = fResponse->GetBinContent(iBin,fCoordinates2N);
85b6bda9 452 Double_t responseError = fResponse->GetBinError (iBin);
c0b10ad4 453 GetCoordinates();
85b6bda9 454 Double_t projValue = fProjResponseInT->GetBinContent(fCoordinatesN_T);
455 Double_t projError = fProjResponseInT->GetBinError (fCoordinatesN_T);
456
457 Double_t fill = responseValue / projValue ;
458 if (fill>0. || fConditional->GetBinContent(fCoordinates2N)>0.) {
459 fConditional->SetBinContent(fCoordinates2N,fill);
460 // gaussian error for the moment
461 Double_t err2 = TMath::Power(responseError*projValue,2) + TMath::Power(responseValue*projError,2) ;
462 Double_t err = TMath::Sqrt(err2);
463 err /= TMath::Power(projValue,2) ;
464 fConditional->SetBinError (fCoordinates2N,err);
465 }
c0b10ad4 466 }
467}
468
469//______________________________________________________________
470
471Double_t AliCFUnfolding::GetChi2() {
472 //
473 // Returns the chi2 between unfolded and a priori spectrum
474 //
475
476 Double_t chi2 = 0. ;
7bc6a013 477 for (Long_t iBin=0; iBin<fPrior->GetNbins(); iBin++) {
c0b10ad4 478 Double_t priorValue = fPrior->GetBinContent(iBin);
7bc6a013 479// chi2 += (priorValue>0. ? TMath::Power(fUnfolded->GetBinContent(iBin) - priorValue,2) / priorValue : 0.) ;
480 chi2 += (fUnfolded->GetBinError(iBin)>0. ? TMath::Power((fUnfolded->GetBinContent(iBin) - priorValue)/fUnfolded->GetBinError(iBin),2) / priorValue : 0.) ;
c0b10ad4 481 }
482 return chi2;
483}
484
485//______________________________________________________________
486
487void AliCFUnfolding::SetMaxChi2PerDOF(Double_t val) {
488 //
489 // Max. chi2 per degree of freedom : user setting
490 //
491
492 Int_t nDOF = 1 ;
493 for (Int_t iDim=0; iDim<fNVariables; iDim++) {
494 nDOF *= fPrior->GetAxis(iDim)->GetNbins();
495 }
496 AliInfo(Form("Number of degrees of freedom = %d",nDOF));
497 fMaxChi2 = val * nDOF ;
498}
499
500//______________________________________________________________
501
85b6bda9 502Short_t AliCFUnfolding::Smooth() {
c0b10ad4 503 //
504 // Smoothes the unfolded spectrum
85b6bda9 505 //
506 // By default each cell content is replaced by the average with the neighbouring bins (but not diagonally-neighbouring bins)
507 // However, if a specific function fcn has been defined in UseSmoothing(fcn), the unfolded will be fit and updated using fcn
c0b10ad4 508 //
509
85b6bda9 510 if (fSmoothFunction) {
7bc6a013 511 AliDebug(2,Form("Smoothing spectrum with fit function %p",fSmoothFunction));
85b6bda9 512 return SmoothUsingFunction();
513 }
514 else return SmoothUsingNeighbours();
515}
516
517//______________________________________________________________
518
519Short_t AliCFUnfolding::SmoothUsingNeighbours() {
520 //
521 // Smoothes the unfolded spectrum using neighouring bins
522 //
523
c0b10ad4 524 Int_t* numBins = new Int_t[fNVariables];
525 for (Int_t iVar=0; iVar<fNVariables; iVar++) numBins[iVar]=fUnfolded->GetAxis(iVar)->GetNbins();
526
527 //need a copy because fUnfolded will be updated during the loop, and this creates problems
528 THnSparse* copy = (THnSparse*)fUnfolded->Clone();
529
7bc6a013 530 for (Long_t iBin=0; iBin<copy->GetNbins(); iBin++) { //loop on non-empty bins
c0b10ad4 531 Double_t content = copy->GetBinContent(iBin,fCoordinatesN_T);
7bc6a013 532 Double_t error2 = TMath::Power(copy->GetBinError(iBin),2);
c0b10ad4 533
534 // skip the under/overflow bins...
535 Bool_t isOutside = kFALSE ;
536 for (Int_t iVar=0; iVar<fNVariables; iVar++) {
537 if (fCoordinatesN_T[iVar]<1 || fCoordinatesN_T[iVar]>numBins[iVar]) {
538 isOutside=kTRUE;
539 break;
540 }
541 }
542 if (isOutside) continue;
543
544 Int_t neighbours = 0; // number of neighbours to average with
545
546 for (Int_t iVar=0; iVar<fNVariables; iVar++) {
547 if (fCoordinatesN_T[iVar] > 1) { // must not be on low edge border
548 fCoordinatesN_T[iVar]-- ; //get lower neighbouring bin
85b6bda9 549 content += copy->GetBinContent(fCoordinatesN_T);
550 error2 += TMath::Power(copy->GetBinError(fCoordinatesN_T),2);
c0b10ad4 551 neighbours++;
552 fCoordinatesN_T[iVar]++ ; //back to initial coordinate
553 }
554 if (fCoordinatesN_T[iVar] < numBins[iVar]) { // must not be on up edge border
555 fCoordinatesN_T[iVar]++ ; //get upper neighbouring bin
85b6bda9 556 content += copy->GetBinContent(fCoordinatesN_T);
557 error2 += TMath::Power(copy->GetBinError(fCoordinatesN_T),2);
c0b10ad4 558 neighbours++;
559 fCoordinatesN_T[iVar]-- ; //back to initial coordinate
560 }
561 }
85b6bda9 562 // make an average
563 fUnfolded->SetBinContent(fCoordinatesN_T,content/(1.+neighbours));
564 fUnfolded->SetBinError (fCoordinatesN_T,TMath::Sqrt(error2)/(1.+neighbours));
c0b10ad4 565 }
566 delete [] numBins;
567 delete copy;
85b6bda9 568 return 0;
c0b10ad4 569}
570
85b6bda9 571//______________________________________________________________
572
573Short_t AliCFUnfolding::SmoothUsingFunction() {
574 //
575 // Fits the unfolded spectrum using the function fSmoothFunction
576 //
577
7bc6a013 578 AliDebug(0,Form("Smooth function is a %s with option \"%s\" and has %d parameters : ",fSmoothFunction->ClassName(),fSmoothOption,fSmoothFunction->GetNpar()));
85b6bda9 579
7bc6a013 580 for (Int_t iPar=0; iPar<fSmoothFunction->GetNpar(); iPar++) AliDebug(0,Form("par[%d]=%e",iPar,fSmoothFunction->GetParameter(iPar)));
581
582 Int_t fitResult = 0;
85b6bda9 583
584 switch (fNVariables) {
7bc6a013 585 case 1 : fitResult = fUnfolded->Projection(0) ->Fit(fSmoothFunction,fSmoothOption); break;
586 case 2 : fitResult = fUnfolded->Projection(1,0) ->Fit(fSmoothFunction,fSmoothOption); break; // (1,0) instead of (0,1) -> TAxis issue
587 case 3 : fitResult = fUnfolded->Projection(0,1,2)->Fit(fSmoothFunction,fSmoothOption); break;
588 default: AliFatal(Form("Cannot handle such fit in %d dimensions",fNVariables)) ; return 1;
589 }
590
591 if (fitResult != 0) {
592 AliWarning(Form("Fit failed with status %d, stopping the loop",fitResult));
593 return 1;
85b6bda9 594 }
595
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] = fUnfolded->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 Double_t x[3] = {0,0,0} ; // value in bin center (max dimension is 3 (TF3))
607
608 // loop on the bins and update of fUnfolded
609 // THnSparse::Multiply(TF1*) doesn't exist, so let's do it bin by bin
610 for (Long_t iBin=0; iBin<nBins; iBin++) {
611 Long_t bin_tmp = iBin ;
612 for (Int_t iVar=0; iVar<nDim; iVar++) {
613 bin[iVar] = 1 + bin_tmp % bins[iVar] ;
614 bin_tmp /= bins[iVar] ;
615 x[iVar] = fUnfolded->GetAxis(iVar)->GetBinCenter(bin[iVar]);
616 }
617 Double_t functionValue = fSmoothFunction->Eval(x[0],x[1],x[2]) ;
618 fUnfolded->SetBinContent(bin,functionValue);
619 fUnfolded->SetBinError (bin,functionValue*fUnfolded->GetBinError(bin));
620 }
621 return 0;
622}
c0b10ad4 623
624//______________________________________________________________
625
626void AliCFUnfolding::CreateFlatPrior() {
627 //
628 // Creates a flat prior distribution
629 //
630
631 AliInfo("Creating a flat a priori distribution");
632
633 // create the frame of the THnSparse given (for example) the one from the efficiency map
634 fPrior = (THnSparse*) fEfficiency->Clone();
635
85b6bda9 636 if (fNVariables != fPrior->GetNdimensions())
637 AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions()));
638
c0b10ad4 639 Int_t nDim = fNVariables;
640 Int_t* bins = new Int_t[nDim]; // number of bins for each variable
641 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
642
643 for (Int_t iVar=0; iVar<nDim; iVar++) {
644 bins[iVar] = fPrior->GetAxis(iVar)->GetNbins();
645 nBins *= bins[iVar];
646 }
647
648 Int_t *bin = new Int_t[nDim]; // bin to fill the THnSparse (holding the bin coordinates)
649
650 // loop that sets 1 in each bin
651 for (Long_t iBin=0; iBin<nBins; iBin++) {
652 Long_t bin_tmp = iBin ;
653 for (Int_t iVar=0; iVar<nDim; iVar++) {
654 bin[iVar] = 1 + bin_tmp % bins[iVar] ;
655 bin_tmp /= bins[iVar] ;
656 }
657 fPrior->SetBinContent(bin,1.); // put 1 everywhere
85b6bda9 658 fPrior->SetBinError (bin,0.); // put 0 everywhere
c0b10ad4 659 }
660
661 fOriginalPrior = (THnSparse*)fPrior->Clone();
662
663 delete [] bin;
664 delete [] bins;
665}