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