1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //---------------------------------------------------------------------//
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 //
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. //
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. //
40 // Currently the similarity criterion is the Chi2 between the a priori //
41 // and the unfolded spectrum (OBSOLETE). //
43 // Chi2 calculation became absolute with the correlated error //
45 // Errors on the unfolded distribution are not known until the end //
46 // Use the convergence criterion instead //
48 // Currently the user has to define the max. number of iterations //
49 // (::SetMaxNumberOfIterations) //
51 // - the chi2 below which the procedure will stop //
52 // (::SetMaxChi2 or ::SetMaxChi2PerDOF) (OBSOLETE) //
53 // - the convergence criterion below which the procedure will stop //
54 // SetMaxConvergencePerDOF(Double_t val); //
56 // Correlated error calculation can be activated by using: //
57 // SetUseCorrelatedErrors(Bool_t b) in combination with convergence //
59 // Documentation about correlated error calculation method can be //
60 // found in AliCFUnfolding::CalculateCorrelatedErrors() //
61 // Author: marta.verweij@cern.ch //
63 // An optional possibility is to smooth the unfolded spectrum at the //
64 // end of each iteration, either using a fit function //
65 // (only if #dimensions <=3) //
66 // or a simple averaging using the neighbouring bins values. //
67 // This is possible calling the function ::UseSmoothing //
68 // If no argument is passed to this function, then the second option //
73 // With this approach, the efficiency map must be calculated //
74 // with *simulated* values only, otherwise the method won't work. //
76 // ex: efficiency(bin_pt) = number_rec(bin_pt) / number_sim(bin_pt) //
78 // the pt bin "bin_pt" must always be the same in both the efficiency //
79 // numerator and denominator. //
80 // This is why the efficiency map has to be created by a method //
81 // from which both reconstructed and simulated values are accessible //
85 //---------------------------------------------------------------------//
86 // Author : renaud.vernet@cern.ch //
87 //---------------------------------------------------------------------//
90 #include "AliCFUnfolding.h"
101 ClassImp(AliCFUnfolding)
103 //______________________________________________________________
105 AliCFUnfolding::AliCFUnfolding() :
112 fMaxNumIterations(20),
114 fUseSmoothing(kFALSE),
115 fSmoothFunction(0x0),
118 fUseCorrelatedErrors(kTRUE),
119 fNRandomIterations(20),
121 fInverseResponse(0x0),
122 fMeasuredEstimate(0x0),
124 fProjResponseInT(0x0),
127 fCoordinatesN_M(0x0),
128 fCoordinatesN_T(0x0),
129 fRandomizedDist(0x0),
131 fRandomUnfolded(0x0),
132 fDeltaUnfoldedP(0x0),
137 // default constructor
141 //______________________________________________________________
143 AliCFUnfolding::AliCFUnfolding(const Char_t* name, const Char_t* title, const Int_t nVar,
144 const THnSparse* response, const THnSparse* efficiency, const THnSparse* measured, const THnSparse* prior) :
146 fResponse((THnSparse*)response->Clone()),
148 fEfficiency((THnSparse*)efficiency->Clone()),
149 fMeasured((THnSparse*)measured->Clone()),
150 fMeasuredOrig((THnSparse*)measured->Clone()),
151 fMaxNumIterations(0),
153 fUseSmoothing(kFALSE),
154 fSmoothFunction(0x0),
157 fUseCorrelatedErrors(kTRUE),
158 fNRandomIterations(20),
160 fInverseResponse(0x0),
161 fMeasuredEstimate(0x0),
163 fProjResponseInT(0x0),
166 fCoordinatesN_M(0x0),
167 fCoordinatesN_T(0x0),
168 fRandomizedDist(0x0),
170 fRandomUnfolded(0x0),
171 fDeltaUnfoldedP(0x0),
179 AliInfo(Form("\n\n--------------------------\nCreating an unfolder :\n--------------------------\nresponse matrix has %d dimension(s)",fResponse->GetNdimensions()));
181 if (!prior) CreateFlatPrior(); // if no prior distribution declared, simply use a flat distribution
183 fPrior = (THnSparse*) prior->Clone();
184 fOriginalPrior = (THnSparse*)fPrior->Clone();
185 if (fPrior->GetNdimensions() != fNVariables)
186 AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions()));
189 if (fEfficiency->GetNdimensions() != fNVariables)
190 AliFatal(Form("The efficiency matrix should have %d dimensions, and it has actually %d",fNVariables,fEfficiency->GetNdimensions()));
191 if (fMeasured->GetNdimensions() != fNVariables)
192 AliFatal(Form("The measured matrix should have %d dimensions, and it has actually %d",fNVariables,fMeasured->GetNdimensions()));
193 if (fResponse->GetNdimensions() != 2*fNVariables)
194 AliFatal(Form("The response matrix should have %d dimensions, and it has actually %d",2*fNVariables,fResponse->GetNdimensions()));
197 for (Int_t iVar=0; iVar<fNVariables; iVar++) {
198 AliInfo(Form("prior matrix has %d bins in dimension %d",fPrior ->GetAxis(iVar)->GetNbins(),iVar));
199 AliInfo(Form("efficiency matrix has %d bins in dimension %d",fEfficiency->GetAxis(iVar)->GetNbins(),iVar));
200 AliInfo(Form("measured matrix has %d bins in dimension %d",fMeasured ->GetAxis(iVar)->GetNbins(),iVar));
207 //______________________________________________________________
209 AliCFUnfolding::AliCFUnfolding(const AliCFUnfolding& c) :
211 fResponse((THnSparse*)c.fResponse->Clone()),
212 fPrior((THnSparse*)c.fPrior->Clone()),
213 fEfficiency((THnSparse*)c.fEfficiency->Clone()),
214 fMeasured((THnSparse*)c.fMeasured->Clone()),
215 fMeasuredOrig((THnSparse*)c.fMeasuredOrig->Clone()),
216 fMaxNumIterations(c.fMaxNumIterations),
217 fNVariables(c.fNVariables),
218 fUseSmoothing(c.fUseSmoothing),
219 fSmoothFunction((TF1*)c.fSmoothFunction->Clone()),
220 fSmoothOption(c.fSmoothOption),
221 fMaxConvergence(c.fMaxConvergence),
222 fUseCorrelatedErrors(c.fUseCorrelatedErrors),
223 fNRandomIterations(c.fNRandomIterations),
224 fOriginalPrior((THnSparse*)c.fOriginalPrior->Clone()),
225 fInverseResponse((THnSparse*)c.fInverseResponse->Clone()),
226 fMeasuredEstimate((THnSparse*)fMeasuredEstimate->Clone()),
227 fConditional((THnSparse*)c.fConditional->Clone()),
228 fProjResponseInT((THnSparse*)c.fProjResponseInT->Clone()),
229 fUnfolded((THnSparse*)c.fUnfolded->Clone()),
230 fCoordinates2N(new Int_t(*c.fCoordinates2N)),
231 fCoordinatesN_M(new Int_t(*c.fCoordinatesN_M)),
232 fCoordinatesN_T(new Int_t(*c.fCoordinatesN_T)),
233 fRandomizedDist((THnSparse*)c.fRandomizedDist->Clone()),
234 fRandom3((TRandom3*)c.fRandom3->Clone()),
235 fRandomUnfolded((THnSparse*)c.fRandomUnfolded->Clone()),
236 fDeltaUnfoldedP((TProfile*)c.fDeltaUnfoldedP),
237 fNCalcCorrErrors(c.fNCalcCorrErrors),
238 fRandomSeed(c.fRandomSeed)
245 //______________________________________________________________
247 AliCFUnfolding& AliCFUnfolding::operator=(const AliCFUnfolding& c) {
249 // assignment operator
253 TNamed::operator=(c);
254 fResponse = (THnSparse*)c.fResponse->Clone() ;
255 fPrior = (THnSparse*)c.fPrior->Clone() ;
256 fEfficiency = (THnSparse*)c.fEfficiency->Clone() ;
257 fMeasured = (THnSparse*)c.fMeasured->Clone() ;
258 fMeasuredOrig = ((THnSparse*)c.fMeasuredOrig->Clone()),
259 fMaxNumIterations = c.fMaxNumIterations ;
260 fNVariables = c.fNVariables ;
261 fMaxConvergence = c.fMaxConvergence ;
262 fUseSmoothing = c.fUseSmoothing ;
263 fSmoothFunction = (TF1*)c.fSmoothFunction->Clone();
264 fSmoothOption = c.fSmoothOption ;
265 fUseCorrelatedErrors = c.fUseCorrelatedErrors;
266 fNRandomIterations = c.fNRandomIterations;
267 fOriginalPrior = (THnSparse*)c.fOriginalPrior->Clone() ;
268 fInverseResponse = (THnSparse*)c.fInverseResponse->Clone() ;
269 fMeasuredEstimate = (THnSparse*)fMeasuredEstimate->Clone() ;
270 fConditional = (THnSparse*)c.fConditional->Clone() ;
271 fProjResponseInT = (THnSparse*)c.fProjResponseInT->Clone() ;
272 fUnfolded = (THnSparse*)c.fUnfolded->Clone() ;
273 fCoordinates2N = new Int_t(*c.fCoordinates2N) ;
274 fCoordinatesN_M = new Int_t(*c.fCoordinatesN_M) ;
275 fCoordinatesN_T = new Int_t(*c.fCoordinatesN_T) ;
276 fRandomizedDist = (THnSparse*)c.fRandomizedDist->Clone();
277 fRandom3 = (TRandom3*)c.fRandom3->Clone();
278 fRandomUnfolded = (THnSparse*)c.fRandomUnfolded->Clone();
279 fDeltaUnfoldedP = (TProfile*)c.fDeltaUnfoldedP;
280 fNCalcCorrErrors = c.fNCalcCorrErrors;
281 fRandomSeed = c.fRandomSeed ;
286 //______________________________________________________________
288 AliCFUnfolding::~AliCFUnfolding() {
292 if (fResponse) delete fResponse;
293 if (fPrior) delete fPrior;
294 if (fEfficiency) delete fEfficiency;
295 if (fMeasured) delete fMeasured;
296 if (fMeasuredOrig) delete fMeasuredOrig;
297 if (fSmoothFunction) delete fSmoothFunction;
298 if (fOriginalPrior) delete fOriginalPrior;
299 if (fInverseResponse) delete fInverseResponse;
300 if (fMeasuredEstimate) delete fMeasuredEstimate;
301 if (fConditional) delete fConditional;
302 if (fProjResponseInT) delete fProjResponseInT;
303 if (fCoordinates2N) delete [] fCoordinates2N;
304 if (fCoordinatesN_M) delete [] fCoordinatesN_M;
305 if (fCoordinatesN_T) delete [] fCoordinatesN_T;
306 if (fRandomizedDist) delete fRandomizedDist;
307 if (fRandom3) delete fRandom3;
308 if (fRandomUnfolded) delete fRandomUnfolded;
309 if (fDeltaUnfoldedP) delete fDeltaUnfoldedP;
312 //______________________________________________________________
314 void AliCFUnfolding::Init() {
316 // initialisation function : creates internal settings
319 fRandom3 = new TRandom3(fRandomSeed);
321 fCoordinates2N = new Int_t[2*fNVariables];
322 fCoordinatesN_M = new Int_t[fNVariables];
323 fCoordinatesN_T = new Int_t[fNVariables];
325 // create the matrix of conditional probabilities P(M|T)
328 // create the frame of the inverse response matrix
329 fInverseResponse = (THnSparse*) fResponse->Clone();
330 // create the frame of the unfolded spectrum
331 fUnfolded = (THnSparse*) fPrior->Clone();
332 // create the frame of the random unfolded spectrum
333 fRandomUnfolded = (THnSparse*) fPrior->Clone();
334 // create the frame of the measurement estimate spectrum
335 fMeasuredEstimate = (THnSparse*) fMeasured->Clone();
336 // create the frame of the original measurement spectrum
337 fMeasuredOrig = (THnSparse*) fMeasured->Clone();
339 InitDeltaUnfoldedProfile();
343 //______________________________________________________________
345 void AliCFUnfolding::CreateEstMeasured() {
347 // This function creates a estimate (M) of the reconstructed spectrum
348 // given the a priori distribution (T), the efficiency (E) and the conditional matrix (COND)
350 // --> P(M) = SUM { P(M|T) * P(T) }
351 // --> M(i) = SUM_k { COND(i,k) * T(k) * E (k)}
353 // This is needed to calculate the inverse response matrix
357 // clean the measured estimate spectrum
358 fMeasuredEstimate->Reset();
360 THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone();
361 priorTimesEff->Multiply(fEfficiency);
364 for (Long_t iBin=0; iBin<fConditional->GetNbins(); iBin++) {
365 Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N);
367 Double_t priorTimesEffValue = priorTimesEff->GetBinContent(fCoordinatesN_T);
368 Double_t fill = conditionalValue * priorTimesEffValue ;
371 fMeasuredEstimate->AddBinContent(fCoordinatesN_M,fill);
372 fMeasuredEstimate->SetBinError(fCoordinatesN_M,0.);
375 delete priorTimesEff ;
378 //______________________________________________________________
380 void AliCFUnfolding::CreateInvResponse() {
382 // Creates the inverse response matrix (INV) with Bayesian method
383 // : uses the conditional matrix (COND), the prior probabilities (T) and the efficiency map (E)
385 // --> P(T|M) = P(M|T) * P(T) * eff(T) / SUM { P(M|T) * P(T) }
386 // --> INV(i,j) = COND(i,j) * T(j) * E(j) / SUM_k { COND(i,k) * T(k) }
389 THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone();
390 priorTimesEff->Multiply(fEfficiency);
392 for (Long_t iBin=0; iBin<fConditional->GetNbins(); iBin++) {
393 Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N);
395 Double_t estMeasuredValue = fMeasuredEstimate->GetBinContent(fCoordinatesN_M);
396 Double_t priorTimesEffValue = priorTimesEff ->GetBinContent(fCoordinatesN_T);
397 Double_t fill = (estMeasuredValue>0. ? conditionalValue * priorTimesEffValue / estMeasuredValue : 0. ) ;
398 if (fill>0. || fInverseResponse->GetBinContent(fCoordinates2N)>0.) {
399 fInverseResponse->SetBinContent(fCoordinates2N,fill);
400 fInverseResponse->SetBinError (fCoordinates2N,0.);
403 delete priorTimesEff ;
406 //______________________________________________________________
408 void AliCFUnfolding::Unfold() {
410 // Main routine called by the user :
411 // it calculates the unfolded spectrum from the response matrix, measured spectrum and efficiency
412 // several iterations are performed until a reasonable chi2 or convergence criterion is reached
415 Int_t iIterBayes = 0 ;
416 Double_t convergence = 0.;
418 for (iIterBayes=0; iIterBayes<fMaxNumIterations; iIterBayes++) { // bayes iterations
422 if (fUseCorrelatedErrors) {
423 convergence = GetConvergence();
424 AliDebug(0,Form("convergence at iteration %d is %e",iIterBayes,convergence));
426 else AliWarning("No errors will be calculated. Activate SetUseCorrelatedErrors(kTRUE)\n");
428 if (fMaxConvergence>0. && convergence<fMaxConvergence && fNCalcCorrErrors<1) {
429 fNRandomIterations=iIterBayes;
433 // update the prior distribution
436 AliError("Couldn't smooth the unfolded spectrum!!");
437 if (fUseCorrelatedErrors) {
438 if (fNCalcCorrErrors>0) {
439 AliInfo(Form("=======================\nUnfolding of randomized distribution finished at iteration %d with convergence %e \n",iIterBayes,convergence));
442 AliInfo(Form("\n\n=======================\nFinished at iteration %d : convergence is %e and you required it to be < %e\n=======================\n\n",iIterBayes,convergence,fMaxConvergence));
448 if(fNCalcCorrErrors>0) {
449 if (fPrior) delete fPrior ;
450 fPrior = (THnSparse*)fRandomUnfolded->Clone() ;
453 if (fPrior) delete fPrior ;
454 fPrior = (THnSparse*)fUnfolded->Clone() ;
458 if (fUseCorrelatedErrors && fNCalcCorrErrors==0) {
460 CalculateCorrelatedErrors();
463 if (fUseCorrelatedErrors) {
464 if (fNCalcCorrErrors>1) {
465 AliInfo(Form("\n\n=======================\nFinished at iteration %d : convergence is %e and you required it to be < %e\n=======================\n\n",iIterBayes,convergence,fMaxConvergence));
467 else if(fNCalcCorrErrors>0) {
468 AliInfo(Form("=======================\nUnfolding of randomized distribution finished at iteration %d with convergence %e \n",iIterBayes,convergence));
473 //______________________________________________________________
475 void AliCFUnfolding::CreateUnfolded() {
477 // Creates the unfolded (T) spectrum from the measured spectrum (M) and the inverse response matrix (INV)
478 // We have P(T) = SUM { P(T|M) * P(M) }
479 // --> T(i) = SUM_k { INV(i,k) * M(k) }
483 // clear the unfolded spectrum
484 if(fNCalcCorrErrors>0) {
485 //unfold randomized dist
486 fRandomUnfolded->Reset();
489 //unfold measured dist
493 for (Long_t iBin=0; iBin<fInverseResponse->GetNbins(); iBin++) {
494 Double_t invResponseValue = fInverseResponse->GetBinContent(iBin,fCoordinates2N);
496 Double_t effValue = fEfficiency->GetBinContent(fCoordinatesN_T);
497 Double_t measuredValue = fMeasured ->GetBinContent(fCoordinatesN_M);
498 Double_t fill = (effValue>0. ? invResponseValue * measuredValue / effValue : 0.) ;
504 if(fNCalcCorrErrors>0) {
505 fRandomUnfolded->SetBinError(fCoordinatesN_T,err);
506 fRandomUnfolded->AddBinContent(fCoordinatesN_T,fill);
509 fUnfolded->SetBinError(fCoordinatesN_T,err);
510 fUnfolded->AddBinContent(fCoordinatesN_T,fill);
516 //______________________________________________________________
518 void AliCFUnfolding::CalculateCorrelatedErrors() {
520 fRandomizedDist = (THnSparse*) fMeasuredOrig->Clone();
521 fPrior = (THnSparse*) fOriginalPrior->Clone();
523 // Step 1: Create randomized distribution (fRandomizedDist) of each bin of the measured spectrum to calculate correlated errors. Poisson statistics: mean = measured value of bin
524 // Step 2: Unfold randomized distribution
525 // Step 3: Store difference of unfolded spectrum from measured distribution and unfolded distribution from randomized distribution -> fDeltaUnfoldedP (TProfile with option "S")
526 // Step 4: Repeat Step 1-3 several times (fNRandomIterations)
527 // Step 5: The spread of fDeltaUnfoldedP for each bin is the error on the unfolded spectrum of that specific bin
529 //Do fNRandomIterations = bayes iterations performed
530 for(int i=0; i<fNRandomIterations; i++) {
531 if (fPrior) delete fPrior ;
532 if (fRandomizedDist) delete fRandomizedDist ;
533 fPrior = (THnSparse*) fOriginalPrior->Clone();
534 fRandomizedDist = (THnSparse*) fMeasuredOrig->Clone();
535 CreateRandomizedDist();
536 if (fMeasured) delete fMeasured ;
537 fMeasured = (THnSparse*) fRandomizedDist->Clone();
538 //Unfold fRandomizedDist
540 FillDeltaUnfoldedProfile();
543 // Get statistical errors for final unfolded spectrum
544 // ie. spread of each pt bin in fDeltaUnfoldedP
547 for (Long_t iBin=0; iBin<fRandomUnfolded->GetNbins(); iBin++) {
548 dummy = fUnfolded->GetBinContent(iBin,fCoordinatesN_M);
549 sigma = fDeltaUnfoldedP->GetBinError(fCoordinatesN_M[0]);
550 fUnfolded->SetBinError(fCoordinatesN_M,sigma);
551 fNCalcCorrErrors = 2;
556 //______________________________________________________________
557 void AliCFUnfolding::InitDeltaUnfoldedProfile() {
559 //Initialize the fDeltaUnfoldedP profile
560 //Errors will be filled with spread between unfolded measured and unfolded randomized spectra
563 Int_t nbinsx = fResponse->GetAxis(0)->GetNbins();
564 Double_t xbins[nbinsx];
565 for(int ix=0; ix<nbinsx; ix++) {
566 xbins[ix] = fResponse->GetAxis(0)->GetBinLowEdge(ix+1);
568 xbins[nbinsx] = fResponse->GetAxis(0)->GetBinUpEdge(nbinsx);
569 fDeltaUnfoldedP = new TProfile("fDeltaUnfoldedP","Profile of pTUnfolded with spread in error",nbinsx,xbins,"S");
571 //______________________________________________________________
572 void AliCFUnfolding::CreateRandomizedDist() {
574 // Create randomized dist from measured distribution
577 Double_t random = 0.;
578 Double_t measuredValue = 0.;
579 Double_t measuredError = 0.;
580 for (Long_t iBin=0; iBin<fRandomizedDist->GetNbins(); iBin++) {
581 measuredValue = fMeasuredOrig->GetBinContent(iBin,fCoordinatesN_M); //used as mean
582 measuredError = fMeasuredOrig->GetBinError(fCoordinatesN_M); //used as sigma
583 // random = fRandom3->PoissonD(measuredValue); //doesn't work for normalized spectra, use Gaus (assuming raw counts in bin is large >10)
584 random = fRandom3->Gaus(measuredValue,measuredError);
585 fRandomizedDist->SetBinContent(iBin,random);
589 //______________________________________________________________
590 void AliCFUnfolding::FillDeltaUnfoldedProfile() {
592 // Store difference of unfolded spectrum from measured distribution and unfolded distribution from randomized distribution
595 for (Long_t iBin2=0; iBin2<fRandomUnfolded->GetNbins(); iBin2++) {
596 Double_t delta = fUnfolded->GetBinContent(iBin2,fCoordinatesN_M)-fRandomUnfolded->GetBinContent(iBin2,fCoordinatesN_M);
597 fDeltaUnfoldedP->Fill(fDeltaUnfoldedP->GetBinCenter(fCoordinatesN_M[0]),delta);
601 //______________________________________________________________
603 void AliCFUnfolding::GetCoordinates() {
605 // assign coordinates in Measured and True spaces (dim=N) from coordinates in global space (dim=2N)
607 for (Int_t i = 0; i<fNVariables ; i++) {
608 fCoordinatesN_M[i] = fCoordinates2N[i];
609 fCoordinatesN_T[i] = fCoordinates2N[i+fNVariables];
613 //______________________________________________________________
615 void AliCFUnfolding::CreateConditional() {
617 // creates the conditional probability matrix (R*) holding the P(M|T), given the reponse matrix R
619 // --> R*(i,j) = R(i,j) / SUM_k{ R(k,j) }
622 fConditional = (THnSparse*) fResponse->Clone(); // output of this function
623 fProjResponseInT = (THnSparse*) fPrior->Clone(); // output denominator :
624 // projection of the response matrix on the TRUE axis
625 Int_t* dim = new Int_t [fNVariables];
626 for (Int_t iDim=0; iDim<fNVariables; iDim++) dim[iDim] = fNVariables+iDim ; //dimensions corresponding to TRUE values (i.e. from N to 2N-1)
627 fProjResponseInT = fConditional->Projection(fNVariables,dim,"E"); //project
630 // fill the conditional probability matrix
631 for (Long_t iBin=0; iBin<fResponse->GetNbins(); iBin++) {
632 Double_t responseValue = fResponse->GetBinContent(iBin,fCoordinates2N);
634 Double_t projValue = fProjResponseInT->GetBinContent(fCoordinatesN_T);
636 Double_t fill = responseValue / projValue ;
637 if (fill>0. || fConditional->GetBinContent(fCoordinates2N)>0.) {
638 fConditional->SetBinContent(fCoordinates2N,fill);
640 fConditional->SetBinError (fCoordinates2N,err);
644 //______________________________________________________________
646 Int_t AliCFUnfolding::GetDOF() {
648 // number of dof = number of bins
652 for (Int_t iDim=0; iDim<fNVariables; iDim++) {
653 nDOF *= fPrior->GetAxis(iDim)->GetNbins();
655 AliDebug(0,Form("Number of degrees of freedom = %d",nDOF));
659 //______________________________________________________________
661 Double_t AliCFUnfolding::GetChi2() {
663 // Returns the chi2 between unfolded and a priori spectrum
664 // This function became absolute with the correlated error calculation.
665 // Errors on the unfolded distribution are not known until the end
666 // Use the convergence criterion instead
670 Double_t error_unf = 0.;
671 for (Long_t iBin=0; iBin<fPrior->GetNbins(); iBin++) {
672 Double_t priorValue = fPrior->GetBinContent(iBin,fCoordinatesN_T);
673 error_unf = fUnfolded->GetBinError(fCoordinatesN_T);
674 chi2 += (error_unf > 0. ? TMath::Power((fUnfolded->GetBinContent(fCoordinatesN_T) - priorValue)/error_unf,2) / priorValue : 0.) ;
679 //______________________________________________________________
681 Double_t AliCFUnfolding::GetConvergence() {
683 // Returns convergence criterion = \sum_t ((U_t^{n-1}-U_t^n)/U_t^{n-1})^2
684 // U is unfolded spectrum, t is the bin, n = current, n-1 = previous
686 Double_t convergence = 0.;
687 Double_t priorValue = 0.;
688 Double_t currentValue = 0.;
689 for (Long_t iBin=0; iBin < fPrior->GetNbins(); iBin++) {
690 priorValue = fPrior->GetBinContent(iBin,fCoordinatesN_T);
691 if (fNCalcCorrErrors > 0)
692 currentValue = fRandomUnfolded->GetBinContent(fCoordinatesN_T);
694 currentValue = fUnfolded->GetBinContent(fCoordinatesN_T);
697 convergence += ((priorValue-currentValue)/priorValue)*((priorValue-currentValue)/priorValue);
699 AliWarning(Form("priorValue = %f. Adding 0 to convergence criterion.",priorValue));
706 //______________________________________________________________
708 void AliCFUnfolding::SetMaxConvergencePerDOF(Double_t val) {
710 // Max. convergence criterion per degree of freedom : user setting
711 // convergence criterion = DOF*val; DOF = number of bins
712 // In Jan-Fiete's multiplicity note: Convergence criterion = DOF*0.001^2
715 fUseCorrelatedErrors = kTRUE;
716 Int_t nDOF = GetDOF() ;
717 fMaxConvergence = val * nDOF ;
718 AliInfo(Form("MaxConvergence = %e. Number of degrees of freedom = %d",fMaxConvergence,nDOF));
721 //______________________________________________________________
723 Short_t AliCFUnfolding::Smooth() {
725 // Smoothes the unfolded spectrum
727 // By default each cell content is replaced by the average with the neighbouring bins (but not diagonally-neighbouring bins)
728 // However, if a specific function fcn has been defined in UseSmoothing(fcn), the unfolded will be fit and updated using fcn
731 if (fSmoothFunction) {
732 AliDebug(2,Form("Smoothing spectrum with fit function %p",fSmoothFunction));
733 return SmoothUsingFunction();
735 else return SmoothUsingNeighbours(fUnfolded);
738 //______________________________________________________________
740 Short_t AliCFUnfolding::SmoothUsingNeighbours(THnSparse* hist) {
742 // Smoothes the unfolded spectrum using neighouring bins
745 Int_t const nDimensions = hist->GetNdimensions() ;
746 Int_t* coordinates = new Int_t[nDimensions];
748 Int_t* numBins = new Int_t[nDimensions];
749 for (Int_t iVar=0; iVar<nDimensions; iVar++) numBins[iVar] = hist->GetAxis(iVar)->GetNbins();
751 //need a copy because hist will be updated during the loop, and this creates problems
752 THnSparse* copy = (THnSparse*)hist->Clone();
754 for (Long_t iBin=0; iBin<copy->GetNbins(); iBin++) { //loop on non-empty bins
755 Double_t content = copy->GetBinContent(iBin,coordinates);
756 Double_t error2 = TMath::Power(copy->GetBinError(iBin),2);
758 // skip the under/overflow bins...
759 Bool_t isOutside = kFALSE ;
760 for (Int_t iVar=0; iVar<nDimensions; iVar++) {
761 if (coordinates[iVar]<1 || coordinates[iVar]>numBins[iVar]) {
766 if (isOutside) continue;
768 Int_t neighbours = 0; // number of neighbours to average with
770 for (Int_t iVar=0; iVar<nDimensions; iVar++) {
771 if (coordinates[iVar] > 1) { // must not be on low edge border
772 coordinates[iVar]-- ; //get lower neighbouring bin
773 content += copy->GetBinContent(coordinates);
774 error2 += TMath::Power(copy->GetBinError(coordinates),2);
776 coordinates[iVar]++ ; //back to initial coordinate
778 if (coordinates[iVar] < numBins[iVar]) { // must not be on up edge border
779 coordinates[iVar]++ ; //get upper neighbouring bin
780 content += copy->GetBinContent(coordinates);
781 error2 += TMath::Power(copy->GetBinError(coordinates),2);
783 coordinates[iVar]-- ; //back to initial coordinate
787 hist->SetBinContent(coordinates,content/(1.+neighbours));
788 hist->SetBinError (coordinates,TMath::Sqrt(error2)/(1.+neighbours));
791 delete [] coordinates ;
796 //______________________________________________________________
798 Short_t AliCFUnfolding::SmoothUsingFunction() {
800 // Fits the unfolded spectrum using the function fSmoothFunction
803 AliDebug(0,Form("Smooth function is a %s with option \"%s\" and has %d parameters : ",fSmoothFunction->ClassName(),fSmoothOption,fSmoothFunction->GetNpar()));
805 for (Int_t iPar=0; iPar<fSmoothFunction->GetNpar(); iPar++) AliDebug(0,Form("par[%d]=%e",iPar,fSmoothFunction->GetParameter(iPar)));
809 switch (fNVariables) {
810 case 1 : fitResult = fUnfolded->Projection(0) ->Fit(fSmoothFunction,fSmoothOption); break;
811 case 2 : fitResult = fUnfolded->Projection(1,0) ->Fit(fSmoothFunction,fSmoothOption); break; // (1,0) instead of (0,1) -> TAxis issue
812 case 3 : fitResult = fUnfolded->Projection(0,1,2)->Fit(fSmoothFunction,fSmoothOption); break;
813 default: AliFatal(Form("Cannot handle such fit in %d dimensions",fNVariables)) ; return 1;
816 if (fitResult != 0) {
817 AliWarning(Form("Fit failed with status %d, stopping the loop",fitResult));
821 Int_t nDim = fNVariables;
822 Int_t* bins = new Int_t[nDim]; // number of bins for each variable
823 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
825 for (Int_t iVar=0; iVar<nDim; iVar++) {
826 bins[iVar] = fUnfolded->GetAxis(iVar)->GetNbins();
830 Int_t *bin = new Int_t[nDim]; // bin to fill the THnSparse (holding the bin coordinates)
831 Double_t x[3] = {0,0,0} ; // value in bin center (max dimension is 3 (TF3))
833 // loop on the bins and update of fUnfolded
834 // THnSparse::Multiply(TF1*) doesn't exist, so let's do it bin by bin
835 for (Long_t iBin=0; iBin<nBins; iBin++) {
836 Long_t bin_tmp = iBin ;
837 for (Int_t iVar=0; iVar<nDim; iVar++) {
838 bin[iVar] = 1 + bin_tmp % bins[iVar] ;
839 bin_tmp /= bins[iVar] ;
840 x[iVar] = fUnfolded->GetAxis(iVar)->GetBinCenter(bin[iVar]);
842 Double_t functionValue = fSmoothFunction->Eval(x[0],x[1],x[2]) ;
843 fUnfolded->SetBinError (bin,fUnfolded->GetBinError(bin)*functionValue/fUnfolded->GetBinContent(bin));
844 fUnfolded->SetBinContent(bin,functionValue);
851 //______________________________________________________________
853 void AliCFUnfolding::CreateFlatPrior() {
855 // Creates a flat prior distribution
858 AliInfo("Creating a flat a priori distribution");
860 // create the frame of the THnSparse given (for example) the one from the efficiency map
861 fPrior = (THnSparse*) fEfficiency->Clone();
863 if (fNVariables != fPrior->GetNdimensions())
864 AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions()));
866 Int_t nDim = fNVariables;
867 Int_t* bins = new Int_t[nDim]; // number of bins for each variable
868 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
870 for (Int_t iVar=0; iVar<nDim; iVar++) {
871 bins[iVar] = fPrior->GetAxis(iVar)->GetNbins();
875 Int_t *bin = new Int_t[nDim]; // bin to fill the THnSparse (holding the bin coordinates)
877 // loop that sets 1 in each bin
878 for (Long_t iBin=0; iBin<nBins; iBin++) {
879 Long_t bin_tmp = iBin ;
880 for (Int_t iVar=0; iVar<nDim; iVar++) {
881 bin[iVar] = 1 + bin_tmp % bins[iVar] ;
882 bin_tmp /= bins[iVar] ;
884 fPrior->SetBinContent(bin,1.); // put 1 everywhere
885 fPrior->SetBinError (bin,0.); // put 0 everywhere
888 fOriginalPrior = (THnSparse*)fPrior->Clone();