]>
Commit | Line | Data |
---|---|---|
1d86b7ca | 1 | /** |
2 | * @file DrawELossPoisson.C | |
3 | * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk> | |
4 | * @date Thu Jul 7 10:24:58 2011 | |
5 | * | |
6 | * @brief A script to draw the Poisson vs Energy Loss correlation | |
7 | * | |
56199f2b | 8 | * |
f0ef1e71 | 9 | * @ingroup pwg2_forward_scripts_qa |
1d86b7ca | 10 | * |
11 | */ | |
d015ecfe | 12 | #ifndef __CINT__ |
13 | # include <TH1.h> | |
14 | # include <TH2.h> | |
15 | # include <TList.h> | |
16 | # include <TFile.h> | |
17 | # include <TString.h> | |
18 | # include <TError.h> | |
19 | # include <TPad.h> | |
20 | # include <TCanvas.h> | |
21 | # include <TMath.h> | |
22 | # include <TF1.h> | |
23 | # include <TLine.h> | |
24 | # include <TLatex.h> | |
25 | # include <TLinearFitter.h> | |
26 | # include <TStyle.h> | |
27 | #else | |
28 | class TList; | |
29 | #endif | |
1d86b7ca | 30 | |
31 | /** | |
32 | * Draw the poisson @f$N_{ch}@f$ estimate against the @f$\Delta@f$ | |
33 | * @f$N_{ch}@f$ estimate and do a regression line between the two | |
34 | * | |
35 | * @param p List | |
36 | * @param d Detector | |
37 | * @param r Ring | |
38 | * @param xmin Minimum | |
39 | * @param xmax Maximum | |
40 | * | |
41 | * @return The regression coefficient | |
56199f2b | 42 | * |
f0ef1e71 | 43 | * @ingroup pwg2_forward_scripts_qa |
1d86b7ca | 44 | */ |
0919aa56 | 45 | Double_t |
797161e8 | 46 | DrawRingELossPoisson(TList* p, UShort_t d, Char_t r, |
47 | Double_t xmin, Double_t xmax) | |
0919aa56 | 48 | { |
d015ecfe | 49 | if (!p) return 0; |
0919aa56 | 50 | |
51 | TList* ring = static_cast<TList*>(p->FindObject(Form("FMD%d%c",d,r))); | |
52 | if (!ring) { | |
53 | Error("DrawELossPoisson", "List FMD%d%c not found in %s",d,r,p->GetName()); | |
d015ecfe | 54 | return 0; |
0919aa56 | 55 | } |
56 | ||
57 | TH2* corr = static_cast<TH2D*>(ring->FindObject("elossVsPoisson")); | |
58 | if (!corr) { | |
59 | Error("DrawRingELossPoisson", "Histogram esdEloss not found in FMD%d%c", | |
60 | d, r); | |
d015ecfe | 61 | return 0; |
0919aa56 | 62 | } |
b5c9a732 | 63 | TPad* pad = static_cast<TPad*>(gPad); |
64 | pad->SetGridy(); | |
65 | pad->SetGridx(); | |
66 | pad->SetLogz(); | |
67 | if (d == 3) { | |
68 | pad->SetPad(pad->GetXlowNDC(), pad->GetYlowNDC(), .99, | |
69 | pad->GetYlowNDC()+pad->GetHNDC()); | |
70 | pad->SetRightMargin(0.15); | |
71 | } | |
72 | pad->SetFillColor(0); | |
797161e8 | 73 | if (xmax < 0) xmax = corr->GetXaxis()->GetXmax(); |
74 | corr->GetXaxis()->SetRangeUser(xmin,xmax); | |
75 | corr->GetYaxis()->SetRangeUser(xmin,xmax); | |
0919aa56 | 76 | corr->SetTitle(Form("FMD%d%c",d,r)); |
d015ecfe | 77 | // Info("", "Entries: %d, integral: %f", |
78 | // int(corr->GetEntries()), corr->Integral()); | |
0919aa56 | 79 | corr->Draw("colz"); |
d015ecfe | 80 | if (corr->GetEntries() <= 0) return 0; |
0919aa56 | 81 | |
797161e8 | 82 | // Calculate the linear regression |
d015ecfe | 83 | // Double_t dx = (xmax-xmin); |
797161e8 | 84 | Double_t rxy = corr->GetCorrelationFactor(); |
85 | Double_t sx = corr->GetRMS(1); | |
86 | Double_t sy = corr->GetRMS(2); | |
87 | Double_t sx2 = sx*sx; | |
88 | Double_t sy2 = sy*sy; | |
89 | Double_t cxy = corr->GetCovariance(); | |
90 | Double_t mx = corr->GetMean(1); | |
91 | Double_t my = corr->GetMean(2); | |
92 | Double_t delta = 1; | |
93 | #if 0 | |
94 | Double_t beta = rxy * sy / sx; | |
95 | #else | |
d015ecfe | 96 | if (TMath::Abs(cxy) < 1e-6) return 0; |
797161e8 | 97 | Double_t beta = ((sy2 - delta*sx2 + |
98 | TMath::Sqrt(TMath::Power(sy2-delta*sx2,2) + | |
99 | 4*delta*cxy*cxy)) | |
100 | / 2 / cxy); | |
101 | #endif | |
102 | Double_t alpha = my - beta * mx; | |
103 | TF1* f = new TF1("f", "[0]+[1]*x", xmin, xmax); | |
104 | f->SetParameters(alpha, beta); | |
105 | f->SetLineWidth(1); | |
106 | f->SetLineStyle(1); | |
107 | f->Draw("same"); | |
108 | ||
109 | TLine* l = new TLine(xmin,xmin,xmax,xmax); | |
0919aa56 | 110 | l->SetLineWidth(1); |
797161e8 | 111 | l->SetLineStyle(2); |
0919aa56 | 112 | l->SetLineColor(kBlack); |
113 | l->Draw(); | |
114 | // corr->GetXaxis()->SetRangeUser(0, 2); | |
115 | ||
d015ecfe | 116 | #if 0 |
797161e8 | 117 | Info("", "FMD%d%c correlation coefficient: %9.5f%% " |
118 | "line y = %f + %f * x", d, r, 100*rxy, alpha, beta); | |
d015ecfe | 119 | #endif |
797161e8 | 120 | |
b5c9a732 | 121 | Double_t x = pad->GetLeftMargin()+.01; |
122 | Double_t y = 1-pad->GetTopMargin()-.01; | |
123 | TLatex* ltx = new TLatex(x, y, Form("FMD%d%c", d, r)); | |
124 | ltx->SetNDC(); | |
797161e8 | 125 | ltx->SetTextAlign(13); |
b5c9a732 | 126 | ltx->SetTextSize(0.08); |
797161e8 | 127 | ltx->Draw(); |
b5c9a732 | 128 | y -= 0.12; |
129 | ltx->SetTextSize(0.06); | |
130 | ltx->DrawLatex(x,y,"Deming regression: y=#alpha+#beta x"); | |
131 | x += .02; | |
132 | y -= .06; | |
797161e8 | 133 | ltx->SetTextSize(0.05); |
797161e8 | 134 | ltx->DrawLatex(x, y, Form("#alpha=%5.3f", alpha)); |
b5c9a732 | 135 | y -= .06; |
797161e8 | 136 | ltx->DrawLatex(x, y, Form("#beta= %5.3f", beta)); |
137 | ||
138 | TLinearFitter* fitter = new TLinearFitter(1); | |
139 | fitter->SetFormula("1 ++ x"); | |
140 | for (Int_t i = 1; i <= corr->GetNbinsX(); i++) { | |
141 | x = corr->GetXaxis()->GetBinCenter(i); | |
142 | if (x < xmin || x > xmax) continue; | |
143 | for (Int_t j = 1; j <= corr->GetNbinsY(); j++) { | |
144 | y = corr->GetYaxis()->GetBinCenter(j); | |
145 | if (y < xmin || y > xmax) continue; | |
146 | Double_t c = corr->GetBinContent(i,j); | |
147 | if (c < .1) continue; | |
148 | fitter->AddPoint(&x, y, c); | |
149 | } | |
150 | } | |
151 | Bool_t robust = false; | |
152 | if (robust) | |
153 | fitter->EvalRobust(); | |
154 | else | |
155 | fitter->Eval(); | |
d015ecfe | 156 | #if 0 |
797161e8 | 157 | for (Int_t i = 0; i < 2; i++) { |
158 | std::cout << i << "\t" | |
159 | << fitter->GetParName(i) << "\t" | |
160 | << fitter->GetParameter(i) << "\t"; | |
161 | if (!robust) | |
162 | std::cout << fitter->GetParError(i); | |
163 | std::cout << std::endl; | |
164 | } | |
165 | Double_t chi2 = fitter->GetChisquare(); | |
166 | Int_t ndf = (fitter->GetNpoints() - | |
167 | fitter->GetNumberFreeParameters() ); | |
b5c9a732 | 168 | std::cout << "chi2/ndf: " << chi2 << '/' << ndf |
169 | << '=' << chi2 / ndf << std::endl; | |
d015ecfe | 170 | #endif |
797161e8 | 171 | |
b5c9a732 | 172 | // corr->Scale(1. / corr->GetMaximum()); |
173 | pad->cd(); | |
0919aa56 | 174 | return corr->GetCorrelationFactor(); |
175 | } | |
176 | ||
1d86b7ca | 177 | /** |
178 | * Draw the correlation between the Poisson @f$N_{ch}@f$ estimate | |
179 | * and the @f$\Delta@f$ @f$N_{ch}@f$ estimate and do a regression | |
180 | * line between the two for each ring | |
181 | * | |
182 | * @param filename File to read | |
183 | * @param xmax Minimum X | |
184 | * @param xmin Maximum X | |
56199f2b | 185 | * |
f0ef1e71 | 186 | * @ingroup pwg2_forward_scripts_qa |
1d86b7ca | 187 | */ |
0919aa56 | 188 | void |
797161e8 | 189 | DrawELossPoisson(const char* filename="forward.root", |
d015ecfe | 190 | const char* folder="ForwardResults", |
797161e8 | 191 | Double_t xmax=-1, |
b5c9a732 | 192 | Double_t xmin=-1) |
0919aa56 | 193 | { |
194 | gStyle->SetPalette(1); | |
195 | gStyle->SetOptFit(0); | |
196 | gStyle->SetOptStat(0); | |
b5c9a732 | 197 | gStyle->SetOptTitle(0); |
0919aa56 | 198 | gStyle->SetTitleW(.4); |
199 | gStyle->SetTitleH(.1); | |
b5c9a732 | 200 | gStyle->SetTitleX(.4); |
201 | // gStyle->SetTitleY(.1); | |
0919aa56 | 202 | gStyle->SetTitleColor(0); |
203 | gStyle->SetTitleStyle(0); | |
204 | gStyle->SetTitleBorderSize(0); | |
0919aa56 | 205 | |
206 | TFile* file = TFile::Open(filename, "READ"); | |
207 | if (!file) { | |
208 | Error("DrawELossPoisson", "failed to open %s", filename); | |
209 | return; | |
210 | } | |
211 | ||
d015ecfe | 212 | TList* forward = static_cast<TList*>(file->Get(folder)); |
0919aa56 | 213 | if (!forward) { |
d015ecfe | 214 | Error("DrawELossPoisson", "List %s not found in %s", folder, filename); |
0919aa56 | 215 | return; |
216 | } | |
217 | ||
218 | TList* dc = static_cast<TList*>(forward->FindObject("fmdDensityCalculator")); | |
219 | if (!dc) { | |
220 | Error("DrawELossPoisson", "List fmdDensityCalculator not found in Forward"); | |
221 | return; | |
222 | } | |
223 | ||
224 | TCanvas* c = new TCanvas("elossVsPoisson", | |
225 | "N_ch from ELoss versus from Poisson", 900, 700); | |
226 | c->SetFillColor(0); | |
227 | c->SetBorderSize(0); | |
228 | c->SetBorderMode(0); | |
229 | c->SetHighLightColor(0); | |
230 | c->Divide(3, 2, 0, 0); | |
231 | ||
232 | Double_t corrs[5]; | |
797161e8 | 233 | c->cd(1); corrs[0] = DrawRingELossPoisson(dc, 1, 'I', xmin, xmax); |
234 | c->cd(2); corrs[1] = DrawRingELossPoisson(dc, 2, 'I', xmin, xmax); | |
235 | c->cd(5); corrs[2] = DrawRingELossPoisson(dc, 2, 'O', xmin, xmax); | |
236 | c->cd(3); corrs[3] = DrawRingELossPoisson(dc, 3, 'I', xmin, xmax); | |
237 | c->cd(6); corrs[4] = DrawRingELossPoisson(dc, 3, 'O', xmin, xmax); | |
0919aa56 | 238 | |
239 | TVirtualPad* p = c->cd(4); | |
240 | p->SetTopMargin(0.05); | |
241 | p->SetRightMargin(0.10); | |
242 | p->SetLeftMargin(0.15); | |
243 | p->SetBottomMargin(0.15); | |
244 | p->SetFillColor(0); | |
245 | ||
246 | TH1D* hc = new TH1D("corrs", "Correlation factors", 5, .5, 5.5); | |
247 | hc->SetFillColor(kRed+1); | |
248 | hc->SetFillStyle(3001); | |
249 | hc->SetMinimum(0.0); | |
797161e8 | 250 | hc->SetMaximum(1.3); |
0919aa56 | 251 | hc->GetXaxis()->SetBinLabel(1,"FMD1i"); hc->SetBinContent(1,corrs[0]); |
252 | hc->GetXaxis()->SetBinLabel(2,"FMD2i"); hc->SetBinContent(2,corrs[1]); | |
253 | hc->GetXaxis()->SetBinLabel(3,"FMD2o"); hc->SetBinContent(3,corrs[2]); | |
254 | hc->GetXaxis()->SetBinLabel(4,"FMD3i"); hc->SetBinContent(4,corrs[3]); | |
255 | hc->GetXaxis()->SetBinLabel(5,"FMD3o"); hc->SetBinContent(5,corrs[4]); | |
256 | hc->GetXaxis()->SetLabelSize(0.08); | |
257 | hc->GetYaxis()->SetTitle("r"); | |
258 | hc->SetMarkerSize(1.5); | |
259 | hc->Draw("text hist"); | |
260 | ||
261 | // TH2D* highCuts = static_cast<TH2D*>(dc->FindObject("highCuts")); | |
262 | // if (highCuts) highCuts->Draw("colz"); | |
263 | c->cd(); | |
797161e8 | 264 | c->SaveAs("elossVsPoisson.png"); |
0919aa56 | 265 | } |
1d86b7ca | 266 | // |
267 | // EOF | |
268 | // |