]> git.uio.no Git - u/mrichter/AliRoot.git/blame - CORRFW/AliCFUnfolding.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
coveritiy
[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 //
769f5114 41// and the unfolded spectrum (OBSOLETE). //
42// //
43// Chi2 calculation became absolute with the correlated error //
44// calculation. //
45// Errors on the unfolded distribution are not known until the end //
46// Use the convergence criterion instead //
7bc6a013 47// //
48// Currently the user has to define the max. number of iterations //
49// (::SetMaxNumberOfIterations) //
769f5114 50// and //
51// - the chi2 below which the procedure will stop //
52// (::SetMaxChi2 or ::SetMaxChi2PerDOF) (OBSOLETE) //
53// - the convergence criterion below which the procedure will stop //
54// SetMaxConvergencePerDOF(Double_t val); //
55// //
56// Correlated error calculation can be activated by using: //
57// SetUseCorrelatedErrors(Bool_t b) in combination with convergence //
58// criterion //
59// Documentation about correlated error calculation method can be //
60// found in AliCFUnfolding::CalculateCorrelatedErrors() //
61// Author: marta.verweij@cern.ch //
7bc6a013 62// //
63// An optional possibility is to smooth the unfolded spectrum at the //
64// end of each iteration, either using a fit function //
65// (only if #dimensions <=3) //
66// or a simple averaging using the neighbouring bins values. //
67// This is possible calling the function ::UseSmoothing //
68// If no argument is passed to this function, then the second option //
69// is used. //
70// //
c7f196bd 71// IMPORTANT: //
72//----------- //
73// With this approach, the efficiency map must be calculated //
74// with *simulated* values only, otherwise the method won't work. //
75// //
76// ex: efficiency(bin_pt) = number_rec(bin_pt) / number_sim(bin_pt) //
77// //
78// the pt bin "bin_pt" must always be the same in both the efficiency //
79// numerator and denominator. //
80// This is why the efficiency map has to be created by a method //
81// from which both reconstructed and simulated values are accessible //
82// simultaneously. //
83// //
84// //
7bc6a013 85//---------------------------------------------------------------------//
86// Author : renaud.vernet@cern.ch //
87//---------------------------------------------------------------------//
c0b10ad4 88
89
90#include "AliCFUnfolding.h"
91#include "TMath.h"
92#include "TAxis.h"
93#include "AliLog.h"
85b6bda9 94#include "TF1.h"
95#include "TH1D.h"
96#include "TH2D.h"
97#include "TH3D.h"
769f5114 98#include "TProfile.h"
99#include "TRandom3.h"
c0b10ad4 100
101ClassImp(AliCFUnfolding)
102
103//______________________________________________________________
104
105AliCFUnfolding::AliCFUnfolding() :
106 TNamed(),
107 fResponse(0x0),
108 fPrior(0x0),
c0b10ad4 109 fEfficiency(0x0),
110 fMeasured(0x0),
769f5114 111 fMeasuredOrig(0x0),
112 fMaxNumIterations(20),
c0b10ad4 113 fNVariables(0),
c0b10ad4 114 fUseSmoothing(kFALSE),
85b6bda9 115 fSmoothFunction(0x0),
116 fSmoothOption(""),
769f5114 117 fMaxConvergence(0),
118 fUseCorrelatedErrors(kTRUE),
119 fNRandomIterations(20),
85b6bda9 120 fOriginalPrior(0x0),
c0b10ad4 121 fInverseResponse(0x0),
122 fMeasuredEstimate(0x0),
123 fConditional(0x0),
124 fProjResponseInT(0x0),
125 fUnfolded(0x0),
126 fCoordinates2N(0x0),
127 fCoordinatesN_M(0x0),
769f5114 128 fCoordinatesN_T(0x0),
129 fRandomizedDist(0x0),
130 fRandom3(0x0),
131 fRandomUnfolded(0x0),
132 fDeltaUnfoldedP(0x0),
133 fNCalcCorrErrors(0),
134 fRandomSeed(0)
c0b10ad4 135{
136 //
137 // default constructor
138 //
139}
140
141//______________________________________________________________
142
143AliCFUnfolding::AliCFUnfolding(const Char_t* name, const Char_t* title, const Int_t nVar,
144 const THnSparse* response, const THnSparse* efficiency, const THnSparse* measured, const THnSparse* prior) :
145 TNamed(name,title),
146 fResponse((THnSparse*)response->Clone()),
147 fPrior(0x0),
c0b10ad4 148 fEfficiency((THnSparse*)efficiency->Clone()),
149 fMeasured((THnSparse*)measured->Clone()),
769f5114 150 fMeasuredOrig((THnSparse*)measured->Clone()),
c0b10ad4 151 fMaxNumIterations(0),
152 fNVariables(nVar),
c0b10ad4 153 fUseSmoothing(kFALSE),
85b6bda9 154 fSmoothFunction(0x0),
155 fSmoothOption(""),
769f5114 156 fMaxConvergence(0),
157 fUseCorrelatedErrors(kTRUE),
158 fNRandomIterations(20),
85b6bda9 159 fOriginalPrior(0x0),
c0b10ad4 160 fInverseResponse(0x0),
161 fMeasuredEstimate(0x0),
162 fConditional(0x0),
163 fProjResponseInT(0x0),
164 fUnfolded(0x0),
165 fCoordinates2N(0x0),
166 fCoordinatesN_M(0x0),
769f5114 167 fCoordinatesN_T(0x0),
168 fRandomizedDist(0x0),
169 fRandom3(0x0),
170 fRandomUnfolded(0x0),
171 fDeltaUnfoldedP(0x0),
172 fNCalcCorrErrors(0),
173 fRandomSeed(0)
c0b10ad4 174{
175 //
176 // named constructor
177 //
178
179 AliInfo(Form("\n\n--------------------------\nCreating an unfolder :\n--------------------------\nresponse matrix has %d dimension(s)",fResponse->GetNdimensions()));
180
181 if (!prior) CreateFlatPrior(); // if no prior distribution declared, simply use a flat distribution
182 else {
183 fPrior = (THnSparse*) prior->Clone();
184 fOriginalPrior = (THnSparse*)fPrior->Clone();
85b6bda9 185 if (fPrior->GetNdimensions() != fNVariables)
186 AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions()));
c0b10ad4 187 }
85b6bda9 188
189 if (fEfficiency->GetNdimensions() != fNVariables)
190 AliFatal(Form("The efficiency matrix should have %d dimensions, and it has actually %d",fNVariables,fEfficiency->GetNdimensions()));
191 if (fMeasured->GetNdimensions() != fNVariables)
192 AliFatal(Form("The measured matrix should have %d dimensions, and it has actually %d",fNVariables,fMeasured->GetNdimensions()));
193 if (fResponse->GetNdimensions() != 2*fNVariables)
194 AliFatal(Form("The response matrix should have %d dimensions, and it has actually %d",2*fNVariables,fResponse->GetNdimensions()));
c0b10ad4 195
85b6bda9 196
c0b10ad4 197 for (Int_t iVar=0; iVar<fNVariables; iVar++) {
198 AliInfo(Form("prior matrix has %d bins in dimension %d",fPrior ->GetAxis(iVar)->GetNbins(),iVar));
199 AliInfo(Form("efficiency matrix has %d bins in dimension %d",fEfficiency->GetAxis(iVar)->GetNbins(),iVar));
200 AliInfo(Form("measured matrix has %d bins in dimension %d",fMeasured ->GetAxis(iVar)->GetNbins(),iVar));
201 }
769f5114 202
c0b10ad4 203 Init();
204}
205
206
207//______________________________________________________________
208
209AliCFUnfolding::AliCFUnfolding(const AliCFUnfolding& c) :
210 TNamed(c),
211 fResponse((THnSparse*)c.fResponse->Clone()),
212 fPrior((THnSparse*)c.fPrior->Clone()),
c0b10ad4 213 fEfficiency((THnSparse*)c.fEfficiency->Clone()),
214 fMeasured((THnSparse*)c.fMeasured->Clone()),
769f5114 215 fMeasuredOrig((THnSparse*)c.fMeasuredOrig->Clone()),
c0b10ad4 216 fMaxNumIterations(c.fMaxNumIterations),
217 fNVariables(c.fNVariables),
c0b10ad4 218 fUseSmoothing(c.fUseSmoothing),
85b6bda9 219 fSmoothFunction((TF1*)c.fSmoothFunction->Clone()),
7c5606f9 220 fSmoothOption(c.fSmoothOption),
769f5114 221 fMaxConvergence(c.fMaxConvergence),
222 fUseCorrelatedErrors(c.fUseCorrelatedErrors),
223 fNRandomIterations(c.fNRandomIterations),
85b6bda9 224 fOriginalPrior((THnSparse*)c.fOriginalPrior->Clone()),
c0b10ad4 225 fInverseResponse((THnSparse*)c.fInverseResponse->Clone()),
226 fMeasuredEstimate((THnSparse*)fMeasuredEstimate->Clone()),
227 fConditional((THnSparse*)c.fConditional->Clone()),
228 fProjResponseInT((THnSparse*)c.fProjResponseInT->Clone()),
229 fUnfolded((THnSparse*)c.fUnfolded->Clone()),
230 fCoordinates2N(new Int_t(*c.fCoordinates2N)),
231 fCoordinatesN_M(new Int_t(*c.fCoordinatesN_M)),
769f5114 232 fCoordinatesN_T(new Int_t(*c.fCoordinatesN_T)),
233 fRandomizedDist((THnSparse*)c.fRandomizedDist->Clone()),
234 fRandom3((TRandom3*)c.fRandom3->Clone()),
235 fRandomUnfolded((THnSparse*)c.fRandomUnfolded->Clone()),
236 fDeltaUnfoldedP((TProfile*)c.fDeltaUnfoldedP),
237 fNCalcCorrErrors(c.fNCalcCorrErrors),
238 fRandomSeed(c.fRandomSeed)
c0b10ad4 239{
240 //
241 // copy constructor
242 //
243}
244
245//______________________________________________________________
246
247AliCFUnfolding& AliCFUnfolding::operator=(const AliCFUnfolding& c) {
248 //
249 // assignment operator
250 //
251
252 if (this!=&c) {
253 TNamed::operator=(c);
254 fResponse = (THnSparse*)c.fResponse->Clone() ;
255 fPrior = (THnSparse*)c.fPrior->Clone() ;
c0b10ad4 256 fEfficiency = (THnSparse*)c.fEfficiency->Clone() ;
257 fMeasured = (THnSparse*)c.fMeasured->Clone() ;
769f5114 258 fMeasuredOrig = ((THnSparse*)c.fMeasuredOrig->Clone()),
c0b10ad4 259 fMaxNumIterations = c.fMaxNumIterations ;
260 fNVariables = c.fNVariables ;
769f5114 261 fMaxConvergence = c.fMaxConvergence ;
c0b10ad4 262 fUseSmoothing = c.fUseSmoothing ;
85b6bda9 263 fSmoothFunction = (TF1*)c.fSmoothFunction->Clone();
264 fSmoothOption = c.fSmoothOption ;
769f5114 265 fUseCorrelatedErrors = c.fUseCorrelatedErrors;
266 fNRandomIterations = c.fNRandomIterations;
85b6bda9 267 fOriginalPrior = (THnSparse*)c.fOriginalPrior->Clone() ;
c0b10ad4 268 fInverseResponse = (THnSparse*)c.fInverseResponse->Clone() ;
269 fMeasuredEstimate = (THnSparse*)fMeasuredEstimate->Clone() ;
270 fConditional = (THnSparse*)c.fConditional->Clone() ;
271 fProjResponseInT = (THnSparse*)c.fProjResponseInT->Clone() ;
272 fUnfolded = (THnSparse*)c.fUnfolded->Clone() ;
273 fCoordinates2N = new Int_t(*c.fCoordinates2N) ;
274 fCoordinatesN_M = new Int_t(*c.fCoordinatesN_M) ;
275 fCoordinatesN_T = new Int_t(*c.fCoordinatesN_T) ;
769f5114 276 fRandomizedDist = (THnSparse*)c.fRandomizedDist->Clone();
277 fRandom3 = (TRandom3*)c.fRandom3->Clone();
278 fRandomUnfolded = (THnSparse*)c.fRandomUnfolded->Clone();
279 fDeltaUnfoldedP = (TProfile*)c.fDeltaUnfoldedP;
280 fNCalcCorrErrors = c.fNCalcCorrErrors;
281 fRandomSeed = c.fRandomSeed ;
c0b10ad4 282 }
283 return *this;
284}
285
286//______________________________________________________________
287
288AliCFUnfolding::~AliCFUnfolding() {
289 //
290 // destructor
291 //
292 if (fResponse) delete fResponse;
293 if (fPrior) delete fPrior;
c0b10ad4 294 if (fEfficiency) delete fEfficiency;
295 if (fMeasured) delete fMeasured;
769f5114 296 if (fMeasuredOrig) delete fMeasuredOrig;
85b6bda9 297 if (fSmoothFunction) delete fSmoothFunction;
298 if (fOriginalPrior) delete fOriginalPrior;
c0b10ad4 299 if (fInverseResponse) delete fInverseResponse;
300 if (fMeasuredEstimate) delete fMeasuredEstimate;
301 if (fConditional) delete fConditional;
302 if (fProjResponseInT) delete fProjResponseInT;
303 if (fCoordinates2N) delete [] fCoordinates2N;
304 if (fCoordinatesN_M) delete [] fCoordinatesN_M;
305 if (fCoordinatesN_T) delete [] fCoordinatesN_T;
769f5114 306 if (fRandomizedDist) delete fRandomizedDist;
307 if (fRandom3) delete fRandom3;
308 if (fRandomUnfolded) delete fRandomUnfolded;
309 if (fDeltaUnfoldedP) delete fDeltaUnfoldedP;
c0b10ad4 310}
311
312//______________________________________________________________
313
314void AliCFUnfolding::Init() {
315 //
316 // initialisation function : creates internal settings
317 //
318
769f5114 319 fRandom3 = new TRandom3(fRandomSeed);
320
c0b10ad4 321 fCoordinates2N = new Int_t[2*fNVariables];
322 fCoordinatesN_M = new Int_t[fNVariables];
323 fCoordinatesN_T = new Int_t[fNVariables];
324
325 // create the matrix of conditional probabilities P(M|T)
326 CreateConditional();
327
328 // create the frame of the inverse response matrix
329 fInverseResponse = (THnSparse*) fResponse->Clone();
330 // create the frame of the unfolded spectrum
331 fUnfolded = (THnSparse*) fPrior->Clone();
769f5114 332 // create the frame of the random unfolded spectrum
333 fRandomUnfolded = (THnSparse*) fPrior->Clone();
c0b10ad4 334 // create the frame of the measurement estimate spectrum
335 fMeasuredEstimate = (THnSparse*) fMeasured->Clone();
769f5114 336 // create the frame of the original measurement spectrum
337 fMeasuredOrig = (THnSparse*) fMeasured->Clone();
338
339 InitDeltaUnfoldedProfile();
340
c0b10ad4 341}
342
343//______________________________________________________________
344
345void AliCFUnfolding::CreateEstMeasured() {
346 //
347 // This function creates a estimate (M) of the reconstructed spectrum
7bc6a013 348 // given the a priori distribution (T), the efficiency (E) and the conditional matrix (COND)
c0b10ad4 349 //
350 // --> P(M) = SUM { P(M|T) * P(T) }
7bc6a013 351 // --> M(i) = SUM_k { COND(i,k) * T(k) * E (k)}
c0b10ad4 352 //
353 // This is needed to calculate the inverse response matrix
354 //
355
356
357 // clean the measured estimate spectrum
85f9f9e1 358 fMeasuredEstimate->Reset();
359
85b6bda9 360 THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone();
361 priorTimesEff->Multiply(fEfficiency);
362
c0b10ad4 363 // fill it
7bc6a013 364 for (Long_t iBin=0; iBin<fConditional->GetNbins(); iBin++) {
c0b10ad4 365 Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N);
366 GetCoordinates();
85b6bda9 367 Double_t priorTimesEffValue = priorTimesEff->GetBinContent(fCoordinatesN_T);
85b6bda9 368 Double_t fill = conditionalValue * priorTimesEffValue ;
369
370 if (fill>0.) {
12e419d5 371 fMeasuredEstimate->AddBinContent(fCoordinatesN_M,fill);
769f5114 372 fMeasuredEstimate->SetBinError(fCoordinatesN_M,0.);
85b6bda9 373 }
c0b10ad4 374 }
85b6bda9 375 delete priorTimesEff ;
c0b10ad4 376}
377
378//______________________________________________________________
379
380void AliCFUnfolding::CreateInvResponse() {
381 //
382 // Creates the inverse response matrix (INV) with Bayesian method
383 // : uses the conditional matrix (COND), the prior probabilities (T) and the efficiency map (E)
384 //
385 // --> P(T|M) = P(M|T) * P(T) * eff(T) / SUM { P(M|T) * P(T) }
386 // --> INV(i,j) = COND(i,j) * T(j) * E(j) / SUM_k { COND(i,k) * T(k) }
387 //
388
85b6bda9 389 THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone();
390 priorTimesEff->Multiply(fEfficiency);
391
7bc6a013 392 for (Long_t iBin=0; iBin<fConditional->GetNbins(); iBin++) {
c0b10ad4 393 Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N);
394 GetCoordinates();
85b6bda9 395 Double_t estMeasuredValue = fMeasuredEstimate->GetBinContent(fCoordinatesN_M);
85b6bda9 396 Double_t priorTimesEffValue = priorTimesEff ->GetBinContent(fCoordinatesN_T);
85b6bda9 397 Double_t fill = (estMeasuredValue>0. ? conditionalValue * priorTimesEffValue / estMeasuredValue : 0. ) ;
85b6bda9 398 if (fill>0. || fInverseResponse->GetBinContent(fCoordinates2N)>0.) {
399 fInverseResponse->SetBinContent(fCoordinates2N,fill);
769f5114 400 fInverseResponse->SetBinError (fCoordinates2N,0.);
85b6bda9 401 }
402 }
403 delete priorTimesEff ;
c0b10ad4 404}
405
406//______________________________________________________________
407
408void AliCFUnfolding::Unfold() {
409 //
410 // Main routine called by the user :
769f5114 411 // it calculates the unfolded spectrum from the response matrix, measured spectrum and efficiency
412 // several iterations are performed until a reasonable chi2 or convergence criterion is reached
c0b10ad4 413 //
414
769f5114 415 Int_t iIterBayes = 0 ;
416 Double_t convergence = 0.;
c0b10ad4 417
418 for (iIterBayes=0; iIterBayes<fMaxNumIterations; iIterBayes++) { // bayes iterations
419 CreateEstMeasured();
420 CreateInvResponse();
421 CreateUnfolded();
769f5114 422 if (fUseCorrelatedErrors) {
423 convergence = GetConvergence();
424 AliDebug(0,Form("convergence at iteration %d is %e",iIterBayes,convergence));
425 }
426 else AliWarning("No errors will be calculated. Activate SetUseCorrelatedErrors(kTRUE)\n");
427
428 if (fMaxConvergence>0. && convergence<fMaxConvergence && fNCalcCorrErrors<1) {
429 fNRandomIterations=iIterBayes;
c0b10ad4 430 break;
431 }
769f5114 432
c0b10ad4 433 // update the prior distribution
85b6bda9 434 if (fUseSmoothing) {
435 if (Smooth()) {
436 AliError("Couldn't smooth the unfolded spectrum!!");
769f5114 437 if (fUseCorrelatedErrors) {
438 if (fNCalcCorrErrors>0) {
439 AliInfo(Form("=======================\nUnfolding of randomized distribution finished at iteration %d with convergence %e \n",iIterBayes,convergence));
440 }
441 else {
442 AliInfo(Form("\n\n=======================\nFinished at iteration %d : convergence is %e and you required it to be < %e\n=======================\n\n",iIterBayes,convergence,fMaxConvergence));
443 }
444 }
85b6bda9 445 return;
446 }
447 }
769f5114 448 if(fNCalcCorrErrors>0) {
449 if (fPrior) delete fPrior ;
450 fPrior = (THnSparse*)fRandomUnfolded->Clone() ;
451 }
452 else {
453 if (fPrior) delete fPrior ;
454 fPrior = (THnSparse*)fUnfolded->Clone() ;
455 }
456 }
457
458 if (fUseCorrelatedErrors && fNCalcCorrErrors==0) {
459 fNCalcCorrErrors=1;
460 CalculateCorrelatedErrors();
461 }
462
463 if (fUseCorrelatedErrors) {
464 if (fNCalcCorrErrors>1) {
465 AliInfo(Form("\n\n=======================\nFinished at iteration %d : convergence is %e and you required it to be < %e\n=======================\n\n",iIterBayes,convergence,fMaxConvergence));
466 }
467 else if(fNCalcCorrErrors>0) {
468 AliInfo(Form("=======================\nUnfolding of randomized distribution finished at iteration %d with convergence %e \n",iIterBayes,convergence));
469 }
c0b10ad4 470 }
c0b10ad4 471}
472
473//______________________________________________________________
474
475void AliCFUnfolding::CreateUnfolded() {
476 //
477 // Creates the unfolded (T) spectrum from the measured spectrum (M) and the inverse response matrix (INV)
478 // We have P(T) = SUM { P(T|M) * P(M) }
479 // --> T(i) = SUM_k { INV(i,k) * M(k) }
480 //
481
482
483 // clear the unfolded spectrum
769f5114 484 if(fNCalcCorrErrors>0) {
485 //unfold randomized dist
486 fRandomUnfolded->Reset();
487 }
488 else {
489 //unfold measured dist
490 fUnfolded->Reset();
491 }
c0b10ad4 492
7bc6a013 493 for (Long_t iBin=0; iBin<fInverseResponse->GetNbins(); iBin++) {
c0b10ad4 494 Double_t invResponseValue = fInverseResponse->GetBinContent(iBin,fCoordinates2N);
495 GetCoordinates();
85b6bda9 496 Double_t effValue = fEfficiency->GetBinContent(fCoordinatesN_T);
85b6bda9 497 Double_t measuredValue = fMeasured ->GetBinContent(fCoordinatesN_M);
85b6bda9 498 Double_t fill = (effValue>0. ? invResponseValue * measuredValue / effValue : 0.) ;
769f5114 499
85b6bda9 500 if (fill>0.) {
769f5114 501
502 Double_t err = 0.;
503
504 if(fNCalcCorrErrors>0) {
505 fRandomUnfolded->SetBinError(fCoordinatesN_T,err);
506 fRandomUnfolded->AddBinContent(fCoordinatesN_T,fill);
507 }
508 else {
509 fUnfolded->SetBinError(fCoordinatesN_T,err);
510 fUnfolded->AddBinContent(fCoordinatesN_T,fill);
511 }
85b6bda9 512 }
c0b10ad4 513 }
514}
769f5114 515
516//______________________________________________________________
517
518void AliCFUnfolding::CalculateCorrelatedErrors() {
519
520 fRandomizedDist = (THnSparse*) fMeasuredOrig->Clone();
521 fPrior = (THnSparse*) fOriginalPrior->Clone();
522
523 // Step 1: Create randomized distribution (fRandomizedDist) of each bin of the measured spectrum to calculate correlated errors. Poisson statistics: mean = measured value of bin
524 // Step 2: Unfold randomized distribution
525 // Step 3: Store difference of unfolded spectrum from measured distribution and unfolded distribution from randomized distribution -> fDeltaUnfoldedP (TProfile with option "S")
526 // Step 4: Repeat Step 1-3 several times (fNRandomIterations)
527 // Step 5: The spread of fDeltaUnfoldedP for each bin is the error on the unfolded spectrum of that specific bin
528
529 //Do fNRandomIterations = bayes iterations performed
530 for(int i=0; i<fNRandomIterations; i++) {
531 if (fPrior) delete fPrior ;
532 if (fRandomizedDist) delete fRandomizedDist ;
533 fPrior = (THnSparse*) fOriginalPrior->Clone();
534 fRandomizedDist = (THnSparse*) fMeasuredOrig->Clone();
535 CreateRandomizedDist();
536 if (fMeasured) delete fMeasured ;
537 fMeasured = (THnSparse*) fRandomizedDist->Clone();
538 //Unfold fRandomizedDist
539 Unfold();
540 FillDeltaUnfoldedProfile();
541 }
542
543 // Get statistical errors for final unfolded spectrum
544 // ie. spread of each pt bin in fDeltaUnfoldedP
545 Double_t sigma = 0.;
546 Double_t dummy = 0.;
547 for (Long_t iBin=0; iBin<fRandomUnfolded->GetNbins(); iBin++) {
548 dummy = fUnfolded->GetBinContent(iBin,fCoordinatesN_M);
549 sigma = fDeltaUnfoldedP->GetBinError(fCoordinatesN_M[0]);
550 fUnfolded->SetBinError(fCoordinatesN_M,sigma);
551 fNCalcCorrErrors = 2;
552 }
553
554}
555
556//______________________________________________________________
557void AliCFUnfolding::InitDeltaUnfoldedProfile() {
558 //
559 //Initialize the fDeltaUnfoldedP profile
560 //Errors will be filled with spread between unfolded measured and unfolded randomized spectra
561 //
562
563 Int_t nbinsx = fResponse->GetAxis(0)->GetNbins();
564 Double_t xbins[nbinsx];
565 for(int ix=0; ix<nbinsx; ix++) {
566 xbins[ix] = fResponse->GetAxis(0)->GetBinLowEdge(ix+1);
567 }
568 xbins[nbinsx] = fResponse->GetAxis(0)->GetBinUpEdge(nbinsx);
569 fDeltaUnfoldedP = new TProfile("fDeltaUnfoldedP","Profile of pTUnfolded with spread in error",nbinsx,xbins,"S");
570}
571//______________________________________________________________
572void AliCFUnfolding::CreateRandomizedDist() {
573 //
574 // Create randomized dist from measured distribution
575 //
576
577 Double_t random = 0.;
578 Double_t measuredValue = 0.;
579 Double_t measuredError = 0.;
580 for (Long_t iBin=0; iBin<fRandomizedDist->GetNbins(); iBin++) {
581 measuredValue = fMeasuredOrig->GetBinContent(iBin,fCoordinatesN_M); //used as mean
582 measuredError = fMeasuredOrig->GetBinError(fCoordinatesN_M); //used as sigma
583 // random = fRandom3->PoissonD(measuredValue); //doesn't work for normalized spectra, use Gaus (assuming raw counts in bin is large >10)
584 random = fRandom3->Gaus(measuredValue,measuredError);
585 fRandomizedDist->SetBinContent(iBin,random);
586 }
587}
588
589//______________________________________________________________
590void AliCFUnfolding::FillDeltaUnfoldedProfile() {
591 //
592 // Store difference of unfolded spectrum from measured distribution and unfolded distribution from randomized distribution
593 //
594
595 for (Long_t iBin2=0; iBin2<fRandomUnfolded->GetNbins(); iBin2++) {
596 Double_t delta = fUnfolded->GetBinContent(iBin2,fCoordinatesN_M)-fRandomUnfolded->GetBinContent(iBin2,fCoordinatesN_M);
597 fDeltaUnfoldedP->Fill(fDeltaUnfoldedP->GetBinCenter(fCoordinatesN_M[0]),delta);
598 }
599}
600
85b6bda9 601//______________________________________________________________
602
c0b10ad4 603void AliCFUnfolding::GetCoordinates() {
604 //
605 // assign coordinates in Measured and True spaces (dim=N) from coordinates in global space (dim=2N)
606 //
607 for (Int_t i = 0; i<fNVariables ; i++) {
608 fCoordinatesN_M[i] = fCoordinates2N[i];
609 fCoordinatesN_T[i] = fCoordinates2N[i+fNVariables];
610 }
611}
612
613//______________________________________________________________
614
615void AliCFUnfolding::CreateConditional() {
616 //
617 // creates the conditional probability matrix (R*) holding the P(M|T), given the reponse matrix R
618 //
619 // --> R*(i,j) = R(i,j) / SUM_k{ R(k,j) }
620 //
621
622 fConditional = (THnSparse*) fResponse->Clone(); // output of this function
623 fProjResponseInT = (THnSparse*) fPrior->Clone(); // output denominator :
624 // projection of the response matrix on the TRUE axis
85b6bda9 625 Int_t* dim = new Int_t [fNVariables];
626 for (Int_t iDim=0; iDim<fNVariables; iDim++) dim[iDim] = fNVariables+iDim ; //dimensions corresponding to TRUE values (i.e. from N to 2N-1)
627 fProjResponseInT = fConditional->Projection(fNVariables,dim,"E"); //project
628 delete [] dim;
769f5114 629
c0b10ad4 630 // fill the conditional probability matrix
7bc6a013 631 for (Long_t iBin=0; iBin<fResponse->GetNbins(); iBin++) {
c0b10ad4 632 Double_t responseValue = fResponse->GetBinContent(iBin,fCoordinates2N);
633 GetCoordinates();
85b6bda9 634 Double_t projValue = fProjResponseInT->GetBinContent(fCoordinatesN_T);
769f5114 635
85b6bda9 636 Double_t fill = responseValue / projValue ;
637 if (fill>0. || fConditional->GetBinContent(fCoordinates2N)>0.) {
638 fConditional->SetBinContent(fCoordinates2N,fill);
769f5114 639 Double_t err = 0.;
85b6bda9 640 fConditional->SetBinError (fCoordinates2N,err);
641 }
c0b10ad4 642 }
643}
769f5114 644//______________________________________________________________
645
646Int_t AliCFUnfolding::GetDOF() {
647 //
648 // number of dof = number of bins
649 //
650
651 Int_t nDOF = 1 ;
652 for (Int_t iDim=0; iDim<fNVariables; iDim++) {
653 nDOF *= fPrior->GetAxis(iDim)->GetNbins();
654 }
655 AliDebug(0,Form("Number of degrees of freedom = %d",nDOF));
656 return nDOF;
657}
c0b10ad4 658
659//______________________________________________________________
660
661Double_t AliCFUnfolding::GetChi2() {
662 //
663 // Returns the chi2 between unfolded and a priori spectrum
769f5114 664 // This function became absolute with the correlated error calculation.
665 // Errors on the unfolded distribution are not known until the end
666 // Use the convergence criterion instead
c0b10ad4 667 //
668
769f5114 669 Double_t chi2 = 0. ;
670 Double_t error_unf = 0.;
7bc6a013 671 for (Long_t iBin=0; iBin<fPrior->GetNbins(); iBin++) {
85f9f9e1 672 Double_t priorValue = fPrior->GetBinContent(iBin,fCoordinatesN_T);
769f5114 673 error_unf = fUnfolded->GetBinError(fCoordinatesN_T);
85f9f9e1 674 chi2 += (error_unf > 0. ? TMath::Power((fUnfolded->GetBinContent(fCoordinatesN_T) - priorValue)/error_unf,2) / priorValue : 0.) ;
c0b10ad4 675 }
676 return chi2;
677}
678
679//______________________________________________________________
680
769f5114 681Double_t AliCFUnfolding::GetConvergence() {
682 //
683 // Returns convergence criterion = \sum_t ((U_t^{n-1}-U_t^n)/U_t^{n-1})^2
684 // U is unfolded spectrum, t is the bin, n = current, n-1 = previous
685 //
686 Double_t convergence = 0.;
687 Double_t priorValue = 0.;
688 Double_t currentValue = 0.;
689 for (Long_t iBin=0; iBin < fPrior->GetNbins(); iBin++) {
690 priorValue = fPrior->GetBinContent(iBin,fCoordinatesN_T);
691 if (fNCalcCorrErrors > 0)
692 currentValue = fRandomUnfolded->GetBinContent(fCoordinatesN_T);
693 else
694 currentValue = fUnfolded->GetBinContent(fCoordinatesN_T);
695
696 if (priorValue > 0.)
697 convergence += ((priorValue-currentValue)/priorValue)*((priorValue-currentValue)/priorValue);
698 else {
699 AliWarning(Form("priorValue = %f. Adding 0 to convergence criterion.",priorValue));
700 convergence += 0.;
701 }
702 }
703 return convergence;
704}
705
706//______________________________________________________________
707
708void AliCFUnfolding::SetMaxConvergencePerDOF(Double_t val) {
c0b10ad4 709 //
769f5114 710 // Max. convergence criterion per degree of freedom : user setting
711 // convergence criterion = DOF*val; DOF = number of bins
712 // In Jan-Fiete's multiplicity note: Convergence criterion = DOF*0.001^2
c0b10ad4 713 //
714
769f5114 715 fUseCorrelatedErrors = kTRUE;
716 Int_t nDOF = GetDOF() ;
717 fMaxConvergence = val * nDOF ;
718 AliInfo(Form("MaxConvergence = %e. Number of degrees of freedom = %d",fMaxConvergence,nDOF));
c0b10ad4 719}
720
721//______________________________________________________________
722
85b6bda9 723Short_t AliCFUnfolding::Smooth() {
c0b10ad4 724 //
725 // Smoothes the unfolded spectrum
85b6bda9 726 //
727 // By default each cell content is replaced by the average with the neighbouring bins (but not diagonally-neighbouring bins)
728 // However, if a specific function fcn has been defined in UseSmoothing(fcn), the unfolded will be fit and updated using fcn
c0b10ad4 729 //
730
85b6bda9 731 if (fSmoothFunction) {
7bc6a013 732 AliDebug(2,Form("Smoothing spectrum with fit function %p",fSmoothFunction));
85b6bda9 733 return SmoothUsingFunction();
734 }
7036630f 735 else return SmoothUsingNeighbours(fUnfolded);
85b6bda9 736}
737
738//______________________________________________________________
739
7036630f 740Short_t AliCFUnfolding::SmoothUsingNeighbours(THnSparse* hist) {
85b6bda9 741 //
742 // Smoothes the unfolded spectrum using neighouring bins
743 //
744
7036630f 745 Int_t const nDimensions = hist->GetNdimensions() ;
746 Int_t* coordinates = new Int_t[nDimensions];
747
748 Int_t* numBins = new Int_t[nDimensions];
749 for (Int_t iVar=0; iVar<nDimensions; iVar++) numBins[iVar] = hist->GetAxis(iVar)->GetNbins();
c0b10ad4 750
7036630f 751 //need a copy because hist will be updated during the loop, and this creates problems
752 THnSparse* copy = (THnSparse*)hist->Clone();
c0b10ad4 753
7bc6a013 754 for (Long_t iBin=0; iBin<copy->GetNbins(); iBin++) { //loop on non-empty bins
7036630f 755 Double_t content = copy->GetBinContent(iBin,coordinates);
7bc6a013 756 Double_t error2 = TMath::Power(copy->GetBinError(iBin),2);
c0b10ad4 757
758 // skip the under/overflow bins...
759 Bool_t isOutside = kFALSE ;
7036630f 760 for (Int_t iVar=0; iVar<nDimensions; iVar++) {
761 if (coordinates[iVar]<1 || coordinates[iVar]>numBins[iVar]) {
c0b10ad4 762 isOutside=kTRUE;
763 break;
764 }
765 }
766 if (isOutside) continue;
767
768 Int_t neighbours = 0; // number of neighbours to average with
769
7036630f 770 for (Int_t iVar=0; iVar<nDimensions; iVar++) {
771 if (coordinates[iVar] > 1) { // must not be on low edge border
772 coordinates[iVar]-- ; //get lower neighbouring bin
773 content += copy->GetBinContent(coordinates);
774 error2 += TMath::Power(copy->GetBinError(coordinates),2);
c0b10ad4 775 neighbours++;
7036630f 776 coordinates[iVar]++ ; //back to initial coordinate
c0b10ad4 777 }
7036630f 778 if (coordinates[iVar] < numBins[iVar]) { // must not be on up edge border
779 coordinates[iVar]++ ; //get upper neighbouring bin
780 content += copy->GetBinContent(coordinates);
781 error2 += TMath::Power(copy->GetBinError(coordinates),2);
c0b10ad4 782 neighbours++;
7036630f 783 coordinates[iVar]-- ; //back to initial coordinate
c0b10ad4 784 }
785 }
85b6bda9 786 // make an average
7036630f 787 hist->SetBinContent(coordinates,content/(1.+neighbours));
788 hist->SetBinError (coordinates,TMath::Sqrt(error2)/(1.+neighbours));
c0b10ad4 789 }
790 delete [] numBins;
7036630f 791 delete [] coordinates ;
c0b10ad4 792 delete copy;
85b6bda9 793 return 0;
c0b10ad4 794}
795
85b6bda9 796//______________________________________________________________
797
798Short_t AliCFUnfolding::SmoothUsingFunction() {
799 //
800 // Fits the unfolded spectrum using the function fSmoothFunction
801 //
802
7bc6a013 803 AliDebug(0,Form("Smooth function is a %s with option \"%s\" and has %d parameters : ",fSmoothFunction->ClassName(),fSmoothOption,fSmoothFunction->GetNpar()));
85b6bda9 804
7bc6a013 805 for (Int_t iPar=0; iPar<fSmoothFunction->GetNpar(); iPar++) AliDebug(0,Form("par[%d]=%e",iPar,fSmoothFunction->GetParameter(iPar)));
806
807 Int_t fitResult = 0;
85b6bda9 808
809 switch (fNVariables) {
7bc6a013 810 case 1 : fitResult = fUnfolded->Projection(0) ->Fit(fSmoothFunction,fSmoothOption); break;
811 case 2 : fitResult = fUnfolded->Projection(1,0) ->Fit(fSmoothFunction,fSmoothOption); break; // (1,0) instead of (0,1) -> TAxis issue
812 case 3 : fitResult = fUnfolded->Projection(0,1,2)->Fit(fSmoothFunction,fSmoothOption); break;
813 default: AliFatal(Form("Cannot handle such fit in %d dimensions",fNVariables)) ; return 1;
814 }
815
816 if (fitResult != 0) {
817 AliWarning(Form("Fit failed with status %d, stopping the loop",fitResult));
818 return 1;
85b6bda9 819 }
820
821 Int_t nDim = fNVariables;
822 Int_t* bins = new Int_t[nDim]; // number of bins for each variable
823 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
824
825 for (Int_t iVar=0; iVar<nDim; iVar++) {
826 bins[iVar] = fUnfolded->GetAxis(iVar)->GetNbins();
827 nBins *= bins[iVar];
828 }
829
830 Int_t *bin = new Int_t[nDim]; // bin to fill the THnSparse (holding the bin coordinates)
831 Double_t x[3] = {0,0,0} ; // value in bin center (max dimension is 3 (TF3))
832
833 // loop on the bins and update of fUnfolded
834 // THnSparse::Multiply(TF1*) doesn't exist, so let's do it bin by bin
835 for (Long_t iBin=0; iBin<nBins; iBin++) {
836 Long_t bin_tmp = iBin ;
837 for (Int_t iVar=0; iVar<nDim; iVar++) {
838 bin[iVar] = 1 + bin_tmp % bins[iVar] ;
839 bin_tmp /= bins[iVar] ;
840 x[iVar] = fUnfolded->GetAxis(iVar)->GetBinCenter(bin[iVar]);
841 }
842 Double_t functionValue = fSmoothFunction->Eval(x[0],x[1],x[2]) ;
12e419d5 843 fUnfolded->SetBinError (bin,fUnfolded->GetBinError(bin)*functionValue/fUnfolded->GetBinContent(bin));
85b6bda9 844 fUnfolded->SetBinContent(bin,functionValue);
85b6bda9 845 }
846 return 0;
847}
c0b10ad4 848
849//______________________________________________________________
850
851void AliCFUnfolding::CreateFlatPrior() {
852 //
853 // Creates a flat prior distribution
854 //
855
856 AliInfo("Creating a flat a priori distribution");
857
858 // create the frame of the THnSparse given (for example) the one from the efficiency map
859 fPrior = (THnSparse*) fEfficiency->Clone();
860
85b6bda9 861 if (fNVariables != fPrior->GetNdimensions())
862 AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions()));
863
c0b10ad4 864 Int_t nDim = fNVariables;
865 Int_t* bins = new Int_t[nDim]; // number of bins for each variable
866 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
867
868 for (Int_t iVar=0; iVar<nDim; iVar++) {
869 bins[iVar] = fPrior->GetAxis(iVar)->GetNbins();
870 nBins *= bins[iVar];
871 }
872
873 Int_t *bin = new Int_t[nDim]; // bin to fill the THnSparse (holding the bin coordinates)
874
875 // loop that sets 1 in each bin
876 for (Long_t iBin=0; iBin<nBins; iBin++) {
877 Long_t bin_tmp = iBin ;
878 for (Int_t iVar=0; iVar<nDim; iVar++) {
879 bin[iVar] = 1 + bin_tmp % bins[iVar] ;
880 bin_tmp /= bins[iVar] ;
881 }
882 fPrior->SetBinContent(bin,1.); // put 1 everywhere
85b6bda9 883 fPrior->SetBinError (bin,0.); // put 0 everywhere
c0b10ad4 884 }
885
886 fOriginalPrior = (THnSparse*)fPrior->Clone();
887
888 delete [] bin;
889 delete [] bins;
890}