Coverity fix
[u/mrichter/AliRoot.git] / CORRFW / AliCFUnfolding.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 //---------------------------------------------------------------------//
17 //                                                                     //
18 // AliCFUnfolding Class                                                //
19 // Class to handle general unfolding procedure                         // 
20 // For the moment only bayesian unfolding is supported                 //
21 // The next steps are to add chi2 minimisation and weighting methods   //
22 //                                                                     //
23 //                                                                     //
24 //                                                                     //
25 // Use :                                                               //
26 // -------                                                             //
27 // The Bayesian unfolding consists of several iterations.              //
28 // At each iteration, an inverse response matrix is calculated, given  //
29 // the measured spectrum, the a priori (guessed) spectrum,             //
30 // the efficiency spectrum and the response matrix.                    //
31 //                                                                     //
32 // Then at each iteration, the unfolded spectrum is calculated using   //
33 // the inverse response : the goal is to get an unfolded spectrum      //
34 // similar (according to some criterion) to the a priori one.          //
35 // If the difference is too big, another iteration is performed :      //
36 // the a priori spectrum is updated to the unfolded one from the       //
37 // previous iteration, and so on so forth, until the maximum number    //
38 // of iterations or the similarity criterion is reached.               //
39 //                                                                     //
40 // Currently the similarity criterion is the Chi2 between the a priori //
41 // and the unfolded spectrum (OBSOLETE).                               //
42 //                                                                     //
43 // Chi2 calculation became absolute with the correlated error          //
44 // calculation.                                                        //
45 // Errors on the unfolded distribution are not known until the end     //
46 // Use the convergence criterion instead                               //
47 //                                                                     //
48 // Currently the user has to define the max. number of iterations      //
49 // (::SetMaxNumberOfIterations)                                        //
50 // and                                                                 //
51 //     - the chi2 below which the procedure will stop                  //
52 // (::SetMaxChi2 or ::SetMaxChi2PerDOF)   (OBSOLETE)                   //
53 //     - the convergence criterion below which the procedure will stop //
54 // SetMaxConvergencePerDOF(Double_t val);                              //
55 //                                                                     //
56 // Correlated error calculation can be activated by using:             //
57 // SetUseCorrelatedErrors(Bool_t b) in combination with convergence    //
58 // criterion                                                           //
59 // Documentation about correlated error calculation method can be      //
60 // found in AliCFUnfolding::CalculateCorrelatedErrors()                //
61 // Author: marta.verweij@cern.ch                                       //
62 //                                                                     //
63 // An optional possibility is to smooth the unfolded spectrum at the   //
64 // end of each iteration, either using a fit function                  //
65 // (only if #dimensions <=3)                                           //
66 // or a simple averaging using the neighbouring bins values.           //
67 // This is possible calling the function ::UseSmoothing                //
68 // If no argument is passed to this function, then the second option   //
69 // is used.                                                            //
70 //                                                                     //
71 // IMPORTANT:                                                          //
72 //-----------                                                          //
73 // With this approach, the efficiency map must be calculated           //
74 // with *simulated* values only, otherwise the method won't work.      //
75 //                                                                     //
76 // ex: efficiency(bin_pt) = number_rec(bin_pt) / number_sim(bin_pt)    //
77 //                                                                     //
78 // the pt bin "bin_pt" must always be the same in both the efficiency  //
79 // numerator and denominator.                                          //
80 // This is why the efficiency map has to be created by a method        //
81 // from which both reconstructed and simulated values are accessible   //
82 // simultaneously.                                                     //
83 //                                                                     //
84 //                                                                     //
85 //---------------------------------------------------------------------//
86 // Author : renaud.vernet@cern.ch                                      //
87 //---------------------------------------------------------------------//
88
89
90 #include "AliCFUnfolding.h"
91 #include "TMath.h"
92 #include "TAxis.h"
93 #include "AliLog.h"
94 #include "TF1.h"
95 #include "TH1D.h"
96 #include "TH2D.h"
97 #include "TH3D.h"
98 #include "TProfile.h"
99 #include "TRandom3.h"
100
101 ClassImp(AliCFUnfolding)
102
103 //______________________________________________________________
104
105 AliCFUnfolding::AliCFUnfolding() :
106   TNamed(),
107   fResponse(0x0),
108   fPrior(0x0),
109   fEfficiency(0x0),
110   fMeasured(0x0),
111   fMeasuredOrig(0x0),
112   fMaxNumIterations(20),
113   fNVariables(0),
114   fUseSmoothing(kFALSE),
115   fSmoothFunction(0x0),
116   fSmoothOption(""),
117   fMaxConvergence(0),
118   fUseCorrelatedErrors(kTRUE),
119   fNRandomIterations(20),
120   fOriginalPrior(0x0),
121   fInverseResponse(0x0),
122   fMeasuredEstimate(0x0),
123   fConditional(0x0),
124   fProjResponseInT(0x0),
125   fUnfolded(0x0),
126   fCoordinates2N(0x0),
127   fCoordinatesN_M(0x0),
128   fCoordinatesN_T(0x0),
129   fRandomizedDist(0x0),
130   fRandom3(0x0),
131   fRandomUnfolded(0x0),
132   fDeltaUnfoldedP(0x0),
133   fNCalcCorrErrors(0),
134   fRandomSeed(0)
135 {
136   //
137   // default constructor
138   //
139 }
140
141 //______________________________________________________________
142
143 AliCFUnfolding::AliCFUnfolding(const Char_t* name, const Char_t* title, const Int_t nVar, 
144                                const THnSparse* response, const THnSparse* efficiency, const THnSparse* measured, const THnSparse* prior) :
145   TNamed(name,title),
146   fResponse((THnSparse*)response->Clone()),
147   fPrior(0x0),
148   fEfficiency((THnSparse*)efficiency->Clone()),
149   fMeasured((THnSparse*)measured->Clone()),
150   fMeasuredOrig((THnSparse*)measured->Clone()),
151   fMaxNumIterations(0),
152   fNVariables(nVar),
153   fUseSmoothing(kFALSE),
154   fSmoothFunction(0x0),
155   fSmoothOption(""),
156   fMaxConvergence(0),
157   fUseCorrelatedErrors(kTRUE),
158   fNRandomIterations(20),
159   fOriginalPrior(0x0),
160   fInverseResponse(0x0),
161   fMeasuredEstimate(0x0),
162   fConditional(0x0),
163   fProjResponseInT(0x0),
164   fUnfolded(0x0),
165   fCoordinates2N(0x0),
166   fCoordinatesN_M(0x0),
167   fCoordinatesN_T(0x0),
168   fRandomizedDist(0x0),
169   fRandom3(0x0),
170   fRandomUnfolded(0x0),
171   fDeltaUnfoldedP(0x0),
172   fNCalcCorrErrors(0),
173   fRandomSeed(0)
174 {
175   //
176   // named constructor
177   //
178
179   AliInfo(Form("\n\n--------------------------\nCreating an unfolder :\n--------------------------\nresponse matrix has %d dimension(s)",fResponse->GetNdimensions()));
180   
181   if (!prior) CreateFlatPrior(); // if no prior distribution declared, simply use a flat distribution
182   else {
183     fPrior = (THnSparse*) prior->Clone();
184     fOriginalPrior = (THnSparse*)fPrior->Clone();
185     if (fPrior->GetNdimensions() != fNVariables) 
186       AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions()));
187   }
188
189   if (fEfficiency->GetNdimensions() != fNVariables) 
190     AliFatal(Form("The efficiency matrix should have %d dimensions, and it has actually %d",fNVariables,fEfficiency->GetNdimensions()));
191   if (fMeasured->GetNdimensions() != fNVariables) 
192     AliFatal(Form("The measured matrix should have %d dimensions, and it has actually %d",fNVariables,fMeasured->GetNdimensions()));
193   if (fResponse->GetNdimensions() != 2*fNVariables) 
194     AliFatal(Form("The response matrix should have %d dimensions, and it has actually %d",2*fNVariables,fResponse->GetNdimensions()));
195   
196
197   for (Int_t iVar=0; iVar<fNVariables; iVar++) {
198     AliInfo(Form("prior      matrix has %d bins in dimension %d",fPrior     ->GetAxis(iVar)->GetNbins(),iVar));
199     AliInfo(Form("efficiency matrix has %d bins in dimension %d",fEfficiency->GetAxis(iVar)->GetNbins(),iVar));
200     AliInfo(Form("measured   matrix has %d bins in dimension %d",fMeasured  ->GetAxis(iVar)->GetNbins(),iVar));
201   }
202
203   Init();
204 }
205
206
207 //______________________________________________________________
208
209 AliCFUnfolding::AliCFUnfolding(const AliCFUnfolding& c) :
210   TNamed(c),
211   fResponse((THnSparse*)c.fResponse->Clone()),
212   fPrior((THnSparse*)c.fPrior->Clone()),
213   fEfficiency((THnSparse*)c.fEfficiency->Clone()),
214   fMeasured((THnSparse*)c.fMeasured->Clone()),
215   fMeasuredOrig((THnSparse*)c.fMeasuredOrig->Clone()),
216   fMaxNumIterations(c.fMaxNumIterations),
217   fNVariables(c.fNVariables),
218   fUseSmoothing(c.fUseSmoothing),
219   fSmoothFunction((TF1*)c.fSmoothFunction->Clone()),
220   fSmoothOption(c.fSmoothOption),
221   fMaxConvergence(c.fMaxConvergence),
222   fUseCorrelatedErrors(c.fUseCorrelatedErrors),
223   fNRandomIterations(c.fNRandomIterations),
224   fOriginalPrior((THnSparse*)c.fOriginalPrior->Clone()),
225   fInverseResponse((THnSparse*)c.fInverseResponse->Clone()),
226   fMeasuredEstimate((THnSparse*)fMeasuredEstimate->Clone()),
227   fConditional((THnSparse*)c.fConditional->Clone()),
228   fProjResponseInT((THnSparse*)c.fProjResponseInT->Clone()),
229   fUnfolded((THnSparse*)c.fUnfolded->Clone()),
230   fCoordinates2N(new Int_t(*c.fCoordinates2N)),
231   fCoordinatesN_M(new Int_t(*c.fCoordinatesN_M)),
232   fCoordinatesN_T(new Int_t(*c.fCoordinatesN_T)),
233   fRandomizedDist((THnSparse*)c.fRandomizedDist->Clone()),
234   fRandom3((TRandom3*)c.fRandom3->Clone()),
235   fRandomUnfolded((THnSparse*)c.fRandomUnfolded->Clone()),
236   fDeltaUnfoldedP((TProfile*)c.fDeltaUnfoldedP),
237   fNCalcCorrErrors(c.fNCalcCorrErrors),
238   fRandomSeed(c.fRandomSeed)
239 {
240   //
241   // copy constructor
242   //
243 }
244
245 //______________________________________________________________
246
247 AliCFUnfolding& AliCFUnfolding::operator=(const AliCFUnfolding& c) {
248   //
249   // assignment operator
250   //
251   
252   if (this!=&c) {
253     TNamed::operator=(c);
254     fResponse = (THnSparse*)c.fResponse->Clone() ;
255     fPrior = (THnSparse*)c.fPrior->Clone() ;
256     fEfficiency = (THnSparse*)c.fEfficiency->Clone() ;
257     fMeasured = (THnSparse*)c.fMeasured->Clone() ;
258     fMeasuredOrig = ((THnSparse*)c.fMeasuredOrig->Clone()),
259     fMaxNumIterations = c.fMaxNumIterations ;
260     fNVariables = c.fNVariables ;
261     fMaxConvergence = c.fMaxConvergence ;
262     fUseSmoothing = c.fUseSmoothing ;
263     fSmoothFunction = (TF1*)c.fSmoothFunction->Clone();
264     fSmoothOption = c.fSmoothOption ;
265     fUseCorrelatedErrors = c.fUseCorrelatedErrors;
266     fNRandomIterations = c.fNRandomIterations;
267     fOriginalPrior = (THnSparse*)c.fOriginalPrior->Clone() ;
268     fInverseResponse = (THnSparse*)c.fInverseResponse->Clone() ;
269     fMeasuredEstimate = (THnSparse*)fMeasuredEstimate->Clone() ;
270     fConditional = (THnSparse*)c.fConditional->Clone() ;
271     fProjResponseInT = (THnSparse*)c.fProjResponseInT->Clone() ;
272     fUnfolded = (THnSparse*)c.fUnfolded->Clone() ;
273     fCoordinates2N  = new Int_t(*c.fCoordinates2N)  ;
274     fCoordinatesN_M = new Int_t(*c.fCoordinatesN_M) ;
275     fCoordinatesN_T = new Int_t(*c.fCoordinatesN_T) ;
276     fRandomizedDist = (THnSparse*)c.fRandomizedDist->Clone();
277     fRandom3 = (TRandom3*)c.fRandom3->Clone();
278     fRandomUnfolded = (THnSparse*)c.fRandomUnfolded->Clone();
279     fDeltaUnfoldedP = (TProfile*)c.fDeltaUnfoldedP;
280     fNCalcCorrErrors = c.fNCalcCorrErrors;
281     fRandomSeed = c.fRandomSeed ;
282   }
283   return *this;
284 }
285
286 //______________________________________________________________
287
288 AliCFUnfolding::~AliCFUnfolding() {
289   //
290   // destructor
291   //
292   if (fResponse)           delete fResponse;
293   if (fPrior)              delete fPrior;
294   if (fEfficiency)         delete fEfficiency;
295   if (fMeasured)           delete fMeasured;
296   if (fMeasuredOrig)       delete fMeasuredOrig;
297   if (fSmoothFunction)     delete fSmoothFunction;
298   if (fOriginalPrior)      delete fOriginalPrior;
299   if (fInverseResponse)    delete fInverseResponse;
300   if (fMeasuredEstimate)   delete fMeasuredEstimate;
301   if (fConditional)        delete fConditional;
302   if (fProjResponseInT)    delete fProjResponseInT;
303   if (fCoordinates2N)      delete [] fCoordinates2N; 
304   if (fCoordinatesN_M)     delete [] fCoordinatesN_M; 
305   if (fCoordinatesN_T)     delete [] fCoordinatesN_T; 
306   if (fRandomizedDist)     delete fRandomizedDist;
307   if (fRandom3)            delete fRandom3;
308   if (fRandomUnfolded)     delete fRandomUnfolded;
309   if (fDeltaUnfoldedP)     delete fDeltaUnfoldedP;
310 }
311
312 //______________________________________________________________
313
314 void AliCFUnfolding::Init() {
315   //
316   // initialisation function : creates internal settings
317   //
318
319   fRandom3 = new TRandom3(fRandomSeed);
320
321   fCoordinates2N  = new Int_t[2*fNVariables];
322   fCoordinatesN_M = new Int_t[fNVariables];
323   fCoordinatesN_T = new Int_t[fNVariables];
324
325   // create the matrix of conditional probabilities P(M|T)
326   CreateConditional();
327   
328   // create the frame of the inverse response matrix
329   fInverseResponse  = (THnSparse*) fResponse->Clone();
330   // create the frame of the unfolded spectrum
331   fUnfolded = (THnSparse*) fPrior->Clone();
332   // create the frame of the random unfolded spectrum
333   fRandomUnfolded = (THnSparse*) fPrior->Clone();
334   // create the frame of the measurement estimate spectrum
335   fMeasuredEstimate = (THnSparse*) fMeasured->Clone();
336   // create the frame of the original measurement spectrum
337   fMeasuredOrig = (THnSparse*) fMeasured->Clone();
338
339   InitDeltaUnfoldedProfile();
340
341 }
342
343 //______________________________________________________________
344
345 void AliCFUnfolding::CreateEstMeasured() {
346   //
347   // This function creates a estimate (M) of the reconstructed spectrum 
348   // given the a priori distribution (T), the efficiency (E) and the conditional matrix (COND)
349   //
350   // --> P(M) = SUM   { P(M|T)    * P(T) }
351   // --> M(i) = SUM_k { COND(i,k) * T(k) * E (k)}
352   //
353   // This is needed to calculate the inverse response matrix
354   //
355
356
357   // clean the measured estimate spectrum
358   fMeasuredEstimate->Reset();
359
360   THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone();
361   priorTimesEff->Multiply(fEfficiency);
362
363   // fill it
364   for (Long_t iBin=0; iBin<fConditional->GetNbins(); iBin++) {
365     Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N);
366     GetCoordinates();
367     Double_t priorTimesEffValue = priorTimesEff->GetBinContent(fCoordinatesN_T);
368     Double_t fill = conditionalValue * priorTimesEffValue ;
369     
370     if (fill>0.) {
371       fMeasuredEstimate->AddBinContent(fCoordinatesN_M,fill);
372       fMeasuredEstimate->SetBinError(fCoordinatesN_M,0.);
373     }
374   }
375   delete priorTimesEff ;
376 }
377
378 //______________________________________________________________
379
380 void AliCFUnfolding::CreateInvResponse() {
381   //
382   // Creates the inverse response matrix (INV) with Bayesian method
383   //  : uses the conditional matrix (COND), the prior probabilities (T) and the efficiency map (E)
384   //
385   // --> P(T|M)   = P(M|T)    * P(T) * eff(T) / SUM   { P(M|T)    * P(T) }
386   // --> INV(i,j) = COND(i,j) * T(j) * E(j)   / SUM_k { COND(i,k) * T(k) }
387   //
388
389   THnSparse* priorTimesEff = (THnSparse*) fPrior->Clone();
390   priorTimesEff->Multiply(fEfficiency);
391
392   for (Long_t iBin=0; iBin<fConditional->GetNbins(); iBin++) {
393     Double_t conditionalValue = fConditional->GetBinContent(iBin,fCoordinates2N);
394     GetCoordinates();
395     Double_t estMeasuredValue   = fMeasuredEstimate->GetBinContent(fCoordinatesN_M);
396     Double_t priorTimesEffValue = priorTimesEff    ->GetBinContent(fCoordinatesN_T);
397     Double_t fill = (estMeasuredValue>0. ? conditionalValue * priorTimesEffValue / estMeasuredValue : 0. ) ;
398     if (fill>0. || fInverseResponse->GetBinContent(fCoordinates2N)>0.) {
399       fInverseResponse->SetBinContent(fCoordinates2N,fill);
400       fInverseResponse->SetBinError  (fCoordinates2N,0.);
401     }
402   } 
403   delete priorTimesEff ;
404 }
405
406 //______________________________________________________________
407
408 void AliCFUnfolding::Unfold() {
409   //
410   // Main routine called by the user : 
411   // it calculates the unfolded spectrum from the response matrix, measured spectrum and efficiency
412   // several iterations are performed until a reasonable chi2 or convergence criterion is reached
413   //
414
415   Int_t iIterBayes     = 0 ;
416   Double_t convergence = 0.;
417
418   for (iIterBayes=0; iIterBayes<fMaxNumIterations; iIterBayes++) { // bayes iterations
419     CreateEstMeasured();
420     CreateInvResponse();
421     CreateUnfolded();
422     if (fUseCorrelatedErrors) {
423       convergence = GetConvergence();
424       AliDebug(0,Form("convergence at iteration %d is %e",iIterBayes,convergence));
425     } 
426     else AliWarning("No errors will be calculated. Activate SetUseCorrelatedErrors(kTRUE)\n");
427
428     if (fMaxConvergence>0. && convergence<fMaxConvergence && fNCalcCorrErrors<1) {
429       fNRandomIterations=iIterBayes;
430       break;
431     }
432
433     // update the prior distribution
434     if (fUseSmoothing) {
435       if (Smooth()) {
436         AliError("Couldn't smooth the unfolded spectrum!!");
437         if (fUseCorrelatedErrors) {
438           if (fNCalcCorrErrors>0) {
439             AliInfo(Form("=======================\nUnfolding of randomized distribution finished at iteration %d with convergence %e \n",iIterBayes,convergence));
440           }
441           else {
442             AliInfo(Form("\n\n=======================\nFinished at iteration %d : convergence is %e and you required it to be < %e\n=======================\n\n",iIterBayes,convergence,fMaxConvergence));
443           }
444         }
445         return;
446       }
447     }
448     if(fNCalcCorrErrors>0) {
449       if (fPrior) delete fPrior ;
450       fPrior = (THnSparse*)fRandomUnfolded->Clone() ;
451     } 
452     else {
453       if (fPrior) delete fPrior ;
454       fPrior = (THnSparse*)fUnfolded->Clone() ;
455     }
456   }
457
458   if (fUseCorrelatedErrors && fNCalcCorrErrors==0) {
459     fNCalcCorrErrors=1;
460     CalculateCorrelatedErrors();
461   }
462
463   if (fUseCorrelatedErrors) {
464     if (fNCalcCorrErrors>1) {
465       AliInfo(Form("\n\n=======================\nFinished at iteration %d : convergence is %e and you required it to be < %e\n=======================\n\n",iIterBayes,convergence,fMaxConvergence));
466     }
467     else if(fNCalcCorrErrors>0) {
468       AliInfo(Form("=======================\nUnfolding of randomized distribution finished at iteration %d with convergence %e \n",iIterBayes,convergence));
469     }
470   }
471 }
472
473 //______________________________________________________________
474
475 void AliCFUnfolding::CreateUnfolded() {
476   //
477   // Creates the unfolded (T) spectrum from the measured spectrum (M) and the inverse response matrix (INV)
478   // We have P(T) = SUM   { P(T|M)   * P(M) } 
479   //   -->   T(i) = SUM_k { INV(i,k) * M(k) }
480   //
481
482
483   // clear the unfolded spectrum
484   if(fNCalcCorrErrors>0) {
485     //unfold randomized dist
486     fRandomUnfolded->Reset();
487   }
488   else {  
489     //unfold measured dist
490     fUnfolded->Reset();
491   }
492   
493   for (Long_t iBin=0; iBin<fInverseResponse->GetNbins(); iBin++) {
494     Double_t invResponseValue = fInverseResponse->GetBinContent(iBin,fCoordinates2N);
495     GetCoordinates();
496     Double_t effValue      = fEfficiency->GetBinContent(fCoordinatesN_T);
497     Double_t measuredValue = fMeasured  ->GetBinContent(fCoordinatesN_M);
498     Double_t fill = (effValue>0. ? invResponseValue * measuredValue / effValue : 0.) ;
499
500     if (fill>0.) {
501         
502       Double_t err = 0.;
503       
504       if(fNCalcCorrErrors>0) {
505         fRandomUnfolded->SetBinError(fCoordinatesN_T,err);
506         fRandomUnfolded->AddBinContent(fCoordinatesN_T,fill);
507       }
508       else {
509         fUnfolded->SetBinError(fCoordinatesN_T,err);
510         fUnfolded->AddBinContent(fCoordinatesN_T,fill);
511       }
512     }
513   }
514 }
515
516 //______________________________________________________________
517
518 void AliCFUnfolding::CalculateCorrelatedErrors() {
519
520   fRandomizedDist     = (THnSparse*) fMeasuredOrig->Clone(); 
521   fPrior              = (THnSparse*) fOriginalPrior->Clone();
522
523   // Step 1: Create randomized distribution (fRandomizedDist) of each bin of the measured spectrum to calculate correlated errors. Poisson statistics: mean = measured value of bin
524   // Step 2: Unfold randomized distribution
525   // Step 3: Store difference of unfolded spectrum from measured distribution and unfolded distribution from randomized distribution -> fDeltaUnfoldedP (TProfile with option "S")
526   // Step 4: Repeat Step 1-3 several times (fNRandomIterations)
527   // Step 5: The spread of fDeltaUnfoldedP for each bin is the error on the unfolded spectrum of that specific bin
528
529     //Do fNRandomIterations = bayes iterations performed
530     for(int i=0; i<fNRandomIterations; i++) {
531       if (fPrior)          delete fPrior ;
532       if (fRandomizedDist) delete fRandomizedDist ;
533       fPrior          = (THnSparse*) fOriginalPrior->Clone();
534       fRandomizedDist = (THnSparse*) fMeasuredOrig->Clone();
535       CreateRandomizedDist();
536       if (fMeasured) delete fMeasured ;
537       fMeasured = (THnSparse*) fRandomizedDist->Clone();
538       //Unfold fRandomizedDist
539       Unfold();
540       FillDeltaUnfoldedProfile();
541     }
542
543     // Get statistical errors for final unfolded spectrum
544     // ie. spread of each pt bin in fDeltaUnfoldedP
545     Double_t sigma       = 0.;
546     Double_t dummy       = 0.;
547     for (Long_t iBin=0; iBin<fRandomUnfolded->GetNbins(); iBin++) {
548       dummy = fUnfolded->GetBinContent(iBin,fCoordinatesN_M);
549       sigma = fDeltaUnfoldedP->GetBinError(fCoordinatesN_M[0]);
550       fUnfolded->SetBinError(fCoordinatesN_M,sigma);
551       fNCalcCorrErrors = 2;
552     }
553
554 }
555
556 //______________________________________________________________
557 void AliCFUnfolding::InitDeltaUnfoldedProfile() {
558   //
559   //Initialize the fDeltaUnfoldedP profile
560   //Errors will be filled with spread between unfolded measured and unfolded randomized spectra
561   //
562
563   Int_t nbinsx = fResponse->GetAxis(0)->GetNbins();
564   Double_t xbins[nbinsx];
565   for(int ix=0; ix<nbinsx; ix++) {
566     xbins[ix] = fResponse->GetAxis(0)->GetBinLowEdge(ix+1);
567   }
568   xbins[nbinsx] = fResponse->GetAxis(0)->GetBinUpEdge(nbinsx);
569   fDeltaUnfoldedP = new TProfile("fDeltaUnfoldedP","Profile of pTUnfolded with spread in error",nbinsx,xbins,"S");
570 }
571 //______________________________________________________________
572 void AliCFUnfolding::CreateRandomizedDist() {
573   //
574   // Create randomized dist from measured distribution
575   //
576
577   Double_t random = 0.;
578   Double_t measuredValue = 0.;
579   Double_t measuredError = 0.;
580   for (Long_t iBin=0; iBin<fRandomizedDist->GetNbins(); iBin++) {
581     measuredValue = fMeasuredOrig->GetBinContent(iBin,fCoordinatesN_M); //used as mean
582     measuredError = fMeasuredOrig->GetBinError(fCoordinatesN_M); //used as sigma
583     // random        = fRandom3->PoissonD(measuredValue); //doesn't work for normalized spectra, use Gaus (assuming raw counts in bin is large >10)
584     random        = fRandom3->Gaus(measuredValue,measuredError);
585     fRandomizedDist->SetBinContent(iBin,random);
586   }
587 }
588
589 //______________________________________________________________
590 void AliCFUnfolding::FillDeltaUnfoldedProfile() {
591   //
592   // Store difference of unfolded spectrum from measured distribution and unfolded distribution from randomized distribution
593   //
594
595   for (Long_t iBin2=0; iBin2<fRandomUnfolded->GetNbins(); iBin2++) {
596     Double_t delta = fUnfolded->GetBinContent(iBin2,fCoordinatesN_M)-fRandomUnfolded->GetBinContent(iBin2,fCoordinatesN_M);
597     fDeltaUnfoldedP->Fill(fDeltaUnfoldedP->GetBinCenter(fCoordinatesN_M[0]),delta);
598   }
599 }
600
601 //______________________________________________________________
602
603 void AliCFUnfolding::GetCoordinates() {
604   //
605   // assign coordinates in Measured and True spaces (dim=N) from coordinates in global space (dim=2N)
606   //
607   for (Int_t i = 0; i<fNVariables ; i++) {
608     fCoordinatesN_M[i] = fCoordinates2N[i];
609     fCoordinatesN_T[i] = fCoordinates2N[i+fNVariables];
610   }
611 }
612
613 //______________________________________________________________
614
615 void AliCFUnfolding::CreateConditional() {
616   //
617   // creates the conditional probability matrix (R*) holding the P(M|T), given the reponse matrix R
618   //
619   //  --> R*(i,j) = R(i,j) / SUM_k{ R(k,j) }
620   //
621
622   fConditional     = (THnSparse*) fResponse->Clone(); // output of this function
623   fProjResponseInT = (THnSparse*) fPrior->Clone();    // output denominator : 
624                                                       // projection of the response matrix on the TRUE axis
625   Int_t* dim = new Int_t [fNVariables];
626   for (Int_t iDim=0; iDim<fNVariables; iDim++) dim[iDim] = fNVariables+iDim ; //dimensions corresponding to TRUE values (i.e. from N to 2N-1)
627   fProjResponseInT = fConditional->Projection(fNVariables,dim,"E"); //project
628   delete [] dim; 
629
630   // fill the conditional probability matrix
631   for (Long_t iBin=0; iBin<fResponse->GetNbins(); iBin++) {
632     Double_t responseValue = fResponse->GetBinContent(iBin,fCoordinates2N);
633     GetCoordinates();
634     Double_t projValue = fProjResponseInT->GetBinContent(fCoordinatesN_T);
635    
636     Double_t fill = responseValue / projValue ;
637     if (fill>0. || fConditional->GetBinContent(fCoordinates2N)>0.) {
638       fConditional->SetBinContent(fCoordinates2N,fill);
639       Double_t err = 0.;
640       fConditional->SetBinError  (fCoordinates2N,err);
641     }
642   }
643 }
644 //______________________________________________________________
645
646 Int_t AliCFUnfolding::GetDOF() {
647   //
648   // number of dof = number of bins
649   //
650
651   Int_t nDOF = 1 ;
652   for (Int_t iDim=0; iDim<fNVariables; iDim++) {
653     nDOF *= fPrior->GetAxis(iDim)->GetNbins();
654   }
655   AliDebug(0,Form("Number of degrees of freedom = %d",nDOF));
656   return nDOF;
657 }
658
659 //______________________________________________________________
660
661 Double_t AliCFUnfolding::GetChi2() {
662   //
663   // Returns the chi2 between unfolded and a priori spectrum
664   // This function became absolute with the correlated error calculation.
665   // Errors on the unfolded distribution are not known until the end
666   // Use the convergence criterion instead
667   //
668
669   Double_t chi2      = 0. ;
670   Double_t error_unf = 0.;
671   for (Long_t iBin=0; iBin<fPrior->GetNbins(); iBin++) {
672     Double_t priorValue = fPrior->GetBinContent(iBin,fCoordinatesN_T);
673     error_unf  = fUnfolded->GetBinError(fCoordinatesN_T);
674     chi2 += (error_unf > 0. ? TMath::Power((fUnfolded->GetBinContent(fCoordinatesN_T) - priorValue)/error_unf,2) / priorValue : 0.) ;
675   }
676   return chi2;
677 }
678
679 //______________________________________________________________
680
681 Double_t AliCFUnfolding::GetConvergence() {
682   //
683   // Returns convergence criterion = \sum_t ((U_t^{n-1}-U_t^n)/U_t^{n-1})^2
684   // U is unfolded spectrum, t is the bin, n = current, n-1 = previous
685   //
686   Double_t convergence = 0.;
687   Double_t priorValue  = 0.;
688   Double_t currentValue = 0.;
689   for (Long_t iBin=0; iBin < fPrior->GetNbins(); iBin++) {
690     priorValue = fPrior->GetBinContent(iBin,fCoordinatesN_T);
691     if (fNCalcCorrErrors > 0) 
692       currentValue = fRandomUnfolded->GetBinContent(fCoordinatesN_T);
693     else 
694       currentValue = fUnfolded->GetBinContent(fCoordinatesN_T);
695
696     if (priorValue > 0.)
697       convergence += ((priorValue-currentValue)/priorValue)*((priorValue-currentValue)/priorValue);
698     else {
699       AliWarning(Form("priorValue = %f. Adding 0 to convergence criterion.",priorValue)); 
700       convergence += 0.;
701     }
702   }
703   return convergence;
704 }
705
706 //______________________________________________________________
707
708 void AliCFUnfolding::SetMaxConvergencePerDOF(Double_t val) {
709   //
710   // Max. convergence criterion per degree of freedom : user setting
711   // convergence criterion = DOF*val; DOF = number of bins
712   // In Jan-Fiete's multiplicity note: Convergence criterion = DOF*0.001^2
713   //
714
715   fUseCorrelatedErrors = kTRUE;
716   Int_t nDOF = GetDOF() ;
717   fMaxConvergence = val * nDOF ;
718   AliInfo(Form("MaxConvergence = %e. Number of degrees of freedom = %d",fMaxConvergence,nDOF));
719 }
720
721 //______________________________________________________________
722
723 Short_t AliCFUnfolding::Smooth() {
724   //
725   // Smoothes the unfolded spectrum
726   //
727   // By default each cell content is replaced by the average with the neighbouring bins (but not diagonally-neighbouring bins)
728   // However, if a specific function fcn has been defined in UseSmoothing(fcn), the unfolded will be fit and updated using fcn 
729   //
730   
731   if (fSmoothFunction) {
732     AliDebug(2,Form("Smoothing spectrum with fit function %p",fSmoothFunction));
733     return SmoothUsingFunction();
734   }
735   else return SmoothUsingNeighbours(fUnfolded);
736 }
737
738 //______________________________________________________________
739
740 Short_t AliCFUnfolding::SmoothUsingNeighbours(THnSparse* hist) {
741   //
742   // Smoothes the unfolded spectrum using neighouring bins
743   //
744
745   Int_t  const nDimensions = hist->GetNdimensions() ;
746   Int_t* coordinates = new Int_t[nDimensions];
747
748   Int_t* numBins = new Int_t[nDimensions];
749   for (Int_t iVar=0; iVar<nDimensions; iVar++) numBins[iVar] = hist->GetAxis(iVar)->GetNbins();
750   
751   //need a copy because hist will be updated during the loop, and this creates problems
752   THnSparse* copy = (THnSparse*)hist->Clone();
753
754   for (Long_t iBin=0; iBin<copy->GetNbins(); iBin++) { //loop on non-empty bins
755     Double_t content = copy->GetBinContent(iBin,coordinates);
756     Double_t error2  = TMath::Power(copy->GetBinError(iBin),2);
757
758     // skip the under/overflow bins...
759     Bool_t isOutside = kFALSE ;
760     for (Int_t iVar=0; iVar<nDimensions; iVar++) {
761       if (coordinates[iVar]<1 || coordinates[iVar]>numBins[iVar]) {
762         isOutside=kTRUE;
763         break;
764       }
765     }
766     if (isOutside) continue;
767     
768     Int_t neighbours = 0; // number of neighbours to average with
769
770     for (Int_t iVar=0; iVar<nDimensions; iVar++) {
771       if (coordinates[iVar] > 1) { // must not be on low edge border
772         coordinates[iVar]-- ; //get lower neighbouring bin 
773         content += copy->GetBinContent(coordinates);
774         error2  += TMath::Power(copy->GetBinError(coordinates),2);
775         neighbours++;
776         coordinates[iVar]++ ; //back to initial coordinate
777       }
778       if (coordinates[iVar] < numBins[iVar]) { // must not be on up edge border
779         coordinates[iVar]++ ; //get upper neighbouring bin
780         content += copy->GetBinContent(coordinates);
781         error2  += TMath::Power(copy->GetBinError(coordinates),2);
782         neighbours++;
783         coordinates[iVar]-- ; //back to initial coordinate
784       }
785     }
786     // make an average
787     hist->SetBinContent(coordinates,content/(1.+neighbours));
788     hist->SetBinError  (coordinates,TMath::Sqrt(error2)/(1.+neighbours));
789   }
790   delete [] numBins;
791   delete [] coordinates ;
792   delete copy;
793   return 0;
794 }
795
796 //______________________________________________________________
797
798 Short_t AliCFUnfolding::SmoothUsingFunction() {
799   //
800   // Fits the unfolded spectrum using the function fSmoothFunction
801   //
802
803   AliDebug(0,Form("Smooth function is a %s with option \"%s\" and has %d parameters : ",fSmoothFunction->ClassName(),fSmoothOption,fSmoothFunction->GetNpar()));
804
805   for (Int_t iPar=0; iPar<fSmoothFunction->GetNpar(); iPar++) AliDebug(0,Form("par[%d]=%e",iPar,fSmoothFunction->GetParameter(iPar)));
806
807   Int_t fitResult = 0;
808
809   switch (fNVariables) {
810   case 1 : fitResult = fUnfolded->Projection(0)    ->Fit(fSmoothFunction,fSmoothOption); break;
811   case 2 : fitResult = fUnfolded->Projection(1,0)  ->Fit(fSmoothFunction,fSmoothOption); break; // (1,0) instead of (0,1) -> TAxis issue
812   case 3 : fitResult = fUnfolded->Projection(0,1,2)->Fit(fSmoothFunction,fSmoothOption); break;
813   default: AliFatal(Form("Cannot handle such fit in %d dimensions",fNVariables)) ; return 1;
814   }
815
816   if (fitResult != 0) {
817     AliWarning(Form("Fit failed with status %d, stopping the loop",fitResult));
818     return 1;
819   }
820
821   Int_t nDim = fNVariables;
822   Int_t* bins = new Int_t[nDim]; // number of bins for each variable
823   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
824
825   for (Int_t iVar=0; iVar<nDim; iVar++) {
826     bins[iVar] = fUnfolded->GetAxis(iVar)->GetNbins();
827     nBins *= bins[iVar];
828   }
829
830   Int_t *bin  = new Int_t[nDim];    // bin to fill the THnSparse (holding the bin coordinates)
831   Double_t x[3] = {0,0,0} ;         // value in bin center (max dimension is 3 (TF3))
832
833   // loop on the bins and update of fUnfolded
834   // THnSparse::Multiply(TF1*) doesn't exist, so let's do it bin by bin
835   for (Long_t iBin=0; iBin<nBins; iBin++) {
836     Long_t bin_tmp = iBin ;
837     for (Int_t iVar=0; iVar<nDim; iVar++) {
838       bin[iVar] = 1 + bin_tmp % bins[iVar] ;
839       bin_tmp /= bins[iVar] ;
840       x[iVar] = fUnfolded->GetAxis(iVar)->GetBinCenter(bin[iVar]);
841     }
842     Double_t functionValue = fSmoothFunction->Eval(x[0],x[1],x[2]) ;
843     fUnfolded->SetBinError  (bin,fUnfolded->GetBinError(bin)*functionValue/fUnfolded->GetBinContent(bin));
844     fUnfolded->SetBinContent(bin,functionValue);
845   }
846   delete [] bins;
847   delete [] bin ;
848   return 0;
849 }
850
851 //______________________________________________________________
852
853 void AliCFUnfolding::CreateFlatPrior() {
854   //
855   // Creates a flat prior distribution
856   // 
857
858   AliInfo("Creating a flat a priori distribution");
859   
860   // create the frame of the THnSparse given (for example) the one from the efficiency map
861   fPrior = (THnSparse*) fEfficiency->Clone();
862
863   if (fNVariables != fPrior->GetNdimensions()) 
864     AliFatal(Form("The prior matrix should have %d dimensions, and it has actually %d",fNVariables,fPrior->GetNdimensions()));
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] = fPrior->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
877   // loop that sets 1 in each bin
878   for (Long_t iBin=0; iBin<nBins; iBin++) {
879     Long_t bin_tmp = iBin ;
880     for (Int_t iVar=0; iVar<nDim; iVar++) {
881       bin[iVar] = 1 + bin_tmp % bins[iVar] ;
882       bin_tmp /= bins[iVar] ;
883     }
884     fPrior->SetBinContent(bin,1.); // put 1 everywhere
885     fPrior->SetBinError  (bin,0.); // put 0 everywhere
886   }
887   
888   fOriginalPrior = (THnSparse*)fPrior->Clone();
889
890   delete [] bin;
891   delete [] bins;
892 }