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