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