reverted change in l. 573
[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   Double_t dummy       = 0.;
486   for (Long_t iBin=0; iBin<fUnfoldedFinal->GetNbins(); iBin++) {
487     dummy = fUnfoldedFinal->GetBinContent(iBin,fCoordinatesN_M);
488     sigma = fDeltaUnfoldedP->GetBinError(fCoordinatesN_M);
489     //AliDebug(2,Form("filling error %e\n",sigma));
490     fUnfoldedFinal->SetBinError(fCoordinatesN_M,sigma);
491   }
492
493   // now errors are calculated
494   fNCalcCorrErrors = 2;
495 }
496
497 //______________________________________________________________
498 void AliCFUnfolding::CreateRandomizedDist() {
499   //
500   // Create randomized dist from original measured distribution
501   // This distribution is created several times, each time with a different random number
502   //
503
504   for (Long_t iBin=0; iBin<fResponseOrig->GetNbins(); iBin++) {
505     Double_t val = fResponseOrig->GetBinContent(iBin,fCoordinatesN_M); //used as mean
506     Double_t err = fResponseOrig->GetBinError(fCoordinatesN_M);        //used as sigma
507     Double_t ran = fRandom3->Gaus(val,err);
508     // random        = fRandom3->PoissonD(measuredValue); //doesn't work for normalized spectra, use Gaus (assuming raw counts in bin is large >10)
509     fRandomResponse->SetBinContent(iBin,ran);
510   }
511   for (Long_t iBin=0; iBin<fEfficiencyOrig->GetNbins(); iBin++) {
512     Double_t val = fEfficiencyOrig->GetBinContent(iBin,fCoordinatesN_M); //used as mean
513     Double_t err = fEfficiencyOrig->GetBinError(fCoordinatesN_M);        //used as sigma
514     Double_t ran = fRandom3->Gaus(val,err);
515     // random        = fRandom3->PoissonD(measuredValue); //doesn't work for normalized spectra, use Gaus (assuming raw counts in bin is large >10)
516     fRandomEfficiency->SetBinContent(iBin,ran);
517   }
518   for (Long_t iBin=0; iBin<fMeasuredOrig->GetNbins(); iBin++) {
519     Double_t val = fMeasuredOrig->GetBinContent(iBin,fCoordinatesN_M); //used as mean
520     Double_t err = fMeasuredOrig->GetBinError(fCoordinatesN_M);        //used as sigma
521     Double_t ran = fRandom3->Gaus(val,err);
522     // random        = fRandom3->PoissonD(measuredValue); //doesn't work for normalized spectra, use Gaus (assuming raw counts in bin is large >10)
523     fRandomMeasured->SetBinContent(iBin,ran);
524   }
525 }
526
527 //______________________________________________________________
528 void AliCFUnfolding::FillDeltaUnfoldedProfile() {
529   //
530   // Store difference of unfolded spectrum from measured distribution and unfolded spectrum from randomized distribution
531   // The delta profile has been set to a THnSparse to handle N dimension
532   // The THnSparse contains in each bin the mean value and spread of the difference 
533   // This function updates the profile wrt to its previous mean and error
534   // The relation between iterations (n+1) and n is as follows :
535   //  mean_{n+1} = (n*mean_n + value_{n+1}) / (n+1)
536   // sigma_{n+1} = sqrt { 1/(n+1) * [ n*sigma_n^2 + (n^2+n)*(mean_{n+1}-mean_n)^2 ] }    (can this be optimized?)
537
538   for (Long_t iBin=0; iBin<fUnfolded->GetNbins(); iBin++) {
539     Double_t deltaInBin   = fUnfoldedFinal->GetBinContent(iBin,fCoordinatesN_M) - fUnfolded->GetBinContent(iBin);
540     Double_t entriesInBin = fDeltaUnfoldedN->GetBinContent(fCoordinatesN_M);
541     //AliDebug(2,Form("%e %e ==> delta = %e\n",fUnfoldedFinal->GetBinContent(iBin,fCoordinatesN_M),fUnfolded->GetBinContent(iBin),deltaInBin));
542
543     Double_t mean_n = fDeltaUnfoldedP->GetBinContent(fCoordinatesN_M) ;
544     Double_t mean_nplus1 = mean_n ;
545     mean_nplus1 *= entriesInBin ;
546     mean_nplus1 += deltaInBin ;
547     mean_nplus1 /= (entriesInBin+1) ;
548
549     Double_t sigma = fDeltaUnfoldedP->GetBinContent(fCoordinatesN_M) ;
550     sigma *= sigma ;
551     sigma *= entriesInBin ;
552     sigma += ( (entriesInBin*entriesInBin+entriesInBin) * TMath::Power(mean_nplus1 - mean_n,2) ) ;
553     sigma /= (entriesInBin+1) ;
554     sigma = TMath::Sqrt(sigma) ;
555
556     //AliDebug(2,Form("sigma = %e\n",sigma));
557
558     fDeltaUnfoldedP->SetBinContent(fCoordinatesN_M,mean_nplus1) ;
559     fDeltaUnfoldedP->SetBinError  (fCoordinatesN_M,sigma) ;
560     fDeltaUnfoldedN->SetBinContent(fCoordinatesN_M,entriesInBin+1);
561   }
562 }
563
564 //______________________________________________________________
565
566 void AliCFUnfolding::GetCoordinates() {
567   //
568   // assign coordinates in Measured and True spaces (dim=N) from coordinates in global space (dim=2N)
569   //
570   for (Int_t i = 0; i<fNVariables ; i++) {
571     fCoordinatesN_M[i] = fCoordinates2N[i];
572     fCoordinatesN_T[i] = fCoordinates2N[i+fNVariables];
573   }
574 }
575
576 //______________________________________________________________
577
578 void AliCFUnfolding::CreateConditional() {
579   //
580   // creates the conditional probability matrix (R*) holding the P(M|T), given the reponse matrix R
581   //
582   //  --> R*(i,j) = R(i,j) / SUM_k{ R(k,j) }
583   //
584
585   fConditional = (THnSparse*) fResponse->Clone();  // output of this function
586
587   Int_t* dim = new Int_t [fNVariables];
588   for (Int_t iDim=0; iDim<fNVariables; iDim++) dim[iDim] = fNVariables+iDim ; //dimensions corresponding to TRUE values (i.e. from N to 2N-1)
589
590   THnSparse* responseInT = fConditional->Projection(fNVariables,dim,"E");     // output denominator : 
591                                                                               // projection of the response matrix on the TRUE axis
592   delete [] dim; 
593
594   // fill the conditional probability matrix
595   for (Long_t iBin=0; iBin<fResponse->GetNbins(); iBin++) {
596     Double_t responseValue = fResponse->GetBinContent(iBin,fCoordinates2N);
597     GetCoordinates();
598     Double_t projValue = responseInT->GetBinContent(fCoordinatesN_T);
599    
600     Double_t fill = responseValue / projValue ;
601     if (fill>0. || fConditional->GetBinContent(fCoordinates2N)>0.) {
602       fConditional->SetBinContent(fCoordinates2N,fill);
603       Double_t err = 0.;
604       fConditional->SetBinError  (fCoordinates2N,err);
605     }
606   }
607   delete responseInT ;
608 }
609 //______________________________________________________________
610
611 Int_t AliCFUnfolding::GetDOF() {
612   //
613   // number of dof = number of bins
614   //
615
616   Int_t nDOF = 1 ;
617   for (Int_t iDim=0; iDim<fNVariables; iDim++) {
618     nDOF *= fPrior->GetAxis(iDim)->GetNbins();
619   }
620   AliDebug(0,Form("Number of degrees of freedom = %d",nDOF));
621   return nDOF;
622 }
623
624 //______________________________________________________________
625
626 Double_t AliCFUnfolding::GetChi2() {
627   //
628   // Returns the chi2 between unfolded and a priori spectrum
629   // This function became absolute with the correlated error calculation.
630   // Errors on the unfolded distribution are not known until the end
631   // Use the convergence criterion instead
632   //
633
634   Double_t chi2      = 0. ;
635   Double_t error_unf = 0.;
636   for (Long_t iBin=0; iBin<fPrior->GetNbins(); iBin++) {
637     Double_t priorValue = fPrior->GetBinContent(iBin,fCoordinatesN_T);
638     error_unf  = fUnfolded->GetBinError(fCoordinatesN_T);
639     chi2 += (error_unf > 0. ? TMath::Power((fUnfolded->GetBinContent(fCoordinatesN_T) - priorValue)/error_unf,2) / priorValue : 0.) ;
640   }
641   return chi2;
642 }
643
644 //______________________________________________________________
645
646 Double_t AliCFUnfolding::GetConvergence() {
647   //
648   // Returns convergence criterion = \sum_t ((U_t^{n-1}-U_t^n)/U_t^{n-1})^2
649   // U is unfolded spectrum, t is the bin, n = current, n-1 = previous
650   //
651   Double_t convergence = 0.;
652   Double_t priorValue  = 0.;
653   Double_t currentValue = 0.;
654   for (Long_t iBin=0; iBin < fPrior->GetNbins(); iBin++) {
655     priorValue = fPrior->GetBinContent(iBin,fCoordinatesN_T);
656     currentValue = fUnfolded->GetBinContent(fCoordinatesN_T);
657
658     if (priorValue > 0.)
659       convergence += ((priorValue-currentValue)/priorValue)*((priorValue-currentValue)/priorValue);
660     else 
661       AliWarning(Form("priorValue = %f. Adding 0 to convergence criterion.",priorValue)); 
662   }
663   return convergence;
664 }
665
666 //______________________________________________________________
667
668 void AliCFUnfolding::SetMaxConvergencePerDOF(Double_t val) {
669   //
670   // Max. convergence criterion per degree of freedom : user setting
671   // convergence criterion = DOF*val; DOF = number of bins
672   // In Jan-Fiete's multiplicity note: Convergence criterion = DOF*0.001^2
673   //
674
675   Int_t nDOF = GetDOF() ;
676   fMaxConvergence = val * nDOF ;
677   AliInfo(Form("MaxConvergence = %e. Number of degrees of freedom = %d",fMaxConvergence,nDOF));
678 }
679
680 //______________________________________________________________
681
682 Short_t AliCFUnfolding::Smooth() {
683   //
684   // Smoothes the unfolded spectrum
685   //
686   // By default each cell content is replaced by the average with the neighbouring bins (but not diagonally-neighbouring bins)
687   // However, if a specific function fcn has been defined in UseSmoothing(fcn), the unfolded will be fit and updated using fcn 
688   //
689   
690   if (fSmoothFunction) {
691     AliDebug(2,Form("Smoothing spectrum with fit function %p",fSmoothFunction));
692     return SmoothUsingFunction();
693   }
694   else return SmoothUsingNeighbours(fUnfolded);
695 }
696
697 //______________________________________________________________
698
699 Short_t AliCFUnfolding::SmoothUsingNeighbours(THnSparse* hist) {
700   //
701   // Smoothes the unfolded spectrum using neighouring bins
702   //
703
704   Int_t  const nDimensions = hist->GetNdimensions() ;
705   Int_t* coordinates = new Int_t[nDimensions];
706
707   Int_t* numBins = new Int_t[nDimensions];
708   for (Int_t iVar=0; iVar<nDimensions; iVar++) numBins[iVar] = hist->GetAxis(iVar)->GetNbins();
709   
710   //need a copy because hist will be updated during the loop, and this creates problems
711   THnSparse* copy = (THnSparse*)hist->Clone();
712
713   for (Long_t iBin=0; iBin<copy->GetNbins(); iBin++) { //loop on non-empty bins
714     Double_t content = copy->GetBinContent(iBin,coordinates);
715     Double_t error2  = TMath::Power(copy->GetBinError(iBin),2);
716
717     // skip the under/overflow bins...
718     Bool_t isOutside = kFALSE ;
719     for (Int_t iVar=0; iVar<nDimensions; iVar++) {
720       if (coordinates[iVar]<1 || coordinates[iVar]>numBins[iVar]) {
721         isOutside=kTRUE;
722         break;
723       }
724     }
725     if (isOutside) continue;
726     
727     Int_t neighbours = 0; // number of neighbours to average with
728
729     for (Int_t iVar=0; iVar<nDimensions; iVar++) {
730       if (coordinates[iVar] > 1) { // must not be on low edge border
731         coordinates[iVar]-- ; //get lower neighbouring bin 
732         content += copy->GetBinContent(coordinates);
733         error2  += TMath::Power(copy->GetBinError(coordinates),2);
734         neighbours++;
735         coordinates[iVar]++ ; //back to initial coordinate
736       }
737       if (coordinates[iVar] < numBins[iVar]) { // must not be on up edge border
738         coordinates[iVar]++ ; //get upper neighbouring bin
739         content += copy->GetBinContent(coordinates);
740         error2  += TMath::Power(copy->GetBinError(coordinates),2);
741         neighbours++;
742         coordinates[iVar]-- ; //back to initial coordinate
743       }
744     }
745     // make an average
746     hist->SetBinContent(coordinates,content/(1.+neighbours));
747     hist->SetBinError  (coordinates,TMath::Sqrt(error2)/(1.+neighbours));
748   }
749   delete [] numBins;
750   delete [] coordinates ;
751   delete copy;
752   return 0;
753 }
754
755 //______________________________________________________________
756
757 Short_t AliCFUnfolding::SmoothUsingFunction() {
758   //
759   // Fits the unfolded spectrum using the function fSmoothFunction
760   //
761
762   AliDebug(0,Form("Smooth function is a %s with option \"%s\" and has %d parameters : ",fSmoothFunction->ClassName(),fSmoothOption,fSmoothFunction->GetNpar()));
763
764   for (Int_t iPar=0; iPar<fSmoothFunction->GetNpar(); iPar++) AliDebug(0,Form("par[%d]=%e",iPar,fSmoothFunction->GetParameter(iPar)));
765
766   Int_t fitResult = 0;
767
768   switch (fNVariables) {
769   case 1 : fitResult = fUnfolded->Projection(0)    ->Fit(fSmoothFunction,fSmoothOption); break;
770   case 2 : fitResult = fUnfolded->Projection(1,0)  ->Fit(fSmoothFunction,fSmoothOption); break; // (1,0) instead of (0,1) -> TAxis issue
771   case 3 : fitResult = fUnfolded->Projection(0,1,2)->Fit(fSmoothFunction,fSmoothOption); break;
772   default: AliFatal(Form("Cannot handle such fit in %d dimensions",fNVariables)) ; return 1;
773   }
774
775   if (fitResult != 0) {
776     AliWarning(Form("Fit failed with status %d, stopping the loop",fitResult));
777     return 1;
778   }
779
780   Int_t nDim = fNVariables;
781   Int_t* bins = new Int_t[nDim]; // number of bins for each variable
782   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
783
784   for (Int_t iVar=0; iVar<nDim; iVar++) {
785     bins[iVar] = fUnfolded->GetAxis(iVar)->GetNbins();
786     nBins *= bins[iVar];
787   }
788
789   Int_t *bin  = new Int_t[nDim];    // bin to fill the THnSparse (holding the bin coordinates)
790   Double_t x[3] = {0,0,0} ;         // value in bin center (max dimension is 3 (TF3))
791
792   // loop on the bins and update of fUnfolded
793   // THnSparse::Multiply(TF1*) doesn't exist, so let's do it bin by bin
794   for (Long_t iBin=0; iBin<nBins; iBin++) {
795     Long_t bin_tmp = iBin ;
796     for (Int_t iVar=0; iVar<nDim; iVar++) {
797       bin[iVar] = 1 + bin_tmp % bins[iVar] ;
798       bin_tmp /= bins[iVar] ;
799       x[iVar] = fUnfolded->GetAxis(iVar)->GetBinCenter(bin[iVar]);
800     }
801     Double_t functionValue = fSmoothFunction->Eval(x[0],x[1],x[2]) ;
802     fUnfolded->SetBinError  (bin,fUnfolded->GetBinError(bin)*functionValue/fUnfolded->GetBinContent(bin));
803     fUnfolded->SetBinContent(bin,functionValue);
804   }
805   delete [] bins;
806   delete [] bin ;
807   return 0;
808 }
809
810 //______________________________________________________________
811
812 void AliCFUnfolding::CreateFlatPrior() {
813   //
814   // Creates a flat prior distribution
815   // 
816
817   AliInfo("Creating a flat a priori distribution");
818   
819   // create the frame of the THnSparse given (for example) the one from the efficiency map
820   fPrior = (THnSparse*) fEfficiency->Clone();
821   fPrior->SetTitle("Prior");
822
823   if (fNVariables != fPrior->GetNdimensions()) 
824     AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions()));
825
826   Int_t nDim = fNVariables;
827   Int_t* bins = new Int_t[nDim]; // number of bins for each variable
828   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
829
830   for (Int_t iVar=0; iVar<nDim; iVar++) {
831     bins[iVar] = fPrior->GetAxis(iVar)->GetNbins();
832     nBins *= bins[iVar];
833   }
834
835   Int_t *bin = new Int_t[nDim]; // bin to fill the THnSparse (holding the bin coordinates)
836
837   // loop that sets 1 in each bin
838   for (Long_t iBin=0; iBin<nBins; iBin++) {
839     Long_t bin_tmp = iBin ;
840     for (Int_t iVar=0; iVar<nDim; iVar++) {
841       bin[iVar] = 1 + bin_tmp % bins[iVar] ;
842       bin_tmp /= bins[iVar] ;
843     }
844     fPrior->SetBinContent(bin,1.); // put 1 everywhere
845     fPrior->SetBinError  (bin,0.); // put 0 everywhere
846   }
847   
848   fPriorOrig = (THnSparse*)fPrior->Clone();
849
850   delete [] bin;
851   delete [] bins;
852 }