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