Update the mult corr histograms
[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 // Currently the similarity criterion is the Chi2 between the a priori //
41 // and the unfolded spectrum.                                          //
42 //                                                                     //
43 // Currently the user has to define the max. number of iterations      //
44 // (::SetMaxNumberOfIterations)                                        //
45 // and the chi2 below which the procedure will stop                    //
46 // (::SetMaxChi2 or ::SetMaxChi2PerDOF)                                //
47 //                                                                     //
48 // An optional possibility is to smooth the unfolded spectrum at the   //
49 // end of each iteration, either using a fit function                  //
50 // (only if #dimensions <=3)                                           //
51 // or a simple averaging using the neighbouring bins values.           //
52 // This is possible calling the function ::UseSmoothing                //
53 // If no argument is passed to this function, then the second option   //
54 // is used.                                                            //
55 //                                                                     //
56 // IMPORTANT:                                                          //
57 //-----------                                                          //
58 // With this approach, the efficiency map must be calculated           //
59 // with *simulated* values only, otherwise the method won't work.      //
60 //                                                                     //
61 // ex: efficiency(bin_pt) = number_rec(bin_pt) / number_sim(bin_pt)    //
62 //                                                                     //
63 // the pt bin "bin_pt" must always be the same in both the efficiency  //
64 // numerator and denominator.                                          //
65 // This is why the efficiency map has to be created by a method        //
66 // from which both reconstructed and simulated values are accessible   //
67 // simultaneously.                                                     //
68 //                                                                     //
69 //                                                                     //
70 //---------------------------------------------------------------------//
71 // Author : renaud.vernet@cern.ch                                      //
72 //---------------------------------------------------------------------//
73
74
75 #include "AliCFUnfolding.h"
76 #include "TMath.h"
77 #include "TAxis.h"
78 #include "AliLog.h"
79 #include "TF1.h"
80 #include "TH1D.h"
81 #include "TH2D.h"
82 #include "TH3D.h"
83
84 ClassImp(AliCFUnfolding)
85
86 //______________________________________________________________
87
88 AliCFUnfolding::AliCFUnfolding() :
89   TNamed(),
90   fResponse(0x0),
91   fPrior(0x0),
92   fEfficiency(0x0),
93   fMeasured(0x0),
94   fMaxNumIterations(0),
95   fNVariables(0),
96   fMaxChi2(0),
97   fUseSmoothing(kFALSE),
98   fSmoothFunction(0x0),
99   fSmoothOption(""),
100   fOriginalPrior(0x0),
101   fInverseResponse(0x0),
102   fMeasuredEstimate(0x0),
103   fConditional(0x0),
104   fProjResponseInT(0x0),
105   fUnfolded(0x0),
106   fCoordinates2N(0x0),
107   fCoordinatesN_M(0x0),
108   fCoordinatesN_T(0x0)
109 {
110   //
111   // default constructor
112   //
113 }
114
115 //______________________________________________________________
116
117 AliCFUnfolding::AliCFUnfolding(const Char_t* name, const Char_t* title, const Int_t nVar, 
118                                const THnSparse* response, const THnSparse* efficiency, const THnSparse* measured, const THnSparse* prior) :
119   TNamed(name,title),
120   fResponse((THnSparse*)response->Clone()),
121   fPrior(0x0),
122   fEfficiency((THnSparse*)efficiency->Clone()),
123   fMeasured((THnSparse*)measured->Clone()),
124   fMaxNumIterations(0),
125   fNVariables(nVar),
126   fMaxChi2(0),
127   fUseSmoothing(kFALSE),
128   fSmoothFunction(0x0),
129   fSmoothOption(""),
130   fOriginalPrior(0x0),
131   fInverseResponse(0x0),
132   fMeasuredEstimate(0x0),
133   fConditional(0x0),
134   fProjResponseInT(0x0),
135   fUnfolded(0x0),
136   fCoordinates2N(0x0),
137   fCoordinatesN_M(0x0),
138   fCoordinatesN_T(0x0)
139 {
140   //
141   // named constructor
142   //
143
144   AliInfo(Form("\n\n--------------------------\nCreating an unfolder :\n--------------------------\nresponse matrix has %d dimension(s)",fResponse->GetNdimensions()));
145   
146   if (!prior) CreateFlatPrior(); // if no prior distribution declared, simply use a flat distribution
147   else {
148     fPrior = (THnSparse*) prior->Clone();
149     fOriginalPrior = (THnSparse*)fPrior->Clone();
150     if (fPrior->GetNdimensions() != fNVariables) 
151       AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions()));
152   }
153
154   if (fEfficiency->GetNdimensions() != fNVariables) 
155     AliFatal(Form("The efficiency matrix should have %d dimensions, and it has actually %d",fNVariables,fEfficiency->GetNdimensions()));
156   if (fMeasured->GetNdimensions() != fNVariables) 
157     AliFatal(Form("The measured matrix should have %d dimensions, and it has actually %d",fNVariables,fMeasured->GetNdimensions()));
158   if (fResponse->GetNdimensions() != 2*fNVariables) 
159     AliFatal(Form("The response matrix should have %d dimensions, and it has actually %d",2*fNVariables,fResponse->GetNdimensions()));
160   
161
162   for (Int_t iVar=0; iVar<fNVariables; iVar++) {
163     AliInfo(Form("prior      matrix has %d bins in dimension %d",fPrior     ->GetAxis(iVar)->GetNbins(),iVar));
164     AliInfo(Form("efficiency matrix has %d bins in dimension %d",fEfficiency->GetAxis(iVar)->GetNbins(),iVar));
165     AliInfo(Form("measured   matrix has %d bins in dimension %d",fMeasured  ->GetAxis(iVar)->GetNbins(),iVar));
166   }
167   Init();
168 }
169
170
171 //______________________________________________________________
172
173 AliCFUnfolding::AliCFUnfolding(const AliCFUnfolding& c) :
174   TNamed(c),
175   fResponse((THnSparse*)c.fResponse->Clone()),
176   fPrior((THnSparse*)c.fPrior->Clone()),
177   fEfficiency((THnSparse*)c.fEfficiency->Clone()),
178   fMeasured((THnSparse*)c.fMeasured->Clone()),
179   fMaxNumIterations(c.fMaxNumIterations),
180   fNVariables(c.fNVariables),
181   fMaxChi2(c.fMaxChi2),
182   fUseSmoothing(c.fUseSmoothing),
183   fSmoothFunction((TF1*)c.fSmoothFunction->Clone()),
184   fSmoothOption(fSmoothOption),
185   fOriginalPrior((THnSparse*)c.fOriginalPrior->Clone()),
186   fInverseResponse((THnSparse*)c.fInverseResponse->Clone()),
187   fMeasuredEstimate((THnSparse*)fMeasuredEstimate->Clone()),
188   fConditional((THnSparse*)c.fConditional->Clone()),
189   fProjResponseInT((THnSparse*)c.fProjResponseInT->Clone()),
190   fUnfolded((THnSparse*)c.fUnfolded->Clone()),
191   fCoordinates2N(new Int_t(*c.fCoordinates2N)),
192   fCoordinatesN_M(new Int_t(*c.fCoordinatesN_M)),
193   fCoordinatesN_T(new Int_t(*c.fCoordinatesN_T))
194 {
195   //
196   // copy constructor
197   //
198 }
199
200 //______________________________________________________________
201
202 AliCFUnfolding& AliCFUnfolding::operator=(const AliCFUnfolding& c) {
203   //
204   // assignment operator
205   //
206   
207   if (this!=&c) {
208     TNamed::operator=(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     fMaxNumIterations = c.fMaxNumIterations ;
214     fNVariables = c.fNVariables ;
215     fMaxChi2 = c.fMaxChi2 ;
216     fUseSmoothing = c.fUseSmoothing ;
217     fSmoothFunction = (TF1*)c.fSmoothFunction->Clone();
218     fSmoothOption = c.fSmoothOption ;
219     fOriginalPrior = (THnSparse*)c.fOriginalPrior->Clone() ;
220     fInverseResponse = (THnSparse*)c.fInverseResponse->Clone() ;
221     fMeasuredEstimate = (THnSparse*)fMeasuredEstimate->Clone() ;
222     fConditional = (THnSparse*)c.fConditional->Clone() ;
223     fProjResponseInT = (THnSparse*)c.fProjResponseInT->Clone() ;
224     fUnfolded = (THnSparse*)c.fUnfolded->Clone() ;
225     fCoordinates2N  = new Int_t(*c.fCoordinates2N)  ;
226     fCoordinatesN_M = new Int_t(*c.fCoordinatesN_M) ;
227     fCoordinatesN_T = new Int_t(*c.fCoordinatesN_T) ;
228   }
229   return *this;
230 }
231
232 //______________________________________________________________
233
234 AliCFUnfolding::~AliCFUnfolding() {
235   //
236   // destructor
237   //
238   if (fResponse)           delete fResponse;
239   if (fPrior)              delete fPrior;
240   if (fEfficiency)         delete fEfficiency;
241   if (fMeasured)           delete fMeasured;
242   if (fSmoothFunction)     delete fSmoothFunction;
243   if (fOriginalPrior)      delete fOriginalPrior;
244   if (fInverseResponse)    delete fInverseResponse;
245   if (fMeasuredEstimate)   delete fMeasuredEstimate;
246   if (fConditional)        delete fConditional;
247   if (fProjResponseInT)    delete fProjResponseInT;
248   if (fCoordinates2N)      delete [] fCoordinates2N; 
249   if (fCoordinatesN_M)     delete [] fCoordinatesN_M; 
250   if (fCoordinatesN_T)     delete [] fCoordinatesN_T; 
251 }
252
253 //______________________________________________________________
254
255 void AliCFUnfolding::Init() {
256   //
257   // initialisation function : creates internal settings
258   //
259
260   fCoordinates2N  = new Int_t[2*fNVariables];
261   fCoordinatesN_M = new Int_t[fNVariables];
262   fCoordinatesN_T = new Int_t[fNVariables];
263
264   // create the matrix of conditional probabilities P(M|T)
265   CreateConditional();
266   
267   // create the frame of the inverse response matrix
268   fInverseResponse  = (THnSparse*) fResponse->Clone();
269   // create the frame of the unfolded spectrum
270   fUnfolded = (THnSparse*) fPrior->Clone();
271   // create the frame of the measurement estimate spectrum
272   fMeasuredEstimate = (THnSparse*) fMeasured->Clone();
273 }
274
275 //______________________________________________________________
276
277 void AliCFUnfolding::CreateEstMeasured() {
278   //
279   // This function creates a estimate (M) of the reconstructed spectrum 
280   // given the a priori distribution (T), the efficiency (E) and the conditional matrix (COND)
281   //
282   // --> P(M) = SUM   { P(M|T)    * P(T) }
283   // --> M(i) = SUM_k { COND(i,k) * T(k) * E (k)}
284   //
285   // This is needed to calculate the inverse response matrix
286   //
287
288
289   // clean the measured estimate spectrum
290   fMeasuredEstimate->Reset();
291
292   THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone();
293   priorTimesEff->Multiply(fEfficiency);
294
295   // fill it
296   for (Long_t iBin=0; iBin<fConditional->GetNbins(); iBin++) {
297     Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N);
298     Double_t conditionalError = fConditional->GetBinError  (iBin);
299     GetCoordinates();
300     Double_t priorTimesEffValue = priorTimesEff->GetBinContent(fCoordinatesN_T);
301     Double_t priorTimesEffError = priorTimesEff->GetBinError  (fCoordinatesN_T);
302     Double_t fill = conditionalValue * priorTimesEffValue ;
303     
304     if (fill>0.) {
305       // error calculation : gaussian error propagation (may be overestimated...)
306       Double_t err2  = TMath::Power(fMeasuredEstimate->GetBinError(fCoordinatesN_M),2) ;
307       err2 += TMath::Power(conditionalValue*priorTimesEffError,2) + TMath::Power(conditionalError*priorTimesEffValue,2) ;
308       Double_t err = TMath::Sqrt(err2);
309       fMeasuredEstimate->SetBinError(fCoordinatesN_M,err);
310
311       fMeasuredEstimate->AddBinContent(fCoordinatesN_M,fill);
312     }
313   }
314   delete priorTimesEff ;
315 }
316
317 //______________________________________________________________
318
319 void AliCFUnfolding::CreateInvResponse() {
320   //
321   // Creates the inverse response matrix (INV) with Bayesian method
322   //  : uses the conditional matrix (COND), the prior probabilities (T) and the efficiency map (E)
323   //
324   // --> P(T|M)   = P(M|T)    * P(T) * eff(T) / SUM   { P(M|T)    * P(T) }
325   // --> INV(i,j) = COND(i,j) * T(j) * E(j)   / SUM_k { COND(i,k) * T(k) }
326   //
327
328   THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone();
329   priorTimesEff->Multiply(fEfficiency);
330
331   for (Long_t iBin=0; iBin<fConditional->GetNbins(); iBin++) {
332     Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N);
333     Double_t conditionalError = fConditional->GetBinError  (iBin);
334     GetCoordinates();
335     Double_t estMeasuredValue   = fMeasuredEstimate->GetBinContent(fCoordinatesN_M);
336     Double_t estMeasuredError   = fMeasuredEstimate->GetBinError  (fCoordinatesN_M);
337     Double_t priorTimesEffValue = priorTimesEff    ->GetBinContent(fCoordinatesN_T);
338     Double_t priorTimesEffError = priorTimesEff    ->GetBinError  (fCoordinatesN_T);
339     Double_t fill = (estMeasuredValue>0. ? conditionalValue * priorTimesEffValue / estMeasuredValue : 0. ) ;
340     // error calculation : gaussian error propagation (may be overestimated...)
341     Double_t err  = 0. ;
342     if (estMeasuredValue>0.) {
343       err = TMath::Sqrt( TMath::Power(conditionalError * priorTimesEffValue * estMeasuredValue ,2) +
344                          TMath::Power(conditionalValue * priorTimesEffError * estMeasuredValue ,2) + 
345                          TMath::Power(conditionalValue * priorTimesEffValue * estMeasuredError ,2) )
346         / TMath::Power(estMeasuredValue,2) ;
347     }
348     if (fill>0. || fInverseResponse->GetBinContent(fCoordinates2N)>0.) {
349       fInverseResponse->SetBinContent(fCoordinates2N,fill);
350       fInverseResponse->SetBinError  (fCoordinates2N,err );
351     }
352   } 
353   delete priorTimesEff ;
354 }
355
356 //______________________________________________________________
357
358 void AliCFUnfolding::Unfold() {
359   //
360   // Main routine called by the user : 
361   // it calculates the unfolded spectrum from the response matrix and the measured spectrum
362   // several iterations are performed until a reasonable chi2 is reached
363   //
364
365   Int_t iIterBayes=0 ;
366   Double_t chi2=0 ;
367
368   for (iIterBayes=0; iIterBayes<fMaxNumIterations; iIterBayes++) { // bayes iterations
369     CreateEstMeasured();
370     CreateInvResponse();
371     CreateUnfolded();
372     chi2 = GetChi2();
373     AliDebug(0,Form("Chi2 at iteration %d is %e",iIterBayes,chi2));
374     if (fMaxChi2>0. && chi2<fMaxChi2) {
375       break;
376     }
377     // update the prior distribution
378     if (fUseSmoothing) {
379       if (Smooth()) {
380         AliError("Couldn't smooth the unfolded spectrum!!");
381         AliInfo(Form("\n\n=======================\nFinished at iteration %d : Chi2 is %e and you required it to be < %e\n=======================\n\n",iIterBayes,chi2,fMaxChi2));
382         return;
383       }
384     }
385     fPrior = (THnSparse*)fUnfolded->Clone() ; // this should be changed (memory)
386   }
387   AliInfo(Form("\n\n=======================\nFinished at iteration %d : Chi2 is %e and you required it to be < %e\n=======================\n\n",iIterBayes,chi2,fMaxChi2));
388 }
389
390 //______________________________________________________________
391
392 void AliCFUnfolding::CreateUnfolded() {
393   //
394   // Creates the unfolded (T) spectrum from the measured spectrum (M) and the inverse response matrix (INV)
395   // We have P(T) = SUM   { P(T|M)   * P(M) } 
396   //   -->   T(i) = SUM_k { INV(i,k) * M(k) }
397   //
398
399
400   // clear the unfolded spectrum
401   fUnfolded->Reset();
402   
403   for (Long_t iBin=0; iBin<fInverseResponse->GetNbins(); iBin++) {
404     Double_t invResponseValue = fInverseResponse->GetBinContent(iBin,fCoordinates2N);
405     Double_t invResponseError = fInverseResponse->GetBinError  (iBin);
406     GetCoordinates();
407     Double_t effValue      = fEfficiency->GetBinContent(fCoordinatesN_T);
408     Double_t effError      = fEfficiency->GetBinError  (fCoordinatesN_T);
409     Double_t measuredValue = fMeasured  ->GetBinContent(fCoordinatesN_M);
410     Double_t measuredError = fMeasured  ->GetBinError  (fCoordinatesN_M);
411     Double_t fill = (effValue>0. ? invResponseValue * measuredValue / effValue : 0.) ;
412     
413     if (fill>0.) {
414       // error calculation : gaussian error propagation (may be overestimated...)
415       Double_t err2 = TMath::Power(fUnfolded->GetBinError(fCoordinatesN_T),2) ;
416       err2 += TMath::Power(invResponseError * measuredValue * effValue,2) / TMath::Power(effValue,4) ;
417       err2 += TMath::Power(invResponseValue * measuredError * effValue,2) / TMath::Power(effValue,4) ;
418       err2 += TMath::Power(invResponseValue * measuredValue * effError,2) / TMath::Power(effValue,4) ;
419       Double_t err = TMath::Sqrt(err2);
420       fUnfolded->SetBinError(fCoordinatesN_T,err);
421
422       fUnfolded->AddBinContent(fCoordinatesN_T,fill);
423     }
424   }
425 }
426     
427 //______________________________________________________________
428
429 void AliCFUnfolding::GetCoordinates() {
430   //
431   // assign coordinates in Measured and True spaces (dim=N) from coordinates in global space (dim=2N)
432   //
433   for (Int_t i = 0; i<fNVariables ; i++) {
434     fCoordinatesN_M[i] = fCoordinates2N[i];
435     fCoordinatesN_T[i] = fCoordinates2N[i+fNVariables];
436   }
437 }
438
439 //______________________________________________________________
440
441 void AliCFUnfolding::CreateConditional() {
442   //
443   // creates the conditional probability matrix (R*) holding the P(M|T), given the reponse matrix R
444   //
445   //  --> R*(i,j) = R(i,j) / SUM_k{ R(k,j) }
446   //
447
448   fConditional     = (THnSparse*) fResponse->Clone(); // output of this function
449   fProjResponseInT = (THnSparse*) fPrior->Clone();    // output denominator : 
450                                                       // projection of the response matrix on the TRUE axis
451   Int_t* dim = new Int_t [fNVariables];
452   for (Int_t iDim=0; iDim<fNVariables; iDim++) dim[iDim] = fNVariables+iDim ; //dimensions corresponding to TRUE values (i.e. from N to 2N-1)
453   fProjResponseInT = fConditional->Projection(fNVariables,dim,"E"); //project
454   delete [] dim; 
455   
456   // fill the conditional probability matrix
457   for (Long_t iBin=0; iBin<fResponse->GetNbins(); iBin++) {
458     Double_t responseValue = fResponse->GetBinContent(iBin,fCoordinates2N);
459     Double_t responseError = fResponse->GetBinError  (iBin);
460     GetCoordinates();
461     Double_t projValue = fProjResponseInT->GetBinContent(fCoordinatesN_T);
462     Double_t projError = fProjResponseInT->GetBinError  (fCoordinatesN_T);
463     
464     Double_t fill = responseValue / projValue ;
465     if (fill>0. || fConditional->GetBinContent(fCoordinates2N)>0.) {
466       fConditional->SetBinContent(fCoordinates2N,fill);
467       // gaussian error for the moment
468       Double_t err2 = TMath::Power(responseError*projValue,2) + TMath::Power(responseValue*projError,2) ;
469       Double_t err = TMath::Sqrt(err2);
470       err /= TMath::Power(projValue,2) ;
471       fConditional->SetBinError  (fCoordinates2N,err);
472     }
473   }
474 }
475
476 //______________________________________________________________
477
478 Double_t AliCFUnfolding::GetChi2() {
479   //
480   // Returns the chi2 between unfolded and a priori spectrum
481   //
482
483   Double_t chi2 = 0. ;
484   for (Long_t iBin=0; iBin<fPrior->GetNbins(); iBin++) {
485     Double_t priorValue = fPrior->GetBinContent(iBin,fCoordinatesN_T);
486     Double_t error_unf  = fUnfolded->GetBinError(fCoordinatesN_T);
487     chi2 += (error_unf > 0. ? TMath::Power((fUnfolded->GetBinContent(fCoordinatesN_T) - priorValue)/error_unf,2) / priorValue : 0.) ;
488   }
489   return chi2;
490 }
491
492 //______________________________________________________________
493
494 void AliCFUnfolding::SetMaxChi2PerDOF(Double_t val) {
495   //
496   // Max. chi2 per degree of freedom : user setting
497   //
498
499   Int_t nDOF = 1 ;
500   for (Int_t iDim=0; iDim<fNVariables; iDim++) {
501     nDOF *= fPrior->GetAxis(iDim)->GetNbins();
502   }
503   AliInfo(Form("Number of degrees of freedom = %d",nDOF));
504   fMaxChi2 = val * nDOF ;
505 }
506
507 //______________________________________________________________
508
509 Short_t AliCFUnfolding::Smooth() {
510   //
511   // Smoothes the unfolded spectrum
512   //
513   // By default each cell content is replaced by the average with the neighbouring bins (but not diagonally-neighbouring bins)
514   // However, if a specific function fcn has been defined in UseSmoothing(fcn), the unfolded will be fit and updated using fcn 
515   //
516   
517   if (fSmoothFunction) {
518     AliDebug(2,Form("Smoothing spectrum with fit function %p",fSmoothFunction));
519     return SmoothUsingFunction();
520   }
521   else return SmoothUsingNeighbours(fUnfolded);
522 }
523
524 //______________________________________________________________
525
526 Short_t AliCFUnfolding::SmoothUsingNeighbours(THnSparse* hist) {
527   //
528   // Smoothes the unfolded spectrum using neighouring bins
529   //
530
531   Int_t  const nDimensions = hist->GetNdimensions() ;
532   Int_t* coordinates = new Int_t[nDimensions];
533
534   Int_t* numBins = new Int_t[nDimensions];
535   for (Int_t iVar=0; iVar<nDimensions; iVar++) numBins[iVar] = hist->GetAxis(iVar)->GetNbins();
536   
537   //need a copy because hist will be updated during the loop, and this creates problems
538   THnSparse* copy = (THnSparse*)hist->Clone();
539
540   for (Long_t iBin=0; iBin<copy->GetNbins(); iBin++) { //loop on non-empty bins
541     Double_t content = copy->GetBinContent(iBin,coordinates);
542     Double_t error2  = TMath::Power(copy->GetBinError(iBin),2);
543
544     printf("coord = [%d,%d]\n",coordinates[0],coordinates[1]);
545
546     // skip the under/overflow bins...
547     Bool_t isOutside = kFALSE ;
548     for (Int_t iVar=0; iVar<nDimensions; iVar++) {
549       if (coordinates[iVar]<1 || coordinates[iVar]>numBins[iVar]) {
550         isOutside=kTRUE;
551         break;
552       }
553     }
554     if (isOutside) continue;
555     
556     Int_t neighbours = 0; // number of neighbours to average with
557
558     for (Int_t iVar=0; iVar<nDimensions; iVar++) {
559       if (coordinates[iVar] > 1) { // must not be on low edge border
560         coordinates[iVar]-- ; //get lower neighbouring bin 
561         content += copy->GetBinContent(coordinates);
562         error2  += TMath::Power(copy->GetBinError(coordinates),2);
563         neighbours++;
564         coordinates[iVar]++ ; //back to initial coordinate
565       }
566       if (coordinates[iVar] < numBins[iVar]) { // must not be on up edge border
567         coordinates[iVar]++ ; //get upper neighbouring bin
568         content += copy->GetBinContent(coordinates);
569         error2  += TMath::Power(copy->GetBinError(coordinates),2);
570         neighbours++;
571         coordinates[iVar]-- ; //back to initial coordinate
572       }
573     }
574     // make an average
575     hist->SetBinContent(coordinates,content/(1.+neighbours));
576     hist->SetBinError  (coordinates,TMath::Sqrt(error2)/(1.+neighbours));
577   }
578   delete [] numBins;
579   delete [] coordinates ;
580   delete copy;
581   return 0;
582 }
583
584 //______________________________________________________________
585
586 Short_t AliCFUnfolding::SmoothUsingFunction() {
587   //
588   // Fits the unfolded spectrum using the function fSmoothFunction
589   //
590
591   AliDebug(0,Form("Smooth function is a %s with option \"%s\" and has %d parameters : ",fSmoothFunction->ClassName(),fSmoothOption,fSmoothFunction->GetNpar()));
592
593   for (Int_t iPar=0; iPar<fSmoothFunction->GetNpar(); iPar++) AliDebug(0,Form("par[%d]=%e",iPar,fSmoothFunction->GetParameter(iPar)));
594
595   Int_t fitResult = 0;
596
597   switch (fNVariables) {
598   case 1 : fitResult = fUnfolded->Projection(0)    ->Fit(fSmoothFunction,fSmoothOption); break;
599   case 2 : fitResult = fUnfolded->Projection(1,0)  ->Fit(fSmoothFunction,fSmoothOption); break; // (1,0) instead of (0,1) -> TAxis issue
600   case 3 : fitResult = fUnfolded->Projection(0,1,2)->Fit(fSmoothFunction,fSmoothOption); break;
601   default: AliFatal(Form("Cannot handle such fit in %d dimensions",fNVariables)) ; return 1;
602   }
603
604   if (fitResult != 0) {
605     AliWarning(Form("Fit failed with status %d, stopping the loop",fitResult));
606     return 1;
607   }
608
609   Int_t nDim = fNVariables;
610   Int_t* bins = new Int_t[nDim]; // number of bins for each variable
611   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
612
613   for (Int_t iVar=0; iVar<nDim; iVar++) {
614     bins[iVar] = fUnfolded->GetAxis(iVar)->GetNbins();
615     nBins *= bins[iVar];
616   }
617
618   Int_t *bin  = new Int_t[nDim];    // bin to fill the THnSparse (holding the bin coordinates)
619   Double_t x[3] = {0,0,0} ;         // value in bin center (max dimension is 3 (TF3))
620
621   // loop on the bins and update of fUnfolded
622   // THnSparse::Multiply(TF1*) doesn't exist, so let's do it bin by bin
623   for (Long_t iBin=0; iBin<nBins; iBin++) {
624     Long_t bin_tmp = iBin ;
625     for (Int_t iVar=0; iVar<nDim; iVar++) {
626       bin[iVar] = 1 + bin_tmp % bins[iVar] ;
627       bin_tmp /= bins[iVar] ;
628       x[iVar] = fUnfolded->GetAxis(iVar)->GetBinCenter(bin[iVar]);
629     }
630     Double_t functionValue = fSmoothFunction->Eval(x[0],x[1],x[2]) ;
631     fUnfolded->SetBinError  (bin,fUnfolded->GetBinError(bin)*functionValue/fUnfolded->GetBinContent(bin));
632     fUnfolded->SetBinContent(bin,functionValue);
633   }
634   return 0;
635 }
636
637 //______________________________________________________________
638
639 void AliCFUnfolding::CreateFlatPrior() {
640   //
641   // Creates a flat prior distribution
642   // 
643
644   AliInfo("Creating a flat a priori distribution");
645   
646   // create the frame of the THnSparse given (for example) the one from the efficiency map
647   fPrior = (THnSparse*) fEfficiency->Clone();
648
649   if (fNVariables != fPrior->GetNdimensions()) 
650     AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions()));
651
652   Int_t nDim = fNVariables;
653   Int_t* bins = new Int_t[nDim]; // number of bins for each variable
654   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
655
656   for (Int_t iVar=0; iVar<nDim; iVar++) {
657     bins[iVar] = fPrior->GetAxis(iVar)->GetNbins();
658     nBins *= bins[iVar];
659   }
660
661   Int_t *bin = new Int_t[nDim]; // bin to fill the THnSparse (holding the bin coordinates)
662
663   // loop that sets 1 in each bin
664   for (Long_t iBin=0; iBin<nBins; iBin++) {
665     Long_t bin_tmp = iBin ;
666     for (Int_t iVar=0; iVar<nDim; iVar++) {
667       bin[iVar] = 1 + bin_tmp % bins[iVar] ;
668       bin_tmp /= bins[iVar] ;
669     }
670     fPrior->SetBinContent(bin,1.); // put 1 everywhere
671     fPrior->SetBinError  (bin,0.); // put 0 everywhere
672   }
673   
674   fOriginalPrior = (THnSparse*)fPrior->Clone();
675
676   delete [] bin;
677   delete [] bins;
678 }