]>
Commit | Line | Data |
---|---|---|
c0b10ad4 | 1 | |
2 | void testUnfolding() { | |
fb494025 | 3 | TBenchmark b; |
4 | ||
5 | b.Start("init"); | |
6 | ||
c0b10ad4 | 7 | Load(); |
8 | ||
fb494025 | 9 | AliLog::SetGlobalDebugLevel(0); |
10 | ||
c0b10ad4 | 11 | // get the essential |
12 | AliCFDataGrid *measured = (AliCFDataGrid*) GetMeasuredSpectrum(); | |
13 | AliCFDataGrid *generated = (AliCFDataGrid*) GetGeneratedSpectrum(); | |
14 | AliCFEffGrid *efficiency = (AliCFEffGrid*) GetEfficiency(); | |
15 | THnSparse *response = (THnSparse*) GetResponseMatrix(); | |
16 | ||
17 | // create a guessed "a priori" distribution using binning of MC | |
18 | THnSparse* guessed = CreateGuessed(((AliCFGridSparse*)generated->GetData())->GetGrid()) ; | |
19 | //---- | |
fb494025 | 20 | |
21 | TF2* fit2D = new TF2("fit2D","[0]*([2]+[3]*y+[4]*y*y)*x*TMath::Exp(-x/[1])",0.1,0.8,-1.5,1.5); | |
22 | fit2D->SetParameter(0,1); | |
23 | fit2D->SetParameter(1,1); | |
24 | fit2D->SetParameter(2,1); | |
25 | fit2D->SetParameter(3,1); | |
26 | fit2D->SetParameter(4,1); | |
27 | fit2D->SetParLimits(1,0,1); | |
28 | ||
29 | AliCFUnfolding unfolding("unfolding","",2,response,efficiency->GetGrid(),((AliCFGridSparse*)measured->GetData())->GetGrid()/*,guessed*/); | |
c0b10ad4 | 30 | unfolding.SetMaxNumberOfIterations(100); |
31 | unfolding.SetMaxChi2PerDOF(1.e-07); | |
fb494025 | 32 | unfolding.UseSmoothing(fit2D,"ren"); |
c0b10ad4 | 33 | //unfolding.UseSmoothing(); |
fb494025 | 34 | |
35 | b.Stop("init"); | |
36 | b.Start("unfolding"); | |
c0b10ad4 | 37 | unfolding.Unfold(); |
fb494025 | 38 | b.Stop("unfolding"); |
39 | b.Start("finish"); | |
c0b10ad4 | 40 | |
fb494025 | 41 | canvas->cd(); |
c0b10ad4 | 42 | |
c0b10ad4 | 43 | TH2D* h_gen = generated->Project(0,1); |
44 | h_gen->SetTitle("generated"); | |
45 | h_gen->Draw("lego2"); | |
fb494025 | 46 | canvas->SaveAs("/tmp/generated.eps"); |
47 | canvas->SaveAs("/tmp/generated.gif"); | |
c0b10ad4 | 48 | |
c0b10ad4 | 49 | TH2D* h_meas = measured->Project(0,1); |
50 | h_meas->SetTitle("measured"); | |
51 | h_meas->Draw("lego2"); | |
fb494025 | 52 | canvas->SaveAs("/tmp/measured.eps"); |
53 | canvas->SaveAs("/tmp/measured.gif"); | |
54 | ||
55 | TH2D* h_guessed = unfolding.GetOriginalPrior()->Projection(1,0); | |
c0b10ad4 | 56 | h_guessed->SetTitle("a priori"); |
57 | h_guessed->Draw("lego2"); | |
fb494025 | 58 | canvas->SaveAs("/tmp/apriori.eps"); |
59 | canvas->SaveAs("/tmp/apriori.gif"); | |
c0b10ad4 | 60 | |
c0b10ad4 | 61 | TH2D* h_eff = efficiency->Project(0,1); |
62 | h_eff->SetTitle("efficiency"); | |
fb494025 | 63 | h_eff->Draw("e"); |
c0b10ad4 | 64 | |
fb494025 | 65 | TH2D* h_unf = unfolding.GetUnfolded()->Projection(1,0); |
c0b10ad4 | 66 | h_unf->SetTitle("unfolded"); |
67 | h_unf->Draw("lego2"); | |
fb494025 | 68 | fit2D->Draw("surf same"); |
69 | return; | |
70 | canvas->SaveAs("/tmp/unfolded.eps"); | |
71 | canvas->SaveAs("/tmp/unfolded.gif"); | |
72 | ||
c0b10ad4 | 73 | |
c0b10ad4 | 74 | TH2D* ratio = (TH2D*)h_unf->Clone(); |
75 | ratio->SetTitle("unfolded / generated"); | |
76 | ratio->Divide(h_unf,h_gen,1,1); | |
fb494025 | 77 | ratio->GetZaxis()->SetRangeUser(0.5,1.5); |
78 | //ratio->Draw("cont4z"); | |
c0b10ad4 | 79 | ratio->Draw("lego2"); |
fb494025 | 80 | //ratio->Draw("e"); |
81 | // ratio->ProjectionY()->Draw(); | |
82 | canvas->SaveAs("/tmp/ratio.eps"); | |
83 | canvas->SaveAs("/tmp/ratio.gif"); | |
84 | ||
c0b10ad4 | 85 | |
86 | canvas->cd(7); | |
87 | TH2* orig = unfolding.GetOriginalPrior()->Projection(1,0); | |
88 | orig->SetTitle("original prior"); | |
89 | orig->Draw("lego2"); | |
90 | ||
91 | canvas->cd(8); | |
fb494025 | 92 | TH2D* h_est = (TH2D*) unfolding.GetEstMeasured()->Projection(1,0); |
93 | h_est->SetTitle("est. measured"); | |
94 | h_est->Draw("e"); | |
c0b10ad4 | 95 | |
96 | canvas->cd(9); | |
fb494025 | 97 | unfolding.GetUnfolded()->Projection(0)->Draw(); |
c0b10ad4 | 98 | |
fb494025 | 99 | canvas->cd(10); |
100 | unfolding.GetUnfolded()->Projection(1)->Draw(); | |
101 | ||
102 | canvas->cd(11); | |
103 | h_gen->ProjectionX()->Draw(); | |
104 | ||
105 | canvas->cd(12); | |
106 | h_gen->ProjectionY()->Draw(); | |
107 | ||
108 | b.Stop("finish"); | |
109 | Float_t x,y; | |
110 | b.Summary(x,y); | |
c0b10ad4 | 111 | } |
112 | ||
113 | // ==================================================================== | |
114 | ||
115 | ||
116 | void Load(){ | |
117 | gSystem->Load("libANALYSIS"); | |
118 | gSystem->Load("libCORRFW"); | |
119 | gStyle->SetPalette(1); | |
120 | gStyle->SetOptStat(0); | |
fb494025 | 121 | TCanvas * canvas = new TCanvas("canvas","",600,400); |
122 | canvas->SetFillColor(10); | |
123 | canvas->SetFrameFillColor(10); | |
c0b10ad4 | 124 | } |
125 | ||
126 | TObject* GetMeasuredSpectrum() { | |
127 | TFile * f = TFile::Open("test/output.root","read"); | |
128 | AliCFContainer* c = (AliCFContainer*)f->Get("container"); | |
129 | AliCFDataGrid* data = new AliCFDataGrid("data","",*c); | |
130 | data->SetMeasured(1); | |
131 | return data; | |
132 | } | |
133 | TObject* GetGeneratedSpectrum() { | |
134 | TFile * f = TFile::Open("test/output.root","read"); | |
135 | AliCFContainer* c = (AliCFContainer*)f->Get("container"); | |
136 | AliCFDataGrid* data = new AliCFDataGrid("data","",*c); | |
137 | data->SetMeasured(0); | |
138 | return data; | |
139 | } | |
140 | TObject* GetEfficiency() { | |
141 | TFile * f = TFile::Open("test/output.root","read"); | |
142 | AliCFContainer* c = (AliCFContainer*)f->Get("container"); | |
143 | AliCFEffGrid* eff = new AliCFEffGrid("eff","",*c); | |
144 | eff->CalculateEfficiency(2,0); | |
145 | return eff; | |
146 | } | |
147 | THnSparse* GetResponseMatrix() { | |
148 | TFile * f = TFile::Open("test/output.root","read"); | |
149 | THnSparse* h = (THnSparse*)f->Get("correlation"); | |
150 | return h; | |
151 | } | |
152 | THnSparse* CreateGuessed(const THnSparse* h) { | |
153 | THnSparse* guessed = (THnSparse*) h->Clone(); | |
154 | return guessed ; | |
155 | } |