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