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