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