]>
Commit | Line | Data |
---|---|---|
6d75bdb8 | 1 | // |
2 | // script with functions to use AliJetSpectrumUnfolding class | |
3 | // | |
4 | ||
5 | ||
6 | void loadlibs(){ | |
7 | // load all the needed libs to run wit root only | |
8 | ||
9 | gSystem->Load("libTree"); | |
10 | gSystem->Load("libPhysics"); | |
11 | gSystem->Load("libHist"); | |
12 | gSystem->Load("libVMC"); | |
13 | gSystem->Load("libSTEERBase"); | |
14 | gSystem->Load("libAOD"); | |
15 | gSystem->Load("libESD"); | |
16 | gSystem->Load("libANALYSIS"); | |
17 | gSystem->Load("libANALYSISalice"); | |
18 | gSystem->Load("libJETAN"); | |
19 | gSystem->Load("libPWG4JetTasks"); | |
20 | ||
21 | ||
22 | } | |
23 | ||
24 | ||
25 | ||
26 | void draw(const char* fileName = "unfolded_pwg4spec.root", const char* folder = "unfolding", Bool_t proj = kTRUE) | |
27 | { | |
28 | ||
29 | loadlibs(); | |
30 | ||
31 | AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding(folder, folder); | |
32 | ||
33 | TFile::Open(fileName); | |
34 | jetSpec->LoadHistograms(); | |
35 | ||
36 | ||
37 | if (proj) | |
38 | { | |
39 | canvas1 = new TCanvas("Response Map Projection", "Response Map Projection", 500, 500); | |
40 | canvas1->Divide(2); | |
41 | ||
42 | Int_t style = 1; | |
43 | const Int_t NRGBs = 5; | |
44 | const Int_t NCont = 500; | |
45 | ||
46 | Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 }; | |
47 | Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 }; | |
48 | Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 }; | |
49 | Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 }; | |
50 | TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont); | |
51 | gStyle->SetNumberContours(NCont); | |
52 | gStyle->SetPalette(style); | |
53 | ||
54 | canvas1->cd(1); | |
55 | gPad->SetLogz(); | |
56 | h2 = (jetSpec->GetCorrelation())->Projection(2,3); | |
57 | h2->SetXTitle("z^{lp}_{gen}"); | |
58 | h2->SetYTitle("z^{lp}_{rec}"); | |
59 | h2->DrawCopy("colz"); | |
60 | ||
61 | canvas1->cd(2); | |
62 | gPad->SetLogz(); | |
63 | h3 = (jetSpec->GetCorrelation())->Projection(0,1); | |
64 | h3->SetXTitle("E^{jet}_{gen} [GeV]"); | |
65 | h3->SetYTitle("E^{jet}_{rec} [GeV]"); | |
66 | h3->DrawCopy("colz"); | |
67 | } | |
68 | jetSpec->DrawComparison("Draw_unfolded_pwg4spec", jetSpec->GetGenSpectrum()); | |
69 | ||
70 | return; | |
71 | } | |
72 | ||
73 | //________________________________________________________________________________________________________________ | |
74 | void unfold(const char* fileNameGen = "gen_pwg4spec.root", const char* folder = "unfolding", const char* fileNameRec = "rec_pwg4spec.root", const char* fileNameUnf = "unfolded_pwg4spec.root") | |
75 | { | |
76 | // function to load jet spectra from the output file of the task AliAnalysisTaskJetSpectrum | |
77 | // to do the unfolding | |
78 | ||
79 | loadlibs(); | |
80 | ||
81 | AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding(folder, folder); | |
82 | ||
df65bddb | 83 | TFile::Open(fileNameGen); |
6d75bdb8 | 84 | jetSpec->LoadHistograms(); |
85 | ||
df65bddb | 86 | TFile::Open(fileNameRec); |
87 | TH2F* hist = (TH2F*) gFile->Get("unfolding/fRecSpectrum"); | |
88 | jetSpec->SetRecSpectrum(hist); | |
4dbfdecc | 89 | hist = (TH2F*) gFile->Get("unfolding/fGenSpectrum"); |
90 | if(hist->GetEntries()>0)jetSpec->SetGenSpectrum(hist); | |
6d75bdb8 | 91 | |
92 | jetSpec->ApplyBayesianMethod(0.3, 20, 0, 0); | |
93 | // last parameter = calculateErrors <- this method to calculate the errors takes a lot of time | |
94 | ||
95 | TFile* file = TFile::Open(fileNameUnf,"RECREATE"); | |
96 | jetSpec->SaveHistograms(); | |
97 | file->Close(); | |
98 | } | |
99 | ||
100 | //___________________________________________________________________________ | |
101 | void buildSpectra(Int_t caseNo, const char* inFile, const char* outFile) | |
102 | { | |
103 | // function to test 2D Bayesian unfolding with other spectra | |
104 | // build from a function | |
105 | ||
106 | loadlibs(); | |
107 | ||
108 | ||
109 | ||
110 | AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding("unfolding", "unfolding"); | |
111 | ||
112 | TFile::Open(inFile); | |
113 | jetSpec->LoadHistograms("unfolding"); | |
114 | ||
115 | TCanvas *c1 = new TCanvas();TH2 *h2 = (jetSpec->GetCorrelation())->Projection(0,3); | |
116 | h2->DrawCopy("colz"); | |
117 | c1->Update(); | |
118 | if(getchar()=='q')return; | |
119 | ||
120 | ||
121 | switch (caseNo) | |
122 | { | |
123 | case 1: func = new TF2("func", "501-x-y"); break; | |
124 | case 2: func = new TF2("func", "1000 * 1/(y+x+1)"); break; | |
125 | case 3: func = new TF2("func", "1./(y*pow(x,5.7))"); break; | |
126 | case 4: func = new TF2("func", "exp(-0.1*x - 0.1*y)"); break; | |
127 | case 5: func = new TF2("func", "x*x + (y**3)/100."); break; | |
128 | case 6: func = new TF2("func", "1000*y*exp(-0.1*x)"); break; | |
129 | case 7: func = new TF2("func", "exp(-((x-100.)/(0.3*x))**2 - ((y-0.6)/(0.8*y))**2)"); break; | |
130 | default: return; | |
131 | } | |
132 | ||
133 | //new TCanvas; func->Draw(); | |
134 | ||
135 | jetSpec->SetGenRecFromFunc(func); | |
136 | ||
137 | TFile* file = TFile::Open(outFile,"RECREATE"); | |
138 | jetSpec->SaveHistograms(); | |
139 | ||
140 | h2 = (jetSpec->GetCorrelation())->Projection(0,3); | |
141 | h2->DrawCopy("colz"); | |
142 | c1->Update(); | |
143 | file->Close(); | |
144 | ||
145 | //new TCanvas; jetSpec->GetRecSpectrum()->DrawCopy("colz"); | |
146 | //new TCanvas; jetSpec->GetGenSpectrum()->DrawCopy("colz"); | |
147 | } | |
148 | ||
149 | //_____________________________________________________________________________________________ | |
150 | void buildResponseMap(const char* fileName = "responseMap.root") | |
151 | { | |
152 | // function to build a Response Map with a gaussian distribution | |
153 | loadlibs(); | |
154 | ||
155 | AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding("unfolding", "unfolding"); | |
156 | ||
157 | TF2* func = new TF2("func", "exp(-((x-[0])/[1])**2 - ((y-[2])/[3])**2)"); | |
158 | ||
159 | bool bPerfect = false; | |
160 | ||
161 | Double_t var[4]; | |
162 | Float_t sigmax, sigmay; | |
163 | for (Int_t tx=1; tx<=jetSpec->GetCorrelation()->GetAxis(0)->GetNbins(); tx++) | |
164 | for (Int_t ty=1; ty<=jetSpec->GetCorrelation()->GetAxis(2)->GetNbins(); ty++) | |
165 | { | |
166 | var[0] = jetSpec->GetCorrelation()->GetAxis(0)->GetBinCenter(tx); | |
167 | var[2] = jetSpec->GetCorrelation()->GetAxis(2)->GetBinCenter(ty); | |
168 | sigmax = 0.2*var[0]; | |
169 | sigmay = 0.2*var[2]; | |
170 | func->SetParameters(var[0],sigmax,var[2],sigmay); | |
171 | for (Int_t mx=1; mx<=jetSpec->GetCorrelation()->GetAxis(1)->GetNbins(); mx++) | |
172 | for (Int_t my=1; my<=jetSpec->GetCorrelation()->GetAxis(3)->GetNbins(); my++) | |
173 | { | |
174 | var[1] = jetSpec->GetCorrelation()->GetAxis(1)->GetBinCenter(mx); | |
175 | var[3] = jetSpec->GetCorrelation()->GetAxis(3)->GetBinCenter(my); | |
176 | ||
177 | ||
178 | if(bPerfect){ | |
179 | if (var[1]==var[0] && var[3]==var[2]) | |
180 | jetSpec->GetCorrelation()->Fill(var,1); | |
181 | } | |
182 | else { | |
183 | // cut at sigma | |
184 | if (TMath::Abs(var[1]-var[0]) < sigmax || TMath::Abs(var[3]-var[2]) < sigmay) | |
185 | jetSpec->GetCorrelation()->Fill(var,func->Eval(var[1],var[3])); | |
186 | } | |
187 | } | |
188 | } | |
189 | ||
190 | ||
191 | TFile* file = TFile::Open(fileName,"RECREATE"); | |
192 | jetSpec->SaveHistograms(); | |
193 | file->Close(); | |
194 | } | |
195 | ||
196 | //_____________________________________________________________________________ | |
197 | void StatisticalUncertainties(const char* fileNameGen = "gen_pwg4spec.root", const char* folder = "unfolding", const char* fileNameRec = "rec_pwg4spec.root", const char* fileNameOut = "Uncertainties2DBayesUnf.root") | |
198 | { | |
199 | // This function gives the statistical uncertainties due to the 2D Bayesian Unfolding | |
200 | ||
201 | gSystem->Load("libANALYSIS"); | |
202 | gSystem->Load("libANALYSISalice"); | |
203 | gSystem->Load("libJETAN"); | |
204 | gSystem->Load("libPWG4JetTasks"); | |
205 | ||
206 | AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding(folder, folder); | |
207 | ||
208 | TFile::Open(fileNameRec); | |
209 | jetSpec->LoadHistograms(); | |
210 | ||
211 | TFile::Open(fileNameGen); | |
212 | TH2F* hist = (TH2F*) gFile->Get("unfolding/fGenSpectrum"); | |
213 | jetSpec->SetGenSpectrum(hist); | |
214 | ||
215 | // create sigma histogram | |
216 | TH2F* sigma = (TH2F*)(jetSpec->GetGenSpectrum())->Clone("sigma"); | |
217 | sigma->Reset(); | |
218 | ||
219 | THnSparseF* hcorr = (THnSparseF*)(jetSpec->GetCorrelation())->Clone(); | |
220 | TH2F* htrue = (TH2F*)(jetSpec->GetGenSpectrum())->Clone(); | |
221 | TH2F* hmeas = (TH2F*)(jetSpec->GetRecSpectrum())->Clone(); | |
222 | TH2F* hunfo = (TH2F*)(jetSpec->GetUnfSpectrum())->Clone(); | |
223 | ||
224 | Int_t nIterations = 1000; | |
225 | Float_t binContent; | |
226 | for(Int_t i=0; i<nIterations; i++) | |
227 | { | |
228 | printf("iteration = %d\n", i); | |
229 | // reset histograms | |
230 | jetSpec->GetRecSpectrum()->Reset(); | |
231 | jetSpec->GetGenSpectrum()->Reset(); | |
232 | jetSpec->GetCorrelation()->Reset(); | |
233 | jetSpec->GetUnfSpectrum()->Reset(); | |
234 | ||
235 | THnSparseF* tmpcorr = (THnSparseF*)hcorr->Clone("tmpcorr"); | |
236 | TH2F* tmptrue = (TH2F*)htrue->Clone("tmptrue"); | |
237 | ||
238 | jetSpec->SetGenSpectrum(tmptrue); | |
239 | jetSpec->SetCorrelation(tmpcorr); | |
240 | ||
241 | // randomize reconstructed distribution | |
242 | for (Int_t me=1; me<=hmeas->GetNbinsX(); me++) | |
243 | for (Int_t mz=1; mz<=hmeas->GetNbinsY(); mz++) | |
244 | { | |
245 | binContent = hmeas->GetBinContent(me,mz); | |
246 | if (binContent>0) | |
247 | { | |
248 | TF1* poisson = new TF1("poisson", "TMath::Poisson(x,[0])",binContent*0.25, binContent*1.25); | |
249 | poisson->SetParameters(binContent,0.); | |
250 | binContent = poisson->GetRandom(); | |
251 | delete poisson; | |
252 | } | |
253 | jetSpec->GetRecSpectrum()->SetBinContent(me,mz, binContent); | |
254 | } | |
255 | ||
256 | // unfold | |
257 | jetSpec->ApplyBayesianMethod(0.2, 20, 0, 0); | |
258 | ||
259 | // calculate sigma^2 | |
260 | for (Int_t te=1; te<=sigma->GetNbinsX(); te++) | |
261 | for (Int_t tz=1; tz<=sigma->GetNbinsY(); tz++) | |
262 | { | |
263 | if (htrue->GetBinContent(te,tz)!=0) | |
264 | { | |
265 | binContent = (jetSpec->GetUnfSpectrum()->GetBinContent(te,tz) - | |
266 | htrue->GetBinContent(te,tz))/htrue->GetBinContent(te,tz); | |
267 | binContent *= binContent; | |
268 | sigma->SetBinContent(te,tz, binContent + sigma->GetBinContent(te,tz)); | |
269 | } | |
270 | } | |
271 | } | |
272 | ||
273 | // calculate sigma | |
274 | for (Int_t te=1; te<=sigma->GetNbinsX(); te++) | |
275 | for (Int_t tz=1; tz<=sigma->GetNbinsY(); tz++) | |
276 | { | |
277 | binContent = sigma->GetBinContent(te,tz); | |
278 | binContent = TMath::Sqrt(binContent/(Float_t)nIterations); | |
279 | sigma->SetBinContent(te,tz, binContent); | |
280 | } | |
281 | ||
282 | const Int_t NRGBs = 5; | |
283 | const Int_t NCont = 500; | |
284 | ||
285 | Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 }; | |
286 | Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 }; | |
287 | Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 }; | |
288 | Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 }; | |
289 | TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont); | |
290 | gStyle->SetNumberContours(NCont); | |
291 | gStyle->SetPalette(1); | |
292 | ||
293 | new TCanvas(); | |
294 | gPad->SetLogz(); | |
295 | sigma->SetTitle("#sigma((U_{R} - U)/U)"); | |
296 | sigma->SetXTitle("E^{jet} [GeV]"); | |
297 | sigma->SetYTitle("z^{lp}"); | |
298 | sigma->DrawCopy(); | |
299 | ||
300 | TFile* file = TFile::Open(fileNameOut,"RECREATE"); | |
301 | sigma->Write(); | |
302 | file->Close(); | |
303 | } | |
304 | ||
305 | //_______________________________________________________________________________________________________________ | |
df65bddb | 306 | void FillSpecFromFiles(const char* fileNameReal = "histos_pwg4spec.root",const char* fileNameSim = "histos_pwg4spec.root") |
6d75bdb8 | 307 | { |
308 | // This functions avoids problems when the number of bins or the bin limits | |
309 | // in the histograms of the AliJetSpectrumUnfolding and AliAnalysisTaskJetSpectrum classes | |
310 | // are different. | |
311 | ||
312 | gSystem->Load("libANALYSIS"); | |
313 | gSystem->Load("libANALYSISalice"); | |
314 | gSystem->Load("libJETAN"); | |
315 | gSystem->Load("libPWG4JetTasks"); | |
316 | ||
df65bddb | 317 | file = new TFile(fileNameSim,"read"); |
6d75bdb8 | 318 | tlist = dynamic_cast<TList*> (file->Get("pwg4spec")); |
6d75bdb8 | 319 | |
df65bddb | 320 | THnSparseF *fhCorrelation = 0; |
321 | for(int ic = 0;ic<3;++ic){ | |
322 | THnSparseF *hTmp = (THnSparseF*)(tlist->FindObject(Form("fhnCorrelation_%d",ic))); | |
323 | if(fhCorrelation==0)fhCorrelation = (THnSparseF*)hTmp->Clone("fhnCorrelation"); | |
324 | else fhCorrelation->Add(hTmp); | |
325 | } | |
326 | TH2F *fhEGenZGen = (TH2F*)(tlist->FindObject("fh2EGenZGen")); | |
6d75bdb8 | 327 | |
328 | AliJetSpectrumUnfolding *jetSpec = new AliJetSpectrumUnfolding("unfolding","unfolding"); | |
329 | ||
330 | // generated jets (true distribution) | |
331 | for (Int_t te=1; te<=fhEGenZGen->GetNbinsX(); te++) | |
332 | for (Int_t tz=1; tz<=fhEGenZGen->GetNbinsY(); tz++) | |
333 | { | |
334 | Float_t ej = fhEGenZGen->GetXaxis()->GetBinCenter(te); | |
335 | Float_t z = fhEGenZGen->GetYaxis()->GetBinCenter(tz); | |
336 | Int_t bine = jetSpec->GetGenSpectrum()->GetXaxis()->FindBin(ej); | |
337 | Int_t binz = jetSpec->GetGenSpectrum()->GetYaxis()->FindBin(z); | |
338 | Float_t cont = jetSpec->GetGenSpectrum()->GetBinContent(bine,binz); | |
339 | Float_t err = jetSpec->GetGenSpectrum()->GetBinError(bine,binz); | |
bebe25a4 | 340 | // merging of bins happens here! |
6d75bdb8 | 341 | jetSpec->GetGenSpectrum()->SetBinContent(bine, binz, cont + fhEGenZGen->GetBinContent(te, tz)); |
342 | jetSpec->GetGenSpectrum()->SetBinError(bine, binz, err + fhEGenZGen->GetBinError(te, tz)); | |
343 | } | |
6d75bdb8 | 344 | |
6d75bdb8 | 345 | |
df65bddb | 346 | Printf("Bins %.0E",jetSpec->GetCorrelation()->GetNbins()); |
6d75bdb8 | 347 | for (Int_t idx=1; idx<=fhCorrelation->GetNbins(); idx++) |
348 | { | |
349 | //printf("idx: %d\n",idx); | |
350 | Double_t var[4]; | |
351 | Int_t bin[4]; | |
352 | Float_t BinContent = fhCorrelation->GetBinContent(idx, bin); | |
353 | var[0] = fhCorrelation->GetAxis(0)->GetBinCenter(bin[0]); | |
354 | var[1] = fhCorrelation->GetAxis(1)->GetBinCenter(bin[1]); | |
355 | var[2] = fhCorrelation->GetAxis(2)->GetBinCenter(bin[2]); | |
356 | var[3] = fhCorrelation->GetAxis(3)->GetBinCenter(bin[3]); | |
357 | bin[0] = jetSpec->GetCorrelation()->GetAxis(0)->FindBin(var[0]); | |
358 | bin[1] = jetSpec->GetCorrelation()->GetAxis(1)->FindBin(var[1]); | |
359 | bin[2] = jetSpec->GetCorrelation()->GetAxis(2)->FindBin(var[2]); | |
360 | bin[3] = jetSpec->GetCorrelation()->GetAxis(3)->FindBin(var[3]); | |
361 | Float_t cont = jetSpec->GetCorrelation()->GetBinContent(bin); | |
362 | Float_t err = jetSpec->GetCorrelation()->GetBinError(bin); | |
bebe25a4 | 363 | // merging of bins happens here! |
6d75bdb8 | 364 | jetSpec->GetCorrelation()->SetBinContent(bin, cont + fhCorrelation->GetBinContent(idx)); |
365 | jetSpec->GetCorrelation()->SetBinError(bin, err + fhCorrelation->GetBinError(idx)); | |
366 | } | |
df65bddb | 367 | // need number of entries for correct drawing |
368 | jetSpec->GetCorrelation()->SetEntries(fhCorrelation->GetEntries()); | |
369 | ||
370 | ||
371 | file = TFile::Open("gen_pwg4spec.root", "RECREATE"); | |
372 | jetSpec->SaveHistograms(); | |
373 | Printf("Bins %.0E",jetSpec->GetCorrelation()->GetNbins()); | |
374 | file->Close(); | |
375 | jetSpec->GetGenSpectrum()->Reset(); | |
376 | printf("true distribution has been set\n"); | |
377 | ||
378 | // reconstructed jets (measured distribution) | |
379 | ||
380 | ||
381 | // Response function | |
382 | jetSpec->GetCorrelation()->Reset(); | |
383 | jetSpec->GetCorrelation()->Sumw2(); | |
384 | jetSpec->GetGenSpectrum()->Reset(); | |
385 | jetSpec->GetRecSpectrum()->Reset(); | |
386 | ||
387 | ||
388 | file = new TFile(fileNameReal,"read"); | |
389 | tlist = dynamic_cast<TList*> (file->Get("pwg4spec")); | |
390 | ||
391 | ||
392 | TH2F *fhERecZRec = (TH2F*)(tlist->FindObject("fh2ERecZRec")); | |
393 | ||
394 | for (Int_t me=1; me<=fhERecZRec->GetNbinsX(); me++) | |
395 | for (Int_t mz=1; mz<=fhERecZRec->GetNbinsY(); mz++) | |
396 | { | |
397 | Float_t erec = fhERecZRec->GetXaxis()->GetBinCenter(me); | |
398 | Float_t zm = fhERecZRec->GetYaxis()->GetBinCenter(mz); | |
399 | Int_t bine = jetSpec->GetRecSpectrum()->GetXaxis()->FindBin(erec); | |
400 | Int_t binz = jetSpec->GetRecSpectrum()->GetYaxis()->FindBin(zm); | |
401 | Float_t cont = jetSpec->GetRecSpectrum()->GetBinContent(bine, binz); | |
402 | Float_t err = jetSpec->GetRecSpectrum()->GetBinError(bine, binz); | |
403 | jetSpec->GetRecSpectrum()->SetBinContent(bine, binz, cont + fhERecZRec->GetBinContent(me, mz)); | |
404 | jetSpec->GetRecSpectrum()->SetBinError(bine, binz, err + fhERecZRec->GetBinError(me, mz)); | |
405 | } | |
6d75bdb8 | 406 | |
4dbfdecc | 407 | |
408 | // for control again, but now from the rec file | |
409 | // generated jets (true distribution) | |
410 | TH2F *fhEGenZGen = (TH2F*)(tlist->FindObject("fh2EGenZGen")); | |
411 | for (Int_t te=1; te<=fhEGenZGen->GetNbinsX(); te++) | |
412 | for (Int_t tz=1; tz<=fhEGenZGen->GetNbinsY(); tz++) | |
413 | { | |
414 | Float_t ej = fhEGenZGen->GetXaxis()->GetBinCenter(te); | |
415 | Float_t z = fhEGenZGen->GetYaxis()->GetBinCenter(tz); | |
416 | Int_t bine = jetSpec->GetGenSpectrum()->GetXaxis()->FindBin(ej); | |
417 | Int_t binz = jetSpec->GetGenSpectrum()->GetYaxis()->FindBin(z); | |
418 | Float_t cont = jetSpec->GetGenSpectrum()->GetBinContent(bine,binz); | |
419 | Float_t err = jetSpec->GetGenSpectrum()->GetBinError(bine,binz); | |
420 | jetSpec->GetGenSpectrum()->SetBinContent(bine, binz, cont + fhEGenZGen->GetBinContent(te, tz)); | |
421 | jetSpec->GetGenSpectrum()->SetBinError(bine, binz, err + fhEGenZGen->GetBinError(te, tz)); | |
422 | } | |
423 | ||
424 | ||
6d75bdb8 | 425 | file = TFile::Open("rec_pwg4spec.root", "RECREATE"); |
426 | jetSpec->SaveHistograms(); | |
427 | file->Close(); | |
428 | ||
429 | printf("reconstructed distribution has been set\n"); | |
430 | printf("response map has been set\n"); | |
431 | ||
432 | } | |
433 | ||
6d75bdb8 | 434 | void correct(){ |
435 | // simple steering to correct a given distribution; | |
436 | loadlibs(); | |
4dbfdecc | 437 | // rec and gen |
438 | // FillSpecFromFiles("pwg4spec_15-50_all.root","pwg4spec_allpt.root"); | |
bebe25a4 | 439 | FillSpecFromFiles("pwg4spec_allpt_16.root","pwg4spec_allpt_16.root"); |
6d75bdb8 | 440 | |
441 | char name[100]; | |
bebe25a4 | 442 | sprintf(name, "unfolded_pwg4spec_16.root"); |
6d75bdb8 | 443 | |
444 | unfold("gen_pwg4spec.root", "unfolding", "rec_pwg4spec.root", name); | |
445 | //draw(name, "unfolding", 1); | |
446 | ||
447 | } | |
4dbfdecc | 448 | |
449 | void mergeJetAnaOutput(){ | |
450 | // This is used to merge the analysis-output from different | |
451 | // data samples/pt_hard bins | |
1e8729f0 | 452 | // in case the eventweigth was set to xsection/ntrials already, this |
453 | // is not needed. Both methods only work in case we do not mix different | |
4dbfdecc | 454 | // pt_hard bins, and do not have overlapping bins |
455 | ||
456 | const Int_t nBins = 2; | |
457 | // LHC08q jetjet100: Mean = 1.42483e-03, RMS = 6.642e-05 | |
458 | // LHC08r jetjet50: Mean = 2.44068e-02, RMS = 1.144e-03 | |
459 | // LHC08v jetjet15-50: Mean = 2.168291 , RMS = 7.119e-02 | |
1e8729f0 | 460 | // const Float_t xsection[nBins] = {2.168291,2.44068e-02}; |
461 | Float_t xsection[nBins] = {0,0}; | |
4dbfdecc | 462 | Float_t nTrials[nBins] = {0,0}; |
463 | Float_t sf[nBins] = {0,0}; | |
464 | ||
1e8729f0 | 465 | const char *cFile[nBins] = {"pwg4spec_0000000-1000000_LHC08v_jetjet15-50.root", |
466 | "pwg4spec_0000000-1000000_LHC08r_jetjet50.root"}; | |
4dbfdecc | 467 | |
468 | ||
469 | TList *lIn[nBins]; | |
470 | TFile *fIn[nBins]; | |
471 | for(int ib = 0;ib < nBins;++ib){ | |
472 | fIn[ib] = TFile::Open(cFile[ib]); | |
473 | lIn[ib] = (TList*)fIn[ib]->Get("pwg4spec"); | |
474 | TH1* hTrials = (TH1F*)lIn[ib]->FindObject("fh1PtHard_Trials"); | |
1e8729f0 | 475 | TProfile* hXsec = (TProfile*)lIn[ib]->FindObject("fh1Xsec"); |
476 | xsection[ib] = hXsec->GetBinContent(1); | |
4dbfdecc | 477 | nTrials[ib] = hTrials->Integral(); |
478 | sf[ib] = xsection[ib]/ nTrials[ib]; | |
479 | } | |
480 | ||
481 | TFile *fOut = new TFile("pwg4spec_allpt.root","RECREATE"); | |
482 | TList *lOut = new TList(); | |
483 | lOut->SetName(lIn[0]->GetName()); | |
484 | // for the start scale all... | |
485 | for(int ie = 0; ie < lIn[0]->GetEntries();++ie){ | |
486 | TH1 *h1Add = 0; | |
487 | THnSparse *hnAdd = 0; | |
488 | Printf("%d: %s",ie, lIn[0]->At(ie)->GetName()); | |
489 | for(int ib = 0;ib < nBins;++ib){ | |
490 | // dynamic cast does not work with cint | |
491 | TObject *h = lIn[ib]->At(ie); | |
492 | if(h->InheritsFrom("TH1")){ | |
493 | Printf("ib %d",ib); | |
494 | TH1 *h1 = (TH1*)h; | |
495 | if(ib==0){ | |
496 | h1Add = (TH1*)h1->Clone(h1->GetName()); | |
497 | h1Add->Scale(sf[ib]); | |
498 | } | |
499 | else{ | |
500 | h1Add->Add(h1,sf[ib]); | |
501 | } | |
502 | } | |
503 | else if(h->InheritsFrom("THnSparse")){ | |
504 | Printf("ib %d",ib); | |
505 | THnSparse *hn = (THnSparse*)h; | |
506 | if(ib==0){ | |
507 | hnAdd = (THnSparse*)hn->Clone(hn->GetName()); | |
508 | hnAdd->Scale(sf[ib]); | |
509 | } | |
510 | else{ | |
511 | hnAdd->Add(hn,sf[ib]); | |
512 | } | |
513 | } | |
514 | ||
515 | ||
516 | }// ib | |
517 | if(h1Add)lOut->Add(h1Add); | |
518 | else if(hnAdd)lOut->Add(hnAdd); | |
519 | } | |
520 | fOut->cd(); | |
521 | lOut->Write(lOut->GetName(),TObject::kSingleKey); | |
522 | fOut->Close(); | |
523 | ||
524 | ||
525 | ||
526 | } |