GraphSysErr.C

Go to the documentation of this file.

# include <TNamed.h>
# include <TAttMarker.h>
# include <TAttLine.h>
# include <TAttFill.h>
#ifndef GraphSysErr_C
# define GraphSysErr_C
# include <TList.h>
# ifndef __CINT__
#  include <TStyle.h>
#  include <TGraph.h>
#  include <TGraphErrors.h>
#  include <TGraphAsymmErrors.h>
#  include <TMultiGraph.h>
#  include <TH1.h>
#  include <TMath.h>
#  include <TFitResultPtr.h>
#  include <TF1.h>
#  include <TROOT.h>
#  include <TSystem.h>
#  include <TDatime.h>
#  include <TList.h>
#  include <TBrowser.h>
#  include <TPRegexp.h>
#  include <TClonesArray.h>
#  include <TArrayL64.h>
#  include <TLine.h>
#  include <TClass.h>
#  include <TObjString.h>
#  include <iostream>
#  include <iomanip>
#  include <fstream>
#  include <sstream>
// #  define TOKEN(C,I) static_cast<TObjString*>(C->At(I))->String()
#  define INCOMPAT_CMN_AS_QUAL 1
# else 
# include <iosfwd>
class TGraph;
class TGraphErrors;
class TGraphAsymmErrors;
class TMultiGraph;
class TFitResultPtr;
class TF1;
class TList;
class TBrowser;
class TArrayD;
class TClonesArray;
class TLine;
# endif


class GraphSysErr : public TNamed, public TAttMarker, 
        public TAttLine, public TAttFill 
{
public:
  enum { 
    kDraw    = 0x1,
    kImport  = 0x2, 
    kExport  = 0x4,
    kRatio   = 0x8,
    kVerbose = 0  // Set to OR of above to enable verbose/debug output 
  };
  typedef TGraphAsymmErrors Graph;
  enum EDrawOption_t { 
    kNormal  = 0, //     - Line with ticks
    kNoTick  = 1, // Z0  - Line with no ticks 
    kArrow   = 2, // >0  - Linw with arrows 
    kRect    = 3, // 20  - Rectangle w/fill
    kBox     = 4, // 50  - Rectangle w/fill & line 
    kFill    = 5, // 30  - Filled area 
    kCurve   = 6, // 40  - Filled smoothed area 
    kHat     = 7, // []0 - Hats 
    kBar     = 8, // ||0 - A bar 
    kNone    = 9, // XP  - No errors
    kLine    = 10, // XC
    kConnect = 11 // XL
  };
  enum EChi2Type {
    kExperimentExperiment, // Ignore errors on both
    kExperimentModel,      // Ignore errors on second
    kModelModel,           // Use errors on both
    kExperimentTruth       // Use errors on first, take second to be exact.
  };
  enum ERatioOption {
    kMax          = 0x00001, // Error is largest relative error
    kMin          = 0x00002, // Error is largest relative error
    kCancel       = 0x00004,
    kDenomRel     = 0x00008,
    kRatioDefault = kCancel
  };
  enum {
    kUsedBit = (1 << 18),
    kOnlyWeightBit = (1 << 19)
  };
  GraphSysErr() 
    : TNamed(), 
      TAttMarker(1,20,1), 
      TAttLine(1,1,1),
      TAttFill(0,0),
      fPoint2Point(), 
      fCommon(), 
      fData(0), 
      fDrawn(0),
      fCounter(0),
      fSumFill(0,0),
      fSumLine(1,1,1),
      fSumTitle(),
      fSumOption(0),
      fCommonSumFill(0,0),
      fCommonSumLine(1,1,1),
      fCommonSumTitle(),
      fCommonSumOption(0),
      fDataOption(0), 
      fXTitle(""),
      fYTitle(""),
      fMap(0),
      fQualifiers(0),
      fStatRelative(false)
  {
    
  }
  GraphSysErr(Int_t n) 
    : TNamed("sysErrGraph","Data"), 
      TAttMarker(1,20,1), 
      TAttLine(1,1,1),
      TAttFill(0,0),
      fPoint2Point(), 
      fCommon(), 
      fData(0), 
      fDrawn(0),
      fCounter(0),
      fSumFill(0,0),
      fSumLine(1,1,1),
      fSumTitle("Uncorr.Errors"),
      fSumOption(0),
      fCommonSumFill(0,0),
      fCommonSumLine(1,1,1),
      fCommonSumTitle("Corr.Errors."),
      fCommonSumOption(0),
      fDataOption(0), 
      fXTitle(""),
      fYTitle(""),
      fMap(0),
      fQualifiers(0),
      fStatRelative(false)
  {
    fPoint2Point.SetName("point2point");
    fPoint2Point.SetOwner();
    fCommon.SetName("common");
    fCommon.SetOwner();
    MakeDataGraph(n);
  }
  GraphSysErr(const char* name, const char* title, Int_t n=0)
    : TNamed(name,title), 
      TAttMarker(1,20,1), 
      TAttLine(1,1,1),
      TAttFill(0,0),
      fPoint2Point(), 
      fCommon(), 
      fData(0), 
      fDrawn(0),
      fCounter(0),
      fSumFill(0,0),
      fSumLine(1,1,1),
      fSumTitle("Errors"),
      fSumOption(0),
      fCommonSumFill(0,0),
      fCommonSumLine(1,1,1),
      fCommonSumTitle("Corr.Errors."),
      fCommonSumOption(0),
      fDataOption(0), 
      fXTitle(""),
      fYTitle(""),
      fMap(0),
      fQualifiers(0),
      fStatRelative(false)
  {
    fPoint2Point.SetOwner();
    fCommon.SetOwner();
    fCommon.SetName("common");
    fPoint2Point("point2point");
    MakeDataGraph(n);
  }
  GraphSysErr(const GraphSysErr& other) 
    : TNamed(other), 
      TAttMarker(other), 
      TAttLine(other),
      TAttFill(other),
      fPoint2Point(), 
      fCommon(), 
      fData(0),
      fDrawn(0),
      fCounter(0),
      fSumFill(other.fSumFill),
      fSumLine(other.fSumLine),
      fSumTitle(other.fSumTitle),
      fSumOption(other.fSumOption),
      fCommonSumFill(other.fCommonSumFill),
      fCommonSumLine(other.fCommonSumLine),
      fCommonSumTitle(other.fCommonSumTitle),
      fCommonSumOption(other.fCommonSumOption),
      fDataOption(other.fDataOption), 
      fXTitle(other.fXTitle),
      fYTitle(other.fYTitle),
      fMap(0),
      fQualifiers(0),
      fStatRelative(false)
  {
    if (other.fData) 
      fData = static_cast<Graph*>(other.fData->Clone());
 
    fPoint2Point.SetOwner();
    fCommon.SetOwner();
    fCommon.SetName(other.GetName());
    fPoint2Point(other.GetName());

    TIter nextC(&other.fCommon);
    HolderCommon* common = 0;
    while ((common = static_cast<HolderCommon*>(nextC()))) {
      fCommon.Add(new HolderCommon(*common));
      fCounter++;
    }

    TIter nextP(&other.fPoint2Point);
    HolderP2P* p2p = 0;
    while ((p2p = static_cast<HolderP2P*>(nextP()))) {
      fPoint2Point.Add(new HolderP2P(*p2p));
      fCounter++;
    }

    CopyKeys(&other, "ar");
    // Print("all");
  }
  virtual ~GraphSysErr()
  {
    fPoint2Point.Delete();
    fCommon.Delete();
    if (fData)         delete fData;
    if (fDrawn)        delete fDrawn;
    if (fMap)          delete fMap;
    if (fQualifiers)   delete fQualifiers;
  }
  GraphSysErr& operator=(const GraphSysErr& other) 
  {
    if (&other == this) return *this;
    
    other.TNamed::Copy(*this);
    other.TAttMarker::Copy(*this);
    other.TAttLine::Copy(*this);
    other.TAttFill::Copy(*this);
    fPoint2Point.Clear();
    fCommon.Clear();
    
    if (fData) delete fData;
    fData = 0;
    if (other.fData)
      fData = static_cast<Graph*>(other.fData->Clone());

    if (fDrawn) { delete fDrawn; fDrawn = 0; }

    fCounter = other.fCounter;

    other.fSumFill.Copy(fSumFill);
    other.fSumLine.Copy(fSumLine);
    fSumTitle     = other.fSumTitle;
    fSumOption    = other.fSumOption;
    
    other.fCommonSumFill.Copy(fCommonSumFill);
    other.fCommonSumLine.Copy(fCommonSumLine);
    fCommonSumTitle     = other.fCommonSumTitle;
    fCommonSumOption    = other.fCommonSumOption;
    
    fXTitle       = other.fXTitle;
    fYTitle       = other.fYTitle;
    fStatRelative = other.fStatRelative;

    TIter nextC(&other.fCommon);
    HolderCommon* common = 0;
    while ((common = static_cast<HolderCommon*>(nextC()))) 
      fCommon.Add(new HolderCommon(*common));

    TIter nextP(&other.fPoint2Point);
    HolderP2P* p2p = 0;
    while ((p2p = static_cast<HolderP2P*>(nextP()))) 
      fPoint2Point.Add(new HolderP2P(*p2p));    

    if (fMap) delete fMap;
    fMap = 0;
    if (other.fMap)
      fMap = static_cast<TList*>(other.fMap->Clone());

    if (fQualifiers) delete fQualifiers;
    fQualifiers = 0;
    if (other.fQualifiers)
      fQualifiers = static_cast<TList*>(other.fQualifiers->Clone());
    
    return *this;
  }
  /* @} */

  virtual void ls(Option_t* option="") const
  {
    Print(option);
  }
  virtual void Print(Option_t* option="R") const
  {
    gROOT->IndentLevel();
    std::cout << GetName() << ": " << GetTitle() << std::endl;

    TString opt(option);    
    opt.ToUpper();
    if (opt.IsNull()) return;

    Bool_t all  = opt.Contains("ALL");
    Bool_t keys = all || opt.Contains("KEY");
    Bool_t qual = all || opt.Contains("QUAL");
    Bool_t sys  = all || opt.Contains("SYS");
    Bool_t cmn  = sys || opt.Contains("COMMON");
    Bool_t p2p  = sys || opt.Contains("P2P");
    Bool_t poi  = all || opt.Contains("XY");
    Bool_t attr = all || opt.Contains("ATTR");
    Bool_t sum  = opt.Contains("SUM");
    
    gROOT->IncreaseDirLevel();
    if (fMap && keys) {
      gROOT->IndentLevel();
      std::cout << "Key/value pairs: " << std::endl;
      gROOT->IncreaseDirLevel();
      TIter nextKV(fMap);
      TObject* kv = 0;
      while ((kv = nextKV())) {
  gROOT->IndentLevel();
  std::cout << '"' << kv->GetName() << '"' << "\t"
      << '"' << kv->GetTitle() << '"' << "\n";
      }
      gROOT->IndentLevel();
      std::cout << "\"XTitle\"\t\"" << fXTitle << "\"\n";
      gROOT->IndentLevel();
      std::cout << "\"YTitle\"\t\"" << fYTitle << "\"\n";      
      // fMap->Print(option);
      gROOT->DecreaseDirLevel();
    }
    if (fQualifiers && qual) {
      gROOT->IndentLevel();
      std::cout << "Qualifier pairs: " << std::endl;
      gROOT->IncreaseDirLevel();
      TIter nextQ(fQualifiers);
      TObject* q = 0;
      while ((q = nextQ())) {
  gROOT->IndentLevel();
  std::cout << '"' << q->GetName() << '"' << "\t"
      << '"' << q->GetTitle() << '"' << "\n";
      }
      // fQualifiers->ls(option);
      gROOT->DecreaseDirLevel();
    }
    if (attr) {
      gROOT->IndentLevel();
      Printf(" [option: %3d line (c/s/w):%3d/%1d/%2d fill (c/s):%3d/%4d marker (c/s/s):%3d/%2d/%f]",
       fDataOption, GetLineColor(), GetLineStyle(), GetLineWidth(),
       GetFillColor(), GetFillStyle(),
       GetMarkerColor(), GetMarkerStyle(), GetMarkerSize());
      gROOT->IndentLevel();
      Printf(" [sum title: %s "
       "option: %3d line (c/s/w):%3d/%1d/%2d fill (c/s):%3d/%4d",
       fSumTitle.Data(), fSumOption,
       fSumLine.GetLineColor(),
       fSumLine.GetLineStyle(),
       fSumLine.GetLineWidth(), 
       fSumFill.GetFillColor(),
       fSumFill.GetFillStyle());
      gROOT->IndentLevel();
      Printf(" [common sum title: %s "
       "option: %3d line (c/s/w):%3d/%1d/%2d fill (c/s):%3d/%4d",
       fCommonSumTitle.Data(),
       fCommonSumOption,
       fCommonSumLine.GetLineColor(),
       fCommonSumLine.GetLineStyle(),
       fCommonSumLine.GetLineWidth(), 
       fCommonSumFill.GetFillColor(),
       fCommonSumFill.GetFillStyle());
    }

    if (cmn) {
      gROOT->IndentLevel();
      std::cout << "Commons: " << std::endl;
      gROOT->IncreaseDirLevel();
      TIter nextC(&fCommon);
      TObject* c = 0;
      while ((c = nextC())) {
  gROOT->IndentLevel();
  c->Print(option);
      }
      // fCommon.Print(option);
      gROOT->DecreaseDirLevel();
    }

    if (p2p) {
      gROOT->IndentLevel();
      std::cout <<  "Point-to-point: " << std::endl;
      gROOT->IncreaseDirLevel();
      TIter nextP(&fPoint2Point);
      TObject* p = 0;
      while ((p = nextP())) {
  gROOT->IndentLevel();
  p->Print(option);
      }
      // fPoint2Point.Print(option);
      gROOT->DecreaseDirLevel();
    }
    if (poi) {
      // gROOT->IndentLevel();
      gROOT->IndentLevel();
      std::cout <<  "Points: " << std::endl;
      gROOT->IncreaseDirLevel();
      TString hs;
      hs.Form(" %3s: %9s %9s %9s -> %9s %9s %9s",
        "#", "X", "-dX", "+dX", "Y", "-stat", "+stat");
      gROOT->IndentLevel();
      std::cout << hs;
      if (sum) {
  hs.Form(" %9s %9s", "-sys", "+sys");
  std::cout << hs;
      }
      else {
  TIter nextP(&fPoint2Point);
  HolderP2P* p = 0;
  while ((p = static_cast<HolderP2P*>(nextP()))) {
    Bool_t      rel  = p->IsRelative();
    const char* post = (rel ? "%" : " ");
    TString nm(p->GetTitle());
    if (nm.Length() > 7) {
      nm.Remove(4,nm.Length()-4);
      nm.Append("...");
    }
    hs.Form(" -%8s%s +%8s%s", nm.Data(), post, nm.Data(), post);
    std::cout << hs;    
  }
      }
      std::cout << std::endl;
      for (Int_t i = 0; i < GetN(); i++) {
  gROOT->IndentLevel();
  Double_t eyl, eyh, wyl, wyh;
  // point,stat,common,quad,nosqrt,eyl,eyh
  Double_t y = GetYandError(i,false,false,true,false,eyl,eyh,wyl,wyh);
  // Double_t y = GetY(i);
  TString s;
  s.Form(" %3d: %+8f -%7f +%7f -> %+8f -%7f +%7f",
         i, GetX(i), GetErrorXLeft(i), GetErrorXRight(i),
         y, GetStatErrorDown(i), GetStatErrorUp(i));
  std::cout << s;
  if (sum) {
    TString sSum;
    sSum.Form(" -%7f +%7f", eyl, eyh);
    std::cout << sSum;
  }
  // if (sum && p2p) std::cout << " [";
  else if (p2p) {
    TIter nextP(&fPoint2Point);
    HolderP2P* p = 0;
    while ((p = static_cast<HolderP2P*>(nextP()))) {
      Bool_t      rel  = p->IsRelative();
      Double_t    fac  = (rel ? (TMath::Abs(y) > 1e-9 ? 100/y : 0) : 1);
      const char* post = (rel ? "%" : " ");
      s.Form(" -%7f%s +%7f%s",
       fac * p->GetYDown(i), post,
       fac * p->GetYUp(i)  , post);
      std::cout << s;
    }
  }
  // if (sum && p2p) std::cout << " ]";
  std::cout << std::endl;
      } // for i 
      gROOT->DecreaseDirLevel();
    } // if poi
    gROOT->DecreaseDirLevel();
  } //*MENU*
  virtual Bool_t IsFolder() const { return true; }
  virtual void Browse(TBrowser* b)
  {
    if (fMap)        b->Add(fMap);
    if (fQualifiers) b->Add(fQualifiers);
    b->Add(&fCommon);
    b->Add(&fPoint2Point);
    if (fData)  b->Add(fData);
    if (fDrawn) b->Add(fDrawn);
  }
  /* @} */
  void Draw(Option_t* option="")
  {
    TString opt(option);
    opt.ToUpper();
    Bool_t   clear   = opt.Contains("CLEAR");
    Bool_t   axis    = opt.Contains("AXIS");
    // --- Optionally clear old stack --------------------------------
    if (clear && fDrawn) { 
      delete fDrawn;
      fDrawn = 0;
    }

    fDrawn = MakeMulti(option);
    if (!fDrawn) return;

    fDrawn->Draw(axis ? "A" : "");
    if (axis) { 
      fDrawn->GetHistogram()->SetXTitle(fXTitle);
      fDrawn->GetHistogram()->SetYTitle(fYTitle);
    }
  }  //*MENU*
  TFitResultPtr Fit(TF1* f1, Option_t* fitOption, Option_t* drawOption, 
        Axis_t min=0, Axis_t max=0)
  {
    TString dOpt(drawOption);
    dOpt.ToUpper();
    Bool_t   clear   = dOpt.Contains("CLEAR");
    Bool_t   axis    = dOpt.Contains("AXIS");

    if (clear && fDrawn) { 
      delete fDrawn;
      fDrawn = 0;
    }

    fDrawn = MakeMulti(drawOption);
    if (!fDrawn || 
  !fDrawn->GetListOfGraphs() || 
  !fDrawn->GetListOfGraphs()->First()) 
      return TFitResultPtr(-1);
    // fDrawn->GetListOfGraphs()->ls();
    Graph* g = static_cast<Graph*>(fDrawn->GetListOfGraphs()->First());
    
    TString fOpt(fitOption);
    fOpt.ToUpper();
    Bool_t noStore = fOpt.Contains("N");
    Bool_t noDraw  = fOpt.Contains("0"); 
    Bool_t range   = fOpt.Contains("R"); 
    fOpt.Append("N0");
    TFitResultPtr r = g->Fit(f1, fOpt.Data(), "", min, max);

    if (!noStore) {
      Int_t status = r;
      if (status == 0) {
  if (min < max && !range) f1->SetRange(min, max);
  fDrawn->GetListOfFunctions()->Add(f1);
      }
    }

    if (!noStore && !noDraw) {
      fDrawn->Draw(axis ? "A" : "");
      if (axis) { 
  fDrawn->GetHistogram()->SetXTitle(fXTitle);
  fDrawn->GetHistogram()->SetYTitle(fYTitle);
      }
    }

    return r;
  }
  TFitResultPtr Fit(const char* formula, 
        Option_t* fitOption, Option_t* drawOption, 
        Axis_t min=0, Axis_t max=0)
  {
    TString fname(formula);
    Bool_t  linear = fname.Contains("++");
    TF1*    f1     = 0;
    if (linear) f1 = new TF1(formula,formula,min,max);
    else { 
      f1 = static_cast<TF1*>(gROOT->GetFunction(formula));
      if (!f1) Warning("Fit", "Unknown function %s", formula);
    }
    if (!f1) return -1;
    
    return Fit(f1, fitOption, drawOption, min, max);
  }
  TMultiGraph* GetMulti(Option_t* option="") 
  {
    if (!fDrawn) fDrawn = MakeMulti(option);
    return fDrawn;
  }
  /* @} */
  void Scale(TF1* f, Double_t s=1)
  {
    HolderCommon* cmn = 0;
    TIter         cNext(&fCommon);
    while ((cmn = static_cast<HolderCommon*>(cNext()))) {
      if (cmn->IsRelative())
  // Do not scale relative errors - they scale by themselves
  continue;
      cmn->Set(cmn->GetYDown()/s, cmn->GetYUp()/s);
    }
    for (Int_t i = 0; i < GetN(); i++) {
      Double_t x = GetX(i);
      Double_t y = GetY(i);
      Double_t g = f->Eval(x);
      SetPoint(i, GetX(i), y/g);
      SetStatError(i, GetStatErrorDown(i)/g,GetStatErrorUp(i)/g);
      
      HolderP2P* p2p = 0;
      TIter      pNext(&fPoint2Point);
      while ((p2p = static_cast<HolderP2P*>(pNext()))) {
  Double_t fac = (p2p->IsRelative() ?
      (TMath::Abs(y) < 1e-9) ? 0 : 1/y : 1/g);
  SetSysError(p2p->GetUniqueID(), i,
        p2p->GetXLeft(i), p2p->GetXRight(i),
        fac*p2p->GetYDown(i), fac*p2p->GetYUp(i));        
      }
    }
  }
  void Scale(Double_t s)
  {
    HolderCommon* cmn = 0;
    TIter         cNext(&fCommon);
    while ((cmn = static_cast<HolderCommon*>(cNext()))) {
      if (cmn->IsRelative())
  // Do not scale relative errors - they scale by themselves
  continue;
      cmn->Set(s * cmn->GetYDown(), s * cmn->GetYUp());
    }
    for (Int_t i = 0; i < GetN(); i++) {
      Double_t y = GetY(i);
      SetPoint(i, GetX(i), s*y);
      SetStatError(i, s*GetStatErrorDown(i),s*GetStatErrorUp(i));
      
      HolderP2P* p2p = 0;
      TIter      pNext(&fPoint2Point);
      while ((p2p = static_cast<HolderP2P*>(pNext()))) {
  Double_t fac = (p2p->IsRelative() ?
      (TMath::Abs(y) < 1e-9) ? 0 : 1/y :
      s);
  SetSysError(p2p->GetUniqueID(), i,
        p2p->GetXLeft(i), p2p->GetXRight(i),
        fac*p2p->GetYDown(i), fac*p2p->GetYUp(i));        
      }
    }
  }
  Int_t Average(const TCollection* others, Bool_t sep=true)
  {
    const Double_t tol = 1e-9;
    TList          sharedC; sharedC.SetOwner();
    TList          sharedP; 
    TIter          oNext(others);
    GraphSysErr*   other = 0;
    Double_t       xMin  = +1e9;
    Double_t       xMax  = -1e9;
    Int_t          nX    = 0;
    Bool_t         first = true;
    Bool_t         nonSh = false;
    Int_t          npS   = 0;
    // Here, we do some initial investigations into the the passed
    // graphs. We should find all common errors that are shared amoung
    // the input graphs, and we need to check for the consistency of
    // those.
    while ((other = static_cast<GraphSysErr*>(oNext()))) {
      Int_t n = other->GetN();
      nX      = TMath::Max(nX,n);
      xMin    = TMath::Min(xMin,other->GetX(0)  -2*other->GetErrorXLeft(0));
      xMax    = TMath::Max(xMax,other->GetX(n-1)+2*other->GetErrorXRight(n-1));

      if (first) { 
  // For the first graph, we simply add the title of all the
  // common systematics.
  HolderCommon* cmn = 0;
  TIter         cNext(&other->fCommon);
  while ((cmn = static_cast<HolderCommon*>(cNext()))) {
    HolderCommon* o = static_cast<HolderCommon*>(cmn->Clone());
    sharedC.Add(o);
  }
  if (sep) {
    // and then the same for the point to point
    HolderP2P* p2p = 0;
    TIter      pNext(&other->fPoint2Point);
    while ((p2p = static_cast<HolderP2P*>(pNext()))) {
      sharedP.Add(p2p->Clone());
    }
  }
  
  // We also copy the keys
  CopyKeys(other, "ar");
      }
      else {
  // Otherwise, we check that this graph has all the common
  // errors in our temporary list, and if not, we remove that
  // error from the list
  TObjLink* cur  = sharedC.FirstLink();
  while (cur) {
    HolderCommon* o  = static_cast<HolderCommon*>(cur->GetObject());
    HolderCommon* c  = other->FindCompat(o, tol, true);
    if (!c) {
      // This error not found in the graph we're looking at.
      // Remove it from the list
      Info("Average", "Common systematic %s not found in %s",
     o->GetTitle(), other->GetName());
      other->Print("sys");
      TObjLink* keep = cur->Next();
      TObject*  obj  = sharedC.Remove(cur);
      cur            = keep;
      nonSh          = true;
      if (obj) delete obj;
      continue;
    }
    cur = cur->Next();
  } // while cur

  // Now check the point-to-point errors 
  cur = sharedP.FirstLink();
  while (cur) {
    HolderP2P* o  = static_cast<HolderP2P*>(cur->GetObject());
    HolderP2P* p  = other->FindCompat(o);
    if (!p) {
      // This error not found in the graph we're looking at.
      // Remove it from the list
      Info("Average", "P2P systematic %s not found in %s",
     o->GetTitle(), other->GetName());
      TObjLink* keep = cur->Next();
      TObject*  obj  = sharedP.Remove(cur);
      cur            = keep;
      npS++;
      if (obj) delete obj;
      continue;
    }
    cur = cur->Next();
  } // while cur    
      } // else
      first = false;
    }   // while cmn
    // Now we have the list of common errors that are shared amoung
    // all the input graphs.  So we define those new common errors on
    // this graph, and disable them on the inputs.
    TIter          nextC(&sharedC);
    HolderCommon*  com = 0;
    while ((com = static_cast<HolderCommon*>(nextC()))) {
      Int_t cId = DefineCommon(com->GetTitle(), com->IsRelative(),
             com->GetYDown(1), com->GetYUp(1), kBox);
      HolderCommon* c = FindCommon(cId);
      c->CopyAttr(com);
      // Info("Average", "Shared common error %s +%7f%s -%7f%s",
      //     c->GetTitle(),
      //     c->GetYDown(1), (c->IsRelative() ? "%" : ""),
      //     c->GetYUp(1),   (c->IsRelative() ? "%" : ""));

      // Now loop over all input graphs, find the shared commons, and
      // mark them as used.
      oNext.Reset();
      while ((other = static_cast<GraphSysErr*>(oNext()))) {
  HolderCommon* oc = other->FindCompat(c, tol, true);
  if (!oc) {
    Error("Average", "This should NOT happen");
    other->Print("sys");
    continue;
  }
  oc->SetBit(kUsedBit);
      } // while other
    } // while com
    sharedC.Clear();
    // If we found some common errors that are not shared, we need to
    // make a point-to-point error for absorbing the remaining common
    // errors.
    Int_t bId = 0;
    if (nonSh || !sep) {
      bId = DeclarePoint2Point("Blob", false, kBox);
      SetSysFillColor(bId, kRed+2);
      SetSysLineColor(bId, kRed+2);
      SetSysFillStyle(bId, 3004);
    }
    
    // At this point, we have found all the common shared errors and
    // added those to ourselves. Those that are not shared go into the
    // blob point-to-point error.
    // 
    // Next step, is to deal with the shared point-to-point errors.
    // If we have shared point-to-point errors, we need to deal with
    // them separately.  To do that, we make a list of identifiers
    // that we should inspect at each X value.  The identifiers point
    // to the p2p errors of each input graph.
    //
    TArrayL64     mp2p;
    Int_t         dOf = 6;
    if (sep) {
      TIter       nextP(&sharedP);
      HolderP2P*  p2p = 0;
      Int_t       j   = 0;
      mp2p.Set(100);
      while ((p2p = static_cast<HolderP2P*>(nextP()))) {
  // Info("Average", "Loop %d shared P2P %s", j,  p2p->GetTitle());
  Int_t      pId = DeclarePoint2Point(p2p->GetTitle(),
            p2p->IsRelative(), kBox);
  HolderP2P* p   = FindP2P(pId);
  p->CopyAttr(p2p);
  Info("Average", "Shared point-to-point error %s", p->GetTitle());
  
  // Now loop over all input graphs, find the shared commons, and
  // mark them as used.
  oNext.Reset();
  Long64_t mask = 0;
  Int_t    off  = 0;
  while ((other = static_cast<GraphSysErr*>(oNext())) && off < 64) {
    HolderP2P* op = other->FindCompat(p, true);
    if (!op) {
      Error("Average", "This should NOT happen");
      continue;
    }
    mask |= ((op->GetUniqueID() & 0x3F) << off);
    off  += dOf;
    // Info("Average", "Mark %s as used", op->GetTitle());
    op->SetBit(kUsedBit); 
  } // while other
  if (off >= 64) {
    Warning("Average",
      "Some shared point-to-point errors could not be encoded "
      "becasue we have too many (%d>%d) input graphs",
      others->GetEntries(), 64/dOf);
  }
  mp2p[pId] = mask;
  j++;
      } // while p2p
      // Info("Average", "Declared new shared P2P - mapping:");
      // for (Int_t j = 1; j < mp2p.GetSize(); j++) {
      //   if (mp2p[j] == 0) continue;
      //   rintf("%3d/%3d", j, mp2p.GetSize());
      //   for (Int_t k = 0; k < 64; k += dOf)
      //     printf(" %3d", (mp2p[j] >> k) & 0x3f);
      //   printf("\n");
      // }
      
      // sharedP.Clear();
      // Now we have all shared point-to-point errors defined, and we
      // have the ids for each graph in the shar array (6bits each).
      //
      // Next, we will declare _all_ other our point-to-point errors
      // separately on this graph.  That means that the point-to-point
      // errors can be 0 in some ranges.
      oNext.Reset();
      Int_t       off  = 0;
      while ((other = static_cast<GraphSysErr*>(oNext()))) {
  Long64_t*   mask = 0;
  TIter       oNextP(&(other->fPoint2Point));
  HolderP2P*  p2p = 0;
  while ((p2p = static_cast<HolderP2P*>(oNextP())) && off < 64) {
    if (p2p->TestBit(kUsedBit)) {
      // Info("Average", "%s marked as used", p2p->GetTitle()); 
      continue; // Already in the shared part
    }
    
    // The point-to-point error should go into the weight, but not
    // in the result.
    p2p->SetBit(kOnlyWeightBit);
    Int_t      pId = 0;
    HolderP2P* p   = FindCompat(p2p);
    if (!p) {
      pId = DeclarePoint2Point(p2p->GetTitle(), p2p->IsRelative(), kBox);
      p   = FindP2P(pId);
      p->CopyAttr(p2p);
      mask =  &(mp2p[pId]);

    }
    *mask |= ((p2p->GetUniqueID() & 0x3F) << off);
  }
  off += dOf;
      } // while other
      if (off >= 64) {
  Warning("Average",
    "Some shared point-to-point errors could not be encoded "
    "becasue we have too many (%d>%d) input graphs",
    others->GetEntries(), 64/dOf);
      }
      // Info("Average", "Declared new shared P2P - mapping:");
      // for (Int_t j = 0; j < mp2p.GetSize(); j++) {
      //  if (mp2p[j] == 0) continue;
      //  printf("%3d/%3d", j, mp2p.GetSize());
      //  for (Int_t k = 0; k < 64; k += dOf)
      //    printf(" %3d", (mp2p[j] >> k) & 0x3f);
      //  printf("\n");
      // }
    }
    
    // Now search for the X values to evaluate at 
    Double_t x    = xMin;
    Double_t oldX = xMin;
    Int_t    cnt  = 0;
    TArrayD  xa(2*nX); // Assume no more than 100 points
    TArrayD  exla(2*nX);
    TArrayD  exha(2*nX);
    do {
      // Find the next X to evaluate at
      oNext.Reset();
      Double_t nextX = xMax;
      Double_t texl  = 0;
      Double_t texh  = 0;
      // Info("Average", "nextX=%f xMax=%f", nextX, xMax);
      while ((other = static_cast<GraphSysErr*>(oNext()))) {
  for (Int_t i = other->GetUniqueID(); i < other->GetN(); i++) {
    if (other->GetX(i) <= x) {
      // Info("Average", "- x_i=%f < x=%f", g->GetX(i), x);
      continue;
    }
    if (other->GetX(i) < nextX) {
      // Info("Average", "Set x to %f (was %f)", g->GetX(i), nextX);
      nextX = other->GetX(i);
      texl  = TMath::Max(texl, other->GetErrorXLeft(i));
      texh  = TMath::Max(texh, other->GetErrorXRight(i));
    }
    else {
      // g->SetUniqueID(i);
      // Info("Average", "No more (%s) here at x=%f (next=%f)",
      //      g->GetName(), g->GetX(i), nextX);
      break;
    }
  }
      }
      if (nextX >= xMax) {
  // Info("Average", "Next x is too large %f>%f", nextX, xMax);
  break;
      }

      if (cnt >= xa.GetSize()) {
  Warning("Average", "increasing size of X cache");
  xa.Set(1.5*cnt);
  exla.Set(1.5*cnt);
  exha.Set(1.5*cnt);
      }
      // Info("Average", "Set next x to %f +/- (%f,%f)", nextX, texl, texh);
      xa[cnt]     = nextX;
      exla[cnt]   = texl;
      exha[cnt]   = texh;
      Double_t dx = texh*2;//(nextX-oldX);
      oldX        = x;
      x           = nextX+dx/4;
      cnt++;
    } while (cnt < 1000); // Do not go for more than 1000 points
    // At this point we have
    //
    //  - All common errors defined, and we have set the shared ones
    //  - All point-to-point errors defined, and we know which input
    //    graph contributes to each.
    //  - A list of X values to evaluate each graph at. 
        
    // Print("sys");
    // Now we loop over the list of X values and evaluate each graph
    // at that X value.
    for (Int_t i = 0; i < cnt; i++) {
      Double_t xi   = xa[i];
      Double_t texl = exla[i];
      Double_t texh = exha[i];
      Double_t nexl = TMath::Abs(xi - (i     == 0   ? xa[i+1] : xa[i-1]))/2;
      Double_t nexh = TMath::Abs(xi - (i + 1 == cnt ? xa[i-1] : xa[i+1]))/2;
      Double_t exl  = TMath::Min(texl, nexl);
      Double_t exh  = TMath::Min(texh, nexh);
      // Info("Average", "@ %3d x=%f, exl=min(%f,%f)=%f exh=min(%f,%f)=%f",
      //      i, xi, texl, nexl, exl, texh, nexh, exh);

      
      // Now loop again - this time calculating our data point
      Double_t sum   = 0;
      Double_t sumW  = 0;
      // Double_t sumWS = 0;
      Double_t sumE  = 0;
      Double_t sumES = 0;
      oNext.Reset();

      // Use our combiner framework for calculating the averaged
      // result and errors.  Note, we define two combiners.  One that
      // include all non-shared common and point-to-point errors
      // (centroid), and one that doesn't have the non-shared common
      // and point-to-point errors as well as the split point-to-point
      // errors (error).  The first one is used to set the centroid of
      // the data, while the second is used to set the combined error
      LinearSigmaCombiner centroid;
      LinearSigmaCombiner error;
      LinearSigmaCombiner stats;
      Int_t               j = 0;
      TArrayI             ai1(others->GetEntries());
      TArrayI             ai2(others->GetEntries());
      TArrayD             ax(others->GetEntries());
      while ((other = static_cast<GraphSysErr*>(oNext()))) {
  Bool_t cmn    = true;   // Include the common errors
  Bool_t stat   = false;  // Include the statistical errors
  Bool_t quad   = true;   // Add in quadrature
  Bool_t nosqrt = false;  // Do not return squared error 

  Double_t y, eyl, eyh, wyl, wyh, syl, syh;

  Int_t  fret   = other->FindPoint(xi, ai1[j], ai2[j]);
  if (fret < -1) { j++; continue; } // Nothing here
  if (fret >= 0) {
    y = other->GetYandError(fret,cmn,stat,quad,nosqrt,eyl, eyh, wyl, wyh);
    syl = other->GetStatErrorDown(fret);
    syh = other->GetStatErrorUp(fret);
  }
  else {
    Double_t eyl1, eyh1, wyl1, wyh1;
    Double_t eyl2, eyh2, wyl2, wyh2;
    Double_t x1  = other->GetX(ai1[j]);
    Double_t x2  = other->GetX(ai2[j]);
    Double_t y1  = other->GetYandError(ai1[j],cmn,stat,quad,nosqrt,
               eyl1,eyh1,wyl1,wyh1);
    Double_t y2  = other->GetYandError(ai2[j],cmn,stat,quad,nosqrt,
               eyl2,eyh2,wyl2,wyh2);
    Double_t syl1 = other->GetStatErrorDown(ai1[j]);
    Double_t syh1 = other->GetStatErrorUp(ai1[j]);
    Double_t syl2 = other->GetStatErrorDown(ai2[j]);
    Double_t syh2 = other->GetStatErrorUp(ai2[j]);
    Double_t dx = (x2-x1);
    ax[j]       = (xi-x1)/dx;
    y           = y1   + ax[j] * (y2-y1);
    eyl         = eyl1 + ax[j] * (eyh1-eyl1);
    eyh         = eyh2 + ax[j] * (eyh2-eyl2);
    wyl         = wyl1 + ax[j] * (wyh1-wyl1);
    wyh         = wyh2 + ax[j] * (wyh2-wyl2);
    syl         = syl1 + ax[j] * (syh1-syl1);
    syh         = syh2 + ax[j] * (syh2-syl2);
  }
  
  centroid.Add(y, wyl, wyh);
  error.Add(y, eyl, eyh);
  stats.Add(y, syl, syh);
  j++;
      } // while (other)
      centroid.Calculate();
      // Info("Average", "Centroid");
      // centroid.Print();
      // centroid.fResult->Print();
      Double_t newY = centroid.fResult->fX;
      SetPoint(i, xi, newY);
      SetPointError(i, exl, exh);

      stats.Calculate();
      // Info("Average", "Stats");
      // stats.Print();
      // stats.fResult->Print();
      stats.Calculate();
      Double_t fac = IsStatRelative() ? 1/newY : 1;
      SetStatError(i, fac*stats.fResult->fEl,fac*stats.fResult->fEh);

      if (bId > 0) {
      // Calculations for error
  error.Calculate();
  // Info("Average", "Error");
  // error.Print();
  // error.fResult->Print();

  fac = IsRelative(bId) ? 1/newY : 1;
  SetSysError(bId,i,0,0,fac*error.fResult->fEl,fac*error.fResult->fEh);
      }

      // Now, we have set the blob point-to-point error.  Next, we
      // need to do the point-to-point erros that should be transfered
      // to the output.  We have stored masks of IDs in mp2p, so we
      // loop over that array first.  The lower 6 bits is the output
      // Id, while the remaining bits (counting from the low side) are
      // ids for each input graph.
      for (Int_t j = 1; j < mp2p.GetSize(); j++) {
  Long64_t mask = mp2p[j];
  if (mask == 0) continue;
  Int_t    tid  = j; // (mask & 0x3f);
  HolderP2P* tp2p = FindP2P(tid);
  if (!tp2p) {
    Warning("Average", "No target point-to-point error at %d for id=%d",
      j, tid);
    continue;
  }
  // Zero by default
  SetSysError(tid, i, 0, 0);
  Long64_t rem  = mask;
  if (rem == 0) 
    // There's no IDs for this, explicitly set the value to 0
    continue;

  // Now we should loop over all the input graphs and get the
  // error and finally combine them into one.  Again, we use a
  // linear sigma approximation from our combiner framework to
  // do the calculations.
  LinearSigmaCombiner sc;
  Int_t k = 0; // Index into index array and for off set in mask
  oNext.Reset();
  while ((other = static_cast<GraphSysErr*>(oNext()))) {
    if (ai1[k] == -1 && ai2[k] == -1) {
      // This graph does not contribute at this point
      // Info("Average", "No contribution at %f (%d:%d,%d) for %s to %s",
      //      xi, k, ai1[k], ai2[k], other->GetName(), tp2p->GetTitle());
      k++;
      continue;
    }
    Int_t sid = ((rem >> dOf*k) & 0x3f);
    if (sid <= 0) {
      // This graph does not contribute to this error
      // Info("Average",
      //      "No contribution for %s to %s, not found [%d:%d,%d]",
      //      other->GetName(), tp2p->GetTitle(), k, ai1[k], ai2[k]);
      k++;
      continue;
    }
    HolderP2P* sp2p = other->FindP2P(sid);
    if (!sp2p) {
      Warning("Average", "Couldn't find point-to-point error %d (%s)"
        "in graph %d (%s) - should not happen",
        sid, tp2p->GetTitle(), k, other->GetName());
      k++;
    }
    Double_t sy   = 0;
    Double_t seyl = 0;
    Double_t seyh = 0;
    if (ai1[k] == ai2[k]) {
      // Exact point
      sy   = other->GetY(ai1[k]);
      seyl = other->GetSysErrorYDown(sid, ai1[k]);
      seyh = other->GetSysErrorYUp(sid, ai1[k]);
    }
    else {
      // Interpolate
      Double_t sy1   = other->GetY(ai1[k]);
      Double_t sy2   = other->GetY(ai2[k]);
      sy             = sy1 + ax[k] * (sy2-sy1);
      seyl           = sp2p->GetYDown(ai1[k],ai2[k],ax[k]);
      seyh           = sp2p->GetYUp(ai1[k],ai2[k],ax[k]);
    }
    //Info("Average", "Contribution (%d:%d-%d:%f,%f) to %s from %s",
    //     k,ai1[k],ai2[k],seyl,seyh,tp2p->GetTitle(),other->GetName());
    sc.Add(sy, seyl, seyh);
    k++;
  } // while other
  sc.Calculate();
  // Info("Average", "%s", tp2p->GetTitle());
  // sc.Print();
  // sc.fResult->Print();
  fac = IsRelative(tid) ? 1/newY : 1;
  SetSysError(tid, i, 0, 0, fac*sc.fResult->fEl, fac*sc.fResult->fEh);
      } // for j (# of p2p)
    } // for i (points)

    oNext.Reset();
    while ((other = static_cast<GraphSysErr*>(oNext()))) 
      other->ClearUsed();
    ClearUsed();
      
    return bId;
  }
  Double_t Integral(Double_t& error,
        Option_t* option="quad sum stat",
        UShort_t  first=0,
        Short_t   last=-1)
  {
    if (last < 0) last = GetN()-1;
    Double_t sum  = 0;
    Double_t err2 = 0;
    TString  opts(option); opts.ToLower();
    Bool_t   cmn  = opts.Contains("comm");
    Bool_t   stat = opts.Contains("stat");
    Bool_t   quad = opts.Contains("quad");
    
    for (Short_t i = first; i <= last; i++) {
      Double_t eyl, eyh;
      Double_t y   = GetYandError(i, cmn, stat, quad, false, eyl, eyh);
      Double_t eym = (eyh+eyl)/2;
      Double_t x   = GetX(i);
      Double_t exl = GetErrorXLeft(i);
      Double_t exr = GetErrorXRight(i);
      Double_t x1  = x - exl;
      Double_t x2  = x + exr;
      Double_t prv = x1;
      Double_t nxt = x2;
      if (i-1 >= 0)     prv = GetX(i-1)+GetErrorXRight(i-1);
      if (i+1 < GetN()) nxt = GetX(i+1)-GetErrorXLeft(i+1);
      if (prv > x1) x1 -= (prv-x1)/2; // Adjust to half-ways 
      if (nxt < x2) x2 -= (x2-nxt)/2; // Adjust to half-ways
      Double_t wdt = (x2-x1);
      if (wdt < 0) continue; // Do not take this point
      sum  += wdt * y;
      err2 += TMath::Power(eym*wdt,2);
    }
    error = TMath::Sqrt(err2);
    return sum;
  }
  static Bool_t NextPoint(Int_t i,
        const GraphSysErr* num,
        const GraphSysErr* denom,
        Double_t& x,
        Double_t& dY,
        Double_t& dEyl,
        Double_t& dEyh,
        Double_t& nY,
        Double_t& nEyl,
        Double_t& nEyh)
  {
    x  = num->GetX(i);
    if (!denom->FindYandError(x, false, true, true, false, dY, dEyl, dEyh))
      return false; 
    if (TMath::Abs(dY) < 1e-9)
      return false;
    
    nY   = num->GetYandError(i,false,true,true,false,nEyl, nEyh);
    return true;
  }
  
  static GraphSysErr* Ratio(const GraphSysErr* num,
          const GraphSysErr* den,
          UInt_t             flags=kRatioDefault,
          Double_t           fac=1)
  {
    // num->Print("sys");
    // den->Print("sys");
    num->ClearUsed();
    den->ClearUsed();
    GraphSysErr* ret = new GraphSysErr(1);
    ret->CopyAttr(num);
    ret->CopyKeys(num, "ar");
    Bool_t   emax    = (flags & kMax);
    Bool_t   cmin    = (flags & kMin);
    Bool_t   cmax    = (flags & kMax);
    Bool_t   denrel  = (flags & kDenomRel);
    Bool_t   cancc   = (flags & kCancel);
    Int_t    id      = 0;
    Int_t    notUsed = 0;
    // Loop over all common errors in the numerator, and see if we
    // have the same common error in the denominator.  If we do, we
    // can (partially) cancel this common error.
    //
    // Technically, we mark each used common error as used, and create
    // a new common error in the output graph.
    //
    // Common errors that are not shared by both the numerator and the
    // denominator will not be cancelled.  Instead, if they are
    // relative, they will be added to the output graph directly.
    // Finally, common systematic errors that are not cancelled or are
    // not relative, will go into the blob of point-to-point errors
    TIter nextNC(&num->fCommon);
    HolderCommon* comN = 0;
    while ((comN = static_cast<HolderCommon*>(nextNC()))) {      
      Bool_t        rel  = comN->IsRelative();
      Double_t      eyl  = comN->GetYDown(1);
      Double_t      eyh  = comN->GetYUp(1);
      Int_t         idD  = den->FindId(comN->GetTitle());      
      HolderCommon* comD = 0;
      if (idD != 0 &&                       // Check it exists
    (comD = den->FindCommon(idD)) &&  // Check it is a common
    (rel == comD->IsRelative())) {    // Same as this
  // We found the same error in the denominator, so we need to
  // cancel errors here.
  Double_t dEyl = comD->GetYDown(1);
  Double_t dEyh = comD->GetYUp(1);
  Double_t nEyl = comN->GetYDown(1);
  Double_t nEyh = comN->GetYUp(1);
  if (cmin) {
    eyl = TMath::Min(nEyl, dEyl);
    eyh = TMath::Min(nEyh, dEyh);
  }
  else if (cmax) {
    eyl = TMath::Max(nEyl, dEyl);
    eyh = TMath::Max(nEyh, dEyh);
  }
  else if (cancc) { 
    eyl  = TMath::Sqrt(TMath::Abs(dEyl*dEyl-nEyl*nEyl));
    eyh  = TMath::Sqrt(TMath::Abs(dEyh*dEyh-nEyh*nEyh));
  }
  const char* post = (rel ? "%" : "");
  Double_t    pfc  = (rel ? 100 : 1);
          
  if (kVerbose & kRatio)
    ::Info("Ratio", "Cancelling the common systematic error %s "
     "between numerator (-%f%s +%f%s) "
     "and denominator (-%f%s +%f%s): -%f%s + %f%s",
     comD->GetTitle(),
     pfc*nEyl, post, pfc*nEyh, post,
     pfc*dEyl, post, pfc*dEyh, post,
     pfc*eyl,  post, pfc*eyh,  post);
         
  // Flag the denominator error as used
  comD->SetBit(kUsedBit);
      }
      else if (!rel) {
  // We cannot cancel, and it is not relative, add to blob
  if (kVerbose & kRatio)
    ::Info("Ratio",
     "Common numerator systematic %s will be added to blob",
     comN->GetTitle());
  notUsed++; // Count how many we're not setting directly
  continue;
      }
      // Now we can define our common error. Which value we set depend
      // on whether the error was cancelled with the denominator
      Int_t        cId  = ret->DefineCommon(comN->GetTitle(),rel,
              eyl, eyh, kBox);
      HolderCommon* c   = ret->FindCommon(cId);
      c->CopyAttr(comN);
      // Flag the numerator error as used
      comN->SetBit(kUsedBit);
    }
    // Loop over all denominator common systematic errors.  If it is
    // not cancelled with the numerator and is relative, we define a
    // new error on the output ratio. Otherwise it will be absorbed in a blob 
    TIter nextDC(&den->fCommon);
    HolderCommon* comD = 0;
    while ((comD = static_cast<HolderCommon*>(nextDC()))) {
      if (comD->TestBit(kUsedBit)) continue;
      if (!comD->IsRelative()) {
  if (kVerbose & kRatio)
    ::Info("Ratio",
     "Common denominator systematic %s will be added to blob",
     comD->GetTitle());
  notUsed++; // Count how many we're not setting directly
  continue;
      }
      if (kVerbose & kRatio)
  ::Info("Ratio", "Propagating common sysmatic %s to ratio",
         comD->GetTitle());
      // ret->Print("sys");
      // Now we can define our common error. Which value we set depend
      // on whether the error was cancelled with the denominator
      Int_t        cId  = ret->DefineCommon(comD->GetTitle(),true,
              comD->GetYDown(1),
              comD->GetYUp(1),
              kBox);
      HolderCommon* c   = ret->FindCommon(cId);
      c->CopyAttr(comD);
      // Flag the denominator error as used
      comD->SetBit(kUsedBit);
      // ret->Print("sys");
    }
    
    // Loop over all point-to-point errors in the numerator, and see
    // if we have the same point-to-point error in the denominator.
    // If we do, we can (partially) cancel this point-to-point error.
    //
    // Technically, we mark each used point-to-point error as used.
    //
    // Point-to-point errors that are not shared by both the numerator
    // and the denominator will not be cancelled.  Instead, they will
    // be added in on the residual point-to-point sum (in quadrature)
    // error.
    //
    // We encode the IDs of the
    // errors in a single mask.
    //
    // 24 -------- 16 -------- 8 -------- 0
    //  | Denom ID  | Num ID   | Ratio ID |
    //  +-----------+----------+----------+
    //
    Int_t      sCnt = 0;
    TArrayI    shared(num->fPoint2Point.GetEntries() +
          den->fPoint2Point.GetEntries());
    TIter      nextNP(&num->fPoint2Point);
    HolderP2P* p2pN = 0;
    while ((p2pN = static_cast<HolderP2P*>(nextNP()))) {      
      Bool_t     rel  = p2pN->IsRelative();
      Int_t      idD  = den->FindId(p2pN->GetTitle());      
      HolderP2P* p2pD = idD == 0 ? 0     : den->FindP2P(idD);
      Bool_t     same = !p2pD    ? false : rel == p2pD->IsRelative();
      Int_t      mask = 0;
      // Printf("Id of denom %s -> %d (%p, %s, %s)", p2pN->GetTitle(), idD,
      //       p2pD, same ? "same" : "diff",
      //       p2pD ? (p2pD->TestBit(kUsedBit) ? "used" : "ready") : "?");
      if (cancc    &&
    p2pD     &&   // Check it is a p2p
    same) {   // Same as this
  // We found the same error in the denominator, so we need to
  // cancel errors when looping.  
  mask = ((p2pD->GetUniqueID() & 0xFF) << 16);
  // Printf("Found %s in both denom (%s,%d) and num (%s,%d)",
  //        p2pN->GetTitle(),
  //        p2pD->GetTitle(), p2pD->GetUniqueID(),
  //        p2pN->GetTitle(), p2pD->GetUniqueID());
  if (kVerbose & kRatio)
    ::Info("Ratio", "Cancelling the p2p systamtic error %s "
     "between numerator and denominator",
     p2pN->GetTitle());
  // Flag the denominator error as used
  p2pD->SetBit(kUsedBit);
      }
      Int_t      pId  = ret->DeclarePoint2Point(p2pN->GetTitle(),rel,
            kBox);
      HolderP2P* p    = ret->FindP2P(pId);
      p->CopyAttr(p2pN);
      // Flag the numerator error as used
      p2pN->SetBit(kUsedBit);
      mask |= (((p2pN->GetUniqueID() & 0xFF) << 8) |
         (pId & 0xFF));
      // Printf("Found shared P2P: %s -> %d",  p2pN->GetTitle(), pId);
      // num->Print("sys");
      // den->Print("sys");
      shared[sCnt] = mask;
      sCnt++;
    }
    // Loop over all denominator point-to-point systematic errors.  If
    // it is not cancelled with the numerator, we define a new error
    // on the output ratio
    TIter      nextDP(&den->fPoint2Point);
    HolderP2P* p2pD = 0;
    while ((p2pD = static_cast<HolderP2P*>(nextDP()))) {      
      if (p2pD->TestBit(kUsedBit)) {
  // Printf("P2P of denom %s already used", p2pD->GetTitle());
  continue;
      }
      Int_t      pId  = ret->DeclarePoint2Point(p2pD->GetTitle(),true,
            kBox);
      HolderP2P* p    = ret->FindP2P(pId);
      p->CopyAttr(p2pD);
      // Flag the numerator error as used
      p2pD->SetBit(kUsedBit);
      Int_t mask = (((p2pD->GetUniqueID() & 0xFF) << 16) |
        (pId & 0xFF));
      // Printf("P2P of %s (%d) denom propagated", p2pD->GetTitle(),
      //        p2pD->GetUniqueID());
      shared[sCnt] = mask;
      sCnt++;
    }    
    Int_t bId = 0;
    if (notUsed > 0) {
      // If we have some error that are not used - i.e., non-relative
      // common errors in the numerator or denominator, we need to
      // make a new point-to-point systematic error for those.
      bId = ret->DeclarePoint2Point("Blob", false, kBox);
    }

    // Now loop over all the points we have in the numerator and
    // figure out what to do for each point.
    Int_t cnt = 0; // Counter of output points
    for (Int_t i = 0; i < num->GetN(); i++) {
      Double_t x    = num->GetX(i);
      Double_t nY   = num->GetY(i);
      Double_t nEyl = 0;
      Double_t nEyh = 0;
      Double_t dY   = 0;
      Double_t dEyl = 0;
      Double_t dEyh = 0;
      Int_t    di1, di2;
      Int_t    di   = den->FindPoint(x, di1, di2, fac);
      Double_t daX  = 0; // Relative distance between interpolated points 
      if (di < -1) {
  if (kVerbose & kRatio)
    ::Warning("Ratio", "Next point %d (%f) not found in denominator",
        i, x);
  continue;
      }
      Bool_t cmn    = true;  // Include the common errors
      Bool_t stat   = false;  // Include the statistical errors
      Bool_t quad   = true;  // Add in quadrature
      Bool_t nosqrt = false; // Do not return squared error 
      if (di >= 0) {
  // Exact point
  // point,common,stat,quadratic,nosqrt,
  dY = den->GetYandError(di,cmn,stat,quad,nosqrt,dEyl,dEyh);
  di1 = di;
      }
      else {
  Double_t eyl1, eyl2, eyh1, eyh2;
  Double_t x1 = den->GetX(di1);
  Double_t x2 = den->GetX(di2);
  Double_t y1 = den->GetYandError(di1,cmn,stat,quad,nosqrt,eyl1,eyh1);
  Double_t y2 = den->GetYandError(di2,cmn,stat,quad,nosqrt,eyl2,eyh2);
  // Linear interpolation
  Double_t dx = (x2-x1);
  daX         = (x-x1)/dx;
  dY          = y1   + daX * (y2 - y1);
  dEyl        = eyl1 + daX * (eyl2 - eyl1);
  dEyh        = eyh1 + daX * (eyh2 - eyh1);
      }
      // At this point, we have the X and Y values, and the non-p2p
      // and non-cancelled and non-relative-common errors calculated
      // for numerator and denominator individually.
      //
      // We start by setting the point value.  Perhaps here, we should
      // also set the statistics.
      Double_t rY    = (dY != 0 ? nY/dY   : 0);
      ret->SetPoint(cnt, x, rY);
      ret->SetPointError(cnt, num->GetErrorXLeft(i), num->GetErrorXRight(i));

      // We have not dealt with the statistical errors yet.  We do
      // that here.
      Double_t sEyl = 0;
      Double_t sEyh = 0;
      if (TMath::Abs(rY) > 1.e-9) {
  Double_t snEyl = num->GetStatErrorDown(cnt);
  Double_t snEyh = num->GetStatErrorUp  (cnt);
  Double_t sdEyl = den->GetStatErrorDown(cnt);
  Double_t sdEyh = den->GetStatErrorUp  (cnt);
  // Statistical errors never cancel, and must be propagted. Note,
  // despite it technically being wrong, we treat the lower and
  // higher errors independently.
  sEyl = TMath::Sqrt(rY*rY*(snEyl*snEyl/nY/nY+sdEyl*sdEyl/dY/dY));
  sEyh = TMath::Sqrt(rY*rY*(snEyh*snEyh/nY/nY+sdEyh*sdEyh/dY/dY));
      }
      ret->SetStatError(cnt, sEyl, sEyh);
      
      // Next, we should update our blob systematic p2p error if it
      // exists.
      if (bId > 0) { 
  Double_t eyl = 0, eyh = 0;
  // Printf("Calculate blob syst.uncer. d=(%f,%f) n=(%f,%f)",
  //        dEyl, dEyh, nEyl, nEyh);
  if (denrel) {
    Double_t rdEyl = (dY != 0 ? dEyl/dY : 0);
    Double_t rdEyh = (dY != 0 ? dEyh/dY : 0);
    eyl            = rdEyl*rY;
    eyh            = rdEyh*rY;
  }
  else if (emax) {
    // Take the largest of the two errors, and divide by the
    // denominator. 
    eyl = TMath::Max(nEyl, dEyl);
    eyh = TMath::Max(nEyh, dEyh);
    eyl          /= dY;
    eyh          /= dY;
  }
  else if (cmin) {
    // Take the largest of the two errors, and divide by the
    // denominator. 
    eyl = TMath::Min(nEyl, dEyl);
    eyh = TMath::Min(nEyh, dEyh);
    eyl          /= dY;
    eyh          /= dY;
  }
  else {
    // Find the relative errors and add in quadrature 
    Double_t rnEyl = (nY != 0 ? nEyl/nY : 0);
    Double_t rnEyh = (nY != 0 ? nEyh/nY : 0);
    Double_t rdEyl = (dY != 0 ? dEyl/dY : 0);
    Double_t rdEyh = (dY != 0 ? dEyh/dY : 0);
    Double_t reyl  = TMath::Sqrt(rnEyl*rnEyl+rdEyl*rdEyl);
    Double_t reyh  = TMath::Sqrt(rnEyh*rnEyh+rdEyh*rdEyh);
    eyl            = reyl*rY;
    eyh            = reyh*rY;
  }
  // Printf("den err (%f%%,%f%%) - blob (%f%%,%f%%)", dEyl*100, dEyh*100,
  //        eyl*100,eyh*100);
  ret->SetSysError(bId, cnt, eyl, eyh);
      }

      // Now we need to update our various point-to-point
      // errors. We've stored a mask for each declared point-to-point
      // error.  Each mask consist of three fields, corresponding to
      // point-to-point errors in the denominator, numerator, and the
      // ratio.  If all three fields are set, we shold cancel between
      // numerator and denominator.  If only the ratio and either of
      // the denominator or numerator fields are set, we simply need
      // to copy that to the ratio point-to-point error - as a
      // relative error of course.
      for (Int_t im = 0; im < shared.GetSize(); im++) {
  Int_t mask = shared[im];
  if (mask <= 0) break;
  Int_t rId = ((mask >>  0) & 0xFF);
  Int_t nId = ((mask >>  8) & 0xFF);
  Int_t dId = ((mask >> 16) & 0xFF);
  if (nId == 0 && dId == 0) {
    ::Warning("Ratio", "Both numerator and denominator IDs are 0");
    break;
  }
  Bool_t   rel   = ret->IsRelative(rId);
  Bool_t   crel  = true; // rel
  Double_t pdEyl = 0;
  Double_t pdEyh = 0;
  Double_t pnEyl = 0;
  Double_t pnEyh = 0;
  if (dId != 0) {
    HolderP2P* pD  =  den->FindP2P(dId);
    Double_t   lfc = (crel ? (TMath::Abs(dY) > 1e-9 ? 1/dY : 0) : 1);
    pdEyl          =  pD->GetYDown(di1, di2, daX) * lfc;
    pdEyh          =  pD->GetYUp  (di1, di2, daX) * lfc;
  }
  if (nId != 0) {
    HolderP2P* pN  =  num->FindP2P(nId);
    Double_t   lfc = (crel ? (TMath::Abs(nY) > 1e-9 ? 1/nY : 0) : 1);
    pnEyl          =  pN->GetYDown(i) * lfc;
    pnEyh          =  pN->GetYUp  (i) * lfc;
  }
  // Printf("Doing error %d from n=%d (%f +%f -%f) d=%d (%f +%f -%f)",
  //        rId, nId, nY, pnEyl, pnEyh, dId, dY, pdEyh, pdEyl);
  Double_t peyl = 0;
  Double_t peyh = 0;
  Double_t lfc  = (!crel ? (TMath::Abs(dY) > 1e-9 ? 1/dY : 0) : 1);
  if (nId == 0) {
    // No numerator value - set directly
    peyl = pdEyl * lfc;
    peyh = pdEyh * lfc;
  } // else if
  else if (dId == 0) {
    // No denomenator value - set directly
    peyl = pnEyl * lfc;
    peyh = pnEyh * lfc;
  } // else if 
  else {
    // Cancel errors
    if (denrel) {
      peyl = pdEyl*lfc;
      peyh = pdEyh*lfc;
      // Printf("den err (%f%%,%f%%)", peyl*100,peyh*100);
    }
    else if (cmin) {
      peyl = TMath::Min(nEyl, dEyl) * lfc;
      peyh = TMath::Min(nEyh, dEyh) * lfc;
    }
    else if (cmax) {
      peyl = TMath::Min(nEyl, dEyl) * lfc;
      peyh = TMath::Min(nEyh, dEyh) * lfc;
    }
    else {
      peyl = TMath::Sqrt(TMath::Abs(pnEyl*pnEyl-pdEyl*pdEyl)) * lfc;
      peyh = TMath::Sqrt(TMath::Abs(pnEyh*pnEyh-pdEyh*pdEyh)) * lfc;
      // ::Info("Ratio", "sqrt(|%8f^2 - %8f^2|)=%8f", pnEyl,pdEyl,peyl);
    }
  } // Else
  // Printf("%f, den err (%f%%,%f%%) - %d (%f%%,%f%%)",
  //        rY, pdEyl*100, pdEyh*100, rId, peyl*100,peyh*100);
  lfc = (rel ? 1 : rY);
  ret->SetSysError(rId, cnt, 0, 0, lfc*peyl, lfc*peyh);
      }

      // Finally, increment counter
      cnt++;
    }

    // Reset the used case bits
    num->ClearUsed();
    den->ClearUsed();

    if (cnt <= 0) {
      ::Warning("", "No common points found");
      delete ret;
      ret = 0;
      return 0;
    }
    return ret;
  }
  void ClearUsed() const
  {
    // Reset the used case bits
    TObject* oe = 0;
    TIter nextC(&fCommon);
    while ((oe = nextC())) {
      oe->ResetBit(kUsedBit);
      oe->ResetBit(kOnlyWeightBit);
    }
    TIter nextP(&fPoint2Point);
    while ((oe = nextP())) {
      oe->ResetBit(kUsedBit);
      oe->ResetBit(kOnlyWeightBit);
    }
  }
  static Double_t KolomogorovTest(const GraphSysErr* g1,
          const GraphSysErr* g2)
  {
    Double_t z;
    return KolomogorovTest(g1, g2, z);
  }
  static Double_t KolomogorovTest(const GraphSysErr* g1,
          const GraphSysErr* g2,
          Double_t&          z)
  {
    TArrayD  a1y;
    TArrayD  a2y;
    TArrayD  a1e2;
    TArrayD  a2e2;
    Int_t    cnt   = CacheGraphs(g1, g2, a1y, a2y, a1e2, a2e2);
    if (cnt <= 0) return -1;

    Double_t s1    = a1y.GetSum();  // Summed content
    Double_t s2    = a2y.GetSum();  // Summed content
    Double_t sw1   = a1e2.GetSum(); // Summed weights
    Double_t sw2   = a2e2.GetSum(); // Summed weights
    if (a1e2.GetSum() <= 0 && a2e2.GetSum()) {
      ::Warning("ChisquareTest", "No errors");
      return -1;
    }

    Double_t sr1 = 0;
    Double_t sr2 = 0;
    Double_t dfMax = 0;
    for (Int_t i = 0; i < cnt; i++) {
      // KolmogorowvTest
      sr1 += a1y[i] / s1;
      sr2 += a2y[i] / s2;
      dfMax =  TMath::Max(dfMax, TMath::Abs(sr1-sr2));
    }
    // K-S test
    Double_t se1 = s1 * s1 / sw1;
    Double_t se2 = s2 * s2 / sw2;
    z   = dfMax * TMath::Sqrt(se1 * se2 / (se1 + se2));
    return TMath::KolmogorovProb(z);
  }
  static Int_t CacheGraphs(const GraphSysErr* g1,
         const GraphSysErr* g2,
         TArrayD&           a1y,
         TArrayD&           a2y,
         TArrayD&           a1e2,
         TArrayD&           a2e2)
  {
    Int_t    cnt    = 0;
    Int_t    n      = g1->GetN();
    a1y.Set(n);
    a2y.Set(n);
    a1e2.Set(n);
    a2e2.Set(n);
    for (Int_t i = 0; i < g1->GetN(); i++) {
      Double_t x    = 0;
      Double_t y1   = 0;
      Double_t e1yl = 0;
      Double_t e1yh = 0;
      Double_t e2yl = 0;
      Double_t e2yh = 0;
      Double_t y2   = 0;
      if (!NextPoint(i, g1, g2, x, y1, e1yl, e1yh, y2, e2yl, e2yh)) {
  ::Warning("ChisquareTest", "Next point - %d %f not found", i, x);
  continue;
      }

      Double_t e1y2   = TMath::Power(TMath::Max(e1yl,e1yh), 2);
      Double_t e2y2   = TMath::Power(TMath::Max(e2yl,e2yh), 2); 

      // Cache the numbers 
      a1y[cnt]  =  y1;   // Store value 
      a2y[cnt]  =  y2;   // Store value 
      a1e2[cnt] =  e1y2; // Store square error
      a2e2[cnt] =  e2y2; // Store square error 
      cnt++;
    }
    if (cnt <= 0) 
      ::Warning("CacheGraphs", "No common points found");
    
    return cnt;
  }
  static Double_t ChisquareTest(const GraphSysErr* g1,
        const GraphSysErr* g2,
        EChi2Type          type=kExperimentTruth)
  {
    Int_t ndf = 0;
    return ChisquareTest(g1, g2, ndf, type);
  }
  static Double_t ChisquareTest(const GraphSysErr* g1,
        const GraphSysErr* g2,
        Int_t&             ndf,
        EChi2Type          type)
  {
    Double_t chi2   = 0;
    ndf             = -1;
    TArrayD  a1y;
    TArrayD  a2y;
    TArrayD  a1e2;
    TArrayD  a2e2;
    Int_t    cnt   = CacheGraphs(g1, g2, a1y, a2y, a1e2, a2e2);
    if (cnt <= 0) return -1;

    Double_t s1    = a1y.GetSum();  // Summed content
    Double_t s2    = a2y.GetSum();  // Summed content
    if (type == kModelModel && (a1e2.GetSum() <= 0 && a2e2.GetSum())) {
      ::Warning("ChisquareTest", "No errors");
      return -1;
    }

    ndf = cnt;
    
    // Now loop over cached points
    Double_t s = s1 + s2;
    for (Int_t i = 0; i < ndf; i++) {
      switch (type) {
      case kExperimentExperiment:  {
  Double_t sum   = a1y[i]+a2y[i];
  Double_t delta = s2 * a1y[i] - s1 * a2y[i];
  chi2           += delta * delta / sum;
      }
  break;
      case kExperimentModel: {
  Double_t v1 =  s2 * a1y[i] - s1 * a2e2[i];
  Double_t v2 =  v1 * v1 + 4 * s2 * s2 * a1y[i] * a2e2[i];
  Double_t pp =  (v1 + v2) / (2 * s2 * s2);
  Double_t p1 =  pp * s1;
  Double_t p2 =  pp * s2;
  Double_t d1 =  a1y[i] - p1;
  Double_t d2 =  a1y[i] - p2;
  chi2        += d1 * d1 / p1;
  if (a2e2[i] > 0) chi2 += d2 * d2 / a2e2[i];
      }
  break;
      case kModelModel: {
  Double_t sigma =  s1 * s1 * a2e2[i] + s2 * s2 * a1e2[i];
  Double_t delta =  s2 * a1y[i] - s1 * a2y[i];
  chi2           += delta * delta / sigma;
      }
  break;
      case kExperimentTruth: {
  Double_t delta =  a2y[i] - a1y[i];
  chi2           += delta * delta / a1e2[i];
      }
  break;
      default:
  ::Warning("ChisquareTest", "Should not happen");
  return -1;
      }
  
    }
    if      (type == kExperimentExperiment) chi2 /= s1 * s2;
    else if (type == kExperimentTruth)      chi2 /= ndf;

    return chi2;
  }
  
  void RemovePoint(Int_t i)
  {
    if (i < 0 || i >= GetN()) return;
    fData->RemovePoint(i);
    TIter next(&fPoint2Point);
    HolderP2P* p2p = 0;
    while ((p2p = static_cast<HolderP2P*>(next()))) {
      p2p->fGraph->RemovePoint(i);
    }
  }
  void SwapPoints(Int_t i, Int_t j, Bool_t reflect=false)
  {
    if (i == j && !reflect) return;
    SwapPoints(fData, i, j, reflect);
    TIter next(&fPoint2Point);
    HolderP2P* p2p = 0;
    while ((p2p = static_cast<HolderP2P*>(next()))) {
      SwapPoints(p2p->fGraph, i, j, reflect);
    }
  }
  GraphSysErr* Reflect(Double_t x0=0) const
  {
    GraphSysErr* cpy = new GraphSysErr(*this);
    for (Int_t i = 0; i < GetN()/2; i++) {
      cpy->SwapPoints(i, GetN()-i-1, true);
    }
    return cpy;
  }
  Bool_t Symmetrize(GraphSysErr* other)
  {
    GraphSysErr*  cpy = new GraphSysErr(*other);

    // Declare "blob" systematic point-to-point error 
    Int_t         id  = DeclarePoint2Point("Syst.unc..", false, kBox);

    // Loop over common errors 
    HolderCommon* cmn = 0;
    TIter         cNext(&fCommon);
    while ((cmn = static_cast<HolderCommon*>(cNext()))) {
      Int_t oId = cpy->FindId(cmn->GetTitle());
      if (oId <= 0) {
  // If the common error ins't found in copy - it should,
  // because we did a straight copy - then give up.
  Error("Symmetrice", "Common error %s not found in %s",
        cmn->GetTitle(), cpy->GetName());
  cpy->Print("all");
  delete cpy;
  cpy = 0;
  break;
  // }
  // }
      } // oId <= 0

      // Zero this error here - we should add it later       
      HolderCommon* oCmn = cpy->FindCommon(oId);
      if (cmn->fEyl <= 0 && cmn->fEyh <= 0) {
  // Copy the error over 
  cmn->fEyl = oCmn->fEyl;
  cmn->fEyh = oCmn->fEyh;
      }
      oCmn->Set(0,0);
    }  // while cmn
    if (!cpy) return false;

    // Loop over all points 
    Int_t   cnt = 0;
    Int_t   n   = other->GetN();
    TArrayI used(n); // book-keeping
    Int_t   nonZero = 0; // how many non-zero sys.uncerr.
    used.Reset(-1);
    for (Int_t i = 0; i < n; i++) {
      if (used[i] >= 0 && used[i] != i) {
  Int_t j = used[i];
  // Info("Symmetrice", "Copy %d to %d instead of %d", j, cnt, i);
  SetPoint(cnt, -GetX(j), GetY(j));
  SetPointError(cnt, GetErrorXRight(j), GetErrorXLeft(j));
  SetStatError(cnt, GetStatErrorDown(j), GetStatErrorUp(j));
  SetSysError(id, cnt, 0, GetSysErrorY(id, j));
  // Info("Symmetrice", "%f -> %f +/- %f",
  //      -GetX(j), GetY(j), GetSysErrorY(id,j));
  cnt++;
  continue;
      }
      used[i] = i;
      Double_t eyl1, eyh1;
      Double_t exl1  = cpy->GetErrorXLeft(i);
      Double_t exh1  = cpy->GetErrorXRight(i);
      Double_t seyl1 = cpy->GetStatErrorDown(i);
      Double_t seyh1 = cpy->GetStatErrorUp(i);
      Double_t y1    = cpy->GetYandError(i,true,false,true,false,eyl1,eyh1);
      Double_t x1    = cpy->GetX(i);
      Double_t y2    = 0;
      Double_t eyl2  = 0;
      Double_t eyh2  = 0;
      Double_t seyl2 = 0;
      Double_t seyh2 = 0;
      Int_t    i1, i2;
      Int_t    j  = cpy->FindPoint(-x1, i1, i2); // Find mirror 
      if (j < -1) {
  // IF there's no mirror point, the just set to current 
  // Info("Symmetrice", "No mirror of %f using %d for %d", x1, i, cnt);
  SetPoint(cnt, x1, y1);
  SetPointError(cnt, exl1, exh1);
  SetStatError(cnt, seyl1, seyh1);
  SetSysError(id, cnt, 0, 0, eyl1, eyh1);
  cnt++;
  continue;
      }
      if (j >= 0) {
  // If we have exact mirror, then use that point and mark as
  // used
  // Info("Symmetrice", "Got exact mirror of %f(%d) at %d", x1, i, j);
  y2    = cpy->GetYandError(j,true,false,true,false,eyl2,eyh2);
  seyl2 = cpy->GetStatErrorDown(j);
  seyh2 = cpy->GetStatErrorUp(j);
  // We should copy values when we get to j
  used[j] = cnt;
      }
      else {
  // Otherwise we use interpolation between two points 
  // Info("Symmetrice", "Interpolate mirror of %f(%d) between %d-%d",
  //      x1, i, i1, i2);
  cpy->FindYandError(-x1,true,true,true,false,y2,eyl2,eyh2,seyl2,seyh2);
      }
      // Now calculate weighted mean and variance
      Double_t newY = (y1+y2) / 2;
      Double_t newV = 0;
      if ((eyh1+eyl1) > 1e-9 && (eyh2+eyl2) > 1e-9) {
  Double_t s1   = 2 * eyl1 * eyh1 / (eyh1 + eyl1);
  Double_t sp1  = (eyh1 - eyl1) / (eyh1 + eyl1);
  Double_t w1   = .5 * TMath::Power(s1+y1*sp1, 3) / s1;
  Double_t s2   = 2 * eyl2 * eyh2 / (eyh2 + eyl2);
  Double_t sp2  = (eyh2 - eyl2) / (eyh2 + eyl2);
  Double_t w2   = .5 * TMath::Power(s2+y2*sp2, 3) / s2;
  // Printf("s1=%f sp1=%f w1=%f s2=%f sp2=%f w2=%f",
  Double_t sumW = (w1+w2);
  newY          = (y1*w1+y2*w2) / sumW;
  newV          = TMath::Sqrt((w1*w1*s1*s1+w2*w2*s2*s2) / sumW / sumW);
  nonZero++;
      }
      Double_t seyl = TMath::Sqrt(seyl1*seyl1+seyl2*seyl2);
      Double_t seyh = TMath::Sqrt(seyh1*seyh1+seyh2*seyh2);
      
      // Info("Symmmetrice", "%f (%f +/- %f) + (%f +/- %f) -> %f +/- %f",
      //      x1, y1, w1, y2, w2, newY, newV);
      SetPoint(cnt, x1, newY);
      SetPointError(cnt, exl1, exh1);
      SetStatError(cnt, seyl, seyh);
      SetSysError(id, cnt, 0, newV);
      cnt++;
    }
    if (nonZero <= 0) RemoveSysError(id);
    // Info("", "After symmetrization:");
    // Print("all");
    return true;
  }
  /* @} */
  void SavePrimitive(std::ostream& out, Option_t* option="")
  {
    TString opt(option);
    opt.ToLower();
    Bool_t load = opt.Contains("load"); opt.ReplaceAll("load", "");
    TString funcName;
    TPRegexp regex("func=([a-zA-z][a-zA-Z0-9_]*)");
    TObjArray* toks = regex.MatchS(opt);
    if (toks) {
      if (toks->GetEntriesFast() > 1) 
  funcName = toks->At(1)->GetName();
      if (toks->GetEntriesFast() > 0) 
  opt.ReplaceAll(toks->At(0)->GetName(), "");
      delete toks;
    }
      
    
    // Save initialization 
    if (!funcName.IsNull()) 
      out << "TObject* " << funcName << "(Option_t* o=\"\")\n";
    out << "{\n";
    if (load)
      out << " // Load class\n"
    << "  if (!gROOT->GetClass(\"GraphSysErr\"))\n"
    << "    gROOT->LoadMacro(\"GraphSysErr.C+\");\n";
    out << "  GraphSysErr* g = new GraphSysErr(\""
  << GetName() << "\",\"" << GetTitle() << "\","
  << GetN() << ");\n";
    // Save attributes 
    out << "  // Point options\n"
  << "  g->SetMarkerStyle("  << GetMarkerStyle() << ");\n"
  << "  g->SetMarkerColor("  << GetMarkerColor() << ");\n"
  << "  g->SetMarkerSize("   << GetMarkerSize() << ");\n"
  << "  g->SetLineStyle("    << GetLineStyle() << ");\n"
  << "  g->SetLineColor("    << GetLineColor() << ");\n"
  << "  g->SetLineWidth("    << GetLineWidth() << ");\n"      
  << "  g->SetFillStyle("    << GetFillStyle() << ");\n"
  << "  g->SetFillColor("    << GetFillColor() << ");\n"
  << "  g->SetXTitle(\""     << fXTitle << "\");\n"
  << "  g->SetYTitle(\""     << fYTitle << "\");\n"
  << "  g->SetDataOption("   << fDataOption << ");\n"
  << "  // Sum options\n"
  << "  g->SetSumOption("    << fSumOption << ");\n"
  << "  g->SetSumTitle(\""   << fSumTitle << "\");\n"
  << "  g->SetSumLineStyle(" << fSumLine.GetLineStyle() << ");\n"
  << "  g->SetSumLineColor(" << fSumLine.GetLineColor() << ");\n"
  << "  g->SetSumLineWidth(" << fSumLine.GetLineWidth() << ");\n"      
  << "  g->SetSumFillStyle(" << fSumFill.GetFillStyle() << ");\n"
  << "  g->SetSumFillColor(" << fSumFill.GetFillColor() << ");\n"
  << "  g->SetCommonSumOption("    << fCommonSumOption << ");\n"
  << "  g->SetCommonSumTitle(\""   << fCommonSumTitle << "\");\n"
  << "  g->SetCommonSumLineStyle(" <<fCommonSumLine.GetLineStyle()<<");\n"
  << "  g->SetCommonSumLineColor(" <<fCommonSumLine.GetLineColor()<<");\n"
  << "  g->SetCommonSumLineWidth(" <<fCommonSumLine.GetLineWidth()<<");\n"
  << "  g->SetCommonSumFillStyle(" <<fCommonSumFill.GetFillStyle()<<");\n"
  << "  g->SetCommonSumFillColor(" <<fCommonSumFill.GetFillColor()<<");\n"
  << "  // Stat options\n"
  << "  g->SetStatRelative(" << fStatRelative << ");\n";
    TIter nextC(&fCommon);
    HolderCommon* cmn = 0;
    while ((cmn = static_cast<HolderCommon*>(nextC())))
      cmn->SavePrimitive(out, "d");
    TIter nextP(&fPoint2Point);
    HolderP2P* p2p = 0;
    while ((p2p = static_cast<HolderP2P*>(nextP())))
      p2p->SavePrimitive(out, "d");
    Int_t n = GetN();
    out << " // " << n << " points\n";
    Bool_t statRel = IsStatRelative();
    for (Int_t i = 0; i < n; i++) {
      Double_t y = GetY(i);
      out << "  g->SetPoint(" << i << ',' << GetX(i) << ',' << y << ");\n"
    << "  g->SetPointError(" << i << ',' << GetErrorXLeft(i) << ','
    << GetErrorXRight(i) << ");\n"
    << "  g->SetStatError(" << i << ','
    << (statRel ? 1/y : 1)*GetStatErrorDown(i) << ','
    << (statRel ? 1/y : 1)*GetStatErrorUp(i) << ");\n";
      nextP.Reset();
      while ((p2p = static_cast<HolderP2P*>(nextP()))) {
  Int_t  id  = p2p->GetUniqueID();
  Bool_t rel = p2p->IsRelative();
  out << "  g->SetSysError(" << id << ',' << i << ','
      << GetSysErrorXLeft(id, i) << ','
      << GetSysErrorXRight(id, i) << ','
      << (rel ? 1/y : 1) * GetSysErrorYDown(id, i) << ','
      << (rel ? 1/y : 1) * GetSysErrorYUp(id, i) << ");\n";
      } // while(p2p)
    } // for (i)
    out << "  if (o && o[0] != '\\0') {\n"
  << "    g->Draw(o);\n";
    if (fDrawn && fDrawn->GetHistogram())
      out << "    if (g->GetMulti() && g->GetMulti()->GetHistogram()) {\n"
    << "      g->GetMulti()->GetHistogram()->SetMinimum("
    << fDrawn->GetHistogram()->GetMinimum() << ");\n"
    << "      g->GetMulti()->GetHistogram()->SetMaximum("
    << fDrawn->GetHistogram()->GetMaximum() << ");"
    << "    }\n";
   
    out << "  }\n";
    if (!funcName.IsNull()) out << "  return g;\n";
    out << "};\n";
  }
  void Save(const char* fileName)
  {
    std::ofstream out(fileName);
    TString funcName(fileName);
    funcName.ReplaceAll(".C","");
    out << "// \n"
  << "// Generated by GraphSysErr.C\n"
  << "// \n"
      // << "class GraphSysErr;\n\n";
  << "\n";
    SavePrimitive(out, Form("load func=%s", funcName.Data()));
    out.close();
  }
  void Export(std::ostream& out=std::cout,
        Option_t*     option="",
        Int_t         nsign=2)
  {
    // Info("Export", "Using %d significant digits on errors", nsign);
    TString opt(option);
    opt.ToLower();
    Bool_t header   = opt.Contains("h");
    Bool_t sysNames = opt.Contains("s");
    Bool_t comment  = opt.Contains("c");
    Bool_t title    = opt.Contains("t");
    
    ExportHeader(out, header, comment, nsign);
    
    // --- Export qualifiers -----------------------------------------
    Bool_t hasTitle = false;
    if (fQualifiers) {
      TIter nextQ(fQualifiers);
      TObject* q = 0;
      while ((q = nextQ())) {
  TString k(q->GetName());
  if (k.EqualTo("RE") || k.EqualTo("title", TString::kIgnoreCase))
    hasTitle = true;
  out << FormatKey("qual") << q->GetName() << " : "
      << q->GetTitle() << std::endl;
      }
    }
    if (!hasTitle && title)
      out << FormatKey("qual") << "RE : " << GetTitle() << std::endl;
    
    // --- Export X/Y titles ----------------------------------------
    const char* fill = "<please fill in>";
    TString xTit = fXTitle; EscapeLtx(xTit, "fill");
    TString yTit = fYTitle; EscapeLtx(yTit, "fill");
    
    out << FormatKey("xheader") << xTit << "\n"
  << FormatKey("yheader") << yTit << std::endl;

    // --- Export common errors --------------------------------------
    TIter nextC(&fCommon);
    HolderCommon* holderCommon = 0;
    while ((holderCommon = static_cast<HolderCommon*>(nextC()))) { 
      Bool_t rel = holderCommon->IsRelative();
      Double_t up = holderCommon->GetYUp(rel ? 100 : 1);
      Double_t down = holderCommon->GetYDown(rel ? 100 : 1);
      out << FormatKey("dserror");
      ExportError(out, down, up, true, rel, nsign);
      out << ":" << holderCommon->GetTitle() << std::endl;
    }

    // --- Export data points ----------------------------------------
    out  << FormatKey("data") << " x : y" << std::endl;
    Int_t n = GetN();
    for (Int_t i = 0; i < n; i++) {
      ExportPoint(out, i, true, sysNames, nsign);
      out << std::endl;
    }
    out << "*dataend:\n" 
  << "# End of dataset\n" << std::endl;
  } //*MENU*
  static void EscapeLtx(TString& val, const TString& fill="")
  {
    if (val.IsNull()) {val = fill; return; }
    if (!val.Contains("#") && !val.Contains("\\")) return;

    if (val[0]             !='$') val.Prepend("$");
    if (val[val.Length()-1]!='$') val.Append("$");
    val.ReplaceAll("#", "\\");
  }
  static void Export(const TSeqCollection* col,
         std::ostream& out,
         Option_t* option="H",
         Int_t nsign=0)
  {
    if (col->GetEntries() < 1) return;

    // --- Deduce options --------------------------------------------
    TString opt(option);
    opt.ToLower();
    Bool_t alsoTop = opt.Contains("h");
    Bool_t alsoCmt = opt.Contains("c");
    Bool_t alsoNme = opt.Contains("s");
    Bool_t alsoTit = opt.Contains("t");

    // --- some variables to use -------------------------------------
    const Double_t tol = 1e-10;
    GraphSysErr*   first = 0;
    GraphSysErr*   gse = 0;
    TList          toExport;
    TList          commons;
    TList          quals;
    TString        data("x ");
    Int_t          nPoints = -1;
    Int_t          idx     = 0;

    // --- Check the passed graphs for compatiblity ------------------
    TIter          nextCheck(col);
    TObject*       o       = 0;
    while ((o = nextCheck())) {
      if (!o->IsA()->InheritsFrom(GraphSysErr::Class())) continue;
      gse = static_cast<GraphSysErr*>(o);
      if (!first) {
  first   = gse;
  nPoints = first->GetN();
      }
      else {
  // --- Check number of points --------------------------------
  if (gse->GetN() != nPoints) {
    Int_t nTmp = TMath::Min(gse->GetN(), nPoints);
    ::Warning("Export", "Incompatible number of points %d in %s"
        "using only %d of them", 
        gse->GetN(), gse->GetName(), nTmp);
    nPoints = nTmp;
  }
  // --- check X values ----------------------------------------
  Bool_t ok = true;
  for (Int_t i = 0; i < nPoints; i++) {
    Double_t x1    = first->GetX(i);
    Double_t exl1  = first->fData->GetErrorXlow(i);
    Double_t exh1  = first->fData->GetErrorXhigh(i);
    Double_t xT    = gse->GetX(i);
    Double_t exlT  = gse->fData->GetErrorXlow(i);
    Double_t exhT  = gse->fData->GetErrorXhigh(i);
    if (TMath::Abs(x1-xT) > tol ||
        (exl1 > tol && TMath::Abs(exl1-exlT) > tol) ||
        (exh1 > tol && TMath::Abs(exh1-exhT) > tol)) {
      ::Warning("Export", "X--coordinate of %s @ point %d: %f (+%f,-%f) "
          "incompatible [%f (+%f,-%f)]", 
          gse->GetTitle(), i, xT, exhT, exlT, x1, exh1, exl1);
      ok = false;
      break;
    }
  } // for i      
  if (!ok) continue;
      } // !first

  // --- Get all possible qualifiers -----------------------------
      toExport.Add(gse);
      TIter    nextQ(gse->fQualifiers);
      TObject* q = 0;
      while ((q = nextQ())) {
  // ::Info("", "qualifier %s=%s", q->GetName(), q->GetTitle());
  StoreQual(quals, idx, q);
      }

      // --- Get all common systematics ------------------------------
      TIter    nextCmn(&(gse->fCommon));
      TObject* oCmn = 0;
      while ((oCmn = nextCmn())) {
  if (commons.FindObject(oCmn->GetTitle())) continue;
  TObjString* cmn = new TObjString(oCmn->GetTitle());
  if (gse == first) cmn->SetUniqueID(gse->FindId(oCmn->GetTitle()));
  commons.Add(cmn);
      }
      
      idx++;
      data += ": y ";
    }
    // --- Now deduce the common common errors -----------------------
    TIter    nextCmn(&(commons));
    TObject* oCmn = 0;
    while ((oCmn = nextCmn())) {
      TString  oNme(oCmn->GetName());
      Bool_t   found = true;
      Double_t ecl1  = -1; 
      Double_t ech1  = -1;
      
      // --- Loop over data ------------------------------------------
      TIter  nextG(&toExport);
      while ((gse = static_cast<GraphSysErr*>(nextG()))) {
  Int_t gseId = first->FindId(oNme);

  if (gseId > 0) {
    // If we do, check compatibility  
    if (ecl1 < 0) {
      // First time we get this error 
      ecl1 = gse->GetCommonErrorYDown(gseId);
      ech1 = gse->GetCommonErrorYUp(gseId);
      continue;
    }
    else {
      // Now check values within tolerance 
      Double_t eclT  = gse->GetCommonErrorYDown(gseId);
      Double_t echT  = gse->GetCommonErrorYUp(gseId);
      if ((ecl1 > tol && TMath::Abs(ecl1-eclT) > tol) ||
    (ech1 > tol && TMath::Abs(ech1-echT) > tol)) {
        found = false;
      } // incompatible
    } // ecl1 >= 0
  } // gseId > 0
  else 
    // This graph does not have this common, flag it 
    found = false;
      } // while(gse)
#if INCOMPAT_CMN_AS_QUAL
#else 
      found = true;
#endif 
      oCmn->SetBit(BIT(14), found);

      // If we found the error in all graphs and they are all
      // compatible, then all is good,
      if (found) continue;

      // If not, we represent the common systematic as a qualifier 
      idx   = 0;      
      nextG.Reset();
      while ((gse = static_cast<GraphSysErr*>(nextG()))) {
  Int_t gseId = first->FindId(oNme);
  TString val;
  if (gseId > 0) {
    // This graph has the value, store as qual
    Bool_t   rel = gse->IsRelative(gseId);
    Double_t ecl = gse->GetCommonErrorYDown(gseId);
    Double_t ech = gse->GetCommonErrorYUp(gseId);
    std::stringstream s;
    ExportError(s, ecl, ech, false, rel, nsign);
    val = s.str().c_str();
  }
  else {
    // This does not have the error, store empty string
    val = "";
  }
  // Now store this in our qualifier table 
  StoreQual(quals, idx, oNme, val);
  idx++;
      } // while(gse)
    } // while common

      // --- Sanity checks ---------------------------------------------
    if (nPoints <= 0) {
      ::Error("Export", "No points to write");
      return;
    }
    if (toExport.GetEntries() <= 0) {
      ::Error("Export", "No graphs to export");
      return;
    }
    if (!first) {
      ::Error("Export", "Didn't get the first graph");
      return;
    }

    // --- Export header --------------------------------------------
    first->ExportHeader(out, alsoTop, alsoCmt, nsign);

    // --- Export qualifiers -----------------------------------------
    // quals.ls();
    Bool_t hasTitle = false;
    TIter nextQ(&quals);
    TList* ql = 0;
    while ((ql = static_cast<TList*>(nextQ()))) {
      TString k(ql->GetName());
      if (k.EqualTo("RE") || k.EqualTo("title", TString::kIgnoreCase))
  hasTitle = true;
      out << FormatKey("qual") << ql->GetName();
      for (Int_t i = 0; i < toExport.GetEntries(); i++) {
  TObject* qv = ql->At(i);
  TString  v  = "";
  if (qv)  v  = qv->GetName();
  out << " : " << v;
      }
      out << std::endl;
    }
    if (!hasTitle && alsoTit) {
      out << FormatKey("qual") << "RE ";
      for (Int_t i = 0; i < toExport.GetEntries(); i++) 
  out << ": " << toExport.At(i)->GetTitle();
      out << std::endl;
    }

    
    // --- Export axis titles ----------------------------------------
    const char*  fill     = "<please fill in>";
    const char*  fields[] = { "xheader", "yheader", 0 };
    const char** pfld     = fields;
    while (*pfld) { 
      out << FormatKey(*pfld);
      TIter nextSpec(&toExport);
      Bool_t one = true;
      while ((gse = static_cast<GraphSysErr*>(nextSpec()))) {
  TString val;
  if      ((*pfld)[0] == 'x') val = gse->fXTitle;
  else if ((*pfld)[0] == 'y') val = gse->fYTitle;
  EscapeLtx(val, fill);
  out << (one ? "" : ":") << val;
  if ((*pfld)[0] == 'x') break;
  one = false;
      }
      pfld++;
      out << std::endl;
    }
      
    // --- Export common systematics ---------------------------------
    TIter    nextC(&commons);
    while ((oCmn = nextC())) {
      TString tit(oCmn->GetName());
      EscapeLtx(tit);
      out << FormatKey("dserror");
#if INCOMPAT_CMN_AS_QUAL
      if (!oCmn->TestBit(BIT(14))) {
  ::Warning("Export",
      "Common systematic error \"%s\" represented by qual",
      oCmn->GetName());
  out << "- : " << tit << std::endl;
  continue;
      }
      UInt_t   id  = first->FindId(oCmn->GetName());
      Bool_t   rel = first->IsRelative(id);
      Double_t ecl = first->GetCommonErrorYDown(id);
      Double_t ech = first->GetCommonErrorYUp(id);
      ExportError(out, ecl, ech, true, rel, nsign);
      out << " : "<< tit << std::endl;
#else 
      out << tit << std::flush;
      TIter  next(&toExport);
      while ((gse = static_cast<GraphSysErr*>(next()))) {
  UInt_t   id  = gse->FindId(oCmn->GetName());
  Bool_t   rel = gse->IsRelative(id);
  Double_t ecl = gse->GetCommonErrorYDown(id);
  Double_t ech = gse->GetCommonErrorYUp(id);
  out << " : ";
  ExportError(out, ecl, ech, true, rel, nsign);
      }
      out << ":" << tit << std::endl;
#endif
    }

    // --- Export points ---------------------------------------------
    out << FormatKey("data") << data << std::endl;
    for (Int_t i = 0; i < nPoints; i++) {
      TIter  next(&toExport);
      Bool_t one = true;
      while ((gse = static_cast<GraphSysErr*>(next()))) {
  gse->ExportPoint(out, i, one, alsoNme, nsign);
  one = false;
      }
      out << std::endl;
    }
    out << "*dataend:\n" 
  << "# End of dataset\n" << std::endl;
    // return out;
  }
  static TSeqCollection* Import(const TString& fileName)
  {
    std::ifstream in(fileName.Data());
    if (!in) {
      ::Warning("Import", "Failed to open \"%s\"", fileName.Data());
      return 0;
    }
    TList* ret = new TList;
    ret->SetOwner();
    GraphSysErr* g     = 0;
    GraphSysErr* first = 0;
    Int_t id = 0;
    do {
      Int_t    sub  = 1;
      UShort_t nIdx = 256;
      int      cur  = in.tellg();
      do {
  in.seekg(cur, in.beg);
  if (!(g = Import(in, sub, &nIdx))) break;
  if (kVerbose & kImport)
    ::Info("Import", "Imported %d of %d", sub, nIdx);
  g->SetName(Form("ds_%d_%d", id, sub-1));
  if (!first) first = g;
  else        g->CopyKeys(first);
  ret->Add(g);
  sub++;
      } while(sub < nIdx);
      id++;
    } while (!in.eof());
    return ret;
  }
  static GraphSysErr* Import(std::istream& in, UShort_t idx=0,
           UShort_t* nIdx=0)
  {
    GraphSysErr* ret    = 0;
    Int_t        n      = 0;
    Bool_t       inSet  = false;
    Bool_t       inData = false;
    TString      xtit   = "";
    TString      ytit   = "";
    TString      tmp    = "";
    TString      last   = "";
    TList        quals;
    TList        keys;
    Int_t        isty   = 0;
    do {
      TString line;
      line.ReadLine(in);
      tmp = line.Strip(TString::kBoth, ' '); line = tmp;
      if (line.IsNull()) continue; 
      if (line[0] == '#') continue;
      if (line[0] == '*') { // We have some sort of header stuff 
  Int_t colon = line.Index(":"); 
  if (colon == kNPOS) continue;
  Int_t   len   = line.Length()-colon-1;
  TString value = line(colon+1,len);
  tmp = value.Strip(TString::kBoth, ' '); value = tmp;
  last = line(1,colon-1);
  if (kVerbose & kImport) 
    ::Info("Import", "Got a key '%s' -> '%s'", last.Data(), value.Data());

  if (last.BeginsWith("dataset", TString::kIgnoreCase)) {
    inSet = true;
    inData = false;
    isty  = 0;
  }
  else if (last.BeginsWith("dataend", TString::kIgnoreCase)) {
    inSet = false;
    inData = false;
    // Always terminate on the end of a data set
    break;
  }
  // else if (last.BeginsWith("dscomment",TString::kIgnoreCase)) {
  // tit = value;
  // }
  else if (last.BeginsWith("qual", TString::kIgnoreCase)) {
    // ::Info("", "Got qual: %s", value.Data());
    // if (!qual.IsNull()) {
    //   ::Info("", "Adding seen qual: %s", qual.Data());
    quals.Add(new TObjString(value));
  }
  else if (last.BeginsWith("xheader", TString::kIgnoreCase))
    xtit = value;
  else if (last.BeginsWith("yheader", TString::kIgnoreCase))
    ytit = value;
  else if (last.BeginsWith("dserror", TString::kIgnoreCase)) {
    // Common systematic error 
    if (kVerbose & kImport) 
      ::Info("Import", "Got common error line: '%s'", line.Data());
    TObjArray* tokens = value.Tokenize(":");
    Double_t el = 0, eh = 0;
    Bool_t   rel = false;
    if (ImportError(Token(tokens, 0), el, eh, rel)) { 
      if (rel) { el /= 100.; eh /= 100.; }
      TString& tnam = Token(tokens, 1);
      TString  nam  = tnam.Strip(TString::kBoth, ' ');
      if (!ret) ret = new GraphSysErr("imported", "");
      Int_t id = ret->DefineCommon(nam, rel, el, eh, kRect);
      ret->SetSysLineColor(id, (isty % 6) + 2);
      ret->SetSysFillColor(id, (isty % 6) + 2);
      ret->SetSysFillStyle(id, 3001 + (isty % 10));
      isty++;
    }
  }
  else if (last.BeginsWith("data", TString::kIgnoreCase)) {
    if (kVerbose & kImport) 
      ::Info("Import", "Got start of data line: '%s'", line.Data());
    // These are the field we can deal with 
    // Let's get the header field value 
    TObjArray* tokens = value.Tokenize(":");
    if (nIdx) {
      *nIdx = tokens->GetEntriesFast();
      // ::Info("Import", "Max index is %d", *nIdx);
    }
    if (tokens->GetEntriesFast() > 2 && idx < 1) { 
      idx = 1;
      ::Warning("Import", "Can only import one data set at a time, "
          "selecting the %d column", idx);
    }
    else if (idx >= 1 && idx >= tokens->GetEntriesFast()) { 
      ::Warning("Import", "column %d not available for this data set",
          idx);
      inSet = false;
    }
    // We do nothing to the tokens here. 
    tokens->Delete();
    inData = true;
  }
  else {   
    // some other key. 
    if (kVerbose & kImport) 
      ::Info("Import", "Got pair line '%s' (%s -> %s)", line.Data(), 
       last.Data(), value.Data());
    keys.Add(new TNamed(last, value));
  }
  continue; 
      }
      if (!inData) { 
  if (kVerbose & kImport) 
    ::Info("Import", "Got a possible contiuation line '%s'", line.Data());
  // Probably a continuation line. 
  if (last.IsNull()) 
    // No last key or no map 
    continue; 

  if (kVerbose & kImport) 
    ::Info("Import", "Got some other line: '%s'", line.Data());

  TNamed* ptr = static_cast<TNamed*>(keys.FindObject(last));
  if (!ptr) 
    // Key is not added, do nothing to this line 
    continue; 

  TString tt(ptr->GetTitle());
  if (!tt.EndsWith(" ") && !tt.EndsWith("\t") && !tt.EndsWith("\n")) 
    // If value does not end in a white-space, add it. 
    tt.Append("\n");

  // Append our line, and set the title
  tt.Append(line);
  ptr->SetTitle(tt);

  continue;
      }
      if (!inSet) { 
  // We're outside a data-set so do nothing
  continue;
      }
      if (idx < 1) idx = 1;

      if (kVerbose & kImport) 
  ::Info("Import", "Got a data line: '%s'", line.Data());
      TObjArray* tokens = line.Tokenize(";");
      if (tokens->GetEntriesFast() < idx+1) {
  ::Warning("Import", "Too few columns %d<%0d in line: %s", 
      tokens->GetEntriesFast(), idx+1, line.Data());
  tokens->Delete();
  continue;
      }
      TString& xCol = static_cast<TObjString*>(tokens->At(0))->String();
      TString& yCol = static_cast<TObjString*>(tokens->At(idx))->String();
      if (kVerbose & kImport) 
  ::Info("Import", "xColumn: '%s', yColumn: '%s'", 
         xCol.Data(), yCol.Data());

      tmp = xCol.Strip(TString::kBoth, ' '); xCol = tmp;
      tmp = yCol.Strip(TString::kBoth, ' '); yCol = tmp;
      TString yFull = yCol;

      Int_t chop = yCol.Last('(');
      if (chop != kNPOS) {
  yCol.Remove(chop, yCol.Length()-chop);
  yCol.Remove(TString::kBoth, ' ');
      }
      if (kVerbose & kImport) 
  ::Info("Import", "xColumn: '%s', yColumn: '%s'", 
         xCol.Data(), yCol.Data());

      if (xCol.IsNull()) { 
  ::Warning("Import", "Empty X column in line: %s", line.Data());
  tokens->Delete();
  continue;
      }
      if (yCol.IsNull()) { 
  ::Warning("Import", "Empty Y column in line: %s", line.Data());
  tokens->Delete();
  continue;
      }
      Bool_t rel = false;
      Double_t x = 0, exl = 0, exh = 0;
      if (!ImportPoint(xCol, x, exl, exh, rel)) { 
  ::Warning("Import", "Failed to import X value in line: %s",line.Data());
  tokens->Delete();
  continue;
      }
      Double_t y = 0, eyl = 0, eyh = 0;
      if (!ImportPoint(yCol, y, eyl, eyh, rel)) {
  ::Warning("Import", "Failed to import X value in line: %s",line.Data());
  tokens->Delete();
  continue;
      }
      if (!ret) ret = new GraphSysErr("imported", "");
      if (rel) ret->SetStatRelative(true);
      ret->SetPoint(n, x, y);
      ret->SetPointError(n, exl, exh);

      if (ret->IsStatRelative()) { eyl /= 100; eyh /= 100; }
      ret->SetStatError(n, eyl, eyh);

      Int_t lparen = yFull.Index("(");
      Int_t rparen = yFull.Index(")");
      if (lparen == rparen) { 
  n++;
  continue;
      }
      TString    rem  = yFull(lparen+1, rparen-lparen-1);
      if (kVerbose & kImport) 
  ::Info("Import", "Got systematic errors: '%s'", rem.Data());

      TObjArray* stok = rem.Tokenize("D");
      Int_t      isys = 0;
      for (Int_t i = 0; i < stok->GetEntriesFast(); i++) {
  TString t = Token(stok,i);
  if (!t.BeginsWith("SYS=",TString::kIgnoreCase)) { 
    ::Warning("Import", "Ignoring unknown spec %s (@ %d) in %s",
        t.Data(), i, rem.Data());
    continue;
  }
  t.Remove(0,4);
  if (t[t.Length()-1] == ',') t.Remove(t.Length()-1,1);
  tmp = t.Strip(TString::kBoth, ' '); t = tmp;

  Double_t el = 0, eh = 0;
  rel = false;
  if (!ImportError(t, el, eh, rel)) continue;

  Int_t colon = t.Index(":");
  TString nam(Form("sys%d", isys));
  if (colon != kNPOS) nam = t(colon+1, t.Length()-colon-1);
    
  UInt_t id = ret->FindId(nam);
  if (id == 0) { 
    id = ret->DeclarePoint2Point(nam, rel, kRect);
    ret->SetSysLineColor(id, (isty % 6) + 2);
    ret->SetSysFillColor(id, (isty % 6) + 2);
    ret->SetSysFillStyle(id, 3001 + (isty % 10));
    isty++;
  }
  HolderP2P* p = ret->FindP2P(id);
  if (p->IsRelative()) { el /= 100.; eh /= 100.; }
  ret->SetSysError(id, n, exl, exh, el, eh);
  isys++;
      }
      stok->Delete();

      n++;      
    } while (!in.eof());
    if (ret) {
      TString rxtit = ExtractField(xtit, idx-1);
      TString rytit = ExtractField(ytit, idx-1);
      if (!rxtit.IsNull()) ret->SetXTitle(rxtit);
      if (!rytit.IsNull()) ret->SetYTitle(rytit);

      Bool_t hasTitle = false;
      TIter nextQ(&quals);
      TObject* qual = 0;
      while ((qual = nextQ())) {
  TString k = ExtractField(qual->GetName(), 0);
  TString q = ExtractField(qual->GetName(), idx);
  /*
    ::Info("LoopQ", "Qualifier string: %s -> %s,%s",
    qual->GetName(), k.Data(), q.Data());
  */
  ret->AddQualifier(k, q);
  if (!hasTitle &&
      (k.EqualTo("RE") ||
       k.EqualTo("title", TString::kIgnoreCase))) {
    hasTitle = true;
    ret->SetTitle(q);
  }
      }
      TIter nextK(&keys);
      TObject* pair = 0;
      while ((pair = nextK())) {
  ret->SetKey(pair->GetName(),pair->GetTitle());
  TString k = pair->GetName();
  if (!hasTitle && (k.EqualTo("location"))) {
    ret->SetTitle(pair->GetTitle());
    hasTitle = true;
  }
      }
      if (!hasTitle) ret->SetTitle("Title");
    }
    keys.SetOwner();
    keys.Delete();
    return ret;
  }
  /* @} */
  //__________________________________________________________________
  UInt_t DefineCommon(const char* title, Bool_t relative, 
          Double_t ey, EDrawOption_t option=kFill)
  {
    return DefineCommon(title, relative, ey, ey, option);
  }
  UInt_t DefineCommon(const char* title, Bool_t relative, 
          Double_t eyl, Double_t eyh, 
          EDrawOption_t option=kRect)
  {
    UInt_t  id = ++fCounter;
    TString name(Form("common_%08x", id));

    HolderCommon* h = new HolderCommon(name, title, relative, option, id);
    h->Set(eyl, eyh);
    
    fCommon.AddLast(h);
    return id;
  }
  UInt_t DeclarePoint2Point(const char* title, Bool_t relative, 
          EDrawOption_t option=kBar)
  {
    UInt_t  id = ++fCounter; 
    TString name(Form("p2p_%08x", id));

    Holder* h = new HolderP2P(name, title, relative, option, id);

    fPoint2Point.AddLast(h);
    return id;
  }
  void RemoveSysError(Int_t id)
  {
    HolderCommon* h = FindCommon(id);
    if (h) {
      fCommon.Remove(h);
      return;
    }
    HolderP2P* p = FindP2P(id);
    if (p) {
      fPoint2Point.Remove(p);
      return;
    }
  }
  UInt_t FindId(const char* title) const
  { 
    TString tit(title);

    TIter nextC(&fCommon);
    Holder* holder = 0;
    while ((holder = static_cast<Holder*>(nextC()))) {
      if (tit.EqualTo(holder->GetTitle())) return holder->GetUniqueID();
    }

    TIter nextP(&fPoint2Point);
    while ((holder = static_cast<Holder*>(nextP()))) {
      if (tit.EqualTo(holder->GetTitle())) return holder->GetUniqueID();
    }
    // fCommon.ls();
    // fPoint2Point.ls();
    return 0;
  }
  Int_t GetNSys() const { return fCounter; }
  /* @} */

  //__________________________________________________________________
  void SetPoint(Int_t i, Double_t x, Double_t y) 
  {
    if (!fData) return;
    fData->SetPoint(i, x, y);
  }
  void SetPointError(Int_t i, Double_t ex)
  {
    SetPointError(i, ex, ex);
  }
  void SetPointError(Int_t i, Double_t exl, Double_t exh)
  {
    if (!fData) return;
    fData->SetPointError(i, exl, exh, 
       fData->GetErrorYlow(i),
       fData->GetErrorYhigh(i));
  }
  void SetStatRelative(Bool_t rel) { fStatRelative = rel; }
  Bool_t IsStatRelative() const { return fStatRelative; }
  void SetStatError(Int_t i, Double_t ey) 
  {
    SetStatError(i, ey, ey);
  }
  /*
   * Set the statistical error on the ith data point
   * 
   * @param i   Point number 
   * @param eyl Lower error on Y
   * @param eyh Higher error on Y
   */
  void SetStatError(Int_t i, Double_t eyl, Double_t eyh) 
  {
    if (!fData) return;
    if (fStatRelative) { 
      eyl *= fData->GetY()[i];
      eyh *= fData->GetY()[i];
    }
    fData->SetPointError(i,
       fData->GetErrorXlow(i),
       fData->GetErrorXhigh(i),
       eyl, eyh);
  }
  void SetSysError(Int_t id, Double_t eyl, Double_t eyh)
  {
    HolderCommon* h = FindCommon(id);
    if (!h) return;
    h->Set(eyl, eyh);
  }
  void SetSysError(Int_t id, Int_t i, Double_t ex, Double_t ey)
  {
    HolderP2P* h = FindP2P(id);
    if (!h) return;
    h->Set(i, fData, ex, ey);
  }
  void SetSysError(Int_t id, Int_t i, Double_t exl, Double_t exh,
       Double_t eyl, Double_t eyh) 
  {
    HolderP2P* h = FindP2P(id);
    if (!h) return;
    h->Set(i, fData, exl, exh, eyl, eyh);
  }
  /* @} */
  //__________________________________________________________________
  void SetTitle(const char* name)
  {
    TNamed::SetTitle(name);
    if (fData) fData->SetTitle(name);
  }
  void SetDataOption(EDrawOption_t opt) { fDataOption = opt; } //*MENU*
  void SetXTitle(const char* title) { fXTitle = title; } //*MENU*
  void SetYTitle(const char* title) { fYTitle = title; } //*MENU*
  /* @} */

  //__________________________________________________________________
  Int_t    GetN() const { return fData ? fData->GetN() : 0; }
  Double_t GetX(Int_t p) const { return fData->GetX()[p]; }
  Double_t GetErrorXLeft(Int_t p) const { return fData->GetErrorXlow(p); }
  Double_t GetErrorXRight(Int_t p) const { return fData->GetErrorXhigh(p); }
  Double_t GetY(Int_t point) const { return fData->GetY()[point]; }
  Double_t GetStatError(Int_t point) const { return fData->GetErrorY(point); }
  Double_t GetStatErrorUp(Int_t point) const
  { 
    return fData->GetErrorYhigh(point); 
  }
  Double_t GetStatErrorDown(Int_t point) const
  {
    return fData->GetErrorYlow(point);
  }
  Bool_t IsCommon(Int_t id) const
  {
    HolderCommon* c = FindCommon(id);
    return c != 0;
  }
  Bool_t IsRelative(Int_t id) const
  {
    HolderP2P* h = FindP2P(id);
    if (h) return h->IsRelative();
    HolderCommon* c = FindCommon(id);
    if (c) return c->IsRelative();
    return false;
  }
  Double_t GetSysErrorX(Int_t id, Int_t point) const 
  { 
    HolderP2P* h = FindP2P(id);
    if (!h) return 0;
    return h->GetX(point);
  }
  Double_t GetSysErrorY(Int_t id, Int_t point) const 
  { 
    HolderP2P* h = FindP2P(id);
    if (!h) return 0;
    return h->GetY(point);
  }
  Double_t GetSysErrorXLeft(Int_t id, Int_t point) const 
  { 
    HolderP2P* h = FindP2P(id);
    if (!h) return 0;
    return h->GetXLeft(point);
  }
  Double_t GetSysErrorXRight(Int_t id, Int_t point) const 
  { 
    HolderP2P* h = FindP2P(id);
    if (!h) return 0;
    return h->GetXRight(point);
  }
  Double_t GetSysErrorYUp(Int_t id, Int_t point) const 
  { 
    HolderP2P* h = FindP2P(id);
    if (h) return h->GetYUp(point);
    HolderCommon* c = FindCommon(id);
    if (c) return c->GetYUp(GetY(point));
    return 0;
  }
  Double_t GetSysErrorYDown(Int_t id, Int_t point) const 
  { 
    HolderP2P* h = FindP2P(id);
    if (h) return h->GetYDown(point);
    HolderCommon* c = FindCommon(id);
    if (c) return c->GetYDown(GetY(point));
    return 0;
  }
  const char* GetSysTitle(Int_t id) const
  {
    Holder* h = Find(id);
    if (!h) return "";
    return h->GetTitle();
  }
  Style_t GetSysFillStyle(Int_t id) const
  {
    Holder* h = Find(id);
    if (!h) return 0;
    return h->GetFillStyle();    
  }
  Style_t GetSysLineStyle(Int_t id) const
  {
    Holder* h = Find(id);
    if (!h) return 0;
    return h->GetLineStyle();    
  }
  Color_t GetSysFillColor(Int_t id) const
  {
    Holder* h = Find(id);
    if (!h) return 0;
    return h->GetFillColor();    
  }
  Color_t GetSysLineColor(Int_t id) const
  {
    Holder* h = Find(id);
    if (!h) return 0;
    return h->GetLineColor();    
  }
  Width_t GetSysLineWidth(Int_t id) const
  {
    Holder* h = Find(id);
    if (!h) return 0;
    return h->GetLineWidth();    
  }
  const char* GetSumTitle() const
  {
    return fSumTitle.Data();
  }
  UInt_t GetSumOption() const { return fSumOption; }
  Style_t GetSumFillStyle() const
  {
    return fSumFill.GetFillStyle();
  }
  Style_t GetSumLineStyle() const
  {
    return fSumLine.GetLineStyle();
  }
  Color_t GetSumFillColor() const
  {
    return fSumFill.GetFillColor();
  }
  Color_t GetSumLineColor() const
  {
    return fSumLine.GetLineColor();
  }
  Width_t GetSumLineWidth() const
  {
    return fSumLine.GetLineWidth();
  }
  const char* GetCommonSumTitle() const
  {
    return fCommonSumTitle.Data();
  }
  UInt_t GetCommonSumOption() const { return fCommonSumOption; }
  Style_t GetCommonSumFillStyle() const
  {
    return fCommonSumFill.GetFillStyle();
  }
  Style_t GetCommonSumLineStyle() const
  {
    return fCommonSumLine.GetLineStyle();
  }
  Color_t GetCommonSumFillColor() const
  {
    return fCommonSumFill.GetFillColor();
  }
  Color_t GetCommonSumLineColor() const
  {
    return fCommonSumLine.GetLineColor();
  }
  Width_t GetCommonSumLineWidth() const
  {
    return fCommonSumLine.GetLineWidth();
  }
  Double_t GetCommonErrorYUp(Int_t id) const
  {
    HolderCommon* c = FindCommon(id);
    if (!c) return 0;
    Bool_t rel = c->IsRelative();
    return c->GetYUp(rel ? 100 : 1);
  }
  Double_t GetCommonErrorYDown(Int_t id) const
  {
    HolderCommon* c = FindCommon(id);
    if (!c) return 0;
    Bool_t rel = c->IsRelative();
    return c->GetYDown(rel ? 100 : 1);
  }
  const char* GetXTitle() const { return fXTitle.Data(); }
  const char* GetYTitle() const { return fYTitle.Data(); }
  void GetMinMax(Option_t* option,
     Double_t& ymin, Double_t& ymax) const
  {
    Double_t xmin, xmax;
    Int_t    imin, imax;
    GetMinMax(option, ymin, ymax, xmin, xmax, imin, imax);
  }
  void GetMinMax(Option_t* option,
     Double_t& ymin, Double_t& ymax,
     Double_t& xmin, Double_t& xmax,
     Int_t&    imin, Int_t&    imax) const
  {
    TString  opt(option); opt.ToUpper();
    Bool_t   cmn     = opt.Contains("COMMON");
    Bool_t   stat    = opt.Contains("STAT");
    Bool_t   quad    = opt.Contains("QUAD");
    Bool_t   noerr   = opt.Contains("Y");
    quad             = !opt.Contains("DIRECT");

    ymin = +1e9;
    ymax = -1e9;
    xmin = -1e9;
    xmax = +1e9;
    imin = -1;
    imax = -1;
    for (Int_t i = 0; i < GetN(); i++) {
      Double_t eyl = 0;
      Double_t eyh = 0;
      Double_t y   = noerr ? GetY(i) : GetYandError(i, cmn, stat, quad, false, eyl, eyh);
      Double_t yl  = y - eyl;
      Double_t yh  = y + eyh;
      if (yl < ymin) {
  ymin = yl;
  xmin = GetX(i);
  imin = i;
      }
      if (yh > ymax) {
  ymax = yh;
  xmax = GetX(i);
  imax = i;
      }
      // Printf("%3d: %f -> %f +%f -%f -> %f %f -> %f %f",
      //        i, GetX(i), y, eyl, eyh, yl, yh, ymin, ymax);
    }
    // Printf("min=(%d,%f,%f) max=(%d,%f,%f)", imin, xmin, ymin, imax, xmax, ymax);
  }
  void FindFwhm(Int_t start, Int_t dir, Double_t ymax,
    Bool_t cmn, Bool_t stat, Bool_t quad,
    Int_t& i1, Int_t& i2) const
  {
    Int_t lim = (dir < 0) ? -1 : GetN();
    i1 = i2   = -1;
    for (Int_t i = start+dir; i != lim; i += dir) {
      Double_t eyl, eyh;
      Double_t y = GetYandError(i, cmn, stat, quad, false, eyl, eyh);
      // Printf("%3d: %2d %f -> %f -%f +%f (%f, %f, %f) -> %2d %2d",
      //        i, dir, GetX(i), y, eyl, eyh, y-eyl, y+eyh, ymax/2,
      //        i1, i2);
      if (y-eyl > ymax/2) continue; // Still above
      if (y+eyh < ymax/2) break; // Found lower limit
      if (i1 < 0)  i1 = i2 = i;
      else         i2 = i;
    }
    // Printf("Found %d,%d", i1, i2);
  }
  Double_t FWHM(Double_t& el, Double_t& eh) const
  {
    Double_t xl, xh;
    return FWHM(el, eh, xl, xh);
  }
  Double_t FWHM(Double_t& el, Double_t& eh, Double_t& xl, Double_t& xh) const
  {
    Double_t ymin, ymax, xmin, xmax;
    Int_t    imin, imax;
    GetMinMax("quad stat", ymin, ymax, xmin, xmax, imin, imax);
    if (ymin > ymax/2) {
      Warning("FWHM", "Half of ymax=%f is out of reach [%f,%f]",
        ymax/2, ymin, ymax);
      return -1;
    }

    // Point bounds of found half-max
    Int_t li1, li2, ri1, ri2; 
    // Search left 
    FindFwhm(imax, -1, ymax, false, true, true, li1, li2);
    // Search right 
    FindFwhm(imax, +1, ymax, false, true, true, ri1, ri2);

    Double_t lsign = 1;
    // Check if we've found no left point
    if (li1 < 0 && li2 < 0) {
      // Check if we've found no right point
      if (ri1 < 0 && ri2 < 0) {
  Warning("FWHM", "No left and right point found");
  return -1;
      }
      // Set sign, and copy right point to left point
      lsign = -1;
      li2   = ri1;
      li1   = ri2;
    }
    Double_t rsign = 1;
    // Check if we've found no right point
    if (ri1 < 0 && ri2 < 0) {
      // Check if we've found no left point
      if (li1 < 0 && li2 < 0) {
  Warning("FWHM", "No right and left point found");
  return -1;
      }
      // Set sign, and copy left point to right point
      rsign = -1;
      ri2   = li1;
      ri1   = li2;
    }
    // Printf("Points low=(%d,%d) high=(%d,%d)", li1, li2, ri1, ri2);

    // Now find best estimate of left X
    LinearSigmaCombiner scl;
    scl.Add(lsign*GetX(li1), GetErrorXLeft(li1), GetErrorXRight(li1));
    scl.Add(lsign*GetX(li2), GetErrorXLeft(li2), GetErrorXRight(li2));
    Combiner::Result* lr = scl.Calculate();
    scl.Print();
    lr->Print();
    
    // Now find best estimate of right X
    LinearSigmaCombiner scr;
    scr.Add(GetX(ri1), GetErrorXLeft(ri1), GetErrorXRight(ri1));
    scr.Add(GetX(ri2), GetErrorXLeft(ri2), GetErrorXRight(ri2));
    Combiner::Result* rr = scr.Calculate();
    scr.Print();
    rr->Print();

    // Get error on left/right hand sides
    el            = lr->fEh+rr->fEl;
    eh            = lr->fEl+rr->fEh;
    if (el < 1e-9 && eh < 1e-9) {
      // If the errors are very small, we take the average
      xl   = (GetX(li2)+GetX(li1))/2;
      xh   = (GetX(ri2)+GetX(ri1))/2;
      el   = (GetX(li2)-GetX(li1))/2;
      eh   = (GetX(ri2)-GetX(ri1))/2;
    }
    else {
      // Otherwise we set left/right hand values to best-fit values
      xl   = lr->fX;
      xh   = rr->fX;
    }
    // Evaluate the full-width half-maximum
    Double_t fwhm = (xh - xl);
    return fwhm; 
  }
    
  Double_t MeanX(Double_t& error,
     Bool_t cmn=false,
     Bool_t stat=true,
     Bool_t quad=true) const
  {
    Double_t meanX=TMath::Infinity(), stdDev, n;
    StatisticsX(meanX, stdDev, n, cmn, stat, quad);
    error = stdDev/n;
    return meanX;
  }
  Double_t MeanX(Bool_t cmn=false,
     Bool_t stat=true,
     Bool_t quad=true) const
  {
    Double_t error;
    return MeanX(error, cmn, stat, quad);
  }
  Double_t StandardDeviationX(Double_t& error,
            Bool_t cmn=false,
            Bool_t stat=true,
            Bool_t quad=true) const
  {
    Double_t meanX=TMath::Infinity(), stdDev, eStdDev, n;
    StatisticsX(meanX, stdDev, n, cmn, stat, quad);
    error = stdDev/TMath::Sqrt(2)/n;
    return stdDev;
  }
  Double_t StandardDeviationX(Bool_t cmn=false,
            Bool_t stat=true,
            Bool_t quad=true) const
  {
    Double_t error;
    return StandardDeviationX(error, cmn, stat, quad);
  }
  Double_t StandardDeviationXMean(Double_t mean,
          Double_t& error,
          Bool_t cmn=false,
          Bool_t stat=true,
          Bool_t quad=true) const
  {
    Double_t meanX, stdDev, eStdDev, n;
    StatisticsX(meanX, stdDev, n, cmn, stat, quad);
    Double_t sd2 = stdDev*stdDev - meanX*meanX;
    stdDev = TMath::Sqrt(sd2 - mean*mean);
    error  = stdDev/TMath::Sqrt(2)/n;
    return stdDev;
  }
    
  void StatisticsX(Double_t& meanX,
       Double_t& stdDevX,
       Double_t& n,
       Bool_t cmn=false,
       Bool_t stat=true,
       Bool_t quad=true) const
  {
    Double_t sumY   = 0;
    Double_t sumXY  = 0;
    Double_t sumXXY = 0;
    Double_t sumE2  = 0;
    for (Int_t i = 0; i < GetN(); i++) {
      Double_t e2yl, e2yh;
      Double_t x = GetX(i);
      Double_t y = GetYandError(i,cmn,stat,quad,true,e2yl,e2yh);

      sumY   += y;
      sumXY  += x*y;
      sumXXY += x*x*y;
      sumE2  += TMath::Max(e2yl,e2yh);
    }
    if (meanX == TMath::Infinity())
      meanX       = sumXY/sumY;
    Double_t sd2  = TMath::Abs(sumXXY/sumY-meanX*meanX);
    stdDevX       = TMath::Sqrt(sd2);
    Double_t eff  = sumE2 > 0 ? sumY*sumY/sumE2 : sumY;
    n             = TMath::Sqrt(eff);
  }  
  /* @} */

  //__________________________________________________________________
  void SetSysLineColor(Int_t id, Color_t color) //*MENU*
  {
    Holder* h = Find(id);
    if (!h) return;
    h->SetLineColor(color);
  }
  void SetSysLineStyle(Int_t id, Style_t style) //*MENU*
  {
    Holder* h = Find(id);
    if (!h) return;
    h->SetLineStyle(style);
  }
  void SetSysLineWidth(Int_t id, Width_t width) //*MENU*
  {
    Holder* h = Find(id);
    if (!h) return;
    h->SetLineWidth(width);
  }
  void SetSysFillColor(Int_t id, Color_t color) //*MENU*
  {
    Holder* h = Find(id);
    if (!h) return;
    h->SetFillColor(color);
  }
  void SetSysFillStyle(Int_t id, Style_t style) //*MENU*
  {
    Holder* h = Find(id);
    if (!h) return;
    h->SetFillStyle(style);
  }
  void SetSysTitle(Int_t id, const char* name) //*MENU*
  {
    Holder* h = Find(id);
    if (!h) return;
    h->SetTitle(name);
  }
  void SetSysOption(Int_t id, EDrawOption_t opt) //*MENU*
  {
    Holder* h = Find(id);
    if (!h) return;
    h->SetDOption(opt);
  }
  /* @} */

  //__________________________________________________________________
  void SetSumOption(EDrawOption_t opt) { fSumOption = opt; } //*MENU*
  void SetSumTitle(const char* title) { fSumTitle = title; } //*MENU*
  void SetSumLineColor(Color_t color) { fSumLine.SetLineColor(color); } //*MENU*
  void SetSumLineStyle(Style_t style){ fSumLine.SetLineStyle(style); } //*MENU*
  void SetSumLineWidth(Width_t width){ fSumLine.SetLineWidth(width); }//*MENU*
  void SetSumFillColor(Color_t color){ fSumFill.SetFillColor(color); }//*MENU*
  void SetSumFillStyle(Style_t style){ fSumFill.SetFillStyle(style); }//*MENU*
  void SetCommonSumOption(EDrawOption_t opt) { fCommonSumOption = opt; } //*MENU*
  void SetCommonSumTitle(const char* title) { fCommonSumTitle = title; } //*MENU*
  void SetCommonSumLineColor(Color_t color) { fCommonSumLine.SetLineColor(color); } //*MENU*
  void SetCommonSumLineStyle(Style_t style){ fCommonSumLine.SetLineStyle(style); } //*MENU*
  void SetCommonSumLineWidth(Width_t width){ fCommonSumLine.SetLineWidth(width); }//*MENU*
  void SetCommonSumFillColor(Color_t color){ fCommonSumFill.SetFillColor(color); }//*MENU*
  void SetCommonSumFillStyle(Style_t style){ fCommonSumFill.SetFillStyle(style); }//*MENU*
  /* @} */
  void SetKey(const char* key, const char* value, Bool_t replace=false) //*MENU*
  {
    if (!fMap) { 
      fMap = new TList();
      fMap->SetOwner();
      fMap->SetName("keys");
      // fMap->SetTitle("key/value pairs");
    }
    TString k(key);
    TString v(value);
    TString t =  v.Strip(TString::kBoth);
    v = t;
    if (k.EqualTo("experiment",TString::kIgnoreCase)) {
      TObjArray* l = v.Tokenize("-");
      t = Token(l, 0); TString lab = t.Strip(TString::kBoth);
      t = Token(l, 1); TString acc = t.Strip(TString::kBoth);
      t = Token(l, 2); TString exp = t.Strip(TString::kBoth);
      SetKey("laboratory",  lab.Data(), replace);
      SetKey("accelerator", acc.Data(), replace);
      SetKey("detector",    exp.Data(), replace);
      l->Delete();
    }
    if (k.EqualTo("comment", TString::kIgnoreCase)) { 
      TString l(GetKey("laboratory"));
      TString a(GetKey("accelerator"));
      if (!l.IsNull() && !a.IsNull())
  v.Remove(0, l.Length()+1+a.Length()+1);
      t = v.Strip(TString::kBoth);
      fMap->Add(new TNamed("abstract", t.Data()));
    }
    else {
      if (replace) {
  TObjLink* cur   = fMap->FirstLink();
  TObjLink* last  = fMap->LastLink();
  while (cur) { //  != last) {
    if (!k.EqualTo(cur->GetObject()->GetName())) {
      cur = cur->Next();
      continue;
    }
    TObjLink* tmp = cur->Next();
    fMap->Remove(cur);
    cur = tmp;
  }
      }
      fMap->Add(new TNamed(k, v));
    }
  }
  const char* GetKey(const char* key) const 
  {
    if (!fMap) return 0;
    TObject* o = fMap->FindObject(key);
    if (!o) return 0;
    return o->GetTitle();
  }
  void CopyKeys(const GraphSysErr* g, Option_t* option="fr")
  {
    if (!g) return;
    
    TString opt(option);
    opt.ToLower();

    Bool_t all    = opt.Contains("a");
    Bool_t file   = opt.Contains("f");
    Bool_t header = opt.Contains("h");
    Bool_t qual   = opt.Contains("q");
    Bool_t repl   = opt.Contains("r");
    
    // With the "all" option we just do a plain copy 
    TList*   map = g->fMap;
    TIter    nextKV(map);
    TObject* kv = 0;
    while ((kv = nextKV())) {
      if (!all) {
  TString key = kv->GetName();
  if (file && !(key.EqualTo("reference")  ||
          key.EqualTo("laboratory") ||
          key.EqualTo("accelerator")||
          key.EqualTo("detector")   ||
          key.EqualTo("abstract")   ||
          key.EqualTo("author")   ||
          key.EqualTo("doi")  ||
          key.EqualTo("inspireId")  ||
          key.EqualTo("cdsId")  ||
          key.EqualTo("durhamId")   ||
          key.EqualTo("title")  ||
          key.EqualTo("status"))) continue;
  if (header && !(key.EqualTo("location") ||
      key.EqualTo("reackey")  ||
      key.EqualTo("obskey")   ||
      key.EqualTo("dscomment"))) continue;  
      }
      SetKey(kv->GetName(), kv->GetTitle(), repl);
    }

    if (!all && !qual) return;
    TList*   quals = g->fQualifiers;
    TIter    nextQ(quals);
    TObject* qv = 0;
    while ((qv = nextQ())) {
      const char* test = GetQualifier(qv->GetName());
      if (!test || test[0] == '\0')
  AddQualifier(qv->GetName(), qv->GetTitle());
    }
  }
  void CopyAttr(const GraphSysErr* f)
  {
    fDataOption = f->fDataOption;
    SetMarkerColor(f->GetMarkerColor());
    SetMarkerStyle(f->GetMarkerStyle());
    SetMarkerSize(f->GetMarkerSize());
    SetLineColor(f->GetLineColor());
    SetLineStyle(f->GetLineStyle());
    SetLineWidth(f->GetLineWidth());
    SetFillColor(f->GetFillColor());
    SetFillStyle(f->GetFillStyle());

    fSumOption = f->fSumOption;
    SetSumTitle(f->GetSumTitle());
    SetSumLineStyle(f->GetSumLineStyle());
    SetSumLineColor(f->GetSumLineColor());
    SetSumLineWidth(f->GetSumLineWidth());
    SetSumFillStyle(f->GetSumFillStyle());
    SetSumFillColor(f->GetSumFillColor());

    fCommonSumOption = f->fCommonSumOption;
    SetCommonSumTitle(f->GetCommonSumTitle());
    SetCommonSumLineStyle(f->GetCommonSumLineStyle());
    SetCommonSumLineColor(f->GetCommonSumLineColor());
    SetCommonSumLineWidth(f->GetCommonSumLineWidth());
    SetCommonSumFillStyle(f->GetCommonSumFillStyle());
    SetCommonSumFillColor(f->GetCommonSumFillColor());

    SetXTitle(f->GetXTitle());
    SetYTitle(f->GetYTitle());
  }

    
    
    
  /* @} */
  //__________________________________________________________________
  void AddQualifier(const TString& key, const TString& value,
        Bool_t replace=false)
  {
    if (!fQualifiers) { 
      fQualifiers = new TList();
      fQualifiers->SetOwner();
      fQualifiers->SetName("qualifiers");
    }
    TString val(value);
    if (key.EqualTo(value)) val = "";

    TString k = key.Strip(TString::kBoth, ' ');
    TObject* o = fQualifiers->FindObject(k);
    if (o) {
      Warning("AddQualifier", "Dataset already has qualifier \"%s\"",
        k.Data());
      if (replace) static_cast<TNamed*>(o)->SetTitle(value);
      return;
    }
    
    fQualifiers->Add(new TNamed(k, value));
  }
  void RemoveQualifier(const TString& key)
  {
    if (!fQualifiers) return;
    TObject* o = fQualifiers->FindObject(key);
    if (!o) return;
    fQualifiers->Remove(o);
    delete o;
  }
  const char* GetQualifier(const char* name) const
  {
    if (!fQualifiers) return "";
    TObject* o = fQualifiers->FindObject(name);
    if (!o) return "";
    return o->GetTitle();
  }
  /* @} */
  Bool_t FindYandError(Double_t    x,
           Bool_t      cmn,
           Bool_t      stat,
           Bool_t      quad,
           Bool_t      nosqrt,
           Double_t&   y,
           Double_t&   eyl,
           Double_t&   eyh,
           Double_t&   seyl,
           Double_t&   seyh) const
  {
    Int_t i1, i2;
    Int_t ret = FindPoint(x, i1, i2);
    seyl = seyh = 0;
    if (ret < -1)
      // No valid point found 
      return false;
    if (ret >= 0) {
      // Got exact point 
      // Info("", "Get point %d (%f)", ret, x);
      y = GetYandError(ret, cmn, stat, quad, nosqrt, eyl, eyh);
      if (!stat) {
  // If stat error not include, get them here 
  seyl = GetStatErrorDown(ret);
  seyh = GetStatErrorUp(ret);
      }
      return true;
    }
    // Interpolate between points 
    Double_t eyl1,  eyl2,  eyh1,  eyh2;
    Double_t x1    = fData->GetX()[i1];
    Double_t x2    = fData->GetX()[i2];
    Double_t y1    = GetYandError(i1, cmn, stat, quad, nosqrt, eyl1, eyh1);
    Double_t y2    = GetYandError(i2, cmn, stat, quad, nosqrt, eyl2, eyh2);
    Double_t seyl1 = GetStatErrorDown(i1);
    Double_t seyh1 = GetStatErrorUp(i1);
    Double_t seyl2 = GetStatErrorDown(i2);
    Double_t seyh2 = GetStatErrorUp(i2);
    
    // Linear interpolation
    Double_t dx = (x2-x1);
    Double_t ax = (x-x1)/dx;
    y   = y1   + ax * (y2 - y1);
    eyl = eyl1 + ax * (eyl2 - eyl1);
    eyh = eyh1 + ax * (eyh2 - eyh1);
    if (!stat) {
      // if stat errors not included, calculate here
      seyl = seyl1 + ax * (seyl2 - seyl1);
      seyh = seyh1 + ax * (seyh2 - seyh1);
    }
    return true;
  }
  Bool_t FindYandError(Double_t    x,
           Bool_t      cmn,
           Bool_t      stat,
           Bool_t      quad,
           Bool_t      nosqrt,
           Double_t&   y,
           Double_t&   eyl,
           Double_t&   eyh) const
  {
    Double_t seyl;
    Double_t seyh;
    return FindYandError(x, cmn, stat, quad, nosqrt, y, eyl, eyh, seyl, seyh);
  }

  //__________________________________________________________________
  struct Combiner : public TObject
  {
    typedef Double_t (*Wrapper_t)(Double_t*,Double_t*);
    struct Observation : public TObject 
    {
      Double_t fX;
      Double_t fEl;
      Double_t fEh;
      Observation(Double_t x=0, Double_t el=0, Double_t eh=0)
  : fX(x), fEl(el), fEh(eh)
      {    
      }
      virtual ~Observation() {}
      Double_t S() const
      {
  if (fEl * fEh == 0) return 0;
  return 2 * fEl * fEh / (fEl + fEh);
      }
      Double_t Sprime() const
      {
  if (TMath::Abs(fEh + fEl) < 1e-9) return 1; 
  return (fEh - fEl) / (fEh + fEl);
      }
      Double_t Svar(Double_t guess=0) const
      {
  return  S() + Sprime() * (guess - fX);
      }
      Double_t V() const  { return fEl * fEh; }
      Double_t Vprime() const { return fEh - fEl;  }
      virtual Double_t Low() const { return fX - 3 * fEl; }
      virtual Double_t High() const { return fX + 3 * fEh; }
      virtual void Print(Option_t* option="") const
      {
  Printf("%10.8f  -%10.8f +%10.8f", fX, fEl, fEh);
      }
      ClassDef(Observation,1); // An experimental observation 
    };
  
    struct Result : public Observation
    {
      Double_t fChi2;
      Double_t fLower;
      Double_t fUpper;
      Result(Double_t x=0, Double_t el=0, Double_t eh=0,
       Double_t chi2=0, Double_t low=0, Double_t high=0)
  : Observation(x, el, eh),
    fChi2(chi2),
    fLower(low),
    fUpper(high)
      {}
      Double_t Low() const { return fLower; }
      Double_t High() const { return fUpper; }
      virtual void Print(Option_t* option="") const
      {
  Printf("%10.8f  -%10.8f +%10.8f  %5.2f", fX, fEl, fEh, fChi2);
      }
      ClassDef(Result,1); // Final result 
    };
    TClonesArray fData;   // Container of observations
    Result*      fResult; // Result of combining

    Combiner()
      : fData("GraphSysErr::Combiner::Observation"), fResult(0)
    {
      fData.SetOwner();
    }
    virtual ~Combiner() { fData.Clear(); if (fResult) delete fResult; }
    void Clear(Option_t* option="")
    {
      fData.Clear();
      if (fResult) delete fResult;
      fResult = 0;
    }
    void Add(const Observation& r)
    {
      Add(r.fX,r.fEl,r.fEh);
    }
    void Add(Double_t x, Double_t el, Double_t eh)
    {
      // if (TMath::Abs(x)  < 1e-10 &&
      //     TMath::Abs(el) < 1e-10 &&
      //     TMath::Abs(eh) < 1e-10) return;
      new (fData[fData.GetEntries()]) Observation(x,el,eh);
    }
    void Print(Option_t* option="") const
    {
      TIter next(&fData);
      Observation* r = 0;
      while ((r = static_cast<Observation*>(next())))
  r->Print(option);
    }
    /* @} */
  
    virtual Double_t W(const Observation& r) const = 0;
    virtual Double_t StepW(Double_t guess, const Observation& r) const = 0;
    virtual Double_t StepOffset(Double_t guess, const Observation& r) const = 0;
    virtual Double_t VarTerm(Double_t guess, const Observation& r) const = 0;
    Double_t ChiTerm(Double_t guess, const Observation& r) const
    {
      Double_t var = VarTerm(guess, r);
      if (var <= 0) return -1000;

      return TMath::Power(guess - r.fX, 2) / var;
    }
    Double_t F(Double_t          guess,
         Double_t          chi2) const
    {
      if (fData.GetEntries() == 1) return 1;
      Double_t s = -chi2;

      TIter next(&fData);
      Observation* r = 0;
      while ((r = static_cast<Observation*>(next())))
  s += ChiTerm(guess, *r);
      return s;
    }
    Double_t E(UShort_t        nIter,
         Int_t           sign,
         Double_t        best,
         Double_t        chi2,
         Double_t        s)
    {
      if (fData.GetEntries() == 1) {
  Observation* o = static_cast<Observation*>(fData[0]);
  return (sign < 0 ? o->fEl : o->fEh);
      }

      // Step size 
      Double_t delta = 0.1 * sign * s;

      // Iterations 
      for (UShort_t i = 0; i < nIter; i++) {
  // Calculate chi^2 with current guess 
  Double_t got = F(best + sign * s, chi2);

  if (TMath::Abs(got-1) < 1e-7)
    // We're close to 1 so get out
    break;

  // The next guess' chi^2 value e
  Double_t guess = F(best + sign * s + delta, chi2);

  // Where to search next 
  if ((got - 1) * (guess - 1) > 0) {
    if ((got - 1) / (guess - 1) < 1)
      delta = -delta;
    else
      s += sign * delta;
    continue;
  }

  // Update error value and decrease step size 
  s     += sign * delta * (1 - got) / (guess - got);
  delta /= 2;
      }
      return s;
    }
    Double_t X(UShort_t      nIter,
         Double_t      lowest,
         Double_t      highest)
    {
      // Starting values
      if (fData.GetEntries() == 1)
  return static_cast<Observation*>(fData[0])->fX;
      Double_t x    = (highest+lowest)/2;
      Double_t oldX = -1e33;
      // Do the iterations 
      for (UShort_t i = 0; i < nIter; i++) {
  Double_t sum    = 0;
  Double_t sumw   = 0;
  Double_t offset = 0;

  // Loop over observations
  TIter next(&fData);
  Observation* r = 0;
  while ((r = static_cast<Observation*>(next()))) {
    Double_t w =  StepW(x, *r);
    offset     += StepOffset(x,  *r);
    sum  += r->fX * w;
    sumw += w;
  }
  x = (TMath::Abs(sumw) > 1e-9) ? (sum - offset) / sumw : 0;

  if (TMath::Abs(x - oldX) < (highest-lowest) * 1e-9) break;
  oldX = x;
      }
      return x;
    }
    Result* Calculate(UShort_t nIter=50)
    {
      Double_t lowest  = +1e33;
      Double_t highest = -1e33;
      Double_t sumLow  = 0;
      Double_t sumHigh = 0;

      // Find boundaries and sum weights
      TIter next(&fData);
      Observation* r = 0;
      while ((r = static_cast<Observation*>(next()))) {
  lowest  = TMath::Min(r->Low(),  lowest);
  highest = TMath::Max(r->High(), highest);
  if (r->fEl > 1e-10)  sumLow  += 1./TMath::Power(r->fEl, 2);
  if (r->fEh > 1e-10)  sumHigh += 1./TMath::Power(r->fEh, 2);
      }
      // Summed weights 
      Double_t sLow  = sumLow  > 1e-10 ? 1. / TMath::Sqrt(sumLow) : 0;
      Double_t sHigh = sumHigh > 1e-10 ? 1. / TMath::Sqrt(sumHigh) : 0;
    
      // Now do the calculations
      Double_t bestX    = X(nIter, lowest, highest);
      Double_t bestChi2 = F(bestX, 0);
      Double_t bestLow  = TMath::Abs(E(nIter,-1,bestX,bestChi2,sLow));
      Double_t bestHigh = TMath::Abs(E(nIter,+1,bestX,bestChi2,sHigh));

      fResult = new Result(bestX, bestLow, bestHigh, bestChi2, lowest, highest);
      return fResult;
    }
    virtual Wrapper_t Wrapper() const = 0;
    
    TF1* MakeF(const Observation& r, Int_t j) const
    {
      TF1* f = new TF1(Form("f%02d", j), Wrapper(), r.Low(), r.High(), 3);
      f->SetParNames("x", "#sigma^{-}", "#sigma^{+}");
      f->SetParameters(r.fX, r.fEl, r.fEh);
      f->SetLineStyle((j%3)+2);
      f->SetLineColor(kBlack);
      f->SetLineWidth(2);
      return f;
    }
    TLine* MakeL(TF1* f) const
    {
      Double_t m = f->GetParameter(0);
      TLine* l = new TLine(m-f->GetParameter(1), 1,
         m+f->GetParameter(2), 1);
      l->SetLineColor(f->GetLineColor());
      l->SetLineStyle(f->GetLineStyle());
      l->SetLineWidth(f->GetLineWidth());
      return l;
    }
    void Draw(Option_t* option="")
    {
      if (!fResult) return;
      
      TList fs; fs.SetOwner(false);
      Int_t j = 0;
      TIter next(&fData);
      Observation* r = 0;
      while ((r = static_cast<Observation*>(next()))) {
  TF1* f = MakeF(*r, j);
  f->SetRange(fResult->Low(), fResult->High());
  fs.Add(f, j == 0 ? "" : "same");
  fs.Add(MakeL(f));
  j++;
      }
      TF1* fr = MakeF(*fResult, j);
      fr->SetLineColor(kRed+2);
      fr->SetLineStyle(1);
      fs.Add(fr);
      fs.Add(MakeL(fr));
      TIter nextD(&fs);
      TObject* o = 0;
      j = 0;
      while((o = nextD())) { o->Draw(j == 0 ? "" : "same"); j++; }
      // fs.Draw();
      static_cast<TF1*>(fs.First())->GetHistogram()
  ->GetYaxis()->SetRangeUser(-.1, 5.1);
    }
    ClassDef(Combiner,1); // Combine systematics 
  };


  struct LinearSigmaCombiner : public Combiner
  {
    Double_t W(const Observation& r) const
    {
      // Double_t s = r.S();
      // if (s < 1e-10) return 0;
      // return .5 * TMath::Power(r.Svar(), 3) / s;
      Double_t w = StepW(0,r);
      return 1/w;
    }
    Double_t StepW(Double_t guess, const Observation& r) const
    {
      Double_t s   = r.S();
      Double_t t   = r.Svar(guess);
      Double_t ret = s / TMath::Power(t,3);
      return ret;
    }
    Double_t StepOffset(Double_t, const Observation&) const
    {
      return 0;
    }
    Double_t VarTerm(Double_t guess, const Observation& r) const
    {
      return TMath::Power(r.Svar(guess),2);
    }
    static Double_t L(Double_t guess, Double_t x, Double_t el, Double_t eh)
    {
      Observation tmp(x,el,eh);
      Double_t d  = (guess-x);
      Double_t t  = tmp.Svar(guess); // (s+sp*d);
      if (TMath::Abs(d) < 1e-10) return 0;
      // if (TMath::Abs(t) < 1e-10) return DBL_MAX;
      Double_t ret = TMath::Power(d/t, 2);
      return ret;
    }
    static Double_t WrapL(Double_t* xp, Double_t* pp)
    {
      return L(xp[0], pp[0], pp[1], pp[2]);
    }
    Wrapper_t Wrapper() const { return WrapL; }
    ClassDef(LinearSigmaCombiner,1); // Sigma combiner 
  };
  struct LinearVarianceCombiner : public Combiner
  {
    Double_t W(const Observation& r) const
    {
      Double_t v  = r.V();
      Double_t vp = r.Vprime();
      return TMath::Power(v + r.fX * vp, 2) / (2 * v + r.fX * vp);
    }
    Double_t StepW(Double_t guess, const Observation& r) const
    {
      Double_t v = r.V();
      return v / TMath::Power(v+r.Vprime()*(guess - r.fX), 2);
    }
    Double_t StepOffset(Double_t guess, const Observation& r) const
    {
      Double_t vp = r.Vprime();
      return 0.5 * vp * TMath::Power((guess-r.fX)/(r.V()+vp*(guess-r.fX)),2);
    }
    Double_t VarTerm(Double_t guess, const Observation& r) const
    {
      return r.V() + r.Vprime() * (guess - r.fX);
    }
    static Double_t L(Double_t guess, Double_t x, Double_t el, Double_t eh)
    {
      Double_t d  = (guess-x);
      Double_t v  = eh * el;
      Double_t vp = eh - el;
      return TMath::Power(d,2) / (v+vp*d);
    }
    static Double_t WrapL(Double_t* xp, Double_t* pp)
    {
      return L(xp[0], pp[0], pp[1], pp[2]);
    }
    Wrapper_t Wrapper() const { return WrapL;  }
    ClassDef(LinearVarianceCombiner,1); // Variance combiner 
  };
  
  //__________________________________________________________________
  struct Holder : public TNamed, TAttLine, TAttFill 
  {
    friend struct GraphSysErr;
    Holder() 
      : TNamed(), 
  TAttLine(),
  TAttFill(),
  fRelative(true),
  fOption(kRect)
    {}
    virtual ~Holder() {}
    void CopyAttr(Holder* h)
    {
      SetLineColor(h->GetLineColor());
      SetLineStyle(h->GetLineStyle());
      SetLineWidth(h->GetLineWidth());
      SetFillColor(h->GetFillColor());
      SetFillStyle(h->GetFillStyle());
      fOption = h->fOption;
    }
  protected:
    Holder(const char* name, const char* title, Bool_t rel,
     UInt_t option, UInt_t id)
      : TNamed(name, title), 
  TAttLine(0,0,0),
  TAttFill(0,0),
  fRelative(rel),
  fOption(option)
    {
      SetUniqueID(id);
    }
    Holder(const Holder& other) 
      : TNamed(other), 
  TAttLine(other),
  TAttFill(other),
  fRelative(other.fRelative),
  fOption(other.fOption)
    {
      SetUniqueID(other.GetUniqueID());
    }
    Holder& operator=(const Holder& other) 
    {
      if (&other == this) return *this;

      other.TNamed::Copy(*this);
      other.TAttLine::Copy(*this);
      other.TAttFill::Copy(*this);

      fRelative = other.fRelative;
      fOption   = other.fOption;
    
      SetUniqueID(other.GetUniqueID());

      return *this;
    }
    
    virtual Graph* StackError(Graph* g, Bool_t ignoreErr, Bool_t quad) const = 0;
    virtual void SumError(Graph* g, Int_t i, Bool_t ignoreErr, 
        Bool_t quad, UInt_t opt) const = 0;
    virtual UInt_t GetDOption() const { return fOption; }
    virtual void SetDOption(EDrawOption_t opt) { fOption = opt; }
    virtual Bool_t IsRelative() const { return fRelative; }

    virtual void Print(Option_t* option="") const
    {
      gROOT->IndentLevel();
      Printf("%s/%s %s(ID # %d%s)", GetName(), GetTitle(),
       IsRelative() ? "PCT " : "",  GetUniqueID(),
       TestBit(kUsedBit) ? " used" : "");
      TString opt(option);
      if (!opt.Contains("ATTR", TString::kIgnoreCase)) return;
      gROOT->IndentLevel();
      Printf(" [option: %3d line (c/s/w):%3d/%1d/%2d fill (c/s):%3d/%4d]",
       fOption, GetLineColor(), GetLineStyle(), GetLineWidth(),
       GetFillColor(), GetFillStyle());
      
    }
    virtual void ls(Option_t* option) const 
    {
      Print(option);
    }
  protected:
    UShort_t XMode(Int_t opt=-1) const
    {
      UShort_t xMode  = 0;
      if (opt < 0) opt = fOption;
      switch (opt) {
      case kNormal: case kNoTick: case kFill: case kCurve: xMode = 0; break;
      case kArrow:  case kHat:    case kBar:               xMode = 1; break;
      case kNone:                                          xMode = 0; break;
      case kRect:   case kBox:                             xMode = 3; break;
      }
      return xMode;
    }
    void DoAdd(UShort_t  xMode,
         Double_t  curExl,
         Double_t  curExh,
         Double_t  curEyl,
         Double_t  curEyh,
         Bool_t    ignoreErr,
         Bool_t    quad, 
         Bool_t    sqOld,
         Double_t& exl, 
         Double_t& exh, 
         Double_t& eyl, 
         Double_t& eyh) const
    {
      // Double_t oexl = exl;
      // Double_t oexh = exh;
      if (quad) { 
  eyl    *= eyl;
  eyh    *= eyh;
  if (sqOld) {
    curEyl *= curEyl;
    curEyh *= curEyh;
  }
      }
  
      if (!ignoreErr) {
  if (kVerbose & kDraw)
    Info("", "old: +%f/-%f this: +%f/-%f -> +%f/-%f",
         curEyl, curEyh, eyl, eyh, eyl + curEyl, eyh + curEyh);
  exl  = TMath::Max(exl, curExl);
  exh  = TMath::Max(exh, curExh);
  eyl  += curEyl;
  eyh  += curEyh;
      }

      switch (xMode) {
      case 0: break;                // kNormal, kNoTick, kFill, kCurve, kNone
      case 1: exl = exh = 0; break; // kArrow, kHat, kBar
      case 2:                       // kRect, kBox
  if (exl <= 0) exl = gStyle->GetErrorX()/2;
  if (exh <= 0) exh = gStyle->GetErrorX()/2;
  break;
      }
    }
    
    void SetAttributes(Graph* g) const
    {
      if (!g) return;
      this->TAttLine::Copy(*g);
      this->TAttFill::Copy(*g);
      g->SetMarkerColor(0);
      g->SetMarkerStyle(0);
      g->SetMarkerSize(0);
    }
    Bool_t  fRelative;
    UInt_t fOption;
    ClassDef(Holder,3);
  };
  //__________________________________________________________________
  struct HolderP2P : public Holder
  {
    friend struct GraphSysErr;
    HolderP2P()
      : Holder(),
  fGraph(0)
    {}
  protected:
    HolderP2P(const char* name, const char* title, Bool_t rel,
        UInt_t  opt, UInt_t id)
      : Holder(name, title, rel, opt, id),
  fGraph(0)
    {
      fGraph = new Graph();
      fGraph->SetName(name);
      fGraph->SetTitle(title);
      SetAttributes(fGraph);
    }
    HolderP2P(const HolderP2P& other)
      : Holder(other),
  fGraph(0)
    {
      if (other.fGraph) 
  fGraph = static_cast<Graph*>(other.fGraph->Clone());

      if (fGraph) SetAttributes(fGraph);
    }
    HolderP2P& operator=(const HolderP2P& other)
    {
      if (&other == this) return *this;
      
      this->Holder::operator=(other);
      if (fGraph) delete fGraph;
      fGraph = 0;
      if (other.fGraph) 
  fGraph = static_cast<Graph*>(other.fGraph->Clone());
      if (fGraph) SetAttributes(fGraph);
      
      return *this;
    }

    void Set(Int_t point, Graph* g, Double_t ex, Double_t ey) 
    {
      Set(point, g, ex, ex, ey, ey);
    }
    void Set(Int_t point, Graph* g, Double_t ex1, Double_t ex2, 
       Double_t ey1, Double_t ey2) 
    {
      if (!fGraph) return;
      if (!g) return;
      if (point >= g->GetN()) return;
      Double_t x = g->GetX()[point];
      Double_t y = g->GetY()[point];
      fGraph->SetPoint(point, x, y);
      if (fRelative) { ey1 *= y; ey2 *= y; }
      fGraph->SetPointError(point, ex1, ex2, ey1, ey2);
    }
    /* @} */

    Double_t GetX(Int_t point) const 
    {
      if (!fGraph) return 0;
      return fGraph->GetErrorX(point);
    }
    Double_t GetXLeft(Int_t point) const 
    {
      if (!fGraph) return 0;
      return fGraph->GetErrorXlow(point);
    }
    Double_t GetXRight(Int_t point) const 
    {
      if (!fGraph) return 0;
      return fGraph->GetErrorXhigh(point);
    }
    Double_t GetY(Int_t point) const 
    {
      if (!fGraph) return 0;
      return fGraph->GetErrorY(point);
    }
    Double_t GetYDown(Int_t i1, Int_t i2, Double_t ax) const 
    {
      if (!fGraph) return 0;
      if (i1 == i2) return GetYDown(i1);
      Double_t e1 = GetYDown(i1);
      Double_t e2 = GetYDown(i2);
      return e1 + ax * (e2-e1); // Linear interpolation
    }
    Double_t GetYDown(Int_t point) const 
    {
      if (!fGraph) return 0;
      return fGraph->GetErrorYlow(point);
    }
    Double_t GetYUp(Int_t i1, Int_t i2, Double_t ax) const 
    {
      if (!fGraph) return 0;
      if (i1 == i2) return GetYUp(i1);
      Double_t e1 = GetYUp(i1);
      Double_t e2 = GetYUp(i2);
      return e1 + ax * (e2-e1); // Linear interpolation
    }
    Double_t GetYUp(Int_t point) const 
    {
      if (!fGraph) return 0;
      return fGraph->GetErrorYhigh(point);
    }
    /* @} */

    void AddError(Int_t     i,
      UShort_t  xMode,
      Bool_t    ignoreErr,
      Bool_t    quad, 
      Bool_t    sqOld,
      Double_t& exl, 
      Double_t& exh, 
      Double_t& eyl, 
      Double_t& eyh) const
    {
      Double_t exlO = fGraph->GetErrorXlow(i);
      Double_t exhO = fGraph->GetErrorXhigh(i);
      Double_t eylO = fGraph->GetErrorYlow(i);
      Double_t eyhO = fGraph->GetErrorYhigh(i);
      Double_t oyl  = eylO;
      Double_t oyh  = eyhO;
      if (exl <= 0) exlO = exl;
      if (exh <= 0) exhO = exh;

      DoAdd(xMode,
      exl, exh, eyl, eyh,
      ignoreErr, quad, sqOld,
      exlO, exhO, eylO, eyhO);
      if (kVerbose & kDraw)
  Info(GetTitle(), "quad: %s this: +%f/-%f, old: +%f/-%f -> +%f/-%f",
       (quad ? "true" : "false"), oyl, oyh, eyl, eyh, eylO, eyhO);

      exl = exlO;
      exh = exhO;
      eyl = eylO;
      eyh = eyhO;
    }
    void StackPointError(Int_t     i,
       UShort_t  xMode,
       Bool_t    ignoreErr, 
       Bool_t    quad,
       Double_t& exl,
       Double_t& exh, 
       Double_t& eyl,
       Double_t& eyh) const
    {
      Double_t oldEyl = eyl;
      Double_t oldEyh = eyh;
      AddError(i, xMode, ignoreErr, quad, true, exl, exh, eyl, eyh);
      if (kVerbose & kDraw)
  Info(GetTitle(), "old= +%f/-%f -> +%f/-%f", oldEyl, oldEyh, eyl, eyh);
    }
    Graph* StackError(Graph* g, Bool_t ignoreErr, Bool_t quad) const
    {
      Graph* ga = new Graph(g->GetN());
      for (Int_t i = 0; i < g->GetN(); i++) { 
  Double_t x    = g->GetX()[i];   
  Double_t y    = g->GetY()[i];
  Double_t exl  = g->GetErrorXlow(i);
  Double_t exh  = g->GetErrorXhigh(i);
  Double_t eyl  = g->GetErrorYlow(i);
  Double_t eyh  = g->GetErrorYhigh(i);
  StackPointError(i, XMode(), ignoreErr, quad, exl, exh, eyl, eyh);
  // AddError(i, g, quad, ignoreErr, true, exl, exh, eyl, eyh);
  ga->SetPoint(i, x, y);
  ga->SetPointError(i, exl,exh,eyl,eyh);
      }
      SetAttributes(ga);
      ga->SetTitle(GetTitle());
      ga->SetName(Form("stack_%08x", GetUniqueID()));
      return ga;
    }
    void SumPointError(Int_t i, UShort_t xMode, Bool_t ignoreErr, Bool_t quad, 
           Double_t& exl,
           Double_t& exh,
           Double_t& eyl,
           Double_t& eyh) const
    {
      AddError(i, xMode, ignoreErr, quad, false, exl, exh, eyl, eyh);
    }
    void SumError(Graph* g, Int_t i, Bool_t ignoreErr,
      Bool_t quad, UInt_t opt) const
    {
      Double_t exl = (g ? g->GetErrorXlow(i)  : 0);
      Double_t exh = (g ? g->GetErrorXhigh(i) : 0);
      Double_t eyl = (g ? g->GetErrorYlow(i)  : 0);
      Double_t eyh = (g ? g->GetErrorYhigh(i) : 0);
      SumPointError(i, XMode(opt), ignoreErr, quad, exl, exh, eyl, eyh);
      g->SetPointError(i, exl,exh,eyl,eyh);
    }
    /* @} */
    void SavePrimitive(std::ostream& out, Option_t* option="")
    {
      if (option[0] == 'd' || option[0] == 'D')
  out << " // Point-2-Point systematic " << GetTitle() << "\n"
      << "  {\n"
      << "    Int_t id = g->DeclarePoint2Point(\"" << GetTitle() << "\","
      << IsRelative() <<  ',' << fOption << ");\n"
      << "    g->SetSysLineColor(id," << GetLineColor() << ");\n"
      << "    g->SetSysLineStyle(id," << GetLineStyle() << ");\n"
      << "    g->SetSysLineWidth(id," << GetLineWidth() << ");\n"
      << "    g->SetSysFillColor(id," << GetFillColor() << ");\n"
      << "    g->SetSysFillStyle(id," << GetFillStyle() << ");\n"
      << "  }\n";
    }
    virtual void Print(Option_t* option="") const
    {
      gROOT->IndentLevel();
      Printf("%s/%s (ID # %d, %d points)", GetName(), GetTitle(), GetUniqueID(),
       (fGraph ? fGraph->GetN() : -1));
      TString opt(option);
      if (!opt.Contains("ATTR", TString::kIgnoreCase)) return;
      gROOT->IndentLevel();
      Printf(" [option: %3d line (c/s/w):%3d/%1d/%2d fill (c/s):%3d/%4d]",
       fOption, GetLineColor(), GetLineStyle(), GetLineWidth(),
       GetFillColor(), GetFillStyle());
      
    }
    // Graph* fGraph;
    TGraphAsymmErrors* fGraph;

    ClassDef(HolderP2P,3);
  };
  //__________________________________________________________________
  struct HolderCommon : public Holder 
  {
    friend struct GraphSysErr;
    HolderCommon()
      : Holder(), fEyl(0), fEyh(0)
    {}
  protected:
    HolderCommon(const char* name, const char* title, Bool_t rel,
     UInt_t opt, UInt_t id)
      : Holder(name, title, rel, opt, id), fEyl(0), fEyh(0)
    {
    }
    HolderCommon(const HolderCommon& other)
      : Holder(other),
  fEyl(other.fEyl), 
  fEyh(other.fEyh)
    {
    }
    HolderCommon& operator=(const HolderCommon& other)
    {
      if (&other == this) return *this;
      this->Holder::operator=(other);
      fEyl = other.fEyl;
      fEyh = other.fEyh;
      return *this;
    }

    void Set(Double_t ey) { Set(ey, ey); }
    void Set(Double_t eyl, Double_t eyh)
    {
      fEyl = eyl;
      fEyh = eyh;
    }
    Double_t GetYDown(Double_t y=0) const 
    {
      return (fRelative ? y : 1) * fEyl;
    }
    Double_t GetYUp(Double_t y=0) const 
    {
      return (fRelative ? y : 1) * fEyh;
    }
    /* @} */

    void AddError(Double_t  y,
      UShort_t  xMode,
      Bool_t    ignoreErr,
      Bool_t    quad, 
      Bool_t    sqOld,
      Double_t& exl, 
      Double_t& exh, 
      Double_t& eyl, 
      Double_t& eyh) const
    {
      //exl         = (i   == 0        ? 0 :(g->GetX()[i]  -g->GetX()[i-1])/2);
      //exh         = (i+1 >= g->GetN()? 0 :(g->GetX()[i+1]-g->GetX()[i])/2);
      Double_t exlO = exl;
      Double_t exhO = exh;
      Double_t eylO = GetYDown(y);
      Double_t eyhO = GetYUp(y);
      if (exlO <= 0) exlO = gStyle->GetErrorX()/2;
      if (exhO <= 0) exhO = gStyle->GetErrorX()/2;

      DoAdd(xMode, exl, exh, eyl, eyh,
      ignoreErr, quad, sqOld,
      exlO, exhO, eylO, eyhO);
      exl = exlO;
      exh = exhO;
      eyl = eylO;
      eyh = eyhO;
    }
    Graph* BarError(Graph* g, Bool_t quad, Double_t x, Double_t y) const
    {
      Double_t exl = (g ? g->GetErrorXlow(0)  : 0);
      Double_t exh = (g ? g->GetErrorXhigh(0) : 0);
      Double_t eyl = (g ? g->GetErrorYlow(0)  : 0);
      Double_t eyh = (g ? g->GetErrorYhigh(0) : 0);
      AddError(y, XMode(), false, quad, quad, exl, exh, eyl, eyh);

      Graph* r = new Graph(1);
      r->SetPoint(0, x, y);
      r->SetPointError(0, exl, exh, eyl, eyh);
      
      SetAttributes(r);
      r->SetTitle(GetTitle());
      r->SetName(Form("bar_%08x", GetUniqueID()));

      return r;
    }
    void StackPointError(Double_t  y,
       UShort_t  xMode,
       Bool_t    ignoreErr, 
       Bool_t    quad,
       Double_t& exl,
       Double_t& exh, 
       Double_t& eyl,
       Double_t& eyh) const
    {
      AddError(y, xMode, ignoreErr, quad, true, exl, exh, eyl, eyh);
    }
    
    Graph* StackError(Graph* g, Bool_t ignoreErr, Bool_t quad) const
    {
      Graph* ga = new Graph(g->GetN());
      for (Int_t i = 0; i < g->GetN(); i++) { 
  Double_t x    = g->GetX()[i];
  Double_t y    = g->GetY()[i]; 
  Double_t exl  = g->GetErrorXlow(i);
  Double_t exh  = g->GetErrorXhigh(i);
  Double_t eyl  = g->GetErrorYlow(i);
  Double_t eyh  = g->GetErrorYhigh(i);
  StackPointError(y, XMode(), ignoreErr, quad, exl, exh, eyl, eyh);
  ga->SetPoint(i, x, y);
  ga->SetPointError(i, exl,exh,eyl,eyh);
      }
      SetAttributes(ga);
      ga->SetTitle(GetTitle());
      ga->SetName(Form("stack_%08x", GetUniqueID()));
      return ga;
    }
    void SumPointError(Double_t y, UShort_t xMode, Bool_t ignoreErr,Bool_t quad,
           Double_t& exl,
           Double_t& exh,
           Double_t& eyl,
           Double_t& eyh) const
    {
      AddError(y, xMode, ignoreErr, quad, false, exl, exh, eyl, eyh);
    }
    void SumError(Graph* g, Int_t i, Bool_t ignoreErr, Bool_t quad,
      UInt_t opt) const
    {
      Double_t y   = (g ? g->GetY()[i]        : 0);
      Double_t exl = (g ? g->GetErrorXlow(i)  : 0);
      Double_t exh = (g ? g->GetErrorXhigh(i) : 0);
      Double_t eyl = (g ? g->GetErrorYlow(i)  : 0);
      Double_t eyh = (g ? g->GetErrorYhigh(i) : 0);
      SumPointError(y, XMode(opt), ignoreErr, quad, exl, exl, eyl, eyh);
      g->SetPointError(i, exl,exh,eyl,eyh);
    }
    /* @} */
    void SavePrimitive(std::ostream& out, Option_t* option="")
    {
      if (option[0] == 'd' || option[0] == 'D')
  out << " // Common systematic " << GetTitle() << "\n"
      << "  {\n"
      << "    Int_t id = g->DefineCommon(\"" << GetTitle() << "\","
      << fRelative << ',' << fEyl << ',' << fEyh << ','
      << fOption << ");\n"
      << "    g->SetSysLineColor(id," << GetLineColor() << ")\n;"
      << "    g->SetSysLineStyle(id," << GetLineStyle() << ")\n;"
      << "    g->SetSysLineWidth(id," << GetLineWidth() << ")\n;"
      << "    g->SetSysFillColor(id," << GetFillColor() << ")\n;"
      << "    g->SetSysFillStyle(id," << GetFillStyle() << ")\n;"
      << "  }\n";
    }
    virtual void Print(Option_t* option="") const
    {
      Bool_t      rel = IsRelative();
      Double_t    fac = (rel ? 100 : 1);
      const char* pst = (rel ? "%" : "");
      gROOT->IndentLevel();
      Printf("%s/%s (ID # %d) -%f%s +%f%s",
       GetName(), GetTitle(), GetUniqueID(),
       fac*fEyl, pst, fac*fEyh, pst);
      TString opt(option);
      if (!opt.Contains("ATTR", TString::kIgnoreCase)) return;
      gROOT->IndentLevel();
      Printf(" [option: %3d line (c/s/w):%3d/%1d/%2d fill (c/s):%3d/%4d]",
       fOption, GetLineColor(), GetLineStyle(), GetLineWidth(),
       GetFillColor(), GetFillStyle());
      
    }
    Double_t fEyl;
    Double_t fEyh;

    ClassDef(HolderCommon,3);
  };
  Double_t GetYandError(Int_t     i,
      Bool_t    cmn,
      Bool_t    stat,
      Bool_t    quad,
      Bool_t    nosqrt,
      Double_t& eyl,
      Double_t& eyh) const
  {
    Double_t wyl = 0;
    Double_t wyh = 0;
    return GetYandError(i, cmn, stat, quad, nosqrt, eyl, eyh, wyl, wyh);
  }
  Double_t GetYandError(Int_t     i,
      Bool_t    cmn,
      Bool_t    stat,
      Bool_t    quad,
      Bool_t    nosqrt,
      Double_t& eyl,
      Double_t& eyh,
      Double_t& wyl,
      Double_t& wyh) const
  {
    // --- Find location for common errors ---------------------------
    Int_t n = fData->GetN();
    if (i >= n) {
      eyl = -1;
      eyh = -1;
      return 0;
    }
    Double_t y   = fData->GetY()[i];
    wyl = GetStatErrorDown(i);
    wyh = GetStatErrorUp(i);
    eyl = (stat ? wyl : 0);
    eyh = (stat ? wyh : 0);
    if (quad) {
      wyl *= wyl;
      wyh *= wyh;
      eyl *= eyl;
      eyh *= eyh;
    }
    Double_t exl = GetErrorXLeft(i);
    Double_t exh = GetErrorXRight(i);

    // Otherwise, we are adding all selected errors togethter 
    // Graph* g = static_cast<Graph*>(fData->Clone("error"));
    // g->SetTitle(fSumTitle);
    Int_t xMode = -1;
    if (cmn) { 
      TIter nextC(&fCommon);
      HolderCommon* hc = 0;
      while ((hc = static_cast<HolderCommon*>(nextC()))) {
  if (hc->TestBit(kUsedBit)) continue;
  if (xMode < 0) xMode = hc->XMode(fSumOption);
  hc->SumPointError(y, xMode, false, quad, exl, exh, wyl, wyh);
  if (hc->TestBit(kOnlyWeightBit)) continue;
  hc->SumPointError(y, xMode, false, quad, exl, exh, eyl, eyh);
      }
    }
    TIter nextP(&fPoint2Point);
    HolderP2P* hp = 0;
    while ((hp = static_cast<HolderP2P*>(nextP()))) {
      if (hp->TestBit(kUsedBit)) continue;
      if (xMode < 0) xMode = hp->XMode(fSumOption); 
      hp->SumPointError(i, xMode, false, quad, exl, exh, wyl, wyh);
      if (hp->TestBit(kOnlyWeightBit)) continue;
      hp->SumPointError(i, xMode, false, quad, exl, exh, eyl, eyh);
    }
    if (quad && !nosqrt) {
      eyl = TMath::Sqrt(eyl);
      eyh = TMath::Sqrt(eyh);
      wyl = TMath::Sqrt(wyl);
      wyh = TMath::Sqrt(wyh);
    }
    return y;
  }
  static Double_t Round(Double_t v, Int_t p, Int_t& rexpo)
  {
    if (v == 0) {
      // ret.Form("%.*f", p-1, 0);
      rexpo = 0;
      return 0; // ret.Data();
    }

    // Get the sign, take absolute value, get exponent, get scalar, scale 
    Bool_t   neg  = (v < 0);
    Double_t tmp  = TMath::Abs(v);
    Int_t    expo = Int_t(TMath::Log10(tmp));
    Double_t tens = TMath::Power(10, expo - p + 1);
    Double_t n    = RoundN(tens, tmp); // TMath::Ceil(tmp/tens);

    if (n < TMath::Power(10, p-1)) {
      // If scaled is less than 10 to our precision, take off one digit
      expo--;
      tens = TMath::Power(10, expo - p + 1);
      n    = RoundN(tens,tmp); // TMath::Ceil(tmp/tens);
    }

    if (TMath::Abs((n+1) * tens) <= TMath::Abs(n * tens)) 
      // If the value is not near original, add one
      n++;

    if (n >= TMath::Power(10,p)) {
      // If value is larger than 10 to our precision, scale down one
      n /= 10;
      expo++;
    }

    rexpo = expo-p+1;
    return (neg ? -1 : 1)*n;
  }
protected:
  static void SwapPoints(Graph* g, Int_t i, Int_t j, Bool_t reflect)
  {
    Double_t s        = (reflect ? -1 : 1);
    Double_t tmpX     = g->GetX()[i];
    Double_t tmpY     = g->GetY()[i];
    Double_t tmpExl   = g->GetErrorXlow(i);
    Double_t tmpExh   = g->GetErrorXhigh(i);
    Double_t tmpEyl   = g->GetErrorYlow(i);
    Double_t tmpEyh   = g->GetErrorYhigh(i);
    if (reflect) {
      // Swap low and high X error when reflecting 
      Double_t tmp = tmpExl;
      tmpExl       = tmpExh;
      tmpExh       = tmp;
    }
    g->SetPoint(i, s*g->GetX()[j], g->GetY()[j]);
    g->SetPointError(i,
         g->GetErrorXlow(j), g->GetErrorXhigh(j),
         g->GetErrorYlow(j), g->GetErrorYhigh(j));
    g->SetPoint(j, s*tmpX, tmpY);
    g->SetPointError(j, tmpExl, tmpExh, tmpEyl, tmpEyh);
  }
  static const char* FormatKey(const char* key)
  {
    TString tmp(Form("*%s:", key));
    return Form("%-12s", tmp.Data());
  }
  static const char* ExtractField(const TString& value, Int_t idx)
  {
    static TString val;
    val = "";
    TObjArray* tokens = value.Tokenize(":");
    Int_t      iVal   = TMath::Min(Int_t(idx),tokens->GetEntriesFast()-1);
    if (iVal < 0) return val;
    TString    tmp    = tokens->At(iVal)->GetName();
    val               = tmp.Strip(TString::kBoth, ' ');
    // tokens->ls();
    // ::Info("", "Selecting %d: %s", iVal, val.Data());
    delete tokens;
    return val.Data();
  }
  static void StoreQual(TList& quals, Int_t idx,
      const char* name, const char* val)
  {
    TObject* oqv = quals.FindObject(name);
    TList*   qv  = 0;
    if (!oqv) {
      // ::Info("StoreQual", "Creating list for qualifier %s", name);
      qv = new TList;
      qv->SetOwner();
      qv->SetName(name);
      quals.Add(qv);
    }
    else qv = static_cast<TList*>(oqv);
    // ::Info("StoreQual", "Storing %d value %s=%s", idx, name, val);
    qv->AddAt(new TObjString(val), idx);
  }
  static void StoreQual(TList& quals, Int_t idx, TObject* q)
  {
    StoreQual(quals, idx, q->GetName(), q->GetTitle());
  }
  std::ostream& ExportKey(std::ostream& out,
        const char*   which,
        Bool_t        alsoKey = true,
        const char*   fill    = "<please fill in>") const
  {
    Bool_t gotit = false;
    if (fMap) { 
      TIter nextK(fMap);
      TObject* obj = 0;
      while ((obj = nextK())) { 
  TString k(obj->GetName());
  if (!k.EqualTo(which)) continue;
  gotit = true;
  if (alsoKey)  out << FormatKey(obj->GetName()); 
  out << obj->GetTitle() << std::endl;
      }
    }
    if (!gotit) {
      if (alsoKey)
  out << FormatKey(which);
      out << fill << std::endl;
    }
    return out;
  }
  std::ostream& ExportHeader(std::ostream& out,
           Bool_t alsoTop=false,
           Bool_t alsoComment=false,
           Int_t  nsign=-1) const
  {
    if (alsoComment) {
      out << "# Generated by GraphSysErr\n"
    << "# See also\n"
    << "# \n"
    << "#  http://hepdata.cedar.ac.uk/resource/sample.input\n"
    << "#  http://hepdata.cedar.ac.uk/submittingdata\n"
    << "# \n"
    << "# Please fill in missing fields\n";
    }
    if (alsoTop) { 
      TDatime      now;
      UserGroup_t* u        = gSystem->GetUserInfo();
      TString      ref      = GetKey("reference");
      const char*  lab      = GetKey("laboratory");
      const char*  accel    = GetKey("accelerator");
      const char*  exper    = GetKey("detector");
      const char*  abs      = GetKey("abstract");
      const char*  months[] = {"JAN","FEB","MAR","APR",
             "MAY","JUN","JUL","AUG",
             "SEP","OCT","NOV","DEC"};
      if (ref.IsNull()) 
  out << "*reference: <journal/archive ref> : <year>" << std::endl;
      else {
  ExportKey(out, "reference", true, "<reference>");
      }
      ExportKey(out, "author",    true, "<Last name of first author in CAPS>");
      ExportKey(out, "doi",       true, "<DOI number>");
      ExportKey(out, "inspireId", true, "<inSpire ID>");
      ExportKey(out, "cdsId",     true, "<CDS ID>");
      ExportKey(out, "durhamId",  true, "<fill in on submission>");
      ExportKey(out, "title",     true, "<Title of paper>");
      
      out << FormatKey("status") << "Encoded " 
    << std::setw(2) << now.GetDay() <<' '<< months[now.GetMonth()-1]<<' '
    << std::setw(4) << now.GetYear() << " by "<< u->fRealName.Data()<<"\n"
    << FormatKey("experiment") << (lab ? lab : "<lab>") << '-' 
    << (accel ? accel : "<accelerator>") << '-' 
    << (exper ? exper : "<experiment>") << '\n'
    << FormatKey("detector") << (exper ? exper : "<experiment>") << '\n'
    << FormatKey("comment")  << (lab ? lab : "<lab>") << '-' 
    << (accel ? accel : "<accelerator>") << ". " 
    << (abs ? abs : "<abstract>") << "\n"
    << std::endl;
    }
    out << "# Start of dataset\n";
    if (nsign >= 0) out << "# Using " << nsign << " significant digits\n";
    out << FormatKey("dataset") << std::endl;
    // out << FormatKey("dscomment") << GetTitle() << std::endl;
    const char* fields[] = { "location",
           "reackey", 
           "obskey",
           "dscomment",
           // "qual",            
           0 };
    const char** pfld = fields;
    while (*pfld) {
      ExportKey(out, *pfld, true);
      pfld++;
    }
    return out;
  }
  static Int_t ExportError(std::ostream& o,
         Double_t low,
         Double_t high,
         Bool_t   nopm,
         Bool_t   rel,
         Int_t    nsign)
  {
    Double_t l = low;
    Double_t h = high;
    Int_t    p = o.precision();
    Int_t    e = 0;
    if (nsign >= 0) {
      Int_t    pp = nsign;
      Int_t    le, he;
      Double_t lm = Round(l, nsign, le);
      Double_t hm = Round(h, nsign, he);
      e           = TMath::Min(le,he);
      l           = lm*TMath::Power(10,le);
      h           = hm*TMath::Power(10,he);
      o.precision(pp);
    }
    // ::Info("ExportError","Before (%2d) -%8f.+%8f After -%8f,+%8f",
    //        nsign, low,high,l,h);
    if (TMath::Abs(h-l) < 1e-12)
      o << (nopm ? "" : "+- ") << TMath::Abs(l) << (rel ? " PCT" : "");
    else
      o << "+" << h << ",-" << l << (rel ? " PCT" : "");
    o.precision(p);
    return e;
  }
  std::ostream& ExportPoint(std::ostream& out,
          Int_t         i,
          Bool_t        alsoX=true,
          Bool_t        sysName=true,
          Int_t         nsign=0) const
  {
    Int_t p = out.precision();
    // out.setf(std::ios::fixed);
    if (alsoX) {
      Double_t x       = GetX(i);
      Double_t exl     = fData->GetErrorXlow(i);
      Double_t exh     = fData->GetErrorXhigh(i);
      Int_t    pp      = p;
      Int_t    ppp     = p;
      if (nsign > 0) {
  Int_t    pxl, pxh;
  Double_t mxl   = Round(exl,nsign,pxl);
  Double_t mxh   = Round(exh,nsign,pxh);
  exl            = mxl*TMath::Power(10, pxl);
  exh            = mxh*TMath::Power(10, pxh);
  Int_t    mp    = TMath::Min(pxl,pxh);
  pp             = nsign;
  ppp            = TMath::Ceil(TMath::Log10(TMath::Abs(x)))-mp;
  // Info("","pxl=%d pxh=%d mp=%d ppp=%d", pxl, pxh, mp, ppp);
      }
      
      if (TMath::Abs(exl) < 1e-10 && TMath::Abs(exh) < 1e-10) {
  if (nsign <= 0) out << ' ' << x;
  else {
    Int_t    px;
    Double_t mx = Round(x, nsign, px);
    Double_t xx = mx*TMath::Power(10,px);
    // Printf("round %f -> %f (%d) -> %.*g", x, mx, px, nsign, xx);
    out.precision(nsign);
    out << ' ' << xx;
  }
      }
      else {
  if (TMath::Abs(exh-exl) < 1e-10) {
    out.precision(ppp);
    out << ' ' << x-exl << " TO " << x+exh;
  }
  else {
    out.precision(ppp);
    out << ' ' << x << ' ';
    out.precision(pp);
    out << '+' << exh << ",-" << exl;
  }
      }
      out.precision(p);
      out << "; ";
    }
    // out.unsetf(std::ios::fixed);
  
    Bool_t   statRel = fStatRelative;
    Double_t y       = GetY(i);
    Double_t fy      = (statRel ? (y == 0 ? 0 : 100./y) : 1);
    Double_t eyl     = GetStatErrorDown(i) * fy;
    Double_t eyh     = GetStatErrorUp(i)   * fy;
    // out << y << ' ';
    std::stringstream tmp;
    // tmp.setf(std::ios::fixed);
    Int_t le = ExportError(tmp, eyl, eyh, false, statRel, nsign);
      
    if (fPoint2Point.GetEntries() > 0) {
      tmp  << " (" << std::flush;
      TIter       nextP(&fPoint2Point);
      HolderP2P*  holderP2P = 0;
      Bool_t      first = true;
      while ((holderP2P = static_cast<HolderP2P*>(nextP()))) { 
  Bool_t rel = holderP2P->IsRelative();
  fy         = (rel ? (y == 0 ? 0 : 100./y) : 1);
  if (!first) tmp << ',';
  Double_t esh = holderP2P->GetYUp(i)*fy;
  Double_t esl = holderP2P->GetYDown(i)*fy;
  tmp << "DSYS=";
  le = TMath::Min(le,ExportError(tmp, esl, esh, true, rel, nsign));
  if (sysName) 
    tmp << ':'  << holderP2P->GetTitle() << std::flush;
  first = false;
      }
      tmp << ')';
    }
    if (nsign) {
      Int_t nY = TMath::Ceil(TMath::Log10(TMath::Abs(y)))-le;
      out.precision(nY);
      out << y << ' ';
      out.precision(p);
    }
    else
      out << y << ' ';
    out << tmp.str() << ';';
    return out;
  }
  static Double_t RoundN(Double_t tens, Double_t tmp)
  {
    // Add small number to not loose least digit if a 1
    Double_t n     = TMath::Floor(100*tmp/tens+.00001);
    Int_t    next  = int(n/10) % 10;
    if (next > 5) 
      // Round up 
      return TMath::Ceil(n/100);
    if (next < 5)
      // Round down
      return TMath::Floor(n/100);
    Int_t nnext = int(n) % 10;
    if (nnext != 0)
      // Round up 
      return TMath::Ceil(n/100);
    Int_t last = int(n/100) % 10;
    if (last % 2 == 0)
      // Round down
      return TMath::Floor(n/100);
    // Round to nearest even
    return TMath::Ceil(n/100);
  }  
  static TString& Token(TObjArray* c, UShort_t idx, Bool_t verbose=true)
  {
    static TString empty("");
    if (idx >= c->GetEntriesFast()) { 
      if (verbose) {
  ::Error("Token", "Token %d does not exist in array", idx);
  c->ls();
      }
      return empty;
    }
    return static_cast<TObjString*>(c->At(idx))->String();
  }
  static Bool_t ImportError(const TString& s, 
          Double_t& el, Double_t& eh, Bool_t& rel)
  {
    TString tmp(s);
    // ::Info("ImportError", "Parsing '%s' for errors", s.Data());
    // Allow for PCT (any case) and % as relative markers 
    if (tmp.Contains("PCT", TString::kIgnoreCase) || tmp.Contains("%")) 
      rel = true;
    if (tmp.BeginsWith("+-") || tmp.BeginsWith("-+")) {
      el = eh = tmp.Atof();
      // ::Info("ImportError", "Got symmetric %f,%f", el, eh);
      return true;
    }
    TObjArray* tokens = tmp.Tokenize(",");
    for (Int_t i = 0; i < tokens->GetEntriesFast(); i++) { 
      TString& t = Token(tokens, i);
      if (t.IsNull()) continue;
      t.Remove(TString::kBoth, ' ');
      Double_t v = t.Atof();
      // ::Info("ImportError", "Got token '%s' -> %f", t.Data(), v);
      if (t[0] != '-' && t[0] != '+') {
  // ::Info("ImportError", "First char not + or -: %c", t.Data()[0]);
  el = eh = TMath::Abs(v);
      }
      else if (v < 0) {
  // ::Info("ImportError", "Value %f < 0 -> setting el=%f", v, -v);
  if (el > 0) { 
    Double_t max = TMath::Max(el, -v);
    ::Warning("ImportError", 
        "Lower error already set to %f, chosing the maximum of "
        "(%f,%f) -> %f", el, el, -v, max);
    v = -max;
  }
  el = -v;
      }
      else {
  // ::Info("ImportError", "Value %f >= 0 -> setting eh=%f", v, v);
  if (eh > 0) { 
    Double_t max = TMath::Max(eh, v);
    ::Warning("ImportError", 
        "Lower error already set to %f, chosing the maximum of "
        "(%f,%f) -> %f", eh, eh, v, max);
    v = max;
  }
  eh = v;
      }
    }
    // ::Info("ImportError", "Got a-symmetric %f,%f", el, eh);
    tokens->Delete();
    delete tokens;
    return true;
  }
      
  static Bool_t ImportPoint(const TString& s, 
          Double_t&      v, 
          Double_t&      el, 
          Double_t&      eh, 
          Bool_t&        rel)
  {
    // ::Info("ImportPoint", "s=%s", s.Data());
    rel            = s.Contains("PCT");
    TObjArray* tok = s.Tokenize(" ,");
    if (tok->GetEntriesFast() >= 4 && 
  Token(tok,2).EqualTo("TO", TString::kIgnoreCase) && 
  Token(tok,1).Contains("BIN", TString::kIgnoreCase)) {
      v = Token(tok,0).Atof();
      TString tmp = Token(tok,1);
      TString s1  = tmp.Strip(TString::kBoth, '(');
      s1.ReplaceAll("BIN=","");
      Double_t v1 = s1.Atof();
      tmp         = Token(tok,3);
      s1          = tmp.Strip(TString::kBoth, ')');
      Double_t v2 = s1.Atof();
      el          = (v-v1);
      eh          = (v2-v);            
    }
    else if (tok->GetEntriesFast() >= 3 && 
       Token(tok,1).EqualTo("TO", TString::kIgnoreCase)) {
      Double_t v1 = Token(tok,0).Atof();
      Double_t v2 = Token(tok,2).Atof();
      v           = (v1+v2)/2;
      el          = (v-v1);
      eh          = (v2-v);
    }
    else {
      v  = Token(tok,0).Atof();
      if (tok->GetEntriesFast() >= 2 && 
    (Token(tok,1).BeginsWith("+-") || 
     Token(tok,1).BeginsWith("-+"))) {
  TString t = Token(tok,1);
  if (t.Length() > 2) {
    t.Remove(0,2);
    // ::Info("ImportPoint", "Extract from %s", t.Data());
    el = eh = t.Atof();
  }
  else {
    // ::Info("ImportPoint","Extract from next %s",Token(tok, 2).Data());
    el = eh = Token(tok,2).Atof();
  }
      }
      else {
  for (Int_t i = 1; i < tok->GetEntriesFast(); i++) {
    TString t = Token(tok, i);
    Char_t  m = t[0];
    // ::Info("ImportPoint", "%s -> %c", t.Data(), m);
    if (m != '+' && m != '-') continue;
    if (t.Length() == 1) t = Token(tok, ++i); // Take next
    else                 t.Remove(0,1);       // Remove sign
    // ::Info("ImportPoint", "%s -> %f", t.Data(), t.Atof());
    if (m == '-') el = t.Atof();
    else          eh = t.Atof(); 
  }
      }
    }
    tok->Delete();
    delete tok;
    return true;
  }
  /* @} */
  void SqrtPoint(Graph* g, Int_t i) 
  {
    if (kVerbose & kDraw) 
      printf("%20s Point %3d %f,%f -> ","",
       i,g->GetEYlow()[i],g->GetEYhigh()[i]);

    g->GetEYlow()[i]  = TMath::Sqrt(g->GetEYlow()[i]);
    g->GetEYhigh()[i] = TMath::Sqrt(g->GetEYhigh()[i]);

    if (kVerbose & kDraw)
      Printf("%f,%f", g->GetEYlow()[i], g->GetEYhigh()[i]);
  }
  void SqrtGraph(Graph* g) 
  {
    for (Int_t i = 0; i < g->GetN(); i++) 
      SqrtPoint(g, i);
  }
  /* @} */
  const char* FormatOption(UInt_t opt)
  {
    static TString ret;
    ret = "";
    switch (opt) {
    case kNormal:  ret = "";     break; // Line with ticks
    case kNoTick:  ret = "Z0";  break; // Line with no ticks 
    case kArrow:   ret = ">0";  break; // Linw with arrows 
    case kRect:    ret = "20";  break; // Rectangle w/fill
    case kBox:     ret = "50";  break; // Rectangle w/fill & line 
    case kFill:    ret = "30";  break; // Filled area 
    case kCurve:   ret = "40";  break; // Filled smoothed area 
    case kHat:     ret = "[]0"; break; // Hats 
    case kBar:     ret = "||0"; break; // A bar 
    case kNone:    ret = "XP";  break; // No errors
    case kLine:    ret = "CX";   break;
    case kConnect: ret = "LX";   break;
    }
    return ret.Data();
  }
  void MakeDataGraph(Int_t n)
  {
    fData = new Graph(n);
    fData->SetName(Form("%s_data", fName.Data()));
    fData->SetTitle(GetTitle());
  }
  Int_t FindPoint(Double_t x, Int_t& i1, Int_t& i2, Double_t fac=1) const
  {
    i1 = -1;
    i2 = -1;
    Int_t n = GetN();
    if (x < GetX(0)-fac*GetErrorXLeft(0) ||
  x > GetX(n-1)+fac*GetErrorXRight(n-1)) {
      // If we're out side the coverage, set return to -1
      // Info("FindPoint", "X=%f Out of range [%f,%f]",
      //      x,GetX(0),GetX(GetN()-1));
      i1 = i2 = -1;
      return -2;
    }
    const Double_t tol = 1e-9;
    for (Int_t i=0; i < n; i++) {
      // if (TMath::Abs(GetX(i)-x) < tol) {
      if (x <= GetX(i)+fac*GetErrorXRight(i) &&
    x >= GetX(i)-fac*GetErrorXLeft(i)) {
  // if ((GetX(i)+GetErrorXRight(i)-x)<tol &&
  //  (x - GetX(i) - GetErrorXLeft(i) <tol)) {
  // Found matching point
  i2 = i1 = i;
  break;
      }
      // Info("FindPoint", "Not @ %d (%f not in %f -%f +%f)",
      //      i, x, GetX(i), GetErrorXLeft(i), GetErrorXRight(i));
      if (x > GetX(i)) {
  // X still larger
  i1 = i;
      }
      else {
  // X is lower
  i2 = i;
  break;
      }
    }
    if (i1 == i2)
      // If equal, return the point 
      return i1;
    if (TMath::Abs((GetX(i1)+fac*GetErrorXRight(i1)) -
       (GetX(i2)-fac*GetErrorXLeft(i2))) > 10*tol) {
      // If the two found points are not adjecent, we are in a hole
      // and we indicate that.
      // Info("FindPoint","Found hole at i1=%d(%f+%f) i2=%d(%f-%f) (%f)",
      //      i1, GetX(i1), GetErrorXRight(i1),
      //      i2, GetX(i2), GetErrorXLeft(i2), x);
      // Print("xy");
      i1 = i2 = -1;
      return -2;
    }
    if (i2 > 0) {
      // Otherwise return negative value to indicate we should look at
      // two points.
      // Info("FindPoint","Found i1=%d i2=%d %f in [%f,%f]", i1, i2,
      //      x, GetX(i1)-GetErrorXLeft(i1), GetX(i2)+GetErrorXRight(i2));
      return -1;
    }
    // Return negative two to indicate nothing was found
    return -2;
  
  }
  TMultiGraph* MakeMulti(Option_t* option)
  {
    if (!fData) {
      Warning("MakeMulti", "No data defined");
      return 0;
    }

    // --- Process options -------------------------------------------
    TString opt(option);
    opt.ToUpper();

    Bool_t   cmn      = opt.Contains("COMMON");
    Bool_t   combine  = opt.Contains("COMBINED");
    Bool_t   stat     = opt.Contains("STAT");
    Bool_t   quad     = opt.Contains("QUAD");
    combine           = !opt.Contains("STACK");
    quad              = !opt.Contains("DIRECT");
    Bool_t   split    = opt.Contains("SPLIT");
    Int_t    xpos     = (opt.Contains("WEST") ? -1 : 1);
    Int_t    ypos     = (opt.Contains("MIN") ? -1 : 
       opt.Contains("MAX") ?  1 : 0);

    // --- Some general information ----------------------------------
    Int_t    n       = fData->GetN();
    Double_t dx       = TMath::Max(gStyle->GetErrorX(),0.0001F);
    if (n <= 0) {
      Warning("MakeMulti", "Nothing to do (n=%d)", n);
      return 0;
    }

    // --- Find location for common errors ---------------------------
    Double_t xBase    = 0;
    if (opt.Contains("XBASE=")) {
      Int_t   idx = opt.Index("XBASE=")+6;
      TString tmp = opt(idx,opt.Length()-idx);
      xBase       = tmp.Atof();
    }
    else {
      xBase   = ((xpos < 0 ? 
      fData->GetX()[0]   - fData->GetEXlow()[0]: 
      fData->GetX()[n-1] + fData->GetEXhigh()[n-1])
     + xpos * 1.2 * dx);
    }
    Double_t yBase   = fData->GetY()[0];
    for (Int_t i = 0; i < n; i++) { 
      Double_t y = fData->GetY()[i];
      yBase      = (ypos < 0 ? TMath::Min(yBase, y) : 
        ypos > 0 ? TMath::Max(yBase, y) : yBase + y);
    }
    if (ypos == 0) yBase /= n;

    // --- Create stack and temp collection --------------------------
    TList drawn;
    drawn.SetOwner(false);

    Graph* prev = 0;
    if (!combine) { 
      // If we're not combining, we basically copy the graphs to the
      // draw list, but modify the errors to show progresive larger
      // errors.
      prev = fData;
      if (cmn) { 
  // If we should add common errors, do that first 
  TIter next(&fCommon);
  Holder* h = 0;
  while ((h = static_cast<Holder*>(next()))) {
    Graph* g = h->StackError(prev, (!stat ? prev == fData : false), quad);
    if (!g) continue;

    if (quad) SqrtGraph(g);
    prev = g;
    drawn.AddFirst(g, FormatOption(h->GetDOption()));
  }
      }
      TIter next(&fPoint2Point);
      Holder* h = 0;
      while ((h = static_cast<Holder*>(next()))) {
  Graph* g = h->StackError(prev, (!stat ? prev == fData : false), quad);
  if (!g) continue;
  
  if (quad) SqrtGraph(g);
  prev = g;
  drawn.AddFirst(g, FormatOption(h->GetDOption()));
      }
    }
    else { 
      // Otherwise, we are adding all selected errors togethter 
      Graph* g = static_cast<Graph*>(fData->Clone("error"));
      g->SetTitle(fSumTitle);
      Holder* h = 0;
      for (Int_t i = 0; i < n; i++) {
  if (!stat) {
    g->SetPointError(i, g->GetEXlow()[i],
         g->GetEXhigh()[i],0,0);
  }
  else if (quad) {
    // g->GetEXlow()[i] *= g->GetEXlow()[i];
    g->GetEYlow()[i] *= g->GetEYlow()[i];
    // g->GetEXhigh()[i] *= g->GetEXhigh()[i];
    g->GetEYhigh()[i] *= g->GetEYhigh()[i];
  }
    

  if (cmn) { 
    TIter nextC(&fCommon);
    while ((h = static_cast<Holder*>(nextC()))) {
      h->SumError(g, i, false, quad, fSumOption);
      
    }
  }
  TIter nextP(&fPoint2Point);
  while ((h = static_cast<Holder*>(nextP()))) 
    h->SumError(g, i, false, quad, fSumOption);

  // Printf("Point # %d", i);
  if (quad) SqrtPoint(g, i);
      }
      fSumLine.Copy(*g);
      fSumFill.Copy(*g);
      g->SetMarkerStyle(0);
      g->SetMarkerSize(0);
      g->SetMarkerColor(0);

      drawn.AddFirst(g, FormatOption(fSumOption));
    }
    // --- Show common errors on the side, if requested --------------
    if (!cmn) {
      TIter nextC(&fCommon);
      HolderCommon* hc = 0;
      if (!combine) {
  prev = 0;
  while ((hc = static_cast<HolderCommon*>(nextC()))) {
    Graph* g = hc->BarError(prev, quad && !split, xBase, yBase);
    if (!g) {
      Warning("MakeMulti", "Got no graph for common error %s", 
        hc->GetTitle());
      continue;
    }
      
    if (quad && !split) SqrtGraph(g);
      
    if (split) xBase += xpos * 1.2 * dx;
    if (!prev)
      drawn.AddLast(g, FormatOption(hc->GetDOption()));
    else {
      // Here, we use the underling linked list of the list so
      // that we can add an object before another including the
      // possible options 
      TObjLink* f = drawn.FirstLink();
      TObjLink* p = 0;
      while (f && f->GetObject()) { 
        if (f->GetObject() == prev) {
    if (p) new TObjOptLink(g, p, FormatOption(hc->GetDOption()));
    else   new TObjOptLink(g, FormatOption(hc->GetDOption()));
    break;
        }
        p = f;
        f = f->Next();
      }
    }
    prev = (split ? 0 : g);
  } // while (hc)
      } // !combine
      else {
  // Combine systematics 
  if (fCommon.GetEntries() > 0) {
    Graph* g = new Graph(1);
    g->SetPoint(0, xBase, yBase);
    Double_t exl =