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