PWGLF/FORWARD/analysis2/qa/QAPlotter.C

   1 /**
   2  * @file   QAPlotter.C
   3  * @author Christian Holm Christensen <cholm@nbi.dk>
   4  * @date   Thu Nov 17 12:08:04 2011
   5  *
   6  * @brief  Class to plot QA trends
   7  *
   8  * @ingroup pwglf_forward_qa_scripts
   9  *
  10  */
  11 #ifndef __CINT__
  12 # include "QABase.h"
  13 # include <TTree.h>
  14 # include <TGraph.h>
  15 # include <TGraphErrors.h>
  16 # include <TGraphAsymmErrors.h>
  17 # include <TMultiGraph.h>
  18 # include <TH1.h>
  19 # include <TLegend.h>
  20 # include <TCanvas.h>
  21 # include <TLine.h>
  22 # include <TArrow.h>
  23 # include <TArrayI.h>
  24 # include <TMath.h>
  25 # include <TChain.h>
  26 #else
  27 class QABase;
  28 class QARing;
  29 class TTree;
  30 class TGraphAsymmErrors;
  31 class TGraph;
  32 class TGraphErrors;
  33 class TMultiGraph;
  34 class RingQuantity;
  35 class TArrayI;
  36 class TH1;
  37 #endif
  38
  39
  40 /**
  41  * Class to plot QA trends
  42  *
  43  * @ingroup pwglf_forward_qa_scripts
  44  */
  45 struct QAPlotter : public QABase
  46 {
  47   /**
  48    * Ring class
  49    */
  50   struct Ring : public QARing
  51   {
  52     /**
  53      * Constuctor
  54      *
  55      * @param d Detector
  56      * @param r Ring
  57      */
  58     Ring(UShort_t d, Char_t r)
  59       : QARing(d, r),
  60         fGChi2(0),
  61         fGC(0),
  62         fGDelta(0),
  63         fGXi(0),
  64         fGSigma(0),
  65         fGLow(0),
  66         fGSingles(0),
  67         fGLoss(0),
  68         fGBeta(0),
  69         fGOccupancy(0)
  70     {
  71       fGChi2           = new TGraphAsymmErrors;
  72       fGC              = new TGraphAsymmErrors;
  73       fGDelta          = new TGraphAsymmErrors;
  74       fGXi             = new TGraphAsymmErrors;
  75       fGSigma          = new TGraphAsymmErrors;
  76       fGLow            = new TGraph;
  77       fGSingles        = new TGraph;
  78       fGLoss           = new TGraph;
  79       fGBeta           = new TGraph;
  80       fGOccupancy      = new TGraphAsymmErrors;
  81
  82       SetAtt(fGChi2,            "chi2",     "#LT#chi^{2}/#nu#GT");
  83       SetAtt(fGC,               "c",        "Constant");
  84       SetAtt(fGDelta,           "delta",    "#LT#Delta_{p}#GT");
  85       SetAtt(fGXi,              "xi",       "#LT#xi#GT");
  86       SetAtt(fGSigma,           "sigma",    "#LT#sigma#GT");
  87       SetAtt(fGLow,             "low",      "# of low statistics bins");
  88       SetAtt(fGSingles,         "singles",  "Fraction of single hits");
  89       SetAtt(fGLoss,            "loss",     "Data lossed due to cuts");
  90       SetAtt(fGOccupancy,       "occupancy","#LTOccupancy#GT");
  91       SetAtt(fGBeta,            "beta",     "Correlation of methods");
  92
  93     }
  94     /**
  95      * Static member function to get the ring color
  96      *
  97      * @param d Detector number
  98      * @param r Ring identifer
  99      *
 100      * @return Color
 101      */
 102     static Color_t RingColor(UShort_t d, Char_t r)
 103     {
 104       return ((d == 1 ? kRed : (d == 2 ? kGreen : kBlue))
 105               + ((r == 'I' || r == 'i') ? 2 : -3));
 106     }
 107     /**
 108      * Set graph attributes
 109      *
 110      * @param g      Graph
 111      * @param name   Name of graph
 112      */
 113     void SetAtt(TGraph* g, const char* name, const char* /*title=""*/)
 114     {
 115       Color_t c = RingColor(fD, fR);
 116       g->SetName(Form("FMD%d%c_%s", fD, fR, name));
 117       // g->SetTitle(Form("FMD%d%c %s", fD, fR,
 118       // !title || title[0] == '\0' ? name : title));
 119       Int_t marker = 20+(fD-1) + (fR == 'I' ? 0 : 4);
 120       g->SetTitle(Form("FMD%d%c", fD, fR));
 121       g->SetLineColor(c);
 122       g->SetFillColor(c);
 123       g->SetMarkerColor(c);
 124       g->SetMarkerStyle(marker);
 125       g->SetLineWidth(2);
 126       switch (marker) {
 127       case 20: g->SetMarkerSize(1.2); break;
 128       case 21: g->SetMarkerSize(1.2); break;
 129       case 22: g->SetMarkerSize(1.3); break;
 130       case 26: g->SetMarkerSize(1.1); break;
 131       }
 132     }
 133     /**
 134      * Update a graph from a RingQuantity
 135      *
 136      * @param g      Graph to update
 137      * @param q      Quantity
 138      * @param runNo  Run number
 139      */
 140     void UpdateGraph(TGraphAsymmErrors* g, RingQuantity& q,
 141                      UInt_t runNo, UInt_t /* n */)
 142     {
 143       Double_t y  = q.mean;
 144       Double_t el = y-q.min;
 145       Double_t eh = q.max-y;
 146       if (TMath::Abs(y) < 1e-6) return;
 147       Int_t    i  = g->GetN();
 148       g->SetPoint(i, runNo, y);
 149       g->SetPointError(i, 0, 0, el, eh);
 150     }
 151     /**
 152      * Update all graphs
 153      *
 154      * @param n     Entry number (not used)
 155      * @param runNo Run number
 156      *
 157      * @return true on success
 158      */
 159     Bool_t Update(UInt_t n, UInt_t runNo)
 160     {
 161       UpdateGraph(fGChi2,      *fChi2,          runNo, n);
 162       UpdateGraph(fGC,         *fC,             runNo, n);
 163       UpdateGraph(fGDelta,     *fDelta,         runNo, n);
 164       UpdateGraph(fGXi,        *fXi,            runNo, n);
 165       UpdateGraph(fGSigma,     *fSigma,         runNo, n);
 166       UpdateGraph(fGOccupancy, *fOccupancy,     runNo, n);
 167
 168       fGLow->SetPoint(n, runNo, fFitStatus->nLow);
 169
 170       if (fMerge->one > 1e-6)
 171         fGSingles->SetPoint(fGSingles->GetN(), runNo, fMerge->one);
 172       if (fCorrelation->beta > 1e-6)
 173         fGBeta->SetPoint(fGBeta->GetN(), runNo, fCorrelation->beta);
 174       if (-fDataLoss->full > 1e-6)
 175         fGLoss->SetPoint(fGLoss->GetN(), runNo, -fDataLoss->full);
 176       return true;
 177     }
 178     TGraphAsymmErrors* fGChi2;     // Graph of ELoss reduced chi-square
 179     TGraphAsymmErrors* fGC;        // Graph of ELoss constant
 180     TGraphAsymmErrors* fGDelta;    // Graph of ELoss MPV
 181     TGraphAsymmErrors* fGXi;       // Graph of ELoss Landau width
 182     TGraphAsymmErrors* fGSigma;    // Graph of ELoss Gaus width
 183     TGraph*            fGLow;      // Graph of bins with low statistics
 184     TGraph*            fGSingles;  // Graph of fraction of singles
 185     TGraph*            fGLoss;     // Graph of 'lost' data
 186     TGraph*            fGBeta;     // Graph of Poisson vs ELoss correlation
 187     TGraphAsymmErrors* fGOccupancy;// Graph of mean occupancy
 188   };
 189   /**
 190    * Constructor
 191    */
 192   QAPlotter()
 193     : QABase(),
 194       fNAccepted(0),
 195       fVz(0)
 196   {
 197     fFMD1i = new Ring(1, 'I');
 198     fFMD2i = new Ring(2, 'I');
 199     fFMD2o = new Ring(2, 'O');
 200     fFMD3i = new Ring(3, 'I');
 201     fFMD3o = new Ring(3, 'O');
 202     fNAccepted = new TGraph;
 203     fNAccepted->SetName("nAccepted");
 204     fNAccepted->SetMarkerStyle(20);
 205     fNAccepted->SetLineWidth(2);
 206
 207     fVz = new TGraphErrors;
 208     fVz->SetName("vz");
 209     fVz->SetMarkerStyle(20);
 210     fVz->SetLineWidth(2);
 211   }
 212   /**
 213    * Add a file to be processed
 214    *
 215    * @param filename Name of file
 216    */
 217   void AddFile(const char* filename)
 218   {
 219     fFiles.Add(new TObjString(filename));
 220   }
 221   /**
 222    * Make a tree
 223    *
 224    * @param read If true, read from file
 225    *
 226    * @return True on success
 227    */
 228   Bool_t MakeTree(bool read)
 229   {
 230     if (fFiles.GetEntriesFast() <= 0) return QABase::MakeTree(read);
 231
 232     TChain* chain = new TChain("T", "T");
 233     if (!chain->AddFileInfoList(&fFiles)) return false;
 234
 235     fTree   = chain;
 236     return true;
 237   }
 238
 239   /**
 240    * Run the job
 241    *
 242    */
 243   void Run()
 244   {
 245     Init(true);
 246     if (!fTree) {
 247       Error("Run", "No input tree");
 248       return;
 249     }
 250
 251     fFirst = 0xFFFFFFFF;
 252     fLast  = 0;
 253
 254     UInt_t nEntries = fTree->GetEntries();
 255     UInt_t j = 0;
 256     fRuns.Set(nEntries);
 257     for (UInt_t i = 0; i < nEntries; i++) {
 258       fTree->GetEntry(i);
 259
 260       UInt_t run = fGlobal->runNo;
 261       UInt_t nev = fGlobal->nAccepted;
 262
 263       fFirst = TMath::Min(run, fFirst);
 264       fLast  = TMath::Max(run, fLast);
 265       fRuns[i] = run;
 266
 267       Info("Run", "Got run %d with %d accepted events", run, nev);
 268       fNAccepted->SetPoint(i, run, nev);
 269       fVz->SetPoint(i, run, fGlobal->meanVz);
 270       fVz->SetPointError(i, 0, fGlobal->sigmaVz);
 271
 272       if (nev <= 100) continue;
 273       static_cast<Ring*>(fFMD1i)->Update(j, run);
 274       static_cast<Ring*>(fFMD2i)->Update(j, run);
 275       static_cast<Ring*>(fFMD2o)->Update(j, run);
 276       static_cast<Ring*>(fFMD3i)->Update(j, run);
 277       static_cast<Ring*>(fFMD3o)->Update(j, run);
 278       j++;
 279     }
 280
 281     Plot();
 282   }
 283   /**
 284    * Plot results
 285    *
 286    */
 287   void Plot()
 288   {
 289
 290     fTeXName = Form("trend_%09d_%09d", fFirst, fLast);
 291     MakeCanvas(Form("QA trends for runs %d --- %d", fFirst, fLast));
 292
 293     CanvasTitle("# of accepted events");
 294     fNAccepted->Draw("apl");
 295     PutCanvasTitle("# of accepted events");
 296     AddRuns(fNAccepted->GetHistogram(), "# of accepted events");
 297
 298     TLine* l = new TLine(fFirst, 100, fLast, 100);
 299     l->SetLineColor(kRed+2);
 300     l->SetLineStyle(2);
 301     l->Draw();
 302     PrintCanvas("nAccepted");
 303
 304     CanvasTitle("#LTv_{z}#GT");
 305     fVz->Draw("apl");
 306     PutCanvasTitle("Mean z coordinate of interaction point");
 307     AddRuns(fVz->GetHistogram(), "#LTv_{z}#GT");
 308     PrintCanvas("vz");
 309
 310     TMultiGraph* chi2           = new TMultiGraph;
 311     TMultiGraph* c              = new TMultiGraph;
 312     TMultiGraph* delta          = new TMultiGraph;
 313     TMultiGraph* xi             = new TMultiGraph;
 314     TMultiGraph* sigma          = new TMultiGraph;
 315     TMultiGraph* low            = new TMultiGraph;
 316     TMultiGraph* singles        = new TMultiGraph;
 317     TMultiGraph* loss           = new TMultiGraph;
 318     TMultiGraph* occ    = new TMultiGraph;
 319     TMultiGraph* beta           = new TMultiGraph;
 320     chi2        ->SetName("chi2");
 321     c           ->SetName("c");
 322     delta       ->SetName("delta");
 323     xi          ->SetName("xi");
 324     sigma       ->SetName("sigma");
 325     low         ->SetName("low");
 326     singles     ->SetName("singles");
 327     loss        ->SetName("loss");
 328     beta        ->SetName("beta");
 329     occ ->SetName("occupancy");
 330
 331
 332     AddToMulti(fFMD1i,chi2, c, delta, xi, sigma, low, singles, loss, beta, occ);
 333     AddToMulti(fFMD2i,chi2, c, delta, xi, sigma, low, singles, loss, beta, occ);
 334     AddToMulti(fFMD2o,chi2, c, delta, xi, sigma, low, singles, loss, beta, occ);
 335     AddToMulti(fFMD3i,chi2, c, delta, xi, sigma, low, singles, loss, beta, occ);
 336     AddToMulti(fFMD3o,chi2, c, delta, xi, sigma, low, singles, loss, beta, occ);
 337
 338     PlotMulti(chi2,     "#LT#chi^{2}/#nu#GT from #Delta fits", true);
 339     PlotMulti(c,        "#LTc#GT from #Delta fits");
 340     PlotMulti(delta,    "#LT#Delta_{p}#GT from #Delta fits");
 341     PlotMulti(xi,       "#LT#xi#GT from #Delta fits");
 342     PlotMulti(sigma,    "#LT#sigma#GT from #Delta fits");
 343     PlotMulti(low,      "Bins with too low statistics");
 344     PlotMulti(singles,  "Fraction of single hits");
 345     PlotMulti(loss,     "% of hits 'lost' due to merging+cuts");
 346     PlotMulti(occ,      "#LTOccupancy#GT [%]", true);
 347     PlotMulti(beta,     "Correlation of methods");
 348
 349     fStore->cd();
 350     fNAccepted->Write();
 351     chi2->Write();
 352     c->Write();
 353     delta->Write();
 354     sigma->Write();
 355     low->Write();
 356     singles->Write();
 357     loss->Write();
 358     occ->Write();
 359     beta->Write();
 360
 361     Close(false); // Do not delete PNGs
 362   }
 363   /**
 364    * Add graphs from a ring to the multi graphs
 365    *
 366    * @param qr        Ring
 367    * @param chi2      ELoss reduced chi-square
 368    * @param c         ELoss constant
 369    * @param delta     ELoss MPV
 370    * @param xi        ELoss Landau width
 371    * @param sigma     ELoss Gaus width
 372    * @param low       bins with low statistics
 373    * @param singles   fraction of singles
 374    * @param loss      'lost' data
 375    * @param beta      Poisson vs ELoss correlation
 376    * @param occupancy mean occupancy
 377    */
 378   void AddToMulti(QARing* qr,
 379                   TMultiGraph* chi2,
 380                   TMultiGraph* c,
 381                   TMultiGraph* delta,
 382                   TMultiGraph* xi,
 383                   TMultiGraph* sigma,
 384                   TMultiGraph* low,
 385                   TMultiGraph* singles,
 386                   TMultiGraph* loss,
 387                   TMultiGraph* beta,
 388                   TMultiGraph* occupancy)
 389   {
 390     Ring* r = static_cast<Ring*>(qr);
 391     chi2        ->Add(r->fGChi2);
 392     c           ->Add(r->fGC);
 393     delta       ->Add(r->fGDelta);
 394     xi          ->Add(r->fGXi);
 395     sigma       ->Add(r->fGSigma);
 396     low         ->Add(r->fGLow);
 397     singles     ->Add(r->fGSingles);
 398     loss        ->Add(r->fGLoss);
 399     occupancy   ->Add(r->fGOccupancy);
 400     beta        ->Add(r->fGBeta);
 401   }
 402   /**
 403    * Plot a multi-graph
 404    *
 405    * @param mg     Multi graph
 406    * @param title  Title
 407    * @param logy   If true, make @f$\log@f$ scale
 408    */
 409   void PlotMulti(TMultiGraph* mg, const char* title, Bool_t logy=false)
 410   {
 411     CanvasTitle(title);
 412     // fCanvas->SetBottomMargin(.15);
 413     fCanvas->SetLeftMargin(.08);
 414     fCanvas->SetTopMargin(.1);
 415     fCanvas->SetLogy(logy);
 416     // fCanvas->SetRightMargin(.2);
 417     mg->Draw("apl");
 418
 419     TH1*     h   = mg->GetHistogram();
 420     Double_t max = h->GetMaximum();
 421     Double_t min = h->GetMinimum();
 422     if (h->GetMinimum() == 0) {
 423       min = min - .1*(max-min);
 424       h->SetMinimum(min);
 425     }
 426     Int_t    x1  = h->GetXaxis()->GetXmin();
 427     Int_t    x2  = h->GetXaxis()->GetXmax();
 428
 429     TLegend* l = new TLegend(.1, .91, .97, .95);
 430     l->SetNColumns(5);
 431     l->SetFillColor(0);
 432     l->SetFillStyle(0);
 433     l->SetBorderSize(0);
 434     l->SetTextFont(42);
 435     TIter next(mg->GetListOfGraphs());
 436     mg->GetListOfGraphs();
 437     TGraph* g = 0;
 438
 439     // Get the runs we have here
 440     TArrayI runs(fRuns.GetSize());
 441     runs.Reset(INT_MAX);
 442     Int_t   j = 0;
 443     while ((g = static_cast<TGraph*>(next()))) {
 444       l->AddEntry(g, g->GetTitle(), "lp");
 445       Double_t* xs = g->GetX();
 446       Int_t     n  = g->GetN();
 447       for (Int_t i = 0; i < n; i++) {
 448         if (FindRun(runs, Int_t(xs[i])) >= 0) continue;
 449         runs.SetAt(xs[i], j++);
 450       }
 451       Double_t ymean = g->GetMean(2);
 452       Double_t xh    = x2 - .03 * Double_t(x2-x1);
 453       TLine* lm = new TLine(x1, ymean, xh, ymean);
 454       lm->SetLineColor(g->GetLineColor());
 455       lm->SetLineStyle(2);
 456       lm->SetLineWidth(1);
 457       lm->Draw();
 458       TLatex* al = new TLatex(xh, ymean, g->GetTitle());
 459       al->SetTextColor(g->GetLineColor());
 460       al->SetTextFont(42);
 461       al->SetTextSize(.02);
 462       al->SetTextAlign(12);
 463       al->Draw();
 464     }
 465     l->Draw();
 466
 467     AddRuns(h, title, &runs);
 468
 469     PrintCanvas(mg->GetName());
 470   }
 471   /**
 472    * Find a run
 473    *
 474    * @param runs List of runs
 475    * @param run  Run to find
 476    *
 477    * @return Index of run in run list
 478    */
 479   Int_t FindRun(const TArrayI& runs, Int_t run)
 480   {
 481     std::sort(&(runs.fArray[0]), &(runs.fArray[runs.GetSize()]));
 482     Int_t idx = TMath::BinarySearch(runs.GetSize(), runs.fArray, run);
 483     if (idx >= runs.GetSize() || idx < 0 || runs[idx] != run) return -1;
 484     return idx;
 485   }
 486   /**
 487    * Add run labels at appropriate places on the plot
 488    *
 489    * @param h      Frame histogram
 490    * @param title  Title
 491    * @param runs   List of runs, if any
 492    */
 493   void AddRuns(TH1* h, const char* title, TArrayI* runs=0)
 494   {
 495     h->GetXaxis()->SetNoExponent();
 496     // h->GetXaxis()->SetTitleOffset(1);
 497     h->SetYTitle(title);
 498     h->SetXTitle("Run #");
 499
 500     Int_t    r1  = h->GetXaxis()->GetXmin();
 501     Int_t    r2  = h->GetXaxis()->GetXmax();
 502     Double_t lx  = 0;
 503     Double_t tx  = .045; // (r2 - r1) / 18;
 504     Double_t wx  = 1 - fCanvas->GetLeftMargin() - fCanvas->GetRightMargin();
 505     Double_t dy  = .025;
 506     Double_t y   = fCanvas->GetBottomMargin()+dy;
 507     for (Int_t i = 0; i < fRuns.GetSize(); i++) {
 508       Int_t    r = fRuns[i];
 509       Double_t x = fCanvas->GetLeftMargin() + wx*Double_t(r-r1)/(r2-r1);
 510
 511       // Skip runs out of range
 512       if (r < r1 || r > r2) continue;
 513
 514       // Skip runs not in the graphs
 515       if (runs) {
 516         if (FindRun(*runs, r) < 0) {
 517           continue;
 518         }
 519       }
 520
 521
 522       if (TMath::Abs(x - lx) < tx) y += dy;
 523       else                         y =  fCanvas->GetBottomMargin() + dy;
 524
 525
 526       // Info(h->GetName(), "%6d (x,y)=(%12f,%12f) |lx-x|=|%f-x|=%f",
 527       //      r, x, y, lx, TMath::Abs(lx-x));
 528       lx = x;
 529
 530       Double_t* xa = fNAccepted->GetX();
 531       Int_t     na = fNAccepted->GetN();
 532       Int_t idx = TMath::BinarySearch(na, xa, Double_t(r));
 533       Color_t color = kBlue+3;
 534       if (idx >= 0  && idx < na  && r == xa[idx] &&
 535           fNAccepted->GetY()[idx] < 10000)
 536         color = kRed+3;
 537
 538       TLatex* ll = new TLatex(x, y, Form("%d", r));
 539       ll->SetNDC();
 540       ll->SetTextAlign(21);
 541       ll->SetTextSize(0.02);
 542       ll->SetTextColor(color);
 543       ll->SetTextFont(42);
 544       // ll->SetTextAngle(90);
 545       ll->Draw();
 546       TLine* tl = new TLine(x, y, x, 1-fCanvas->GetTopMargin());
 547       tl->SetBit(TLine::kLineNDC);
 548       tl->SetLineStyle(3);
 549       tl->SetLineColor(color);
 550       tl->Draw();
 551     }
 552   }
 553   /**
 554    * Output list of runs
 555    *
 556    * @param o Output stream
 557    */
 558   void WriteRuns(std::ostream& o)
 559   {
 560     o << "<div class='jobid'><!--JOBID--></div>\n"
 561       << "<div class='runs'>\n"
 562       << "  Runs:<br> ";
 563     for (Int_t i = 0; i < fRuns.GetSize(); i++) {
 564       o << "<a href='qa_" << Form("%09d", fRuns[i]) << ".html'>"
 565         << fRuns[i] << "</a> " << std::flush;
 566     }
 567     o << "\n"
 568       << "</div>" << std::endl;
 569   }
 570   /**
 571    * Write page footer
 572    *
 573    */
 574   void WriteFooter()
 575   {
 576     WriteRuns(*fHtml);
 577     QABase::WriteFooter();
 578   }
 579   /**
 580    * Write out image footer
 581    *
 582    * @param o Output stream
 583    * @param pngName Name of the PNG file
 584    */
 585   void WriteImageFooter(std::ostream& o, const char* pngName)
 586   {
 587     WriteRuns(o);
 588     QABase::WriteImageFooter(o, pngName);
 589   }
 590   TGraph*       fNAccepted; // Graph of number of accepted events
 591   TGraphErrors* fVz;        // Graph of mean vertex
 592   UInt_t        fFirst;     // First run
 593   UInt_t        fLast;      // Last run
 594   TArrayI       fRuns;      // Seen runs
 595   TObjArray     fFiles;
 596 };
 597
 598 //
 599 // EOF
 600 //
 601
 602
 603