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