e25996e265faf1d10f9345ae3bc4bba6383d7ec0
[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 #include <TLegend.h>
38 #include <TLine.h>
39 #include <TRandom.h>
40 #include <TProfile.h>
41 #include <TProfile2D.h>
42
43 ClassImp(AliMultiplicityCorrection)
44
45 const Int_t AliMultiplicityCorrection::fgMaxInput  = 250;  // bins in measured histogram
46 const Int_t AliMultiplicityCorrection::fgMaxParams = 250;  // bins in unfolded histogram = number of fit params
47
48 TH1* AliMultiplicityCorrection::fCurrentESD = 0;
49 TH1* AliMultiplicityCorrection::fCurrentCorrelation = 0;
50 TH1* AliMultiplicityCorrection::fCurrentEfficiency = 0;
51 TMatrixD* AliMultiplicityCorrection::fCorrelationMatrix = 0;
52 TMatrixD* AliMultiplicityCorrection::fCorrelationCovarianceMatrix = 0;
53 TVectorD* AliMultiplicityCorrection::fCurrentESDVector = 0;
54 TVectorD* AliMultiplicityCorrection::fEntropyAPriori = 0;
55 AliMultiplicityCorrection::RegularizationType AliMultiplicityCorrection::fRegularizationType = AliMultiplicityCorrection::kPol1;
56 Float_t AliMultiplicityCorrection::fRegularizationWeight = 5000;
57 TF1* AliMultiplicityCorrection::fNBD = 0;
58
59 // These are the areas where the quality of the unfolding results are evaluated
60 // Default defined here, call SetQualityRegions to change them
61 // unit is in multiplicity (not in bin!)
62
63 // SPD:   peak area - flat area - low stat area
64 Int_t AliMultiplicityCorrection::fgQualityRegionsB[kQualityRegions] = {4,  60,  180};
65 Int_t AliMultiplicityCorrection::fgQualityRegionsE[kQualityRegions] = {20, 140, 210};
66
67 //____________________________________________________________________
68 void AliMultiplicityCorrection::SetQualityRegions(Bool_t SPDStudy)
69 {
70   //
71   // sets the quality region definition to TPC or SPD
72   //
73
74   if (SPDStudy)
75   {
76     // SPD:   peak area - flat area - low stat area
77     fgQualityRegionsB[0] = 4;
78     fgQualityRegionsE[0] = 20;
79
80     fgQualityRegionsB[1] = 60;
81     fgQualityRegionsE[1] = 140;
82
83     fgQualityRegionsB[2] = 180;
84     fgQualityRegionsE[2] = 210;
85   }
86   else
87   {
88     // TPC:   peak area - flat area - low stat area
89     fgQualityRegionsB[0] = 4;
90     fgQualityRegionsE[0] = 12;
91
92     fgQualityRegionsB[1] = 25;
93     fgQualityRegionsE[1] = 55;
94
95     fgQualityRegionsB[2] = 88;
96     fgQualityRegionsE[2] = 108;
97   }
98 }
99
100 //____________________________________________________________________
101 AliMultiplicityCorrection::AliMultiplicityCorrection() :
102   TNamed(), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
103 {
104   //
105   // default constructor
106   //
107
108   for (Int_t i = 0; i < kESDHists; ++i)
109     fMultiplicityESD[i] = 0;
110
111   for (Int_t i = 0; i < kMCHists; ++i)
112   {
113     fMultiplicityVtx[i] = 0;
114     fMultiplicityMB[i] = 0;
115     fMultiplicityINEL[i] = 0;
116   }
117
118   for (Int_t i = 0; i < kCorrHists; ++i)
119   {
120     fCorrelation[i] = 0;
121     fMultiplicityESDCorrected[i] = 0;
122   }
123 }
124
125 //____________________________________________________________________
126 AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const Char_t* title) :
127   TNamed(name, title), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
128 {
129   //
130   // named constructor
131   //
132
133   // do not add this hists to the directory
134   Bool_t oldStatus = TH1::AddDirectoryStatus();
135   TH1::AddDirectory(kFALSE);
136
137   /*Float_t binLimitsVtx[] = {-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10};
138   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,
139                           10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
140                           20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
141                           30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
142                           40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5,
143                           50.5, 55.5, 60.5, 65.5, 70.5, 75.5, 80.5, 85.5, 90.5, 95.5,
144                           100.5, 105.5, 110.5, 115.5, 120.5, 125.5, 130.5, 135.5, 140.5, 145.5,
145                           150.5, 160.5, 170.5, 180.5, 190.5, 200.5, 210.5, 220.5, 230.5, 240.5,
146                           250.5, 275.5, 300.5, 325.5, 350.5, 375.5, 400.5, 425.5, 450.5, 475.5,
147                           500.5 };
148                                 //525.5, 550.5, 575.5, 600.5, 625.5, 650.5, 675.5, 700.5, 725.5,
149                           //750.5, 775.5, 800.5, 825.5, 850.5, 875.5, 900.5, 925.5, 950.5, 975.5,
150                           //1000.5 };
151
152   #define VTXBINNING 10, binLimitsVtx
153   #define NBINNING fgMaxParams, binLimitsN*/
154
155   #define NBINNING 501, -0.5, 500.5
156   #define VTXBINNING 1, -10, 10
157
158   for (Int_t i = 0; i < kESDHists; ++i)
159     fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", VTXBINNING, NBINNING);
160
161   for (Int_t i = 0; i < kMCHists; ++i)
162   {
163     fMultiplicityVtx[i] = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("fMultiplicityVtx%d", i)));
164     fMultiplicityVtx[i]->SetTitle("fMultiplicityVtx;vtx-z;Npart");
165
166     fMultiplicityMB[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityMB%d", i)));
167     fMultiplicityMB[i]->SetTitle("fMultiplicityMB");
168
169     fMultiplicityINEL[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityINEL%d", i)));
170     fMultiplicityINEL[i]->SetTitle("fMultiplicityINEL");
171   }
172
173   for (Int_t i = 0; i < kCorrHists; ++i)
174   {
175     fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;Npart;Ntracks", VTXBINNING, NBINNING, NBINNING);
176     fMultiplicityESDCorrected[i] = new TH1F(Form("fMultiplicityESDCorrected%d", i), "fMultiplicityESDCorrected;Npart;dN/dN", NBINNING);
177   }
178
179   TH1::AddDirectory(oldStatus);
180 }
181
182 //____________________________________________________________________
183 AliMultiplicityCorrection::~AliMultiplicityCorrection()
184 {
185   //
186   // Destructor
187   //
188 }
189
190 //____________________________________________________________________
191 Long64_t AliMultiplicityCorrection::Merge(TCollection* list)
192 {
193   // Merge a list of AliMultiplicityCorrection objects with this (needed for
194   // PROOF).
195   // Returns the number of merged objects (including this).
196
197   if (!list)
198     return 0;
199
200   if (list->IsEmpty())
201     return 1;
202
203   TIterator* iter = list->MakeIterator();
204   TObject* obj;
205
206   // collections of all histograms
207   TList collections[kESDHists+kMCHists*3+kCorrHists*2];
208
209   Int_t count = 0;
210   while ((obj = iter->Next())) {
211
212     AliMultiplicityCorrection* entry = dynamic_cast<AliMultiplicityCorrection*> (obj);
213     if (entry == 0) 
214       continue;
215
216     for (Int_t i = 0; i < kESDHists; ++i)
217       collections[i].Add(entry->fMultiplicityESD[i]);
218
219     for (Int_t i = 0; i < kMCHists; ++i)
220     {
221       collections[kESDHists+i].Add(entry->fMultiplicityVtx[i]);
222       collections[kESDHists+kMCHists+i].Add(entry->fMultiplicityMB[i]);
223       collections[kESDHists+kMCHists*2+i].Add(entry->fMultiplicityINEL[i]);
224     }
225
226     for (Int_t i = 0; i < kCorrHists; ++i)
227       collections[kESDHists+kMCHists*3+i].Add(entry->fCorrelation[i]);
228
229     for (Int_t i = 0; i < kCorrHists; ++i)
230       collections[kESDHists+kMCHists*3+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]);
231
232     count++;
233   }
234
235   for (Int_t i = 0; i < kESDHists; ++i)
236     fMultiplicityESD[i]->Merge(&collections[i]);
237
238   for (Int_t i = 0; i < kMCHists; ++i)
239   {
240     fMultiplicityVtx[i]->Merge(&collections[kESDHists+i]);
241     fMultiplicityMB[i]->Merge(&collections[kESDHists+kMCHists+i]);
242     fMultiplicityINEL[i]->Merge(&collections[kESDHists+kMCHists*2+i]);
243   }
244
245   for (Int_t i = 0; i < kCorrHists; ++i)
246     fCorrelation[i]->Merge(&collections[kESDHists+kMCHists*3+i]);
247
248   for (Int_t i = 0; i < kCorrHists; ++i)
249     fMultiplicityESDCorrected[i]->Merge(&collections[kESDHists+kMCHists*3+kCorrHists+i]);
250
251   delete iter;
252
253   return count+1;
254 }
255
256 //____________________________________________________________________
257 Bool_t AliMultiplicityCorrection::LoadHistograms(const Char_t* dir)
258 {
259   //
260   // loads the histograms from a file
261   // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
262   //
263
264   if (!dir)
265     dir = GetName();
266
267   if (!gDirectory->cd(dir))
268     return kFALSE;
269
270   // TODO memory leak. old histograms needs to be deleted.
271
272   Bool_t success = kTRUE;
273
274   for (Int_t i = 0; i < kESDHists; ++i)
275   {
276     fMultiplicityESD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityESD[i]->GetName()));
277     if (!fMultiplicityESD[i])
278       success = kFALSE;
279   }
280
281   for (Int_t i = 0; i < kMCHists; ++i)
282   {
283     fMultiplicityVtx[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityVtx[i]->GetName()));
284     fMultiplicityMB[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityMB[i]->GetName()));
285     fMultiplicityINEL[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityINEL[i]->GetName()));
286     if (!fMultiplicityVtx[i] || !fMultiplicityMB[i] || !fMultiplicityINEL[i])
287       success = kFALSE;
288   }
289
290   for (Int_t i = 0; i < kCorrHists; ++i)
291   {
292     fCorrelation[i] = dynamic_cast<TH3F*> (gDirectory->Get(fCorrelation[i]->GetName()));
293     if (!fCorrelation[i])
294       success = kFALSE;
295     fMultiplicityESDCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fMultiplicityESDCorrected[i]->GetName()));
296     if (!fMultiplicityESDCorrected[i])
297       success = kFALSE;
298   }
299
300   gDirectory->cd("..");
301
302   return success;
303 }
304
305 //____________________________________________________________________
306 void AliMultiplicityCorrection::SaveHistograms()
307 {
308   //
309   // saves the histograms
310   //
311
312   gDirectory->mkdir(GetName());
313   gDirectory->cd(GetName());
314
315   for (Int_t i = 0; i < kESDHists; ++i)
316     if (fMultiplicityESD[i])
317       fMultiplicityESD[i]->Write();
318
319   for (Int_t i = 0; i < kMCHists; ++i)
320   {
321     if (fMultiplicityVtx[i])
322       fMultiplicityVtx[i]->Write();
323     if (fMultiplicityMB[i])
324       fMultiplicityMB[i]->Write();
325     if (fMultiplicityINEL[i])
326       fMultiplicityINEL[i]->Write();
327   }
328
329   for (Int_t i = 0; i < kCorrHists; ++i)
330   {
331     if (fCorrelation[i])
332       fCorrelation[i]->Write();
333     if (fMultiplicityESDCorrected[i])
334       fMultiplicityESDCorrected[i]->Write();
335   }
336
337   gDirectory->cd("..");
338 }
339
340 //____________________________________________________________________
341 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)
342 {
343   //
344   // Fills an event from MC
345   //
346
347   if (triggered)
348   {
349     fMultiplicityMB[0]->Fill(vtx, generated05);
350     fMultiplicityMB[1]->Fill(vtx, generated10);
351     fMultiplicityMB[2]->Fill(vtx, generated15);
352     fMultiplicityMB[3]->Fill(vtx, generated20);
353     fMultiplicityMB[4]->Fill(vtx, generatedAll);
354
355     if (vertex)
356     {
357       fMultiplicityVtx[0]->Fill(vtx, generated05);
358       fMultiplicityVtx[1]->Fill(vtx, generated10);
359       fMultiplicityVtx[2]->Fill(vtx, generated15);
360       fMultiplicityVtx[3]->Fill(vtx, generated20);
361       fMultiplicityVtx[4]->Fill(vtx, generatedAll);
362     }
363   }
364
365   fMultiplicityINEL[0]->Fill(vtx, generated05);
366   fMultiplicityINEL[1]->Fill(vtx, generated10);
367   fMultiplicityINEL[2]->Fill(vtx, generated15);
368   fMultiplicityINEL[3]->Fill(vtx, generated20);
369   fMultiplicityINEL[4]->Fill(vtx, generatedAll);
370 }
371
372 //____________________________________________________________________
373 void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20)
374 {
375   //
376   // Fills an event from ESD
377   //
378
379   fMultiplicityESD[0]->Fill(vtx, measured05);
380   fMultiplicityESD[1]->Fill(vtx, measured10);
381   fMultiplicityESD[2]->Fill(vtx, measured15);
382   fMultiplicityESD[3]->Fill(vtx, measured20);
383 }
384
385 //____________________________________________________________________
386 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)
387 {
388   //
389   // Fills an event into the correlation map with the information from MC and ESD
390   //
391
392   fCorrelation[0]->Fill(vtx, generated05, measured05);
393   fCorrelation[1]->Fill(vtx, generated10, measured10);
394   fCorrelation[2]->Fill(vtx, generated15, measured15);
395   fCorrelation[3]->Fill(vtx, generated20, measured20);
396
397   fCorrelation[4]->Fill(vtx, generatedAll, measured05);
398   fCorrelation[5]->Fill(vtx, generatedAll, measured10);
399   fCorrelation[6]->Fill(vtx, generatedAll, measured15);
400   fCorrelation[7]->Fill(vtx, generatedAll, measured20);
401 }
402
403 //____________________________________________________________________
404 Double_t AliMultiplicityCorrection::RegularizationPol0(TVectorD& params)
405 {
406   // homogenity term for minuit fitting
407   // pure function of the parameters
408   // prefers constant function (pol0)
409
410   Double_t chi2 = 0;
411
412   // ignore the first bin here. on purpose...
413   for (Int_t i=2; i<fgMaxParams; ++i)
414   {
415     Double_t right  = params[i];
416     Double_t left   = params[i-1];
417
418     if (right != 0)
419     {
420       Double_t diff = 1 - left / right;
421       chi2 += diff * diff;
422     }
423   }
424
425   return chi2 / 100.0;
426 }
427
428 //____________________________________________________________________
429 Double_t AliMultiplicityCorrection::RegularizationPol1(TVectorD& params)
430 {
431   // homogenity term for minuit fitting
432   // pure function of the parameters
433   // prefers linear function (pol1)
434
435   Double_t chi2 = 0;
436
437   // ignore the first bin here. on purpose...
438   for (Int_t i=2+1; i<fgMaxParams; ++i)
439   {
440     if (params[i-1] == 0)
441       continue;
442
443     Double_t right  = params[i];
444     Double_t middle = params[i-1];
445     Double_t left   = params[i-2];
446
447     Double_t der1 = (right - middle);
448     Double_t der2 = (middle - left);
449
450     Double_t diff = (der1 - der2) / middle;
451
452     chi2 += diff * diff;
453   }
454
455   return chi2;
456 }
457
458 //____________________________________________________________________
459 Double_t AliMultiplicityCorrection::RegularizationLog(TVectorD& params)
460 {
461   // homogenity term for minuit fitting
462   // pure function of the parameters
463   // prefers linear function (pol1)
464
465   Double_t chi2 = 0;
466
467   /*Float_t der[fgMaxParams];
468
469   for (Int_t i=0; i<fgMaxParams-1; ++i)
470     der[i] = params[i+1] - params[i];
471
472   for (Int_t i=0; i<fgMaxParams-6; ++i)
473   {
474     for (Int_t j=1; j<=5; ++j)
475     {
476       Double_t diff = der[i+j] - der[i];
477       chi2 += diff * diff;
478     }
479   }*/
480
481   // ignore the first bin here. on purpose...
482   for (Int_t i=2+1; i<fgMaxParams; ++i)
483   {
484     if (params[i-1] == 0)
485       continue;
486
487     Double_t right  = log(params[i]);
488     Double_t middle = log(params[i-1]);
489     Double_t left   = log(params[i-2]);
490
491     Double_t der1 = (right - middle);
492     Double_t der2 = (middle - left);
493
494     Double_t diff = (der1 - der2) / middle;
495
496     chi2 += diff * diff;
497   }
498
499   return chi2 * 100;
500 }
501
502 //____________________________________________________________________
503 Double_t AliMultiplicityCorrection::RegularizationTotalCurvature(TVectorD& params)
504 {
505   // homogenity term for minuit fitting
506   // pure function of the parameters
507   // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments,
508   // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp.
509
510   Double_t chi2 = 0;
511
512   // ignore the first bin here. on purpose...
513   for (Int_t i=3; i<fgMaxParams; ++i)
514   {
515     Double_t right  = params[i];
516     Double_t middle = params[i-1];
517     Double_t left   = params[i-2];
518
519     Double_t der1 = (right - middle);
520     Double_t der2 = (middle - left);
521
522     Double_t diff = (der1 - der2);
523
524     chi2 += diff * diff;
525   }
526
527   return chi2 * 1e4;
528 }
529
530 //____________________________________________________________________
531 Double_t AliMultiplicityCorrection::RegularizationEntropy(TVectorD& params)
532 {
533   // homogenity term for minuit fitting
534   // pure function of the parameters
535   // calculates entropy, from
536   // The method of reduced cross-entropy (M. Schmelling 1993)
537
538   Double_t paramSum = 0;
539   // ignore the first bin here. on purpose...
540   for (Int_t i=1; i<fgMaxParams; ++i)
541     paramSum += params[i];
542
543   Double_t chi2 = 0;
544   for (Int_t i=1; i<fgMaxParams; ++i)
545   {
546     //Double_t tmp = params[i] / paramSum;
547     Double_t tmp = params[i];
548     if (tmp > 0 && (*fEntropyAPriori)[i] > 0)
549       chi2 += tmp * TMath::Log(tmp / (*fEntropyAPriori)[i]);
550   }
551
552   return 10.0 + chi2;
553 }
554
555 //____________________________________________________________________
556 void AliMultiplicityCorrection::MinuitNBD(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3)
557 {
558   //
559   // fit function for minuit
560   // does: nbd
561   // func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 50);
562   // func->SetParNames("scaling", "averagen", "k");
563   // func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
564   // func->SetParLimits(1, 0.001, 1000);
565   // func->SetParLimits(2, 0.001, 1000);
566   //
567
568   fNBD->SetParameters(params[0], params[1], params[2]);
569
570   Double_t params2[fgMaxParams];
571
572   for (Int_t i=0; i<fgMaxParams; ++i)
573     params2[i] = fNBD->Eval(i);
574
575   MinuitFitFunction(unused1, unused2, chi2, params2, unused3);
576
577   printf("%f %f %f --> %f\n", params[0], params[1], params[2], chi2);
578 }
579
580 //____________________________________________________________________
581 void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t)
582 {
583   //
584   // fit function for minuit
585   // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix
586   //
587
588   // d
589   static TVectorD paramsVector(fgMaxParams);
590   for (Int_t i=0; i<fgMaxParams; ++i)
591     paramsVector[i] = params[i] * params[i];
592
593   // calculate penalty factor
594   Double_t penaltyVal = 0;
595   switch (fRegularizationType)
596   {
597     case kNone:       break;
598     case kPol0:       penaltyVal = RegularizationPol0(paramsVector); break;
599     case kPol1:       penaltyVal = RegularizationPol1(paramsVector); break;
600     case kCurvature:  penaltyVal = RegularizationTotalCurvature(paramsVector); break;
601     case kEntropy:    penaltyVal = RegularizationEntropy(paramsVector); break;
602     case kLog:        penaltyVal = RegularizationLog(paramsVector); break;
603   }
604
605   //if (penaltyVal > 1e10)
606   //  paramsVector2.Print();
607
608   penaltyVal *= fRegularizationWeight;
609
610   // Ad
611   TVectorD measGuessVector(fgMaxInput);
612   measGuessVector = (*fCorrelationMatrix) * paramsVector;
613
614   // Ad - m
615   measGuessVector -= (*fCurrentESDVector);
616
617   TVectorD copy(measGuessVector);
618
619   // (Ad - m) W
620   // this step can be optimized because currently only the diagonal elements of fCorrelationCovarianceMatrix are used
621   // normal way is like this:
622   // measGuessVector *= (*fCorrelationCovarianceMatrix);
623   // optimized way like this:
624   for (Int_t i=0; i<fgMaxParams; ++i)
625     measGuessVector[i] *= (*fCorrelationCovarianceMatrix)(i, i);
626
627   //measGuessVector.Print();
628
629   // (Ad - m) W (Ad - m)
630   // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
631   // big number ((*fCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
632   Double_t chi2FromFit = measGuessVector * copy * 1e6;
633
634   chi2 = chi2FromFit + penaltyVal;
635
636   static Int_t callCount = 0;
637   if ((callCount++ % 10000) == 0)
638     printf("%d) %f %f --> %f\n", callCount, chi2FromFit, penaltyVal, chi2);
639 }
640
641 //____________________________________________________________________
642 void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t createBigBin)
643 {
644   //
645   // fills fCurrentESD, fCurrentCorrelation
646   // resets fMultiplicityESDCorrected
647   //
648
649   Bool_t display = kFALSE;
650
651   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
652
653   fMultiplicityESDCorrected[correlationID]->Reset();
654   fMultiplicityESDCorrected[correlationID]->Sumw2();
655
656   fCurrentESD = fMultiplicityESD[inputRange]->ProjectionY("fCurrentESD"); //"fCurrentESD", -1, -1, "e");
657   fCurrentESD->Sumw2();
658
659   // empty under/overflow bins in x, otherwise Project3D takes them into account
660   TH3* hist = fCorrelation[correlationID];
661   for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
662   {
663     for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
664     {
665       hist->SetBinContent(0, y, z, 0);
666       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
667     }
668   }
669
670   fCurrentCorrelation = hist->Project3D("zy");
671   //((TH2*) fCurrentCorrelation)->Rebin2D(2, 1);
672   //fMultiplicityESDCorrected[correlationID]->Rebin(2);
673   fCurrentCorrelation->Sumw2();
674
675   if (createBigBin)
676   {
677     Int_t maxBin = 0;
678     for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i)
679     {
680       if (fCurrentESD->GetBinContent(i) <= 5)
681       {
682         maxBin = i;
683         break;
684       }
685     }
686
687     if (maxBin > 0)
688     {
689       TCanvas* canvas = 0;
690
691       if (display)
692       {
693         canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800);
694         canvas->Divide(2, 2);
695
696         canvas->cd(1);
697         fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
698         fCurrentESD->SetStats(kFALSE);
699         fCurrentESD->SetTitle(";measured multiplicity");
700         fCurrentESD->DrawCopy();
701         gPad->SetLogy();
702
703         canvas->cd(2);
704         fCurrentCorrelation->GetXaxis()->SetRangeUser(0, 250);
705         fCurrentCorrelation->GetYaxis()->SetRangeUser(0, 250);
706         fCurrentCorrelation->SetTitle(";true multiplicity;measured multiplicity");
707         fCurrentCorrelation->SetStats(kFALSE);
708         fCurrentCorrelation->DrawCopy("COLZ");
709       }
710
711       printf("Bin limit in measured spectrum is %d.\n", maxBin);
712       fCurrentESD->SetBinContent(maxBin, fCurrentESD->Integral(maxBin, fCurrentESD->GetNbinsX()));
713       for (Int_t i=maxBin+1; i<=fCurrentESD->GetNbinsX(); ++i)
714       {
715         fCurrentESD->SetBinContent(i, 0);
716         fCurrentESD->SetBinError(i, 0);
717       }
718       // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
719       fCurrentESD->SetBinError(maxBin, TMath::Sqrt(fCurrentESD->GetBinContent(maxBin)));
720
721       printf("This bin has now %f +- %f entries\n", fCurrentESD->GetBinContent(maxBin), fCurrentESD->GetBinError(maxBin));
722
723       for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
724       {
725         fCurrentCorrelation->SetBinContent(i, maxBin, fCurrentCorrelation->Integral(i, i, maxBin, fCurrentCorrelation->GetNbinsY()));
726         // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
727         fCurrentCorrelation->SetBinError(i, maxBin, TMath::Sqrt(fCurrentCorrelation->GetBinContent(i, maxBin)));
728
729         for (Int_t j=maxBin+1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
730         {
731           fCurrentCorrelation->SetBinContent(i, j, 0);
732           fCurrentCorrelation->SetBinError(i, j, 0);
733         }
734       }
735
736       printf("Adjusted correlation matrix!\n");
737
738       if (canvas)
739       {
740         canvas->cd(3);
741         fCurrentESD->DrawCopy();
742         gPad->SetLogy();
743
744         canvas->cd(4);
745         fCurrentCorrelation->DrawCopy("COLZ");
746       }
747     }
748   }
749
750 #if 0 // does not help
751   // null bins with one entry
752   Int_t nNulledBins = 0;
753   for (Int_t x=1; x<=fCurrentCorrelation->GetXaxis()->GetNbins(); ++x)
754     for (Int_t y=1; y<=fCurrentCorrelation->GetYaxis()->GetNbins(); ++y)
755     {
756       if (fCurrentCorrelation->GetBinContent(x, y) == 1)
757       {
758         fCurrentCorrelation->SetBinContent(x, y, 0);
759         fCurrentCorrelation->SetBinError(x, y, 0);
760
761         ++nNulledBins;
762       }
763     }
764   Printf("Nulled %d bins", nNulledBins);
765 #endif
766
767   fCurrentEfficiency = GetEfficiency(inputRange, eventType);
768   //fCurrentEfficiency->Rebin(2);
769   //fCurrentEfficiency->Scale(0.5);
770 }
771
772 //____________________________________________________________________
773 TH1* AliMultiplicityCorrection::GetEfficiency(Int_t inputRange, EventType eventType)
774 {
775   //
776   // calculates efficiency for given event type
777   //
778
779   TH1* divisor = 0;
780   switch (eventType)
781   {
782     case kTrVtx : divisor = fMultiplicityVtx[inputRange]->ProjectionY("divisor", -1, -1, "e"); break;
783     case kMB: divisor = fMultiplicityMB[inputRange]->ProjectionY("divisor", -1, -1, "e"); break;
784     case kINEL: divisor = fMultiplicityINEL[inputRange]->ProjectionY("divisor", -1, -1, "e"); break;
785   }
786   TH1* eff = fMultiplicityVtx[inputRange]->ProjectionY("CurrentEfficiency", -1, -1, "e");
787   eff->Divide(eff, divisor, 1, 1, "B");
788   return eff;
789 }
790
791 //____________________________________________________________________
792 void AliMultiplicityCorrection::SetRegularizationParameters(RegularizationType type, Float_t weight)
793 {
794   fRegularizationType = type;
795   fRegularizationWeight = weight;
796
797   printf("AliMultiplicityCorrection::SetRegularizationParameters --> Regularization set to %d with weight %f\n", (Int_t) type, weight);
798 }
799
800 //____________________________________________________________________
801 Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* inputDist)
802 {
803   //
804   // correct spectrum using minuit chi2 method
805   //
806   // if check is kTRUE the input MC solution (by definition the right solution) is used
807   // no fit is made and just the chi2 is calculated
808   //
809
810   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
811   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
812
813   SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kTRUE);
814   //normalize ESD
815   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
816
817   fCorrelationMatrix = new TMatrixD(fgMaxInput, fgMaxParams);
818   fCorrelationCovarianceMatrix = new TMatrixD(fgMaxInput, fgMaxInput);
819   fCurrentESDVector = new TVectorD(fgMaxInput);
820   fEntropyAPriori = new TVectorD(fgMaxParams);
821
822   /*new TCanvas; fCurrentESD->DrawCopy();
823   fCurrentESD = ((TH2*) fCurrentCorrelation)->ProjectionY("check-proj2");
824   fCurrentESD->Sumw2();
825   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
826   fCurrentESD->SetLineColor(2);
827   fCurrentESD->DrawCopy("SAME");*/
828
829   // normalize correction for given nPart
830   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
831   {
832     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
833     if (sum <= 0)
834       continue;
835     Float_t maxValue = 0;
836     Int_t maxBin = -1;
837     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
838     {
839       // find most probably value
840       if (maxValue < fCurrentCorrelation->GetBinContent(i, j))
841       {
842         maxValue = fCurrentCorrelation->GetBinContent(i, j);
843         maxBin = j;
844       }
845
846       // npart sum to 1
847       fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
848       fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
849
850       if (i <= fgMaxParams && j <= fgMaxInput)
851         (*fCorrelationMatrix)(j-1, i-1) = fCurrentCorrelation->GetBinContent(i, j);
852     }
853
854     //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
855   }
856
857   // Initialize TMinuit via generic fitter interface
858   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgMaxParams);
859   Double_t arglist[100];
860
861   // disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find out why...
862   arglist[0] = 0;
863   minuit->ExecuteCommand("SET PRINT", arglist, 1);
864
865   // however, enable warnings
866   //minuit->ExecuteCommand("SET WAR", arglist, 0);
867
868   // set minimization function
869   minuit->SetFCN(MinuitFitFunction);
870
871   for (Int_t i=0; i<fgMaxParams; i++)
872     (*fEntropyAPriori)[i] = 1;
873
874   if (inputDist)
875   {
876     printf("Using different starting conditions...\n");
877     new TCanvas;
878     inputDist->DrawCopy();
879
880     inputDist->Scale(1.0 / inputDist->Integral());
881     for (Int_t i=0; i<fgMaxParams; i++)
882       if (inputDist->GetBinContent(i+1) > 0)
883         (*fEntropyAPriori)[i] = inputDist->GetBinContent(i+1);
884   }
885   else
886     inputDist = fCurrentESD;
887
888
889   //Float_t minStartValue = 0; //1e-3;
890
891   //new TCanvas; fMultiplicityVtx[mcTarget]->Draw("COLZ");
892   TH1* proj = fMultiplicityVtx[mcTarget]->ProjectionY("check-proj");
893   //proj->Rebin(2);
894   proj->Scale(1.0 / proj->Integral());
895
896   Double_t results[fgMaxParams+1];
897   for (Int_t i=0; i<fgMaxParams; ++i)
898   {
899     results[i] = inputDist->GetBinContent(i+1);
900
901     if (check)
902       results[i] = proj->GetBinContent(i+1);
903
904     // minimum value
905     results[i] = TMath::Max(results[i], 1e-3);
906
907     Float_t stepSize = 0.1;
908
909     // minuit sees squared values to prevent it from going negative...
910     results[i] = TMath::Sqrt(results[i]);
911     //stepSize /= results[i]; // keep relative step size
912
913     minuit->SetParameter(i, Form("param%d", i), results[i], stepSize, 0, 0);
914   }
915   //minuit->SetParameter(fgMaxParams, "alpha", 1e4, 1, 0, 0);
916   // bin 0 is filled with value from bin 1 (otherwise it's 0)
917   //minuit->SetParameter(0, "param0", (results[1] > minStartValue) ? results[1] : minStartValue, 0.1, 0, 0);
918   //results[0] = 0;
919   //minuit->SetParameter(0, "param0", 0, 0, 0, 0);
920
921   for (Int_t i=0; i<fgMaxInput; ++i)
922   {
923     //printf("%f %f %f\n", inputDist->GetBinContent(i+1), TMath::Sqrt(inputDist->GetBinContent(i+1)), inputDist->GetBinError(i+1));
924
925     (*fCurrentESDVector)[i] = fCurrentESD->GetBinContent(i+1);
926     if (fCurrentESD->GetBinError(i+1) > 0)
927       (*fCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / fCurrentESD->GetBinError(i+1) / fCurrentESD->GetBinError(i+1);
928
929     if ((*fCorrelationCovarianceMatrix)(i, i) > 1e7)
930       (*fCorrelationCovarianceMatrix)(i, i) = 0;
931
932     //printf("%d --> %e\n", i, (*fCorrelationCovarianceMatrix)(i, i));
933   }
934
935   Int_t dummy;
936   Double_t chi2;
937   MinuitFitFunction(dummy, 0, chi2, results, 0);
938   printf("Chi2 of initial parameters is = %f\n", chi2);
939
940   if (check)
941     return -1;
942
943   // first param is number of iterations, second is precision....
944   arglist[0] = 1e6;
945   //arglist[1] = 1e-5;
946   //minuit->ExecuteCommand("SCAN", arglist, 0);
947   Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1);
948   printf("MINUIT status is %d\n", status);
949   //minuit->ExecuteCommand("MIGRAD", arglist, 1);
950   //minuit->ExecuteCommand("MIGRAD", arglist, 1);
951   //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!");
952   //minuit->ExecuteCommand("MINOS", arglist, 0);
953
954   for (Int_t i=0; i<fgMaxParams; ++i)
955   {
956     if (fCurrentEfficiency->GetBinContent(i+1) > 0)
957     {
958       fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, minuit->GetParameter(i) * minuit->GetParameter(i) / fCurrentEfficiency->GetBinContent(i+1));
959       // error is : (relError) * (value) = (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i))
960       fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, minuit->GetParError(i) * minuit->GetParameter(i) /  fCurrentEfficiency->GetBinContent(i+1));
961     }
962     else
963     {
964       fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, 0);
965       fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0);
966     }
967   }
968
969   return status;
970 }
971
972 //____________________________________________________________________
973 void AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace)
974 {
975   //
976   // correct spectrum using minuit chi2 method applying a NBD fit
977   //
978
979   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
980
981   SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE);
982   //normalize ESD
983   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
984
985   fCorrelationMatrix = new TMatrixD(fgMaxParams, fgMaxParams);
986   fCorrelationCovarianceMatrix = new TMatrixD(fgMaxParams, fgMaxParams);
987   fCurrentESDVector = new TVectorD(fgMaxParams);
988
989   // normalize correction for given nPart
990   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
991   {
992     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
993     if (sum <= 0)
994       continue;
995     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
996     {
997       // npart sum to 1
998       fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
999       fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1000
1001       if (i <= fgMaxParams && j <= fgMaxParams)
1002         (*fCorrelationMatrix)(j-1, i-1) = fCurrentCorrelation->GetBinContent(i, j);
1003     }
1004   }
1005
1006   for (Int_t i=0; i<fgMaxParams; ++i)
1007   {
1008     (*fCurrentESDVector)[i] = fCurrentESD->GetBinContent(i+1);
1009     if (fCurrentESD->GetBinError(i+1) > 0)
1010       (*fCorrelationCovarianceMatrix)(i, i) = 1.0 / (fCurrentESD->GetBinError(i+1) * fCurrentESD->GetBinError(i+1));
1011   }
1012
1013   // Create NBD function
1014   if (!fNBD)
1015   {
1016     fNBD = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 250);
1017     fNBD->SetParNames("scaling", "averagen", "k");
1018   }
1019
1020   // Initialize TMinuit via generic fitter interface
1021   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, 3);
1022
1023   minuit->SetFCN(MinuitNBD);
1024   minuit->SetParameter(0, "param0", 1, 0.1, 0.001, fCurrentESD->GetMaximum() * 1000);
1025   minuit->SetParameter(1, "param1", 10, 1, 0.001, 1000);
1026   minuit->SetParameter(2, "param2", 2, 1, 0.001, 1000);
1027
1028   Double_t arglist[100];
1029   arglist[0] = 0;
1030   minuit->ExecuteCommand("SET PRINT", arglist, 1);
1031   minuit->ExecuteCommand("MIGRAD", arglist, 0);
1032   //minuit->ExecuteCommand("MINOS", arglist, 0);
1033
1034   fNBD->SetParameters(minuit->GetParameter(0), minuit->GetParameter(1), minuit->GetParameter(2));
1035
1036   for (Int_t i=0; i<fgMaxParams; ++i)
1037   {
1038     printf("%d : %f\n", i, fNBD->Eval(i));
1039     if (fNBD->Eval(i) > 0)
1040     {
1041       fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, fNBD->Eval(i));
1042       fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0);
1043     }
1044   }
1045 }
1046
1047 //____________________________________________________________________
1048 void AliMultiplicityCorrection::NormalizeToBinWidth(TH1* hist)
1049 {
1050   //
1051   // normalizes a 1-d histogram to its bin width
1052   //
1053
1054   for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
1055   {
1056     hist->SetBinContent(i, hist->GetBinContent(i) / hist->GetBinWidth(i));
1057     hist->SetBinError(i, hist->GetBinError(i) / hist->GetBinWidth(i));
1058   }
1059 }
1060
1061 //____________________________________________________________________
1062 void AliMultiplicityCorrection::NormalizeToBinWidth(TH2* hist)
1063 {
1064   //
1065   // normalizes a 2-d histogram to its bin width (x width * y width)
1066   //
1067
1068   for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
1069     for (Int_t j=1; j<=hist->GetNbinsY(); ++j)
1070     {
1071       Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j);
1072       hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor);
1073       hist->SetBinError(i, j, hist->GetBinError(i, j) / factor);
1074     }
1075 }
1076
1077 //____________________________________________________________________
1078 void AliMultiplicityCorrection::DrawHistograms()
1079 {
1080   printf("ESD:\n");
1081
1082   TCanvas* canvas1 = new TCanvas("fMultiplicityESD", "fMultiplicityESD", 900, 600);
1083   canvas1->Divide(3, 2);
1084   for (Int_t i = 0; i < kESDHists; ++i)
1085   {
1086     canvas1->cd(i+1);
1087     fMultiplicityESD[i]->DrawCopy("COLZ");
1088     printf("%d --> %f\n", i, (Float_t) fMultiplicityESD[i]->ProjectionY()->GetMean());
1089   }
1090
1091   printf("Vtx:\n");
1092
1093   TCanvas* canvas2 = new TCanvas("fMultiplicityMC", "fMultiplicityMC", 900, 600);
1094   canvas2->Divide(3, 2);
1095   for (Int_t i = 0; i < kMCHists; ++i)
1096   {
1097     canvas2->cd(i+1);
1098     fMultiplicityVtx[i]->DrawCopy("COLZ");
1099     printf("%d --> %f\n", i, fMultiplicityVtx[i]->ProjectionY()->GetMean());
1100   }
1101
1102   TCanvas* canvas3 = new TCanvas("fCorrelation", "fCorrelation", 900, 900);
1103   canvas3->Divide(3, 3);
1104   for (Int_t i = 0; i < kCorrHists; ++i)
1105   {
1106     canvas3->cd(i+1);
1107     TH3* hist = dynamic_cast<TH3*> (fCorrelation[i]->Clone());
1108     for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
1109     {
1110       for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
1111       {
1112         hist->SetBinContent(0, y, z, 0);
1113         hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
1114       }
1115     }
1116     TH1* proj = hist->Project3D("zy");
1117     proj->DrawCopy("COLZ");
1118   }
1119 }
1120
1121 //____________________________________________________________________
1122 void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t normalizeESD, TH1* mcHist, Bool_t simple)
1123 {
1124   //mcHist->Rebin(2);
1125
1126   Int_t esdCorrId = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1127
1128   TString tmpStr;
1129   tmpStr.Form("%s_DrawComparison_%d", name, esdCorrId);
1130
1131   if (fMultiplicityESDCorrected[esdCorrId]->Integral() == 0)
1132   {
1133     printf("ERROR. Unfolded histogram is empty\n");
1134     return;
1135   }
1136
1137   //regain measured distribution used for unfolding, because the bins at high mult. were modified in SetupCurrentHists
1138   fCurrentESD = fMultiplicityESD[esdCorrId]->ProjectionY();
1139   fCurrentESD->Sumw2();
1140   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
1141
1142   // normalize unfolded result to 1
1143   fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral());
1144
1145   //fCurrentESD->Scale(mcHist->Integral(2, 200));
1146
1147   //new TCanvas;
1148   /*TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
1149   ratio->Divide(mcHist);
1150   ratio->Draw("HIST");
1151   ratio->Fit("pol0", "W0", "", 20, 120);
1152   Float_t scalingFactor = ratio->GetFunction("pol0")->GetParameter(0);
1153   delete ratio;
1154   mcHist->Scale(scalingFactor);*/
1155
1156   // find bin with <= 100 entries. this is used as limit for the chi2 calculation
1157   Int_t mcBinLimit = 0;
1158   for (Int_t i=20; i<mcHist->GetNbinsX(); ++i)
1159   {
1160     if (mcHist->GetBinContent(i) > 50)
1161     {
1162       mcBinLimit = i;
1163     }
1164     else
1165       break;
1166   }
1167   Printf("AliMultiplicityCorrection::DrawComparison: MC bin limit is %d", mcBinLimit);
1168
1169   // scale to 1
1170   mcHist->Sumw2();
1171   mcHist->Scale(1.0 / mcHist->Integral());
1172
1173   // calculate residual
1174
1175   // for that we convolute the response matrix with the gathered result
1176   // if normalizeESD is not set, the histogram is already normalized, this needs to be passed to CalculateMultiplicityESD
1177   TH1* tmpESDEfficiencyRecorrected = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("tmpESDEfficiencyRecorrected");
1178   tmpESDEfficiencyRecorrected->Multiply(fCurrentEfficiency);
1179   TH2* convoluted = CalculateMultiplicityESD(tmpESDEfficiencyRecorrected, esdCorrId);
1180   TH1* convolutedProj = convoluted->ProjectionY("convolutedProj", -1, -1, "e");
1181   if (convolutedProj->Integral() > 0)
1182   {
1183     convolutedProj->Scale(1.0 / convolutedProj->Integral());
1184   }
1185   else
1186     printf("ERROR: convolutedProj is empty. Something went wrong calculating the convoluted histogram.\n");
1187
1188   //NormalizeToBinWidth(proj2);
1189
1190   TH1* residual = (TH1*) convolutedProj->Clone("residual");
1191   residual->SetTitle("Residuals;Ntracks;(folded unfolded measured - measured) / e");
1192
1193   residual->Add(fCurrentESD, -1);
1194   //residual->Divide(residual, fCurrentESD, 1, 1, "B");
1195
1196   TH1* residualHist = new TH1F("residualHist", "residualHist", 51, -5, 5);
1197
1198   // TODO fix errors
1199   Float_t chi2 = 0;
1200   for (Int_t i=1; i<=residual->GetNbinsX(); ++i)
1201   {
1202     if (fCurrentESD->GetBinError(i) > 0)
1203     {
1204       Float_t value = residual->GetBinContent(i) / fCurrentESD->GetBinError(i);
1205       if (i > 1)
1206         chi2 += value * value;
1207       residual->SetBinContent(i, value);
1208       residualHist->Fill(value);
1209     }
1210     else
1211     {
1212       //printf("Residual bin %d set to 0\n", i);
1213       residual->SetBinContent(i, 0);
1214     }
1215     convolutedProj->SetBinError(i, 0);
1216     residual->SetBinError(i, 0);
1217
1218     if (i == 200)
1219       fLastChi2Residuals = chi2;
1220   }
1221
1222   new TCanvas; residualHist->DrawCopy();
1223
1224   //residualHist->Fit("gaus", "N");
1225   //delete residualHist;
1226
1227   printf("Difference (Residuals) is %f for bin 2-200\n", chi2);
1228
1229   TCanvas* canvas1 = 0;
1230   if (simple)
1231   {
1232     canvas1 = new TCanvas(tmpStr, tmpStr, 900, 400);
1233     canvas1->Divide(2, 1);
1234   }
1235   else
1236   {
1237     canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 1200);
1238     canvas1->Divide(2, 3);
1239   }
1240
1241   canvas1->cd(1);
1242   TH1* proj = (TH1*) mcHist->Clone("proj");
1243   NormalizeToBinWidth(proj);
1244
1245   if (normalizeESD)
1246     NormalizeToBinWidth(fMultiplicityESDCorrected[esdCorrId]);
1247
1248   proj->GetXaxis()->SetRangeUser(0, 200);
1249   proj->SetTitle(";true multiplicity;Entries");
1250   proj->SetStats(kFALSE);
1251   proj->DrawCopy("HIST");
1252   gPad->SetLogy();
1253
1254   //fMultiplicityESDCorrected[esdCorrId]->SetMarkerStyle(3);
1255   fMultiplicityESDCorrected[esdCorrId]->SetLineColor(2);
1256   fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAME HIST E");
1257
1258   TLegend* legend = new TLegend(0.3, 0.8, 0.93, 0.93);
1259   legend->AddEntry(proj, "true distribution");
1260   legend->AddEntry(fMultiplicityESDCorrected[esdCorrId], "unfolded distribution");
1261   legend->SetFillColor(0);
1262   legend->Draw();
1263   // unfortunately does not work. maybe a bug? --> legend->SetTextSizePixels(14);
1264
1265   canvas1->cd(2);
1266
1267   gPad->SetLogy();
1268   fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
1269   //fCurrentESD->SetLineColor(2);
1270   fCurrentESD->SetTitle(";measured multiplicity;Entries");
1271   fCurrentESD->SetStats(kFALSE);
1272   fCurrentESD->DrawCopy("HIST E");
1273
1274   convolutedProj->SetLineColor(2);
1275   //proj2->SetMarkerColor(2);
1276   //proj2->SetMarkerStyle(5);
1277   convolutedProj->DrawCopy("HIST SAME");
1278
1279   legend = new TLegend(0.3, 0.8, 0.93, 0.93);
1280   legend->AddEntry(fCurrentESD, "measured distribution");
1281   legend->AddEntry(convolutedProj, "R #otimes unfolded distribution");
1282   legend->SetFillColor(0);
1283   legend->Draw();
1284
1285   //TH1* diffMCUnfolded = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded"));
1286   //diffMCUnfolded->Add(fMultiplicityESDCorrected[esdCorrId], -1);
1287
1288   /*Float_t chi2 = 0;
1289   Float_t chi = 0;
1290   fLastChi2MCLimit = 0;
1291   Int_t limit = 0;
1292   for (Int_t i=2; i<=200; ++i)
1293   {
1294     if (proj->GetBinContent(i) != 0)
1295     {
1296       Float_t value = (proj->GetBinContent(i) - fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i)) / proj->GetBinContent(i);
1297       chi2 += value * value;
1298       chi += TMath::Abs(value);
1299
1300       //printf("%d %f\n", i, chi);
1301
1302       if (chi2 < 0.2)
1303         fLastChi2MCLimit = i;
1304
1305       if (chi < 3)
1306         limit = i;
1307
1308     }
1309   }*/
1310
1311   /*chi2 = 0;
1312   Float_t chi = 0;
1313   Int_t limit = 0;
1314   for (Int_t i=1; i<=diffMCUnfolded->GetNbinsX(); ++i)
1315   {
1316     if (fMultiplicityESDCorrected[esdCorrId]->GetBinError(i) > 0)
1317     {
1318       Double_t value = diffMCUnfolded->GetBinContent(i) / fMultiplicityESDCorrected[esdCorrId]->GetBinError(i);
1319       if (value > 1e8)
1320         value = 1e8; //prevent arithmetic exception
1321       else if (value < -1e8)
1322         value = -1e8;
1323       if (i > 1)
1324       {
1325         chi2 += value * value;
1326         chi += TMath::Abs(value);
1327       }
1328       diffMCUnfolded->SetBinContent(i, value);
1329     }
1330     else
1331     {
1332       //printf("diffMCUnfolded bin %d set to 0\n", i);
1333       diffMCUnfolded->SetBinContent(i, 0);
1334     }
1335     if (chi2 < 1000)
1336       fLastChi2MCLimit = i;
1337     if (chi < 1000)
1338       limit = i;
1339     if (i == 150)
1340       fLastChi2MC = chi2;
1341   }
1342
1343   printf("limits %d %d\n", fLastChi2MCLimit, limit);
1344   fLastChi2MCLimit = limit;
1345
1346   printf("Difference (from MC) is %f for bin 2-150. Limit is %d.\n", fLastChi2MC, fLastChi2MCLimit);*/
1347
1348   if (!simple)
1349   {
1350     /*canvas1->cd(3);
1351
1352     diffMCUnfolded->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e(unfolded)");
1353     //diffMCUnfolded->GetYaxis()->SetRangeUser(-20, 20);
1354     diffMCUnfolded->GetXaxis()->SetRangeUser(0, 200);
1355     diffMCUnfolded->DrawCopy("HIST");
1356
1357     TH1F* fluctuation = new TH1F("fluctuation", "fluctuation", 20, -5, 5);
1358     for (Int_t i=20; i<=diffMCUnfolded->GetNbinsX(); ++i)
1359       fluctuation->Fill(diffMCUnfolded->GetBinContent(i));
1360
1361     //new TCanvas; fluctuation->DrawCopy();
1362     delete fluctuation;*/
1363
1364     /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1365     legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
1366     legend->AddEntry(fMultiplicityMC, "MC");
1367     legend->AddEntry(fMultiplicityESD, "ESD");
1368     legend->Draw();*/
1369
1370     canvas1->cd(4);
1371     //residual->GetYaxis()->SetRangeUser(-0.2, 0.2);
1372     residual->GetXaxis()->SetRangeUser(0, 200);
1373     residual->DrawCopy();
1374
1375     canvas1->cd(5);
1376
1377     TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
1378     ratio->Divide(mcHist);
1379     ratio->SetTitle("Ratio;true multiplicity;Unfolded / MC");
1380     ratio->GetYaxis()->SetRangeUser(0.5, 1.5);
1381     ratio->GetXaxis()->SetRangeUser(0, 200);
1382     ratio->SetStats(kFALSE);
1383     ratio->Draw("HIST");
1384
1385     Double_t ratioChi2 = 0;
1386     fRatioAverage = 0;
1387     fLastChi2MCLimit = 0;
1388     for (Int_t i=2; i<=150; ++i)
1389     {
1390       Float_t value = ratio->GetBinContent(i) - 1;
1391       if (value > 1e8)
1392         value = 1e8; //prevent arithmetic exception
1393       else if (value < -1e8)
1394         value = -1e8;
1395
1396       ratioChi2 += value * value;
1397       fRatioAverage += TMath::Abs(value);
1398
1399       if (ratioChi2 < 0.1)
1400         fLastChi2MCLimit = i;
1401     }
1402     fRatioAverage /= 149;
1403
1404     printf("Sum over (ratio-1)^2 (2..150) is %f; average of |ratio-1| is %f\n", ratioChi2, fRatioAverage);
1405     // TODO FAKE
1406     fLastChi2MC = ratioChi2;
1407
1408     // FFT of ratio
1409     canvas1->cd(6);
1410     const Int_t kFFT = 128;
1411     Double_t fftReal[kFFT];
1412     Double_t fftImag[kFFT];
1413
1414     for (Int_t i=0; i<kFFT; ++i)
1415     {
1416       fftReal[i] = ratio->GetBinContent(i+1+10);
1417       // test: ;-)
1418       //fftReal[i] = cos(TMath::Pi() * 5 * 2 * i / 128);
1419       fftImag[i] = 0;
1420     }
1421
1422     FFT(-1, TMath::Nint(log(kFFT) / log(2)), fftReal, fftImag);
1423
1424     TH1* fftResult = (TH1*) ratio->Clone("fftResult");
1425     fftResult->SetTitle("FFT;true multiplicity;coeff. (10...137)");
1426     TH1* fftResult2 = (TH1*) ratio->Clone("fftResult2");
1427     TH1* fftResult3 = (TH1*) ratio->Clone("fftResult3");
1428     fftResult->Reset();
1429     fftResult2->Reset();
1430     fftResult3->Reset();
1431     fftResult->GetYaxis()->UnZoom();
1432     fftResult2->GetYaxis()->UnZoom();
1433     fftResult3->GetYaxis()->UnZoom();
1434
1435     Printf("AFTER FFT");
1436     for (Int_t i=0; i<kFFT; ++i)
1437     {
1438       Printf("%d: %f", i, fftReal[i]);
1439       fftResult->SetBinContent(i+1, fftReal[i]);
1440       /*if (i != 0 && TMath::Abs(fftReal[i]) > 0.5)
1441       {
1442         Printf("Nulled %d", i);
1443         fftReal[i] = 0;
1444       }*/
1445     }
1446
1447     fftResult->SetLineColor(1);
1448     fftResult->DrawCopy();
1449
1450
1451     Printf("IMAG");
1452     for (Int_t i=0; i<kFFT; ++i)
1453     {
1454       Printf("%d: %f", i, fftImag[i]);
1455       fftResult2->SetBinContent(i+1, fftImag[i]);
1456
1457       fftResult3->SetBinContent(i+1, TMath::Sqrt(fftReal[i] * fftReal[i] + fftImag[i] * fftImag[i]));
1458     }
1459
1460     fftResult2->SetLineColor(2);
1461     fftResult2->DrawCopy("SAME");
1462
1463     fftResult3->SetLineColor(4);
1464     fftResult3->DrawCopy("SAME");
1465
1466     for (Int_t i=1; i<kFFT - 1; ++i)
1467     {
1468       if (TMath::Sqrt(fftReal[i] * fftReal[i] + fftImag[i] * fftImag[i]) > 3)
1469       {
1470         fftReal[i-1] = 0;
1471         fftReal[i] = 0;
1472         fftReal[i+1] = 0;
1473         fftImag[i-1] = 0;
1474         fftImag[i] = 0;
1475         fftImag[i+1] = 0;
1476         //fftReal[i] = (fftReal[i-1] + fftReal[i+1]) / 2;
1477         //fftImag[i]  = (fftImag[i-1] + fftImag[i+1]) / 2;
1478         Printf("Nulled %d to %f %f", i, fftReal[i], fftImag[i]);
1479       }
1480     }
1481
1482     FFT(1, TMath::Nint(log(kFFT) / log(2)), fftReal, fftImag);
1483
1484     TH1* fftResult4 = (TH1*) fftResult3->Clone("fftResult4");
1485     fftResult4->Reset();
1486
1487     Printf("AFTER BACK-TRAFO");
1488     for (Int_t i=0; i<kFFT; ++i)
1489     {
1490       Printf("%d: %f", i, fftReal[i]);
1491       fftResult4->SetBinContent(i+1+10, fftReal[i]);
1492     }
1493
1494     canvas1->cd(5);
1495     fftResult4->SetLineColor(4);
1496     fftResult4->DrawCopy("SAME");
1497
1498     // plot (MC - Unfolded) / error (MC)
1499     canvas1->cd(3);
1500
1501     TH1* diffMCUnfolded2 = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded2"));
1502     diffMCUnfolded2->Add(fMultiplicityESDCorrected[esdCorrId], -1);
1503
1504     Int_t ndfQual[kQualityRegions];
1505     for (Int_t region=0; region<kQualityRegions; ++region)
1506     {
1507       fQuality[region] = 0;
1508       ndfQual[region] = 0;
1509     }
1510
1511
1512     Double_t newChi2 = 0;
1513     Double_t newChi2_150 = 0;
1514     Int_t ndf = 0;
1515     for (Int_t i=1; i<=TMath::Min(diffMCUnfolded2->GetNbinsX(), fgMaxParams+1); ++i)
1516     {
1517       Double_t value = 0;
1518       if (proj->GetBinError(i) > 0)
1519       {
1520         value = diffMCUnfolded2->GetBinContent(i) / proj->GetBinError(i);
1521         newChi2 += value * value;
1522         if (i > 1 && i <= mcBinLimit)
1523           newChi2_150 += value * value;
1524         ++ndf;
1525
1526         for (Int_t region=0; region<kQualityRegions; ++region)
1527           if (diffMCUnfolded2->GetXaxis()->GetBinCenter(i) >= fgQualityRegionsB[region] - 0.1 && diffMCUnfolded2->GetXaxis()->GetBinCenter(i) <= fgQualityRegionsE[region] + 0.1) // 0.1 to avoid e.g. 3.9999 < 4 problem
1528           {
1529             fQuality[region] += TMath::Abs(value);
1530             ++ndfQual[region];
1531           }
1532       }
1533
1534       diffMCUnfolded2->SetBinContent(i, value);
1535     }
1536
1537     // normalize region to the number of entries
1538     for (Int_t region=0; region<kQualityRegions; ++region)
1539     {
1540       if (ndfQual[region] > 0)
1541         fQuality[region] /= ndfQual[region];
1542       Printf("Quality parameter %d (%d <= mult <= %d) is %f with %d df", region, fgQualityRegionsB[region], fgQualityRegionsE[region], fQuality[region], ndfQual[region]);
1543     }
1544
1545     if (mcBinLimit > 1)
1546     {
1547       // TODO another FAKE
1548       fLastChi2MC = newChi2_150 / (mcBinLimit - 1);
1549       Printf("Chi2 (2..%d) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", mcBinLimit, newChi2_150, mcBinLimit - 1, fLastChi2MC);
1550     }
1551     else
1552       fLastChi2MC = -1;
1553
1554     Printf("Chi2 (full range) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", newChi2, ndf, newChi2 / ndf);
1555
1556     diffMCUnfolded2->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e(MC)");
1557     //diffMCUnfolded->GetYaxis()->SetRangeUser(-20, 20);
1558     diffMCUnfolded2->GetXaxis()->SetRangeUser(0, 200);
1559     diffMCUnfolded2->DrawCopy("HIST");
1560   }
1561
1562   canvas1->SaveAs(Form("%s.gif", canvas1->GetName()));
1563 }
1564
1565 //____________________________________________________________________
1566 void AliMultiplicityCorrection::FFT(Int_t dir, Int_t m, Double_t *x, Double_t *y)
1567 {
1568 /*-------------------------------------------------------------------------
1569    This computes an in-place complex-to-complex FFT
1570    x and y are the real and imaginary arrays of 2^m points.
1571    dir =  1 gives forward transform
1572    dir = -1 gives reverse transform
1573
1574      Formula: forward
1575                   N-1
1576                   ---
1577               1   \          - j k 2 pi n / N
1578       X(n) = ---   >   x(k) e                    = forward transform
1579               N   /                                n=0..N-1
1580                   ---
1581                   k=0
1582
1583       Formula: reverse
1584                   N-1
1585                   ---
1586                   \          j k 2 pi n / N
1587       X(n) =       >   x(k) e                    = forward transform
1588                   /                                n=0..N-1
1589                   ---
1590                   k=0
1591 */
1592
1593    Long_t   nn, i, i1, j, k, i2, l, l1, l2;
1594    Double_t c1, c2, tx, ty, t1, t2, u1, u2, z;
1595
1596    /* Calculate the number of points */
1597    nn = 1;
1598    for (i = 0; i < m; i++)
1599        nn *= 2;
1600
1601    /* Do the bit reversal */
1602    i2 = nn >> 1;
1603    j = 0;
1604    for (i= 0; i < nn - 1; i++) {
1605        if (i < j) {
1606            tx = x[i];
1607            ty = y[i];
1608            x[i] = x[j];
1609            y[i] = y[j];
1610            x[j] = tx;
1611            y[j] = ty;
1612        }
1613        k = i2;
1614        while (k <= j) {
1615            j -= k;
1616            k >>= 1;
1617        }
1618        j += k;
1619    }
1620
1621    /* Compute the FFT */
1622    c1 = -1.0;
1623    c2 = 0.0;
1624    l2 = 1;
1625    for (l = 0; l < m; l++) {
1626        l1 = l2;
1627        l2 <<= 1;
1628        u1 = 1.0;
1629        u2 = 0.0;
1630        for (j = 0;j < l1; j++) {
1631            for (i = j; i < nn; i += l2) {
1632                i1 = i + l1;
1633                t1 = u1 * x[i1] - u2 * y[i1];
1634                t2 = u1 * y[i1] + u2 * x[i1];
1635                x[i1] = x[i] - t1;
1636                y[i1] = y[i] - t2;
1637                x[i] += t1;
1638                y[i] += t2;
1639            }
1640            z =  u1 * c1 - u2 * c2;
1641            u2 = u1 * c2 + u2 * c1;
1642            u1 = z;
1643        }
1644        c2 = TMath::Sqrt((1.0 - c1) / 2.0);
1645        if (dir == 1)
1646            c2 = -c2;
1647        c1 = TMath::Sqrt((1.0 + c1) / 2.0);
1648    }
1649
1650    /* Scaling for forward transform */
1651    if (dir == 1) {
1652        for (i=0;i<nn;i++) {
1653            x[i] /= (Double_t)nn;
1654            y[i] /= (Double_t)nn;
1655        }
1656    }
1657 }
1658
1659 //____________________________________________________________________
1660 void AliMultiplicityCorrection::GetComparisonResults(Float_t* mc, Int_t* mcLimit, Float_t* residuals, Float_t* ratioAverage)
1661 {
1662   // Returns the chi2 between the MC and the unfolded ESD as well as between the ESD and the folded unfolded ESD
1663   // These values are computed during DrawComparison, thus this function picks up the
1664   // last calculation
1665
1666   if (mc)
1667     *mc = fLastChi2MC;
1668   if (mcLimit)
1669     *mcLimit = fLastChi2MCLimit;
1670   if (residuals)
1671     *residuals = fLastChi2Residuals;
1672   if (ratioAverage)
1673     *ratioAverage = fRatioAverage;
1674 }
1675
1676 //____________________________________________________________________
1677 TH2F* AliMultiplicityCorrection::GetMultiplicityMC(Int_t i, EventType eventType)
1678 {
1679   //
1680   // returns the corresponding MC spectrum
1681   //
1682
1683   switch (eventType)
1684   {
1685     case kTrVtx : return fMultiplicityVtx[i]; break;
1686     case kMB : return fMultiplicityMB[i]; break;
1687     case kINEL : return fMultiplicityINEL[i]; break;
1688   }
1689
1690   return 0;
1691 }
1692
1693 //____________________________________________________________________
1694 TH1* AliMultiplicityCorrection::BayesianStatisticsEffect(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType,  Bool_t randomizeMeasured, Bool_t randomizeResponse, Float_t regPar, Int_t nIterations, TH1* compareTo)
1695 {
1696   //
1697   // evaluates the uncertainty that arises from the non-infinite statistics in the response matrix
1698   // the function unfolds the spectrum using the default response matrix and several modified ones
1699   // the modified ones are created by randomizing each cell using poisson statistics with the mean = bin value
1700   // these unfolded results are compared to the first result gained with the default response OR to the histogram given
1701   // in <compareTo> (optional)
1702   //
1703   // returns the error assigned to the measurement
1704   //
1705
1706   // initialize seed with current time
1707   gRandom->SetSeed(0);
1708
1709   const Int_t kErrorIterations = 150;
1710
1711   TH1* maxError = 0;
1712   TH1* firstResult = 0;
1713
1714   TH1** results = new TH1*[kErrorIterations];
1715
1716   for (Int_t n=0; n<kErrorIterations; ++n)
1717   {
1718     Printf("Iteration %d of %d...", n, kErrorIterations);
1719
1720     SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
1721
1722     TH1* measured = (TH1*) fCurrentESD->Clone("measured");
1723
1724     if (n > 0)
1725     {
1726       if (randomizeResponse)
1727       {
1728         // randomize response matrix
1729         for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1730           for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1731             fCurrentCorrelation->SetBinContent(i, j, gRandom->Poisson(fCurrentCorrelation->GetBinContent(i, j)));
1732       }
1733
1734       if (randomizeMeasured)
1735       {
1736         // randomize measured spectrum
1737         for (Int_t x=1; x<=measured->GetNbinsX(); x++) // mult. axis
1738         {
1739           Int_t randomValue = gRandom->Poisson(fCurrentESD->GetBinContent(x));
1740           measured->SetBinContent(x, randomValue);
1741           measured->SetBinError(x, TMath::Sqrt(randomValue));
1742         }
1743       }
1744     }
1745
1746     for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1747     {
1748       // with this it is normalized to 1
1749       Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1750
1751       // with this normalized to the given efficiency
1752       if (fCurrentEfficiency->GetBinContent(i) > 0)
1753         sum /= fCurrentEfficiency->GetBinContent(i);
1754       else
1755         sum = 0;
1756
1757       for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1758       {
1759         if (sum > 0)
1760         {
1761           fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1762           fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1763         }
1764         else
1765         {
1766           fCurrentCorrelation->SetBinContent(i, j, 0);
1767           fCurrentCorrelation->SetBinError(i, j, 0);
1768         }
1769       }
1770     }
1771
1772     TH1* result = 0;
1773     if (n == 0 && compareTo)
1774     {
1775       // in this case we just store the histogram we want to compare to
1776       result = (TH1*) compareTo->Clone("compareTo");
1777       result->Sumw2();
1778     }
1779     else
1780       result = UnfoldWithBayesian(measured, regPar, nIterations, 0);
1781
1782     if (!result)
1783       return 0;
1784
1785     // normalize
1786     result->Scale(1.0 / result->Integral());
1787
1788     if (n == 0)
1789     {
1790       firstResult = (TH1*) result->Clone("firstResult");
1791
1792       maxError = (TH1*) result->Clone("maxError");
1793       maxError->Reset();
1794     }
1795     else
1796     {
1797       // calculate ratio
1798       TH1* ratio = (TH1*) firstResult->Clone("ratio");
1799       ratio->Divide(result);
1800
1801       // find max. deviation
1802       for (Int_t x=1; x<=ratio->GetNbinsX(); x++)
1803         maxError->SetBinContent(x, TMath::Max(maxError->GetBinContent(x), TMath::Abs(1 - ratio->GetBinContent(x))));
1804
1805       delete ratio;
1806     }
1807
1808     results[n] = result;
1809   }
1810
1811   // find covariance matrix
1812   // results[n] is X_x
1813   // cov. matrix is M_xy = E ( (X_x - E(X_x)) * (X_y - E(X_y))), with E() = expectation value
1814
1815   Int_t nBins = results[0]->GetNbinsX();
1816   Float_t lowEdge = results[0]->GetXaxis()->GetBinLowEdge(1);
1817   Float_t upEdge = results[0]->GetXaxis()->GetBinUpEdge(nBins);
1818
1819   // find average, E(X)
1820   TProfile* average = new TProfile("average", "average", nBins, lowEdge, upEdge);
1821   for (Int_t n=1; n<kErrorIterations; ++n)
1822     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1823       average->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x));
1824   //new TCanvas; average->DrawClone();
1825
1826   // find cov. matrix
1827   TProfile2D* covMatrix = new TProfile2D("covMatrix", "covMatrix", nBins, lowEdge, upEdge, nBins, lowEdge, upEdge);
1828
1829   for (Int_t n=1; n<kErrorIterations; ++n)
1830     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1831       for (Int_t y=1; y<=results[n]->GetNbinsX(); y++)
1832       {
1833         // (X_x - E(X_x)) * (X_y - E(X_y)
1834         Float_t cov = (results[n]->GetBinContent(x) - average->GetBinContent(x)) * (results[n]->GetBinContent(y) - average->GetBinContent(y));
1835         covMatrix->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetXaxis()->GetBinCenter(y), cov);
1836       }
1837   new TCanvas; covMatrix->DrawClone("COLZ");
1838
1839   // fill 2D histogram that contains deviation from first
1840   TH2F* deviations = new TH2F("deviations", "deviations", nBins, lowEdge, upEdge, 1000, -0.01, 0.01);
1841   for (Int_t n=1; n<kErrorIterations; ++n)
1842     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1843       deviations->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x) - results[0]->GetBinContent(x));
1844   new TCanvas; deviations->DrawCopy("COLZ");
1845
1846   TH1* standardDeviation = (TH1*) maxError->Clone("standardDeviation");
1847   standardDeviation->Reset();
1848
1849   // get standard deviation "by hand"
1850   for (Int_t x=1; x<=nBins; x++)
1851   {
1852     if (results[0]->GetBinContent(x) > 0)
1853     {
1854       TH1* proj = deviations->ProjectionY("projRMS", x, x, "e");
1855       Float_t standardDev = proj->GetRMS(); // this is standard deviation in fact
1856       standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x));
1857     }
1858   }
1859
1860   // compare maxError to RMS of profile (variable name: average)
1861   // first: calculate rms in percent of value
1862   TH1* rmsError = (TH1*) maxError->Clone("rmsError");
1863   rmsError->Reset();
1864
1865   // enable error to be standard deviation (see http://root.cern.ch/root/html/TProfile.html#TProfile:SetErrorOption)
1866   average->SetErrorOption("s");
1867   for (Int_t x=1; x<=rmsError->GetNbinsX(); x++)
1868     if (average->GetBinContent(x) > 0)
1869       rmsError->SetBinContent(x, average->GetBinError(x) / average->GetBinContent(x));
1870
1871   // find maxError deviation from average (not from "first result"/mc as above)
1872   TH1* maxError2 = (TH1*) maxError->Clone("maxError2");
1873   maxError2->Reset();
1874   for (Int_t n=1; n<kErrorIterations; ++n)
1875     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1876       if (average->GetBinContent(x) > 0)
1877         maxError2->SetBinContent(x, TMath::Max(maxError2->GetBinContent(x), TMath::Abs((results[n]->GetBinContent(x) - average->GetBinContent(x)) / average->GetBinContent(x))));
1878
1879   new TCanvas; maxError2->DrawCopy(); rmsError->SetLineColor(2); rmsError->DrawCopy("SAME"); standardDeviation->SetLineColor(3); standardDeviation->DrawCopy("SAME");
1880
1881   // plot difference between average and MC/first result
1882   TH1* averageFirstRatio = (TH1*) results[0]->Clone("averageFirstRatio");
1883   averageFirstRatio->Reset();
1884   averageFirstRatio->Divide(results[0], average);
1885
1886   new TCanvas; results[0]->DrawCopy(); average->SetLineColor(2); average->DrawClone("SAME");
1887   new TCanvas; averageFirstRatio->DrawCopy();
1888
1889   // clean up
1890   for (Int_t n=0; n<kErrorIterations; ++n)
1891     delete results[n];
1892   delete[] results;
1893
1894   // fill into result histogram
1895   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1896
1897   for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
1898     fMultiplicityESDCorrected[correlationID]->SetBinContent(i, firstResult->GetBinContent(i));
1899
1900   for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
1901     fMultiplicityESDCorrected[correlationID]->SetBinError(i, maxError->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
1902
1903   return standardDeviation;
1904 }
1905
1906 //____________________________________________________________________
1907 void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar, Int_t nIterations, TH1* inputDist, Bool_t determineError)
1908 {
1909   //
1910   // correct spectrum using bayesian method
1911   //
1912
1913   // initialize seed with current time
1914   gRandom->SetSeed(0);
1915
1916   SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
1917
1918   // normalize correction for given nPart
1919   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1920   {
1921     // with this it is normalized to 1
1922     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1923
1924     // with this normalized to the given efficiency
1925     if (fCurrentEfficiency->GetBinContent(i) > 0)
1926       sum /= fCurrentEfficiency->GetBinContent(i);
1927     else
1928       sum = 0;
1929
1930     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1931     {
1932       if (sum > 0)
1933       {
1934         fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1935         fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1936       }
1937       else
1938       {
1939         fCurrentCorrelation->SetBinContent(i, j, 0);
1940         fCurrentCorrelation->SetBinError(i, j, 0);
1941       }
1942     }
1943   }
1944
1945   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1946
1947   TH1* result = UnfoldWithBayesian(fCurrentESD, regPar, nIterations, inputDist);
1948   if (!result)
1949     return;
1950
1951   for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
1952     fMultiplicityESDCorrected[correlationID]->SetBinContent(i, result->GetBinContent(i));
1953
1954   if (!determineError)
1955   {
1956     Printf("AliMultiplicityCorrection::ApplyBayesianMethod: WARNING: No errors calculated.");
1957     return;
1958   }
1959
1960   // evaluate errors, this is done by randomizing the measured spectrum following Poission statistics
1961   // this (new) measured spectrum is then unfolded and the different to the result from the "real" measured
1962   // spectrum calculated. This is performed N times and the maximum difference is taken as the systematic
1963   // error of the unfolding method itself.
1964
1965   TH1* maxError = (TH1*) result->Clone("maxError");
1966   maxError->Reset();
1967
1968   TH1* normalizedResult = (TH1*) result->Clone("normalizedResult");
1969   normalizedResult->Scale(1.0 / normalizedResult->Integral());
1970
1971   const Int_t kErrorIterations = 20;
1972
1973   printf("Spectrum unfolded. Determining error (%d iterations)...\n", kErrorIterations);
1974
1975   TH1* randomized = (TH1*) fCurrentESD->Clone("randomized");
1976   for (Int_t n=0; n<kErrorIterations; ++n)
1977   {
1978     // randomize the content of clone following a poisson with the mean = the value of that bin
1979     for (Int_t x=1; x<=randomized->GetNbinsX(); x++) // mult. axis
1980     {
1981       Int_t randomValue = gRandom->Poisson(fCurrentESD->GetBinContent(x));
1982       //printf("%d --> %d\n", fCurrentESD->GetBinContent(x), randomValue);
1983       randomized->SetBinContent(x, randomValue);
1984       randomized->SetBinError(x, TMath::Sqrt(randomValue));
1985     }
1986
1987     TH1* result2 = UnfoldWithBayesian(randomized, regPar, nIterations, inputDist);
1988     if (!result2)
1989       return;
1990
1991     result2->Scale(1.0 / result2->Integral());
1992
1993     // calculate ratio
1994     TH1* ratio = (TH1*) result2->Clone("ratio");
1995     ratio->Divide(normalizedResult);
1996
1997     //new TCanvas; ratio->DrawCopy("HIST");
1998
1999     // find max. deviation
2000     for (Int_t x=1; x<=ratio->GetNbinsX(); x++)
2001       maxError->SetBinContent(x, TMath::Max(maxError->GetBinContent(x), TMath::Abs(1 - ratio->GetBinContent(x))));
2002
2003     delete ratio;
2004     delete result2;
2005   }
2006
2007   for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
2008     fMultiplicityESDCorrected[correlationID]->SetBinError(i, maxError->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
2009
2010   delete maxError;
2011   delete normalizedResult;
2012 }
2013
2014 //____________________________________________________________________
2015 TH1* AliMultiplicityCorrection::UnfoldWithBayesian(TH1* measured, Float_t regPar, Int_t nIterations, TH1* inputDist)
2016 {
2017   //
2018   // unfolds a spectrum
2019   //
2020
2021   if (measured->Integral() <= 0)
2022   {
2023     Printf("AliMultiplicityCorrection::UnfoldWithBayesian: ERROR: The measured spectrum is empty");
2024     return 0;
2025   }
2026
2027   const Int_t kStartBin = 0;
2028
2029   const Int_t kMaxM = fgMaxInput; //<= fCurrentCorrelation->GetNbinsY(); // max measured axis
2030   const Int_t kMaxT = fgMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis
2031
2032   // store information in arrays, to increase processing speed (~ factor 5)
2033   Double_t measuredCopy[kMaxM];
2034   Double_t prior[kMaxT];
2035   Double_t result[kMaxT];
2036   Double_t efficiency[kMaxT];
2037
2038   Double_t response[kMaxT][kMaxM];
2039   Double_t inverseResponse[kMaxT][kMaxM];
2040
2041   // for normalization
2042   Float_t measuredIntegral = measured->Integral();
2043   for (Int_t m=0; m<kMaxM; m++)
2044   {
2045     measuredCopy[m] = measured->GetBinContent(m+1) / measuredIntegral;
2046
2047     for (Int_t t=0; t<kMaxT; t++)
2048     {
2049       response[t][m] = fCurrentCorrelation->GetBinContent(t+1, m+1);
2050       inverseResponse[t][m] = 0;
2051     }
2052   }
2053
2054   for (Int_t t=0; t<kMaxT; t++)
2055   {
2056     efficiency[t] = fCurrentEfficiency->GetBinContent(t+1);
2057     prior[t] = measuredCopy[t];
2058     result[t] = 0;
2059   }
2060
2061   // pick prior distribution
2062   if (inputDist)
2063   {
2064     printf("Using different starting conditions...\n");
2065     // for normalization
2066     Float_t inputDistIntegral = inputDist->Integral();
2067     for (Int_t i=0; i<kMaxT; i++)
2068       prior[i] = inputDist->GetBinContent(i+1) / inputDistIntegral;
2069   }
2070
2071   //TH1F* convergence = new TH1F("convergence", "convergence", 50, 0.5, 50.5);
2072
2073   // unfold...
2074   for (Int_t i=0; i<nIterations; i++)
2075   {
2076     //printf(" iteration %i \n", i);
2077
2078     // calculate IR from Beyes theorem
2079     // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
2080
2081     for (Int_t m=0; m<kMaxM; m++)
2082     {
2083       Float_t norm = 0;
2084       for (Int_t t = kStartBin; t<kMaxT; t++)
2085         norm += response[t][m] * prior[t];
2086
2087       if (norm > 0)
2088       {
2089         for (Int_t t = kStartBin; t<kMaxT; t++)
2090           inverseResponse[t][m] = response[t][m] * prior[t] / norm;
2091       }
2092       else
2093       {
2094         for (Int_t t = kStartBin; t<kMaxT; t++)
2095           inverseResponse[t][m] = 0;
2096       }
2097     }
2098
2099     for (Int_t t = kStartBin; t<kMaxT; t++)
2100     {
2101       Float_t value = 0;
2102       for (Int_t m=0; m<kMaxM; m++)
2103         value += inverseResponse[t][m] * measuredCopy[m];
2104
2105       if (efficiency[t] > 0)
2106         result[t] = value / efficiency[t];
2107       else
2108         result[t] = 0;
2109     }
2110
2111     // regularization (simple smoothing)
2112     for (Int_t t=kStartBin; t<kMaxT; t++)
2113     {
2114       Float_t newValue = 0;
2115       // 0 bin excluded from smoothing
2116       if (t > kStartBin+1 && t<kMaxT-1)
2117       {
2118         Float_t average = (result[t-1] + result[t] + result[t+1]) / 3;
2119
2120         // weight the average with the regularization parameter
2121         newValue = (1 - regPar) * result[t] + regPar * average;
2122       }
2123       else
2124         newValue = result[t];
2125
2126       prior[t] = newValue;
2127     }
2128
2129     // calculate chi2 (change from last iteration)
2130     //Double_t chi2 = 0;
2131     /*for (Int_t t=kStartBin; t<=hTrueSmoothed->GetNbinsX(); t++)
2132     {
2133       Float_t newValue = hTrueSmoothed->GetBinContent(t);
2134       //Float_t diff = hPrior->GetBinContent(t) - newValue;
2135       //chi2 += (Double_t) diff * diff;
2136
2137       hPrior->SetBinContent(t, newValue);
2138     }
2139     //printf("Chi2 of %d iteration = %.10f\n", i, chi2);
2140     //convergence->Fill(i+1, chi2);
2141
2142     //if (i % 5 == 0)
2143     //  DrawComparison(Form("Bayesian_%d", i), mcTarget, correlationID, kTRUE, eventType);
2144
2145     delete hTrueSmoothed;*/
2146   } // end of iterations
2147
2148   //new TCanvas;
2149   //convergence->DrawCopy();
2150   //gPad->SetLogy();
2151
2152   //delete convergence;
2153
2154   // TODO this does not work when the number of bins is differnt
2155   TH1* resultHist = (TH1*) measured->Clone("resultHist");
2156   resultHist->Reset();
2157   for (Int_t t=0; t<kMaxT; t++)
2158     resultHist->SetBinContent(t+1, result[t]);
2159
2160   return resultHist;
2161
2162   // ********
2163   // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995
2164
2165   /*printf("Calculating covariance matrix. This may take some time...\n");
2166
2167   // TODO check if this is the right one...
2168   TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
2169
2170   Int_t xBins = hInverseResponseBayes->GetNbinsX();
2171   Int_t yBins = hInverseResponseBayes->GetNbinsY();
2172
2173   // calculate "unfolding matrix" Mij
2174   Float_t matrixM[251][251];
2175   for (Int_t i=1; i<=xBins; i++)
2176   {
2177     for (Int_t j=1; j<=yBins; j++)
2178     {
2179       if (fCurrentEfficiency->GetBinContent(i) > 0)
2180         matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i);
2181       else
2182         matrixM[i-1][j-1] = 0;
2183     }
2184   }
2185
2186   Float_t* vectorn = new Float_t[yBins];
2187   for (Int_t j=1; j<=yBins; j++)
2188     vectorn[j-1] = fCurrentESD->GetBinContent(j);
2189
2190   // first part of covariance matrix, depends on input distribution n(E)
2191   Float_t cov1[251][251];
2192
2193   Float_t nEvents = fCurrentESD->Integral(); // N
2194
2195   xBins = 20;
2196   yBins = 20;
2197
2198   for (Int_t k=0; k<xBins; k++)
2199   {
2200     printf("In Cov1: %d\n", k);
2201     for (Int_t l=0; l<yBins; l++)
2202     {
2203       cov1[k][l] = 0;
2204
2205       // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N)
2206       for (Int_t j=0; j<yBins; j++)
2207         cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j]
2208           * (1.0 - vectorn[j] / nEvents);
2209
2210       // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N
2211       for (Int_t i=0; i<yBins; i++)
2212         for (Int_t j=0; j<yBins; j++)
2213         {
2214           if (i == j)
2215             continue;
2216           cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i]
2217             * vectorn[j] / nEvents;
2218          }
2219     }
2220   }
2221
2222   printf("Cov1 finished\n");
2223
2224   TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov");
2225   cov->Reset();
2226
2227   for (Int_t i=1; i<=xBins; i++)
2228     for (Int_t j=1; j<=yBins; j++)
2229       cov->SetBinContent(i, j, cov1[i-1][j-1]);
2230
2231   new TCanvas;
2232   cov->Draw("COLZ");
2233
2234   // second part of covariance matrix, depends on response matrix
2235   Float_t cov2[251][251];
2236
2237   // Cov[P(Er|Cu), P(Es|Cu)] term
2238   Float_t covTerm[100][100][100];
2239   for (Int_t r=0; r<yBins; r++)
2240     for (Int_t u=0; u<xBins; u++)
2241       for (Int_t s=0; s<yBins; s++)
2242       {
2243         if (r == s)
2244           covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
2245             * (1.0 - hResponse->GetBinContent(u+1, r+1));
2246         else
2247           covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
2248             * hResponse->GetBinContent(u+1, s+1);
2249       }
2250
2251   for (Int_t k=0; k<xBins; k++)
2252     for (Int_t l=0; l<yBins; l++)
2253     {
2254       cov2[k][l] = 0;
2255       printf("In Cov2: %d %d\n", k, l);
2256       for (Int_t i=0; i<yBins; i++)
2257         for (Int_t j=0; j<yBins; j++)
2258         {
2259           //printf("In Cov2: %d %d %d %d\n", k, l, i, j);
2260           // calculate Cov(Mki, Mlj) = sum{ru},{su} ...
2261           Float_t tmpCov = 0;
2262           for (Int_t r=0; r<yBins; r++)
2263             for (Int_t u=0; u<xBins; u++)
2264               for (Int_t s=0; s<yBins; s++)
2265               {
2266                 if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0
2267                   || hResponse->GetBinContent(u+1, i+1) == 0)
2268                   continue;
2269
2270                 tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u)
2271                   * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u)
2272                   * covTerm[r][u][s];
2273               }
2274
2275           cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov;
2276         }
2277     }
2278
2279   printf("Cov2 finished\n");
2280
2281   for (Int_t i=1; i<=xBins; i++)
2282     for (Int_t j=1; j<=yBins; j++)
2283       cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]);
2284
2285   new TCanvas;
2286   cov->Draw("COLZ");*/
2287 }
2288
2289 //____________________________________________________________________
2290 Float_t AliMultiplicityCorrection::BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, Int_t k, Int_t i, Int_t r, Int_t u)
2291 {
2292   //
2293   // helper function for the covariance matrix of the bayesian method
2294   //
2295
2296   Float_t result = 0;
2297
2298   if (k == u && r == i)
2299     result += 1.0 / hResponse->GetBinContent(u+1, r+1);
2300
2301   if (k == u)
2302     result -= 1.0 / fCurrentEfficiency->GetBinContent(u+1);
2303
2304   if (r == i)
2305     result -= matrixM[u][i] * fCurrentEfficiency->GetBinContent(u+1) / hResponse->GetBinContent(u+1, i+1);
2306
2307   result *= matrixM[k][i];
2308
2309   return result;
2310 }
2311
2312 //____________________________________________________________________
2313 void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
2314 {
2315   //
2316   // correct spectrum using bayesian method
2317   //
2318
2319   Float_t regPar = 0;
2320
2321   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
2322   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
2323
2324   SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
2325   //normalize ESD
2326   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
2327
2328   // TODO should be taken from correlation map
2329   //TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
2330
2331   // normalize correction for given nPart
2332   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2333   {
2334     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
2335     //Double_t sum = sumHist->GetBinContent(i);
2336     if (sum <= 0)
2337       continue;
2338     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
2339     {
2340       // npart sum to 1
2341       fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
2342       fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
2343     }
2344   }
2345
2346   new TCanvas;
2347   fCurrentCorrelation->Draw("COLZ");
2348
2349   // FAKE
2350   fCurrentEfficiency = ((TH2*) fCurrentCorrelation)->ProjectionX("eff");
2351
2352   // pick prior distribution
2353   TH1F* hPrior = (TH1F*)fCurrentESD->Clone("prior");
2354   Float_t norm = 1; //hPrior->Integral();
2355   for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
2356     hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm);
2357
2358   // zero distribution
2359   TH1F* zero =  (TH1F*)hPrior->Clone("zero");
2360
2361   // define temp hist
2362   TH1F* hTemp = (TH1F*)fCurrentESD->Clone("temp");
2363   hTemp->Reset();
2364
2365   // just a shortcut
2366   TH2F* hResponse = (TH2F*) fCurrentCorrelation;
2367
2368   // unfold...
2369   Int_t iterations = 25;
2370   for (Int_t i=0; i<iterations; i++)
2371   {
2372     //printf(" iteration %i \n", i);
2373
2374     for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
2375     {
2376       Float_t value = 0;
2377       for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
2378         value += hResponse->GetBinContent(t, m) * hPrior->GetBinContent(t);
2379       hTemp->SetBinContent(m, value);
2380       //printf("%d %f %f %f\n", m, zero->GetBinContent(m), hPrior->GetBinContent(m), value);
2381     }
2382
2383     // regularization (simple smoothing)
2384     TH1F* hTrueSmoothed = (TH1F*) hTemp->Clone("truesmoothed");
2385
2386     for (Int_t t=2; t<hTrueSmoothed->GetNbinsX(); t++)
2387     {
2388       Float_t average = (hTemp->GetBinContent(t-1) / hTemp->GetBinWidth(t-1)
2389                          + hTemp->GetBinContent(t) / hTemp->GetBinWidth(t)
2390                          + hTemp->GetBinContent(t+1) / hTemp->GetBinWidth(t+1)) / 3.;
2391       average *= hTrueSmoothed->GetBinWidth(t);
2392
2393       // weight the average with the regularization parameter
2394       hTrueSmoothed->SetBinContent(t, (1-regPar) * hTemp->GetBinContent(t) + regPar * average);
2395     }
2396
2397     for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
2398       hTemp->SetBinContent(m, zero->GetBinContent(m) + hPrior->GetBinContent(m) - hTrueSmoothed->GetBinContent(m));
2399
2400     // fill guess
2401     for (Int_t t=1; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++)
2402     {
2403       fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTemp->GetBinContent(t));
2404       fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t)); // TODO
2405
2406       //printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
2407     }
2408
2409
2410     // calculate chi2 (change from last iteration)
2411     Double_t chi2 = 0;
2412
2413     // use smoothed true (normalized) as new prior
2414     Float_t norm = 1; //hTrueSmoothed->Integral();
2415
2416     for (Int_t t=1; t<hTrueSmoothed->GetNbinsX(); t++)
2417     {
2418       Float_t newValue = hTemp->GetBinContent(t)/norm;
2419       Float_t diff = hPrior->GetBinContent(t) - newValue;
2420       chi2 += (Double_t) diff * diff;
2421
2422       hPrior->SetBinContent(t, newValue);
2423     }
2424
2425     printf("Chi2 of %d iteration = %.10f\n", i, chi2);
2426
2427     //if (i % 5 == 0)
2428       DrawComparison(Form("Laszlo_%d", i), inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
2429
2430     delete hTrueSmoothed;
2431   } // end of iterations
2432
2433   DrawComparison("Laszlo", inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
2434 }
2435
2436 //____________________________________________________________________
2437 void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace)
2438 {
2439   //
2440   // correct spectrum using a simple Gaussian approach, that is model-dependent
2441   //
2442
2443   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
2444   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
2445
2446   SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE);
2447   //normalize ESD
2448   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
2449
2450   NormalizeToBinWidth((TH2*) fCurrentCorrelation);
2451
2452   TH1D* correction = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianMean"));
2453   correction->SetTitle("GaussianMean");
2454
2455   TH1D* correctionWidth = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianWidth"));
2456   correction->SetTitle("GaussianWidth");
2457
2458   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2459   {
2460     TH1D* proj = (dynamic_cast<TH2*> (fCurrentCorrelation))->ProjectionX("_px", i, i+1);
2461     proj->Fit("gaus", "0Q");
2462     correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
2463     correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
2464
2465     continue;
2466
2467     // draw for debugging
2468     new TCanvas;
2469     proj->DrawCopy();
2470     proj->GetFunction("gaus")->DrawCopy("SAME");
2471   }
2472
2473   TH1* target = fMultiplicityESDCorrected[correlationID];
2474
2475   Int_t nBins = target->GetNbinsX()*10+1;
2476   Float_t* binning = new Float_t[nBins];
2477   for (Int_t i=1; i<=target->GetNbinsX(); ++i)
2478     for (Int_t j=0; j<10; ++j)
2479       binning[(i-1)*10 + j] = target->GetXaxis()->GetBinLowEdge(i) + target->GetXaxis()->GetBinWidth(i) / 10 * j;
2480
2481   binning[nBins-1] = target->GetXaxis()->GetBinUpEdge(target->GetNbinsX());
2482
2483   TH1F* fineBinned = new TH1F("targetFineBinned", "targetFineBinned", nBins-1, binning);
2484
2485   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2486   {
2487     Float_t mean = correction->GetBinContent(i);
2488     Float_t width = correctionWidth->GetBinContent(i);
2489
2490     Int_t fillBegin = fineBinned->FindBin(mean - width * 5);
2491     Int_t fillEnd   = fineBinned->FindBin(mean + width * 5);
2492     //printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
2493
2494     for (Int_t j=fillBegin; j <= fillEnd; ++j)
2495     {
2496       fineBinned->AddBinContent(j, TMath::Gaus(fineBinned->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fCurrentESD->GetBinContent(i));
2497     }
2498   }
2499
2500   for (Int_t i=1; i<=target->GetNbinsX(); ++i)
2501   {
2502     Float_t sum = 0;
2503     for (Int_t j=1; j<=10; ++j)
2504       sum += fineBinned->GetBinContent((i-1)*10 + j);
2505     target->SetBinContent(i, sum / 10);
2506   }
2507
2508   delete[] binning;
2509
2510   DrawComparison("Gaussian", inputRange, fullPhaseSpace, kFALSE, GetMultiplicityMC(mcTarget, kTrVtx)->ProjectionY());
2511 }
2512
2513 //____________________________________________________________________
2514 TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap)
2515 {
2516   // runs the distribution given in inputMC through the response matrix identified by
2517   // correlationMap and produces a measured distribution
2518   // although it is a TH2F the vertex axis is not used at the moment and all entries are filled in mid-vertex
2519   // if normalized is set, inputMC is expected to be normalized to the bin width
2520
2521   if (!inputMC)
2522     return 0;
2523
2524   if (correlationMap < 0 || correlationMap >= kCorrHists)
2525     return 0;
2526
2527   // empty under/overflow bins in x, otherwise Project3D takes them into account
2528   TH3* hist = fCorrelation[correlationMap];
2529   for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
2530   {
2531     for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
2532     {
2533       hist->SetBinContent(0, y, z, 0);
2534       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
2535     }
2536   }
2537
2538   TH2* corr = (TH2*) hist->Project3D("zy");
2539   //corr->Rebin2D(2, 1);
2540   corr->Sumw2();
2541
2542   // normalize correction for given nPart
2543   for (Int_t i=1; i<=corr->GetNbinsX(); ++i)
2544   {
2545     Double_t sum = corr->Integral(i, i, 1, corr->GetNbinsY());
2546     if (sum <= 0)
2547       continue;
2548
2549     for (Int_t j=1; j<=corr->GetNbinsY(); ++j)
2550     {
2551       // npart sum to 1
2552       corr->SetBinContent(i, j, corr->GetBinContent(i, j) / sum);
2553       corr->SetBinError(i, j, corr->GetBinError(i, j) / sum);
2554     }
2555   }
2556
2557   TH2F* target = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("%s_measured", inputMC->GetName())));
2558   target->Reset();
2559
2560   for (Int_t meas=1; meas<=corr->GetNbinsY(); ++meas)
2561   {
2562     Float_t measured = 0;
2563     Float_t error = 0;
2564
2565     for (Int_t gen=1; gen<=corr->GetNbinsX(); ++gen)
2566     {
2567       Int_t mcGenBin = inputMC->GetXaxis()->FindBin(corr->GetXaxis()->GetBinCenter(gen));
2568
2569       measured += inputMC->GetBinContent(mcGenBin) * corr->GetBinContent(gen, meas);
2570       error += inputMC->GetBinError(mcGenBin) * corr->GetBinContent(gen, meas);
2571     }
2572
2573     //printf("%f +- %f ; %f +- %f \n", inputMC->GetBinContent(meas), inputMC->GetBinError(meas), measured, error);
2574
2575     target->SetBinContent(1 + target->GetNbinsX() / 2, meas, measured);
2576     target->SetBinError(1 + target->GetNbinsX() / 2, meas, error);
2577   }
2578
2579   return target;
2580 }
2581
2582 //____________________________________________________________________
2583 void AliMultiplicityCorrection::SetGenMeasFromFunc(TF1* inputMC, Int_t id)
2584 {
2585   // uses the given function to fill the input MC histogram and generates from that
2586   // the measured histogram by applying the response matrix
2587   // this can be used to evaluate if the methods work indepedently of the input
2588   // distribution
2589   // WARNING does not respect the vertex distribution, just fills central vertex bin
2590
2591   if (!inputMC)
2592     return;
2593
2594   if (id < 0 || id >= kESDHists)
2595     return;
2596
2597   TH2F* mc = fMultiplicityVtx[id];
2598
2599   mc->Reset();
2600
2601   for (Int_t i=1; i<=mc->GetNbinsY(); ++i)
2602   {
2603     mc->SetBinContent(mc->GetNbinsX() / 2 + 1, i, inputMC->Eval(mc->GetYaxis()->GetBinCenter(i)) * mc->GetYaxis()->GetBinWidth(i));
2604     mc->SetBinError(mc->GetNbinsX() / 2 + 1, i, 0);
2605   }
2606
2607   //new TCanvas;
2608   //mc->Draw("COLZ");
2609
2610   TH1* proj = mc->ProjectionY();
2611
2612   TString tmp(fMultiplicityESD[id]->GetName());
2613   delete fMultiplicityESD[id];
2614   fMultiplicityESD[id] = CalculateMultiplicityESD(proj, id);
2615   fMultiplicityESD[id]->SetName(tmp);
2616 }