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