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