]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/macros/jetspectrum/runUnfoldingCF.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
fix bug from last commit - Marta V.
[u/mrichter/AliRoot.git] / PWGJE / macros / jetspectrum / runUnfoldingCF.C
CommitLineData
e2afec0c 1void Load();
2TList *GetResults(Char_t *testfile,char* listName);
3TObject* GetMeasuredSpectrum(TList *fContainer);
4TObject* GetGeneratedSpectrum(TList *fContainer);
5TObject* GetEfficiency(TList *fContainer,Int_t iCase = 44);
6THnSparse* GetResponseMatrix(TList *fContainer);
4f01f301 7THnSparse* CreateEmptyResponseMatrix();
e2afec0c 8Float_t GetNTrials(TList *fContainer);
9THnSparse* CreateGuessed(const THnSparse* h);
10void PrintErrors(THnSparse *h);
11
12Int_t fDebug = 10;
13
14int iRebin = 8; //
15
16// *************************************************************************************************************
17// Macro for unfolding jet spectra:
18// * First Unfolding is always the MC itself
19// * if bMCtest we take a second MC which is unfolded with the first MC otherwise real data is unfolded with first MC
20// * bMC2 controls wether we correct to charged particle jets or not
21// ############### CHANGELOG ################
22// 23.11. added method to create an empty response matrix
23
24void runUnfoldingCF( int gIterations = 10, bool bMC2 = true,bool bNoPrior = true,bool bMCtest = true, int iJetF = 0) {
25 Load();
26
27 const bool bPreSmooth = false;
28
29 // can be different depeding on the prior
30 Int_t fIterations = gIterations;
31 Int_t fIterations2 = gIterations;
32 Int_t fIterations3 = gIterations;
33
34 // the first one is usually MC only
35 TString cJetF;
36 TString listName;
37 TString testfile;
38 if(bMCtest) testfile = "allpt_lhc10e14_100903_Chunk1.root"; // the higher stats chunck
39 else testfile = "allpt_lhc10e14_100903.root";
40 TList *results = 0;
41 if(iJetF == 0){
42 cJetF = "anti-k_{T}";
43 listName = "spec2_jetsAOD_FASTJET04_jetsAODMC%s_FASTJET04_0000000000";
44 }
45 else if(iJetF == 1){
46 cJetF = "k_{T}";
47 listName = "spec2_jetsAOD_FASTKT04_jetsAODMC%s_FASTKT04_0000000000";
48 }
49 else if(iJetF == 2){
50 cJetF = "UA1";
51 listName = "spec2_jets_jetsAODMC%s_UA104_0000000000";
52 }
53 TList *results = GetResults(testfile.Data(),Form(listName.Data(),(bMC2?"2":"")));
54
55 // here comes the real data
56 TString cJetF2;
57 TString listName2;
58 TString testfile2;
59
60 if(bMCtest)testfile2 = "allpt_lhc10e14_100903_Chunk0.root"; // the smalle stats
61 else testfile2 = "/Users/kleinb/alice/jets/train/100910/PWG4_JetTasksOutput_Merge_bcd.root";
62 if(iJetF == 0){
63 cJetF = "anti-k_{T}";
64 if(bMCtest)listName2 = "spec2_jetsAOD_FASTJET04_jetsAODMC%s_FASTJET04_0000000000";
65 else listName2 = "spec2_jetsAOD_FASTJET04__0000000000";
66
67 }
68 else if(iJetF == 1){
69 cJetF = "k_{T}";
70 if(bMCtest)listName2 = "spec2_jetsAOD_FASTKT04_jetsAODMC%s_FASTKT04_0000000000";
71 else listName2 = "spec2_jetsAOD_FASTKT04__0000000000";
72
73 }
74 else if(iJetF == 2){
75 cJetF = "k_{T}";
76 if(bMCtest)listName2 = "spec2_jets_jetsAODMC%s_UA104_0000000000";
77 else listName2 = "spec2_jets__0000000000";
78
79 }
80 TList *results2 = GetResults(testfile2.Data(),Form(listName2.Data(),(bMC2?"2":"")));bool bScale = false;
81
82
83
84 if(!results){
85 Error("No output objects: Calculation will terminate here");
86 return;
87 }
88
89 bool b1D = true;
90 const Float_t kMaxPt = 140; // only for drawin...
91 const Float_t kMaxPt2 = 140;
92
93
94 char cString[] = Form("unfolded/101010_unfolding%s_MC%s%s%%s_%sIter%04d_Rebin%02d_jet%d",(bNoPrior?"_NoPrior":""),(bMC2?"2":""),(bPreSmooth?"_PreSmooth":""),(bMCtest?"MCtest_":""),gIterations,iRebin,iJetF);
95 char cPrintMask[] = Form("%s.png",cString);
96
97 TF3 *fSmooth = new TF3("fSmooth","[0]*pow(x,[1])+[2]*y+[3]*z",5.,150,-10,10,-10,10);
98 fSmooth = 0;
99 TF3 *fSmooth2 = 0;
100 TF3 *fSmooth3 = 0;
101
102 // fSmooth = 0;
103 AliLog::SetGlobalLogLevel(AliLog::kDebug);
104 // get the essential
105 if(fSmooth){
106 fSmooth->SetParameters(1,-8,0,0);
107 fSmooth2 = (TF3*)fSmooth->Clone("fSmooth2");
108 fSmooth3 = (TF3*)fSmooth->Clone("fSmooth3");
109 }
110
111
112
113
114 THnSparseF *efficiency = (THnSparseF*) GetEfficiency(results,41);
115 THnSparseF *efficiencyRecMatch = (THnSparseF*) GetEfficiency(results,131);
116 THnSparseF *response = (THnSparseF*) GetResponseMatrix(results);
117
118
119
120 THnSparseF *measuredIn2 = (THnSparseF*) GetMeasuredSpectrum(results2);
121 measuredIn2->Multiply(efficiencyRecMatch);
122 THnSparseF *measuredIn = (THnSparseF*) GetMeasuredSpectrum(results);
123 measuredIn->Multiply(efficiencyRecMatch);
124 THnSparseF *generatedOut = (THnSparseF*) GetGeneratedSpectrum(results);
125 THnSparseF *generatedOut2 = (THnSparseF*) GetGeneratedSpectrum(results2);
126
127 Float_t fTrials = GetNTrials(results);
128 Float_t fTrials2 = GetNTrials(results2);
129 if(bScale){
130 Printf("Trials %.3E %.3E",fTrials,fTrials2);
131 measuredIn2->Scale(1./fTrials2);
132 measuredIn->Scale(1./fTrials);
133 }
134
135 // set the content to zero above threshold
136 /*
137 const float lastPt = 250;
138 Int_t ibxyz[3];
139 for(int ibx = 1;ibx<=measuredIn2->GetAxis(0)->GetNbins();ibx++){
140 float pt = measuredIn2->GetAxis(0)->GetBinCenter(ibx);
141
142 if(pt<lastPt)continue;
143 ibxyz[0] = ibx;
144 for(int iby = 1;iby<=measuredIn2->GetAxis(1)->GetNbins();iby++){
145 ibxyz[1] = iby;
146 for(int ibz = 1;ibz<=measuredIn2->GetAxis(2)->GetNbins();ibz++){
147 ibxyz[2] = ibz;
148 measuredIn2->SetBinContent(ibxyz,0);
149 measuredIn2->SetBinError(ibxyz,0);
150 }
151 }
152 }
153 */
154
155
156 // for testing
157 const int nDim = 3;
158 const int dimrec[nDim] = {0,1,2};
159 const int dimgen[nDim] = {3,4,5};
160
161
162 /*
163 THnSparseF *generated = response->Projection(nDim,dimgen);
164 THnSparseF *measured = response->Projection(nDim,dimrec);
165 */
166
167 // create a guessed "a priori" distribution using binning of MC
168 THnSparse* guessed2 = CreateGuessed(generatedOut) ; // can at best take the measured?
169 THnSparse* guessed = CreateGuessed(generatedOut) ; // can at best take the measured?
170 Printf("%s:%d %d",(char*)__FILE__,__LINE__,guessed2->GetNbins());
171 Printf("%s:%d %d",(char*)__FILE__,__LINE__,guessed->GetNbins());
172
173 //---- Dbug show the errrors
174 // PrintErrors(guessed);
175 // PrintErrors(response);
176 // PrintErrors(efficiency);
177 // PrintErrors(measuredIn);
178 // PrintErrors(generatedOut);
179
180 Bool_t bCorrelatedErrors = true;
181 AliCFUnfolding unfolding("unfolding","title",3,response,efficiency,measuredIn,(bNoPrior?0:guessed));
182 unfolding.SetMaxNumberOfIterations(fIterations); // regulate flutuations...
183 unfolding.SetMaxChi2PerDOF(0);
184 if(fSmooth)unfolding.UseSmoothing(fSmooth);
185 if(bCorrelatedErrors){
186 unfolding.SetUseCorrelatedErrors(kTRUE);
187 unfolding.SetMaxConvergencePerDOF(0.01*0.01);
188 }
189
190 unfolding.Unfold();
191
192 THnSparse* result = unfolding.GetUnfolded();
193 THnSparse* estMeasured = unfolding.GetEstMeasured();
194
195
196 //----
197 AliCFUnfolding unfolding2("unfolding2","title",3,response,efficiency,measuredIn2,(bNoPrior?0:guessed));
198 // AliCFUnfolding unfolding2("unfolding2","title",3,response,efficiency,measuredIn2,0);
199 unfolding2.SetMaxNumberOfIterations(fIterations2);
200 unfolding2.SetMaxChi2PerDOF(0);
201 // carefull with 0x0 pointer neighbouring bin smoothing is switched on...
202 if(fSmooth2)unfolding2.UseSmoothing(fSmooth2);
203 if(bCorrelatedErrors){
204 unfolding2.SetUseCorrelatedErrors(kTRUE);
205 unfolding2.SetMaxConvergencePerDOF(0.01*0.01);
206 }
207 unfolding2.Unfold();
208
209
210
211 THnSparse* estMeasuredTmp = unfolding2.GetEstMeasured();
212 THnSparse* result2 = 0;
213 THnSparse* estMeasured2 = 0;
214 if(bPreSmooth){
215 // use the result of the previous as smoothed input...
216 AliCFUnfolding unfolding3("unfolding3","title",3,response,efficiency,estMeasuredTmp,(bNoPrior?0:guessed));
217 unfolding3.SetMaxNumberOfIterations(fIterations3);
218 unfolding3.SetMaxChi2PerDOF(0);
219 unfolding3.UseSmoothing(fSmooth3);
220 unfolding3.Unfold();
221 result2 = unfolding3.GetUnfolded();
222 estMeasured2 = unfolding3.GetEstMeasured();
223 }
224 else{
225 result2 = unfolding2.GetUnfolded();
226 estMeasured2 = estMeasuredTmp;
227 }
228
229
230
231 TCanvas * canvas = new TCanvas("canvas","",1200,900);
232 canvas->Divide(2,3);
233 TCanvas * cPrint = new TCanvas("cPrint","");
234
235
236
237 TFile *f1 = new TFile(Form("%s.root",Form(cString,"")),"RECREATE");
238 TParameter<Long_t>* fIterationsPara = new TParameter<Long_t> ("fIterations", 0);
239 fIterationsPara->SetVal((Long_t)gIterations);
240 fIterationsPara->Write();
241
242 TDirectory *dMC = f1->mkdir("unfoldMC");
243 TDirectory *dReal = f1->mkdir("unfoldReal");
244
245 gROOT->cd();
246 // color code black is unfolded
247 // blue is measured
248 // red is generated
249 // kGreen is guessed
250
251
252
253 if(b1D){
254
255 TH1* h_gen = generatedOut->Projection(0);
256
257 h_gen->SetMarkerColor(kRed);
258 h_gen->SetMarkerStyle(kFullSquare);
259 h_gen->SetName("generated");
260 h_gen->SetTitle("generated");
261 TH1* h_meas = measuredIn->Projection(0);
262 h_meas->SetMarkerColor(kBlue);
263 h_meas->SetMarkerStyle(kFullSquare);
264 h_meas->SetName("measured");
265 h_meas->SetTitle("measured");
266 TH1* h_guess = guessed->Projection(0);
267 h_guess->SetMarkerColor(kGreen);
268 h_guess->SetMarkerStyle(kFullSquare);
269 h_guess->SetName("guessed");
270 h_guess->SetTitle("guesse");
271 TH1* h_unf = result->Projection(0);
272 h_unf->SetMarkerColor(kBlack);
273 h_unf->SetMarkerStyle(kFullSquare);
274 h_unf->SetName("unfolded");
275 h_unf->SetTitle("unfolded");
276 TH1* h_estmeas = estMeasured->Projection(0);
277 h_estmeas->SetMarkerColor(kGray);
278 h_estmeas->SetMarkerStyle(kFullSquare);
279 h_estmeas->SetName("estmeas");
280 h_estmeas->SetTitle("estmeas");
281
282 TLegend *leg = new TLegend(0.6,0.6,0.85,0.8);
283 leg->SetFillColor(0);
284 leg->SetTextFont(gStyle->GetTextFont());
285 leg->SetBorderSize(0);
286 leg->AddEntry(h_meas,"measured","P");
287 leg->AddEntry(h_gen,"generated","P");
288 leg->AddEntry(h_unf,"unfolded","P");
289 leg->AddEntry(h_estmeas,"R * u","P");
290
291
292
293 canvas->cd(1)->SetLogy();
294 h_gen->SetAxisRange(0,kMaxPt);
295 h_meas->SetAxisRange(0,kMaxPt);
296 h_meas->SetXTitle("p_{T} (GeV)");
297 h_meas->SetYTitle("yield (arb. units)");
298 h_meas->DrawCopy("P");
299 h_gen->DrawCopy("Psame");
300 h_unf->DrawCopy("Psame");
301 h_estmeas->DrawCopy("Psame");
302 leg->Draw("");
303 cPrint->cd()->SetLogy();
304 h_meas->DrawCopy("P");
305 h_gen->DrawCopy("Psame");
306 h_unf->DrawCopy("Psame");
307 h_estmeas->DrawCopy("Psame");
308 leg->Draw("");
309 canvas->Update();
310 cPrint->Update();
311 cPrint->SaveAs(Form(cPrintMask,"spectrum_mc"));
312 if(!gROOT->IsBatch()){
313 if(getchar()=='q')return;
314 }
315 // the same for the real data
316 TH1* h_gen2 = generatedOut2->Projection(0);
317
318 h_gen2->SetMarkerColor(kRed);
319 h_gen2->SetMarkerStyle(kFullCircle);
320 h_gen2->SetName("generated2");
321 h_gen2->SetTitle("generated2");
322 TH1* h_meas2 = measuredIn2->Projection(0);
323 h_meas2->SetMarkerColor(kBlue);
324 h_meas2->SetMarkerStyle(kFullCircle);
325 h_meas2->SetName("measured2");
326 h_meas2->SetTitle("measured2");
327 TH1* h_guess2 = guessed2->Projection(0);
328 h_guess2->SetMarkerColor(kGreen);
329 h_guess2->SetMarkerStyle(kFullCircle);
330 h_guess2->SetName("guessed2");
331 h_guess2->SetTitle("guesse2");
332 TH1* h_unf2 = result2->Projection(0);
333 PrintErrors(result2);
334 h_unf2->SetMarkerColor(kBlack);
335 h_unf2->SetMarkerStyle(kFullCircle);
336 h_unf2->SetName("unfolded2");
337 h_unf2->SetTitle("unfolded2");
338 TH1* h_estmeas2 = estMeasured2->Projection(0);
339 h_estmeas2->SetMarkerColor(kGray);
340 h_estmeas2->SetMarkerStyle(kFullCircle);
341 h_estmeas2->SetName("estmeas2");
342 h_estmeas2->SetTitle("estmeas2");
343
344 TLegend *leg2 = new TLegend(0.6,0.6,0.85,0.8);
345 leg2->SetFillColor(0);
346 leg2->SetTextFont(gStyle->GetTextFont());
347 leg2->SetBorderSize(0);
348 leg2->AddEntry(h_meas2,"measured","P");
349 leg2->AddEntry(h_gen2,"generated","P");
350 leg2->AddEntry(h_unf2,"unfolded","P");
351 leg2->AddEntry(h_estmeas2,"R * u","P");
352
353 Printf("Integral Measured: %E Unfolded %E R*U %E",h_meas2->Integral(),h_unf2->Integral(),h_estmeas2->Integral());
354
355
356
357 canvas->cd(2)->SetLogy();
358 //
359
360 h_meas2->SetXTitle("p_{T} (GeV)");
361 h_meas2->SetYTitle("yield (arb. units)");
362 // h_gen2->DrawCopy("P");
363 h_unf2->SetAxisRange(0,kMaxPt2);
364 h_meas2->SetAxisRange(0,kMaxPt2);
365 h_meas2->DrawCopy("P");
366 h_unf2->DrawCopy("Psame");
367 h_estmeas2->DrawCopy("Psame");
368 leg2->Draw();
369 cPrint->cd()->SetLogy();
370 h_unf2->SetAxisRange(0,kMaxPt2);
371 h_meas2->SetAxisRange(0,kMaxPt2);
372 h_meas2->DrawCopy("P");
373 h_unf2->DrawCopy("Psame");
374 h_estmeas2->DrawCopy("Psame");
375 canvas->Update();
376 cPrint->Update();
377 leg2->Draw();
378 cPrint->SaveAs(Form(cPrintMask,"spectrum_real"));
379 if(!gROOT->IsBatch()){
380 if(getchar()=='q')return;
381 }
382
383
384
385 // residuals
386
387 TH1D *hResiduals = (TH1D*) h_meas->Clone("hResiduals");
388 hResiduals->Reset();
389 for(int ib = 1;ib < hResiduals->GetNbinsX();ib++){
390 float val1 = h_meas->GetBinContent(ib);
391 float err1 = h_meas->GetBinError(ib);
392 float val2 = h_estmeas->GetBinContent(ib);
393 if(err1>0){
394 Float_t res = (val1-val2)/err1;
395 Float_t res_err = (val1-val2)/err1*0.01; // error bars of 1%
396 hResiduals->SetBinContent(ib,res);
397 hResiduals->SetBinError(ib,0.01);
398 }
399
400 }
401 hResiduals->SetXTitle("p_{T} (GeV)");
402 hResiduals->SetYTitle("#frac{m - R * u}{e_{m}}");
403 hResiduals->SetMaximum(4.5);
404 hResiduals->SetMinimum(-4.5);
405
406
407 canvas->cd(3);
408 // hRatioGenUnf->DrawCopy("p");
409 hResiduals->SetAxisRange(0,kMaxPt);
410 hResiduals->DrawCopy("P");
411 cPrint->cd();
412 gPad->SetLogy(0);
413 hResiduals->DrawCopy("P");
414 canvas->Update();
415 cPrint->Update();
416 cPrint->SaveAs(Form(cPrintMask,"residuals_mc"));
417 if(!gROOT->IsBatch()){
418 if(getchar()=='q')return;
419 }
420
421 TH1D *hResiduals2 = (TH1D*)h_meas2->Clone("hResiduals2");
422 hResiduals2->Reset();
423 for(int ib = 1;ib < hResiduals2->GetNbinsX();ib++){
424 float val1 = h_meas2->GetBinContent(ib);
425 float err1 = h_meas2->GetBinError(ib);
426 float val2 = h_estmeas2->GetBinContent(ib);
427 if(err1>0){
428 Float_t res = (val1-val2)/err1;
429 Float_t res_err = (val1-val2)/err1*0.01; // error bars of 1%
430 hResiduals2->SetBinContent(ib,res);
431 hResiduals2->SetBinError(ib,0.01);
432 }
433
434 }
435 hResiduals2->SetMaximum(4.5);
436 hResiduals2->SetMinimum(-4.5);
437 hResiduals2->SetYTitle("#frac{m - R * u}{e_{m}}");
438 hResiduals2->SetAxisRange(0,kMaxPt2);
439 cPrint->cd();
440 gPad->SetLogy(0);
441 hResiduals2->DrawCopy("P");
442 canvas->Update();
443 cPrint->Update();
444 cPrint->SaveAs(Form(cPrintMask,"residuals_real"));
445 if(!gROOT->IsBatch()){
446 if(getchar()=='q')return;
447 }
448
449 canvas->Update();
450 cPrint->Update();
451
452 if(!gROOT->IsBatch()){
453 if(getchar()=='q')return;
454 }
455 hResiduals2->DrawCopy("P");
456 canvas->Update();
457 cPrint->Update();
458 if(!gROOT->IsBatch()){
459 if(getchar()=='q')return;
460 }
461
462 TH1D *hRatioGenUnf2 = (TH1D*) h_gen2->Clone("hRatioGenUnf2");
463 hRatioGenUnf2->Divide(h_unf2);
464 canvas->cd(4);
465 // hRatioGenUnf2->DrawCopy("p");
466
467 hResiduals2->DrawCopy("P");
468
469 TH1D *hRatioUnfMeas = (TH1D*) h_unf->Clone("hRatioUnfMeas");
470 hRatioUnfMeas->SetXTitle("p_{T} (GeV)");
471 hRatioUnfMeas->SetYTitle("unfolded/meas");
472 hRatioUnfMeas->Divide(h_meas);
473 hRatioUnfMeas->SetMaximum(15);
474 if(bMC2)hRatioUnfMeas->SetMaximum(5);
475 hRatioUnfMeas->SetMinimum(0);
476 TH1D *hRatioGenMeas = (TH1D*)h_gen->Clone("hRatioGenMeas");
477 hRatioGenMeas->Divide(h_meas);
478
479 canvas->cd(5);
480
481 hRatioUnfMeas->SetAxisRange(0.,kMaxPt);
482 hRatioUnfMeas->DrawCopy("p");
483 hRatioGenMeas->DrawCopy("psame");
484
485
486
487
488 TH1D *hRatioUnfMeas2 = (TH1D*) h_unf2->Clone("hRatioUnfMeas2");
489 hRatioUnfMeas2->Divide(h_meas2);
490 canvas->cd(6);
491 hRatioUnfMeas2->SetXTitle("p_{T} (GeV)");
492 hRatioUnfMeas2->SetYTitle("unfolded/meas");
493 hRatioUnfMeas2->SetAxisRange(0.,kMaxPt2);
494 hRatioUnfMeas2->SetMaximum(15);
495 if(bMC2)hRatioUnfMeas2->SetMaximum(5);
496 hRatioUnfMeas2->SetMinimum(0);
497 hRatioUnfMeas2->DrawCopy("p");
498 hRatioUnfMeas->DrawCopy("psame");
499
500 // plot the effieciencies
501 cPrint->cd();
502 TH1* hEffGen = efficiency->Projection(0);
503 hEffGen->SetName("hEffGen");
504 TH1* hEffRec = efficiencyRecMatch->Projection(0);
505 hEffRec->SetName("hEffRec");
506
507 hEffGen->SetAxisRange(0,kMaxPt2);
508 hEffRec->SetAxisRange(0,kMaxPt2);
509
510 hEffRec->SetXTitle("p_{T,gen/rec} (GeV/c)");
511 hEffRec->SetYTitle("Efficiency");
512 hEffGen->SetXTitle("p_{T,gen/rec} (GeV/c)");
513 hEffGen->SetYTitle("Efficiency");
514
515 hEffGen->SetMarkerColor(kBlue);
516 hEffGen->SetMarkerStyle(kFullCircle);
517 hEffRec->SetMarkerColor(kBlue);
518 hEffRec->SetMarkerStyle(kOpenCircle);
519 hEffGen->SetMaximum(1.2);
520 hEffGen->SetMinimum(0.3);
521 hEffGen->DrawCopy("p");
522 hEffRec->DrawCopy("psame");
523 cPrint->Update();
524 cPrint->SaveAs(Form(cPrintMask,"effs"));
525
526 h_unf->SetDirectory(dMC);
527 h_meas->SetDirectory(dMC);
528 hResiduals->SetDirectory(dMC);
529 h_estmeas->SetDirectory(dMC);
530 h_gen->SetDirectory(dMC);
531
532
533 // put the used effs and response to the unfolded of the MC
534
535 hEffGen->SetDirectory(dMC);
536 hEffRec->SetDirectory(dMC);
537 // response->SetDirectory(dMC);
538 dMC->cd();
539 guessed->Write();
540 response->Write();
541 measuredIn->Write();
542 gROOT->cd();
543
544
545 h_unf2->SetDirectory(dReal);
546 h_meas2->SetDirectory(dReal);
547 hResiduals2->SetDirectory(dReal);
548 h_estmeas2->SetDirectory(dReal);
549 h_gen2->SetDirectory(dReal); // will be empty if we use real data, but can check with different data set as well...
550
551 dReal->cd();
552 guessed2->Write();
553 measuredIn2->Write();
554 gROOT->cd();
555 // store the parameters of the unfloding with tparamter...
556
557
558 }
559 canvas->cd();
560 canvas->Modified();
561 canvas->SaveAs(Form(cPrintMask,"all"));
562
563 f1->Write();
564 f1->Close();
565
566
567 return;
568 if(fDebug)Printf("Line:%d %s Entries %f",__LINE__,h_gen->GetName(),h_gen->GetEntries());
569 //printf("c1\n");
570 if(b1D)h_gen->DrawCopy("P");
571 else h_gen->DrawCopy("lego2");
572
573
574
575
576 if(b1D) canvas->cd(2)->SetLogy();
577 else canvas->cd(2)->SetLogz();
578 TH1* h_meas = 0;
579 if(b1D) h_meas = measuredIn->Projection(0);
580 else h_guessed = measuredIn->Projection(0,1);
581 h_meas->SetName("measured");
582 h_meas->SetTitle("measured");
583 if(fDebug)Printf("Line:%d %s Entries %f",__LINE__,h_meas->GetName(),h_meas->GetEntries());
584
585 //printf("c2\n");
586 if(b1D) h_meas->Draw("hist");
587 else h_meas->Draw("lego2");
588
589
590 if(b1D) canvas->cd(3)->SetLogy();
591 else canvas->cd(3)->SetLogz();
592
593 TH1* h_guessed = 0;
594 if(b1D) h_guessed = guessed2->Projection(0);
595 else h_guessed = guessed2->Projection(0,1);
596
597 h_guessed->SetTitle("a priori");
598 h_guessed->SetName("apriori");
599 if(fDebug)Printf("Line:%d %s Entries %f",__LINE__,h_guessed->GetName(),h_guessed->GetEntries());
600 //printf("c3\n");
601 if(b1D) h_guessed->Draw("hist");
602 else h_guessed->Draw("lego2");
603
604 canvas->cd(4);
605 TH1* h_eff = 0;
606 if(b1D) h_eff = efficiency->Projection(0);
607 else h_eff = efficiency->Projection(0,1);
608
609 h_eff->SetTitle("efficiency");
610 h_eff->SetName("efficiency");
611 if(fDebug)Printf("Line:%d %s Entries %f",__LINE__,h_eff->GetName(),h_eff->GetEntries());
612 //printf("c4\n");
613 if(b1D) h_eff->Draw("hist");
614 else h_eff->Draw("lego2");
615
616
617 if(b1D) canvas->cd(5)->SetLogy();
618 else canvas->cd(5)->SetLogz();
619 TH1* h_unf = 0;
620 if(b1D) h_unf = result->Projection(0);
621 else h_unf = result->Projection(0,1);
622
623 h_unf->SetName("unfolded");
624 h_unf->SetTitle("unfolded");
625 //printf("c5\n");
626 if(b1D) h_unf->Draw("hist");
627 else h_unf->Draw("lego2");
628
629 canvas->cd(6);
630 TH1* ratio = 0;
631 if(b1D){
632 ratio = (TH1D*)h_unf->Clone("ratio");
633 ratio->SetTitle("unfolded / generated");
634 ratio->Divide(h_unf,h_gen,1,1);
635 }
636 else{
637 ratio = (TH2D*)h_unf->Clone("ratio");
638 ratio->SetTitle("unfolded / generated");
639 ratio->Divide(h_unf,h_gen,1,1);
640 }
641// ratio->Draw("cont4z");
642// ratio->Draw("surf2");
643 //printf("c6\n");
644 if(b1D)ratio->Draw("hist");
645 else ratio->Draw("colz");
646
647 return;
648
649 canvas->cd(7);
650 TH2* orig = unfolding.GetOriginalPrior()->Projection(1,0);
651 orig->SetName("originalprior");
652 orig->SetTitle("original prior");
653 //printf("c7\n");
654 orig->Draw("lego2");
655
656 canvas->cd(8);
657 THnSparseF* corrected = (THnSparseF*)measured->Clone("corrected");
658 //corrected->ApplyEffCorrection(*efficiency);
659 TH2D* corr = corrected->Projection(0,1);
660 corr->SetTitle("simple correction");
661 corr->SetName("simplecorrection");
662 //printf("c8\n");
663 corr->Draw("lego2");
664
665 canvas->cd(9);
666 TH2D* ratio2 = (TH2D*) corr->Clone();
667 ratio2->Divide(corr,h_gen,1,1);
668 ratio2->SetTitle("simple correction / generated");
669 //printf("c9\n");
670 ratio2->Draw("lego2");
671
672 return;
673}
674
675// ====================================================================
676
677
678void Load(){
679 gSystem->Load("libTree");
680 gSystem->Load("libPhysics");
681 gSystem->Load("libHist");
682 gSystem->Load("libVMC");
683 gSystem->Load("libSTEERBase");
684 gSystem->Load("libAOD");
685 gSystem->Load("libESD");
686 gSystem->Load("libANALYSIS");
687 gSystem->Load("libANALYSISalice");
688 gSystem->Load("libJETAN");
689 gSystem->Load("libCORRFW");
690 gSystem->Load("libPWG4JetTasks");
691 gStyle->SetPalette(1);
692 gStyle->SetOptStat(0);
693}
694
695TList *GetResults( Char_t *testfile, Char_t *cName){
696 //
697 // read output
698 //
699 TFile *f = TFile::Open(testfile);
700 if(!f || f->IsZombie()){
701 Error("File not readable");
702 return 0x0;
703 }
704 TDirectory *dir = dynamic_cast<TDirectory*>(f->Get(Form("PWG4_%s",cName)));
705 if(!dir){
706 Printf("%s:%d: Output directory %s not found",(char*)__FILE__,__LINE__,Form("PWG4_%s",cName));
707 f->ls();
708 f->Close(); delete f;
709 return 0x0;
710 }
711 TList *l = dynamic_cast<TList *>(dir->Get(Form("pwg4%s",cName)));
712 if(!l){
713 Printf("%s:%d: Output list %s not found",(char*)__FILE__,__LINE__,Form("pwg4%s",cName));
714 dir->ls();
715 f->Close(); delete f;
716 return 0x0;
717 }
718 // TList *returnlist = dynamic_cast<TList *>(l->Clone());
719 // f->Close(); delete f;
720 return l;
721}
722
723
724TObject* GetMeasuredSpectrum(TList *fContainer) {
725 THnSparseF* htmp = (THnSparseF*)fContainer->FindObject(Form("fhnJetContainer%d",AliAnalysisTaskJetSpectrum2::kMaxStep+AliAnalysisTaskJetSpectrum2::kStep1));
726 THnSparseF* hm = (THnSparseF*)htmp->Clone("measured");
727 hm->SetName("measured");
728 hm->SetTitle("measured");
729 if(iRebin==1)return hm;
730 Int_t fRebin[3];
731 fRebin[0] = iRebin;
732 fRebin[1] = 1;
733 fRebin[2] = 1;
734 THnSparse *hm2 = hm->Rebin(fRebin); // this creates a new histo...
735 // hm2->Scale(iRebin);
736 if(fDebug)Printf("Line: %d %s Entries %f",__LINE__,hm->GetName(),hm->GetEntries());
737 return hm2;
738}
739
740Float_t GetNTrials(TList *fContainer){
741 TH1F* htmp = (TH1F*)fContainer->FindObject("fh1Trials");
742 return htmp->Integral();
743}
744
745TObject* GetGeneratedSpectrum(TList *fContainer) {
746 // the generated jet spectrum within eta < 0.5
747 THnSparseF* htmp = (THnSparseF*)fContainer->FindObject(Form("fhnJetContainer%d",AliAnalysisTaskJetSpectrum2::kStep1));
748 THnSparseF* hg = (THnSparseF*)htmp->Clone("generated");
749 hg->SetName("generated");
750 hg->SetTitle("generated");
751
752 if(iRebin==1)return hg;
753 Int_t fRebin[3];
754 fRebin[0] = iRebin;
755 fRebin[1] = 1;
756 fRebin[2] = 1;
757 THnSparse* hg2 = hg->Rebin(fRebin);
758 // hg2->Scale(iRebin);
759 if(fDebug)Printf("Line: %d %s Entries %f",__LINE__,hg->GetName(),hg->GetEntries());
760 return hg2;
761}
762
763TObject* GetEfficiency(TList *fContainer,Int_t iCase) {
764
765 static int count = 0;
766 THnSparseF* hn1 = 0;
767 THnSparseF* hn2 = 0;
768 THnSparseF* he = 0;
769 if(iCase==44){
770 // Unitiy efficiency...
771 hn1 = (THnSparseF*)fContainer->FindObject(Form("fhnJetContainer%d",AliAnalysisTaskJetSpectrum2::kStep4));
772 hn2 = (THnSparseF*)fContainer->FindObject(Form("fhnJetContainer%d",AliAnalysisTaskJetSpectrum2::kStep4));
773 }
774 else if(iCase==14){
775 // eficiency Generated with gen < 0.5 --> gen with Rec partner in 0.5
776 hn1 = (THnSparseF*)fContainer->FindObject(Form("fhnJetContainer%d",AliAnalysisTaskJetSpectrum2::kStep1));
777 hn2 = (THnSparseF*)fContainer->FindObject(Form("fhnJetContainer%d",AliAnalysisTaskJetSpectrum2::kStep4));
778 }
779 else if(iCase==41){
780 // eficiency Generated with gen < 0.5 --> gen with Rec partner in 0.5
781 hn1 = (THnSparseF*)fContainer->FindObject(Form("fhnJetContainer%d",AliAnalysisTaskJetSpectrum2::kStep4));
782 hn2 = (THnSparseF*)fContainer->FindObject(Form("fhnJetContainer%d",AliAnalysisTaskJetSpectrum2::kStep1));
783 }
784 else if(iCase==131){
785 // eficiency reconstruted with < 0.5 --> with partner
786 hn1 = (THnSparseF*)fContainer->FindObject(Form("fhnJetContainer%d",AliAnalysisTaskJetSpectrum2::kMaxStep+AliAnalysisTaskJetSpectrum2::kStep3));
787 hn2 = (THnSparseF*)fContainer->FindObject(Form("fhnJetContainer%d",AliAnalysisTaskJetSpectrum2::kMaxStep+AliAnalysisTaskJetSpectrum2::kStep1));
788 }
789
790 if(iRebin==1){
791 he = (THnSparseF*)hn1->Clone("efficiency");
792 he->Divide(hn1,hn2,1.,1.,"B");
793 he->SetName(Form("efficiency_%d_%d",iCase,count));
794 he->SetTitle(Form("efficiency_%d_%d",iCase,count));
795 count++;
796 return he;
797 }
798 Int_t fRebin[3];
799 fRebin[0] = iRebin;
800 fRebin[1] = 1;
801 fRebin[2] = 1;
802
803 THnSparseF* hn1_2 = hn1->Rebin(fRebin);
804 THnSparseF* hn2_2 = hn2->Rebin(fRebin);
805 he = (THnSparseF*)hn1_2->Clone("efficiency");
806 he->Divide(hn1_2,hn2_2,1.,1.,"B");
807 he->SetName(Form("efficiency_%d_%d",iCase,count));
808 he->SetTitle(Form("efficiency_%d_%d",iCase,count));
809 return he;
810}
811
4f01f301 812THnSparse* CreateEmptyResponseMatrix(){
e2afec0c 813 // Creates an empty response matrix as in the Spectrum2 Task
814
815 const Int_t kNvar = 3 ; //number of variables on the grid:pt,eta, phi
816 const Double_t kPtmin = 0.0, kPtmax = 320.;
817 const Double_t kEtamin = -3.0, kEtamax = 3.0;
818 const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
819 const Double_t kZmin = 0., kZmax = 1;
820
821 //arrays for the number of bins in each dimension
822 Int_t iBin[kNvar];
823 iBin[0] = 320; //bins in pt
824 iBin[1] = 1; //bins in eta
825 iBin[2] = 1; // bins in phi
826
827 //arrays for lower bounds :
828 Double_t* binEdges[kNvar];
829 for(Int_t ivar = 0; ivar < kNvar; ivar++)
830 binEdges[ivar] = new Double_t[iBin[ivar] + 1];
831
832 //values for bin lower bounds
833 // for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)TMath::Power(10,TMath::Log10(kPtmin) + (TMath::Log10(kPtmax)-TMath::Log10(kPtmin))/iBin[0]*(Double_t)i);
834 for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
835 for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
836 for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
837
838 Int_t thnDim[2*kNvar];
839 for (int k=0; k<kNvar; k++) {
840 //first half : reconstructed
841 //second half : MC
842 thnDim[k] = iBin[k];
843 thnDim[k+kNvar] = iBin[k];
844 }
845
846 THnSparseF* fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
847 for (int k=0; k<kNvar; k++) {
848 fhnCorrelation->SetBinEdges(k,binEdges[k]);
849 fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
850 }
851 fhnCorrelation->Sumw2();
852 return fhnCorrelation;
853
854}
855
856THnSparse* GetResponseMatrix(TList *fContainer) {
857 THnSparse* h = (THnSparse*)fContainer->FindObject("fhnCorrelation");
858 if(fDebug)Printf("Line: %d %s Entries %f",__LINE__,h->GetName(),h->GetEntries());
859 if(iRebin==1)return h;
860 Int_t fRebin[6];
861 fRebin[0] = iRebin;
862 fRebin[1] = 1;
863 fRebin[2] = 1;
864 fRebin[3] = iRebin;
865 fRebin[4] = 1;
866 fRebin[5] = 1;
867
868 THnSparse* h2 = h->Rebin(fRebin);
869 return h2;
870}
871
872THnSparse* CreateGuessed(const THnSparse* h) {
873
874 // is already rebinned
875 THnSparse* guessed = (THnSparse*) h->Clone();
876 static int icount = 0;
877 guessed->SetName(Form("%s_guess_%d",h->GetName(),icount++));
878
879
880 return guessed;
881
882}
883
884void PrintErrors(THnSparse *h){
885
886 for(Long_t i = 0; i< h->GetNbins();i++){
887 Double_t val = h->GetBinContent(i);
888 Double_t err = h->GetBinError (i);
889 if(val>0){
890 Printf("%s %ld %1.3E +- %1.3E rel Error: %1.3E",h->GetName(),i,val,err,err/val);
891 }
892 }
893
894}