Adding draw and Print function to AliESDtrackCuts
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AliMultiplicityCorrection.cxx
CommitLineData
0a173978 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>
9ca6be09 34#include <TF1.h>
d560b581 35#include <TMath.h>
3602328d 36#include <TCollection.h>
0a173978 37
38ClassImp(AliMultiplicityCorrection)
39
b477f4f2 40const Int_t AliMultiplicityCorrection::fgMaxParams = 251;
9ca6be09 41TH1* AliMultiplicityCorrection::fCurrentESD = 0;
42TH1* AliMultiplicityCorrection::fCurrentCorrelation = 0;
0a173978 43
44//____________________________________________________________________
45AliMultiplicityCorrection::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//____________________________________________________________________
66AliMultiplicityCorrection::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
b477f4f2 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};
0a173978 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,
9ca6be09 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 };
0a173978 91
b477f4f2 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
0a173978 97
98 for (Int_t i = 0; i < kESDHists; ++i)
b477f4f2 99 fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", VTXBINNING, NBINNING);
0a173978 100
101 for (Int_t i = 0; i < kMCHists; ++i)
102 {
103 fMultiplicityMC[i] = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("fMultiplicityMC%d", i)));
b477f4f2 104 fMultiplicityMC[i]->SetTitle("fMultiplicityMC;vtx-z;Npart");
0a173978 105 }
106
107 for (Int_t i = 0; i < kCorrHists; ++i)
108 {
b477f4f2 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);
0a173978 111 }
112
113 TH1::AddDirectory(oldStatus);
114}
115
116//____________________________________________________________________
117AliMultiplicityCorrection::~AliMultiplicityCorrection()
118{
119 //
120 // Destructor
121 //
122}
123
124//____________________________________________________________________
125Long64_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//____________________________________________________________________
183Bool_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//____________________________________________________________________
230void 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//____________________________________________________________________
259void 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//____________________________________________________________________
273void 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//____________________________________________________________________
286void 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//____________________________________________________________________
9ca6be09 304Double_t AliMultiplicityCorrection::RegularizationPol0(Double_t *params)
0a173978 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
9ca6be09 317 Double_t right = params[i] / fCurrentESD->GetBinWidth(i+1);
318 Double_t left = params[i-1] / fCurrentESD->GetBinWidth(i);
0a173978 319
320 Double_t diff = (right - left) / right;
321 chi2 += diff * diff;
322 }
323
324 return chi2;
325}
326
327//____________________________________________________________________
9ca6be09 328Double_t AliMultiplicityCorrection::RegularizationPol1(Double_t *params)
0a173978 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
9ca6be09 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);
0a173978 344
b477f4f2 345 Double_t der1 = (right - middle) / middle / fCurrentESD->GetBinWidth(i);
346 Double_t der2 = (middle - left) / middle / fCurrentESD->GetBinWidth(i-1);
0a173978 347
3602328d 348 Double_t diff = (der1 - der2); /// fCurrentESD->GetBinWidth(i);
0a173978 349
350 // TODO give additonal weight to big bins
9ca6be09 351 chi2 += diff * diff * fCurrentESD->GetBinWidth(i) * fCurrentESD->GetBinWidth(i-1);
0a173978 352
353 //printf("%d --> %f\n", i, diff);
354 }
355
b477f4f2 356 return chi2 * 1000;
0a173978 357}
358
9ca6be09 359//____________________________________________________________________
360Double_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
3602328d 389 return chi2 * 10;
390}
391
392//____________________________________________________________________
393Double_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;
9ca6be09 425}
426
0a173978 427//____________________________________________________________________
428void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t)
429{
430 //
431 // fit function for minuit
432 //
433
0a173978 434 static Int_t callCount = 0;
435
436 Double_t chi2FromFit = 0;
437
438 // loop over multiplicity (x axis of fMultiplicityESD)
9ca6be09 439 for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i)
0a173978 440 {
9ca6be09 441 if (i > fCurrentCorrelation->GetNbinsY())
0a173978 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
9ca6be09 447 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsX(); ++j)
0a173978 448 {
449 if (j > fgMaxParams)
450 break;
451
9ca6be09 452 sum += fCurrentCorrelation->GetBinContent(j, i) * params[j-1];
0a173978 453
454 //if (params[j-1] > 0)
9ca6be09 455 // error += fCurrentCorrelation->GetBinError(j, i) * fCurrentCorrelation->GetBinError(j, i) * params[j-1];
0a173978 456 //printf("%f ", sum);
457 }
458
9ca6be09 459 Double_t diff = fCurrentESD->GetBinContent(i) - sum;
3602328d 460
461 // include error
462 if (fCurrentESD->GetBinError(i) > 0)
463 diff /= fCurrentESD->GetBinError(i);
464
465 /*if (fCurrentESD->GetBinContent(i) > 0)
9ca6be09 466 diff /= fCurrentESD->GetBinContent(i);
0a173978 467 else
3602328d 468 diff /= fCurrentESD->Integral();*/
9ca6be09 469
470 // weight with bin width
3602328d 471 diff *= fCurrentESD->GetBinWidth(i);
9ca6be09 472
473 //error = TMath::Sqrt(error) + fCurrentESD->GetBinError(i);
0a173978 474 //if (error <= 0)
475 // error = 1;
476
477 //Double_t tmp = diff / error;
478 //chi2 += tmp * tmp;
479 chi2FromFit += diff * diff;
480
9ca6be09 481 //printf("\nExpected sum = %f; Diff for bin %d is %f\n**********************************\n", fCurrentESD->GetBinContent(i), i, diff);
0a173978 482 //printf("Diff for bin %d is %f\n", i, diff);
483 }
484
b477f4f2 485 Double_t penaltyVal = RegularizationPol1(params);
0a173978 486
3602328d 487 chi2 = chi2FromFit + penaltyVal;
488 //chi2 = penaltyVal;
0a173978 489
490 if ((callCount++ % 100) == 0)
491 printf("%f %f --> %f\n", chi2FromFit, penaltyVal, chi2);
492}
493
494//____________________________________________________________________
9ca6be09 495void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace)
0a173978 496{
497 //
9ca6be09 498 // fills fCurrentESD, fCurrentCorrelation
499 // resets fMultiplicityESDCorrected
0a173978 500 //
501
502 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
9ca6be09 503 fMultiplicityESDCorrected[correlationID]->Reset();
0a173978 504
9ca6be09 505 fCurrentESD = fMultiplicityESD[inputRange]->ProjectionY();
506 fCurrentESD->Sumw2();
0a173978 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
9ca6be09 519 fCurrentCorrelation = hist->Project3D("zy");
520 fCurrentCorrelation->Sumw2();
521}
522
523//____________________________________________________________________
524void 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);
0a173978 537
538 // normalize correction for given nPart
9ca6be09 539 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
0a173978 540 {
9ca6be09 541 Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
0a173978 542 if (sum <= 0)
543 continue;
9ca6be09 544 Float_t maxValue = 0;
545 Int_t maxBin = -1;
546 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
0a173978 547 {
9ca6be09 548 // find most probably value
549 if (maxValue < fCurrentCorrelation->GetBinContent(i, j))
550 {
551 maxValue = fCurrentCorrelation->GetBinContent(i, j);
552 maxBin = j;
553 }
554
0a173978 555 // npart sum to 1
9ca6be09 556 fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
557 fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
0a173978 558 }
9ca6be09 559
560 printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
0a173978 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)
9ca6be09 567 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
568 fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i-1, j));
0a173978 569 }*/
570
571 TCanvas* canvas1 = new TCanvas("ApplyMinuitFit", "ApplyMinuitFit", 800, 400);
572 canvas1->Divide(2, 1);
573 canvas1->cd(1);
9ca6be09 574 fCurrentESD->DrawCopy();
0a173978 575 canvas1->cd(2);
9ca6be09 576 fCurrentCorrelation->DrawCopy("COLZ");
0a173978 577
578 // Initialize TMinuit via generic fitter interface
579 TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgMaxParams);
580
581 minuit->SetFCN(MinuitFitFunction);
582
3602328d 583 static Double_t results[fgMaxParams];
584 printf("%x\n", (void*) results);
0a173978 585
9ca6be09 586 TH1* proj = fMultiplicityMC[mcTarget]->ProjectionY();
0a173978 587 for (Int_t i=0; i<fgMaxParams; ++i)
588 {
589 //results[i] = 1;
9ca6be09 590 results[i] = fCurrentESD->GetBinContent(i+1);
b477f4f2 591 if (results[i] < 1e-2)
592 results[i] = 1e-2;
9ca6be09 593 if (check)
594 results[i] = proj->GetBinContent(i+1);
b477f4f2 595 minuit->SetParameter(i, Form("param%d", i), results[i], results[i] / 10, 0.001, fCurrentESD->GetMaximum() * 10);
3602328d 596 //minuit->SetParameter(i, Form("param%d", i), fCurrentESD->GetBinWidth(i+1), 0.01, 0.001, 50);
0a173978 597 }
3602328d 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);
0a173978 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
9ca6be09 606 if (check)
607 return;
0a173978 608
609 Double_t arglist[100];
3602328d 610 arglist[0] = 0;
611 minuit->ExecuteCommand("SET PRINT", arglist, 1);
0a173978 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
3602328d 622 //printf("Penalty is %f\n", RegularizationTotalCurvature(results));
0a173978 623
9ca6be09 624 DrawComparison("MinuitChi2", mcTarget, correlationID);
0a173978 625}
626
627//____________________________________________________________________
628void 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//____________________________________________________________________
642void 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//____________________________________________________________________
658void 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//____________________________________________________________________
672void 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);
b477f4f2 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");
0a173978 711 proj->DrawCopy("COLZ");
712 }
713}
714
715//____________________________________________________________________
9ca6be09 716void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t mcID, Int_t esdCorrId, Bool_t normalizeESD)
0a173978 717{
9ca6be09 718 TString tmpStr;
719 tmpStr.Form("%s_DrawComparison_%d_%d", name, mcID, esdCorrId);
0a173978 720
b477f4f2 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);
0a173978 743
744 canvas1->cd(1);
745 TH1* proj = fMultiplicityMC[mcID]->ProjectionY();
746 NormalizeToBinWidth(proj);
9ca6be09 747
748 if (normalizeESD)
749 NormalizeToBinWidth(fMultiplicityESDCorrected[esdCorrId]);
750
0a173978 751 proj->GetXaxis()->SetRangeUser(0, 200);
b477f4f2 752 proj->DrawCopy("HIST E");
0a173978 753 gPad->SetLogy();
754
9ca6be09 755 NormalizeToBinWidth(fCurrentESD);
756 fCurrentESD->SetLineColor(2);
b477f4f2 757 fCurrentESD->DrawCopy("HIST SAME E");
758
759 proj2->SetLineColor(2);
760 proj2->SetMarkerColor(2);
761 proj2->SetMarkerStyle(5);
762 proj2->DrawCopy("P SAME");
9ca6be09 763
0a173978 764 fMultiplicityESDCorrected[esdCorrId]->SetMarkerStyle(3);
b477f4f2 765 fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAME P E");
0a173978 766
767 canvas1->cd(2);
9ca6be09 768 fMultiplicityESDCorrected[esdCorrId]->GetXaxis()->SetRangeUser(0, 200);
769 fMultiplicityESDCorrected[esdCorrId]->DrawCopy("HIST");
0a173978 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();*/
3602328d 784
b477f4f2 785 canvas1->cd(4);
786
787 residual->GetYaxis()->SetRangeUser(-0.2, 0.2);
788 residual->GetXaxis()->SetRangeUser(0, 200);
789 residual->DrawCopy();
3602328d 790
b477f4f2 791 canvas1->SaveAs(Form("%s.gif", canvas1->GetName()));
0a173978 792}
793
794//____________________________________________________________________
795void 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
9ca6be09 804 SetupCurrentHists(inputRange, fullPhaseSpace);
0a173978 805
806 // normalize correction for given nPart
9ca6be09 807 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
0a173978 808 {
9ca6be09 809 Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
0a173978 810 if (sum <= 0)
811 continue;
9ca6be09 812 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
0a173978 813 {
814 // npart sum to 1
9ca6be09 815 fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
816 fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
0a173978 817 }
818 }
819
6127aca6 820 Float_t regPar = 0.01;
821
822 Float_t norm = 0;
823 // pick prior distribution
9ca6be09 824 TH1F* hPrior = (TH1F*)fCurrentESD->Clone("prior");
6127aca6 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
9ca6be09 835 TH1F* hTrue = (TH1F*)fCurrentESD->Clone("prior");
6127aca6 836 // hTrue->Reset();
837
838 // clone response R
9ca6be09 839 TH2F* hResponse = (TH2F*)fCurrentCorrelation->Clone("pij");
6127aca6 840
9ca6be09 841 // histogram to store the inverse response calculated from bayes theorem (from prior and response)
6127aca6 842 // IR
9ca6be09 843 //TAxis* xAxis = hResponse->GetYaxis();
844 //TAxis* yAxis = hResponse->GetXaxis();
6127aca6 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();
b477f4f2 851
6127aca6 852 // unfold...
3602328d 853 Int_t iterations = 20;
6127aca6 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)
9ca6be09 859 for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++) {
6127aca6 860 for (Int_t m=1; m<=hResponse->GetNbinsY(); m++) {
861
862 norm = 0;
863 for(Int_t t2 = 1; t2<=hResponse->GetNbinsX(); t2++)
9ca6be09 864 norm += (hResponse->GetBinContent(t2,m) * hPrior->GetBinContent(t2));
865
6127aca6 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 }
9ca6be09 872 // calculate true assuming IR is the correct inverse response
873 for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++) {
6127aca6 874 hTrue ->SetBinContent(t, 0);
875 for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
9ca6be09 876 hTrue->SetBinContent(t, hTrue->GetBinContent(t) + fCurrentESD->GetBinContent(m)*hInverseResponseBayes->GetBinContent(t,m));
6127aca6 877 }
9ca6be09 878
b477f4f2 879 // regularization (simple smoothing)
6127aca6 880 TH1F* hTrueSmoothed = (TH1F*)hTrue->Clone("truesmoothed");
9ca6be09 881
6127aca6 882 for (Int_t t=2; t<hTrue->GetNbinsX()-1; t++) {
9ca6be09 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
6127aca6 888 // weight the average with the regularization parameter
889 hTrueSmoothed->SetBinContent(t, (1-regPar)*hTrue->GetBinContent(t) + (regPar*average));
890 }
9ca6be09 891
892 // use smoothed true (normalized) as new prior
6127aca6 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 }
9ca6be09 900
901 delete hTrueSmoothed;
6127aca6 902 } // end of iterations
903
9ca6be09 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));
6127aca6 907
9ca6be09 908 printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
6127aca6 909 }
910
9ca6be09 911 DrawComparison("Bayesian", mcTarget, correlationID);
912}
913
914//____________________________________________________________________
915void 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
3602328d 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);
9ca6be09 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);
0a173978 987}
3602328d 988
989//____________________________________________________________________
b477f4f2 990TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap, Bool_t normalized)
3602328d 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
b477f4f2 995 // if normalized is set, inputMC is expected to be normalized to the bin width
3602328d 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
b477f4f2 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);
3602328d 1050 }
1051
1052 // bin width of the <measured> axis has to be taken into account
b477f4f2 1053 //measured /= target->GetYaxis()->GetBinWidth(meas);
1054 //error /= target->GetYaxis()->GetBinWidth(meas);
3602328d 1055
b477f4f2 1056 printf("%f +- %f ; %f +- %f \n", inputMC->GetBinContent(meas), inputMC->GetBinError(meas), measured, error);
3602328d 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//____________________________________________________________________
1066void 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)
b477f4f2 1085 {
3602328d 1086 mc->SetBinContent(mc->GetNbinsX() / 2, i, inputMC->Eval(mc->GetYaxis()->GetBinCenter(i)) * mc->GetYaxis()->GetBinWidth(i));
b477f4f2 1087 mc->SetBinError(mc->GetNbinsX() / 2, i, 0);
1088 }
3602328d 1089
1090 //new TCanvas;
1091 //mc->Draw("COLZ");
1092
1093 TH1* proj = mc->ProjectionY();
1094
1095 fMultiplicityESD[id] = CalculateMultiplicityESD(proj, id);
1096}