Coding rules (Markus)
[u/mrichter/AliRoot.git] / PWG3 / hfe / testUnfolding.C
CommitLineData
259c3296 1void Load();
2TList *GetResults(const Char_t *testfile);
3TObject* GetMeasuredSpectrum(TList *fContainer);
4TObject* GetGeneratedSpectrum(TList *fContainer);
5TObject* GetEfficiency(TList *fContainer);
6THnSparse* GetResponseMatrix(TList *fContainer);
7THnSparse* CreateGuessed(const THnSparse* h);
8
9void testUnfolding(const char *testfile = "HFEtask.root") {
10 Load();
11 TList *results = GetResults(testfile);
12 if(!results){
13 Error("No output objects: Calculation will terminate here");
14 return;
15 }
16 // get the essential
17 AliCFDataGrid *measured = (AliCFDataGrid*) GetMeasuredSpectrum(testfile);
18 AliCFDataGrid *generated = (AliCFDataGrid*) GetGeneratedSpectrum(testfile);
19 AliCFEffGrid *efficiency = (AliCFEffGrid*) GetEfficiency(testfile);
20 THnSparse *response = (THnSparse*) GetResponseMatrix(testfile);
21
22 // create a guessed "a priori" distribution using binning of MC
23 THnSparse* guessed = CreateGuessed(((AliCFGridSparse*)generated->GetData())->GetGrid()) ;
24 //----
25 AliCFUnfolding unfolding("unfolding","",3,response,efficiency->GetGrid(),((AliCFGridSparse*)measured->GetData())->GetGrid(),guessed);
26 unfolding.SetMaxNumberOfIterations(100);
27 unfolding.SetMaxChi2PerDOF(1.e-07);
28 //unfolding.UseSmoothing();
29 unfolding.Unfold();
30
31 THnSparse* result = unfolding.GetUnfolded();
32 //----
33
34 TCanvas * canvas = new TCanvas("canvas","",1000,700);
35 canvas->Divide(3,3);
36
37 canvas->cd(1);
38 TH2D* h_gen = generated->Project(0,1);
39 h_gen->SetTitle("generated");
40 //printf("c1\n");
41 h_gen->Draw("lego2");
42
43 canvas->cd(2);
44 TH2D* h_meas = measured->Project(0,1);
45 h_meas->SetTitle("measured");
46 //printf("c2\n");
47 h_meas->Draw("lego2");
48
49 canvas->cd(3);
50 TH2D* h_guessed = guessed->Projection(1,0);
51 h_guessed->SetTitle("a priori");
52 //printf("c3\n");
53 h_guessed->Draw("lego2");
54
55 canvas->cd(4);
56 TH2D* h_eff = efficiency->Project(0,1);
57 h_eff->SetTitle("efficiency");
58 //printf("c4\n");
59 h_eff->Draw("lego2");
60
61 canvas->cd(5);
62 TH2D* h_unf = result->Projection(1,0);
63 h_unf->SetTitle("unfolded");
64 //printf("c5\n");
65 h_unf->Draw("lego2");
66
67 canvas->cd(6);
68 TH2D* ratio = (TH2D*)h_unf->Clone();
69 ratio->SetTitle("unfolded / generated");
70 ratio->Divide(h_unf,h_gen,1,1);
71// ratio->Draw("cont4z");
72// ratio->Draw("surf2");
73 //printf("c6\n");
74 ratio->Draw("lego2");
75
76 canvas->cd(7);
77 TH2* orig = unfolding.GetOriginalPrior()->Projection(1,0);
78 orig->SetTitle("original prior");
79 //printf("c7\n");
80 orig->Draw("lego2");
81
82 canvas->cd(8);
83 AliCFDataGrid* corrected = (AliCFDataGrid*)measured->Clone();
84 //corrected->ApplyEffCorrection(*efficiency);
85 TH2D* corr = corrected->Project(0,1);
86 corr->SetTitle("simple correction");
87 //printf("c8\n");
88 corr->Draw("lego2");
89
90 canvas->cd(9);
91 TH2D* ratio2 = (TH2D*) corr->Clone();
92 ratio2->Divide(corr,h_gen,1,1);
93 ratio2->SetTitle("simple correction / generated");
94 //printf("c9\n");
95 ratio2->Draw("lego2");
96
97 return;
98}
99
100// ====================================================================
101
102
103void Load(){
104 gSystem->Load("libANALYSIS");
105 gSystem->Load("libCORRFW");
106 gStyle->SetPalette(1);
107 gStyle->SetOptStat(0);
108}
109
110TList *GetResults(const Char_t *testfile){
111 //
112 // read output
113 //
114 TFile *f = TFile::Open(testfile);
115 if(!f || f->IsZombie()){
116 Error("File not readable");
117 return 0x0;
118 }
119 TList *l = dynamic_cast<TList *>(f->Get("Results"));
120 if(l){
121 Error("Output container not found");
122 f->Close(); delete f;
123 return 0x0;
124 }
125 TList *returnlist = dynamic_cast<TList *>(l->Clone());
126 f->Close(); delete f;
127 return returnlist;
128}
129
130TObject* GetMeasuredSpectrum(TList *fContainer) {
131 AliCFContainer* c = (AliCFContainer*)fContainer->FindObject("container");
132 AliCFDataGrid* data = new AliCFDataGrid("data","",*c);
133 //data->SetMeasured(AliHFEcuts::kStepHFEcuts + 1);
134 data->SetMeasured(5);
135 return data;
136}
137
138TObject* GetGeneratedSpectrum(TList *fContainer) {
139 AliCFContainer* c = (AliCFContainer*)fContainer->FindObject("container");
140 AliCFDataGrid* data = new AliCFDataGrid("data","",*c);
141 //data->SetMeasured(AliHFEcuts::kStepMCGenerated);
142 data->SetMeasured(0);
143 return data;
144}
145
146TObject* GetEfficiency(TList *fContainer) {
147 AliCFContainer* c = (AliCFContainer*)fContainer->FindObject("container");
148 AliCFEffGrid* eff = new AliCFEffGrid("eff","",*c);
149 //eff->CalculateEfficiency(AliHFEcuts::kStepHFEcuts + 2,AliHFEcuts::kStepMCGenerated);
150 eff->CalculateEfficiency(6,0);
151 return eff;
152}
153
154THnSparse* GetResponseMatrix(TList *fContainer) {
155 THnSparse* h = (THnSparse*)f->Get("correlation");
156 return h;
157}
158
159THnSparse* CreateGuessed(const THnSparse* h) {
160 THnSparse* guessed = (THnSparse*) h->Clone();
161 return guessed ;
162}