Small update to distinguish, in the HV plot, channels which are OFF, from channels...
[u/mrichter/AliRoot.git] / PWG3 / hfe / testUnfolding.C
1 void Load();
2 TList *GetResults(const Char_t *testfile);
3 TObject* GetMeasuredSpectrum(TList *fContainer);
4 TObject* GetGeneratedSpectrum(TList *fContainer);
5 TObject* GetEfficiency(TList *fContainer);
6 THnSparse* GetResponseMatrix(TList *fContainer);
7 THnSparse* CreateGuessed(const THnSparse* h);
8
9 void 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
103 void Load(){
104   gSystem->Load("libANALYSIS");
105   gSystem->Load("libCORRFW");
106   gStyle->SetPalette(1);
107   gStyle->SetOptStat(0);
108 }
109
110 TList *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
130 TObject* 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
138 TObject* 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
146 TObject* 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
154 THnSparse* GetResponseMatrix(TList *fContainer) {
155   THnSparse* h = (THnSparse*)f->Get("correlation");
156   return h;
157 }
158
159 THnSparse* CreateGuessed(const THnSparse* h) {
160   THnSparse* guessed = (THnSparse*) h->Clone();
161   return guessed ;
162 }