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