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