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);
9 void testUnfolding(const char *testfile = "HFEtask.root") {
11 TList *results = GetResults(testfile);
13 Error("No output objects: Calculation will terminate here");
17 AliCFDataGrid *measured = (AliCFDataGrid*) GetMeasuredSpectrum(testfile);
18 AliCFDataGrid *generated = (AliCFDataGrid*) GetGeneratedSpectrum(testfile);
19 AliCFEffGrid *efficiency = (AliCFEffGrid*) GetEfficiency(testfile);
20 THnSparse *response = (THnSparse*) GetResponseMatrix(testfile);
22 // create a guessed "a priori" distribution using binning of MC
23 THnSparse* guessed = CreateGuessed(((AliCFGridSparse*)generated->GetData())->GetGrid()) ;
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();
31 THnSparse* result = unfolding.GetUnfolded();
34 TCanvas * canvas = new TCanvas("canvas","",1000,700);
38 TH2D* h_gen = generated->Project(0,1);
39 h_gen->SetTitle("generated");
44 TH2D* h_meas = measured->Project(0,1);
45 h_meas->SetTitle("measured");
47 h_meas->Draw("lego2");
50 TH2D* h_guessed = guessed->Projection(1,0);
51 h_guessed->SetTitle("a priori");
53 h_guessed->Draw("lego2");
56 TH2D* h_eff = efficiency->Project(0,1);
57 h_eff->SetTitle("efficiency");
62 TH2D* h_unf = result->Projection(1,0);
63 h_unf->SetTitle("unfolded");
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");
77 TH2* orig = unfolding.GetOriginalPrior()->Projection(1,0);
78 orig->SetTitle("original prior");
83 AliCFDataGrid* corrected = (AliCFDataGrid*)measured->Clone();
84 //corrected->ApplyEffCorrection(*efficiency);
85 TH2D* corr = corrected->Project(0,1);
86 corr->SetTitle("simple correction");
91 TH2D* ratio2 = (TH2D*) corr->Clone();
92 ratio2->Divide(corr,h_gen,1,1);
93 ratio2->SetTitle("simple correction / generated");
95 ratio2->Draw("lego2");
100 // ====================================================================
104 gSystem->Load("libANALYSIS");
105 gSystem->Load("libCORRFW");
106 gStyle->SetPalette(1);
107 gStyle->SetOptStat(0);
110 TList *GetResults(const Char_t *testfile){
114 TFile *f = TFile::Open(testfile);
115 if(!f || f->IsZombie()){
116 Error("File not readable");
119 TList *l = dynamic_cast<TList *>(f->Get("Results"));
121 Error("Output container not found");
122 f->Close(); delete f;
125 TList *returnlist = dynamic_cast<TList *>(l->Clone());
126 f->Close(); delete f;
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);
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);
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);
154 THnSparse* GetResponseMatrix(TList *fContainer) {
155 THnSparse* h = (THnSparse*)f->Get("correlation");
159 THnSparse* CreateGuessed(const THnSparse* h) {
160 THnSparse* guessed = (THnSparse*) h->Clone();