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