PWGLF/FORWARD/analysis2/qa/QATrender.C

   1 /**
   2  * @file   QATrender.C
   3  * @author Christian Holm Christensen <cholm@nbi.dk>
   4  * @date   Thu Nov 17 12:28:17 2011
   5  *
   6  * @brief  Class to make the QA trending tree
   7  *
   8  * @ingroup pwglf_forward_qa_scripts
   9  */
  10 #ifndef __CINT__
  11 # include <TFile.h>
  12 # include <TTree.h>
  13 # include <TH1.h>
  14 # include <TH2.h>
  15 # include <TList.h>
  16 # include <TStyle.h>
  17 # include <TCanvas.h>
  18 # include <TPad.h>
  19 # include <TLatex.h>
  20 # include <TMath.h>
  21 # include <THStack.h>
  22 # include <TLegend.h>
  23 # include <TLegendEntry.h>
  24 # include <TLine.h>
  25 # include <TLinearFitter.h>
  26 # include <TF1.h>
  27 # include <TSystem.h>
  28 # include <fstream>
  29 # include <TFile.h>
  30 # include "QAStructs.h"
  31 # include "QARing.h"
  32 # include "QABase.h"
  33 #else
  34 class TList;
  35 class TFile;
  36 class TH1;
  37 class Quantity;
  38 class RingQuantity;
  39 class Merge;
  40 class FitStatus;
  41 class Correlation;
  42 class QARing;
  43 class QABase;
  44 class Global;
  45 class THStack;
  46 class TLatex;
  47 class TVirtualPad;
  48 #endif
  49
  50
  51 // --- A quantity ----------------------------------------------------
  52 /**
  53  * Class to make the QA trending tree
  54  *
  55  * @ingroup pwglf_forward_qa_scripts
  56  */
  57 struct QATrender : public QABase
  58 {
  59 public:
  60   /******************************************************************/
  61   /**
  62    * A ring object
  63    *
  64    */
  65   struct Ring : public QARing
  66   {
  67     /**
  68      * Constructor
  69      *
  70      * @param d Detector
  71      * @param r Ring
  72      */
  73     Ring(UShort_t d, Char_t r) : QARing(d, r)  {}
  74     /**
  75      * Get the detector specific list
  76      *
  77      * @param parent Parent list
  78      *
  79      * @return Found list or null
  80      */
  81     TList* GetDetectorList(const TList* parent) const
  82     {
  83       TObject* o = QATrender::GetSubList(parent, Form("FMD%d%c", fD, fR));
  84       return static_cast<TList*>(o);
  85     }
  86     /**
  87      * Calculate numbers from read histograms
  88      *
  89      * @param parent Parent list
  90      * @param name   Name of histogram
  91      * @param q      Qauntity to store in
  92      *
  93      * @return true on success
  94      */
  95     Bool_t ExtractYQuantity(const TList* parent, const char* name,
  96                             Quantity* q)
  97     {
  98       q->mean  = 0;
  99       q->var   = 0;
 100       q->min   = 0;
 101       q->max   = 0;
 102
 103       TH1* h = QATrender::GetHistogram(parent, name);
 104       if (!h) {
 105         // Warning("ExtractYQuantity", "Histogram %s not found", name);
 106         return false;
 107       }
 108
 109       Double_t sum   = 0;
 110       Double_t sumw  = 0;
 111       Double_t min   = 1e12;
 112       Double_t max   = 0;
 113       for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
 114         Double_t y = h->GetBinContent(i);
 115         Double_t e = h->GetBinError(i);
 116         Double_t w = 1;
 117         if (y <= 1e-12) continue;
 118         if (e != 0) w = 1 / (e*e);
 119         sum   += w * y;
 120         sumw  += w;
 121         min  =  TMath::Min(min, y);
 122         max  =  TMath::Max(max, y);
 123         // Info("", "  %s %3d: y=%f,\t e=%f,\t w=%f", name, i, y, e, w);
 124       }
 125       q->min  = min;
 126       q->max  = max;
 127       if (sumw <= 1e-6) {
 128         // Warning("ExtractYQuantity",
 129         //         "Sum of weights for %s too small %g<1e-6",
 130         //         name, sumw);
 131         return true;
 132       }
 133       q->mean = sum / sumw;
 134       q->var  = TMath::Sqrt(1/sumw);
 135       // Info(name, " mean=%f, var=%f", q->mean, q->var);
 136       return true;
 137     }
 138     /**
 139      * Process data from the energy loss fits
 140      *
 141      * @param parent Parent list
 142      *
 143      * @return true on success
 144      */
 145     Bool_t ProcessEnergyLoss(const TList* parent)
 146     {
 147       TList* det = GetDetectorList(parent);
 148       if (!det) return false;
 149
 150       TList* res = QATrender::GetSubList(det, "FitResults");
 151       if (!res) return false;
 152
 153       ExtractYQuantity(res, Form("FMD%d%c_chi2",  fD, fR), fChi2);
 154       ExtractYQuantity(res, Form("FMD%d%c_c",     fD, fR), fC);
 155       ExtractYQuantity(res, Form("FMD%d%c_delta", fD, fR), fDelta);
 156       ExtractYQuantity(res, Form("FMD%d%c_xi",    fD, fR), fXi);
 157       ExtractYQuantity(res, Form("FMD%d%c_sigma", fD, fR), fSigma);
 158
 159       TH1* status = QATrender::GetHistogram(res, "status");
 160       if (!status) return false;
 161       fFitStatus->nLow        = UShort_t(status->GetBinContent(3));
 162       fFitStatus->nCandidates = UShort_t(status->GetBinContent(4));
 163       fFitStatus->nFitted     = UShort_t(status->GetBinContent(5));
 164
 165       return true;
 166     }
 167     /**
 168      * Process data on neighbors
 169      *
 170      * @param parent Parent list
 171      * @param p      (optional) Pad to draw
 172      *
 173      * @return true on success
 174      */
 175     Bool_t ProcessNeighbors(const TList* parent, TVirtualPad* p)
 176     {
 177       if (!p) return true;  // Case of no drawing
 178
 179       TList* det = GetDetectorList(parent);
 180       if (!det) return false;
 181
 182       TH1* before = QATrender::GetHistogram(det, "neighborsBefore");
 183       TH1* after  = QATrender::GetHistogram(det, "neighborsAfter");
 184       if (!before || !after) return false;
 185
 186       if (before->GetMaximum() > 0) p->SetLogz();
 187       p->SetFillColor(0);
 188       if (fD == 3) p->SetRightMargin(0.15);
 189
 190       before->SetTitle(Form("FMD%d%c",fD,fR));
 191       before->Draw("colz");
 192       after->Draw("same box");
 193
 194       before->GetXaxis()->SetRangeUser(-.5, 2);
 195       before->GetYaxis()->SetRangeUser(-.5, 2);
 196
 197       TLatex* ltx = new TLatex(p->GetLeftMargin()+.01,
 198                                p->GetBottomMargin()+.01,
 199                                before->GetTitle());
 200       ltx->SetNDC();
 201       ltx->SetTextSize(.07);
 202       ltx->Draw();
 203
 204       return true;
 205     }
 206     /**
 207      * Process data on single, double, and triple hits
 208      *
 209      * @param parent Parent list
 210      * @param p      (optional) Pad to draw
 211      *
 212      * @return true on success
 213      */
 214     Bool_t Process123(const TList* parent, TVirtualPad* p)
 215     {
 216       TList* det = GetDetectorList(parent);
 217       if (!det) return false;
 218
 219       fMerge->one     = 0;
 220       fMerge->two     = 0;
 221       fMerge->three   = 0;
 222
 223       TH1* one   = QATrender::GetHistogram(det, "singleEloss");
 224       TH1* two   = QATrender::GetHistogram(det, "doubleEloss");
 225       TH1* three = QATrender::GetHistogram(det, "tripleEloss");
 226       if (!one || !two || !three) return false;
 227
 228       Int_t    nOne   = one->GetEntries();
 229       Int_t    nTwo   = two->GetEntries();
 230       Int_t    nThree = three->GetEntries();
 231       Int_t    total  = nOne + nTwo + nThree;
 232       if (total > 0) {
 233         fMerge->one     = Double_t(nOne) / total;
 234         fMerge->two     = Double_t(nTwo) / total;
 235         fMerge->three   = Double_t(nThree) / total;
 236       }
 237
 238       if (!p) return true;
 239
 240       one->SetStats(0);
 241       one->SetTitle(Form("FMD%d%c", fD, fR));
 242       one->GetXaxis()->SetRangeUser(0, 8);
 243
 244       if (one->GetMaximum() > 0) p->SetLogy();
 245       p->SetFillColor(0);
 246       if ((p->GetNumber() % 3) != 1) p->SetLeftMargin(.1);
 247
 248       one->SetLineStyle(1);
 249       two->SetLineStyle(2);
 250       three->SetLineStyle(3);
 251
 252       one->Draw();
 253       two->Draw("same");
 254       three->Draw("same");
 255
 256       // Double_t ts  = 0.03;
 257       Double_t y = 1-p->GetTopMargin()-1.2*gStyle->GetTitleH(); // +ts;
 258       Double_t x = p->GetLeftMargin() + .25;
 259       TLatex* ltx = new TLatex(x, y, "Fraction of ...");
 260       ltx->SetNDC();
 261       ltx->SetTextColor(kBlue+3);
 262       ltx->SetTextAlign(13);
 263       // ltx->SetTextAlign(33);
 264       ltx->SetTextSize(.07);
 265       ltx->Draw();
 266
 267       // y -= .1;
 268       DrawText(ltx, x, y, "Singles:", Form("%5.2f%%", fMerge->one * 100));
 269       DrawText(ltx, x, y, "Doubles:", Form("%5.2f%%", fMerge->two * 100));
 270       DrawText(ltx, x, y, "Triples:", Form("%5.2f%%", fMerge->three * 100));
 271
 272       p->cd();
 273
 274       return true;
 275     }
 276     /**
 277      * Process energy loss distributions
 278      *
 279      * @param p1     First parent list (sharing filter)
 280      * @param p2     Second parent list (density calculator)
 281      * @param p      (optional) Pad to draw
 282      *
 283      * @return true on success
 284      */
 285     Bool_t ProcessELoss(TList* p1, TList* p2, TVirtualPad* p)
 286     {
 287       TList* det1 = GetDetectorList(p1);
 288       if (!det1) return false;
 289
 290       TList* det2 = GetDetectorList(p2);
 291       if (!det2) return false;
 292
 293       TH1* before    =  QATrender::GetHistogram(det1, "esdEloss");
 294       TH1* after     =  QATrender::GetHistogram(det1, "anaEloss");
 295       TH1* presented =  QATrender::GetHistogram(det2, "eloss");
 296       TH1* used      =  QATrender::GetHistogram(det2, "elossUsed");
 297
 298       if (!before || !after || !presented || !used) return false;
 299
 300       const Double_t lowCut = 0.15;
 301
 302       Int_t low = before->GetXaxis()->FindBin(lowCut);
 303       Int_t ib  = Int_t(before->Integral(low,before->GetNbinsX()));
 304       Int_t ia  = Int_t(after->Integral(low,after->GetNbinsX()));
 305       Int_t ip  = Int_t(presented->Integral(low,presented->GetNbinsX()));
 306       Int_t iu  = Int_t(used->Integral(low,used->GetNbinsX()));
 307
 308       Double_t dba = ib > 0 ? (100.*(ia-ib))/ib : 0;
 309       Double_t dbp = ib > 0 ? (100.*(ip-ib))/ib : 0;
 310       Double_t dbu = ib > 0 ? (100.*(iu-ib))/ib : 0;
 311       Double_t dap = ia > 0 ? (100.*(ip-ia))/ia : 0;
 312       Double_t dau = ia > 0 ? (100.*(iu-ia))/ia : 0;
 313       Double_t dpu = ip > 0 ? (100.*(iu-ip))/ip : 0;
 314
 315       fDataLoss->merge   = dba;
 316       fDataLoss->density = dau;
 317       fDataLoss->full    = dbu;
 318
 319       if (!p) return true;
 320
 321       if (before->GetMaximum() > 0) p->SetLogy();
 322       p->SetFillColor(0);
 323       before->SetTitle(Form("FMD%d%c",fD,fR));
 324
 325       before->SetLineStyle(1);
 326       after->SetLineStyle(2);
 327       used->SetLineStyle(3);
 328
 329       before->Draw("");
 330       after->Draw("same");
 331       presented->Draw("same");
 332       used->Draw("same");
 333
 334       Double_t ts  = 0.03;
 335       Double_t x   = p->GetLeftMargin() + .25;
 336       Double_t y   = 1-p->GetTopMargin()-gStyle->GetTitleH()+ts;
 337       TLatex*  ltx = new TLatex(x, y, Form("FMD%d%c", fD, fR));
 338       ltx->SetNDC();
 339       ltx->SetTextAlign(13);
 340       ltx->SetTextSize(ts);
 341       ltx->SetTextColor(kBlue+3);
 342       // ltx->Draw();
 343       // ltx->SetTextSize(.05);
 344       TString inte(Form("Integral [%4.2f,#infty]", lowCut));
 345       DrawText(ltx, x, y, Form("%s before:", inte.Data()), Form("%9d", ib));
 346       DrawText(ltx, x, y, Form("%s after:",  inte.Data()), Form("%9d", ia));
 347       DrawText(ltx, x, y, Form("%s input:",  inte.Data()), Form("%9d", ip));
 348       DrawText(ltx, x, y, Form("%s user:",   inte.Data()), Form("%9d", iu));
 349       TLine* l = new TLine;
 350       l->SetLineWidth(1);
 351       l->DrawLineNDC(x, y-0.9*ts, 1-p->GetRightMargin()-0.01, y-0.9*ts);
 352       if (ib != 0 && ia != 0) {
 353         DrawText(ltx, x, y, "Change (merging)", Form("%5.1f%%", dba));
 354         DrawText(ltx, x, y, "Change (input)",   Form("%5.1f%% (%5.1f%%)",
 355                                                      dap, dbp));
 356         DrawText(ltx, x, y, "Change (use)",     Form("%5.1f%% (%5.1f%%)",
 357                                                      dpu, dbu));
 358       }
 359       before->GetXaxis()->SetRangeUser(0, 4);
 360       p->cd();
 361
 362       return true;
 363     }
 364     /**
 365      * Process occupancy information
 366      *
 367      * @param parent Parent list
 368      * @param p      (optional) Pad to draw
 369      *
 370      * @return true on success
 371      */
 372     Bool_t ProcessOccupancy(const TList* parent, TVirtualPad* p)
 373     {
 374       TList* det = GetDetectorList(parent);
 375       if (!det) return false;
 376
 377       TH1* occ   = QATrender::GetHistogram(det, "occupancy");
 378       if (!occ) return false;
 379
 380       fOccupancy->mean = occ->GetMean();
 381       fOccupancy->var  = occ->GetRMS();
 382       fOccupancy->min  = 1;
 383       fOccupancy->max  = 0;
 384
 385       for (Int_t i = occ->GetNbinsX(); i >= 1; i--) {
 386         Double_t y = occ->GetBinContent(i);
 387         if (y < 1e-6) continue;
 388         fOccupancy->max = occ->GetXaxis()->GetBinUpEdge(i);
 389         break;
 390       }
 391
 392       if (!p) return true;
 393       p->SetGridy();
 394       p->SetGridx();
 395       if (occ->GetMaximum() > 0) p->SetLogy();
 396       p->SetFillColor(0);
 397       p->SetRightMargin(0.01);
 398       occ->Rebin(4);
 399       occ->SetTitle(Form("FMD%d%c", fD, fR));
 400       occ->Draw("hist");
 401
 402       // Double_t ts  = 0.03;
 403       Double_t y = 1-p->GetTopMargin()-1.2*gStyle->GetTitleH(); // +ts;
 404       Double_t x = p->GetLeftMargin() + .25;
 405       TLatex* ltx = new TLatex(x, y, "");
 406       ltx->SetNDC();
 407       ltx->SetTextColor(kBlue+3);
 408       ltx->SetTextAlign(13);
 409       // ltx->SetTextAlign(33);
 410       ltx->SetTextSize(.07);
 411       ltx->Draw();
 412
 413       // y -= .1;
 414       DrawText(ltx, x, y, "Mean:",     Form("%5.2f%%", fOccupancy->mean));
 415       DrawText(ltx, x, y, "Variance:", Form("%5.2f%%", fOccupancy->var));
 416       DrawText(ltx, x, y, "Max:",      Form("%5.2f%%", fOccupancy->max));
 417
 418       p->cd();
 419
 420       return true;
 421     }
 422     /**
 423      * Process method correlations
 424      *
 425      * @param parent Parent list
 426      * @param p      (optional) Pad to draw
 427      *
 428      * @return true on success
 429      */
 430     Bool_t ProcessCorrelation(const TList* parent, TVirtualPad* p)
 431     {
 432       TList* det = GetDetectorList(parent);
 433       if (!det) return false;
 434
 435       TH1* co   = QATrender::GetHistogram(det, "elossVsPoisson");
 436       if (!co) return false;
 437       TH2* corr = static_cast<TH2*>(co);
 438
 439       fCorrelation->alpha = 0;
 440       fCorrelation->beta  = 0;
 441       fCorrelation->a     = 0;
 442       fCorrelation->ea    = 0;
 443       fCorrelation->b     = 0;
 444       fCorrelation->eb    = 0;
 445       fCorrelation->chi2  = 0;
 446
 447       Double_t xmin = -1;
 448       Double_t xmax = corr->GetXaxis()->GetXmax();
 449
 450       if (corr->GetEntries() > 0) {
 451         // Calculate the linear regression
 452         // Double_t dx    = (xmax-xmin);
 453         // Double_t rxy   = corr->GetCorrelationFactor();
 454         Double_t sx    = corr->GetRMS(1);
 455         Double_t sy    = corr->GetRMS(2);
 456         Double_t sx2   = sx*sx;
 457         Double_t sy2   = sy*sy;
 458         Double_t cxy   = corr->GetCovariance();
 459         Double_t mx    = corr->GetMean(1);
 460         Double_t my    = corr->GetMean(2);
 461         Double_t delta = 1;
 462         if (TMath::Abs(cxy) > 1e-6) {
 463           fCorrelation->beta  = ((sy2 - delta*sx2 +
 464                                   TMath::Sqrt(TMath::Power(sy2-delta*sx2,2) +
 465                                               4*delta*cxy*cxy)) / 2 / cxy);
 466           fCorrelation->alpha = my - fCorrelation->beta * mx;
 467         }
 468
 469
 470 #if 0
 471         TLinearFitter* fitter = new TLinearFitter(1);
 472         fitter->SetFormula("1 ++ x");
 473         for (Int_t i = 1; i <= corr->GetNbinsX(); i++) {
 474           Double_t x = corr->GetXaxis()->GetBinCenter(i);
 475           if (x < -1 || x > xmax) continue;
 476           for (Int_t j = 1; j <= corr->GetNbinsY(); j++) {
 477             Double_t y = corr->GetYaxis()->GetBinCenter(j);
 478             if (y < -1 || y > xmax) continue;
 479             Double_t c = corr->GetBinContent(i,j);
 480             if (c < .1) continue;
 481             fitter->AddPoint(&x, y, c);
 482           }
 483         }
 484         fitter->Eval();
 485         fCorrelation->a    = fitter->GetParameter(0);
 486         fCorrelation->ea   = fitter->GetParError(0);
 487         fCorrelation->b    = fitter->GetParameter(1);
 488         fCorrelation->eb   = fitter->GetParError(1);
 489         Double_t chi2 = fitter->GetChisquare();
 490         Int_t    ndf  = (fitter->GetNpoints() -
 491                          fitter->GetNumberFreeParameters() );
 492         fCorrelation->chi2 = chi2 / ndf;
 493 #endif
 494       }
 495       if (!p) return true;
 496
 497       p->SetGridy();
 498       p->SetGridx();
 499       if (corr->GetMaximum() > 0) p->SetLogz();
 500       // if (fD == 3) p->SetRightMargin(0.15);
 501       p->SetFillColor(0);
 502       corr->GetXaxis()->SetRangeUser(-1,xmax);
 503       corr->GetYaxis()->SetRangeUser(-1,xmax);
 504       corr->SetTitle(Form("FMD%d%c",fD,fR));
 505       corr->Draw("colz");
 506
 507       TF1* f = new TF1("f", "[0]+[1]*x", xmin, xmax);
 508       f->SetParameters(fCorrelation->alpha, fCorrelation->beta);
 509       f->SetLineWidth(1);
 510       f->SetLineStyle(1);
 511       f->Draw("same");
 512
 513       TLine* l = new TLine(-1,-1,xmax,xmax);
 514       l->SetLineWidth(1);
 515       l->SetLineStyle(2);
 516       l->SetLineColor(kBlack);
 517       l->Draw();
 518
 519       Double_t x   = p->GetLeftMargin() + .05;
 520       Double_t y   = 1-p->GetTopMargin()-gStyle->GetTitleH(); // +ts;
 521       TLatex* ltx = new TLatex(x, y, "Deming regression: y=#alpha+#beta x");
 522       ltx->SetNDC();
 523       ltx->SetTextAlign(13);
 524       ltx->SetTextSize(0.06);
 525       ltx->SetTextColor(kBlue+3);
 526       ltx->Draw();
 527
 528       DrawText(ltx, x, y, "#alpha:", Form("%5.3f", fCorrelation->alpha));
 529       DrawText(ltx, x, y, "#beta:",  Form("%5.3f", fCorrelation->beta));
 530
 531 #if 0
 532       DrawText(ltx, x, y, "A:",      Form("%5.3f#pm%5.3f",
 533                                           fCorrelation->a,
 534                                           fCorrelation->ea));
 535       DrawText(ltx, x, y, "B:",      Form("%5.3f#pm%5.3f",
 536                                           fCorrelation->b,
 537                                           fCorrelation->eb));
 538       DrawText(ltx, x, y, "#chi^{2}/#nu:", Form("%5.3f", fCorrelation->chi2));
 539 #endif
 540       p->cd();
 541
 542       return true;
 543     }
 544   };
 545
 546   /******************************************************************/
 547   /**
 548    * CTOR
 549    *
 550    * @param keep   Whehter to keep all info
 551    * @param single Single run input
 552    */
 553   QATrender(Bool_t keep=false, Bool_t single=false)
 554     : QABase(single),
 555       fCurrentFile(0),
 556       fSharingFilter(0),
 557       fEventInspector(0),
 558       fDensityCalculator(0),
 559       fEnergyFitter(0),
 560       fFiles(0),
 561       fKeep(keep)
 562   {
 563     fFMD1i = new Ring(1, 'I');
 564     fFMD2i = new Ring(2, 'I');
 565     fFMD2o = new Ring(2, 'O');
 566     fFMD3i = new Ring(3, 'I');
 567     fFMD3o = new Ring(3, 'O');
 568   }
 569   /**
 570    * DTOR
 571    */
 572   virtual ~QATrender() {}
 573   /**
 574    * Copy CTOR
 575    *
 576    * @param o Object to copy from
 577    */
 578   QATrender(const QATrender& o)
 579     : QABase(o),
 580       fCurrentFile(o.fCurrentFile),
 581       fSharingFilter(o.fSharingFilter),
 582       fEventInspector(o.fEventInspector),
 583       fDensityCalculator(o.fDensityCalculator),
 584       fEnergyFitter(o.fEnergyFitter),
 585       fFiles(0),
 586       fKeep(o.fKeep)
 587   {}
 588   /**
 589    * Assignment operator
 590    *
 591    * @return Reference to this
 592    */
 593   QATrender operator=(const QATrender&) { return *this; }
 594
 595
 596   // --- Interface ---------------------------------------------------
 597   /**
 598    * Add a file to be processed
 599    *
 600    * @param filename Name of file
 601    */
 602   void AddFile(const char* filename)
 603   {
 604     fFiles.Add(new TObjString(filename));
 605   }
 606   /**
 607    * Run the job
 608    *
 609    * @return true on success
 610    */
 611   Bool_t Run()
 612   {
 613     Init(false);
 614     TIter next(&fFiles);
 615     TObject* o = 0;
 616     Bool_t ret = true;
 617     while ((o = next())) {
 618       if (!ProcessOne(o->GetName())) ret = false;
 619     }
 620     Finish();
 621     return ret;
 622   }
 623   /**
 624    * Finish the job
 625    *
 626    */
 627   void Finish()
 628   {
 629     if (!fOutput) return;
 630     fOutput->Write();
 631     fOutput->Close();
 632     fOutput = 0;
 633     fTree   = 0;
 634     gSystem->Exec(Form("chmod g+rw %s", OutputName()));
 635   }
 636   /**
 637    * Process a single file
 638    *
 639    * @param filename File to open.
 640    *
 641    * @return true on success
 642    */
 643   Bool_t ProcessOne(const char* filename)
 644   {
 645     if (fCurrentFile) {
 646       fCurrentFile->Close();
 647       fCurrentFile = 0;
 648     }
 649
 650     fCurrentFile = TFile::Open(filename, "READ");
 651     if (!fCurrentFile) {
 652       Error("ProcessOne", "Failed to open %s", filename);
 653       return false;
 654     }
 655
 656     if (!GetLists()) {
 657       // Error("ProcessOne", "Failed to get lists from %s", filename);
 658       return false;
 659     }
 660
 661     if (!ProcessGlobal()) {
 662       Error("ProcessOne", "Failed to get global stuff from %s", filename);
 663       return false;
 664     }
 665     fTeXName = Form("qa_%09d", fGlobal->runNo);
 666     MakeCanvas(Form("QA plots for run %d", fGlobal->runNo));
 667     Bool_t eloss = ProcessELossFitter();
 668     Bool_t merge = ProcessSharingFilter();
 669     Bool_t dense = ProcessDensityCalculator();
 670     if (fTree) fTree->Fill();
 671     CloseCurrent();
 672
 673     return eloss && merge && dense;
 674   }
 675   // --- Processing member functions ---------------------------------
 676   /**
 677    * Get global stuff
 678    *
 679    * @return true on success
 680    */
 681   Bool_t ProcessGlobal()
 682   {
 683     TObject* oRun = GetObject(fEventInspector, "runNo");
 684     if (!oRun) return false;
 685
 686     fGlobal->runNo = oRun->GetUniqueID();
 687     TH1* oAcc = GetHistogram(fEventInspector,"nEventsAccepted");
 688     if (!oAcc) return false;
 689
 690     fGlobal->nAccepted = oAcc->GetEntries();
 691     fGlobal->meanVz    = oAcc->GetMean();
 692     fGlobal->sigmaVz   = oAcc->GetRMS();
 693
 694     return true;
 695   }
 696   /**
 697    * Clean a stack for histograms that match @a what
 698    *
 699    * @param stack Stack to clean
 700    * @param what  Pattern
 701    */
 702   void CleanStack(THStack* stack, const char* what)
 703   {
 704     TList*   l = stack->GetHists();
 705
 706     // Clean up list of histogram.  Histograms with no entries or
 707     // no functions are deleted.  We have to do this using the TObjLink
 708     // objects stored in the list since ROOT cannot guaranty the validity
 709     // of iterators when removing from a list - tsck.  Should just implement
 710     // TIter::Remove().
 711     TObjLink* lnk = l->FirstLink();
 712     while (lnk) {
 713       TObject* o = lnk->GetObject();
 714       TString  s(o->GetName());
 715       if (s.Contains(what)) {
 716         TObjLink* keep = lnk->Next();
 717         l->Remove(lnk);
 718         lnk = keep;
 719         continue;
 720       }
 721       lnk = lnk->Next();
 722     }
 723   }
 724   /**
 725    * Process the information from the energy loss fitter
 726    *
 727    *
 728    * @return true on succes
 729    */
 730   Bool_t ProcessELossFitter()
 731   {
 732     MakeCanvasTitle("Summary of energy loss fits");
 733
 734     THStack* chi2  = static_cast<THStack*>(GetObject(fEnergyFitter, "chi2"));
 735     THStack* c     = static_cast<THStack*>(GetObject(fEnergyFitter, "c"));
 736     THStack* delta = static_cast<THStack*>(GetObject(fEnergyFitter, "delta"));
 737     THStack* xi    = static_cast<THStack*>(GetObject(fEnergyFitter, "xi"));
 738     THStack* sigma = static_cast<THStack*>(GetObject(fEnergyFitter, "sigma"));
 739
 740     if (!chi2)  return false;
 741     if (!c)     return false;
 742     if (!delta) return false;
 743     if (!xi)    return false;
 744     if (!sigma) return false;
 745
 746     CleanStack(chi2, "_n");
 747     CleanStack(c,    "status");
 748     if (fCanvas) {
 749       const int nL = 3;
 750       THStack* stacks[] = { chi2, c, delta, xi, sigma, 0 };
 751       for (int i = 0; i < 5; i++) {
 752         THStack*     stack = stacks[i];
 753         TVirtualPad* p     = GetPad(i+1);
 754         // stack->GetHists()->ls();
 755
 756         p->SetLeftMargin(.6/nL);
 757         p->SetTopMargin(.01);
 758         p->SetRightMargin(.01);
 759         p->SetFillColor(0);
 760         p->SetFillStyle(0);
 761         p->SetGridx();
 762         stack->Draw("nostack");
 763         stack->GetHistogram()->SetYTitle(stack->GetTitle());
 764         stack->GetHistogram()->SetXTitle("#eta");
 765
 766         TAxis* yaxis = stack->GetHistogram()->GetYaxis();
 767         if (i == 0) yaxis->SetRangeUser(0,20); // Chi2
 768         if (i == 1) stack->SetMaximum(1);      // c
 769         if (i == 2) stack->SetMaximum(1);      // delta
 770         if (i == 3) stack->SetMaximum(0.1);   // xi
 771         if (i == 4) stack->SetMaximum(0.5);    // sigma{,n}
 772         // if (i == 0) p->SetLogy();
 773         yaxis->SetTitleSize(0.3/nL);
 774         yaxis->SetLabelSize(0.08);
 775         yaxis->SetTitleOffset(3/nL);
 776         yaxis->SetNdivisions(5);
 777         yaxis->SetTitleFont(42);
 778         yaxis->SetLabelFont(42);
 779         yaxis->SetDecimals();
 780
 781         TAxis* xaxis = stack->GetHistogram()->GetXaxis();
 782         xaxis->SetTitleSize(0.3/nL);
 783         xaxis->SetLabelSize(0.08);
 784         xaxis->SetTitleOffset(2./nL);
 785         xaxis->SetNdivisions(10);
 786         xaxis->SetTitleFont(42);
 787         xaxis->SetLabelFont(42);
 788         xaxis->SetDecimals();
 789
 790         stack->Draw("nostack");
 791         p->cd();
 792       }
 793       TVirtualPad* p     = GetPad(6);
 794       p->SetFillColor(kWhite);
 795       Double_t x = .3;
 796       Double_t y = .8;
 797       TLatex* l = new TLatex(x, y, "Fits to #Delta (energy loss) spectra");
 798       l->SetTextColor(kBlue+3);
 799       l->SetNDC();
 800       l->Draw();
 801       x = .05;
 802       y -= 2 * 1.2*l->GetTextSize();
 803       l->DrawLatex(x, y, "F(#Delta;c,#Delta_{p},#xi,#sigma)="
 804                    "#frac{c}{#sqrt{2#pi}#sigma}#int_{-#infty}^{#infty}d#Delta'"
 805                    "L(#Delta;#Delta',#xi) G(#Delta_{p};#Delta',#sigma^{2})");
 806       y -= 1.2*l->GetTextSize();
 807       x += .1;
 808       DrawText(l, x, y, "#chi^{2}/#nu", "Goodness of fit",         .2);
 809       DrawText(l, x, y, "c",            "Overall constant",        .2);
 810       DrawText(l, x, y, "#Delta_{p}",   "Most probable value",     .2);
 811       DrawText(l, x, y, "#xi",          "'Width' of Landau (L)",   .2);
 812       DrawText(l, x, y, "#sigma",       "'Width' of Gaussian (G)", .2);
 813
 814       // stack->GetHists()->ls();
 815
 816       PrintCanvas("fitResults", fGlobal->runNo);
 817     }
 818
 819     static_cast<Ring*>(fFMD1i)->ProcessEnergyLoss(fEnergyFitter);
 820     static_cast<Ring*>(fFMD2i)->ProcessEnergyLoss(fEnergyFitter);
 821     static_cast<Ring*>(fFMD2o)->ProcessEnergyLoss(fEnergyFitter);
 822     static_cast<Ring*>(fFMD3i)->ProcessEnergyLoss(fEnergyFitter);
 823     static_cast<Ring*>(fFMD3o)->ProcessEnergyLoss(fEnergyFitter);
 824
 825     return true;
 826   }
 827   /**
 828    * Process the information from the sharing filter
 829    *
 830    *
 831    * @return true on success
 832    */
 833   Bool_t ProcessSharingFilter()
 834   {
 835     // --- Neighbors -------------------------------------------------
 836     MakeCanvasTitle("Correlation of neighboring strips");
 837
 838     for (Int_t i = 1; i <= 6; i++) {
 839       TVirtualPad* p = GetPad(i);
 840
 841       Ring* r = GetRing(i);
 842       if (!r) continue;
 843
 844       r->ProcessNeighbors(fSharingFilter, p);
 845     }
 846     TVirtualPad* p = 0;
 847     if ((p = GetPad(4))) {
 848       p->SetFillColor(kWhite);
 849
 850       TLatex* l = new TLatex(.2, .7, "Gradient: before merging");
 851       l->SetNDC();
 852       l->SetTextColor(kBlue+3);
 853       l->Draw();
 854       l->DrawText(.2, .6, "Boxes: after merging");
 855
 856       fCanvas->cd();
 857       PrintCanvas("neighbors", fGlobal->runNo);
 858     }
 859
 860     // --- 123 -------------------------------------------------------
 861     MakeCanvasTitle("#Delta for singles, doubles, and triples");
 862
 863     for (Int_t i = 1; i <= 6; i++) {
 864       p = GetPad(i);
 865
 866       Ring* r = GetRing(i);
 867       if (!r) continue;
 868
 869       r->Process123(fSharingFilter, p);
 870     }
 871     if ((p = GetPad(4))) {
 872       TLegend* ll = new TLegend(.2, .2, .8, .8);
 873       ll->SetFillColor(0);
 874       ll->SetBorderSize(0);
 875       TLegendEntry* e = ll->AddEntry("dummy", "Singles", "l");
 876       e->SetLineStyle(1);
 877       e = ll->AddEntry("dummy", "Doubles", "l");
 878       e->SetLineStyle(2);
 879       e = ll->AddEntry("dummy", "Triples", "l");
 880       e->SetLineStyle(3);
 881       ll->Draw();
 882
 883       PrintCanvas("123", fGlobal->runNo);
 884     }
 885     return true;
 886   }
 887   /**
 888    * Process the information from the density calculator
 889    *
 890    *
 891    * @return true on success
 892    */
 893   Bool_t ProcessDensityCalculator()
 894   {
 895     // --- ELoss -----------------------------------------------------
 896     MakeCanvasTitle("Energy loss from ESD, after merging, used");
 897
 898     for (Int_t i = 1; i <= 6; i++) {
 899       TVirtualPad* p = GetPad(i);
 900
 901       Ring* r = GetRing(i);
 902       if (!r) continue;
 903
 904       r->ProcessELoss(fSharingFilter, fDensityCalculator, p);
 905     }
 906     TVirtualPad* p = 0;
 907     if ((p = GetPad(4))) {
 908       TLegend* ll = new TLegend(.2, .2, .8, .8);
 909       ll->SetFillColor(0);
 910       ll->SetBorderSize(0);
 911       TLegendEntry* e = ll->AddEntry("dummy", "From ESDs", "l");
 912       e->SetLineStyle(1);
 913       e = ll->AddEntry("dummy", "After Merging", "l");
 914       e->SetLineStyle(2);
 915       e = ll->AddEntry("dummy", "Used", "l");
 916       e->SetLineStyle(3);
 917       ll->Draw();
 918
 919       PrintCanvas("recAna", fGlobal->runNo);
 920     }
 921
 922     // --- Occupancy -------------------------------------------------
 923     MakeCanvasTitle("Occupancy");
 924
 925     for (Int_t i = 1; i <= 6; i++) {
 926       p = GetPad(i);
 927
 928       Ring* r = GetRing(i);
 929       if (!r) continue;
 930
 931       r->ProcessOccupancy(fDensityCalculator, p);
 932     }
 933     if ((p = GetPad(4))) {
 934       TLatex* ltx = new TLatex(.2, .8, "Calculated assuming Poisson stat.");
 935       ltx->SetNDC();
 936       ltx->SetTextColor(kBlue+3);
 937       ltx->Draw();
 938
 939       TObject* etaL = GetObject(fDensityCalculator, "etaLumping");
 940       TObject* phiL = GetObject(fDensityCalculator, "phiLumping");
 941       if (etaL && phiL)
 942         ltx->DrawLatex(.2, .7, Form("Regions of %s strips #times %s sectors",
 943                                     etaL->GetTitle(), phiL->GetTitle()));
 944
 945       PrintCanvas("occupancy", fGlobal->runNo);
 946     }
 947
 948     // --- Correlation of methods ------------------------------------
 949     MakeCanvasTitle("Correlation of N_{ch} methods");
 950
 951     for (Int_t i = 1; i <= 6; i++) {
 952       p = GetPad(i);
 953
 954       Ring* r = GetRing(i);
 955       if (!r) continue;
 956
 957       r->ProcessCorrelation(fDensityCalculator, p);
 958     }
 959     if ((p = GetPad(4))) {
 960       TLatex* ltx = new TLatex(.2, .8, "Correlation of N_{ch} methods");
 961       ltx->SetNDC();
 962       ltx->SetTextColor(kBlue+3);
 963       ltx->Draw();
 964       ltx->DrawLatex(.24, .7, "From #DeltaE fits along x");
 965       ltx->DrawLatex(.24, .6, "From Poisson assumption along y");
 966       ltx->DrawLatex(.24, .4, "Solid line: regression");
 967       ltx->DrawLatex(.24, .3, "Dashed line: x=y to guide the eye");
 968
 969       PrintCanvas("elossVsPoisson", fGlobal->runNo);
 970     }
 971
 972     return true;
 973   }
 974
 975   // --- Utilities ---------------------------------------------------
 976   /**
 977    * Close current file if open, and reset other pointers
 978    *
 979    */
 980   void CloseCurrent()
 981   {
 982     if (fCurrentFile) {
 983       fCurrentFile->Close();
 984       fCurrentFile       = 0;
 985       fSharingFilter     = 0;
 986       fEventInspector    = 0;
 987       fDensityCalculator = 0;
 988       fEnergyFitter      = 0;
 989     }
 990     bool keep = fKeep;
 991     if (fSingle) keep = true;
 992
 993     Close(!keep);
 994   }
 995   /**
 996    * Write full job footer
 997    *
 998    */
 999   virtual void WriteLinks()
1000   {
1001     const char* lbUrl = "https://alice-logbook.cern.ch/logbook/date_online.php";
1002     const char* crUrl = "http://alimonitor.cern.ch/configuration/index.jsp";
1003     const char* rpUrl = "http://alimonitor.cern.ch/runview/?run=";
1004     QABase::WriteLinks();
1005
1006     TString resName(fTeXName);
1007     resName.ReplaceAll("qa", "QAresults");
1008     resName.Append(".root");
1009     TFile* results = TFile::Open(resName, "READ");
1010     if (results) {
1011       *fHtml << "<h3>QA results</h3>\n"
1012              << "<ul>\n"
1013              << "<li><a href='" << resName << "'>ROOT file</a></li>\n"
1014              << "</ul>"
1015              << std::endl;
1016       results->Close();
1017     }
1018
1019     *fHtml << "<h3>External links</h3>\n"
1020            << "<ul>\n"
1021            << "  <li><a target='_blank' href='" << lbUrl
1022            << "?p_cont=rund&p_run=" <<  fGlobal->runNo
1023            << "&p_tab=dqm&p_dqmagent=FMDQAshifter'>DQM plots</a>"
1024            << " (restricted)</li>\n"
1025            << "  <li><a target='_blank' href='" << crUrl
1026            << "?raw_run=" << fGlobal->runNo
1027            << "'>Condition Table</a></li>\n"
1028            << "  <li><a target='_blank' href='" << rpUrl << fGlobal->runNo
1029            << "'>Processing Details</a></li>\n"
1030            << "</ul>" << std::endl;
1031   }
1032
1033   /**
1034    * Get the ring corresponding to a pad
1035    *
1036    * @param padNo Pad number
1037    *
1038    * @return Pointer to ring
1039    */
1040   Ring* GetRing(Int_t padNo)
1041   {
1042     QARing* r = 0;
1043     switch(padNo) {
1044     case 1: r = fFMD1i; break;
1045     case 2: r = fFMD2i; break;
1046     case 3: r = fFMD3i; break;
1047     case 5: r = fFMD2o; break;
1048     case 6: r = fFMD3o; break;
1049     }
1050     if (r) return static_cast<Ring*>(r);
1051     return 0;
1052   }
1053   /**
1054    * Get the pad corresponding to a number.  Also clears the pad
1055    *
1056    * @param padNo Pad number to get
1057    *
1058    * @return Pointer to pad or null
1059    */
1060   TVirtualPad* GetPad(Int_t padNo)
1061   {
1062     if (!fCanvas) return 0;
1063     TVirtualPad* p = fCanvas->cd(padNo);
1064     if (!p) return 0;
1065     p->Clear();
1066     p->SetFillColor(kWhite);
1067     return p;
1068   }
1069
1070   /**
1071    * Make canvas title. Canvas is divided into 3x2 pads
1072    *
1073    * @param what Title on canvas
1074    */
1075   void MakeCanvasTitle(const char* what)
1076   {
1077     if (!fCanvas) return;
1078     fCanvas->cd();
1079
1080     CanvasTitle(what);
1081     fCanvas->Divide(3,2,0,0);
1082   }
1083   // --- List utilities ----------------------------------------------
1084   /**
1085    * Get a sub-list from parent list
1086    *
1087    * @param parent Parent list
1088    * @param name   Name of sub-list
1089    *
1090    * @return Pointer to the list
1091    */
1092   static TList* GetSubList(const TList* parent, const char* name)
1093   {
1094     TList* tmp = static_cast<TList*>(parent->FindObject(name));
1095     if (!tmp) {
1096       Error("GetLists", "List %s not found in %s", name,
1097             parent->GetName());
1098       return 0;
1099     }
1100     return tmp;
1101   }
1102   /**
1103    * Get the lists from the file
1104    *
1105    *
1106    * @return true on success
1107    */
1108   Bool_t GetLists()
1109   {
1110     if (!fCurrentFile) return 0;
1111
1112     const char* folder = "ForwardResults";
1113     TList* forward = static_cast<TList*>(fCurrentFile->Get(folder));
1114     if (!forward) {
1115       Error("GetLists", "List %s not found in %s", folder,
1116             fCurrentFile->GetName());
1117       return false;
1118     }
1119
1120     fSharingFilter     = GetSubList(forward, "fmdSharingFilter");
1121     fEventInspector    = GetSubList(forward, "fmdEventInspector");
1122     fDensityCalculator = GetSubList(forward, "fmdDensityCalculator");
1123     fEnergyFitter      = GetSubList(forward, "fmdEnergyFitter");
1124
1125     if (!fSharingFilter)     return false;
1126     if (!fEventInspector)    return false;
1127     if (!fDensityCalculator) return false;
1128     if (!fEnergyFitter)     return false;
1129
1130     return true;
1131   }
1132   /**
1133    * Find and object in a list
1134    *
1135    * @param list List to look in
1136    * @param name Name of object
1137    *
1138    * @return Pointer to object or null
1139    */
1140   static TObject* GetObject(const TList* list, const char* name)
1141   {
1142     if (!list) return 0;
1143     TObject* o = list->FindObject(name);
1144     if (!o) {
1145       Error("GetObject", "Failed to find object %s in %s",
1146             name, list->GetName());
1147       return 0;
1148     }
1149     return o;
1150   }
1151   /**
1152    * Get a histogram from a list
1153    *
1154    * @param list List
1155    * @param name Name of histogram
1156    *
1157    * @return Pointer to object or null
1158    */
1159   static TH1* GetHistogram(const TList* list, const char* name)
1160   {
1161     return static_cast<TH1*>(GetObject(list, name));
1162   }
1163   /**
1164    * Draw some text on a pad
1165    *
1166    * @param l   LaTeX object to use
1167    * @param x   x coordinate
1168    * @param y   y coordinate (incremented on return)
1169    * @param c1  First text
1170    * @param c2  Second text
1171    * @param dx  Distance between c1 start and c2 start
1172    */
1173   static void
1174   DrawText(TLatex* l, Double_t x, Double_t& y, const char* c1, const char* c2,
1175            Double_t dx=.4)
1176   {
1177     y -= 1.2*l->GetTextSize();
1178     l->DrawLatex(x,    y, c1);
1179     l->DrawLatex(x+dx, y, c2);
1180   }
1181
1182
1183   // --- Members -----------------------------------------------------
1184   TFile*  fCurrentFile;          // Current input file
1185   TList*  fSharingFilter;        // Sharing filter list
1186   TList*  fEventInspector;       // Event inspector list
1187   TList*  fDensityCalculator;    // Density calculator list
1188   TList*  fEnergyFitter;         // Energy fitter list
1189   TList   fFiles;                // List of files to process
1190   Bool_t  fKeep;                 // Keep PNGs
1191 };
1192 //
1193 // EOF
1194 //
1195
1196
1197