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