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_analysis_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_analysis_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_analysis_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_analysis_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 /**
 135  * Draw summary
 136  *
 137  * @param fname
 138  *
 139  * @ingroup pwg2_forward_analysis_scripts_corr
 140  */
 141 void DrawSummary(const char* fname="forward_eloss.root")
 142 {
 143   if (!CheckFitter(fname)) {
 144     Error("DrawFits", "File not opened");
 145     return;
 146   }
 147   if (!CheckCanvas()) {
 148     Error("DrawFits", "No canvas");
 149     return;
 150   }
 151   canvas->Clear();
 152
 153   THStack* chi2nu;
 154   THStack* c;
 155   THStack* delta;
 156   THStack* xi;
 157   THStack* sigma;
 158   THStack* sigman;
 159   THStack* n;
 160   TList stacks;
 161   stacks.Add(chi2nu = static_cast<THStack*>(fitter->FindObject("chi2")));
 162   stacks.Add(c      = static_cast<THStack*>(fitter->FindObject("c")));
 163   stacks.Add(delta  = static_cast<THStack*>(fitter->FindObject("delta")));
 164   stacks.Add(xi     = static_cast<THStack*>(fitter->FindObject("xi")));
 165   stacks.Add(sigma  = static_cast<THStack*>(fitter->FindObject("sigma")));
 166   stacks.Add(sigman = static_cast<THStack*>(fitter->FindObject("sigman")));
 167   stacks.Add(n      = static_cast<THStack*>(fitter->FindObject("n")));
 168   Int_t baseA = stacks.GetEntries()+1;
 169   Int_t i=2;
 170   while (true) {
 171     TObject* o = fitter->FindObject(Form("a%d",i++));
 172     if (!o) break;
 173     Info("DrawFits", "Adding %s", o->GetName());
 174     stacks.Add(o);
 175   }
 176   // stacks.ls();
 177   Int_t nMax = stacks.GetEntries();
 178   for (Int_t i = nMax-1; i >= baseA; i--) {
 179     THStack* stack   = static_cast<THStack*>(stacks.At(i));
 180     TIter    nextH(stack->GetHists());
 181     TH1*     hist    = 0;
 182     Bool_t   hasData = kFALSE;
 183     while ((hist = static_cast<TH1*>(nextH())))
 184       if (hist->Integral() > 0) hasData = kTRUE;
 185     if (!hasData) nMax--;
 186   }
 187
 188   canvas->SetRightMargin(0.05);
 189   canvas->SetTopMargin(0.05);
 190   canvas->Divide(2, (nMax+1)/2, 0.1, 0, 0);
 191
 192   TIter next(&stacks);
 193   THStack* stack = 0;
 194   i = 0;
 195   Int_t b = 1;
 196   while ((stack = static_cast<THStack*>(next()))) {
 197     if (i > nMax) break;
 198     TVirtualPad* p = canvas->cd(1+i/5 + 2*(i%5));
 199     p->SetLeftMargin(.15);
 200     p->SetFillColor(0);
 201     p->SetFillStyle(0);
 202     p->SetGridx();
 203     stack->Draw("nostack");
 204     stack->GetHistogram()->SetYTitle(stack->GetTitle());
 205     stack->GetHistogram()->SetXTitle("#eta");
 206
 207     TAxis* yaxis = stack->GetHistogram()->GetYaxis();
 208     if (i == 0)                     yaxis->SetRangeUser(0,20); // Chi2
 209     if (i == 1)                     stack->SetMaximum(1);      // c
 210     if (i == 2)                     stack->SetMaximum(1);      // delta
 211     if (i == 3)                     stack->SetMaximum(0.1);   // xi
 212     if (i == 4 || i == 5)           stack->SetMaximum(0.5);    // sigma{,n}
 213     if (i == 7)                     stack->SetMaximum(0.5);    // a
 214     yaxis->SetTitleSize(0.15);
 215     yaxis->SetLabelSize(0.08);
 216     yaxis->SetTitleOffset(0.35);
 217     yaxis->SetNdivisions(5);
 218
 219     TAxis* xaxis = stack->GetHistogram()->GetXaxis();
 220     xaxis->SetTitleSize(0.15);
 221     xaxis->SetLabelSize(0.08);
 222     xaxis->SetTitleOffset(0.35);
 223     xaxis->SetNdivisions(10);
 224
 225     // Redraw
 226     stack->Draw("nostack");
 227     i++;
 228     if (i >= 5) b = 2;
 229     p->cd();
 230   }
 231   canvas->Print(pdfName, "Title:Fit summary");
 232 }
 233
 234 //____________________________________________________________________
 235 /**
 236  * Draw ring fits
 237  *
 238  * @param fname
 239  *
 240  * @ingroup pwg2_forward_analysis_scripts_corr
 241  */
 242 void DrawRings(const char* fname="AnalysisResults.root")
 243 {
 244   if (!CheckFitter(fname)) {
 245     Error("DrawFits", "File not opened");
 246     return;
 247   }
 248   if (!CheckCanvas()) {
 249     Error("DrawFits", "No canvas");
 250     return;
 251   }
 252   canvas->Clear();
 253
 254   canvas->Clear();
 255   canvas->Divide(1, 5,0,0);
 256
 257   const char* dets[] = { "FMD1I", "FMD2I", "FMD2O", "FMD3I", "FMD3O", 0 };
 258   for (Int_t i = 0; i < 5; i++) {
 259     TVirtualPad* p = canvas->cd(i+1);
 260     p->SetGridx();
 261     p->SetFillColor(0);
 262     p->SetFillStyle(0);
 263     p->SetLogy();
 264     TList* d = static_cast<TList*>(fitter->FindObject(dets[i]));
 265     if (!d) {
 266       Warning("DrawFits", "List %s not found", dets[i]);
 267       continue;
 268     }
 269     TH1* edist = static_cast<TH1*>(d->FindObject(Form("%s_edist", dets[i])));
 270     if (!edist) {
 271       Warning("DrawFits", "Histogram %s_edist not found", dets[i]);
 272       continue;
 273     }
 274     edist->Draw();
 275     TF1*   f = 0;
 276     TIter  nextF(edist->GetListOfFunctions());
 277     while ((f = static_cast<TF1*>(nextF()))) {
 278       Double_t chi2 = f->GetChisquare();
 279       Int_t    ndf  = f->GetNDF();
 280       Printf("%s %s:\n  Range: %f-%f\n"
 281              "chi^2/ndf= %f / %d = %f",
 282              edist->GetName(), f->GetName(),
 283              f->GetXmin(), f->GetXmax(), chi2, ndf,
 284              (ndf > 0) ? chi2/ndf : 0);
 285       for (Int_t j = 0; j < f->GetNpar(); j++) {
 286         Printf("  %-20s : %9.4f +/- %9.4f",
 287                f->GetParName(j), f->GetParameter(j), f->GetParError(j));
 288       }
 289     }
 290     p->cd();
 291   }
 292   canvas->cd();
 293   canvas->Print(pdfName, "Title:Fit to rings");
 294 }
 295
 296 //____________________________________________________________________
 297 /**
 298  * Draw fits in eta bins
 299  *
 300  * @param fname
 301  *
 302  * @ingroup pwg2_forward_analysis_scripts_corr
 303  */
 304 void DrawEtaBins(const char* fname="AnalysisResults.root")
 305 {
 306   if (!CheckFitter(fname)) {
 307     Error("DrawFits", "File not opened");
 308     return;
 309   }
 310   if (!CheckCanvas()) {
 311     Error("DrawFits", "No canvas");
 312     return;
 313   }
 314   canvas->Clear();
 315   canvas->Divide(2,2,0,0);
 316
 317   for (UShort_t d=1; d<=3; d++) {
 318     UShort_t nQ = (d == 1 ? 1 : 2);
 319     for (UShort_t q = 0; q < nQ; q++) {
 320       Char_t r = (q == 0 ? 'I' : 'O');
 321
 322       TList* ring =
 323         static_cast<TList*>(fitter->FindObject(Form("FMD%d%c",d,r)));
 324       if (!ring) {
 325         Error("PrintFits", "Couldn't get ring FMD%d%c", d,r);
 326         continue;
 327       }
 328       TList* edists = static_cast<TList*>(ring->FindObject("EDists"));
 329       if (!edists) {
 330         Error("PrintFits", "Couldn't get EDists list for FMD%d%c", d,r);
 331         continue;
 332       }
 333       TIter next(edists);
 334       TH1*  dist = 0;
 335       Int_t i    = 0;
 336       Int_t j    = 1;
 337       while ((dist = static_cast<TH1*>(next()))) {
 338         if (i == 4) {
 339           i = 0;
 340           j++;
 341           canvas->Print(pdfName, Form("Title:FMD%d%c page %2d", d,r,j));
 342         }
 343         TVirtualPad* p = canvas->cd(++i);
 344         p->SetFillColor(kWhite);
 345         p->SetFillStyle(0);
 346         p->SetBorderSize(0);
 347         p->SetLogy();
 348         dist->SetMaximum(15);
 349         dist->Draw();
 350
 351       }
 352       if (i != 0) {
 353         i++;
 354         for (; i <= 4; i++) {
 355           TVirtualPad* p = canvas->cd(i);
 356           p->Clear();
 357           p->SetFillColor(kMagenta-3);
 358           p->SetFillStyle(0);
 359           p->SetBorderSize(0);
 360         }
 361         canvas->Print(pdfName, Form("FMD%d%c page %2d", d,r,j++));
 362       }
 363     }
 364   }
 365 }
 366
 367 //____________________________________________________________________
 368 /**
 369  * Draw energy loss fits to a multi-page PDF
 370  *
 371  * The input file is the result of running AliFMDELossFitter -
 372  * either directly via AliFMDELossFitterTask or as part of a larger
 373  * train (AliForwardMultiplicityTask or AliForwardMCMultiplicityTask).
 374  *
 375  * @verbatim
 376  *  file
 377  *   +- ForwardResults
 378  *       +- fmdEnergyFitter
 379  *           +- chi2   (THStack)
 380  *           +- c      (THStack)
 381  *           +- delta  (THStack)
 382  *           +- xi     (THStack)
 383  *           +- sigma  (THStack)
 384  *           +- sigman (THStack)
 385  *           +- n      (THStack)
 386  *           +- a2     (THStack)
 387  *           +- ...    (THStack)
 388  *           +- an     (THStack)
 389  *           +- FMD1I (TList)
 390  *           |   +- FMD1I_edist (TH1)
 391  *           |   +- EDists      (TList)
 392  *           ...
 393  * @endverbatim
 394  *
 395  * @param fname
 396  *
 397  * @ingroup pwg2_forward_analysis_scripts_corr
 398  */
 399 void
 400 DrawAnaELoss(const char* fname="forward_eloss.root")
 401 {
 402   if (!CheckCanvas()) {
 403     Error("DrawFits", "No canvas");
 404     return;
 405   }
 406   canvas->Print(Form("%s[", pdfName));
 407   DrawSummary(fname);
 408   DrawRings(fname);
 409   DrawEtaBins(fname);
 410   canvas->Print(Form("%s]", pdfName));
 411 }
 412 //
 413 // EOF
 414 //