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