]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG2/FORWARD/analysis2/DrawdNdeta.C
Mior fixes
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / DrawdNdeta.C
index 05e83e1648c2667133c3064cc4ce8892f03e80ae..7883d71739605ddd1dc30f5c0077ef1035323835 100644 (file)
@@ -3,15 +3,18 @@
  * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
  * @date   Wed Mar 23 14:07:10 2011
  * 
- * @brief  Script to visualise the dN/deta 
+ * @brief  Script to visualise the dN/deta for pp and PbPb
  *
  * This script is independent of any AliROOT code - and should stay
  * that way.
  * 
+ * The script is <emph>very</emph> long - sigh - the joy of drawing
+ * things nicely in ROOT
  * 
  * @ingroup pwg2_forward_dndeta
  */
 #include <TH1.h>
+#include <TColor.h>
 #include <THStack.h>
 #include <TGraphErrors.h>
 #include <TGraphAsymmErrors.h>
 #include <TLegendEntry.h>
 #include <TLatex.h>
 #include <TImage.h>
+#include <TRandom.h>
+#include <fstream>
+#define SYSERR_COLOR kBlue-10
+#define SYSERR_STYLE 1001
+
+Double_t myFunc(Double_t* xp, Double_t* pp);
 
 /**
  * Class to draw dN/deta results 
@@ -56,21 +65,35 @@ struct dNdetaDrawer
    * 
    */
   dNdetaDrawer()
-    : fShowOthers(false),      // Bool_t 
-      fShowAlice(false),        // Bool_t
-      fShowRatios(false),      // Bool_t 
-      fShowLeftRight(false),   // Bool_t
-      fRebin(5),               // UShort_t
-      fCutEdges(false),                // Bool_t
-      fTitle(""),              // TString
-      fTrigString(0),          // TNamed*
-      fSNNString(0),           // TNamed*
-      fSysString(0),           // TNamed*
-      fVtxAxis(0),             // TAxis*
-      fCentAxis(0),             // TAxis*
-      fTriggers(0),            // TH1*
-      fRangeParam(0)
-
+    : // Options 
+    fShowRatios(false),    // Show ratios 
+    fShowLeftRight(false), // Show asymmetry 
+    fShowRings(false),     // Show rings too
+    fExport(false),        // Export data to script
+    fCutEdges(false),      // Whether to cut edges
+    fRemoveOuters(false),  // Whether to remove outers
+    fShowOthers(0),        // Show other data
+    // Settings 
+    fRebin(0),             // Rebinning factor 
+    fFwdSysErr(0.076),     // Systematic error in forward range
+    fCenSysErr(0),         // Systematic error in central range 
+    fTitle(""),            // Title on plot
+    fClusterScale(""),     // Scaling of clusters to tracklets      
+    // Read (or set) information 
+    fTrigString(0),        // Trigger string (read, or set)
+    fNormString(0),        // Normalisation string (read, or set)
+    fSNNString(0),         // Energy string (read, or set)
+    fSysString(0),         // Collision system string (read or set)
+    fVtxAxis(0),           // Vertex cuts (read or set)
+    fCentAxis(0),          // Centrality axis
+    // Resulting plots 
+    fResults(0),           // Stack of results 
+    fRatios(0),            // Stack of ratios 
+    fLeftRight(0),         // Left-right asymmetry
+    fOthers(0),            // Older data 
+    fTriggers(0),          // Number of triggers
+    fTruth(0),             // Pointer to truth 
+    fRangeParam(0)         // Parameter object for range zoom 
   {
     fRangeParam = new RangeParam;
     fRangeParam->fMasterAxis = 0;
@@ -79,6 +102,9 @@ struct dNdetaDrawer
     fRangeParam->fSlave2Axis = 0;
     fRangeParam->fSlave2Pad  = 0;
   }
+  dNdetaDrawer(const dNdetaDrawer&) {}
+  dNdetaDrawer& operator=(const dNdetaDrawer&) { return *this; }
+
   //__________________________________________________________________
   virtual ~dNdetaDrawer()
   {
@@ -107,14 +133,7 @@ struct dNdetaDrawer
    * 
    * @param x Whether to show or not 
    */
-  void SetShowOthers(Bool_t x)    { fShowOthers = x; }
-  //__________________________________________________________________
-  /** 
-   * Show ALICE published data 
-   * 
-   * @param x Wheter to show or not 
-   */
-  void SetShowAlice(Bool_t x)     { fShowAlice = x; }
+  void SetShowOthers(UShort_t x)    { fShowOthers = x; }
   //__________________________________________________________________
   /** 
    * Whether to show ratios or not.  If there's nothing to compare to,
@@ -132,6 +151,20 @@ struct dNdetaDrawer
    */
   void SetShowLeftRight(Bool_t x) { fShowLeftRight = x; }
   //__________________________________________________________________
+  /** 
+   * Whether to show rings 
+   * 
+   * @param x To show or not 
+   */
+  void SetShowRings(Bool_t x) { fShowRings = x; }
+  //__________________________________________________________________
+  /** 
+   * Whether to export results to a script 
+   *
+   * @param x Wheter to export results to a script
+   */
+  void SetExport(Bool_t x)     { fExport = x; }
+  //__________________________________________________________________
   /** 
    * Set the rebinning factor 
    * 
@@ -152,6 +185,20 @@ struct dNdetaDrawer
    * @param x Title
    */
   void SetTitle(TString x)        { fTitle = x; }
+  //__________________________________________________________________
+  /** 
+   * Set the systematic error in the forward region
+   * 
+   * @param e Systematic error in the forward region 
+   */
+  void SetForwardSysError(Double_t e=0) { fFwdSysErr = e; }
+  //__________________________________________________________________
+  /** 
+   * Set the systematic error in the forward region
+   * 
+   * @param e Systematic error in the forward region 
+   */
+  void SetCentralSysError(Double_t e=0) { fCenSysErr = e; }
   /* @} */
   //==================================================================  
   /** 
@@ -245,6 +292,10 @@ struct dNdetaDrawer
     TList* clusters = static_cast<TList*>(file->Get("CentralResults"));
     if (!clusters) Warning("Open", "Couldn't find list CentralResults");
 
+    // --- Get the central results -----------------------------------
+    TList* mcTruth = static_cast<TList*>(file->Get("MCTruthResults"));
+    if (!mcTruth) Warning("Open", "Couldn't find list MCTruthResults");
+
     // --- Make our containtes ---------------------------------------
     fResults   = new THStack("results", "Results");
     fRatios    = new THStack("ratios",  "Ratios");
@@ -252,8 +303,9 @@ struct dNdetaDrawer
     fOthers    = new TMultiGraph();
     
     // --- Loop over input data --------------------------------------
-    FetchResults(forward,  "Forward", max, rmax, amax);
-    FetchResults(clusters, "Central", max, rmax, amax);
+    /*TObjArray* mcA =*/ FetchResults(mcTruth,  "MCTruth", max, rmax, amax);
+    TObjArray* fwdA = FetchResults(forward,  "Forward", max, rmax, amax);
+    TObjArray* cenA = FetchResults(clusters, "Central", max, rmax, amax);
 
     // --- Get trigger information -----------------------------------
     TList* sums = static_cast<TList*>(file->Get("ForwardSums"));
@@ -282,7 +334,7 @@ struct dNdetaDrawer
     if (!fOthers->GetListOfGraphs() || 
        fOthers->GetListOfGraphs()->GetEntries() <= 0) { 
       Warning("Run", "No other data found - disabling that");
-      fShowOthers = false;
+      fShowOthers = 0;
     }
     if (!fRatios->GetHists() || 
        fRatios->GetHists()->GetEntries() <= 0) { 
@@ -296,6 +348,32 @@ struct dNdetaDrawer
       // fLeftRight->ls();
       fShowLeftRight = false;
     }
+    if (fFwdSysErr > 0) { 
+      if (fCenSysErr <= 0) fCenSysErr = fFwdSysErr;
+      for (Int_t i = 0; i < fwdA->GetEntriesFast(); i++) {
+       TH1* fwd = static_cast<TH1*>(fwdA->At(i));
+       TH1* cen = static_cast<TH1*>(cenA->At(i));
+       CorrectForward(fwd);
+       CorrectCentral(cen);
+       Double_t low, high;
+       TH1* tmp = Merge(cen, fwd, low, high);
+       TF1* f   = FitMerged(tmp, low, high);
+       MakeSysError(tmp, cen, fwd, f);
+       delete f;
+       Info("", "Adding systematic error histogram %s", 
+            tmp->GetName());
+       fResults->GetHists()->AddFirst(tmp, "e5");
+       TH1* tmp2 = Symmetrice(tmp);
+       tmp2->SetFillColor(tmp->GetFillColor());
+       tmp2->SetFillStyle(tmp->GetFillStyle());
+       tmp2->SetMarkerStyle(tmp->GetMarkerStyle());
+       tmp2->SetLineWidth(tmp->GetLineWidth());
+       fResults->GetHists()->AddFirst(tmp2, "e5");
+       fResults->Modified();
+      }
+    }
+    delete fwdA;
+    delete cenA;
     
     // --- Close the input file --------------------------------------
     file->Close();
@@ -315,7 +393,9 @@ struct dNdetaDrawer
   void FetchInformation(const TList* results)
   {
     if (!fTrigString) 
-      fTrigString = static_cast<TNamed*>(results->FindObject("trigString"));
+      fTrigString = static_cast<TNamed*>(results->FindObject("trigger"));
+    if (!fNormString) 
+      fNormString = static_cast<TNamed*>(results->FindObject("scheme"));
     if (!fSNNString) 
       fSNNString  = static_cast<TNamed*>(results->FindObject("sNN"));
     if (!fSysString) 
@@ -324,7 +404,10 @@ struct dNdetaDrawer
       fVtxAxis    = static_cast<TAxis*>(results->FindObject("vtxAxis"));
     if (!fCentAxis) 
       fCentAxis   = static_cast<TAxis*>(results->FindObject("centAxis"));
-    if (!fTrigString) fTrigString = new TNamed("trigString", "unknown");
+
+    TNamed* options = static_cast<TAxis*>(results->FindObject("options"));
+    if (!fTrigString) fTrigString = new TNamed("trigger", "unknown");
+    if (!fNormString) fNormString = new TNamed("scheme", "unknown");
     if (!fSNNString)  fSNNString  = new TNamed("sNN", "unknown");
     if (!fSysString)  fSysString  = new TNamed("sys", "unknown");
     if (!fVtxAxis) { 
@@ -333,48 +416,46 @@ struct dNdetaDrawer
       fVtxAxis->SetTitle("v_{z} range unspecified");
     }
 
-    TString centTxt;
+    TString centTxt("none");
     if (fCentAxis) { 
       Int_t nCent = fCentAxis->GetNbins();
-      centTxt = Form("\n   Centrality: %d bins", nCent);
+      centTxt = Form("%d bins", nCent);
       for (Int_t i = 0; i <= nCent; i++) 
-       centTxt.Append(Form("%c%d", i == 0 ? ' ' : ',', 
+       centTxt.Append(Form("%c%d", i == 0 ? ' ' : '-', 
                            int(fCentAxis->GetXbins()->At(i))));
     }
     Info("FetchInformation", 
         "Initialized for\n"
-        "   Trigger:    %s  (%d)\n"
-        "   sqrt(sNN):  %s  (%dGeV)\n"
-        "   System:     %s  (%d)\n"
-        "   Vz range:   %s  (%f,%f)%s",
+        "   Trigger:       %-30s  (%d)\n"
+        "   sqrt(sNN):     %-30s  (%dGeV)\n"
+        "   System:        %-30s  (%d)\n"
+        "   Vz range:      %-30s  (%f,%f)\n"
+        "   Normalization: %-30s  (%d)\n"
+        "   Centrality:    %s\n"
+        "   Options:       %s",
         fTrigString->GetTitle(), fTrigString->GetUniqueID(), 
         fSNNString->GetTitle(),  fSNNString->GetUniqueID(), 
         fSysString->GetTitle(),  fSysString->GetUniqueID(), 
         fVtxAxis->GetTitle(), fVtxAxis->GetXmin(), fVtxAxis->GetXmax(),
-        centTxt.Data());
+        fNormString->GetTitle(), fNormString->GetUniqueID(),
+        centTxt.Data(), (options ? options->GetTitle() : "none"));
   }
   //__________________________________________________________________
   TMultiGraph* FetchOthers(UShort_t centLow, UShort_t centHigh)
   {
     TMultiGraph* thisOther = 0;
-    if (!fShowOthers && !fShowAlice) return 0;
+    if (fShowOthers == 0) return 0;
 
-    Bool_t   onlya = (fShowOthers ? false : true);
     UShort_t sys   = (fSysString  ? fSysString->GetUniqueID() : 0);
     UShort_t trg   = (fTrigString ? fTrigString->GetUniqueID() : 0);
     UShort_t snn   = (fSNNString  ? fSNNString->GetUniqueID() : 0);
     Long_t   ret   = gROOT->ProcessLine(Form("GetData(%d,%d,%d,%d,%d,%d);",
                                             sys,snn,trg,
-                                            centLow,centHigh,onlya));
-    if (!ret) { 
-      Warning("FetchResults", "No other data found for sys=%d, sNN=%d, "
-             "trigger=%d %d%%-%d%% central %s", 
-             sys, snn, trg, centLow, centHigh, 
-             onlya ? " (ALICE results)" : "all");
-      return 0;
-    }
-    thisOther = reinterpret_cast<TMultiGraph*>(ret);
-    
+                                            centLow,centHigh,
+                                            fShowOthers));
+    if (!ret) return 0;
+
+    thisOther = reinterpret_cast<TMultiGraph*>(ret);    
     return thisOther;
   }
   //__________________________________________________________________
@@ -387,46 +468,61 @@ struct dNdetaDrawer
    * @param rmax  On return, maximum of ratios
    * @param amax  On return, maximum of left-right comparisons
    */
-  void FetchResults(const TList* list, 
-                   const char*  name, 
-                   Double_t&    max,
-                   Double_t&    rmax,
-                   Double_t&    amax)
+  TObjArray* 
+  FetchResults(const TList* list, 
+              const char*  name, 
+              Double_t&    max,
+              Double_t&    rmax,
+              Double_t&    amax)
   {
-    UShort_t n = fCentAxis ? fCentAxis->GetNbins() : 0;
+    UShort_t   n = fCentAxis ? fCentAxis->GetNbins() : 0;
     if (n == 0) {
       TList* all = static_cast<TList*>(list->FindObject("all"));
-      if (!all)
+      if (!all) {
        Error("FetchResults", "Couldn't find list 'all' in %s", 
              list->GetName());
-      else 
-       FetchResults(all, name, FetchOthers(0,0), -1, 0, max, rmax, amax);
-      return;
+       return 0;
+      }
+      TObjArray* a = new TObjArray;
+      TH1*       h = FetchResults(all, name, FetchOthers(0,0), 
+                                 -1000, 0, max, rmax, amax);
+      a->AddAt(h, 0);
+      return a;
     }
     
-    Int_t   nCol = gStyle->GetNumberOfColors();
+    TObjArray* a = new TObjArray;
     for (UShort_t i = 0; i < n; i++) { 
-      UShort_t centLow  = fCentAxis->GetXbins()->At(i);
-      UShort_t centHigh = fCentAxis->GetXbins()->At(i+1);
+      UShort_t centLow  = fCentAxis->GetBinLowEdge(i+1);
+      UShort_t centHigh = fCentAxis->GetBinUpEdge(i+1);
       TString  lname    = Form("cent%03d_%03d", centLow, centHigh);
       TList*   thisCent = static_cast<TList*>(list->FindObject(lname));
-
-      Float_t fc   = (centLow+double(centHigh-centLow)/2) / 100;
-      Int_t   icol = TMath::Min(nCol-1,int(fc * nCol + .5));
-      Int_t   col  = gStyle->GetColorPalette(icol);
-      Info("FetchResults","Centrality %d-%d color index %d (=%d*%f) -> %d", 
-          centLow, centHigh, icol, nCol, fc, col);
+      Int_t    col      = GetCentralityColor(i+1);
 
       TString centTxt = Form("%3d%%-%3d%% central", centLow, centHigh);
-      if (!thisCent)
+      if (!thisCent) {
        Error("FetchResults", "Couldn't find list '%s' in %s", 
              lname.Data(), list->GetName());
-      else 
-       FetchResults(thisCent, name, FetchOthers(centLow, centHigh), 
-                    col, centTxt.Data(), max, rmax, amax);
+       continue;
+      }
+      TH1* h = FetchResults(thisCent, name, FetchOthers(centLow, centHigh), 
+                           col, centTxt.Data(), max, rmax, amax);
+      a->AddAt(h, i);
     }
+    return a;
   } 
   //__________________________________________________________________
+  Int_t GetCentralityColor(Int_t bin) const
+  {
+    UShort_t centLow  = fCentAxis->GetBinLowEdge(bin);
+    UShort_t centHigh = fCentAxis->GetBinUpEdge(bin);
+    Float_t  fc       = (centLow+double(centHigh-centLow)/2) / 100;
+    Int_t    nCol     = gStyle->GetNumberOfColors();
+    Int_t    icol     = TMath::Min(nCol-1,int(fc * nCol + .5));
+    Int_t    col      = gStyle->GetColorPalette(icol);
+    //Info("GetCentralityColor","%3d: %3d-%3d -> %3d",bin,centLow,centHigh,col);
+    return col;
+  }
+  //__________________________________________________________________
   void SetAttributes(TH1* h, Int_t color)
   {
     if (!h) return;
@@ -445,10 +541,11 @@ struct dNdetaDrawer
     // g->SetFillColor(color);
   }
   //__________________________________________________________________
-  void ModifyTitle(TNamed* h, const char* centTxt)
+  void ModifyTitle(TNamed* /*h*/, const char* /*centTxt*/)
   {
-    if (!centTxt || !h) return;
-    h->SetTitle(Form("%s, %s", h->GetTitle(), centTxt));
+    return;
+    // if (!centTxt || !h) return;
+    // h->SetTitle(Form("%s, %s", h->GetTitle(), centTxt));
   }
 
   //__________________________________________________________________
@@ -464,7 +561,7 @@ struct dNdetaDrawer
    * @param rmax       On return, ratio maximum 
    * @param amax       On return, left-right maximum 
    */
-  void FetchResults(const TList* list, 
+  TH1* FetchResults(const TList* list, 
                    const char*  name, 
                    TMultiGraph* thisOther,
                    Int_t        color,
@@ -478,64 +575,77 @@ struct dNdetaDrawer
     TH1* dndetaTruth = FetchResult(list, "dndetaTruth");
     TH1* dndetaSym   = 0;
     TH1* dndetaMCSym = 0;
-    TH1* tracks      = FetchResult(list, "tracks");
-    if (tracks) tracks->SetTitle("ALICE Tracks");
     SetAttributes(dndeta,     color);
-    SetAttributes(dndetaMC,   color+2);
+    SetAttributes(dndetaMC,   fCentAxis ? color : color+2);
     SetAttributes(dndetaTruth,color);
     SetAttributes(dndetaSym,  color);
-    SetAttributes(dndetaMCSym,color+2);
-    SetAttributes(tracks,     color+3);
+    SetAttributes(dndetaMCSym,fCentAxis ? color : color+2);
+    if (dndetaMC && fCentAxis) 
+      dndetaMC->SetMarkerStyle(dndetaMC->GetMarkerStyle()+2);
+    if (dndetaMCSym && fCentAxis) 
+      dndetaMCSym->SetMarkerStyle(dndetaMCSym->GetMarkerStyle()+2);
+    if (dndetaTruth && fCentAxis) {
+      dndetaTruth->SetMarkerStyle(34);
+      dndetaTruth->SetMarkerColor(kYellow-1);
+    }
     ModifyTitle(dndeta,     centTxt);
     ModifyTitle(dndetaMC,   centTxt);
     ModifyTitle(dndetaTruth,centTxt);
     ModifyTitle(dndetaSym,  centTxt);
     ModifyTitle(dndetaMCSym,centTxt);
-    ModifyTitle(tracks,     centTxt);
       
 
     max = TMath::Max(max, AddHistogram(fResults, dndetaTruth, "e5 p"));
     max = TMath::Max(max, AddHistogram(fResults, dndetaMC,    dndetaMCSym));
     max = TMath::Max(max, AddHistogram(fResults, dndeta,      dndetaSym));
-    max = TMath::Max(max, AddHistogram(fResults, tracks));
-    
-    // Info("FetchResults", "Got %p, %p, %p from %s with name %s, max=%f", 
-    //      dndeta, dndetaMC, dndetaTruth, list->GetName(), name, max);
 
-    if (fShowLeftRight) {
-      fLeftRight->Add(Asymmetry(dndeta,    amax));
-      fLeftRight->Add(Asymmetry(dndetaMC,  amax));
-      fLeftRight->Add(Asymmetry(tracks,    amax));
+    if (dndetaTruth) {
+      fTruth = dndetaTruth;
     }
-
-    if (thisOther) {
-      TIter next(thisOther->GetListOfGraphs());
-      TGraph* g = 0;
-      while ((g = static_cast<TGraph*>(next()))) {
-       fRatios->Add(Ratio(dndeta,    g, rmax));
-       fRatios->Add(Ratio(dndetaSym, g, rmax));
-       SetAttributes(g, color);
-       ModifyTitle(g, centTxt);
-       if (!fOthers->GetListOfGraphs() || 
-           !fOthers->GetListOfGraphs()->FindObject(g->GetName())) {
-         max = TMath::Max(max,TMath::MaxElement(g->GetN(), g->GetY()));
-         fOthers->Add(g);
+    else {
+      if (fShowRings) {
+       THStack* rings = static_cast<THStack*>(list->FindObject("dndetaRings"));
+       if (rings) { 
+         TIter next(rings->GetHists());
+         TH1*  hist = 0;
+         while ((hist = static_cast<TH1*>(next()))) 
+           max = TMath::Max(max, AddHistogram(fResults, hist));
+       }
+      }
+      // Info("FetchResults", "Got %p, %p, %p from %s with name %s, max=%f", 
+      //      dndeta, dndetaMC, dndetaTruth, list->GetName(), name, max);
+      
+      if (fShowLeftRight) {
+       fLeftRight->Add(Asymmetry(dndeta,    amax));
+       fLeftRight->Add(Asymmetry(dndetaMC,  amax));
+      }
+      
+      if (thisOther) {
+       TIter next(thisOther->GetListOfGraphs());
+       TGraph* g = 0;
+       while ((g = static_cast<TGraph*>(next()))) {
+         fRatios->Add(Ratio(dndeta,    g, rmax));
+         fRatios->Add(Ratio(dndetaSym, g, rmax));
+         SetAttributes(g, color);
+         ModifyTitle(g, centTxt);
+         if (!fOthers->GetListOfGraphs() || 
+             !fOthers->GetListOfGraphs()->FindObject(g->GetName())) {
+           max = TMath::Max(max,TMath::MaxElement(g->GetN(), g->GetY()));
+           fOthers->Add(g);
+         }
        }
+       // fOthers->Add(thisOther);
       }
-      // fOthers->Add(thisOther);
     }
-    if (tracks) { 
-      if (!fRatios->GetHists() || 
-         !fRatios->GetHists()->FindObject(tracks->GetName()))
-       fRatios->Add(Ratio(dndeta, tracks, rmax));
+    if (dndetaMC) { 
+      fRatios->Add(Ratio(dndeta,    dndetaMC,    rmax));
+      fRatios->Add(Ratio(dndetaSym, dndetaMCSym, rmax));
     }
-
-    if (dndetaTruth) { 
-      fRatios->Add(Ratio(dndeta,      dndetaTruth, rmax));
-      fRatios->Add(Ratio(dndetaSym,   dndetaTruth, rmax));
-      fRatios->Add(Ratio(dndetaMC,    dndetaTruth, rmax));
-      fRatios->Add(Ratio(dndetaMCSym, dndetaTruth, rmax));
+    if (fTruth) {
+      fRatios->Add(Ratio(dndeta,      fTruth, rmax));
+      fRatios->Add(Ratio(dndetaSym,   fTruth, rmax));
     }
+    return dndeta;
   }
   //__________________________________________________________________
   /** 
@@ -557,9 +667,9 @@ struct dNdetaDrawer
     Double_t y1 = 0;
     Double_t y2 = 0;
     Double_t y3 = 0;
-    if (!fShowRatios)    w  *= 1.4;
+    if (!fShowRatios)    w  *= 1.3;
     else                 y1 =  0.3;
-    if (!fShowLeftRight) w  *= 1.4;
+    if (!fShowLeftRight) w  *= 1.3;
     else { 
       Double_t y11 = y1;
       y1 = (y11 > 0.0001 ? 0.4 : 0.2);
@@ -570,12 +680,6 @@ struct dNdetaDrawer
     c->SetBorderSize(0);
     c->SetBorderMode(0);
 
-#if 1
-    Info("Plot", "y1=%f, y2=%f, y3=%f extra: %s %s", 
-        y1, y2, y2, (fShowRatios ? "ratios" : ""), 
-        (fShowLeftRight ? "right/left" : ""));
-    Info("Plot", "Maximum is %f", max);
-#endif
     PlotResults(max, y1);
     c->cd();
 
@@ -592,6 +696,8 @@ struct dNdetaDrawer
     Int_t   nev  = 0;
     if (fTriggers) nev = fTriggers->GetBinContent(1);
     trg          = trg.Strip(TString::kBoth);
+    trg.ReplaceAll(" ", "_");
+    trg.ReplaceAll(">", "Gt");
     TString base(Form("dndeta_%s_%s_%s_%c%02d%c%02dcm_%09dev",
                      fSysString->GetTitle(), 
                      fSNNString->GetTitle(), 
@@ -602,6 +708,8 @@ struct dNdetaDrawer
     c->SaveAs(Form("%s.png",  base.Data()));
     c->SaveAs(Form("%s.root", base.Data()));
     c->SaveAs(Form("%s.C",    base.Data()));
+    base.ReplaceAll("dndeta", "export");
+    Export(base);
   }
   //__________________________________________________________________
   /** 
@@ -618,31 +726,91 @@ struct dNdetaDrawer
                   Double_t x1, Double_t y1, Double_t x2, Double_t y2)
   {
     TLegend* l = new TLegend(x1,y1,x2,y2);
-    l->SetNColumns(2);
+    l->SetNColumns(fCentAxis ? 1 : 2);
     l->SetFillColor(0);
     l->SetFillStyle(0);
     l->SetBorderSize(0);
     l->SetTextFont(132);
 
+    // Loop over items in stack and get unique items, while ignoring
+    // mirrored data and systematic error bands 
     TIter    next(stack->GetHists());
-    TObject* hist = 0;
-    while ((hist = next())) { 
-      TString n(hist->GetTitle());
-      if (n.Contains("mirrored")) continue;
-      l->AddEntry(hist, hist->GetTitle(), "pl");
+    TH1*     hist = 0;
+    TObjArray unique;
+    unique.SetOwner();
+    Bool_t   sysErrSeen = false;
+    while ((hist = static_cast<TH1*>(next()))) { 
+      TString t(hist->GetTitle());
+      TString n(hist->GetName());
+      n.ToLower();
+      if (t.Contains("mirrored")) continue;
+      if (n.Contains("syserror")) { sysErrSeen = true; continue; }
+      if (unique.FindObject(t.Data())) continue;
+      TObjString* s = new TObjString(hist->GetTitle());
+      s->SetUniqueID(((hist->GetMarkerStyle() & 0xFFFF) << 16) |
+                    ((hist->GetMarkerColor() & 0xFFFF) <<  0));
+      unique.Add(s);
+      // l->AddEntry(hist, hist->GetTitle(), "pl");
     }
     if (mg) {
+      // If we have other stuff, scan for unique names 
       TIter nexto(mg->GetListOfGraphs());
-      while ((hist = nexto())) { 
-       TString n(hist->GetTitle());
+      TGraph* g = 0;
+      while ((g = static_cast<TGraph*>(nexto()))) { 
+       TString n(g->GetTitle());
        if (n.Contains("mirrored")) continue;
-       l->AddEntry(hist, hist->GetTitle(), "pl");
+       if (unique.FindObject(n.Data())) continue;
+       TObjString* s = new TObjString(n);
+       s->SetUniqueID(((g->GetMarkerStyle() & 0xFFFF) << 16) |
+                      ((g->GetMarkerColor() & 0xFFFF) <<  0));
+       unique.Add(s);
+       // l->AddEntry(hist, hist->GetTitle(), "pl");
       }
     }
+
+    // Add legend entries for unique items only
+    TIter nextu(&unique);
+    TObject* s = 0;
+    Int_t i = 0;
+    while ((s = nextu())) { 
+      TLegendEntry* dd = l->AddEntry(Form("data%2d", i++), 
+                                    s->GetName(), "lp");
+      Int_t style = (s->GetUniqueID() >> 16) & 0xFFFF;
+      Int_t color = (s->GetUniqueID() >>  0) & 0xFFFF;
+      dd->SetLineColor(kBlack);
+      if (fCentAxis) dd->SetMarkerColor(kBlack);
+      else           dd->SetMarkerColor(color);
+      dd->SetMarkerStyle(style);
+    }
+    if (sysErrSeen) {
+      // Add entry for systematic errors 
+      TLegendEntry* d0 = l->AddEntry("d0", Form("%4.1f%% Systematic error", 
+                                               100*fFwdSysErr), "f");
+      d0->SetLineColor(SYSERR_COLOR);
+      d0->SetMarkerColor(SYSERR_COLOR);
+      d0->SetFillColor(SYSERR_COLOR);
+      d0->SetFillStyle(SYSERR_STYLE);
+      d0->SetMarkerStyle(0);
+      d0->SetLineWidth(0);
+      i++;
+    }
+    if (!fCentAxis && i % 2 == 1)  {
+      // To make sure the 'data' and 'mirrored' entries are on a line
+      // by themselves 
+      TLegendEntry* dd = l->AddEntry("dd", "   ", "");
+      dd->SetTextSize(0);
+      dd->SetFillStyle(0);
+      dd->SetFillColor(0);
+      dd->SetLineWidth(0);
+      dd->SetLineColor(0);
+      dd->SetMarkerSize(0);
+    }
+    // Add entry for 'data'
     TLegendEntry* d1 = l->AddEntry("d1", "Data", "lp");
     d1->SetLineColor(kBlack);
     d1->SetMarkerColor(kBlack);
     d1->SetMarkerStyle(20);
+    // Add entry for 'mirrored data'
     TLegendEntry* d2 = l->AddEntry("d2", "Mirrored data", "lp");
     d2->SetLineColor(kBlack);
     d2->SetMarkerColor(kBlack);
@@ -651,6 +819,38 @@ struct dNdetaDrawer
     l->Draw();
   }
   //__________________________________________________________________
+  /** 
+   * Build centrality legend 
+   * 
+   * @param axis    Stack to include 
+   * @param x1      Lower X coordinate in the range [0,1]
+   * @param y1      Lower Y coordinate in the range [0,1]
+   * @param x2      Upper X coordinate in the range [0,1]
+   * @param y2             Upper Y coordinate in the range [0,1]
+   */
+  void BuildCentLegend(Double_t x1, Double_t y1, Double_t x2, Double_t y2)
+  {
+    if (!fCentAxis) return;
+
+    TLegend* l = new TLegend(x1,y1,x2,y2);
+    l->SetNColumns(1);
+    l->SetFillColor(0);
+    l->SetFillStyle(0);
+    l->SetBorderSize(0);
+    l->SetTextFont(132);
+
+    Int_t n = fCentAxis->GetNbins();
+    for (Int_t i = 1; i <= n; i++) { 
+      Double_t low = fCentAxis->GetBinLowEdge(i);
+      Double_t upp = fCentAxis->GetBinUpEdge(i);
+      TLegendEntry* e = l->AddEntry(Form("dummy%02d", i),
+                                   Form("%3d%% - %3d%%", 
+                                        int(low), int(upp)), "pl");
+      e->SetMarkerColor(GetCentralityColor(i));
+    }
+    l->Draw();
+  }
+  //__________________________________________________________________
   /** 
    * Plot the results
    *    
@@ -665,18 +865,20 @@ struct dNdetaDrawer
     p1->SetBorderSize(0);
     p1->SetBorderMode(0);
     p1->SetBottomMargin(yd > 0.001 ? 0.001 : 0.1);
-    p1->SetRightMargin(0.05);
-    p1->SetGridx();
+    p1->SetRightMargin(0.03);
+    if (fShowLeftRight || fShowRatios) p1->SetGridx();
     p1->SetTicks(1,1);
     p1->SetNumber(1);
     p1->Draw();
     p1->cd();
-    
+
     // Info("PlotResults", "Plotting results with max=%f", max);
     fResults->SetMaximum(1.15*max);
     fResults->SetMinimum(yd > 0.00001 ? -0.1 : 0);
 
-    FixAxis(fResults, 1/(1-yd)/1.7, "#frac{1}{N} #frac{dN_{ch}}{d#eta}");
+    FixAxis(fResults, (1-yd)*(yd > .001 ? 1 : .9 / 1.2), 
+           "#font[12]{#frac{1}{N} "
+           "#frac{dN_{#font[132]{ch}}}{d#font[152]{#eta}}}");
 
     p1->Clear();
     fResults->DrawClone("nostack e1");
@@ -685,7 +887,7 @@ struct dNdetaDrawer
     fRangeParam->fSlave1Pad  = p1;
 
     // Draw other data
-    if (fShowOthers || fShowAlice) {
+    if (fShowOthers != 0) {
       TGraphAsymmErrors* o      = 0;
       TIter              nextG(fOthers->GetListOfGraphs());
       while ((o = static_cast<TGraphAsymmErrors*>(nextG())))
@@ -693,15 +895,14 @@ struct dNdetaDrawer
     }
 
     // Make a legend in the main result pad
-    BuildLegend(fResults, fOthers, .15,p1->GetBottomMargin()+.01,.90,.35);
-#if 0
-    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);
-#endif
+    BuildCentLegend(.12, 1-p1->GetTopMargin()-.01-.5,  
+                   .35, 1-p1->GetTopMargin()-.01-.1);
+    Double_t x1 = (fCentAxis ? .7 : .15); 
+    Double_t x2 = (fCentAxis ? 1-p1->GetRightMargin()-.01: .90);
+    Double_t y1 = (fCentAxis ? .5: p1->GetBottomMargin()+.01); 
+    Double_t y2 = (fCentAxis ? 1-p1->GetTopMargin()-.01-.15 : .35);
+                  
+    BuildLegend(fResults, fOthers, x1, y1, x2, y2);
 
     // Put a title on top
     fTitle.ReplaceAll("@", " ");
@@ -711,16 +912,21 @@ struct dNdetaDrawer
     tit->SetTextSize(0.05);
     tit->Draw();
 
+    Int_t aliceBlue = TColor::GetColor(41,73,156);
     // Put a nice label in the plot
     TString     eS;
     UShort_t    snn = fSNNString->GetUniqueID();
     const char* sys = fSysString->GetTitle();
+    if (snn == 2750) snn = 2760;
     if (snn > 1000) eS = Form("%4.2fTeV", float(snn)/1000);
     else            eS = Form("%3dGeV", snn);
-    TLatex* tt = new TLatex(.93, .93, Form("%s #sqrt{s}=%s, %s", 
+    TLatex* tt = new TLatex(.93, .93, Form("%s #sqrt{s%s}=%s, %s", 
                                           sys, 
+                                          (fCentAxis ? "_{NN}" : ""),
                                           eS.Data(), 
+                                          fCentAxis ? "by centrality" : 
                                           fTrigString->GetTitle()));
+    tt->SetTextColor(aliceBlue);
     tt->SetNDC();
     tt->SetTextFont(132);
     tt->SetTextAlign(33);
@@ -730,6 +936,7 @@ struct dNdetaDrawer
     Int_t nev = 0;
     if (fTriggers) nev = fTriggers->GetBinContent(1);
     TLatex* et = new TLatex(.93, .83, Form("%d events", nev));
+    et->SetTextColor(aliceBlue);
     et->SetNDC();
     et->SetTextFont(132);
     et->SetTextAlign(33);
@@ -741,6 +948,7 @@ struct dNdetaDrawer
       vt->SetNDC();
       vt->SetTextFont(132);
       vt->SetTextAlign(33);
+      vt->SetTextColor(aliceBlue);
       vt->Draw();
     }
     // results->Draw("nostack e1 same");
@@ -750,22 +958,29 @@ struct dNdetaDrawer
 
 
     // Mark the plot as preliminary
-    TLatex* pt = new TLatex(.12, .93, "Preliminary");
+    TLatex* pt = new TLatex(.12, .93, "Work in progress");
     pt->SetNDC();
     pt->SetTextFont(22);
-    pt->SetTextSize(0.07);
-    pt->SetTextColor(kRed+1);
+    // pt->SetTextSize();
+    pt->SetTextColor(TColor::GetColor(234,26,46));
     pt->SetTextAlign(13);
     pt->Draw();
 
-    if (!gSystem->AccessPathName("ALICE.png")) { 
-      TPad* logo = new TPad("logo", "logo", .12, .65, .25, .85, 0, 0, 0);
-      logo->SetFillStyle(0);
-      logo->Draw();
-      logo->cd();
+    const char*  logos[] = { "ALICE.png", "FMD.png", 0 };
+    const char** logo    = logos;
+    while (*logo) {
+      if (gSystem->AccessPathName(*logo)) {
+       logo++;
+       continue;
+      }
+      TPad* pad = new TPad("logo", "logo", .12, .7, .25, .9, 0, 0, 0);
+      pad->SetFillStyle(0);
+      pad->Draw();
+      pad->cd();
       TImage* i = TImage::Create();
-      i->ReadImage("ALICE.png");
+      i->ReadImage(*logo);
       i->Draw();
+      break;
     }
     p1->cd();
   }
@@ -779,14 +994,14 @@ struct dNdetaDrawer
    */
   void PlotRatios(Double_t max, Double_t y1, Double_t y2) 
   {
-    if (!fShowRatios) return;
+    if (fShowRatios == 0) return;
 
     bool isBottom = (y1 < 0.0001);
     Double_t yd = y2 - y1;
     // Make a sub-pad for the result itself
     TPad* p2 = new TPad("p2", "p2", 0, y1, 1.0, y2, 0, 0, 0);
     p2->SetTopMargin(0.001);
-    p2->SetRightMargin(0.05);
+    p2->SetRightMargin(0.03);
     p2->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001);
     p2->SetGridx();
     p2->SetTicks(1,1);
@@ -795,7 +1010,7 @@ struct dNdetaDrawer
     p2->cd();
 
     // Fix up axis
-    FixAxis(fRatios, 1/yd/1.7, "Ratios", 7);
+    FixAxis(fRatios, yd, "Ratios", 7);
 
     fRatios->SetMaximum(1+TMath::Max(.22,1.05*max));
     fRatios->SetMinimum(1-TMath::Max(.32,1.05*max));
@@ -858,7 +1073,7 @@ struct dNdetaDrawer
     // Make a sub-pad for the result itself
     TPad* p3 = new TPad("p3", "p3", 0, y1, 1.0, y2, 0, 0, 0);
     p3->SetTopMargin(0.001);
-    p3->SetRightMargin(0.05);
+    p3->SetRightMargin(0.03);
     p3->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001);
     p3->SetGridx();
     p3->SetTicks(1,1);
@@ -867,7 +1082,7 @@ struct dNdetaDrawer
     p3->cd();
 
     // Fix up axis
-    FixAxis(fLeftRight, 1/yd/1.7, "Right/Left", 4);
+    FixAxis(fLeftRight, yd, "Right/Left", 4);
 
     fLeftRight->SetMaximum(1+TMath::Max(.12,1.05*max));
     fLeftRight->SetMinimum(1-TMath::Max(.15,1.05*max));
@@ -876,12 +1091,17 @@ struct dNdetaDrawer
 
     
     // Make a legend
-    TLegend* l2 = p3->BuildLegend(.15,p3->GetBottomMargin()+.01,.9,.5);
-    l2->SetNColumns(2);
-    l2->SetFillColor(0);
-    l2->SetFillStyle(0);
-    l2->SetBorderSize(0);
-    l2->SetTextFont(132);
+    Double_t xx1 = (fCentAxis ? .7                           : .15); 
+    Double_t xx2 = (fCentAxis ? 1-p3->GetRightMargin()-.01   : .90);
+    Double_t yy1 = p3->GetBottomMargin()+.01;
+    Double_t yy2 = (fCentAxis ? 1-p3->GetTopMargin()-.01-.15 : .5);
+    BuildLegend(fLeftRight, 0, xx1, yy1, xx2, yy2);
+    // TLegend* l2 = p3->BuildLegend(.15,p3->GetBottomMargin()+.01,.9,.5);
+    // l2->SetNColumns(2);
+    // l2->SetFillColor(0);
+    // l2->SetFillStyle(0);
+    // l2->SetBorderSize(0);
+    // l2->SetTextFont(132);
 
     // Make a nice band from 0.9 to 1.1
     TGraphErrors* band = new TGraphErrors(2);
@@ -926,10 +1146,12 @@ struct dNdetaDrawer
     if (!list) return 0;
     
     TH1* ret = static_cast<TH1*>(list->FindObject(name));
+#if 0
     if (!ret) {
       // all->ls();
       Warning("GetResult", "Histogram %s not found", name);
     }
+#endif
 
     return ret;
   }
@@ -951,10 +1173,7 @@ struct dNdetaDrawer
     // Rebin if needed 
     Rebin(hist);
 
-    // Info("AddHistogram", "Adding %s to %s", 
-    //      hist->GetName(), stack->GetName());
     stack->Add(hist, option);
-    // stack->ls();
     return hist->GetMaximum();
   }
   //__________________________________________________________________
@@ -982,8 +1201,6 @@ struct dNdetaDrawer
     sym = Symmetrice(hist);
     stack->Add(sym, option);
 
-    // Info("AddHistogram", "Adding %s and %s to %s", 
-    //      hist->GetName(), sym->GetName(), stack->GetName());
     return hist->GetMaximum();
   }
 
@@ -1227,32 +1444,36 @@ struct dNdetaDrawer
    */
   TH1* Ratio(const TObject* o1, const TObject* o2, Double_t& max) const
   {
-    TH1* r = 0;
+    if (!o1 || !o2) return 0;
+
+    TH1*        r  = 0;
+    const TAttMarker* m1 = 0;
+    const TAttMarker* m2 = 0;
     const TH1* h1 = dynamic_cast<const TH1*>(o1); 
     if (h1) { 
-      // Info("Ratio", "First is a TH1");
+      m1 = h1;
       const TH1* h2 = dynamic_cast<const TH1*>(o2); 
       if (h2) { 
-       // Info("Ratio", "Second is a TH1");
-       r = RatioHH(h1,h2,max);
+       m2 = h2;
+       r = RatioHH(h1,h2);
       }
       else {
        const TGraph* g2 = dynamic_cast<const TGraph*>(o2);
        if (g2) { 
-         // Info("Ratio", "Second os a TGraph");
-         r = RatioHG(h1,g2,max);      
+         m2 = g2;
+         r = RatioHG(h1,g2);      
        }
       }
     }
     else {
       const TGraphAsymmErrors* g1 = dynamic_cast<const TGraphAsymmErrors*>(o1);
       if (g1) { 
-       // Info("Ratio", "First is a TGraphAsymmErrors");
+       m1 = g1;
        const TGraphAsymmErrors* g2 = 
          dynamic_cast<const TGraphAsymmErrors*>(o2);
        if (g2) {
-         // Info("Ratio", "Second is a TGraphAsymmErrors");
-         r = RatioGG(g1, g2, max);
+         m2 = g2;
+         r = RatioGG(g1, g2);
        }
       }
     }
@@ -1266,9 +1487,16 @@ struct dNdetaDrawer
       delete r; 
       r = 0; 
     }
-    // for (Int_t bin = 1; bin <= r->GetNbinsX(); bin++) 
-    //   if (r->GetBinContent(bin) != 0) return r;
-      
+    if (r) {
+      r->SetMarkerStyle(m2->GetMarkerStyle());
+      r->SetMarkerColor(m1->GetMarkerColor());
+      r->SetMarkerSize(0.9*m1->GetMarkerSize());
+      r->SetName(Form("%s_over_%s", o1->GetName(), o2->GetName()));
+      r->SetTitle(Form("%s / %s", o1->GetTitle(), o2->GetTitle()));
+      r->SetDirectory(0);
+      max = TMath::Max(RatioMax(r), max);
+    }
+
     return r;
   }
   //__________________________________________________________________
@@ -1278,25 +1506,15 @@ struct dNdetaDrawer
    * 
    * @param h  Numerator 
    * @param g  Divisor 
-   * @param max Maximum diviation from 1 
    * 
    * @return h/g 
    */
-  TH1* RatioHG(const TH1* h, const TGraph* g, Double_t& max) const 
+  TH1* RatioHG(const TH1* h, const TGraph* g) const 
   {
     if (!h || !g) return 0;
 
     TH1* ret = static_cast<TH1*>(h->Clone("tmp"));
-    ret->SetName(Form("%s_over_%s", h->GetName(), g->GetName()));
-    ret->SetTitle(Form("%s / %s", h->GetTitle(), g->GetTitle()));
     ret->Reset();
-    ret->SetMarkerStyle(h->GetMarkerStyle());
-    ret->SetMarkerColor(g->GetMarkerColor());
-    ret->SetMarkerSize(0.9*h->GetMarkerSize());
-    // ret->SetMarkerStyle(g->GetMarkerStyle());
-    // ret->SetMarkerColor(h->GetMarkerColor());
-    // ret->SetMarkerSize(0.9*g->GetMarkerSize());
-    ret->SetDirectory(0);
     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; }
@@ -1314,8 +1532,6 @@ struct dNdetaDrawer
       ret->SetBinContent(i, c / f);
       ret->SetBinError(i, h->GetBinError(i) / f);
     }
-    if (ret->GetEntries() > 0) 
-      max = TMath::Max(RatioMax(ret), max);
     return ret;
   }
   //__________________________________________________________________
@@ -1324,25 +1540,14 @@ struct dNdetaDrawer
    * 
    * @param h1 First histogram (numerator) 
    * @param h2 Second histogram (denominator)
-   * @param max Maximum diviation from 1 
    * 
    * @return h1 / h2
    */
-  TH1* RatioHH(const TH1* h1, const TH1* h2, Double_t& max) const
+  TH1* RatioHH(const TH1* h1, const TH1* h2) const
   {
     if (!h1 || !h2) return 0;
-    TH1* t1 = static_cast<TH1*>(h1->Clone(Form("%s_%s", 
-                                              h1->GetName(), 
-                                              h2->GetName())));
-    t1->SetTitle(Form("%s / %s", h1->GetTitle(), h2->GetTitle()));
+    TH1* t1 = static_cast<TH1*>(h1->Clone("tmp"));
     t1->Divide(h2);
-    // t1->SetMarkerColor(h1->GetMarkerColor());
-    // t1->SetMarkerStyle(h2->GetMarkerStyle());
-    t1->SetMarkerColor(h2->GetMarkerColor());
-    t1->SetMarkerStyle(h1->GetMarkerStyle());
-    t1->SetMarkerSize(0.9*h1->GetMarkerSize());
-    t1->SetDirectory(0);
-    max = TMath::Max(RatioMax(t1), max);
     return t1;
   }
   //__________________________________________________________________
@@ -1351,12 +1556,11 @@ struct dNdetaDrawer
    * 
    * @param g1 Numerator 
    * @param g2 Denominator
-   * @param max Maximum diviation from 1 
    * 
    * @return g1 / g2 in a histogram 
    */
   TH1* RatioGG(const TGraphAsymmErrors* g1, 
-              const TGraphAsymmErrors* g2, Double_t& max) const
+              const TGraphAsymmErrors* g2) const
   {
     Int_t    nBins = g1->GetN();
     TArrayF  bins(nBins+1);
@@ -1371,22 +1575,13 @@ struct dNdetaDrawer
       if (i == 0) dx  = dxi;
       else if (dxi != dx) dx = 0;
     }
-    TString name(Form("%s_%s", g1->GetName(), g2->GetName()));
-    TString title(Form("%s / %s", g1->GetTitle(), g2->GetTitle()));
     TH1* h = 0;
     if (dx != 0) {
-      h = new TH1F(name.Data(), title.Data(), nBins, bins[0], bins[nBins]);
+      h = new TH1F("tmp", "tmp", nBins, bins[0], bins[nBins]);
     }
     else {
-      h = new TH1F(name.Data(), title.Data(), nBins, bins.fArray);
+      h = new TH1F("tmp", "tmp", nBins, bins.fArray);
     }
-    h->SetMarkerStyle(g1->GetMarkerStyle());
-    h->SetMarkerColor(g2->GetMarkerColor());
-    h->SetMarkerSize(0.9*g1->GetMarkerSize());
-    // h->SetMarkerStyle(g2->GetMarkerStyle());
-    // h->SetMarkerColor(g1->GetMarkerColor());
-    // h->SetMarkerSize(0.9*g2->GetMarkerSize());
-    h->SetDirectory(0);
 
     Double_t low  = g2->GetX()[0];
     Double_t high = g2->GetX()[g2->GetN()-1];
@@ -1401,7 +1596,6 @@ struct dNdetaDrawer
       h->SetBinContent(i+1, c1 / c2);
       h->SetBinError(i+1, e1 / c2);
     }
-    max = TMath::Max(RatioMax(h), max);
     return h;
   }  
   /* @} */
@@ -1448,7 +1642,7 @@ struct dNdetaDrawer
    * @param force  Whether to draw the stack first or not
    * @param ynDiv  Divisions on Y axis
    */
-  void FixAxis(THStack* stack, Double_t s, const char* ytitle,
+  void FixAxis(THStack* stack, Double_t yd, const char* ytitle,
                Int_t ynDiv=210, Bool_t force=true)
   {
     if (!stack) { 
@@ -1462,17 +1656,25 @@ struct dNdetaDrawer
       Warning("FixAxis", "Stack %s has no histogram", stack->GetName());
       return;
     }
-    
-    h->SetXTitle("#eta");
+    Double_t s = 1/yd/1.2;
+    // Info("FixAxis", "for %s, s=1/%f=%f", stack->GetName(), yd, s);
+
+    h->SetXTitle("#font[152]{#eta}");
     h->SetYTitle(ytitle);
     TAxis* xa = h->GetXaxis();
     TAxis* ya = h->GetYaxis();
+
+    // Int_t   npixels = 20
+    // Float_t dy = gPad->AbsPixeltoY(0) - gPad->AbsPixeltoY(npixels);
+    // Float_t ts = dy/(gPad->GetY2() - gPad->GetY1());
+
     if (xa) {
-      xa->SetTitle("#eta");
+      // xa->SetTitle(h->GetXTitle());
       // xa->SetTicks("+-");
       xa->SetTitleSize(s*xa->GetTitleSize());
       xa->SetLabelSize(s*xa->GetLabelSize());
       xa->SetTickLength(s*xa->GetTickLength());
+      // xa->SetTitleOffset(xa->GetTitleOffset()/s);
 
       if (stack != fResults) {
        TAxis* rxa = fResults->GetXaxis();
@@ -1480,7 +1682,7 @@ struct dNdetaDrawer
       }
     }
     if (ya) {
-      ya->SetTitle(ytitle);
+      // ya->SetTitle(h->GetYTitle());
       ya->SetDecimals();
       // ya->SetTicks("+-");
       ya->SetNdivisions(ynDiv);
@@ -1489,32 +1691,279 @@ struct dNdetaDrawer
       ya->SetLabelSize(s*ya->GetLabelSize());
     }
   }
+  //__________________________________________________________________
+  /** 
+   * Merge two histograms into one 
+   * 
+   * @param cen    Central part
+   * @param fwd    Forward part
+   * @param xlow   On return, lower eta bound
+   * @param xhigh  On return, upper eta bound
+   * 
+   * @return Newly allocated histogram or null
+   */
+  TH1* 
+  Merge(const TH1* cen, const TH1* fwd, Double_t& xlow, Double_t& xhigh)
+  {
+    TH1* tmp = static_cast<TH1*>(fwd->Clone("tmp"));
+    TString name(fwd->GetName());
+    name.ReplaceAll("Forward", "Merged");
+    tmp->SetName(name);
+
+    // tmp->SetMarkerStyle(28);
+    // tmp->SetMarkerColor(kBlack);
+    tmp->SetDirectory(0);
+    xlow  = 100;
+    xhigh = -100;
+    for (Int_t i = 1; i <= tmp->GetNbinsX(); i++) {
+      Double_t cc = cen->GetBinContent(i);
+      Double_t cf = fwd->GetBinContent(i);
+      Double_t ec = cen->GetBinError(i);
+      Double_t ef = fwd->GetBinError(i);
+      Double_t nc = cf;
+      Double_t ne = ef;
+      if (cc < 0.001 && cf < 0.01) continue;
+      xlow  = TMath::Min(tmp->GetXaxis()->GetBinLowEdge(i),xlow);
+      xhigh = TMath::Max(tmp->GetXaxis()->GetBinUpEdge(i),xhigh);
+      if (cc > 0.001) {
+       nc = cc;
+       ne = ec;
+       if (cf > 0.001) {
+         nc  = (cf + cc) / 2;
+         ne  = TMath::Sqrt(ec*ec + ef*ef);
+       }
+      }
+      tmp->SetBinContent(i, nc);
+      tmp->SetBinError(i, ne);
+    }
+    return tmp;
+  }
+  //____________________________________________________________________
+  /** 
+   * Fit  @f$g(x;A_1,A_2,\sigma_1,\sigma_2)@f$ to histogram data 
+   * 
+   * @param tmp    Histogram
+   * @param xlow   Lower x bound
+   * @param xhigh  Upper x bound 
+   *
+   * @return Fitted function 
+   */
+  TF1* 
+  FitMerged(TH1* tmp, Double_t xlow, Double_t xhigh)
+  {
+    TF1* tmpf  = new TF1("tmpf", "gaus", xlow, xhigh);
+    tmp->Fit(tmpf, "NQ", "");
+    tmp->SetDirectory(0);
+
+    TF1* fit = new TF1("f", myFunc, xlow, xhigh, 4);
+    fit->SetParNames("a_{1}", "a_{2}", "#sigma_{1}", "#sigma_{2}");
+    fit->SetParameters(tmpf->GetParameter(0), 
+                      .2, 
+                      tmpf->GetParameter(2), 
+                      tmpf->GetParameter(2)/4);
+    fit->SetParLimits(3, 0, 100);
+    fit->SetParLimits(4, 0, 100);
+    tmp->Fit(fit,"0W","");
+
+    delete tmpf;
+    return fit;
+  }
+  //____________________________________________________________________
+  /** 
+   * Make band of systematic errors 
+   * 
+   * @param tmp Histogram
+   * @param fit Fit 
+   */
+  void
+  MakeSysError(TH1* tmp, TH1* cen, TH1* fwd, TF1* fit)
+  {
+    for (Int_t i = 1; i <= tmp->GetNbinsX(); i++) {
+      Double_t tc = tmp->GetBinContent(i);
+      if (tc < 0.01) continue;
+      Double_t fc = fwd->GetBinContent(i);
+      Double_t cc = cen->GetBinContent(i);
+      Double_t sysErr = fFwdSysErr;
+      if (cc > .01 && fc > 0.01) 
+       sysErr = (fFwdSysErr+fCenSysErr) / 2;
+      else if (cc > .01) 
+       sysErr = fCenSysErr;
+      Double_t x = tmp->GetXaxis()->GetBinCenter(i);
+      Double_t y = fit->Eval(x);
+      tmp->SetBinContent(i, y);
+      tmp->SetBinError(i,sysErr*y);
+    }
+    TString name(tmp->GetName());
+    name.ReplaceAll("Merged", "SysError");
+    tmp->SetName(name);
+    tmp->SetMarkerColor(SYSERR_COLOR);
+    tmp->SetLineColor(SYSERR_COLOR);
+    tmp->SetFillColor(SYSERR_COLOR);
+    tmp->SetFillStyle(SYSERR_STYLE);
+    tmp->SetMarkerStyle(0);
+    tmp->SetLineWidth(0);
+  }
+  void CorrectForward(TH1* h) const
+  {
+    if (!fRemoveOuters) return;
+    
+    for (Int_t i = 1; i <= h->GetNbinsX(); i++) { 
+      Double_t eta = h->GetBinCenter(i);
+      if (TMath::Abs(eta) < 2.3) { 
+       h->SetBinContent(i, 0);
+       h->SetBinError(i, 0);
+      }
+    }
+  }
+  void CorrectCentral(TH1* h) const 
+  {
+    if (fClusterScale.IsNull()) return;
+    TString t(h->GetTitle());
+    Info("CorrectCentral", "Replacing Central with Tracklets in %s", t.Data());
+    t.ReplaceAll("Central", "Tracklets");
+    h->SetTitle(t);
+
+    TF1* cf = new TF1("clusterScale", fClusterScale);
+    for (Int_t i = 1; i <= h->GetNbinsX(); i++) { 
+      Double_t eta = h->GetBinCenter(i);
+      Double_t f   = cf->Eval(eta);
+      Double_t c   = h->GetBinContent(i);
+      if (f < .1) f = 1;
+      h->SetBinContent(i, c / f);
+    }
+    delete cf;
+  }
+  //____________________________________________________________________
+  void Export(const char* basename)
+  {
+    TString bname(basename);
+    bname.ReplaceAll(" ", "_");
+    bname.ReplaceAll("-", "_");
+    TString fname(Form("%s.C", bname.Data()));
+
+    std::ofstream out(fname.Data());
+    if (!out) { 
+      Error("Export", "Failed to open output file %s", fname.Data());
+      return;
+    }
+    out << "// Create by dNdetaDrawer\n"
+       << "void " << bname << "(THStack* stack, TLegend* l, Int_t m)\n"
+       << "{"
+       << "   Int_t ma[] = { 24, 25, 26, 32,\n"
+       << "                  20, 21, 22, 33,\n"
+       << "                  34, 30, 29, 0, \n"
+       << "                  23, 27 };\n"
+       << "   Int_t mm = ((m < 20 || m > 34) ? 0 : ma[m-20]);\n\n";
+    TList* hists = fResults->GetHists();
+    TIter  next(hists);
+    TH1*   hist = 0;
+    while ((hist = static_cast<TH1*>(next()))) { 
+      TString hname = hist->GetName();
+      hname.Append(Form("_%04x", (gRandom->Integer(0xffff) & 0xffff)));
+      hist->SetName(hname);
+      hist->GetListOfFunctions()->Clear();
+      hist->SavePrimitive(out, "nodraw");
+      bool mirror = hname.Contains("mirror");
+      bool syserr = hname.Contains("SysError");
+      if (!syserr) 
+       out << "   " << hname << "->SetMarkerStyle(" 
+           << (mirror ? "mm" : "m") << ");\n";
+      else 
+       out << "   " << hname << "->SetMarkerStyle(1);\n";
+      out << "   stack->Add(" << hname 
+         << (syserr ? ",\"e5\"" : "") << ");\n\n";
+    }
+    UShort_t    snn = fSNNString->GetUniqueID();
+    // const char* sys = fSysString->GetTitle();
+    TString eS;
+    if      (snn == 2750)     snn = 2760;
+    if      (snn < 1000)      eS = Form("%3dGeV", snn);
+    else if (snn % 1000 == 0) eS = Form("%dTeV", snn/1000);
+    else                      eS = Form("%4.2fTeV", float(snn)/1000);
+    out << "  if (l) {\n"
+       << "    TLegendEntry* e = l->AddEntry(\"\",\"" << eS << "\",\"pl\");\n"
+       << "    e->SetMarkerStyle(m);\n"
+       << "    e->SetMarkerColor(kBlack);\n"
+       << "  }\n"
+       << "}\n" << std::endl;
+  }
+             
   /* @} */
 
 
 
   //__________________________________________________________________
-  Bool_t       fShowOthers;   // Show other data
-  Bool_t       fShowAlice;    // Show ALICE published data
+  /** 
+   * @{ 
+   * @name Options 
+   */
   Bool_t       fShowRatios;   // Show ratios 
   Bool_t       fShowLeftRight;// Show asymmetry 
-  UShort_t     fRebin;        // Rebinning factor 
+  Bool_t       fShowRings;    // Show rings too
+  Bool_t       fExport;       // Export results to file
   Bool_t       fCutEdges;     // Whether to cut edges
+  Bool_t       fRemoveOuters; // Whether to remove outers
+  UShort_t     fShowOthers;   // Show other data
+  /* @} */
+  /** 
+   * @{ 
+   * @name Settings 
+   */
+  UShort_t     fRebin;        // Rebinning factor 
+  Double_t     fFwdSysErr;    // Systematic error in forward range
+  Double_t     fCenSysErr;    // Systematic error in central range 
   TString      fTitle;        // Title on plot
+  TString      fClusterScale; // Scaling of clusters to tracklets      
+  /* @} */
+  /** 
+   * @{ 
+   * @name Read (or set) information 
+   */
   TNamed*      fTrigString;   // Trigger string (read, or set)
+  TNamed*      fNormString;   // Normalisation string (read, or set)
   TNamed*      fSNNString;    // Energy string (read, or set)
   TNamed*      fSysString;    // Collision system string (read or set)
   TAxis*       fVtxAxis;      // Vertex cuts (read or set)
   TAxis*       fCentAxis;     // Centrality axis
+  /* @} */
+  /** 
+   * @{ 
+   * @name Resulting plots 
+   */
   THStack*     fResults;      // Stack of results 
   THStack*     fRatios;       // Stack of ratios 
   THStack*     fLeftRight;    // Left-right asymmetry
   TMultiGraph* fOthers;       // Older data 
   TH1*         fTriggers;     // Number of triggers
+  TH1*         fTruth;        // Pointer to truth 
+  /* @} */
   RangeParam*  fRangeParam;   // Parameter object for range zoom 
-  
 };
 
+//____________________________________________________________________
+/** 
+ * Function to calculate 
+ * @f[
+ *  g(x;A_1,A_2,\sigma_1,\sigma_2) = 
+ *       A_1\left(\frac{1}{2\pi\sigma_1}e^{(x/\sigma_1)^2} - 
+ *           A_2\frac{1}{2\pi\sigma_2}e^{(x/\sigma_2)^2}\right)
+ * @f]
+ * 
+ * @param xp Pointer to x array
+ * @param pp Pointer to parameter array 
+ * 
+ * @return @f$g(x;A_1,A_2,\sigma_1,\sigma_2)@f$
+ */
+Double_t myFunc(Double_t* xp, Double_t* pp)
+{
+  Double_t x  = xp[0];
+  Double_t a1 = pp[0];
+  Double_t a2 = pp[1];
+  Double_t s1 = pp[2];
+  Double_t s2 = pp[3];
+  return a1*(TMath::Gaus(x, 0, s1) - a2 * TMath::Gaus(x, 0, s2));
+}
+
 //=== Stuff for auto zooming =========================================
 void UpdateRange(dNdetaDrawer::RangeParam* p)
 {
@@ -1530,7 +1979,6 @@ void UpdateRange(dNdetaDrawer::RangeParam* p)
   Int_t    last  = p->fMasterAxis->GetLast();
   Double_t x1    = p->fMasterAxis->GetBinCenter(first);
   Double_t x2    = p->fMasterAxis->GetBinCenter(last);
-  //Info("UpdateRange", "Range set to [%3d,%3d]->[%f,%f]", first, last, x1,x2);
 
   if (p->fSlave1Axis) { 
     Int_t i1 = p->fSlave1Axis->FindBin(x1);
@@ -1571,7 +2019,34 @@ void RangeExec(dNdetaDrawer::RangeParam* p)
   UpdateRange(p);
 }
 
-//=== Steering function ==============================================  
+//=== Steering functions
+//==============================================  
+void
+Usage()
+{
+  Info("DrawdNdeta", "Usage: DrawdNdeta(FILE,TITLE,REBIN,OTHERS,FLAGS)\n\n"
+       "  const char* FILE   File name to open (\"forward_root\")\n"
+       "  const char* TITLE  Title to put on plot (\"\")\n"
+       "  UShort_t    REBIN  Rebinning factor (1)\n"
+       "  UShort_t    OTHERS Other data to draw - more below (0x7)\n"
+       "  UShort_t    FLAGS  Visualisation flags - more below (0x7)\n\n"
+       " OTHERS is a bit mask of\n\n"
+       "  0x1   Show UA5 data (INEL,NSD, ppbar, 900GeV)\n"
+       "  0x2   Show CMS data (NSD, pp)\n"
+       "  0x4   Show published ALICE data (INEL,INEL>0,NSD, pp)\n"
+       "  0x8   Show event genertor data\n\n"
+       " FLAGS is a bit mask of\n\n"
+       "  0x1   Show ratios of data to other data and possibly MC\n"
+       "  0x2   Show left-right asymmetry\n"
+       "  0x4   Show systematic error band\n"
+       "  0x8   Show individual ring results (INEL only)\n"
+       "  0x10  Cut edges when rebinning\n"
+       "  0x20  Remove FMDxO points\n"
+       "  0x40  Do not make our own canvas\n"
+       );
+}
+
+//____________________________________________________________________
 /** 
  * Draw @f$ dN/d\eta@f$ 
  * 
@@ -1590,25 +2065,36 @@ void RangeExec(dNdetaDrawer::RangeParam* p)
  */
 void
 DrawdNdeta(const char* filename="forward_dndeta.root", 
-          Int_t       flags=0xf,
           const char* title="",
           UShort_t    rebin=5, 
-          Bool_t      cutEdges=false,
+          UShort_t    others=0x7,
+          UShort_t    flags=0x7,
           UShort_t    sNN=0, 
           UShort_t    sys=0,
           UShort_t    trg=0,
           Float_t     vzMin=999, 
           Float_t     vzMax=-999)
 {
+  TString fname(filename);
+  fname.ToLower();
+  if (fname.CompareTo("help") == 0 || 
+      fname.CompareTo("--help") == 0) { 
+    Usage();
+    return;
+  }
   dNdetaDrawer* pd = new dNdetaDrawer;
   dNdetaDrawer& d = *pd;
   d.SetRebin(rebin);
-  d.SetCutEdges(cutEdges);
   d.SetTitle(title);
-  d.SetShowOthers(flags & 0x1);
-  d.SetShowAlice(flags & 0x2);
-  d.SetShowRatios(flags & 0x4);
-  d.SetShowLeftRight(flags & 0x8);
+  d.SetShowOthers(others);
+  d.SetShowRatios(flags & 0x1);
+  d.SetShowLeftRight(flags & 0x2);
+  d.SetForwardSysError(flags & 0x4 ? 0.076 : 0);
+  d.SetShowRings(flags & 0x8);
+  d.SetCutEdges(flags & 0x10);
+  d.fRemoveOuters = (flags & 0x20);
+  d.SetExport(flags & 0x40);
+  // d.fClusterScale = "1.06 -0.003*x +0.0119*x*x";
   // Do the below if your input data does not contain these settings 
   if (sNN > 0) d.SetSNN(sNN);     // Collision energy per nucleon pair (GeV)
   if (sys > 0) d.SetSys(sys);     // Collision system (1:pp, 2:PbPB)