]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - CORRFW/AliCFUnfolding.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
fix finding of pad neighbours; remove methods to write them in OCDB
[u/mrichter/AliRoot.git] / CORRFW / AliCFUnfolding.cxx
... / ...
CommitLineData
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
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. //
31// //
32// Then at each iteration, the unfolded spectrum is calculated using //
33// the inverse response : the goal is to get an unfolded spectrum //
34// similar (according to some criterion) to the a priori one. //
35// If the difference is too big, another iteration is performed : //
36// the a priori spectrum is updated to the unfolded one from the //
37// previous iteration, and so on so forth, until the maximum number //
38// of iterations or the similarity criterion is reached. //
39// //
40// 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 //
44// //
45// Currently the user has to define the max. number of iterations //
46// (::SetMaxNumberOfIterations) //
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 //
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// //
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// //
82//---------------------------------------------------------------------//
83// Author : renaud.vernet@cern.ch //
84//---------------------------------------------------------------------//
85
86
87#include "AliCFUnfolding.h"
88#include "TMath.h"
89#include "TAxis.h"
90#include "TF1.h"
91#include "TH1D.h"
92#include "TH2D.h"
93#include "TH3D.h"
94#include "TRandom3.h"
95
96
97ClassImp(AliCFUnfolding)
98
99//______________________________________________________________
100
101AliCFUnfolding::AliCFUnfolding() :
102 TNamed(),
103 fResponseOrig(0x0),
104 fPriorOrig(0x0),
105 fEfficiencyOrig(0x0),
106 fMeasuredOrig(0x0),
107 fMaxNumIterations(0),
108 fNVariables(0),
109 fUseSmoothing(kFALSE),
110 fSmoothFunction(0x0),
111 fSmoothOption("iremn"),
112 fMaxConvergence(0),
113 fNRandomIterations(0),
114 fResponse(0x0),
115 fPrior(0x0),
116 fEfficiency(0x0),
117 fMeasured(0x0),
118 fInverseResponse(0x0),
119 fMeasuredEstimate(0x0),
120 fConditional(0x0),
121 fUnfolded(0x0),
122 fUnfoldedFinal(0x0),
123 fCoordinates2N(0x0),
124 fCoordinatesN_M(0x0),
125 fCoordinatesN_T(0x0),
126 fRandomResponse(0x0),
127 fRandomEfficiency(0x0),
128 fRandomMeasured(0x0),
129 fRandom3(0x0),
130 fDeltaUnfoldedP(0x0),
131 fDeltaUnfoldedN(0x0),
132 fNCalcCorrErrors(0),
133 fRandomSeed(0)
134{
135 //
136 // default constructor
137 //
138}
139
140//______________________________________________________________
141
142AliCFUnfolding::AliCFUnfolding(const Char_t* name, const Char_t* title, const Int_t nVar,
143 const THnSparse* response, const THnSparse* efficiency, const THnSparse* measured, const THnSparse* prior ,
144 Double_t maxConvergencePerDOF, UInt_t randomSeed, Int_t maxNumIterations
145 ) :
146 TNamed(name,title),
147 fResponseOrig((THnSparse*)response->Clone()),
148 fPriorOrig(0x0),
149 fEfficiencyOrig((THnSparse*)efficiency->Clone()),
150 fMeasuredOrig((THnSparse*)measured->Clone()),
151 fMaxNumIterations(maxNumIterations),
152 fNVariables(nVar),
153 fUseSmoothing(kFALSE),
154 fSmoothFunction(0x0),
155 fSmoothOption("iremn"),
156 fMaxConvergence(0),
157 fNRandomIterations(maxNumIterations),
158 fResponse((THnSparse*)response->Clone()),
159 fPrior(0x0),
160 fEfficiency((THnSparse*)efficiency->Clone()),
161 fMeasured((THnSparse*)measured->Clone()),
162 fInverseResponse(0x0),
163 fMeasuredEstimate(0x0),
164 fConditional(0x0),
165 fUnfolded(0x0),
166 fUnfoldedFinal(0x0),
167 fCoordinates2N(0x0),
168 fCoordinatesN_M(0x0),
169 fCoordinatesN_T(0x0),
170 fRandomResponse((THnSparse*)response->Clone()),
171 fRandomEfficiency((THnSparse*)efficiency->Clone()),
172 fRandomMeasured((THnSparse*)measured->Clone()),
173 fRandom3(0x0),
174 fDeltaUnfoldedP(0x0),
175 fDeltaUnfoldedN(0x0),
176 fNCalcCorrErrors(0),
177 fRandomSeed(randomSeed)
178{
179 //
180 // named constructor
181 //
182
183 AliInfo(Form("\n\n--------------------------\nCreating an unfolder :\n--------------------------\nresponse matrix has %d dimension(s)",fResponse->GetNdimensions()));
184
185 if (!prior) CreateFlatPrior(); // if no prior distribution declared, simply use a flat distribution
186 else {
187 fPrior = (THnSparse*) prior->Clone();
188 fPriorOrig = (THnSparse*)fPrior->Clone();
189 if (fPrior->GetNdimensions() != fNVariables)
190 AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions()));
191 }
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()));
199
200
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 }
206
207 fRandomResponse ->SetTitle("Randomized response matrix");
208 fRandomEfficiency->SetTitle("Randomized efficiency");
209 fRandomMeasured ->SetTitle("Randomized measured");
210 SetMaxConvergencePerDOF(maxConvergencePerDOF) ;
211 Init();
212}
213
214//______________________________________________________________
215
216AliCFUnfolding::~AliCFUnfolding() {
217 //
218 // destructor
219 //
220 if (fResponse) delete fResponse;
221 if (fPrior) delete fPrior;
222 if (fEfficiency) delete fEfficiency;
223 if (fEfficiencyOrig) delete fEfficiencyOrig;
224 if (fMeasured) delete fMeasured;
225 if (fMeasuredOrig) delete fMeasuredOrig;
226 if (fSmoothFunction) delete fSmoothFunction;
227 if (fPriorOrig) delete fPriorOrig;
228 if (fInverseResponse) delete fInverseResponse;
229 if (fMeasuredEstimate) delete fMeasuredEstimate;
230 if (fConditional) delete fConditional;
231 if (fUnfolded) delete fUnfolded;
232 if (fUnfoldedFinal) delete fUnfoldedFinal;
233 if (fCoordinates2N) delete [] fCoordinates2N;
234 if (fCoordinatesN_M) delete [] fCoordinatesN_M;
235 if (fCoordinatesN_T) delete [] fCoordinatesN_T;
236 if (fRandomResponse) delete fRandomResponse;
237 if (fRandomEfficiency) delete fRandomEfficiency;
238 if (fRandomMeasured) delete fRandomMeasured;
239 if (fRandom3) delete fRandom3;
240 if (fDeltaUnfoldedP) delete fDeltaUnfoldedP;
241 if (fDeltaUnfoldedN) delete fDeltaUnfoldedN;
242
243}
244
245//______________________________________________________________
246
247void AliCFUnfolding::Init() {
248 //
249 // initialisation function : creates internal settings
250 //
251
252 fRandom3 = new TRandom3(fRandomSeed);
253
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)
259 CreateConditional(); //done only once at initialization
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();
265 fUnfolded->SetTitle("Unfolded");
266 // create the frame of the measurement estimate spectrum
267 fMeasuredEstimate = (THnSparse*) fMeasured->Clone();
268
269 // create the frame of the delta profiles
270 fDeltaUnfoldedP = (THnSparse*)fPrior->Clone();
271 fDeltaUnfoldedP->SetTitle("#Delta unfolded");
272 fDeltaUnfoldedP->Reset();
273 fDeltaUnfoldedN = (THnSparse*)fPrior->Clone();
274 fDeltaUnfoldedN->SetTitle("");
275 fDeltaUnfoldedN->Reset();
276
277
278}
279
280
281//______________________________________________________________
282
283void AliCFUnfolding::CreateEstMeasured() {
284 //
285 // This function creates a estimate (M) of the reconstructed spectrum
286 // given the a priori distribution (T), the efficiency (E) and the conditional matrix (COND)
287 //
288 // --> P(M) = SUM { P(M|T) * P(T) }
289 // --> M(i) = SUM_k { COND(i,k) * T(k) * E (k)}
290 //
291 // This is needed to calculate the inverse response matrix
292 //
293
294
295 // clean the measured estimate spectrum
296 fMeasuredEstimate->Reset();
297
298 THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone();
299 priorTimesEff->Multiply(fEfficiency);
300
301 // fill it
302 for (Long_t iBin=0; iBin<fConditional->GetNbins(); iBin++) {
303 Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N);
304 GetCoordinates();
305 Double_t priorTimesEffValue = priorTimesEff->GetBinContent(fCoordinatesN_T);
306 Double_t fill = conditionalValue * priorTimesEffValue ;
307
308 if (fill>0.) {
309 fMeasuredEstimate->AddBinContent(fCoordinatesN_M,fill);
310 fMeasuredEstimate->SetBinError(fCoordinatesN_M,0.);
311 }
312 }
313 delete priorTimesEff ;
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
327 THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone();
328 priorTimesEff->Multiply(fEfficiency);
329
330 for (Long_t iBin=0; iBin<fConditional->GetNbins(); iBin++) {
331 Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N);
332 GetCoordinates();
333 Double_t estMeasuredValue = fMeasuredEstimate->GetBinContent(fCoordinatesN_M);
334 Double_t priorTimesEffValue = priorTimesEff ->GetBinContent(fCoordinatesN_T);
335 Double_t fill = (estMeasuredValue>0. ? conditionalValue * priorTimesEffValue / estMeasuredValue : 0. ) ;
336 if (fill>0. || fInverseResponse->GetBinContent(fCoordinates2N)>0.) {
337 fInverseResponse->SetBinContent(fCoordinates2N,fill);
338 fInverseResponse->SetBinError (fCoordinates2N,0.);
339 }
340 }
341 delete priorTimesEff ;
342}
343
344//______________________________________________________________
345
346void AliCFUnfolding::Unfold() {
347 //
348 // Main routine called by the user :
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
351 //
352
353 Int_t iIterBayes = 0 ;
354 Double_t convergence = 0.;
355
356 for (iIterBayes=0; iIterBayes<fMaxNumIterations; iIterBayes++) { // bayes iterations
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));
368 break;
369 }
370
371 if (fUseSmoothing) {
372 if (Smooth()) {
373 AliError("Couldn't smooth the unfolded spectrum!!");
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));
379 }
380 return;
381 }
382 }
383
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;
400 CalculateCorrelatedErrors();
401 }
402
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));
408 }
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
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();
426
427 for (Long_t iBin=0; iBin<fInverseResponse->GetNbins(); iBin++) {
428 Double_t invResponseValue = fInverseResponse->GetBinContent(iBin,fCoordinates2N);
429 GetCoordinates();
430 Double_t effValue = fEfficiency->GetBinContent(fCoordinatesN_T);
431 Double_t measuredValue = fMeasured ->GetBinContent(fCoordinatesN_M);
432 Double_t fill = (effValue>0. ? invResponseValue * measuredValue / effValue : 0.) ;
433
434 if (fill>0.) {
435 // set errors to zero
436 // true errors will be filled afterwards
437 Double_t err = 0.;
438 fUnfolded->SetBinError (fCoordinatesN_T,err);
439 fUnfolded->AddBinContent(fCoordinatesN_T,fill);
440 }
441 }
442}
443
444//______________________________________________________________
445
446void AliCFUnfolding::CalculateCorrelatedErrors() {
447
448 // Step 1: Create randomized distribution (fRandomXXXX) of each bin of
449 // the measured spectrum to calculate correlated errors.
450 // Poisson statistics: mean = measured value of bin
451 // Step 2: Unfold randomized distribution
452 // Step 3: Store difference of unfolded spectrum from measured distribution and
453 // unfolded distribution from randomized distribution
454 // -> fDeltaUnfoldedP (TProfile with option "S")
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
458
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 ;
464 fPrior = (THnSparse*) fPriorOrig->Clone();
465
466 // create randomized distribution and stick measured spectrum to it
467 CreateRandomizedDist();
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();
479 fMeasured->SetTitle("Measured");
480
481 //unfold with randomized distributions
482 Unfold();
483 FillDeltaUnfoldedProfile();
484 }
485
486 // Get statistical errors for final unfolded spectrum
487 // ie. spread of each pt bin in fDeltaUnfoldedP
488 Double_t meanx2 = 0.;
489 Double_t mean = 0.;
490 Double_t checksigma = 0.;
491 Double_t entriesInBin = 0.;
492 for (Long_t iBin=0; iBin<fUnfoldedFinal->GetNbins(); iBin++) {
493 fUnfoldedFinal->GetBinContent(iBin,fCoordinatesN_M);
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);
499 //AliDebug(2,Form("filling error %e\n",sigma));
500 fUnfoldedFinal->SetBinError(fCoordinatesN_M,checksigma);
501 }
502
503 // now errors are calculated
504 fNCalcCorrErrors = 2;
505}
506
507//______________________________________________________________
508void AliCFUnfolding::CreateRandomizedDist() {
509 //
510 // Create randomized dist from original measured distribution
511 // This distribution is created several times, each time with a different random number
512 //
513
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);
532 // random = fRandom3->PoissonD(measuredValue); //doesn't work for normalized spectra, use Gaus (assuming raw counts in bin is large >10)
533 fRandomMeasured->SetBinContent(iBin,ran);
534 }
535}
536
537//______________________________________________________________
538void AliCFUnfolding::FillDeltaUnfoldedProfile() {
539 //
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
548 for (Long_t iBin=0; iBin<fUnfolded->GetNbins(); iBin++) {
549 Double_t deltaInBin = fUnfoldedFinal->GetBinContent(iBin,fCoordinatesN_M) - fUnfolded->GetBinContent(iBin);
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
553 //printf("deltaInBin %f\n",deltaInBin);
554 //printf("pt %f\n",ptaxis->GetBinCenter(iBin+1));
555
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
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) ;
567
568 //AliDebug(2,Form("sigma = %e\n",sigma));
569
570 fDeltaUnfoldedP->SetBinError(fCoordinatesN_M,meanx2_nplus1) ;
571 fDeltaUnfoldedP->SetBinContent(fCoordinatesN_M,mean_nplus1) ;
572 fDeltaUnfoldedN->SetBinContent(fCoordinatesN_M,entriesInBin+1);
573 }
574}
575
576//______________________________________________________________
577
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
597 fConditional = (THnSparse*) fResponse->Clone(); // output of this function
598
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)
601
602 THnSparse* responseInT = fConditional->Projection(fNVariables,dim,"E"); // output denominator :
603 // projection of the response matrix on the TRUE axis
604 delete [] dim;
605
606 // fill the conditional probability matrix
607 for (Long_t iBin=0; iBin<fResponse->GetNbins(); iBin++) {
608 Double_t responseValue = fResponse->GetBinContent(iBin,fCoordinates2N);
609 GetCoordinates();
610 Double_t projValue = responseInT->GetBinContent(fCoordinatesN_T);
611
612 Double_t fill = responseValue / projValue ;
613 if (fill>0. || fConditional->GetBinContent(fCoordinates2N)>0.) {
614 fConditional->SetBinContent(fCoordinates2N,fill);
615 Double_t err = 0.;
616 fConditional->SetBinError (fCoordinates2N,err);
617 }
618 }
619 delete responseInT ;
620}
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}
635
636//______________________________________________________________
637
638Double_t AliCFUnfolding::GetChi2() {
639 //
640 // Returns the chi2 between unfolded and a priori spectrum
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
644 //
645
646 Double_t chi2 = 0. ;
647 Double_t error_unf = 0.;
648 for (Long_t iBin=0; iBin<fPrior->GetNbins(); iBin++) {
649 Double_t priorValue = fPrior->GetBinContent(iBin,fCoordinatesN_T);
650 error_unf = fUnfolded->GetBinError(fCoordinatesN_T);
651 chi2 += (error_unf > 0. ? TMath::Power((fUnfolded->GetBinContent(fCoordinatesN_T) - priorValue)/error_unf,2) / priorValue : 0.) ;
652 }
653 return chi2;
654}
655
656//______________________________________________________________
657
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);
668 currentValue = fUnfolded->GetBinContent(fCoordinatesN_T);
669
670 if (priorValue > 0.)
671 convergence += ((priorValue-currentValue)/priorValue)*((priorValue-currentValue)/priorValue);
672 else
673 AliWarning(Form("priorValue = %f. Adding 0 to convergence criterion.",priorValue));
674 }
675 return convergence;
676}
677
678//______________________________________________________________
679
680void AliCFUnfolding::SetMaxConvergencePerDOF(Double_t val) {
681 //
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
685 //
686
687 Int_t nDOF = GetDOF() ;
688 fMaxConvergence = val * nDOF ;
689 AliInfo(Form("MaxConvergence = %e. Number of degrees of freedom = %d",fMaxConvergence,nDOF));
690}
691
692//______________________________________________________________
693
694Short_t AliCFUnfolding::Smooth() {
695 //
696 // Smoothes the unfolded spectrum
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
700 //
701
702 if (fSmoothFunction) {
703 AliDebug(2,Form("Smoothing spectrum with fit function %p",fSmoothFunction));
704 return SmoothUsingFunction();
705 }
706 else return SmoothUsingNeighbours(fUnfolded);
707}
708
709//______________________________________________________________
710
711Short_t AliCFUnfolding::SmoothUsingNeighbours(THnSparse* hist) {
712 //
713 // Smoothes the unfolded spectrum using neighouring bins
714 //
715
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();
721
722 //need a copy because hist will be updated during the loop, and this creates problems
723 THnSparse* copy = (THnSparse*)hist->Clone();
724
725 for (Long_t iBin=0; iBin<copy->GetNbins(); iBin++) { //loop on non-empty bins
726 Double_t content = copy->GetBinContent(iBin,coordinates);
727 Double_t error2 = TMath::Power(copy->GetBinError(iBin),2);
728
729 // skip the under/overflow bins...
730 Bool_t isOutside = kFALSE ;
731 for (Int_t iVar=0; iVar<nDimensions; iVar++) {
732 if (coordinates[iVar]<1 || coordinates[iVar]>numBins[iVar]) {
733 isOutside=kTRUE;
734 break;
735 }
736 }
737 if (isOutside) continue;
738
739 Int_t neighbours = 0; // number of neighbours to average with
740
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);
746 neighbours++;
747 coordinates[iVar]++ ; //back to initial coordinate
748 }
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);
753 neighbours++;
754 coordinates[iVar]-- ; //back to initial coordinate
755 }
756 }
757 // make an average
758 hist->SetBinContent(coordinates,content/(1.+neighbours));
759 hist->SetBinError (coordinates,TMath::Sqrt(error2)/(1.+neighbours));
760 }
761 delete [] numBins;
762 delete [] coordinates ;
763 delete copy;
764 return 0;
765}
766
767//______________________________________________________________
768
769Short_t AliCFUnfolding::SmoothUsingFunction() {
770 //
771 // Fits the unfolded spectrum using the function fSmoothFunction
772 //
773
774 AliDebug(0,Form("Smooth function is a %s with option \"%s\" and has %d parameters : ",fSmoothFunction->ClassName(),fSmoothOption,fSmoothFunction->GetNpar()));
775
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;
779
780 switch (fNVariables) {
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;
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]) ;
814 fUnfolded->SetBinError (bin,fUnfolded->GetBinError(bin)*functionValue/fUnfolded->GetBinContent(bin));
815 fUnfolded->SetBinContent(bin,functionValue);
816 }
817 delete [] bins;
818 delete [] bin ;
819 return 0;
820}
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();
833 fPrior->SetTitle("Prior");
834
835 if (fNVariables != fPrior->GetNdimensions())
836 AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions()));
837
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
857 fPrior->SetBinError (bin,0.); // put 0 everywhere
858 }
859
860 fPriorOrig = (THnSparse*)fPrior->Clone();
861
862 delete [] bin;
863 delete [] bins;
864}