AliPhysics  cdeda5a (cdeda5a)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DrawELossPoisson.C
Go to the documentation of this file.
1 
13 #ifndef __CINT__
14 # include <TH1.h>
15 # include <TH2.h>
16 # include <TList.h>
17 # include <TFile.h>
18 # include <TString.h>
19 # include <TError.h>
20 # include <TPad.h>
21 # include <TCanvas.h>
22 # include <TMath.h>
23 # include <TF1.h>
24 # include <TLine.h>
25 # include <TLatex.h>
26 # include <TLinearFitter.h>
27 # include <TStyle.h>
28 #else
29 class TList;
30 #endif
31 
49  Double_t xmin, Double_t xmax)
50 {
51  if (!p) return 0;
52 
53  TList* ring = static_cast<TList*>(p->FindObject(Form("FMD%d%c",d,r)));
54  if (!ring) {
55  Error("DrawELossPoisson", "List FMD%d%c not found in %s",d,r,p->GetName());
56  return 0;
57  }
58 
59  TH2* corr = static_cast<TH2D*>(ring->FindObject("elossVsPoisson"));
60  if (!corr) {
61  Error("DrawRingELossPoisson", "Histogram esdEloss not found in FMD%d%c",
62  d, r);
63  return 0;
64  }
65  TPad* pad = static_cast<TPad*>(gPad);
66  pad->SetGridy();
67  pad->SetGridx();
68  pad->SetLogz();
69  if (d == 3) {
70  pad->SetPad(pad->GetXlowNDC(), pad->GetYlowNDC(), .99,
71  pad->GetYlowNDC()+pad->GetHNDC());
72  pad->SetRightMargin(0.15);
73  }
74  pad->SetFillColor(0);
75  if (xmax < 0) xmax = corr->GetXaxis()->GetXmax();
76  corr->GetXaxis()->SetRangeUser(xmin,xmax);
77  corr->GetYaxis()->SetRangeUser(xmin,xmax);
78  corr->SetTitle(Form("FMD%d%c",d,r));
79  // Info("", "Entries: %d, integral: %f",
80  // int(corr->GetEntries()), corr->Integral());
81  corr->Draw("colz");
82  if (corr->GetEntries() <= 0) return 0;
83 
84  // Calculate the linear regression
85  // Double_t dx = (xmax-xmin);
86  Double_t rxy = corr->GetCorrelationFactor();
87  Double_t sx = corr->GetRMS(1);
88  Double_t sy = corr->GetRMS(2);
89  Double_t sx2 = sx*sx;
90  Double_t sy2 = sy*sy;
91  Double_t cxy = corr->GetCovariance();
92  Double_t mx = corr->GetMean(1);
93  Double_t my = corr->GetMean(2);
94  Double_t delta = 1;
95 #if 0
96  Double_t beta = rxy * sy / sx;
97 #else
98  if (TMath::Abs(cxy) < 1e-6) return 0;
99  Double_t beta = ((sy2 - delta*sx2 +
100  TMath::Sqrt(TMath::Power(sy2-delta*sx2,2) +
101  4*delta*cxy*cxy))
102  / 2 / cxy);
103 #endif
104  Double_t alpha = my - beta * mx;
105  TF1* f = new TF1("f", "[0]+[1]*x", xmin, xmax);
106  f->SetParameters(alpha, beta);
107  f->SetLineWidth(1);
108  f->SetLineStyle(1);
109  f->Draw("same");
110 
111  TLine* l = new TLine(xmin,xmin,xmax,xmax);
112  l->SetLineWidth(1);
113  l->SetLineStyle(2);
114  l->SetLineColor(kBlack);
115  l->Draw();
116  // corr->GetXaxis()->SetRangeUser(0, 2);
117 
118 #if 0
119  Info("", "FMD%d%c correlation coefficient: %9.5f%% "
120  "line y = %f + %f * x", d, r, 100*rxy, alpha, beta);
121 #endif
122 
123  Double_t x = pad->GetLeftMargin()+.01;
124  Double_t y = 1-pad->GetTopMargin()-.01;
125  TLatex* ltx = new TLatex(x, y, Form("FMD%d%c", d, r));
126  ltx->SetNDC();
127  ltx->SetTextAlign(13);
128  ltx->SetTextSize(0.08);
129  ltx->Draw();
130  y -= 0.12;
131  ltx->SetTextSize(0.06);
132  ltx->DrawLatex(x,y,"Deming regression: y=#alpha+#beta x");
133  x += .02;
134  y -= .06;
135  ltx->SetTextSize(0.05);
136  ltx->DrawLatex(x, y, Form("#alpha=%5.3f", alpha));
137  y -= .06;
138  ltx->DrawLatex(x, y, Form("#beta= %5.3f", beta));
139 
140  TLinearFitter* fitter = new TLinearFitter(1);
141  fitter->SetFormula("1 ++ x");
142  for (Int_t i = 1; i <= corr->GetNbinsX(); i++) {
143  x = corr->GetXaxis()->GetBinCenter(i);
144  if (x < xmin || x > xmax) continue;
145  for (Int_t j = 1; j <= corr->GetNbinsY(); j++) {
146  y = corr->GetYaxis()->GetBinCenter(j);
147  if (y < xmin || y > xmax) continue;
148  Double_t c = corr->GetBinContent(i,j);
149  if (c < .1) continue;
150  fitter->AddPoint(&x, y, c);
151  }
152  }
153  Bool_t robust = false;
154  if (robust)
155  fitter->EvalRobust();
156  else
157  fitter->Eval();
158 #if 0
159  for (Int_t i = 0; i < 2; i++) {
160  std::cout << i << "\t"
161  << fitter->GetParName(i) << "\t"
162  << fitter->GetParameter(i) << "\t";
163  if (!robust)
164  std::cout << fitter->GetParError(i);
165  std::cout << std::endl;
166  }
167  Double_t chi2 = fitter->GetChisquare();
168  Int_t ndf = (fitter->GetNpoints() -
169  fitter->GetNumberFreeParameters() );
170  std::cout << "chi2/ndf: " << chi2 << '/' << ndf
171  << '=' << chi2 / ndf << std::endl;
172 #endif
173 
174  // corr->Scale(1. / corr->GetMaximum());
175  pad->cd();
176  return corr->GetCorrelationFactor();
177 }
178 
192 void
193 DrawELossPoisson(const char* filename="forward.root",
194  const char* folder="ForwardResults",
195  Double_t xmax=-1,
196  Double_t xmin=-1)
197 {
198  gStyle->SetPalette(1);
199  gStyle->SetOptFit(0);
200  gStyle->SetOptStat(0);
201  gStyle->SetOptTitle(0);
202  gStyle->SetTitleW(.4);
203  gStyle->SetTitleH(.1);
204  gStyle->SetTitleX(.4);
205  // gStyle->SetTitleY(.1);
206  gStyle->SetTitleColor(0);
207  gStyle->SetTitleStyle(0);
208  gStyle->SetTitleBorderSize(0);
209 
210  TFile* file = TFile::Open(filename, "READ");
211  if (!file) {
212  Error("DrawELossPoisson", "failed to open %s", filename);
213  return;
214  }
215 
216  TList* forward = static_cast<TList*>(file->Get(folder));
217  if (!forward) {
218  Error("DrawELossPoisson", "List %s not found in %s", folder, filename);
219  return;
220  }
221 
222  TList* dc = static_cast<TList*>(forward->FindObject("fmdDensityCalculator"));
223  if (!dc) {
224  Error("DrawELossPoisson", "List fmdDensityCalculator not found in Forward");
225  return;
226  }
227 
228  TCanvas* c = new TCanvas("elossVsPoisson",
229  "N_ch from ELoss versus from Poisson", 900, 700);
230  c->SetFillColor(0);
231  c->SetBorderSize(0);
232  c->SetBorderMode(0);
233  c->SetHighLightColor(0);
234  c->Divide(3, 2, 0, 0);
235 
236  Double_t corrs[5];
237  c->cd(1); corrs[0] = DrawRingELossPoisson(dc, 1, 'I', xmin, xmax);
238  c->cd(2); corrs[1] = DrawRingELossPoisson(dc, 2, 'I', xmin, xmax);
239  c->cd(5); corrs[2] = DrawRingELossPoisson(dc, 2, 'O', xmin, xmax);
240  c->cd(3); corrs[3] = DrawRingELossPoisson(dc, 3, 'I', xmin, xmax);
241  c->cd(6); corrs[4] = DrawRingELossPoisson(dc, 3, 'O', xmin, xmax);
242 
243  TVirtualPad* p = c->cd(4);
244  p->SetTopMargin(0.05);
245  p->SetRightMargin(0.10);
246  p->SetLeftMargin(0.15);
247  p->SetBottomMargin(0.15);
248  p->SetFillColor(0);
249 
250  TH1D* hc = new TH1D("corrs", "Correlation factors", 5, .5, 5.5);
251  hc->SetFillColor(kRed+1);
252  hc->SetFillStyle(3001);
253  hc->SetMinimum(0.0);
254  hc->SetMaximum(1.3);
255  hc->GetXaxis()->SetBinLabel(1,"FMD1i"); hc->SetBinContent(1,corrs[0]);
256  hc->GetXaxis()->SetBinLabel(2,"FMD2i"); hc->SetBinContent(2,corrs[1]);
257  hc->GetXaxis()->SetBinLabel(3,"FMD2o"); hc->SetBinContent(3,corrs[2]);
258  hc->GetXaxis()->SetBinLabel(4,"FMD3i"); hc->SetBinContent(4,corrs[3]);
259  hc->GetXaxis()->SetBinLabel(5,"FMD3o"); hc->SetBinContent(5,corrs[4]);
260  hc->GetXaxis()->SetLabelSize(0.08);
261  hc->GetYaxis()->SetTitle("r");
262  hc->SetMarkerSize(1.5);
263  hc->Draw("text hist");
264 
265  // TH2D* highCuts = static_cast<TH2D*>(dc->FindObject("highCuts"));
266  // if (highCuts) highCuts->Draw("colz");
267  c->cd();
268  c->SaveAs("elossVsPoisson.png");
269 }
270 //
271 // EOF
272 //
const char * filename
Definition: TestFCM.C:1
double Double_t
Definition: External.C:58
char Char_t
Definition: External.C:18
TCanvas * c
Definition: TestFitELoss.C:172
int Int_t
Definition: External.C:63
Definition: External.C:228
Definition: External.C:212
Double_t DrawRingELossPoisson(TList *p, UShort_t d, Char_t r, Double_t xmin, Double_t xmax)
TList * fitter
Definition: DrawAnaELoss.C:26
void DrawELossPoisson(const char *filename="forward.root", const char *folder="ForwardResults", Double_t xmax=-1, Double_t xmin=-1)
Definition: External.C:220
TFile * file
unsigned short UShort_t
Definition: External.C:28
bool Bool_t
Definition: External.C:53