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