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