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