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