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