Minor fixes
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / DrawRes.C
index 8b540f49d80689a9b9341d8e314aa7d0257b3042..359f98da92b51cf4713f3e94fd3f92514733a56e 100644 (file)
@@ -13,6 +13,7 @@
 #include <TGraph.h>
 #include <TGraphErrors.h>
 #include <TLatex.h>
+#include <TSystem.h>
 #include "AliAODForwardMult.h"
 #include "OtherData.C"
 
@@ -49,6 +50,8 @@ public:
   Double_t           fVtxEff;
   /** Title to put on the plot */
   TString            fTitle;
+  /** Do HHD comparison */
+  Bool_t             fDoHHD;
 
   //__________________________________________________________________
   /** 
@@ -61,7 +64,8 @@ public:
       fOut(0), 
       fSum(0),
       fVtxEff(0),
-      fTitle("")
+      fTitle(""),
+      fDoHHD(kTRUE)
   {}
   //__________________________________________________________________
   /** 
@@ -102,7 +106,7 @@ public:
   Bool_t Run(const char* file="AliAODs.root", 
             Double_t vzMin=-10, Double_t vzMax=10, Int_t rebin=1,
             Int_t mask=AliAODForwardMult::kInel, Int_t energy=900,
-            const char* title="")
+            const char* title="", Bool_t doHHD=false, Bool_t doComp=false)
   {
     TString trgName;
     if    (mask & AliAODForwardMult::kInel)    trgName.Append("inel-");
@@ -119,7 +123,7 @@ public:
     fTitle = title;
     if (!Open(file, outName)) return kFALSE;
     if (!Process(vzMin,vzMax,mask)) return kFALSE;
-    if (!Finish(rebin, mask, energy)) return kFALSE;
+    if (!Finish(rebin, mask, energy,doHHD,doComp)) return kFALSE;
 
     return kTRUE;
   }
@@ -179,6 +183,11 @@ public:
  
     for (int i = 0; i < nAvailable; i++) { 
       fTree->GetEntry(i);
+      if (((i+1) % 100) == 0) {
+       fprintf(stdout,"Event # %9d of %9d, %9d accepted so far\r", 
+              i+1, nAvailable, nAccepted);
+       fflush(stdout);
+      }
 
       // Create sum histogram on first event - to match binning to input
       if (!fSum) {
@@ -208,6 +217,7 @@ public:
       // Add contribution from this event
       fSum->Add(&(fAOD->GetHistogram()));
     }
+    printf("\n");
     fVtxEff = Double_t(nWithVertex)/nTriggered;
 
     Info("Process", "Total of %9d events\n"
@@ -228,7 +238,8 @@ public:
    * 
    * @return true on success, false otherwise 
    */
-  Bool_t Finish(Int_t rebin, Int_t mask, Int_t energy)
+  Bool_t Finish(Int_t rebin, Int_t mask, Int_t energy, 
+               Bool_t doHHD, Bool_t doComp)
   {
     fOut->cd();
 
@@ -245,95 +256,136 @@ public:
     dndeta->SetMarkerStyle(20);
     dndeta->SetMarkerSize(1);
     dndeta->SetFillStyle(0);
+    Rebin(dndeta, rebin);
+    DrawIt(dndeta, mask, energy, doHHD, doComp);
 
+    return kTRUE;
+  }
+  //__________________________________________________________________
+  /** 
+   */
+  void DrawIt(TH1* dndeta, Int_t mask, Int_t energy, 
+             Bool_t doHHD, Bool_t doComp)
+  {
+    // --- 1st part - prepare data -----------------------------------
     TH1* sym = Symmetrice(dndeta);
-    
-    Rebin(dndeta, rebin);
-    Rebin(sym, rebin);
+    // Rebin(sym, rebin);
 
     THStack* stack = new THStack("results", "Results");
     stack->Add(dndeta);
     stack->Add(sym);
 
-    TString hhdName(fOut->GetName());
-    hhdName.ReplaceAll("dndeta", "hhd");    
-    TH1* hhd    = GetHHD(hhdName.Data(), mask & AliAODForwardMult::kNSD);
+    // Get the data from HHD's analysis - if available 
+    TH1* hhd    = 0;
     TH1* hhdsym = 0;
-    if (hhd) { 
-      hhdsym = Symmetrice(hhd);
-      stack->Add(hhd);
-      stack->Add(hhdsym);
+    Info("DrawIt", "Will %sdraw HHD results", (doHHD ? "" : "not "));
+    if (doHHD) {
+      TString hhdName(fOut->GetName());
+      hhdName.ReplaceAll("dndeta", "hhd");    
+      hhd    = GetHHD(hhdName.Data(), mask & AliAODForwardMult::kNSD);
+      hhdsym = 0;
+      if (hhd) { 
+       // Symmetrice and add to stack 
+       hhdsym = Symmetrice(hhd);
+       stack->Add(hhd);
+       stack->Add(hhdsym);
+      }
+      else 
+       Warning("DrawIt", "No HHD data found");
     }
 
-    TMultiGraph* mg     = GetData(energy, mask);
+    // Get graph of 'other' data - e.g., UA5, CMS, ... - and check if
+    // there's any graphs.  Set the pad division based on that.
+    Info("DrawIt", "Will %sdraw other results", (doComp ? "" : "not "));
+    TMultiGraph* other    = (doComp ? GetData(energy, mask) : 0);
+    THStack*     ratios   = MakeRatios(dndeta, sym, hhd, hhdsym, other);
+
+    // Check if we have ratios 
 
+    // --- 2nd part - Draw in canvas ---------------------------------
+    // Make a canvas 
     gStyle->SetOptTitle(0);
-    Double_t yd = (mg->GetListOfGraphs()->GetEntries() ? 0.3 : 0);
+    gStyle->SetTitleFont(132, "xyz");
+    gStyle->SetLabelFont(132, "xyz");
     TCanvas* c = new TCanvas("c", "C", 900, 800);
     c->SetFillColor(0);
     c->SetBorderMode(0);
     c->SetBorderSize(0);
     c->cd();
 
+    Double_t yd = (ratios ? 0.3 : 0);
+
+    // Make a sub-pad for the result itself
     TPad* p1 = new TPad("p1", "p1", 0, yd, 1.0, 1.0, 0, 0);
     p1->SetTopMargin(0.05);
-    p1->SetBottomMargin(0.001);
+    p1->SetBottomMargin(ratios ? 0.001 : 0.1);
     p1->SetRightMargin(0.05);
     p1->SetGridx();
     p1->SetTicks(1,1);
     p1->Draw();
     p1->cd();
-    stack->SetMinimum(-0.1);
+
+    // Fix the apperance of the stack and redraw. 
+    stack->SetMinimum(ratios ? -0.1 : 0);
     FixAxis(stack, 1/(1-yd)/1.5, "#frac{1}{N} #frac{dN_{ch}}{#eta}");
     p1->Clear();
     stack->DrawClone("nostack e1");
 
-
-    THStack* ratios = new THStack("ratios", "Ratios");
-    TGraphAsymmErrors* o = 0;
-    TIter nextG(mg->GetListOfGraphs());
-    while ((o = static_cast<TGraphAsymmErrors*>(nextG()))) {
-      o->Draw("same p");
-      ratios->Add(Ratio(dndeta, o));
-      ratios->Add(Ratio(sym,    o));
-      ratios->Add(Ratio(hhd,    o));
-      ratios->Add(Ratio(hhdsym, o));
-    }
-    if (hhd && hhdsym) { 
-      TH1F* t1 = static_cast<TH1F*>(dndeta->Clone(Form("%s_%s", 
-                                                      dndeta->GetName(), 
-                                                      hhd->GetName())));
-      t1->SetTitle(Form("%s / %s", dndeta->GetTitle(), hhd->GetTitle()));
-      TH1F* t2 = static_cast<TH1F*>(sym->Clone(Form("%s_%s", 
-                                                   sym->GetName(), 
-                                                   hhdsym->GetName())));
-      t2->SetTitle(Form("%s / %s", sym->GetTitle(), hhdsym->GetTitle()));
-      t1->Divide(hhd);
-      t2->Divide(hhdsym);
-      ratios->Add(t1);
-      ratios->Add(t2);
+    // Draw other data 
+    if (other) {
+      TGraphAsymmErrors* o      = 0;
+      TIter              nextG(other->GetListOfGraphs());
+      while ((o = static_cast<TGraphAsymmErrors*>(nextG()))) 
+       o->Draw("same p");
     }
 
-    // Make a legend
+
+    // Make a legend in the main result pad 
     TString trigs(AliAODForwardMult::GetTriggerString(mask));
-    TLegend* l = p1->BuildLegend(.3,p1->GetBottomMargin()+.01,.7,.4,
-                                Form("#sqrt{s}=%dGeV, %s", 
-                                     energy, trigs.Data()));
+    TLegend* l = p1->BuildLegend(.15,p1->GetBottomMargin()+.01,.90,.35);
+    l->SetNColumns(2);
     l->SetFillColor(0);
     l->SetFillStyle(0);
     l->SetBorderSize(0);
+    l->SetTextFont(132);
     p1->cd();
+
+    // Put a title on top 
+    TLatex* tit = new TLatex(0.10, 0.95, fTitle.Data());
+    tit->SetNDC();
+    tit->SetTextFont(132);
+    tit->SetTextSize(0.05);
+    tit->Draw();
+
+    // Put a nice label in the plot 
+    TLatex* tt = new TLatex(.93, .93, 
+                           Form("#sqrt{s}=%dGeV, %s", energy,
+                                AliAODForwardMult::GetTriggerString(mask)));
+    tt->SetNDC();
+    tt->SetTextFont(132);
+    tt->SetTextAlign(33);
+    tt->Draw();
+
+    // Mark the plot as preliminary 
+    TLatex* pt = new TLatex(.93, .88, "Preliminary");
+    pt->SetNDC();
+    pt->SetTextFont(22);
+    pt->SetTextSize(0.07);
+    pt->SetTextColor(kRed+1);
+    pt->SetTextAlign(33);
+    pt->Draw();
     c->cd();
 
-    if (!ratios->GetHists() || ratios->GetHists()->GetEntries() <= 0) {
-      p1->SetPad(0, 0, 1, 1);
-      p1->SetBottomMargin(0.1);
-      l->SetY1(0.11);
-      FixAxis(stack, (1-yd)/1,  "#frac{1}{N} #frac{dN_{ch}}{#eta}",false);
-      p1->cd();
+    // If we do not have the ratios, fix up the display 
+    // p1->SetPad(0, 0, 1, 1);
+    // p1->SetBottomMargin(0.1);
+    // l->SetY1(0.11);
+    // stack->SetMinimum(0);
+    // FixAxis(stack, (1-yd)/1,  "#frac{1}{N} #frac{dN_{ch}}{#eta}",10,false);
+    if (ratios) {
+      // If we do have the ratios, then make a new pad and draw the 
+      // ratios there 
       c->cd();
-    }
-    else {
       TPad* p2 = new TPad("p2", "p2", 0, 0.0, 1.0, yd, 0, 0);
       p2->SetTopMargin(0.001);
       p2->SetRightMargin(0.05);
@@ -342,16 +394,25 @@ public:
       p2->SetTicks(1,1);
       p2->Draw();
       p2->cd();
-      FixAxis(ratios, 1/yd/1.5, "Ratios");
+
+      // Fix up axis 
+      FixAxis(ratios, 1/yd/1.5, "Ratios", 5);
+
+      // Fix up y range and redraw 
+      ratios->SetMinimum(.58);
+      ratios->SetMaximum(1.22);
       p2->Clear();
       ratios->DrawClone("nostack e1");
       
-      TLegend* l2 = p2->BuildLegend(.3,p2->GetBottomMargin()+.01,.7,.6);
+      // Make a legend 
+      TLegend* l2 = p2->BuildLegend(.15,p2->GetBottomMargin()+.01,.9,.6);
+      l2->SetNColumns(2);
       l2->SetFillColor(0);
       l2->SetFillStyle(0);
       l2->SetBorderSize(0);
-      
-      p2->cd();
+      l2->SetTextFont(132);
+
+      // Make a nice band from 0.9 to 1.1 
       TGraphErrors* band = new TGraphErrors(2);
       band->SetPoint(0, sym->GetXaxis()->GetXmin(), 1);
       band->SetPoint(1, dndeta->GetXaxis()->GetXmax(), 1);
@@ -360,27 +421,26 @@ public:
       band->SetFillColor(kYellow+2);
       band->SetFillStyle(3002);
       band->Draw("3 same");
+
+      // Replot the ratios on top 
       ratios->DrawClone("nostack e1 same");
 
       c->cd();
     }
-    p1->cd();
-    TLatex* tit = new TLatex(0.10, 0.95, fTitle.Data());
-    tit->SetNDC();
-    tit->SetTextSize(0.05);
-    tit->Draw();
     
+    // Plot to disk
     TString imgName(fOut->GetName());
     imgName.ReplaceAll(".root", ".png");
     c->SaveAs(imgName.Data());
+    imgName.ReplaceAll(".png", ".C");
+    c->SaveAs(imgName.Data());
 
     stack->Write();
-    mg->Write();
-    ratios->Write();
+    if (other)  other->Write();
+    if (ratios) ratios->Write();
 
+    // Close our file 
     fOut->Close();
-
-    return kTRUE;
   }
   //__________________________________________________________________
   /** 
@@ -393,6 +453,10 @@ public:
   TH1* GetHHD(const char* fn="fmd_dNdeta_mult.root", bool nsd=false) 
   {
     TDirectory* savdir = gDirectory;
+    if (gSystem->AccessPathName(fn)) { 
+      Warning("GetHHD", "Output of HHD analysis (%s) not available", fn);
+      return 0;
+    }
     TFile* file = TFile::Open(fn, "READ");
     if (!file) { 
       Warning("GetHHD", "couldn't open HHD file %s", fn);
@@ -404,7 +468,6 @@ public:
     if (!h) { 
       Warning("GetHHD", "Couldn't find HHD histogram %s in %s", 
              hist.Data(), fn);
-      file->ls();
       file->Close();
       savdir->cd();
       return 0;
@@ -421,6 +484,58 @@ public:
     savdir->cd();
     return r;
   }
+  //__________________________________________________________________
+  /** 
+   */ 
+  THStack* MakeRatios(const TH1* dndeta, const TH1* sym, 
+                     const TH1* hhd,    const TH1* hhdsym, 
+                     TMultiGraph* other) const 
+  {
+    // If we have 'other' data, then do the ratio of the results to that
+    Bool_t hasOther = (other && other->GetListOfGraphs() && 
+                      other->GetListOfGraphs()->GetEntries() > 0);
+    Bool_t hasHhd   = (hhd && hhdsym);
+    if (!hasOther || !hasHhd) return 0;
+
+    THStack* ratios = new THStack("ratios", "Ratios");
+    if (hasOther) {
+      TGraphAsymmErrors* o      = 0;
+      TIter              nextG(other->GetListOfGraphs());
+      while ((o = static_cast<TGraphAsymmErrors*>(nextG()))) {
+       ratios->Add(Ratio(dndeta, o));
+       ratios->Add(Ratio(sym,    o));
+       ratios->Add(Ratio(hhd,    o));
+       ratios->Add(Ratio(hhdsym, o));
+      }
+    }
+
+    // If we have data from HHD's analysis, then do the ratio of 
+    // our result to that data. 
+    if (hasHhd) { 
+      TH1F* t1 = static_cast<TH1F*>(dndeta->Clone(Form("%s_%s", 
+                                                      dndeta->GetName(), 
+                                                      hhd->GetName())));
+      t1->SetTitle(Form("%s / %s", dndeta->GetTitle(), hhd->GetTitle()));
+      TH1F* t2 = static_cast<TH1F*>(sym->Clone(Form("%s_%s", 
+                                                   sym->GetName(), 
+                                                   hhdsym->GetName())));
+      t2->SetTitle(Form("%s / %s", sym->GetTitle(), hhdsym->GetTitle()));
+      t1->Divide(hhd);
+      t2->Divide(hhdsym);
+      ratios->Add(t1);
+      ratios->Add(t2);
+    }
+
+    // Check if we have ratios 
+    Bool_t   hasRatios = (ratios->GetHists() && 
+                         (ratios->GetHists()->GetEntries() > 0));
+    Info("MakeRatios", "Got a total of %d ratios", !hasRatios ? 0 :
+        ratios->GetHists()->GetEntries());
+
+    if (!hasRatios) { delete ratios; ratios = 0; }
+    return ratios;
+  }
+
   //__________________________________________________________________
   /** 
    * Fix the apperance of the axis in a stack 
@@ -431,7 +546,7 @@ public:
    * @param force  Whether to draw the stack first or not 
    */
   void FixAxis(THStack* stack, Double_t s, const char* ytitle, 
-              Bool_t force=true) 
+              Int_t ynDiv=10, Bool_t force=true) 
   {
     if (force) stack->Draw("nostack e1");
 
@@ -451,8 +566,9 @@ public:
     }
     if (ya) { 
       ya->SetTitle(ytitle);
+      ya->SetDecimals();
       // ya->SetTicks("+-");
-      ya->SetNdivisions(10);
+      ya->SetNdivisions(ynDiv);
       ya->SetTitleSize(s*ya->GetTitleSize());
       ya->SetLabelSize(s*ya->GetLabelSize());
     }      
@@ -477,7 +593,7 @@ public:
     ret->Reset();
     ret->SetMarkerStyle(h->GetMarkerStyle());
     ret->SetMarkerColor(g->GetMarkerColor());
-    ret->SetMarkerSize(0.7*g->GetMarkerSize());
+    ret->SetMarkerSize(0.9*g->GetMarkerSize());
     Double_t xlow  = g->GetX()[0];
     Double_t xhigh = g->GetX()[g->GetN()-1];
     if (xlow > xhigh) { Double_t t = xhigh; xhigh = xlow; xlow = t; }