1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 // This class is used to store correction maps, raw input and results of the multiplicity
19 // measurement with the ITS or TPC
20 // It also contains functions to correct the spectrum using different methods.
21 // e.g. chi2 minimization and bayesian unfolding
23 // Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch
25 #include "AliMultiplicityCorrection.h"
31 #include <TDirectory.h>
36 #include <TCollection.h>
41 #include <TProfile2D.h>
44 ClassImp(AliMultiplicityCorrection)
46 // Defined where the efficiency drops below 1/3
47 // |eta| < 1.4 --> -0.3 ... 0.8
48 // |eta| < 1.3 --> -1.9 ... 2.4
49 // |eta| < 1.0 --> -5.6 ... 6.1
50 //Double_t AliMultiplicityCorrection::fgVtxRangeBegin[kESDHists] = { -10.0, -5.6, -1.9 };
51 //Double_t AliMultiplicityCorrection::fgVtxRangeEnd[kESDHists] = { 10.0, 6.1, 2.4 };
52 Double_t AliMultiplicityCorrection::fgVtxRangeBegin[kESDHists] = { -10.0, -5.5, -1.9 };
53 Double_t AliMultiplicityCorrection::fgVtxRangeEnd[kESDHists] = { 10.0, 5.5, 2.4 };
55 // These are the areas where the quality of the unfolding results are evaluated
56 // Default defined here, call SetQualityRegions to change them
57 // unit is in multiplicity (not in bin!)
58 // SPD: peak area - flat area - low stat area
59 Int_t AliMultiplicityCorrection::fgQualityRegionsB[kQualityRegions] = {1, 20, 70};
60 Int_t AliMultiplicityCorrection::fgQualityRegionsE[kQualityRegions] = {10, 65, 80};
62 //____________________________________________________________________
63 void AliMultiplicityCorrection::SetQualityRegions(Bool_t SPDStudy)
66 // sets the quality region definition to TPC or SPD
71 // SPD: peak area - flat area - low stat area
72 fgQualityRegionsB[0] = 1;
73 fgQualityRegionsE[0] = 10;
75 fgQualityRegionsB[1] = 20;
76 fgQualityRegionsE[1] = 65;
78 fgQualityRegionsB[2] = 70;
79 fgQualityRegionsE[2] = 80;
81 Printf("AliMultiplicityCorrection::SetQualityRegions --> Enabled quality regions for SPD");
85 // TPC: peak area - flat area - low stat area
86 fgQualityRegionsB[0] = 4;
87 fgQualityRegionsE[0] = 12;
89 fgQualityRegionsB[1] = 25;
90 fgQualityRegionsE[1] = 55;
92 fgQualityRegionsB[2] = 88;
93 fgQualityRegionsE[2] = 108;
95 Printf("AliMultiplicityCorrection::SetQualityRegions --> Enabled quality regions for TPC");
99 //____________________________________________________________________
100 AliMultiplicityCorrection::AliMultiplicityCorrection() :
101 TNamed(), fCurrentESD(0), fCurrentCorrelation(0), fCurrentEfficiency(0), fLastBinLimit(0), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0), fVtxBegin(0), fVtxEnd(0)
104 // default constructor
107 for (Int_t i = 0; i < kESDHists; ++i)
109 fMultiplicityESD[i] = 0;
110 fTriggeredEvents[i] = 0;
111 fNoVertexEvents[i] = 0;
114 for (Int_t i = 0; i < kMCHists; ++i)
116 fMultiplicityVtx[i] = 0;
117 fMultiplicityMB[i] = 0;
118 fMultiplicityINEL[i] = 0;
119 fMultiplicityNSD[i] = 0;
122 for (Int_t i = 0; i < kCorrHists; ++i)
125 fMultiplicityESDCorrected[i] = 0;
128 for (Int_t i = 0; i < kQualityRegions; ++i)
132 //____________________________________________________________________
133 AliMultiplicityCorrection* AliMultiplicityCorrection::Open(const char* fileName, const char* folderName)
135 // opens the given file, reads the multiplicity from the given folder and returns the object
137 TFile* file = TFile::Open(fileName);
140 Printf("ERROR: Could not open %s", fileName);
144 Printf("AliMultiplicityCorrection::Open: Reading file %s", fileName);
146 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folderName, folderName);
147 mult->LoadHistograms();
149 // TODO closing the file does not work here, because the histograms cannot be read anymore. LoadHistograms need to be adapted
154 //____________________________________________________________________
155 AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const Char_t* title) :
158 fCurrentCorrelation(0),
159 fCurrentEfficiency(0),
163 fLastChi2Residuals(0),
172 // do not add this hists to the directory
173 Bool_t oldStatus = TH1::AddDirectoryStatus();
174 TH1::AddDirectory(kFALSE);
176 /*Float_t binLimitsVtx[] = {-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10};
177 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,
178 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
179 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
180 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
181 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5,
182 50.5, 55.5, 60.5, 65.5, 70.5, 75.5, 80.5, 85.5, 90.5, 95.5,
183 100.5, 105.5, 110.5, 115.5, 120.5, 125.5, 130.5, 135.5, 140.5, 145.5,
184 150.5, 160.5, 170.5, 180.5, 190.5, 200.5, 210.5, 220.5, 230.5, 240.5,
185 250.5, 275.5, 300.5, 325.5, 350.5, 375.5, 400.5, 425.5, 450.5, 475.5,
187 //525.5, 550.5, 575.5, 600.5, 625.5, 650.5, 675.5, 700.5, 725.5,
188 //750.5, 775.5, 800.5, 825.5, 850.5, 875.5, 900.5, 925.5, 950.5, 975.5,
191 #define VTXBINNING 10, binLimitsVtx
192 #define NBINNING fgkMaxParams, binLimitsN*/
194 #define NBINNING 201, -0.5, 200.5
196 for (Int_t i = 0; i < kESDHists; ++i)
198 fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;measured multiplicity;Count", 1, fgVtxRangeBegin[i], fgVtxRangeEnd[i], NBINNING);
199 fTriggeredEvents[i] = new TH1F(Form("fTriggeredEvents%d", i), "fTriggeredEvents;measured multiplicity;Count", NBINNING);
200 fNoVertexEvents[i] = new TH1F(Form("fNoVertexEvents%d", i), "fNoVertexEvents;generated multiplicity;Count", NBINNING);
203 for (Int_t i = 0; i < kMCHists; ++i)
205 fMultiplicityVtx[i] = dynamic_cast<TH2F*> (fMultiplicityESD[i%3]->Clone(Form("fMultiplicityVtx%d", i)));
206 fMultiplicityVtx[i]->SetTitle("fMultiplicityVtx;vtx-z;Npart");
208 fMultiplicityMB[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[i]->Clone(Form("fMultiplicityMB%d", i)));
209 fMultiplicityMB[i]->SetTitle("fMultiplicityMB");
211 fMultiplicityINEL[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[i]->Clone(Form("fMultiplicityINEL%d", i)));
212 fMultiplicityINEL[i]->SetTitle("fMultiplicityINEL");
214 fMultiplicityNSD[i] = dynamic_cast<TH2F*> (fMultiplicityVtx[i]->Clone(Form("fMultiplicityNSD%d", i)));
215 fMultiplicityNSD[i]->SetTitle("fMultiplicityNSD");
218 for (Int_t i = 0; i < kCorrHists; ++i)
220 fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;true multiplicity;measured multiplicity", 1, fgVtxRangeBegin[i%3], fgVtxRangeEnd[i%3], NBINNING, NBINNING);
221 fMultiplicityESDCorrected[i] = new TH1F(Form("fMultiplicityESDCorrected%d", i), "fMultiplicityESDCorrected;true multiplicity;Count", NBINNING);
224 for (Int_t i = 0; i < kQualityRegions; ++i)
227 TH1::AddDirectory(oldStatus);
229 AliUnfolding::SetNbins(120, 120);
230 AliUnfolding::SetSkipBinsBegin(1);
231 //AliUnfolding::SetNormalizeInput(kTRUE);
234 //____________________________________________________________________
235 void AliMultiplicityCorrection::Rebin2DY(TH2F*& hist, Int_t nBins, Double_t* newBins) const
238 // rebins the y axis of a two-dimensional histogram giving variable size binning (missing in ROOT v5/25/02)
241 TH2F* temp = new TH2F(hist->GetName(), Form("%s;%s;%s", hist->GetTitle(), hist->GetXaxis()->GetTitle(), hist->GetYaxis()->GetTitle()), hist->GetNbinsX(), hist->GetXaxis()->GetXmin(), hist->GetXaxis()->GetXmax(), nBins, newBins);
243 for (Int_t x=0; x<=hist->GetNbinsX()+1; x++)
244 for (Int_t y=0; y<=hist->GetNbinsY()+1; y++)
245 temp->Fill(hist->GetXaxis()->GetBinCenter(x), hist->GetYaxis()->GetBinCenter(y), hist->GetBinContent(x, y));
247 for (Int_t x=0; x<=temp->GetNbinsX()+1; x++)
248 for (Int_t y=0; y<=temp->GetNbinsY()+1; y++)
249 temp->SetBinError(x, y, TMath::Sqrt(temp->GetBinContent(x, y)));
255 //____________________________________________________________________
256 void AliMultiplicityCorrection::Rebin3DY(TH3F*& hist, Int_t nBins, Double_t* newBins) const
259 // rebins the y axis of a three-dimensional histogram giving variable size binning (missing in ROOT v5/25/02)
260 // this function is a mess - and it should have been Fons who should have gone through the pain of writing it! (JF)
263 // construct variable size arrays for fixed size binning axes because TH3 lacks some constructors
264 Double_t* xBins = new Double_t[hist->GetNbinsX()+1];
265 Double_t* zBins = new Double_t[hist->GetNbinsZ()+1];
267 for (Int_t x=1; x<=hist->GetNbinsX()+1; x++)
269 xBins[x-1] = hist->GetXaxis()->GetBinLowEdge(x);
270 //Printf("%d %f", x, xBins[x-1]);
273 for (Int_t z=1; z<=hist->GetNbinsZ()+1; z++)
275 zBins[z-1] = hist->GetZaxis()->GetBinLowEdge(z);
276 //Printf("%d %f", y, yBins[y-1]);
279 TH3F* temp = new TH3F(hist->GetName(), Form("%s;%s;%s;%s", hist->GetTitle(), hist->GetXaxis()->GetTitle(), hist->GetYaxis()->GetTitle(), hist->GetZaxis()->GetTitle()), hist->GetNbinsX(), xBins, nBins, newBins, hist->GetNbinsZ(), zBins);
281 for (Int_t x=0; x<=hist->GetNbinsX()+1; x++)
282 for (Int_t y=0; y<=hist->GetNbinsY()+1; y++)
283 for (Int_t z=0; z<=hist->GetNbinsZ()+1; z++)
284 temp->Fill(hist->GetXaxis()->GetBinCenter(x), hist->GetYaxis()->GetBinCenter(y), hist->GetZaxis()->GetBinCenter(z), hist->GetBinContent(x, y, z));
286 for (Int_t x=0; x<=temp->GetNbinsX()+1; x++)
287 for (Int_t y=0; y<=temp->GetNbinsY()+1; y++)
288 for (Int_t z=0; z<=hist->GetNbinsZ()+1; z++)
289 temp->SetBinError(x, y, z, TMath::Sqrt(temp->GetBinContent(x, y, z)));
298 //____________________________________________________________________
299 void AliMultiplicityCorrection::RebinGenerated(Int_t nBins05, Double_t* newBins05, Int_t nBins10, Double_t* newBins10, Int_t nBins13, Double_t* newBins13)
302 // Rebins the (and only the) generated multiplicity axis
305 Printf("Rebinning generated-multiplicity axis...");
307 // do not add this hists to the directory
308 Bool_t oldStatus = TH1::AddDirectoryStatus();
309 TH1::AddDirectory(kFALSE);
312 AliFatal("This function only works for three ESD hists!");
314 for (Int_t i = 0; i < kESDHists; ++i)
317 Double_t* newBins = 0;
337 fNoVertexEvents[i] = (TH1F*) fNoVertexEvents[i]->Rebin(nBins, fNoVertexEvents[i]->GetName(), newBins);
338 fMultiplicityESDCorrected[i] = (TH1F*) fMultiplicityESDCorrected[i]->Rebin(nBins, fMultiplicityESDCorrected[i]->GetName(), newBins);
341 Rebin2DY(fMultiplicityVtx[i], nBins, newBins);
342 Rebin2DY(fMultiplicityMB[i], nBins, newBins);
343 Rebin2DY(fMultiplicityINEL[i], nBins, newBins);
344 Rebin2DY(fMultiplicityNSD[i], nBins, newBins);
347 Rebin3DY(fCorrelation[i], nBins, newBins);
350 TH1::AddDirectory(oldStatus);
353 //____________________________________________________________________
354 AliMultiplicityCorrection::~AliMultiplicityCorrection()
360 Printf("AliMultiplicityCorrection::~AliMultiplicityCorrection called");
362 for (Int_t i = 0; i < kESDHists; ++i)
364 if (fMultiplicityESD[i])
365 delete fMultiplicityESD[i];
366 fMultiplicityESD[i] = 0;
368 if (fTriggeredEvents[i])
369 delete fTriggeredEvents[i];
370 fTriggeredEvents[i]= 0;
372 if (fNoVertexEvents[i])
373 delete fNoVertexEvents[i];
374 fNoVertexEvents[i]= 0;
377 for (Int_t i = 0; i < kMCHists; ++i)
379 if (fMultiplicityVtx[i])
380 delete fMultiplicityVtx[i];
381 fMultiplicityVtx[i] = 0;
383 if (fMultiplicityMB[i])
384 delete fMultiplicityMB[i];
385 fMultiplicityMB[i] = 0;
387 if (fMultiplicityINEL[i])
388 delete fMultiplicityINEL[i];
389 fMultiplicityINEL[i] = 0;
391 if (fMultiplicityNSD[i])
392 delete fMultiplicityNSD[i];
393 fMultiplicityNSD[i] = 0;
396 for (Int_t i = 0; i < kCorrHists; ++i)
399 delete fCorrelation[i];
402 if (fMultiplicityESDCorrected[i])
403 delete fMultiplicityESDCorrected[i];
404 fMultiplicityESDCorrected[i] = 0;
408 //____________________________________________________________________
409 Long64_t AliMultiplicityCorrection::Merge(const TCollection* list)
411 // Merge a list of AliMultiplicityCorrection objects with this (needed for
413 // Returns the number of merged objects (including this).
421 TIterator* iter = list->MakeIterator();
424 // collections of all histograms
425 TList collections[3*kESDHists+kMCHists*4+kCorrHists*2];
428 while ((obj = iter->Next())) {
430 AliMultiplicityCorrection* entry = dynamic_cast<AliMultiplicityCorrection*> (obj);
434 for (Int_t i = 0; i < kESDHists; ++i)
436 collections[i].Add(entry->fMultiplicityESD[i]);
437 collections[kESDHists+i].Add(entry->fTriggeredEvents[i]);
438 collections[kESDHists*2+i].Add(entry->fNoVertexEvents[i]);
441 for (Int_t i = 0; i < kMCHists; ++i)
443 collections[3*kESDHists+i].Add(entry->fMultiplicityVtx[i]);
444 collections[3*kESDHists+kMCHists+i].Add(entry->fMultiplicityMB[i]);
445 collections[3*kESDHists+kMCHists*2+i].Add(entry->fMultiplicityINEL[i]);
446 collections[3*kESDHists+kMCHists*3+i].Add(entry->fMultiplicityNSD[i]);
449 for (Int_t i = 0; i < kCorrHists; ++i)
450 collections[3*kESDHists+kMCHists*4+i].Add(entry->fCorrelation[i]);
452 for (Int_t i = 0; i < kCorrHists; ++i)
453 collections[3*kESDHists+kMCHists*4+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]);
458 for (Int_t i = 0; i < kESDHists; ++i)
460 fMultiplicityESD[i]->Merge(&collections[i]);
461 fTriggeredEvents[i]->Merge(&collections[kESDHists+i]);
462 fNoVertexEvents[i]->Merge(&collections[2*kESDHists+i]);
465 for (Int_t i = 0; i < kMCHists; ++i)
467 fMultiplicityVtx[i]->Merge(&collections[3*kESDHists+i]);
468 fMultiplicityMB[i]->Merge(&collections[3*kESDHists+kMCHists+i]);
469 fMultiplicityINEL[i]->Merge(&collections[3*kESDHists+kMCHists*2+i]);
470 fMultiplicityNSD[i]->Merge(&collections[3*kESDHists+kMCHists*3+i]);
473 for (Int_t i = 0; i < kCorrHists; ++i)
474 fCorrelation[i]->Merge(&collections[3*kESDHists+kMCHists*4+i]);
476 for (Int_t i = 0; i < kCorrHists; ++i)
477 fMultiplicityESDCorrected[i]->Merge(&collections[3*kESDHists+kMCHists*4+kCorrHists+i]);
484 //____________________________________________________________________
485 Bool_t AliMultiplicityCorrection::LoadHistograms(const Char_t* dir)
488 // loads the histograms from a file
489 // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
495 if (!gDirectory->cd(dir))
498 // store old hists to delete them later
500 oldObjects.SetOwner(1);
501 for (Int_t i = 0; i < kESDHists; ++i)
503 if (fMultiplicityESD[i])
504 oldObjects.Add(fMultiplicityESD[i]);
505 if (fTriggeredEvents[i])
506 oldObjects.Add(fTriggeredEvents[i]);
507 if (fNoVertexEvents[i])
508 oldObjects.Add(fNoVertexEvents[i]);
511 for (Int_t i = 0; i < kMCHists; ++i)
513 if (fMultiplicityVtx[i])
514 oldObjects.Add(fMultiplicityVtx[i]);
515 if (fMultiplicityMB[i])
516 oldObjects.Add(fMultiplicityMB[i]);
517 if (fMultiplicityINEL[i])
518 oldObjects.Add(fMultiplicityINEL[i]);
519 if (fMultiplicityNSD[i])
520 oldObjects.Add(fMultiplicityNSD[i]);
523 for (Int_t i = 0; i < kCorrHists; ++i)
525 oldObjects.Add(fCorrelation[i]);
529 Bool_t success = kTRUE;
531 for (Int_t i = 0; i < kESDHists; ++i)
533 fMultiplicityESD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityESD[i]->GetName()));
534 fTriggeredEvents[i] = dynamic_cast<TH1F*> (gDirectory->Get(fTriggeredEvents[i]->GetName()));
535 fNoVertexEvents[i] = dynamic_cast<TH1F*> (gDirectory->Get(fNoVertexEvents[i]->GetName()));
536 if (!fMultiplicityESD[i] || !fTriggeredEvents[i] || !fNoVertexEvents[i])
540 for (Int_t i = 0; i < kMCHists; ++i)
542 fMultiplicityVtx[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityVtx[i]->GetName()));
543 fMultiplicityMB[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityMB[i]->GetName()));
544 fMultiplicityINEL[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityINEL[i]->GetName()));
545 fMultiplicityNSD[i] = dynamic_cast<TH2F*> (gDirectory->Get(fMultiplicityNSD[i]->GetName()));
546 if (!fMultiplicityVtx[i] || !fMultiplicityMB[i] || !fMultiplicityINEL[i])
550 for (Int_t i = 0; i < kCorrHists; ++i)
552 fCorrelation[i] = dynamic_cast<TH3F*> (gDirectory->Get(fCorrelation[i]->GetName()));
553 if (!fCorrelation[i])
555 fMultiplicityESDCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fMultiplicityESDCorrected[i]->GetName()));
556 if (!fMultiplicityESDCorrected[i])
560 gDirectory->cd("..");
568 //____________________________________________________________________
569 void AliMultiplicityCorrection::SaveHistograms(const char* dir)
572 // saves the histograms
578 gDirectory->mkdir(dir);
581 for (Int_t i = 0; i < kESDHists; ++i)
583 if (fMultiplicityESD[i])
585 fMultiplicityESD[i]->Write();
586 fMultiplicityESD[i]->ProjectionY(Form("%s_px", fMultiplicityESD[i]->GetName()), 1, fMultiplicityESD[i]->GetNbinsX())->Write();
588 if (fTriggeredEvents[i])
589 fTriggeredEvents[i]->Write();
590 if (fNoVertexEvents[i])
591 fNoVertexEvents[i]->Write();
594 for (Int_t i = 0; i < kMCHists; ++i)
596 if (fMultiplicityVtx[i])
598 fMultiplicityVtx[i]->Write();
599 fMultiplicityVtx[i]->ProjectionY(Form("%s_px", fMultiplicityVtx[i]->GetName()), 1, fMultiplicityVtx[i]->GetNbinsX())->Write();
601 if (fMultiplicityMB[i])
603 fMultiplicityMB[i]->Write();
604 fMultiplicityMB[i]->ProjectionY(Form("%s_px", fMultiplicityMB[i]->GetName()), 1, fMultiplicityMB[i]->GetNbinsX())->Write();
606 if (fMultiplicityINEL[i])
608 fMultiplicityINEL[i]->Write();
609 fMultiplicityINEL[i]->ProjectionY(Form("%s_px", fMultiplicityINEL[i]->GetName()), 1, fMultiplicityINEL[i]->GetNbinsX())->Write();
611 if (fMultiplicityNSD[i])
613 fMultiplicityNSD[i]->Write();
614 fMultiplicityNSD[i]->ProjectionY(Form("%s_px", fMultiplicityNSD[i]->GetName()), 1, fMultiplicityNSD[i]->GetNbinsX())->Write();
618 for (Int_t i = 0; i < kCorrHists; ++i)
621 fCorrelation[i]->Write();
622 if (fMultiplicityESDCorrected[i])
623 fMultiplicityESDCorrected[i]->Write();
626 gDirectory->cd("..");
629 //____________________________________________________________________
630 void 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)
633 // Fills an event from MC
638 fMultiplicityMB[0]->Fill(vtx, generated05);
639 fMultiplicityMB[1]->Fill(vtx, generated10);
640 fMultiplicityMB[2]->Fill(vtx, generated14);
641 fMultiplicityMB[3]->Fill(vtx, generatedAll);
645 fMultiplicityVtx[0]->Fill(vtx, generated05);
646 fMultiplicityVtx[1]->Fill(vtx, generated10);
647 fMultiplicityVtx[2]->Fill(vtx, generated14);
648 fMultiplicityVtx[3]->Fill(vtx, generatedAll);
652 fMultiplicityINEL[0]->Fill(vtx, generated05);
653 fMultiplicityINEL[1]->Fill(vtx, generated10);
654 fMultiplicityINEL[2]->Fill(vtx, generated14);
655 fMultiplicityINEL[3]->Fill(vtx, generatedAll);
657 if (processType != AliPWG0Helper::kSD)
659 fMultiplicityNSD[0]->Fill(vtx, generated05);
660 fMultiplicityNSD[1]->Fill(vtx, generated10);
661 fMultiplicityNSD[2]->Fill(vtx, generated14);
662 fMultiplicityNSD[3]->Fill(vtx, generatedAll);
666 //____________________________________________________________________
667 void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured14)
670 // Fills an event from ESD
673 fMultiplicityESD[0]->Fill(vtx, measured05);
674 fMultiplicityESD[1]->Fill(vtx, measured10);
675 fMultiplicityESD[2]->Fill(vtx, measured14);
678 //____________________________________________________________________
679 void AliMultiplicityCorrection::FillTriggeredEvent(Int_t measured05, Int_t measured10, Int_t measured14)
682 // fills raw distribution of triggered events
685 fTriggeredEvents[0]->Fill(measured05);
686 fTriggeredEvents[1]->Fill(measured10);
687 fTriggeredEvents[2]->Fill(measured14);
690 //____________________________________________________________________
691 void AliMultiplicityCorrection::FillNoVertexEvent(Float_t vtx, Bool_t vertexReconstructed, Int_t generated05, Int_t generated10, Int_t generated14, Int_t measured05, Int_t measured10, Int_t measured14)
694 // fills raw distribution of triggered events
697 if (vtx > fgVtxRangeBegin[0] && vtx < fgVtxRangeEnd[0] && (!vertexReconstructed || measured05 == 0))
698 fNoVertexEvents[0]->Fill(generated05);
700 if (vtx > fgVtxRangeBegin[1] && vtx < fgVtxRangeEnd[1] && (!vertexReconstructed || measured10 == 0))
701 fNoVertexEvents[1]->Fill(generated10);
703 if (vtx > fgVtxRangeBegin[2] && vtx < fgVtxRangeEnd[2] && (!vertexReconstructed || measured14 == 0))
704 fNoVertexEvents[2]->Fill(generated14);
707 //____________________________________________________________________
708 void 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)
711 // Fills an event into the correlation map with the information from MC and ESD
714 fCorrelation[0]->Fill(vtx, generated05, measured05);
715 fCorrelation[1]->Fill(vtx, generated10, measured10);
716 fCorrelation[2]->Fill(vtx, generated14, measured14);
718 fCorrelation[3]->Fill(vtx, generatedAll, measured05);
719 fCorrelation[4]->Fill(vtx, generatedAll, measured10);
720 fCorrelation[5]->Fill(vtx, generatedAll, measured14);
723 //____________________________________________________________________
724 void AliMultiplicityCorrection::SetupCurrentHists(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
727 // fills fCurrentESD, fCurrentCorrelation
728 // resets fMultiplicityESDCorrected
731 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
733 fMultiplicityESDCorrected[correlationID]->Reset();
734 fMultiplicityESDCorrected[correlationID]->Sumw2();
736 // project without under/overflow bins
738 Int_t end = fMultiplicityESD[inputRange]->GetXaxis()->GetNbins();
739 if (fVtxEnd > fVtxBegin)
744 fCurrentESD = fMultiplicityESD[inputRange]->ProjectionY("fCurrentESD", begin, end);
745 fCurrentESD->Sumw2();
747 // empty under/overflow bins in x, otherwise Project3D takes them into account
748 TH3* hist = fCorrelation[correlationID];
749 for (Int_t y=0; y<=hist->GetYaxis()->GetNbins()+1; ++y)
751 for (Int_t z=0; z<=hist->GetZaxis()->GetNbins()+1; ++z)
753 hist->SetBinContent(0, y, z, 0);
754 hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
758 if (fVtxEnd > fVtxBegin)
759 hist->GetXaxis()->SetRange(fVtxBegin, fVtxEnd);
761 fCurrentCorrelation = (TH2*) hist->Project3D("zy");
762 fCurrentCorrelation->Sumw2();
764 Printf("AliMultiplicityCorrection::SetupCurrentHists: Statistics information: %.f entries in correlation map; %.f entries in measured spectrum", fCurrentCorrelation->Integral(), fCurrentESD->Integral());
766 #if 0 // does not help
767 // null bins with one entry
768 Int_t nNulledBins = 0;
769 for (Int_t x=1; x<=fCurrentCorrelation->GetXaxis()->GetNbins(); ++x)
770 for (Int_t y=1; y<=fCurrentCorrelation->GetYaxis()->GetNbins(); ++y)
772 if (fCurrentCorrelation->GetBinContent(x, y) == 1)
774 fCurrentCorrelation->SetBinContent(x, y, 0);
775 fCurrentCorrelation->SetBinError(x, y, 0);
780 Printf("Nulled %d bins", nNulledBins);
783 fCurrentEfficiency = GetEfficiency(inputRange, eventType);
784 //fCurrentEfficiency->Rebin(2);
785 //fCurrentEfficiency->Scale(0.5);
788 //____________________________________________________________________
789 TH1* AliMultiplicityCorrection::GetEfficiency(Int_t inputRange, EventType eventType)
792 // calculates efficiency for given event type
796 name1.Form("divisor%d", inputRange);
799 name2.Form("CurrentEfficiency%d", inputRange);
805 case kMB: divisor = fMultiplicityMB[inputRange]->ProjectionY(name1, 1, fMultiplicityMB[inputRange]->GetNbinsX(), "e"); break;
806 case kINEL: divisor = fMultiplicityINEL[inputRange]->ProjectionY(name1, 1, fMultiplicityINEL[inputRange]->GetNbinsX(), "e"); break;
807 case kNSD: divisor = fMultiplicityNSD[inputRange]->ProjectionY(name1, 1, fMultiplicityNSD[inputRange]->GetNbinsX(), "e"); break;
809 TH1* eff = fMultiplicityVtx[inputRange]->ProjectionY(name2, 1, fMultiplicityVtx[inputRange]->GetNbinsX(), "e");
811 if (eventType == kTrVtx)
813 for (Int_t i=0; i<= eff->GetNbinsX()+1; i++)
814 eff->SetBinContent(i, 1);
817 eff->Divide(eff, divisor, 1, 1, "B");
822 //____________________________________________________________________
823 TH1* AliMultiplicityCorrection::GetTriggerEfficiency(Int_t inputRange, EventType eventType)
826 // calculates efficiency for given event type
830 name1.Form("divisor%d", inputRange);
833 name2.Form("CurrentEfficiency%d", inputRange);
838 case kTrVtx : AliFatal("Not supported!"); break;
839 case kMB: divisor = fMultiplicityMB[inputRange]->ProjectionY (name1, 1, fMultiplicityMB[inputRange]->GetNbinsX(), "e"); break;
840 case kINEL: divisor = fMultiplicityINEL[inputRange]->ProjectionY(name1, 1, fMultiplicityINEL[inputRange]->GetNbinsX(), "e"); break;
841 case kNSD: divisor = fMultiplicityNSD[inputRange]->ProjectionY (name1, 1, fMultiplicityNSD[inputRange]->GetNbinsX(), "e"); break;
843 TH1* eff = fMultiplicityMB[inputRange]->ProjectionY(name2, 1, fMultiplicityMB[inputRange]->GetNbinsX(), "e");
845 eff->Divide(eff, divisor, 1, 1, "B");
849 //____________________________________________________________________
850 Int_t AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Int_t zeroBinEvents, Bool_t check, TH1* initialConditions, Bool_t errorAsBias)
853 // correct spectrum using minuit chi2 method
855 // for description of parameters, see AliUnfolding::Unfold
858 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
860 //AliUnfolding::SetCreateOverflowBin(5);
861 AliUnfolding::SetUnfoldingMethod(AliUnfolding::kChi2Minimization);
862 AliUnfolding::SetMinimumInitialValue(kTRUE, 0.1);
864 // use here only vtx efficiency (to MB sample) which is always needed if we use the 0 bin
865 SetupCurrentHists(inputRange, fullPhaseSpace, (eventType == kTrVtx) ? kTrVtx : kMB);
867 // TODO set errors on measured with 0.5 * TMath::ChisquareQuantile(0.1, 20000)
868 // see PDG: Statistics / Poission or binomial data / Eq. 32.49a/b in 2004 edition
870 Calculate0Bin(inputRange, eventType, zeroBinEvents);
872 Int_t resultCode = -1;
873 if (errorAsBias == kFALSE)
875 resultCode = AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID], check);
879 resultCode = AliUnfolding::UnfoldGetBias(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID]);
882 // HACK store new vertex reco efficiency for bin 0, changing number of events with trigger and vertex in MC map
883 if (zeroBinEvents > 0)
885 Printf("WARNING: Stored vertex reco efficiency from unfolding for bin 0.");
886 fMultiplicityVtx[inputRange]->SetBinContent(1, 1, fMultiplicityMB[inputRange]->GetBinContent(1, 1) * fCurrentEfficiency->GetBinContent(1));
889 // correct for the trigger bias if requested
892 Printf("Applying trigger efficiency");
893 TH1* eff = GetTriggerEfficiency(inputRange, eventType);
894 for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); i++)
896 fMultiplicityESDCorrected[correlationID]->SetBinContent(i, fMultiplicityESDCorrected[correlationID]->GetBinContent(i) / eff->GetBinContent(i));
897 fMultiplicityESDCorrected[correlationID]->SetBinError(i, fMultiplicityESDCorrected[correlationID]->GetBinError(i) / eff->GetBinContent(i));
904 //____________________________________________________________________
905 void AliMultiplicityCorrection::Calculate0Bin(Int_t inputRange, EventType eventType, Int_t zeroBinEvents)
909 if (eventType == kTrVtx)
912 Double_t fractionEventsInVertexRange = fMultiplicityESD[inputRange]->Integral(1, fMultiplicityESD[inputRange]->GetXaxis()->GetNbins()) / fMultiplicityESD[inputRange]->Integral(0, fMultiplicityESD[inputRange]->GetXaxis()->GetNbins() + 1);
914 // difference of fraction that is inside the considered range between triggered events and events with vertex
915 // Extension to NSD not needed, INEL and NSD vertex distributions are nature-given and unbiased!
916 Double_t differenceVtxDist = (fMultiplicityINEL[inputRange]->Integral(1, fMultiplicityINEL[inputRange]->GetXaxis()->GetNbins()) / fMultiplicityINEL[inputRange]->Integral(0, fMultiplicityINEL[inputRange]->GetXaxis()->GetNbins() + 1)) / (fMultiplicityVtx[inputRange]->Integral(1, fMultiplicityVtx[inputRange]->GetXaxis()->GetNbins()) / fMultiplicityVtx[inputRange]->Integral(0, fMultiplicityVtx[inputRange]->GetXaxis()->GetNbins() + 1));
918 Printf("Enabling 0 bin estimate for triggered events without vertex.");
919 Printf(" Events in 0 bin: %d", zeroBinEvents);
920 Printf(" Fraction in range: %.1f%%", fractionEventsInVertexRange * 100);
921 Printf(" Difference Vtx Dist: %f", differenceVtxDist);
923 AliUnfolding::SetNotFoundEvents(zeroBinEvents * fractionEventsInVertexRange * differenceVtxDist);
926 //____________________________________________________________________
927 void AliMultiplicityCorrection::FixTriggerEfficiencies(Int_t start)
930 // sets trigger and vertex efficiencies to 1 for large multiplicities where no event was simulated
933 for (Int_t etaRange = 0; etaRange < kMCHists; etaRange++)
935 for (Int_t i = 1; i <= fMultiplicityINEL[etaRange]->GetNbinsY(); i++)
937 if (fMultiplicityINEL[etaRange]->GetYaxis()->GetBinCenter(i) < start)
940 if (fMultiplicityINEL[etaRange]->GetBinContent(1, i) > 0)
943 fMultiplicityINEL[etaRange]->SetBinContent(1, i, 1);
944 fMultiplicityNSD[etaRange]->SetBinContent(1, i, 1);
945 fMultiplicityMB[etaRange]->SetBinContent(1, i, 1);
946 fMultiplicityVtx[etaRange]->SetBinContent(1, i, 1);
951 //____________________________________________________________________
952 Float_t AliMultiplicityCorrection::GetFraction0Generated(Int_t inputRange)
955 // returns the fraction of events that have 0 generated particles in the given range among all events without vertex OR 0 tracklets
958 TH1* multMB = GetNoVertexEvents(inputRange);
959 return multMB->GetBinContent(1) / multMB->Integral();
962 //____________________________________________________________________
963 Int_t AliMultiplicityCorrection::ApplyNBDFit(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
966 // correct spectrum using minuit chi2 method with a NBD function
968 // for description of parameters, see AliUnfolding::Unfold
971 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
973 AliUnfolding::SetUnfoldingMethod(AliUnfolding::kFunction);
974 SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
976 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])");
977 func->SetParNames("scaling", "averagen", "k");
978 func->SetParLimits(0, 0, 1000);
979 func->SetParLimits(1, 1, 50);
980 func->SetParLimits(2, 1, 10);
981 func->SetParameters(1, 10, 2);
982 AliUnfolding::SetFunction(func);
984 return AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, 0, fMultiplicityESDCorrected[correlationID]);
987 //____________________________________________________________________
988 void AliMultiplicityCorrection::DrawHistograms()
991 // draws the histograms of this class
996 TCanvas* canvas1 = new TCanvas("fMultiplicityESD", "fMultiplicityESD", 900, 600);
997 canvas1->Divide(3, 2);
998 for (Int_t i = 0; i < kESDHists; ++i)
1001 fMultiplicityESD[i]->DrawCopy("COLZ");
1002 printf("%d --> %f\n", i, (Float_t) fMultiplicityESD[i]->ProjectionY()->GetMean());
1007 TCanvas* canvas2 = new TCanvas("fMultiplicityMC", "fMultiplicityMC", 900, 600);
1008 canvas2->Divide(3, 2);
1009 for (Int_t i = 0; i < kMCHists; ++i)
1012 fMultiplicityVtx[i]->DrawCopy("COLZ");
1013 printf("%d --> %f\n", i, fMultiplicityVtx[i]->ProjectionY()->GetMean());
1016 TCanvas* canvas3 = new TCanvas("fCorrelation", "fCorrelation", 900, 900);
1017 canvas3->Divide(3, 3);
1018 for (Int_t i = 0; i < kCorrHists; ++i)
1021 TH3* hist = static_cast<TH3*> (fCorrelation[i]->Clone());
1022 for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y)
1024 for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z)
1026 hist->SetBinContent(0, y, z, 0);
1027 hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
1030 TH1* proj = hist->Project3D("zy");
1031 proj->DrawCopy("COLZ");
1035 //____________________________________________________________________
1036 void AliMultiplicityCorrection::DrawComparison(const char* name, Int_t inputRange, Bool_t fullPhaseSpace, Bool_t /*normalizeESD*/, TH1* mcHist, Bool_t simple, EventType eventType)
1038 // draw comparison plots
1040 Int_t esdCorrId = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1043 tmpStr.Form("%s_DrawComparison_%d", name, esdCorrId);
1045 if (fMultiplicityESDCorrected[esdCorrId]->Integral() == 0)
1047 printf("ERROR. Unfolded histogram is empty\n");
1052 Int_t end = fMultiplicityESD[inputRange]->GetXaxis()->GetNbins();
1053 if (fVtxEnd > fVtxBegin)
1058 fCurrentESD = fMultiplicityESD[esdCorrId]->ProjectionY("fCurrentESD", begin, end);
1059 fCurrentESD->Sumw2();
1063 for (Int_t i=5; i<=mcHist->GetNbinsX(); ++i)
1065 if (mcHist->GetBinContent(i) > 0)
1066 mcMax = (Int_t) mcHist->GetXaxis()->GetBinCenter(i) + 2;
1070 for (Int_t i=5; i<=fMultiplicityESDCorrected[esdCorrId]->GetNbinsX(); ++i)
1071 if (fMultiplicityESDCorrected[esdCorrId]->GetBinContent(i) > 1)
1072 mcMax = (Int_t) fMultiplicityESDCorrected[esdCorrId]->GetXaxis()->GetBinCenter(i) + 2;
1074 Printf("AliMultiplicityCorrection::DrawComparison: MC bin limit is %d", mcMax);
1075 // calculate residual
1077 TH1* convolutedProj = (TH1*) GetConvoluted(esdCorrId, eventType)->Clone("convolutedProj");
1078 TH1* residual = GetResiduals(esdCorrId, eventType, tmp);
1080 TH1* residualHist = new TH1F("residualHist", "residualHist", 51, -5, 5);
1083 for (Int_t i=1; i<=TMath::Min(residual->GetNbinsX(), 75); ++i)
1085 Float_t value = residual->GetBinContent(i);
1086 // TODO has to get a parameter (used in Chi2Function, GetResiduals, and here)
1088 chi2 += value * value;
1089 Printf("%d --> %f (%f)", i, value * value, chi2);
1090 residualHist->Fill(value);
1091 convolutedProj->SetBinError(i, 0);
1093 fLastChi2Residuals = chi2;
1095 //new TCanvas; residualHist->DrawCopy();
1097 printf("Difference (Residuals) is %f\n", fLastChi2Residuals);
1099 TCanvas* canvas1 = 0;
1102 canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 600);
1103 canvas1->Divide(2, 1);
1107 canvas1 = new TCanvas(tmpStr, tmpStr, 1200, 1200);
1108 canvas1->Divide(2, 3);
1112 canvas1->cd(1)->SetGridx();
1113 canvas1->cd(1)->SetGridy();
1114 canvas1->cd(1)->SetTopMargin(0.05);
1115 canvas1->cd(1)->SetRightMargin(0.05);
1116 canvas1->cd(1)->SetLeftMargin(0.12);
1117 canvas1->cd(1)->SetBottomMargin(0.12);
1118 TH1* proj = (TH1*) mcHist->Clone("proj");
1119 if (proj->GetEntries() > 0)
1120 AliPWG0Helper::NormalizeToBinWidth(proj);
1122 proj->GetXaxis()->SetRangeUser(0, mcMax);
1123 proj->GetYaxis()->SetTitleOffset(1.4);
1124 proj->SetTitle(Form(";True multiplicity in |#eta| < %.1f;Entries", (inputRange+1)*0.5));
1125 proj->SetStats(kFALSE);
1127 fMultiplicityESDCorrected[esdCorrId]->GetXaxis()->SetRangeUser(0, mcMax);
1128 fMultiplicityESDCorrected[esdCorrId]->GetYaxis()->SetRangeUser(0.1, fMultiplicityESDCorrected[esdCorrId]->GetMaximum() * 1.5);
1129 fMultiplicityESDCorrected[esdCorrId]->GetYaxis()->SetTitleOffset(1.4);
1130 fMultiplicityESDCorrected[esdCorrId]->SetTitle(Form(";True multiplicity in |#eta| < %.1f;Entries", (inputRange+1)*0.5));
1131 fMultiplicityESDCorrected[esdCorrId]->SetStats(kFALSE);
1133 fMultiplicityESDCorrected[esdCorrId]->SetLineColor(2);
1134 fMultiplicityESDCorrected[esdCorrId]->SetMarkerColor(2);
1136 TH1* esdCorrected = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("esdCorrected");
1137 AliPWG0Helper::NormalizeToBinWidth(esdCorrected);
1139 if (proj->GetEntries() > 0) {
1140 proj->DrawCopy("HIST");
1141 esdCorrected->DrawCopy("SAME HIST E");
1144 esdCorrected->DrawCopy("HIST E");
1148 TLegend* legend = new TLegend(0.3, 0.8, 0.93, 0.93);
1149 legend->AddEntry(proj, "True distribution");
1150 legend->AddEntry(fMultiplicityESDCorrected[esdCorrId], "Unfolded distribution");
1151 legend->SetFillColor(0);
1152 legend->SetTextSize(0.04);
1156 canvas1->cd(2)->SetGridx();
1157 canvas1->cd(2)->SetGridy();
1158 canvas1->cd(2)->SetTopMargin(0.05);
1159 canvas1->cd(2)->SetRightMargin(0.05);
1160 canvas1->cd(2)->SetLeftMargin(0.12);
1161 canvas1->cd(2)->SetBottomMargin(0.12);
1164 fCurrentESD->GetXaxis()->SetRangeUser(0, mcMax);
1165 fCurrentESD->SetTitle(Form(";Measured multiplicity in |#eta| < %.1f;Entries", (inputRange+1)*0.5));
1166 fCurrentESD->SetStats(kFALSE);
1167 fCurrentESD->GetYaxis()->SetTitleOffset(1.4);
1168 fCurrentESD->DrawCopy("HIST E");
1170 convolutedProj->SetLineColor(2);
1171 convolutedProj->SetMarkerColor(2);
1172 convolutedProj->SetMarkerStyle(5);
1173 convolutedProj->DrawCopy("HIST SAME P");
1175 legend = new TLegend(0.3, 0.8, 0.93, 0.93);
1176 legend->AddEntry(fCurrentESD, "Measured distribution");
1177 legend->AddEntry(convolutedProj, "R #otimes unfolded distribution", "P");
1178 legend->SetFillColor(0);
1179 legend->SetTextSize(0.04);
1185 residual->GetYaxis()->SetRangeUser(-5, 5);
1186 residual->GetXaxis()->SetRangeUser(0, mcMax);
1187 residual->SetStats(kFALSE);
1188 residual->DrawCopy();
1191 TH1* ratio = (TH1*) fMultiplicityESDCorrected[esdCorrId]->Clone("ratio");
1192 ratio->Divide(mcHist);
1193 ratio->SetTitle("Ratio;true multiplicity;Unfolded / MC");
1194 ratio->GetYaxis()->SetRangeUser(0.5, 1.5);
1195 ratio->GetXaxis()->SetRangeUser(0, mcMax);
1196 ratio->SetStats(kFALSE);
1197 ratio->Draw("HIST");
1199 // plot (MC - Unfolded) / error (MC)
1202 TH1* diffMCUnfolded2 = static_cast<TH1*> (proj->Clone("diffMCUnfolded2"));
1203 diffMCUnfolded2->Add(esdCorrected, -1);
1205 Int_t ndfQual[kQualityRegions];
1206 for (Int_t region=0; region<kQualityRegions; ++region)
1208 fQuality[region] = 0;
1209 ndfQual[region] = 0;
1212 Double_t newChi2 = 0;
1213 Double_t newChi2Limit150 = 0;
1215 for (Int_t i=1; i<=diffMCUnfolded2->GetNbinsX(); ++i)
1218 if (proj->GetBinError(i) > 0)
1220 value = diffMCUnfolded2->GetBinContent(i) / proj->GetBinError(i);
1221 newChi2 += value * value;
1222 if (i > 1 && i <= mcMax)
1223 newChi2Limit150 += value * value;
1226 for (Int_t region=0; region<kQualityRegions; ++region)
1227 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
1229 fQuality[region] += TMath::Abs(value);
1234 diffMCUnfolded2->SetBinContent(i, value);
1237 // normalize region to the number of entries
1238 for (Int_t region=0; region<kQualityRegions; ++region)
1240 if (ndfQual[region] > 0)
1241 fQuality[region] /= ndfQual[region];
1242 Printf("Quality parameter %d (%d <= mult <= %d) is %f with %d df", region, fgQualityRegionsB[region], fgQualityRegionsE[region], fQuality[region], ndfQual[region]);
1247 fLastChi2MC = newChi2Limit150 / (mcMax - 1);
1248 Printf("Chi2 (2..%d) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", mcMax, newChi2Limit150, mcMax - 1, fLastChi2MC);
1253 Printf("Chi2 (full range) from (MC - Unfolded) / e(MC) is: %.2f ndf is %d --> chi2 / ndf = %.2f", newChi2, ndf, ((ndf > 0) ? newChi2 / ndf : -1));
1255 diffMCUnfolded2->SetTitle("#chi^{2};true multiplicity;(MC - Unfolded) / e(MC)");
1256 diffMCUnfolded2->GetYaxis()->SetRangeUser(-5, 5);
1257 diffMCUnfolded2->GetXaxis()->SetRangeUser(0, mcMax);
1258 diffMCUnfolded2->DrawCopy("HIST");
1261 // draw penalty factor
1263 TH1* penalty = AliUnfolding::GetPenaltyPlot(fMultiplicityESDCorrected[esdCorrId]);
1264 penalty->SetStats(0);
1265 penalty->GetXaxis()->SetRangeUser(0, mcMax);
1266 penalty->DrawCopy("HIST");
1270 //____________________________________________________________________
1271 void AliMultiplicityCorrection::FFT(Int_t dir, Int_t m, Double_t *x, Double_t *y) const
1273 /*-------------------------------------------------------------------------
1274 This computes an in-place complex-to-complex FFT
1275 x and y are the real and imaginary arrays of 2^m points.
1276 dir = 1 gives forward transform
1277 dir = -1 gives reverse transform
1282 1 \ - j k 2 pi n / N
1283 X(n) = --- > x(k) e = forward transform
1292 X(n) = > x(k) e = forward transform
1298 Long_t nn, i, i1, j, k, i2, l, l1, l2;
1299 Double_t c1, c2, tx, ty, t1, t2, u1, u2, z;
1301 /* Calculate the number of points */
1303 for (i = 0; i < m; i++)
1306 /* Do the bit reversal */
1309 for (i= 0; i < nn - 1; i++) {
1326 /* Compute the FFT */
1330 for (l = 0; l < m; l++) {
1335 for (j = 0;j < l1; j++) {
1336 for (i = j; i < nn; i += l2) {
1338 t1 = u1 * x[i1] - u2 * y[i1];
1339 t2 = u1 * y[i1] + u2 * x[i1];
1345 z = u1 * c1 - u2 * c2;
1346 u2 = u1 * c2 + u2 * c1;
1349 c2 = TMath::Sqrt((1.0 - c1) / 2.0);
1352 c1 = TMath::Sqrt((1.0 + c1) / 2.0);
1355 /* Scaling for forward transform */
1357 for (i=0;i<nn;i++) {
1358 x[i] /= (Double_t)nn;
1359 y[i] /= (Double_t)nn;
1364 //____________________________________________________________________
1365 void AliMultiplicityCorrection::GetComparisonResults(Float_t* const mc, Int_t* const mcLimit, Float_t* const residuals, Float_t* const ratioAverage) const
1367 // Returns the chi2 between the MC and the unfolded ESD as well as between the ESD and the folded unfolded ESD
1368 // These values are computed during DrawComparison, thus this function picks up the
1374 *mcLimit = fLastChi2MCLimit;
1376 *residuals = fLastChi2Residuals;
1378 *ratioAverage = fRatioAverage;
1381 //____________________________________________________________________
1382 TH2F* AliMultiplicityCorrection::GetMultiplicityMC(Int_t i, EventType eventType) const
1385 // returns the corresponding MC spectrum
1390 case kTrVtx : return fMultiplicityVtx[i]; break;
1391 case kMB : return fMultiplicityMB[i]; break;
1392 case kINEL : return fMultiplicityINEL[i]; break;
1393 case kNSD : return fMultiplicityNSD[i]; break;
1399 //____________________________________________________________________
1400 void AliMultiplicityCorrection::SetMultiplicityMC(Int_t i, EventType eventType, TH2F* const hist)
1403 // returns the corresponding MC spectrum
1408 case kTrVtx : fMultiplicityVtx[i] = hist; break;
1409 case kMB : fMultiplicityMB[i] = hist; break;
1410 case kINEL : fMultiplicityINEL[i] = hist; break;
1411 case kNSD : fMultiplicityNSD[i] = hist; break;
1415 //____________________________________________________________________
1416 TH1* AliMultiplicityCorrection::CalculateStdDev(TH1** results, Int_t max)
1418 // calculate standard deviation of (results[0] - results[k]) k=1...max-1
1419 // per bin one gets: sigma(r[0] - r[n]) / r[0]
1421 TH1* standardDeviation = (TH1*) results[0]->Clone("standardDeviation");
1422 standardDeviation->Reset();
1424 for (Int_t x=1; x<=results[0]->GetNbinsX(); x++)
1426 if (results[0]->GetBinContent(x) > 0)
1428 Double_t average = 0;
1429 for (Int_t n=1; n<max; ++n)
1430 average += results[n]->GetBinContent(x) - results[0]->GetBinContent(x);
1433 Double_t variance = 0;
1434 for (Int_t n=1; n<max; ++n)
1436 Double_t value = results[n]->GetBinContent(x) - results[0]->GetBinContent(x) - average;
1437 variance += value * value;
1441 Double_t standardDev = TMath::Sqrt(variance);
1442 standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x));
1443 //Printf("sigma_%d is %f value %f --> error %f", x, standardDev, results[0]->GetBinContent(x), standardDev / results[0]->GetBinContent(x));
1447 return standardDeviation;
1450 //____________________________________________________________________
1451 TH1* AliMultiplicityCorrection::StatisticalUncertainty(AliUnfolding::MethodType methodType, Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Int_t zeroBinEvents, Bool_t randomizeMeasured, Bool_t randomizeResponse, const TH1* compareTo)
1454 // evaluates the uncertainty that arises from the non-infinite statistics in the response matrix
1455 // the function unfolds the spectrum using the default response matrix and several modified ones
1456 // the modified ones are created by randomizing each cell using poisson statistics with the mean = bin value
1457 // these unfolded results are compared to the first result gained with the default response OR to the histogram given
1458 // in <compareTo> (optional)
1460 // returns the error assigned to the measurement
1463 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1465 // initialize seed with current time
1466 gRandom->SetSeed(0);
1468 if (methodType == AliUnfolding::kChi2Minimization)
1470 Calculate0Bin(inputRange, eventType, zeroBinEvents);
1471 AliUnfolding::SetMinimumInitialValue(kTRUE, 0.1);
1474 AliUnfolding::SetUnfoldingMethod(methodType);
1476 const Int_t kErrorIterations = 20;
1479 TH1* firstResult = 0;
1481 TH1** results = new TH1*[kErrorIterations];
1483 for (Int_t n=0; n<kErrorIterations; ++n)
1485 Printf("Iteration %d of %d...", n, kErrorIterations);
1487 SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
1489 TH1* measured = (TH1*) fCurrentESD->Clone("measured");
1493 if (randomizeResponse)
1495 // randomize response matrix
1496 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1497 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1498 fCurrentCorrelation->SetBinContent(i, j, gRandom->Poisson(fCurrentCorrelation->GetBinContent(i, j)));
1501 if (randomizeMeasured)
1503 // randomize measured spectrum
1504 for (Int_t x=1; x<=measured->GetNbinsX(); x++) // mult. axis
1506 Int_t randomValue = gRandom->Poisson(fCurrentESD->GetBinContent(x));
1507 measured->SetBinContent(x, randomValue);
1508 measured->SetBinError(x, TMath::Sqrt(randomValue));
1513 // only for bayesian method we have to do it before the call to Unfold...
1514 if (methodType == AliUnfolding::kBayesian)
1516 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1518 // with this it is normalized to 1
1519 Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1521 // with this normalized to the given efficiency
1522 if (fCurrentEfficiency->GetBinContent(i) > 0)
1523 sum /= fCurrentEfficiency->GetBinContent(i);
1527 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1531 fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1532 fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1536 fCurrentCorrelation->SetBinContent(i, j, 0);
1537 fCurrentCorrelation->SetBinError(i, j, 0);
1544 if (n == 0 && compareTo)
1546 // in this case we just store the histogram we want to compare to
1547 result = (TH1*) compareTo->Clone("compareTo");
1552 result = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone(Form("result_%d", n));
1554 Int_t returnCode = AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, measured, 0, result);
1556 if (returnCode != 0)
1564 result->Scale(1.0 / result->Integral());
1568 firstResult = (TH1*) result->Clone("firstResult");
1570 maxError = (TH1*) result->Clone("maxError");
1576 TH1* ratio = (TH1*) firstResult->Clone("ratio");
1577 ratio->Divide(result);
1579 // find max. deviation
1580 for (Int_t x=1; x<=ratio->GetNbinsX(); x++)
1581 maxError->SetBinContent(x, TMath::Max(maxError->GetBinContent(x), TMath::Abs(1 - ratio->GetBinContent(x))));
1586 results[n] = result;
1589 // find covariance matrix
1590 // results[n] is X_x
1591 // cov. matrix is M_xy = E ( (X_x - E(X_x)) * (X_y - E(X_y))), with E() = expectation value
1593 Int_t nBins = results[0]->GetNbinsX();
1594 Float_t lowEdge = results[0]->GetXaxis()->GetBinLowEdge(1);
1595 Float_t upEdge = results[0]->GetXaxis()->GetBinUpEdge(nBins);
1597 // find average, E(X)
1598 TProfile* average = new TProfile("average", "average", nBins, lowEdge, upEdge);
1599 for (Int_t n=1; n<kErrorIterations; ++n)
1600 for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1601 average->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x));
1602 //new TCanvas; average->DrawClone();
1605 TProfile2D* covMatrix = new TProfile2D("covMatrix", "covMatrix", nBins, lowEdge, upEdge, nBins, lowEdge, upEdge);
1607 for (Int_t n=1; n<kErrorIterations; ++n)
1608 for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1609 for (Int_t y=1; y<=results[n]->GetNbinsX(); y++)
1611 // (X_x - E(X_x)) * (X_y - E(X_y)
1612 Float_t cov = (results[n]->GetBinContent(x) - average->GetBinContent(x)) * (results[n]->GetBinContent(y) - average->GetBinContent(y));
1613 covMatrix->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetXaxis()->GetBinCenter(y), cov);
1615 TCanvas* c = new TCanvas; c->cd(); covMatrix->DrawCopy("COLZ");
1617 // // fill 2D histogram that contains deviation from first
1618 // TH2F* deviations = new TH2F("deviations", "deviations", nBins, lowEdge, upEdge, 1000, -0.01, 0.01);
1619 // for (Int_t n=1; n<kErrorIterations; ++n)
1620 // for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1621 // deviations->Fill(results[n]->GetXaxis()->GetBinCenter(x), results[n]->GetBinContent(x) - results[0]->GetBinContent(x));
1622 // //new TCanvas; deviations->DrawCopy("COLZ");
1624 // // get standard deviation "by hand"
1625 // for (Int_t x=1; x<=nBins; x++)
1627 // if (results[0]->GetBinContent(x) > 0)
1629 // TH1* proj = deviations->ProjectionY("projRMS", x, x, "e");
1630 // Float_t standardDev = proj->GetRMS(); // this is standard deviation in fact
1631 // //standardDeviation->SetBinContent(x, standardDev / results[0]->GetBinContent(x));
1632 // Printf("sigma_%d is %f value %f --> error %f", x, standardDev, results[0]->GetBinContent(x), standardDev / results[0]->GetBinContent(x));
1636 TH1* standardDeviation = CalculateStdDev(results, kErrorIterations);
1638 // compare maxError to RMS of profile (variable name: average)
1639 // first: calculate rms in percent of value
1640 TH1* rmsError = (TH1*) maxError->Clone("rmsError");
1643 // enable error to be standard deviation (see http://root.cern.ch/root/html/TProfile.html#TProfile:SetErrorOption)
1644 average->SetErrorOption("s");
1645 for (Int_t x=1; x<=rmsError->GetNbinsX(); x++)
1646 if (average->GetBinContent(x) > 0)
1647 rmsError->SetBinContent(x, average->GetBinError(x) / average->GetBinContent(x));
1649 // find maxError deviation from average (not from "first result"/mc as above)
1650 TH1* maxError2 = (TH1*) maxError->Clone("maxError2");
1652 for (Int_t n=1; n<kErrorIterations; ++n)
1653 for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1654 if (average->GetBinContent(x) > 0)
1655 maxError2->SetBinContent(x, TMath::Max(maxError2->GetBinContent(x), TMath::Abs((results[n]->GetBinContent(x) - average->GetBinContent(x)) / average->GetBinContent(x))));
1657 //new TCanvas; maxError2->DrawCopy(); rmsError->SetLineColor(2); rmsError->DrawCopy("SAME"); standardDeviation->SetLineColor(3); standardDeviation->DrawCopy("SAME");
1659 // plot difference between average and MC/first result
1660 TH1* averageFirstRatio = (TH1*) results[0]->Clone("averageFirstRatio");
1661 averageFirstRatio->Reset();
1662 averageFirstRatio->Divide(results[0], average);
1664 //new TCanvas; results[0]->DrawCopy(); average->SetLineColor(2); average->DrawClone("SAME");
1665 //new TCanvas; averageFirstRatio->DrawCopy();
1667 static TH1* temp = 0;
1670 temp = (TH1*) standardDeviation->Clone("temp");
1671 for (Int_t x=1; x<=results[0]->GetNbinsX(); x++)
1672 temp->SetBinContent(x, temp->GetBinContent(x) * results[0]->GetBinContent(x));
1676 // find difference from result[0] as TH2
1677 TH2F* pulls = new TH2F("pulls", "pulls;multiplicity;difference", nBins, lowEdge, upEdge, 1000, -10, 10);
1678 for (Int_t n=1; n<kErrorIterations; ++n)
1679 for (Int_t x=1; x<=results[n]->GetNbinsX(); x++)
1680 if (temp->GetBinContent(x) > 0)
1681 pulls->Fill(results[n]->GetXaxis()->GetBinCenter(x), (results[0]->GetBinContent(x) - results[n]->GetBinContent(x)) / temp->GetBinContent(x));
1682 new TCanvas("pulls", "pulls", 800, 600); pulls->DrawCopy(); pulls->FitSlicesY();
1686 for (Int_t n=0; n<kErrorIterations; ++n)
1690 // fill into result histogram
1691 for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
1692 fMultiplicityESDCorrected[correlationID]->SetBinContent(i, firstResult->GetBinContent(i));
1694 for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
1695 fMultiplicityESDCorrected[correlationID]->SetBinError(i, standardDeviation->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
1697 return standardDeviation;
1700 //____________________________________________________________________
1701 void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType, Float_t regPar, Int_t nIterations, TH1* initialConditions, Int_t determineError)
1704 // correct spectrum using bayesian method
1708 // 1 = from randomizing
1709 // 2 = with UnfoldGetBias
1711 // initialize seed with current time
1712 gRandom->SetSeed(0);
1714 SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
1716 // normalize correction for given nPart
1717 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1719 // with this it is normalized to 1
1720 Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1722 // with this normalized to the given efficiency
1723 if (fCurrentEfficiency->GetBinContent(i) > 0)
1724 sum /= fCurrentEfficiency->GetBinContent(i);
1728 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1732 fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1733 fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1737 fCurrentCorrelation->SetBinContent(i, j, 0);
1738 fCurrentCorrelation->SetBinError(i, j, 0);
1743 Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1745 AliUnfolding::SetBayesianParameters(regPar, nIterations);
1746 AliUnfolding::SetUnfoldingMethod(AliUnfolding::kBayesian);
1748 if (determineError <= 1)
1750 if (AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID]) != 0)
1753 else if (determineError == 2)
1755 AliUnfolding::UnfoldGetBias(fCurrentCorrelation, fCurrentEfficiency, fCurrentESD, initialConditions, fMultiplicityESDCorrected[correlationID]);
1759 if (determineError == 0)
1761 Printf("AliMultiplicityCorrection::ApplyBayesianMethod: WARNING: No errors calculated.");
1765 // evaluate errors, this is done by randomizing the measured spectrum following Poission statistics
1766 // this (new) measured spectrum is then unfolded and the different to the result from the "real" measured
1767 // spectrum calculated. This is performed N times and the sigma is taken as the statistical
1768 // error of the unfolding method itself.
1770 const Int_t kErrorIterations = 20;
1772 Printf("Spectrum unfolded. Determining error (%d iterations)...", kErrorIterations);
1774 TH1* randomized = (TH1*) fCurrentESD->Clone("randomized");
1775 TH1* resultArray[kErrorIterations+1];
1776 for (Int_t n=0; n<kErrorIterations; ++n)
1778 // randomize the content of clone following a poisson with the mean = the value of that bin
1779 for (Int_t x=1; x<=randomized->GetNbinsX(); x++) // mult. axis
1781 Int_t randomValue = gRandom->Poisson(fCurrentESD->GetBinContent(x));
1782 //printf("%d --> %d\n", fCurrentESD->GetBinContent(x), randomValue);
1783 randomized->SetBinContent(x, randomValue);
1784 randomized->SetBinError(x, TMath::Sqrt(randomValue));
1787 TH1* result2 = (TH1*) fMultiplicityESDCorrected[correlationID]->Clone("result2");
1789 if (AliUnfolding::Unfold(fCurrentCorrelation, fCurrentEfficiency, randomized, initialConditions, result2) != 0)
1795 resultArray[n+1] = result2;
1799 resultArray[0] = fMultiplicityESDCorrected[correlationID];
1800 TH1* error = CalculateStdDev(resultArray, kErrorIterations+1);
1802 for (Int_t n=0; n<kErrorIterations; ++n)
1803 delete resultArray[n+1];
1805 Printf("Comparing bias and error:");
1806 for (Int_t i=1; i<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); ++i)
1808 Printf("Bin %d: Content: %f Error: %f Bias: %f", i, fMultiplicityESDCorrected[correlationID]->GetBinContent(i), error->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i), fMultiplicityESDCorrected[correlationID]->GetBinError(i));
1809 fMultiplicityESDCorrected[correlationID]->SetBinError(i, error->GetBinContent(i) * fMultiplicityESDCorrected[correlationID]->GetBinContent(i));
1815 //____________________________________________________________________
1816 Float_t AliMultiplicityCorrection::BayesCovarianceDerivate(Float_t matrixM[251][251], const TH2* hResponse, Int_t k, Int_t i, Int_t r, Int_t u)
1819 // helper function for the covariance matrix of the bayesian method
1824 if (k == u && r == i)
1825 result += 1.0 / hResponse->GetBinContent(u+1, r+1);
1828 result -= 1.0 / fCurrentEfficiency->GetBinContent(u+1);
1831 result -= matrixM[u][i] * fCurrentEfficiency->GetBinContent(u+1) / hResponse->GetBinContent(u+1, i+1);
1833 result *= matrixM[k][i];
1838 //____________________________________________________________________
1839 void AliMultiplicityCorrection::ApplyLaszloMethod(Int_t inputRange, Bool_t fullPhaseSpace, EventType eventType)
1842 // correct spectrum using bayesian method
1847 Int_t correlationID = inputRange; // + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1848 Int_t mcTarget = inputRange; //((fullPhaseSpace == kFALSE) ? inputRange : 4);
1850 SetupCurrentHists(inputRange, fullPhaseSpace, eventType);
1852 fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
1854 // TODO should be taken from correlation map
1855 //TH1* sumHist = GetMultiplicityMC(inputRange, eventType)->ProjectionY("sumHist", 1, GetMultiplicityMC(inputRange, eventType)->GetNbinsX());
1857 // normalize correction for given nPart
1858 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1860 Double_t sum = fCurrentCorrelation->Integral(i, i, 1, fCurrentCorrelation->GetNbinsY());
1861 //Double_t sum = sumHist->GetBinContent(i);
1864 for (Int_t j=1; j<=fCurrentCorrelation->GetNbinsY(); ++j)
1867 fCurrentCorrelation->SetBinContent(i, j, fCurrentCorrelation->GetBinContent(i, j) / sum);
1868 fCurrentCorrelation->SetBinError(i, j, fCurrentCorrelation->GetBinError(i, j) / sum);
1873 fCurrentCorrelation->Draw("COLZ");
1876 fCurrentEfficiency = ((TH2*) fCurrentCorrelation)->ProjectionX("eff");
1878 // pick prior distribution
1879 TH1F* hPrior = (TH1F*)fCurrentESD->Clone("prior");
1880 Float_t norm = 1; //hPrior->Integral();
1881 for (Int_t t=1; t<=hPrior->GetNbinsX(); t++)
1882 hPrior->SetBinContent(t, hPrior->GetBinContent(t)/norm);
1884 // zero distribution
1885 TH1F* zero = (TH1F*)hPrior->Clone("zero");
1888 TH1F* hTemp = (TH1F*)fCurrentESD->Clone("temp");
1892 TH2F* hResponse = (TH2F*) fCurrentCorrelation;
1895 Int_t iterations = 25;
1896 for (Int_t i=0; i<iterations; i++)
1898 //printf(" iteration %i \n", i);
1900 for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
1903 for (Int_t t = 1; t<=hResponse->GetNbinsX(); t++)
1904 value += hResponse->GetBinContent(t, m) * hPrior->GetBinContent(t);
1905 hTemp->SetBinContent(m, value);
1906 //printf("%d %f %f %f\n", m, zero->GetBinContent(m), hPrior->GetBinContent(m), value);
1909 // regularization (simple smoothing)
1910 TH1F* hTrueSmoothed = (TH1F*) hTemp->Clone("truesmoothed");
1912 for (Int_t t=2; t<hTrueSmoothed->GetNbinsX(); t++)
1914 Float_t average = (hTemp->GetBinContent(t-1) / hTemp->GetBinWidth(t-1)
1915 + hTemp->GetBinContent(t) / hTemp->GetBinWidth(t)
1916 + hTemp->GetBinContent(t+1) / hTemp->GetBinWidth(t+1)) / 3.;
1917 average *= hTrueSmoothed->GetBinWidth(t);
1919 // weight the average with the regularization parameter
1920 hTrueSmoothed->SetBinContent(t, (1-regPar) * hTemp->GetBinContent(t) + regPar * average);
1923 for (Int_t m=1; m<=hResponse->GetNbinsY(); m++)
1924 hTemp->SetBinContent(m, zero->GetBinContent(m) + hPrior->GetBinContent(m) - hTrueSmoothed->GetBinContent(m));
1927 for (Int_t t=1; t<=fMultiplicityESDCorrected[correlationID]->GetNbinsX(); t++)
1929 fMultiplicityESDCorrected[correlationID]->SetBinContent(t, hTemp->GetBinContent(t));
1930 fMultiplicityESDCorrected[correlationID]->SetBinError(t, 0.05 * hTemp->GetBinContent(t)); // TODO
1932 //printf(" bin %d content %f \n", t, fMultiplicityESDCorrected[correlationID]->GetBinContent(t));
1936 // calculate chi2 (change from last iteration)
1939 // use smoothed true (normalized) as new prior
1940 norm = 1; //hTrueSmoothed->Integral();
1942 for (Int_t t=1; t<hTrueSmoothed->GetNbinsX(); t++)
1944 Float_t newValue = hTemp->GetBinContent(t)/norm;
1945 Float_t diff = hPrior->GetBinContent(t) - newValue;
1946 chi2 += (Double_t) diff * diff;
1948 hPrior->SetBinContent(t, newValue);
1951 printf("Chi2 of %d iteration = %.10f\n", i, chi2);
1954 DrawComparison(Form("Laszlo_%d", i), inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
1956 delete hTrueSmoothed;
1957 } // end of iterations
1959 DrawComparison("Laszlo", inputRange, fullPhaseSpace, kTRUE, GetMultiplicityMC(mcTarget, eventType)->ProjectionY());
1962 //____________________________________________________________________
1963 void AliMultiplicityCorrection::ApplyGaussianMethod(Int_t inputRange, Bool_t fullPhaseSpace)
1966 // correct spectrum using a simple Gaussian approach, that is model-dependent
1969 Int_t correlationID = inputRange; // + ((fullPhaseSpace == kFALSE) ? 0 : 4);
1970 Int_t mcTarget = inputRange; //((fullPhaseSpace == kFALSE) ? inputRange : 4);
1972 SetupCurrentHists(inputRange, fullPhaseSpace, kTrVtx);
1974 fCurrentESD->Scale(1.0 / fCurrentESD->Integral());
1976 TH1D* correction = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianMean"));
1977 correction->SetTitle("GaussianMean");
1979 TH1D* correctionWidth = dynamic_cast<TH1D*> (fCurrentESD->Clone("GaussianWidth"));
1980 correction->SetTitle("GaussianWidth");
1982 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
1984 TH1D* proj = (dynamic_cast<TH2*> (fCurrentCorrelation))->ProjectionX("_px", i, i+1);
1985 proj->Fit("gaus", "0Q");
1986 correction->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(1));
1987 correctionWidth->SetBinContent(i, proj->GetFunction("gaus")->GetParameter(2));
1989 // draw for debugging
1992 proj->GetFunction("gaus")->DrawCopy("SAME");
1996 TH1* target = fMultiplicityESDCorrected[correlationID];
1998 Int_t nBins = target->GetNbinsX()*10+1;
1999 Float_t* binning = new Float_t[nBins];
2000 for (Int_t i=1; i<=target->GetNbinsX(); ++i)
2001 for (Int_t j=0; j<10; ++j)
2002 binning[(i-1)*10 + j] = target->GetXaxis()->GetBinLowEdge(i) + target->GetXaxis()->GetBinWidth(i) / 10 * j;
2004 binning[nBins-1] = target->GetXaxis()->GetBinUpEdge(target->GetNbinsX());
2006 TH1F* fineBinned = new TH1F("targetFineBinned", "targetFineBinned", nBins-1, binning);
2008 for (Int_t i=1; i<=fCurrentCorrelation->GetNbinsX(); ++i)
2010 Float_t mean = correction->GetBinContent(i);
2011 Float_t width = correctionWidth->GetBinContent(i);
2013 Int_t fillBegin = fineBinned->FindBin(mean - width * 5);
2014 Int_t fillEnd = fineBinned->FindBin(mean + width * 5);
2015 //printf("bin %d mean %f width %f, filling from %d to %d\n", i, mean, width, fillBegin, fillEnd);
2017 for (Int_t j=fillBegin; j <= fillEnd; ++j)
2019 fineBinned->AddBinContent(j, TMath::Gaus(fineBinned->GetXaxis()->GetBinCenter(j), mean, width, kTRUE) * fCurrentESD->GetBinContent(i));
2023 for (Int_t i=1; i<=target->GetNbinsX(); ++i)
2026 for (Int_t j=1; j<=10; ++j)
2027 sum += fineBinned->GetBinContent((i-1)*10 + j);
2028 target->SetBinContent(i, sum / 10);
2033 DrawComparison("Gaussian", inputRange, fullPhaseSpace, kFALSE, GetMultiplicityMC(mcTarget, kTrVtx)->ProjectionY());
2036 //____________________________________________________________________
2037 TH1* AliMultiplicityCorrection::GetConvoluted(Int_t i, EventType eventType)
2039 // convolutes the corrected histogram i with the response matrix
2041 TH1* corrected = (TH1*) fMultiplicityESDCorrected[i]->Clone("corrected");
2043 // undo efficiency correction (hist must not be deleted, is reused)
2044 TH1* efficiency = GetEfficiency(i, eventType);
2045 //new TCanvas; efficiency->DrawCopy();
2046 corrected->Multiply(efficiency);
2048 TH2* convoluted = CalculateMultiplicityESD(corrected, i);
2049 TH1* convolutedProj = convoluted->ProjectionY("GetConvoluted_convolutedProj", 1, convoluted->GetNbinsX());
2054 return convolutedProj;
2057 //____________________________________________________________________
2058 TH1* AliMultiplicityCorrection::GetResiduals(Int_t i, EventType eventType, Float_t& residualSum)
2060 // creates the residual histogram from the corrected histogram i corresponding to an eventType event sample using the corresponding correlation matrix
2061 // residual is : M - UT / eM
2062 // residualSum contains the squared sum of the residuals
2064 TH1* corrected = (TH1*) fMultiplicityESDCorrected[i]->Clone("corrected");
2065 TH1* convolutedProj = GetConvoluted(i, eventType);
2068 Int_t end = fMultiplicityESD[i]->GetNbinsX();
2069 if (fVtxEnd > fVtxBegin)
2074 TH1* measuredProj = fMultiplicityESD[i]->ProjectionY("measuredProj", begin, end);
2076 TH1* residuals = (TH1*) measuredProj->Clone("GetResiduals_residuals");
2077 residuals->SetTitle(";measured multiplicity;residuals (M-Ut)/e");
2078 residuals->Add(convolutedProj, -1);
2081 for (Int_t j=1; j<=residuals->GetNbinsX(); j++)
2083 if (measuredProj->GetBinContent(j) > 0)
2084 residuals->SetBinContent(j, residuals->GetBinContent(j) / TMath::Sqrt(measuredProj->GetBinContent(j)));
2085 residuals->SetBinError(j, 0);
2088 residualSum += residuals->GetBinContent(j) * residuals->GetBinContent(j);
2092 delete convolutedProj;
2093 delete measuredProj;
2098 //____________________________________________________________________
2099 TH2F* AliMultiplicityCorrection::CalculateMultiplicityESD(TH1* inputMC, Int_t correlationMap)
2101 // runs the distribution given in inputMC through the response matrix identified by
2102 // correlationMap and produces a measured distribution
2103 // although it is a TH2F the vertex axis is not used at the moment and all entries are filled in mid-vertex
2108 if (correlationMap < 0 || correlationMap >= kCorrHists)
2111 // empty under/overflow bins in x, otherwise Project3D takes them into account
2112 TH3* hist = fCorrelation[correlationMap];
2113 for (Int_t y=0; y<=hist->GetYaxis()->GetNbins()+1; ++y)
2115 for (Int_t z=0; z<=hist->GetZaxis()->GetNbins()+1; ++z)
2117 hist->SetBinContent(0, y, z, 0);
2118 hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0);
2122 if (fVtxEnd > fVtxBegin)
2123 hist->GetXaxis()->SetRange(fVtxBegin, fVtxEnd);
2125 TH2* corr = (TH2*) hist->Project3D("zy");
2126 //corr->Rebin2D(2, 1);
2128 Printf("Correction histogram used for convolution has %f entries", corr->Integral());
2130 // normalize correction for given nPart
2131 for (Int_t i=1; i<=corr->GetNbinsX(); ++i)
2133 Double_t sum = corr->Integral(i, i, 1, corr->GetNbinsY());
2137 for (Int_t j=1; j<=corr->GetNbinsY(); ++j)
2140 corr->SetBinContent(i, j, corr->GetBinContent(i, j) / sum);
2141 corr->SetBinError(i, j, corr->GetBinError(i, j) / sum);
2145 TH2F* target = static_cast<TH2F*> (fMultiplicityESD[0]->Clone(Form("%s_measured", inputMC->GetName())));
2148 for (Int_t meas=1; meas<=corr->GetNbinsY(); ++meas)
2150 Float_t measured = 0;
2153 for (Int_t gen=1; gen<=corr->GetNbinsX(); ++gen)
2155 Int_t mcGenBin = inputMC->GetXaxis()->FindBin(corr->GetXaxis()->GetBinCenter(gen));
2157 measured += inputMC->GetBinContent(mcGenBin) * corr->GetBinContent(gen, meas);
2159 error += inputMC->GetBinError(mcGenBin) * corr->GetBinContent(gen, meas);
2162 //printf("%f +- %f ; %f +- %f \n", inputMC->GetBinContent(meas), inputMC->GetBinError(meas), measured, error);
2164 target->SetBinContent(1 + target->GetNbinsX() / 2, meas, measured);
2165 target->SetBinError(1 + target->GetNbinsX() / 2, meas, error);
2171 //____________________________________________________________________
2172 void AliMultiplicityCorrection::SetGenMeasFromFunc(const TF1* inputMC, Int_t id)
2174 // uses the given function to fill the input MC histogram and generates from that
2175 // the measured histogram by applying the response matrix
2176 // this can be used to evaluate if the methods work indepedently of the input
2178 // WARNING does not respect the vertex distribution, just fills central vertex bin
2183 if (id < 0 || id >= kESDHists)
2186 // fill histogram used for random generation
2187 TH1* tmp = fMultiplicityVtx[id]->ProjectionY("tmp");
2190 for (Int_t i=1; i<=tmp->GetNbinsX(); ++i)
2191 tmp->SetBinContent(i, inputMC->Eval(tmp->GetXaxis()->GetBinCenter(i)) * tmp->GetXaxis()->GetBinWidth(i));
2193 TH1* mcRnd = fMultiplicityVtx[id]->ProjectionY("mcRnd");
2195 mcRnd->FillRandom(tmp, (Int_t) tmp->Integral());
2197 //new TCanvas; tmp->Draw();
2198 //new TCanvas; mcRnd->Draw();
2200 // and move into 2d histogram
2201 TH1* mc = fMultiplicityVtx[id];
2203 for (Int_t i=1; i<=mc->GetNbinsY(); ++i)
2205 mc->SetBinContent(mc->GetNbinsX() / 2 + 1, i, mcRnd->GetBinContent(i));
2206 mc->SetBinError(mc->GetNbinsX() / 2 + 1, i, TMath::Sqrt(mcRnd->GetBinContent(i)));
2209 //new TCanvas; mc->Draw("COLZ");
2211 // now randomize the measured histogram; funcMeasured is used as pilot function to generated the measured entries
2212 TH1* funcMeasured = CalculateMultiplicityESD(tmp, id)->ProjectionY("funcMeasured");
2214 //new TCanvas; funcMeasured->Draw();
2216 fMultiplicityESD[id]->Reset();
2218 TH1* measRnd = fMultiplicityESD[id]->ProjectionY("measRnd");
2219 measRnd->FillRandom(funcMeasured, (Int_t) tmp->Integral());
2221 //new TCanvas; measRnd->Draw();
2223 fMultiplicityESD[id]->Reset();
2224 for (Int_t i=1; i<=fMultiplicityESD[id]->GetNbinsY(); ++i)
2226 fMultiplicityESD[id]->SetBinContent(fMultiplicityESD[id]->GetNbinsX() / 2 + 1, i, measRnd->GetBinContent(i));
2227 fMultiplicityESD[id]->SetBinError(fMultiplicityESD[id]->GetNbinsX() / 2 + 1, i, TMath::Sqrt(measRnd->GetBinContent(i)));