Major overhaul of the QA code.
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / qa / DrawELossPoisson.C
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  * 
8  *
9  * @deprecated Use QATrender instead 
10  * @ingroup pwg2_forward_scripts_qa
11  * 
12  */
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
32 /** 
33  * Draw the poisson @f$N_{ch}@f$ estimate against the @f$\Delta@f$
34  * @f$N_{ch}@f$ estimate and do a regression line between the two 
35  * 
36  * @param p            List 
37  * @param d            Detector
38  * @param r            Ring
39  * @param xmin         Minimum
40  * @param xmax         Maximum
41  * 
42  * @return The regression coefficient 
43  *
44  * @deprecated Use QATrender instead 
45  * @ingroup pwg2_forward_scripts_qa
46  */
47 Double_t
48 DrawRingELossPoisson(TList* p, UShort_t d, Char_t r, 
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
179 /** 
180  * Draw the correlation between the Poisson @f$N_{ch}@f$ estimate
181  * and the @f$\Delta@f$ @f$N_{ch}@f$ estimate and do a regression
182  * line between the two for each ring
183  * 
184  * @param filename File to read
185  * @param xmax     Minimum X
186  * @param xmin     Maximum X 
187  *
188  * @deprecated Use QATrender instead 
189  * @ingroup pwg2_forward_scripts_qa
190  */
191 void
192 DrawELossPoisson(const char* filename="forward.root", 
193                  const char* folder="ForwardResults",
194                  Double_t xmax=-1,
195                  Double_t xmin=-1)
196 {
197   gStyle->SetPalette(1);
198   gStyle->SetOptFit(0);
199   gStyle->SetOptStat(0);
200   gStyle->SetOptTitle(0);
201   gStyle->SetTitleW(.4);
202   gStyle->SetTitleH(.1);
203   gStyle->SetTitleX(.4);
204   // gStyle->SetTitleY(.1);
205   gStyle->SetTitleColor(0);
206   gStyle->SetTitleStyle(0);
207   gStyle->SetTitleBorderSize(0);
208   
209   TFile* file = TFile::Open(filename, "READ");
210   if (!file) { 
211     Error("DrawELossPoisson", "failed to open %s", filename);
212     return;
213   }
214
215   TList* forward = static_cast<TList*>(file->Get(folder));
216   if (!forward) { 
217     Error("DrawELossPoisson", "List %s not found in %s", folder, filename);
218     return;
219   }
220
221   TList* dc = static_cast<TList*>(forward->FindObject("fmdDensityCalculator"));
222   if (!dc) { 
223     Error("DrawELossPoisson", "List fmdDensityCalculator not found in Forward");
224     return;
225   }
226   
227   TCanvas* c = new TCanvas("elossVsPoisson", 
228                            "N_ch from ELoss versus from Poisson", 900, 700);
229   c->SetFillColor(0);
230   c->SetBorderSize(0);
231   c->SetBorderMode(0);
232   c->SetHighLightColor(0);
233   c->Divide(3, 2, 0, 0);
234
235   Double_t corrs[5];
236   c->cd(1); corrs[0] = DrawRingELossPoisson(dc, 1, 'I', xmin, xmax);
237   c->cd(2); corrs[1] = DrawRingELossPoisson(dc, 2, 'I', xmin, xmax);
238   c->cd(5); corrs[2] = DrawRingELossPoisson(dc, 2, 'O', xmin, xmax);
239   c->cd(3); corrs[3] = DrawRingELossPoisson(dc, 3, 'I', xmin, xmax);
240   c->cd(6); corrs[4] = DrawRingELossPoisson(dc, 3, 'O', xmin, xmax);
241
242   TVirtualPad* p = c->cd(4);
243   p->SetTopMargin(0.05);
244   p->SetRightMargin(0.10);
245   p->SetLeftMargin(0.15);
246   p->SetBottomMargin(0.15);
247   p->SetFillColor(0);
248
249   TH1D* hc = new TH1D("corrs", "Correlation factors", 5, .5, 5.5);
250   hc->SetFillColor(kRed+1);
251   hc->SetFillStyle(3001);
252   hc->SetMinimum(0.0);
253   hc->SetMaximum(1.3);
254   hc->GetXaxis()->SetBinLabel(1,"FMD1i"); hc->SetBinContent(1,corrs[0]);
255   hc->GetXaxis()->SetBinLabel(2,"FMD2i"); hc->SetBinContent(2,corrs[1]);
256   hc->GetXaxis()->SetBinLabel(3,"FMD2o"); hc->SetBinContent(3,corrs[2]);
257   hc->GetXaxis()->SetBinLabel(4,"FMD3i"); hc->SetBinContent(4,corrs[3]);
258   hc->GetXaxis()->SetBinLabel(5,"FMD3o"); hc->SetBinContent(5,corrs[4]);
259   hc->GetXaxis()->SetLabelSize(0.08);
260   hc->GetYaxis()->SetTitle("r");
261   hc->SetMarkerSize(1.5);
262   hc->Draw("text hist");
263
264   // TH2D* highCuts = static_cast<TH2D*>(dc->FindObject("highCuts"));
265   // if (highCuts) highCuts->Draw("colz");
266   c->cd();
267   c->SaveAs("elossVsPoisson.png");
268 }
269 //
270 // EOF
271 //