adding calculation covariance matrix to bayesian method
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AliMultiplicityCorrection.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 /* $Id$ */
17
18 // This class is used to store correction maps, raw input and results of the multiplicity
19 // measurement with the ITS or TPC
20 // It also contains functions to correct the spectrum using different methods.
21 //
22 //  Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch
23
24 #include "AliMultiplicityCorrection.h"
25
26 #include <TFile.h>
27 #include <TH1F.h>
28 #include <TH2F.h>
29 #include <TH3F.h>
30 #include <TDirectory.h>
31 #include <TVirtualFitter.h>
32 #include <TCanvas.h>
33 #include <TString.h>
34 #include <TF1.h>
35 #include <TMath.h>
36 #include <TCollection.h>
37
38 ClassImp(AliMultiplicityCorrection)
39
40 const Int_t AliMultiplicityCorrection::fgMaxParams = 200;
41 TH1* AliMultiplicityCorrection::fCurrentESD = 0;
42 TH1* AliMultiplicityCorrection::fCurrentCorrelation = 0;
43 TH1* AliMultiplicityCorrection::fCurrentEfficiency = 0;
44 TMatrixF* AliMultiplicityCorrection::fCorrelationMatrix = 0;
45 TMatrixF* AliMultiplicityCorrection::fCorrelationCovarianceMatrix = 0;
46 TVectorF* AliMultiplicityCorrection::fCurrentESDVector = 0;
47 AliMultiplicityCorrection::RegularizationType AliMultiplicityCorrection::fRegularizationType = AliMultiplicityCorrection::kPol1;
48 Float_t AliMultiplicityCorrection::fRegularizationWeight = 1e4;
49
50 //____________________________________________________________________
51 AliMultiplicityCorrection::AliMultiplicityCorrection() :
52   TNamed(), fLastChi2MC(0), fLastChi2Residuals(0)
53 {
54   //
55   // default constructor
56   //
57
58   for (Int_t i = 0; i < kESDHists; ++i)
59     fMultiplicityESD[i] = 0;
60
61   for (Int_t i = 0; i < kMCHists; ++i)
62   {
63     fMultiplicityVtx[i] = 0;
64     fMultiplicityMB[i] = 0;
65     fMultiplicityINEL[i] = 0;
66   }
67
68   for (Int_t i = 0; i < kCorrHists; ++i)
69   {
70     fCorrelation[i] = 0;
71     fMultiplicityESDCorrected[i] = 0;
72   }
73 }
74
75 //____________________________________________________________________
76 AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const Char_t* title) :
77   TNamed(name, title), fLastChi2MC(0), fLastChi2Residuals(0)
78 {
79   //
80   // named constructor
81   //
82
83   // do not add this hists to the directory
84   Bool_t oldStatus = TH1::AddDirectoryStatus();
85   TH1::AddDirectory(kFALSE);
86
87   /*Float_t binLimitsVtx[] = {-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10};
88   Float_t binLimitsN[] = {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
89                           10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
90                           20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
91                           30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
92                           40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5,
93                           50.5, 55.5, 60.5, 65.5, 70.5, 75.5, 80.5, 85.5, 90.5, 95.5,
94                           100.5, 105.5, 110.5, 115.5, 120.5, 125.5, 130.5, 135.5, 140.5, 145.5,
95                           150.5, 160.5, 170.5, 180.5, 190.5, 200.5, 210.5, 220.5, 230.5, 240.5,
96                           250.5, 275.5, 300.5, 325.5, 350.5, 375.5, 400.5, 425.5, 450.5, 475.5,
97                           500.5 };
98                                 //525.5, 550.5, 575.5, 600.5, 625.5, 650.5, 675.5, 700.5, 725.5,
99                           //750.5, 775.5, 800.5, 825.5, 850.5, 875.5, 900.5, 925.5, 950.5, 975.5,
100                           //1000.5 };
101
102   #define VTXBINNING 10, binLimitsVtx
103   #define NBINNING fgMaxParams, binLimitsN*/
104
105   #define NBINNING 251, -0.5, 250.5
106   #define VTXBINNING 10, -10, 10
107
108   for (Int_t i = 0; i < kESDHists; ++i)
109     fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", VTXBINNING, NBINNING);
110
111   for (Int_t i = 0; i < kMCHists; ++i)
112   {
113     fMultiplicityVtx[i] = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("fMultiplicityVtx%d", i)));
114     fMultiplicityVtx[i]->SetTitle("fMultiplicityVtx;vtx-z;Npart");
115
116     fMultiplicityMB[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityMB%d", i)));
117     fMultiplicityMB[i]->SetTitle("fMultiplicityMB");
118
119     fMultiplicityINEL[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityINEL%d", i)));
120     fMultiplicityINEL[i]->SetTitle("fMultiplicityINEL");
121   }
122
123   for (Int_t i = 0; i < kCorrHists; ++i)
124   {
125     fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;Npart;Ntracks", VTXBINNING, NBINNING, NBINNING);
126     fMultiplicityESDCorrected[i] = new TH1F(Form("fMultiplicityESDCorrected%d", i), "fMultiplicityESDCorrected;Npart;dN/dN", NBINNING);
127   }
128
129   TH1::AddDirectory(oldStatus);
130 }
131
132 //____________________________________________________________________
133 AliMultiplicityCorrection::~AliMultiplicityCorrection()
134 {
135   //
136   // Destructor
137   //
138 }
139
140 //____________________________________________________________________
141 Long64_t AliMultiplicityCorrection::Merge(TCollection* list)
142 {
143   // Merge a list of AliMultiplicityCorrection objects with this (needed for
144   // PROOF).
145   // Returns the number of merged objects (including this).
146
147   if (!list)
148     return 0;
149
150   if (list->IsEmpty())
151     return 1;
152
153   TIterator* iter = list->MakeIterator();
154   TObject* obj;
155
156   // collections of all histograms
157   TList collections[kESDHists+kMCHists*3+kCorrHists*2];
158
159   Int_t count = 0;
160   while ((obj = iter->Next())) {
161
162     AliMultiplicityCorrection* entry = dynamic_cast<AliMultiplicityCorrection*> (obj);
163     if (entry == 0) 
164       continue;
165
166     for (Int_t i = 0; i < kESDHists; ++i)
167       collections[i].Add(entry->fMultiplicityESD[i]);
168
169     for (Int_t i = 0; i < kMCHists; ++i)
170     {
171       collections[kESDHists+i].Add(entry->fMultiplicityVtx[i]);
172       collections[kESDHists+kMCHists+i].Add(entry->fMultiplicityMB[i]);
173       collections[kESDHists+kMCHists*2+i].Add(entry->fMultiplicityINEL[i]);
174     }
175
176     for (Int_t i = 0; i < kCorrHists; ++i)
177       collections[kESDHists+kMCHists*3+i].Add(entry->fCorrelation[i]);
178
179     for (Int_t i = 0; i < kCorrHists; ++i)
180       collections[kESDHists+kMCHists*3+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]);
181
182     count++;
183   }
184
185   for (Int_t i = 0; i < kESDHists; ++i)
186     fMultiplicityESD[i]->Merge(&collections[i]);
187
188   for (Int_t i = 0; i < kMCHists; ++i)
189   {
190     fMultiplicityVtx[i]->Merge(&collections[kESDHists+i]);
191     fMultiplicityMB[i]->Merge(&collections[kESDHists+kMCHists+i]);
192     fMultiplicityINEL[i]->Merge(&collections[kESDHists+kMCHists*2+i]);
193   }
194
195   for (Int_t i = 0; i < kCorrHists; ++i)
196     fCorrelation[i]->Merge(&collections[kESDHists+kMCHists*3+i]);
197
198   for (Int_t i = 0; i < kCorrHists; ++i)
199     fMultiplicityESDCorrected[i]->Merge(&collections[kESDHists+kMCHists*3+kCorrHists+i]);
200
201   delete iter;
202
203   return count+1;
204 }
205
206 //____________________________________________________________________
207 Bool_t AliMultiplicityCorrection::LoadHistograms(const Char_t* dir)
208 {
209   //
210   // loads the histograms from a file
211   // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
212   //
213
214   if (!dir)
215     dir = GetName();
216
217   if (!gDirectory->cd(dir))
218     return kFALSE;
219
220   // TODO memory leak. old histograms needs to be deleted.
221
222   Bool_t success = kTRUE;
223
224   for (Int_t i = 0; i < kESDHists; ++i)
225   {
226     fMultiplicityESD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityESD[i]->GetName()));
227     if (!fMultiplicityESD[i])
228       success = kFALSE;
229   }
230
231   for (Int_t i = 0; i < kMCHists; ++i)
232   {
233     fMultiplicityVtx[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityVtx[i]->GetName()));
234     fMultiplicityMB[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityMB[i]->GetName()));
235     fMultiplicityINEL[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityINEL[i]->GetName()));
236     if (!fMultiplicityVtx[i] || !fMultiplicityMB[i] || !fMultiplicityINEL[i])
237       success = kFALSE;
238   }
239
240   for (Int_t i = 0; i < kCorrHists; ++i)
241   {
242     fCorrelation[i] = dynamic_cast<TH3F*> (gDirectory->Get(fCorrelation[i]->GetName()));
243     if (!fCorrelation[i])
244       success = kFALSE;
245     fMultiplicityESDCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fMultiplicityESDCorrected[i]->GetName()));
246     if (!fMultiplicityESDCorrected[i])
247       success = kFALSE;
248   }
249
250   gDirectory->cd("..");
251
252   return success;
253 }
254
255 //____________________________________________________________________
256 void AliMultiplicityCorrection::SaveHistograms()
257 {
258   //
259   // saves the histograms
260   //
261
262   gDirectory->mkdir(GetName());
263   gDirectory->cd(GetName());
264
265   for (Int_t i = 0; i < kESDHists; ++i)
266     if (fMultiplicityESD[i])
267       fMultiplicityESD[i]->Write();
268
269   for (Int_t i = 0; i < kMCHists; ++i)
270   {
271     if (fMultiplicityVtx[i])
272       fMultiplicityVtx[i]->Write();
273     if (fMultiplicityMB[i])
274       fMultiplicityMB[i]->Write();
275     if (fMultiplicityINEL[i])
276       fMultiplicityINEL[i]->Write();
277   }
278
279   for (Int_t i = 0; i < kCorrHists; ++i)
280   {
281     if (fCorrelation[i])
282       fCorrelation[i]->Write();
283     if (fMultiplicityESDCorrected[i])
284       fMultiplicityESDCorrected[i]->Write();
285   }
286
287   gDirectory->cd("..");
288 }
289
290 //____________________________________________________________________
291 void AliMultiplicityCorrection::FillGenerated(Float_t vtx, Bool_t triggered, Bool_t vertex, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll)
292 {
293   //
294   // Fills an event from MC
295   //
296
297   if (triggered)
298   {
299     fMultiplicityMB[0]->Fill(vtx, generated05);
300     fMultiplicityMB[1]->Fill(vtx, generated10);
301     fMultiplicityMB[2]->Fill(vtx, generated15);
302     fMultiplicityMB[3]->Fill(vtx, generated20);
303     fMultiplicityMB[4]->Fill(vtx, generatedAll);
304
305     if (vertex)
306     {
307       fMultiplicityVtx[0]->Fill(vtx, generated05);
308       fMultiplicityVtx[1]->Fill(vtx, generated10);
309       fMultiplicityVtx[2]->Fill(vtx, generated15);
310       fMultiplicityVtx[3]->Fill(vtx, generated20);
311       fMultiplicityVtx[4]->Fill(vtx, generatedAll);
312     }
313   }
314
315   fMultiplicityINEL[0]->Fill(vtx, generated05);
316   fMultiplicityINEL[1]->Fill(vtx, generated10);
317   fMultiplicityINEL[2]->Fill(vtx, generated15);
318   fMultiplicityINEL[3]->Fill(vtx, generated20);
319   fMultiplicityINEL[4]->Fill(vtx, generatedAll);
320 }
321
322 //____________________________________________________________________
323 void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20)
324 {
325   //
326   // Fills an event from ESD
327   //
328
329   fMultiplicityESD[0]->Fill(vtx, measured05);
330   fMultiplicityESD[1]->Fill(vtx, measured10);
331   fMultiplicityESD[2]->Fill(vtx, measured15);
332   fMultiplicityESD[3]->Fill(vtx, measured20);
333 }
334
335 //____________________________________________________________________
336 void AliMultiplicityCorrection::FillCorrection(Float_t vtx, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20)
337 {
338   //
339   // Fills an event into the correlation map with the information from MC and ESD
340   //
341
342   fCorrelation[0]->Fill(vtx, generated05, measured05);
343   fCorrelation[1]->Fill(vtx, generated10, measured10);
344   fCorrelation[2]->Fill(vtx, generated15, measured15);
345   fCorrelation[3]->Fill(vtx, generated20, measured20);
346
347   fCorrelation[4]->Fill(vtx, generatedAll, measured05);
348   fCorrelation[5]->Fill(vtx, generatedAll, measured10);
349   fCorrelation[6]->Fill(vtx, generatedAll, measured15);
350   fCorrelation[7]->Fill(vtx, generatedAll, measured20);
351 }
352
353 //____________________________________________________________________
354 Double_t AliMultiplicityCorrection::RegularizationPol0(Double_t *params)
355 {
356   // homogenity term for minuit fitting
357   // pure function of the parameters
358   // prefers constant function (pol0)
359
360   Double_t chi2 = 0;
361
362   for (Int_t i=1; i<fgMaxParams; ++i)
363   {
364     if (params[i] == 0)
365       continue;
366
367     Double_t right  = params[i]; // / fCurrentESD->GetBinWidth(i+1);
368     Double_t left   = params[i-1]; // / fCurrentESD->GetBinWidth(i);
369
370     Double_t diff = (right - left) / right;
371     chi2 += diff * diff;
372   }
373
374   return chi2;
375 }
376
377 //____________________________________________________________________
378 Double_t AliMultiplicityCorrection::RegularizationPol1(Double_t *params)
379 {
380   // homogenity term for minuit fitting
381   // pure function of the parameters
382   // prefers linear function (pol1)
383
384   Double_t chi2 = 0;
385
386   for (Int_t i=2; i<fgMaxParams; ++i)
387   {
388     if (params[i] == 0 || params[i-1] == 0)
389       continue;
390
391     Double_t right  = params[i]; // / fCurrentESD->GetBinWidth(i+1);
392     Double_t middle = params[i-1]; // / fCurrentESD->GetBinWidth(i);
393     Double_t left   = params[i-2]; // / fCurrentESD->GetBinWidth(i-1);
394
395     Double_t der1 = (right - middle); // / fCurrentESD->GetBinWidth(i);
396     Double_t der2 = (middle - left); //  / fCurrentESD->GetBinWidth(i-1);
397
398     Double_t diff = (der1 - der2) / middle; /// fCurrentESD->GetBinWidth(i);
399
400     // TODO give additonal weight to big bins
401     chi2 += diff * diff; // * fCurrentESD->GetBinWidth(i) * fCurrentESD->GetBinWidth(i-1);
402
403     //printf("%d --> %f\n", i, diff);
404   }
405
406   return chi2; // 10000
407 }
408
409 //____________________________________________________________________
410 Double_t AliMultiplicityCorrection::RegularizationTotalCurvature(Double_t *params)
411 {
412   // homogenity term for minuit fitting
413   // pure function of the parameters
414   // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments,
415   // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp.
416
417   Double_t chi2 = 0;
418
419   for (Int_t i=2; i<fgMaxParams; ++i)
420   {
421     if (params[i] == 0 || params[i-1] == 0)
422       continue;
423
424     Double_t right  = params[i]; // / fCurrentESD->GetBinWidth(i+1);
425     Double_t middle = params[i-1]; // / fCurrentESD->GetBinWidth(i);
426     Double_t left   = params[i-2]; // / fCurrentESD->GetBinWidth(i-1);
427
428     Double_t der1 = (right - middle); // / fCurrentESD->GetBinWidth(i);
429     Double_t der2 = (middle - left); //  / fCurrentESD->GetBinWidth(i-1);
430
431     Double_t secDer = (der1 - der2); // / fCurrentESD->GetBinWidth(i);
432
433     // square and weight with the bin width
434     chi2 += secDer * secDer; // * fCurrentESD->GetBinWidth(i) * fCurrentESD->GetBinWidth(i-1);
435
436     //printf("%d --> %f\n", i, secDer);
437   }
438
439   return chi2; // 10
440 }
441
442 //____________________________________________________________________
443 Double_t AliMultiplicityCorrection::RegularizationEntropy(Double_t *params)
444 {
445   // homogenity term for minuit fitting
446   // pure function of the parameters
447   // calculates entropy, from
448   // The method of reduced cross-entropy (M. Schmelling 1993)
449
450   //static Int_t callCount = 0;
451
452   Double_t paramSum = 0;
453   for (Int_t i=0; i<fgMaxParams; ++i)
454     paramSum += params[i]; // / fCurrentESD->GetBinWidth(i+1);
455
456   Double_t chi2 = 0;
457   for (Int_t i=0; i<fgMaxParams; ++i)
458   {
459     Double_t tmp = params[i] / paramSum; // / fCurrentESD->GetBinWidth(i+1);
460     if (tmp > 0)
461     {
462       chi2 += tmp * log(tmp); /* * fCurrentESD->GetBinWidth(i+1);
463       //                                   /\
464       //                                    ------------------------------------
465       // TODO WARNING the absolute fitting values seem to depend on that value |
466       // this should be impossible...
467       //if ((callCount % 100) == 0)
468       //  printf("%f => %f\n", params[i], tmp * log(tmp)); */
469     }
470   }
471   //if ((callCount++ % 100) == 0)
472   //  printf("\n");
473
474   return chi2; // 1000
475 }
476
477 //____________________________________________________________________
478 void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t)
479 {
480   //
481   // fit function for minuit
482   // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix
483   //
484
485   static Int_t callCount = 0;
486
487   Double_t chi2FromFit = 0;
488
489   // d
490   TVectorF paramsVector(fgMaxParams);
491   for (Int_t i=0; i<fgMaxParams; ++i)
492     paramsVector[i] = params[i];
493
494   // Ad
495   paramsVector = (*fCorrelationMatrix) * paramsVector;
496
497   // Ad - m
498   paramsVector -= (*fCurrentESDVector);
499
500   TVectorF copy(paramsVector);
501
502   // (Ad - m) W
503   paramsVector *= (*fCorrelationCovarianceMatrix);
504
505   // (Ad - m) W (Ad - m)
506   chi2FromFit = paramsVector * copy;
507
508   /*printf("chi2FromFit = %f\n", chi2FromFit);
509   chi2FromFit = 0;
510
511   // loop over multiplicity (x axis of fMultiplicityESD)
512   for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i)
513   {
514     if (i > fCurrentCorrelation->GetNbinsY() || i > fgMaxParams)
515       break;
516
517     Double_t sum = 0;
518     //Double_t error = 0;
519     // loop over generated particles (x axis of fCorrelation) that resulted in reconstruction of i tracks
520     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsX(); ++j)
521     {
522       if (j > fgMaxParams)
523         break;
524
525       sum += fCurrentCorrelation->GetBinContent(j, i) * params[j-1];
526
527       //if (params[j-1] > 0)
528       //  error += fCurrentCorrelation->GetBinError(j, i) * fCurrentCorrelation->GetBinError(j, i) * params[j-1];
529     }
530
531     //printf("%f\n", sum);
532
533     Double_t diff = fCurrentESD->GetBinContent(i) - sum;
534
535     // include error
536     if (fCurrentESD->GetBinError(i) > 0)
537       diff /= fCurrentESD->GetBinError(i);
538
539     //if (fCurrentESD->GetBinContent(i) > 0)
540     //  diff /= fCurrentESD->GetBinContent(i);
541     //else
542     //  diff /= fCurrentESD->Integral();
543
544     // weight with bin width
545     //diff *= fCurrentESD->GetBinWidth(i);
546
547     //error = TMath::Sqrt(error) + fCurrentESD->GetBinError(i);
548     //if (error <= 0)
549     // error = 1;
550
551     //Double_t tmp = diff / error;
552     //chi2 += tmp * tmp;
553     chi2FromFit += diff * diff;
554
555     //printf("\nExpected sum = %f; Diff for bin %d is %f\n**********************************\n", fCurrentESD->GetBinContent(i), i, diff);
556     //printf("Diff for bin %d is %f\n", i, diff);
557   }
558
559   printf("chi2FromFit = %f\n", chi2FromFit);*/
560
561   Double_t penaltyVal = 0;
562
563   switch (fRegularizationType)
564   {
565     case kNone:      break;
566     case kPol0:      penaltyVal = RegularizationPol0(params); break;
567     case kPol1:      penaltyVal = RegularizationPol1(params); break;
568     case kCurvature: penaltyVal = RegularizationTotalCurvature(params); break;
569     case kEntropy:   penaltyVal = RegularizationEntropy(params); break;
570   }
571
572   penaltyVal *= fRegularizationWeight;
573
574   chi2 = chi2FromFit + penaltyVal;
575
576   if ((callCount++ % 1000) == 0)
577     printf("%f %f --> %f\n", chi2FromFit, penaltyVal, chi2);
578 }
579
580 //____________________________________________________________________
581 void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
582 {
583   //
584   // fills fCurrentESD, fCurrentCorrelation
585   // resets fMultiplicityESDCorrected
586   //
587
588   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
589   fMultiplicityESDCorrected[correlationID]->Reset();
590
591   fCurrentESD = fMultiplicityESD[inputRange]->ProjectionY();
592   fCurrentESD->Sumw2();
593
594   // empty under/overflow bins in x, otherwise Project3D takes them into account
595   TH3* hist = fCorrelation[correlationID];
596   for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
597   {
598     for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
599     {
600       hist->SetBinContent(0, y, z, 0);
601       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
602     }
603   }
604
605   fCurrentCorrelation = hist->Project3D("zy");
606   fCurrentCorrelation->Sumw2();
607
608   fCurrentEfficiency = fMultiplicityVtx[inputRange]->ProjectionY("CurrentEfficiency");
609   TH2* divisor = 0;
610   switch (eventType)
611   {
612     case kTrVtx : divisor = fMultiplicityVtx[inputRange]; break;
613     case kMB: divisor = fMultiplicityMB[inputRange]; break;
614     case kINEL: divisor = fMultiplicityINEL[inputRange]; break;
615   }
616   fCurrentEfficiency->Divide(divisor->ProjectionY());
617 }
618
619 //____________________________________________________________________
620 void AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, Bool_t check)
621 {
622   //
623   // correct spectrum using minuit chi2 method
624   //
625   // if check is kTRUE the input MC solution (by definition the right solution) is used
626   // no fit is made and just the chi2 is calculated
627   //
628
629   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
630   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
631
632   SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx);
633
634   fCorrelationMatrix = new TMatrixF(fgMaxParams, fgMaxParams);
635   fCorrelationCovarianceMatrix = new TMatrixF(fgMaxParams, fgMaxParams);
636   fCurrentESDVector = new TVectorF(fgMaxParams);
637
638   // normalize correction for given nPart
639   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
640   {
641     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
642     if (sum <= 0)
643       continue;
644     Float_t maxValue = 0;
645     Int_t maxBin = -1;
646     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
647     {
648       // find most probably value
649       if (maxValue < fCurrentCorrelation->GetBinContent(i, j))
650       {
651         maxValue = fCurrentCorrelation->GetBinContent(i, j);
652         maxBin = j;
653       }
654
655       // npart sum to 1
656       fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
657       fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
658
659       if (i <= fgMaxParams && j <= fgMaxParams)
660         (*fCorrelationMatrix)(j-1, i-1) = fCurrentCorrelation->GetBinContent(i, j);
661     }
662
663     //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
664   }
665
666   // small hack to get around charge conservation for full phase space ;-)
667   /*if (fullPhaseSpace)
668   {
669     for (Int_t i=2; i<=50; i+=2)
670       for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
671         fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i-1, j));
672   }*/
673
674   /*TCanvas* canvas1 = new TCanvas("ApplyMinuitFit", "ApplyMinuitFit", 800, 400);
675   canvas1->Divide(2, 1);
676   canvas1->cd(1);
677   fCurrentESD->DrawCopy();
678   canvas1->cd(2);
679   fCurrentCorrelation->DrawCopy("COLZ");*/
680
681   // Initialize TMinuit via generic fitter interface
682   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgMaxParams);
683
684   minuit->SetFCN(MinuitFitFunction);
685
686   static Double_t results[fgMaxParams];
687   //printf("%x\n", (void*) results);
688
689   TH1* proj = fMultiplicityVtx[mcTarget]->ProjectionY();
690   for (Int_t i=0; i<fgMaxParams; ++i)
691   {
692     //results[i] = 1;
693     results[i] = fCurrentESD->GetBinContent(i+1);
694     if (results[i] < 1e-2)
695       results[i] = 1e-2;
696     if (check)
697       results[i] = proj->GetBinContent(i+1);
698     minuit->SetParameter(i, Form("param%d", i), results[i], results[i] / 10, 0.001, fCurrentESD->GetMaximum() * 10);
699
700     (*fCurrentESDVector)[i] = fCurrentESD->GetBinContent(i+1);
701     if (fCurrentESD->GetBinError(i+1) > 0)
702       (*fCorrelationCovarianceMatrix)(i, i) = 1.0 / (fCurrentESD->GetBinError(i+1) * fCurrentESD->GetBinError(i+1));
703
704     //minuit->SetParameter(i, Form("param%d", i), fCurrentESD->GetBinWidth(i+1), 0.01, 0.001, 50);
705   }
706   minuit->SetParameter(0, "param0", results[1], ((results[1] > 1) ? (results[1] / 10) : 0), 0.001, fCurrentESD->GetMaximum() * 10);
707   //minuit->SetParameter(0, "param0", 1, 0, 1, 1);
708
709   Int_t dummy;
710   Double_t chi2;
711   MinuitFitFunction(dummy, 0, chi2, results, 0);
712   printf("Chi2 of right solution is = %f\n", chi2);
713
714   if (check)
715     return;
716
717   Double_t arglist[100];
718   arglist[0] = 0;
719   minuit->ExecuteCommand("SET PRINT", arglist, 1);
720   minuit->ExecuteCommand("MIGRAD", arglist, 1);
721   //minuit->ExecuteCommand("MINOS", arglist, 0);
722
723   for (Int_t i=0; i<fgMaxParams; ++i)
724   {
725     results[i] = minuit->GetParameter(i);
726     fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, results[i]);
727     fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0);
728   }
729
730   //printf("Penalty is %f\n", RegularizationTotalCurvature(results));
731 }
732
733 //____________________________________________________________________
734 void AliMultiplicityCorrection::NormalizeToBinWidth(TH1* hist)
735 {
736   //
737   // normalizes a 1-d histogram to its bin width
738   //
739
740   for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
741   {
742     hist->SetBinContent(i, hist->GetBinContent(i) / hist->GetBinWidth(i));
743     hist->SetBinError(i, hist->GetBinError(i) / hist->GetBinWidth(i));
744   }
745 }
746
747 //____________________________________________________________________
748 void AliMultiplicityCorrection::NormalizeToBinWidth(TH2* hist)
749 {
750   //
751   // normalizes a 2-d histogram to its bin width (x width * y width)
752   //
753
754   for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
755     for (Int_t j=1; j<=hist->GetNbinsY(); ++j)
756     {
757       Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j);
758       hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor);
759       hist->SetBinError(i, j, hist->GetBinError(i, j) / factor);
760     }
761 }
762
763 //____________________________________________________________________
764 void AliMultiplicityCorrection::ApplyMinuitFitAll()
765 {
766   //
767   // fit all eta ranges
768   //
769
770   for (Int_t i=0; i<kESDHists; ++i)
771   {
772     ApplyMinuitFit(i, kFALSE);
773     ApplyMinuitFit(i, kTRUE);
774   }
775 }
776
777 //____________________________________________________________________
778 void AliMultiplicityCorrection::DrawHistograms()
779 {
780   printf("ESD:\n");
781
782   TCanvas* canvas1 = new TCanvas("fMultiplicityESD", "fMultiplicityESD", 900, 600);
783   canvas1->Divide(3, 2);
784   for (Int_t i = 0; i < kESDHists; ++i)
785   {
786     canvas1->cd(i+1);
787     fMultiplicityESD[i]->DrawCopy("COLZ");
788     printf("%d --> %f\n", i, (Float_t) fMultiplicityESD[i]->ProjectionY()->GetMean());
789   }
790
791   printf("Vtx:\n");
792
793   TCanvas* canvas2 = new TCanvas("fMultiplicityMC", "fMultiplicityMC", 900, 600);
794   canvas2->Divide(3, 2);
795   for (Int_t i = 0; i < kMCHists; ++i)
796   {
797     canvas2->cd(i+1);
798     fMultiplicityVtx[i]->DrawCopy("COLZ");
799     printf("%d --> %f\n", i, fMultiplicityVtx[i]->ProjectionY()->GetMean());
800   }
801
802   TCanvas* canvas3 = new TCanvas("fCorrelation", "fCorrelation", 900, 900);
803   canvas3->Divide(3, 3);
804   for (Int_t i = 0; i < kCorrHists; ++i)
805   {
806     canvas3->cd(i+1);
807     TH3* hist = dynamic_cast<TH3*> (fCorrelation[i]->Clone());
808     for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
809     {
810       for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
811       {
812         hist->SetBinContent(0, y, z, 0);
813         hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
814       }
815     }
816     TH1* proj = hist->Project3D("zy");
817     proj->DrawCopy("COLZ");
818   }
819 }
820
821 //____________________________________________________________________
822 void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t normalizeESD, TH1* mcHist)
823 {
824   Int_t esdCorrId = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
825
826   TString tmpStr;
827   tmpStr.Form("%s_DrawComparison_%d", name, esdCorrId);
828
829   // calculate residual
830
831   // for that we convolute the response matrix with the gathered result
832   // if normalizeESD is not set, the histogram is already normalized, this needs to be passed to CalculateMultiplicityESD
833   TH2* convoluted = CalculateMultiplicityESD(fMultiplicityESDCorrected[esdCorrId], esdCorrId, !normalizeESD);
834   TH1* proj2 = convoluted->ProjectionY("proj2", -1, -1, "e");
835   NormalizeToBinWidth(proj2);
836   TH1* residual = convoluted->ProjectionY("residual", -1, -1, "e");
837   residual->SetTitle("Residuals");
838
839   residual->Add(fCurrentESD, -1);
840   residual->Divide(residual, fCurrentESD, 1, 1, "B");
841
842   // TODO fix errors
843   for (Int_t i=1; i<residual->GetNbinsX(); ++i)
844   {
845     proj2->SetBinError(i, 0);
846     residual->SetBinError(i, 0);
847   }
848
849   TCanvas* canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 800);
850   canvas1->Divide(2, 2);
851
852   canvas1->cd(1);
853   TH1* proj = (TH1*) mcHist->Clone("proj");
854   NormalizeToBinWidth(proj);
855
856   if (normalizeESD)
857     NormalizeToBinWidth(fMultiplicityESDCorrected[esdCorrId]);
858
859   proj->GetXaxis()->SetRangeUser(0, 150);
860   proj->DrawCopy("HIST");
861   gPad->SetLogy();
862
863   //fMultiplicityESDCorrected[esdCorrId]->SetMarkerStyle(3);
864   fMultiplicityESDCorrected[esdCorrId]->SetLineColor(2);
865   fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAME HIST");
866
867   Float_t chi2 = 0;
868   for (Int_t i=1; i<=100; ++i)
869   {
870     if (fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i) != 0)
871     {
872       Float_t value = (proj->GetBinContent(i) - fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i)) / fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i);
873       chi2 += value * value;
874     }
875   }
876
877   printf("Difference (from MC) is %f for bin 1-100\n", chi2);
878   fLastChi2MC = chi2;
879
880   canvas1->cd(2);
881
882   NormalizeToBinWidth(fCurrentESD);
883   gPad->SetLogy();
884   fCurrentESD->GetXaxis()->SetRangeUser(0, 150);
885   //fCurrentESD->SetLineColor(2);
886   fCurrentESD->DrawCopy("HIST");
887
888   proj2->SetLineColor(2);
889   //proj2->SetMarkerColor(2);
890   //proj2->SetMarkerStyle(5);
891   proj2->DrawCopy("HIST SAME");
892
893   chi2 = 0;
894   for (Int_t i=1; i<=100; ++i)
895   {
896     if (fCurrentESD->GetBinContent(i) != 0)
897     {
898       Float_t value = (proj2->GetBinContent(i) - fCurrentESD->GetBinContent(i)) / fCurrentESD->GetBinContent(i);
899       chi2 += value * value;
900     }
901   }
902
903   printf("Difference (Residuals) is %f for bin 1-100\n", chi2);
904   fLastChi2Residuals = chi2;
905
906   canvas1->cd(3);
907   TH1* clone = dynamic_cast<TH1*> (proj->Clone("clone"));
908   clone->Divide(fMultiplicityESDCorrected[esdCorrId]);
909   clone->SetTitle("Ratio;Npart;MC/ESD");
910   clone->GetYaxis()->SetRangeUser(0.8, 1.2);
911   clone->GetXaxis()->SetRangeUser(0, 150);
912   clone->DrawCopy("HIST");
913
914   /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
915   legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
916   legend->AddEntry(fMultiplicityMC, "MC");
917   legend->AddEntry(fMultiplicityESD, "ESD");
918   legend->Draw();*/
919
920   canvas1->cd(4);
921
922   residual->GetYaxis()->SetRangeUser(-0.2, 0.2);
923   residual->GetXaxis()->SetRangeUser(0, 150);
924   residual->DrawCopy();
925
926   canvas1->SaveAs(Form("%s.gif", canvas1->GetName()));
927 }
928
929 //____________________________________________________________________
930 void AliMultiplicityCorrection::GetComparisonResults(Float_t& mc, Float_t& residuals)
931 {
932   // Returns the chi2 between the MC and the unfolded ESD as well as between the ESD and the folded unfolded ESD
933   // These values are computed during DrawComparison, thus this function picks up the
934   // last calculation
935
936   mc = fLastChi2MC;
937   residuals = fLastChi2Residuals;
938 }
939
940
941 //____________________________________________________________________
942 TH2F* AliMultiplicityCorrection::GetMultiplicityMC(Int_t i, EventType eventType)
943 {
944   //
945   // returns the corresponding MC spectrum
946   //
947
948   switch (eventType)
949   {
950     case kTrVtx : return fMultiplicityVtx[i]; break;
951     case kMB : return fMultiplicityMB[i]; break;
952     case kINEL : return fMultiplicityINEL[i]; break;
953   }
954
955   return 0;
956 }
957
958 //____________________________________________________________________
959 void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar)
960 {
961   //
962   // correct spectrum using bayesian method
963   //
964
965   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
966
967   SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
968
969   // TODO should be taken from correlation map
970   TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
971
972   //new TCanvas;
973   //sumHist->Draw();
974
975   // normalize correction for given nPart
976   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
977   {
978     //Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
979     Double_t sum = sumHist->GetBinContent(i);
980     if (sum <= 0)
981       continue;
982     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
983     {
984       // npart sum to 1
985       fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
986       fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
987     }
988   }
989
990   //new TCanvas;
991   //fCurrentCorrelation->Draw("COLZ");
992   //return;
993
994   // normalize correction for given nTrack
995   /*for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
996   {
997     Double_t sum = fCurrentCorrelation->Integral(1, fCurrentCorrelation->GetNbinsX(), j, j);
998     if (sum <= 0)
999       continue;
1000     for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsY(); ++i)
1001     {
1002       // ntrack sum to 1
1003       fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1004       fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1005     }
1006   }*/
1007
1008   // FAKE
1009   fCurrentEfficiency = ((TH2*) fCurrentCorrelation)->ProjectionX("eff");
1010   //new TCanvas;
1011   //fCurrentEfficiency->Draw();
1012
1013   // pick prior distribution
1014   TH1F* hPrior = (TH1F*)fCurrentESD->Clone("prior");
1015   Float_t norm = hPrior->Integral();
1016   for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
1017     hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm);
1018
1019   // define temp hist
1020   TH1F* hTemp = (TH1F*)fCurrentESD->Clone("temp");
1021   hTemp->Reset();
1022
1023   // just a shortcut
1024   TH2F* hResponse = (TH2F*) fCurrentCorrelation;
1025
1026   // histogram to store the inverse response calculated from bayes theorem (from prior and response) IR
1027   TH2F* hInverseResponseBayes = (TH2F*)hResponse->Clone("pji");
1028   hInverseResponseBayes->Reset();
1029
1030   // unfold...
1031   Int_t iterations = 20;
1032   for (Int_t i=0; i<iterations; i++)
1033   {
1034     //printf(" iteration %i \n", i);
1035
1036     // calculate IR from Beyes theorem
1037     // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
1038
1039     for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
1040     {
1041       Float_t norm = 0;
1042       for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
1043         norm += hResponse->GetBinContent(t,m) * hPrior->GetBinContent(t);
1044       for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
1045       {
1046         if (norm==0)
1047           hInverseResponseBayes->SetBinContent(t,m,0);
1048         else
1049           hInverseResponseBayes->SetBinContent(t,m, hResponse->GetBinContent(t,m) * hPrior->GetBinContent(t)/norm);
1050       }
1051     }
1052
1053     for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
1054     {
1055       Float_t value = 0;
1056       for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
1057         value += fCurrentESD->GetBinContent(m) * hInverseResponseBayes->GetBinContent(t,m);
1058
1059       /*if (eventType == kTrVtx)
1060       {
1061         hTemp->SetBinContent(t, value);
1062       }
1063       else*/
1064       {
1065         if (fCurrentEfficiency->GetBinContent(t) > 0)
1066           hTemp->SetBinContent(t, value / fCurrentEfficiency->GetBinContent(t));
1067         else
1068           hTemp->SetBinContent(t, 0);
1069       }
1070     }
1071
1072     // this is the last guess, fill before (!) smoothing
1073     for (Int_t t=1; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++)
1074     {
1075       fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTemp->GetBinContent(t));
1076       fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t)); // TODO
1077
1078       //printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
1079     }
1080
1081     // regularization (simple smoothing)
1082     TH1F* hTrueSmoothed = (TH1F*) hTemp->Clone("truesmoothed");
1083
1084     for (Int_t t=2; t<hTrueSmoothed->GetNbinsX(); t++)
1085     {
1086       Float_t average = (hTemp->GetBinContent(t-1) / hTemp->GetBinWidth(t-1)
1087                          + hTemp->GetBinContent(t) / hTemp->GetBinWidth(t)
1088                          + hTemp->GetBinContent(t+1) / hTemp->GetBinWidth(t+1)) / 3.;
1089       average *= hTrueSmoothed->GetBinWidth(t);
1090
1091       // weight the average with the regularization parameter
1092       hTrueSmoothed->SetBinContent(t, (1-regPar) * hTemp->GetBinContent(t) + regPar * average);
1093     }
1094
1095     // calculate chi2 (change from last iteration)
1096     Double_t chi2 = 0;
1097
1098     // use smoothed true (normalized) as new prior
1099     Float_t norm = 1;
1100     //for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
1101     //  norm = norm + hTrueSmoothed->GetBinContent(t);
1102
1103     for (Int_t t=1; t<hTrueSmoothed->GetNbinsX(); t++)
1104     {
1105       Float_t newValue = hTrueSmoothed->GetBinContent(t)/norm;
1106       Float_t diff = hPrior->GetBinContent(t) - newValue;
1107       chi2 += (Double_t) diff * diff;
1108
1109       hPrior->SetBinContent(t, newValue);
1110     }
1111
1112     //printf("Chi2 of %d iteration = %.10f\n", i, chi2);
1113
1114     //if (i % 5 == 0)
1115     //  DrawComparison(Form("Bayesian_%d", i), mcTarget, correlationID, kTRUE, eventType);
1116
1117     delete hTrueSmoothed;
1118   } // end of iterations
1119
1120
1121   return;
1122
1123   // ********
1124   // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995
1125
1126   printf("Calculating covariance matrix. This may take some time...\n");
1127
1128   Int_t xBins = hInverseResponseBayes->GetNbinsX();
1129   Int_t yBins = hInverseResponseBayes->GetNbinsY();
1130
1131   // calculate "unfolding matrix" Mij
1132   Float_t matrixM[251][251];
1133   for (Int_t i=1; i<=xBins; i++)
1134   {
1135     for (Int_t j=1; j<=yBins; j++)
1136     {
1137       if (fCurrentEfficiency->GetBinContent(i) > 0)
1138         matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i);
1139       else
1140         matrixM[i-1][j-1] = 0;
1141     }
1142   }
1143
1144   Float_t* vectorn = new Float_t[yBins];
1145   for (Int_t j=1; j<=yBins; j++)
1146     vectorn[j-1] = fCurrentESD->GetBinContent(j);
1147
1148   // first part of covariance matrix, depends on input distribution n(E)
1149   Float_t cov1[251][251];
1150
1151   Float_t nEvents = fCurrentESD->Integral(); // N
1152
1153   xBins = 20;
1154   yBins = 20;
1155
1156   for (Int_t k=0; k<xBins; k++)
1157   {
1158     printf("In Cov1: %d\n", k);
1159     for (Int_t l=0; l<yBins; l++)
1160     {
1161       cov1[k][l] = 0;
1162
1163       // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N)
1164       for (Int_t j=0; j<yBins; j++)
1165         cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j]
1166           * (1.0 - vectorn[j] / nEvents);
1167
1168       // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N
1169       for (Int_t i=0; i<yBins; i++)
1170         for (Int_t j=0; j<yBins; j++)
1171         {
1172           if (i == j)
1173             continue;
1174           cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i]
1175             * vectorn[j] / nEvents;
1176          }
1177     }
1178   }
1179
1180   printf("Cov1 finished\n");
1181
1182   TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov");
1183   cov->Reset();
1184
1185   for (Int_t i=1; i<=xBins; i++)
1186     for (Int_t j=1; j<=yBins; j++)
1187       cov->SetBinContent(i, j, cov1[i-1][j-1]);
1188
1189   new TCanvas;
1190   cov->Draw("COLZ");
1191
1192   // second part of covariance matrix, depends on response matrix
1193   Float_t cov2[251][251];
1194
1195   // Cov[P(Er|Cu), P(Es|Cu)] term
1196   Float_t covTerm[100][100][100];
1197   for (Int_t r=0; r<yBins; r++)
1198     for (Int_t u=0; u<xBins; u++)
1199       for (Int_t s=0; s<yBins; s++)
1200       {
1201         if (r == s)
1202           covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
1203             * (1.0 - hResponse->GetBinContent(u+1, r+1));
1204         else
1205           covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
1206             * hResponse->GetBinContent(u+1, s+1);
1207       }
1208
1209   for (Int_t k=0; k<xBins; k++)
1210     for (Int_t l=0; l<yBins; l++)
1211     {
1212       cov2[k][l] = 0;
1213       printf("In Cov2: %d %d\n", k, l);
1214       for (Int_t i=0; i<yBins; i++)
1215         for (Int_t j=0; j<yBins; j++)
1216         {
1217           //printf("In Cov2: %d %d %d %d\n", k, l, i, j);
1218           // calculate Cov(Mki, Mlj) = sum{ru},{su} ...
1219           Float_t tmpCov = 0;
1220           for (Int_t r=0; r<yBins; r++)
1221             for (Int_t u=0; u<xBins; u++)
1222               for (Int_t s=0; s<yBins; s++)
1223               {
1224                 if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0
1225                   || hResponse->GetBinContent(u+1, i+1) == 0)
1226                   continue;
1227
1228                 tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u)
1229                   * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u)
1230                   * covTerm[r][u][s];
1231               }
1232
1233           cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov;
1234         }
1235     }
1236
1237   printf("Cov2 finished\n");
1238
1239   for (Int_t i=1; i<=xBins; i++)
1240     for (Int_t j=1; j<=yBins; j++)
1241       cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]);
1242
1243   new TCanvas;
1244   cov->Draw("COLZ");
1245 }
1246
1247 //____________________________________________________________________
1248 void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
1249 {
1250   //
1251   // correct spectrum using bayesian method
1252   //
1253
1254   Float_t regPar = 0.05;
1255
1256   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1257   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
1258
1259   SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
1260
1261   // TODO should be taken from correlation map
1262   //TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
1263
1264   // normalize correction for given nPart
1265   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1266   {
1267     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1268     //Double_t sum = sumHist->GetBinContent(i);
1269     if (sum <= 0)
1270       continue;
1271     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1272     {
1273       // npart sum to 1
1274       fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1275       fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1276     }
1277   }
1278
1279   new TCanvas;
1280   fCurrentCorrelation->Draw("COLZ");
1281
1282   // FAKE
1283   fCurrentEfficiency = ((TH2*) fCurrentCorrelation)->ProjectionX("eff");
1284
1285   // pick prior distribution
1286   TH1F* hPrior = (TH1F*)fCurrentESD->Clone("prior");
1287   Float_t norm = 1; //hPrior->Integral();
1288   for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
1289     hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm);
1290
1291   // zero distribution
1292   TH1F* zero =  (TH1F*)hPrior->Clone("zero");
1293
1294   // define temp hist
1295   TH1F* hTemp = (TH1F*)fCurrentESD->Clone("temp");
1296   hTemp->Reset();
1297
1298   // just a shortcut
1299   TH2F* hResponse = (TH2F*) fCurrentCorrelation;
1300
1301   // unfold...
1302   Int_t iterations = 20;
1303   for (Int_t i=0; i<iterations; i++)
1304   {
1305     //printf(" iteration %i \n", i);
1306
1307     for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
1308     {
1309       Float_t value = 0;
1310       for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
1311         value += hResponse->GetBinContent(t, m) * hPrior->GetBinContent(t);
1312       hTemp->SetBinContent(m, value);
1313       //printf("%d %f %f %f\n", m, zero->GetBinContent(m), hPrior->GetBinContent(m), value);
1314     }
1315
1316     // regularization (simple smoothing)
1317     TH1F* hTrueSmoothed = (TH1F*) hTemp->Clone("truesmoothed");
1318
1319     for (Int_t t=2; t<hTrueSmoothed->GetNbinsX(); t++)
1320     {
1321       Float_t average = (hTemp->GetBinContent(t-1) / hTemp->GetBinWidth(t-1)
1322                          + hTemp->GetBinContent(t) / hTemp->GetBinWidth(t)
1323                          + hTemp->GetBinContent(t+1) / hTemp->GetBinWidth(t+1)) / 3.;
1324       average *= hTrueSmoothed->GetBinWidth(t);
1325
1326       // weight the average with the regularization parameter
1327       hTrueSmoothed->SetBinContent(t, (1-regPar) * hTemp->GetBinContent(t) + regPar * average);
1328     }
1329
1330     for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
1331       hTemp->SetBinContent(m, zero->GetBinContent(m) + hPrior->GetBinContent(m) - hTrueSmoothed->GetBinContent(m));
1332
1333     // fill guess
1334     for (Int_t t=1; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++)
1335     {
1336       fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTemp->GetBinContent(t));
1337       fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t)); // TODO
1338
1339       //printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
1340     }
1341
1342
1343     // calculate chi2 (change from last iteration)
1344     Double_t chi2 = 0;
1345
1346     // use smoothed true (normalized) as new prior
1347     Float_t norm = 1; //hTrueSmoothed->Integral();
1348
1349     for (Int_t t=1; t<hTrueSmoothed->GetNbinsX(); t++)
1350     {
1351       Float_t newValue = hTrueSmoothed->GetBinContent(t)/norm;
1352       Float_t diff = hPrior->GetBinContent(t) - newValue;
1353       chi2 += (Double_t) diff * diff;
1354
1355       hPrior->SetBinContent(t, newValue);
1356     }
1357
1358     printf("Chi2 of %d iteration = %.10f\n", i, chi2);
1359
1360     if (i % 5 == 0)
1361       DrawComparison(Form("Laszlo_%d", i), inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
1362
1363     delete hTrueSmoothed;
1364   } // end of iterations
1365
1366   DrawComparison("Laszlo", inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
1367 }
1368
1369 //____________________________________________________________________
1370 Float_t AliMultiplicityCorrection::BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse,
1371   TH1* fCurrentEfficiency, Int_t k, Int_t i, Int_t r, Int_t u)
1372 {
1373   //
1374   //
1375   //
1376
1377   Float_t result = 0;
1378
1379   if (k == u && r == i)
1380     result += 1.0 / hResponse->GetBinContent(u+1, r+1);
1381
1382   if (k == u)
1383     result -= 1.0 / fCurrentEfficiency->GetBinContent(u+1);
1384
1385   if (r == i)
1386     result -= matrixM[u][i] * fCurrentEfficiency->GetBinContent(u+1) / hResponse->GetBinContent(u+1, i+1);
1387
1388   result *= matrixM[k][i];
1389
1390   return result;
1391 }
1392
1393 //____________________________________________________________________
1394 void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace)
1395 {
1396   //
1397   // correct spectrum using a simple Gaussian approach, that is model-dependent
1398   //
1399
1400   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1401   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
1402
1403   SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx);
1404
1405   NormalizeToBinWidth((TH2*) fCurrentCorrelation);
1406
1407   TH1D* correction = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianMean"));
1408   correction->SetTitle("GaussianMean");
1409
1410   TH1D* correctionWidth = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianWidth"));
1411   correction->SetTitle("GaussianWidth");
1412
1413   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1414   {
1415     TH1D* proj = (dynamic_cast<TH2*> (fCurrentCorrelation))->ProjectionX("_px", i, i+1);
1416     proj->Fit("gaus", "0Q");
1417     correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
1418     correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
1419
1420     continue;
1421
1422     // draw for debugging
1423     new TCanvas;
1424     proj->DrawCopy();
1425     proj->GetFunction("gaus")->DrawCopy("SAME");
1426   }
1427
1428   TH1* target = fMultiplicityESDCorrected[correlationID];
1429
1430   Int_t nBins = target->GetNbinsX()*10+1;
1431   Float_t* binning = new Float_t[nBins];
1432   for (Int_t i=1; i<=target->GetNbinsX(); ++i)
1433     for (Int_t j=0; j<10; ++j)
1434       binning[(i-1)*10 + j] = target->GetXaxis()->GetBinLowEdge(i) + target->GetXaxis()->GetBinWidth(i) / 10 * j;
1435
1436   binning[nBins-1] = target->GetXaxis()->GetBinUpEdge(target->GetNbinsX());
1437
1438   TH1F* fineBinned = new TH1F("targetFineBinned", "targetFineBinned", nBins-1, binning);
1439
1440   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1441   {
1442     Float_t mean = correction->GetBinContent(i);
1443     Float_t width = correctionWidth->GetBinContent(i);
1444
1445     Int_t fillBegin = fineBinned->FindBin(mean - width * 5);
1446     Int_t fillEnd   = fineBinned->FindBin(mean + width * 5);
1447     //printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
1448
1449     for (Int_t j=fillBegin; j <= fillEnd; ++j)
1450     {
1451       fineBinned->AddBinContent(j, TMath::Gaus(fineBinned->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fCurrentESD->GetBinContent(i));
1452     }
1453   }
1454
1455   for (Int_t i=1; i<=target->GetNbinsX(); ++i)
1456   {
1457     Float_t sum = 0;
1458     for (Int_t j=1; j<=10; ++j)
1459       sum += fineBinned->GetBinContent((i-1)*10 + j);
1460     target->SetBinContent(i, sum / 10);
1461   }
1462
1463   delete[] binning;
1464
1465   DrawComparison("Gaussian", inputRange, fullPhaseSpace, kFALSE, GetMultiplicityMC(mcTarget, kTrVtx)->ProjectionY());
1466 }
1467
1468 //____________________________________________________________________
1469 TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap, Bool_t normalized)
1470 {
1471   // runs the distribution given in inputMC through the response matrix identified by
1472   // correlationMap and produces a measured distribution
1473   // although it is a TH2F the vertex axis is not used at the moment and all entries are filled in mid-vertex
1474   // if normalized is set, inputMC is expected to be normalized to the bin width
1475
1476   if (!inputMC)
1477     return 0;
1478
1479   if (correlationMap < 0 || correlationMap >= kCorrHists)
1480     return 0;
1481
1482   // empty under/overflow bins in x, otherwise Project3D takes them into account
1483   TH3* hist = fCorrelation[correlationMap];
1484   for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
1485   {
1486     for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
1487     {
1488       hist->SetBinContent(0, y, z, 0);
1489       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
1490     }
1491   }
1492
1493   TH1* corr = hist->Project3D("zy");
1494   corr->Sumw2();
1495
1496   // normalize correction for given nPart
1497   for (Int_t i=1; i<=corr->GetNbinsX(); ++i)
1498   {
1499     Double_t sum = corr->Integral(i, i, 1, corr->GetNbinsY());
1500     if (sum <= 0)
1501       continue;
1502
1503     for (Int_t j=1; j<=corr->GetNbinsY(); ++j)
1504     {
1505       // npart sum to 1
1506       corr->SetBinContent(i, j, corr->GetBinContent(i, j) / sum);
1507       corr->SetBinError(i, j, corr->GetBinError(i, j) / sum);
1508     }
1509   }
1510
1511   TH2F* target = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("%s_measured", inputMC->GetName())));
1512   target->Reset();
1513
1514   for (Int_t meas=1; meas<=corr->GetNbinsY(); ++meas)
1515   {
1516     Float_t measured = 0;
1517     Float_t error = 0;
1518
1519     for (Int_t gen=1; gen<=target->GetNbinsY(); ++gen)
1520     {
1521       Int_t mcGenBin = inputMC->GetXaxis()->FindBin(target->GetYaxis()->GetBinCenter(gen));
1522
1523       Float_t factor = 1;
1524       if (normalized)
1525         factor = inputMC->GetBinWidth(mcGenBin);
1526
1527       measured += inputMC->GetBinContent(mcGenBin) * factor * corr->GetBinContent(gen, meas);
1528       error += inputMC->GetBinError(mcGenBin) * factor * corr->GetBinContent(gen, meas);
1529     }
1530
1531     // bin width of the <measured> axis has to be taken into account
1532     //measured /= target->GetYaxis()->GetBinWidth(meas);
1533     //error /= target->GetYaxis()->GetBinWidth(meas);
1534
1535     //printf("%f +- %f ; %f +- %f \n", inputMC->GetBinContent(meas), inputMC->GetBinError(meas), measured, error);
1536
1537     target->SetBinContent(target->GetNbinsX() / 2, meas, measured);
1538     target->SetBinError(target->GetNbinsX() / 2, meas, error);
1539   }
1540
1541   return target;
1542 }
1543
1544 //____________________________________________________________________
1545 void AliMultiplicityCorrection::SetGenMeasFromFunc(TF1* inputMC, Int_t id)
1546 {
1547   // uses the given function to fill the input MC histogram and generates from that
1548   // the measured histogram by applying the response matrix
1549   // this can be used to evaluate if the methods work indepedently of the input
1550   // distribution
1551   // WARNING does not respect the vertex distribution, just fills central vertex bin
1552
1553   if (!inputMC)
1554     return;
1555
1556   if (id < 0 || id >= kESDHists)
1557     return;
1558
1559   TH2F* mc = fMultiplicityINEL[id];
1560
1561   mc->Reset();
1562
1563   for (Int_t i=1; i<=mc->GetNbinsY(); ++i)
1564   {
1565     mc->SetBinContent(mc->GetNbinsX() / 2, i, inputMC->Eval(mc->GetYaxis()->GetBinCenter(i)) * mc->GetYaxis()->GetBinWidth(i));
1566     mc->SetBinError(mc->GetNbinsX() / 2, i, 0);
1567   }
1568
1569   //new TCanvas;
1570   //mc->Draw("COLZ");
1571
1572   TH1* proj = mc->ProjectionY();
1573
1574   fMultiplicityESD[id] = CalculateMultiplicityESD(proj, id);
1575 }