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