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