]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG2/FORWARD/analysis2/DrawdNdeta.C
Added script to test the energy loss assumptions and
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / DrawdNdeta.C
index 33086534ac670b46bb2be934354eb3bc298d6ff0..7883d71739605ddd1dc30f5c0077ef1035323835 100644 (file)
@@ -14,6 +14,7 @@
  * @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 
@@ -58,22 +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*
-      fNormString(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;
@@ -82,6 +102,9 @@ struct dNdetaDrawer
     fRangeParam->fSlave2Axis = 0;
     fRangeParam->fSlave2Pad  = 0;
   }
+  dNdetaDrawer(const dNdetaDrawer&) {}
+  dNdetaDrawer& operator=(const dNdetaDrawer&) { return *this; }
+
   //__________________________________________________________________
   virtual ~dNdetaDrawer()
   {
@@ -110,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,
@@ -135,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 
    * 
@@ -155,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; }
   /* @} */
   //==================================================================  
   /** 
@@ -248,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");
@@ -255,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"));
@@ -285,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) { 
@@ -299,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();
@@ -369,26 +444,18 @@ struct dNdetaDrawer
   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) { 
-#if 0
-      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");
-#endif
-      return 0;
-    }
-    thisOther = reinterpret_cast<TMultiGraph*>(ret);
-    
+                                            centLow,centHigh,
+                                            fShowOthers));
+    if (!ret) return 0;
+
+    thisOther = reinterpret_cast<TMultiGraph*>(ret);    
     return thisOther;
   }
   //__________________________________________________________________
@@ -401,23 +468,29 @@ 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), -1000, 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;
     }
     
+    TObjArray* a = new TObjArray;
     for (UShort_t i = 0; i < n; i++) { 
       UShort_t centLow  = fCentAxis->GetBinLowEdge(i+1);
       UShort_t centHigh = fCentAxis->GetBinUpEdge(i+1);
@@ -426,13 +499,16 @@ struct dNdetaDrawer
       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
@@ -465,7 +541,7 @@ struct dNdetaDrawer
     // g->SetFillColor(color);
   }
   //__________________________________________________________________
-  void ModifyTitle(TNamed* h, const char* centTxt)
+  void ModifyTitle(TNamed* /*h*/, const char* /*centTxt*/)
   {
     return;
     // if (!centTxt || !h) return;
@@ -485,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,
@@ -499,79 +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,   fCentAxis ? color : color+2);
     SetAttributes(dndetaTruth,color);
     SetAttributes(dndetaSym,  color);
     SetAttributes(dndetaMCSym,fCentAxis ? color : color+2);
-    SetAttributes(tracks,     fCentAxis ? color : color+3);
     if (dndetaMC && fCentAxis) 
       dndetaMC->SetMarkerStyle(dndetaMC->GetMarkerStyle()+2);
     if (dndetaMCSym && fCentAxis) 
       dndetaMCSym->SetMarkerStyle(dndetaMCSym->GetMarkerStyle()+2);
-    if (dndetaTruth && fCentAxis) 
+    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));
-
-    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));
-      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));
        }
       }
-      // fOthers->Add(thisOther);
-    }
-    if (tracks) { 
-      if (!fRatios->GetHists() || 
-         !fRatios->GetHists()->FindObject(tracks->GetName()))
-       fRatios->Add(Ratio(dndeta, tracks, rmax));
+      // 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);
+      }
     }
-
     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));
+    if (fTruth) {
+      fRatios->Add(Ratio(dndeta,      fTruth, rmax));
+      fRatios->Add(Ratio(dndetaSym,   fTruth, rmax));
     }
+    return dndeta;
   }
   //__________________________________________________________________
   /** 
@@ -593,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);
@@ -622,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(), 
@@ -632,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);
   }
   //__________________________________________________________________
   /** 
@@ -654,21 +732,28 @@ struct dNdetaDrawer
     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());
     TH1*     hist = 0;
     TObjArray unique;
     unique.SetOwner();
+    Bool_t   sysErrSeen = false;
     while ((hist = static_cast<TH1*>(next()))) { 
-      TString n(hist->GetTitle());
-      if (n.Contains("mirrored")) continue;
-      if (unique.FindObject(n.Data())) continue;
-      TObjString* s = new TObjString(n);
+      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());
       TGraph* g = 0;
       while ((g = static_cast<TGraph*>(nexto()))) { 
@@ -682,7 +767,8 @@ struct dNdetaDrawer
        // l->AddEntry(hist, hist->GetTitle(), "pl");
       }
     }
-    unique.ls();
+
+    // Add legend entries for unique items only
     TIter nextu(&unique);
     TObject* s = 0;
     Int_t i = 0;
@@ -696,10 +782,35 @@ struct dNdetaDrawer
       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);
@@ -754,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");
@@ -774,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())))
@@ -799,6 +912,7 @@ 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();
@@ -812,6 +926,7 @@ struct dNdetaDrawer
                                           eS.Data(), 
                                           fCentAxis ? "by centrality" : 
                                           fTrigString->GetTitle()));
+    tt->SetTextColor(aliceBlue);
     tt->SetNDC();
     tt->SetTextFont(132);
     tt->SetTextAlign(33);
@@ -821,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);
@@ -832,6 +948,7 @@ struct dNdetaDrawer
       vt->SetNDC();
       vt->SetTextFont(132);
       vt->SetTextAlign(33);
+      vt->SetTextColor(aliceBlue);
       vt->Draw();
     }
     // results->Draw("nostack e1 same");
@@ -841,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();
   }
@@ -870,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);
@@ -886,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));
@@ -949,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);
@@ -958,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));
@@ -1364,8 +1488,8 @@ struct dNdetaDrawer
       r = 0; 
     }
     if (r) {
-      r->SetMarkerStyle(m1->GetMarkerStyle());
-      r->SetMarkerColor(m2->GetMarkerColor());
+      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()));
@@ -1518,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) { 
@@ -1532,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();
@@ -1550,7 +1682,7 @@ struct dNdetaDrawer
       }
     }
     if (ya) {
-      ya->SetTitle(ytitle);
+      // ya->SetTitle(h->GetYTitle());
       ya->SetDecimals();
       // ya->SetTicks("+-");
       ya->SetNdivisions(ynDiv);
@@ -1559,33 +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)
 {
@@ -1641,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$ 
  * 
@@ -1660,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)