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