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