PWG2/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 pwg2_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 pwg2_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   QATrender(Bool_t keep=false, Bool_t single=false)
 551     : QABase(single),
 552       fCurrentFile(0),
 553       fSharingFilter(0),
 554       fEventInspector(0),
 555       fDensityCalculator(0),
 556       fEnergyFitter(0),
 557       fFiles(0),
 558       fKeep(keep)
 559   {
 560     fFMD1i = new Ring(1, 'I');
 561     fFMD2i = new Ring(2, 'I');
 562     fFMD2o = new Ring(2, 'O');
 563     fFMD3i = new Ring(3, 'I');
 564     fFMD3o = new Ring(3, 'O');
 565   }
 566   /**
 567    * DTOR
 568    */
 569   virtual ~QATrender() {}
 570   /**
 571    * Copy CTOR
 572    *
 573    */
 574   QATrender(const QATrender& o)
 575     : QABase(o),
 576       fCurrentFile(o.fCurrentFile),
 577       fSharingFilter(o.fSharingFilter),
 578       fEventInspector(o.fEventInspector),
 579       fDensityCalculator(o.fDensityCalculator),
 580       fEnergyFitter(o.fEnergyFitter),
 581       fFiles(0),
 582       fKeep(o.fKeep)
 583   {}
 584   /**
 585    * Assignment operator
 586    *
 587    * @return Reference to this
 588    */
 589   QATrender operator=(const QATrender&) { return *this; }
 590
 591
 592   // --- Interface ---------------------------------------------------
 593   /**
 594    * Add a file to be processed
 595    *
 596    * @param filename Name of file
 597    */
 598   void AddFile(const char* filename)
 599   {
 600     fFiles.Add(new TObjString(filename));
 601   }
 602   /**
 603    * Run the job
 604    *
 605    */
 606   void Run()
 607   {
 608     Init(false);
 609     TIter next(&fFiles);
 610     TObject* o = 0;
 611     while ((o = next())) {
 612       ProcessOne(o->GetName());
 613     }
 614     Finish();
 615   }
 616   /**
 617    * Finish the job
 618    *
 619    */
 620   void Finish()
 621   {
 622     if (!fOutput) return;
 623     fOutput->Write();
 624     fOutput->Close();
 625     fOutput = 0;
 626     fTree   = 0;
 627     gSystem->Exec(Form("chmod g+rw %s", OutputName()));
 628   }
 629   /**
 630    * Process a single file
 631    *
 632    * @param filename File to open.
 633    *
 634    * @return true on success
 635    */
 636   Bool_t ProcessOne(const char* filename)
 637   {
 638     if (fCurrentFile) {
 639       fCurrentFile->Close();
 640       fCurrentFile = 0;
 641     }
 642
 643     fCurrentFile = TFile::Open(filename, "READ");
 644     if (!fCurrentFile) {
 645       Error("ProcessOne", "Failed to open %s", filename);
 646       return false;
 647     }
 648
 649     if (!GetLists()) {
 650       // Error("ProcessOne", "Failed to get lists from %s", filename);
 651       return false;
 652     }
 653
 654     if (!ProcessGlobal()) {
 655       Error("ProcessOne", "Failed to get global stuff from %s", filename);
 656       return false;
 657     }
 658     fTeXName = Form("qa_%09d", fGlobal->runNo);
 659     MakeCanvas(Form("QA plots for run %d", fGlobal->runNo));
 660     Bool_t eloss = ProcessELossFitter();
 661     Bool_t merge = ProcessSharingFilter();
 662     Bool_t dense = ProcessDensityCalculator();
 663     if (fTree) fTree->Fill();
 664     CloseCurrent();
 665
 666     return eloss && merge && dense;
 667   }
 668   // --- Processing member functions ---------------------------------
 669   /**
 670    * Get global stuff
 671    *
 672    * @return true on success
 673    */
 674   Bool_t ProcessGlobal()
 675   {
 676     TObject* oRun = GetObject(fEventInspector, "runNo");
 677     if (!oRun) return false;
 678
 679     fGlobal->runNo = oRun->GetUniqueID();
 680     TH1* oAcc = GetHistogram(fEventInspector,"nEventsAccepted");
 681     if (!oAcc) return false;
 682
 683     fGlobal->nAccepted = oAcc->GetEntries();
 684     fGlobal->meanVz    = oAcc->GetMean();
 685     fGlobal->sigmaVz   = oAcc->GetRMS();
 686
 687     return true;
 688   }
 689   /**
 690    * Clean a stack for histograms that match @a what
 691    *
 692    * @param stack Stack to clean
 693    * @param what  Pattern
 694    */
 695   void CleanStack(THStack* stack, const char* what)
 696   {
 697     TList*   l = stack->GetHists();
 698
 699     // Clean up list of histogram.  Histograms with no entries or
 700     // no functions are deleted.  We have to do this using the TObjLink
 701     // objects stored in the list since ROOT cannot guaranty the validity
 702     // of iterators when removing from a list - tsck.  Should just implement
 703     // TIter::Remove().
 704     TObjLink* lnk = l->FirstLink();
 705     while (lnk) {
 706       TObject* o = lnk->GetObject();
 707       TString  s(o->GetName());
 708       if (s.Contains(what)) {
 709         TObjLink* keep = lnk->Next();
 710         l->Remove(lnk);
 711         lnk = keep;
 712         continue;
 713       }
 714       lnk = lnk->Next();
 715     }
 716   }
 717   /**
 718    * Process the information from the energy loss fitter
 719    *
 720    *
 721    * @return true on succes
 722    */
 723   Bool_t ProcessELossFitter()
 724   {
 725     MakeCanvasTitle("Summary of energy loss fits");
 726
 727     THStack* chi2  = static_cast<THStack*>(GetObject(fEnergyFitter, "chi2"));
 728     THStack* c     = static_cast<THStack*>(GetObject(fEnergyFitter, "c"));
 729     THStack* delta = static_cast<THStack*>(GetObject(fEnergyFitter, "delta"));
 730     THStack* xi    = static_cast<THStack*>(GetObject(fEnergyFitter, "xi"));
 731     THStack* sigma = static_cast<THStack*>(GetObject(fEnergyFitter, "sigma"));
 732
 733     if (!chi2)  return false;
 734     if (!c)     return false;
 735     if (!delta) return false;
 736     if (!xi)    return false;
 737     if (!sigma) return false;
 738
 739     CleanStack(chi2, "_n");
 740     CleanStack(c,    "status");
 741     if (fCanvas) {
 742       const int nL = 3;
 743       THStack* stacks[] = { chi2, c, delta, xi, sigma, 0 };
 744       for (int i = 0; i < 5; i++) {
 745         THStack*     stack = stacks[i];
 746         TVirtualPad* p     = GetPad(i+1);
 747         // stack->GetHists()->ls();
 748
 749         p->SetLeftMargin(.6/nL);
 750         p->SetTopMargin(.01);
 751         p->SetRightMargin(.01);
 752         p->SetFillColor(0);
 753         p->SetFillStyle(0);
 754         p->SetGridx();
 755         stack->Draw("nostack");
 756         stack->GetHistogram()->SetYTitle(stack->GetTitle());
 757         stack->GetHistogram()->SetXTitle("#eta");
 758
 759         TAxis* yaxis = stack->GetHistogram()->GetYaxis();
 760         if (i == 0) yaxis->SetRangeUser(0,20); // Chi2
 761         if (i == 1) stack->SetMaximum(1);      // c
 762         if (i == 2) stack->SetMaximum(1);      // delta
 763         if (i == 3) stack->SetMaximum(0.1);   // xi
 764         if (i == 4) stack->SetMaximum(0.5);    // sigma{,n}
 765         // if (i == 0) p->SetLogy();
 766         yaxis->SetTitleSize(0.3/nL);
 767         yaxis->SetLabelSize(0.08);
 768         yaxis->SetTitleOffset(3/nL);
 769         yaxis->SetNdivisions(5);
 770         yaxis->SetTitleFont(42);
 771         yaxis->SetLabelFont(42);
 772         yaxis->SetDecimals();
 773
 774         TAxis* xaxis = stack->GetHistogram()->GetXaxis();
 775         xaxis->SetTitleSize(0.3/nL);
 776         xaxis->SetLabelSize(0.08);
 777         xaxis->SetTitleOffset(2./nL);
 778         xaxis->SetNdivisions(10);
 779         xaxis->SetTitleFont(42);
 780         xaxis->SetLabelFont(42);
 781         xaxis->SetDecimals();
 782
 783         stack->Draw("nostack");
 784         p->cd();
 785       }
 786       TVirtualPad* p     = GetPad(6);
 787       p->SetFillColor(kWhite);
 788       Double_t x = .3;
 789       Double_t y = .8;
 790       TLatex* l = new TLatex(x, y, "Fits to #Delta (energy loss) spectra");
 791       l->SetTextColor(kBlue+3);
 792       l->SetNDC();
 793       l->Draw();
 794       x = .05;
 795       y -= 2 * 1.2*l->GetTextSize();
 796       l->DrawLatex(x, y, "F(#Delta;c,#Delta_{p},#xi,#sigma)="
 797                    "#frac{c}{#sqrt{2#pi}#sigma}#int_{-#infty}^{#infty}d#Delta'"
 798                    "L(#Delta;#Delta',#xi) G(#Delta_{p};#Delta',#sigma^{2})");
 799       y -= 1.2*l->GetTextSize();
 800       x += .1;
 801       DrawText(l, x, y, "#chi^{2}/#nu", "Goodness of fit",         .2);
 802       DrawText(l, x, y, "c",            "Overall constant",        .2);
 803       DrawText(l, x, y, "#Delta_{p}",   "Most probable value",     .2);
 804       DrawText(l, x, y, "#xi",          "'Width' of Landau (L)",   .2);
 805       DrawText(l, x, y, "#sigma",       "'Width' of Gaussian (G)", .2);
 806
 807       // stack->GetHists()->ls();
 808
 809       PrintCanvas("fitResults", fGlobal->runNo);
 810     }
 811
 812     static_cast<Ring*>(fFMD1i)->ProcessEnergyLoss(fEnergyFitter);
 813     static_cast<Ring*>(fFMD2i)->ProcessEnergyLoss(fEnergyFitter);
 814     static_cast<Ring*>(fFMD2o)->ProcessEnergyLoss(fEnergyFitter);
 815     static_cast<Ring*>(fFMD3i)->ProcessEnergyLoss(fEnergyFitter);
 816     static_cast<Ring*>(fFMD3o)->ProcessEnergyLoss(fEnergyFitter);
 817
 818     return true;
 819   }
 820   /**
 821    * Process the information from the sharing filter
 822    *
 823    *
 824    * @return true on success
 825    */
 826   Bool_t ProcessSharingFilter()
 827   {
 828     // --- Neighbors -------------------------------------------------
 829     MakeCanvasTitle("Correlation of neighboring strips");
 830
 831     for (Int_t i = 1; i <= 6; i++) {
 832       TVirtualPad* p = GetPad(i);
 833
 834       Ring* r = GetRing(i);
 835       if (!r) continue;
 836
 837       r->ProcessNeighbors(fSharingFilter, p);
 838     }
 839     TVirtualPad* p = 0;
 840     if ((p = GetPad(4))) {
 841       p->SetFillColor(kWhite);
 842
 843       TLatex* l = new TLatex(.2, .7, "Gradient: before merging");
 844       l->SetNDC();
 845       l->SetTextColor(kBlue+3);
 846       l->Draw();
 847       l->DrawText(.2, .6, "Boxes: after merging");
 848
 849       fCanvas->cd();
 850       PrintCanvas("neighbors", fGlobal->runNo);
 851     }
 852
 853     // --- 123 -------------------------------------------------------
 854     MakeCanvasTitle("#Delta for singles, doubles, and triples");
 855
 856     for (Int_t i = 1; i <= 6; i++) {
 857       p = GetPad(i);
 858
 859       Ring* r = GetRing(i);
 860       if (!r) continue;
 861
 862       r->Process123(fSharingFilter, p);
 863     }
 864     if ((p = GetPad(4))) {
 865       TLegend* ll = new TLegend(.2, .2, .8, .8);
 866       ll->SetFillColor(0);
 867       ll->SetBorderSize(0);
 868       TLegendEntry* e = ll->AddEntry("dummy", "Singles", "l");
 869       e->SetLineStyle(1);
 870       e = ll->AddEntry("dummy", "Doubles", "l");
 871       e->SetLineStyle(2);
 872       e = ll->AddEntry("dummy", "Triples", "l");
 873       e->SetLineStyle(3);
 874       ll->Draw();
 875
 876       PrintCanvas("123", fGlobal->runNo);
 877     }
 878     return true;
 879   }
 880   /**
 881    * Process the information from the density calculator
 882    *
 883    *
 884    * @return true on success
 885    */
 886   Bool_t ProcessDensityCalculator()
 887   {
 888     // --- ELoss -----------------------------------------------------
 889     MakeCanvasTitle("Energy loss from ESD, after merging, used");
 890
 891     for (Int_t i = 1; i <= 6; i++) {
 892       TVirtualPad* p = GetPad(i);
 893
 894       Ring* r = GetRing(i);
 895       if (!r) continue;
 896
 897       r->ProcessELoss(fSharingFilter, fDensityCalculator, p);
 898     }
 899     TVirtualPad* p = 0;
 900     if ((p = GetPad(4))) {
 901       TLegend* ll = new TLegend(.2, .2, .8, .8);
 902       ll->SetFillColor(0);
 903       ll->SetBorderSize(0);
 904       TLegendEntry* e = ll->AddEntry("dummy", "From ESDs", "l");
 905       e->SetLineStyle(1);
 906       e = ll->AddEntry("dummy", "After Merging", "l");
 907       e->SetLineStyle(2);
 908       e = ll->AddEntry("dummy", "Used", "l");
 909       e->SetLineStyle(3);
 910       ll->Draw();
 911
 912       PrintCanvas("recAna", fGlobal->runNo);
 913     }
 914
 915     // --- Occupancy -------------------------------------------------
 916     MakeCanvasTitle("Occupancy");
 917
 918     for (Int_t i = 1; i <= 6; i++) {
 919       p = GetPad(i);
 920
 921       Ring* r = GetRing(i);
 922       if (!r) continue;
 923
 924       r->ProcessOccupancy(fDensityCalculator, p);
 925     }
 926     if ((p = GetPad(4))) {
 927       TLatex* ltx = new TLatex(.2, .8, "Calculated assuming Poisson stat.");
 928       ltx->SetNDC();
 929       ltx->SetTextColor(kBlue+3);
 930       ltx->Draw();
 931
 932       TObject* etaL = GetObject(fDensityCalculator, "etaLumping");
 933       TObject* phiL = GetObject(fDensityCalculator, "phiLumping");
 934       if (etaL && phiL)
 935         ltx->DrawLatex(.2, .7, Form("Regions of %s strips #times %s sectors",
 936                                     etaL->GetTitle(), phiL->GetTitle()));
 937
 938       PrintCanvas("occupancy", fGlobal->runNo);
 939     }
 940
 941     // --- Correlation of methods ------------------------------------
 942     MakeCanvasTitle("Correlation of N_{ch} methods");
 943
 944     for (Int_t i = 1; i <= 6; i++) {
 945       p = GetPad(i);
 946
 947       Ring* r = GetRing(i);
 948       if (!r) continue;
 949
 950       r->ProcessCorrelation(fDensityCalculator, p);
 951     }
 952     if ((p = GetPad(4))) {
 953       TLatex* ltx = new TLatex(.2, .8, "Correlation of N_{ch} methods");
 954       ltx->SetNDC();
 955       ltx->SetTextColor(kBlue+3);
 956       ltx->Draw();
 957       ltx->DrawLatex(.24, .7, "From #DeltaE fits along x");
 958       ltx->DrawLatex(.24, .6, "From Poisson assumption along y");
 959       ltx->DrawLatex(.24, .4, "Solid line: regression");
 960       ltx->DrawLatex(.24, .3, "Dashed line: x=y to guide the eye");
 961
 962       PrintCanvas("elossVsPoisson", fGlobal->runNo);
 963     }
 964
 965     return true;
 966   }
 967
 968   // --- Utilities ---------------------------------------------------
 969   /**
 970    * Close current file if open, and reset other pointers
 971    *
 972    */
 973   void CloseCurrent()
 974   {
 975     if (fCurrentFile) {
 976       fCurrentFile->Close();
 977       fCurrentFile       = 0;
 978       fSharingFilter     = 0;
 979       fEventInspector    = 0;
 980       fDensityCalculator = 0;
 981       fEnergyFitter      = 0;
 982     }
 983     bool keep = fKeep;
 984     if (fSingle) keep = true;
 985     Close(!keep);
 986   }
 987   /**
 988    * Get the ring corresponding to a pad
 989    *
 990    * @param padNo Pad number
 991    *
 992    * @return Pointer to ring
 993    */
 994   Ring* GetRing(Int_t padNo)
 995   {
 996     QARing* r = 0;
 997     switch(padNo) {
 998     case 1: r = fFMD1i; break;
 999     case 2: r = fFMD2i; break;
1000     case 3: r = fFMD3i; break;
1001     case 5: r = fFMD2o; break;
1002     case 6: r = fFMD3o; break;
1003     }
1004     if (r) return static_cast<Ring*>(r);
1005     return 0;
1006   }
1007   /**
1008    * Get the pad corresponding to a number.  Also clears the pad
1009    *
1010    * @param padNo Pad number to get
1011    *
1012    * @return Pointer to pad or null
1013    */
1014   TVirtualPad* GetPad(Int_t padNo)
1015   {
1016     if (!fCanvas) return 0;
1017     TVirtualPad* p = fCanvas->cd(padNo);
1018     if (!p) return 0;
1019     p->Clear();
1020     p->SetFillColor(kWhite);
1021     return p;
1022   }
1023
1024   /**
1025    * Make canvas title. Canvas is divided into 3x2 pads
1026    *
1027    * @param what Title on canvas
1028    */
1029   void MakeCanvasTitle(const char* what)
1030   {
1031     if (!fCanvas) return;
1032     fCanvas->cd();
1033
1034     CanvasTitle(what);
1035     fCanvas->Divide(3,2,0,0);
1036   }
1037   // --- List utilities ----------------------------------------------
1038   /**
1039    * Get a sub-list from parent list
1040    *
1041    * @param parent Parent list
1042    * @param name   Name of sub-list
1043    *
1044    * @return Pointer to the list
1045    */
1046   static TList* GetSubList(const TList* parent, const char* name)
1047   {
1048     TList* tmp = static_cast<TList*>(parent->FindObject(name));
1049     if (!tmp) {
1050       Error("GetLists", "List %s not found in %s", name,
1051             parent->GetName());
1052       return 0;
1053     }
1054     return tmp;
1055   }
1056   /**
1057    * Get the lists from the file
1058    *
1059    *
1060    * @return true on success
1061    */
1062   Bool_t GetLists()
1063   {
1064     if (!fCurrentFile) return 0;
1065
1066     const char* folder = "ForwardResults";
1067     TList* forward = static_cast<TList*>(fCurrentFile->Get(folder));
1068     if (!forward) {
1069       Error("GetLists", "List %s not found in %s", folder,
1070             fCurrentFile->GetName());
1071       return false;
1072     }
1073
1074     fSharingFilter     = GetSubList(forward, "fmdSharingFilter");
1075     fEventInspector    = GetSubList(forward, "fmdEventInspector");
1076     fDensityCalculator = GetSubList(forward, "fmdDensityCalculator");
1077     fEnergyFitter      = GetSubList(forward, "fmdEnergyFitter");
1078
1079     if (!fSharingFilter)     return false;
1080     if (!fEventInspector)    return false;
1081     if (!fDensityCalculator) return false;
1082     if (!fEnergyFitter)     return false;
1083
1084     return true;
1085   }
1086   /**
1087    * Find and object in a list
1088    *
1089    * @param list List to look in
1090    * @param name Name of object
1091    *
1092    * @return Pointer to object or null
1093    */
1094   static TObject* GetObject(const TList* list, const char* name)
1095   {
1096     if (!list) return 0;
1097     TObject* o = list->FindObject(name);
1098     if (!o) {
1099       Error("GetObject", "Failed to find object %s in %s",
1100             name, list->GetName());
1101       return 0;
1102     }
1103     return o;
1104   }
1105   /**
1106    * Get a histogram from a list
1107    *
1108    * @param list List
1109    * @param name Name of histogram
1110    *
1111    * @return Pointer to object or null
1112    */
1113   static TH1* GetHistogram(const TList* list, const char* name)
1114   {
1115     return static_cast<TH1*>(GetObject(list, name));
1116   }
1117   /**
1118    * Draw some text on a pad
1119    *
1120    * @param l   LaTeX object to use
1121    * @param x   x coordinate
1122    * @param y   y coordinate (incremented on return)
1123    * @param c1  First text
1124    * @param c2  Second text
1125    * @param dx  Distance between c1 start and c2 start
1126    */
1127   static void
1128   DrawText(TLatex* l, Double_t x, Double_t& y, const char* c1, const char* c2,
1129            Double_t dx=.4)
1130   {
1131     y -= 1.2*l->GetTextSize();
1132     l->DrawLatex(x,    y, c1);
1133     l->DrawLatex(x+dx, y, c2);
1134   }
1135
1136
1137   // --- Members -----------------------------------------------------
1138   TFile*  fCurrentFile;          // Current input file
1139   TList*  fSharingFilter;        // Sharing filter list
1140   TList*  fEventInspector;       // Event inspector list
1141   TList*  fDensityCalculator;    // Density calculator list
1142   TList*  fEnergyFitter;         // Energy fitter list
1143   TList   fFiles;                // List of files to process
1144   Bool_t  fKeep;                 // Keep PNGs
1145 };
1146 //
1147 // EOF
1148 //
1149
1150
1151