22b47d2055538508b409b25359f3a4ddbecafd99
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AliMultiplicityCorrection.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 // This class is used to store correction maps, raw input and results of the multiplicity
19 // measurement with the ITS or TPC
20 // It also contains functions to correct the spectrum using different methods.
21 //
22 //  Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch
23
24 #include "AliMultiplicityCorrection.h"
25
26 #include <TFile.h>
27 #include <TH1F.h>
28 #include <TH2F.h>
29 #include <TH3F.h>
30 #include <TDirectory.h>
31 #include <TVirtualFitter.h>
32 #include <TCanvas.h>
33 #include <TString.h>
34 #include <TF1.h>
35 #include <TMath.h>
36
37 ClassImp(AliMultiplicityCorrection)
38
39 const Int_t AliMultiplicityCorrection::fgMaxParams = 90;
40 TH1* AliMultiplicityCorrection::fCurrentESD = 0;
41 TH1* AliMultiplicityCorrection::fCurrentCorrelation = 0;
42
43 //____________________________________________________________________
44 AliMultiplicityCorrection::AliMultiplicityCorrection() :
45   TNamed()
46 {
47   //
48   // default constructor
49   //
50
51   for (Int_t i = 0; i < kESDHists; ++i)
52     fMultiplicityESD[i] = 0;
53
54   for (Int_t i = 0; i < kMCHists; ++i)
55     fMultiplicityMC[i] = 0;
56
57   for (Int_t i = 0; i < kCorrHists; ++i)
58   {
59     fCorrelation[i] = 0;
60     fMultiplicityESDCorrected[i] = 0;
61   }
62 }
63
64 //____________________________________________________________________
65 AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const Char_t* title) :
66   TNamed(name, title)
67 {
68   //
69   // named constructor
70   //
71
72   // do not add this hists to the directory
73   Bool_t oldStatus = TH1::AddDirectoryStatus();
74   TH1::AddDirectory(kFALSE);
75
76   Float_t binLimitsVtx[] = {-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10};
77   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,
78                           10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
79                           20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
80                           30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
81                           40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5,
82                           50.5, 55.5, 60.5, 65.5, 70.5, 75.5, 80.5, 85.5, 90.5, 95.5,
83                           100.5, 105.5, 110.5, 115.5, 120.5, 125.5, 130.5, 135.5, 140.5, 145.5,
84                           150.5, 160.5, 170.5, 180.5, 190.5, 200.5, 210.5, 220.5, 230.5, 240.5,
85                           250.5, 275.5, 300.5, 325.5, 350.5, 375.5, 400.5, 425.5, 450.5, 475.5,
86                           500.5 };
87                                 //525.5, 550.5, 575.5, 600.5, 625.5, 650.5, 675.5, 700.5, 725.5,
88                           //750.5, 775.5, 800.5, 825.5, 850.5, 875.5, 900.5, 925.5, 950.5, 975.5,
89                           //1000.5 };
90
91   #define BINNING fgMaxParams, binLimitsN
92
93   for (Int_t i = 0; i < kESDHists; ++i)
94     fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", 10, binLimitsVtx, BINNING);
95
96   for (Int_t i = 0; i < kMCHists; ++i)
97   {
98     fMultiplicityMC[i] = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("fMultiplicityMC%d", i)));
99     fMultiplicityMC[i]->SetTitle("fMultiplicityMC");
100   }
101
102   for (Int_t i = 0; i < kCorrHists; ++i)
103   {
104     fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;Npart;Ntracks", 10, binLimitsVtx, BINNING, BINNING);
105     fMultiplicityESDCorrected[i] = new TH1F(Form("fMultiplicityESDCorrected%d", i), "fMultiplicityESDCorrected;Npart;dN/dN", BINNING);
106   }
107
108   TH1::AddDirectory(oldStatus);
109 }
110
111 //____________________________________________________________________
112 AliMultiplicityCorrection::~AliMultiplicityCorrection()
113 {
114   //
115   // Destructor
116   //
117 }
118
119 //____________________________________________________________________
120 Long64_t AliMultiplicityCorrection::Merge(TCollection* list)
121 {
122   // Merge a list of AliMultiplicityCorrection objects with this (needed for
123   // PROOF).
124   // Returns the number of merged objects (including this).
125
126   if (!list)
127     return 0;
128
129   if (list->IsEmpty())
130     return 1;
131
132   TIterator* iter = list->MakeIterator();
133   TObject* obj;
134
135   // collections of all histograms
136   TList collections[kESDHists+kMCHists+kCorrHists*2];
137
138   Int_t count = 0;
139   while ((obj = iter->Next())) {
140
141     AliMultiplicityCorrection* entry = dynamic_cast<AliMultiplicityCorrection*> (obj);
142     if (entry == 0) 
143       continue;
144
145     for (Int_t i = 0; i < kESDHists; ++i)
146       collections[i].Add(entry->fMultiplicityESD[i]);
147
148     for (Int_t i = 0; i < kMCHists; ++i)
149       collections[kESDHists+i].Add(entry->fMultiplicityMC[i]);
150
151     for (Int_t i = 0; i < kCorrHists; ++i)
152       collections[kESDHists+kMCHists+i].Add(entry->fCorrelation[i]);
153
154     for (Int_t i = 0; i < kCorrHists; ++i)
155       collections[kESDHists+kMCHists+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]);
156
157     count++;
158   }
159
160   for (Int_t i = 0; i < kESDHists; ++i)
161     fMultiplicityESD[i]->Merge(&collections[i]);
162
163   for (Int_t i = 0; i < kMCHists; ++i)
164     fMultiplicityMC[i]->Merge(&collections[kESDHists+i]);
165
166   for (Int_t i = 0; i < kCorrHists; ++i)
167     fCorrelation[i]->Merge(&collections[kESDHists+kMCHists+i]);
168
169   for (Int_t i = 0; i < kCorrHists; ++i)
170     fMultiplicityESDCorrected[i]->Merge(&collections[kESDHists+kMCHists+kCorrHists+i]);
171
172   delete iter;
173
174   return count+1;
175 }
176
177 //____________________________________________________________________
178 Bool_t AliMultiplicityCorrection::LoadHistograms(const Char_t* dir)
179 {
180   //
181   // loads the histograms from a file
182   // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
183   //
184
185   if (!dir)
186     dir = GetName();
187
188   if (!gDirectory->cd(dir))
189     return kFALSE;
190
191   // TODO memory leak. old histograms needs to be deleted.
192
193   Bool_t success = kTRUE;
194
195   for (Int_t i = 0; i < kESDHists; ++i)
196   {
197     fMultiplicityESD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityESD[i]->GetName()));
198     if (!fMultiplicityESD[i])
199       success = kFALSE;
200   }
201
202   for (Int_t i = 0; i < kMCHists; ++i)
203   {
204     fMultiplicityMC[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityMC[i]->GetName()));
205     if (!fMultiplicityMC[i])
206       success = kFALSE;
207   }
208
209   for (Int_t i = 0; i < kCorrHists; ++i)
210   {
211     fCorrelation[i] = dynamic_cast<TH3F*> (gDirectory->Get(fCorrelation[i]->GetName()));
212     if (!fCorrelation[i])
213       success = kFALSE;
214     fMultiplicityESDCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fMultiplicityESDCorrected[i]->GetName()));
215     if (!fMultiplicityESDCorrected[i])
216       success = kFALSE;
217   }
218
219   gDirectory->cd("..");
220
221   return success;
222 }
223
224 //____________________________________________________________________
225 void AliMultiplicityCorrection::SaveHistograms()
226 {
227   //
228   // saves the histograms
229   //
230
231   gDirectory->mkdir(GetName());
232   gDirectory->cd(GetName());
233
234   for (Int_t i = 0; i < kESDHists; ++i)
235     if (fMultiplicityESD[i])
236       fMultiplicityESD[i]->Write();
237
238   for (Int_t i = 0; i < kMCHists; ++i)
239     if (fMultiplicityMC[i])
240       fMultiplicityMC[i]->Write();
241
242   for (Int_t i = 0; i < kCorrHists; ++i)
243   {
244     if (fCorrelation[i])
245       fCorrelation[i]->Write();
246     if (fMultiplicityESDCorrected[i])
247       fMultiplicityESDCorrected[i]->Write();
248   }
249
250   gDirectory->cd("..");
251 }
252
253 //____________________________________________________________________
254 void AliMultiplicityCorrection::FillGenerated(Float_t vtx, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll)
255 {
256   //
257   // Fills an event from MC
258   //
259
260   fMultiplicityMC[0]->Fill(vtx, generated05);
261   fMultiplicityMC[1]->Fill(vtx, generated10);
262   fMultiplicityMC[2]->Fill(vtx, generated15);
263   fMultiplicityMC[3]->Fill(vtx, generated20);
264   fMultiplicityMC[4]->Fill(vtx, generatedAll);
265 }
266
267 //____________________________________________________________________
268 void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20)
269 {
270   //
271   // Fills an event from ESD
272   //
273
274   fMultiplicityESD[0]->Fill(vtx, measured05);
275   fMultiplicityESD[1]->Fill(vtx, measured10);
276   fMultiplicityESD[2]->Fill(vtx, measured15);
277   fMultiplicityESD[3]->Fill(vtx, measured20);
278 }
279
280 //____________________________________________________________________
281 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)
282 {
283   //
284   // Fills an event into the correlation map with the information from MC and ESD
285   //
286
287   fCorrelation[0]->Fill(vtx, generated05, measured05);
288   fCorrelation[1]->Fill(vtx, generated10, measured10);
289   fCorrelation[2]->Fill(vtx, generated15, measured15);
290   fCorrelation[3]->Fill(vtx, generated20, measured20);
291
292   fCorrelation[4]->Fill(vtx, generatedAll, measured05);
293   fCorrelation[5]->Fill(vtx, generatedAll, measured10);
294   fCorrelation[6]->Fill(vtx, generatedAll, measured15);
295   fCorrelation[7]->Fill(vtx, generatedAll, measured20);
296 }
297
298 //____________________________________________________________________
299 Double_t AliMultiplicityCorrection::RegularizationPol0(Double_t *params)
300 {
301   // homogenity term for minuit fitting
302   // pure function of the parameters
303   // prefers constant function (pol0)
304
305   Double_t chi2 = 0;
306
307   for (Int_t i=1; i<fgMaxParams; ++i)
308   {
309     if (params[i] == 0)
310       continue;
311
312     Double_t right  = params[i] / fCurrentESD->GetBinWidth(i+1);
313     Double_t left   = params[i-1] / fCurrentESD->GetBinWidth(i);
314
315     Double_t diff = (right - left) / right;
316     chi2 += diff * diff;
317   }
318
319   return chi2;
320 }
321
322 //____________________________________________________________________
323 Double_t AliMultiplicityCorrection::RegularizationPol1(Double_t *params)
324 {
325   // homogenity term for minuit fitting
326   // pure function of the parameters
327   // prefers linear function (pol1)
328
329   Double_t chi2 = 0;
330
331   for (Int_t i=2; i<fgMaxParams; ++i)
332   {
333     if (params[i] == 0 || params[i-1] == 0)
334       continue;
335
336     Double_t right  = params[i] / fCurrentESD->GetBinWidth(i+1);
337     Double_t middle = params[i-1] / fCurrentESD->GetBinWidth(i);
338     Double_t left   = params[i-2] / fCurrentESD->GetBinWidth(i-1);
339
340     Double_t der1 = (right - middle) / fCurrentESD->GetBinWidth(i);
341     Double_t der2 = (middle - left)  / fCurrentESD->GetBinWidth(i-1);
342
343     Double_t diff = der1 - der2;
344
345     // TODO give additonal weight to big bins
346     chi2 += diff * diff * fCurrentESD->GetBinWidth(i) * fCurrentESD->GetBinWidth(i-1);
347
348     //printf("%d --> %f\n", i, diff);
349   }
350
351   return chi2 / 1e5 / 2;
352 }
353
354 //____________________________________________________________________
355 Double_t AliMultiplicityCorrection::RegularizationTotalCurvature(Double_t *params)
356 {
357   // homogenity term for minuit fitting
358   // pure function of the parameters
359   // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments,
360   // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp.
361
362   Double_t chi2 = 0;
363
364   for (Int_t i=2; i<fgMaxParams; ++i)
365   {
366     if (params[i] == 0 || params[i-1] == 0)
367       continue;
368
369     Double_t right  = params[i] / fCurrentESD->GetBinWidth(i+1);
370     Double_t middle = params[i-1] / fCurrentESD->GetBinWidth(i);
371     Double_t left   = params[i-2] / fCurrentESD->GetBinWidth(i-1);
372
373     Double_t der1 = (right - middle) / fCurrentESD->GetBinWidth(i);
374     Double_t der2 = (middle - left)  / fCurrentESD->GetBinWidth(i-1);
375
376     Double_t secDer = (der1 - der2) / fCurrentESD->GetBinWidth(i);
377
378     // square and weight with the bin width
379     chi2 += secDer * secDer * fCurrentESD->GetBinWidth(i) * fCurrentESD->GetBinWidth(i-1);
380
381     //printf("%d --> %f\n", i, secDer);
382   }
383
384   return chi2 / 1e5;
385 }
386
387 //____________________________________________________________________
388 void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t)
389 {
390   //
391   // fit function for minuit
392   //
393
394   // TODO take errors into account
395
396   static Int_t callCount = 0;
397
398   Double_t chi2FromFit = 0;
399
400   // loop over multiplicity (x axis of fMultiplicityESD)
401   for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i)
402   {
403     if (i > fCurrentCorrelation->GetNbinsY())
404       break;
405
406     Double_t sum = 0;
407     //Double_t error = 0;
408     // loop over generated particles (x axis of fCorrelation) that resulted in reconstruction of i tracks
409     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsX(); ++j)
410     {
411       if (j > fgMaxParams)
412         break;
413
414       sum += fCurrentCorrelation->GetBinContent(j, i) * params[j-1];
415
416       //if (params[j-1] > 0)
417       //  error += fCurrentCorrelation->GetBinError(j, i) * fCurrentCorrelation->GetBinError(j, i) * params[j-1];
418       //printf("%f  ", sum);
419     }
420
421     Double_t diff = fCurrentESD->GetBinContent(i) - sum;
422     if (fCurrentESD->GetBinContent(i) > 0)
423       diff /= fCurrentESD->GetBinContent(i);
424     else
425       diff /= fCurrentESD->Integral();
426
427     // weight with bin width
428     //diff *= fCurrentESD->GetBinWidth(i);
429
430     //error = TMath::Sqrt(error) + fCurrentESD->GetBinError(i);
431     //if (error <= 0)
432     // error = 1;
433
434     //Double_t tmp = diff / error;
435     //chi2 += tmp * tmp;
436     chi2FromFit += diff * diff;
437
438     //printf("\nExpected sum = %f; Diff for bin %d is %f\n**********************************\n", fCurrentESD->GetBinContent(i), i, diff);
439     //printf("Diff for bin %d is %f\n", i, diff);
440   }
441
442   Double_t penaltyVal = RegularizationTotalCurvature(params);
443
444   chi2 = chi2FromFit * chi2FromFit + penaltyVal * penaltyVal;
445
446   if ((callCount++ % 100) == 0)
447     printf("%f %f --> %f\n", chi2FromFit, penaltyVal, chi2);
448 }
449
450 //____________________________________________________________________
451 void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace)
452 {
453   //
454   // fills fCurrentESD, fCurrentCorrelation
455   // resets fMultiplicityESDCorrected
456   //
457
458   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
459   fMultiplicityESDCorrected[correlationID]->Reset();
460
461   fCurrentESD = fMultiplicityESD[inputRange]->ProjectionY();
462   fCurrentESD->Sumw2();
463
464   // empty under/overflow bins in x, otherwise Project3D takes them into account
465   TH3* hist = fCorrelation[correlationID];
466   for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
467   {
468     for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
469     {
470       hist->SetBinContent(0, y, z, 0);
471       hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
472     }
473   }
474
475   fCurrentCorrelation = hist->Project3D("zy");
476   fCurrentCorrelation->Sumw2();
477 }
478
479 //____________________________________________________________________
480 void AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, Bool_t check)
481 {
482   //
483   // correct spectrum using minuit chi2 method
484   //
485   // if check is kTRUE the input MC solution (by definition the right solution) is used
486   // no fit is made and just the chi2 is calculated
487   //
488
489   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
490   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
491
492   SetupCurrentHists(inputRange, fullPhaseSpace);
493
494   // normalize correction for given nPart
495   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
496   {
497     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
498     if (sum <= 0)
499       continue;
500     Float_t maxValue = 0;
501     Int_t maxBin = -1;
502     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
503     {
504       // find most probably value
505       if (maxValue < fCurrentCorrelation->GetBinContent(i, j))
506       {
507         maxValue = fCurrentCorrelation->GetBinContent(i, j);
508         maxBin = j;
509       }
510
511       // npart sum to 1
512       fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
513       fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
514     }
515
516     printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
517   }
518
519   // small hack to get around charge conservation for full phase space ;-)
520   /*if (fullPhaseSpace)
521   {
522     for (Int_t i=2; i<=50; i+=2)
523       for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
524         fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i-1, j));
525   }*/
526
527   TCanvas* canvas1 = new TCanvas("ApplyMinuitFit", "ApplyMinuitFit", 800, 400);
528   canvas1->Divide(2, 1);
529   canvas1->cd(1);
530   fCurrentESD->DrawCopy();
531   canvas1->cd(2);
532   fCurrentCorrelation->DrawCopy("COLZ");
533
534   // Initialize TMinuit via generic fitter interface
535   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgMaxParams);
536
537   minuit->SetFCN(MinuitFitFunction);
538
539   Double_t results[fgMaxParams];
540
541   TH1* proj = fMultiplicityMC[mcTarget]->ProjectionY();
542   for (Int_t i=0; i<fgMaxParams; ++i)
543   {
544     //results[i] = 1;
545     results[i] = fCurrentESD->GetBinContent(i+1);
546     if (check)
547       results[i] = proj->GetBinContent(i+1);
548     minuit->SetParameter(i, Form("param%d", i), results[i], ((results[i] > 1) ? (results[i] / 10) : 0), 0, fCurrentESD->GetMaximum() * 100);
549   }
550
551   Int_t dummy;
552   Double_t chi2;
553   MinuitFitFunction(dummy, 0, chi2, results, 0);
554   printf("Chi2 of right solution is = %f\n", chi2);
555
556   if (check)
557     return;
558
559   Double_t arglist[100];
560   arglist[0] = 100000;
561   //minuit->ExecuteCommand("SET PRINT", arglist, 1);
562   minuit->ExecuteCommand("MIGRAD", arglist, 1);
563   //minuit->ExecuteCommand("MINOS", arglist, 0);
564
565   for (Int_t i=0; i<fgMaxParams; ++i)
566   {
567     results[i] = minuit->GetParameter(i);
568     fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, results[i]);
569     fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0);
570   }
571
572   printf("Penalty is %f\n", RegularizationTotalCurvature(results));
573
574   DrawComparison("MinuitChi2", mcTarget, correlationID);
575 }
576
577 //____________________________________________________________________
578 void AliMultiplicityCorrection::NormalizeToBinWidth(TH1* hist)
579 {
580   //
581   // normalizes a 1-d histogram to its bin width
582   //
583
584   for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
585   {
586     hist->SetBinContent(i, hist->GetBinContent(i) / hist->GetBinWidth(i));
587     hist->SetBinError(i, hist->GetBinError(i) / hist->GetBinWidth(i));
588   }
589 }
590
591 //____________________________________________________________________
592 void AliMultiplicityCorrection::NormalizeToBinWidth(TH2* hist)
593 {
594   //
595   // normalizes a 2-d histogram to its bin width (x width * y width)
596   //
597
598   for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
599     for (Int_t j=1; j<=hist->GetNbinsY(); ++j)
600     {
601       Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j);
602       hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor);
603       hist->SetBinError(i, j, hist->GetBinError(i, j) / factor);
604     }
605 }
606
607 //____________________________________________________________________
608 void AliMultiplicityCorrection::ApplyMinuitFitAll()
609 {
610   //
611   // fit all eta ranges
612   //
613
614   for (Int_t i=0; i<kESDHists; ++i)
615   {
616     ApplyMinuitFit(i, kFALSE);
617     ApplyMinuitFit(i, kTRUE);
618   }
619 }
620
621 //____________________________________________________________________
622 void AliMultiplicityCorrection::DrawHistograms()
623 {
624   printf("ESD:\n");
625
626   TCanvas* canvas1 = new TCanvas("fMultiplicityESD", "fMultiplicityESD", 900, 600);
627   canvas1->Divide(3, 2);
628   for (Int_t i = 0; i < kESDHists; ++i)
629   {
630     canvas1->cd(i+1);
631     fMultiplicityESD[i]->DrawCopy("COLZ");
632     printf("%d --> %f\n", i, (Float_t) fMultiplicityESD[i]->ProjectionY()->GetMean());
633   }
634
635   printf("MC:\n");
636
637   TCanvas* canvas2 = new TCanvas("fMultiplicityMC", "fMultiplicityMC", 900, 600);
638   canvas2->Divide(3, 2);
639   for (Int_t i = 0; i < kMCHists; ++i)
640   {
641     canvas2->cd(i+1);
642     fMultiplicityMC[i]->DrawCopy("COLZ");
643     printf("%d --> %f\n", i, fMultiplicityMC[i]->ProjectionY()->GetMean());
644   }
645
646   TCanvas* canvas3 = new TCanvas("fCorrelation", "fCorrelation", 900, 900);
647   canvas3->Divide(3, 3);
648   for (Int_t i = 0; i < kCorrHists; ++i)
649   {
650     canvas3->cd(i+1);
651     TH1* proj = fCorrelation[i]->Project3D("zy");
652     proj->DrawCopy("COLZ");
653   }
654 }
655
656 //____________________________________________________________________
657 void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t mcID, Int_t esdCorrId, Bool_t normalizeESD)
658 {
659   TString tmpStr;
660   tmpStr.Form("%s_DrawComparison_%d_%d", name, mcID, esdCorrId);
661
662   TCanvas* canvas1 = new TCanvas(tmpStr, tmpStr, 900, 300);
663   canvas1->Divide(3, 1);
664
665   canvas1->cd(1);
666   TH1* proj = fMultiplicityMC[mcID]->ProjectionY();
667   NormalizeToBinWidth(proj);
668
669   if (normalizeESD)
670     NormalizeToBinWidth(fMultiplicityESDCorrected[esdCorrId]);
671
672   proj->GetXaxis()->SetRangeUser(0, 200);
673   proj->DrawCopy("HIST");
674   gPad->SetLogy();
675
676   NormalizeToBinWidth(fCurrentESD);
677   fCurrentESD->SetLineColor(2);
678   fCurrentESD->DrawCopy("HISTSAME");
679
680   fMultiplicityESDCorrected[esdCorrId]->SetMarkerStyle(3);
681   fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAMEP");
682
683   canvas1->cd(2);
684   fMultiplicityESDCorrected[esdCorrId]->GetXaxis()->SetRangeUser(0, 200);
685   fMultiplicityESDCorrected[esdCorrId]->DrawCopy("HIST");
686
687   canvas1->cd(3);
688   TH1* clone = dynamic_cast<TH1*> (proj->Clone("clone"));
689   clone->Divide(fMultiplicityESDCorrected[esdCorrId]);
690   clone->SetTitle("Ratio;Npart;MC/ESD");
691   clone->GetYaxis()->SetRangeUser(0.8, 1.2);
692   clone->GetXaxis()->SetRangeUser(0, 200);
693   clone->DrawCopy("HIST");
694
695   /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
696   legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
697   legend->AddEntry(fMultiplicityMC, "MC");
698   legend->AddEntry(fMultiplicityESD, "ESD");
699   legend->Draw();*/
700 }
701
702 //____________________________________________________________________
703 void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace)
704 {
705   //
706   // correct spectrum using bayesian method
707   //
708
709   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
710   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
711
712   SetupCurrentHists(inputRange, fullPhaseSpace);
713
714   // normalize correction for given nPart
715   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
716   {
717     Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
718     if (sum <= 0)
719       continue;
720     for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
721     {
722       // npart sum to 1
723       fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
724       fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
725     }
726   }
727
728   Float_t regPar = 0.01;
729
730   Float_t norm = 0;
731   // pick prior distribution
732   TH1F* hPrior = (TH1F*)fCurrentESD->Clone("prior");
733   for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
734     norm = norm + hPrior->GetBinContent(t);
735   for (Int_t t=1; t<=hPrior->GetNbinsX(); t++) {
736     hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm);
737
738     printf(" bin %d content %f \n", t, hPrior->GetBinContent(t));
739
740   }
741   
742   // define hist to store guess of true
743   TH1F* hTrue  = (TH1F*)fCurrentESD->Clone("prior");
744   //  hTrue->Reset();
745
746   // clone response R
747   TH2F* hResponse = (TH2F*)fCurrentCorrelation->Clone("pij");
748
749   // histogram to store the inverse response calculated from bayes theorem (from prior and response)
750   // IR
751   //TAxis* xAxis = hResponse->GetYaxis();
752   //TAxis* yAxis = hResponse->GetXaxis();
753
754   TH2F* hInverseResponseBayes = (TH2F*)hResponse->Clone("pji");
755   //new TH2F("pji","pji",
756   //                                     xAxis->GetNbins(),xAxis->GetXbins()->GetArray(),
757   //                                     yAxis->GetNbins(),yAxis->GetXbins()->GetArray());
758   hInverseResponseBayes->Reset();
759   
760   // unfold...
761   Int_t iterations = 50;
762   for (Int_t i=0; i<iterations; i++) {
763     printf(" iteration %i \n", i);
764     
765     // calculate IR from Beyes theorem
766     // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
767     for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++) {
768       for (Int_t m=1; m<=hResponse->GetNbinsY(); m++) {
769
770         norm = 0;
771         for(Int_t t2 = 1; t2<=hResponse->GetNbinsX(); t2++)
772           norm += (hResponse->GetBinContent(t2,m) * hPrior->GetBinContent(t2));
773
774         if (norm==0)
775           hInverseResponseBayes->SetBinContent(t,m,0);
776         else
777           hInverseResponseBayes->SetBinContent(t,m, hResponse->GetBinContent(t,m) * hPrior->GetBinContent(t)/norm);
778       }
779     }
780     // calculate true assuming IR is the correct inverse response
781     for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++) {
782       hTrue ->SetBinContent(t, 0);
783       for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
784         hTrue->SetBinContent(t, hTrue->GetBinContent(t) + fCurrentESD->GetBinContent(m)*hInverseResponseBayes->GetBinContent(t,m));
785     }
786
787     // regularization (simple smoothing) NB : not good bin width!!!
788     TH1F* hTrueSmoothed = (TH1F*)hTrue->Clone("truesmoothed");
789
790     for (Int_t t=2; t<hTrue->GetNbinsX()-1; t++) {
791       Float_t average = (hTrue->GetBinContent(t-1) / hTrue->GetBinWidth(t-1)
792                          + hTrue->GetBinContent(t) / hTrue->GetBinWidth(t)
793                          + hTrue->GetBinContent(t+1) / hTrue->GetBinWidth(t+1)) / 3.;
794       average *= hTrueSmoothed->GetBinWidth(t);
795
796       // weight the average with the regularization parameter
797       hTrueSmoothed->SetBinContent(t, (1-regPar)*hTrue->GetBinContent(t) + (regPar*average));
798     }
799
800     // use smoothed true (normalized) as new prior
801     norm = 0;
802     for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
803       norm = norm + hTrueSmoothed->GetBinContent(t);
804     for (Int_t t=1; t<=hPrior->GetNbinsX(); t++) {
805       hPrior->SetBinContent(t, hTrueSmoothed->GetBinContent(t)/norm);
806       hTrue->SetBinContent(t, hTrueSmoothed->GetBinContent(t));
807     }
808
809     delete hTrueSmoothed;
810   } // end of iterations
811
812   for (Int_t t=1; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++) {
813     fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTrue->GetBinContent(t));
814     fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05*hTrue->GetBinContent(t));
815
816     printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
817   }
818   
819   DrawComparison("Bayesian", mcTarget, correlationID);
820 }
821
822 //____________________________________________________________________
823 void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace)
824 {
825   //
826   // correct spectrum using a simple Gaussian approach, that is model-dependent
827   //
828
829   Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
830   Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
831
832   SetupCurrentHists(inputRange, fullPhaseSpace);
833
834   NormalizeToBinWidth((TH2*) fCurrentCorrelation);
835
836   TH1D* correction = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianMean"));
837   correction->SetTitle("GaussianMean");
838
839   TH1D* correctionWidth = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianWidth"));
840   correction->SetTitle("GaussianWidth");
841
842   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
843   {
844     TH1D* proj = (dynamic_cast<TH2*> (fCurrentCorrelation))->ProjectionX("_px", i, i+1);
845     proj->Fit("gaus", "0Q");
846     correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
847     correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
848
849     continue;
850
851     // draw for debugging
852     new TCanvas;
853     proj->DrawCopy();
854     proj->GetFunction("gaus")->DrawCopy("SAME");
855   }
856
857   TH1* target = fMultiplicityESDCorrected[correlationID];
858
859   Int_t nBins = target->GetNbinsX()*10+1;
860   Float_t* binning = new Float_t[nBins];
861   for (Int_t i=1; i<=target->GetNbinsX(); ++i)
862     for (Int_t j=0; j<10; ++j)
863       binning[(i-1)*10 + j] = target->GetXaxis()->GetBinLowEdge(i) + target->GetXaxis()->GetBinWidth(i) / 10 * j;
864
865   binning[nBins-1] = target->GetXaxis()->GetBinUpEdge(target->GetNbinsX());
866
867   TH1F* fineBinned = new TH1F("targetFineBinned", "targetFineBinned", nBins-1, binning);
868
869   for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
870   {
871     Float_t mean = correction->GetBinContent(i);
872     Float_t width = correctionWidth->GetBinContent(i);
873
874     Int_t fillBegin = fineBinned->FindBin(mean - width * 3);
875     Int_t fillEnd   = fineBinned->FindBin(mean + width * 3);
876     printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
877
878     for (Int_t j=fillBegin; j <= fillEnd; ++j)
879     {
880       fineBinned->AddBinContent(j, TMath::Gaus(fineBinned->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fCurrentESD->GetBinContent(i));
881     }
882   }
883
884   for (Int_t i=1; i<=target->GetNbinsX(); ++i)
885   {
886     Float_t sum = 0;
887     for (Int_t j=1; j<=10; ++j)
888       sum += fineBinned->GetBinContent((i-1)*10 + j);
889     target->SetBinContent(i, sum / 10);
890   }
891
892   delete[] binning;
893
894   DrawComparison("Gaussian", mcTarget, correlationID, kFALSE);
895 }