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