]>
Commit | Line | Data |
---|---|---|
6d75bdb8 | 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 | // This class is used to store correlation maps, generated and reconstructed data of the jet spectrum | |
17 | // It also contains functions to correct the spectrum using the bayesian unfolding | |
18 | // | |
19 | ||
20 | #include "AliJetSpectrumUnfolding.h" | |
21 | ||
22 | #include <TFile.h> | |
23 | #include <TH1F.h> | |
24 | #include <TH2F.h> | |
25 | #include <TH3F.h> | |
22ead7c4 | 26 | #include <THnSparse.h> |
6d75bdb8 | 27 | #include <TDirectory.h> |
28 | #include <TVirtualFitter.h> | |
29 | #include <TCanvas.h> | |
30 | #include <TString.h> | |
31 | #include <TF1.h> | |
32 | #include <TF2.h> | |
33 | #include <TMath.h> | |
34 | #include <TCollection.h> | |
35 | #include <TLegend.h> | |
36 | #include <TLine.h> | |
37 | #include <TRandom.h> | |
38 | #include <TProfile.h> | |
39 | #include <TProfile2D.h> | |
40 | #include <TStyle.h> | |
41 | #include <TColor.h> | |
8de9091a | 42 | #include <TVectorD.h> |
6d75bdb8 | 43 | |
44 | #include <THnSparse.h> | |
45 | ||
46 | ClassImp(AliJetSpectrumUnfolding) | |
47 | ||
8de9091a | 48 | const Int_t AliJetSpectrumUnfolding::fgkNBINSE = 50; |
49 | const Int_t AliJetSpectrumUnfolding::fgkNBINSZ = 50; | |
50 | const Int_t AliJetSpectrumUnfolding::fgkNEVENTS = 500000; | |
51 | const Double_t AliJetSpectrumUnfolding::fgkaxisLowerLimitE = 0.; | |
52 | const Double_t AliJetSpectrumUnfolding::fgkaxisLowerLimitZ = 0.; | |
53 | const Double_t AliJetSpectrumUnfolding::fgkaxisUpperLimitE = 250.; | |
54 | const Double_t AliJetSpectrumUnfolding::fgkaxisUpperLimitZ = 1.; | |
6d75bdb8 | 55 | |
56 | Float_t AliJetSpectrumUnfolding::fgBayesianSmoothing = 1; // smoothing parameter (0 = no smoothing) | |
57 | Int_t AliJetSpectrumUnfolding::fgBayesianIterations = 100; // number of iterations in Bayesian method | |
58 | ||
59 | //____________________________________________________________________ | |
60 | ||
61 | AliJetSpectrumUnfolding::AliJetSpectrumUnfolding() : | |
8de9091a | 62 | TNamed(), fCurrentRec(0), fCurrentCorrelation(0), fRecSpectrum(0), fGenSpectrum(0), |
63 | fUnfSpectrum(0), fCorrelation(0), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0) | |
6d75bdb8 | 64 | { |
65 | // | |
66 | // default constructor | |
67 | // | |
68 | ||
69 | fGenSpectrum = 0; | |
70 | fRecSpectrum = 0; | |
71 | fUnfSpectrum = 0; | |
72 | fCorrelation = 0; | |
73 | } | |
74 | ||
75 | //____________________________________________________________________ | |
76 | AliJetSpectrumUnfolding::AliJetSpectrumUnfolding(const Char_t* name, const Char_t* title) : | |
8de9091a | 77 | TNamed(name, title), fCurrentRec(0), fCurrentCorrelation(0), fRecSpectrum(0), |
78 | fGenSpectrum(0), fUnfSpectrum(0), fCorrelation(0), fLastChi2MC(0), fLastChi2MCLimit(0), fLastChi2Residuals(0), fRatioAverage(0) | |
6d75bdb8 | 79 | { |
80 | // | |
81 | // named constructor | |
82 | // | |
83 | ||
84 | // do not add this hists to the directory | |
85 | Bool_t oldStatus = TH1::AddDirectoryStatus(); | |
86 | TH1::AddDirectory(kFALSE); | |
8de9091a | 87 | fRecSpectrum = new TH2F("fRecSpectrum", "Reconstructed Spectrum;E^{jet}_{rec} [GeV];z^{lp}_{rec}", |
88 | fgkNBINSE, fgkaxisLowerLimitE, fgkaxisUpperLimitE, | |
89 | fgkNBINSZ, fgkaxisLowerLimitZ, fgkaxisUpperLimitZ); | |
90 | fGenSpectrum = new TH2F("fGenSpectrum", "Generated Spectrum;E^{jet}_{gen} [GeV];z^{lp}_{gen}", | |
91 | fgkNBINSE, fgkaxisLowerLimitE, fgkaxisUpperLimitE, | |
92 | fgkNBINSZ, fgkaxisLowerLimitZ, fgkaxisUpperLimitZ); | |
93 | fUnfSpectrum = new TH2F("fUnfSpectrum", "Unfolded Spectrum;E^{jet} [GeV];z^{lp}", | |
94 | fgkNBINSE, fgkaxisLowerLimitE, fgkaxisUpperLimitE, | |
95 | fgkNBINSZ, fgkaxisLowerLimitZ, fgkaxisUpperLimitZ); | |
96 | ||
97 | const Int_t nbin[4]={fgkNBINSE, fgkNBINSE, fgkNBINSZ, fgkNBINSZ}; | |
6d75bdb8 | 98 | //arrays for bin limits |
8de9091a | 99 | Double_t lowEdge[4] = {fgkaxisLowerLimitE, fgkaxisLowerLimitE, fgkaxisLowerLimitZ, fgkaxisLowerLimitZ}; |
100 | Double_t upEdge[4] = {fgkaxisUpperLimitE, fgkaxisUpperLimitE, fgkaxisUpperLimitZ, fgkaxisUpperLimitZ}; | |
6d75bdb8 | 101 | |
8de9091a | 102 | fCorrelation = new THnSparseF("fCorrelation", "Correlation Function", 4, nbin, lowEdge, upEdge); |
6d75bdb8 | 103 | |
104 | TH1::AddDirectory(oldStatus); | |
105 | } | |
106 | ||
107 | //____________________________________________________________________ | |
108 | AliJetSpectrumUnfolding::~AliJetSpectrumUnfolding() | |
109 | { | |
110 | // | |
111 | // Destructor | |
112 | // | |
113 | ||
114 | if (fGenSpectrum) | |
115 | delete fGenSpectrum; | |
116 | fGenSpectrum = 0; | |
117 | ||
118 | if (fRecSpectrum) | |
119 | delete fRecSpectrum; | |
120 | fRecSpectrum = 0; | |
121 | ||
122 | if (fUnfSpectrum) | |
123 | delete fUnfSpectrum; | |
124 | fUnfSpectrum = 0; | |
125 | ||
126 | if (fCorrelation) | |
127 | delete fCorrelation; | |
128 | fCorrelation = 0; | |
129 | ||
130 | } | |
131 | ||
132 | //____________________________________________________________________ | |
133 | Long64_t AliJetSpectrumUnfolding::Merge(TCollection* list) | |
134 | { | |
135 | // Merge a list of AliJetSpectrumUnfolding objects with this (needed for | |
136 | // PROOF). | |
137 | // Returns the number of merged objects (including this). | |
138 | ||
139 | if (!list) | |
140 | return 0; | |
141 | ||
142 | if (list->IsEmpty()) | |
143 | return 1; | |
144 | ||
145 | TIterator* iter = list->MakeIterator(); | |
146 | TObject* obj; | |
147 | ||
148 | // collections of all histograms | |
149 | TList collections[4]; | |
150 | ||
151 | Int_t count = 0; | |
152 | while ((obj = iter->Next())) { | |
153 | ||
154 | AliJetSpectrumUnfolding* entry = dynamic_cast<AliJetSpectrumUnfolding*> (obj); | |
155 | if (entry == 0) | |
156 | continue; | |
157 | ||
158 | collections[0].Add(entry->fGenSpectrum); | |
159 | collections[1].Add(entry->fRecSpectrum); | |
160 | collections[2].Add(entry->fUnfSpectrum); | |
161 | collections[3].Add(entry->fCorrelation); | |
162 | ||
163 | count++; | |
164 | } | |
165 | ||
166 | fGenSpectrum->Merge(&collections[0]); | |
167 | fRecSpectrum->Merge(&collections[1]); | |
168 | fUnfSpectrum->Merge(&collections[2]); | |
169 | fCorrelation->Merge(&collections[3]); | |
170 | ||
171 | delete iter; | |
172 | ||
173 | return count+1; | |
174 | } | |
175 | ||
176 | //____________________________________________________________________ | |
177 | Bool_t AliJetSpectrumUnfolding::LoadHistograms(const Char_t* dir) | |
178 | { | |
179 | // | |
180 | // loads the histograms from a file | |
181 | // if dir is empty a directory with the name of this object is taken (like in SaveHistogram) | |
182 | // | |
183 | ||
184 | if (!dir) | |
185 | dir = GetName(); | |
186 | ||
187 | if (!gDirectory->cd(dir)) | |
188 | return kFALSE; | |
189 | ||
190 | Bool_t success = kTRUE; | |
191 | ||
192 | // store old histograms to delete them later | |
193 | TList oldHistograms; | |
194 | oldHistograms.SetOwner(1); | |
195 | ||
196 | if (fGenSpectrum) oldHistograms.Add(fGenSpectrum); | |
b85b355c | 197 | else return kFALSE; |
6d75bdb8 | 198 | if (fRecSpectrum) oldHistograms.Add(fRecSpectrum); |
b85b355c | 199 | else return kFALSE; |
6d75bdb8 | 200 | if (fUnfSpectrum) oldHistograms.Add(fUnfSpectrum); |
b85b355c | 201 | else return kFALSE; |
6d75bdb8 | 202 | if (fCorrelation) oldHistograms.Add(fCorrelation); |
b85b355c | 203 | else return kFALSE; |
6d75bdb8 | 204 | |
205 | // load new histograms | |
206 | fGenSpectrum = dynamic_cast<TH2F*> (gDirectory->Get(fGenSpectrum->GetName())); | |
207 | if (!fGenSpectrum) | |
208 | success = kFALSE; | |
209 | ||
210 | fRecSpectrum = dynamic_cast<TH2F*> (gDirectory->Get(fRecSpectrum->GetName())); | |
211 | if (!fRecSpectrum) | |
212 | success = kFALSE; | |
213 | ||
214 | fUnfSpectrum = dynamic_cast<TH2F*> (gDirectory->Get(fUnfSpectrum->GetName())); | |
215 | if (!fUnfSpectrum) | |
216 | success = kFALSE; | |
217 | ||
218 | fCorrelation = dynamic_cast<THnSparseF*> (gDirectory->Get(fCorrelation->GetName())); | |
219 | if (!fCorrelation) | |
220 | success = kFALSE; | |
221 | ||
222 | gDirectory->cd(".."); | |
223 | ||
224 | // delete old histograms | |
225 | oldHistograms.Delete(); | |
226 | ||
227 | return success; | |
228 | } | |
229 | ||
230 | //____________________________________________________________________ | |
231 | void AliJetSpectrumUnfolding::SaveHistograms() | |
232 | { | |
233 | // | |
234 | // saves the histograms | |
235 | // | |
236 | ||
237 | gDirectory->mkdir(GetName()); | |
238 | gDirectory->cd(GetName()); | |
239 | ||
240 | if (fGenSpectrum) | |
241 | fGenSpectrum->Write(); | |
242 | ||
243 | if (fRecSpectrum) | |
244 | fRecSpectrum->Write(); | |
245 | ||
246 | if (fUnfSpectrum) | |
247 | fUnfSpectrum->Write(); | |
248 | ||
249 | if (fCorrelation) | |
250 | fCorrelation->Write(); | |
251 | ||
252 | gDirectory->cd(".."); | |
253 | } | |
254 | ||
255 | //____________________________________________________________________ | |
256 | void AliJetSpectrumUnfolding::SetupCurrentHists(Bool_t createBigBin) | |
257 | { | |
258 | // | |
259 | // resets fUnfSpectrum | |
260 | // | |
261 | ||
262 | fUnfSpectrum->Reset(); | |
263 | fUnfSpectrum->Sumw2(); | |
264 | ||
265 | fCurrentRec = (TH2F*)fRecSpectrum->Clone("fCurrentRec"); | |
266 | fCurrentRec->Sumw2(); | |
267 | ||
268 | fCurrentCorrelation = (THnSparseF*)fCorrelation->Clone("fCurrentCorrelation"); | |
269 | fCurrentCorrelation->Sumw2(); | |
270 | ||
fd5db301 | 271 | Printf("Correlation Matrix has %d filled bins", (Int_t)fCurrentCorrelation->GetNbins()); |
df65bddb | 272 | |
6d75bdb8 | 273 | if (createBigBin) |
274 | { | |
275 | Int_t maxBinE = 0, maxBinZ = 0; | |
276 | Float_t maxE = 0, maxZ = 0; | |
277 | for (Int_t me=1; me<=fCurrentRec->GetNbinsX(); me++) | |
278 | for (Int_t mz=1; mz<=fCurrentRec->GetNbinsY(); mz++) | |
279 | { | |
8de9091a | 280 | if (fCurrentRec->GetBinContent(me,mz) <= 5 && me>fgkNBINSE/2 && mz>fgkNBINSZ/2) |
6d75bdb8 | 281 | { |
282 | maxBinE = me; | |
283 | maxBinZ = mz; | |
284 | maxE = fCurrentRec->GetXaxis()->GetBinCenter(me); | |
285 | maxZ = fCurrentRec->GetYaxis()->GetBinCenter(mz); | |
286 | break; | |
287 | } | |
288 | } | |
289 | ||
290 | if (maxBinE > 0 || maxBinZ > 0) | |
291 | { | |
292 | printf("Bin limit in measured spectrum is e = %d and z = %d.\n", maxBinE, maxBinZ); | |
293 | fCurrentRec->SetBinContent(maxBinE, maxBinZ, fCurrentRec->Integral(maxBinE, fCurrentRec->GetNbinsX(), maxBinZ, fCurrentRec->GetNbinsY())); | |
294 | for (Int_t me=maxBinE+1; me<=fCurrentRec->GetNbinsX(); me++) | |
295 | for (Int_t mz=maxBinZ+1; mz<=fCurrentRec->GetNbinsY(); mz++) | |
296 | { | |
297 | fCurrentRec->SetBinContent(me, mz, 0); | |
298 | fCurrentRec->SetBinError(me, mz, 0); | |
299 | } | |
300 | // the error is set to sqrt(N), better solution possible?, sum of relative errors of all contributions??? | |
301 | fCurrentRec->SetBinError(maxBinE, maxBinZ, TMath::Sqrt(fCurrentRec->GetBinContent(maxBinE, maxBinZ))); | |
302 | ||
303 | printf("This bin has now %f +- %f entries\n", fCurrentRec->GetBinContent(maxBinE, maxBinZ), fCurrentRec->GetBinError(maxBinE, maxBinZ)); | |
304 | ||
305 | /* for (Int_t te=1; te<=NBINSE; te++) | |
306 | { | |
307 | for (Int_t tz=1; tz<=NBINSZ; tz++) | |
308 | { | |
309 | Int_t binMin[4] = {te, maxBinE, tz, maxBinZ}; | |
310 | Int_t binMax[4] = {NBINSE, NBINSE, NBINSZ, NBINSZ}; | |
311 | Float_t sum=0; | |
312 | for (Int_t ite=te; ite<=NBINSE; ite++) | |
313 | for (Int_t itz=tz; itz<=NBINSZ; itz++) | |
314 | for (Int_t ime=maxBinE; ime<=NBINSE; ime++) | |
315 | for (Int_t imz=maxBinZ; imz<=NBINSZ; imz++) | |
316 | { | |
317 | Int_t bin[4] = {ite, ime, itz, imz}; | |
318 | sum += fCurrentCorrelation->GetBinContent(bin); | |
319 | } | |
320 | fCurrentCorrelation->SetBinContent(binMin, sum); | |
321 | fCurrentCorrelation->SetBinError(binMin, TMath::Sqrt(fCurrentCorrelation->GetBinContent(binMin))); | |
322 | printf("create big bin1, nbins = %d, te = %d, tz = %d\n", NBINSE, te, tz); | |
323 | for (Int_t me=maxBinE; me<=NBINSE; me++) | |
324 | { | |
325 | for (Int_t mz=maxBinZ; mz<=NBINSZ; mz++) | |
326 | { | |
327 | Int_t bin[4] = {te, me, tz, mz}; | |
328 | fCurrentCorrelation->SetBinContent(bin, 0.); | |
329 | fCurrentCorrelation->SetBinError(bin, 0.); | |
330 | printf("create big bin2\n"); | |
331 | } | |
332 | } | |
333 | } | |
334 | }*/ | |
335 | ||
336 | for(Int_t idx = 0; idx<=fCurrentCorrelation->GetNbins(); idx++) | |
337 | { | |
338 | Int_t bin[4]; | |
339 | Float_t binContent = fCurrentCorrelation->GetBinContent(idx,bin); | |
340 | Float_t binError = fCurrentCorrelation->GetBinError(idx); | |
341 | Int_t binMin[4] = {bin[0], maxBinE, bin[2], maxBinZ}; | |
8de9091a | 342 | if ( (bin[1]>maxBinE && bin[1]<=fgkNBINSE) && (bin[3]>maxBinZ && bin[3]<=fgkNBINSZ) ) |
6d75bdb8 | 343 | { |
344 | fCurrentCorrelation->SetBinContent(binMin, binContent + fCurrentCorrelation->GetBinContent(binMin)); | |
345 | fCurrentCorrelation->SetBinError(binMin, binError + TMath::Sqrt(fCurrentCorrelation->GetBinContent(binMin))); | |
346 | fCurrentCorrelation->SetBinContent(bin, 0.); | |
347 | fCurrentCorrelation->SetBinError(bin, 0.); | |
348 | } | |
8de9091a | 349 | printf("create big bin1, nbins = %d, te = %d, tz = %d\n", fgkNBINSE, bin[0], bin[1]); |
6d75bdb8 | 350 | } |
351 | ||
352 | printf("Adjusted correlation matrix!\n"); | |
353 | } | |
354 | } // end Create Big Bin | |
355 | ||
356 | } | |
357 | ||
358 | //____________________________________________________________________ | |
359 | void AliJetSpectrumUnfolding::SetBayesianParameters(Float_t smoothing, Int_t nIterations) | |
360 | { | |
361 | // | |
362 | // sets the parameters for Bayesian unfolding | |
363 | // | |
364 | ||
365 | fgBayesianSmoothing = smoothing; | |
366 | fgBayesianIterations = nIterations; | |
367 | ||
368 | printf("AliJetSpectrumUnfolding::SetBayesianParameters --> Paramaters set to %d iterations with smoothing %f\n", fgBayesianIterations, fgBayesianSmoothing); | |
369 | } | |
370 | ||
371 | //____________________________________________________________________ | |
8de9091a | 372 | void AliJetSpectrumUnfolding::NormalizeToBinWidth(TH2* const hist) |
6d75bdb8 | 373 | { |
374 | // | |
375 | // normalizes a 2-d histogram to its bin width (x width * y width) | |
376 | // | |
377 | ||
378 | for (Int_t i=1; i<=hist->GetNbinsX(); i++) | |
379 | for (Int_t j=1; j<=hist->GetNbinsY(); j++) | |
380 | { | |
381 | Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j); | |
382 | hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor); | |
383 | hist->SetBinError(i, j, hist->GetBinError(i, j) / factor); | |
384 | } | |
385 | } | |
386 | ||
387 | //____________________________________________________________________ | |
388 | void AliJetSpectrumUnfolding::DrawHistograms() | |
389 | { | |
390 | // | |
391 | // draws the histograms of this class | |
392 | // | |
393 | ||
394 | gStyle->SetPalette(1); | |
395 | ||
396 | TCanvas* canvas1 = new TCanvas("fRecSpectrum", "fRecSpectrum", 900, 600); | |
397 | gPad->SetLogz(); | |
398 | fRecSpectrum->DrawCopy("COLZ"); | |
399 | ||
400 | TCanvas* canvas2 = new TCanvas("fGenSpectrum", "fGenSpectrum", 900, 600); | |
8de9091a | 401 | canvas2->cd(); |
6d75bdb8 | 402 | gPad->SetLogz(); |
403 | fGenSpectrum->DrawCopy("COLZ"); | |
404 | ||
405 | TCanvas* canvas3 = new TCanvas("fUnfSpectrum", "fUnfSpectrum", 900, 600); | |
8de9091a | 406 | canvas3->cd(); |
6d75bdb8 | 407 | gPad->SetLogz(); |
408 | fUnfSpectrum->DrawCopy("COLZ"); | |
409 | ||
410 | TCanvas* canvas4 = new TCanvas("fCorrelation", "fCorrelation", 500, 500); | |
411 | canvas1->Divide(2); | |
412 | ||
413 | canvas4->cd(1); | |
414 | gPad->SetLogz(); | |
415 | TH2D* h0 = fCorrelation->Projection(1,0); | |
416 | h0->SetXTitle("E^{jet}_{gen} [GeV]"); | |
417 | h0->SetYTitle("E^{jet}_{rec} [GeV]"); | |
418 | h0->SetTitle("Projection: Jet Energy"); | |
419 | h0->DrawCopy("colz"); | |
420 | ||
421 | canvas1->cd(2); | |
422 | gPad->SetLogz(); | |
423 | TH2D* h1 = fCorrelation->Projection(3,2); | |
424 | h1->SetXTitle("z^{lp}_{gen}"); | |
425 | h1->SetYTitle("z^{lp}_{rec}"); | |
426 | h1->SetTitle("Projection: Leading Particle Fragmentation"); | |
427 | h1->DrawCopy("colz"); | |
428 | ||
429 | } | |
430 | ||
431 | //____________________________________________________________________ | |
8de9091a | 432 | void AliJetSpectrumUnfolding::DrawComparison(const char* name, TH2* const genHist) |
6d75bdb8 | 433 | { |
8de9091a | 434 | // |
435 | // Draws the copmparison plot (gen,rec and unfolded distributions | |
436 | // | |
6d75bdb8 | 437 | |
438 | if (fUnfSpectrum->Integral() == 0) | |
439 | { | |
440 | printf("ERROR. Unfolded histogram is empty\n"); | |
441 | return; | |
442 | } | |
443 | ||
444 | //regain measured distribution used for unfolding, because the bins were modified in SetupCurrentHists | |
445 | //in create big bin | |
446 | fCurrentRec = (TH2F*)fRecSpectrum->Clone(); | |
447 | fCurrentRec->Sumw2(); | |
448 | fCurrentRec->Scale(1.0 / fCurrentRec->Integral()); | |
449 | ||
450 | // normalize unfolded result to 1 | |
451 | fUnfSpectrum->Scale(1.0 / fUnfSpectrum->Integral()); | |
452 | ||
453 | // find bin with <= 100 entries. this is used as limit for the chi2 calculation | |
454 | Int_t mcBinLimitE = 0, mcBinLimitZ = 0; | |
455 | for (Int_t i=0; i<genHist->GetNbinsX(); ++i) | |
456 | for (Int_t j=0; j<genHist->GetNbinsY(); ++j) | |
457 | { | |
458 | if (genHist->GetBinContent(i,j) > 100) | |
459 | { | |
460 | mcBinLimitE = i; | |
461 | mcBinLimitZ = j; | |
462 | } | |
463 | else | |
464 | break; | |
465 | } | |
466 | Printf("AliJetSpectrumUnfolding::DrawComparison: Gen bin limit is (x,y) = (%d,%d)", mcBinLimitE,mcBinLimitZ); | |
467 | ||
468 | // scale to 1 true spectrum | |
469 | genHist->Sumw2(); | |
470 | genHist->Scale(1.0 / genHist->Integral()); | |
471 | ||
472 | // calculate residual | |
473 | // for that we convolute the response matrix with the gathered result | |
474 | TH2* tmpRecRecorrected = (TH2*) fUnfSpectrum->Clone("tmpRecRecorrected"); | |
475 | TH2* convoluted = CalculateRecSpectrum(tmpRecRecorrected); | |
476 | if (convoluted->Integral() > 0) | |
477 | convoluted->Scale(1.0 / convoluted->Integral()); | |
478 | else | |
479 | printf("ERROR: convoluted is empty. Something went wrong calculating the convoluted histogram.\n"); | |
480 | ||
481 | TH2* residual = (TH2*) convoluted->Clone("residual"); | |
482 | residual->SetTitle("(R#otimesUnfolded - Reconstructed)/Reconstructed;E^{jet} [GeV]; z^{lp}"); | |
483 | ||
484 | fCurrentRec->Scale(1./fCurrentRec->Integral()); | |
485 | residual->Add(fCurrentRec, -1); | |
486 | //residual->Divide(residual, fCurrentRec, 1, 1, "B"); | |
487 | ||
488 | // draw canvas | |
489 | TCanvas* canvas1 = new TCanvas(name, name, 1000, 1000); | |
490 | canvas1->Divide(2, 3); | |
491 | ||
492 | Int_t style = 1; | |
8de9091a | 493 | const Int_t nRGBs = 5; |
494 | const Int_t nCont = 500; | |
6d75bdb8 | 495 | |
8de9091a | 496 | Double_t stops[nRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 }; |
497 | Double_t red[nRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 }; | |
498 | Double_t green[nRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 }; | |
499 | Double_t blue[nRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 }; | |
500 | TColor::CreateGradientColorTable(nRGBs, stops, red, green, blue, nCont); | |
501 | gStyle->SetNumberContours(nCont); | |
6d75bdb8 | 502 | |
503 | canvas1->cd(1); | |
504 | gStyle->SetPalette(style); | |
505 | gPad->SetLogz(); | |
506 | genHist->SetTitle("Generated Spectrum;E^{jet}_{gen} [GeV];z^{lp}"); | |
507 | genHist->SetStats(0); | |
508 | genHist->DrawCopy("colz"); | |
509 | ||
510 | canvas1->cd(2); | |
511 | gStyle->SetPalette(style); | |
512 | gPad->SetLogz(); | |
513 | fUnfSpectrum->SetStats(0); | |
514 | fUnfSpectrum->DrawCopy("colz"); | |
515 | ||
516 | canvas1->cd(3); | |
517 | gStyle->SetPalette(style); | |
518 | gPad->SetLogz(); | |
519 | fCurrentRec->SetTitle(fRecSpectrum->GetTitle()); | |
520 | fCurrentRec->SetStats(0); | |
521 | fCurrentRec->DrawCopy("colz"); | |
522 | ||
523 | canvas1->cd(4); | |
524 | gStyle->SetPalette(style); | |
525 | gPad->SetLogy(); | |
526 | TH1D* projGenX = genHist->ProjectionX(); | |
527 | projGenX->SetTitle("Projection: Jet Energy; E^{jet} [GeV]"); | |
528 | TH1D* projUnfX = fUnfSpectrum->ProjectionX(); | |
529 | TH1D* projRecX = fCurrentRec->ProjectionX(); | |
530 | projGenX->SetStats(0); | |
531 | projRecX->SetStats(0); | |
532 | projUnfX->SetStats(0); | |
533 | projGenX->SetLineColor(8); | |
534 | projRecX->SetLineColor(2); | |
535 | projGenX->DrawCopy(); | |
536 | projUnfX->DrawCopy("same"); | |
537 | projRecX->DrawCopy("same"); | |
538 | ||
539 | TLegend* legend = new TLegend(0.6, 0.85, 0.98, 0.98); | |
540 | legend->AddEntry(projGenX, "Generated Spectrum"); | |
541 | legend->AddEntry(projUnfX, "Unfolded Spectrum"); | |
542 | legend->AddEntry(projRecX, "Reconstructed Spectrum"); | |
543 | //legend->SetFillColor(0); | |
544 | legend->Draw("same"); | |
545 | ||
546 | canvas1->cd(5); | |
547 | gPad->SetLogy(); | |
548 | gStyle->SetPalette(style); | |
549 | TH1D* projGenY = genHist->ProjectionY(); | |
550 | projGenY->SetTitle("Projection: Leading Particle Fragmentation; z^{lp}"); | |
551 | TH1D* projUnfY = fUnfSpectrum->ProjectionY(); | |
552 | TH1D* projRecY = fCurrentRec->ProjectionY(); | |
553 | projGenY->SetStats(0); | |
554 | projRecY->SetStats(0); | |
555 | projUnfY->SetStats(0); | |
556 | projGenY->SetLineColor(8); | |
557 | projRecY->SetLineColor(2); | |
558 | projGenY->DrawCopy(); | |
559 | projUnfY->DrawCopy("same"); | |
560 | projRecY->DrawCopy("same"); | |
561 | ||
562 | TLegend* legend1 = new TLegend(0.6, 0.85, 0.98, 0.98); | |
563 | legend1->AddEntry(projGenY, "Generated Spectrum"); | |
564 | legend1->AddEntry(projUnfY, "Unfolded Spectrum"); | |
565 | legend1->AddEntry(projRecY, "Recontructed Spectrum"); | |
566 | //legend1->SetFillColor(0); | |
567 | legend1->Draw("same"); | |
568 | ||
569 | // Draw residuals | |
570 | canvas1->cd(6); | |
571 | gStyle->SetPalette(style); | |
572 | gPad->SetLogz(); | |
573 | residual->SetStats(0); | |
574 | residual->DrawCopy("colz"); | |
575 | ||
576 | canvas1->SaveAs(Form("%s.png", canvas1->GetName())); | |
577 | } | |
578 | ||
579 | ||
6d75bdb8 | 580 | //____________________________________________________________________ |
8de9091a | 581 | void AliJetSpectrumUnfolding::ApplyBayesianMethod(Float_t regPar, Int_t nIterations, TH2* const initialConditions, Bool_t determineError) |
6d75bdb8 | 582 | { |
583 | // | |
584 | // correct spectrum using bayesian unfolding | |
585 | // | |
586 | ||
587 | // initialize seed with current time | |
588 | gRandom->SetSeed(0); | |
589 | ||
590 | printf("seting up current arrays and histograms...\n"); | |
591 | SetupCurrentHists(kFALSE); // kFALSE to not create big bin | |
592 | ||
593 | // normalize Correlation Map to convert number of events into probabilities | |
594 | /*for (Int_t te=1; te<=NBINSE; te++) | |
595 | for (Int_t tz=1; tz<=NBINSZ; tz++) | |
596 | { | |
597 | Int_t bin[4]; | |
598 | Float_t sum=0.; | |
599 | for (Int_t me = 1; me<=NBINSE; me++) | |
600 | for (Int_t mz = 1; mz<=NBINSZ; mz++) | |
601 | { | |
602 | bin[0] = te; bin[1] = me; | |
603 | bin[2] = tz; bin[3] = mz; | |
604 | sum += fCurrentCorrelation->GetBinContent(bin); | |
605 | } | |
606 | if (sum > 0.) | |
607 | for (Int_t me = 1; me<=NBINSE; me++) | |
608 | for (Int_t mz = 1; mz<=NBINSZ; mz++) | |
609 | { | |
610 | bin[0] = te; bin[1] = me; | |
611 | bin[2] = tz; bin[3] = mz; | |
612 | fCurrentCorrelation->SetBinContent(bin, fCurrentCorrelation->GetBinContent(bin)/sum); | |
613 | fCurrentCorrelation->SetBinError(bin, fCurrentCorrelation->GetBinError(bin)/sum); | |
614 | } | |
615 | }*/ | |
8de9091a | 616 | Float_t sum[fgkNBINSE+2][fgkNBINSZ+2]; |
617 | memset(sum,0,sizeof(Float_t)*(fgkNBINSE+2)*(fgkNBINSZ+2)); | |
6d75bdb8 | 618 | |
619 | for (Int_t idx=0; idx<=fCurrentCorrelation->GetNbins(); idx++) | |
620 | { | |
621 | Int_t bin[4]; | |
622 | Float_t binContent = fCurrentCorrelation->GetBinContent(idx, bin); | |
8de9091a | 623 | if ( (bin[1]>0 && bin[1]<=fgkNBINSE) && (bin[3]>0 && bin[3]<=fgkNBINSZ) ) |
6d75bdb8 | 624 | sum[bin[0]][bin[2]] += binContent; |
625 | } | |
626 | ||
627 | for (Int_t idx=0; idx<=fCurrentCorrelation->GetNbins(); idx++) | |
628 | { | |
629 | Int_t bin[4]; | |
630 | Float_t binContent = fCurrentCorrelation->GetBinContent(idx, bin); | |
631 | Float_t binError = fCurrentCorrelation->GetBinError(bin); | |
8de9091a | 632 | if (sum[bin[0]][bin[2]]>0 && (bin[1]>0 && bin[1]<=fgkNBINSE) && |
633 | (bin[3]>0 && bin[3]<=fgkNBINSZ) && (bin[0]>0 && bin[0]<=fgkNBINSE) && (bin[2]>0 && bin[2]<=fgkNBINSZ) ) | |
6d75bdb8 | 634 | { |
635 | fCurrentCorrelation->SetBinContent(bin, binContent/sum[bin[0]][bin[2]]); | |
636 | fCurrentCorrelation->SetBinError(bin, binError/sum[bin[0]][bin[2]]); | |
637 | } | |
638 | } | |
639 | ||
640 | printf("calling UnfoldWithBayesian\n"); | |
641 | Int_t success = UnfoldWithBayesian(fCurrentCorrelation, fCurrentRec, initialConditions, fUnfSpectrum, regPar, nIterations, kFALSE); | |
642 | ||
643 | if ( success != 0) | |
644 | return; | |
645 | ||
646 | if (!determineError) | |
647 | { | |
648 | Printf("AliJetSpectrumUnfolding::ApplyBayesianMethod: WARNING: No errors calculated."); | |
649 | return; | |
650 | } | |
651 | ||
652 | // evaluate errors, this is done by randomizing the measured spectrum following Poission statistics | |
653 | // this (new) measured spectrum is then unfolded and the different to the result from the "real" measured | |
654 | // spectrum calculated. This is performed N times and the maximum difference is taken as the systematic | |
655 | // error of the unfolding method itself. | |
656 | ||
657 | TH2* maxError = (TH2*) fUnfSpectrum->Clone("maxError"); | |
658 | maxError->Reset(); | |
659 | ||
660 | TH2* normalizedResult = (TH2*) fUnfSpectrum->Clone("normalizedResult"); | |
661 | normalizedResult->Scale(1.0 / normalizedResult->Integral()); | |
662 | ||
663 | const Int_t kErrorIterations = 20; | |
664 | ||
665 | printf("Spectrum unfolded. Determining error (%d iterations)...\n", kErrorIterations); | |
666 | ||
667 | TH2* randomized = (TH2*) fCurrentRec->Clone("randomized"); | |
668 | TH2* result2 = (TH2*) fUnfSpectrum->Clone("result2"); | |
669 | for (Int_t n=0; n<kErrorIterations; ++n) | |
670 | { | |
671 | // randomize the content of clone following a poisson with the mean = the value of that bin | |
672 | for (Int_t x=1; x<=randomized->GetNbinsX(); x++) | |
673 | for (Int_t y=1; y<=randomized->GetNbinsY(); y++) | |
674 | { | |
675 | Float_t randomValue = fCurrentRec->GetBinContent(x,y); | |
676 | TF1* poisson = new TF1("poisson", "TMath::Poisson(x,[0])",randomValue*0.25, randomValue*1.25); | |
677 | poisson->SetParameters(randomValue,0.); | |
678 | randomValue = poisson->GetRandom(); | |
679 | //printf("%e --> %e\n", fCurrentRec->GetBinContent(x,y), (Double_t)randomValue); | |
680 | randomized->SetBinContent(x, y, randomValue); | |
681 | delete poisson; | |
682 | } | |
683 | ||
684 | result2->Reset(); | |
685 | if (UnfoldWithBayesian(fCurrentCorrelation, randomized, initialConditions, result2, regPar, nIterations) != 0) | |
686 | return; | |
687 | ||
688 | result2->Scale(1.0 / result2->Integral()); | |
689 | ||
690 | // calculate ratio | |
691 | result2->Divide(normalizedResult); | |
692 | ||
693 | //new TCanvas; result2->DrawCopy("HIST"); | |
694 | ||
695 | // find max. deviation | |
696 | for (Int_t i=1; i<=result2->GetNbinsX(); i++) | |
697 | for (Int_t j=1; j<=result2->GetNbinsY(); j++) | |
698 | maxError->SetBinContent(i, j, TMath::Max(maxError->GetBinContent(i,j), TMath::Abs(1 - result2->GetBinContent(i,j)))); | |
699 | } | |
700 | delete randomized; | |
701 | delete result2; | |
702 | ||
703 | for (Int_t i=1; i<=fUnfSpectrum->GetNbinsX(); i++) | |
704 | for (Int_t j=1; j<=fUnfSpectrum->GetNbinsY(); j++) | |
705 | fUnfSpectrum->SetBinError(i, j, fUnfSpectrum->GetBinError(i,j) + maxError->GetBinContent(i,j)*fUnfSpectrum->GetBinContent(i,j)); | |
706 | ||
707 | delete maxError; | |
708 | delete normalizedResult; | |
709 | } | |
710 | ||
711 | //____________________________________________________________________ | |
8de9091a | 712 | Int_t AliJetSpectrumUnfolding::UnfoldWithBayesian(THnSparseF* const correlation, TH2* const measured, TH2* const initialConditions, TH2* const aResult, Float_t regPar, Int_t nIterations, Bool_t calculateErrors) |
6d75bdb8 | 713 | { |
714 | // | |
715 | // unfolds a spectrum | |
716 | // | |
717 | ||
718 | if (measured->Integral() <= 0) | |
719 | { | |
720 | Printf("AliJetSpectrumUnfolding::UnfoldWithBayesian: ERROR: The measured spectrum is empty"); | |
721 | return 1; | |
722 | } | |
8de9091a | 723 | const Int_t nFillesBins = correlation->GetNbins(); |
6d75bdb8 | 724 | const Int_t kStartBin = 1; |
725 | ||
8de9091a | 726 | const Int_t kMaxTZ = fgkNBINSZ; // max true axis fragmentation function |
727 | const Int_t kMaxMZ = fgkNBINSZ; // max measured axis fragmentation function | |
728 | const Int_t kMaxTE = fgkNBINSE; // max true axis energy | |
729 | const Int_t kMaxME = fgkNBINSE; // max measured axis energy | |
6d75bdb8 | 730 | |
8de9091a | 731 | printf("NbinsE=%d - NbinsZ=%d\n", fgkNBINSE, fgkNBINSZ); |
6d75bdb8 | 732 | |
733 | // store information in arrays, to increase processing speed | |
734 | Double_t measuredCopy[kMaxME+1][kMaxMZ+1]; | |
735 | Double_t prior[kMaxTE+1][kMaxTZ+1]; | |
736 | Double_t errors[kMaxTE+1][kMaxTZ+1]; | |
737 | Double_t result[kMaxTE+1][kMaxTZ+1]; | |
738 | ||
739 | THnSparseF *inverseCorrelation; | |
740 | inverseCorrelation = (THnSparseF*)correlation->Clone("inverseCorrelation"); | |
741 | inverseCorrelation->Reset(); | |
742 | ||
743 | Float_t inputDistIntegral = 1; | |
744 | if (initialConditions) | |
745 | { | |
746 | printf("Using different starting conditions...\n"); | |
747 | inputDistIntegral = initialConditions->Integral(); | |
748 | } | |
749 | Float_t measuredIntegral = measured->Integral(); | |
750 | for (Int_t me=1; me<=kMaxME; me++) | |
751 | for (Int_t mz=1; mz<=kMaxMZ; mz++) | |
752 | { | |
753 | // normalization of the measured spectrum | |
754 | measuredCopy[me][mz] = measured->GetBinContent(me,mz) / measuredIntegral; | |
755 | errors[me][mz] = measured->GetBinError(me, mz) / measuredIntegral; | |
756 | // pick prior distribution and normalize it | |
757 | if (initialConditions) | |
758 | prior[me][mz] = initialConditions->GetBinContent(me,mz) / inputDistIntegral; | |
759 | else | |
760 | prior[me][mz] = measured->GetBinContent(me,mz) / measuredIntegral; | |
761 | } | |
762 | ||
763 | // unfold... | |
764 | for (Int_t i=0; i<nIterations; i++) | |
765 | { | |
766 | // calculate Inverse Correlation Map from Bayes theorem: | |
767 | // IR_ji = R_ij * prior_i / sum_k(R_kj * prior_k) | |
768 | /*Float_t norm = 0; | |
769 | for (Int_t me=1; me<=kMaxME; me++) | |
770 | for (Int_t mz=1; mz<=kMaxMZ; mz++) | |
771 | { | |
772 | norm = 0; | |
773 | for (Int_t te=kStartBin; te<=kMaxTE; te++) | |
774 | for (Int_t tz=kStartBin; tz<=kMaxTZ; tz++) | |
775 | { | |
776 | Int_t bin[4] = {te, me, tz, mz}; | |
777 | norm += correlation->GetBinContent(bin)*prior[te][tz]; | |
778 | } | |
779 | if (norm > 0) | |
780 | for (Int_t te = kStartBin; te <= kMaxTE; te++) | |
781 | for (Int_t tz = kStartBin; tz <= kMaxTZ; tz++) | |
782 | { | |
783 | Int_t bin[4] = {te, me, tz, mz}; | |
784 | inverseCorrelation->SetBinContent(bin, correlation->GetBinContent(bin)*prior[te][tz]/norm ); | |
785 | } | |
786 | //else | |
787 | // inverse response set to '0' wich has been already done in line 2069 | |
788 | }*/ | |
789 | inverseCorrelation->Reset(); | |
790 | Float_t norm[kMaxTE+2][kMaxTZ+2]; | |
791 | for (Int_t te=0; te<(kMaxTE+2); te++) | |
792 | for (Int_t tz=0; tz<(kMaxTZ+2); tz++) | |
793 | norm[te][tz]=0; | |
794 | for (Int_t idx=0; idx<=correlation->GetNbins(); idx++) | |
795 | { | |
796 | Int_t bin[4]; | |
797 | Float_t binContent = correlation->GetBinContent(idx, bin); | |
8de9091a | 798 | if (bin[1]>0 && bin[1]<=fgkNBINSE && bin[3]>0 && bin[3]<=fgkNBINSZ && |
799 | bin[0]>0 && bin[0]<=fgkNBINSE && bin[2]>0 && bin[2]<=fgkNBINSZ) | |
6d75bdb8 | 800 | norm[bin[1]][bin[3]] += binContent*prior[bin[0]][bin[2]]; |
801 | } | |
802 | Float_t chi2Measured=0, diff; | |
803 | for (Int_t idx=0; idx<=correlation->GetNbins(); idx++) | |
804 | { | |
805 | Int_t bin[4]; | |
806 | Float_t binContent = correlation->GetBinContent(idx, bin); | |
8de9091a | 807 | if (norm[bin[1]][bin[3]]>0 && bin[1]>0 && bin[1]<=fgkNBINSE && |
808 | bin[3]>0 && bin[3]<=fgkNBINSZ && bin[0]>0 && bin[2]>0 && bin[0]<=fgkNBINSE && bin[2]<=fgkNBINSZ) | |
6d75bdb8 | 809 | { |
810 | inverseCorrelation->SetBinContent(bin, binContent*prior[bin[0]][bin[2]]/norm[bin[1]][bin[3]]); | |
811 | if (errors[bin[1]][bin[3]]>0) | |
812 | { | |
813 | diff = ((measuredCopy[bin[1]][bin[3]]-norm[bin[1]][bin[3]])/(errors[bin[1]][bin[3]])); | |
814 | chi2Measured += diff*diff; | |
815 | } | |
816 | } | |
817 | } | |
818 | ||
819 | // calculate "generated" spectrum | |
820 | for (Int_t te = kStartBin; te<=kMaxTE; te++) | |
821 | for (Int_t tz = kStartBin; tz<=kMaxTZ; tz++) | |
822 | { | |
823 | Float_t value = 0; | |
824 | for (Int_t me=1; me<=kMaxME; me++) | |
825 | for (Int_t mz=1; mz<=kMaxMZ; mz++) | |
826 | { | |
827 | Int_t bin[4] = {te, me, tz, mz}; | |
828 | value += inverseCorrelation->GetBinContent(bin)*measuredCopy[me][mz]; | |
829 | } | |
830 | result[te][tz] = value; | |
831 | //printf("%e\n", result[te][tz]); | |
832 | } | |
833 | ||
834 | // regularization (simple smoothing) | |
835 | Float_t chi2LastIter = 0; | |
836 | for (Int_t te=kStartBin; te<=kMaxTE; te++) | |
837 | for (Int_t tz=kStartBin; tz<=kMaxTZ; tz++) | |
838 | { | |
839 | Float_t newValue = 0; | |
840 | // 0 bin excluded from smoothing | |
841 | if (( te >(kStartBin+1) && te<(kMaxTE-1) ) && ( tz > (kStartBin+1) && tz<(kMaxTZ-1) )) | |
842 | { | |
843 | Float_t average = ((result[te-1][tz-1] + result[te-1][tz] + result[te-1][tz+1])+(result[te][tz-1] + result[te][tz] + result[te][tz+1])+(result[te+1][tz-1] + result[te+1][tz] + result[te+1][tz+1]))/9.; | |
844 | ||
845 | // weight the average with the regularization parameter | |
846 | newValue = (1 - regPar) * result[te][tz] + regPar * average; | |
847 | } | |
848 | else | |
849 | newValue = result[te][tz]; | |
850 | if (prior[te][tz]>1.e-5) | |
851 | { | |
852 | diff = ((prior[te][tz]-newValue)/prior[te][tz]); | |
853 | chi2LastIter = diff*diff; | |
854 | } | |
855 | prior[te][tz] = newValue; | |
856 | } | |
857 | //printf(" iteration %d - chi2LastIter = %e - chi2Measured = %e \n", i, chi2LastIter/((Float_t)kMaxTE*(Float_t)kMaxTZ), chi2Measured/((Float_t)kMaxTE*(Float_t)kMaxTZ)); | |
858 | if (chi2LastIter/((Float_t)kMaxTE*(Float_t)kMaxTZ)<5.e-6 && chi2Measured/((Float_t)kMaxTE*(Float_t)kMaxTZ)<5.e-3) | |
859 | break; | |
860 | } // end of iterations | |
861 | ||
862 | // propagate errors of the reconstructed distribution through the unfolding | |
863 | for (Int_t te = kStartBin; te<=kMaxTE; te++) | |
864 | for (Int_t tz = kStartBin; tz<=kMaxTZ; tz++) | |
865 | { | |
866 | Float_t valueError = 0; | |
8de9091a | 867 | // Float_t binError = 0; |
6d75bdb8 | 868 | for (Int_t me=1; me<=kMaxME; me++) |
869 | for (Int_t mz=1; mz<=kMaxMZ; mz++) | |
870 | { | |
871 | Int_t bin[4] = {te, me, tz, mz}; | |
872 | valueError += inverseCorrelation->GetBinContent(bin)*inverseCorrelation->GetBinContent(bin)*errors[me][mz]*errors[me][mz]; | |
873 | } | |
874 | //if (errors[te][tz]!=0)printf("errors[%d][%d]=%e\n", te, tz, valueError); | |
875 | aResult->SetBinContent(te, tz, prior[te][tz]); | |
876 | aResult->SetBinError(te, tz, TMath::Sqrt(valueError)); | |
877 | } | |
878 | ||
879 | // *********************************************************************************************************** | |
880 | // Calculate the covariance matrix, all arguments are taken from G. D'Agostini (p.6-8) | |
881 | if (calculateErrors) | |
882 | { | |
883 | printf("Covariance matrix will be calculated... this will take a lot of time (>1 day) ;)\n"); | |
884 | ||
885 | //Variables and Matrices that will be use along the calculation | |
8de9091a | 886 | const Int_t binsV[4] = {fgkNBINSE,fgkNBINSE, fgkNBINSZ, fgkNBINSZ}; |
887 | const Double_t lowEdgeV[4] = {fgkaxisLowerLimitE, fgkaxisLowerLimitE, fgkaxisLowerLimitZ, fgkaxisLowerLimitZ}; | |
888 | const Double_t upEdgeV[4] = {fgkaxisUpperLimitE, fgkaxisUpperLimitE, fgkaxisUpperLimitZ, fgkaxisUpperLimitZ}; | |
6d75bdb8 | 889 | |
8de9091a | 890 | const Double_t nTrue = (Double_t)measured->Integral(); |
6d75bdb8 | 891 | |
8de9091a | 892 | THnSparseF *v = new THnSparseF("V","",4, binsV, lowEdgeV, upEdgeV); |
893 | v->Reset(); | |
894 | Double_t invCorrContent1, nt; | |
895 | Double_t invCorrContent2, v11,v12, v2; | |
6d75bdb8 | 896 | // calculate V1 and V2 |
8de9091a | 897 | for (Int_t idx1=0; idx1<=nFillesBins; idx1++) |
6d75bdb8 | 898 | { |
8de9091a | 899 | printf("Covariance Matrix calculation: iteration idx1=%d of %d\n", idx1, nFillesBins); |
900 | for (Int_t idx2=0; idx2<=nFillesBins; idx2++) | |
6d75bdb8 | 901 | { |
902 | Int_t bin1[4]; | |
903 | Int_t bin2[4]; | |
904 | invCorrContent1 = inverseCorrelation->GetBinContent(idx1, bin1); | |
905 | invCorrContent2 = inverseCorrelation->GetBinContent(idx2, bin2); | |
906 | v11=0; v12=0; v2=0; | |
8de9091a | 907 | if(bin1[0]>0 && bin1[0]<=fgkNBINSE && bin1[1]>0 && bin1[1]<=fgkNBINSE && |
908 | bin1[2]>0 && bin1[2]<=fgkNBINSZ && bin1[3]>0 && bin1[3]<=fgkNBINSZ && | |
909 | bin2[0]>0 && bin2[0]<=fgkNBINSE && bin2[1]>0 && bin2[1]<=fgkNBINSE && | |
910 | bin2[2]>0 && bin2[2]<=fgkNBINSZ && bin2[3]>0 && bin2[3]<=fgkNBINSZ) | |
6d75bdb8 | 911 | { |
912 | if (bin1[1]==bin2[1] && bin1[3]==bin2[3]) | |
913 | v11 = invCorrContent1*invCorrContent2*measuredCopy[bin1[1]][bin1[3]] | |
8de9091a | 914 | *(1. - measuredCopy[bin2[1]][bin2[3]]/nTrue); |
6d75bdb8 | 915 | else |
916 | v12 = invCorrContent1*invCorrContent2*measuredCopy[bin1[1]][bin1[3]]* | |
8de9091a | 917 | measuredCopy[bin2[1]][bin2[3]]/nTrue; |
918 | nt = (Double_t)prior[bin2[0]][bin2[2]]; | |
6d75bdb8 | 919 | v2 = measuredCopy[bin1[1]][bin1[3]]*measuredCopy[bin2[1]][bin2[3]]* |
920 | invCorrContent1*invCorrContent2* | |
8de9091a | 921 | BayesUncertaintyTerms(inverseCorrelation, correlation, bin1, bin2, nt); |
6d75bdb8 | 922 | Int_t binV[4] = {bin1[0],bin2[0],bin1[2],bin2[2]}; |
8de9091a | 923 | v->SetBinContent(binV,v11-v12 + v2); |
6d75bdb8 | 924 | } |
925 | } | |
926 | } | |
927 | ||
8de9091a | 928 | for(Int_t te = 1; te<=fgkNBINSE; te++) |
929 | for(Int_t tz = 1; tz<=fgkNBINSZ; tz++) | |
6d75bdb8 | 930 | { |
931 | Int_t binV[4] = {te,te,tz,tz}; | |
8de9091a | 932 | aResult->SetBinError( te, tz, v->GetBinContent(binV) ); |
6d75bdb8 | 933 | } |
934 | ||
935 | TFile* f = new TFile("Covariance_UnfSpectrum.root"); | |
936 | f->Open("RECREATE"); | |
8de9091a | 937 | v->Write(); |
6d75bdb8 | 938 | f->Close(); |
939 | } | |
940 | ||
941 | return 0; | |
942 | ||
943 | } | |
944 | ||
945 | //____________________________________________________________________ | |
22ead7c4 | 946 | Double_t AliJetSpectrumUnfolding::BayesUncertaintyTerms(THnSparseF* const M, THnSparseF* const C,const Int_t* const binTM, const Int_t* const binTM1, Double_t nt) |
6d75bdb8 | 947 | { |
948 | // | |
949 | // helper function for the covariance matrix of the bayesian method | |
950 | // | |
951 | ||
952 | Double_t result = 0; | |
85587ff2 | 953 | Float_t term[9] = {0.}; |
954 | Int_t tmpBin[4] = {0}, tmpBin1[4] = {0}; | |
6d75bdb8 | 955 | const Int_t nFilledBins = C->GetNbins(); |
22ead7c4 | 956 | if (!(nt>0&&nt<0)) |
6d75bdb8 | 957 | return 0; |
958 | ||
8de9091a | 959 | Float_t corrContent; |
960 | Float_t invCorrContent; | |
6d75bdb8 | 961 | |
962 | tmpBin[0] =binTM[0]; tmpBin[1] =binTM[1]; tmpBin[2] =binTM[2]; tmpBin[3] =binTM[3]; | |
963 | tmpBin1[0]=binTM[0]; tmpBin1[1]=binTM1[1]; tmpBin1[2]=binTM[2]; tmpBin1[3]=binTM1[3]; | |
964 | if (C->GetBinContent(tmpBin)!=0 && C->GetBinContent(tmpBin1)!=0) | |
965 | { | |
966 | if (binTM[0]==binTM1[0] && binTM[2]==binTM1[2]) | |
967 | term[0] = BayesCov(M, C, tmpBin, tmpBin1)/ | |
968 | (C->GetBinContent(tmpBin)*C->GetBinContent(tmpBin1)); | |
969 | term[2] = term[0]*M->GetBinContent(tmpBin1); | |
970 | } | |
971 | else | |
972 | { | |
973 | term[0] = 0; | |
974 | term[2] = 0; | |
975 | } | |
976 | ||
977 | tmpBin[0]=binTM1[0]; tmpBin[1]=binTM[1]; tmpBin[2]=binTM1[2]; tmpBin[3]=binTM[3]; | |
978 | tmpBin1[0]=binTM1[0]; tmpBin1[1]=binTM1[1]; tmpBin1[2]=binTM1[2]; tmpBin1[3]=binTM1[3]; | |
979 | if (C->GetBinContent(tmpBin)!=0 && C->GetBinContent(tmpBin1)!=0) | |
980 | term[6] = BayesCov(M, C, tmpBin, tmpBin1)* | |
981 | M->GetBinContent(tmpBin)/ | |
982 | (C->GetBinContent(tmpBin)*C->GetBinContent(tmpBin1)); | |
983 | else | |
984 | term[6] = 0; | |
985 | ||
986 | for(Int_t idx1=0; idx1<=nFilledBins; idx1++) | |
987 | { | |
988 | Int_t bin1[4]; | |
8de9091a | 989 | corrContent = C->GetBinContent(idx1, bin1); |
990 | invCorrContent = M->GetBinContent(idx1, bin1); | |
991 | if(bin1[0]>0 && bin1[0]<=fgkNBINSE && bin1[1]>0 && bin1[1]<=fgkNBINSE && | |
992 | bin1[2]>0 && bin1[2]<=fgkNBINSZ && bin1[3]>0 && bin1[3]<=fgkNBINSZ) | |
6d75bdb8 | 993 | { |
994 | tmpBin[0] =binTM[0]; tmpBin[1] =binTM[1]; tmpBin[2] =binTM[2]; tmpBin[3] =binTM[3]; | |
995 | tmpBin1[0]=binTM[0]; tmpBin1[1]=bin1[1]; tmpBin1[2]=binTM[2]; tmpBin1[3]=bin1[3]; | |
996 | if (C->GetBinContent(tmpBin)!=0 && | |
997 | binTM[0]==binTM1[0] && binTM[2]==binTM1[2]) | |
998 | term[1] = BayesCov(M, C, tmpBin, tmpBin1)/C->GetBinContent(tmpBin); | |
999 | else | |
1000 | term[1] = 0; | |
1001 | ||
1002 | tmpBin[0] =binTM[0]; tmpBin[1] =bin1[1]; tmpBin[2] =binTM[2]; tmpBin[3] =bin1[3]; | |
1003 | tmpBin1[0]=binTM[0]; tmpBin1[1]=binTM1[1]; tmpBin1[2]=binTM[2]; tmpBin1[3]=binTM1[3]; | |
1004 | if (C->GetBinContent(tmpBin1)!=0) | |
1005 | { | |
1006 | if (binTM[0]==binTM1[0] && binTM[2]==binTM1[2]) | |
1007 | term[3] = BayesCov(M, C, tmpBin, tmpBin1)/ | |
1008 | C->GetBinContent(tmpBin1); | |
1009 | term[5] = BayesCov(M, C, tmpBin, tmpBin1)*M->GetBinContent(tmpBin1)/ | |
1010 | C->GetBinContent(tmpBin1); | |
1011 | } | |
1012 | else | |
1013 | { | |
1014 | term[3] = 0; | |
1015 | term[5] = 0; | |
1016 | } | |
1017 | ||
1018 | tmpBin[0] =binTM1[0]; tmpBin[1] =binTM[1]; tmpBin[2] =binTM1[2]; tmpBin[3] =binTM[3]; | |
1019 | tmpBin1[0]=binTM1[0]; tmpBin1[1]=bin1[1]; tmpBin1[2]=binTM1[2]; tmpBin1[3]=bin1[3]; | |
1020 | if (C->GetBinContent(tmpBin)!=0) | |
1021 | term[7] = BayesCov(M, C, tmpBin, tmpBin1)*M->GetBinContent(tmpBin)/ | |
1022 | C->GetBinContent(tmpBin); | |
1023 | else | |
1024 | term[7] = 0; | |
1025 | ||
1026 | tmpBin[0] =bin1[0]; tmpBin[1] =binTM[1]; tmpBin[2] =bin1[2]; tmpBin[3] =binTM[3]; | |
1027 | tmpBin1[0]=bin1[0]; tmpBin1[1]=binTM1[1]; tmpBin1[2]=bin1[2]; tmpBin1[3]=binTM1[3]; | |
1028 | if (C->GetBinContent(tmpBin)!=0 && C->GetBinContent(tmpBin1)!=0) | |
1029 | term[8] = BayesCov(M, C, tmpBin, tmpBin1)* | |
1030 | M->GetBinContent(tmpBin)*M->GetBinContent(tmpBin)/ | |
1031 | (C->GetBinContent(tmpBin)*C->GetBinContent(tmpBin1)); | |
1032 | else | |
1033 | term[8] = 0; | |
1034 | ||
1035 | for (Int_t i=0; i<9; i++) | |
8de9091a | 1036 | result += term[i]/nt; |
6d75bdb8 | 1037 | } |
1038 | } | |
1039 | ||
1040 | return result; | |
1041 | } | |
1042 | ||
1043 | //____________________________________________________________________ | |
22ead7c4 | 1044 | Double_t AliJetSpectrumUnfolding::BayesCov(THnSparseF* const M, THnSparseF* const correlation,const Int_t* const binTM,const Int_t* const bin1) |
6d75bdb8 | 1045 | { |
8de9091a | 1046 | |
1047 | // | |
1048 | // get the covariance matrix | |
1049 | // | |
1050 | ||
1051 | ||
6d75bdb8 | 1052 | Double_t result, result1, result2, result3; |
1053 | ||
1054 | if (binTM[0]==bin1[0] && binTM[2]==bin1[2]) | |
1055 | { | |
1056 | if (correlation->GetBinContent(bin1)!=0) | |
1057 | result1 = 1./correlation->GetBinContent(bin1); | |
1058 | else | |
1059 | result1 = 0; | |
1060 | result2 = 1.; | |
1061 | } | |
1062 | else | |
1063 | { | |
1064 | result1 = 0; | |
1065 | result2 = 0; | |
1066 | } | |
1067 | ||
1068 | if (binTM[1]==bin1[1] && binTM[3]==bin1[3]) | |
1069 | { | |
1070 | Int_t tmpbin[4] = {bin1[0], binTM[1], bin1[2], binTM[3]}; | |
1071 | if(correlation->GetBinContent(tmpbin)!=0) | |
1072 | result3 = M->GetBinContent(tmpbin)/correlation->GetBinContent(tmpbin); | |
1073 | else | |
1074 | result3 = 0; | |
1075 | } | |
1076 | else | |
1077 | { | |
1078 | result1 = 0; | |
1079 | result3 = 0; | |
1080 | } | |
1081 | ||
1082 | return result = result1 + result2 + result3; | |
1083 | } | |
1084 | ||
1085 | //____________________________________________________________________ | |
8de9091a | 1086 | TH2F* AliJetSpectrumUnfolding::CalculateRecSpectrum(TH2* const inputGen) |
6d75bdb8 | 1087 | { |
1088 | // runs the distribution given in inputGen through the correlation histogram identified by | |
1089 | // fCorrelation and produces a reconstructed spectrum | |
1090 | ||
1091 | if (!inputGen) | |
1092 | return 0; | |
1093 | ||
1094 | // normalize to convert number of events into probability | |
1095 | /*for (Int_t te=1; te<=NBINSE; te++) | |
1096 | for (Int_t tz=1; tz<=NBINSZ; tz++) | |
1097 | { | |
1098 | Int_t bin[4]; | |
1099 | Float_t sum=0.; | |
1100 | for (Int_t me = 1; me<=NBINSE; me++) | |
1101 | for (Int_t mz = 1; mz<=NBINSZ; mz++) | |
1102 | { | |
1103 | bin[0] = te; bin[1] = me; | |
1104 | bin[2] = tz; bin[3] = mz; | |
1105 | sum += fCorrelation[correlationMap]->GetBinContent(bin); | |
1106 | } | |
1107 | if (sum > 0.) | |
1108 | for (Int_t me = 1; me<=NBINSE; me++) | |
1109 | for (Int_t mz = 1; mz<=NBINSZ; mz++) | |
1110 | { | |
1111 | bin[0] = te; bin[1] = me; | |
1112 | bin[2] = tz; bin[3] = mz; | |
1113 | fCorrelation[correlationMap]->SetBinContent(bin, fCorrelation[correlationMap]->GetBinContent(bin)/sum); | |
1114 | fCorrelation[correlationMap]->SetBinError(bin, fCorrelation[correlationMap]->GetBinError(bin)/sum); | |
1115 | } | |
1116 | }*/ | |
1117 | // normalize to convert number of events into probability (the following loop is much faster) | |
8de9091a | 1118 | Float_t sum[fgkNBINSE+2][fgkNBINSZ+2]; |
1119 | memset(sum,0,sizeof(Float_t)*(fgkNBINSE+2)*(fgkNBINSZ+2)); | |
6d75bdb8 | 1120 | |
1121 | for (Int_t idx=0; idx<fCorrelation->GetNbins(); idx++) | |
1122 | { | |
1123 | Int_t bin[4]; | |
1124 | Float_t binContent = fCorrelation->GetBinContent(idx, bin); | |
8de9091a | 1125 | if (bin[1]>0 && bin[1]<=fgkNBINSE && bin[3]>0 && bin[3]<=fgkNBINSZ){ |
6d75bdb8 | 1126 | sum[bin[0]][bin[2]] += binContent; |
1127 | } | |
1128 | } | |
1129 | ||
1130 | for (Int_t idx=0; idx<fCorrelation->GetNbins(); idx++) | |
1131 | { | |
1132 | Int_t bin[4]; | |
1133 | Float_t binContent = fCorrelation->GetBinContent(idx, bin); | |
1134 | Float_t binError = fCorrelation->GetBinError(bin); | |
8de9091a | 1135 | if (sum[bin[0]][bin[2]]>0 && bin[1]>0 && bin[1]<=fgkNBINSE && |
1136 | bin[3]>0 && bin[3]<=fgkNBINSZ && bin[0]>0 && bin[2]>0 && bin[0]<=fgkNBINSE && bin[2]<=fgkNBINSZ) | |
6d75bdb8 | 1137 | { |
1138 | fCorrelation->SetBinContent(bin, binContent/sum[bin[0]][bin[2]]); | |
1139 | fCorrelation->SetBinError(bin, binError/sum[bin[0]][bin[2]]); | |
1140 | } | |
1141 | } | |
1142 | ||
1143 | TH2F* target = dynamic_cast<TH2F*> (fRecSpectrum->Clone(Form("reconstructed_%s", inputGen->GetName()))); | |
b85b355c | 1144 | if(!targer)return 0; |
1145 | ||
6d75bdb8 | 1146 | target->Reset(); |
1147 | ||
8de9091a | 1148 | for (Int_t me=1; me<=fgkNBINSE; ++me) |
1149 | for (Int_t mz=1; mz<=fgkNBINSZ; ++mz) | |
6d75bdb8 | 1150 | { |
1151 | Float_t measured = 0; | |
1152 | Float_t error = 0; | |
1153 | ||
8de9091a | 1154 | for (Int_t te=1; te<=fgkNBINSE; ++te) |
1155 | for (Int_t tz=1; tz<=fgkNBINSZ; ++tz) | |
6d75bdb8 | 1156 | { |
1157 | Int_t bin[4] = {te, me, tz, mz}; | |
1158 | measured += inputGen->GetBinContent(te,tz) * fCorrelation->GetBinContent(bin); | |
1159 | error += inputGen->GetBinError(te,tz) * fCorrelation->GetBinContent(bin); | |
1160 | } | |
1161 | target->SetBinContent(me, mz, measured); | |
1162 | target->SetBinError(me, mz, error); | |
1163 | } | |
1164 | ||
1165 | return target; | |
1166 | } | |
1167 | ||
1168 | //__________________________________________________________________________________________________ | |
22ead7c4 | 1169 | void AliJetSpectrumUnfolding::SetGenRecFromFunc(const TF2* const inputGen) |
6d75bdb8 | 1170 | { |
1171 | // uses the given function to fill the input Generated histogram and generates from that | |
1172 | // the reconstructed histogram by applying the response histogram | |
1173 | // this can be used to evaluate if the methods work indepedently of the input | |
1174 | // distribution | |
1175 | ||
1176 | if (!inputGen) | |
1177 | return; | |
1178 | ||
8de9091a | 1179 | TH2F* histtmp = new TH2F("histtmp", "tmp", |
1180 | fgkNBINSE, fgkaxisLowerLimitE, fgkaxisUpperLimitE, | |
1181 | fgkNBINSZ, fgkaxisLowerLimitZ, fgkaxisUpperLimitZ); | |
1182 | ||
6d75bdb8 | 1183 | TH2F* gen = fGenSpectrum; |
1184 | ||
1185 | histtmp->Reset(); | |
1186 | gen->Reset(); | |
1187 | ||
8de9091a | 1188 | histtmp->FillRandom(inputGen->GetName(), fgkNEVENTS); |
6d75bdb8 | 1189 | |
1190 | for (Int_t i=1; i<=gen->GetNbinsX(); ++i) | |
1191 | for (Int_t j=1; j<=gen->GetNbinsY(); ++j) | |
1192 | { | |
1193 | gen->SetBinContent(i, j, histtmp->GetBinContent(i,j)); | |
1194 | gen->SetBinError(i, j, histtmp->GetBinError(i,j)); | |
1195 | } | |
1196 | ||
1197 | delete histtmp; | |
1198 | ||
1199 | //new TCanvas; | |
1200 | //gStyle->SetPalette(1); | |
1201 | //gPad->SetLogz(); | |
1202 | //gen->Draw("COLZ"); | |
1203 | ||
1204 | ||
1205 | TH2 *recsave = fRecSpectrum; | |
1206 | ||
1207 | fRecSpectrum = CalculateRecSpectrum(gen); | |
1208 | fRecSpectrum->SetName(recsave->GetName()); | |
1209 | delete recsave; | |
1210 | ||
1211 | return; | |
1212 | } | |
1213 | //________________________________________________________________________________________ |