some fixes/improvements while moving to new ESD format
[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::NormalizeToBinWidth(TH1* hist)
1229 {
1230   //
1231   // normalizes a 1-d histogram to its bin width
1232   //
1233
1234   for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
1235   {
1236     hist->SetBinContent(i, hist->GetBinContent(i) / hist->GetBinWidth(i));
1237     hist->SetBinError(i, hist->GetBinError(i) / hist->GetBinWidth(i));
1238   }
1239 }
1240
1241 //____________________________________________________________________
1242 void AliMultiplicityCorrection::NormalizeToBinWidth(TH2* hist)
1243 {
1244   //
1245   // normalizes a 2-d histogram to its bin width (x width * y width)
1246   //
1247
1248   for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
1249     for (Int_t j=1; j<=hist->GetNbinsY(); ++j)
1250     {
1251       Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j);
1252       hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor);
1253       hist->SetBinError(i, j, hist->GetBinError(i, j) / factor);
1254     }
1255 }
1256
1257 //____________________________________________________________________
1258 void AliMultiplicityCorrection::DrawHistograms()
1259 {
1260   //
1261   // draws the histograms of this class
1262   //
1263
1264   printf("ESD:\n");
1265
1266   TCanvas* canvas1 = new TCanvas("fMultiplicityESD", "fMultiplicityESD", 900, 600);
1267   canvas1->Divide(3, 2);
1268   for (Int_t i = 0; i < kESDHists; ++i)
1269   {
1270     canvas1->cd(i+1);
1271     fMultiplicityESD[i]->DrawCopy("COLZ");
1272     printf("%d --> %f\n", i, (Float_t) fMultiplicityESD[i]->ProjectionY()->GetMean());
1273   }
1274
1275   printf("Vtx:\n");
1276
1277   TCanvas* canvas2 = new TCanvas("fMultiplicityMC", "fMultiplicityMC", 900, 600);
1278   canvas2->Divide(3, 2);
1279   for (Int_t i = 0; i < kMCHists; ++i)
1280   {
1281     canvas2->cd(i+1);
1282     fMultiplicityVtx[i]->DrawCopy("COLZ");
1283     printf("%d --> %f\n", i, fMultiplicityVtx[i]->ProjectionY()->GetMean());
1284   }
1285
1286   TCanvas* canvas3 = new TCanvas("fCorrelation", "fCorrelation", 900, 900);
1287   canvas3->Divide(3, 3);
1288   for (Int_t i = 0; i < kCorrHists; ++i)
1289   {
1290     canvas3->cd(i+1);
1291     TH3* hist = dynamic_cast<TH3*> (fCorrelation[i]->Clone());
1292     for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
1293     {
1294       for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
1295       {
1296         hist->SetBinContent(0, y, z, 0);
1297         hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
1298       }
1299     }
1300     TH1* proj = hist->Project3D("zy");
1301     proj->DrawCopy("COLZ");
1302   }
1303 }
1304
1305 //____________________________________________________________________
1306 void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t normalizeESD, TH1* mcHist, Bool_t simple)
1307 {
1308   //mcHist->Rebin(2);
1309
1310   Int_t esdCorrId = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1311
1312   TString tmpStr;
1313   tmpStr.Form("%s_DrawComparison_%d", name, esdCorrId);
1314
1315   if (fMultiplicityESDCorrected[esdCorrId]->Integral() == 0)
1316   {
1317     printf("ERROR. Unfolded histogram is empty\n");
1318     return;
1319   }
1320
1321   //regain measured distribution used for unfolding, because the bins at high mult. were modified in SetupCurrentHists
1322   fCurrentESD = fMultiplicityESD[esdCorrId]->ProjectionY();
1323   fCurrentESD->Sumw2();
1324   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
1325
1326   // normalize unfolded result to 1
1327   fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral());
1328
1329   //fCurrentESD->Scale(mcHist->Integral(2, 200));
1330
1331   //new TCanvas;
1332   /*TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
1333   ratio->Divide(mcHist);
1334   ratio->Draw("HIST");
1335   ratio->Fit("pol0", "W0", "", 20, 120);
1336   Float_t scalingFactor = ratio->GetFunction("pol0")->GetParameter(0);
1337   delete ratio;
1338   mcHist->Scale(scalingFactor);*/
1339
1340   // find bin with <= 50 entries. this is used as limit for the chi2 calculation
1341   Int_t mcBinLimit = 0;
1342   for (Int_t i=20; i<mcHist->GetNbinsX(); ++i)
1343   {
1344     if (mcHist->GetBinContent(i) > 50)
1345     {
1346       mcBinLimit = i;
1347     }
1348     else
1349       break;
1350   }
1351   Printf("AliMultiplicityCorrection::DrawComparison: MC bin limit is %d", mcBinLimit);
1352
1353   // scale to 1
1354   mcHist->Sumw2();
1355   if (mcHist->Integral() > 0)
1356     mcHist->Scale(1.0 / mcHist->Integral());
1357
1358   // calculate residual
1359
1360   // for that we convolute the response matrix with the gathered result
1361   // if normalizeESD is not set, the histogram is already normalized, this needs to be passed to CalculateMultiplicityESD
1362   TH1* tmpESDEfficiencyRecorrected = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("tmpESDEfficiencyRecorrected");
1363   tmpESDEfficiencyRecorrected->Multiply(fCurrentEfficiency);
1364   TH2* convoluted = CalculateMultiplicityESD(tmpESDEfficiencyRecorrected, esdCorrId);
1365   TH1* convolutedProj = convoluted->ProjectionY("convolutedProj", -1, -1, "e");
1366   if (convolutedProj->Integral() > 0)
1367   {
1368     convolutedProj->Scale(1.0 / convolutedProj->Integral());
1369   }
1370   else
1371     printf("ERROR: convolutedProj is empty. Something went wrong calculating the convoluted histogram.\n");
1372
1373   TH1* residual = (TH1*) convolutedProj->Clone("residual");
1374   residual->SetTitle("Residuals;Ntracks;(folded unfolded measured - measured) / e");
1375
1376   residual->Add(fCurrentESD, -1);
1377   //residual->Divide(residual, fCurrentESD, 1, 1, "B");
1378
1379   TH1* residualHist = new TH1F("residualHist", "residualHist", 51, -5, 5);
1380
1381   // TODO fix errors
1382   Float_t chi2 = 0;
1383   for (Int_t i=1; i<=residual->GetNbinsX(); ++i)
1384   {
1385     if (fCurrentESD->GetBinError(i) > 0)
1386     {
1387       Float_t value = residual->GetBinContent(i) / fCurrentESD->GetBinError(i);
1388       if (i > 1 && i <= fLastBinLimit)
1389         chi2 += value * value;
1390       Printf("%d --> %f (%f)", i, value * value, chi2);
1391       residual->SetBinContent(i, value);
1392       residualHist->Fill(value);
1393     }
1394     else
1395     {
1396       //printf("Residual bin %d set to 0\n", i);
1397       residual->SetBinContent(i, 0);
1398     }
1399     convolutedProj->SetBinError(i, 0);
1400     residual->SetBinError(i, 0);
1401   }
1402   fLastChi2Residuals = chi2;
1403
1404   new TCanvas; residualHist->DrawCopy();
1405
1406   //residualHist->Fit("gaus", "N");
1407   //delete residualHist;
1408
1409   printf("Difference (Residuals) is %f for bin 2-%d\n", fLastChi2Residuals, fLastBinLimit);
1410
1411   TCanvas* canvas1 = 0;
1412   if (simple)
1413   {
1414     canvas1 = new TCanvas(tmpStr, tmpStr, 900, 400);
1415     canvas1->Divide(2, 1);
1416   }
1417   else
1418   {
1419     canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 1200);
1420     canvas1->Divide(2, 3);
1421   }
1422
1423   canvas1->cd(1);
1424   TH1* proj = (TH1*) mcHist->Clone("proj");
1425
1426   proj->GetXaxis()->SetRangeUser(0, 200);
1427   proj->SetTitle(";true multiplicity;Entries");
1428   proj->SetStats(kFALSE);
1429
1430   fMultiplicityESDCorrected[esdCorrId]->SetLineColor(2);
1431
1432   if (proj->GetEntries() > 0) {
1433     proj->DrawCopy("HIST");
1434     fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAME HIST E");
1435   }
1436   else
1437     fMultiplicityESDCorrected[esdCorrId]->DrawCopy("HIST E");
1438
1439   gPad->SetLogy();
1440
1441   TLegend* legend = new TLegend(0.3, 0.8, 0.93, 0.93);
1442   legend->AddEntry(proj, "true distribution");
1443   legend->AddEntry(fMultiplicityESDCorrected[esdCorrId], "unfolded distribution");
1444   legend->SetFillColor(0);
1445   legend->Draw();
1446   // unfortunately does not work. maybe a bug? --> legend->SetTextSizePixels(14);
1447
1448   canvas1->cd(2);
1449
1450   gPad->SetLogy();
1451   fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
1452   //fCurrentESD->SetLineColor(2);
1453   fCurrentESD->SetTitle(";measured multiplicity;Entries");
1454   fCurrentESD->SetStats(kFALSE);
1455   fCurrentESD->DrawCopy("HIST E");
1456
1457   convolutedProj->SetLineColor(2);
1458   //proj2->SetMarkerColor(2);
1459   //proj2->SetMarkerStyle(5);
1460   convolutedProj->DrawCopy("HIST SAME");
1461
1462   legend = new TLegend(0.3, 0.8, 0.93, 0.93);
1463   legend->AddEntry(fCurrentESD, "measured distribution");
1464   legend->AddEntry(convolutedProj, "R #otimes unfolded distribution");
1465   legend->SetFillColor(0);
1466   legend->Draw();
1467
1468   //TH1* diffMCUnfolded = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded"));
1469   //diffMCUnfolded->Add(fMultiplicityESDCorrected[esdCorrId], -1);
1470
1471   /*Float_t chi2 = 0;
1472   Float_t chi = 0;
1473   fLastChi2MCLimit = 0;
1474   Int_t limit = 0;
1475   for (Int_t i=2; i<=200; ++i)
1476   {
1477     if (proj->GetBinContent(i) != 0)
1478     {
1479       Float_t value = (proj->GetBinContent(i) - fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i)) / proj->GetBinContent(i);
1480       chi2 += value * value;
1481       chi += TMath::Abs(value);
1482
1483       //printf("%d %f\n", i, chi);
1484
1485       if (chi2 < 0.2)
1486         fLastChi2MCLimit = i;
1487
1488       if (chi < 3)
1489         limit = i;
1490
1491     }
1492   }*/
1493
1494   /*chi2 = 0;
1495   Float_t chi = 0;
1496   Int_t limit = 0;
1497   for (Int_t i=1; i<=diffMCUnfolded->GetNbinsX(); ++i)
1498   {
1499     if (fMultiplicityESDCorrected[esdCorrId]->GetBinError(i) > 0)
1500     {
1501       Double_t value = diffMCUnfolded->GetBinContent(i) / fMultiplicityESDCorrected[esdCorrId]->GetBinError(i);
1502       if (value > 1e8)
1503         value = 1e8; //prevent arithmetic exception
1504       else if (value < -1e8)
1505         value = -1e8;
1506       if (i > 1)
1507       {
1508         chi2 += value * value;
1509         chi += TMath::Abs(value);
1510       }
1511       diffMCUnfolded->SetBinContent(i, value);
1512     }
1513     else
1514     {
1515       //printf("diffMCUnfolded bin %d set to 0\n", i);
1516       diffMCUnfolded->SetBinContent(i, 0);
1517     }
1518     if (chi2 < 1000)
1519       fLastChi2MCLimit = i;
1520     if (chi < 1000)
1521       limit = i;
1522     if (i == 150)
1523       fLastChi2MC = chi2;
1524   }
1525
1526   printf("limits %d %d\n", fLastChi2MCLimit, limit);
1527   fLastChi2MCLimit = limit;
1528
1529   printf("Difference (from MC) is %f for bin 2-150. Limit is %d.\n", fLastChi2MC, fLastChi2MCLimit);*/
1530
1531   if (!simple)
1532   {
1533     /*canvas1->cd(3);
1534
1535     diffMCUnfolded->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e(unfolded)");
1536     //diffMCUnfolded->GetYaxis()->SetRangeUser(-20, 20);
1537     diffMCUnfolded->GetXaxis()->SetRangeUser(0, 200);
1538     diffMCUnfolded->DrawCopy("HIST");
1539
1540     TH1F* fluctuation = new TH1F("fluctuation", "fluctuation", 20, -5, 5);
1541     for (Int_t i=20; i<=diffMCUnfolded->GetNbinsX(); ++i)
1542       fluctuation->Fill(diffMCUnfolded->GetBinContent(i));
1543
1544     //new TCanvas; fluctuation->DrawCopy();
1545     delete fluctuation;*/
1546
1547     /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1548     legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
1549     legend->AddEntry(fMultiplicityMC, "MC");
1550     legend->AddEntry(fMultiplicityESD, "ESD");
1551     legend->Draw();*/
1552
1553     canvas1->cd(4);
1554     //residual->GetYaxis()->SetRangeUser(-0.2, 0.2);
1555     residual->GetXaxis()->SetRangeUser(0, 200);
1556     residual->DrawCopy();
1557
1558     canvas1->cd(5);
1559
1560     TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
1561     ratio->Divide(mcHist);
1562     ratio->SetTitle("Ratio;true multiplicity;Unfolded / MC");
1563     ratio->GetYaxis()->SetRangeUser(0.5, 1.5);
1564     ratio->GetXaxis()->SetRangeUser(0, 200);
1565     ratio->SetStats(kFALSE);
1566     ratio->Draw("HIST");
1567
1568     Double_t ratioChi2 = 0;
1569     fRatioAverage = 0;
1570     fLastChi2MCLimit = 0;
1571     for (Int_t i=2; i<=150; ++i)
1572     {
1573       Float_t value = ratio->GetBinContent(i) - 1;
1574       if (value > 1e8)
1575         value = 1e8; //prevent arithmetic exception
1576       else if (value < -1e8)
1577         value = -1e8;
1578
1579       ratioChi2 += value * value;
1580       fRatioAverage += TMath::Abs(value);
1581
1582       if (ratioChi2 < 0.1)
1583         fLastChi2MCLimit = i;
1584     }
1585     fRatioAverage /= 149;
1586
1587     printf("Sum over (ratio-1)^2 (2..150) is %f; average of |ratio-1| is %f\n", ratioChi2, fRatioAverage);
1588     // TODO FAKE
1589     fLastChi2MC = ratioChi2;
1590
1591     // FFT of ratio
1592     canvas1->cd(6);
1593     const Int_t kFFT = 128;
1594     Double_t fftReal[kFFT];
1595     Double_t fftImag[kFFT];
1596
1597     for (Int_t i=0; i<kFFT; ++i)
1598     {
1599       fftReal[i] = ratio->GetBinContent(i+1+10);
1600       // test: ;-)
1601       //fftReal[i] = cos(TMath::Pi() * 5 * 2 * i / 128);
1602       fftImag[i] = 0;
1603     }
1604
1605     FFT(-1, TMath::Nint(log(kFFT) / log(2)), fftReal, fftImag);
1606
1607     TH1* fftResult = (TH1*) ratio->Clone("fftResult");
1608     fftResult->SetTitle("FFT;true multiplicity;coeff. (10...137)");
1609     TH1* fftResult2 = (TH1*) ratio->Clone("fftResult2");
1610     TH1* fftResult3 = (TH1*) ratio->Clone("fftResult3");
1611     fftResult->Reset();
1612     fftResult2->Reset();
1613     fftResult3->Reset();
1614     fftResult->GetYaxis()->UnZoom();
1615     fftResult2->GetYaxis()->UnZoom();
1616     fftResult3->GetYaxis()->UnZoom();
1617
1618     //Printf("AFTER FFT");
1619     for (Int_t i=0; i<kFFT; ++i)
1620     {
1621       //Printf("%d: %f", i, fftReal[i]);
1622       fftResult->SetBinContent(i+1, fftReal[i]);
1623       /*if (i != 0 && TMath::Abs(fftReal[i]) > 0.5)
1624       {
1625         Printf("Nulled %d", i);
1626         fftReal[i] = 0;
1627       }*/
1628     }
1629
1630     fftResult->SetLineColor(1);
1631     fftResult->DrawCopy();
1632
1633
1634     //Printf("IMAG");
1635     for (Int_t i=0; i<kFFT; ++i)
1636     {
1637       //Printf("%d: %f", i, fftImag[i]);
1638       fftResult2->SetBinContent(i+1, fftImag[i]);
1639
1640       fftResult3->SetBinContent(i+1, TMath::Sqrt(fftReal[i] * fftReal[i] + fftImag[i] * fftImag[i]));
1641     }
1642
1643     fftResult2->SetLineColor(2);
1644     fftResult2->DrawCopy("SAME");
1645
1646     fftResult3->SetLineColor(4);
1647     fftResult3->DrawCopy("SAME");
1648
1649     for (Int_t i=1; i<kFFT - 1; ++i)
1650     {
1651       if (TMath::Sqrt(fftReal[i] * fftReal[i] + fftImag[i] * fftImag[i]) > 3)
1652       {
1653         fftReal[i-1] = 0;
1654         fftReal[i] = 0;
1655         fftReal[i+1] = 0;
1656         fftImag[i-1] = 0;
1657         fftImag[i] = 0;
1658         fftImag[i+1] = 0;
1659         //fftReal[i] = (fftReal[i-1] + fftReal[i+1]) / 2;
1660         //fftImag[i]  = (fftImag[i-1] + fftImag[i+1]) / 2;
1661         //Printf("Nulled %d to %f %f", i, fftReal[i], fftImag[i]);
1662       }
1663     }
1664
1665     FFT(1, TMath::Nint(log(kFFT) / log(2)), fftReal, fftImag);
1666
1667     TH1* fftResult4 = (TH1*) fftResult3->Clone("fftResult4");
1668     fftResult4->Reset();
1669
1670     //Printf("AFTER BACK-TRAFO");
1671     for (Int_t i=0; i<kFFT; ++i)
1672     {
1673       //Printf("%d: %f", i, fftReal[i]);
1674       fftResult4->SetBinContent(i+1+10, fftReal[i]);
1675     }
1676
1677     canvas1->cd(5);
1678     fftResult4->SetLineColor(4);
1679     fftResult4->DrawCopy("SAME");
1680
1681     // plot (MC - Unfolded) / error (MC)
1682     canvas1->cd(3);
1683
1684     TH1* diffMCUnfolded2 = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded2"));
1685     diffMCUnfolded2->Add(fMultiplicityESDCorrected[esdCorrId], -1);
1686
1687     Int_t ndfQual[kQualityRegions];
1688     for (Int_t region=0; region<kQualityRegions; ++region)
1689     {
1690       fQuality[region] = 0;
1691       ndfQual[region] = 0;
1692     }
1693
1694
1695     Double_t newChi2 = 0;
1696     Double_t newChi2Limit150 = 0;
1697     Int_t ndf = 0;
1698     for (Int_t i=1; i<=TMath::Min(diffMCUnfolded2->GetNbinsX(), fgkMaxParams+1); ++i)
1699     {
1700       Double_t value = 0;
1701       if (proj->GetBinError(i) > 0)
1702       {
1703         value = diffMCUnfolded2->GetBinContent(i) / proj->GetBinError(i);
1704         newChi2 += value * value;
1705         if (i > 1 && i <= mcBinLimit)
1706           newChi2Limit150 += value * value;
1707         ++ndf;
1708
1709         for (Int_t region=0; region<kQualityRegions; ++region)
1710           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
1711           {
1712             fQuality[region] += TMath::Abs(value);
1713             ++ndfQual[region];
1714           }
1715       }
1716
1717       diffMCUnfolded2->SetBinContent(i, value);
1718     }
1719
1720     // normalize region to the number of entries
1721     for (Int_t region=0; region<kQualityRegions; ++region)
1722     {
1723       if (ndfQual[region] > 0)
1724         fQuality[region] /= ndfQual[region];
1725       Printf("Quality parameter %d (%d <= mult <= %d) is %f with %d df", region, fgQualityRegionsB[region], fgQualityRegionsE[region], fQuality[region], ndfQual[region]);
1726     }
1727
1728     if (mcBinLimit > 1)
1729     {
1730       // TODO another FAKE
1731       fLastChi2MC = newChi2Limit150 / (mcBinLimit - 1);
1732       Printf("Chi2 (2..%d) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", mcBinLimit, newChi2Limit150, mcBinLimit - 1, fLastChi2MC);
1733     }
1734     else
1735       fLastChi2MC = -1;
1736
1737     Printf("Chi2 (full range) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", newChi2, ndf, ((ndf > 0) ? newChi2 / ndf : -1));
1738
1739     diffMCUnfolded2->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e(MC)");
1740     //diffMCUnfolded->GetYaxis()->SetRangeUser(-20, 20);
1741     diffMCUnfolded2->GetXaxis()->SetRangeUser(0, 200);
1742     diffMCUnfolded2->DrawCopy("HIST");
1743   }
1744
1745   canvas1->SaveAs(Form("%s.gif", canvas1->GetName()));
1746 }
1747
1748 //____________________________________________________________________
1749 void AliMultiplicityCorrection::FFT(Int_t dir, Int_t m, Double_t *x, Double_t *y)
1750 {
1751 /*-------------------------------------------------------------------------
1752    This computes an in-place complex-to-complex FFT
1753    x and y are the real and imaginary arrays of 2^m points.
1754    dir =  1 gives forward transform
1755    dir = -1 gives reverse transform
1756
1757      Formula: forward
1758                   N-1
1759                   ---
1760               1   \          - j k 2 pi n / N
1761       X(n) = ---   >   x(k) e                    = forward transform
1762               N   /                                n=0..N-1
1763                   ---
1764                   k=0
1765
1766       Formula: reverse
1767                   N-1
1768                   ---
1769                   \          j k 2 pi n / N
1770       X(n) =       >   x(k) e                    = forward transform
1771                   /                                n=0..N-1
1772                   ---
1773                   k=0
1774 */
1775
1776    Long_t   nn, i, i1, j, k, i2, l, l1, l2;
1777    Double_t c1, c2, tx, ty, t1, t2, u1, u2, z;
1778
1779    /* Calculate the number of points */
1780    nn = 1;
1781    for (i = 0; i < m; i++)
1782        nn *= 2;
1783
1784    /* Do the bit reversal */
1785    i2 = nn >> 1;
1786    j = 0;
1787    for (i= 0; i < nn - 1; i++) {
1788        if (i < j) {
1789            tx = x[i];
1790            ty = y[i];
1791            x[i] = x[j];
1792            y[i] = y[j];
1793            x[j] = tx;
1794            y[j] = ty;
1795        }
1796        k = i2;
1797        while (k <= j) {
1798            j -= k;
1799            k >>= 1;
1800        }
1801        j += k;
1802    }
1803
1804    /* Compute the FFT */
1805    c1 = -1.0;
1806    c2 = 0.0;
1807    l2 = 1;
1808    for (l = 0; l < m; l++) {
1809        l1 = l2;
1810        l2 <<= 1;
1811        u1 = 1.0;
1812        u2 = 0.0;
1813        for (j = 0;j < l1; j++) {
1814            for (i = j; i < nn; i += l2) {
1815                i1 = i + l1;
1816                t1 = u1 * x[i1] - u2 * y[i1];
1817                t2 = u1 * y[i1] + u2 * x[i1];
1818                x[i1] = x[i] - t1;
1819                y[i1] = y[i] - t2;
1820                x[i] += t1;
1821                y[i] += t2;
1822            }
1823            z =  u1 * c1 - u2 * c2;
1824            u2 = u1 * c2 + u2 * c1;
1825            u1 = z;
1826        }
1827        c2 = TMath::Sqrt((1.0 - c1) / 2.0);
1828        if (dir == 1)
1829            c2 = -c2;
1830        c1 = TMath::Sqrt((1.0 + c1) / 2.0);
1831    }
1832
1833    /* Scaling for forward transform */
1834    if (dir == 1) {
1835        for (i=0;i<nn;i++) {
1836            x[i] /= (Double_t)nn;
1837            y[i] /= (Double_t)nn;
1838        }
1839    }
1840 }
1841
1842 //____________________________________________________________________
1843 void AliMultiplicityCorrection::GetComparisonResults(Float_t* mc, Int_t* mcLimit, Float_t* residuals, Float_t* ratioAverage) const
1844 {
1845   // Returns the chi2 between the MC and the unfolded ESD as well as between the ESD and the folded unfolded ESD
1846   // These values are computed during DrawComparison, thus this function picks up the
1847   // last calculation
1848
1849   if (mc)
1850     *mc = fLastChi2MC;
1851   if (mcLimit)
1852     *mcLimit = fLastChi2MCLimit;
1853   if (residuals)
1854     *residuals = fLastChi2Residuals;
1855   if (ratioAverage)
1856     *ratioAverage = fRatioAverage;
1857 }
1858
1859 //____________________________________________________________________
1860 TH2F* AliMultiplicityCorrection::GetMultiplicityMC(Int_t i, EventType eventType)
1861 {
1862   //
1863   // returns the corresponding MC spectrum
1864   //
1865
1866   switch (eventType)
1867   {
1868     case kTrVtx : return fMultiplicityVtx[i]; break;
1869     case kMB : return fMultiplicityMB[i]; break;
1870     case kINEL : return fMultiplicityINEL[i]; break;
1871   }
1872
1873   return 0;
1874 }
1875
1876 //____________________________________________________________________
1877 TH1* AliMultiplicityCorrection::CalculateStdDev(TH1** results, Int_t max)
1878 {
1879   // calculate standard deviation of (results[0] - results[k]) k=1...max-1
1880   // per bin one gets: sigma(r[0] - r[n]) / r[0]
1881
1882   TH1* standardDeviation = (TH1*) results[0]->Clone("standardDeviation");
1883   standardDeviation->Reset();
1884
1885   for (Int_t x=1; x<=results[0]->GetNbinsX(); x++)
1886   {
1887     if (results[0]->GetBinContent(x) > 0)
1888     {
1889       Double_t average = 0;
1890       for (Int_t n=1; n<max; ++n)
1891         average += results[n]->GetBinContent(x) - results[0]->GetBinContent(x);
1892       average /= max-1;
1893
1894       Double_t variance = 0;
1895       for (Int_t n=1; n<max; ++n)
1896       {
1897         Double_t value = results[n]->GetBinContent(x) - results[0]->GetBinContent(x) - average;
1898         variance += value * value;
1899       }
1900       variance /= max-1;
1901
1902       Double_t standardDev = TMath::Sqrt(variance);
1903       standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x));
1904       //Printf("sigma_%d is %f value %f --> error %f", x, standardDev, results[0]->GetBinContent(x), standardDev / results[0]->GetBinContent(x));
1905     }
1906   }
1907
1908   return standardDeviation;
1909 }
1910
1911 //____________________________________________________________________
1912 TH1* AliMultiplicityCorrection::StatisticalUncertainty(MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, TH1* compareTo)
1913 {
1914   //
1915   // evaluates the uncertainty that arises from the non-infinite statistics in the response matrix
1916   // the function unfolds the spectrum using the default response matrix and several modified ones
1917   // the modified ones are created by randomizing each cell using poisson statistics with the mean = bin value
1918   // these unfolded results are compared to the first result gained with the default response OR to the histogram given
1919   // in <compareTo> (optional)
1920   //
1921   // returns the error assigned to the measurement
1922   //
1923
1924   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1925
1926   // initialize seed with current time
1927   gRandom->SetSeed(0);
1928
1929   const Int_t kErrorIterations = 150;
1930
1931   TH1* maxError = 0;
1932   TH1* firstResult = 0;
1933
1934   TH1** results = new TH1*[kErrorIterations];
1935
1936   for (Int_t n=0; n<kErrorIterations; ++n)
1937   {
1938     Printf("Iteration %d of %d...", n, kErrorIterations);
1939
1940     // only done for chi2 minimization
1941     Bool_t createBigBin = kFALSE;
1942     if (methodType == kChi2Minimization)
1943       createBigBin = kTRUE;
1944
1945     SetupCurrentHists(inputRange, fullPhaseSpace, eventType, createBigBin);
1946
1947     TH1* measured = (TH1*) fCurrentESD->Clone("measured");
1948
1949     if (n > 0)
1950     {
1951       if (randomizeResponse)
1952       {
1953         // randomize response matrix
1954         for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1955           for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1956             fCurrentCorrelation->SetBinContent(i, j, gRandom->Poisson(fCurrentCorrelation->GetBinContent(i, j)));
1957       }
1958
1959       if (randomizeMeasured)
1960       {
1961         // randomize measured spectrum
1962         for (Int_t x=1; x<=measured->GetNbinsX(); x++) // mult. axis
1963         {
1964           Int_t randomValue = gRandom->Poisson(fCurrentESD->GetBinContent(x));
1965           measured->SetBinContent(x, randomValue);
1966           measured->SetBinError(x, TMath::Sqrt(randomValue));
1967         }
1968       }
1969     }
1970
1971     // only for bayesian method we have to do it before the call to Unfold...
1972     if (methodType == kBayesian)
1973     {
1974       for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1975       {
1976         // with this it is normalized to 1
1977         Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1978
1979         // with this normalized to the given efficiency
1980         if (fCurrentEfficiency->GetBinContent(i) > 0)
1981           sum /= fCurrentEfficiency->GetBinContent(i);
1982         else
1983           sum = 0;
1984
1985         for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1986         {
1987           if (sum > 0)
1988           {
1989             fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1990             fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1991           }
1992           else
1993           {
1994             fCurrentCorrelation->SetBinContent(i, j, 0);
1995             fCurrentCorrelation->SetBinError(i, j, 0);
1996           }
1997         }
1998       }
1999     }
2000
2001     TH1* result = 0;
2002     if (n == 0 && compareTo)
2003     {
2004       // in this case we just store the histogram we want to compare to
2005       result = (TH1*) compareTo->Clone("compareTo");
2006       result->Sumw2();
2007     }
2008     else
2009     {
2010       result = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone(Form("result_%d", n));
2011
2012       Int_t returnCode = -1;
2013       if (methodType == kChi2Minimization)
2014       {
2015         returnCode = UnfoldWithMinuit(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result, kFALSE);
2016       }
2017       else if (methodType == kBayesian)
2018         returnCode = UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result, fgBayesianSmoothing, fgBayesianIterations);
2019
2020       if (returnCode != 0)
2021         return 0;
2022     }
2023
2024     // normalize
2025     result->Scale(1.0 / result->Integral());
2026
2027     if (n == 0)
2028     {
2029       firstResult = (TH1*) result->Clone("firstResult");
2030
2031       maxError = (TH1*) result->Clone("maxError");
2032       maxError->Reset();
2033     }
2034     else
2035     {
2036       // calculate ratio
2037       TH1* ratio = (TH1*) firstResult->Clone("ratio");
2038       ratio->Divide(result);
2039
2040       // find max. deviation
2041       for (Int_t x=1; x<=ratio->GetNbinsX(); x++)
2042         maxError->SetBinContent(x, TMath::Max(maxError->GetBinContent(x), TMath::Abs(1 - ratio->GetBinContent(x))));
2043
2044       delete ratio;
2045     }
2046
2047     results[n] = result;
2048   }
2049
2050   // find covariance matrix
2051   // results[n] is X_x
2052   // cov. matrix is M_xy = E ( (X_x - E(X_x)) * (X_y - E(X_y))), with E() = expectation value
2053
2054   Int_t nBins = results[0]->GetNbinsX();
2055   Float_t lowEdge = results[0]->GetXaxis()->GetBinLowEdge(1);
2056   Float_t upEdge = results[0]->GetXaxis()->GetBinUpEdge(nBins);
2057
2058   // find average, E(X)
2059   TProfile* average = new TProfile("average", "average", nBins, lowEdge, upEdge);
2060   for (Int_t n=1; n<kErrorIterations; ++n)
2061     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
2062       average->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x));
2063   //new TCanvas; average->DrawClone();
2064
2065   // find cov. matrix
2066   TProfile2D* covMatrix = new TProfile2D("covMatrix", "covMatrix", nBins, lowEdge, upEdge, nBins, lowEdge, upEdge);
2067
2068   for (Int_t n=1; n<kErrorIterations; ++n)
2069     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
2070       for (Int_t y=1; y<=results[n]->GetNbinsX(); y++)
2071       {
2072         // (X_x - E(X_x)) * (X_y - E(X_y)
2073         Float_t cov = (results[n]->GetBinContent(x) - average->GetBinContent(x)) * (results[n]->GetBinContent(y) - average->GetBinContent(y));
2074         covMatrix->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetXaxis()->GetBinCenter(y), cov);
2075       }
2076   new TCanvas; covMatrix->DrawClone("COLZ");
2077
2078 //   // fill 2D histogram that contains deviation from first
2079 //   TH2F* deviations = new TH2F("deviations", "deviations", nBins, lowEdge, upEdge, 1000, -0.01, 0.01);
2080 //   for (Int_t n=1; n<kErrorIterations; ++n)
2081 //     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
2082 //       deviations->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x) - results[0]->GetBinContent(x));
2083 //   //new TCanvas; deviations->DrawCopy("COLZ");
2084 // 
2085 //   // get standard deviation "by hand"
2086 //   for (Int_t x=1; x<=nBins; x++)
2087 //   {
2088 //     if (results[0]->GetBinContent(x) > 0)
2089 //     {
2090 //       TH1* proj = deviations->ProjectionY("projRMS", x, x, "e");
2091 //       Float_t standardDev = proj->GetRMS(); // this is standard deviation in fact
2092 //       //standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x));
2093 //       Printf("sigma_%d is %f value %f --> error %f", x, standardDev, results[0]->GetBinContent(x), standardDev / results[0]->GetBinContent(x));
2094 //     }
2095 //   }
2096
2097   TH1* standardDeviation = CalculateStdDev(results, kErrorIterations);
2098
2099   // compare maxError to RMS of profile (variable name: average)
2100   // first: calculate rms in percent of value
2101   TH1* rmsError = (TH1*) maxError->Clone("rmsError");
2102   rmsError->Reset();
2103
2104   // enable error to be standard deviation (see http://root.cern.ch/root/html/TProfile.html#TProfile:SetErrorOption)
2105   average->SetErrorOption("s");
2106   for (Int_t x=1; x<=rmsError->GetNbinsX(); x++)
2107     if (average->GetBinContent(x) > 0)
2108       rmsError->SetBinContent(x, average->GetBinError(x) / average->GetBinContent(x));
2109
2110   // find maxError deviation from average (not from "first result"/mc as above)
2111   TH1* maxError2 = (TH1*) maxError->Clone("maxError2");
2112   maxError2->Reset();
2113   for (Int_t n=1; n<kErrorIterations; ++n)
2114     for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
2115       if (average->GetBinContent(x) > 0)
2116         maxError2->SetBinContent(x, TMath::Max(maxError2->GetBinContent(x), TMath::Abs((results[n]->GetBinContent(x) - average->GetBinContent(x)) / average->GetBinContent(x))));
2117
2118   //new TCanvas; maxError2->DrawCopy(); rmsError->SetLineColor(2); rmsError->DrawCopy("SAME"); standardDeviation->SetLineColor(3); standardDeviation->DrawCopy("SAME");
2119
2120   // plot difference between average and MC/first result
2121   TH1* averageFirstRatio = (TH1*) results[0]->Clone("averageFirstRatio");
2122   averageFirstRatio->Reset();
2123   averageFirstRatio->Divide(results[0], average);
2124
2125   //new TCanvas; results[0]->DrawCopy(); average->SetLineColor(2); average->DrawClone("SAME");
2126   //new TCanvas; averageFirstRatio->DrawCopy();
2127
2128   // clean up
2129   for (Int_t n=0; n<kErrorIterations; ++n)
2130     delete results[n];
2131   delete[] results;
2132
2133   // fill into result histogram
2134   for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
2135     fMultiplicityESDCorrected[correlationID]->SetBinContent(i, firstResult->GetBinContent(i));
2136
2137   for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
2138     fMultiplicityESDCorrected[correlationID]->SetBinError(i, maxError->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
2139
2140   return standardDeviation;
2141 }
2142
2143 //____________________________________________________________________
2144 void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar, Int_t nIterations, TH1* initialConditions, Bool_t determineError)
2145 {
2146   //
2147   // correct spectrum using bayesian method
2148   //
2149
2150   // initialize seed with current time
2151   gRandom->SetSeed(0);
2152
2153   SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
2154
2155   // normalize correction for given nPart
2156   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2157   {
2158     // with this it is normalized to 1
2159     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
2160
2161     // with this normalized to the given efficiency
2162     if (fCurrentEfficiency->GetBinContent(i) > 0)
2163       sum /= fCurrentEfficiency->GetBinContent(i);
2164     else
2165       sum = 0;
2166
2167     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
2168     {
2169       if (sum > 0)
2170       {
2171         fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
2172         fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
2173       }
2174       else
2175       {
2176         fCurrentCorrelation->SetBinContent(i, j, 0);
2177         fCurrentCorrelation->SetBinError(i, j, 0);
2178       }
2179     }
2180   }
2181
2182   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
2183
2184   if (UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], regPar, nIterations) != 0)
2185     return;
2186
2187   if (!determineError)
2188   {
2189     Printf("AliMultiplicityCorrection::ApplyBayesianMethod: WARNING: No errors calculated.");
2190     return;
2191   }
2192
2193   // evaluate errors, this is done by randomizing the measured spectrum following Poission statistics
2194   // this (new) measured spectrum is then unfolded and the different to the result from the "real" measured
2195   // spectrum calculated. This is performed N times and the sigma is taken as the statistical
2196   // error of the unfolding method itself.
2197
2198   const Int_t kErrorIterations = 20;
2199
2200   printf("Spectrum unfolded. Determining error (%d iterations)...\n", kErrorIterations);
2201
2202   TH1* randomized = (TH1*) fCurrentESD->Clone("randomized");
2203   TH1* resultArray[kErrorIterations+1];
2204   for (Int_t n=0; n<kErrorIterations; ++n)
2205   {
2206     // randomize the content of clone following a poisson with the mean = the value of that bin
2207     for (Int_t x=1; x<=randomized->GetNbinsX(); x++) // mult. axis
2208     {
2209       Int_t randomValue = gRandom->Poisson(fCurrentESD->GetBinContent(x));
2210       //printf("%d --> %d\n", fCurrentESD->GetBinContent(x), randomValue);
2211       randomized->SetBinContent(x, randomValue);
2212       randomized->SetBinError(x, TMath::Sqrt(randomValue));
2213     }
2214
2215     TH1* result2 = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone("result2");
2216     result2->Reset();
2217     if (UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, randomized, initialConditions, result2, regPar, nIterations) != 0)
2218       return;
2219
2220     result2->Scale(1.0 / result2->Integral());
2221
2222     resultArray[n+1] = result2;
2223   }
2224   delete randomized;
2225
2226   resultArray[0] = fMultiplicityESDCorrected[correlationID];
2227   TH1* error = CalculateStdDev(resultArray, kErrorIterations+1);
2228
2229   for (Int_t n=0; n<kErrorIterations; ++n)
2230     delete resultArray[n+1];
2231
2232   for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
2233     fMultiplicityESDCorrected[correlationID]->SetBinError(i, error->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
2234
2235   delete error;
2236 }
2237
2238 //____________________________________________________________________
2239 Int_t AliMultiplicityCorrection::UnfoldWithBayesian(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult, Float_t regPar, Int_t nIterations)
2240 {
2241   //
2242   // unfolds a spectrum
2243   //
2244
2245   if (measured->Integral() <= 0)
2246   {
2247     Printf("AliMultiplicityCorrection::UnfoldWithBayesian: ERROR: The measured spectrum is empty");
2248     return -1;
2249   }
2250
2251   const Int_t kStartBin = 0;
2252
2253   const Int_t kMaxM = fgkMaxInput;  //<= fCurrentCorrelation->GetNbinsY(); // max measured axis
2254   const Int_t kMaxT = fgkMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis
2255
2256   // convergence limit: kMaxT * 0.001^2 = kMaxT * 1e-6 (e.g. 250 bins --> 2.5 e-4)
2257   const Double_t kConvergenceLimit = kMaxT * 1e-6;
2258
2259   // store information in arrays, to increase processing speed (~ factor 5)
2260   Double_t measuredCopy[kMaxM];
2261   Double_t measuredError[kMaxM];
2262   Double_t prior[kMaxT];
2263   Double_t result[kMaxT];
2264   Double_t efficiency[kMaxT];
2265
2266   Double_t response[kMaxT][kMaxM];
2267   Double_t inverseResponse[kMaxT][kMaxM];
2268
2269   // for normalization
2270   Float_t measuredIntegral = measured->Integral();
2271   for (Int_t m=0; m<kMaxM; m++)
2272   {
2273     measuredCopy[m] = measured->GetBinContent(m+1) / measuredIntegral;
2274     measuredError[m] = measured->GetBinError(m+1) / measuredIntegral;
2275
2276     for (Int_t t=0; t<kMaxT; t++)
2277     {
2278       response[t][m] = correlation->GetBinContent(t+1, m+1);
2279       inverseResponse[t][m] = 0;
2280     }
2281   }
2282
2283   for (Int_t t=0; t<kMaxT; t++)
2284   {
2285     efficiency[t] = aEfficiency->GetBinContent(t+1);
2286     prior[t] = measuredCopy[t];
2287     result[t] = 0;
2288   }
2289
2290   // pick prior distribution
2291   if (initialConditions)
2292   {
2293     printf("Using different starting conditions...\n");
2294     // for normalization
2295     Float_t inputDistIntegral = initialConditions->Integral();
2296     for (Int_t i=0; i<kMaxT; i++)
2297       prior[i] = initialConditions->GetBinContent(i+1) / inputDistIntegral;
2298   }
2299
2300   //TH1F* convergence = new TH1F("convergence", "convergence", 200, 0.5, 200.5);
2301
2302   // unfold...
2303   for (Int_t i=0; i<nIterations || nIterations < 0; i++)
2304   {
2305     //printf(" iteration %i \n", i);
2306
2307     // calculate IR from Beyes theorem
2308     // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
2309
2310     Double_t chi2Measured = 0;
2311     for (Int_t m=0; m<kMaxM; m++)
2312     {
2313       Float_t norm = 0;
2314       for (Int_t t = kStartBin; t<kMaxT; t++)
2315         norm += response[t][m] * prior[t];
2316
2317       // calc. chi2: (measured - response * prior) / error
2318       if (measuredError[m] > 0)
2319       {
2320         Double_t value = (measuredCopy[m] - norm) / measuredError[m];
2321         chi2Measured += value * value;
2322       }
2323
2324       if (norm > 0)
2325       {
2326         for (Int_t t = kStartBin; t<kMaxT; t++)
2327           inverseResponse[t][m] = response[t][m] * prior[t] / norm;
2328       }
2329       else
2330       {
2331         for (Int_t t = kStartBin; t<kMaxT; t++)
2332           inverseResponse[t][m] = 0;
2333       }
2334     }
2335     //Printf("chi2Measured of the last prior is %e", chi2Measured);
2336
2337     for (Int_t t = kStartBin; t<kMaxT; t++)
2338     {
2339       Float_t value = 0;
2340       for (Int_t m=0; m<kMaxM; m++)
2341         value += inverseResponse[t][m] * measuredCopy[m];
2342
2343       if (efficiency[t] > 0)
2344         result[t] = value / efficiency[t];
2345       else
2346         result[t] = 0;
2347     }
2348
2349     Double_t chi2LastIter = 0;
2350     // regularization (simple smoothing)
2351     for (Int_t t=kStartBin; t<kMaxT; t++)
2352     {
2353       Float_t newValue = 0;
2354       // 0 bin excluded from smoothing
2355       if (t > kStartBin+1 && t<kMaxT-1)
2356       {
2357         Float_t average = (result[t-1] + result[t] + result[t+1]) / 3;
2358
2359         // weight the average with the regularization parameter
2360         newValue = (1 - regPar) * result[t] + regPar * average;
2361       }
2362       else
2363         newValue = result[t];
2364
2365       // calculate chi2 (change from last iteration)
2366       if (prior[t] > 1e-5)
2367       {
2368         Double_t diff = (prior[t] - newValue) / prior[t];
2369         chi2LastIter += diff * diff;
2370       }
2371
2372       prior[t] = newValue;
2373     }
2374     //printf("Chi2 of %d iteration = %e\n", i, chi2LastIter);
2375     //convergence->Fill(i+1, chi2LastIter);
2376
2377     if (nIterations < 0 && chi2LastIter < kConvergenceLimit)
2378     {
2379       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);
2380       break;
2381     }
2382
2383     //if (i % 5 == 0)
2384     //  DrawComparison(Form("Bayesian_%d", i), mcTarget, correlationID, kTRUE, eventType);
2385   } // end of iterations
2386
2387   //new TCanvas; convergence->DrawCopy(); gPad->SetLogy();
2388   //delete convergence;
2389
2390   for (Int_t t=0; t<kMaxT; t++)
2391     aResult->SetBinContent(t+1, result[t]);
2392
2393   return 0;
2394
2395   // ********
2396   // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995
2397
2398   /*printf("Calculating covariance matrix. This may take some time...\n");
2399
2400   // TODO check if this is the right one...
2401   TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
2402
2403   Int_t xBins = hInverseResponseBayes->GetNbinsX();
2404   Int_t yBins = hInverseResponseBayes->GetNbinsY();
2405
2406   // calculate "unfolding matrix" Mij
2407   Float_t matrixM[251][251];
2408   for (Int_t i=1; i<=xBins; i++)
2409   {
2410     for (Int_t j=1; j<=yBins; j++)
2411     {
2412       if (fCurrentEfficiency->GetBinContent(i) > 0)
2413         matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i);
2414       else
2415         matrixM[i-1][j-1] = 0;
2416     }
2417   }
2418
2419   Float_t* vectorn = new Float_t[yBins];
2420   for (Int_t j=1; j<=yBins; j++)
2421     vectorn[j-1] = fCurrentESD->GetBinContent(j);
2422
2423   // first part of covariance matrix, depends on input distribution n(E)
2424   Float_t cov1[251][251];
2425
2426   Float_t nEvents = fCurrentESD->Integral(); // N
2427
2428   xBins = 20;
2429   yBins = 20;
2430
2431   for (Int_t k=0; k<xBins; k++)
2432   {
2433     printf("In Cov1: %d\n", k);
2434     for (Int_t l=0; l<yBins; l++)
2435     {
2436       cov1[k][l] = 0;
2437
2438       // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N)
2439       for (Int_t j=0; j<yBins; j++)
2440         cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j]
2441           * (1.0 - vectorn[j] / nEvents);
2442
2443       // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N
2444       for (Int_t i=0; i<yBins; i++)
2445         for (Int_t j=0; j<yBins; j++)
2446         {
2447           if (i == j)
2448             continue;
2449           cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i]
2450             * vectorn[j] / nEvents;
2451          }
2452     }
2453   }
2454
2455   printf("Cov1 finished\n");
2456
2457   TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov");
2458   cov->Reset();
2459
2460   for (Int_t i=1; i<=xBins; i++)
2461     for (Int_t j=1; j<=yBins; j++)
2462       cov->SetBinContent(i, j, cov1[i-1][j-1]);
2463
2464   new TCanvas;
2465   cov->Draw("COLZ");
2466
2467   // second part of covariance matrix, depends on response matrix
2468   Float_t cov2[251][251];
2469
2470   // Cov[P(Er|Cu), P(Es|Cu)] term
2471   Float_t covTerm[100][100][100];
2472   for (Int_t r=0; r<yBins; r++)
2473     for (Int_t u=0; u<xBins; u++)
2474       for (Int_t s=0; s<yBins; s++)
2475       {
2476         if (r == s)
2477           covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
2478             * (1.0 - hResponse->GetBinContent(u+1, r+1));
2479         else
2480           covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
2481             * hResponse->GetBinContent(u+1, s+1);
2482       }
2483
2484   for (Int_t k=0; k<xBins; k++)
2485     for (Int_t l=0; l<yBins; l++)
2486     {
2487       cov2[k][l] = 0;
2488       printf("In Cov2: %d %d\n", k, l);
2489       for (Int_t i=0; i<yBins; i++)
2490         for (Int_t j=0; j<yBins; j++)
2491         {
2492           //printf("In Cov2: %d %d %d %d\n", k, l, i, j);
2493           // calculate Cov(Mki, Mlj) = sum{ru},{su} ...
2494           Float_t tmpCov = 0;
2495           for (Int_t r=0; r<yBins; r++)
2496             for (Int_t u=0; u<xBins; u++)
2497               for (Int_t s=0; s<yBins; s++)
2498               {
2499                 if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0
2500                   || hResponse->GetBinContent(u+1, i+1) == 0)
2501                   continue;
2502
2503                 tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u)
2504                   * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u)
2505                   * covTerm[r][u][s];
2506               }
2507
2508           cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov;
2509         }
2510     }
2511
2512   printf("Cov2 finished\n");
2513
2514   for (Int_t i=1; i<=xBins; i++)
2515     for (Int_t j=1; j<=yBins; j++)
2516       cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]);
2517
2518   new TCanvas;
2519   cov->Draw("COLZ");*/
2520 }
2521
2522 //____________________________________________________________________
2523 Float_t AliMultiplicityCorrection::BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, Int_t k, Int_t i, Int_t r, Int_t u)
2524 {
2525   //
2526   // helper function for the covariance matrix of the bayesian method
2527   //
2528
2529   Float_t result = 0;
2530
2531   if (k == u && r == i)
2532     result += 1.0 / hResponse->GetBinContent(u+1, r+1);
2533
2534   if (k == u)
2535     result -= 1.0 / fCurrentEfficiency->GetBinContent(u+1);
2536
2537   if (r == i)
2538     result -= matrixM[u][i] * fCurrentEfficiency->GetBinContent(u+1) / hResponse->GetBinContent(u+1, i+1);
2539
2540   result *= matrixM[k][i];
2541
2542   return result;
2543 }
2544
2545 //____________________________________________________________________
2546 void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
2547 {
2548   //
2549   // correct spectrum using bayesian method
2550   //
2551
2552   Float_t regPar = 0;
2553
2554   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
2555   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
2556
2557   SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
2558   //normalize ESD
2559   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
2560
2561   // TODO should be taken from correlation map
2562   //TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
2563
2564   // normalize correction for given nPart
2565   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2566   {
2567     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
2568     //Double_t sum = sumHist->GetBinContent(i);
2569     if (sum <= 0)
2570       continue;
2571     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
2572     {
2573       // npart sum to 1
2574       fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
2575       fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
2576     }
2577   }
2578
2579   new TCanvas;
2580   fCurrentCorrelation->Draw("COLZ");
2581
2582   // FAKE
2583   fCurrentEfficiency = ((TH2*) fCurrentCorrelation)->ProjectionX("eff");
2584
2585   // pick prior distribution
2586   TH1F* hPrior = (TH1F*)fCurrentESD->Clone("prior");
2587   Float_t norm = 1; //hPrior->Integral();
2588   for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
2589     hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm);
2590
2591   // zero distribution
2592   TH1F* zero =  (TH1F*)hPrior->Clone("zero");
2593
2594   // define temp hist
2595   TH1F* hTemp = (TH1F*)fCurrentESD->Clone("temp");
2596   hTemp->Reset();
2597
2598   // just a shortcut
2599   TH2F* hResponse = (TH2F*) fCurrentCorrelation;
2600
2601   // unfold...
2602   Int_t iterations = 25;
2603   for (Int_t i=0; i<iterations; i++)
2604   {
2605     //printf(" iteration %i \n", i);
2606
2607     for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
2608     {
2609       Float_t value = 0;
2610       for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
2611         value += hResponse->GetBinContent(t, m) * hPrior->GetBinContent(t);
2612       hTemp->SetBinContent(m, value);
2613       //printf("%d %f %f %f\n", m, zero->GetBinContent(m), hPrior->GetBinContent(m), value);
2614     }
2615
2616     // regularization (simple smoothing)
2617     TH1F* hTrueSmoothed = (TH1F*) hTemp->Clone("truesmoothed");
2618
2619     for (Int_t t=2; t<hTrueSmoothed->GetNbinsX(); t++)
2620     {
2621       Float_t average = (hTemp->GetBinContent(t-1) / hTemp->GetBinWidth(t-1)
2622                          + hTemp->GetBinContent(t) / hTemp->GetBinWidth(t)
2623                          + hTemp->GetBinContent(t+1) / hTemp->GetBinWidth(t+1)) / 3.;
2624       average *= hTrueSmoothed->GetBinWidth(t);
2625
2626       // weight the average with the regularization parameter
2627       hTrueSmoothed->SetBinContent(t, (1-regPar) * hTemp->GetBinContent(t) + regPar * average);
2628     }
2629
2630     for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
2631       hTemp->SetBinContent(m, zero->GetBinContent(m) + hPrior->GetBinContent(m) - hTrueSmoothed->GetBinContent(m));
2632
2633     // fill guess
2634     for (Int_t t=1; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++)
2635     {
2636       fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTemp->GetBinContent(t));
2637       fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t)); // TODO
2638
2639       //printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
2640     }
2641
2642
2643     // calculate chi2 (change from last iteration)
2644     Double_t chi2 = 0;
2645
2646     // use smoothed true (normalized) as new prior
2647     Float_t norm = 1; //hTrueSmoothed->Integral();
2648
2649     for (Int_t t=1; t<hTrueSmoothed->GetNbinsX(); t++)
2650     {
2651       Float_t newValue = hTemp->GetBinContent(t)/norm;
2652       Float_t diff = hPrior->GetBinContent(t) - newValue;
2653       chi2 += (Double_t) diff * diff;
2654
2655       hPrior->SetBinContent(t, newValue);
2656     }
2657
2658     printf("Chi2 of %d iteration = %.10f\n", i, chi2);
2659
2660     //if (i % 5 == 0)
2661       DrawComparison(Form("Laszlo_%d", i), inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
2662
2663     delete hTrueSmoothed;
2664   } // end of iterations
2665
2666   DrawComparison("Laszlo", inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
2667 }
2668
2669 //____________________________________________________________________
2670 void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace)
2671 {
2672   //
2673   // correct spectrum using a simple Gaussian approach, that is model-dependent
2674   //
2675
2676   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
2677   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
2678
2679   SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE);
2680   //normalize ESD
2681   fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
2682
2683   TH1D* correction = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianMean"));
2684   correction->SetTitle("GaussianMean");
2685
2686   TH1D* correctionWidth = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianWidth"));
2687   correction->SetTitle("GaussianWidth");
2688
2689   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2690   {
2691     TH1D* proj = (dynamic_cast<TH2*> (fCurrentCorrelation))->ProjectionX("_px", i, i+1);
2692     proj->Fit("gaus", "0Q");
2693     correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
2694     correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
2695
2696     continue;
2697
2698     // draw for debugging
2699     new TCanvas;
2700     proj->DrawCopy();
2701     proj->GetFunction("gaus")->DrawCopy("SAME");
2702   }
2703
2704   TH1* target = fMultiplicityESDCorrected[correlationID];
2705
2706   Int_t nBins = target->GetNbinsX()*10+1;
2707   Float_t* binning = new Float_t[nBins];
2708   for (Int_t i=1; i<=target->GetNbinsX(); ++i)
2709     for (Int_t j=0; j<10; ++j)
2710       binning[(i-1)*10 + j] = target->GetXaxis()->GetBinLowEdge(i) + target->GetXaxis()->GetBinWidth(i) / 10 * j;
2711
2712   binning[nBins-1] = target->GetXaxis()->GetBinUpEdge(target->GetNbinsX());
2713
2714   TH1F* fineBinned = new TH1F("targetFineBinned", "targetFineBinned", nBins-1, binning);
2715
2716   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2717   {
2718     Float_t mean = correction->GetBinContent(i);
2719     Float_t width = correctionWidth->GetBinContent(i);
2720
2721     Int_t fillBegin = fineBinned->FindBin(mean - width * 5);
2722     Int_t fillEnd   = fineBinned->FindBin(mean + width * 5);
2723     //printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
2724
2725     for (Int_t j=fillBegin; j <= fillEnd; ++j)
2726     {
2727       fineBinned->AddBinContent(j, TMath::Gaus(fineBinned->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fCurrentESD->GetBinContent(i));
2728     }
2729   }
2730
2731   for (Int_t i=1; i<=target->GetNbinsX(); ++i)
2732   {
2733     Float_t sum = 0;
2734     for (Int_t j=1; j<=10; ++j)
2735       sum += fineBinned->GetBinContent((i-1)*10 + j);
2736     target->SetBinContent(i, sum / 10);
2737   }
2738
2739   delete[] binning;
2740
2741   DrawComparison("Gaussian", inputRange, fullPhaseSpace, kFALSE, GetMultiplicityMC(mcTarget, kTrVtx)->ProjectionY());
2742 }
2743
2744 //____________________________________________________________________
2745 TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap)
2746 {
2747   // runs the distribution given in inputMC through the response matrix identified by
2748   // correlationMap and produces a measured distribution
2749   // although it is a TH2F the vertex axis is not used at the moment and all entries are filled in mid-vertex
2750   // if normalized is set, inputMC is expected to be normalized to the bin width
2751
2752   if (!inputMC)
2753     return 0;
2754
2755   if (correlationMap < 0 || correlationMap >= kCorrHists)
2756     return 0;
2757
2758   // empty under/overflow bins in x, otherwise Project3D takes them into account
2759   TH3* hist = fCorrelation[correlationMap];
2760   for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
2761   {
2762     for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
2763     {
2764       hist->SetBinContent(0, y, z, 0);
2765       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
2766     }
2767   }
2768
2769   TH2* corr = (TH2*) hist->Project3D("zy");
2770   //corr->Rebin2D(2, 1);
2771   corr->Sumw2();
2772
2773   // normalize correction for given nPart
2774   for (Int_t i=1; i<=corr->GetNbinsX(); ++i)
2775   {
2776     Double_t sum = corr->Integral(i, i, 1, corr->GetNbinsY());
2777     if (sum <= 0)
2778       continue;
2779
2780     for (Int_t j=1; j<=corr->GetNbinsY(); ++j)
2781     {
2782       // npart sum to 1
2783       corr->SetBinContent(i, j, corr->GetBinContent(i, j) / sum);
2784       corr->SetBinError(i, j, corr->GetBinError(i, j) / sum);
2785     }
2786   }
2787
2788   TH2F* target = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("%s_measured", inputMC->GetName())));
2789   target->Reset();
2790
2791   for (Int_t meas=1; meas<=corr->GetNbinsY(); ++meas)
2792   {
2793     Float_t measured = 0;
2794     Float_t error = 0;
2795
2796     for (Int_t gen=1; gen<=corr->GetNbinsX(); ++gen)
2797     {
2798       Int_t mcGenBin = inputMC->GetXaxis()->FindBin(corr->GetXaxis()->GetBinCenter(gen));
2799
2800       measured += inputMC->GetBinContent(mcGenBin) * corr->GetBinContent(gen, meas);
2801       error += inputMC->GetBinError(mcGenBin) * corr->GetBinContent(gen, meas);
2802     }
2803
2804     //printf("%f +- %f ; %f +- %f \n", inputMC->GetBinContent(meas), inputMC->GetBinError(meas), measured, error);
2805
2806     target->SetBinContent(1 + target->GetNbinsX() / 2, meas, measured);
2807     target->SetBinError(1 + target->GetNbinsX() / 2, meas, error);
2808   }
2809
2810   return target;
2811 }
2812
2813 //____________________________________________________________________
2814 void AliMultiplicityCorrection::SetGenMeasFromFunc(TF1* inputMC, Int_t id)
2815 {
2816   // uses the given function to fill the input MC histogram and generates from that
2817   // the measured histogram by applying the response matrix
2818   // this can be used to evaluate if the methods work indepedently of the input
2819   // distribution
2820   // WARNING does not respect the vertex distribution, just fills central vertex bin
2821
2822   if (!inputMC)
2823     return;
2824
2825   if (id < 0 || id >= kESDHists)
2826     return;
2827
2828   TH2F* mc = fMultiplicityVtx[id];
2829
2830   mc->Reset();
2831
2832   for (Int_t i=1; i<=mc->GetNbinsY(); ++i)
2833   {
2834     mc->SetBinContent(mc->GetNbinsX() / 2 + 1, i, inputMC->Eval(mc->GetYaxis()->GetBinCenter(i)) * mc->GetYaxis()->GetBinWidth(i));
2835     mc->SetBinError(mc->GetNbinsX() / 2 + 1, i, 0);
2836   }
2837
2838   //new TCanvas;
2839   //mc->Draw("COLZ");
2840
2841   TH1* proj = mc->ProjectionY();
2842
2843   TString tmp(fMultiplicityESD[id]->GetName());
2844   delete fMultiplicityESD[id];
2845   fMultiplicityESD[id] = CalculateMultiplicityESD(proj, id);
2846   fMultiplicityESD[id]->SetName(tmp);
2847 }