PWGLF/FORWARD/analysis2/corrs/DrawAnaELoss.C

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