coding conventions
[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.
6d81c2de 21// e.g. chi2 minimization and bayesian unfolding
0a173978 22//
23// Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch
24
25#include "AliMultiplicityCorrection.h"
26
27#include <TFile.h>
28#include <TH1F.h>
29#include <TH2F.h>
30#include <TH3F.h>
31#include <TDirectory.h>
32#include <TVirtualFitter.h>
33#include <TCanvas.h>
34#include <TString.h>
9ca6be09 35#include <TF1.h>
d560b581 36#include <TMath.h>
3602328d 37#include <TCollection.h>
447c325d 38#include <TLegend.h>
39#include <TLine.h>
0b4bfd98 40#include <TRandom.h>
41#include <TProfile.h>
42#include <TProfile2D.h>
0a173978 43
44ClassImp(AliMultiplicityCorrection)
45
6d81c2de 46const Int_t AliMultiplicityCorrection::fgkMaxInput = 250; // bins in measured histogram
47const Int_t AliMultiplicityCorrection::fgkMaxParams = 250; // bins in unfolded histogram = number of fit params
447c325d 48
6d81c2de 49TMatrixD* AliMultiplicityCorrection::fgCorrelationMatrix = 0;
50TMatrixD* AliMultiplicityCorrection::fgCorrelationCovarianceMatrix = 0;
51TVectorD* AliMultiplicityCorrection::fgCurrentESDVector = 0;
52TVectorD* AliMultiplicityCorrection::fgEntropyAPriori = 0;
53
54AliMultiplicityCorrection::RegularizationType AliMultiplicityCorrection::fgRegularizationType = AliMultiplicityCorrection::kPol1;
55Float_t AliMultiplicityCorrection::fgRegularizationWeight = 10000;
56
57TF1* AliMultiplicityCorrection::fgNBD = 0;
58
59Float_t AliMultiplicityCorrection::fgBayesianSmoothing = 1; // smoothing parameter (0 = no smoothing)
60Int_t AliMultiplicityCorrection::fgBayesianIterations = 100; // number of iterations in Bayesian method
0a173978 61
0b4bfd98 62// These are the areas where the quality of the unfolding results are evaluated
63// Default defined here, call SetQualityRegions to change them
64// unit is in multiplicity (not in bin!)
65
66// SPD: peak area - flat area - low stat area
67Int_t AliMultiplicityCorrection::fgQualityRegionsB[kQualityRegions] = {4, 60, 180};
68Int_t AliMultiplicityCorrection::fgQualityRegionsE[kQualityRegions] = {20, 140, 210};
69
70//____________________________________________________________________
71void AliMultiplicityCorrection::SetQualityRegions(Bool_t SPDStudy)
72{
73 //
74 // sets the quality region definition to TPC or SPD
75 //
76
77 if (SPDStudy)
78 {
79 // SPD: peak area - flat area - low stat area
80 fgQualityRegionsB[0] = 4;
81 fgQualityRegionsE[0] = 20;
82
83 fgQualityRegionsB[1] = 60;
84 fgQualityRegionsE[1] = 140;
85
86 fgQualityRegionsB[2] = 180;
87 fgQualityRegionsE[2] = 210;
6d81c2de 88
89 Printf("AliMultiplicityCorrection::SetQualityRegions --> Enabled quality regions for SPD");
0b4bfd98 90 }
91 else
92 {
93 // TPC: peak area - flat area - low stat area
94 fgQualityRegionsB[0] = 4;
95 fgQualityRegionsE[0] = 12;
96
97 fgQualityRegionsB[1] = 25;
98 fgQualityRegionsE[1] = 55;
99
100 fgQualityRegionsB[2] = 88;
101 fgQualityRegionsE[2] = 108;
6d81c2de 102
103 Printf("AliMultiplicityCorrection::SetQualityRegions --> Enabled quality regions for TPC");
0b4bfd98 104 }
105}
106
0a173978 107//____________________________________________________________________
108AliMultiplicityCorrection::AliMultiplicityCorrection() :
6d81c2de 109 TNamed(), fCurrentESD(0), fCurrentCorrelation(0), fCurrentEfficiency(0), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
0a173978 110{
111 //
112 // default constructor
113 //
114
115 for (Int_t i = 0; i < kESDHists; ++i)
116 fMultiplicityESD[i] = 0;
117
118 for (Int_t i = 0; i < kMCHists; ++i)
cfc19dd5 119 {
120 fMultiplicityVtx[i] = 0;
121 fMultiplicityMB[i] = 0;
122 fMultiplicityINEL[i] = 0;
123 }
0a173978 124
125 for (Int_t i = 0; i < kCorrHists; ++i)
126 {
127 fCorrelation[i] = 0;
128 fMultiplicityESDCorrected[i] = 0;
129 }
6d81c2de 130
131 for (Int_t i = 0; i < kQualityRegions; ++i)
132 fQuality[i] = 0;
0a173978 133}
134
135//____________________________________________________________________
136AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const Char_t* title) :
0b4bfd98 137 TNamed(name, title), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0)
0a173978 138{
139 //
140 // named constructor
141 //
142
143 // do not add this hists to the directory
144 Bool_t oldStatus = TH1::AddDirectoryStatus();
145 TH1::AddDirectory(kFALSE);
146
b477f4f2 147 /*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 148 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,
149 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
150 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
151 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
152 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5,
153 50.5, 55.5, 60.5, 65.5, 70.5, 75.5, 80.5, 85.5, 90.5, 95.5,
154 100.5, 105.5, 110.5, 115.5, 120.5, 125.5, 130.5, 135.5, 140.5, 145.5,
155 150.5, 160.5, 170.5, 180.5, 190.5, 200.5, 210.5, 220.5, 230.5, 240.5,
156 250.5, 275.5, 300.5, 325.5, 350.5, 375.5, 400.5, 425.5, 450.5, 475.5,
9ca6be09 157 500.5 };
158 //525.5, 550.5, 575.5, 600.5, 625.5, 650.5, 675.5, 700.5, 725.5,
159 //750.5, 775.5, 800.5, 825.5, 850.5, 875.5, 900.5, 925.5, 950.5, 975.5,
160 //1000.5 };
0a173978 161
b477f4f2 162 #define VTXBINNING 10, binLimitsVtx
6d81c2de 163 #define NBINNING fgkMaxParams, binLimitsN*/
b477f4f2 164
447c325d 165 #define NBINNING 501, -0.5, 500.5
166 #define VTXBINNING 1, -10, 10
0a173978 167
168 for (Int_t i = 0; i < kESDHists; ++i)
b477f4f2 169 fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", VTXBINNING, NBINNING);
0a173978 170
171 for (Int_t i = 0; i < kMCHists; ++i)
172 {
cfc19dd5 173 fMultiplicityVtx[i] = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("fMultiplicityVtx%d", i)));
174 fMultiplicityVtx[i]->SetTitle("fMultiplicityVtx;vtx-z;Npart");
175
176 fMultiplicityMB[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityMB%d", i)));
177 fMultiplicityMB[i]->SetTitle("fMultiplicityMB");
178
179 fMultiplicityINEL[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityINEL%d", i)));
180 fMultiplicityINEL[i]->SetTitle("fMultiplicityINEL");
0a173978 181 }
182
183 for (Int_t i = 0; i < kCorrHists; ++i)
184 {
b477f4f2 185 fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;Npart;Ntracks", VTXBINNING, NBINNING, NBINNING);
186 fMultiplicityESDCorrected[i] = new TH1F(Form("fMultiplicityESDCorrected%d", i), "fMultiplicityESDCorrected;Npart;dN/dN", NBINNING);
0a173978 187 }
188
189 TH1::AddDirectory(oldStatus);
190}
191
192//____________________________________________________________________
193AliMultiplicityCorrection::~AliMultiplicityCorrection()
194{
195 //
196 // Destructor
197 //
6d81c2de 198
199 for (Int_t i = 0; i < kESDHists; ++i)
200 {
201 if (fMultiplicityESD[i])
202 delete fMultiplicityESD[i];
203 fMultiplicityESD[i] = 0;
204 }
205
206 for (Int_t i = 0; i < kMCHists; ++i)
207 {
208 if (fMultiplicityVtx[i])
209 delete fMultiplicityVtx[i];
210 fMultiplicityVtx[i] = 0;
211
212 if (fMultiplicityMB[i])
213 delete fMultiplicityMB[i];
214 fMultiplicityMB[i] = 0;
215
216 if (fMultiplicityINEL[i])
217 delete fMultiplicityINEL[i];
218 fMultiplicityINEL[i] = 0;
219 }
220
221 for (Int_t i = 0; i < kCorrHists; ++i)
222 {
223 if (fCorrelation[i])
224 delete fCorrelation[i];
225 fCorrelation[i] = 0;
226
227 if (fMultiplicityESDCorrected[i])
228 delete fMultiplicityESDCorrected[i];
229 fMultiplicityESDCorrected[i] = 0;
230 }
0a173978 231}
232
233//____________________________________________________________________
234Long64_t AliMultiplicityCorrection::Merge(TCollection* list)
235{
236 // Merge a list of AliMultiplicityCorrection objects with this (needed for
237 // PROOF).
238 // Returns the number of merged objects (including this).
239
240 if (!list)
241 return 0;
242
243 if (list->IsEmpty())
244 return 1;
245
246 TIterator* iter = list->MakeIterator();
247 TObject* obj;
248
249 // collections of all histograms
cfc19dd5 250 TList collections[kESDHists+kMCHists*3+kCorrHists*2];
0a173978 251
252 Int_t count = 0;
253 while ((obj = iter->Next())) {
254
255 AliMultiplicityCorrection* entry = dynamic_cast<AliMultiplicityCorrection*> (obj);
256 if (entry == 0)
257 continue;
258
259 for (Int_t i = 0; i < kESDHists; ++i)
260 collections[i].Add(entry->fMultiplicityESD[i]);
261
262 for (Int_t i = 0; i < kMCHists; ++i)
cfc19dd5 263 {
264 collections[kESDHists+i].Add(entry->fMultiplicityVtx[i]);
265 collections[kESDHists+kMCHists+i].Add(entry->fMultiplicityMB[i]);
266 collections[kESDHists+kMCHists*2+i].Add(entry->fMultiplicityINEL[i]);
267 }
0a173978 268
269 for (Int_t i = 0; i < kCorrHists; ++i)
cfc19dd5 270 collections[kESDHists+kMCHists*3+i].Add(entry->fCorrelation[i]);
0a173978 271
272 for (Int_t i = 0; i < kCorrHists; ++i)
cfc19dd5 273 collections[kESDHists+kMCHists*3+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]);
0a173978 274
275 count++;
276 }
277
278 for (Int_t i = 0; i < kESDHists; ++i)
279 fMultiplicityESD[i]->Merge(&collections[i]);
280
281 for (Int_t i = 0; i < kMCHists; ++i)
cfc19dd5 282 {
283 fMultiplicityVtx[i]->Merge(&collections[kESDHists+i]);
284 fMultiplicityMB[i]->Merge(&collections[kESDHists+kMCHists+i]);
285 fMultiplicityINEL[i]->Merge(&collections[kESDHists+kMCHists*2+i]);
286 }
0a173978 287
288 for (Int_t i = 0; i < kCorrHists; ++i)
cfc19dd5 289 fCorrelation[i]->Merge(&collections[kESDHists+kMCHists*3+i]);
0a173978 290
291 for (Int_t i = 0; i < kCorrHists; ++i)
cfc19dd5 292 fMultiplicityESDCorrected[i]->Merge(&collections[kESDHists+kMCHists*3+kCorrHists+i]);
0a173978 293
294 delete iter;
295
296 return count+1;
297}
298
299//____________________________________________________________________
300Bool_t AliMultiplicityCorrection::LoadHistograms(const Char_t* dir)
301{
302 //
303 // loads the histograms from a file
304 // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
305 //
306
307 if (!dir)
308 dir = GetName();
309
310 if (!gDirectory->cd(dir))
311 return kFALSE;
312
313 // TODO memory leak. old histograms needs to be deleted.
314
315 Bool_t success = kTRUE;
316
317 for (Int_t i = 0; i < kESDHists; ++i)
318 {
319 fMultiplicityESD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityESD[i]->GetName()));
320 if (!fMultiplicityESD[i])
321 success = kFALSE;
322 }
323
324 for (Int_t i = 0; i < kMCHists; ++i)
325 {
cfc19dd5 326 fMultiplicityVtx[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityVtx[i]->GetName()));
327 fMultiplicityMB[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityMB[i]->GetName()));
328 fMultiplicityINEL[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityINEL[i]->GetName()));
329 if (!fMultiplicityVtx[i] || !fMultiplicityMB[i] || !fMultiplicityINEL[i])
0a173978 330 success = kFALSE;
331 }
332
333 for (Int_t i = 0; i < kCorrHists; ++i)
334 {
335 fCorrelation[i] = dynamic_cast<TH3F*> (gDirectory->Get(fCorrelation[i]->GetName()));
336 if (!fCorrelation[i])
337 success = kFALSE;
338 fMultiplicityESDCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fMultiplicityESDCorrected[i]->GetName()));
339 if (!fMultiplicityESDCorrected[i])
340 success = kFALSE;
341 }
342
343 gDirectory->cd("..");
344
345 return success;
346}
347
348//____________________________________________________________________
349void AliMultiplicityCorrection::SaveHistograms()
350{
351 //
352 // saves the histograms
353 //
354
355 gDirectory->mkdir(GetName());
356 gDirectory->cd(GetName());
357
358 for (Int_t i = 0; i < kESDHists; ++i)
359 if (fMultiplicityESD[i])
360 fMultiplicityESD[i]->Write();
361
362 for (Int_t i = 0; i < kMCHists; ++i)
cfc19dd5 363 {
364 if (fMultiplicityVtx[i])
365 fMultiplicityVtx[i]->Write();
366 if (fMultiplicityMB[i])
367 fMultiplicityMB[i]->Write();
368 if (fMultiplicityINEL[i])
369 fMultiplicityINEL[i]->Write();
370 }
0a173978 371
372 for (Int_t i = 0; i < kCorrHists; ++i)
373 {
374 if (fCorrelation[i])
375 fCorrelation[i]->Write();
376 if (fMultiplicityESDCorrected[i])
377 fMultiplicityESDCorrected[i]->Write();
378 }
379
380 gDirectory->cd("..");
381}
382
383//____________________________________________________________________
cfc19dd5 384void AliMultiplicityCorrection::FillGenerated(Float_t vtx, Bool_t triggered, Bool_t vertex, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll)
0a173978 385{
386 //
387 // Fills an event from MC
388 //
389
cfc19dd5 390 if (triggered)
391 {
392 fMultiplicityMB[0]->Fill(vtx, generated05);
393 fMultiplicityMB[1]->Fill(vtx, generated10);
394 fMultiplicityMB[2]->Fill(vtx, generated15);
395 fMultiplicityMB[3]->Fill(vtx, generated20);
396 fMultiplicityMB[4]->Fill(vtx, generatedAll);
397
398 if (vertex)
399 {
400 fMultiplicityVtx[0]->Fill(vtx, generated05);
401 fMultiplicityVtx[1]->Fill(vtx, generated10);
402 fMultiplicityVtx[2]->Fill(vtx, generated15);
403 fMultiplicityVtx[3]->Fill(vtx, generated20);
404 fMultiplicityVtx[4]->Fill(vtx, generatedAll);
405 }
406 }
407
408 fMultiplicityINEL[0]->Fill(vtx, generated05);
409 fMultiplicityINEL[1]->Fill(vtx, generated10);
410 fMultiplicityINEL[2]->Fill(vtx, generated15);
411 fMultiplicityINEL[3]->Fill(vtx, generated20);
412 fMultiplicityINEL[4]->Fill(vtx, generatedAll);
0a173978 413}
414
415//____________________________________________________________________
416void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20)
417{
418 //
419 // Fills an event from ESD
420 //
421
422 fMultiplicityESD[0]->Fill(vtx, measured05);
423 fMultiplicityESD[1]->Fill(vtx, measured10);
424 fMultiplicityESD[2]->Fill(vtx, measured15);
425 fMultiplicityESD[3]->Fill(vtx, measured20);
426}
427
428//____________________________________________________________________
429void 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)
430{
431 //
432 // Fills an event into the correlation map with the information from MC and ESD
433 //
434
435 fCorrelation[0]->Fill(vtx, generated05, measured05);
436 fCorrelation[1]->Fill(vtx, generated10, measured10);
437 fCorrelation[2]->Fill(vtx, generated15, measured15);
438 fCorrelation[3]->Fill(vtx, generated20, measured20);
439
440 fCorrelation[4]->Fill(vtx, generatedAll, measured05);
441 fCorrelation[5]->Fill(vtx, generatedAll, measured10);
442 fCorrelation[6]->Fill(vtx, generatedAll, measured15);
443 fCorrelation[7]->Fill(vtx, generatedAll, measured20);
444}
445
446//____________________________________________________________________
447c325d 447Double_t AliMultiplicityCorrection::RegularizationPol0(TVectorD& params)
0a173978 448{
449 // homogenity term for minuit fitting
450 // pure function of the parameters
451 // prefers constant function (pol0)
452
453 Double_t chi2 = 0;
454
447c325d 455 // ignore the first bin here. on purpose...
6d81c2de 456 for (Int_t i=2; i<fgkMaxParams; ++i)
0a173978 457 {
dd701109 458 Double_t right = params[i];
459 Double_t left = params[i-1];
0a173978 460
447c325d 461 if (right != 0)
462 {
463 Double_t diff = 1 - left / right;
464 chi2 += diff * diff;
465 }
0a173978 466 }
467
447c325d 468 return chi2 / 100.0;
0a173978 469}
470
471//____________________________________________________________________
447c325d 472Double_t AliMultiplicityCorrection::RegularizationPol1(TVectorD& params)
0a173978 473{
474 // homogenity term for minuit fitting
475 // pure function of the parameters
476 // prefers linear function (pol1)
477
478 Double_t chi2 = 0;
479
447c325d 480 // ignore the first bin here. on purpose...
6d81c2de 481 for (Int_t i=2+1; i<fgkMaxParams; ++i)
0a173978 482 {
dd701109 483 if (params[i-1] == 0)
0a173978 484 continue;
485
dd701109 486 Double_t right = params[i];
487 Double_t middle = params[i-1];
488 Double_t left = params[i-2];
489
490 Double_t der1 = (right - middle);
491 Double_t der2 = (middle - left);
492
493 Double_t diff = (der1 - der2) / middle;
494
495 chi2 += diff * diff;
496 }
0a173978 497
dd701109 498 return chi2;
499}
0a173978 500
dd701109 501//____________________________________________________________________
0b4bfd98 502Double_t AliMultiplicityCorrection::RegularizationLog(TVectorD& params)
dd701109 503{
504 // homogenity term for minuit fitting
505 // pure function of the parameters
506 // prefers linear function (pol1)
0a173978 507
dd701109 508 Double_t chi2 = 0;
0a173978 509
6d81c2de 510 /*Float_t der[fgkMaxParams];
dd701109 511
6d81c2de 512 for (Int_t i=0; i<fgkMaxParams-1; ++i)
dd701109 513 der[i] = params[i+1] - params[i];
514
6d81c2de 515 for (Int_t i=0; i<fgkMaxParams-6; ++i)
dd701109 516 {
517 for (Int_t j=1; j<=5; ++j)
518 {
519 Double_t diff = der[i+j] - der[i];
520 chi2 += diff * diff;
521 }
0b4bfd98 522 }*/
523
524 // ignore the first bin here. on purpose...
6d81c2de 525 for (Int_t i=2+1; i<fgkMaxParams; ++i)
0b4bfd98 526 {
527 if (params[i-1] == 0)
528 continue;
529
530 Double_t right = log(params[i]);
531 Double_t middle = log(params[i-1]);
532 Double_t left = log(params[i-2]);
533
534 Double_t der1 = (right - middle);
535 Double_t der2 = (middle - left);
536
537 Double_t diff = (der1 - der2) / middle;
538
539 chi2 += diff * diff;
0a173978 540 }
541
0b4bfd98 542 return chi2 * 100;
0a173978 543}
544
9ca6be09 545//____________________________________________________________________
447c325d 546Double_t AliMultiplicityCorrection::RegularizationTotalCurvature(TVectorD& params)
9ca6be09 547{
548 // homogenity term for minuit fitting
549 // pure function of the parameters
550 // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments,
551 // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp.
552
553 Double_t chi2 = 0;
554
447c325d 555 // ignore the first bin here. on purpose...
6d81c2de 556 for (Int_t i=3; i<fgkMaxParams; ++i)
9ca6be09 557 {
dd701109 558 Double_t right = params[i];
559 Double_t middle = params[i-1];
560 Double_t left = params[i-2];
9ca6be09 561
dd701109 562 Double_t der1 = (right - middle);
563 Double_t der2 = (middle - left);
9ca6be09 564
447c325d 565 Double_t diff = (der1 - der2);
9ca6be09 566
447c325d 567 chi2 += diff * diff;
9ca6be09 568 }
569
447c325d 570 return chi2 * 1e4;
3602328d 571}
572
573//____________________________________________________________________
447c325d 574Double_t AliMultiplicityCorrection::RegularizationEntropy(TVectorD& params)
3602328d 575{
576 // homogenity term for minuit fitting
577 // pure function of the parameters
578 // calculates entropy, from
579 // The method of reduced cross-entropy (M. Schmelling 1993)
580
3602328d 581 Double_t paramSum = 0;
447c325d 582 // ignore the first bin here. on purpose...
6d81c2de 583 for (Int_t i=1; i<fgkMaxParams; ++i)
dd701109 584 paramSum += params[i];
3602328d 585
586 Double_t chi2 = 0;
6d81c2de 587 for (Int_t i=1; i<fgkMaxParams; ++i)
3602328d 588 {
0b4bfd98 589 //Double_t tmp = params[i] / paramSum;
590 Double_t tmp = params[i];
6d81c2de 591 if (tmp > 0 && (*fgEntropyAPriori)[i] > 0)
592 chi2 += tmp * TMath::Log(tmp / (*fgEntropyAPriori)[i]);
3602328d 593 }
3602328d 594
447c325d 595 return 10.0 + chi2;
dd701109 596}
597
598//____________________________________________________________________
599void AliMultiplicityCorrection::MinuitNBD(Int_t& unused1, Double_t* unused2, Double_t& chi2, Double_t *params, Int_t unused3)
600{
601 //
602 // fit function for minuit
603 // does: nbd
604 // func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 50);
605 // func->SetParNames("scaling", "averagen", "k");
447c325d 606 // func->SetParLimits(0, 0.001, fCurrentESD->GetMaximum() * 1000);
607 // func->SetParLimits(1, 0.001, 1000);
608 // func->SetParLimits(2, 0.001, 1000);
dd701109 609 //
610
6d81c2de 611 fgNBD->SetParameters(params[0], params[1], params[2]);
dd701109 612
6d81c2de 613 Double_t params2[fgkMaxParams];
dd701109 614
6d81c2de 615 for (Int_t i=0; i<fgkMaxParams; ++i)
616 params2[i] = fgNBD->Eval(i);
dd701109 617
618 MinuitFitFunction(unused1, unused2, chi2, params2, unused3);
619
620 printf("%f %f %f --> %f\n", params[0], params[1], params[2], chi2);
9ca6be09 621}
622
0a173978 623//____________________________________________________________________
624void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t)
625{
626 //
627 // fit function for minuit
cfc19dd5 628 // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix
0a173978 629 //
630
cfc19dd5 631 // d
6d81c2de 632 static TVectorD paramsVector(fgkMaxParams);
633 for (Int_t i=0; i<fgkMaxParams; ++i)
447c325d 634 paramsVector[i] = params[i] * params[i];
635
636 // calculate penalty factor
637 Double_t penaltyVal = 0;
6d81c2de 638 switch (fgRegularizationType)
447c325d 639 {
640 case kNone: break;
641 case kPol0: penaltyVal = RegularizationPol0(paramsVector); break;
642 case kPol1: penaltyVal = RegularizationPol1(paramsVector); break;
643 case kCurvature: penaltyVal = RegularizationTotalCurvature(paramsVector); break;
644 case kEntropy: penaltyVal = RegularizationEntropy(paramsVector); break;
0b4bfd98 645 case kLog: penaltyVal = RegularizationLog(paramsVector); break;
447c325d 646 }
647
648 //if (penaltyVal > 1e10)
649 // paramsVector2.Print();
650
6d81c2de 651 penaltyVal *= fgRegularizationWeight;
cfc19dd5 652
653 // Ad
6d81c2de 654 TVectorD measGuessVector(fgkMaxInput);
655 measGuessVector = (*fgCorrelationMatrix) * paramsVector;
cfc19dd5 656
657 // Ad - m
6d81c2de 658 measGuessVector -= (*fgCurrentESDVector);
cfc19dd5 659
447c325d 660 TVectorD copy(measGuessVector);
cfc19dd5 661
662 // (Ad - m) W
6d81c2de 663 // this step can be optimized because currently only the diagonal elements of fgCorrelationCovarianceMatrix are used
0b4bfd98 664 // normal way is like this:
6d81c2de 665 // measGuessVector *= (*fgCorrelationCovarianceMatrix);
0b4bfd98 666 // optimized way like this:
6d81c2de 667 for (Int_t i=0; i<fgkMaxParams; ++i)
668 measGuessVector[i] *= (*fgCorrelationCovarianceMatrix)(i, i);
cfc19dd5 669
447c325d 670 //measGuessVector.Print();
cfc19dd5 671
447c325d 672 // (Ad - m) W (Ad - m)
0b4bfd98 673 // the factor 1e6 prevents that a very small number (measGuessVector[i]) is multiplied with a very
6d81c2de 674 // big number ((*fgCorrelationCovarianceMatrix)(i, i)) (see ApplyMinuitFit)
447c325d 675 Double_t chi2FromFit = measGuessVector * copy * 1e6;
0a173978 676
3602328d 677 chi2 = chi2FromFit + penaltyVal;
0a173978 678
447c325d 679 static Int_t callCount = 0;
680 if ((callCount++ % 10000) == 0)
681 printf("%d) %f %f --> %f\n", callCount, chi2FromFit, penaltyVal, chi2);
0a173978 682}
683
684//____________________________________________________________________
447c325d 685void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t createBigBin)
0a173978 686{
687 //
9ca6be09 688 // fills fCurrentESD, fCurrentCorrelation
689 // resets fMultiplicityESDCorrected
0a173978 690 //
691
0b4bfd98 692 Bool_t display = kFALSE;
693
0a173978 694 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
0b4bfd98 695
9ca6be09 696 fMultiplicityESDCorrected[correlationID]->Reset();
0b4bfd98 697 fMultiplicityESDCorrected[correlationID]->Sumw2();
0a173978 698
0b4bfd98 699 fCurrentESD = fMultiplicityESD[inputRange]->ProjectionY("fCurrentESD"); //"fCurrentESD", -1, -1, "e");
9ca6be09 700 fCurrentESD->Sumw2();
0a173978 701
702 // empty under/overflow bins in x, otherwise Project3D takes them into account
703 TH3* hist = fCorrelation[correlationID];
704 for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
705 {
706 for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
707 {
708 hist->SetBinContent(0, y, z, 0);
709 hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
710 }
711 }
712
9ca6be09 713 fCurrentCorrelation = hist->Project3D("zy");
447c325d 714 //((TH2*) fCurrentCorrelation)->Rebin2D(2, 1);
715 //fMultiplicityESDCorrected[correlationID]->Rebin(2);
9ca6be09 716 fCurrentCorrelation->Sumw2();
cfc19dd5 717
447c325d 718 if (createBigBin)
719 {
720 Int_t maxBin = 0;
721 for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i)
722 {
723 if (fCurrentESD->GetBinContent(i) <= 5)
724 {
725 maxBin = i;
726 break;
727 }
728 }
729
730 if (maxBin > 0)
731 {
0b4bfd98 732 TCanvas* canvas = 0;
447c325d 733
0b4bfd98 734 if (display)
735 {
736 canvas = new TCanvas("StatSolution", "StatSolution", 1000, 800);
737 canvas->Divide(2, 2);
738
739 canvas->cd(1);
740 fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
741 fCurrentESD->SetStats(kFALSE);
742 fCurrentESD->SetTitle(";measured multiplicity");
743 fCurrentESD->DrawCopy();
744 gPad->SetLogy();
745
746 canvas->cd(2);
747 fCurrentCorrelation->GetXaxis()->SetRangeUser(0, 250);
748 fCurrentCorrelation->GetYaxis()->SetRangeUser(0, 250);
749 fCurrentCorrelation->SetTitle(";true multiplicity;measured multiplicity");
750 fCurrentCorrelation->SetStats(kFALSE);
751 fCurrentCorrelation->DrawCopy("COLZ");
752 }
447c325d 753
754 printf("Bin limit in measured spectrum is %d.\n", maxBin);
755 fCurrentESD->SetBinContent(maxBin, fCurrentESD->Integral(maxBin, fCurrentESD->GetNbinsX()));
756 for (Int_t i=maxBin+1; i<=fCurrentESD->GetNbinsX(); ++i)
757 {
758 fCurrentESD->SetBinContent(i, 0);
759 fCurrentESD->SetBinError(i, 0);
760 }
761 // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
762 fCurrentESD->SetBinError(maxBin, TMath::Sqrt(fCurrentESD->GetBinContent(maxBin)));
763
764 printf("This bin has now %f +- %f entries\n", fCurrentESD->GetBinContent(maxBin), fCurrentESD->GetBinError(maxBin));
765
766 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
767 {
768 fCurrentCorrelation->SetBinContent(i, maxBin, fCurrentCorrelation->Integral(i, i, maxBin, fCurrentCorrelation->GetNbinsY()));
769 // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions???
770 fCurrentCorrelation->SetBinError(i, maxBin, TMath::Sqrt(fCurrentCorrelation->GetBinContent(i, maxBin)));
771
772 for (Int_t j=maxBin+1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
773 {
774 fCurrentCorrelation->SetBinContent(i, j, 0);
775 fCurrentCorrelation->SetBinError(i, j, 0);
776 }
777 }
778
779 printf("Adjusted correlation matrix!\n");
780
0b4bfd98 781 if (canvas)
782 {
783 canvas->cd(3);
784 fCurrentESD->DrawCopy();
785 gPad->SetLogy();
447c325d 786
0b4bfd98 787 canvas->cd(4);
788 fCurrentCorrelation->DrawCopy("COLZ");
789 }
447c325d 790 }
791 }
792
0b4bfd98 793#if 0 // does not help
794 // null bins with one entry
795 Int_t nNulledBins = 0;
796 for (Int_t x=1; x<=fCurrentCorrelation->GetXaxis()->GetNbins(); ++x)
797 for (Int_t y=1; y<=fCurrentCorrelation->GetYaxis()->GetNbins(); ++y)
798 {
799 if (fCurrentCorrelation->GetBinContent(x, y) == 1)
800 {
801 fCurrentCorrelation->SetBinContent(x, y, 0);
802 fCurrentCorrelation->SetBinError(x, y, 0);
447c325d 803
0b4bfd98 804 ++nNulledBins;
805 }
806 }
807 Printf("Nulled %d bins", nNulledBins);
808#endif
809
810 fCurrentEfficiency = GetEfficiency(inputRange, eventType);
811 //fCurrentEfficiency->Rebin(2);
812 //fCurrentEfficiency->Scale(0.5);
813}
814
815//____________________________________________________________________
816TH1* AliMultiplicityCorrection::GetEfficiency(Int_t inputRange, EventType eventType)
817{
818 //
819 // calculates efficiency for given event type
820 //
821
822 TH1* divisor = 0;
cfc19dd5 823 switch (eventType)
824 {
0b4bfd98 825 case kTrVtx : divisor = fMultiplicityVtx[inputRange]->ProjectionY("divisor", -1, -1, "e"); break;
826 case kMB: divisor = fMultiplicityMB[inputRange]->ProjectionY("divisor", -1, -1, "e"); break;
827 case kINEL: divisor = fMultiplicityINEL[inputRange]->ProjectionY("divisor", -1, -1, "e"); break;
cfc19dd5 828 }
0b4bfd98 829 TH1* eff = fMultiplicityVtx[inputRange]->ProjectionY("CurrentEfficiency", -1, -1, "e");
830 eff->Divide(eff, divisor, 1, 1, "B");
831 return eff;
9ca6be09 832}
833
834//____________________________________________________________________
447c325d 835void AliMultiplicityCorrection::SetRegularizationParameters(RegularizationType type, Float_t weight)
836{
6d81c2de 837 //
838 // sets the parameters for chi2 minimization
839 //
840
841 fgRegularizationType = type;
842 fgRegularizationWeight = weight;
447c325d 843
844 printf("AliMultiplicityCorrection::SetRegularizationParameters --> Regularization set to %d with weight %f\n", (Int_t) type, weight);
845}
846
847//____________________________________________________________________
6d81c2de 848void AliMultiplicityCorrection::SetBayesianParameters(Float_t smoothing, Int_t nIterations)
9ca6be09 849{
850 //
6d81c2de 851 // sets the parameters for Bayesian unfolding
9ca6be09 852 //
853
6d81c2de 854 fgBayesianSmoothing = smoothing;
855 fgBayesianIterations = nIterations;
9ca6be09 856
6d81c2de 857 printf("AliMultiplicityCorrection::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f\n", fgBayesianIterations, fgBayesianSmoothing);
858}
cfc19dd5 859
6d81c2de 860//____________________________________________________________________
861Int_t AliMultiplicityCorrection::UnfoldWithMinuit(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* result, Bool_t check)
862{
863 //
864 // implemenation of unfolding (internal function)
865 //
866 // unfolds <measured> using response from <correlation> and effiency <aEfficiency>
867 // output is in <result>
868 // <initialConditions> set the initial values for the minimization, if 0 <measured> is used
869 // if <check> is true no unfolding is made, instead only the chi2 without unfolding is printed
870 //
871 // returns minuit status (0 = success), (-1 when check was set)
872 //
447c325d 873
6d81c2de 874 //normalize measured
875 measured->Scale(1.0 / measured->Integral());
876
877 if (!fgCorrelationMatrix)
878 fgCorrelationMatrix = new TMatrixD(fgkMaxInput, fgkMaxParams);
879 if (!fgCorrelationCovarianceMatrix)
880 fgCorrelationCovarianceMatrix = new TMatrixD(fgkMaxInput, fgkMaxInput);
881 if (!fgCurrentESDVector)
882 fgCurrentESDVector = new TVectorD(fgkMaxInput);
883 if (!fgEntropyAPriori)
884 fgEntropyAPriori = new TVectorD(fgkMaxParams);
885
886 fgCorrelationMatrix->Zero();
887 fgCorrelationCovarianceMatrix->Zero();
888 fgCurrentESDVector->Zero();
889 fgEntropyAPriori->Zero();
0a173978 890
891 // normalize correction for given nPart
6d81c2de 892 for (Int_t i=1; i<=correlation->GetNbinsX(); ++i)
0a173978 893 {
6d81c2de 894 Double_t sum = correlation->Integral(i, i, 1, correlation->GetNbinsY());
0a173978 895 if (sum <= 0)
896 continue;
9ca6be09 897 Float_t maxValue = 0;
898 Int_t maxBin = -1;
6d81c2de 899 for (Int_t j=1; j<=correlation->GetNbinsY(); ++j)
0a173978 900 {
9ca6be09 901 // find most probably value
6d81c2de 902 if (maxValue < correlation->GetBinContent(i, j))
9ca6be09 903 {
6d81c2de 904 maxValue = correlation->GetBinContent(i, j);
9ca6be09 905 maxBin = j;
906 }
907
0a173978 908 // npart sum to 1
6d81c2de 909 correlation->SetBinContent(i, j, correlation->GetBinContent(i, j) / sum);
910 correlation->SetBinError(i, j, correlation->GetBinError(i, j) / sum);
cfc19dd5 911
6d81c2de 912 if (i <= fgkMaxParams && j <= fgkMaxInput)
913 (*fgCorrelationMatrix)(j-1, i-1) = correlation->GetBinContent(i, j);
0a173978 914 }
9ca6be09 915
cfc19dd5 916 //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
0a173978 917 }
918
0a173978 919 // Initialize TMinuit via generic fitter interface
6d81c2de 920 TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgkMaxParams);
447c325d 921 Double_t arglist[100];
0b4bfd98 922
6d81c2de 923 // disable any output (-1), unfortuantly we do not see warnings anymore then. Have to find another way...
447c325d 924 arglist[0] = 0;
925 minuit->ExecuteCommand("SET PRINT", arglist, 1);
0a173978 926
0b4bfd98 927 // however, enable warnings
928 //minuit->ExecuteCommand("SET WAR", arglist, 0);
929
930 // set minimization function
0a173978 931 minuit->SetFCN(MinuitFitFunction);
932
6d81c2de 933 for (Int_t i=0; i<fgkMaxParams; i++)
934 (*fgEntropyAPriori)[i] = 1;
0a173978 935
6d81c2de 936 if (initialConditions)
dd701109 937 {
938 printf("Using different starting conditions...\n");
6d81c2de 939 //new TCanvas; initialConditions->DrawCopy();
447c325d 940
6d81c2de 941 initialConditions->Scale(1.0 / initialConditions->Integral());
942 for (Int_t i=0; i<fgkMaxParams; i++)
943 if (initialConditions->GetBinContent(i+1) > 0)
944 (*fgEntropyAPriori)[i] = initialConditions->GetBinContent(i+1);
dd701109 945 }
946 else
6d81c2de 947 initialConditions = measured;
447c325d 948
6d81c2de 949 Double_t results[fgkMaxParams+1];
950 for (Int_t i=0; i<fgkMaxParams; ++i)
0a173978 951 {
6d81c2de 952 results[i] = initialConditions->GetBinContent(i+1);
dd701109 953
447c325d 954 // minimum value
955 results[i] = TMath::Max(results[i], 1e-3);
956
957 Float_t stepSize = 0.1;
958
959 // minuit sees squared values to prevent it from going negative...
960 results[i] = TMath::Sqrt(results[i]);
961 //stepSize /= results[i]; // keep relative step size
962
963 minuit->SetParameter(i, Form("param%d", i), results[i], stepSize, 0, 0);
964 }
6d81c2de 965 //minuit->SetParameter(fgkMaxParams, "alpha", 1e4, 1, 0, 0);
447c325d 966 // bin 0 is filled with value from bin 1 (otherwise it's 0)
967 //minuit->SetParameter(0, "param0", (results[1] > minStartValue) ? results[1] : minStartValue, 0.1, 0, 0);
968 //results[0] = 0;
969 //minuit->SetParameter(0, "param0", 0, 0, 0, 0);
970
6d81c2de 971 for (Int_t i=0; i<fgkMaxInput; ++i)
447c325d 972 {
6d81c2de 973 (*fgCurrentESDVector)[i] = measured->GetBinContent(i+1);
974 if (measured->GetBinError(i+1) > 0)
975 (*fgCorrelationCovarianceMatrix)(i, i) = (Double_t) 1e-6 / measured->GetBinError(i+1) / measured->GetBinError(i+1);
cfc19dd5 976
6d81c2de 977 if ((*fgCorrelationCovarianceMatrix)(i, i) > 1e7)
978 (*fgCorrelationCovarianceMatrix)(i, i) = 0;
0a173978 979 }
980
981 Int_t dummy;
982 Double_t chi2;
983 MinuitFitFunction(dummy, 0, chi2, results, 0);
dd701109 984 printf("Chi2 of initial parameters is = %f\n", chi2);
0a173978 985
9ca6be09 986 if (check)
447c325d 987 return -1;
988
989 // first param is number of iterations, second is precision....
990 arglist[0] = 1e6;
991 //arglist[1] = 1e-5;
992 //minuit->ExecuteCommand("SCAN", arglist, 0);
993 Int_t status = minuit->ExecuteCommand("MIGRAD", arglist, 1);
994 printf("MINUIT status is %d\n", status);
995 //minuit->ExecuteCommand("MIGRAD", arglist, 1);
996 //minuit->ExecuteCommand("MIGRAD", arglist, 1);
dd701109 997 //printf("!!!!!!!!!!!!!! MIGRAD finished: Starting MINOS !!!!!!!!!!!!!!");
0a173978 998 //minuit->ExecuteCommand("MINOS", arglist, 0);
999
6d81c2de 1000 for (Int_t i=0; i<fgkMaxParams; ++i)
0a173978 1001 {
6d81c2de 1002 if (aEfficiency->GetBinContent(i+1) > 0)
dd701109 1003 {
6d81c2de 1004 result->SetBinContent(i+1, minuit->GetParameter(i) * minuit->GetParameter(i) / aEfficiency->GetBinContent(i+1));
447c325d 1005 // error is : (relError) * (value) = (minuit->GetParError(i) / minuit->GetParameter(i)) * (minuit->GetParameter(i) * minuit->GetParameter(i))
6d81c2de 1006 result->SetBinError(i+1, minuit->GetParError(i) * minuit->GetParameter(i) / aEfficiency->GetBinContent(i+1));
dd701109 1007 }
1008 else
1009 {
6d81c2de 1010 result->SetBinContent(i+1, 0);
1011 result->SetBinError(i+1, 0);
dd701109 1012 }
1013 }
447c325d 1014
1015 return status;
dd701109 1016}
1017
6d81c2de 1018//____________________________________________________________________
1019Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t check, TH1* initialConditions)
1020{
1021 //
1022 // correct spectrum using minuit chi2 method
1023 //
1024 // for description of parameters, see UnfoldWithMinuit
1025 //
1026
1027 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1028 SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kTRUE);
1029
1030 return UnfoldWithMinuit(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], check);
1031}
1032
dd701109 1033//____________________________________________________________________
1034void AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace)
1035{
1036 //
1037 // correct spectrum using minuit chi2 method applying a NBD fit
1038 //
1039
1040 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1041
447c325d 1042 SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE);
0b4bfd98 1043 //normalize ESD
1044 fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
dd701109 1045
6d81c2de 1046 fgCorrelationMatrix = new TMatrixD(fgkMaxParams, fgkMaxParams);
1047 fgCorrelationCovarianceMatrix = new TMatrixD(fgkMaxParams, fgkMaxParams);
1048 fgCurrentESDVector = new TVectorD(fgkMaxParams);
dd701109 1049
1050 // normalize correction for given nPart
1051 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1052 {
1053 Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1054 if (sum <= 0)
1055 continue;
1056 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1057 {
1058 // npart sum to 1
1059 fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1060 fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1061
6d81c2de 1062 if (i <= fgkMaxParams && j <= fgkMaxParams)
1063 (*fgCorrelationMatrix)(j-1, i-1) = fCurrentCorrelation->GetBinContent(i, j);
dd701109 1064 }
1065 }
1066
6d81c2de 1067 for (Int_t i=0; i<fgkMaxParams; ++i)
dd701109 1068 {
6d81c2de 1069 (*fgCurrentESDVector)[i] = fCurrentESD->GetBinContent(i+1);
dd701109 1070 if (fCurrentESD->GetBinError(i+1) > 0)
6d81c2de 1071 (*fgCorrelationCovarianceMatrix)(i, i) = 1.0 / (fCurrentESD->GetBinError(i+1) * fCurrentESD->GetBinError(i+1));
dd701109 1072 }
1073
1074 // Create NBD function
6d81c2de 1075 if (!fgNBD)
dd701109 1076 {
6d81c2de 1077 fgNBD = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 250);
1078 fgNBD->SetParNames("scaling", "averagen", "k");
0a173978 1079 }
1080
dd701109 1081 // Initialize TMinuit via generic fitter interface
1082 TVirtualFitter *minuit = TVirtualFitter::Fitter(0, 3);
1083
1084 minuit->SetFCN(MinuitNBD);
1085 minuit->SetParameter(0, "param0", 1, 0.1, 0.001, fCurrentESD->GetMaximum() * 1000);
1086 minuit->SetParameter(1, "param1", 10, 1, 0.001, 1000);
1087 minuit->SetParameter(2, "param2", 2, 1, 0.001, 1000);
1088
1089 Double_t arglist[100];
1090 arglist[0] = 0;
1091 minuit->ExecuteCommand("SET PRINT", arglist, 1);
447c325d 1092 minuit->ExecuteCommand("MIGRAD", arglist, 0);
dd701109 1093 //minuit->ExecuteCommand("MINOS", arglist, 0);
1094
6d81c2de 1095 fgNBD->SetParameters(minuit->GetParameter(0), minuit->GetParameter(1), minuit->GetParameter(2));
dd701109 1096
6d81c2de 1097 for (Int_t i=0; i<fgkMaxParams; ++i)
dd701109 1098 {
6d81c2de 1099 printf("%d : %f\n", i, fgNBD->Eval(i));
1100 if (fgNBD->Eval(i) > 0)
dd701109 1101 {
6d81c2de 1102 fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, fgNBD->Eval(i));
dd701109 1103 fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0);
1104 }
1105 }
0a173978 1106}
1107
1108//____________________________________________________________________
1109void AliMultiplicityCorrection::NormalizeToBinWidth(TH1* hist)
1110{
1111 //
1112 // normalizes a 1-d histogram to its bin width
1113 //
1114
1115 for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
1116 {
1117 hist->SetBinContent(i, hist->GetBinContent(i) / hist->GetBinWidth(i));
1118 hist->SetBinError(i, hist->GetBinError(i) / hist->GetBinWidth(i));
1119 }
1120}
1121
1122//____________________________________________________________________
1123void AliMultiplicityCorrection::NormalizeToBinWidth(TH2* hist)
1124{
1125 //
1126 // normalizes a 2-d histogram to its bin width (x width * y width)
1127 //
1128
1129 for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
1130 for (Int_t j=1; j<=hist->GetNbinsY(); ++j)
1131 {
1132 Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j);
1133 hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor);
1134 hist->SetBinError(i, j, hist->GetBinError(i, j) / factor);
1135 }
1136}
1137
0a173978 1138//____________________________________________________________________
1139void AliMultiplicityCorrection::DrawHistograms()
1140{
6d81c2de 1141 //
1142 // draws the histograms of this class
1143 //
1144
0a173978 1145 printf("ESD:\n");
1146
1147 TCanvas* canvas1 = new TCanvas("fMultiplicityESD", "fMultiplicityESD", 900, 600);
1148 canvas1->Divide(3, 2);
1149 for (Int_t i = 0; i < kESDHists; ++i)
1150 {
1151 canvas1->cd(i+1);
1152 fMultiplicityESD[i]->DrawCopy("COLZ");
1153 printf("%d --> %f\n", i, (Float_t) fMultiplicityESD[i]->ProjectionY()->GetMean());
1154 }
1155
cfc19dd5 1156 printf("Vtx:\n");
0a173978 1157
1158 TCanvas* canvas2 = new TCanvas("fMultiplicityMC", "fMultiplicityMC", 900, 600);
1159 canvas2->Divide(3, 2);
1160 for (Int_t i = 0; i < kMCHists; ++i)
1161 {
1162 canvas2->cd(i+1);
cfc19dd5 1163 fMultiplicityVtx[i]->DrawCopy("COLZ");
1164 printf("%d --> %f\n", i, fMultiplicityVtx[i]->ProjectionY()->GetMean());
0a173978 1165 }
1166
1167 TCanvas* canvas3 = new TCanvas("fCorrelation", "fCorrelation", 900, 900);
1168 canvas3->Divide(3, 3);
1169 for (Int_t i = 0; i < kCorrHists; ++i)
1170 {
1171 canvas3->cd(i+1);
b477f4f2 1172 TH3* hist = dynamic_cast<TH3*> (fCorrelation[i]->Clone());
1173 for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
1174 {
1175 for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
1176 {
1177 hist->SetBinContent(0, y, z, 0);
1178 hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
1179 }
1180 }
1181 TH1* proj = hist->Project3D("zy");
0a173978 1182 proj->DrawCopy("COLZ");
1183 }
1184}
1185
1186//____________________________________________________________________
447c325d 1187void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t normalizeESD, TH1* mcHist, Bool_t simple)
0a173978 1188{
447c325d 1189 //mcHist->Rebin(2);
1190
cfc19dd5 1191 Int_t esdCorrId = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1192
9ca6be09 1193 TString tmpStr;
cfc19dd5 1194 tmpStr.Form("%s_DrawComparison_%d", name, esdCorrId);
0a173978 1195
447c325d 1196 if (fMultiplicityESDCorrected[esdCorrId]->Integral() == 0)
1197 {
1198 printf("ERROR. Unfolded histogram is empty\n");
1199 return;
1200 }
1201
1202 //regain measured distribution used for unfolding, because the bins at high mult. were modified in SetupCurrentHists
1203 fCurrentESD = fMultiplicityESD[esdCorrId]->ProjectionY();
1204 fCurrentESD->Sumw2();
1205 fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
1206
1207 // normalize unfolded result to 1
1208 fMultiplicityESDCorrected[esdCorrId]->Scale(1.0 / fMultiplicityESDCorrected[esdCorrId]->Integral());
1209
1210 //fCurrentESD->Scale(mcHist->Integral(2, 200));
1211
1212 //new TCanvas;
1213 /*TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
1214 ratio->Divide(mcHist);
1215 ratio->Draw("HIST");
1216 ratio->Fit("pol0", "W0", "", 20, 120);
1217 Float_t scalingFactor = ratio->GetFunction("pol0")->GetParameter(0);
1218 delete ratio;
1219 mcHist->Scale(scalingFactor);*/
1220
0b4bfd98 1221 // find bin with <= 100 entries. this is used as limit for the chi2 calculation
1222 Int_t mcBinLimit = 0;
1223 for (Int_t i=20; i<mcHist->GetNbinsX(); ++i)
1224 {
1225 if (mcHist->GetBinContent(i) > 50)
1226 {
1227 mcBinLimit = i;
1228 }
1229 else
1230 break;
1231 }
1232 Printf("AliMultiplicityCorrection::DrawComparison: MC bin limit is %d", mcBinLimit);
1233
1234 // scale to 1
1235 mcHist->Sumw2();
447c325d 1236 mcHist->Scale(1.0 / mcHist->Integral());
1237
b477f4f2 1238 // calculate residual
1239
1240 // for that we convolute the response matrix with the gathered result
1241 // if normalizeESD is not set, the histogram is already normalized, this needs to be passed to CalculateMultiplicityESD
dd701109 1242 TH1* tmpESDEfficiencyRecorrected = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("tmpESDEfficiencyRecorrected");
1243 tmpESDEfficiencyRecorrected->Multiply(fCurrentEfficiency);
447c325d 1244 TH2* convoluted = CalculateMultiplicityESD(tmpESDEfficiencyRecorrected, esdCorrId);
1245 TH1* convolutedProj = convoluted->ProjectionY("convolutedProj", -1, -1, "e");
1246 if (convolutedProj->Integral() > 0)
1247 {
1248 convolutedProj->Scale(1.0 / convolutedProj->Integral());
1249 }
1250 else
1251 printf("ERROR: convolutedProj is empty. Something went wrong calculating the convoluted histogram.\n");
1252
447c325d 1253 TH1* residual = (TH1*) convolutedProj->Clone("residual");
1254 residual->SetTitle("Residuals;Ntracks;(folded unfolded measured - measured) / e");
b477f4f2 1255
1256 residual->Add(fCurrentESD, -1);
447c325d 1257 //residual->Divide(residual, fCurrentESD, 1, 1, "B");
1258
0b4bfd98 1259 TH1* residualHist = new TH1F("residualHist", "residualHist", 51, -5, 5);
b477f4f2 1260
1261 // TODO fix errors
447c325d 1262 Float_t chi2 = 0;
1263 for (Int_t i=1; i<=residual->GetNbinsX(); ++i)
b477f4f2 1264 {
447c325d 1265 if (fCurrentESD->GetBinError(i) > 0)
1266 {
1267 Float_t value = residual->GetBinContent(i) / fCurrentESD->GetBinError(i);
1268 if (i > 1)
1269 chi2 += value * value;
1270 residual->SetBinContent(i, value);
1271 residualHist->Fill(value);
1272 }
1273 else
1274 {
1275 //printf("Residual bin %d set to 0\n", i);
1276 residual->SetBinContent(i, 0);
1277 }
1278 convolutedProj->SetBinError(i, 0);
b477f4f2 1279 residual->SetBinError(i, 0);
b477f4f2 1280
447c325d 1281 if (i == 200)
1282 fLastChi2Residuals = chi2;
1283 }
dd701109 1284
0b4bfd98 1285 new TCanvas; residualHist->DrawCopy();
1286
1287 //residualHist->Fit("gaus", "N");
1288 //delete residualHist;
dd701109 1289
447c325d 1290 printf("Difference (Residuals) is %f for bin 2-200\n", chi2);
dd701109 1291
447c325d 1292 TCanvas* canvas1 = 0;
1293 if (simple)
1294 {
1295 canvas1 = new TCanvas(tmpStr, tmpStr, 900, 400);
1296 canvas1->Divide(2, 1);
1297 }
1298 else
1299 {
1300 canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 1200);
1301 canvas1->Divide(2, 3);
1302 }
0a173978 1303
1304 canvas1->cd(1);
cfc19dd5 1305 TH1* proj = (TH1*) mcHist->Clone("proj");
9ca6be09 1306
dd701109 1307 proj->GetXaxis()->SetRangeUser(0, 200);
447c325d 1308 proj->SetTitle(";true multiplicity;Entries");
1309 proj->SetStats(kFALSE);
cfc19dd5 1310 proj->DrawCopy("HIST");
0a173978 1311 gPad->SetLogy();
1312
cfc19dd5 1313 //fMultiplicityESDCorrected[esdCorrId]->SetMarkerStyle(3);
1314 fMultiplicityESDCorrected[esdCorrId]->SetLineColor(2);
447c325d 1315 fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAME HIST E");
cfc19dd5 1316
447c325d 1317 TLegend* legend = new TLegend(0.3, 0.8, 0.93, 0.93);
1318 legend->AddEntry(proj, "true distribution");
1319 legend->AddEntry(fMultiplicityESDCorrected[esdCorrId], "unfolded distribution");
1320 legend->SetFillColor(0);
1321 legend->Draw();
1322 // unfortunately does not work. maybe a bug? --> legend->SetTextSizePixels(14);
1323
1324 canvas1->cd(2);
1325
1326 gPad->SetLogy();
1327 fCurrentESD->GetXaxis()->SetRangeUser(0, 200);
1328 //fCurrentESD->SetLineColor(2);
1329 fCurrentESD->SetTitle(";measured multiplicity;Entries");
1330 fCurrentESD->SetStats(kFALSE);
1331 fCurrentESD->DrawCopy("HIST E");
1332
1333 convolutedProj->SetLineColor(2);
1334 //proj2->SetMarkerColor(2);
1335 //proj2->SetMarkerStyle(5);
1336 convolutedProj->DrawCopy("HIST SAME");
1337
1338 legend = new TLegend(0.3, 0.8, 0.93, 0.93);
1339 legend->AddEntry(fCurrentESD, "measured distribution");
1340 legend->AddEntry(convolutedProj, "R #otimes unfolded distribution");
1341 legend->SetFillColor(0);
1342 legend->Draw();
1343
0b4bfd98 1344 //TH1* diffMCUnfolded = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded"));
1345 //diffMCUnfolded->Add(fMultiplicityESDCorrected[esdCorrId], -1);
447c325d 1346
1347 /*Float_t chi2 = 0;
1348 Float_t chi = 0;
dd701109 1349 fLastChi2MCLimit = 0;
447c325d 1350 Int_t limit = 0;
dd701109 1351 for (Int_t i=2; i<=200; ++i)
cfc19dd5 1352 {
dd701109 1353 if (proj->GetBinContent(i) != 0)
cfc19dd5 1354 {
dd701109 1355 Float_t value = (proj->GetBinContent(i) - fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i)) / proj->GetBinContent(i);
cfc19dd5 1356 chi2 += value * value;
447c325d 1357 chi += TMath::Abs(value);
dd701109 1358
447c325d 1359 //printf("%d %f\n", i, chi);
cfc19dd5 1360
447c325d 1361 if (chi2 < 0.2)
1362 fLastChi2MCLimit = i;
cfc19dd5 1363
447c325d 1364 if (chi < 3)
1365 limit = i;
cfc19dd5 1366
447c325d 1367 }
1368 }*/
9ca6be09 1369
0b4bfd98 1370 /*chi2 = 0;
447c325d 1371 Float_t chi = 0;
1372 Int_t limit = 0;
1373 for (Int_t i=1; i<=diffMCUnfolded->GetNbinsX(); ++i)
cfc19dd5 1374 {
447c325d 1375 if (fMultiplicityESDCorrected[esdCorrId]->GetBinError(i) > 0)
cfc19dd5 1376 {
447c325d 1377 Double_t value = diffMCUnfolded->GetBinContent(i) / fMultiplicityESDCorrected[esdCorrId]->GetBinError(i);
1378 if (value > 1e8)
1379 value = 1e8; //prevent arithmetic exception
1380 else if (value < -1e8)
1381 value = -1e8;
1382 if (i > 1)
1383 {
1384 chi2 += value * value;
1385 chi += TMath::Abs(value);
1386 }
1387 diffMCUnfolded->SetBinContent(i, value);
cfc19dd5 1388 }
447c325d 1389 else
1390 {
1391 //printf("diffMCUnfolded bin %d set to 0\n", i);
1392 diffMCUnfolded->SetBinContent(i, 0);
1393 }
1394 if (chi2 < 1000)
1395 fLastChi2MCLimit = i;
1396 if (chi < 1000)
1397 limit = i;
1398 if (i == 150)
1399 fLastChi2MC = chi2;
cfc19dd5 1400 }
0a173978 1401
447c325d 1402 printf("limits %d %d\n", fLastChi2MCLimit, limit);
1403 fLastChi2MCLimit = limit;
1404
0b4bfd98 1405 printf("Difference (from MC) is %f for bin 2-150. Limit is %d.\n", fLastChi2MC, fLastChi2MCLimit);*/
0a173978 1406
447c325d 1407 if (!simple)
1408 {
0b4bfd98 1409 /*canvas1->cd(3);
447c325d 1410
0b4bfd98 1411 diffMCUnfolded->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e(unfolded)");
447c325d 1412 //diffMCUnfolded->GetYaxis()->SetRangeUser(-20, 20);
1413 diffMCUnfolded->GetXaxis()->SetRangeUser(0, 200);
1414 diffMCUnfolded->DrawCopy("HIST");
0a173978 1415
447c325d 1416 TH1F* fluctuation = new TH1F("fluctuation", "fluctuation", 20, -5, 5);
1417 for (Int_t i=20; i<=diffMCUnfolded->GetNbinsX(); ++i)
1418 fluctuation->Fill(diffMCUnfolded->GetBinContent(i));
3602328d 1419
0b4bfd98 1420 //new TCanvas; fluctuation->DrawCopy();
1421 delete fluctuation;*/
447c325d 1422
1423 /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
1424 legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
1425 legend->AddEntry(fMultiplicityMC, "MC");
1426 legend->AddEntry(fMultiplicityESD, "ESD");
1427 legend->Draw();*/
1428
1429 canvas1->cd(4);
1430 //residual->GetYaxis()->SetRangeUser(-0.2, 0.2);
1431 residual->GetXaxis()->SetRangeUser(0, 200);
1432 residual->DrawCopy();
1433
1434 canvas1->cd(5);
1435
1436 TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
1437 ratio->Divide(mcHist);
1438 ratio->SetTitle("Ratio;true multiplicity;Unfolded / MC");
1439 ratio->GetYaxis()->SetRangeUser(0.5, 1.5);
1440 ratio->GetXaxis()->SetRangeUser(0, 200);
1441 ratio->SetStats(kFALSE);
1442 ratio->Draw("HIST");
1443
1444 Double_t ratioChi2 = 0;
0b4bfd98 1445 fRatioAverage = 0;
447c325d 1446 fLastChi2MCLimit = 0;
1447 for (Int_t i=2; i<=150; ++i)
1448 {
1449 Float_t value = ratio->GetBinContent(i) - 1;
1450 if (value > 1e8)
1451 value = 1e8; //prevent arithmetic exception
1452 else if (value < -1e8)
1453 value = -1e8;
1454
1455 ratioChi2 += value * value;
0b4bfd98 1456 fRatioAverage += TMath::Abs(value);
447c325d 1457
1458 if (ratioChi2 < 0.1)
1459 fLastChi2MCLimit = i;
1460 }
0b4bfd98 1461 fRatioAverage /= 149;
447c325d 1462
0b4bfd98 1463 printf("Sum over (ratio-1)^2 (2..150) is %f; average of |ratio-1| is %f\n", ratioChi2, fRatioAverage);
447c325d 1464 // TODO FAKE
1465 fLastChi2MC = ratioChi2;
0b4bfd98 1466
1467 // FFT of ratio
1468 canvas1->cd(6);
1469 const Int_t kFFT = 128;
1470 Double_t fftReal[kFFT];
1471 Double_t fftImag[kFFT];
1472
1473 for (Int_t i=0; i<kFFT; ++i)
1474 {
1475 fftReal[i] = ratio->GetBinContent(i+1+10);
1476 // test: ;-)
1477 //fftReal[i] = cos(TMath::Pi() * 5 * 2 * i / 128);
1478 fftImag[i] = 0;
1479 }
1480
1481 FFT(-1, TMath::Nint(log(kFFT) / log(2)), fftReal, fftImag);
1482
1483 TH1* fftResult = (TH1*) ratio->Clone("fftResult");
1484 fftResult->SetTitle("FFT;true multiplicity;coeff. (10...137)");
1485 TH1* fftResult2 = (TH1*) ratio->Clone("fftResult2");
1486 TH1* fftResult3 = (TH1*) ratio->Clone("fftResult3");
1487 fftResult->Reset();
1488 fftResult2->Reset();
1489 fftResult3->Reset();
1490 fftResult->GetYaxis()->UnZoom();
1491 fftResult2->GetYaxis()->UnZoom();
1492 fftResult3->GetYaxis()->UnZoom();
1493
1494 Printf("AFTER FFT");
1495 for (Int_t i=0; i<kFFT; ++i)
1496 {
1497 Printf("%d: %f", i, fftReal[i]);
1498 fftResult->SetBinContent(i+1, fftReal[i]);
1499 /*if (i != 0 && TMath::Abs(fftReal[i]) > 0.5)
1500 {
1501 Printf("Nulled %d", i);
1502 fftReal[i] = 0;
1503 }*/
1504 }
1505
1506 fftResult->SetLineColor(1);
1507 fftResult->DrawCopy();
1508
1509
1510 Printf("IMAG");
1511 for (Int_t i=0; i<kFFT; ++i)
1512 {
1513 Printf("%d: %f", i, fftImag[i]);
1514 fftResult2->SetBinContent(i+1, fftImag[i]);
1515
1516 fftResult3->SetBinContent(i+1, TMath::Sqrt(fftReal[i] * fftReal[i] + fftImag[i] * fftImag[i]));
1517 }
1518
1519 fftResult2->SetLineColor(2);
1520 fftResult2->DrawCopy("SAME");
1521
1522 fftResult3->SetLineColor(4);
1523 fftResult3->DrawCopy("SAME");
1524
1525 for (Int_t i=1; i<kFFT - 1; ++i)
1526 {
1527 if (TMath::Sqrt(fftReal[i] * fftReal[i] + fftImag[i] * fftImag[i]) > 3)
1528 {
1529 fftReal[i-1] = 0;
1530 fftReal[i] = 0;
1531 fftReal[i+1] = 0;
1532 fftImag[i-1] = 0;
1533 fftImag[i] = 0;
1534 fftImag[i+1] = 0;
1535 //fftReal[i] = (fftReal[i-1] + fftReal[i+1]) / 2;
1536 //fftImag[i] = (fftImag[i-1] + fftImag[i+1]) / 2;
1537 Printf("Nulled %d to %f %f", i, fftReal[i], fftImag[i]);
1538 }
1539 }
1540
1541 FFT(1, TMath::Nint(log(kFFT) / log(2)), fftReal, fftImag);
1542
1543 TH1* fftResult4 = (TH1*) fftResult3->Clone("fftResult4");
1544 fftResult4->Reset();
1545
1546 Printf("AFTER BACK-TRAFO");
1547 for (Int_t i=0; i<kFFT; ++i)
1548 {
1549 Printf("%d: %f", i, fftReal[i]);
1550 fftResult4->SetBinContent(i+1+10, fftReal[i]);
1551 }
1552
1553 canvas1->cd(5);
1554 fftResult4->SetLineColor(4);
1555 fftResult4->DrawCopy("SAME");
1556
1557 // plot (MC - Unfolded) / error (MC)
1558 canvas1->cd(3);
1559
1560 TH1* diffMCUnfolded2 = dynamic_cast<TH1*> (proj->Clone("diffMCUnfolded2"));
1561 diffMCUnfolded2->Add(fMultiplicityESDCorrected[esdCorrId], -1);
1562
1563 Int_t ndfQual[kQualityRegions];
1564 for (Int_t region=0; region<kQualityRegions; ++region)
1565 {
1566 fQuality[region] = 0;
1567 ndfQual[region] = 0;
1568 }
1569
1570
1571 Double_t newChi2 = 0;
6d81c2de 1572 Double_t newChi2Limit150 = 0;
0b4bfd98 1573 Int_t ndf = 0;
6d81c2de 1574 for (Int_t i=1; i<=TMath::Min(diffMCUnfolded2->GetNbinsX(), fgkMaxParams+1); ++i)
0b4bfd98 1575 {
1576 Double_t value = 0;
1577 if (proj->GetBinError(i) > 0)
1578 {
1579 value = diffMCUnfolded2->GetBinContent(i) / proj->GetBinError(i);
1580 newChi2 += value * value;
1581 if (i > 1 && i <= mcBinLimit)
6d81c2de 1582 newChi2Limit150 += value * value;
0b4bfd98 1583 ++ndf;
1584
1585 for (Int_t region=0; region<kQualityRegions; ++region)
1586 if (diffMCUnfolded2->GetXaxis()->GetBinCenter(i) >= fgQualityRegionsB[region] - 0.1 && diffMCUnfolded2->GetXaxis()->GetBinCenter(i) <= fgQualityRegionsE[region] + 0.1) // 0.1 to avoid e.g. 3.9999 < 4 problem
1587 {
1588 fQuality[region] += TMath::Abs(value);
1589 ++ndfQual[region];
1590 }
1591 }
1592
1593 diffMCUnfolded2->SetBinContent(i, value);
1594 }
1595
1596 // normalize region to the number of entries
1597 for (Int_t region=0; region<kQualityRegions; ++region)
1598 {
1599 if (ndfQual[region] > 0)
1600 fQuality[region] /= ndfQual[region];
1601 Printf("Quality parameter %d (%d <= mult <= %d) is %f with %d df", region, fgQualityRegionsB[region], fgQualityRegionsE[region], fQuality[region], ndfQual[region]);
1602 }
1603
1604 if (mcBinLimit > 1)
1605 {
1606 // TODO another FAKE
6d81c2de 1607 fLastChi2MC = newChi2Limit150 / (mcBinLimit - 1);
1608 Printf("Chi2 (2..%d) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", mcBinLimit, newChi2Limit150, mcBinLimit - 1, fLastChi2MC);
0b4bfd98 1609 }
1610 else
1611 fLastChi2MC = -1;
1612
1613 Printf("Chi2 (full range) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", newChi2, ndf, newChi2 / ndf);
1614
1615 diffMCUnfolded2->SetTitle("#chi^{2};Npart;(MC - Unfolded) / e(MC)");
1616 //diffMCUnfolded->GetYaxis()->SetRangeUser(-20, 20);
1617 diffMCUnfolded2->GetXaxis()->SetRangeUser(0, 200);
1618 diffMCUnfolded2->DrawCopy("HIST");
447c325d 1619 }
3602328d 1620
b477f4f2 1621 canvas1->SaveAs(Form("%s.gif", canvas1->GetName()));
0a173978 1622}
1623
1624//____________________________________________________________________
0b4bfd98 1625void AliMultiplicityCorrection::FFT(Int_t dir, Int_t m, Double_t *x, Double_t *y)
1626{
1627/*-------------------------------------------------------------------------
1628 This computes an in-place complex-to-complex FFT
1629 x and y are the real and imaginary arrays of 2^m points.
1630 dir = 1 gives forward transform
1631 dir = -1 gives reverse transform
1632
1633 Formula: forward
1634 N-1
1635 ---
1636 1 \ - j k 2 pi n / N
1637 X(n) = --- > x(k) e = forward transform
1638 N / n=0..N-1
1639 ---
1640 k=0
1641
1642 Formula: reverse
1643 N-1
1644 ---
1645 \ j k 2 pi n / N
1646 X(n) = > x(k) e = forward transform
1647 / n=0..N-1
1648 ---
1649 k=0
1650*/
1651
1652 Long_t nn, i, i1, j, k, i2, l, l1, l2;
1653 Double_t c1, c2, tx, ty, t1, t2, u1, u2, z;
1654
1655 /* Calculate the number of points */
1656 nn = 1;
1657 for (i = 0; i < m; i++)
1658 nn *= 2;
1659
1660 /* Do the bit reversal */
1661 i2 = nn >> 1;
1662 j = 0;
1663 for (i= 0; i < nn - 1; i++) {
1664 if (i < j) {
1665 tx = x[i];
1666 ty = y[i];
1667 x[i] = x[j];
1668 y[i] = y[j];
1669 x[j] = tx;
1670 y[j] = ty;
1671 }
1672 k = i2;
1673 while (k <= j) {
1674 j -= k;
1675 k >>= 1;
1676 }
1677 j += k;
1678 }
1679
1680 /* Compute the FFT */
1681 c1 = -1.0;
1682 c2 = 0.0;
1683 l2 = 1;
1684 for (l = 0; l < m; l++) {
1685 l1 = l2;
1686 l2 <<= 1;
1687 u1 = 1.0;
1688 u2 = 0.0;
1689 for (j = 0;j < l1; j++) {
1690 for (i = j; i < nn; i += l2) {
1691 i1 = i + l1;
1692 t1 = u1 * x[i1] - u2 * y[i1];
1693 t2 = u1 * y[i1] + u2 * x[i1];
1694 x[i1] = x[i] - t1;
1695 y[i1] = y[i] - t2;
1696 x[i] += t1;
1697 y[i] += t2;
1698 }
1699 z = u1 * c1 - u2 * c2;
1700 u2 = u1 * c2 + u2 * c1;
1701 u1 = z;
1702 }
1703 c2 = TMath::Sqrt((1.0 - c1) / 2.0);
1704 if (dir == 1)
1705 c2 = -c2;
1706 c1 = TMath::Sqrt((1.0 + c1) / 2.0);
1707 }
1708
1709 /* Scaling for forward transform */
1710 if (dir == 1) {
1711 for (i=0;i<nn;i++) {
1712 x[i] /= (Double_t)nn;
1713 y[i] /= (Double_t)nn;
1714 }
1715 }
1716}
1717
1718//____________________________________________________________________
6d81c2de 1719void AliMultiplicityCorrection::GetComparisonResults(Float_t* mc, Int_t* mcLimit, Float_t* residuals, Float_t* ratioAverage) const
cfc19dd5 1720{
1721 // Returns the chi2 between the MC and the unfolded ESD as well as between the ESD and the folded unfolded ESD
1722 // These values are computed during DrawComparison, thus this function picks up the
1723 // last calculation
1724
dd701109 1725 if (mc)
1726 *mc = fLastChi2MC;
1727 if (mcLimit)
1728 *mcLimit = fLastChi2MCLimit;
1729 if (residuals)
1730 *residuals = fLastChi2Residuals;
0b4bfd98 1731 if (ratioAverage)
1732 *ratioAverage = fRatioAverage;
cfc19dd5 1733}
1734
cfc19dd5 1735//____________________________________________________________________
1736TH2F* AliMultiplicityCorrection::GetMultiplicityMC(Int_t i, EventType eventType)
1737{
1738 //
1739 // returns the corresponding MC spectrum
1740 //
1741
1742 switch (eventType)
1743 {
1744 case kTrVtx : return fMultiplicityVtx[i]; break;
1745 case kMB : return fMultiplicityMB[i]; break;
1746 case kINEL : return fMultiplicityINEL[i]; break;
1747 }
1748
1749 return 0;
1750}
1751
1752//____________________________________________________________________
6d81c2de 1753TH1* AliMultiplicityCorrection::StatisticalUncertainty(MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Bool_t randomizeMeasured, Bool_t randomizeResponse, TH1* compareTo)
0a173978 1754{
1755 //
0b4bfd98 1756 // evaluates the uncertainty that arises from the non-infinite statistics in the response matrix
1757 // the function unfolds the spectrum using the default response matrix and several modified ones
1758 // the modified ones are created by randomizing each cell using poisson statistics with the mean = bin value
1759 // these unfolded results are compared to the first result gained with the default response OR to the histogram given
1760 // in <compareTo> (optional)
0a173978 1761 //
0b4bfd98 1762 // returns the error assigned to the measurement
1763 //
1764
6d81c2de 1765 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1766
0b4bfd98 1767 // initialize seed with current time
1768 gRandom->SetSeed(0);
1769
1770 const Int_t kErrorIterations = 150;
1771
1772 TH1* maxError = 0;
1773 TH1* firstResult = 0;
1774
1775 TH1** results = new TH1*[kErrorIterations];
1776
1777 for (Int_t n=0; n<kErrorIterations; ++n)
1778 {
1779 Printf("Iteration %d of %d...", n, kErrorIterations);
1780
6d81c2de 1781 // only done for chi2 minimization
1782 Bool_t createBigBin = kFALSE;
1783 if (methodType == kChi2Minimization)
1784 createBigBin = kTRUE;
1785
1786 SetupCurrentHists(inputRange, fullPhaseSpace, eventType, createBigBin);
0b4bfd98 1787
1788 TH1* measured = (TH1*) fCurrentESD->Clone("measured");
1789
1790 if (n > 0)
1791 {
1792 if (randomizeResponse)
1793 {
1794 // randomize response matrix
1795 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1796 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1797 fCurrentCorrelation->SetBinContent(i, j, gRandom->Poisson(fCurrentCorrelation->GetBinContent(i, j)));
1798 }
1799
1800 if (randomizeMeasured)
1801 {
1802 // randomize measured spectrum
1803 for (Int_t x=1; x<=measured->GetNbinsX(); x++) // mult. axis
1804 {
1805 Int_t randomValue = gRandom->Poisson(fCurrentESD->GetBinContent(x));
1806 measured->SetBinContent(x, randomValue);
1807 measured->SetBinError(x, TMath::Sqrt(randomValue));
1808 }
1809 }
1810 }
1811
6d81c2de 1812 // only for bayesian method we have to do it before the call to Unfold...
1813 if (methodType == kBayesian)
0b4bfd98 1814 {
6d81c2de 1815 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
0b4bfd98 1816 {
6d81c2de 1817 // with this it is normalized to 1
1818 Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1819
1820 // with this normalized to the given efficiency
1821 if (fCurrentEfficiency->GetBinContent(i) > 0)
1822 sum /= fCurrentEfficiency->GetBinContent(i);
0b4bfd98 1823 else
6d81c2de 1824 sum = 0;
1825
1826 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
0b4bfd98 1827 {
6d81c2de 1828 if (sum > 0)
1829 {
1830 fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1831 fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1832 }
1833 else
1834 {
1835 fCurrentCorrelation->SetBinContent(i, j, 0);
1836 fCurrentCorrelation->SetBinError(i, j, 0);
1837 }
0b4bfd98 1838 }
1839 }
1840 }
1841
1842 TH1* result = 0;
1843 if (n == 0 && compareTo)
1844 {
1845 // in this case we just store the histogram we want to compare to
1846 result = (TH1*) compareTo->Clone("compareTo");
1847 result->Sumw2();
1848 }
1849 else
6d81c2de 1850 {
1851 result = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone(Form("result_%d", n));
0b4bfd98 1852
6d81c2de 1853 Int_t returnCode = -1;
1854 if (methodType == kChi2Minimization)
1855 {
1856 returnCode = UnfoldWithMinuit(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result, kFALSE);
1857 }
1858 else if (methodType == kBayesian)
1859 returnCode = UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result, fgBayesianSmoothing, fgBayesianIterations);
1860
1861 if (returnCode != 0)
1862 return 0;
1863 }
0b4bfd98 1864
1865 // normalize
1866 result->Scale(1.0 / result->Integral());
1867
1868 if (n == 0)
1869 {
1870 firstResult = (TH1*) result->Clone("firstResult");
1871
1872 maxError = (TH1*) result->Clone("maxError");
1873 maxError->Reset();
1874 }
1875 else
1876 {
1877 // calculate ratio
1878 TH1* ratio = (TH1*) firstResult->Clone("ratio");
1879 ratio->Divide(result);
1880
1881 // find max. deviation
1882 for (Int_t x=1; x<=ratio->GetNbinsX(); x++)
1883 maxError->SetBinContent(x, TMath::Max(maxError->GetBinContent(x), TMath::Abs(1 - ratio->GetBinContent(x))));
1884
1885 delete ratio;
1886 }
1887
1888 results[n] = result;
1889 }
1890
1891 // find covariance matrix
1892 // results[n] is X_x
1893 // cov. matrix is M_xy = E ( (X_x - E(X_x)) * (X_y - E(X_y))), with E() = expectation value
1894
1895 Int_t nBins = results[0]->GetNbinsX();
1896 Float_t lowEdge = results[0]->GetXaxis()->GetBinLowEdge(1);
1897 Float_t upEdge = results[0]->GetXaxis()->GetBinUpEdge(nBins);
1898
1899 // find average, E(X)
1900 TProfile* average = new TProfile("average", "average", nBins, lowEdge, upEdge);
1901 for (Int_t n=1; n<kErrorIterations; ++n)
1902 for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1903 average->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x));
1904 //new TCanvas; average->DrawClone();
1905
1906 // find cov. matrix
1907 TProfile2D* covMatrix = new TProfile2D("covMatrix", "covMatrix", nBins, lowEdge, upEdge, nBins, lowEdge, upEdge);
1908
1909 for (Int_t n=1; n<kErrorIterations; ++n)
1910 for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1911 for (Int_t y=1; y<=results[n]->GetNbinsX(); y++)
1912 {
1913 // (X_x - E(X_x)) * (X_y - E(X_y)
1914 Float_t cov = (results[n]->GetBinContent(x) - average->GetBinContent(x)) * (results[n]->GetBinContent(y) - average->GetBinContent(y));
1915 covMatrix->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetXaxis()->GetBinCenter(y), cov);
1916 }
1917 new TCanvas; covMatrix->DrawClone("COLZ");
1918
6d81c2de 1919// // fill 2D histogram that contains deviation from first
1920// TH2F* deviations = new TH2F("deviations", "deviations", nBins, lowEdge, upEdge, 1000, -0.01, 0.01);
1921// for (Int_t n=1; n<kErrorIterations; ++n)
1922// for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1923// deviations->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x) - results[0]->GetBinContent(x));
1924// //new TCanvas; deviations->DrawCopy("COLZ");
1925//
1926// // get standard deviation "by hand"
1927// for (Int_t x=1; x<=nBins; x++)
1928// {
1929// if (results[0]->GetBinContent(x) > 0)
1930// {
1931// TH1* proj = deviations->ProjectionY("projRMS", x, x, "e");
1932// Float_t standardDev = proj->GetRMS(); // this is standard deviation in fact
1933// //standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x));
1934// Printf("sigma_%d is %f value %f --> error %f", x, standardDev, results[0]->GetBinContent(x), standardDev / results[0]->GetBinContent(x));
1935// }
1936// }
1937
1938 // calculate standard deviation of (results[0] - results[n])
0b4bfd98 1939 TH1* standardDeviation = (TH1*) maxError->Clone("standardDeviation");
1940 standardDeviation->Reset();
0a173978 1941
6d81c2de 1942 for (Int_t x=1; x<=results[0]->GetNbinsX(); x++)
0b4bfd98 1943 {
1944 if (results[0]->GetBinContent(x) > 0)
1945 {
6d81c2de 1946 Double_t average = 0;
1947 for (Int_t n=1; n<kErrorIterations; ++n)
1948 average += results[n]->GetBinContent(x) - results[0]->GetBinContent(x);
1949 average /= kErrorIterations-1;
1950
1951 Double_t variance = 0;
1952 for (Int_t n=1; n<kErrorIterations; ++n)
1953 {
1954 Double_t value = results[n]->GetBinContent(x) - results[0]->GetBinContent(x) - average;
1955 variance += value * value;
1956 }
1957 variance /= kErrorIterations-1;
1958
1959 Double_t standardDev = TMath::Sqrt(variance);
0b4bfd98 1960 standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x));
6d81c2de 1961 Printf("sigma_%d is %f value %f --> error %f", x, standardDev, results[0]->GetBinContent(x), standardDev / results[0]->GetBinContent(x));
0b4bfd98 1962 }
1963 }
1964
1965 // compare maxError to RMS of profile (variable name: average)
1966 // first: calculate rms in percent of value
1967 TH1* rmsError = (TH1*) maxError->Clone("rmsError");
1968 rmsError->Reset();
1969
1970 // enable error to be standard deviation (see http://root.cern.ch/root/html/TProfile.html#TProfile:SetErrorOption)
1971 average->SetErrorOption("s");
1972 for (Int_t x=1; x<=rmsError->GetNbinsX(); x++)
1973 if (average->GetBinContent(x) > 0)
1974 rmsError->SetBinContent(x, average->GetBinError(x) / average->GetBinContent(x));
1975
1976 // find maxError deviation from average (not from "first result"/mc as above)
1977 TH1* maxError2 = (TH1*) maxError->Clone("maxError2");
1978 maxError2->Reset();
1979 for (Int_t n=1; n<kErrorIterations; ++n)
1980 for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1981 if (average->GetBinContent(x) > 0)
1982 maxError2->SetBinContent(x, TMath::Max(maxError2->GetBinContent(x), TMath::Abs((results[n]->GetBinContent(x) - average->GetBinContent(x)) / average->GetBinContent(x))));
1983
6d81c2de 1984 //new TCanvas; maxError2->DrawCopy(); rmsError->SetLineColor(2); rmsError->DrawCopy("SAME"); standardDeviation->SetLineColor(3); standardDeviation->DrawCopy("SAME");
0b4bfd98 1985
1986 // plot difference between average and MC/first result
1987 TH1* averageFirstRatio = (TH1*) results[0]->Clone("averageFirstRatio");
1988 averageFirstRatio->Reset();
1989 averageFirstRatio->Divide(results[0], average);
1990
6d81c2de 1991 //new TCanvas; results[0]->DrawCopy(); average->SetLineColor(2); average->DrawClone("SAME");
1992 //new TCanvas; averageFirstRatio->DrawCopy();
0b4bfd98 1993
1994 // clean up
1995 for (Int_t n=0; n<kErrorIterations; ++n)
1996 delete results[n];
1997 delete[] results;
1998
1999 // fill into result histogram
0b4bfd98 2000 for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
2001 fMultiplicityESDCorrected[correlationID]->SetBinContent(i, firstResult->GetBinContent(i));
cfc19dd5 2002
0b4bfd98 2003 for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
2004 fMultiplicityESDCorrected[correlationID]->SetBinError(i, maxError->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
447c325d 2005
0b4bfd98 2006 return standardDeviation;
2007}
2008
2009//____________________________________________________________________
6d81c2de 2010void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar, Int_t nIterations, TH1* initialConditions, Bool_t determineError)
0b4bfd98 2011{
2012 //
2013 // correct spectrum using bayesian method
2014 //
2015
2016 // initialize seed with current time
2017 gRandom->SetSeed(0);
2018
2019 SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
0a173978 2020
2021 // normalize correction for given nPart
9ca6be09 2022 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
0a173978 2023 {
dd701109 2024 // with this it is normalized to 1
2025 Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
2026
2027 // with this normalized to the given efficiency
2028 if (fCurrentEfficiency->GetBinContent(i) > 0)
2029 sum /= fCurrentEfficiency->GetBinContent(i);
2030 else
2031 sum = 0;
2032
9ca6be09 2033 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
0a173978 2034 {
dd701109 2035 if (sum > 0)
2036 {
2037 fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
2038 fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
2039 }
2040 else
2041 {
2042 fCurrentCorrelation->SetBinContent(i, j, 0);
2043 fCurrentCorrelation->SetBinError(i, j, 0);
2044 }
0a173978 2045 }
2046 }
2047
0b4bfd98 2048 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
2049
6d81c2de 2050 if (UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], regPar, nIterations) != 0)
0b4bfd98 2051 return;
2052
0b4bfd98 2053 if (!determineError)
447c325d 2054 {
0b4bfd98 2055 Printf("AliMultiplicityCorrection::ApplyBayesianMethod: WARNING: No errors calculated.");
2056 return;
2057 }
447c325d 2058
0b4bfd98 2059 // evaluate errors, this is done by randomizing the measured spectrum following Poission statistics
2060 // this (new) measured spectrum is then unfolded and the different to the result from the "real" measured
2061 // spectrum calculated. This is performed N times and the maximum difference is taken as the systematic
2062 // error of the unfolding method itself.
2063
6d81c2de 2064 TH1* maxError = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone("maxError");
0b4bfd98 2065 maxError->Reset();
2066
6d81c2de 2067 TH1* normalizedResult = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone("normalizedResult");
0b4bfd98 2068 normalizedResult->Scale(1.0 / normalizedResult->Integral());
2069
2070 const Int_t kErrorIterations = 20;
2071
2072 printf("Spectrum unfolded. Determining error (%d iterations)...\n", kErrorIterations);
2073
2074 TH1* randomized = (TH1*) fCurrentESD->Clone("randomized");
6d81c2de 2075 TH1* result2 = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone("result2");
0b4bfd98 2076 for (Int_t n=0; n<kErrorIterations; ++n)
2077 {
2078 // randomize the content of clone following a poisson with the mean = the value of that bin
2079 for (Int_t x=1; x<=randomized->GetNbinsX(); x++) // mult. axis
447c325d 2080 {
0b4bfd98 2081 Int_t randomValue = gRandom->Poisson(fCurrentESD->GetBinContent(x));
2082 //printf("%d --> %d\n", fCurrentESD->GetBinContent(x), randomValue);
2083 randomized->SetBinContent(x, randomValue);
2084 randomized->SetBinError(x, TMath::Sqrt(randomValue));
447c325d 2085 }
447c325d 2086
6d81c2de 2087 result2->Reset();
2088 if (UnfoldWithBayesian(fCurrentCorrelation, fCurrentEfficiency, randomized, initialConditions, result2, regPar, nIterations) != 0)
0b4bfd98 2089 return;
cfc19dd5 2090
0b4bfd98 2091 result2->Scale(1.0 / result2->Integral());
cfc19dd5 2092
0b4bfd98 2093 // calculate ratio
6d81c2de 2094 result2->Divide(normalizedResult);
dd701109 2095
0b4bfd98 2096 //new TCanvas; ratio->DrawCopy("HIST");
6127aca6 2097
0b4bfd98 2098 // find max. deviation
6d81c2de 2099 for (Int_t x=1; x<=result2->GetNbinsX(); x++)
2100 maxError->SetBinContent(x, TMath::Max(maxError->GetBinContent(x), TMath::Abs(1 - result2->GetBinContent(x))));
0b4bfd98 2101 }
6d81c2de 2102 delete randomized;
2103 delete result2;
447c325d 2104
0b4bfd98 2105 for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
2106 fMultiplicityESDCorrected[correlationID]->SetBinError(i, maxError->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
447c325d 2107
0b4bfd98 2108 delete maxError;
2109 delete normalizedResult;
2110}
447c325d 2111
0b4bfd98 2112//____________________________________________________________________
6d81c2de 2113Int_t AliMultiplicityCorrection::UnfoldWithBayesian(TH1* correlation, TH1* aEfficiency, TH1* measured, TH1* initialConditions, TH1* aResult, Float_t regPar, Int_t nIterations)
0b4bfd98 2114{
2115 //
2116 // unfolds a spectrum
2117 //
2118
2119 if (measured->Integral() <= 0)
dd701109 2120 {
0b4bfd98 2121 Printf("AliMultiplicityCorrection::UnfoldWithBayesian: ERROR: The measured spectrum is empty");
6d81c2de 2122 return -1;
dd701109 2123 }
6127aca6 2124
0b4bfd98 2125 const Int_t kStartBin = 0;
6127aca6 2126
6d81c2de 2127 const Int_t kMaxM = fgkMaxInput; //<= fCurrentCorrelation->GetNbinsY(); // max measured axis
2128 const Int_t kMaxT = fgkMaxParams; //<= fCurrentCorrelation->GetNbinsX(); // max true axis
6127aca6 2129
0b4bfd98 2130 // store information in arrays, to increase processing speed (~ factor 5)
2131 Double_t measuredCopy[kMaxM];
2132 Double_t prior[kMaxT];
2133 Double_t result[kMaxT];
2134 Double_t efficiency[kMaxT];
b477f4f2 2135
0b4bfd98 2136 Double_t response[kMaxT][kMaxM];
2137 Double_t inverseResponse[kMaxT][kMaxM];
dd701109 2138
0b4bfd98 2139 // for normalization
2140 Float_t measuredIntegral = measured->Integral();
2141 for (Int_t m=0; m<kMaxM; m++)
2142 {
2143 measuredCopy[m] = measured->GetBinContent(m+1) / measuredIntegral;
2144
2145 for (Int_t t=0; t<kMaxT; t++)
2146 {
6d81c2de 2147 response[t][m] = correlation->GetBinContent(t+1, m+1);
0b4bfd98 2148 inverseResponse[t][m] = 0;
2149 }
2150 }
2151
2152 for (Int_t t=0; t<kMaxT; t++)
2153 {
6d81c2de 2154 efficiency[t] = aEfficiency->GetBinContent(t+1);
0b4bfd98 2155 prior[t] = measuredCopy[t];
2156 result[t] = 0;
2157 }
2158
2159 // pick prior distribution
6d81c2de 2160 if (initialConditions)
0b4bfd98 2161 {
2162 printf("Using different starting conditions...\n");
2163 // for normalization
6d81c2de 2164 Float_t inputDistIntegral = initialConditions->Integral();
0b4bfd98 2165 for (Int_t i=0; i<kMaxT; i++)
6d81c2de 2166 prior[i] = initialConditions->GetBinContent(i+1) / inputDistIntegral;
0b4bfd98 2167 }
2168
2169 //TH1F* convergence = new TH1F("convergence", "convergence", 50, 0.5, 50.5);
447c325d 2170
6127aca6 2171 // unfold...
dd701109 2172 for (Int_t i=0; i<nIterations; i++)
cfc19dd5 2173 {
2174 //printf(" iteration %i \n", i);
2175
6127aca6 2176 // calculate IR from Beyes theorem
2177 // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
9ca6be09 2178
0b4bfd98 2179 for (Int_t m=0; m<kMaxM; m++)
cfc19dd5 2180 {
2181 Float_t norm = 0;
0b4bfd98 2182 for (Int_t t = kStartBin; t<kMaxT; t++)
2183 norm += response[t][m] * prior[t];
2184
2185 if (norm > 0)
cfc19dd5 2186 {
0b4bfd98 2187 for (Int_t t = kStartBin; t<kMaxT; t++)
2188 inverseResponse[t][m] = response[t][m] * prior[t] / norm;
2189 }
2190 else
2191 {
2192 for (Int_t t = kStartBin; t<kMaxT; t++)
2193 inverseResponse[t][m] = 0;
6127aca6 2194 }
2195 }
cfc19dd5 2196
0b4bfd98 2197 for (Int_t t = kStartBin; t<kMaxT; t++)
cfc19dd5 2198 {
2199 Float_t value = 0;
0b4bfd98 2200 for (Int_t m=0; m<kMaxM; m++)
2201 value += inverseResponse[t][m] * measuredCopy[m];
cfc19dd5 2202
0b4bfd98 2203 if (efficiency[t] > 0)
2204 result[t] = value / efficiency[t];
dd701109 2205 else
0b4bfd98 2206 result[t] = 0;
cfc19dd5 2207 }
2208
b477f4f2 2209 // regularization (simple smoothing)
0b4bfd98 2210 for (Int_t t=kStartBin; t<kMaxT; t++)
cfc19dd5 2211 {
0b4bfd98 2212 Float_t newValue = 0;
2213 // 0 bin excluded from smoothing
2214 if (t > kStartBin+1 && t<kMaxT-1)
2215 {
2216 Float_t average = (result[t-1] + result[t] + result[t+1]) / 3;
9ca6be09 2217
0b4bfd98 2218 // weight the average with the regularization parameter
2219 newValue = (1 - regPar) * result[t] + regPar * average;
2220 }
2221 else
2222 newValue = result[t];
2223
2224 prior[t] = newValue;
6127aca6 2225 }
9ca6be09 2226
cfc19dd5 2227 // calculate chi2 (change from last iteration)
0b4bfd98 2228 //Double_t chi2 = 0;
2229 /*for (Int_t t=kStartBin; t<=hTrueSmoothed->GetNbinsX(); t++)
cfc19dd5 2230 {
0b4bfd98 2231 Float_t newValue = hTrueSmoothed->GetBinContent(t);
2232 //Float_t diff = hPrior->GetBinContent(t) - newValue;
2233 //chi2 += (Double_t) diff * diff;
cfc19dd5 2234
2235 hPrior->SetBinContent(t, newValue);
6127aca6 2236 }
cfc19dd5 2237 //printf("Chi2 of %d iteration = %.10f\n", i, chi2);
0b4bfd98 2238 //convergence->Fill(i+1, chi2);
dd701109 2239
cfc19dd5 2240 //if (i % 5 == 0)
2241 // DrawComparison(Form("Bayesian_%d", i), mcTarget, correlationID, kTRUE, eventType);
2242
0b4bfd98 2243 delete hTrueSmoothed;*/
6127aca6 2244 } // end of iterations
2245
dd701109 2246 //new TCanvas;
2247 //convergence->DrawCopy();
2248 //gPad->SetLogy();
2249
0b4bfd98 2250 //delete convergence;
2251
0b4bfd98 2252 for (Int_t t=0; t<kMaxT; t++)
6d81c2de 2253 aResult->SetBinContent(t+1, result[t]);
6127aca6 2254
6d81c2de 2255 return 0;
cfc19dd5 2256
2257 // ********
2258 // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995
2259
dd701109 2260 /*printf("Calculating covariance matrix. This may take some time...\n");
2261
2262 // TODO check if this is the right one...
2263 TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
cfc19dd5 2264
2265 Int_t xBins = hInverseResponseBayes->GetNbinsX();
2266 Int_t yBins = hInverseResponseBayes->GetNbinsY();
2267
2268 // calculate "unfolding matrix" Mij
2269 Float_t matrixM[251][251];
2270 for (Int_t i=1; i<=xBins; i++)
2271 {
2272 for (Int_t j=1; j<=yBins; j++)
2273 {
2274 if (fCurrentEfficiency->GetBinContent(i) > 0)
2275 matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i);
2276 else
2277 matrixM[i-1][j-1] = 0;
2278 }
2279 }
2280
2281 Float_t* vectorn = new Float_t[yBins];
2282 for (Int_t j=1; j<=yBins; j++)
2283 vectorn[j-1] = fCurrentESD->GetBinContent(j);
2284
2285 // first part of covariance matrix, depends on input distribution n(E)
2286 Float_t cov1[251][251];
2287
2288 Float_t nEvents = fCurrentESD->Integral(); // N
2289
2290 xBins = 20;
2291 yBins = 20;
2292
2293 for (Int_t k=0; k<xBins; k++)
2294 {
2295 printf("In Cov1: %d\n", k);
2296 for (Int_t l=0; l<yBins; l++)
2297 {
2298 cov1[k][l] = 0;
2299
2300 // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N)
2301 for (Int_t j=0; j<yBins; j++)
2302 cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j]
2303 * (1.0 - vectorn[j] / nEvents);
2304
2305 // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N
2306 for (Int_t i=0; i<yBins; i++)
2307 for (Int_t j=0; j<yBins; j++)
2308 {
2309 if (i == j)
2310 continue;
2311 cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i]
2312 * vectorn[j] / nEvents;
2313 }
2314 }
2315 }
2316
2317 printf("Cov1 finished\n");
2318
2319 TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov");
2320 cov->Reset();
2321
2322 for (Int_t i=1; i<=xBins; i++)
2323 for (Int_t j=1; j<=yBins; j++)
2324 cov->SetBinContent(i, j, cov1[i-1][j-1]);
2325
2326 new TCanvas;
2327 cov->Draw("COLZ");
2328
2329 // second part of covariance matrix, depends on response matrix
2330 Float_t cov2[251][251];
2331
2332 // Cov[P(Er|Cu), P(Es|Cu)] term
2333 Float_t covTerm[100][100][100];
2334 for (Int_t r=0; r<yBins; r++)
2335 for (Int_t u=0; u<xBins; u++)
2336 for (Int_t s=0; s<yBins; s++)
2337 {
2338 if (r == s)
2339 covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
2340 * (1.0 - hResponse->GetBinContent(u+1, r+1));
2341 else
2342 covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
2343 * hResponse->GetBinContent(u+1, s+1);
2344 }
2345
2346 for (Int_t k=0; k<xBins; k++)
2347 for (Int_t l=0; l<yBins; l++)
2348 {
2349 cov2[k][l] = 0;
2350 printf("In Cov2: %d %d\n", k, l);
2351 for (Int_t i=0; i<yBins; i++)
2352 for (Int_t j=0; j<yBins; j++)
2353 {
2354 //printf("In Cov2: %d %d %d %d\n", k, l, i, j);
2355 // calculate Cov(Mki, Mlj) = sum{ru},{su} ...
2356 Float_t tmpCov = 0;
2357 for (Int_t r=0; r<yBins; r++)
2358 for (Int_t u=0; u<xBins; u++)
2359 for (Int_t s=0; s<yBins; s++)
2360 {
2361 if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0
2362 || hResponse->GetBinContent(u+1, i+1) == 0)
2363 continue;
2364
2365 tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u)
2366 * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u)
2367 * covTerm[r][u][s];
2368 }
2369
2370 cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov;
2371 }
2372 }
2373
2374 printf("Cov2 finished\n");
2375
2376 for (Int_t i=1; i<=xBins; i++)
2377 for (Int_t j=1; j<=yBins; j++)
2378 cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]);
2379
2380 new TCanvas;
dd701109 2381 cov->Draw("COLZ");*/
2382}
2383
2384//____________________________________________________________________
0b4bfd98 2385Float_t AliMultiplicityCorrection::BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse, Int_t k, Int_t i, Int_t r, Int_t u)
dd701109 2386{
2387 //
2388 // helper function for the covariance matrix of the bayesian method
2389 //
2390
2391 Float_t result = 0;
2392
2393 if (k == u && r == i)
2394 result += 1.0 / hResponse->GetBinContent(u+1, r+1);
2395
2396 if (k == u)
2397 result -= 1.0 / fCurrentEfficiency->GetBinContent(u+1);
2398
2399 if (r == i)
2400 result -= matrixM[u][i] * fCurrentEfficiency->GetBinContent(u+1) / hResponse->GetBinContent(u+1, i+1);
2401
2402 result *= matrixM[k][i];
2403
2404 return result;
cfc19dd5 2405}
2406
2407//____________________________________________________________________
2408void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
2409{
2410 //
2411 // correct spectrum using bayesian method
2412 //
2413
dd701109 2414 Float_t regPar = 0;
cfc19dd5 2415
2416 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
2417 Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
2418
447c325d 2419 SetupCurrentHists(inputRange, fullPhaseSpace, eventType, kFALSE);
0b4bfd98 2420 //normalize ESD
2421 fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
cfc19dd5 2422
2423 // TODO should be taken from correlation map
2424 //TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
2425
2426 // normalize correction for given nPart
2427 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2428 {
2429 Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
2430 //Double_t sum = sumHist->GetBinContent(i);
2431 if (sum <= 0)
2432 continue;
2433 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
2434 {
2435 // npart sum to 1
2436 fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
2437 fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
2438 }
6127aca6 2439 }
cfc19dd5 2440
2441 new TCanvas;
2442 fCurrentCorrelation->Draw("COLZ");
2443
2444 // FAKE
2445 fCurrentEfficiency = ((TH2*) fCurrentCorrelation)->ProjectionX("eff");
2446
2447 // pick prior distribution
2448 TH1F* hPrior = (TH1F*)fCurrentESD->Clone("prior");
2449 Float_t norm = 1; //hPrior->Integral();
2450 for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
2451 hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm);
2452
2453 // zero distribution
2454 TH1F* zero = (TH1F*)hPrior->Clone("zero");
2455
2456 // define temp hist
2457 TH1F* hTemp = (TH1F*)fCurrentESD->Clone("temp");
2458 hTemp->Reset();
2459
2460 // just a shortcut
2461 TH2F* hResponse = (TH2F*) fCurrentCorrelation;
2462
2463 // unfold...
dd701109 2464 Int_t iterations = 25;
cfc19dd5 2465 for (Int_t i=0; i<iterations; i++)
2466 {
2467 //printf(" iteration %i \n", i);
2468
2469 for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
2470 {
2471 Float_t value = 0;
2472 for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
2473 value += hResponse->GetBinContent(t, m) * hPrior->GetBinContent(t);
2474 hTemp->SetBinContent(m, value);
2475 //printf("%d %f %f %f\n", m, zero->GetBinContent(m), hPrior->GetBinContent(m), value);
2476 }
2477
2478 // regularization (simple smoothing)
2479 TH1F* hTrueSmoothed = (TH1F*) hTemp->Clone("truesmoothed");
2480
2481 for (Int_t t=2; t<hTrueSmoothed->GetNbinsX(); t++)
2482 {
2483 Float_t average = (hTemp->GetBinContent(t-1) / hTemp->GetBinWidth(t-1)
2484 + hTemp->GetBinContent(t) / hTemp->GetBinWidth(t)
2485 + hTemp->GetBinContent(t+1) / hTemp->GetBinWidth(t+1)) / 3.;
2486 average *= hTrueSmoothed->GetBinWidth(t);
2487
2488 // weight the average with the regularization parameter
2489 hTrueSmoothed->SetBinContent(t, (1-regPar) * hTemp->GetBinContent(t) + regPar * average);
2490 }
2491
2492 for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
2493 hTemp->SetBinContent(m, zero->GetBinContent(m) + hPrior->GetBinContent(m) - hTrueSmoothed->GetBinContent(m));
2494
2495 // fill guess
2496 for (Int_t t=1; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++)
2497 {
2498 fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTemp->GetBinContent(t));
2499 fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t)); // TODO
2500
2501 //printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
2502 }
2503
2504
2505 // calculate chi2 (change from last iteration)
2506 Double_t chi2 = 0;
2507
2508 // use smoothed true (normalized) as new prior
2509 Float_t norm = 1; //hTrueSmoothed->Integral();
2510
2511 for (Int_t t=1; t<hTrueSmoothed->GetNbinsX(); t++)
2512 {
dd701109 2513 Float_t newValue = hTemp->GetBinContent(t)/norm;
cfc19dd5 2514 Float_t diff = hPrior->GetBinContent(t) - newValue;
2515 chi2 += (Double_t) diff * diff;
2516
2517 hPrior->SetBinContent(t, newValue);
2518 }
2519
2520 printf("Chi2 of %d iteration = %.10f\n", i, chi2);
2521
dd701109 2522 //if (i % 5 == 0)
cfc19dd5 2523 DrawComparison(Form("Laszlo_%d", i), inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
2524
2525 delete hTrueSmoothed;
2526 } // end of iterations
2527
2528 DrawComparison("Laszlo", inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
2529}
2530
9ca6be09 2531//____________________________________________________________________
2532void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace)
2533{
2534 //
2535 // correct spectrum using a simple Gaussian approach, that is model-dependent
2536 //
2537
2538 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
2539 Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
2540
447c325d 2541 SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx, kFALSE);
0b4bfd98 2542 //normalize ESD
2543 fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
9ca6be09 2544
9ca6be09 2545 TH1D* correction = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianMean"));
2546 correction->SetTitle("GaussianMean");
2547
2548 TH1D* correctionWidth = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianWidth"));
2549 correction->SetTitle("GaussianWidth");
2550
2551 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2552 {
2553 TH1D* proj = (dynamic_cast<TH2*> (fCurrentCorrelation))->ProjectionX("_px", i, i+1);
2554 proj->Fit("gaus", "0Q");
2555 correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
2556 correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
2557
2558 continue;
2559
2560 // draw for debugging
2561 new TCanvas;
2562 proj->DrawCopy();
2563 proj->GetFunction("gaus")->DrawCopy("SAME");
2564 }
2565
2566 TH1* target = fMultiplicityESDCorrected[correlationID];
2567
2568 Int_t nBins = target->GetNbinsX()*10+1;
2569 Float_t* binning = new Float_t[nBins];
2570 for (Int_t i=1; i<=target->GetNbinsX(); ++i)
2571 for (Int_t j=0; j<10; ++j)
2572 binning[(i-1)*10 + j] = target->GetXaxis()->GetBinLowEdge(i) + target->GetXaxis()->GetBinWidth(i) / 10 * j;
2573
2574 binning[nBins-1] = target->GetXaxis()->GetBinUpEdge(target->GetNbinsX());
2575
2576 TH1F* fineBinned = new TH1F("targetFineBinned", "targetFineBinned", nBins-1, binning);
2577
2578 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2579 {
2580 Float_t mean = correction->GetBinContent(i);
2581 Float_t width = correctionWidth->GetBinContent(i);
2582
3602328d 2583 Int_t fillBegin = fineBinned->FindBin(mean - width * 5);
2584 Int_t fillEnd = fineBinned->FindBin(mean + width * 5);
2585 //printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
9ca6be09 2586
2587 for (Int_t j=fillBegin; j <= fillEnd; ++j)
2588 {
2589 fineBinned->AddBinContent(j, TMath::Gaus(fineBinned->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fCurrentESD->GetBinContent(i));
2590 }
2591 }
2592
2593 for (Int_t i=1; i<=target->GetNbinsX(); ++i)
2594 {
2595 Float_t sum = 0;
2596 for (Int_t j=1; j<=10; ++j)
2597 sum += fineBinned->GetBinContent((i-1)*10 + j);
2598 target->SetBinContent(i, sum / 10);
2599 }
2600
2601 delete[] binning;
2602
cfc19dd5 2603 DrawComparison("Gaussian", inputRange, fullPhaseSpace, kFALSE, GetMultiplicityMC(mcTarget, kTrVtx)->ProjectionY());
0a173978 2604}
3602328d 2605
2606//____________________________________________________________________
447c325d 2607TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap)
3602328d 2608{
2609 // runs the distribution given in inputMC through the response matrix identified by
2610 // correlationMap and produces a measured distribution
2611 // although it is a TH2F the vertex axis is not used at the moment and all entries are filled in mid-vertex
b477f4f2 2612 // if normalized is set, inputMC is expected to be normalized to the bin width
3602328d 2613
2614 if (!inputMC)
2615 return 0;
2616
2617 if (correlationMap < 0 || correlationMap >= kCorrHists)
2618 return 0;
2619
2620 // empty under/overflow bins in x, otherwise Project3D takes them into account
2621 TH3* hist = fCorrelation[correlationMap];
2622 for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
2623 {
2624 for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
2625 {
2626 hist->SetBinContent(0, y, z, 0);
2627 hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
2628 }
2629 }
2630
447c325d 2631 TH2* corr = (TH2*) hist->Project3D("zy");
2632 //corr->Rebin2D(2, 1);
3602328d 2633 corr->Sumw2();
2634
2635 // normalize correction for given nPart
2636 for (Int_t i=1; i<=corr->GetNbinsX(); ++i)
2637 {
2638 Double_t sum = corr->Integral(i, i, 1, corr->GetNbinsY());
2639 if (sum <= 0)
2640 continue;
2641
2642 for (Int_t j=1; j<=corr->GetNbinsY(); ++j)
2643 {
2644 // npart sum to 1
2645 corr->SetBinContent(i, j, corr->GetBinContent(i, j) / sum);
2646 corr->SetBinError(i, j, corr->GetBinError(i, j) / sum);
2647 }
2648 }
2649
2650 TH2F* target = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("%s_measured", inputMC->GetName())));
2651 target->Reset();
2652
2653 for (Int_t meas=1; meas<=corr->GetNbinsY(); ++meas)
2654 {
2655 Float_t measured = 0;
2656 Float_t error = 0;
2657
447c325d 2658 for (Int_t gen=1; gen<=corr->GetNbinsX(); ++gen)
3602328d 2659 {
447c325d 2660 Int_t mcGenBin = inputMC->GetXaxis()->FindBin(corr->GetXaxis()->GetBinCenter(gen));
b477f4f2 2661
447c325d 2662 measured += inputMC->GetBinContent(mcGenBin) * corr->GetBinContent(gen, meas);
2663 error += inputMC->GetBinError(mcGenBin) * corr->GetBinContent(gen, meas);
3602328d 2664 }
2665
cfc19dd5 2666 //printf("%f +- %f ; %f +- %f \n", inputMC->GetBinContent(meas), inputMC->GetBinError(meas), measured, error);
3602328d 2667
447c325d 2668 target->SetBinContent(1 + target->GetNbinsX() / 2, meas, measured);
2669 target->SetBinError(1 + target->GetNbinsX() / 2, meas, error);
3602328d 2670 }
2671
2672 return target;
2673}
2674
2675//____________________________________________________________________
2676void AliMultiplicityCorrection::SetGenMeasFromFunc(TF1* inputMC, Int_t id)
2677{
2678 // uses the given function to fill the input MC histogram and generates from that
2679 // the measured histogram by applying the response matrix
2680 // this can be used to evaluate if the methods work indepedently of the input
2681 // distribution
2682 // WARNING does not respect the vertex distribution, just fills central vertex bin
2683
2684 if (!inputMC)
2685 return;
2686
2687 if (id < 0 || id >= kESDHists)
2688 return;
2689
dd701109 2690 TH2F* mc = fMultiplicityVtx[id];
3602328d 2691
2692 mc->Reset();
2693
2694 for (Int_t i=1; i<=mc->GetNbinsY(); ++i)
b477f4f2 2695 {
0b4bfd98 2696 mc->SetBinContent(mc->GetNbinsX() / 2 + 1, i, inputMC->Eval(mc->GetYaxis()->GetBinCenter(i)) * mc->GetYaxis()->GetBinWidth(i));
2697 mc->SetBinError(mc->GetNbinsX() / 2 + 1, i, 0);
b477f4f2 2698 }
3602328d 2699
2700 //new TCanvas;
2701 //mc->Draw("COLZ");
2702
2703 TH1* proj = mc->ProjectionY();
2704
dd701109 2705 TString tmp(fMultiplicityESD[id]->GetName());
2706 delete fMultiplicityESD[id];
3602328d 2707 fMultiplicityESD[id] = CalculateMultiplicityESD(proj, id);
dd701109 2708 fMultiplicityESD[id]->SetName(tmp);
3602328d 2709}