git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Minor fixes (Andrei)
[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>
34
35ClassImp(AliMultiplicityCorrection)
36
37const Int_t AliMultiplicityCorrection::fgMaxParams = 110;
38TH1* AliMultiplicityCorrection::fCurrentMinuitESD = 0;
39TH1* AliMultiplicityCorrection::fCurrentMinuitCorrelation = 0;
40
41//____________________________________________________________________
42AliMultiplicityCorrection::AliMultiplicityCorrection() :
43 TNamed()
44{
45 //
46 // default constructor
47 //
48
49 for (Int_t i = 0; i < kESDHists; ++i)
50 fMultiplicityESD[i] = 0;
51
52 for (Int_t i = 0; i < kMCHists; ++i)
53 fMultiplicityMC[i] = 0;
54
55 for (Int_t i = 0; i < kCorrHists; ++i)
56 {
57 fCorrelation[i] = 0;
58 fMultiplicityESDCorrected[i] = 0;
59 }
60}
61
62//____________________________________________________________________
63AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const Char_t* title) :
64 TNamed(name, title)
65{
66 //
67 // named constructor
68 //
69
70 // do not add this hists to the directory
71 Bool_t oldStatus = TH1::AddDirectoryStatus();
72 TH1::AddDirectory(kFALSE);
73
74 Float_t binLimitsVtx[] = {-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10};
75 Float_t binLimitsN[] = {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
76 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
77 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
78 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
79 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5,
80 50.5, 55.5, 60.5, 65.5, 70.5, 75.5, 80.5, 85.5, 90.5, 95.5,
81 100.5, 105.5, 110.5, 115.5, 120.5, 125.5, 130.5, 135.5, 140.5, 145.5,
82 150.5, 160.5, 170.5, 180.5, 190.5, 200.5, 210.5, 220.5, 230.5, 240.5,
83 250.5, 275.5, 300.5, 325.5, 350.5, 375.5, 400.5, 425.5, 450.5, 475.5,
84 500.5, 525.5, 550.5, 575.5, 600.5, 625.5, 650.5, 675.5, 700.5, 725.5,
85 750.5, 775.5, 800.5, 825.5, 850.5, 875.5, 900.5, 925.5, 950.5, 975.5,
86 1000.5 };
87
88 #define BINNING 110, binLimitsN
89
90 for (Int_t i = 0; i < kESDHists; ++i)
91 fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", 10, binLimitsVtx, BINNING);
92
93 for (Int_t i = 0; i < kMCHists; ++i)
94 {
95 fMultiplicityMC[i] = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("fMultiplicityMC%d", i)));
96 fMultiplicityMC[i]->SetTitle("fMultiplicityMC");
97 }
98
99 for (Int_t i = 0; i < kCorrHists; ++i)
100 {
101 fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;Npart;Ntracks", 10, binLimitsVtx, BINNING, BINNING);
102 fMultiplicityESDCorrected[i] = new TH1F(Form("fMultiplicityESDCorrected%d", i), "fMultiplicityESDCorrected;Npart;dN/dN", BINNING);
103 }
104
105 TH1::AddDirectory(oldStatus);
106}
107
108//____________________________________________________________________
109AliMultiplicityCorrection::~AliMultiplicityCorrection()
110{
111 //
112 // Destructor
113 //
114}
115
116//____________________________________________________________________
117Long64_t AliMultiplicityCorrection::Merge(TCollection* list)
118{
119 // Merge a list of AliMultiplicityCorrection objects with this (needed for
120 // PROOF).
121 // Returns the number of merged objects (including this).
122
123 if (!list)
124 return 0;
125
126 if (list->IsEmpty())
127 return 1;
128
129 TIterator* iter = list->MakeIterator();
130 TObject* obj;
131
132 // collections of all histograms
133 TList collections[kESDHists+kMCHists+kCorrHists*2];
134
135 Int_t count = 0;
136 while ((obj = iter->Next())) {
137
138 AliMultiplicityCorrection* entry = dynamic_cast<AliMultiplicityCorrection*> (obj);
139 if (entry == 0)
140 continue;
141
142 for (Int_t i = 0; i < kESDHists; ++i)
143 collections[i].Add(entry->fMultiplicityESD[i]);
144
145 for (Int_t i = 0; i < kMCHists; ++i)
146 collections[kESDHists+i].Add(entry->fMultiplicityMC[i]);
147
148 for (Int_t i = 0; i < kCorrHists; ++i)
149 collections[kESDHists+kMCHists+i].Add(entry->fCorrelation[i]);
150
151 for (Int_t i = 0; i < kCorrHists; ++i)
152 collections[kESDHists+kMCHists+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]);
153
154 count++;
155 }
156
157 for (Int_t i = 0; i < kESDHists; ++i)
158 fMultiplicityESD[i]->Merge(&collections[i]);
159
160 for (Int_t i = 0; i < kMCHists; ++i)
161 fMultiplicityMC[i]->Merge(&collections[kESDHists+i]);
162
163 for (Int_t i = 0; i < kCorrHists; ++i)
164 fCorrelation[i]->Merge(&collections[kESDHists+kMCHists+i]);
165
166 for (Int_t i = 0; i < kCorrHists; ++i)
167 fMultiplicityESDCorrected[i]->Merge(&collections[kESDHists+kMCHists+kCorrHists+i]);
168
169 delete iter;
170
171 return count+1;
172}
173
174//____________________________________________________________________
175Bool_t AliMultiplicityCorrection::LoadHistograms(const Char_t* dir)
176{
177 //
178 // loads the histograms from a file
179 // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
180 //
181
182 if (!dir)
183 dir = GetName();
184
185 if (!gDirectory->cd(dir))
186 return kFALSE;
187
188 // TODO memory leak. old histograms needs to be deleted.
189
190 Bool_t success = kTRUE;
191
192 for (Int_t i = 0; i < kESDHists; ++i)
193 {
194 fMultiplicityESD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityESD[i]->GetName()));
195 if (!fMultiplicityESD[i])
196 success = kFALSE;
197 }
198
199 for (Int_t i = 0; i < kMCHists; ++i)
200 {
201 fMultiplicityMC[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityMC[i]->GetName()));
202 if (!fMultiplicityMC[i])
203 success = kFALSE;
204 }
205
206 for (Int_t i = 0; i < kCorrHists; ++i)
207 {
208 fCorrelation[i] = dynamic_cast<TH3F*> (gDirectory->Get(fCorrelation[i]->GetName()));
209 if (!fCorrelation[i])
210 success = kFALSE;
211 fMultiplicityESDCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fMultiplicityESDCorrected[i]->GetName()));
212 if (!fMultiplicityESDCorrected[i])
213 success = kFALSE;
214 }
215
216 gDirectory->cd("..");
217
218 return success;
219}
220
221//____________________________________________________________________
222void AliMultiplicityCorrection::SaveHistograms()
223{
224 //
225 // saves the histograms
226 //
227
228 gDirectory->mkdir(GetName());
229 gDirectory->cd(GetName());
230
231 for (Int_t i = 0; i < kESDHists; ++i)
232 if (fMultiplicityESD[i])
233 fMultiplicityESD[i]->Write();
234
235 for (Int_t i = 0; i < kMCHists; ++i)
236 if (fMultiplicityMC[i])
237 fMultiplicityMC[i]->Write();
238
239 for (Int_t i = 0; i < kCorrHists; ++i)
240 {
241 if (fCorrelation[i])
242 fCorrelation[i]->Write();
243 if (fMultiplicityESDCorrected[i])
244 fMultiplicityESDCorrected[i]->Write();
245 }
246
247 gDirectory->cd("..");
248}
249
250//____________________________________________________________________
251void AliMultiplicityCorrection::FillGenerated(Float_t vtx, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll)
252{
253 //
254 // Fills an event from MC
255 //
256
257 fMultiplicityMC[0]->Fill(vtx, generated05);
258 fMultiplicityMC[1]->Fill(vtx, generated10);
259 fMultiplicityMC[2]->Fill(vtx, generated15);
260 fMultiplicityMC[3]->Fill(vtx, generated20);
261 fMultiplicityMC[4]->Fill(vtx, generatedAll);
262}
263
264//____________________________________________________________________
265void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20)
266{
267 //
268 // Fills an event from ESD
269 //
270
271 fMultiplicityESD[0]->Fill(vtx, measured05);
272 fMultiplicityESD[1]->Fill(vtx, measured10);
273 fMultiplicityESD[2]->Fill(vtx, measured15);
274 fMultiplicityESD[3]->Fill(vtx, measured20);
275}
276
277//____________________________________________________________________
278void AliMultiplicityCorrection::FillCorrection(Float_t vtx, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20)
279{
280 //
281 // Fills an event into the correlation map with the information from MC and ESD
282 //
283
284 fCorrelation[0]->Fill(vtx, generated05, measured05);
285 fCorrelation[1]->Fill(vtx, generated10, measured10);
286 fCorrelation[2]->Fill(vtx, generated15, measured15);
287 fCorrelation[3]->Fill(vtx, generated20, measured20);
288
289 fCorrelation[4]->Fill(vtx, generatedAll, measured05);
290 fCorrelation[5]->Fill(vtx, generatedAll, measured10);
291 fCorrelation[6]->Fill(vtx, generatedAll, measured15);
292 fCorrelation[7]->Fill(vtx, generatedAll, measured20);
293}
294
295//____________________________________________________________________
296Double_t AliMultiplicityCorrection::MinuitHomogenityPol0(Double_t *params)
297{
298 // homogenity term for minuit fitting
299 // pure function of the parameters
300 // prefers constant function (pol0)
301
302 Double_t chi2 = 0;
303
304 for (Int_t i=1; i<fgMaxParams; ++i)
305 {
306 if (params[i] == 0)
307 continue;
308
309 Double_t right = params[i] / fCurrentMinuitESD->GetBinWidth(i+1);
310 Double_t left = params[i-1] / fCurrentMinuitESD->GetBinWidth(i);
311
312 Double_t diff = (right - left) / right;
313 chi2 += diff * diff;
314 }
315
316 return chi2;
317}
318
319//____________________________________________________________________
320Double_t AliMultiplicityCorrection::MinuitHomogenityPol1(Double_t *params)
321{
322 // homogenity term for minuit fitting
323 // pure function of the parameters
324 // prefers linear function (pol1)
325
326 Double_t chi2 = 0;
327
328 for (Int_t i=2; i<fgMaxParams; ++i)
329 {
330 if (params[i] == 0 || params[i-1] == 0)
331 continue;
332
333 Double_t right = params[i] / fCurrentMinuitESD->GetBinWidth(i+1);
334 Double_t middle = params[i-1] / fCurrentMinuitESD->GetBinWidth(i);
335 Double_t left = params[i-2] / fCurrentMinuitESD->GetBinWidth(i-1);
336
337 Double_t der1 = (right - middle) / fCurrentMinuitESD->GetBinWidth(i);
338 Double_t der2 = (middle - left) / fCurrentMinuitESD->GetBinWidth(i-1);
339
340 Double_t diff = der1 - der2;
341
342 // TODO give additonal weight to big bins
343 chi2 += diff * diff * fCurrentMinuitESD->GetBinWidth(i) * fCurrentMinuitESD->GetBinWidth(i-1);
344
345 //printf("%d --> %f\n", i, diff);
346 }
347
348 return chi2 / 1e5 / 2;
349}
350
351//____________________________________________________________________
352void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t)
353{
354 //
355 // fit function for minuit
356 //
357
358 // TODO take errors into account
359
360 static Int_t callCount = 0;
361
362 Double_t chi2FromFit = 0;
363
364 // loop over multiplicity (x axis of fMultiplicityESD)
365 for (Int_t i=1; i<=fCurrentMinuitESD->GetNbinsX(); ++i)
366 {
367 if (i > fCurrentMinuitCorrelation->GetNbinsY())
368 break;
369
370 Double_t sum = 0;
371 //Double_t error = 0;
372 // loop over generated particles (x axis of fCorrelation) that resulted in reconstruction of i tracks
373 for (Int_t j=1; j<=fCurrentMinuitCorrelation->GetNbinsX(); ++j)
374 {
375 if (j > fgMaxParams)
376 break;
377
378 sum += fCurrentMinuitCorrelation->GetBinContent(j, i) * params[j-1];
379
380 //if (params[j-1] > 0)
381 // error += fCurrentMinuitCorrelation->GetBinError(j, i) * fCurrentMinuitCorrelation->GetBinError(j, i) * params[j-1];
382 //printf("%f ", sum);
383 }
384
385 Double_t diff = fCurrentMinuitESD->GetBinContent(i) - sum;
386 if (fCurrentMinuitESD->GetBinContent(i) > 0)
387 diff /= fCurrentMinuitESD->GetBinContent(i);
388 else
389 diff /= fCurrentMinuitESD->Integral();
390 //error = TMath::Sqrt(error) + fCurrentMinuitESD->GetBinError(i);
391 //if (error <= 0)
392 // error = 1;
393
394 //Double_t tmp = diff / error;
395 //chi2 += tmp * tmp;
396 chi2FromFit += diff * diff;
397
398 //printf("\nExpected sum = %f; Diff for bin %d is %f\n**********************************\n", fCurrentMinuitESD->GetBinContent(i), i, diff);
399 //printf("Diff for bin %d is %f\n", i, diff);
400 }
401
402 Double_t penaltyVal = MinuitHomogenityPol1(params);
403
404 chi2 = chi2FromFit * chi2FromFit + penaltyVal * penaltyVal;
405
406 if ((callCount++ % 100) == 0)
407 printf("%f %f --> %f\n", chi2FromFit, penaltyVal, chi2);
408}
409
410//____________________________________________________________________
411void AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace)
412{
413 //
414 // correct spectrum using minuit chi2 method
415 //
416
417 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
418 Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
419
420 fCurrentMinuitESD = fMultiplicityESD[inputRange]->ProjectionY();
421 fCurrentMinuitESD->Sumw2();
422
423 // empty under/overflow bins in x, otherwise Project3D takes them into account
424 TH3* hist = fCorrelation[correlationID];
425 for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
426 {
427 for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
428 {
429 hist->SetBinContent(0, y, z, 0);
430 hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
431 }
432 }
433
434 fCurrentMinuitCorrelation = hist->Project3D("zy");
435 fCurrentMinuitCorrelation->Sumw2();
436
437 // normalize correction for given nPart
438 for (Int_t i=1; i<=fCurrentMinuitCorrelation->GetNbinsX(); ++i)
439 {
440 Double_t sum = fCurrentMinuitCorrelation->Integral(i, i, 1, fCurrentMinuitCorrelation->GetNbinsY());
441 if (sum <= 0)
442 continue;
443 for (Int_t j=1; j<=fCurrentMinuitCorrelation->GetNbinsY(); ++j)
444 {
445 // npart sum to 1
446 fCurrentMinuitCorrelation->SetBinContent(i, j, fCurrentMinuitCorrelation->GetBinContent(i, j) / sum);
447 fCurrentMinuitCorrelation->SetBinError(i, j, fCurrentMinuitCorrelation->GetBinError(i, j) / sum);
448 }
449 }
450
451 // small hack to get around charge conservation for full phase space ;-)
452 /*if (fullPhaseSpace)
453 {
454 for (Int_t i=2; i<=50; i+=2)
455 for (Int_t j=1; j<=fCurrentMinuitCorrelation->GetNbinsY(); ++j)
456 fCurrentMinuitCorrelation->SetBinContent(i, j, fCurrentMinuitCorrelation->GetBinContent(i-1, j));
457 }*/
458
459 TCanvas* canvas1 = new TCanvas("ApplyMinuitFit", "ApplyMinuitFit", 800, 400);
460 canvas1->Divide(2, 1);
461 canvas1->cd(1);
462 fCurrentMinuitESD->DrawCopy();
463 canvas1->cd(2);
464 fCurrentMinuitCorrelation->DrawCopy("COLZ");
465
466 // Initialize TMinuit via generic fitter interface
467 TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgMaxParams);
468
469 minuit->SetFCN(MinuitFitFunction);
470
471 Double_t results[fgMaxParams];
472
473 //TH1* proj = fMultiplicityMC[mcTarget]->ProjectionY();
474 for (Int_t i=0; i<fgMaxParams; ++i)
475 {
476 //results[i] = 1;
477 results[i] = fCurrentMinuitESD->GetBinContent(i+1);
478 //results[i] = proj->GetBinContent(i+1);
479 minuit->SetParameter(i, Form("param%d", i), results[i], ((results[i] > 1) ? (results[i] / 10) : 0), 0, fCurrentMinuitESD->GetMaximum() * 100);
480 }
481
482 Int_t dummy;
483 Double_t chi2;
484 MinuitFitFunction(dummy, 0, chi2, results, 0);
485 printf("Chi2 of right solution is = %f\n", chi2);
486
487 //return;
488
489 Double_t arglist[100];
490 arglist[0] = 100000;
491 //minuit->ExecuteCommand("SET PRINT", arglist, 1);
492 minuit->ExecuteCommand("MIGRAD", arglist, 1);
493 //minuit->ExecuteCommand("MINOS", arglist, 0);
494
495 for (Int_t i=0; i<fgMaxParams; ++i)
496 {
497 results[i] = minuit->GetParameter(i);
498 fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, results[i]);
499 fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0);
500 }
501
502 printf("Penalty is %f\n", MinuitHomogenityPol1(results));
503
504 fMultiplicityESDCorrected[correlationID]->GetXaxis()->SetRangeUser(0, fgMaxParams);
505
506 DrawComparison(mcTarget, correlationID);
507}
508
509//____________________________________________________________________
510void AliMultiplicityCorrection::NormalizeToBinWidth(TH1* hist)
511{
512 //
513 // normalizes a 1-d histogram to its bin width
514 //
515
516 for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
517 {
518 hist->SetBinContent(i, hist->GetBinContent(i) / hist->GetBinWidth(i));
519 hist->SetBinError(i, hist->GetBinError(i) / hist->GetBinWidth(i));
520 }
521}
522
523//____________________________________________________________________
524void AliMultiplicityCorrection::NormalizeToBinWidth(TH2* hist)
525{
526 //
527 // normalizes a 2-d histogram to its bin width (x width * y width)
528 //
529
530 for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
531 for (Int_t j=1; j<=hist->GetNbinsY(); ++j)
532 {
533 Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j);
534 hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor);
535 hist->SetBinError(i, j, hist->GetBinError(i, j) / factor);
536 }
537}
538
539//____________________________________________________________________
540void AliMultiplicityCorrection::ApplyMinuitFitAll()
541{
542 //
543 // fit all eta ranges
544 //
545
546 for (Int_t i=0; i<kESDHists; ++i)
547 {
548 ApplyMinuitFit(i, kFALSE);
549 ApplyMinuitFit(i, kTRUE);
550 }
551}
552
553//____________________________________________________________________
554void AliMultiplicityCorrection::DrawHistograms()
555{
556 printf("ESD:\n");
557
558 TCanvas* canvas1 = new TCanvas("fMultiplicityESD", "fMultiplicityESD", 900, 600);
559 canvas1->Divide(3, 2);
560 for (Int_t i = 0; i < kESDHists; ++i)
561 {
562 canvas1->cd(i+1);
563 fMultiplicityESD[i]->DrawCopy("COLZ");
564 printf("%d --> %f\n", i, (Float_t) fMultiplicityESD[i]->ProjectionY()->GetMean());
565 }
566
567 printf("MC:\n");
568
569 TCanvas* canvas2 = new TCanvas("fMultiplicityMC", "fMultiplicityMC", 900, 600);
570 canvas2->Divide(3, 2);
571 for (Int_t i = 0; i < kMCHists; ++i)
572 {
573 canvas2->cd(i+1);
574 fMultiplicityMC[i]->DrawCopy("COLZ");
575 printf("%d --> %f\n", i, fMultiplicityMC[i]->ProjectionY()->GetMean());
576 }
577
578 TCanvas* canvas3 = new TCanvas("fCorrelation", "fCorrelation", 900, 900);
579 canvas3->Divide(3, 3);
580 for (Int_t i = 0; i < kCorrHists; ++i)
581 {
582 canvas3->cd(i+1);
583 TH1* proj = fCorrelation[i]->Project3D("zy");
584 proj->DrawCopy("COLZ");
585 }
586}
587
588//____________________________________________________________________
589void AliMultiplicityCorrection::DrawComparison(Int_t mcID, Int_t esdCorrId)
590{
591 TString name;
592 name.Form("DrawComparison_%d_%d", mcID, esdCorrId);
593
594 TCanvas* canvas1 = new TCanvas(name, name, 900, 300);
595 canvas1->Divide(3, 1);
596
597 canvas1->cd(1);
598 TH1* proj = fMultiplicityMC[mcID]->ProjectionY();
599 NormalizeToBinWidth(proj);
600 NormalizeToBinWidth(fMultiplicityESDCorrected[esdCorrId]);
601 proj->GetXaxis()->SetRangeUser(0, 200);
602 proj->DrawCopy("HIST");
603 gPad->SetLogy();
604
605 fMultiplicityESDCorrected[esdCorrId]->SetMarkerStyle(3);
606 fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAMEP");
607
608 canvas1->cd(2);
609 fMultiplicityESDCorrected[esdCorrId]->DrawCopy("");
610
611 canvas1->cd(3);
612 TH1* clone = dynamic_cast<TH1*> (proj->Clone("clone"));
613 clone->Divide(fMultiplicityESDCorrected[esdCorrId]);
614 clone->SetTitle("Ratio;Npart;MC/ESD");
615 clone->GetYaxis()->SetRangeUser(0.8, 1.2);
616 clone->GetXaxis()->SetRangeUser(0, 200);
617 clone->DrawCopy("HIST");
618
619 /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
620 legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
621 legend->AddEntry(fMultiplicityMC, "MC");
622 legend->AddEntry(fMultiplicityESD, "ESD");
623 legend->Draw();*/
624}
625
626//____________________________________________________________________
627void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace)
628{
629 //
630 // correct spectrum using bayesian method
631 //
632
633 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
634 Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
635
636 fCurrentMinuitESD = fMultiplicityESD[inputRange]->ProjectionY();
637 fCurrentMinuitESD->Sumw2();
638
639 // empty under/overflow bins in x, otherwise Project3D takes them into account
640 TH3* hist = fCorrelation[correlationID];
641 for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
642 {
643 for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
644 {
645 hist->SetBinContent(0, y, z, 0);
646 hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
647 }
648 }
649
650 fCurrentMinuitCorrelation = hist->Project3D("zy");
651 fCurrentMinuitCorrelation->Sumw2();
652
653 // normalize correction for given nPart
654 for (Int_t i=1; i<=fCurrentMinuitCorrelation->GetNbinsX(); ++i)
655 {
656 Double_t sum = fCurrentMinuitCorrelation->Integral(i, i, 1, fCurrentMinuitCorrelation->GetNbinsY());
657 if (sum <= 0)
658 continue;
659 for (Int_t j=1; j<=fCurrentMinuitCorrelation->GetNbinsY(); ++j)
660 {
661 // npart sum to 1
662 fCurrentMinuitCorrelation->SetBinContent(i, j, fCurrentMinuitCorrelation->GetBinContent(i, j) / sum);
663 fCurrentMinuitCorrelation->SetBinError(i, j, fCurrentMinuitCorrelation->GetBinError(i, j) / sum);
664 }
665 }
666
6127aca6 667 Float_t regPar = 0.01;
668
669 Float_t norm = 0;
670 // pick prior distribution
671 TH1F* hPrior = (TH1F*)fCurrentMinuitESD->Clone("prior");
672 for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
673 norm = norm + hPrior->GetBinContent(t);
674 for (Int_t t=1; t<=hPrior->GetNbinsX(); t++) {
675 hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm);
676
677 printf(" bin %d content %f \n", t, hPrior->GetBinContent(t));
678
679 }
680
681 // define hist to store guess of true
682 TH1F* hTrue = (TH1F*)fCurrentMinuitESD->Clone("prior");
683 // hTrue->Reset();
684
685 // clone response R
686 TH2F* hResponse = (TH2F*)fCurrentMinuitCorrelation->Clone("pij");
687
688 // histogram to store the inverse response calculated from bayes theorem (from prior and response)
689 // IR
690 TAxis* xAxis = hResponse->GetYaxis();
691 TAxis* yAxis = hResponse->GetXaxis();
692
693 TH2F* hInverseResponseBayes = (TH2F*)hResponse->Clone("pji");
694 //new TH2F("pji","pji",
695 // xAxis->GetNbins(),xAxis->GetXbins()->GetArray(),
696 // yAxis->GetNbins(),yAxis->GetXbins()->GetArray());
697 hInverseResponseBayes->Reset();
698
699 // unfold...
700 Int_t iterations = 20;
701 for (Int_t i=0; i<iterations; i++) {
702 printf(" iteration %i \n", i);
703
704 // calculate IR from Beyes theorem
705 // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
706 for(Int_t t = 1; t<=hResponse->GetNbinsX(); t++) {
707 for (Int_t m=1; m<=hResponse->GetNbinsY(); m++) {
708
709 norm = 0;
710 for(Int_t t2 = 1; t2<=hResponse->GetNbinsX(); t2++)
711 norm = norm + (hResponse->GetBinContent(t2,m) * hPrior->GetBinContent(t2));
712
713 if (norm==0)
714 hInverseResponseBayes->SetBinContent(t,m,0);
715 else
716 hInverseResponseBayes->SetBinContent(t,m, hResponse->GetBinContent(t,m) * hPrior->GetBinContent(t)/norm);
717 }
718 }
719 // calculate true assuming IR is the correct inverse response
720 for(Int_t t = 1; t<=hResponse->GetNbinsX(); t++) {
721 hTrue ->SetBinContent(t, 0);
722 for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
723 hTrue->SetBinContent(t, hTrue->GetBinContent(t) + fCurrentMinuitESD->GetBinContent(m)*hInverseResponseBayes->GetBinContent(t,m));
724 }
725
726 // regularization (simple smoothing) NB : not good bin width!!!
727 TH1F* hTrueSmoothed = (TH1F*)hTrue->Clone("truesmoothed");
728
729 for (Int_t t=2; t<hTrue->GetNbinsX()-1; t++) {
730 Float_t average = (hTrue->GetBinContent(t-1) + hTrue->GetBinContent(t) + hTrue->GetBinContent(t+1))/3.;
731 // weight the average with the regularization parameter
732 hTrueSmoothed->SetBinContent(t, (1-regPar)*hTrue->GetBinContent(t) + (regPar*average));
733 }
734
735 // use smoothed true (normalized) as new prior
736 norm = 0;
737 for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
738 norm = norm + hTrueSmoothed->GetBinContent(t);
739 for (Int_t t=1; t<=hPrior->GetNbinsX(); t++) {
740 hPrior->SetBinContent(t, hTrueSmoothed->GetBinContent(t)/norm);
741 hTrue->SetBinContent(t, hTrueSmoothed->GetBinContent(t));
742 }
743 } // end of iterations
744
745 for (Int_t t=1; t<=fMultiplicityESDCorrected[inputRange]->GetNbinsX(); t++) {
746 fMultiplicityESDCorrected[inputRange]->SetBinContent(t, hTrue->GetBinContent(t));
747 fMultiplicityESDCorrected[inputRange]->SetBinError(t, 0.05*hTrue->GetBinContent(t));
748
749 printf(" bin %d content %f \n", t, hPrior->GetBinContent(t));
750 }
751
0a173978 752}