Added transformation function (Marian)
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AliMultiplicityCorrection.cxx
CommitLineData
0a173978 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
18// This class is used to store correction maps, raw input and results of the multiplicity
19// measurement with the ITS or TPC
20// It also contains functions to correct the spectrum using different methods.
21//
22// Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch
23
24#include "AliMultiplicityCorrection.h"
25
26#include <TFile.h>
27#include <TH1F.h>
28#include <TH2F.h>
29#include <TH3F.h>
30#include <TDirectory.h>
31#include <TVirtualFitter.h>
32#include <TCanvas.h>
33#include <TString.h>
9ca6be09 34#include <TF1.h>
d560b581 35#include <TMath.h>
3602328d 36#include <TCollection.h>
0a173978 37
38ClassImp(AliMultiplicityCorrection)
39
cfc19dd5 40const Int_t AliMultiplicityCorrection::fgMaxParams = 200;
9ca6be09 41TH1* AliMultiplicityCorrection::fCurrentESD = 0;
42TH1* AliMultiplicityCorrection::fCurrentCorrelation = 0;
cfc19dd5 43TH1* AliMultiplicityCorrection::fCurrentEfficiency = 0;
44TMatrixF* AliMultiplicityCorrection::fCorrelationMatrix = 0;
45TMatrixF* AliMultiplicityCorrection::fCorrelationCovarianceMatrix = 0;
46TVectorF* AliMultiplicityCorrection::fCurrentESDVector = 0;
47AliMultiplicityCorrection::RegularizationType AliMultiplicityCorrection::fRegularizationType = AliMultiplicityCorrection::kPol1;
48Float_t AliMultiplicityCorrection::fRegularizationWeight = 1e4;
0a173978 49
50//____________________________________________________________________
51AliMultiplicityCorrection::AliMultiplicityCorrection() :
cfc19dd5 52 TNamed(), fLastChi2MC(0), fLastChi2Residuals(0)
0a173978 53{
54 //
55 // default constructor
56 //
57
58 for (Int_t i = 0; i < kESDHists; ++i)
59 fMultiplicityESD[i] = 0;
60
61 for (Int_t i = 0; i < kMCHists; ++i)
cfc19dd5 62 {
63 fMultiplicityVtx[i] = 0;
64 fMultiplicityMB[i] = 0;
65 fMultiplicityINEL[i] = 0;
66 }
0a173978 67
68 for (Int_t i = 0; i < kCorrHists; ++i)
69 {
70 fCorrelation[i] = 0;
71 fMultiplicityESDCorrected[i] = 0;
72 }
73}
74
75//____________________________________________________________________
76AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const Char_t* title) :
cfc19dd5 77 TNamed(name, title), fLastChi2MC(0), fLastChi2Residuals(0)
0a173978 78{
79 //
80 // named constructor
81 //
82
83 // do not add this hists to the directory
84 Bool_t oldStatus = TH1::AddDirectoryStatus();
85 TH1::AddDirectory(kFALSE);
86
b477f4f2 87 /*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 88 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,
89 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
90 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
91 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
92 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5,
93 50.5, 55.5, 60.5, 65.5, 70.5, 75.5, 80.5, 85.5, 90.5, 95.5,
94 100.5, 105.5, 110.5, 115.5, 120.5, 125.5, 130.5, 135.5, 140.5, 145.5,
95 150.5, 160.5, 170.5, 180.5, 190.5, 200.5, 210.5, 220.5, 230.5, 240.5,
96 250.5, 275.5, 300.5, 325.5, 350.5, 375.5, 400.5, 425.5, 450.5, 475.5,
9ca6be09 97 500.5 };
98 //525.5, 550.5, 575.5, 600.5, 625.5, 650.5, 675.5, 700.5, 725.5,
99 //750.5, 775.5, 800.5, 825.5, 850.5, 875.5, 900.5, 925.5, 950.5, 975.5,
100 //1000.5 };
0a173978 101
b477f4f2 102 #define VTXBINNING 10, binLimitsVtx
103 #define NBINNING fgMaxParams, binLimitsN*/
104
105 #define NBINNING 251, -0.5, 250.5
106 #define VTXBINNING 10, -10, 10
0a173978 107
108 for (Int_t i = 0; i < kESDHists; ++i)
b477f4f2 109 fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", VTXBINNING, NBINNING);
0a173978 110
111 for (Int_t i = 0; i < kMCHists; ++i)
112 {
cfc19dd5 113 fMultiplicityVtx[i] = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("fMultiplicityVtx%d", i)));
114 fMultiplicityVtx[i]->SetTitle("fMultiplicityVtx;vtx-z;Npart");
115
116 fMultiplicityMB[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityMB%d", i)));
117 fMultiplicityMB[i]->SetTitle("fMultiplicityMB");
118
119 fMultiplicityINEL[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[0]->Clone(Form("fMultiplicityINEL%d", i)));
120 fMultiplicityINEL[i]->SetTitle("fMultiplicityINEL");
0a173978 121 }
122
123 for (Int_t i = 0; i < kCorrHists; ++i)
124 {
b477f4f2 125 fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;Npart;Ntracks", VTXBINNING, NBINNING, NBINNING);
126 fMultiplicityESDCorrected[i] = new TH1F(Form("fMultiplicityESDCorrected%d", i), "fMultiplicityESDCorrected;Npart;dN/dN", NBINNING);
0a173978 127 }
128
129 TH1::AddDirectory(oldStatus);
130}
131
132//____________________________________________________________________
133AliMultiplicityCorrection::~AliMultiplicityCorrection()
134{
135 //
136 // Destructor
137 //
138}
139
140//____________________________________________________________________
141Long64_t AliMultiplicityCorrection::Merge(TCollection* list)
142{
143 // Merge a list of AliMultiplicityCorrection objects with this (needed for
144 // PROOF).
145 // Returns the number of merged objects (including this).
146
147 if (!list)
148 return 0;
149
150 if (list->IsEmpty())
151 return 1;
152
153 TIterator* iter = list->MakeIterator();
154 TObject* obj;
155
156 // collections of all histograms
cfc19dd5 157 TList collections[kESDHists+kMCHists*3+kCorrHists*2];
0a173978 158
159 Int_t count = 0;
160 while ((obj = iter->Next())) {
161
162 AliMultiplicityCorrection* entry = dynamic_cast<AliMultiplicityCorrection*> (obj);
163 if (entry == 0)
164 continue;
165
166 for (Int_t i = 0; i < kESDHists; ++i)
167 collections[i].Add(entry->fMultiplicityESD[i]);
168
169 for (Int_t i = 0; i < kMCHists; ++i)
cfc19dd5 170 {
171 collections[kESDHists+i].Add(entry->fMultiplicityVtx[i]);
172 collections[kESDHists+kMCHists+i].Add(entry->fMultiplicityMB[i]);
173 collections[kESDHists+kMCHists*2+i].Add(entry->fMultiplicityINEL[i]);
174 }
0a173978 175
176 for (Int_t i = 0; i < kCorrHists; ++i)
cfc19dd5 177 collections[kESDHists+kMCHists*3+i].Add(entry->fCorrelation[i]);
0a173978 178
179 for (Int_t i = 0; i < kCorrHists; ++i)
cfc19dd5 180 collections[kESDHists+kMCHists*3+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]);
0a173978 181
182 count++;
183 }
184
185 for (Int_t i = 0; i < kESDHists; ++i)
186 fMultiplicityESD[i]->Merge(&collections[i]);
187
188 for (Int_t i = 0; i < kMCHists; ++i)
cfc19dd5 189 {
190 fMultiplicityVtx[i]->Merge(&collections[kESDHists+i]);
191 fMultiplicityMB[i]->Merge(&collections[kESDHists+kMCHists+i]);
192 fMultiplicityINEL[i]->Merge(&collections[kESDHists+kMCHists*2+i]);
193 }
0a173978 194
195 for (Int_t i = 0; i < kCorrHists; ++i)
cfc19dd5 196 fCorrelation[i]->Merge(&collections[kESDHists+kMCHists*3+i]);
0a173978 197
198 for (Int_t i = 0; i < kCorrHists; ++i)
cfc19dd5 199 fMultiplicityESDCorrected[i]->Merge(&collections[kESDHists+kMCHists*3+kCorrHists+i]);
0a173978 200
201 delete iter;
202
203 return count+1;
204}
205
206//____________________________________________________________________
207Bool_t AliMultiplicityCorrection::LoadHistograms(const Char_t* dir)
208{
209 //
210 // loads the histograms from a file
211 // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
212 //
213
214 if (!dir)
215 dir = GetName();
216
217 if (!gDirectory->cd(dir))
218 return kFALSE;
219
220 // TODO memory leak. old histograms needs to be deleted.
221
222 Bool_t success = kTRUE;
223
224 for (Int_t i = 0; i < kESDHists; ++i)
225 {
226 fMultiplicityESD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityESD[i]->GetName()));
227 if (!fMultiplicityESD[i])
228 success = kFALSE;
229 }
230
231 for (Int_t i = 0; i < kMCHists; ++i)
232 {
cfc19dd5 233 fMultiplicityVtx[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityVtx[i]->GetName()));
234 fMultiplicityMB[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityMB[i]->GetName()));
235 fMultiplicityINEL[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityINEL[i]->GetName()));
236 if (!fMultiplicityVtx[i] || !fMultiplicityMB[i] || !fMultiplicityINEL[i])
0a173978 237 success = kFALSE;
238 }
239
240 for (Int_t i = 0; i < kCorrHists; ++i)
241 {
242 fCorrelation[i] = dynamic_cast<TH3F*> (gDirectory->Get(fCorrelation[i]->GetName()));
243 if (!fCorrelation[i])
244 success = kFALSE;
245 fMultiplicityESDCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fMultiplicityESDCorrected[i]->GetName()));
246 if (!fMultiplicityESDCorrected[i])
247 success = kFALSE;
248 }
249
250 gDirectory->cd("..");
251
252 return success;
253}
254
255//____________________________________________________________________
256void AliMultiplicityCorrection::SaveHistograms()
257{
258 //
259 // saves the histograms
260 //
261
262 gDirectory->mkdir(GetName());
263 gDirectory->cd(GetName());
264
265 for (Int_t i = 0; i < kESDHists; ++i)
266 if (fMultiplicityESD[i])
267 fMultiplicityESD[i]->Write();
268
269 for (Int_t i = 0; i < kMCHists; ++i)
cfc19dd5 270 {
271 if (fMultiplicityVtx[i])
272 fMultiplicityVtx[i]->Write();
273 if (fMultiplicityMB[i])
274 fMultiplicityMB[i]->Write();
275 if (fMultiplicityINEL[i])
276 fMultiplicityINEL[i]->Write();
277 }
0a173978 278
279 for (Int_t i = 0; i < kCorrHists; ++i)
280 {
281 if (fCorrelation[i])
282 fCorrelation[i]->Write();
283 if (fMultiplicityESDCorrected[i])
284 fMultiplicityESDCorrected[i]->Write();
285 }
286
287 gDirectory->cd("..");
288}
289
290//____________________________________________________________________
cfc19dd5 291void 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 292{
293 //
294 // Fills an event from MC
295 //
296
cfc19dd5 297 if (triggered)
298 {
299 fMultiplicityMB[0]->Fill(vtx, generated05);
300 fMultiplicityMB[1]->Fill(vtx, generated10);
301 fMultiplicityMB[2]->Fill(vtx, generated15);
302 fMultiplicityMB[3]->Fill(vtx, generated20);
303 fMultiplicityMB[4]->Fill(vtx, generatedAll);
304
305 if (vertex)
306 {
307 fMultiplicityVtx[0]->Fill(vtx, generated05);
308 fMultiplicityVtx[1]->Fill(vtx, generated10);
309 fMultiplicityVtx[2]->Fill(vtx, generated15);
310 fMultiplicityVtx[3]->Fill(vtx, generated20);
311 fMultiplicityVtx[4]->Fill(vtx, generatedAll);
312 }
313 }
314
315 fMultiplicityINEL[0]->Fill(vtx, generated05);
316 fMultiplicityINEL[1]->Fill(vtx, generated10);
317 fMultiplicityINEL[2]->Fill(vtx, generated15);
318 fMultiplicityINEL[3]->Fill(vtx, generated20);
319 fMultiplicityINEL[4]->Fill(vtx, generatedAll);
0a173978 320}
321
322//____________________________________________________________________
323void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20)
324{
325 //
326 // Fills an event from ESD
327 //
328
329 fMultiplicityESD[0]->Fill(vtx, measured05);
330 fMultiplicityESD[1]->Fill(vtx, measured10);
331 fMultiplicityESD[2]->Fill(vtx, measured15);
332 fMultiplicityESD[3]->Fill(vtx, measured20);
333}
334
335//____________________________________________________________________
336void 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)
337{
338 //
339 // Fills an event into the correlation map with the information from MC and ESD
340 //
341
342 fCorrelation[0]->Fill(vtx, generated05, measured05);
343 fCorrelation[1]->Fill(vtx, generated10, measured10);
344 fCorrelation[2]->Fill(vtx, generated15, measured15);
345 fCorrelation[3]->Fill(vtx, generated20, measured20);
346
347 fCorrelation[4]->Fill(vtx, generatedAll, measured05);
348 fCorrelation[5]->Fill(vtx, generatedAll, measured10);
349 fCorrelation[6]->Fill(vtx, generatedAll, measured15);
350 fCorrelation[7]->Fill(vtx, generatedAll, measured20);
351}
352
353//____________________________________________________________________
9ca6be09 354Double_t AliMultiplicityCorrection::RegularizationPol0(Double_t *params)
0a173978 355{
356 // homogenity term for minuit fitting
357 // pure function of the parameters
358 // prefers constant function (pol0)
359
360 Double_t chi2 = 0;
361
362 for (Int_t i=1; i<fgMaxParams; ++i)
363 {
364 if (params[i] == 0)
365 continue;
366
cfc19dd5 367 Double_t right = params[i]; // / fCurrentESD->GetBinWidth(i+1);
368 Double_t left = params[i-1]; // / fCurrentESD->GetBinWidth(i);
0a173978 369
370 Double_t diff = (right - left) / right;
371 chi2 += diff * diff;
372 }
373
374 return chi2;
375}
376
377//____________________________________________________________________
9ca6be09 378Double_t AliMultiplicityCorrection::RegularizationPol1(Double_t *params)
0a173978 379{
380 // homogenity term for minuit fitting
381 // pure function of the parameters
382 // prefers linear function (pol1)
383
384 Double_t chi2 = 0;
385
386 for (Int_t i=2; i<fgMaxParams; ++i)
387 {
388 if (params[i] == 0 || params[i-1] == 0)
389 continue;
390
cfc19dd5 391 Double_t right = params[i]; // / fCurrentESD->GetBinWidth(i+1);
392 Double_t middle = params[i-1]; // / fCurrentESD->GetBinWidth(i);
393 Double_t left = params[i-2]; // / fCurrentESD->GetBinWidth(i-1);
0a173978 394
cfc19dd5 395 Double_t der1 = (right - middle); // / fCurrentESD->GetBinWidth(i);
396 Double_t der2 = (middle - left); // / fCurrentESD->GetBinWidth(i-1);
0a173978 397
cfc19dd5 398 Double_t diff = (der1 - der2) / middle; /// fCurrentESD->GetBinWidth(i);
0a173978 399
400 // TODO give additonal weight to big bins
cfc19dd5 401 chi2 += diff * diff; // * fCurrentESD->GetBinWidth(i) * fCurrentESD->GetBinWidth(i-1);
0a173978 402
403 //printf("%d --> %f\n", i, diff);
404 }
405
cfc19dd5 406 return chi2; // 10000
0a173978 407}
408
9ca6be09 409//____________________________________________________________________
410Double_t AliMultiplicityCorrection::RegularizationTotalCurvature(Double_t *params)
411{
412 // homogenity term for minuit fitting
413 // pure function of the parameters
414 // minimizes the total curvature (from Unfolding Methods In High-Energy Physics Experiments,
415 // V. Blobel (Hamburg U.) . DESY 84/118, Dec 1984. 40pp.
416
417 Double_t chi2 = 0;
418
419 for (Int_t i=2; i<fgMaxParams; ++i)
420 {
421 if (params[i] == 0 || params[i-1] == 0)
422 continue;
423
cfc19dd5 424 Double_t right = params[i]; // / fCurrentESD->GetBinWidth(i+1);
425 Double_t middle = params[i-1]; // / fCurrentESD->GetBinWidth(i);
426 Double_t left = params[i-2]; // / fCurrentESD->GetBinWidth(i-1);
9ca6be09 427
cfc19dd5 428 Double_t der1 = (right - middle); // / fCurrentESD->GetBinWidth(i);
429 Double_t der2 = (middle - left); // / fCurrentESD->GetBinWidth(i-1);
9ca6be09 430
cfc19dd5 431 Double_t secDer = (der1 - der2); // / fCurrentESD->GetBinWidth(i);
9ca6be09 432
433 // square and weight with the bin width
cfc19dd5 434 chi2 += secDer * secDer; // * fCurrentESD->GetBinWidth(i) * fCurrentESD->GetBinWidth(i-1);
9ca6be09 435
436 //printf("%d --> %f\n", i, secDer);
437 }
438
cfc19dd5 439 return chi2; // 10
3602328d 440}
441
442//____________________________________________________________________
443Double_t AliMultiplicityCorrection::RegularizationEntropy(Double_t *params)
444{
445 // homogenity term for minuit fitting
446 // pure function of the parameters
447 // calculates entropy, from
448 // The method of reduced cross-entropy (M. Schmelling 1993)
449
450 //static Int_t callCount = 0;
451
452 Double_t paramSum = 0;
453 for (Int_t i=0; i<fgMaxParams; ++i)
cfc19dd5 454 paramSum += params[i]; // / fCurrentESD->GetBinWidth(i+1);
3602328d 455
456 Double_t chi2 = 0;
457 for (Int_t i=0; i<fgMaxParams; ++i)
458 {
cfc19dd5 459 Double_t tmp = params[i] / paramSum; // / fCurrentESD->GetBinWidth(i+1);
3602328d 460 if (tmp > 0)
461 {
462 chi2 += tmp * log(tmp); /* * fCurrentESD->GetBinWidth(i+1);
463 // /\
464 // ------------------------------------
465 // TODO WARNING the absolute fitting values seem to depend on that value |
466 // this should be impossible...
467 //if ((callCount % 100) == 0)
468 // printf("%f => %f\n", params[i], tmp * log(tmp)); */
469 }
470 }
471 //if ((callCount++ % 100) == 0)
472 // printf("\n");
473
cfc19dd5 474 return chi2; // 1000
9ca6be09 475}
476
0a173978 477//____________________________________________________________________
478void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t)
479{
480 //
481 // fit function for minuit
cfc19dd5 482 // does: (m - Ad)W(m - Ad) where m = measured, A correlation matrix, d = guess, W = covariance matrix
0a173978 483 //
484
0a173978 485 static Int_t callCount = 0;
486
487 Double_t chi2FromFit = 0;
488
cfc19dd5 489 // d
490 TVectorF paramsVector(fgMaxParams);
491 for (Int_t i=0; i<fgMaxParams; ++i)
492 paramsVector[i] = params[i];
493
494 // Ad
495 paramsVector = (*fCorrelationMatrix) * paramsVector;
496
497 // Ad - m
498 paramsVector -= (*fCurrentESDVector);
499
500 TVectorF copy(paramsVector);
501
502 // (Ad - m) W
503 paramsVector *= (*fCorrelationCovarianceMatrix);
504
505 // (Ad - m) W (Ad - m)
506 chi2FromFit = paramsVector * copy;
507
508 /*printf("chi2FromFit = %f\n", chi2FromFit);
509 chi2FromFit = 0;
510
0a173978 511 // loop over multiplicity (x axis of fMultiplicityESD)
9ca6be09 512 for (Int_t i=1; i<=fCurrentESD->GetNbinsX(); ++i)
0a173978 513 {
cfc19dd5 514 if (i > fCurrentCorrelation->GetNbinsY() || i > fgMaxParams)
0a173978 515 break;
516
517 Double_t sum = 0;
518 //Double_t error = 0;
519 // loop over generated particles (x axis of fCorrelation) that resulted in reconstruction of i tracks
9ca6be09 520 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsX(); ++j)
0a173978 521 {
522 if (j > fgMaxParams)
523 break;
524
9ca6be09 525 sum += fCurrentCorrelation->GetBinContent(j, i) * params[j-1];
0a173978 526
527 //if (params[j-1] > 0)
9ca6be09 528 // error += fCurrentCorrelation->GetBinError(j, i) * fCurrentCorrelation->GetBinError(j, i) * params[j-1];
0a173978 529 }
530
cfc19dd5 531 //printf("%f\n", sum);
532
9ca6be09 533 Double_t diff = fCurrentESD->GetBinContent(i) - sum;
3602328d 534
535 // include error
536 if (fCurrentESD->GetBinError(i) > 0)
537 diff /= fCurrentESD->GetBinError(i);
538
cfc19dd5 539 //if (fCurrentESD->GetBinContent(i) > 0)
540 // diff /= fCurrentESD->GetBinContent(i);
541 //else
542 // diff /= fCurrentESD->Integral();
9ca6be09 543
544 // weight with bin width
cfc19dd5 545 //diff *= fCurrentESD->GetBinWidth(i);
9ca6be09 546
547 //error = TMath::Sqrt(error) + fCurrentESD->GetBinError(i);
0a173978 548 //if (error <= 0)
549 // error = 1;
550
551 //Double_t tmp = diff / error;
552 //chi2 += tmp * tmp;
553 chi2FromFit += diff * diff;
554
9ca6be09 555 //printf("\nExpected sum = %f; Diff for bin %d is %f\n**********************************\n", fCurrentESD->GetBinContent(i), i, diff);
0a173978 556 //printf("Diff for bin %d is %f\n", i, diff);
557 }
558
cfc19dd5 559 printf("chi2FromFit = %f\n", chi2FromFit);*/
560
561 Double_t penaltyVal = 0;
562
563 switch (fRegularizationType)
564 {
565 case kNone: break;
566 case kPol0: penaltyVal = RegularizationPol0(params); break;
567 case kPol1: penaltyVal = RegularizationPol1(params); break;
568 case kCurvature: penaltyVal = RegularizationTotalCurvature(params); break;
569 case kEntropy: penaltyVal = RegularizationEntropy(params); break;
570 }
571
572 penaltyVal *= fRegularizationWeight;
0a173978 573
3602328d 574 chi2 = chi2FromFit + penaltyVal;
0a173978 575
cfc19dd5 576 if ((callCount++ % 1000) == 0)
0a173978 577 printf("%f %f --> %f\n", chi2FromFit, penaltyVal, chi2);
578}
579
580//____________________________________________________________________
cfc19dd5 581void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
0a173978 582{
583 //
9ca6be09 584 // fills fCurrentESD, fCurrentCorrelation
585 // resets fMultiplicityESDCorrected
0a173978 586 //
587
588 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
9ca6be09 589 fMultiplicityESDCorrected[correlationID]->Reset();
0a173978 590
9ca6be09 591 fCurrentESD = fMultiplicityESD[inputRange]->ProjectionY();
592 fCurrentESD->Sumw2();
0a173978 593
594 // empty under/overflow bins in x, otherwise Project3D takes them into account
595 TH3* hist = fCorrelation[correlationID];
596 for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
597 {
598 for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
599 {
600 hist->SetBinContent(0, y, z, 0);
601 hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
602 }
603 }
604
9ca6be09 605 fCurrentCorrelation = hist->Project3D("zy");
606 fCurrentCorrelation->Sumw2();
cfc19dd5 607
608 fCurrentEfficiency = fMultiplicityVtx[inputRange]->ProjectionY("CurrentEfficiency");
609 TH2* divisor = 0;
610 switch (eventType)
611 {
612 case kTrVtx : divisor = fMultiplicityVtx[inputRange]; break;
613 case kMB: divisor = fMultiplicityMB[inputRange]; break;
614 case kINEL: divisor = fMultiplicityINEL[inputRange]; break;
615 }
616 fCurrentEfficiency->Divide(divisor->ProjectionY());
9ca6be09 617}
618
619//____________________________________________________________________
620void AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, Bool_t check)
621{
622 //
623 // correct spectrum using minuit chi2 method
624 //
625 // if check is kTRUE the input MC solution (by definition the right solution) is used
626 // no fit is made and just the chi2 is calculated
627 //
628
629 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
630 Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
631
cfc19dd5 632 SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx);
633
634 fCorrelationMatrix = new TMatrixF(fgMaxParams, fgMaxParams);
635 fCorrelationCovarianceMatrix = new TMatrixF(fgMaxParams, fgMaxParams);
636 fCurrentESDVector = new TVectorF(fgMaxParams);
0a173978 637
638 // normalize correction for given nPart
9ca6be09 639 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
0a173978 640 {
9ca6be09 641 Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
0a173978 642 if (sum <= 0)
643 continue;
9ca6be09 644 Float_t maxValue = 0;
645 Int_t maxBin = -1;
646 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
0a173978 647 {
9ca6be09 648 // find most probably value
649 if (maxValue < fCurrentCorrelation->GetBinContent(i, j))
650 {
651 maxValue = fCurrentCorrelation->GetBinContent(i, j);
652 maxBin = j;
653 }
654
0a173978 655 // npart sum to 1
9ca6be09 656 fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
657 fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
cfc19dd5 658
659 if (i <= fgMaxParams && j <= fgMaxParams)
660 (*fCorrelationMatrix)(j-1, i-1) = fCurrentCorrelation->GetBinContent(i, j);
0a173978 661 }
9ca6be09 662
cfc19dd5 663 //printf("MPV for Ntrue = %f is %f\n", fCurrentCorrelation->GetXaxis()->GetBinCenter(i), fCurrentCorrelation->GetYaxis()->GetBinCenter(maxBin));
0a173978 664 }
665
666 // small hack to get around charge conservation for full phase space ;-)
667 /*if (fullPhaseSpace)
668 {
669 for (Int_t i=2; i<=50; i+=2)
9ca6be09 670 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
671 fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i-1, j));
0a173978 672 }*/
673
cfc19dd5 674 /*TCanvas* canvas1 = new TCanvas("ApplyMinuitFit", "ApplyMinuitFit", 800, 400);
0a173978 675 canvas1->Divide(2, 1);
676 canvas1->cd(1);
9ca6be09 677 fCurrentESD->DrawCopy();
0a173978 678 canvas1->cd(2);
cfc19dd5 679 fCurrentCorrelation->DrawCopy("COLZ");*/
0a173978 680
681 // Initialize TMinuit via generic fitter interface
682 TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgMaxParams);
683
684 minuit->SetFCN(MinuitFitFunction);
685
3602328d 686 static Double_t results[fgMaxParams];
cfc19dd5 687 //printf("%x\n", (void*) results);
0a173978 688
cfc19dd5 689 TH1* proj = fMultiplicityVtx[mcTarget]->ProjectionY();
0a173978 690 for (Int_t i=0; i<fgMaxParams; ++i)
691 {
692 //results[i] = 1;
9ca6be09 693 results[i] = fCurrentESD->GetBinContent(i+1);
b477f4f2 694 if (results[i] < 1e-2)
695 results[i] = 1e-2;
9ca6be09 696 if (check)
697 results[i] = proj->GetBinContent(i+1);
b477f4f2 698 minuit->SetParameter(i, Form("param%d", i), results[i], results[i] / 10, 0.001, fCurrentESD->GetMaximum() * 10);
cfc19dd5 699
700 (*fCurrentESDVector)[i] = fCurrentESD->GetBinContent(i+1);
701 if (fCurrentESD->GetBinError(i+1) > 0)
702 (*fCorrelationCovarianceMatrix)(i, i) = 1.0 / (fCurrentESD->GetBinError(i+1) * fCurrentESD->GetBinError(i+1));
703
3602328d 704 //minuit->SetParameter(i, Form("param%d", i), fCurrentESD->GetBinWidth(i+1), 0.01, 0.001, 50);
0a173978 705 }
3602328d 706 minuit->SetParameter(0, "param0", results[1], ((results[1] > 1) ? (results[1] / 10) : 0), 0.001, fCurrentESD->GetMaximum() * 10);
707 //minuit->SetParameter(0, "param0", 1, 0, 1, 1);
0a173978 708
709 Int_t dummy;
710 Double_t chi2;
711 MinuitFitFunction(dummy, 0, chi2, results, 0);
712 printf("Chi2 of right solution is = %f\n", chi2);
713
9ca6be09 714 if (check)
715 return;
0a173978 716
717 Double_t arglist[100];
3602328d 718 arglist[0] = 0;
719 minuit->ExecuteCommand("SET PRINT", arglist, 1);
0a173978 720 minuit->ExecuteCommand("MIGRAD", arglist, 1);
721 //minuit->ExecuteCommand("MINOS", arglist, 0);
722
723 for (Int_t i=0; i<fgMaxParams; ++i)
724 {
725 results[i] = minuit->GetParameter(i);
726 fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, results[i]);
727 fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0);
728 }
729
3602328d 730 //printf("Penalty is %f\n", RegularizationTotalCurvature(results));
0a173978 731}
732
733//____________________________________________________________________
734void AliMultiplicityCorrection::NormalizeToBinWidth(TH1* hist)
735{
736 //
737 // normalizes a 1-d histogram to its bin width
738 //
739
740 for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
741 {
742 hist->SetBinContent(i, hist->GetBinContent(i) / hist->GetBinWidth(i));
743 hist->SetBinError(i, hist->GetBinError(i) / hist->GetBinWidth(i));
744 }
745}
746
747//____________________________________________________________________
748void AliMultiplicityCorrection::NormalizeToBinWidth(TH2* hist)
749{
750 //
751 // normalizes a 2-d histogram to its bin width (x width * y width)
752 //
753
754 for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
755 for (Int_t j=1; j<=hist->GetNbinsY(); ++j)
756 {
757 Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j);
758 hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor);
759 hist->SetBinError(i, j, hist->GetBinError(i, j) / factor);
760 }
761}
762
763//____________________________________________________________________
764void AliMultiplicityCorrection::ApplyMinuitFitAll()
765{
766 //
767 // fit all eta ranges
768 //
769
770 for (Int_t i=0; i<kESDHists; ++i)
771 {
772 ApplyMinuitFit(i, kFALSE);
773 ApplyMinuitFit(i, kTRUE);
774 }
775}
776
777//____________________________________________________________________
778void AliMultiplicityCorrection::DrawHistograms()
779{
780 printf("ESD:\n");
781
782 TCanvas* canvas1 = new TCanvas("fMultiplicityESD", "fMultiplicityESD", 900, 600);
783 canvas1->Divide(3, 2);
784 for (Int_t i = 0; i < kESDHists; ++i)
785 {
786 canvas1->cd(i+1);
787 fMultiplicityESD[i]->DrawCopy("COLZ");
788 printf("%d --> %f\n", i, (Float_t) fMultiplicityESD[i]->ProjectionY()->GetMean());
789 }
790
cfc19dd5 791 printf("Vtx:\n");
0a173978 792
793 TCanvas* canvas2 = new TCanvas("fMultiplicityMC", "fMultiplicityMC", 900, 600);
794 canvas2->Divide(3, 2);
795 for (Int_t i = 0; i < kMCHists; ++i)
796 {
797 canvas2->cd(i+1);
cfc19dd5 798 fMultiplicityVtx[i]->DrawCopy("COLZ");
799 printf("%d --> %f\n", i, fMultiplicityVtx[i]->ProjectionY()->GetMean());
0a173978 800 }
801
802 TCanvas* canvas3 = new TCanvas("fCorrelation", "fCorrelation", 900, 900);
803 canvas3->Divide(3, 3);
804 for (Int_t i = 0; i < kCorrHists; ++i)
805 {
806 canvas3->cd(i+1);
b477f4f2 807 TH3* hist = dynamic_cast<TH3*> (fCorrelation[i]->Clone());
808 for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
809 {
810 for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
811 {
812 hist->SetBinContent(0, y, z, 0);
813 hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
814 }
815 }
816 TH1* proj = hist->Project3D("zy");
0a173978 817 proj->DrawCopy("COLZ");
818 }
819}
820
821//____________________________________________________________________
cfc19dd5 822void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t normalizeESD, TH1* mcHist)
0a173978 823{
cfc19dd5 824 Int_t esdCorrId = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
825
9ca6be09 826 TString tmpStr;
cfc19dd5 827 tmpStr.Form("%s_DrawComparison_%d", name, esdCorrId);
0a173978 828
b477f4f2 829 // calculate residual
830
831 // for that we convolute the response matrix with the gathered result
832 // if normalizeESD is not set, the histogram is already normalized, this needs to be passed to CalculateMultiplicityESD
833 TH2* convoluted = CalculateMultiplicityESD(fMultiplicityESDCorrected[esdCorrId], esdCorrId, !normalizeESD);
834 TH1* proj2 = convoluted->ProjectionY("proj2", -1, -1, "e");
835 NormalizeToBinWidth(proj2);
836 TH1* residual = convoluted->ProjectionY("residual", -1, -1, "e");
837 residual->SetTitle("Residuals");
838
839 residual->Add(fCurrentESD, -1);
840 residual->Divide(residual, fCurrentESD, 1, 1, "B");
841
842 // TODO fix errors
843 for (Int_t i=1; i<residual->GetNbinsX(); ++i)
844 {
845 proj2->SetBinError(i, 0);
846 residual->SetBinError(i, 0);
847 }
848
cfc19dd5 849 TCanvas* canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 800);
b477f4f2 850 canvas1->Divide(2, 2);
0a173978 851
852 canvas1->cd(1);
cfc19dd5 853 TH1* proj = (TH1*) mcHist->Clone("proj");
0a173978 854 NormalizeToBinWidth(proj);
9ca6be09 855
856 if (normalizeESD)
857 NormalizeToBinWidth(fMultiplicityESDCorrected[esdCorrId]);
858
cfc19dd5 859 proj->GetXaxis()->SetRangeUser(0, 150);
860 proj->DrawCopy("HIST");
0a173978 861 gPad->SetLogy();
862
cfc19dd5 863 //fMultiplicityESDCorrected[esdCorrId]->SetMarkerStyle(3);
864 fMultiplicityESDCorrected[esdCorrId]->SetLineColor(2);
865 fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAME HIST");
866
867 Float_t chi2 = 0;
868 for (Int_t i=1; i<=100; ++i)
869 {
870 if (fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i) != 0)
871 {
872 Float_t value = (proj->GetBinContent(i) - fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i)) / fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i);
873 chi2 += value * value;
874 }
875 }
876
877 printf("Difference (from MC) is %f for bin 1-100\n", chi2);
878 fLastChi2MC = chi2;
879
880 canvas1->cd(2);
881
9ca6be09 882 NormalizeToBinWidth(fCurrentESD);
cfc19dd5 883 gPad->SetLogy();
884 fCurrentESD->GetXaxis()->SetRangeUser(0, 150);
885 //fCurrentESD->SetLineColor(2);
886 fCurrentESD->DrawCopy("HIST");
b477f4f2 887
888 proj2->SetLineColor(2);
cfc19dd5 889 //proj2->SetMarkerColor(2);
890 //proj2->SetMarkerStyle(5);
891 proj2->DrawCopy("HIST SAME");
9ca6be09 892
cfc19dd5 893 chi2 = 0;
894 for (Int_t i=1; i<=100; ++i)
895 {
896 if (fCurrentESD->GetBinContent(i) != 0)
897 {
898 Float_t value = (proj2->GetBinContent(i) - fCurrentESD->GetBinContent(i)) / fCurrentESD->GetBinContent(i);
899 chi2 += value * value;
900 }
901 }
0a173978 902
cfc19dd5 903 printf("Difference (Residuals) is %f for bin 1-100\n", chi2);
904 fLastChi2Residuals = chi2;
0a173978 905
906 canvas1->cd(3);
907 TH1* clone = dynamic_cast<TH1*> (proj->Clone("clone"));
908 clone->Divide(fMultiplicityESDCorrected[esdCorrId]);
909 clone->SetTitle("Ratio;Npart;MC/ESD");
910 clone->GetYaxis()->SetRangeUser(0.8, 1.2);
cfc19dd5 911 clone->GetXaxis()->SetRangeUser(0, 150);
0a173978 912 clone->DrawCopy("HIST");
913
914 /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85);
915 legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected");
916 legend->AddEntry(fMultiplicityMC, "MC");
917 legend->AddEntry(fMultiplicityESD, "ESD");
918 legend->Draw();*/
3602328d 919
b477f4f2 920 canvas1->cd(4);
921
922 residual->GetYaxis()->SetRangeUser(-0.2, 0.2);
cfc19dd5 923 residual->GetXaxis()->SetRangeUser(0, 150);
b477f4f2 924 residual->DrawCopy();
3602328d 925
b477f4f2 926 canvas1->SaveAs(Form("%s.gif", canvas1->GetName()));
0a173978 927}
928
929//____________________________________________________________________
cfc19dd5 930void AliMultiplicityCorrection::GetComparisonResults(Float_t& mc, Float_t& residuals)
931{
932 // Returns the chi2 between the MC and the unfolded ESD as well as between the ESD and the folded unfolded ESD
933 // These values are computed during DrawComparison, thus this function picks up the
934 // last calculation
935
936 mc = fLastChi2MC;
937 residuals = fLastChi2Residuals;
938}
939
940
941//____________________________________________________________________
942TH2F* AliMultiplicityCorrection::GetMultiplicityMC(Int_t i, EventType eventType)
943{
944 //
945 // returns the corresponding MC spectrum
946 //
947
948 switch (eventType)
949 {
950 case kTrVtx : return fMultiplicityVtx[i]; break;
951 case kMB : return fMultiplicityMB[i]; break;
952 case kINEL : return fMultiplicityINEL[i]; break;
953 }
954
955 return 0;
956}
957
958//____________________________________________________________________
959void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar)
0a173978 960{
961 //
962 // correct spectrum using bayesian method
963 //
964
965 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
0a173978 966
cfc19dd5 967 SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
968
969 // TODO should be taken from correlation map
970 TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
971
972 //new TCanvas;
973 //sumHist->Draw();
0a173978 974
975 // normalize correction for given nPart
9ca6be09 976 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
0a173978 977 {
cfc19dd5 978 //Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
979 Double_t sum = sumHist->GetBinContent(i);
0a173978 980 if (sum <= 0)
981 continue;
9ca6be09 982 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
0a173978 983 {
984 // npart sum to 1
9ca6be09 985 fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
986 fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
0a173978 987 }
988 }
989
cfc19dd5 990 //new TCanvas;
991 //fCurrentCorrelation->Draw("COLZ");
992 //return;
993
994 // normalize correction for given nTrack
995 /*for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
996 {
997 Double_t sum = fCurrentCorrelation->Integral(1, fCurrentCorrelation->GetNbinsX(), j, j);
998 if (sum <= 0)
999 continue;
1000 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsY(); ++i)
1001 {
1002 // ntrack sum to 1
1003 fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1004 fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1005 }
1006 }*/
1007
1008 // FAKE
1009 fCurrentEfficiency = ((TH2*) fCurrentCorrelation)->ProjectionX("eff");
1010 //new TCanvas;
1011 //fCurrentEfficiency->Draw();
6127aca6 1012
6127aca6 1013 // pick prior distribution
9ca6be09 1014 TH1F* hPrior = (TH1F*)fCurrentESD->Clone("prior");
cfc19dd5 1015 Float_t norm = hPrior->Integral();
6127aca6 1016 for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
6127aca6 1017 hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm);
1018
cfc19dd5 1019 // define temp hist
1020 TH1F* hTemp = (TH1F*)fCurrentESD->Clone("temp");
1021 hTemp->Reset();
6127aca6 1022
cfc19dd5 1023 // just a shortcut
1024 TH2F* hResponse = (TH2F*) fCurrentCorrelation;
6127aca6 1025
cfc19dd5 1026 // histogram to store the inverse response calculated from bayes theorem (from prior and response) IR
6127aca6 1027 TH2F* hInverseResponseBayes = (TH2F*)hResponse->Clone("pji");
6127aca6 1028 hInverseResponseBayes->Reset();
b477f4f2 1029
6127aca6 1030 // unfold...
3602328d 1031 Int_t iterations = 20;
cfc19dd5 1032 for (Int_t i=0; i<iterations; i++)
1033 {
1034 //printf(" iteration %i \n", i);
1035
6127aca6 1036 // calculate IR from Beyes theorem
1037 // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k)
9ca6be09 1038
cfc19dd5 1039 for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
1040 {
1041 Float_t norm = 0;
1042 for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
1043 norm += hResponse->GetBinContent(t,m) * hPrior->GetBinContent(t);
1044 for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
1045 {
6127aca6 1046 if (norm==0)
cfc19dd5 1047 hInverseResponseBayes->SetBinContent(t,m,0);
6127aca6 1048 else
cfc19dd5 1049 hInverseResponseBayes->SetBinContent(t,m, hResponse->GetBinContent(t,m) * hPrior->GetBinContent(t)/norm);
6127aca6 1050 }
1051 }
cfc19dd5 1052
1053 for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
1054 {
1055 Float_t value = 0;
6127aca6 1056 for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
cfc19dd5 1057 value += fCurrentESD->GetBinContent(m) * hInverseResponseBayes->GetBinContent(t,m);
1058
1059 /*if (eventType == kTrVtx)
1060 {
1061 hTemp->SetBinContent(t, value);
1062 }
1063 else*/
1064 {
1065 if (fCurrentEfficiency->GetBinContent(t) > 0)
1066 hTemp->SetBinContent(t, value / fCurrentEfficiency->GetBinContent(t));
1067 else
1068 hTemp->SetBinContent(t, 0);
1069 }
1070 }
1071
1072 // this is the last guess, fill before (!) smoothing
1073 for (Int_t t=1; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++)
1074 {
1075 fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTemp->GetBinContent(t));
1076 fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t)); // TODO
1077
1078 //printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
6127aca6 1079 }
9ca6be09 1080
b477f4f2 1081 // regularization (simple smoothing)
cfc19dd5 1082 TH1F* hTrueSmoothed = (TH1F*) hTemp->Clone("truesmoothed");
9ca6be09 1083
cfc19dd5 1084 for (Int_t t=2; t<hTrueSmoothed->GetNbinsX(); t++)
1085 {
1086 Float_t average = (hTemp->GetBinContent(t-1) / hTemp->GetBinWidth(t-1)
1087 + hTemp->GetBinContent(t) / hTemp->GetBinWidth(t)
1088 + hTemp->GetBinContent(t+1) / hTemp->GetBinWidth(t+1)) / 3.;
9ca6be09 1089 average *= hTrueSmoothed->GetBinWidth(t);
1090
6127aca6 1091 // weight the average with the regularization parameter
cfc19dd5 1092 hTrueSmoothed->SetBinContent(t, (1-regPar) * hTemp->GetBinContent(t) + regPar * average);
6127aca6 1093 }
9ca6be09 1094
cfc19dd5 1095 // calculate chi2 (change from last iteration)
1096 Double_t chi2 = 0;
1097
9ca6be09 1098 // use smoothed true (normalized) as new prior
cfc19dd5 1099 Float_t norm = 1;
1100 //for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
1101 // norm = norm + hTrueSmoothed->GetBinContent(t);
1102
1103 for (Int_t t=1; t<hTrueSmoothed->GetNbinsX(); t++)
1104 {
1105 Float_t newValue = hTrueSmoothed->GetBinContent(t)/norm;
1106 Float_t diff = hPrior->GetBinContent(t) - newValue;
1107 chi2 += (Double_t) diff * diff;
1108
1109 hPrior->SetBinContent(t, newValue);
6127aca6 1110 }
9ca6be09 1111
cfc19dd5 1112 //printf("Chi2 of %d iteration = %.10f\n", i, chi2);
1113
1114 //if (i % 5 == 0)
1115 // DrawComparison(Form("Bayesian_%d", i), mcTarget, correlationID, kTRUE, eventType);
1116
9ca6be09 1117 delete hTrueSmoothed;
6127aca6 1118 } // end of iterations
1119
6127aca6 1120
cfc19dd5 1121 return;
1122
1123 // ********
1124 // Calculate the covariance matrix, all arguments are taken from NIM,A362,487-498,1995
1125
1126 printf("Calculating covariance matrix. This may take some time...\n");
1127
1128 Int_t xBins = hInverseResponseBayes->GetNbinsX();
1129 Int_t yBins = hInverseResponseBayes->GetNbinsY();
1130
1131 // calculate "unfolding matrix" Mij
1132 Float_t matrixM[251][251];
1133 for (Int_t i=1; i<=xBins; i++)
1134 {
1135 for (Int_t j=1; j<=yBins; j++)
1136 {
1137 if (fCurrentEfficiency->GetBinContent(i) > 0)
1138 matrixM[i-1][j-1] = hInverseResponseBayes->GetBinContent(i, j) / fCurrentEfficiency->GetBinContent(i);
1139 else
1140 matrixM[i-1][j-1] = 0;
1141 }
1142 }
1143
1144 Float_t* vectorn = new Float_t[yBins];
1145 for (Int_t j=1; j<=yBins; j++)
1146 vectorn[j-1] = fCurrentESD->GetBinContent(j);
1147
1148 // first part of covariance matrix, depends on input distribution n(E)
1149 Float_t cov1[251][251];
1150
1151 Float_t nEvents = fCurrentESD->Integral(); // N
1152
1153 xBins = 20;
1154 yBins = 20;
1155
1156 for (Int_t k=0; k<xBins; k++)
1157 {
1158 printf("In Cov1: %d\n", k);
1159 for (Int_t l=0; l<yBins; l++)
1160 {
1161 cov1[k][l] = 0;
1162
1163 // sum_j Mkj Mlj n(Ej) * (1 - n(Ej) / N)
1164 for (Int_t j=0; j<yBins; j++)
1165 cov1[k][l] += matrixM[k][j] * matrixM[l][j] * vectorn[j]
1166 * (1.0 - vectorn[j] / nEvents);
1167
1168 // - sum_i,j (i != j) Mki Mlj n(Ei) n(Ej) / N
1169 for (Int_t i=0; i<yBins; i++)
1170 for (Int_t j=0; j<yBins; j++)
1171 {
1172 if (i == j)
1173 continue;
1174 cov1[k][l] -= matrixM[k][i] * matrixM[l][j] * vectorn[i]
1175 * vectorn[j] / nEvents;
1176 }
1177 }
1178 }
1179
1180 printf("Cov1 finished\n");
1181
1182 TH2F* cov = (TH2F*) hInverseResponseBayes->Clone("cov");
1183 cov->Reset();
1184
1185 for (Int_t i=1; i<=xBins; i++)
1186 for (Int_t j=1; j<=yBins; j++)
1187 cov->SetBinContent(i, j, cov1[i-1][j-1]);
1188
1189 new TCanvas;
1190 cov->Draw("COLZ");
1191
1192 // second part of covariance matrix, depends on response matrix
1193 Float_t cov2[251][251];
1194
1195 // Cov[P(Er|Cu), P(Es|Cu)] term
1196 Float_t covTerm[100][100][100];
1197 for (Int_t r=0; r<yBins; r++)
1198 for (Int_t u=0; u<xBins; u++)
1199 for (Int_t s=0; s<yBins; s++)
1200 {
1201 if (r == s)
1202 covTerm[r][u][s] = 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
1203 * (1.0 - hResponse->GetBinContent(u+1, r+1));
1204 else
1205 covTerm[r][u][s] = - 1.0 / sumHist->GetBinContent(u+1) * hResponse->GetBinContent(u+1, r+1)
1206 * hResponse->GetBinContent(u+1, s+1);
1207 }
1208
1209 for (Int_t k=0; k<xBins; k++)
1210 for (Int_t l=0; l<yBins; l++)
1211 {
1212 cov2[k][l] = 0;
1213 printf("In Cov2: %d %d\n", k, l);
1214 for (Int_t i=0; i<yBins; i++)
1215 for (Int_t j=0; j<yBins; j++)
1216 {
1217 //printf("In Cov2: %d %d %d %d\n", k, l, i, j);
1218 // calculate Cov(Mki, Mlj) = sum{ru},{su} ...
1219 Float_t tmpCov = 0;
1220 for (Int_t r=0; r<yBins; r++)
1221 for (Int_t u=0; u<xBins; u++)
1222 for (Int_t s=0; s<yBins; s++)
1223 {
1224 if (hResponse->GetBinContent(u+1, r+1) == 0 || hResponse->GetBinContent(u+1, s+1) == 0
1225 || hResponse->GetBinContent(u+1, i+1) == 0)
1226 continue;
1227
1228 tmpCov += BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, k, i, r, u)
1229 * BayesCovarianceDerivate(matrixM, hResponse, fCurrentEfficiency, l, j, s, u)
1230 * covTerm[r][u][s];
1231 }
1232
1233 cov2[k][l] += fCurrentESD->GetBinContent(i+1) * fCurrentESD->GetBinContent(j+1) * tmpCov;
1234 }
1235 }
1236
1237 printf("Cov2 finished\n");
1238
1239 for (Int_t i=1; i<=xBins; i++)
1240 for (Int_t j=1; j<=yBins; j++)
1241 cov->SetBinContent(i, j, cov1[i-1][j-1] + cov2[i-1][j-1]);
1242
1243 new TCanvas;
1244 cov->Draw("COLZ");
1245}
1246
1247//____________________________________________________________________
1248void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
1249{
1250 //
1251 // correct spectrum using bayesian method
1252 //
1253
1254 Float_t regPar = 0.05;
1255
1256 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1257 Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
1258
1259 SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
1260
1261 // TODO should be taken from correlation map
1262 //TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
1263
1264 // normalize correction for given nPart
1265 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1266 {
1267 Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1268 //Double_t sum = sumHist->GetBinContent(i);
1269 if (sum <= 0)
1270 continue;
1271 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1272 {
1273 // npart sum to 1
1274 fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1275 fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1276 }
6127aca6 1277 }
cfc19dd5 1278
1279 new TCanvas;
1280 fCurrentCorrelation->Draw("COLZ");
1281
1282 // FAKE
1283 fCurrentEfficiency = ((TH2*) fCurrentCorrelation)->ProjectionX("eff");
1284
1285 // pick prior distribution
1286 TH1F* hPrior = (TH1F*)fCurrentESD->Clone("prior");
1287 Float_t norm = 1; //hPrior->Integral();
1288 for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
1289 hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm);
1290
1291 // zero distribution
1292 TH1F* zero = (TH1F*)hPrior->Clone("zero");
1293
1294 // define temp hist
1295 TH1F* hTemp = (TH1F*)fCurrentESD->Clone("temp");
1296 hTemp->Reset();
1297
1298 // just a shortcut
1299 TH2F* hResponse = (TH2F*) fCurrentCorrelation;
1300
1301 // unfold...
1302 Int_t iterations = 20;
1303 for (Int_t i=0; i<iterations; i++)
1304 {
1305 //printf(" iteration %i \n", i);
1306
1307 for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
1308 {
1309 Float_t value = 0;
1310 for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
1311 value += hResponse->GetBinContent(t, m) * hPrior->GetBinContent(t);
1312 hTemp->SetBinContent(m, value);
1313 //printf("%d %f %f %f\n", m, zero->GetBinContent(m), hPrior->GetBinContent(m), value);
1314 }
1315
1316 // regularization (simple smoothing)
1317 TH1F* hTrueSmoothed = (TH1F*) hTemp->Clone("truesmoothed");
1318
1319 for (Int_t t=2; t<hTrueSmoothed->GetNbinsX(); t++)
1320 {
1321 Float_t average = (hTemp->GetBinContent(t-1) / hTemp->GetBinWidth(t-1)
1322 + hTemp->GetBinContent(t) / hTemp->GetBinWidth(t)
1323 + hTemp->GetBinContent(t+1) / hTemp->GetBinWidth(t+1)) / 3.;
1324 average *= hTrueSmoothed->GetBinWidth(t);
1325
1326 // weight the average with the regularization parameter
1327 hTrueSmoothed->SetBinContent(t, (1-regPar) * hTemp->GetBinContent(t) + regPar * average);
1328 }
1329
1330 for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
1331 hTemp->SetBinContent(m, zero->GetBinContent(m) + hPrior->GetBinContent(m) - hTrueSmoothed->GetBinContent(m));
1332
1333 // fill guess
1334 for (Int_t t=1; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++)
1335 {
1336 fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTemp->GetBinContent(t));
1337 fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t)); // TODO
1338
1339 //printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
1340 }
1341
1342
1343 // calculate chi2 (change from last iteration)
1344 Double_t chi2 = 0;
1345
1346 // use smoothed true (normalized) as new prior
1347 Float_t norm = 1; //hTrueSmoothed->Integral();
1348
1349 for (Int_t t=1; t<hTrueSmoothed->GetNbinsX(); t++)
1350 {
1351 Float_t newValue = hTrueSmoothed->GetBinContent(t)/norm;
1352 Float_t diff = hPrior->GetBinContent(t) - newValue;
1353 chi2 += (Double_t) diff * diff;
1354
1355 hPrior->SetBinContent(t, newValue);
1356 }
1357
1358 printf("Chi2 of %d iteration = %.10f\n", i, chi2);
1359
1360 if (i % 5 == 0)
1361 DrawComparison(Form("Laszlo_%d", i), inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
1362
1363 delete hTrueSmoothed;
1364 } // end of iterations
1365
1366 DrawComparison("Laszlo", inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
1367}
1368
1369//____________________________________________________________________
1370Float_t AliMultiplicityCorrection::BayesCovarianceDerivate(Float_t matrixM[251][251], TH2* hResponse,
1371 TH1* fCurrentEfficiency, Int_t k, Int_t i, Int_t r, Int_t u)
1372{
1373 //
1374 //
1375 //
1376
1377 Float_t result = 0;
1378
1379 if (k == u && r == i)
1380 result += 1.0 / hResponse->GetBinContent(u+1, r+1);
1381
1382 if (k == u)
1383 result -= 1.0 / fCurrentEfficiency->GetBinContent(u+1);
1384
1385 if (r == i)
1386 result -= matrixM[u][i] * fCurrentEfficiency->GetBinContent(u+1) / hResponse->GetBinContent(u+1, i+1);
1387
1388 result *= matrixM[k][i];
1389
1390 return result;
9ca6be09 1391}
1392
1393//____________________________________________________________________
1394void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace)
1395{
1396 //
1397 // correct spectrum using a simple Gaussian approach, that is model-dependent
1398 //
1399
1400 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1401 Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4);
1402
cfc19dd5 1403 SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx);
9ca6be09 1404
1405 NormalizeToBinWidth((TH2*) fCurrentCorrelation);
1406
1407 TH1D* correction = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianMean"));
1408 correction->SetTitle("GaussianMean");
1409
1410 TH1D* correctionWidth = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianWidth"));
1411 correction->SetTitle("GaussianWidth");
1412
1413 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1414 {
1415 TH1D* proj = (dynamic_cast<TH2*> (fCurrentCorrelation))->ProjectionX("_px", i, i+1);
1416 proj->Fit("gaus", "0Q");
1417 correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
1418 correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
1419
1420 continue;
1421
1422 // draw for debugging
1423 new TCanvas;
1424 proj->DrawCopy();
1425 proj->GetFunction("gaus")->DrawCopy("SAME");
1426 }
1427
1428 TH1* target = fMultiplicityESDCorrected[correlationID];
1429
1430 Int_t nBins = target->GetNbinsX()*10+1;
1431 Float_t* binning = new Float_t[nBins];
1432 for (Int_t i=1; i<=target->GetNbinsX(); ++i)
1433 for (Int_t j=0; j<10; ++j)
1434 binning[(i-1)*10 + j] = target->GetXaxis()->GetBinLowEdge(i) + target->GetXaxis()->GetBinWidth(i) / 10 * j;
1435
1436 binning[nBins-1] = target->GetXaxis()->GetBinUpEdge(target->GetNbinsX());
1437
1438 TH1F* fineBinned = new TH1F("targetFineBinned", "targetFineBinned", nBins-1, binning);
1439
1440 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1441 {
1442 Float_t mean = correction->GetBinContent(i);
1443 Float_t width = correctionWidth->GetBinContent(i);
1444
3602328d 1445 Int_t fillBegin = fineBinned->FindBin(mean - width * 5);
1446 Int_t fillEnd = fineBinned->FindBin(mean + width * 5);
1447 //printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
9ca6be09 1448
1449 for (Int_t j=fillBegin; j <= fillEnd; ++j)
1450 {
1451 fineBinned->AddBinContent(j, TMath::Gaus(fineBinned->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fCurrentESD->GetBinContent(i));
1452 }
1453 }
1454
1455 for (Int_t i=1; i<=target->GetNbinsX(); ++i)
1456 {
1457 Float_t sum = 0;
1458 for (Int_t j=1; j<=10; ++j)
1459 sum += fineBinned->GetBinContent((i-1)*10 + j);
1460 target->SetBinContent(i, sum / 10);
1461 }
1462
1463 delete[] binning;
1464
cfc19dd5 1465 DrawComparison("Gaussian", inputRange, fullPhaseSpace, kFALSE, GetMultiplicityMC(mcTarget, kTrVtx)->ProjectionY());
0a173978 1466}
3602328d 1467
1468//____________________________________________________________________
b477f4f2 1469TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap, Bool_t normalized)
3602328d 1470{
1471 // runs the distribution given in inputMC through the response matrix identified by
1472 // correlationMap and produces a measured distribution
1473 // although it is a TH2F the vertex axis is not used at the moment and all entries are filled in mid-vertex
b477f4f2 1474 // if normalized is set, inputMC is expected to be normalized to the bin width
3602328d 1475
1476 if (!inputMC)
1477 return 0;
1478
1479 if (correlationMap < 0 || correlationMap >= kCorrHists)
1480 return 0;
1481
1482 // empty under/overflow bins in x, otherwise Project3D takes them into account
1483 TH3* hist = fCorrelation[correlationMap];
1484 for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
1485 {
1486 for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
1487 {
1488 hist->SetBinContent(0, y, z, 0);
1489 hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
1490 }
1491 }
1492
1493 TH1* corr = hist->Project3D("zy");
1494 corr->Sumw2();
1495
1496 // normalize correction for given nPart
1497 for (Int_t i=1; i<=corr->GetNbinsX(); ++i)
1498 {
1499 Double_t sum = corr->Integral(i, i, 1, corr->GetNbinsY());
1500 if (sum <= 0)
1501 continue;
1502
1503 for (Int_t j=1; j<=corr->GetNbinsY(); ++j)
1504 {
1505 // npart sum to 1
1506 corr->SetBinContent(i, j, corr->GetBinContent(i, j) / sum);
1507 corr->SetBinError(i, j, corr->GetBinError(i, j) / sum);
1508 }
1509 }
1510
1511 TH2F* target = dynamic_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("%s_measured", inputMC->GetName())));
1512 target->Reset();
1513
1514 for (Int_t meas=1; meas<=corr->GetNbinsY(); ++meas)
1515 {
1516 Float_t measured = 0;
1517 Float_t error = 0;
1518
1519 for (Int_t gen=1; gen<=target->GetNbinsY(); ++gen)
1520 {
1521 Int_t mcGenBin = inputMC->GetXaxis()->FindBin(target->GetYaxis()->GetBinCenter(gen));
1522
b477f4f2 1523 Float_t factor = 1;
1524 if (normalized)
1525 factor = inputMC->GetBinWidth(mcGenBin);
1526
1527 measured += inputMC->GetBinContent(mcGenBin) * factor * corr->GetBinContent(gen, meas);
1528 error += inputMC->GetBinError(mcGenBin) * factor * corr->GetBinContent(gen, meas);
3602328d 1529 }
1530
1531 // bin width of the <measured> axis has to be taken into account
b477f4f2 1532 //measured /= target->GetYaxis()->GetBinWidth(meas);
1533 //error /= target->GetYaxis()->GetBinWidth(meas);
3602328d 1534
cfc19dd5 1535 //printf("%f +- %f ; %f +- %f \n", inputMC->GetBinContent(meas), inputMC->GetBinError(meas), measured, error);
3602328d 1536
1537 target->SetBinContent(target->GetNbinsX() / 2, meas, measured);
1538 target->SetBinError(target->GetNbinsX() / 2, meas, error);
1539 }
1540
1541 return target;
1542}
1543
1544//____________________________________________________________________
1545void AliMultiplicityCorrection::SetGenMeasFromFunc(TF1* inputMC, Int_t id)
1546{
1547 // uses the given function to fill the input MC histogram and generates from that
1548 // the measured histogram by applying the response matrix
1549 // this can be used to evaluate if the methods work indepedently of the input
1550 // distribution
1551 // WARNING does not respect the vertex distribution, just fills central vertex bin
1552
1553 if (!inputMC)
1554 return;
1555
1556 if (id < 0 || id >= kESDHists)
1557 return;
1558
cfc19dd5 1559 TH2F* mc = fMultiplicityINEL[id];
3602328d 1560
1561 mc->Reset();
1562
1563 for (Int_t i=1; i<=mc->GetNbinsY(); ++i)
b477f4f2 1564 {
3602328d 1565 mc->SetBinContent(mc->GetNbinsX() / 2, i, inputMC->Eval(mc->GetYaxis()->GetBinCenter(i)) * mc->GetYaxis()->GetBinWidth(i));
b477f4f2 1566 mc->SetBinError(mc->GetNbinsX() / 2, i, 0);
1567 }
3602328d 1568
1569 //new TCanvas;
1570 //mc->Draw("COLZ");
1571
1572 TH1* proj = mc->ProjectionY();
1573
1574 fMultiplicityESD[id] = CalculateMultiplicityESD(proj, id);
1575}