introducing multiplicity measurement with the ITS or TPC
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AliMultiplicityCorrection.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 // This class is used to store correction maps, raw input and results of the multiplicity
19 // measurement with the ITS or TPC
20 // It also contains functions to correct the spectrum using different methods.
21 //
22 //  Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch
23
24 #include "AliMultiplicityCorrection.h"
25
26 #include <TFile.h>
27 #include <TH1F.h>
28 #include <TH2F.h>
29 #include <TH3F.h>
30 #include <TDirectory.h>
31 #include <TVirtualFitter.h>
32 #include <TCanvas.h>
33 #include <TString.h>
34
35 ClassImp(AliMultiplicityCorrection)
36
37 const Int_t AliMultiplicityCorrection::fgMaxParams = 110;
38 TH1* AliMultiplicityCorrection::fCurrentMinuitESD = 0;
39 TH1* AliMultiplicityCorrection::fCurrentMinuitCorrelation = 0;
40
41 //____________________________________________________________________
42 AliMultiplicityCorrection::AliMultiplicityCorrection() :
43   TNamed()
44 {
45   //
46   // default constructor
47   //
48
49   for (Int_t i = 0; i < kESDHists; ++i)
50     fMultiplicityESD[i] = 0;
51
52   for (Int_t i = 0; i < kMCHists; ++i)
53     fMultiplicityMC[i] = 0;
54
55   for (Int_t i = 0; i < kCorrHists; ++i)
56   {
57     fCorrelation[i] = 0;
58     fMultiplicityESDCorrected[i] = 0;
59   }
60 }
61
62 //____________________________________________________________________
63 AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const Char_t* title) :
64   TNamed(name, title)
65 {
66   //
67   // named constructor
68   //
69
70   // do not add this hists to the directory
71   Bool_t oldStatus = TH1::AddDirectoryStatus();
72   TH1::AddDirectory(kFALSE);
73
74   Float_t binLimitsVtx[] = {-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10};
75   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,
76                           10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
77                           20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
78                           30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
79                           40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5,
80                           50.5, 55.5, 60.5, 65.5, 70.5, 75.5, 80.5, 85.5, 90.5, 95.5,
81                           100.5, 105.5, 110.5, 115.5, 120.5, 125.5, 130.5, 135.5, 140.5, 145.5,
82                           150.5, 160.5, 170.5, 180.5, 190.5, 200.5, 210.5, 220.5, 230.5, 240.5,
83                           250.5, 275.5, 300.5, 325.5, 350.5, 375.5, 400.5, 425.5, 450.5, 475.5,
84                           500.5, 525.5, 550.5, 575.5, 600.5, 625.5, 650.5, 675.5, 700.5, 725.5,
85                           750.5, 775.5, 800.5, 825.5, 850.5, 875.5, 900.5, 925.5, 950.5, 975.5,
86                           1000.5 };
87
88   #define BINNING 110, binLimitsN
89
90   for (Int_t i = 0; i < kESDHists; ++i)
91     fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", 10, binLimitsVtx, BINNING);
92
93   for (Int_t i = 0; i < kMCHists; ++i)
94   {
95     fMultiplicityMC[i] = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("fMultiplicityMC%d", i)));
96     fMultiplicityMC[i]->SetTitle("fMultiplicityMC");
97   }
98
99   for (Int_t i = 0; i < kCorrHists; ++i)
100   {
101     fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;Npart;Ntracks", 10, binLimitsVtx, BINNING, BINNING);
102     fMultiplicityESDCorrected[i] = new TH1F(Form("fMultiplicityESDCorrected%d", i), "fMultiplicityESDCorrected;Npart;dN/dN", BINNING);
103   }
104
105   TH1::AddDirectory(oldStatus);
106 }
107
108 //____________________________________________________________________
109 AliMultiplicityCorrection::~AliMultiplicityCorrection()
110 {
111   //
112   // Destructor
113   //
114 }
115
116 //____________________________________________________________________
117 Long64_t AliMultiplicityCorrection::Merge(TCollection* list)
118 {
119   // Merge a list of AliMultiplicityCorrection objects with this (needed for
120   // PROOF).
121   // Returns the number of merged objects (including this).
122
123   if (!list)
124     return 0;
125
126   if (list->IsEmpty())
127     return 1;
128
129   TIterator* iter = list->MakeIterator();
130   TObject* obj;
131
132   // collections of all histograms
133   TList collections[kESDHists+kMCHists+kCorrHists*2];
134
135   Int_t count = 0;
136   while ((obj = iter->Next())) {
137
138     AliMultiplicityCorrection* entry = dynamic_cast<AliMultiplicityCorrection*> (obj);
139     if (entry == 0) 
140       continue;
141
142     for (Int_t i = 0; i < kESDHists; ++i)
143       collections[i].Add(entry->fMultiplicityESD[i]);
144
145     for (Int_t i = 0; i < kMCHists; ++i)
146       collections[kESDHists+i].Add(entry->fMultiplicityMC[i]);
147
148     for (Int_t i = 0; i < kCorrHists; ++i)
149       collections[kESDHists+kMCHists+i].Add(entry->fCorrelation[i]);
150
151     for (Int_t i = 0; i < kCorrHists; ++i)
152       collections[kESDHists+kMCHists+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]);
153
154     count++;
155   }
156
157   for (Int_t i = 0; i < kESDHists; ++i)
158     fMultiplicityESD[i]->Merge(&collections[i]);
159
160   for (Int_t i = 0; i < kMCHists; ++i)
161     fMultiplicityMC[i]->Merge(&collections[kESDHists+i]);
162
163   for (Int_t i = 0; i < kCorrHists; ++i)
164     fCorrelation[i]->Merge(&collections[kESDHists+kMCHists+i]);
165
166   for (Int_t i = 0; i < kCorrHists; ++i)
167     fMultiplicityESDCorrected[i]->Merge(&collections[kESDHists+kMCHists+kCorrHists+i]);
168
169   delete iter;
170
171   return count+1;
172 }
173
174 //____________________________________________________________________
175 Bool_t AliMultiplicityCorrection::LoadHistograms(const Char_t* dir)
176 {
177   //
178   // loads the histograms from a file
179   // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
180   //
181
182   if (!dir)
183     dir = GetName();
184
185   if (!gDirectory->cd(dir))
186     return kFALSE;
187
188   // TODO memory leak. old histograms needs to be deleted.
189
190   Bool_t success = kTRUE;
191
192   for (Int_t i = 0; i < kESDHists; ++i)
193   {
194     fMultiplicityESD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityESD[i]->GetName()));
195     if (!fMultiplicityESD[i])
196       success = kFALSE;
197   }
198
199   for (Int_t i = 0; i < kMCHists; ++i)
200   {
201     fMultiplicityMC[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityMC[i]->GetName()));
202     if (!fMultiplicityMC[i])
203       success = kFALSE;
204   }
205
206   for (Int_t i = 0; i < kCorrHists; ++i)
207   {
208     fCorrelation[i] = dynamic_cast<TH3F*> (gDirectory->Get(fCorrelation[i]->GetName()));
209     if (!fCorrelation[i])
210       success = kFALSE;
211     fMultiplicityESDCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fMultiplicityESDCorrected[i]->GetName()));
212     if (!fMultiplicityESDCorrected[i])
213       success = kFALSE;
214   }
215
216   gDirectory->cd("..");
217
218   return success;
219 }
220
221 //____________________________________________________________________
222 void AliMultiplicityCorrection::SaveHistograms()
223 {
224   //
225   // saves the histograms
226   //
227
228   gDirectory->mkdir(GetName());
229   gDirectory->cd(GetName());
230
231   for (Int_t i = 0; i < kESDHists; ++i)
232     if (fMultiplicityESD[i])
233       fMultiplicityESD[i]->Write();
234
235   for (Int_t i = 0; i < kMCHists; ++i)
236     if (fMultiplicityMC[i])
237       fMultiplicityMC[i]->Write();
238
239   for (Int_t i = 0; i < kCorrHists; ++i)
240   {
241     if (fCorrelation[i])
242       fCorrelation[i]->Write();
243     if (fMultiplicityESDCorrected[i])
244       fMultiplicityESDCorrected[i]->Write();
245   }
246
247   gDirectory->cd("..");
248 }
249
250 //____________________________________________________________________
251 void AliMultiplicityCorrection::FillGenerated(Float_t vtx, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll)
252 {
253   //
254   // Fills an event from MC
255   //
256
257   fMultiplicityMC[0]->Fill(vtx, generated05);
258   fMultiplicityMC[1]->Fill(vtx, generated10);
259   fMultiplicityMC[2]->Fill(vtx, generated15);
260   fMultiplicityMC[3]->Fill(vtx, generated20);
261   fMultiplicityMC[4]->Fill(vtx, generatedAll);
262 }
263
264 //____________________________________________________________________
265 void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20)
266 {
267   //
268   // Fills an event from ESD
269   //
270
271   fMultiplicityESD[0]->Fill(vtx, measured05);
272   fMultiplicityESD[1]->Fill(vtx, measured10);
273   fMultiplicityESD[2]->Fill(vtx, measured15);
274   fMultiplicityESD[3]->Fill(vtx, measured20);
275 }
276
277 //____________________________________________________________________
278 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)
279 {
280   //
281   // Fills an event into the correlation map with the information from MC and ESD
282   //
283
284   fCorrelation[0]->Fill(vtx, generated05, measured05);
285   fCorrelation[1]->Fill(vtx, generated10, measured10);
286   fCorrelation[2]->Fill(vtx, generated15, measured15);
287   fCorrelation[3]->Fill(vtx, generated20, measured20);
288
289   fCorrelation[4]->Fill(vtx, generatedAll, measured05);
290   fCorrelation[5]->Fill(vtx, generatedAll, measured10);
291   fCorrelation[6]->Fill(vtx, generatedAll, measured15);
292   fCorrelation[7]->Fill(vtx, generatedAll, measured20);
293 }
294
295 //____________________________________________________________________
296 Double_t AliMultiplicityCorrection::MinuitHomogenityPol0(Double_t *params)
297 {
298   // homogenity term for minuit fitting
299   // pure function of the parameters
300   // prefers constant function (pol0)
301
302   Double_t chi2 = 0;
303
304   for (Int_t i=1; i<fgMaxParams; ++i)
305   {
306     if (params[i] == 0)
307       continue;
308
309     Double_t right  = params[i] / fCurrentMinuitESD->GetBinWidth(i+1);
310     Double_t left   = params[i-1] / fCurrentMinuitESD->GetBinWidth(i);
311
312     Double_t diff = (right - left) / right;
313     chi2 += diff * diff;
314   }
315
316   return chi2;
317 }
318
319 //____________________________________________________________________
320 Double_t AliMultiplicityCorrection::MinuitHomogenityPol1(Double_t *params)
321 {
322   // homogenity term for minuit fitting
323   // pure function of the parameters
324   // prefers linear function (pol1)
325
326   Double_t chi2 = 0;
327
328   for (Int_t i=2; i<fgMaxParams; ++i)
329   {
330     if (params[i] == 0 || params[i-1] == 0)
331       continue;
332
333     Double_t right  = params[i] / fCurrentMinuitESD->GetBinWidth(i+1);
334     Double_t middle = params[i-1] / fCurrentMinuitESD->GetBinWidth(i);
335     Double_t left   = params[i-2] / fCurrentMinuitESD->GetBinWidth(i-1);
336
337     Double_t der1 = (right - middle) / fCurrentMinuitESD->GetBinWidth(i);
338     Double_t der2 = (middle - left)  / fCurrentMinuitESD->GetBinWidth(i-1);
339
340     Double_t diff = der1 - der2;
341
342     // TODO give additonal weight to big bins
343     chi2 += diff * diff * fCurrentMinuitESD->GetBinWidth(i) * fCurrentMinuitESD->GetBinWidth(i-1);
344
345     //printf("%d --> %f\n", i, diff);
346   }
347
348   return chi2 / 1e5 / 2;
349 }
350
351 //____________________________________________________________________
352 void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t)
353 {
354   //
355   // fit function for minuit
356   //
357
358   // TODO take errors into account
359
360   static Int_t callCount = 0;
361
362   Double_t chi2FromFit = 0;
363
364   // loop over multiplicity (x axis of fMultiplicityESD)
365   for (Int_t i=1; i<=fCurrentMinuitESD->GetNbinsX(); ++i)
366   {
367     if (i > fCurrentMinuitCorrelation->GetNbinsY())
368       break;
369
370     Double_t sum = 0;
371     //Double_t error = 0;
372     // loop over generated particles (x axis of fCorrelation) that resulted in reconstruction of i tracks
373     for (Int_t j=1; j<=fCurrentMinuitCorrelation->GetNbinsX(); ++j)
374     {
375       if (j > fgMaxParams)
376         break;
377
378       sum += fCurrentMinuitCorrelation->GetBinContent(j, i) * params[j-1];
379
380       //if (params[j-1] > 0)
381       //  error += fCurrentMinuitCorrelation->GetBinError(j, i) * fCurrentMinuitCorrelation->GetBinError(j, i) * params[j-1];
382       //printf("%f  ", sum);
383     }
384
385     Double_t diff = fCurrentMinuitESD->GetBinContent(i) - sum;
386     if (fCurrentMinuitESD->GetBinContent(i) > 0)
387       diff /= fCurrentMinuitESD->GetBinContent(i);
388     else
389       diff /= fCurrentMinuitESD->Integral();
390     //error = TMath::Sqrt(error) + fCurrentMinuitESD->GetBinError(i);
391     //if (error <= 0)
392     // error = 1;
393
394     //Double_t tmp = diff / error;
395     //chi2 += tmp * tmp;
396     chi2FromFit += diff * diff;
397
398     //printf("\nExpected sum = %f; Diff for bin %d is %f\n**********************************\n", fCurrentMinuitESD->GetBinContent(i), i, diff);
399     //printf("Diff for bin %d is %f\n", i, diff);
400   }
401
402   Double_t penaltyVal = MinuitHomogenityPol1(params);
403
404   chi2 = chi2FromFit * chi2FromFit + penaltyVal * penaltyVal;
405
406   if ((callCount++ % 100) == 0)
407     printf("%f %f --> %f\n", chi2FromFit, penaltyVal, chi2);
408 }
409
410 //____________________________________________________________________
411 void AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace)
412 {
413   //
414   // correct spectrum using minuit chi2 method
415   //
416
417   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
418   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
419
420   fCurrentMinuitESD = fMultiplicityESD[inputRange]->ProjectionY();
421   fCurrentMinuitESD->Sumw2();
422
423   // empty under/overflow bins in x, otherwise Project3D takes them into account
424   TH3* hist = fCorrelation[correlationID];
425   for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
426   {
427     for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
428     {
429       hist->SetBinContent(0, y, z, 0);
430       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
431     }
432   }
433
434   fCurrentMinuitCorrelation = hist->Project3D("zy");
435   fCurrentMinuitCorrelation->Sumw2();
436
437   // normalize correction for given nPart
438   for (Int_t i=1; i<=fCurrentMinuitCorrelation->GetNbinsX(); ++i)
439   {
440     Double_t sum = fCurrentMinuitCorrelation->Integral(i, i, 1, fCurrentMinuitCorrelation->GetNbinsY());
441     if (sum <= 0)
442       continue;
443     for (Int_t j=1; j<=fCurrentMinuitCorrelation->GetNbinsY(); ++j)
444     {
445       // npart sum to 1
446       fCurrentMinuitCorrelation->SetBinContent(i, j, fCurrentMinuitCorrelation->GetBinContent(i, j) / sum);
447       fCurrentMinuitCorrelation->SetBinError(i, j, fCurrentMinuitCorrelation->GetBinError(i, j) / sum);
448     }
449   }
450
451   // small hack to get around charge conservation for full phase space ;-)
452   /*if (fullPhaseSpace)
453   {
454     for (Int_t i=2; i<=50; i+=2)
455       for (Int_t j=1; j<=fCurrentMinuitCorrelation->GetNbinsY(); ++j)
456         fCurrentMinuitCorrelation->SetBinContent(i, j, fCurrentMinuitCorrelation->GetBinContent(i-1, j));
457   }*/
458
459   TCanvas* canvas1 = new TCanvas("ApplyMinuitFit", "ApplyMinuitFit", 800, 400);
460   canvas1->Divide(2, 1);
461   canvas1->cd(1);
462   fCurrentMinuitESD->DrawCopy();
463   canvas1->cd(2);
464   fCurrentMinuitCorrelation->DrawCopy("COLZ");
465
466   // Initialize TMinuit via generic fitter interface
467   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgMaxParams);
468
469   minuit->SetFCN(MinuitFitFunction);
470
471   Double_t results[fgMaxParams];
472
473   //TH1* proj = fMultiplicityMC[mcTarget]->ProjectionY();
474   for (Int_t i=0; i<fgMaxParams; ++i)
475   {
476     //results[i] = 1;
477     results[i] = fCurrentMinuitESD->GetBinContent(i+1);
478     //results[i] = proj->GetBinContent(i+1);
479     minuit->SetParameter(i, Form("param%d", i), results[i], ((results[i] > 1) ? (results[i] / 10) : 0), 0, fCurrentMinuitESD->GetMaximum() * 100);
480   }
481
482   Int_t dummy;
483   Double_t chi2;
484   MinuitFitFunction(dummy, 0, chi2, results, 0);
485   printf("Chi2 of right solution is = %f\n", chi2);
486
487   //return;
488
489   Double_t arglist[100];
490   arglist[0] = 100000;
491   //minuit->ExecuteCommand("SET PRINT", arglist, 1);
492   minuit->ExecuteCommand("MIGRAD", arglist, 1);
493   //minuit->ExecuteCommand("MINOS", arglist, 0);
494
495   for (Int_t i=0; i<fgMaxParams; ++i)
496   {
497     results[i] = minuit->GetParameter(i);
498     fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, results[i]);
499     fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0);
500   }
501
502   printf("Penalty is %f\n", MinuitHomogenityPol1(results));
503
504   fMultiplicityESDCorrected[correlationID]->GetXaxis()->SetRangeUser(0, fgMaxParams);
505
506   DrawComparison(mcTarget, correlationID);
507 }
508
509 //____________________________________________________________________
510 void AliMultiplicityCorrection::NormalizeToBinWidth(TH1* hist)
511 {
512   //
513   // normalizes a 1-d histogram to its bin width
514   //
515
516   for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
517   {
518     hist->SetBinContent(i, hist->GetBinContent(i) / hist->GetBinWidth(i));
519     hist->SetBinError(i, hist->GetBinError(i) / hist->GetBinWidth(i));
520   }
521 }
522
523 //____________________________________________________________________
524 void AliMultiplicityCorrection::NormalizeToBinWidth(TH2* hist)
525 {
526   //
527   // normalizes a 2-d histogram to its bin width (x width * y width)
528   //
529
530   for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
531     for (Int_t j=1; j<=hist->GetNbinsY(); ++j)
532     {
533       Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j);
534       hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor);
535       hist->SetBinError(i, j, hist->GetBinError(i, j) / factor);
536     }
537 }
538
539 //____________________________________________________________________
540 void AliMultiplicityCorrection::ApplyMinuitFitAll()
541 {
542   //
543   // fit all eta ranges
544   //
545
546   for (Int_t i=0; i<kESDHists; ++i)
547   {
548     ApplyMinuitFit(i, kFALSE);
549     ApplyMinuitFit(i, kTRUE);
550   }
551 }
552
553 //____________________________________________________________________
554 void AliMultiplicityCorrection::DrawHistograms()
555 {
556   printf("ESD:\n");
557
558   TCanvas* canvas1 = new TCanvas("fMultiplicityESD", "fMultiplicityESD", 900, 600);
559   canvas1->Divide(3, 2);
560   for (Int_t i = 0; i < kESDHists; ++i)
561   {
562     canvas1->cd(i+1);
563     fMultiplicityESD[i]->DrawCopy("COLZ");
564     printf("%d --> %f\n", i, (Float_t) fMultiplicityESD[i]->ProjectionY()->GetMean());
565   }
566
567   printf("MC:\n");
568
569   TCanvas* canvas2 = new TCanvas("fMultiplicityMC", "fMultiplicityMC", 900, 600);
570   canvas2->Divide(3, 2);
571   for (Int_t i = 0; i < kMCHists; ++i)
572   {
573     canvas2->cd(i+1);
574     fMultiplicityMC[i]->DrawCopy("COLZ");
575     printf("%d --> %f\n", i, fMultiplicityMC[i]->ProjectionY()->GetMean());
576   }
577
578   TCanvas* canvas3 = new TCanvas("fCorrelation", "fCorrelation", 900, 900);
579   canvas3->Divide(3, 3);
580   for (Int_t i = 0; i < kCorrHists; ++i)
581   {
582     canvas3->cd(i+1);
583     TH1* proj = fCorrelation[i]->Project3D("zy");
584     proj->DrawCopy("COLZ");
585   }
586 }
587
588 //____________________________________________________________________
589 void AliMultiplicityCorrection::DrawComparison(Int_t mcID, Int_t esdCorrId)
590 {
591   TString name;
592   name.Form("DrawComparison_%d_%d", mcID, esdCorrId);
593
594   TCanvas* canvas1 = new TCanvas(name, name, 900, 300);
595   canvas1->Divide(3, 1);
596
597   canvas1->cd(1);
598   TH1* proj = fMultiplicityMC[mcID]->ProjectionY();
599   NormalizeToBinWidth(proj);
600   NormalizeToBinWidth(fMultiplicityESDCorrected[esdCorrId]);
601   proj->GetXaxis()->SetRangeUser(0, 200);
602   proj->DrawCopy("HIST");
603   gPad->SetLogy();
604
605   fMultiplicityESDCorrected[esdCorrId]->SetMarkerStyle(3);
606   fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAMEP");
607
608   canvas1->cd(2);
609   fMultiplicityESDCorrected[esdCorrId]->DrawCopy("");
610
611   canvas1->cd(3);
612   TH1* clone = dynamic_cast<TH1*> (proj->Clone("clone"));
613   clone->Divide(fMultiplicityESDCorrected[esdCorrId]);
614   clone->SetTitle("Ratio;Npart;MC/ESD");
615   clone->GetYaxis()->SetRangeUser(0.8, 1.2);
616   clone->GetXaxis()->SetRangeUser(0, 200);
617   clone->DrawCopy("HIST");
618
619   /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
620   legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
621   legend->AddEntry(fMultiplicityMC, "MC");
622   legend->AddEntry(fMultiplicityESD, "ESD");
623   legend->Draw();*/
624 }
625
626 //____________________________________________________________________
627 void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace)
628 {
629   //
630   // correct spectrum using bayesian method
631   //
632
633   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
634   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
635
636   fCurrentMinuitESD = fMultiplicityESD[inputRange]->ProjectionY();
637   fCurrentMinuitESD->Sumw2();
638
639   // empty under/overflow bins in x, otherwise Project3D takes them into account
640   TH3* hist = fCorrelation[correlationID];
641   for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
642   {
643     for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
644     {
645       hist->SetBinContent(0, y, z, 0);
646       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
647     }
648   }
649
650   fCurrentMinuitCorrelation = hist->Project3D("zy");
651   fCurrentMinuitCorrelation->Sumw2();
652
653   // normalize correction for given nPart
654   for (Int_t i=1; i<=fCurrentMinuitCorrelation->GetNbinsX(); ++i)
655   {
656     Double_t sum = fCurrentMinuitCorrelation->Integral(i, i, 1, fCurrentMinuitCorrelation->GetNbinsY());
657     if (sum <= 0)
658       continue;
659     for (Int_t j=1; j<=fCurrentMinuitCorrelation->GetNbinsY(); ++j)
660     {
661       // npart sum to 1
662       fCurrentMinuitCorrelation->SetBinContent(i, j, fCurrentMinuitCorrelation->GetBinContent(i, j) / sum);
663       fCurrentMinuitCorrelation->SetBinError(i, j, fCurrentMinuitCorrelation->GetBinError(i, j) / sum);
664     }
665   }
666
667   // TODO to be implemented
668 }