PWGLF/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 pwglf_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 pwglf_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 pwglf_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 //