PWG2/FORWARD/analysis2/corrs/DrawAnaELoss.C

   1 /**
   2  * Script to draw the energy loss fits from the output file of
   3  * AliFMDELossFitter(Task).
   4  *
   5  *
   6  * @ingroup pwg2_forward_scripts_corr
   7  */
   8 #include <TFile.h>
   9 #include <THStack.h>
  10 #include <TList.h>
  11 #include <TError.h>
  12 #include <TCanvas.h>
  13 #include <TPad.h>
  14 #include <TStyle.h>
  15 #include <TF1.h>
  16 #include <TLegend.h>
  17 #include <TMath.h>
  18
  19 //____________________________________________________________________
  20 // Global
  21 TList* fitter = 0;
  22 TCanvas* canvas = 0;
  23 const char* pdfName = "FitResults.pdf";
  24
  25 //____________________________________________________________________
  26 /**
  27  * Open a file.  The file is expected to contain the directory
  28  * structure
  29  *
  30  * @verbatim
  31  *  file
  32  *   +- ForwardResults
  33  *       +- fmdEnergyFitter
  34  *           +- chi2   (THStack)
  35  *           +- c      (THStack)
  36  *           +- delta  (THStack)
  37  *           +- xi     (THStack)
  38  *           +- sigma  (THStack)
  39  *           +- sigman (THStack)
  40  *           +- n      (THStack)
  41  *           +- a2     (THStack)
  42  *           +- ...    (THStack)
  43  *           +- an     (THStack)
  44  *           +- FMD1I (TList)
  45  *           |   +- FMD1I_edist (TH1)
  46  *           |   +- EDists      (TList)
  47  *           ...
  48  * @endverbatim
  49  *
  50  * @param fname File to open
  51  *
  52  * @return Pointer to the list of objects
  53  *
  54  * @ingroup pwg2_forward_scripts_corr
  55  */
  56 TList* OpenFile(const char* fname)
  57 {
  58   TFile* file = TFile::Open(fname, "READ");
  59   if (!file) {
  60     Error("DrawFits", "Couldn't open %s", fname);
  61     return 0;
  62   }
  63
  64   TList* forward = static_cast<TList*>(file->Get("ForwardResults"));
  65   // static_cast<TList*>(file->Get("PWG2forwardDnDeta/Forward"));
  66   if (!forward) {
  67     Error("DrawFits", "Couldn't get forward list from %s", fname);
  68     return 0;
  69   }
  70
  71   fitter = static_cast<TList*>(forward->FindObject("fmdEnergyFitter"));
  72   if (!fitter) {
  73     Error("DrawFits", "Couldn't get fitter folder");
  74     return 0;
  75   }
  76   return fitter;
  77 }
  78 //____________________________________________________________________
  79 /**
  80  * Open file if not done already
  81  *
  82  * @param fname File to open
  83  *
  84  * @return List of fits
  85  *
  86  * @ingroup pwg2_forward_scripts_corr
  87  */
  88 TList* CheckFitter(const char* fname="AnalysisResults.root")
  89 {
  90   if (!fitter) return OpenFile(fname);
  91   return fitter;
  92 }
  93
  94 //____________________________________________________________________
  95 /**
  96  * Make canvas if not done already
  97  *
  98  *
  99  * @return Canvas
 100  *
 101  * @ingroup pwg2_forward_scripts_corr
 102  */
 103 TCanvas* CheckCanvas()
 104 {
 105   if (canvas) return canvas;
 106
 107   gStyle->SetOptFit(111111);
 108   gStyle->SetTitleFillColor(0);
 109   gStyle->SetTitleBorderSize(0);
 110   gStyle->SetTitleX(0);
 111   gStyle->SetTitleY(.9);
 112   gStyle->SetTitleW(.4);
 113   gStyle->SetTitleH(.1);
 114   gStyle->SetStatW(0.2);
 115   gStyle->SetStatH(0.2);
 116   gStyle->SetStatColor(0);
 117   gStyle->SetStatBorderSize(1);
 118   gStyle->SetOptTitle(0);
 119   gStyle->SetOptFit(1111);
 120   gStyle->SetStatW(0.3);
 121   gStyle->SetStatH(0.15);
 122   gStyle->SetStatColor(0);
 123   gStyle->SetStatBorderSize(1);
 124
 125   canvas = new TCanvas("c", "C", Int_t(800 / TMath::Sqrt(2)), 800);
 126   canvas->SetFillColor(0);
 127   canvas->SetBorderMode(0);
 128   canvas->SetBorderSize(0);
 129   canvas->SetBottomMargin(0.15);
 130   return canvas;
 131 }
 132
 133 //____________________________________________________________________
 134 void CleanStack(THStack* stack)
 135 {
 136   TIter next(stack->GetHists());
 137   TObject* o = 0;
 138   while ((o = next())) {
 139     TString name(o->GetName());
 140     if (name.Contains("_t_"))
 141       stack->RecursiveRemove(o);
 142   }
 143 }
 144
 145
 146 //____________________________________________________________________
 147 THStack*
 148 AddToStack(TList* stacks, TList* fitter, const char* name)
 149 {
 150   TObject* o = fitter->FindObject(name);
 151   if (!o) {
 152     Warning("AddToStack", "Object %s not found in %s", name,
 153             fitter->GetName());
 154     // fitter->ls();
 155     return 0;
 156   }
 157   THStack* toAdd = static_cast<THStack*>(o);
 158   CleanStack(toAdd);
 159   Info("AddToStack", "Adding %s to stacks", name);
 160   stacks->Add(toAdd);
 161   return toAdd;
 162 }
 163
 164
 165 //____________________________________________________________________
 166 /**
 167  * Draw summary
 168  *
 169  * @param fname
 170  *
 171  * @ingroup pwg2_forward_scripts_corr
 172  */
 173 void DrawSummary(const char* fname="forward_eloss.root")
 174 {
 175   if (!CheckFitter(fname)) {
 176     Error("DrawFits", "File not opened");
 177     return;
 178   }
 179   if (!CheckCanvas()) {
 180     Error("DrawFits", "No canvas");
 181     return;
 182   }
 183   canvas->Clear();
 184
 185   TList    stacks;
 186   THStack* chi2nu  = AddToStack(&stacks, fitter, "chi2");
 187   THStack* c       = AddToStack(&stacks, fitter, "c");
 188   THStack* delta   = AddToStack(&stacks, fitter, "delta");
 189   THStack* xi      = AddToStack(&stacks, fitter, "xi");
 190   THStack* sigma   = AddToStack(&stacks, fitter, "sigma");
 191   THStack* sigman  = AddToStack(&stacks, fitter, "sigman");
 192   THStack* n       = AddToStack(&stacks, fitter, "n");
 193   Int_t    baseA   = stacks.GetEntries()+1;
 194   Int_t    i       = 2;
 195   while (true) {
 196     if (!AddToStack(&stacks, fitter, Form("a%d",i++)))
 197       break;
 198   }
 199   // stacks.ls();
 200   Int_t nMax = stacks.GetEntries();
 201   for (Int_t i = nMax-1; i >= baseA; i--) {
 202     THStack* stack   = static_cast<THStack*>(stacks.At(i));
 203     TIter    nextH(stack->GetHists());
 204     TH1*     hist    = 0;
 205     Bool_t   hasData = kFALSE;
 206     while ((hist = static_cast<TH1*>(nextH())))
 207       if (hist->Integral() > 0) hasData = kTRUE;
 208     if (!hasData) nMax--;
 209   }
 210
 211   canvas->SetRightMargin(0.01);
 212   canvas->SetTopMargin(0.01);
 213   Int_t nX = 2;
 214   Int_t nY = (nMax+1) / 2;
 215   canvas->Divide(nX, nY, 0.1, 0, 0);
 216
 217   TIter next(&stacks);
 218   THStack* stack = 0;
 219   i = 0;
 220   Int_t b = 1;
 221   while ((stack = static_cast<THStack*>(next()))) {
 222     if (i > nMax) break;
 223     Int_t ipad = 1+i/nY + 2 * (i % nY);
 224     Info("DrawSummary", "cd'ing to canvas %d for %s", ipad,
 225          stack->GetName());
 226     TVirtualPad* p = canvas->cd(ipad);
 227     p->SetLeftMargin(.6/nY);
 228     p->SetTopMargin(.01);
 229     p->SetRightMargin(.01);
 230     p->SetFillColor(0);
 231     p->SetFillStyle(0);
 232     p->SetGridx();
 233     stack->Draw("nostack");
 234     stack->GetHistogram()->SetYTitle(stack->GetTitle());
 235     stack->GetHistogram()->SetXTitle("#eta");
 236
 237     TAxis* yaxis = stack->GetHistogram()->GetYaxis();
 238     if (i == 0)                     yaxis->SetRangeUser(0,20); // Chi2
 239     if (i == 1)                     stack->SetMaximum(1);      // c
 240     if (i == 2)                     stack->SetMaximum(1);      // delta
 241     if (i == 3)                     stack->SetMaximum(0.1);   // xi
 242     if (i == 4 || i == 5)           stack->SetMaximum(0.5);    // sigma{,n}
 243     if (i == 7)                     stack->SetMaximum(0.5);    // a
 244     if (i == 0) p->SetLogy();
 245     yaxis->SetTitleSize(0.3/nY);
 246     yaxis->SetLabelSize(0.08);
 247     yaxis->SetTitleOffset(2.5/nY);
 248     yaxis->SetNdivisions(5);
 249     yaxis->SetTitleFont(42);
 250     yaxis->SetLabelFont(42);
 251     yaxis->SetDecimals();
 252
 253     TAxis* xaxis = stack->GetHistogram()->GetXaxis();
 254     xaxis->SetTitleSize(0.3/nY);
 255     xaxis->SetLabelSize(0.08);
 256     xaxis->SetTitleOffset(2./nY);
 257     xaxis->SetNdivisions(10);
 258     xaxis->SetTitleFont(42);
 259     xaxis->SetLabelFont(42);
 260     xaxis->SetDecimals();
 261
 262     // Redraw
 263     stack->Draw("nostack");
 264     i++;
 265     if (i >= 5) b = 2;
 266     p->cd();
 267   }
 268   canvas->SaveAs("fit_results.png");
 269   canvas->Print(pdfName, "Title:Fit summary");
 270 }
 271
 272 //____________________________________________________________________
 273 /**
 274  * Draw ring fits
 275  *
 276  * @param fname
 277  *
 278  * @ingroup pwg2_forward_scripts_corr
 279  */
 280 void DrawRings(const char* fname="AnalysisResults.root")
 281 {
 282   if (!CheckFitter(fname)) {
 283     Error("DrawFits", "File not opened");
 284     return;
 285   }
 286   if (!CheckCanvas()) {
 287     Error("DrawFits", "No canvas");
 288     return;
 289   }
 290   canvas->Clear();
 291
 292   canvas->Clear();
 293   canvas->Divide(1, 5,0,0);
 294
 295   const char* dets[] = { "FMD1I", "FMD2I", "FMD2O", "FMD3I", "FMD3O", 0 };
 296   for (Int_t i = 0; i < 5; i++) {
 297     TVirtualPad* p = canvas->cd(i+1);
 298     p->SetGridx();
 299     p->SetFillColor(0);
 300     p->SetFillStyle(0);
 301     p->SetLogy();
 302     TList* d = static_cast<TList*>(fitter->FindObject(dets[i]));
 303     if (!d) {
 304       Warning("DrawFits", "List %s not found", dets[i]);
 305       continue;
 306     }
 307     TH1* edist = static_cast<TH1*>(d->FindObject(Form("%s_edist", dets[i])));
 308     if (!edist) {
 309       Warning("DrawFits", "Histogram %s_edist not found", dets[i]);
 310       continue;
 311     }
 312     edist->Draw();
 313     TF1*   f = 0;
 314     TIter  nextF(edist->GetListOfFunctions());
 315     while ((f = static_cast<TF1*>(nextF()))) {
 316       Double_t chi2 = f->GetChisquare();
 317       Int_t    ndf  = f->GetNDF();
 318       Printf("%s %s:\n  Range: %f-%f\n"
 319              "chi^2/ndf= %f / %d = %f",
 320              edist->GetName(), f->GetName(),
 321              f->GetXmin(), f->GetXmax(), chi2, ndf,
 322              (ndf > 0) ? chi2/ndf : 0);
 323       for (Int_t j = 0; j < f->GetNpar(); j++) {
 324         Printf("  %-20s : %9.4f +/- %9.4f",
 325                f->GetParName(j), f->GetParameter(j), f->GetParError(j));
 326       }
 327     }
 328     p->cd();
 329   }
 330   canvas->cd();
 331   canvas->Print(pdfName, "Title:Fit to rings");
 332 }
 333
 334 //____________________________________________________________________
 335 /**
 336  * Draw fits in eta bins
 337  *
 338  * @param fname
 339  *
 340  * @ingroup pwg2_forward_scripts_corr
 341  */
 342 void DrawEtaBins(const char* fname="AnalysisResults.root")
 343 {
 344   if (!CheckFitter(fname)) {
 345     Error("DrawFits", "File not opened");
 346     return;
 347   }
 348   if (!CheckCanvas()) {
 349     Error("DrawFits", "No canvas");
 350     return;
 351   }
 352   canvas->Clear();
 353   canvas->Divide(2,2,0,0);
 354
 355   for (UShort_t d=1; d<=3; d++) {
 356     UShort_t nQ = (d == 1 ? 1 : 2);
 357     for (UShort_t q = 0; q < nQ; q++) {
 358       Char_t r = (q == 0 ? 'I' : 'O');
 359
 360       TList* ring =
 361         static_cast<TList*>(fitter->FindObject(Form("FMD%d%c",d,r)));
 362       if (!ring) {
 363         Error("PrintFits", "Couldn't get ring FMD%d%c", d,r);
 364         continue;
 365       }
 366       TList* edists = static_cast<TList*>(ring->FindObject("EDists"));
 367       if (!edists) {
 368         Error("PrintFits", "Couldn't get EDists list for FMD%d%c", d,r);
 369         continue;
 370       }
 371
 372       Info("DrawEtaBins", "Drawing for FMD%d%c", d, r);
 373       TIter next(edists);
 374       TH1*  dist = 0;
 375       Int_t i    = 0;
 376       Int_t j    = 1;
 377       while ((dist = static_cast<TH1*>(next()))) {
 378         Info("DrawEtaBins", "FMD%d%c: %s", d, r, dist->GetName());
 379         if (i == 4) {
 380           i = 0;
 381           j++;
 382           canvas->Print(pdfName, Form("Title:FMD%d%c page %2d", d,r,j));
 383         }
 384         TVirtualPad* p = canvas->cd(++i);
 385         p->SetFillColor(kWhite);
 386         p->SetFillStyle(0);
 387         p->SetBorderSize(0);
 388         p->SetLogy();
 389         dist->SetMaximum(15);
 390         dist->Draw();
 391
 392       }
 393       if (i != 0) {
 394         i++;
 395         for (; i <= 4; i++) {
 396           TVirtualPad* p = canvas->cd(i);
 397           p->Clear();
 398           p->SetFillColor(kMagenta-3);
 399           p->SetFillStyle(0);
 400           p->SetBorderSize(0);
 401         }
 402         canvas->Print(pdfName, Form("FMD%d%c page %2d", d,r,j++));
 403       }
 404     }
 405   }
 406 }
 407
 408 //____________________________________________________________________
 409 /**
 410  * Draw energy loss fits to a multi-page PDF
 411  *
 412  * The input file is the result of running AliFMDELossFitter -
 413  * either directly via AliFMDELossFitterTask or as part of a larger
 414  * train (AliForwardMultiplicityTask or AliForwardMCMultiplicityTask).
 415  *
 416  * @verbatim
 417  *  file
 418  *   +- ForwardResults
 419  *       +- fmdEnergyFitter
 420  *           +- chi2   (THStack)
 421  *           +- c      (THStack)
 422  *           +- delta  (THStack)
 423  *           +- xi     (THStack)
 424  *           +- sigma  (THStack)
 425  *           +- sigman (THStack)
 426  *           +- n      (THStack)
 427  *           +- a2     (THStack)
 428  *           +- ...    (THStack)
 429  *           +- an     (THStack)
 430  *           +- FMD1I (TList)
 431  *           |   +- FMD1I_edist (TH1)
 432  *           |   +- EDists      (TList)
 433  *           ...
 434  * @endverbatim
 435  *
 436  * @param fname
 437  *
 438  * @ingroup pwg2_forward_scripts_corr
 439  */
 440 void
 441 DrawAnaELoss(const char* fname="forward_eloss.root", bool onlySummary=true)
 442 {
 443   if (!CheckCanvas()) {
 444     Error("DrawFits", "No canvas");
 445     return;
 446   }
 447   if (!onlySummary) canvas->Print(Form("%s[", pdfName));
 448   DrawSummary(fname);
 449   if (onlySummary) return;
 450   DrawRings(fname);
 451   DrawEtaBins(fname);
 452   canvas->Print(Form("%s]", pdfName));
 453 }
 454 //
 455 // EOF
 456 //