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 // Chi2 calculation became absolute with the correlated error //
42 // Errors on the unfolded distribution are not known until the end //
43 // Use the convergence criterion instead //
45 // Currently the user has to define the max. number of iterations //
46 // (::SetMaxNumberOfIterations) //
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); //
53 // Correlated error calculation can be activated by using: //
54 // SetUseCorrelatedErrors(Bool_t b) in combination with convergence //
56 // Documentation about correlated error calculation method can be //
57 // found in AliCFUnfolding::CalculateCorrelatedErrors() //
58 // Author: marta.verweij@cern.ch //
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 //
70 // With this approach, the efficiency map must be calculated //
71 // with *simulated* values only, otherwise the method won't work. //
73 // ex: efficiency(bin_pt) = number_rec(bin_pt) / number_sim(bin_pt) //
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 //
82 //---------------------------------------------------------------------//
83 // Author : renaud.vernet@cern.ch //
84 //---------------------------------------------------------------------//
87 #include "AliCFUnfolding.h"
96 ClassImp(AliCFUnfolding)
98 //______________________________________________________________
100 AliCFUnfolding::AliCFUnfolding() :
107 fMaxNumIterations(0),
109 fUseSmoothing(kFALSE),
110 fSmoothFunction(0x0),
111 fSmoothOption("iremn"),
113 fNRandomIterations(0),
115 fInverseResponse(0x0),
116 fMeasuredEstimate(0x0),
121 fCoordinatesN_M(0x0),
122 fCoordinatesN_T(0x0),
123 fRandomizedDist(0x0),
125 fDeltaUnfoldedP(0x0),
126 fDeltaUnfoldedN(0x0),
131 // default constructor
135 //______________________________________________________________
137 AliCFUnfolding::AliCFUnfolding(const Char_t* name, const Char_t* title, const Int_t nVar,
138 const THnSparse* response, const THnSparse* efficiency, const THnSparse* measured, const THnSparse* prior ,
139 Double_t maxConvergencePerDOF, UInt_t randomSeed, Int_t maxNumIterations
142 fResponse((THnSparse*)response->Clone()),
144 fEfficiency((THnSparse*)efficiency->Clone()),
145 fMeasured((THnSparse*)measured->Clone()),
146 fMeasuredOrig((THnSparse*)measured->Clone()),
147 fMaxNumIterations(maxNumIterations),
149 fUseSmoothing(kFALSE),
150 fSmoothFunction(0x0),
151 fSmoothOption("iremn"),
153 fNRandomIterations(maxNumIterations),
155 fInverseResponse(0x0),
156 fMeasuredEstimate(0x0),
161 fCoordinatesN_M(0x0),
162 fCoordinatesN_T(0x0),
163 fRandomizedDist(0x0),
165 fDeltaUnfoldedP(0x0),
166 fDeltaUnfoldedN(0x0),
168 fRandomSeed(randomSeed)
174 AliInfo(Form("\n\n--------------------------\nCreating an unfolder :\n--------------------------\nresponse matrix has %d dimension(s)",fResponse->GetNdimensions()));
176 if (!prior) CreateFlatPrior(); // if no prior distribution declared, simply use a flat distribution
178 fPrior = (THnSparse*) prior->Clone();
179 fOriginalPrior = (THnSparse*)fPrior->Clone();
180 if (fPrior->GetNdimensions() != fNVariables)
181 AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions()));
184 if (fEfficiency->GetNdimensions() != fNVariables)
185 AliFatal(Form("The efficiency matrix should have %d dimensions, and it has actually %d",fNVariables,fEfficiency->GetNdimensions()));
186 if (fMeasured->GetNdimensions() != fNVariables)
187 AliFatal(Form("The measured matrix should have %d dimensions, and it has actually %d",fNVariables,fMeasured->GetNdimensions()));
188 if (fResponse->GetNdimensions() != 2*fNVariables)
189 AliFatal(Form("The response matrix should have %d dimensions, and it has actually %d",2*fNVariables,fResponse->GetNdimensions()));
192 for (Int_t iVar=0; iVar<fNVariables; iVar++) {
193 AliInfo(Form("prior matrix has %d bins in dimension %d",fPrior ->GetAxis(iVar)->GetNbins(),iVar));
194 AliInfo(Form("efficiency matrix has %d bins in dimension %d",fEfficiency->GetAxis(iVar)->GetNbins(),iVar));
195 AliInfo(Form("measured matrix has %d bins in dimension %d",fMeasured ->GetAxis(iVar)->GetNbins(),iVar));
198 fRandomizedDist = (THnSparse*) fMeasuredOrig->Clone();
199 fRandomizedDist->SetTitle("Randomized");
200 SetMaxConvergencePerDOF(maxConvergencePerDOF) ;
205 //______________________________________________________________
207 AliCFUnfolding::AliCFUnfolding(const AliCFUnfolding& c) :
209 fResponse((THnSparse*)c.fResponse->Clone()),
210 fPrior((THnSparse*)c.fPrior->Clone()),
211 fEfficiency((THnSparse*)c.fEfficiency->Clone()),
212 fMeasured((THnSparse*)c.fMeasured->Clone()),
213 fMeasuredOrig((THnSparse*)c.fMeasuredOrig->Clone()),
214 fMaxNumIterations(c.fMaxNumIterations),
215 fNVariables(c.fNVariables),
216 fUseSmoothing(c.fUseSmoothing),
217 fSmoothFunction((TF1*)c.fSmoothFunction->Clone()),
218 fSmoothOption(c.fSmoothOption),
219 fMaxConvergence(c.fMaxConvergence),
220 fNRandomIterations(c.fNRandomIterations),
221 fOriginalPrior((THnSparse*)c.fOriginalPrior->Clone()),
222 fInverseResponse((THnSparse*)c.fInverseResponse->Clone()),
223 fMeasuredEstimate((THnSparse*)fMeasuredEstimate->Clone()),
224 fConditional((THnSparse*)c.fConditional->Clone()),
225 fUnfolded((THnSparse*)c.fUnfolded->Clone()),
226 fUnfoldedFinal((THnSparse*)c.fUnfoldedFinal->Clone()),
227 fCoordinates2N(new Int_t(*c.fCoordinates2N)),
228 fCoordinatesN_M(new Int_t(*c.fCoordinatesN_M)),
229 fCoordinatesN_T(new Int_t(*c.fCoordinatesN_T)),
230 fRandomizedDist((THnSparse*)c.fRandomizedDist->Clone()),
231 fRandom3((TRandom3*)c.fRandom3->Clone()),
232 fDeltaUnfoldedP((THnSparse*)c.fDeltaUnfoldedP),
233 fDeltaUnfoldedN((THnSparse*)c.fDeltaUnfoldedN),
234 fNCalcCorrErrors(c.fNCalcCorrErrors),
235 fRandomSeed(c.fRandomSeed)
242 //______________________________________________________________
244 AliCFUnfolding& AliCFUnfolding::operator=(const AliCFUnfolding& c) {
246 // assignment operator
250 TNamed::operator=(c);
251 fResponse = (THnSparse*)c.fResponse->Clone() ;
252 fPrior = (THnSparse*)c.fPrior->Clone() ;
253 fEfficiency = (THnSparse*)c.fEfficiency->Clone() ;
254 fMeasured = (THnSparse*)c.fMeasured->Clone() ;
255 fMeasuredOrig = ((THnSparse*)c.fMeasuredOrig->Clone()),
256 fMaxNumIterations = c.fMaxNumIterations ;
257 fNVariables = c.fNVariables ;
258 fMaxConvergence = c.fMaxConvergence ;
259 fUseSmoothing = c.fUseSmoothing ;
260 fSmoothFunction = (TF1*)c.fSmoothFunction->Clone();
261 fSmoothOption = c.fSmoothOption ;
262 fNRandomIterations = c.fNRandomIterations;
263 fOriginalPrior = (THnSparse*)c.fOriginalPrior->Clone() ;
264 fInverseResponse = (THnSparse*)c.fInverseResponse->Clone() ;
265 fMeasuredEstimate = (THnSparse*)fMeasuredEstimate->Clone() ;
266 fConditional = (THnSparse*)c.fConditional->Clone() ;
267 fUnfolded = (THnSparse*)c.fUnfolded->Clone() ;
268 fUnfoldedFinal = (THnSparse*)c.fUnfoldedFinal->Clone() ;
269 fCoordinates2N = new Int_t(*c.fCoordinates2N) ;
270 fCoordinatesN_M = new Int_t(*c.fCoordinatesN_M) ;
271 fCoordinatesN_T = new Int_t(*c.fCoordinatesN_T) ;
272 fRandomizedDist = (THnSparse*)c.fRandomizedDist->Clone();
273 fRandom3 = (TRandom3*)c.fRandom3->Clone();
274 fDeltaUnfoldedP = (THnSparse*)c.fDeltaUnfoldedP;
275 fDeltaUnfoldedN = (THnSparse*)c.fDeltaUnfoldedN;
276 fNCalcCorrErrors = c.fNCalcCorrErrors ;
277 fRandomSeed = c.fRandomSeed ;
282 //______________________________________________________________
284 AliCFUnfolding::~AliCFUnfolding() {
288 if (fResponse) delete fResponse;
289 if (fPrior) delete fPrior;
290 if (fEfficiency) delete fEfficiency;
291 if (fMeasured) delete fMeasured;
292 if (fMeasuredOrig) delete fMeasuredOrig;
293 if (fSmoothFunction) delete fSmoothFunction;
294 if (fOriginalPrior) delete fOriginalPrior;
295 if (fInverseResponse) delete fInverseResponse;
296 if (fMeasuredEstimate) delete fMeasuredEstimate;
297 if (fConditional) delete fConditional;
298 if (fUnfolded) delete fUnfolded;
299 if (fUnfoldedFinal) delete fUnfoldedFinal;
300 if (fCoordinates2N) delete [] fCoordinates2N;
301 if (fCoordinatesN_M) delete [] fCoordinatesN_M;
302 if (fCoordinatesN_T) delete [] fCoordinatesN_T;
303 if (fRandomizedDist) delete fRandomizedDist;
304 if (fRandom3) delete fRandom3;
305 if (fDeltaUnfoldedP) delete fDeltaUnfoldedP;
306 if (fDeltaUnfoldedN) delete fDeltaUnfoldedN;
309 //______________________________________________________________
311 void AliCFUnfolding::Init() {
313 // initialisation function : creates internal settings
316 fRandom3 = new TRandom3(fRandomSeed);
318 fCoordinates2N = new Int_t[2*fNVariables];
319 fCoordinatesN_M = new Int_t[fNVariables];
320 fCoordinatesN_T = new Int_t[fNVariables];
322 // create the matrix of conditional probabilities P(M|T)
323 CreateConditional(); //done only once at initialization
325 // create the frame of the inverse response matrix
326 fInverseResponse = (THnSparse*) fResponse->Clone();
327 // create the frame of the unfolded spectrum
328 fUnfolded = (THnSparse*) fPrior->Clone();
329 fUnfolded->SetTitle("Unfolded");
330 // create the frame of the measurement estimate spectrum
331 fMeasuredEstimate = (THnSparse*) fMeasured->Clone();
332 // create the frame of the original measurement spectrum
333 fMeasuredOrig = (THnSparse*) fMeasured->Clone();
335 fDeltaUnfoldedP = (THnSparse*)fPrior->Clone();
336 fDeltaUnfoldedP->SetTitle("#Delta unfolded");
337 fDeltaUnfoldedP->Reset();
338 fDeltaUnfoldedN = (THnSparse*)fPrior->Clone();
339 fDeltaUnfoldedP->SetTitle("");
340 fDeltaUnfoldedN->Reset();
344 //______________________________________________________________
346 void AliCFUnfolding::CreateEstMeasured() {
348 // This function creates a estimate (M) of the reconstructed spectrum
349 // given the a priori distribution (T), the efficiency (E) and the conditional matrix (COND)
351 // --> P(M) = SUM { P(M|T) * P(T) }
352 // --> M(i) = SUM_k { COND(i,k) * T(k) * E (k)}
354 // This is needed to calculate the inverse response matrix
358 // clean the measured estimate spectrum
359 fMeasuredEstimate->Reset();
361 THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone();
362 priorTimesEff->Multiply(fEfficiency);
365 for (Long_t iBin=0; iBin<fConditional->GetNbins(); iBin++) {
366 Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N);
368 Double_t priorTimesEffValue = priorTimesEff->GetBinContent(fCoordinatesN_T);
369 Double_t fill = conditionalValue * priorTimesEffValue ;
372 fMeasuredEstimate->AddBinContent(fCoordinatesN_M,fill);
373 fMeasuredEstimate->SetBinError(fCoordinatesN_M,0.);
376 delete priorTimesEff ;
379 //______________________________________________________________
381 void AliCFUnfolding::CreateInvResponse() {
383 // Creates the inverse response matrix (INV) with Bayesian method
384 // : uses the conditional matrix (COND), the prior probabilities (T) and the efficiency map (E)
386 // --> P(T|M) = P(M|T) * P(T) * eff(T) / SUM { P(M|T) * P(T) }
387 // --> INV(i,j) = COND(i,j) * T(j) * E(j) / SUM_k { COND(i,k) * T(k) }
390 THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone();
391 priorTimesEff->Multiply(fEfficiency);
393 for (Long_t iBin=0; iBin<fConditional->GetNbins(); iBin++) {
394 Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N);
396 Double_t estMeasuredValue = fMeasuredEstimate->GetBinContent(fCoordinatesN_M);
397 Double_t priorTimesEffValue = priorTimesEff ->GetBinContent(fCoordinatesN_T);
398 Double_t fill = (estMeasuredValue>0. ? conditionalValue * priorTimesEffValue / estMeasuredValue : 0. ) ;
399 if (fill>0. || fInverseResponse->GetBinContent(fCoordinates2N)>0.) {
400 fInverseResponse->SetBinContent(fCoordinates2N,fill);
401 fInverseResponse->SetBinError (fCoordinates2N,0.);
404 delete priorTimesEff ;
407 //______________________________________________________________
409 void AliCFUnfolding::Unfold() {
411 // Main routine called by the user :
412 // it calculates the unfolded spectrum from the response matrix, measured spectrum and efficiency
413 // several iterations are performed until a reasonable chi2 or convergence criterion is reached
416 Int_t iIterBayes = 0 ;
417 Double_t convergence = 0.;
419 for (iIterBayes=0; iIterBayes<fMaxNumIterations; iIterBayes++) { // bayes iterations
421 CreateEstMeasured(); // create measured estimate from prior
422 CreateInvResponse(); // create inverse response from prior
423 CreateUnfolded(); // create unfoled spectrum from measured and inverse response
425 convergence = GetConvergence();
426 AliDebug(0,Form("convergence at iteration %d is %e",iIterBayes,convergence));
428 if (fMaxConvergence>0. && convergence<fMaxConvergence && fNCalcCorrErrors == 0) {
429 fNRandomIterations = iIterBayes;
430 AliDebug(0,Form("convergence is met at iteration %d",iIterBayes));
436 AliError("Couldn't smooth the unfolded spectrum!!");
437 if (fNCalcCorrErrors>0) {
438 AliInfo(Form("=======================\nUnfold of randomized distribution finished at iteration %d with convergence %e \n",iIterBayes,convergence));
441 AliInfo(Form("\n\n=======================\nFinish at iteration %d : convergence is %e and you required it to be < %e\n=======================\n\n",iIterBayes,convergence,fMaxConvergence));
447 // update the prior distribution
448 if (fPrior) delete fPrior ;
449 fPrior = (THnSparse*)fUnfolded->Clone() ;
450 fPrior->SetTitle("Prior");
452 } // end bayes iteration
454 if (fNCalcCorrErrors==0) fUnfoldedFinal = (THnSparse*) fUnfolded->Clone() ;
457 //for (Long_t iBin=0; iBin<fUnfoldedFinal->GetNbins(); iBin++) AliDebug(2,Form("%e\n",fUnfoldedFinal->GetBinError(iBin)));
460 if (fNCalcCorrErrors == 0) {
461 AliInfo("\n================================================\nFinished bayes iteration, now calculating errors...\n================================================\n");
462 fNCalcCorrErrors = 1;
463 CalculateCorrelatedErrors();
466 if (fNCalcCorrErrors >1 ) {
467 AliInfo(Form("\n\n=======================\nFinished at iteration %d : convergence is %e and you required it to be < %e\n=======================\n\n",iIterBayes,convergence,fMaxConvergence));
469 else if(fNCalcCorrErrors>0) {
470 AliInfo(Form("=======================\nUnfolding of randomized distribution finished at iteration %d with convergence %e \n",iIterBayes,convergence));
474 //______________________________________________________________
476 void AliCFUnfolding::CreateUnfolded() {
478 // Creates the unfolded (T) spectrum from the measured spectrum (M) and the inverse response matrix (INV)
479 // We have P(T) = SUM { P(T|M) * P(M) }
480 // --> T(i) = SUM_k { INV(i,k) * M(k) }
484 // clear the unfolded spectrum
485 // if in the process of error calculation, the random unfolded spectrum is created
486 // otherwise the normal unfolded spectrum is created
490 for (Long_t iBin=0; iBin<fInverseResponse->GetNbins(); iBin++) {
491 Double_t invResponseValue = fInverseResponse->GetBinContent(iBin,fCoordinates2N);
493 Double_t effValue = fEfficiency->GetBinContent(fCoordinatesN_T);
494 Double_t measuredValue = fMeasured ->GetBinContent(fCoordinatesN_M);
495 Double_t fill = (effValue>0. ? invResponseValue * measuredValue / effValue : 0.) ;
498 // set errors to zero
499 // true errors will be filled afterwards
501 fUnfolded->SetBinError (fCoordinatesN_T,err);
502 fUnfolded->AddBinContent(fCoordinatesN_T,fill);
507 //______________________________________________________________
509 void AliCFUnfolding::CalculateCorrelatedErrors() {
511 // Step 1: Create randomized distribution (fRandomizedDist) of each bin of
512 // the measured spectrum to calculate correlated errors.
513 // Poisson statistics: mean = measured value of bin
514 // Step 2: Unfold randomized distribution
515 // Step 3: Store difference of unfolded spectrum from measured distribution and
516 // unfolded distribution from randomized distribution
517 // -> fDeltaUnfoldedP (TProfile with option "S")
518 // Step 4: Repeat Step 1-3 several times (fNRandomIterations)
519 // Step 5: The spread of fDeltaUnfoldedP for each bin is the error on the unfolded spectrum of that specific bin
522 //Do fNRandomIterations = bayes iterations performed
523 for (int i=0; i<fNRandomIterations; i++) {
525 // reset prior to original one
526 if (fPrior) delete fPrior ;
527 fPrior = (THnSparse*) fOriginalPrior->Clone();
529 // create randomized distribution and stick measured spectrum to it
530 CreateRandomizedDist();
531 if (fMeasured) delete fMeasured ;
532 fMeasured = (THnSparse*) fRandomizedDist->Clone();
533 fMeasured->SetTitle("Measured");
535 //unfold fRandomizedDist
537 FillDeltaUnfoldedProfile();
540 // Get statistical errors for final unfolded spectrum
541 // ie. spread of each pt bin in fDeltaUnfoldedP
544 for (Long_t iBin=0; iBin<fUnfoldedFinal->GetNbins(); iBin++) {
545 dummy = fUnfoldedFinal->GetBinContent(iBin,fCoordinatesN_M);
546 sigma = fDeltaUnfoldedP->GetBinError(fCoordinatesN_M);
547 //AliDebug(2,Form("filling error %e\n",sigma));
548 fUnfoldedFinal->SetBinError(fCoordinatesN_M,sigma);
551 // now errors are calculated
552 fNCalcCorrErrors = 2;
555 //______________________________________________________________
556 void AliCFUnfolding::CreateRandomizedDist() {
558 // Create randomized dist from original measured distribution
559 // This distribution is created several times, each time with a different random number
562 Double_t random = 0.;
563 Double_t measuredValue = 0.;
564 Double_t measuredError = 0.;
566 for (Long_t iBin=0; iBin<fRandomizedDist->GetNbins(); iBin++) {
567 measuredValue = fMeasuredOrig->GetBinContent(iBin,fCoordinatesN_M); //used as mean
568 measuredError = fMeasuredOrig->GetBinError(fCoordinatesN_M); //used as sigma
569 // random = fRandom3->PoissonD(measuredValue); //doesn't work for normalized spectra, use Gaus (assuming raw counts in bin is large >10)
570 random = fRandom3->Gaus(measuredValue,measuredError);
571 fRandomizedDist->SetBinContent(iBin,random);
575 //______________________________________________________________
576 void AliCFUnfolding::FillDeltaUnfoldedProfile() {
578 // Store difference of unfolded spectrum from measured distribution and unfolded spectrum from randomized distribution
579 // The delta profile has been set to a THnSparse to handle N dimension
580 // The THnSparse contains in each bin the mean value and spread of the difference
581 // This function updates the profile wrt to its previous mean and error
582 // The relation between iterations (n+1) and n is as follows :
583 // mean_{n+1} = (n*mean_n + value_{n+1}) / (n+1)
584 // sigma_{n+1} = sqrt { 1/(n+1) * [ n*sigma_n^2 + (n^2+n)*(mean_{n+1}-mean_n)^2 ] } (can this be optimized?)
586 for (Long_t iBin=0; iBin<fUnfolded->GetNbins(); iBin++) {
587 Double_t deltaInBin = fUnfoldedFinal->GetBinContent(iBin,fCoordinatesN_M) - fUnfolded->GetBinContent(iBin);
588 Double_t entriesInBin = fDeltaUnfoldedN->GetBinContent(fCoordinatesN_M);
589 //AliDebug(2,Form("%e %e ==> delta = %e\n",fUnfoldedFinal->GetBinContent(iBin,fCoordinatesN_M),fUnfolded->GetBinContent(iBin),deltaInBin));
591 Double_t mean_n = fDeltaUnfoldedP->GetBinContent(fCoordinatesN_M) ;
592 Double_t mean_nplus1 = mean_n ;
593 mean_nplus1 *= entriesInBin ;
594 mean_nplus1 += deltaInBin ;
595 mean_nplus1 /= (entriesInBin+1) ;
597 Double_t sigma = fDeltaUnfoldedP->GetBinContent(fCoordinatesN_M) ;
599 sigma *= entriesInBin ;
600 sigma += ( (entriesInBin*entriesInBin+entriesInBin) * TMath::Power(mean_nplus1 - mean_n,2) ) ;
601 sigma /= (entriesInBin+1) ;
602 sigma = TMath::Sqrt(sigma) ;
604 //AliDebug(2,Form("sigma = %e\n",sigma));
606 fDeltaUnfoldedP->SetBinContent(fCoordinatesN_M,mean_nplus1) ;
607 fDeltaUnfoldedP->SetBinError (fCoordinatesN_M,sigma) ;
608 fDeltaUnfoldedN->SetBinContent(fCoordinatesN_M,entriesInBin+1);
612 //______________________________________________________________
614 void AliCFUnfolding::GetCoordinates() {
616 // assign coordinates in Measured and True spaces (dim=N) from coordinates in global space (dim=2N)
618 for (Int_t i = 0; i<fNVariables ; i++) {
619 fCoordinatesN_M[i] = fCoordinates2N[i];
620 fCoordinatesN_T[i] = fCoordinates2N[i+fNVariables];
624 //______________________________________________________________
626 void AliCFUnfolding::CreateConditional() {
628 // creates the conditional probability matrix (R*) holding the P(M|T), given the reponse matrix R
630 // --> R*(i,j) = R(i,j) / SUM_k{ R(k,j) }
633 fConditional = (THnSparse*) fResponse->Clone(); // output of this function
635 Int_t* dim = new Int_t [fNVariables];
636 for (Int_t iDim=0; iDim<fNVariables; iDim++) dim[iDim] = fNVariables+iDim ; //dimensions corresponding to TRUE values (i.e. from N to 2N-1)
638 THnSparse* responseInT = fConditional->Projection(fNVariables,dim,"E"); // output denominator :
639 // projection of the response matrix on the TRUE axis
642 // fill the conditional probability matrix
643 for (Long_t iBin=0; iBin<fResponse->GetNbins(); iBin++) {
644 Double_t responseValue = fResponse->GetBinContent(iBin,fCoordinates2N);
646 Double_t projValue = responseInT->GetBinContent(fCoordinatesN_T);
648 Double_t fill = responseValue / projValue ;
649 if (fill>0. || fConditional->GetBinContent(fCoordinates2N)>0.) {
650 fConditional->SetBinContent(fCoordinates2N,fill);
652 fConditional->SetBinError (fCoordinates2N,err);
657 //______________________________________________________________
659 Int_t AliCFUnfolding::GetDOF() {
661 // number of dof = number of bins
665 for (Int_t iDim=0; iDim<fNVariables; iDim++) {
666 nDOF *= fPrior->GetAxis(iDim)->GetNbins();
668 AliDebug(0,Form("Number of degrees of freedom = %d",nDOF));
672 //______________________________________________________________
674 Double_t AliCFUnfolding::GetChi2() {
676 // Returns the chi2 between unfolded and a priori spectrum
677 // This function became absolute with the correlated error calculation.
678 // Errors on the unfolded distribution are not known until the end
679 // Use the convergence criterion instead
683 Double_t error_unf = 0.;
684 for (Long_t iBin=0; iBin<fPrior->GetNbins(); iBin++) {
685 Double_t priorValue = fPrior->GetBinContent(iBin,fCoordinatesN_T);
686 error_unf = fUnfolded->GetBinError(fCoordinatesN_T);
687 chi2 += (error_unf > 0. ? TMath::Power((fUnfolded->GetBinContent(fCoordinatesN_T) - priorValue)/error_unf,2) / priorValue : 0.) ;
692 //______________________________________________________________
694 Double_t AliCFUnfolding::GetConvergence() {
696 // Returns convergence criterion = \sum_t ((U_t^{n-1}-U_t^n)/U_t^{n-1})^2
697 // U is unfolded spectrum, t is the bin, n = current, n-1 = previous
699 Double_t convergence = 0.;
700 Double_t priorValue = 0.;
701 Double_t currentValue = 0.;
702 for (Long_t iBin=0; iBin < fPrior->GetNbins(); iBin++) {
703 priorValue = fPrior->GetBinContent(iBin,fCoordinatesN_T);
704 currentValue = fUnfolded->GetBinContent(fCoordinatesN_T);
707 convergence += ((priorValue-currentValue)/priorValue)*((priorValue-currentValue)/priorValue);
709 AliWarning(Form("priorValue = %f. Adding 0 to convergence criterion.",priorValue));
714 //______________________________________________________________
716 void AliCFUnfolding::SetMaxConvergencePerDOF(Double_t val) {
718 // Max. convergence criterion per degree of freedom : user setting
719 // convergence criterion = DOF*val; DOF = number of bins
720 // In Jan-Fiete's multiplicity note: Convergence criterion = DOF*0.001^2
723 Int_t nDOF = GetDOF() ;
724 fMaxConvergence = val * nDOF ;
725 AliInfo(Form("MaxConvergence = %e. Number of degrees of freedom = %d",fMaxConvergence,nDOF));
728 //______________________________________________________________
730 Short_t AliCFUnfolding::Smooth() {
732 // Smoothes the unfolded spectrum
734 // By default each cell content is replaced by the average with the neighbouring bins (but not diagonally-neighbouring bins)
735 // However, if a specific function fcn has been defined in UseSmoothing(fcn), the unfolded will be fit and updated using fcn
738 if (fSmoothFunction) {
739 AliDebug(2,Form("Smoothing spectrum with fit function %p",fSmoothFunction));
740 return SmoothUsingFunction();
742 else return SmoothUsingNeighbours(fUnfolded);
745 //______________________________________________________________
747 Short_t AliCFUnfolding::SmoothUsingNeighbours(THnSparse* hist) {
749 // Smoothes the unfolded spectrum using neighouring bins
752 Int_t const nDimensions = hist->GetNdimensions() ;
753 Int_t* coordinates = new Int_t[nDimensions];
755 Int_t* numBins = new Int_t[nDimensions];
756 for (Int_t iVar=0; iVar<nDimensions; iVar++) numBins[iVar] = hist->GetAxis(iVar)->GetNbins();
758 //need a copy because hist will be updated during the loop, and this creates problems
759 THnSparse* copy = (THnSparse*)hist->Clone();
761 for (Long_t iBin=0; iBin<copy->GetNbins(); iBin++) { //loop on non-empty bins
762 Double_t content = copy->GetBinContent(iBin,coordinates);
763 Double_t error2 = TMath::Power(copy->GetBinError(iBin),2);
765 // skip the under/overflow bins...
766 Bool_t isOutside = kFALSE ;
767 for (Int_t iVar=0; iVar<nDimensions; iVar++) {
768 if (coordinates[iVar]<1 || coordinates[iVar]>numBins[iVar]) {
773 if (isOutside) continue;
775 Int_t neighbours = 0; // number of neighbours to average with
777 for (Int_t iVar=0; iVar<nDimensions; iVar++) {
778 if (coordinates[iVar] > 1) { // must not be on low edge border
779 coordinates[iVar]-- ; //get lower neighbouring bin
780 content += copy->GetBinContent(coordinates);
781 error2 += TMath::Power(copy->GetBinError(coordinates),2);
783 coordinates[iVar]++ ; //back to initial coordinate
785 if (coordinates[iVar] < numBins[iVar]) { // must not be on up edge border
786 coordinates[iVar]++ ; //get upper neighbouring bin
787 content += copy->GetBinContent(coordinates);
788 error2 += TMath::Power(copy->GetBinError(coordinates),2);
790 coordinates[iVar]-- ; //back to initial coordinate
794 hist->SetBinContent(coordinates,content/(1.+neighbours));
795 hist->SetBinError (coordinates,TMath::Sqrt(error2)/(1.+neighbours));
798 delete [] coordinates ;
803 //______________________________________________________________
805 Short_t AliCFUnfolding::SmoothUsingFunction() {
807 // Fits the unfolded spectrum using the function fSmoothFunction
810 AliDebug(0,Form("Smooth function is a %s with option \"%s\" and has %d parameters : ",fSmoothFunction->ClassName(),fSmoothOption,fSmoothFunction->GetNpar()));
812 for (Int_t iPar=0; iPar<fSmoothFunction->GetNpar(); iPar++) AliDebug(0,Form("par[%d]=%e",iPar,fSmoothFunction->GetParameter(iPar)));
816 switch (fNVariables) {
817 case 1 : fitResult = fUnfolded->Projection(0) ->Fit(fSmoothFunction,fSmoothOption); break;
818 case 2 : fitResult = fUnfolded->Projection(1,0) ->Fit(fSmoothFunction,fSmoothOption); break; // (1,0) instead of (0,1) -> TAxis issue
819 case 3 : fitResult = fUnfolded->Projection(0,1,2)->Fit(fSmoothFunction,fSmoothOption); break;
820 default: AliFatal(Form("Cannot handle such fit in %d dimensions",fNVariables)) ; return 1;
823 if (fitResult != 0) {
824 AliWarning(Form("Fit failed with status %d, stopping the loop",fitResult));
828 Int_t nDim = fNVariables;
829 Int_t* bins = new Int_t[nDim]; // number of bins for each variable
830 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
832 for (Int_t iVar=0; iVar<nDim; iVar++) {
833 bins[iVar] = fUnfolded->GetAxis(iVar)->GetNbins();
837 Int_t *bin = new Int_t[nDim]; // bin to fill the THnSparse (holding the bin coordinates)
838 Double_t x[3] = {0,0,0} ; // value in bin center (max dimension is 3 (TF3))
840 // loop on the bins and update of fUnfolded
841 // THnSparse::Multiply(TF1*) doesn't exist, so let's do it bin by bin
842 for (Long_t iBin=0; iBin<nBins; iBin++) {
843 Long_t bin_tmp = iBin ;
844 for (Int_t iVar=0; iVar<nDim; iVar++) {
845 bin[iVar] = 1 + bin_tmp % bins[iVar] ;
846 bin_tmp /= bins[iVar] ;
847 x[iVar] = fUnfolded->GetAxis(iVar)->GetBinCenter(bin[iVar]);
849 Double_t functionValue = fSmoothFunction->Eval(x[0],x[1],x[2]) ;
850 fUnfolded->SetBinError (bin,fUnfolded->GetBinError(bin)*functionValue/fUnfolded->GetBinContent(bin));
851 fUnfolded->SetBinContent(bin,functionValue);
858 //______________________________________________________________
860 void AliCFUnfolding::CreateFlatPrior() {
862 // Creates a flat prior distribution
865 AliInfo("Creating a flat a priori distribution");
867 // create the frame of the THnSparse given (for example) the one from the efficiency map
868 fPrior = (THnSparse*) fEfficiency->Clone();
869 fPrior->SetTitle("Prior");
871 if (fNVariables != fPrior->GetNdimensions())
872 AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions()));
874 Int_t nDim = fNVariables;
875 Int_t* bins = new Int_t[nDim]; // number of bins for each variable
876 Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
878 for (Int_t iVar=0; iVar<nDim; iVar++) {
879 bins[iVar] = fPrior->GetAxis(iVar)->GetNbins();
883 Int_t *bin = new Int_t[nDim]; // bin to fill the THnSparse (holding the bin coordinates)
885 // loop that sets 1 in each bin
886 for (Long_t iBin=0; iBin<nBins; iBin++) {
887 Long_t bin_tmp = iBin ;
888 for (Int_t iVar=0; iVar<nDim; iVar++) {
889 bin[iVar] = 1 + bin_tmp % bins[iVar] ;
890 bin_tmp /= bins[iVar] ;
892 fPrior->SetBinContent(bin,1.); // put 1 everywhere
893 fPrior->SetBinError (bin,0.); // put 0 everywhere
896 fOriginalPrior = (THnSparse*)fPrior->Clone();