AliPhysics  c7b8e89 (c7b8e89)
GraphSysErr.C
Go to the documentation of this file.
1 
12 # include <TNamed.h>
13 # include <TAttMarker.h>
14 # include <TAttLine.h>
15 # include <TAttFill.h>
16 #ifndef GraphSysErr_C
17 # define GraphSysErr_C
18 # include <TList.h>
19 # ifndef __CINT__
20 # include <TStyle.h>
21 # include <TGraph.h>
22 # include <TGraphErrors.h>
23 # include <TGraphAsymmErrors.h>
24 # include <TMultiGraph.h>
25 # include <TH1.h>
26 # include <TMath.h>
27 # include <TFitResultPtr.h>
28 # include <TF1.h>
29 # include <TROOT.h>
30 # include <TSystem.h>
31 # include <TDatime.h>
32 # include <TList.h>
33 # include <TBrowser.h>
34 # include <TPRegexp.h>
35 # include <TClonesArray.h>
36 # include <TArrayL64.h>
37 # include <TLine.h>
38 # include <TClass.h>
39 # include <TObjString.h>
40 # include <iostream>
41 # include <iomanip>
42 # include <fstream>
43 # include <sstream>
44 // # define TOKEN(C,I) static_cast<TObjString*>(C->At(I))->String()
45 # define INCOMPAT_CMN_AS_QUAL 1
46 # else
47 # include <iosfwd>
48 class TGraph;
49 class TGraphErrors;
50 class TGraphAsymmErrors;
51 class TMultiGraph;
52 class TFitResultPtr;
53 class TF1;
54 class TList;
55 class TBrowser;
56 class TArrayD;
57 class TClonesArray;
58 class TLine;
59 # endif
60 
61 
184 class GraphSysErr : public TNamed, public TAttMarker,
185  public TAttLine, public TAttFill
186 {
187 public:
188  enum {
189  kDraw = 0x1,
190  kImport = 0x2,
191  kExport = 0x4,
192  kRatio = 0x8,
193  kVerbose = 0 // Set to OR of above to enable verbose/debug output
194  };
201  kNormal = 0, // - Line with ticks
202  kNoTick = 1, // Z0 - Line with no ticks
203  kArrow = 2, // >0 - Linw with arrows
204  kRect = 3, // 20 - Rectangle w/fill
205  kBox = 4, // 50 - Rectangle w/fill & line
206  kFill = 5, // 30 - Filled area
207  kCurve = 6, // 40 - Filled smoothed area
208  kHat = 7, // []0 - Hats
209  kBar = 8, // ||0 - A bar
210  kNone = 9, // XP - No errors
211  kLine = 10, // XC
212  kConnect = 11 // XL
213  };
219  enum EChi2Type {
220  kExperimentExperiment, // Ignore errors on both
221  kExperimentModel, // Ignore errors on second
222  kModelModel, // Use errors on both
223  kExperimentTruth // Use errors on first, take second to be exact.
224  };
234  kMax = 0x00001, // Error is largest relative error
235  kMin = 0x00002, // Error is largest relative error
236  kCancel = 0x00004,
237  kDenomRel = 0x00008,
239  };
240  enum {
241  kUsedBit = (1 << 18),
242  kOnlyWeightBit = (1 << 19)
243  };
252  : TNamed(),
253  TAttMarker(1,20,1),
254  TAttLine(1,1,1),
255  TAttFill(0,0),
256  fPoint2Point(),
257  fCommon(),
258  fData(0),
259  fDrawn(0),
260  fCounter(0),
261  fSumFill(0,0),
262  fSumLine(1,1,1),
263  fSumTitle(),
264  fSumOption(0),
265  fCommonSumFill(0,0),
266  fCommonSumLine(1,1,1),
267  fCommonSumTitle(),
268  fCommonSumOption(0),
269  fDataOption(0),
270  fXTitle(""),
271  fYTitle(""),
272  fMap(0),
273  fQualifiers(0),
274  fStatRelative(false)
275  {
276 
277  }
284  : TNamed("sysErrGraph","Data"),
285  TAttMarker(1,20,1),
286  TAttLine(1,1,1),
287  TAttFill(0,0),
288  fPoint2Point(),
289  fCommon(),
290  fData(0),
291  fDrawn(0),
292  fCounter(0),
293  fSumFill(0,0),
294  fSumLine(1,1,1),
295  fSumTitle("Uncorr.Errors"),
296  fSumOption(0),
297  fCommonSumFill(0,0),
298  fCommonSumLine(1,1,1),
299  fCommonSumTitle("Corr.Errors."),
300  fCommonSumOption(0),
301  fDataOption(0),
302  fXTitle(""),
303  fYTitle(""),
304  fMap(0),
305  fQualifiers(0),
306  fStatRelative(false)
307  {
308  fPoint2Point.SetName("point2point");
309  fPoint2Point.SetOwner();
310  fCommon.SetName("common");
311  fCommon.SetOwner();
312  MakeDataGraph(n);
313  }
321  GraphSysErr(const char* name, const char* title, Int_t n=0)
322  : TNamed(name,title),
323  TAttMarker(1,20,1),
324  TAttLine(1,1,1),
325  TAttFill(0,0),
326  fPoint2Point(),
327  fCommon(),
328  fData(0),
329  fDrawn(0),
330  fCounter(0),
331  fSumFill(0,0),
332  fSumLine(1,1,1),
333  fSumTitle("Errors"),
334  fSumOption(0),
335  fCommonSumFill(0,0),
336  fCommonSumLine(1,1,1),
337  fCommonSumTitle("Corr.Errors."),
338  fCommonSumOption(0),
339  fDataOption(0),
340  fXTitle(""),
341  fYTitle(""),
342  fMap(0),
343  fQualifiers(0),
344  fStatRelative(false)
345  {
346  fPoint2Point.SetOwner();
347  fCommon.SetOwner();
348  fCommon.SetName("common");
349  fPoint2Point("point2point");
350  MakeDataGraph(n);
351  }
357  GraphSysErr(const GraphSysErr& other)
358  : TNamed(other),
359  TAttMarker(other),
360  TAttLine(other),
361  TAttFill(other),
362  fPoint2Point(),
363  fCommon(),
364  fData(0),
365  fDrawn(0),
366  fCounter(0),
367  fSumFill(other.fSumFill),
368  fSumLine(other.fSumLine),
369  fSumTitle(other.fSumTitle),
370  fSumOption(other.fSumOption),
375  fDataOption(other.fDataOption),
376  fXTitle(other.fXTitle),
377  fYTitle(other.fYTitle),
378  fMap(0),
379  fQualifiers(0),
380  fStatRelative(false)
381  {
382  if (other.fData)
383  fData = static_cast<Graph*>(other.fData->Clone());
384 
385  fPoint2Point.SetOwner();
386  fCommon.SetOwner();
387  fCommon.SetName(other.GetName());
388  fPoint2Point(other.GetName());
389 
390  TIter nextC(&other.fCommon);
391  HolderCommon* common = 0;
392  while ((common = static_cast<HolderCommon*>(nextC()))) {
393  fCommon.Add(new HolderCommon(*common));
394  fCounter++;
395  }
396 
397  TIter nextP(&other.fPoint2Point);
398  HolderP2P* p2p = 0;
399  while ((p2p = static_cast<HolderP2P*>(nextP()))) {
400  fPoint2Point.Add(new HolderP2P(*p2p));
401  fCounter++;
402  }
403 
404  CopyKeys(&other, "ar");
405  // Print("all");
406  }
410  virtual ~GraphSysErr()
411  {
412  fPoint2Point.Delete();
413  fCommon.Delete();
414  if (fData) delete fData;
415  if (fDrawn) delete fDrawn;
416  if (fMap) delete fMap;
417  if (fQualifiers) delete fQualifiers;
418  }
427  {
428  if (&other == this) return *this;
429 
430  other.TNamed::Copy(*this);
431  other.TAttMarker::Copy(*this);
432  other.TAttLine::Copy(*this);
433  other.TAttFill::Copy(*this);
434  fPoint2Point.Clear();
435  fCommon.Clear();
436 
437  if (fData) delete fData;
438  fData = 0;
439  if (other.fData)
440  fData = static_cast<Graph*>(other.fData->Clone());
441 
442  if (fDrawn) { delete fDrawn; fDrawn = 0; }
443 
444  fCounter = other.fCounter;
445 
446  other.fSumFill.Copy(fSumFill);
447  other.fSumLine.Copy(fSumLine);
448  fSumTitle = other.fSumTitle;
449  fSumOption = other.fSumOption;
450 
451  other.fCommonSumFill.Copy(fCommonSumFill);
452  other.fCommonSumLine.Copy(fCommonSumLine);
455 
456  fXTitle = other.fXTitle;
457  fYTitle = other.fYTitle;
459 
460  TIter nextC(&other.fCommon);
461  HolderCommon* common = 0;
462  while ((common = static_cast<HolderCommon*>(nextC())))
463  fCommon.Add(new HolderCommon(*common));
464 
465  TIter nextP(&other.fPoint2Point);
466  HolderP2P* p2p = 0;
467  while ((p2p = static_cast<HolderP2P*>(nextP())))
468  fPoint2Point.Add(new HolderP2P(*p2p));
469 
470  if (fMap) delete fMap;
471  fMap = 0;
472  if (other.fMap)
473  fMap = static_cast<TList*>(other.fMap->Clone());
474 
475  if (fQualifiers) delete fQualifiers;
476  fQualifiers = 0;
477  if (other.fQualifiers)
478  fQualifiers = static_cast<TList*>(other.fQualifiers->Clone());
479 
480  return *this;
481  }
482  /* @} */
483 
493  virtual void ls(Option_t* option="") const
494  {
495  Print(option);
496  }
502  virtual void Print(Option_t* option="R") const
503  {
504  gROOT->IndentLevel();
505  std::cout << GetName() << ": " << GetTitle() << std::endl;
506 
507  TString opt(option);
508  opt.ToUpper();
509  if (opt.IsNull()) return;
510 
511  Bool_t all = opt.Contains("ALL");
512  Bool_t keys = all || opt.Contains("KEY");
513  Bool_t qual = all || opt.Contains("QUAL");
514  Bool_t sys = all || opt.Contains("SYS");
515  Bool_t cmn = sys || opt.Contains("COMMON");
516  Bool_t p2p = sys || opt.Contains("P2P");
517  Bool_t poi = all || opt.Contains("XY");
518  Bool_t attr = all || opt.Contains("ATTR");
519  Bool_t sum = opt.Contains("SUM");
520 
521  gROOT->IncreaseDirLevel();
522  if (fMap && keys) {
523  gROOT->IndentLevel();
524  std::cout << "Key/value pairs: " << std::endl;
525  gROOT->IncreaseDirLevel();
526  TIter nextKV(fMap);
527  TObject* kv = 0;
528  while ((kv = nextKV())) {
529  gROOT->IndentLevel();
530  std::cout << '"' << kv->GetName() << '"' << "\t"
531  << '"' << kv->GetTitle() << '"' << "\n";
532  }
533  gROOT->IndentLevel();
534  std::cout << "\"XTitle\"\t\"" << fXTitle << "\"\n";
535  gROOT->IndentLevel();
536  std::cout << "\"YTitle\"\t\"" << fYTitle << "\"\n";
537  // fMap->Print(option);
538  gROOT->DecreaseDirLevel();
539  }
540  if (fQualifiers && qual) {
541  gROOT->IndentLevel();
542  std::cout << "Qualifier pairs: " << std::endl;
543  gROOT->IncreaseDirLevel();
544  TIter nextQ(fQualifiers);
545  TObject* q = 0;
546  while ((q = nextQ())) {
547  gROOT->IndentLevel();
548  std::cout << '"' << q->GetName() << '"' << "\t"
549  << '"' << q->GetTitle() << '"' << "\n";
550  }
551  // fQualifiers->ls(option);
552  gROOT->DecreaseDirLevel();
553  }
554  if (attr) {
555  gROOT->IndentLevel();
556  Printf(" [option: %3d line (c/s/w):%3d/%1d/%2d fill (c/s):%3d/%4d marker (c/s/s):%3d/%2d/%f]",
557  fDataOption, GetLineColor(), GetLineStyle(), GetLineWidth(),
558  GetFillColor(), GetFillStyle(),
559  GetMarkerColor(), GetMarkerStyle(), GetMarkerSize());
560  gROOT->IndentLevel();
561  Printf(" [sum title: %s "
562  "option: %3d line (c/s/w):%3d/%1d/%2d fill (c/s):%3d/%4d",
563  fSumTitle.Data(), fSumOption,
564  fSumLine.GetLineColor(),
565  fSumLine.GetLineStyle(),
566  fSumLine.GetLineWidth(),
567  fSumFill.GetFillColor(),
568  fSumFill.GetFillStyle());
569  gROOT->IndentLevel();
570  Printf(" [common sum title: %s "
571  "option: %3d line (c/s/w):%3d/%1d/%2d fill (c/s):%3d/%4d",
572  fCommonSumTitle.Data(),
574  fCommonSumLine.GetLineColor(),
575  fCommonSumLine.GetLineStyle(),
576  fCommonSumLine.GetLineWidth(),
577  fCommonSumFill.GetFillColor(),
578  fCommonSumFill.GetFillStyle());
579  }
580 
581  if (cmn) {
582  gROOT->IndentLevel();
583  std::cout << "Commons: " << std::endl;
584  gROOT->IncreaseDirLevel();
585  TIter nextC(&fCommon);
586  TObject* c = 0;
587  while ((c = nextC())) {
588  gROOT->IndentLevel();
589  c->Print(option);
590  }
591  // fCommon.Print(option);
592  gROOT->DecreaseDirLevel();
593  }
594 
595  if (p2p) {
596  gROOT->IndentLevel();
597  std::cout << "Point-to-point: " << std::endl;
598  gROOT->IncreaseDirLevel();
599  TIter nextP(&fPoint2Point);
600  TObject* p = 0;
601  while ((p = nextP())) {
602  gROOT->IndentLevel();
603  p->Print(option);
604  }
605  // fPoint2Point.Print(option);
606  gROOT->DecreaseDirLevel();
607  }
608  if (poi) {
609  // gROOT->IndentLevel();
610  gROOT->IndentLevel();
611  std::cout << "Points: " << std::endl;
612  gROOT->IncreaseDirLevel();
613  TString hs;
614  hs.Form(" %3s: %9s %9s %9s -> %9s %9s %9s",
615  "#", "X", "-dX", "+dX", "Y", "-stat", "+stat");
616  gROOT->IndentLevel();
617  std::cout << hs;
618  if (sum) {
619  hs.Form(" %9s %9s", "-sys", "+sys");
620  std::cout << hs;
621  }
622  else {
623  TIter nextP(&fPoint2Point);
624  HolderP2P* p = 0;
625  while ((p = static_cast<HolderP2P*>(nextP()))) {
626  Bool_t rel = p->IsRelative();
627  const char* post = (rel ? "%" : " ");
628  TString nm(p->GetTitle());
629  if (nm.Length() > 7) {
630  nm.Remove(4,nm.Length()-4);
631  nm.Append("...");
632  }
633  hs.Form(" -%8s%s +%8s%s", nm.Data(), post, nm.Data(), post);
634  std::cout << hs;
635  }
636  }
637  std::cout << std::endl;
638  for (Int_t i = 0; i < GetN(); i++) {
639  gROOT->IndentLevel();
640  Double_t eyl, eyh, wyl, wyh;
641  // point,stat,common,quad,nosqrt,eyl,eyh
642  Double_t y = GetYandError(i,false,false,true,false,eyl,eyh,wyl,wyh);
643  // Double_t y = GetY(i);
644  TString s;
645  s.Form(" %3d: %+8f -%7f +%7f -> %+8f -%7f +%7f",
646  i, GetX(i), GetErrorXLeft(i), GetErrorXRight(i),
647  y, GetStatErrorDown(i), GetStatErrorUp(i));
648  std::cout << s;
649  if (sum) {
650  TString sSum;
651  sSum.Form(" -%7f +%7f", eyl, eyh);
652  std::cout << sSum;
653  }
654  // if (sum && p2p) std::cout << " [";
655  else if (p2p) {
656  TIter nextP(&fPoint2Point);
657  HolderP2P* p = 0;
658  while ((p = static_cast<HolderP2P*>(nextP()))) {
659  Bool_t rel = p->IsRelative();
660  Double_t fac = (rel ? (TMath::Abs(y) > 1e-9 ? 100/y : 0) : 1);
661  const char* post = (rel ? "%" : " ");
662  s.Form(" -%7f%s +%7f%s",
663  fac * p->GetYDown(i), post,
664  fac * p->GetYUp(i) , post);
665  std::cout << s;
666  }
667  }
668  // if (sum && p2p) std::cout << " ]";
669  std::cout << std::endl;
670  } // for i
671  gROOT->DecreaseDirLevel();
672  } // if poi
673  gROOT->DecreaseDirLevel();
674  } //*MENU*
680  virtual Bool_t IsFolder() const { return true; }
686  virtual void Browse(TBrowser* b)
687  {
688  if (fMap) b->Add(fMap);
689  if (fQualifiers) b->Add(fQualifiers);
690  b->Add(&fCommon);
691  b->Add(&fPoint2Point);
692  if (fData) b->Add(fData);
693  if (fDrawn) b->Add(fDrawn);
694  }
695  /* @} */
791  void Draw(Option_t* option="")
792  {
793  TString opt(option);
794  opt.ToUpper();
795  Bool_t clear = opt.Contains("CLEAR");
796  Bool_t axis = opt.Contains("AXIS");
797  // --- Optionally clear old stack --------------------------------
798  if (clear && fDrawn) {
799  delete fDrawn;
800  fDrawn = 0;
801  }
802 
803  fDrawn = MakeMulti(option);
804  if (!fDrawn) return;
805 
806  fDrawn->Draw(axis ? "A" : "");
807  if (axis) {
808  fDrawn->GetHistogram()->SetXTitle(fXTitle);
809  fDrawn->GetHistogram()->SetYTitle(fYTitle);
810  }
811  } //*MENU*
833  TFitResultPtr Fit(TF1* f1, Option_t* fitOption, Option_t* drawOption,
834  Axis_t min=0, Axis_t max=0)
835  {
836  TString dOpt(drawOption);
837  dOpt.ToUpper();
838  Bool_t clear = dOpt.Contains("CLEAR");
839  Bool_t axis = dOpt.Contains("AXIS");
840 
841  if (clear && fDrawn) {
842  delete fDrawn;
843  fDrawn = 0;
844  }
845 
846  fDrawn = MakeMulti(drawOption);
847  if (!fDrawn ||
848  !fDrawn->GetListOfGraphs() ||
849  !fDrawn->GetListOfGraphs()->First())
850  return TFitResultPtr(-1);
851  // fDrawn->GetListOfGraphs()->ls();
852  Graph* g = static_cast<Graph*>(fDrawn->GetListOfGraphs()->First());
853 
854  TString fOpt(fitOption);
855  fOpt.ToUpper();
856  Bool_t noStore = fOpt.Contains("N");
857  Bool_t noDraw = fOpt.Contains("0");
858  Bool_t range = fOpt.Contains("R");
859  fOpt.Append("N0");
860  TFitResultPtr r = g->Fit(f1, fOpt.Data(), "", min, max);
861 
862  if (!noStore) {
863  Int_t status = r;
864  if (status == 0) {
865  if (min < max && !range) f1->SetRange(min, max);
866  fDrawn->GetListOfFunctions()->Add(f1);
867  }
868  }
869 
870  if (!noStore && !noDraw) {
871  fDrawn->Draw(axis ? "A" : "");
872  if (axis) {
873  fDrawn->GetHistogram()->SetXTitle(fXTitle);
874  fDrawn->GetHistogram()->SetYTitle(fYTitle);
875  }
876  }
877 
878  return r;
879  }
901  TFitResultPtr Fit(const char* formula,
902  Option_t* fitOption, Option_t* drawOption,
903  Axis_t min=0, Axis_t max=0)
904  {
905  TString fname(formula);
906  Bool_t linear = fname.Contains("++");
907  TF1* f1 = 0;
908  if (linear) f1 = new TF1(formula,formula,min,max);
909  else {
910  f1 = static_cast<TF1*>(gROOT->GetFunction(formula));
911  if (!f1) Warning("Fit", "Unknown function %s", formula);
912  }
913  if (!f1) return -1;
914 
915  return Fit(f1, fitOption, drawOption, min, max);
916  }
924  TMultiGraph* GetMulti(Option_t* option="")
925  {
926  if (!fDrawn) fDrawn = MakeMulti(option);
927  return fDrawn;
928  }
929  /* @} */
940  void Scale(TF1* f, Double_t s=1)
941  {
942  HolderCommon* cmn = 0;
943  TIter cNext(&fCommon);
944  while ((cmn = static_cast<HolderCommon*>(cNext()))) {
945  if (cmn->IsRelative())
946  // Do not scale relative errors - they scale by themselves
947  continue;
948  cmn->Set(cmn->GetYDown()/s, cmn->GetYUp()/s);
949  }
950  for (Int_t i = 0; i < GetN(); i++) {
951  Double_t x = GetX(i);
952  Double_t y = GetY(i);
953  Double_t g = f->Eval(x);
954  SetPoint(i, GetX(i), y/g);
956 
957  HolderP2P* p2p = 0;
958  TIter pNext(&fPoint2Point);
959  while ((p2p = static_cast<HolderP2P*>(pNext()))) {
960  Double_t fac = (p2p->IsRelative() ?
961  (TMath::Abs(y) < 1e-9) ? 0 : 1/y : 1/g);
962  SetSysError(p2p->GetUniqueID(), i,
963  p2p->GetXLeft(i), p2p->GetXRight(i),
964  fac*p2p->GetYDown(i), fac*p2p->GetYUp(i));
965  }
966  }
967  }
973  void Scale(Double_t s)
974  {
975  HolderCommon* cmn = 0;
976  TIter cNext(&fCommon);
977  while ((cmn = static_cast<HolderCommon*>(cNext()))) {
978  if (cmn->IsRelative())
979  // Do not scale relative errors - they scale by themselves
980  continue;
981  cmn->Set(s * cmn->GetYDown(), s * cmn->GetYUp());
982  }
983  for (Int_t i = 0; i < GetN(); i++) {
984  Double_t y = GetY(i);
985  SetPoint(i, GetX(i), s*y);
987 
988  HolderP2P* p2p = 0;
989  TIter pNext(&fPoint2Point);
990  while ((p2p = static_cast<HolderP2P*>(pNext()))) {
991  Double_t fac = (p2p->IsRelative() ?
992  (TMath::Abs(y) < 1e-9) ? 0 : 1/y :
993  s);
994  SetSysError(p2p->GetUniqueID(), i,
995  p2p->GetXLeft(i), p2p->GetXRight(i),
996  fac*p2p->GetYDown(i), fac*p2p->GetYUp(i));
997  }
998  }
999  }
1013  Int_t Average(const TCollection* others, Bool_t sep=true)
1014  {
1015  const Double_t tol = 1e-9;
1016  TList sharedC; sharedC.SetOwner();
1017  TList sharedP;
1018  TIter oNext(others);
1019  GraphSysErr* other = 0;
1020  Double_t xMin = +1e9;
1021  Double_t xMax = -1e9;
1022  Int_t nX = 0;
1023  Bool_t first = true;
1024  Bool_t nonSh = false;
1025  Int_t npS = 0;
1026  // Here, we do some initial investigations into the the passed
1027  // graphs. We should find all common errors that are shared amoung
1028  // the input graphs, and we need to check for the consistency of
1029  // those.
1030  while ((other = static_cast<GraphSysErr*>(oNext()))) {
1031  Int_t n = other->GetN();
1032  nX = TMath::Max(nX,n);
1033  xMin = TMath::Min(xMin,other->GetX(0) -2*other->GetErrorXLeft(0));
1034  xMax = TMath::Max(xMax,other->GetX(n-1)+2*other->GetErrorXRight(n-1));
1035 
1036  if (first) {
1037  // For the first graph, we simply add the title of all the
1038  // common systematics.
1039  HolderCommon* cmn = 0;
1040  TIter cNext(&other->fCommon);
1041  while ((cmn = static_cast<HolderCommon*>(cNext()))) {
1042  HolderCommon* o = static_cast<HolderCommon*>(cmn->Clone());
1043  sharedC.Add(o);
1044  }
1045  if (sep) {
1046  // and then the same for the point to point
1047  HolderP2P* p2p = 0;
1048  TIter pNext(&other->fPoint2Point);
1049  while ((p2p = static_cast<HolderP2P*>(pNext()))) {
1050  sharedP.Add(p2p->Clone());
1051  }
1052  }
1053 
1054  // We also copy the keys
1055  CopyKeys(other, "ar");
1056  }
1057  else {
1058  // Otherwise, we check that this graph has all the common
1059  // errors in our temporary list, and if not, we remove that
1060  // error from the list
1061  TObjLink* cur = sharedC.FirstLink();
1062  while (cur) {
1063  HolderCommon* o = static_cast<HolderCommon*>(cur->GetObject());
1064  HolderCommon* c = other->FindCompat(o, tol, true);
1065  if (!c) {
1066  // This error not found in the graph we're looking at.
1067  // Remove it from the list
1068  Info("Average", "Common systematic %s not found in %s",
1069  o->GetTitle(), other->GetName());
1070  other->Print("sys");
1071  TObjLink* keep = cur->Next();
1072  TObject* obj = sharedC.Remove(cur);
1073  cur = keep;
1074  nonSh = true;
1075  if (obj) delete obj;
1076  continue;
1077  }
1078  cur = cur->Next();
1079  } // while cur
1080 
1081  // Now check the point-to-point errors
1082  cur = sharedP.FirstLink();
1083  while (cur) {
1084  HolderP2P* o = static_cast<HolderP2P*>(cur->GetObject());
1085  HolderP2P* p = other->FindCompat(o);
1086  if (!p) {
1087  // This error not found in the graph we're looking at.
1088  // Remove it from the list
1089  Info("Average", "P2P systematic %s not found in %s",
1090  o->GetTitle(), other->GetName());
1091  TObjLink* keep = cur->Next();
1092  TObject* obj = sharedP.Remove(cur);
1093  cur = keep;
1094  npS++;
1095  if (obj) delete obj;
1096  continue;
1097  }
1098  cur = cur->Next();
1099  } // while cur
1100  } // else
1101  first = false;
1102  } // while cmn
1103  // Now we have the list of common errors that are shared amoung
1104  // all the input graphs. So we define those new common errors on
1105  // this graph, and disable them on the inputs.
1106  TIter nextC(&sharedC);
1107  HolderCommon* com = 0;
1108  while ((com = static_cast<HolderCommon*>(nextC()))) {
1109  Int_t cId = DefineCommon(com->GetTitle(), com->IsRelative(),
1110  com->GetYDown(1), com->GetYUp(1), kBox);
1111  HolderCommon* c = FindCommon(cId);
1112  c->CopyAttr(com);
1113  // Info("Average", "Shared common error %s +%7f%s -%7f%s",
1114  // c->GetTitle(),
1115  // c->GetYDown(1), (c->IsRelative() ? "%" : ""),
1116  // c->GetYUp(1), (c->IsRelative() ? "%" : ""));
1117 
1118  // Now loop over all input graphs, find the shared commons, and
1119  // mark them as used.
1120  oNext.Reset();
1121  while ((other = static_cast<GraphSysErr*>(oNext()))) {
1122  HolderCommon* oc = other->FindCompat(c, tol, true);
1123  if (!oc) {
1124  Error("Average", "This should NOT happen");
1125  other->Print("sys");
1126  continue;
1127  }
1128  oc->SetBit(kUsedBit);
1129  } // while other
1130  } // while com
1131  sharedC.Clear();
1132  // If we found some common errors that are not shared, we need to
1133  // make a point-to-point error for absorbing the remaining common
1134  // errors.
1135  Int_t bId = 0;
1136  if (nonSh || !sep) {
1137  bId = DeclarePoint2Point("Blob", false, kBox);
1138  SetSysFillColor(bId, kRed+2);
1139  SetSysLineColor(bId, kRed+2);
1140  SetSysFillStyle(bId, 3004);
1141  }
1142 
1143  // At this point, we have found all the common shared errors and
1144  // added those to ourselves. Those that are not shared go into the
1145  // blob point-to-point error.
1146  //
1147  // Next step, is to deal with the shared point-to-point errors.
1148  // If we have shared point-to-point errors, we need to deal with
1149  // them separately. To do that, we make a list of identifiers
1150  // that we should inspect at each X value. The identifiers point
1151  // to the p2p errors of each input graph.
1152  //
1153  TArrayL64 mp2p;
1154  Int_t dOf = 6;
1155  if (sep) {
1156  TIter nextP(&sharedP);
1157  HolderP2P* p2p = 0;
1158  Int_t j = 0;
1159  mp2p.Set(100);
1160  while ((p2p = static_cast<HolderP2P*>(nextP()))) {
1161  // Info("Average", "Loop %d shared P2P %s", j, p2p->GetTitle());
1162  Int_t pId = DeclarePoint2Point(p2p->GetTitle(),
1163  p2p->IsRelative(), kBox);
1164  HolderP2P* p = FindP2P(pId);
1165  p->CopyAttr(p2p);
1166  Info("Average", "Shared point-to-point error %s", p->GetTitle());
1167 
1168  // Now loop over all input graphs, find the shared commons, and
1169  // mark them as used.
1170  oNext.Reset();
1171  Long64_t mask = 0;
1172  Int_t off = 0;
1173  while ((other = static_cast<GraphSysErr*>(oNext())) && off < 64) {
1174  HolderP2P* op = other->FindCompat(p, true);
1175  if (!op) {
1176  Error("Average", "This should NOT happen");
1177  continue;
1178  }
1179  mask |= ((op->GetUniqueID() & 0x3F) << off);
1180  off += dOf;
1181  // Info("Average", "Mark %s as used", op->GetTitle());
1182  op->SetBit(kUsedBit);
1183  } // while other
1184  if (off >= 64) {
1185  Warning("Average",
1186  "Some shared point-to-point errors could not be encoded "
1187  "becasue we have too many (%d>%d) input graphs",
1188  others->GetEntries(), 64/dOf);
1189  }
1190  mp2p[pId] = mask;
1191  j++;
1192  } // while p2p
1193  // Info("Average", "Declared new shared P2P - mapping:");
1194  // for (Int_t j = 1; j < mp2p.GetSize(); j++) {
1195  // if (mp2p[j] == 0) continue;
1196  // rintf("%3d/%3d", j, mp2p.GetSize());
1197  // for (Int_t k = 0; k < 64; k += dOf)
1198  // printf(" %3d", (mp2p[j] >> k) & 0x3f);
1199  // printf("\n");
1200  // }
1201 
1202  // sharedP.Clear();
1203  // Now we have all shared point-to-point errors defined, and we
1204  // have the ids for each graph in the shar array (6bits each).
1205  //
1206  // Next, we will declare _all_ other our point-to-point errors
1207  // separately on this graph. That means that the point-to-point
1208  // errors can be 0 in some ranges.
1209  oNext.Reset();
1210  Int_t off = 0;
1211  while ((other = static_cast<GraphSysErr*>(oNext()))) {
1212  Long64_t* mask = 0;
1213  TIter oNextP(&(other->fPoint2Point));
1214  HolderP2P* p2p = 0;
1215  while ((p2p = static_cast<HolderP2P*>(oNextP())) && off < 64) {
1216  if (p2p->TestBit(kUsedBit)) {
1217  // Info("Average", "%s marked as used", p2p->GetTitle());
1218  continue; // Already in the shared part
1219  }
1220 
1221  // The point-to-point error should go into the weight, but not
1222  // in the result.
1223  p2p->SetBit(kOnlyWeightBit);
1224  Int_t pId = 0;
1225  HolderP2P* p = FindCompat(p2p);
1226  if (!p) {
1227  pId = DeclarePoint2Point(p2p->GetTitle(), p2p->IsRelative(), kBox);
1228  p = FindP2P(pId);
1229  p->CopyAttr(p2p);
1230  mask = &(mp2p[pId]);
1231 
1232  }
1233  *mask |= ((p2p->GetUniqueID() & 0x3F) << off);
1234  }
1235  off += dOf;
1236  } // while other
1237  if (off >= 64) {
1238  Warning("Average",
1239  "Some shared point-to-point errors could not be encoded "
1240  "becasue we have too many (%d>%d) input graphs",
1241  others->GetEntries(), 64/dOf);
1242  }
1243  // Info("Average", "Declared new shared P2P - mapping:");
1244  // for (Int_t j = 0; j < mp2p.GetSize(); j++) {
1245  // if (mp2p[j] == 0) continue;
1246  // printf("%3d/%3d", j, mp2p.GetSize());
1247  // for (Int_t k = 0; k < 64; k += dOf)
1248  // printf(" %3d", (mp2p[j] >> k) & 0x3f);
1249  // printf("\n");
1250  // }
1251  }
1252 
1253  // Now search for the X values to evaluate at
1254  Double_t x = xMin;
1255  Double_t oldX = xMin;
1256  Int_t cnt = 0;
1257  TArrayD xa(2*nX); // Assume no more than 100 points
1258  TArrayD exla(2*nX);
1259  TArrayD exha(2*nX);
1260  do {
1261  // Find the next X to evaluate at
1262  oNext.Reset();
1263  Double_t nextX = xMax;
1264  Double_t texl = 0;
1265  Double_t texh = 0;
1266  // Info("Average", "nextX=%f xMax=%f", nextX, xMax);
1267  while ((other = static_cast<GraphSysErr*>(oNext()))) {
1268  for (Int_t i = other->GetUniqueID(); i < other->GetN(); i++) {
1269  if (other->GetX(i) <= x) {
1270  // Info("Average", "- x_i=%f < x=%f", g->GetX(i), x);
1271  continue;
1272  }
1273  if (other->GetX(i) < nextX) {
1274  // Info("Average", "Set x to %f (was %f)", g->GetX(i), nextX);
1275  nextX = other->GetX(i);
1276  texl = TMath::Max(texl, other->GetErrorXLeft(i));
1277  texh = TMath::Max(texh, other->GetErrorXRight(i));
1278  }
1279  else {
1280  // g->SetUniqueID(i);
1281  // Info("Average", "No more (%s) here at x=%f (next=%f)",
1282  // g->GetName(), g->GetX(i), nextX);
1283  break;
1284  }
1285  }
1286  }
1287  if (nextX >= xMax) {
1288  // Info("Average", "Next x is too large %f>%f", nextX, xMax);
1289  break;
1290  }
1291 
1292  if (cnt >= xa.GetSize()) {
1293  Warning("Average", "increasing size of X cache");
1294  xa.Set(1.5*cnt);
1295  exla.Set(1.5*cnt);
1296  exha.Set(1.5*cnt);
1297  }
1298  // Info("Average", "Set next x to %f +/- (%f,%f)", nextX, texl, texh);
1299  xa[cnt] = nextX;
1300  exla[cnt] = texl;
1301  exha[cnt] = texh;
1302  Double_t dx = texh*2;//(nextX-oldX);
1303  oldX = x;
1304  x = nextX+dx/4;
1305  cnt++;
1306  } while (cnt < 1000); // Do not go for more than 1000 points
1307  // At this point we have
1308  //
1309  // - All common errors defined, and we have set the shared ones
1310  // - All point-to-point errors defined, and we know which input
1311  // graph contributes to each.
1312  // - A list of X values to evaluate each graph at.
1313 
1314  // Print("sys");
1315  // Now we loop over the list of X values and evaluate each graph
1316  // at that X value.
1317  for (Int_t i = 0; i < cnt; i++) {
1318  Double_t xi = xa[i];
1319  Double_t texl = exla[i];
1320  Double_t texh = exha[i];
1321  Double_t nexl = TMath::Abs(xi - (i == 0 ? xa[i+1] : xa[i-1]))/2;
1322  Double_t nexh = TMath::Abs(xi - (i + 1 == cnt ? xa[i-1] : xa[i+1]))/2;
1323  Double_t exl = TMath::Min(texl, nexl);
1324  Double_t exh = TMath::Min(texh, nexh);
1325  // Info("Average", "@ %3d x=%f, exl=min(%f,%f)=%f exh=min(%f,%f)=%f",
1326  // i, xi, texl, nexl, exl, texh, nexh, exh);
1327 
1328 
1329  // Now loop again - this time calculating our data point
1330  Double_t sum = 0;
1331  Double_t sumW = 0;
1332  // Double_t sumWS = 0;
1333  Double_t sumE = 0;
1334  Double_t sumES = 0;
1335  oNext.Reset();
1336 
1337  // Use our combiner framework for calculating the averaged
1338  // result and errors. Note, we define two combiners. One that
1339  // include all non-shared common and point-to-point errors
1340  // (centroid), and one that doesn't have the non-shared common
1341  // and point-to-point errors as well as the split point-to-point
1342  // errors (error). The first one is used to set the centroid of
1343  // the data, while the second is used to set the combined error
1344  LinearSigmaCombiner centroid;
1345  LinearSigmaCombiner error;
1346  LinearSigmaCombiner stats;
1347  Int_t j = 0;
1348  TArrayI ai1(others->GetEntries());
1349  TArrayI ai2(others->GetEntries());
1350  TArrayD ax(others->GetEntries());
1351  while ((other = static_cast<GraphSysErr*>(oNext()))) {
1352  Bool_t cmn = true; // Include the common errors
1353  Bool_t stat = false; // Include the statistical errors
1354  Bool_t quad = true; // Add in quadrature
1355  Bool_t nosqrt = false; // Do not return squared error
1356 
1357  Double_t y, eyl, eyh, wyl, wyh, syl, syh;
1358 
1359  Int_t fret = other->FindPoint(xi, ai1[j], ai2[j]);
1360  if (fret < -1) { j++; continue; } // Nothing here
1361  if (fret >= 0) {
1362  y = other->GetYandError(fret,cmn,stat,quad,nosqrt,eyl, eyh, wyl, wyh);
1363  syl = other->GetStatErrorDown(fret);
1364  syh = other->GetStatErrorUp(fret);
1365  }
1366  else {
1367  Double_t eyl1, eyh1, wyl1, wyh1;
1368  Double_t eyl2, eyh2, wyl2, wyh2;
1369  Double_t x1 = other->GetX(ai1[j]);
1370  Double_t x2 = other->GetX(ai2[j]);
1371  Double_t y1 = other->GetYandError(ai1[j],cmn,stat,quad,nosqrt,
1372  eyl1,eyh1,wyl1,wyh1);
1373  Double_t y2 = other->GetYandError(ai2[j],cmn,stat,quad,nosqrt,
1374  eyl2,eyh2,wyl2,wyh2);
1375  Double_t syl1 = other->GetStatErrorDown(ai1[j]);
1376  Double_t syh1 = other->GetStatErrorUp(ai1[j]);
1377  Double_t syl2 = other->GetStatErrorDown(ai2[j]);
1378  Double_t syh2 = other->GetStatErrorUp(ai2[j]);
1379  Double_t dx = (x2-x1);
1380  ax[j] = (xi-x1)/dx;
1381  y = y1 + ax[j] * (y2-y1);
1382  eyl = eyl1 + ax[j] * (eyh1-eyl1);
1383  eyh = eyh2 + ax[j] * (eyh2-eyl2);
1384  wyl = wyl1 + ax[j] * (wyh1-wyl1);
1385  wyh = wyh2 + ax[j] * (wyh2-wyl2);
1386  syl = syl1 + ax[j] * (syh1-syl1);
1387  syh = syh2 + ax[j] * (syh2-syl2);
1388  }
1389 
1390  centroid.Add(y, wyl, wyh);
1391  error.Add(y, eyl, eyh);
1392  stats.Add(y, syl, syh);
1393  j++;
1394  } // while (other)
1395  centroid.Calculate();
1396  // Info("Average", "Centroid");
1397  // centroid.Print();
1398  // centroid.fResult->Print();
1399  Double_t newY = centroid.fResult->fX;
1400  SetPoint(i, xi, newY);
1401  SetPointError(i, exl, exh);
1402 
1403  stats.Calculate();
1404  // Info("Average", "Stats");
1405  // stats.Print();
1406  // stats.fResult->Print();
1407  stats.Calculate();
1408  Double_t fac = IsStatRelative() ? 1/newY : 1;
1409  SetStatError(i, fac*stats.fResult->fEl,fac*stats.fResult->fEh);
1410 
1411  if (bId > 0) {
1412  // Calculations for error
1413  error.Calculate();
1414  // Info("Average", "Error");
1415  // error.Print();
1416  // error.fResult->Print();
1417 
1418  fac = IsRelative(bId) ? 1/newY : 1;
1419  SetSysError(bId,i,0,0,fac*error.fResult->fEl,fac*error.fResult->fEh);
1420  }
1421 
1422  // Now, we have set the blob point-to-point error. Next, we
1423  // need to do the point-to-point erros that should be transfered
1424  // to the output. We have stored masks of IDs in mp2p, so we
1425  // loop over that array first. The lower 6 bits is the output
1426  // Id, while the remaining bits (counting from the low side) are
1427  // ids for each input graph.
1428  for (Int_t j = 1; j < mp2p.GetSize(); j++) {
1429  Long64_t mask = mp2p[j];
1430  if (mask == 0) continue;
1431  Int_t tid = j; // (mask & 0x3f);
1432  HolderP2P* tp2p = FindP2P(tid);
1433  if (!tp2p) {
1434  Warning("Average", "No target point-to-point error at %d for id=%d",
1435  j, tid);
1436  continue;
1437  }
1438  // Zero by default
1439  SetSysError(tid, i, 0, 0);
1440  Long64_t rem = mask;
1441  if (rem == 0)
1442  // There's no IDs for this, explicitly set the value to 0
1443  continue;
1444 
1445  // Now we should loop over all the input graphs and get the
1446  // error and finally combine them into one. Again, we use a
1447  // linear sigma approximation from our combiner framework to
1448  // do the calculations.
1450  Int_t k = 0; // Index into index array and for off set in mask
1451  oNext.Reset();
1452  while ((other = static_cast<GraphSysErr*>(oNext()))) {
1453  if (ai1[k] == -1 && ai2[k] == -1) {
1454  // This graph does not contribute at this point
1455  // Info("Average", "No contribution at %f (%d:%d,%d) for %s to %s",
1456  // xi, k, ai1[k], ai2[k], other->GetName(), tp2p->GetTitle());
1457  k++;
1458  continue;
1459  }
1460  Int_t sid = ((rem >> dOf*k) & 0x3f);
1461  if (sid <= 0) {
1462  // This graph does not contribute to this error
1463  // Info("Average",
1464  // "No contribution for %s to %s, not found [%d:%d,%d]",
1465  // other->GetName(), tp2p->GetTitle(), k, ai1[k], ai2[k]);
1466  k++;
1467  continue;
1468  }
1469  HolderP2P* sp2p = other->FindP2P(sid);
1470  if (!sp2p) {
1471  Warning("Average", "Couldn't find point-to-point error %d (%s)"
1472  "in graph %d (%s) - should not happen",
1473  sid, tp2p->GetTitle(), k, other->GetName());
1474  k++;
1475  }
1476  Double_t sy = 0;
1477  Double_t seyl = 0;
1478  Double_t seyh = 0;
1479  if (ai1[k] == ai2[k]) {
1480  // Exact point
1481  sy = other->GetY(ai1[k]);
1482  seyl = other->GetSysErrorYDown(sid, ai1[k]);
1483  seyh = other->GetSysErrorYUp(sid, ai1[k]);
1484  }
1485  else {
1486  // Interpolate
1487  Double_t sy1 = other->GetY(ai1[k]);
1488  Double_t sy2 = other->GetY(ai2[k]);
1489  sy = sy1 + ax[k] * (sy2-sy1);
1490  seyl = sp2p->GetYDown(ai1[k],ai2[k],ax[k]);
1491  seyh = sp2p->GetYUp(ai1[k],ai2[k],ax[k]);
1492  }
1493  //Info("Average", "Contribution (%d:%d-%d:%f,%f) to %s from %s",
1494  // k,ai1[k],ai2[k],seyl,seyh,tp2p->GetTitle(),other->GetName());
1495  sc.Add(sy, seyl, seyh);
1496  k++;
1497  } // while other
1498  sc.Calculate();
1499  // Info("Average", "%s", tp2p->GetTitle());
1500  // sc.Print();
1501  // sc.fResult->Print();
1502  fac = IsRelative(tid) ? 1/newY : 1;
1503  SetSysError(tid, i, 0, 0, fac*sc.fResult->fEl, fac*sc.fResult->fEh);
1504  } // for j (# of p2p)
1505  } // for i (points)
1506 
1507  oNext.Reset();
1508  while ((other = static_cast<GraphSysErr*>(oNext())))
1509  other->ClearUsed();
1510  ClearUsed();
1511 
1512  return bId;
1513  }
1525  Option_t* option="quad sum stat",
1526  UShort_t first=0,
1527  Short_t last=-1)
1528  {
1529  if (last < 0) last = GetN()-1;
1530  Double_t sum = 0;
1531  Double_t err2 = 0;
1532  TString opts(option); opts.ToLower();
1533  Bool_t cmn = opts.Contains("comm");
1534  Bool_t stat = opts.Contains("stat");
1535  Bool_t quad = opts.Contains("quad");
1536 
1537  for (Short_t i = first; i <= last; i++) {
1538  Double_t eyl, eyh;
1539  Double_t y = GetYandError(i, cmn, stat, quad, false, eyl, eyh);
1540  Double_t eym = (eyh+eyl)/2;
1541  Double_t x = GetX(i);
1542  Double_t exl = GetErrorXLeft(i);
1543  Double_t exr = GetErrorXRight(i);
1544  Double_t x1 = x - exl;
1545  Double_t x2 = x + exr;
1546  Double_t prv = x1;
1547  Double_t nxt = x2;
1548  if (i-1 >= 0) prv = GetX(i-1)+GetErrorXRight(i-1);
1549  if (i+1 < GetN()) nxt = GetX(i+1)-GetErrorXLeft(i+1);
1550  if (prv > x1) x1 -= (prv-x1)/2; // Adjust to half-ways
1551  if (nxt < x2) x2 -= (x2-nxt)/2; // Adjust to half-ways
1552  Double_t wdt = (x2-x1);
1553  if (wdt < 0) continue; // Do not take this point
1554  sum += wdt * y;
1555  err2 += TMath::Power(eym*wdt,2);
1556  }
1557  error = TMath::Sqrt(err2);
1558  return sum;
1559  }
1577  const GraphSysErr* num,
1578  const GraphSysErr* denom,
1579  Double_t& x,
1580  Double_t& dY,
1581  Double_t& dEyl,
1582  Double_t& dEyh,
1583  Double_t& nY,
1584  Double_t& nEyl,
1585  Double_t& nEyh)
1586  {
1587  x = num->GetX(i);
1588  if (!denom->FindYandError(x, false, true, true, false, dY, dEyl, dEyh))
1589  return false;
1590  if (TMath::Abs(dY) < 1e-9)
1591  return false;
1592 
1593  nY = num->GetYandError(i,false,true,true,false,nEyl, nEyh);
1594  return true;
1595  }
1596 
1608  static GraphSysErr* Ratio(const GraphSysErr* num,
1609  const GraphSysErr* den,
1610  UInt_t flags=kRatioDefault,
1611  Double_t fac=1)
1612  {
1613  // num->Print("sys");
1614  // den->Print("sys");
1615  num->ClearUsed();
1616  den->ClearUsed();
1617  GraphSysErr* ret = new GraphSysErr(1);
1618  ret->CopyAttr(num);
1619  ret->CopyKeys(num, "ar");
1620  Bool_t emax = (flags & kMax);
1621  Bool_t cmin = (flags & kMin);
1622  Bool_t cmax = (flags & kMax);
1623  Bool_t denrel = (flags & kDenomRel);
1624  Bool_t cancc = (flags & kCancel);
1625  Int_t id = 0;
1626  Int_t notUsed = 0;
1627  // Loop over all common errors in the numerator, and see if we
1628  // have the same common error in the denominator. If we do, we
1629  // can (partially) cancel this common error.
1630  //
1631  // Technically, we mark each used common error as used, and create
1632  // a new common error in the output graph.
1633  //
1634  // Common errors that are not shared by both the numerator and the
1635  // denominator will not be cancelled. Instead, if they are
1636  // relative, they will be added to the output graph directly.
1637  // Finally, common systematic errors that are not cancelled or are
1638  // not relative, will go into the blob of point-to-point errors
1639  TIter nextNC(&num->fCommon);
1640  HolderCommon* comN = 0;
1641  while ((comN = static_cast<HolderCommon*>(nextNC()))) {
1642  Bool_t rel = comN->IsRelative();
1643  Double_t eyl = comN->GetYDown(1);
1644  Double_t eyh = comN->GetYUp(1);
1645  Int_t idD = den->FindId(comN->GetTitle());
1646  HolderCommon* comD = 0;
1647  if (idD != 0 && // Check it exists
1648  (comD = den->FindCommon(idD)) && // Check it is a common
1649  (rel == comD->IsRelative())) { // Same as this
1650  // We found the same error in the denominator, so we need to
1651  // cancel errors here.
1652  Double_t dEyl = comD->GetYDown(1);
1653  Double_t dEyh = comD->GetYUp(1);
1654  Double_t nEyl = comN->GetYDown(1);
1655  Double_t nEyh = comN->GetYUp(1);
1656  if (cmin) {
1657  eyl = TMath::Min(nEyl, dEyl);
1658  eyh = TMath::Min(nEyh, dEyh);
1659  }
1660  else if (cmax) {
1661  eyl = TMath::Max(nEyl, dEyl);
1662  eyh = TMath::Max(nEyh, dEyh);
1663  }
1664  else if (cancc) {
1665  eyl = TMath::Sqrt(TMath::Abs(dEyl*dEyl-nEyl*nEyl));
1666  eyh = TMath::Sqrt(TMath::Abs(dEyh*dEyh-nEyh*nEyh));
1667  }
1668  const char* post = (rel ? "%" : "");
1669  Double_t pfc = (rel ? 100 : 1);
1670 
1671  if (kVerbose & kRatio)
1672  ::Info("Ratio", "Cancelling the common systematic error %s "
1673  "between numerator (-%f%s +%f%s) "
1674  "and denominator (-%f%s +%f%s): -%f%s + %f%s",
1675  comD->GetTitle(),
1676  pfc*nEyl, post, pfc*nEyh, post,
1677  pfc*dEyl, post, pfc*dEyh, post,
1678  pfc*eyl, post, pfc*eyh, post);
1679 
1680  // Flag the denominator error as used
1681  comD->SetBit(kUsedBit);
1682  }
1683  else if (!rel) {
1684  // We cannot cancel, and it is not relative, add to blob
1685  if (kVerbose & kRatio)
1686  ::Info("Ratio",
1687  "Common numerator systematic %s will be added to blob",
1688  comN->GetTitle());
1689  notUsed++; // Count how many we're not setting directly
1690  continue;
1691  }
1692  // Now we can define our common error. Which value we set depend
1693  // on whether the error was cancelled with the denominator
1694  Int_t cId = ret->DefineCommon(comN->GetTitle(),rel,
1695  eyl, eyh, kBox);
1696  HolderCommon* c = ret->FindCommon(cId);
1697  c->CopyAttr(comN);
1698  // Flag the numerator error as used
1699  comN->SetBit(kUsedBit);
1700  }
1701  // Loop over all denominator common systematic errors. If it is
1702  // not cancelled with the numerator and is relative, we define a
1703  // new error on the output ratio. Otherwise it will be absorbed in a blob
1704  TIter nextDC(&den->fCommon);
1705  HolderCommon* comD = 0;
1706  while ((comD = static_cast<HolderCommon*>(nextDC()))) {
1707  if (comD->TestBit(kUsedBit)) continue;
1708  if (!comD->IsRelative()) {
1709  if (kVerbose & kRatio)
1710  ::Info("Ratio",
1711  "Common denominator systematic %s will be added to blob",
1712  comD->GetTitle());
1713  notUsed++; // Count how many we're not setting directly
1714  continue;
1715  }
1716  if (kVerbose & kRatio)
1717  ::Info("Ratio", "Propagating common sysmatic %s to ratio",
1718  comD->GetTitle());
1719  // ret->Print("sys");
1720  // Now we can define our common error. Which value we set depend
1721  // on whether the error was cancelled with the denominator
1722  Int_t cId = ret->DefineCommon(comD->GetTitle(),true,
1723  comD->GetYDown(1),
1724  comD->GetYUp(1),
1725  kBox);
1726  HolderCommon* c = ret->FindCommon(cId);
1727  c->CopyAttr(comD);
1728  // Flag the denominator error as used
1729  comD->SetBit(kUsedBit);
1730  // ret->Print("sys");
1731  }
1732 
1733  // Loop over all point-to-point errors in the numerator, and see
1734  // if we have the same point-to-point error in the denominator.
1735  // If we do, we can (partially) cancel this point-to-point error.
1736  //
1737  // Technically, we mark each used point-to-point error as used.
1738  //
1739  // Point-to-point errors that are not shared by both the numerator
1740  // and the denominator will not be cancelled. Instead, they will
1741  // be added in on the residual point-to-point sum (in quadrature)
1742  // error.
1743  //
1744  // We encode the IDs of the
1745  // errors in a single mask.
1746  //
1747  // 24 -------- 16 -------- 8 -------- 0
1748  // | Denom ID | Num ID | Ratio ID |
1749  // +-----------+----------+----------+
1750  //
1751  Int_t sCnt = 0;
1752  TArrayI shared(num->fPoint2Point.GetEntries() +
1753  den->fPoint2Point.GetEntries());
1754  TIter nextNP(&num->fPoint2Point);
1755  HolderP2P* p2pN = 0;
1756  while ((p2pN = static_cast<HolderP2P*>(nextNP()))) {
1757  Bool_t rel = p2pN->IsRelative();
1758  Int_t idD = den->FindId(p2pN->GetTitle());
1759  HolderP2P* p2pD = idD == 0 ? 0 : den->FindP2P(idD);
1760  Bool_t same = !p2pD ? false : rel == p2pD->IsRelative();
1761  Int_t mask = 0;
1762  // Printf("Id of denom %s -> %d (%p, %s, %s)", p2pN->GetTitle(), idD,
1763  // p2pD, same ? "same" : "diff",
1764  // p2pD ? (p2pD->TestBit(kUsedBit) ? "used" : "ready") : "?");
1765  if (cancc &&
1766  p2pD && // Check it is a p2p
1767  same) { // Same as this
1768  // We found the same error in the denominator, so we need to
1769  // cancel errors when looping.
1770  mask = ((p2pD->GetUniqueID() & 0xFF) << 16);
1771  // Printf("Found %s in both denom (%s,%d) and num (%s,%d)",
1772  // p2pN->GetTitle(),
1773  // p2pD->GetTitle(), p2pD->GetUniqueID(),
1774  // p2pN->GetTitle(), p2pD->GetUniqueID());
1775  if (kVerbose & kRatio)
1776  ::Info("Ratio", "Cancelling the p2p systamtic error %s "
1777  "between numerator and denominator",
1778  p2pN->GetTitle());
1779  // Flag the denominator error as used
1780  p2pD->SetBit(kUsedBit);
1781  }
1782  Int_t pId = ret->DeclarePoint2Point(p2pN->GetTitle(),rel,
1783  kBox);
1784  HolderP2P* p = ret->FindP2P(pId);
1785  p->CopyAttr(p2pN);
1786  // Flag the numerator error as used
1787  p2pN->SetBit(kUsedBit);
1788  mask |= (((p2pN->GetUniqueID() & 0xFF) << 8) |
1789  (pId & 0xFF));
1790  // Printf("Found shared P2P: %s -> %d", p2pN->GetTitle(), pId);
1791  // num->Print("sys");
1792  // den->Print("sys");
1793  shared[sCnt] = mask;
1794  sCnt++;
1795  }
1796  // Loop over all denominator point-to-point systematic errors. If
1797  // it is not cancelled with the numerator, we define a new error
1798  // on the output ratio
1799  TIter nextDP(&den->fPoint2Point);
1800  HolderP2P* p2pD = 0;
1801  while ((p2pD = static_cast<HolderP2P*>(nextDP()))) {
1802  if (p2pD->TestBit(kUsedBit)) {
1803  // Printf("P2P of denom %s already used", p2pD->GetTitle());
1804  continue;
1805  }
1806  Int_t pId = ret->DeclarePoint2Point(p2pD->GetTitle(),true,
1807  kBox);
1808  HolderP2P* p = ret->FindP2P(pId);
1809  p->CopyAttr(p2pD);
1810  // Flag the numerator error as used
1811  p2pD->SetBit(kUsedBit);
1812  Int_t mask = (((p2pD->GetUniqueID() & 0xFF) << 16) |
1813  (pId & 0xFF));
1814  // Printf("P2P of %s (%d) denom propagated", p2pD->GetTitle(),
1815  // p2pD->GetUniqueID());
1816  shared[sCnt] = mask;
1817  sCnt++;
1818  }
1819  Int_t bId = 0;
1820  if (notUsed > 0) {
1821  // If we have some error that are not used - i.e., non-relative
1822  // common errors in the numerator or denominator, we need to
1823  // make a new point-to-point systematic error for those.
1824  bId = ret->DeclarePoint2Point("Blob", false, kBox);
1825  }
1826 
1827  // Now loop over all the points we have in the numerator and
1828  // figure out what to do for each point.
1829  Int_t cnt = 0; // Counter of output points
1830  for (Int_t i = 0; i < num->GetN(); i++) {
1831  Double_t x = num->GetX(i);
1832  Double_t nY = num->GetY(i);
1833  Double_t nEyl = 0;
1834  Double_t nEyh = 0;
1835  Double_t dY = 0;
1836  Double_t dEyl = 0;
1837  Double_t dEyh = 0;
1838  Int_t di1, di2;
1839  Int_t di = den->FindPoint(x, di1, di2, fac);
1840  Double_t daX = 0; // Relative distance between interpolated points
1841  if (di < -1) {
1842  if (kVerbose & kRatio)
1843  ::Warning("Ratio", "Next point %d (%f) not found in denominator",
1844  i, x);
1845  continue;
1846  }
1847  Bool_t cmn = true; // Include the common errors
1848  Bool_t stat = false; // Include the statistical errors
1849  Bool_t quad = true; // Add in quadrature
1850  Bool_t nosqrt = false; // Do not return squared error
1851  if (di >= 0) {
1852  // Exact point
1853  // point,common,stat,quadratic,nosqrt,
1854  dY = den->GetYandError(di,cmn,stat,quad,nosqrt,dEyl,dEyh);
1855  di1 = di;
1856  }
1857  else {
1858  Double_t eyl1, eyl2, eyh1, eyh2;
1859  Double_t x1 = den->GetX(di1);
1860  Double_t x2 = den->GetX(di2);
1861  Double_t y1 = den->GetYandError(di1,cmn,stat,quad,nosqrt,eyl1,eyh1);
1862  Double_t y2 = den->GetYandError(di2,cmn,stat,quad,nosqrt,eyl2,eyh2);
1863  // Linear interpolation
1864  Double_t dx = (x2-x1);
1865  daX = (x-x1)/dx;
1866  dY = y1 + daX * (y2 - y1);
1867  dEyl = eyl1 + daX * (eyl2 - eyl1);
1868  dEyh = eyh1 + daX * (eyh2 - eyh1);
1869  }
1870  // At this point, we have the X and Y values, and the non-p2p
1871  // and non-cancelled and non-relative-common errors calculated
1872  // for numerator and denominator individually.
1873  //
1874  // We start by setting the point value. Perhaps here, we should
1875  // also set the statistics.
1876  Double_t rY = (dY != 0 ? nY/dY : 0);
1877  ret->SetPoint(cnt, x, rY);
1878  ret->SetPointError(cnt, num->GetErrorXLeft(i), num->GetErrorXRight(i));
1879 
1880  // We have not dealt with the statistical errors yet. We do
1881  // that here.
1882  Double_t sEyl = 0;
1883  Double_t sEyh = 0;
1884  if (TMath::Abs(rY) > 1.e-9) {
1885  Double_t snEyl = num->GetStatErrorDown(cnt);
1886  Double_t snEyh = num->GetStatErrorUp (cnt);
1887  Double_t sdEyl = den->GetStatErrorDown(cnt);
1888  Double_t sdEyh = den->GetStatErrorUp (cnt);
1889  // Statistical errors never cancel, and must be propagted. Note,
1890  // despite it technically being wrong, we treat the lower and
1891  // higher errors independently.
1892  sEyl = TMath::Sqrt(rY*rY*(snEyl*snEyl/nY/nY+sdEyl*sdEyl/dY/dY));
1893  sEyh = TMath::Sqrt(rY*rY*(snEyh*snEyh/nY/nY+sdEyh*sdEyh/dY/dY));
1894  }
1895  ret->SetStatError(cnt, sEyl, sEyh);
1896 
1897  // Next, we should update our blob systematic p2p error if it
1898  // exists.
1899  if (bId > 0) {
1900  Double_t eyl = 0, eyh = 0;
1901  // Printf("Calculate blob syst.uncer. d=(%f,%f) n=(%f,%f)",
1902  // dEyl, dEyh, nEyl, nEyh);
1903  if (denrel) {
1904  Double_t rdEyl = (dY != 0 ? dEyl/dY : 0);
1905  Double_t rdEyh = (dY != 0 ? dEyh/dY : 0);
1906  eyl = rdEyl*rY;
1907  eyh = rdEyh*rY;
1908  }
1909  else if (emax) {
1910  // Take the largest of the two errors, and divide by the
1911  // denominator.
1912  eyl = TMath::Max(nEyl, dEyl);
1913  eyh = TMath::Max(nEyh, dEyh);
1914  eyl /= dY;
1915  eyh /= dY;
1916  }
1917  else if (cmin) {
1918  // Take the largest of the two errors, and divide by the
1919  // denominator.
1920  eyl = TMath::Min(nEyl, dEyl);
1921  eyh = TMath::Min(nEyh, dEyh);
1922  eyl /= dY;
1923  eyh /= dY;
1924  }
1925  else {
1926  // Find the relative errors and add in quadrature
1927  Double_t rnEyl = (nY != 0 ? nEyl/nY : 0);
1928  Double_t rnEyh = (nY != 0 ? nEyh/nY : 0);
1929  Double_t rdEyl = (dY != 0 ? dEyl/dY : 0);
1930  Double_t rdEyh = (dY != 0 ? dEyh/dY : 0);
1931  Double_t reyl = TMath::Sqrt(rnEyl*rnEyl+rdEyl*rdEyl);
1932  Double_t reyh = TMath::Sqrt(rnEyh*rnEyh+rdEyh*rdEyh);
1933  eyl = reyl*rY;
1934  eyh = reyh*rY;
1935  }
1936  // Printf("den err (%f%%,%f%%) - blob (%f%%,%f%%)", dEyl*100, dEyh*100,
1937  // eyl*100,eyh*100);
1938  ret->SetSysError(bId, cnt, eyl, eyh);
1939  }
1940 
1941  // Now we need to update our various point-to-point
1942  // errors. We've stored a mask for each declared point-to-point
1943  // error. Each mask consist of three fields, corresponding to
1944  // point-to-point errors in the denominator, numerator, and the
1945  // ratio. If all three fields are set, we shold cancel between
1946  // numerator and denominator. If only the ratio and either of
1947  // the denominator or numerator fields are set, we simply need
1948  // to copy that to the ratio point-to-point error - as a
1949  // relative error of course.
1950  for (Int_t im = 0; im < shared.GetSize(); im++) {
1951  Int_t mask = shared[im];
1952  if (mask <= 0) break;
1953  Int_t rId = ((mask >> 0) & 0xFF);
1954  Int_t nId = ((mask >> 8) & 0xFF);
1955  Int_t dId = ((mask >> 16) & 0xFF);
1956  if (nId == 0 && dId == 0) {
1957  ::Warning("Ratio", "Both numerator and denominator IDs are 0");
1958  break;
1959  }
1960  Bool_t rel = ret->IsRelative(rId);
1961  Bool_t crel = true; // rel
1962  Double_t pdEyl = 0;
1963  Double_t pdEyh = 0;
1964  Double_t pnEyl = 0;
1965  Double_t pnEyh = 0;
1966  if (dId != 0) {
1967  HolderP2P* pD = den->FindP2P(dId);
1968  Double_t lfc = (crel ? (TMath::Abs(dY) > 1e-9 ? 1/dY : 0) : 1);
1969  pdEyl = pD->GetYDown(di1, di2, daX) * lfc;
1970  pdEyh = pD->GetYUp (di1, di2, daX) * lfc;
1971  }
1972  if (nId != 0) {
1973  HolderP2P* pN = num->FindP2P(nId);
1974  Double_t lfc = (crel ? (TMath::Abs(nY) > 1e-9 ? 1/nY : 0) : 1);
1975  pnEyl = pN->GetYDown(i) * lfc;
1976  pnEyh = pN->GetYUp (i) * lfc;
1977  }
1978  // Printf("Doing error %d from n=%d (%f +%f -%f) d=%d (%f +%f -%f)",
1979  // rId, nId, nY, pnEyl, pnEyh, dId, dY, pdEyh, pdEyl);
1980  Double_t peyl = 0;
1981  Double_t peyh = 0;
1982  Double_t lfc = (!crel ? (TMath::Abs(dY) > 1e-9 ? 1/dY : 0) : 1);
1983  if (nId == 0) {
1984  // No numerator value - set directly
1985  peyl = pdEyl * lfc;
1986  peyh = pdEyh * lfc;
1987  } // else if
1988  else if (dId == 0) {
1989  // No denomenator value - set directly
1990  peyl = pnEyl * lfc;
1991  peyh = pnEyh * lfc;
1992  } // else if
1993  else {
1994  // Cancel errors
1995  if (denrel) {
1996  peyl = pdEyl*lfc;
1997  peyh = pdEyh*lfc;
1998  // Printf("den err (%f%%,%f%%)", peyl*100,peyh*100);
1999  }
2000  else if (cmin) {
2001  peyl = TMath::Min(nEyl, dEyl) * lfc;
2002  peyh = TMath::Min(nEyh, dEyh) * lfc;
2003  }
2004  else if (cmax) {
2005  peyl = TMath::Min(nEyl, dEyl) * lfc;
2006  peyh = TMath::Min(nEyh, dEyh) * lfc;
2007  }
2008  else {
2009  peyl = TMath::Sqrt(TMath::Abs(pnEyl*pnEyl-pdEyl*pdEyl)) * lfc;
2010  peyh = TMath::Sqrt(TMath::Abs(pnEyh*pnEyh-pdEyh*pdEyh)) * lfc;
2011  // ::Info("Ratio", "sqrt(|%8f^2 - %8f^2|)=%8f", pnEyl,pdEyl,peyl);
2012  }
2013  } // Else
2014  // Printf("%f, den err (%f%%,%f%%) - %d (%f%%,%f%%)",
2015  // rY, pdEyl*100, pdEyh*100, rId, peyl*100,peyh*100);
2016  lfc = (rel ? 1 : rY);
2017  ret->SetSysError(rId, cnt, 0, 0, lfc*peyl, lfc*peyh);
2018  }
2019 
2020  // Finally, increment counter
2021  cnt++;
2022  }
2023 
2024  // Reset the used case bits
2025  num->ClearUsed();
2026  den->ClearUsed();
2027 
2028  if (cnt <= 0) {
2029  ::Warning("", "No common points found");
2030  delete ret;
2031  ret = 0;
2032  return 0;
2033  }
2034  return ret;
2035  }
2040  void ClearUsed() const
2041  {
2042  // Reset the used case bits
2043  TObject* oe = 0;
2044  TIter nextC(&fCommon);
2045  while ((oe = nextC())) {
2046  oe->ResetBit(kUsedBit);
2047  oe->ResetBit(kOnlyWeightBit);
2048  }
2049  TIter nextP(&fPoint2Point);
2050  while ((oe = nextP())) {
2051  oe->ResetBit(kUsedBit);
2052  oe->ResetBit(kOnlyWeightBit);
2053  }
2054  }
2066  const GraphSysErr* g2)
2067  {
2068  Double_t z;
2069  return KolomogorovTest(g1, g2, z);
2070  }
2083  const GraphSysErr* g2,
2084  Double_t& z)
2085  {
2086  TArrayD a1y;
2087  TArrayD a2y;
2088  TArrayD a1e2;
2089  TArrayD a2e2;
2090  Int_t cnt = CacheGraphs(g1, g2, a1y, a2y, a1e2, a2e2);
2091  if (cnt <= 0) return -1;
2092 
2093  Double_t s1 = a1y.GetSum(); // Summed content
2094  Double_t s2 = a2y.GetSum(); // Summed content
2095  Double_t sw1 = a1e2.GetSum(); // Summed weights
2096  Double_t sw2 = a2e2.GetSum(); // Summed weights
2097  if (a1e2.GetSum() <= 0 && a2e2.GetSum()) {
2098  ::Warning("ChisquareTest", "No errors");
2099  return -1;
2100  }
2101 
2102  Double_t sr1 = 0;
2103  Double_t sr2 = 0;
2104  Double_t dfMax = 0;
2105  for (Int_t i = 0; i < cnt; i++) {
2106  // KolmogorowvTest
2107  sr1 += a1y[i] / s1;
2108  sr2 += a2y[i] / s2;
2109  dfMax = TMath::Max(dfMax, TMath::Abs(sr1-sr2));
2110  }
2111  // K-S test
2112  Double_t se1 = s1 * s1 / sw1;
2113  Double_t se2 = s2 * s2 / sw2;
2114  z = dfMax * TMath::Sqrt(se1 * se2 / (se1 + se2));
2115  return TMath::KolmogorovProb(z);
2116  }
2131  static Int_t CacheGraphs(const GraphSysErr* g1,
2132  const GraphSysErr* g2,
2133  TArrayD& a1y,
2134  TArrayD& a2y,
2135  TArrayD& a1e2,
2136  TArrayD& a2e2)
2137  {
2138  Int_t cnt = 0;
2139  Int_t n = g1->GetN();
2140  a1y.Set(n);
2141  a2y.Set(n);
2142  a1e2.Set(n);
2143  a2e2.Set(n);
2144  for (Int_t i = 0; i < g1->GetN(); i++) {
2145  Double_t x = 0;
2146  Double_t y1 = 0;
2147  Double_t e1yl = 0;
2148  Double_t e1yh = 0;
2149  Double_t e2yl = 0;
2150  Double_t e2yh = 0;
2151  Double_t y2 = 0;
2152  if (!NextPoint(i, g1, g2, x, y1, e1yl, e1yh, y2, e2yl, e2yh)) {
2153  ::Warning("ChisquareTest", "Next point - %d %f not found", i, x);
2154  continue;
2155  }
2156 
2157  Double_t e1y2 = TMath::Power(TMath::Max(e1yl,e1yh), 2);
2158  Double_t e2y2 = TMath::Power(TMath::Max(e2yl,e2yh), 2);
2159 
2160  // Cache the numbers
2161  a1y[cnt] = y1; // Store value
2162  a2y[cnt] = y2; // Store value
2163  a1e2[cnt] = e1y2; // Store square error
2164  a2e2[cnt] = e2y2; // Store square error
2165  cnt++;
2166  }
2167  if (cnt <= 0)
2168  ::Warning("CacheGraphs", "No common points found");
2169 
2170  return cnt;
2171  }
2185  const GraphSysErr* g2,
2187  {
2188  Int_t ndf = 0;
2189  return ChisquareTest(g1, g2, ndf, type);
2190  }
2205  const GraphSysErr* g2,
2206  Int_t& ndf,
2207  EChi2Type type)
2208  {
2209  Double_t chi2 = 0;
2210  ndf = -1;
2211  TArrayD a1y;
2212  TArrayD a2y;
2213  TArrayD a1e2;
2214  TArrayD a2e2;
2215  Int_t cnt = CacheGraphs(g1, g2, a1y, a2y, a1e2, a2e2);
2216  if (cnt <= 0) return -1;
2217 
2218  Double_t s1 = a1y.GetSum(); // Summed content
2219  Double_t s2 = a2y.GetSum(); // Summed content
2220  if (type == kModelModel && (a1e2.GetSum() <= 0 && a2e2.GetSum())) {
2221  ::Warning("ChisquareTest", "No errors");
2222  return -1;
2223  }
2224 
2225  ndf = cnt;
2226 
2227  // Now loop over cached points
2228  Double_t s = s1 + s2;
2229  for (Int_t i = 0; i < ndf; i++) {
2230  switch (type) {
2231  case kExperimentExperiment: {
2232  Double_t sum = a1y[i]+a2y[i];
2233  Double_t delta = s2 * a1y[i] - s1 * a2y[i];
2234  chi2 += delta * delta / sum;
2235  }
2236  break;
2237  case kExperimentModel: {
2238  Double_t v1 = s2 * a1y[i] - s1 * a2e2[i];
2239  Double_t v2 = v1 * v1 + 4 * s2 * s2 * a1y[i] * a2e2[i];
2240  Double_t pp = (v1 + v2) / (2 * s2 * s2);
2241  Double_t p1 = pp * s1;
2242  Double_t p2 = pp * s2;
2243  Double_t d1 = a1y[i] - p1;
2244  Double_t d2 = a1y[i] - p2;
2245  chi2 += d1 * d1 / p1;
2246  if (a2e2[i] > 0) chi2 += d2 * d2 / a2e2[i];
2247  }
2248  break;
2249  case kModelModel: {
2250  Double_t sigma = s1 * s1 * a2e2[i] + s2 * s2 * a1e2[i];
2251  Double_t delta = s2 * a1y[i] - s1 * a2y[i];
2252  chi2 += delta * delta / sigma;
2253  }
2254  break;
2255  case kExperimentTruth: {
2256  Double_t delta = a2y[i] - a1y[i];
2257  chi2 += delta * delta / a1e2[i];
2258  }
2259  break;
2260  default:
2261  ::Warning("ChisquareTest", "Should not happen");
2262  return -1;
2263  }
2264 
2265  }
2266  if (type == kExperimentExperiment) chi2 /= s1 * s2;
2267  else if (type == kExperimentTruth) chi2 /= ndf;
2268 
2269  return chi2;
2270  }
2271 
2278  {
2279  if (i < 0 || i >= GetN()) return;
2280  fData->RemovePoint(i);
2281  TIter next(&fPoint2Point);
2282  HolderP2P* p2p = 0;
2283  while ((p2p = static_cast<HolderP2P*>(next()))) {
2284  p2p->fGraph->RemovePoint(i);
2285  }
2286  }
2294  void SwapPoints(Int_t i, Int_t j, Bool_t reflect=false)
2295  {
2296  if (i == j && !reflect) return;
2297  SwapPoints(fData, i, j, reflect);
2298  TIter next(&fPoint2Point);
2299  HolderP2P* p2p = 0;
2300  while ((p2p = static_cast<HolderP2P*>(next()))) {
2301  SwapPoints(p2p->fGraph, i, j, reflect);
2302  }
2303  }
2312  {
2313  GraphSysErr* cpy = new GraphSysErr(*this);
2314  for (Int_t i = 0; i < GetN()/2; i++) {
2315  cpy->SwapPoints(i, GetN()-i-1, true);
2316  }
2317  return cpy;
2318  }
2327  {
2328  GraphSysErr* cpy = new GraphSysErr(*other);
2329 
2330  // Declare "blob" systematic point-to-point error
2331  Int_t id = DeclarePoint2Point("Syst.unc..", false, kBox);
2332 
2333  // Loop over common errors
2334  HolderCommon* cmn = 0;
2335  TIter cNext(&fCommon);
2336  while ((cmn = static_cast<HolderCommon*>(cNext()))) {
2337  Int_t oId = cpy->FindId(cmn->GetTitle());
2338  if (oId <= 0) {
2339  // If the common error ins't found in copy - it should,
2340  // because we did a straight copy - then give up.
2341  Error("Symmetrice", "Common error %s not found in %s",
2342  cmn->GetTitle(), cpy->GetName());
2343  cpy->Print("all");
2344  delete cpy;
2345  cpy = 0;
2346  break;
2347  // }
2348  // }
2349  } // oId <= 0
2350 
2351  // Zero this error here - we should add it later
2352  HolderCommon* oCmn = cpy->FindCommon(oId);
2353  if (cmn->fEyl <= 0 && cmn->fEyh <= 0) {
2354  // Copy the error over
2355  cmn->fEyl = oCmn->fEyl;
2356  cmn->fEyh = oCmn->fEyh;
2357  }
2358  oCmn->Set(0,0);
2359  } // while cmn
2360  if (!cpy) return false;
2361 
2362  // Loop over all points
2363  Int_t cnt = 0;
2364  Int_t n = other->GetN();
2365  TArrayI used(n); // book-keeping
2366  Int_t nonZero = 0; // how many non-zero sys.uncerr.
2367  used.Reset(-1);
2368  for (Int_t i = 0; i < n; i++) {
2369  if (used[i] >= 0 && used[i] != i) {
2370  Int_t j = used[i];
2371  // Info("Symmetrice", "Copy %d to %d instead of %d", j, cnt, i);
2372  SetPoint(cnt, -GetX(j), GetY(j));
2375  SetSysError(id, cnt, 0, GetSysErrorY(id, j));
2376  // Info("Symmetrice", "%f -> %f +/- %f",
2377  // -GetX(j), GetY(j), GetSysErrorY(id,j));
2378  cnt++;
2379  continue;
2380  }
2381  used[i] = i;
2382  Double_t eyl1, eyh1;
2383  Double_t exl1 = cpy->GetErrorXLeft(i);
2384  Double_t exh1 = cpy->GetErrorXRight(i);
2385  Double_t seyl1 = cpy->GetStatErrorDown(i);
2386  Double_t seyh1 = cpy->GetStatErrorUp(i);
2387  Double_t y1 = cpy->GetYandError(i,true,false,true,false,eyl1,eyh1);
2388  Double_t x1 = cpy->GetX(i);
2389  Double_t y2 = 0;
2390  Double_t eyl2 = 0;
2391  Double_t eyh2 = 0;
2392  Double_t seyl2 = 0;
2393  Double_t seyh2 = 0;
2394  Int_t i1, i2;
2395  Int_t j = cpy->FindPoint(-x1, i1, i2); // Find mirror
2396  if (j < -1) {
2397  // IF there's no mirror point, the just set to current
2398  // Info("Symmetrice", "No mirror of %f using %d for %d", x1, i, cnt);
2399  SetPoint(cnt, x1, y1);
2400  SetPointError(cnt, exl1, exh1);
2401  SetStatError(cnt, seyl1, seyh1);
2402  SetSysError(id, cnt, 0, 0, eyl1, eyh1);
2403  cnt++;
2404  continue;
2405  }
2406  if (j >= 0) {
2407  // If we have exact mirror, then use that point and mark as
2408  // used
2409  // Info("Symmetrice", "Got exact mirror of %f(%d) at %d", x1, i, j);
2410  y2 = cpy->GetYandError(j,true,false,true,false,eyl2,eyh2);
2411  seyl2 = cpy->GetStatErrorDown(j);
2412  seyh2 = cpy->GetStatErrorUp(j);
2413  // We should copy values when we get to j
2414  used[j] = cnt;
2415  }
2416  else {
2417  // Otherwise we use interpolation between two points
2418  // Info("Symmetrice", "Interpolate mirror of %f(%d) between %d-%d",
2419  // x1, i, i1, i2);
2420  cpy->FindYandError(-x1,true,true,true,false,y2,eyl2,eyh2,seyl2,seyh2);
2421  }
2422  // Now calculate weighted mean and variance
2423  Double_t newY = (y1+y2) / 2;
2424  Double_t newV = 0;
2425  if ((eyh1+eyl1) > 1e-9 && (eyh2+eyl2) > 1e-9) {
2426  Double_t s1 = 2 * eyl1 * eyh1 / (eyh1 + eyl1);
2427  Double_t sp1 = (eyh1 - eyl1) / (eyh1 + eyl1);
2428  Double_t w1 = .5 * TMath::Power(s1+y1*sp1, 3) / s1;
2429  Double_t s2 = 2 * eyl2 * eyh2 / (eyh2 + eyl2);
2430  Double_t sp2 = (eyh2 - eyl2) / (eyh2 + eyl2);
2431  Double_t w2 = .5 * TMath::Power(s2+y2*sp2, 3) / s2;
2432  // Printf("s1=%f sp1=%f w1=%f s2=%f sp2=%f w2=%f",
2434  Double_t sumW = (w1+w2);
2435  newY = (y1*w1+y2*w2) / sumW;
2436  newV = TMath::Sqrt((w1*w1*s1*s1+w2*w2*s2*s2) / sumW / sumW);
2437  nonZero++;
2438  }
2439  Double_t seyl = TMath::Sqrt(seyl1*seyl1+seyl2*seyl2);
2440  Double_t seyh = TMath::Sqrt(seyh1*seyh1+seyh2*seyh2);
2441 
2442  // Info("Symmmetrice", "%f (%f +/- %f) + (%f +/- %f) -> %f +/- %f",
2443  // x1, y1, w1, y2, w2, newY, newV);
2444  SetPoint(cnt, x1, newY);
2445  SetPointError(cnt, exl1, exh1);
2446  SetStatError(cnt, seyl, seyh);
2447  SetSysError(id, cnt, 0, newV);
2448  cnt++;
2449  }
2450  if (nonZero <= 0) RemoveSysError(id);
2451  // Info("", "After symmetrization:");
2452  // Print("all");
2453  return true;
2454  }
2455  /* @} */
2460  void SavePrimitive(std::ostream& out, Option_t* option="")
2461  {
2462  TString opt(option);
2463  opt.ToLower();
2464  Bool_t load = opt.Contains("load"); opt.ReplaceAll("load", "");
2465  TString funcName;
2466  TPRegexp regex("func=([a-zA-z][a-zA-Z0-9_]*)");
2467  TObjArray* toks = regex.MatchS(opt);
2468  if (toks) {
2469  if (toks->GetEntriesFast() > 1)
2470  funcName = toks->At(1)->GetName();
2471  if (toks->GetEntriesFast() > 0)
2472  opt.ReplaceAll(toks->At(0)->GetName(), "");
2473  delete toks;
2474  }
2475 
2476 
2477  // Save initialization
2478  if (!funcName.IsNull())
2479  out << "TObject* " << funcName << "(Option_t* o=\"\")\n";
2480  out << "{\n";
2481  if (load)
2482  out << " // Load class\n"
2483  << " if (!gROOT->GetClass(\"GraphSysErr\"))\n"
2484  << " gROOT->LoadMacro(\"GraphSysErr.C+\");\n";
2485  out << " GraphSysErr* g = new GraphSysErr(\""
2486  << GetName() << "\",\"" << GetTitle() << "\","
2487  << GetN() << ");\n";
2488  // Save attributes
2489  out << " // Point options\n"
2490  << " g->SetMarkerStyle(" << GetMarkerStyle() << ");\n"
2491  << " g->SetMarkerColor(" << GetMarkerColor() << ");\n"
2492  << " g->SetMarkerSize(" << GetMarkerSize() << ");\n"
2493  << " g->SetLineStyle(" << GetLineStyle() << ");\n"
2494  << " g->SetLineColor(" << GetLineColor() << ");\n"
2495  << " g->SetLineWidth(" << GetLineWidth() << ");\n"
2496  << " g->SetFillStyle(" << GetFillStyle() << ");\n"
2497  << " g->SetFillColor(" << GetFillColor() << ");\n"
2498  << " g->SetXTitle(\"" << fXTitle << "\");\n"
2499  << " g->SetYTitle(\"" << fYTitle << "\");\n"
2500  << " g->SetDataOption(" << fDataOption << ");\n"
2501  << " // Sum options\n"
2502  << " g->SetSumOption(" << fSumOption << ");\n"
2503  << " g->SetSumTitle(\"" << fSumTitle << "\");\n"
2504  << " g->SetSumLineStyle(" << fSumLine.GetLineStyle() << ");\n"
2505  << " g->SetSumLineColor(" << fSumLine.GetLineColor() << ");\n"
2506  << " g->SetSumLineWidth(" << fSumLine.GetLineWidth() << ");\n"
2507  << " g->SetSumFillStyle(" << fSumFill.GetFillStyle() << ");\n"
2508  << " g->SetSumFillColor(" << fSumFill.GetFillColor() << ");\n"
2509  << " g->SetCommonSumOption(" << fCommonSumOption << ");\n"
2510  << " g->SetCommonSumTitle(\"" << fCommonSumTitle << "\");\n"
2511  << " g->SetCommonSumLineStyle(" <<fCommonSumLine.GetLineStyle()<<");\n"
2512  << " g->SetCommonSumLineColor(" <<fCommonSumLine.GetLineColor()<<");\n"
2513  << " g->SetCommonSumLineWidth(" <<fCommonSumLine.GetLineWidth()<<");\n"
2514  << " g->SetCommonSumFillStyle(" <<fCommonSumFill.GetFillStyle()<<");\n"
2515  << " g->SetCommonSumFillColor(" <<fCommonSumFill.GetFillColor()<<");\n"
2516  << " // Stat options\n"
2517  << " g->SetStatRelative(" << fStatRelative << ");\n";
2518  TIter nextC(&fCommon);
2519  HolderCommon* cmn = 0;
2520  while ((cmn = static_cast<HolderCommon*>(nextC())))
2521  cmn->SavePrimitive(out, "d");
2522  TIter nextP(&fPoint2Point);
2523  HolderP2P* p2p = 0;
2524  while ((p2p = static_cast<HolderP2P*>(nextP())))
2525  p2p->SavePrimitive(out, "d");
2526  Int_t n = GetN();
2527  out << " // " << n << " points\n";
2528  Bool_t statRel = IsStatRelative();
2529  for (Int_t i = 0; i < n; i++) {
2530  Double_t y = GetY(i);
2531  out << " g->SetPoint(" << i << ',' << GetX(i) << ',' << y << ");\n"
2532  << " g->SetPointError(" << i << ',' << GetErrorXLeft(i) << ','
2533  << GetErrorXRight(i) << ");\n"
2534  << " g->SetStatError(" << i << ','
2535  << (statRel ? 1/y : 1)*GetStatErrorDown(i) << ','
2536  << (statRel ? 1/y : 1)*GetStatErrorUp(i) << ");\n";
2537  nextP.Reset();
2538  while ((p2p = static_cast<HolderP2P*>(nextP()))) {
2539  Int_t id = p2p->GetUniqueID();
2540  Bool_t rel = p2p->IsRelative();
2541  out << " g->SetSysError(" << id << ',' << i << ','
2542  << GetSysErrorXLeft(id, i) << ','
2543  << GetSysErrorXRight(id, i) << ','
2544  << (rel ? 1/y : 1) * GetSysErrorYDown(id, i) << ','
2545  << (rel ? 1/y : 1) * GetSysErrorYUp(id, i) << ");\n";
2546  } // while(p2p)
2547  } // for (i)
2548  out << " if (o && o[0] != '\\0') {\n"
2549  << " g->Draw(o);\n";
2550  if (fDrawn && fDrawn->GetHistogram())
2551  out << " if (g->GetMulti() && g->GetMulti()->GetHistogram()) {\n"
2552  << " g->GetMulti()->GetHistogram()->SetMinimum("
2553  << fDrawn->GetHistogram()->GetMinimum() << ");\n"
2554  << " g->GetMulti()->GetHistogram()->SetMaximum("
2555  << fDrawn->GetHistogram()->GetMaximum() << ");"
2556  << " }\n";
2557 
2558  out << " }\n";
2559  if (!funcName.IsNull()) out << " return g;\n";
2560  out << "};\n";
2561  }
2567  void Save(const char* fileName)
2568  {
2569  std::ofstream out(fileName);
2570  TString funcName(fileName);
2571  funcName.ReplaceAll(".C","");
2572  out << "// \n"
2573  << "// Generated by GraphSysErr.C\n"
2574  << "// \n"
2575  // << "class GraphSysErr;\n\n";
2576  << "\n";
2577  SavePrimitive(out, Form("load func=%s", funcName.Data()));
2578  out.close();
2579  }
2612  void Export(std::ostream& out=std::cout,
2613  Option_t* option="",
2614  Int_t nsign=2)
2615  {
2616  // Info("Export", "Using %d significant digits on errors", nsign);
2617  TString opt(option);
2618  opt.ToLower();
2619  Bool_t header = opt.Contains("h");
2620  Bool_t sysNames = opt.Contains("s");
2621  Bool_t comment = opt.Contains("c");
2622  Bool_t title = opt.Contains("t");
2623 
2624  ExportHeader(out, header, comment, nsign);
2625 
2626  // --- Export qualifiers -----------------------------------------
2627  Bool_t hasTitle = false;
2628  if (fQualifiers) {
2629  TIter nextQ(fQualifiers);
2630  TObject* q = 0;
2631  while ((q = nextQ())) {
2632  TString k(q->GetName());
2633  if (k.EqualTo("RE") || k.EqualTo("title", TString::kIgnoreCase))
2634  hasTitle = true;
2635  out << FormatKey("qual") << q->GetName() << " : "
2636  << q->GetTitle() << std::endl;
2637  }
2638  }
2639  if (!hasTitle && title)
2640  out << FormatKey("qual") << "RE : " << GetTitle() << std::endl;
2641 
2642  // --- Export X/Y titles ----------------------------------------
2643  const char* fill = "<please fill in>";
2644  TString xTit = fXTitle; EscapeLtx(xTit, "fill");
2645  TString yTit = fYTitle; EscapeLtx(yTit, "fill");
2646 
2647  out << FormatKey("xheader") << xTit << "\n"
2648  << FormatKey("yheader") << yTit << std::endl;
2649 
2650  // --- Export common errors --------------------------------------
2651  TIter nextC(&fCommon);
2652  HolderCommon* holderCommon = 0;
2653  while ((holderCommon = static_cast<HolderCommon*>(nextC()))) {
2654  Bool_t rel = holderCommon->IsRelative();
2655  Double_t up = holderCommon->GetYUp(rel ? 100 : 1);
2656  Double_t down = holderCommon->GetYDown(rel ? 100 : 1);
2657  out << FormatKey("dserror");
2658  ExportError(out, down, up, true, rel, nsign);
2659  out << ":" << holderCommon->GetTitle() << std::endl;
2660  }
2661 
2662  // --- Export data points ----------------------------------------
2663  out << FormatKey("data") << " x : y" << std::endl;
2664  Int_t n = GetN();
2665  for (Int_t i = 0; i < n; i++) {
2666  ExportPoint(out, i, true, sysNames, nsign);
2667  out << std::endl;
2668  }
2669  out << "*dataend:\n"
2670  << "# End of dataset\n" << std::endl;
2671  } //*MENU*
2678  static void EscapeLtx(TString& val, const TString& fill="")
2679  {
2680  if (val.IsNull()) {val = fill; return; }
2681  if (!val.Contains("#") && !val.Contains("\\")) return;
2682 
2683  if (val[0] !='$') val.Prepend("$");
2684  if (val[val.Length()-1]!='$') val.Append("$");
2685  val.ReplaceAll("#", "\\");
2686  }
2701  static void Export(const TSeqCollection* col,
2702  std::ostream& out,
2703  Option_t* option="H",
2704  Int_t nsign=0)
2705  {
2706  if (col->GetEntries() < 1) return;
2707 
2708  // --- Deduce options --------------------------------------------
2709  TString opt(option);
2710  opt.ToLower();
2711  Bool_t alsoTop = opt.Contains("h");
2712  Bool_t alsoCmt = opt.Contains("c");
2713  Bool_t alsoNme = opt.Contains("s");
2714  Bool_t alsoTit = opt.Contains("t");
2715 
2716  // --- some variables to use -------------------------------------
2717  const Double_t tol = 1e-10;
2718  GraphSysErr* first = 0;
2719  GraphSysErr* gse = 0;
2720  TList toExport;
2721  TList commons;
2722  TList quals;
2723  TString data("x ");
2724  Int_t nPoints = -1;
2725  Int_t idx = 0;
2726 
2727  // --- Check the passed graphs for compatiblity ------------------
2728  TIter nextCheck(col);
2729  TObject* o = 0;
2730  while ((o = nextCheck())) {
2731  if (!o->IsA()->InheritsFrom(GraphSysErr::Class())) continue;
2732  gse = static_cast<GraphSysErr*>(o);
2733  if (!first) {
2734  first = gse;
2735  nPoints = first->GetN();
2736  }
2737  else {
2738  // --- Check number of points --------------------------------
2739  if (gse->GetN() != nPoints) {
2740  Int_t nTmp = TMath::Min(gse->GetN(), nPoints);
2741  ::Warning("Export", "Incompatible number of points %d in %s"
2742  "using only %d of them",
2743  gse->GetN(), gse->GetName(), nTmp);
2744  nPoints = nTmp;
2745  }
2746  // --- check X values ----------------------------------------
2747  Bool_t ok = true;
2748  for (Int_t i = 0; i < nPoints; i++) {
2749  Double_t x1 = first->GetX(i);
2750  Double_t exl1 = first->fData->GetErrorXlow(i);
2751  Double_t exh1 = first->fData->GetErrorXhigh(i);
2752  Double_t xT = gse->GetX(i);
2753  Double_t exlT = gse->fData->GetErrorXlow(i);
2754  Double_t exhT = gse->fData->GetErrorXhigh(i);
2755  if (TMath::Abs(x1-xT) > tol ||
2756  (exl1 > tol && TMath::Abs(exl1-exlT) > tol) ||
2757  (exh1 > tol && TMath::Abs(exh1-exhT) > tol)) {
2758  ::Warning("Export", "X--coordinate of %s @ point %d: %f (+%f,-%f) "
2759  "incompatible [%f (+%f,-%f)]",
2760  gse->GetTitle(), i, xT, exhT, exlT, x1, exh1, exl1);
2761  ok = false;
2762  break;
2763  }
2764  } // for i
2765  if (!ok) continue;
2766  } // !first
2767 
2768  // --- Get all possible qualifiers -----------------------------
2769  toExport.Add(gse);
2770  TIter nextQ(gse->fQualifiers);
2771  TObject* q = 0;
2772  while ((q = nextQ())) {
2773  // ::Info("", "qualifier %s=%s", q->GetName(), q->GetTitle());
2774  StoreQual(quals, idx, q);
2775  }
2776 
2777  // --- Get all common systematics ------------------------------
2778  TIter nextCmn(&(gse->fCommon));
2779  TObject* oCmn = 0;
2780  while ((oCmn = nextCmn())) {
2781  if (commons.FindObject(oCmn->GetTitle())) continue;
2782  TObjString* cmn = new TObjString(oCmn->GetTitle());
2783  if (gse == first) cmn->SetUniqueID(gse->FindId(oCmn->GetTitle()));
2784  commons.Add(cmn);
2785  }
2786 
2787  idx++;
2788  data += ": y ";
2789  }
2790  // --- Now deduce the common common errors -----------------------
2791  TIter nextCmn(&(commons));
2792  TObject* oCmn = 0;
2793  while ((oCmn = nextCmn())) {
2794  TString oNme(oCmn->GetName());
2795  Bool_t found = true;
2796  Double_t ecl1 = -1;
2797  Double_t ech1 = -1;
2798 
2799  // --- Loop over data ------------------------------------------
2800  TIter nextG(&toExport);
2801  while ((gse = static_cast<GraphSysErr*>(nextG()))) {
2803  Int_t gseId = first->FindId(oNme);
2804 
2805  if (gseId > 0) {
2806  // If we do, check compatibility
2807  if (ecl1 < 0) {
2808  // First time we get this error
2809  ecl1 = gse->GetCommonErrorYDown(gseId);
2810  ech1 = gse->GetCommonErrorYUp(gseId);
2811  continue;
2812  }
2813  else {
2814  // Now check values within tolerance
2815  Double_t eclT = gse->GetCommonErrorYDown(gseId);
2816  Double_t echT = gse->GetCommonErrorYUp(gseId);
2817  if ((ecl1 > tol && TMath::Abs(ecl1-eclT) > tol) ||
2818  (ech1 > tol && TMath::Abs(ech1-echT) > tol)) {
2819  found = false;
2820  } // incompatible
2821  } // ecl1 >= 0
2822  } // gseId > 0
2823  else
2824  // This graph does not have this common, flag it
2825  found = false;
2826  } // while(gse)
2827 #if INCOMPAT_CMN_AS_QUAL
2828 #else
2829  found = true;
2830 #endif
2831  oCmn->SetBit(BIT(14), found);
2832 
2833  // If we found the error in all graphs and they are all
2834  // compatible, then all is good,
2835  if (found) continue;
2836 
2837  // If not, we represent the common systematic as a qualifier
2838  idx = 0;
2839  nextG.Reset();
2840  while ((gse = static_cast<GraphSysErr*>(nextG()))) {
2841  Int_t gseId = first->FindId(oNme);
2842  TString val;
2843  if (gseId > 0) {
2844  // This graph has the value, store as qual
2845  Bool_t rel = gse->IsRelative(gseId);
2846  Double_t ecl = gse->GetCommonErrorYDown(gseId);
2847  Double_t ech = gse->GetCommonErrorYUp(gseId);
2848  std::stringstream s;
2849  ExportError(s, ecl, ech, false, rel, nsign);
2850  val = s.str().c_str();
2851  }
2852  else {
2853  // This does not have the error, store empty string
2854  val = "";
2855  }
2856  // Now store this in our qualifier table
2857  StoreQual(quals, idx, oNme, val);
2858  idx++;
2859  } // while(gse)
2860  } // while common
2861 
2862  // --- Sanity checks ---------------------------------------------
2863  if (nPoints <= 0) {
2864  ::Error("Export", "No points to write");
2865  return;
2866  }
2867  if (toExport.GetEntries() <= 0) {
2868  ::Error("Export", "No graphs to export");
2869  return;
2870  }
2871  if (!first) {
2872  ::Error("Export", "Didn't get the first graph");
2873  return;
2874  }
2875 
2876  // --- Export header --------------------------------------------
2877  first->ExportHeader(out, alsoTop, alsoCmt, nsign);
2878 
2879  // --- Export qualifiers -----------------------------------------
2880  // quals.ls();
2881  Bool_t hasTitle = false;
2882  TIter nextQ(&quals);
2883  TList* ql = 0;
2884  while ((ql = static_cast<TList*>(nextQ()))) {
2885  TString k(ql->GetName());
2886  if (k.EqualTo("RE") || k.EqualTo("title", TString::kIgnoreCase))
2887  hasTitle = true;
2888  out << FormatKey("qual") << ql->GetName();
2889  for (Int_t i = 0; i < toExport.GetEntries(); i++) {
2890  TObject* qv = ql->At(i);
2891  TString v = "";
2892  if (qv) v = qv->GetName();
2893  out << " : " << v;
2894  }
2895  out << std::endl;
2896  }
2897  if (!hasTitle && alsoTit) {
2898  out << FormatKey("qual") << "RE ";
2899  for (Int_t i = 0; i < toExport.GetEntries(); i++)
2900  out << ": " << toExport.At(i)->GetTitle();
2901  out << std::endl;
2902  }
2903 
2904 
2905  // --- Export axis titles ----------------------------------------
2906  const char* fill = "<please fill in>";
2907  const char* fields[] = { "xheader", "yheader", 0 };
2908  const char** pfld = fields;
2909  while (*pfld) {
2910  out << FormatKey(*pfld);
2911  TIter nextSpec(&toExport);
2912  Bool_t one = true;
2913  while ((gse = static_cast<GraphSysErr*>(nextSpec()))) {
2914  TString val;
2915  if ((*pfld)[0] == 'x') val = gse->fXTitle;
2916  else if ((*pfld)[0] == 'y') val = gse->fYTitle;
2917  EscapeLtx(val, fill);
2918  out << (one ? "" : ":") << val;
2919  if ((*pfld)[0] == 'x') break;
2920  one = false;
2921  }
2922  pfld++;
2923  out << std::endl;
2924  }
2925 
2926  // --- Export common systematics ---------------------------------
2927  TIter nextC(&commons);
2928  while ((oCmn = nextC())) {
2929  TString tit(oCmn->GetName());
2930  EscapeLtx(tit);
2931  out << FormatKey("dserror");
2932 #if INCOMPAT_CMN_AS_QUAL
2933  if (!oCmn->TestBit(BIT(14))) {
2934  ::Warning("Export",
2935  "Common systematic error \"%s\" represented by qual",
2936  oCmn->GetName());
2937  out << "- : " << tit << std::endl;
2938  continue;
2939  }
2940  UInt_t id = first->FindId(oCmn->GetName());
2941  Bool_t rel = first->IsRelative(id);
2942  Double_t ecl = first->GetCommonErrorYDown(id);
2943  Double_t ech = first->GetCommonErrorYUp(id);
2944  ExportError(out, ecl, ech, true, rel, nsign);
2945  out << " : "<< tit << std::endl;
2946 #else
2947  out << tit << std::flush;
2948  TIter next(&toExport);
2949  while ((gse = static_cast<GraphSysErr*>(next()))) {
2950  UInt_t id = gse->FindId(oCmn->GetName());
2951  Bool_t rel = gse->IsRelative(id);
2952  Double_t ecl = gse->GetCommonErrorYDown(id);
2953  Double_t ech = gse->GetCommonErrorYUp(id);
2954  out << " : ";
2955  ExportError(out, ecl, ech, true, rel, nsign);
2956  }
2957  out << ":" << tit << std::endl;
2958 #endif
2959  }
2960 
2961  // --- Export points ---------------------------------------------
2962  out << FormatKey("data") << data << std::endl;
2963  for (Int_t i = 0; i < nPoints; i++) {
2964  TIter next(&toExport);
2965  Bool_t one = true;
2966  while ((gse = static_cast<GraphSysErr*>(next()))) {
2967  gse->ExportPoint(out, i, one, alsoNme, nsign);
2968  one = false;
2969  }
2970  out << std::endl;
2971  }
2972  out << "*dataend:\n"
2973  << "# End of dataset\n" << std::endl;
2974  // return out;
2975  }
2993  static TSeqCollection* Import(const TString& fileName)
2994  {
2995  std::ifstream in(fileName.Data());
2996  if (!in) {
2997  ::Warning("Import", "Failed to open \"%s\"", fileName.Data());
2998  return 0;
2999  }
3000  TList* ret = new TList;
3001  ret->SetOwner();
3002  GraphSysErr* g = 0;
3003  GraphSysErr* first = 0;
3004  Int_t id = 0;
3005  do {
3006  Int_t sub = 1;
3007  UShort_t nIdx = 256;
3008  int cur = in.tellg();
3009  do {
3010  in.seekg(cur, in.beg);
3011  if (!(g = Import(in, sub, &nIdx))) break;
3012  if (kVerbose & kImport)
3013  ::Info("Import", "Imported %d of %d", sub, nIdx);
3014  g->SetName(Form("ds_%d_%d", id, sub-1));
3015  if (!first) first = g;
3016  else g->CopyKeys(first);
3017  ret->Add(g);
3018  sub++;
3019  } while(sub < nIdx);
3020  id++;
3021  } while (!in.eof());
3022  return ret;
3023  }
3043  static GraphSysErr* Import(std::istream& in, UShort_t idx=0,
3044  UShort_t* nIdx=0)
3045  {
3046  GraphSysErr* ret = 0;
3047  Int_t n = 0;
3048  Bool_t inSet = false;
3049  Bool_t inData = false;
3050  TString xtit = "";
3051  TString ytit = "";
3052  TString tmp = "";
3053  TString last = "";
3054  TList quals;
3055  TList keys;
3056  Int_t isty = 0;
3057  do {
3058  TString line;
3059  line.ReadLine(in);
3060  tmp = line.Strip(TString::kBoth, ' '); line = tmp;
3061  if (line.IsNull()) continue;
3062  if (line[0] == '#') continue;
3063  if (line[0] == '*') { // We have some sort of header stuff
3064  Int_t colon = line.Index(":");
3065  if (colon == kNPOS) continue;
3066  Int_t len = line.Length()-colon-1;
3067  TString value = line(colon+1,len);
3068  tmp = value.Strip(TString::kBoth, ' '); value = tmp;
3069  last = line(1,colon-1);
3070  if (kVerbose & kImport)
3071  ::Info("Import", "Got a key '%s' -> '%s'", last.Data(), value.Data());
3072 
3073  if (last.BeginsWith("dataset", TString::kIgnoreCase)) {
3074  inSet = true;
3075  inData = false;
3076  isty = 0;
3077  }
3078  else if (last.BeginsWith("dataend", TString::kIgnoreCase)) {
3079  inSet = false;
3080  inData = false;
3081  // Always terminate on the end of a data set
3082  break;
3083  }
3084  // else if (last.BeginsWith("dscomment",TString::kIgnoreCase)) {
3085  // tit = value;
3086  // }
3087  else if (last.BeginsWith("qual", TString::kIgnoreCase)) {
3088  // ::Info("", "Got qual: %s", value.Data());
3089  // if (!qual.IsNull()) {
3090  // ::Info("", "Adding seen qual: %s", qual.Data());
3091  quals.Add(new TObjString(value));
3092  }
3093  else if (last.BeginsWith("xheader", TString::kIgnoreCase))
3094  xtit = value;
3095  else if (last.BeginsWith("yheader", TString::kIgnoreCase))
3096  ytit = value;
3097  else if (last.BeginsWith("dserror", TString::kIgnoreCase)) {
3098  // Common systematic error
3099  if (kVerbose & kImport)
3100  ::Info("Import", "Got common error line: '%s'", line.Data());
3101  TObjArray* tokens = value.Tokenize(":");
3102  Double_t el = 0, eh = 0;
3103  Bool_t rel = false;
3104  if (ImportError(Token(tokens, 0), el, eh, rel)) {
3105  if (rel) { el /= 100.; eh /= 100.; }
3106  TString& tnam = Token(tokens, 1);
3107  TString nam = tnam.Strip(TString::kBoth, ' ');
3108  if (!ret) ret = new GraphSysErr("imported", "");
3109  Int_t id = ret->DefineCommon(nam, rel, el, eh, kRect);
3110  ret->SetSysLineColor(id, (isty % 6) + 2);
3111  ret->SetSysFillColor(id, (isty % 6) + 2);
3112  ret->SetSysFillStyle(id, 3001 + (isty % 10));
3113  isty++;
3114  }
3115  }
3116  else if (last.BeginsWith("data", TString::kIgnoreCase)) {
3117  if (kVerbose & kImport)
3118  ::Info("Import", "Got start of data line: '%s'", line.Data());
3119  // These are the field we can deal with
3120  // Let's get the header field value
3121  TObjArray* tokens = value.Tokenize(":");
3122  if (nIdx) {
3123  *nIdx = tokens->GetEntriesFast();
3124  // ::Info("Import", "Max index is %d", *nIdx);
3125  }
3126  if (tokens->GetEntriesFast() > 2 && idx < 1) {
3127  idx = 1;
3128  ::Warning("Import", "Can only import one data set at a time, "
3129  "selecting the %d column", idx);
3130  }
3131  else if (idx >= 1 && idx >= tokens->GetEntriesFast()) {
3132  ::Warning("Import", "column %d not available for this data set",
3133  idx);
3134  inSet = false;
3135  }
3136  // We do nothing to the tokens here.
3137  tokens->Delete();
3138  inData = true;
3139  }
3140  else {
3141  // some other key.
3142  if (kVerbose & kImport)
3143  ::Info("Import", "Got pair line '%s' (%s -> %s)", line.Data(),
3144  last.Data(), value.Data());
3145  keys.Add(new TNamed(last, value));
3146  }
3147  continue;
3148  }
3149  if (!inData) {
3150  if (kVerbose & kImport)
3151  ::Info("Import", "Got a possible contiuation line '%s'", line.Data());
3152  // Probably a continuation line.
3153  if (last.IsNull())
3154  // No last key or no map
3155  continue;
3156 
3157  if (kVerbose & kImport)
3158  ::Info("Import", "Got some other line: '%s'", line.Data());
3159 
3160  TNamed* ptr = static_cast<TNamed*>(keys.FindObject(last));
3161  if (!ptr)
3162  // Key is not added, do nothing to this line
3163  continue;
3164 
3165  TString tt(ptr->GetTitle());
3166  if (!tt.EndsWith(" ") && !tt.EndsWith("\t") && !tt.EndsWith("\n"))
3167  // If value does not end in a white-space, add it.
3168  tt.Append("\n");
3169 
3170  // Append our line, and set the title
3171  tt.Append(line);
3172  ptr->SetTitle(tt);
3173 
3174  continue;
3175  }
3176  if (!inSet) {
3177  // We're outside a data-set so do nothing
3178  continue;
3179  }
3180  if (idx < 1) idx = 1;
3181 
3182  if (kVerbose & kImport)
3183  ::Info("Import", "Got a data line: '%s'", line.Data());
3184  TObjArray* tokens = line.Tokenize(";");
3185  if (tokens->GetEntriesFast() < idx+1) {
3186  ::Warning("Import", "Too few columns %d<%0d in line: %s",
3187  tokens->GetEntriesFast(), idx+1, line.Data());
3188  tokens->Delete();
3189  continue;
3190  }
3191  TString& xCol = static_cast<TObjString*>(tokens->At(0))->String();
3192  TString& yCol = static_cast<TObjString*>(tokens->At(idx))->String();
3193  if (kVerbose & kImport)
3194  ::Info("Import", "xColumn: '%s', yColumn: '%s'",
3195  xCol.Data(), yCol.Data());
3196 
3197  tmp = xCol.Strip(TString::kBoth, ' '); xCol = tmp;
3198  tmp = yCol.Strip(TString::kBoth, ' '); yCol = tmp;
3199  TString yFull = yCol;
3200 
3201  Int_t chop = yCol.Last('(');
3202  if (chop != kNPOS) {
3203  yCol.Remove(chop, yCol.Length()-chop);
3204  yCol.Remove(TString::kBoth, ' ');
3205  }
3206  if (kVerbose & kImport)
3207  ::Info("Import", "xColumn: '%s', yColumn: '%s'",
3208  xCol.Data(), yCol.Data());
3209 
3210  if (xCol.IsNull()) {
3211  ::Warning("Import", "Empty X column in line: %s", line.Data());
3212  tokens->Delete();
3213  continue;
3214  }
3215  if (yCol.IsNull()) {
3216  ::Warning("Import", "Empty Y column in line: %s", line.Data());
3217  tokens->Delete();
3218  continue;
3219  }
3220  Bool_t rel = false;
3221  Double_t x = 0, exl = 0, exh = 0;
3222  if (!ImportPoint(xCol, x, exl, exh, rel)) {
3223  ::Warning("Import", "Failed to import X value in line: %s",line.Data());
3224  tokens->Delete();
3225  continue;
3226  }
3227  Double_t y = 0, eyl = 0, eyh = 0;
3228  if (!ImportPoint(yCol, y, eyl, eyh, rel)) {
3229  ::Warning("Import", "Failed to import X value in line: %s",line.Data());
3230  tokens->Delete();
3231  continue;
3232  }
3233  if (!ret) ret = new GraphSysErr("imported", "");
3234  if (rel) ret->SetStatRelative(true);
3235  ret->SetPoint(n, x, y);
3236  ret->SetPointError(n, exl, exh);
3237 
3238  if (ret->IsStatRelative()) { eyl /= 100; eyh /= 100; }
3239  ret->SetStatError(n, eyl, eyh);
3240 
3241  Int_t lparen = yFull.Index("(");
3242  Int_t rparen = yFull.Index(")");
3243  if (lparen == rparen) {
3244  n++;
3245  continue;
3246  }
3247  TString rem = yFull(lparen+1, rparen-lparen-1);
3248  if (kVerbose & kImport)
3249  ::Info("Import", "Got systematic errors: '%s'", rem.Data());
3250 
3251  TObjArray* stok = rem.Tokenize("D");
3252  Int_t isys = 0;
3253  for (Int_t i = 0; i < stok->GetEntriesFast(); i++) {
3254  TString t = Token(stok,i);
3255  if (!t.BeginsWith("SYS=",TString::kIgnoreCase)) {
3256  ::Warning("Import", "Ignoring unknown spec %s (@ %d) in %s",
3257  t.Data(), i, rem.Data());
3258  continue;
3259  }
3260  t.Remove(0,4);
3261  if (t[t.Length()-1] == ',') t.Remove(t.Length()-1,1);
3262  tmp = t.Strip(TString::kBoth, ' '); t = tmp;
3263 
3264  Double_t el = 0, eh = 0;
3265  rel = false;
3266  if (!ImportError(t, el, eh, rel)) continue;
3267 
3268  Int_t colon = t.Index(":");
3269  TString nam(Form("sys%d", isys));
3270  if (colon != kNPOS) nam = t(colon+1, t.Length()-colon-1);
3271 
3272  UInt_t id = ret->FindId(nam);
3273  if (id == 0) {
3274  id = ret->DeclarePoint2Point(nam, rel, kRect);
3275  ret->SetSysLineColor(id, (isty % 6) + 2);
3276  ret->SetSysFillColor(id, (isty % 6) + 2);
3277  ret->SetSysFillStyle(id, 3001 + (isty % 10));
3278  isty++;
3279  }
3280  HolderP2P* p = ret->FindP2P(id);
3281  if (p->IsRelative()) { el /= 100.; eh /= 100.; }
3282  ret->SetSysError(id, n, exl, exh, el, eh);
3283  isys++;
3284  }
3285  stok->Delete();
3286 
3287  n++;
3288  } while (!in.eof());
3289  if (ret) {
3290  TString rxtit = ExtractField(xtit, idx-1);
3291  TString rytit = ExtractField(ytit, idx-1);
3292  if (!rxtit.IsNull()) ret->SetXTitle(rxtit);
3293  if (!rytit.IsNull()) ret->SetYTitle(rytit);
3294 
3295  Bool_t hasTitle = false;
3296  TIter nextQ(&quals);
3297  TObject* qual = 0;
3298  while ((qual = nextQ())) {
3299  TString k = ExtractField(qual->GetName(), 0);
3300  TString q = ExtractField(qual->GetName(), idx);
3301  /*
3302  ::Info("LoopQ", "Qualifier string: %s -> %s,%s",
3303  qual->GetName(), k.Data(), q.Data());
3304  */
3305  ret->AddQualifier(k, q);
3306  if (!hasTitle &&
3307  (k.EqualTo("RE") ||
3308  k.EqualTo("title", TString::kIgnoreCase))) {
3309  hasTitle = true;
3310  ret->SetTitle(q);
3311  }
3312  }
3313  TIter nextK(&keys);
3314  TObject* pair = 0;
3315  while ((pair = nextK())) {
3316  ret->SetKey(pair->GetName(),pair->GetTitle());
3317  TString k = pair->GetName();
3318  if (!hasTitle && (k.EqualTo("location"))) {
3319  ret->SetTitle(pair->GetTitle());
3320  hasTitle = true;
3321  }
3322  }
3323  if (!hasTitle) ret->SetTitle("Title");
3324  }
3325  keys.SetOwner();
3326  keys.Delete();
3327  return ret;
3328  }
3329  /* @} */
3330  //__________________________________________________________________
3350  UInt_t DefineCommon(const char* title, Bool_t relative,
3351  Double_t ey, EDrawOption_t option=kFill)
3352  {
3353  return DefineCommon(title, relative, ey, ey, option);
3354  }
3366  UInt_t DefineCommon(const char* title, Bool_t relative,
3367  Double_t eyl, Double_t eyh,
3368  EDrawOption_t option=kRect)
3369  {
3370  UInt_t id = ++fCounter;
3371  TString name(Form("common_%08x", id));
3372 
3373  HolderCommon* h = new HolderCommon(name, title, relative, option, id);
3374  h->Set(eyl, eyh);
3375 
3376  fCommon.AddLast(h);
3377  return id;
3378  }
3393  UInt_t DeclarePoint2Point(const char* title, Bool_t relative,
3394  EDrawOption_t option=kBar)
3395  {
3396  UInt_t id = ++fCounter;
3397  TString name(Form("p2p_%08x", id));
3398 
3399  Holder* h = new HolderP2P(name, title, relative, option, id);
3400 
3401  fPoint2Point.AddLast(h);
3402  return id;
3403  }
3405  {
3406  HolderCommon* h = FindCommon(id);
3407  if (h) {
3408  fCommon.Remove(h);
3409  return;
3410  }
3411  HolderP2P* p = FindP2P(id);
3412  if (p) {
3413  fPoint2Point.Remove(p);
3414  return;
3415  }
3416  }
3424  UInt_t FindId(const char* title) const
3425  {
3426  TString tit(title);
3427 
3428  TIter nextC(&fCommon);
3429  Holder* holder = 0;
3430  while ((holder = static_cast<Holder*>(nextC()))) {
3431  if (tit.EqualTo(holder->GetTitle())) return holder->GetUniqueID();
3432  }
3433 
3434  TIter nextP(&fPoint2Point);
3435  while ((holder = static_cast<Holder*>(nextP()))) {
3436  if (tit.EqualTo(holder->GetTitle())) return holder->GetUniqueID();
3437  }
3438  // fCommon.ls();
3439  // fPoint2Point.ls();
3440  return 0;
3441  }
3442  Int_t GetNSys() const { return fCounter; }
3443  /* @} */
3444 
3445  //__________________________________________________________________
3458  {
3459  if (!fData) return;
3460  fData->SetPoint(i, x, y);
3461  }
3469  {
3470  SetPointError(i, ex, ex);
3471  }
3480  {
3481  if (!fData) return;
3482  fData->SetPointError(i, exl, exh,
3483  fData->GetErrorYlow(i),
3484  fData->GetErrorYhigh(i));
3485  }
3492  void SetStatRelative(Bool_t rel) { fStatRelative = rel; }
3506  {
3507  SetStatError(i, ey, ey);
3508  }
3509  /*
3510  * Set the statistical error on the ith data point
3511  *
3512  * @param i Point number
3513  * @param eyl Lower error on Y
3514  * @param eyh Higher error on Y
3515  */
3517  {
3518  if (!fData) return;
3519  if (fStatRelative) {
3520  eyl *= fData->GetY()[i];
3521  eyh *= fData->GetY()[i];
3522  }
3523  fData->SetPointError(i,
3524  fData->GetErrorXlow(i),
3525  fData->GetErrorXhigh(i),
3526  eyl, eyh);
3527  }
3528  void SetSysError(Int_t id, Double_t eyl, Double_t eyh)
3529  {
3530  HolderCommon* h = FindCommon(id);
3531  if (!h) return;
3532  h->Set(eyl, eyh);
3533  }
3543  {
3544  HolderP2P* h = FindP2P(id);
3545  if (!h) return;
3546  h->Set(i, fData, ex, ey);
3547  }
3558  void SetSysError(Int_t id, Int_t i, Double_t exl, Double_t exh,
3559  Double_t eyl, Double_t eyh)
3560  {
3561  HolderP2P* h = FindP2P(id);
3562  if (!h) return;
3563  h->Set(i, fData, exl, exh, eyl, eyh);
3564  }
3565  /* @} */
3566  //__________________________________________________________________
3576  void SetTitle(const char* name)
3577  {
3578  TNamed::SetTitle(name);
3579  if (fData) fData->SetTitle(name);
3580  }
3592  void SetXTitle(const char* title) { fXTitle = title; } //*MENU*
3598  void SetYTitle(const char* title) { fYTitle = title; } //*MENU*
3599  /* @} */
3600 
3601  //__________________________________________________________________
3609  Int_t GetN() const { return fData ? fData->GetN() : 0; }
3615  Double_t GetX(Int_t p) const { return fData->GetX()[p]; }
3621  Double_t GetErrorXLeft(Int_t p) const { return fData->GetErrorXlow(p); }
3627  Double_t GetErrorXRight(Int_t p) const { return fData->GetErrorXhigh(p); }
3633  Double_t GetY(Int_t point) const { return fData->GetY()[point]; }
3639  Double_t GetStatError(Int_t point) const { return fData->GetErrorY(point); }
3646  {
3647  return fData->GetErrorYhigh(point);
3648  }
3655  {
3656  return fData->GetErrorYlow(point);
3657  }
3659  {
3660  HolderCommon* c = FindCommon(id);
3661  return c != 0;
3662  }
3671  {
3672  HolderP2P* h = FindP2P(id);
3673  if (h) return h->IsRelative();
3674  HolderCommon* c = FindCommon(id);
3675  if (c) return c->IsRelative();
3676  return false;
3677  }
3684  Double_t GetSysErrorX(Int_t id, Int_t point) const
3685  {
3686  HolderP2P* h = FindP2P(id);
3687  if (!h) return 0;
3688  return h->GetX(point);
3689  }
3696  Double_t GetSysErrorY(Int_t id, Int_t point) const
3697  {
3698  HolderP2P* h = FindP2P(id);
3699  if (!h) return 0;
3700  return h->GetY(point);
3701  }
3709  {
3710  HolderP2P* h = FindP2P(id);
3711  if (!h) return 0;
3712  return h->GetXLeft(point);
3713  }
3721  {
3722  HolderP2P* h = FindP2P(id);
3723  if (!h) return 0;
3724  return h->GetXRight(point);
3725  }
3733  {
3734  HolderP2P* h = FindP2P(id);
3735  if (h) return h->GetYUp(point);
3736  HolderCommon* c = FindCommon(id);
3737  if (c) return c->GetYUp(GetY(point));
3738  return 0;
3739  }
3747  {
3748  HolderP2P* h = FindP2P(id);
3749  if (h) return h->GetYDown(point);
3750  HolderCommon* c = FindCommon(id);
3751  if (c) return c->GetYDown(GetY(point));
3752  return 0;
3753  }
3761  const char* GetSysTitle(Int_t id) const
3762  {
3763  Holder* h = Find(id);
3764  if (!h) return "";
3765  return h->GetTitle();
3766  }
3774  Style_t GetSysFillStyle(Int_t id) const
3775  {
3776  Holder* h = Find(id);
3777  if (!h) return 0;
3778  return h->GetFillStyle();
3779  }
3787  Style_t GetSysLineStyle(Int_t id) const
3788  {
3789  Holder* h = Find(id);
3790  if (!h) return 0;
3791  return h->GetLineStyle();
3792  }
3800  Color_t GetSysFillColor(Int_t id) const
3801  {
3802  Holder* h = Find(id);
3803  if (!h) return 0;
3804  return h->GetFillColor();
3805  }
3813  Color_t GetSysLineColor(Int_t id) const
3814  {
3815  Holder* h = Find(id);
3816  if (!h) return 0;
3817  return h->GetLineColor();
3818  }
3826  Width_t GetSysLineWidth(Int_t id) const
3827  {
3828  Holder* h = Find(id);
3829  if (!h) return 0;
3830  return h->GetLineWidth();
3831  }
3837  const char* GetSumTitle() const
3838  {
3839  return fSumTitle.Data();
3840  }
3841  UInt_t GetSumOption() const { return fSumOption; }
3847  Style_t GetSumFillStyle() const
3848  {
3849  return fSumFill.GetFillStyle();
3850  }
3856  Style_t GetSumLineStyle() const
3857  {
3858  return fSumLine.GetLineStyle();
3859  }
3865  Color_t GetSumFillColor() const
3866  {
3867  return fSumFill.GetFillColor();
3868  }
3874  Color_t GetSumLineColor() const
3875  {
3876  return fSumLine.GetLineColor();
3877  }
3883  Width_t GetSumLineWidth() const
3884  {
3885  return fSumLine.GetLineWidth();
3886  }
3892  const char* GetCommonSumTitle() const
3893  {
3894  return fCommonSumTitle.Data();
3895  }
3902  Style_t GetCommonSumFillStyle() const
3903  {
3904  return fCommonSumFill.GetFillStyle();
3905  }
3911  Style_t GetCommonSumLineStyle() const
3912  {
3913  return fCommonSumLine.GetLineStyle();
3914  }
3920  Color_t GetCommonSumFillColor() const
3921  {
3922  return fCommonSumFill.GetFillColor();
3923  }
3929  Color_t GetCommonSumLineColor() const
3930  {
3931  return fCommonSumLine.GetLineColor();
3932  }
3938  Width_t GetCommonSumLineWidth() const
3939  {
3940  return fCommonSumLine.GetLineWidth();
3941  }
3948  {
3949  HolderCommon* c = FindCommon(id);
3950  if (!c) return 0;
3951  Bool_t rel = c->IsRelative();
3952  return c->GetYUp(rel ? 100 : 1);
3953  }
3960  {
3961  HolderCommon* c = FindCommon(id);
3962  if (!c) return 0;
3963  Bool_t rel = c->IsRelative();
3964  return c->GetYDown(rel ? 100 : 1);
3965  }
3971  const char* GetXTitle() const { return fXTitle.Data(); }
3977  const char* GetYTitle() const { return fYTitle.Data(); }
3985  void GetMinMax(Option_t* option,
3986  Double_t& ymin, Double_t& ymax) const
3987  {
3988  Double_t xmin, xmax;
3989  Int_t imin, imax;
3990  GetMinMax(option, ymin, ymax, xmin, xmax, imin, imax);
3991  }
4003  void GetMinMax(Option_t* option,
4005  Double_t& xmin, Double_t& xmax,
4006  Int_t& imin, Int_t& imax) const
4007  {
4008  TString opt(option); opt.ToUpper();
4009  Bool_t cmn = opt.Contains("COMMON");
4010  Bool_t stat = opt.Contains("STAT");
4011  Bool_t quad = opt.Contains("QUAD");
4012  Bool_t noerr = opt.Contains("Y");
4013  quad = !opt.Contains("DIRECT");
4014 
4015  ymin = +1e9;
4016  ymax = -1e9;
4017  xmin = -1e9;
4018  xmax = +1e9;
4019  imin = -1;
4020  imax = -1;
4021  for (Int_t i = 0; i < GetN(); i++) {
4022  Double_t eyl = 0;
4023  Double_t eyh = 0;
4024  Double_t y = noerr ? GetY(i) : GetYandError(i, cmn, stat, quad, false, eyl, eyh);
4025  Double_t yl = y - eyl;
4026  Double_t yh = y + eyh;
4027  if (yl < ymin) {
4028  ymin = yl;
4029  xmin = GetX(i);
4030  imin = i;
4031  }
4032  if (yh > ymax) {
4033  ymax = yh;
4034  xmax = GetX(i);
4035  imax = i;
4036  }
4037  // Printf("%3d: %f -> %f +%f -%f -> %f %f -> %f %f",
4038  // i, GetX(i), y, eyl, eyh, yl, yh, ymin, ymax);
4039  }
4040  // Printf("min=(%d,%f,%f) max=(%d,%f,%f)", imin, xmin, ymin, imax, xmax, ymax);
4041  }
4055  Bool_t cmn, Bool_t stat, Bool_t quad,
4056  Int_t& i1, Int_t& i2) const
4057  {
4058  Int_t lim = (dir < 0) ? -1 : GetN();
4059  i1 = i2 = -1;
4060  for (Int_t i = start+dir; i != lim; i += dir) {
4061  Double_t eyl, eyh;
4062  Double_t y = GetYandError(i, cmn, stat, quad, false, eyl, eyh);
4063  // Printf("%3d: %2d %f -> %f -%f +%f (%f, %f, %f) -> %2d %2d",
4064  // i, dir, GetX(i), y, eyl, eyh, y-eyl, y+eyh, ymax/2,
4065  // i1, i2);
4066  if (y-eyl > ymax/2) continue; // Still above
4067  if (y+eyh < ymax/2) break; // Found lower limit
4068  if (i1 < 0) i1 = i2 = i;
4069  else i2 = i;
4070  }
4071  // Printf("Found %d,%d", i1, i2);
4072  }
4081  Double_t FWHM(Double_t& el, Double_t& eh) const
4082  {
4083  Double_t xl, xh;
4084  return FWHM(el, eh, xl, xh);
4085  }
4096  Double_t FWHM(Double_t& el, Double_t& eh, Double_t& xl, Double_t& xh) const
4097  {
4098  Double_t ymin, ymax, xmin, xmax;
4099  Int_t imin, imax;
4100  GetMinMax("quad stat", ymin, ymax, xmin, xmax, imin, imax);
4101  if (ymin > ymax/2) {
4102  Warning("FWHM", "Half of ymax=%f is out of reach [%f,%f]",
4103  ymax/2, ymin, ymax);
4104  return -1;
4105  }
4106 
4107  // Point bounds of found half-max
4108  Int_t li1, li2, ri1, ri2;
4109  // Search left
4110  FindFwhm(imax, -1, ymax, false, true, true, li1, li2);
4111  // Search right
4112  FindFwhm(imax, +1, ymax, false, true, true, ri1, ri2);
4113 
4114  Double_t lsign = 1;
4115  // Check if we've found no left point
4116  if (li1 < 0 && li2 < 0) {
4117  // Check if we've found no right point
4118  if (ri1 < 0 && ri2 < 0) {
4119  Warning("FWHM", "No left and right point found");
4120  return -1;
4121  }
4122  // Set sign, and copy right point to left point
4123  lsign = -1;
4124  li2 = ri1;
4125  li1 = ri2;
4126  }
4127  Double_t rsign = 1;
4128  // Check if we've found no right point
4129  if (ri1 < 0 && ri2 < 0) {
4130  // Check if we've found no left point
4131  if (li1 < 0 && li2 < 0) {
4132  Warning("FWHM", "No right and left point found");
4133  return -1;
4134  }
4135  // Set sign, and copy left point to right point
4136  rsign = -1;
4137  ri2 = li1;
4138  ri1 = li2;
4139  }
4140  // Printf("Points low=(%d,%d) high=(%d,%d)", li1, li2, ri1, ri2);
4141 
4142  // Now find best estimate of left X
4143  LinearSigmaCombiner scl;
4144  scl.Add(lsign*GetX(li1), GetErrorXLeft(li1), GetErrorXRight(li1));
4145  scl.Add(lsign*GetX(li2), GetErrorXLeft(li2), GetErrorXRight(li2));
4146  Combiner::Result* lr = scl.Calculate();
4147  scl.Print();
4148  lr->Print();
4149 
4150  // Now find best estimate of right X
4151  LinearSigmaCombiner scr;
4152  scr.Add(GetX(ri1), GetErrorXLeft(ri1), GetErrorXRight(ri1));
4153  scr.Add(GetX(ri2), GetErrorXLeft(ri2), GetErrorXRight(ri2));
4154  Combiner::Result* rr = scr.Calculate();
4155  scr.Print();
4156  rr->Print();
4157 
4158  // Get error on left/right hand sides
4159  el = lr->fEh+rr->fEl;
4160  eh = lr->fEl+rr->fEh;
4161  if (el < 1e-9 && eh < 1e-9) {
4162  // If the errors are very small, we take the average
4163  xl = (GetX(li2)+GetX(li1))/2;
4164  xh = (GetX(ri2)+GetX(ri1))/2;
4165  el = (GetX(li2)-GetX(li1))/2;
4166  eh = (GetX(ri2)-GetX(ri1))/2;
4167  }
4168  else {
4169  // Otherwise we set left/right hand values to best-fit values
4170  xl = lr->fX;
4171  xh = rr->fX;
4172  }
4173  // Evaluate the full-width half-maximum
4174  Double_t fwhm = (xh - xl);
4175  return fwhm;
4176  }
4177 
4195  Bool_t cmn=false,
4196  Bool_t stat=true,
4197  Bool_t quad=true) const
4198  {
4199  Double_t meanX=TMath::Infinity(), stdDev, n;
4200  StatisticsX(meanX, stdDev, n, cmn, stat, quad);
4201  error = stdDev/n;
4202  return meanX;
4203  }
4219  Double_t MeanX(Bool_t cmn=false,
4220  Bool_t stat=true,
4221  Bool_t quad=true) const
4222  {
4223  Double_t error;
4224  return MeanX(error, cmn, stat, quad);
4225  }
4243  Bool_t cmn=false,
4244  Bool_t stat=true,
4245  Bool_t quad=true) const
4246  {
4247  Double_t meanX=TMath::Infinity(), stdDev, eStdDev, n;
4248  StatisticsX(meanX, stdDev, n, cmn, stat, quad);
4249  error = stdDev/TMath::Sqrt(2)/n;
4250  return stdDev;
4251  }
4268  Bool_t stat=true,
4269  Bool_t quad=true) const
4270  {
4271  Double_t error;
4272  return StandardDeviationX(error, cmn, stat, quad);
4273  }
4275  Double_t& error,
4276  Bool_t cmn=false,
4277  Bool_t stat=true,
4278  Bool_t quad=true) const
4279  {
4280  Double_t meanX, stdDev, eStdDev, n;
4281  StatisticsX(meanX, stdDev, n, cmn, stat, quad);
4282  Double_t sd2 = stdDev*stdDev - meanX*meanX;
4283  stdDev = TMath::Sqrt(sd2 - mean*mean);
4284  error = stdDev/TMath::Sqrt(2)/n;
4285  return stdDev;
4286  }
4287 
4318  void StatisticsX(Double_t& meanX,
4319  Double_t& stdDevX,
4320  Double_t& n,
4321  Bool_t cmn=false,
4322  Bool_t stat=true,
4323  Bool_t quad=true) const
4324  {
4325  Double_t sumY = 0;
4326  Double_t sumXY = 0;
4327  Double_t sumXXY = 0;
4328  Double_t sumE2 = 0;
4329  for (Int_t i = 0; i < GetN(); i++) {
4330  Double_t e2yl, e2yh;
4331  Double_t x = GetX(i);
4332  Double_t y = GetYandError(i,cmn,stat,quad,true,e2yl,e2yh);
4333 
4334  sumY += y;
4335  sumXY += x*y;
4336  sumXXY += x*x*y;
4337  sumE2 += TMath::Max(e2yl,e2yh);
4338  }
4339  if (meanX == TMath::Infinity())
4340  meanX = sumXY/sumY;
4341  Double_t sd2 = TMath::Abs(sumXXY/sumY-meanX*meanX);
4342  stdDevX = TMath::Sqrt(sd2);
4343  Double_t eff = sumE2 > 0 ? sumY*sumY/sumE2 : sumY;
4344  n = TMath::Sqrt(eff);
4345  }
4346  /* @} */
4347 
4348  //__________________________________________________________________
4359  void SetSysLineColor(Int_t id, Color_t color) //*MENU*
4360  {
4361  Holder* h = Find(id);
4362  if (!h) return;
4363  h->SetLineColor(color);
4364  }
4371  void SetSysLineStyle(Int_t id, Style_t style) //*MENU*
4372  {
4373  Holder* h = Find(id);
4374  if (!h) return;
4375  h->SetLineStyle(style);
4376  }
4383  void SetSysLineWidth(Int_t id, Width_t width) //*MENU*
4384  {
4385  Holder* h = Find(id);
4386  if (!h) return;
4387  h->SetLineWidth(width);
4388  }
4395  void SetSysFillColor(Int_t id, Color_t color) //*MENU*
4396  {
4397  Holder* h = Find(id);
4398  if (!h) return;
4399  h->SetFillColor(color);
4400  }
4407  void SetSysFillStyle(Int_t id, Style_t style) //*MENU*
4408  {
4409  Holder* h = Find(id);
4410  if (!h) return;
4411  h->SetFillStyle(style);
4412  }
4419  void SetSysTitle(Int_t id, const char* name) //*MENU*
4420  {
4421  Holder* h = Find(id);
4422  if (!h) return;
4423  h->SetTitle(name);
4424  }
4432  {
4433  Holder* h = Find(id);
4434  if (!h) return;
4435  h->SetDOption(opt);
4436  }
4437  /* @} */
4438 
4439  //__________________________________________________________________
4455  void SetSumTitle(const char* title) { fSumTitle = title; } //*MENU*
4461  void SetSumLineColor(Color_t color) { fSumLine.SetLineColor(color); } //*MENU*
4467  void SetSumLineStyle(Style_t style){ fSumLine.SetLineStyle(style); } //*MENU*
4473  void SetSumLineWidth(Width_t width){ fSumLine.SetLineWidth(width); }//*MENU*
4479  void SetSumFillColor(Color_t color){ fSumFill.SetFillColor(color); }//*MENU*
4485  void SetSumFillStyle(Style_t style){ fSumFill.SetFillStyle(style); }//*MENU*
4497  void SetCommonSumTitle(const char* title) { fCommonSumTitle = title; } //*MENU*
4503  void SetCommonSumLineColor(Color_t color) { fCommonSumLine.SetLineColor(color); } //*MENU*
4509  void SetCommonSumLineStyle(Style_t style){ fCommonSumLine.SetLineStyle(style); } //*MENU*
4515  void SetCommonSumLineWidth(Width_t width){ fCommonSumLine.SetLineWidth(width); }//*MENU*
4521  void SetCommonSumFillColor(Color_t color){ fCommonSumFill.SetFillColor(color); }//*MENU*
4527  void SetCommonSumFillStyle(Style_t style){ fCommonSumFill.SetFillStyle(style); }//*MENU*
4528  /* @} */
4587  void SetKey(const char* key, const char* value, Bool_t replace=false) //*MENU*
4588  {
4589  if (!fMap) {
4590  fMap = new TList();
4591  fMap->SetOwner();
4592  fMap->SetName("keys");
4593  // fMap->SetTitle("key/value pairs");
4594  }
4595  TString k(key);
4596  TString v(value);
4597  TString t = v.Strip(TString::kBoth);
4598  v = t;
4599  if (k.EqualTo("experiment",TString::kIgnoreCase)) {
4600  TObjArray* l = v.Tokenize("-");
4601  t = Token(l, 0); TString lab = t.Strip(TString::kBoth);
4602  t = Token(l, 1); TString acc = t.Strip(TString::kBoth);
4603  t = Token(l, 2); TString exp = t.Strip(TString::kBoth);
4604  SetKey("laboratory", lab.Data(), replace);
4605  SetKey("accelerator", acc.Data(), replace);
4606  SetKey("detector", exp.Data(), replace);
4607  l->Delete();
4608  }
4609  if (k.EqualTo("comment", TString::kIgnoreCase)) {
4610  TString l(GetKey("laboratory"));
4611  TString a(GetKey("accelerator"));
4612  if (!l.IsNull() && !a.IsNull())
4613  v.Remove(0, l.Length()+1+a.Length()+1);
4614  t = v.Strip(TString::kBoth);
4615  fMap->Add(new TNamed("abstract", t.Data()));
4616  }
4617  else {
4618  if (replace) {
4619  TObjLink* cur = fMap->FirstLink();
4620  TObjLink* last = fMap->LastLink();
4621  while (cur) { // != last) {
4622  if (!k.EqualTo(cur->GetObject()->GetName())) {
4623  cur = cur->Next();
4624  continue;
4625  }
4626  TObjLink* tmp = cur->Next();
4627  fMap->Remove(cur);
4628  cur = tmp;
4629  }
4630  }
4631  fMap->Add(new TNamed(k, v));
4632  }
4633  }
4641  const char* GetKey(const char* key) const
4642  {
4643  if (!fMap) return 0;
4644  TObject* o = fMap->FindObject(key);
4645  if (!o) return 0;
4646  return o->GetTitle();
4647  }
4661  void CopyKeys(const GraphSysErr* g, Option_t* option="fr")
4662  {
4663  if (!g) return;
4664 
4665  TString opt(option);
4666  opt.ToLower();
4667 
4668  Bool_t all = opt.Contains("a");
4669  Bool_t file = opt.Contains("f");
4670  Bool_t header = opt.Contains("h");
4671  Bool_t qual = opt.Contains("q");
4672  Bool_t repl = opt.Contains("r");
4673 
4674  // With the "all" option we just do a plain copy
4675  TList* map = g->fMap;
4676  TIter nextKV(map);
4677  TObject* kv = 0;
4678  while ((kv = nextKV())) {
4679  if (!all) {
4680  TString key = kv->GetName();
4681  if (file && !(key.EqualTo("reference") ||
4682  key.EqualTo("laboratory") ||
4683  key.EqualTo("accelerator")||
4684  key.EqualTo("detector") ||
4685  key.EqualTo("abstract") ||
4686  key.EqualTo("author") ||
4687  key.EqualTo("doi") ||
4688  key.EqualTo("inspireId") ||
4689  key.EqualTo("cdsId") ||
4690  key.EqualTo("durhamId") ||
4691  key.EqualTo("title") ||
4692  key.EqualTo("status"))) continue;
4693  if (header && !(key.EqualTo("location") ||
4694  key.EqualTo("reackey") ||
4695  key.EqualTo("obskey") ||
4696  key.EqualTo("dscomment"))) continue;
4697  }
4698  SetKey(kv->GetName(), kv->GetTitle(), repl);
4699  }
4700 
4701  if (!all && !qual) return;
4702  TList* quals = g->fQualifiers;
4703  TIter nextQ(quals);
4704  TObject* qv = 0;
4705  while ((qv = nextQ())) {
4706  const char* test = GetQualifier(qv->GetName());
4707  if (!test || test[0] == '\0')
4708  AddQualifier(qv->GetName(), qv->GetTitle());
4709  }
4710  }
4711  void CopyAttr(const GraphSysErr* f)
4712  {
4713  fDataOption = f->fDataOption;
4714  SetMarkerColor(f->GetMarkerColor());
4715  SetMarkerStyle(f->GetMarkerStyle());
4716  SetMarkerSize(f->GetMarkerSize());
4717  SetLineColor(f->GetLineColor());
4718  SetLineStyle(f->GetLineStyle());
4719  SetLineWidth(f->GetLineWidth());
4720  SetFillColor(f->GetFillColor());
4721  SetFillStyle(f->GetFillStyle());
4722 
4723  fSumOption = f->fSumOption;
4724  SetSumTitle(f->GetSumTitle());
4730 
4738 
4739  SetXTitle(f->GetXTitle());
4740  SetYTitle(f->GetYTitle());
4741  }
4742 
4743 
4744 
4745 
4746  /* @} */
4747  //__________________________________________________________________
4759  void AddQualifier(const TString& key, const TString& value,
4760  Bool_t replace=false)
4761  {
4762  if (!fQualifiers) {
4763  fQualifiers = new TList();
4764  fQualifiers->SetOwner();
4765  fQualifiers->SetName("qualifiers");
4766  }
4767  TString val(value);
4768  if (key.EqualTo(value)) val = "";
4769 
4770  TString k = key.Strip(TString::kBoth, ' ');
4771  TObject* o = fQualifiers->FindObject(k);
4772  if (o) {
4773  Warning("AddQualifier", "Dataset already has qualifier \"%s\"",
4774  k.Data());
4775  if (replace) static_cast<TNamed*>(o)->SetTitle(value);
4776  return;
4777  }
4778 
4779  fQualifiers->Add(new TNamed(k, value));
4780  }
4786  void RemoveQualifier(const TString& key)
4787  {
4788  if (!fQualifiers) return;
4789  TObject* o = fQualifiers->FindObject(key);
4790  if (!o) return;
4791  fQualifiers->Remove(o);
4792  delete o;
4793  }
4801  const char* GetQualifier(const char* name) const
4802  {
4803  if (!fQualifiers) return "";
4804  TObject* o = fQualifiers->FindObject(name);
4805  if (!o) return "";
4806  return o->GetTitle();
4807  }
4808  /* @} */
4826  Bool_t cmn,
4827  Bool_t stat,
4828  Bool_t quad,
4829  Bool_t nosqrt,
4830  Double_t& y,
4831  Double_t& eyl,
4832  Double_t& eyh,
4833  Double_t& seyl,
4834  Double_t& seyh) const
4835  {
4836  Int_t i1, i2;
4837  Int_t ret = FindPoint(x, i1, i2);
4838  seyl = seyh = 0;
4839  if (ret < -1)
4840  // No valid point found
4841  return false;
4842  if (ret >= 0) {
4843  // Got exact point
4844  // Info("", "Get point %d (%f)", ret, x);
4845  y = GetYandError(ret, cmn, stat, quad, nosqrt, eyl, eyh);
4846  if (!stat) {
4847  // If stat error not include, get them here
4848  seyl = GetStatErrorDown(ret);
4849  seyh = GetStatErrorUp(ret);
4850  }
4851  return true;
4852  }
4853  // Interpolate between points
4854  Double_t eyl1, eyl2, eyh1, eyh2;
4855  Double_t x1 = fData->GetX()[i1];
4856  Double_t x2 = fData->GetX()[i2];
4857  Double_t y1 = GetYandError(i1, cmn, stat, quad, nosqrt, eyl1, eyh1);
4858  Double_t y2 = GetYandError(i2, cmn, stat, quad, nosqrt, eyl2, eyh2);
4859  Double_t seyl1 = GetStatErrorDown(i1);
4860  Double_t seyh1 = GetStatErrorUp(i1);
4861  Double_t seyl2 = GetStatErrorDown(i2);
4862  Double_t seyh2 = GetStatErrorUp(i2);
4863 
4864  // Linear interpolation
4865  Double_t dx = (x2-x1);
4866  Double_t ax = (x-x1)/dx;
4867  y = y1 + ax * (y2 - y1);
4868  eyl = eyl1 + ax * (eyl2 - eyl1);
4869  eyh = eyh1 + ax * (eyh2 - eyh1);
4870  if (!stat) {
4871  // if stat errors not included, calculate here
4872  seyl = seyl1 + ax * (seyl2 - seyl1);
4873  seyh = seyh1 + ax * (seyh2 - seyh1);
4874  }
4875  return true;
4876  }
4892  Bool_t cmn,
4893  Bool_t stat,
4894  Bool_t quad,
4895  Bool_t nosqrt,
4896  Double_t& y,
4897  Double_t& eyl,
4898  Double_t& eyh) const
4899  {
4900  Double_t seyl;
4901  Double_t seyh;
4902  return FindYandError(x, cmn, stat, quad, nosqrt, y, eyl, eyh, seyl, seyh);
4903  }
4904 
4905  //__________________________________________________________________
4911  struct Combiner : public TObject
4912  {
4913  typedef Double_t (*Wrapper_t)(Double_t*,Double_t*);
4917  struct Observation : public TObject
4918  {
4920  Double_t fX;
4922  Double_t fEl;
4924  Double_t fEh;
4932  Observation(Double_t x=0, Double_t el=0, Double_t eh=0)
4933  : fX(x), fEl(el), fEh(eh)
4934  {
4935  }
4939  virtual ~Observation() {}
4950  Double_t S() const
4951  {
4952  if (fEl * fEh == 0) return 0;
4953  return 2 * fEl * fEh / (fEl + fEh);
4954  }
4965  Double_t Sprime() const
4966  {
4967  if (TMath::Abs(fEh + fEl) < 1e-9) return 1;
4968  return (fEh - fEl) / (fEh + fEl);
4969  }
4970  Double_t Svar(Double_t guess=0) const
4971  {
4972  return S() + Sprime() * (guess - fX);
4973  }
4984  Double_t V() const { return fEl * fEh; }
4995  Double_t Vprime() const { return fEh - fEl; }
5001  virtual Double_t Low() const { return fX - 3 * fEl; }
5007  virtual Double_t High() const { return fX + 3 * fEh; }
5013  virtual void Print(Option_t* option="") const
5014  {
5015  Printf("%10.8f -%10.8f +%10.8f", fX, fEl, fEh);
5016  }
5017  ClassDef(Observation,1); // An experimental observation
5018  };
5019 
5023  struct Result : public Observation
5024  {
5026  Double_t fChi2;
5028  Double_t fLower;
5030  Double_t fUpper;
5041  Result(Double_t x=0, Double_t el=0, Double_t eh=0,
5042  Double_t chi2=0, Double_t low=0, Double_t high=0)
5043  : Observation(x, el, eh),
5044  fChi2(chi2),
5045  fLower(low),
5046  fUpper(high)
5047  {}
5053  Double_t Low() const { return fLower; }
5059  Double_t High() const { return fUpper; }
5065  virtual void Print(Option_t* option="") const
5066  {
5067  Printf("%10.8f -%10.8f +%10.8f %5.2f", fX, fEl, fEh, fChi2);
5068  }
5069  ClassDef(Result,1); // Final result
5070  };
5071  TClonesArray fData; // Container of observations
5072  Result* fResult; // Result of combining
5073 
5078  : fData("GraphSysErr::Combiner::Observation"), fResult(0)
5079  {
5080  fData.SetOwner();
5081  }
5085  virtual ~Combiner() { fData.Clear(); if (fResult) delete fResult; }
5095  void Clear(Option_t* option="")
5096  {
5097  fData.Clear();
5098  if (fResult) delete fResult;
5099  fResult = 0;
5100  }
5106  void Add(const Observation& r)
5107  {
5108  Add(r.fX,r.fEl,r.fEh);
5109  }
5117  void Add(Double_t x, Double_t el, Double_t eh)
5118  {
5119  // if (TMath::Abs(x) < 1e-10 &&
5120  // TMath::Abs(el) < 1e-10 &&
5121  // TMath::Abs(eh) < 1e-10) return;
5122  new (fData[fData.GetEntries()]) Observation(x,el,eh);
5123  }
5129  void Print(Option_t* option="") const
5130  {
5131  TIter next(&fData);
5132  Observation* r = 0;
5133  while ((r = static_cast<Observation*>(next())))
5134  r->Print(option);
5135  }
5136  /* @} */
5137 
5145  virtual Double_t W(const Observation& r) const = 0;
5154  virtual Double_t StepW(Double_t guess, const Observation& r) const = 0;
5160  virtual Double_t StepOffset(Double_t guess, const Observation& r) const = 0;
5167  virtual Double_t VarTerm(Double_t guess, const Observation& r) const = 0;
5183  Double_t ChiTerm(Double_t guess, const Observation& r) const
5184  {
5185  Double_t var = VarTerm(guess, r);
5186  if (var <= 0) return -1000;
5187 
5188  return TMath::Power(guess - r.fX, 2) / var;
5189  }
5199  Double_t F(Double_t guess,
5200  Double_t chi2) const
5201  {
5202  if (fData.GetEntries() == 1) return 1;
5203  Double_t s = -chi2;
5204 
5205  TIter next(&fData);
5206  Observation* r = 0;
5207  while ((r = static_cast<Observation*>(next())))
5208  s += ChiTerm(guess, *r);
5209  return s;
5210  }
5222  Double_t E(UShort_t nIter,
5223  Int_t sign,
5224  Double_t best,
5225  Double_t chi2,
5226  Double_t s)
5227  {
5228  if (fData.GetEntries() == 1) {
5229  Observation* o = static_cast<Observation*>(fData[0]);
5230  return (sign < 0 ? o->fEl : o->fEh);
5231  }
5232 
5233  // Step size
5234  Double_t delta = 0.1 * sign * s;
5235 
5236  // Iterations
5237  for (UShort_t i = 0; i < nIter; i++) {
5238  // Calculate chi^2 with current guess
5239  Double_t got = F(best + sign * s, chi2);
5240 
5241  if (TMath::Abs(got-1) < 1e-7)
5242  // We're close to 1 so get out
5243  break;
5244 
5245  // The next guess' chi^2 value e
5246  Double_t guess = F(best + sign * s + delta, chi2);
5247 
5248  // Where to search next
5249  if ((got - 1) * (guess - 1) > 0) {
5250  if ((got - 1) / (guess - 1) < 1)
5251  delta = -delta;
5252  else
5253  s += sign * delta;
5254  continue;
5255  }
5256 
5257  // Update error value and decrease step size
5258  s += sign * delta * (1 - got) / (guess - got);
5259  delta /= 2;
5260  }
5261  return s;
5262  }
5272  Double_t X(UShort_t nIter,
5273  Double_t lowest,
5274  Double_t highest)
5275  {
5276  // Starting values
5277  if (fData.GetEntries() == 1)
5278  return static_cast<Observation*>(fData[0])->fX;
5279  Double_t x = (highest+lowest)/2;
5280  Double_t oldX = -1e33;
5281  // Do the iterations
5282  for (UShort_t i = 0; i < nIter; i++) {
5283  Double_t sum = 0;
5284  Double_t sumw = 0;
5285  Double_t offset = 0;
5286 
5287  // Loop over observations
5288  TIter next(&fData);
5289  Observation* r = 0;
5290  while ((r = static_cast<Observation*>(next()))) {
5291  Double_t w = StepW(x, *r);
5292  offset += StepOffset(x, *r);
5293  sum += r->fX * w;
5294  sumw += w;
5295  }
5296  x = (TMath::Abs(sumw) > 1e-9) ? (sum - offset) / sumw : 0;
5297 
5298  if (TMath::Abs(x - oldX) < (highest-lowest) * 1e-9) break;
5299  oldX = x;
5300  }
5301  return x;
5302  }
5311  {
5312  Double_t lowest = +1e33;
5313  Double_t highest = -1e33;
5314  Double_t sumLow = 0;
5315  Double_t sumHigh = 0;
5316 
5317  // Find boundaries and sum weights
5318  TIter next(&fData);
5319  Observation* r = 0;
5320  while ((r = static_cast<Observation*>(next()))) {
5321  lowest = TMath::Min(r->Low(), lowest);
5322  highest = TMath::Max(r->High(), highest);
5323  if (r->fEl > 1e-10) sumLow += 1./TMath::Power(r->fEl, 2);
5324  if (r->fEh > 1e-10) sumHigh += 1./TMath::Power(r->fEh, 2);
5325  }
5326  // Summed weights
5327  Double_t sLow = sumLow > 1e-10 ? 1. / TMath::Sqrt(sumLow) : 0;
5328  Double_t sHigh = sumHigh > 1e-10 ? 1. / TMath::Sqrt(sumHigh) : 0;
5329 
5330  // Now do the calculations
5331  Double_t bestX = X(nIter, lowest, highest);
5332  Double_t bestChi2 = F(bestX, 0);
5333  Double_t bestLow = TMath::Abs(E(nIter,-1,bestX,bestChi2,sLow));
5334  Double_t bestHigh = TMath::Abs(E(nIter,+1,bestX,bestChi2,sHigh));
5335 
5336  fResult = new Result(bestX, bestLow, bestHigh, bestChi2, lowest, highest);
5337  return fResult;
5338  }
5344  virtual Wrapper_t Wrapper() const = 0;
5345 
5355  TF1* MakeF(const Observation& r, Int_t j) const
5356  {
5357  TF1* f = new TF1(Form("f%02d", j), Wrapper(), r.Low(), r.High(), 3);
5358  f->SetParNames("x", "#sigma^{-}", "#sigma^{+}");
5359  f->SetParameters(r.fX, r.fEl, r.fEh);
5360  f->SetLineStyle((j%3)+2);
5361  f->SetLineColor(kBlack);
5362  f->SetLineWidth(2);
5363  return f;
5364  }
5372  TLine* MakeL(TF1* f) const
5373  {
5374  Double_t m = f->GetParameter(0);
5375  TLine* l = new TLine(m-f->GetParameter(1), 1,
5376  m+f->GetParameter(2), 1);
5377  l->SetLineColor(f->GetLineColor());
5378  l->SetLineStyle(f->GetLineStyle());
5379  l->SetLineWidth(f->GetLineWidth());
5380  return l;
5381  }
5382  void Draw(Option_t* option="")
5383  {
5384  if (!fResult) return;
5385 
5386  TList fs; fs.SetOwner(false);
5387  Int_t j = 0;
5388  TIter next(&fData);
5389  Observation* r = 0;
5390  while ((r = static_cast<Observation*>(next()))) {
5391  TF1* f = MakeF(*r, j);
5392  f->SetRange(fResult->Low(), fResult->High());
5393  fs.Add(f, j == 0 ? "" : "same");
5394  fs.Add(MakeL(f));
5395  j++;
5396  }
5397  TF1* fr = MakeF(*fResult, j);
5398  fr->SetLineColor(kRed+2);
5399  fr->SetLineStyle(1);
5400  fs.Add(fr);
5401  fs.Add(MakeL(fr));
5402  TIter nextD(&fs);
5403  TObject* o = 0;
5404  j = 0;
5405  while((o = nextD())) { o->Draw(j == 0 ? "" : "same"); j++; }
5406  // fs.Draw();
5407  static_cast<TF1*>(fs.First())->GetHistogram()
5408  ->GetYaxis()->SetRangeUser(-.1, 5.1);
5409  }
5410  ClassDef(Combiner,1); // Combine systematics
5411  };
5412 
5413 
5418  {
5430  Double_t W(const Observation& r) const
5431  {
5432  // Double_t s = r.S();
5433  // if (s < 1e-10) return 0;
5434  // return .5 * TMath::Power(r.Svar(), 3) / s;
5435  Double_t w = StepW(0,r);
5436  return 1/w;
5437  }
5450  Double_t StepW(Double_t guess, const Observation& r) const
5451  {
5452  Double_t s = r.S();
5453  Double_t t = r.Svar(guess);
5454  Double_t ret = s / TMath::Power(t,3);
5455  return ret;
5456  }
5463  {
5464  return 0;
5465  }
5479  Double_t VarTerm(Double_t guess, const Observation& r) const
5480  {
5481  return TMath::Power(r.Svar(guess),2);
5482  }
5505  static Double_t L(Double_t guess, Double_t x, Double_t el, Double_t eh)
5506  {
5507  Observation tmp(x,el,eh);
5508  Double_t d = (guess-x);
5509  Double_t t = tmp.Svar(guess); // (s+sp*d);
5510  if (TMath::Abs(d) < 1e-10) return 0;
5511  // if (TMath::Abs(t) < 1e-10) return DBL_MAX;
5512  Double_t ret = TMath::Power(d/t, 2);
5513  return ret;
5514  }
5523  static Double_t WrapL(Double_t* xp, Double_t* pp)
5524  {
5525  return L(xp[0], pp[0], pp[1], pp[2]);
5526  }
5532  Wrapper_t Wrapper() const { return WrapL; }
5533  ClassDef(LinearSigmaCombiner,1); // Sigma combiner
5534  };
5539  {
5551  Double_t W(const Observation& r) const
5552  {
5553  Double_t v = r.V();
5554  Double_t vp = r.Vprime();
5555  return TMath::Power(v + r.fX * vp, 2) / (2 * v + r.fX * vp);
5556  }
5569  Double_t StepW(Double_t guess, const Observation& r) const
5570  {
5571  Double_t v = r.V();
5572  return v / TMath::Power(v+r.Vprime()*(guess - r.fX), 2);
5573  }
5586  Double_t StepOffset(Double_t guess, const Observation& r) const
5587  {
5588  Double_t vp = r.Vprime();
5589  return 0.5 * vp * TMath::Power((guess-r.fX)/(r.V()+vp*(guess-r.fX)),2);
5590  }
5604  Double_t VarTerm(Double_t guess, const Observation& r) const
5605  {
5606  return r.V() + r.Vprime() * (guess - r.fX);
5607  }
5627  static Double_t L(Double_t guess, Double_t x, Double_t el, Double_t eh)
5628  {
5629  Double_t d = (guess-x);
5630  Double_t v = eh * el;
5631  Double_t vp = eh - el;
5632  return TMath::Power(d,2) / (v+vp*d);
5633  }
5642  static Double_t WrapL(Double_t* xp, Double_t* pp)
5643  {
5644  return L(xp[0], pp[0], pp[1], pp[2]);
5645  }
5651  Wrapper_t Wrapper() const { return WrapL; }
5652  ClassDef(LinearVarianceCombiner,1); // Variance combiner
5653  };
5654 
5655  //__________________________________________________________________
5659  struct Holder : public TNamed, TAttLine, TAttFill
5660  {
5662  friend struct GraphSysErr;
5667  : TNamed(),
5668  TAttLine(),
5669  TAttFill(),
5670  fRelative(true),
5671  fOption(kRect)
5672  {}
5676  virtual ~Holder() {}
5677  void CopyAttr(Holder* h)
5678  {
5679  SetLineColor(h->GetLineColor());
5680  SetLineStyle(h->GetLineStyle());
5681  SetLineWidth(h->GetLineWidth());
5682  SetFillColor(h->GetFillColor());
5683  SetFillStyle(h->GetFillStyle());
5684  fOption = h->fOption;
5685  }
5686  protected:
5696  Holder(const char* name, const char* title, Bool_t rel,
5697  UInt_t option, UInt_t id)
5698  : TNamed(name, title),
5699  TAttLine(0,0,0),
5700  TAttFill(0,0),
5701  fRelative(rel),
5702  fOption(option)
5703  {
5704  SetUniqueID(id);
5705  }
5711  Holder(const Holder& other)
5712  : TNamed(other),
5713  TAttLine(other),
5714  TAttFill(other),
5715  fRelative(other.fRelative),
5716  fOption(other.fOption)
5717  {
5718  SetUniqueID(other.GetUniqueID());
5719  }
5727  Holder& operator=(const Holder& other)
5728  {
5729  if (&other == this) return *this;
5730 
5731  other.TNamed::Copy(*this);
5732  other.TAttLine::Copy(*this);
5733  other.TAttFill::Copy(*this);
5734 
5735  fRelative = other.fRelative;
5736  fOption = other.fOption;
5737 
5738  SetUniqueID(other.GetUniqueID());
5739 
5740  return *this;
5741  }
5742 
5752  virtual Graph* StackError(Graph* g, Bool_t ignoreErr, Bool_t quad) const = 0;
5762  virtual void SumError(Graph* g, Int_t i, Bool_t ignoreErr,
5763  Bool_t quad, UInt_t opt) const = 0;
5767  virtual UInt_t GetDOption() const { return fOption; }
5773  virtual void SetDOption(EDrawOption_t opt) { fOption = opt; }
5779  virtual Bool_t IsRelative() const { return fRelative; }
5780 
5781  virtual void Print(Option_t* option="") const
5782  {
5783  gROOT->IndentLevel();
5784  Printf("%s/%s %s(ID # %d%s)", GetName(), GetTitle(),
5785  IsRelative() ? "PCT " : "", GetUniqueID(),
5786  TestBit(kUsedBit) ? " used" : "");
5787  TString opt(option);
5788  if (!opt.Contains("ATTR", TString::kIgnoreCase)) return;
5789  gROOT->IndentLevel();
5790  Printf(" [option: %3d line (c/s/w):%3d/%1d/%2d fill (c/s):%3d/%4d]",
5791  fOption, GetLineColor(), GetLineStyle(), GetLineWidth(),
5792  GetFillColor(), GetFillStyle());
5793 
5794  }
5795  virtual void ls(Option_t* option) const
5796  {
5797  Print(option);
5798  }
5799  protected:
5800  UShort_t XMode(Int_t opt=-1) const
5801  {
5802  UShort_t xMode = 0;
5803  if (opt < 0) opt = fOption;
5804  switch (opt) {
5805  case kNormal: case kNoTick: case kFill: case kCurve: xMode = 0; break;
5806  case kArrow: case kHat: case kBar: xMode = 1; break;
5807  case kNone: xMode = 0; break;
5808  case kRect: case kBox: xMode = 3; break;
5809  }
5810  return xMode;
5811  }
5828  void DoAdd(UShort_t xMode,
5829  Double_t curExl,
5830  Double_t curExh,
5831  Double_t curEyl,
5832  Double_t curEyh,
5833  Bool_t ignoreErr,
5834  Bool_t quad,
5835  Bool_t sqOld,
5836  Double_t& exl,
5837  Double_t& exh,
5838  Double_t& eyl,
5839  Double_t& eyh) const
5840  {
5841  // Double_t oexl = exl;
5842  // Double_t oexh = exh;
5843  if (quad) {
5844  eyl *= eyl;
5845  eyh *= eyh;
5846  if (sqOld) {
5847  curEyl *= curEyl;
5848  curEyh *= curEyh;
5849  }
5850  }
5851 
5852  if (!ignoreErr) {
5853  if (kVerbose & kDraw)
5854  Info("", "old: +%f/-%f this: +%f/-%f -> +%f/-%f",
5855  curEyl, curEyh, eyl, eyh, eyl + curEyl, eyh + curEyh);
5856  exl = TMath::Max(exl, curExl);
5857  exh = TMath::Max(exh, curExh);
5858  eyl += curEyl;
5859  eyh += curEyh;
5860  }
5861 
5862  switch (xMode) {
5863  case 0: break; // kNormal, kNoTick, kFill, kCurve, kNone
5864  case 1: exl = exh = 0; break; // kArrow, kHat, kBar
5865  case 2: // kRect, kBox
5866  if (exl <= 0) exl = gStyle->GetErrorX()/2;
5867  if (exh <= 0) exh = gStyle->GetErrorX()/2;
5868  break;
5869  }
5870  }
5871 
5877  void SetAttributes(Graph* g) const
5878  {
5879  if (!g) return;
5880  this->TAttLine::Copy(*g);
5881  this->TAttFill::Copy(*g);
5882  g->SetMarkerColor(0);
5883  g->SetMarkerStyle(0);
5884  g->SetMarkerSize(0);
5885  }
5890  ClassDef(Holder,3);
5891  };
5892  //__________________________________________________________________
5896  struct HolderP2P : public Holder
5897  {
5899  friend struct GraphSysErr;
5904  : Holder(),
5905  fGraph(0)
5906  {}
5907  protected:
5917  HolderP2P(const char* name, const char* title, Bool_t rel,
5918  UInt_t opt, UInt_t id)
5919  : Holder(name, title, rel, opt, id),
5920  fGraph(0)
5921  {
5922  fGraph = new Graph();
5923  fGraph->SetName(name);
5924  fGraph->SetTitle(title);
5925  SetAttributes(fGraph);
5926  }
5932  HolderP2P(const HolderP2P& other)
5933  : Holder(other),
5934  fGraph(0)
5935  {
5936  if (other.fGraph)
5937  fGraph = static_cast<Graph*>(other.fGraph->Clone());
5938 
5939  if (fGraph) SetAttributes(fGraph);
5940  }
5949  {
5950  if (&other == this) return *this;
5951 
5952  this->Holder::operator=(other);
5953  if (fGraph) delete fGraph;
5954  fGraph = 0;
5955  if (other.fGraph)
5956  fGraph = static_cast<Graph*>(other.fGraph->Clone());
5957  if (fGraph) SetAttributes(fGraph);
5958 
5959  return *this;
5960  }
5961 
5974  void Set(Int_t point, Graph* g, Double_t ex, Double_t ey)
5975  {
5976  Set(point, g, ex, ex, ey, ey);
5977  }
5988  void Set(Int_t point, Graph* g, Double_t ex1, Double_t ex2,
5989  Double_t ey1, Double_t ey2)
5990  {
5991  if (!fGraph) return;
5992  if (!g) return;
5993  if (point >= g->GetN()) return;
5994  Double_t x = g->GetX()[point];
5995  Double_t y = g->GetY()[point];
5996  fGraph->SetPoint(point, x, y);
5997  if (fRelative) { ey1 *= y; ey2 *= y; }
5998  fGraph->SetPointError(point, ex1, ex2, ey1, ey2);
5999  }
6000  /* @} */
6001 
6013  Double_t GetX(Int_t point) const
6014  {
6015  if (!fGraph) return 0;
6016  return fGraph->GetErrorX(point);
6017  }
6025  Double_t GetXLeft(Int_t point) const
6026  {
6027  if (!fGraph) return 0;
6028  return fGraph->GetErrorXlow(point);
6029  }
6037  Double_t GetXRight(Int_t point) const
6038  {
6039  if (!fGraph) return 0;
6040  return fGraph->GetErrorXhigh(point);
6041  }
6049  Double_t GetY(Int_t point) const
6050  {
6051  if (!fGraph) return 0;
6052  return fGraph->GetErrorY(point);
6053  }
6063  Double_t GetYDown(Int_t i1, Int_t i2, Double_t ax) const
6064  {
6065  if (!fGraph) return 0;
6066  if (i1 == i2) return GetYDown(i1);
6067  Double_t e1 = GetYDown(i1);
6068  Double_t e2 = GetYDown(i2);
6069  return e1 + ax * (e2-e1); // Linear interpolation
6070  }
6078  Double_t GetYDown(Int_t point) const
6079  {
6080  if (!fGraph) return 0;
6081  return fGraph->GetErrorYlow(point);
6082  }
6092  Double_t GetYUp(Int_t i1, Int_t i2, Double_t ax) const
6093  {
6094  if (!fGraph) return 0;
6095  if (i1 == i2) return GetYUp(i1);
6096  Double_t e1 = GetYUp(i1);
6097  Double_t e2 = GetYUp(i2);
6098  return e1 + ax * (e2-e1); // Linear interpolation
6099  }
6107  Double_t GetYUp(Int_t point) const
6108  {
6109  if (!fGraph) return 0;
6110  return fGraph->GetErrorYhigh(point);
6111  }
6112  /* @} */
6113 
6131  void AddError(Int_t i,
6132  UShort_t xMode,
6133  Bool_t ignoreErr,
6134  Bool_t quad,
6135  Bool_t sqOld,
6136  Double_t& exl,
6137  Double_t& exh,
6138  Double_t& eyl,
6139  Double_t& eyh) const
6140  {
6141  Double_t exlO = fGraph->GetErrorXlow(i);
6142  Double_t exhO = fGraph->GetErrorXhigh(i);
6143  Double_t eylO = fGraph->GetErrorYlow(i);
6144  Double_t eyhO = fGraph->GetErrorYhigh(i);
6145  Double_t oyl = eylO;
6146  Double_t oyh = eyhO;
6147  if (exl <= 0) exlO = exl;
6148  if (exh <= 0) exhO = exh;
6149 
6150  DoAdd(xMode,
6151  exl, exh, eyl, eyh,
6152  ignoreErr, quad, sqOld,
6153  exlO, exhO, eylO, eyhO);
6154  if (kVerbose & kDraw)
6155  Info(GetTitle(), "quad: %s this: +%f/-%f, old: +%f/-%f -> +%f/-%f",
6156  (quad ? "true" : "false"), oyl, oyh, eyl, eyh, eylO, eyhO);
6157 
6158  exl = exlO;
6159  exh = exhO;
6160  eyl = eylO;
6161  eyh = eyhO;
6162  }
6176  UShort_t xMode,
6177  Bool_t ignoreErr,
6178  Bool_t quad,
6179  Double_t& exl,
6180  Double_t& exh,
6181  Double_t& eyl,
6182  Double_t& eyh) const
6183  {
6184  Double_t oldEyl = eyl;
6185  Double_t oldEyh = eyh;
6186  AddError(i, xMode, ignoreErr, quad, true, exl, exh, eyl, eyh);
6187  if (kVerbose & kDraw)
6188  Info(GetTitle(), "old= +%f/-%f -> +%f/-%f", oldEyl, oldEyh, eyl, eyh);
6189  }
6199  Graph* StackError(Graph* g, Bool_t ignoreErr, Bool_t quad) const
6200  {
6201  Graph* ga = new Graph(g->GetN());
6202  for (Int_t i = 0; i < g->GetN(); i++) {
6203  Double_t x = g->GetX()[i];
6204  Double_t y = g->GetY()[i];
6205  Double_t exl = g->GetErrorXlow(i);
6206  Double_t exh = g->GetErrorXhigh(i);
6207  Double_t eyl = g->GetErrorYlow(i);
6208  Double_t eyh = g->GetErrorYhigh(i);
6209  StackPointError(i, XMode(), ignoreErr, quad, exl, exh, eyl, eyh);
6210  // AddError(i, g, quad, ignoreErr, true, exl, exh, eyl, eyh);
6211  ga->SetPoint(i, x, y);
6212  ga->SetPointError(i, exl,exh,eyl,eyh);
6213  }
6214  SetAttributes(ga);
6215  ga->SetTitle(GetTitle());
6216  ga->SetName(Form("stack_%08x", GetUniqueID()));
6217  return ga;
6218  }
6231  void SumPointError(Int_t i, UShort_t xMode, Bool_t ignoreErr, Bool_t quad,
6232  Double_t& exl,
6233  Double_t& exh,
6234  Double_t& eyl,
6235  Double_t& eyh) const
6236  {
6237  AddError(i, xMode, ignoreErr, quad, false, exl, exh, eyl, eyh);
6238  }
6248  void SumError(Graph* g, Int_t i, Bool_t ignoreErr,
6249  Bool_t quad, UInt_t opt) const
6250  {
6251  Double_t exl = (g ? g->GetErrorXlow(i) : 0);
6252  Double_t exh = (g ? g->GetErrorXhigh(i) : 0);
6253  Double_t eyl = (g ? g->GetErrorYlow(i) : 0);
6254  Double_t eyh = (g ? g->GetErrorYhigh(i) : 0);
6255  SumPointError(i, XMode(opt), ignoreErr, quad, exl, exh, eyl, eyh);
6256  g->SetPointError(i, exl,exh,eyl,eyh);
6257  }
6258  /* @} */
6259  void SavePrimitive(std::ostream& out, Option_t* option="")
6260  {
6261  if (option[0] == 'd' || option[0] == 'D')
6262  out << " // Point-2-Point systematic " << GetTitle() << "\n"
6263  << " {\n"
6264  << " Int_t id = g->DeclarePoint2Point(\"" << GetTitle() << "\","
6265  << IsRelative() << ',' << fOption << ");\n"
6266  << " g->SetSysLineColor(id," << GetLineColor() << ");\n"
6267  << " g->SetSysLineStyle(id," << GetLineStyle() << ");\n"
6268  << " g->SetSysLineWidth(id," << GetLineWidth() << ");\n"
6269  << " g->SetSysFillColor(id," << GetFillColor() << ");\n"
6270  << " g->SetSysFillStyle(id," << GetFillStyle() << ");\n"
6271  << " }\n";
6272  }
6273  virtual void Print(Option_t* option="") const
6274  {
6275  gROOT->IndentLevel();
6276  Printf("%s/%s (ID # %d, %d points)", GetName(), GetTitle(), GetUniqueID(),
6277  (fGraph ? fGraph->GetN() : -1));
6278  TString opt(option);
6279  if (!opt.Contains("ATTR", TString::kIgnoreCase)) return;
6280  gROOT->IndentLevel();
6281  Printf(" [option: %3d line (c/s/w):%3d/%1d/%2d fill (c/s):%3d/%4d]",
6282  fOption, GetLineColor(), GetLineStyle(), GetLineWidth(),
6283  GetFillColor(), GetFillStyle());
6284 
6285  }
6287  // Graph* fGraph;
6289 
6290  ClassDef(HolderP2P,3);
6291  };
6292  //__________________________________________________________________
6297  struct HolderCommon : public Holder
6298  {
6300  friend struct GraphSysErr;
6305  : Holder(), fEyl(0), fEyh(0)
6306  {}
6307  protected:
6317  HolderCommon(const char* name, const char* title, Bool_t rel,
6318  UInt_t opt, UInt_t id)
6319  : Holder(name, title, rel, opt, id), fEyl(0), fEyh(0)
6320  {
6321  }
6328  : Holder(other),
6329  fEyl(other.fEyl),
6330  fEyh(other.fEyh)
6331  {
6332  }
6341  {
6342  if (&other == this) return *this;
6343  this->Holder::operator=(other);
6344  fEyl = other.fEyl;
6345  fEyh = other.fEyh;
6346  return *this;
6347  }
6348 
6358  void Set(Double_t ey) { Set(ey, ey); }
6365  void Set(Double_t eyl, Double_t eyh)
6366  {
6367  fEyl = eyl;
6368  fEyh = eyh;
6369  }
6378  {
6379  return (fRelative ? y : 1) * fEyl;
6380  }
6389  {
6390  return (fRelative ? y : 1) * fEyh;
6391  }
6392  /* @} */
6393 
6413  UShort_t xMode,
6414  Bool_t ignoreErr,
6415  Bool_t quad,
6416  Bool_t sqOld,
6417  Double_t& exl,
6418  Double_t& exh,
6419  Double_t& eyl,
6420  Double_t& eyh) const
6421  {
6422  //exl = (i == 0 ? 0 :(g->GetX()[i] -g->GetX()[i-1])/2);
6423  //exh = (i+1 >= g->GetN()? 0 :(g->GetX()[i+1]-g->GetX()[i])/2);
6424  Double_t exlO = exl;
6425  Double_t exhO = exh;
6426  Double_t eylO = GetYDown(y);
6427  Double_t eyhO = GetYUp(y);
6428  if (exlO <= 0) exlO = gStyle->GetErrorX()/2;
6429  if (exhO <= 0) exhO = gStyle->GetErrorX()/2;
6430 
6431  DoAdd(xMode, exl, exh, eyl, eyh,
6432  ignoreErr, quad, sqOld,
6433  exlO, exhO, eylO, eyhO);
6434  exl = exlO;
6435  exh = exhO;
6436  eyl = eylO;
6437  eyh = eyhO;
6438  }
6449  Graph* BarError(Graph* g, Bool_t quad, Double_t x, Double_t y) const
6450  {
6451  Double_t exl = (g ? g->GetErrorXlow(0) : 0);
6452  Double_t exh = (g ? g->GetErrorXhigh(0) : 0);
6453  Double_t eyl = (g ? g->GetErrorYlow(0) : 0);
6454  Double_t eyh = (g ? g->GetErrorYhigh(0) : 0);
6455  AddError(y, XMode(), false, quad, quad, exl, exh, eyl, eyh);
6456 
6457  Graph* r = new Graph(1);
6458  r->SetPoint(0, x, y);
6459  r->SetPointError(0, exl, exh, eyl, eyh);
6460 
6461  SetAttributes(r);
6462  r->SetTitle(GetTitle());
6463  r->SetName(Form("bar_%08x", GetUniqueID()));
6464 
6465  return r;
6466  }
6480  UShort_t xMode,
6481  Bool_t ignoreErr,
6482  Bool_t quad,
6483  Double_t& exl,
6484  Double_t& exh,
6485  Double_t& eyl,
6486  Double_t& eyh) const
6487  {
6488  AddError(y, xMode, ignoreErr, quad, true, exl, exh, eyl, eyh);
6489  }
6490 
6500  Graph* StackError(Graph* g, Bool_t ignoreErr, Bool_t quad) const
6501  {
6502  Graph* ga = new Graph(g->GetN());
6503  for (Int_t i = 0; i < g->GetN(); i++) {
6504  Double_t x = g->GetX()[i];
6505  Double_t y = g->GetY()[i];
6506  Double_t exl = g->GetErrorXlow(i);
6507  Double_t exh = g->GetErrorXhigh(i);
6508  Double_t eyl = g->GetErrorYlow(i);
6509  Double_t eyh = g->GetErrorYhigh(i);
6510  StackPointError(y, XMode(), ignoreErr, quad, exl, exh, eyl, eyh);
6511  ga->SetPoint(i, x, y);
6512  ga->SetPointError(i, exl,exh,eyl,eyh);
6513  }
6514  SetAttributes(ga);
6515  ga->SetTitle(GetTitle());
6516  ga->SetName(Form("stack_%08x", GetUniqueID()));
6517  return ga;
6518  }
6531  void SumPointError(Double_t y, UShort_t xMode, Bool_t ignoreErr,Bool_t quad,
6532  Double_t& exl,
6533  Double_t& exh,
6534  Double_t& eyl,
6535  Double_t& eyh) const
6536  {
6537  AddError(y, xMode, ignoreErr, quad, false, exl, exh, eyl, eyh);
6538  }
6548  void SumError(Graph* g, Int_t i, Bool_t ignoreErr, Bool_t quad,
6549  UInt_t opt) const
6550  {
6551  Double_t y = (g ? g->GetY()[i] : 0);
6552  Double_t exl = (g ? g->GetErrorXlow(i) : 0);
6553  Double_t exh = (g ? g->GetErrorXhigh(i) : 0);
6554  Double_t eyl = (g ? g->GetErrorYlow(i) : 0);
6555  Double_t eyh = (g ? g->GetErrorYhigh(i) : 0);
6556  SumPointError(y, XMode(opt), ignoreErr, quad, exl, exl, eyl, eyh);
6557  g->SetPointError(i, exl,exh,eyl,eyh);
6558  }
6559  /* @} */
6560  void SavePrimitive(std::ostream& out, Option_t* option="")
6561  {
6562  if (option[0] == 'd' || option[0] == 'D')
6563  out << " // Common systematic " << GetTitle() << "\n"
6564  << " {\n"
6565  << " Int_t id = g->DefineCommon(\"" << GetTitle() << "\","
6566  << fRelative << ',' << fEyl << ',' << fEyh << ','
6567  << fOption << ");\n"
6568  << " g->SetSysLineColor(id," << GetLineColor() << ")\n;"
6569  << " g->SetSysLineStyle(id," << GetLineStyle() << ")\n;"
6570  << " g->SetSysLineWidth(id," << GetLineWidth() << ")\n;"
6571  << " g->SetSysFillColor(id," << GetFillColor() << ")\n;"
6572  << " g->SetSysFillStyle(id," << GetFillStyle() << ")\n;"
6573  << " }\n";
6574  }
6575  virtual void Print(Option_t* option="") const
6576  {
6577  Bool_t rel = IsRelative();
6578  Double_t fac = (rel ? 100 : 1);
6579  const char* pst = (rel ? "%" : "");
6580  gROOT->IndentLevel();
6581  Printf("%s/%s (ID # %d) -%f%s +%f%s",
6582  GetName(), GetTitle(), GetUniqueID(),
6583  fac*fEyl, pst, fac*fEyh, pst);
6584  TString opt(option);
6585  if (!opt.Contains("ATTR", TString::kIgnoreCase)) return;
6586  gROOT->IndentLevel();
6587  Printf(" [option: %3d line (c/s/w):%3d/%1d/%2d fill (c/s):%3d/%4d]",
6588  fOption, GetLineColor(), GetLineStyle(), GetLineWidth(),
6589  GetFillColor(), GetFillStyle());
6590 
6591  }
6596 
6597  ClassDef(HolderCommon,3);
6598  };
6613  Bool_t cmn,
6614  Bool_t stat,
6615  Bool_t quad,
6616  Bool_t nosqrt,
6617  Double_t& eyl,
6618  Double_t& eyh) const
6619  {
6620  Double_t wyl = 0;
6621  Double_t wyh = 0;
6622  return GetYandError(i, cmn, stat, quad, nosqrt, eyl, eyh, wyl, wyh);
6623  }
6644  Bool_t cmn,
6645  Bool_t stat,
6646  Bool_t quad,
6647  Bool_t nosqrt,
6648  Double_t& eyl,
6649  Double_t& eyh,
6650  Double_t& wyl,
6651  Double_t& wyh) const
6652  {
6653  // --- Find location for common errors ---------------------------
6654  Int_t n = fData->GetN();
6655  if (i >= n) {
6656  eyl = -1;
6657  eyh = -1;
6658  return 0;
6659  }
6660  Double_t y = fData->GetY()[i];
6661  wyl = GetStatErrorDown(i);
6662  wyh = GetStatErrorUp(i);
6663  eyl = (stat ? wyl : 0);
6664  eyh = (stat ? wyh : 0);
6665  if (quad) {
6666  wyl *= wyl;
6667  wyh *= wyh;
6668  eyl *= eyl;
6669  eyh *= eyh;
6670  }
6671  Double_t exl = GetErrorXLeft(i);
6672  Double_t exh = GetErrorXRight(i);
6673 
6674  // Otherwise, we are adding all selected errors togethter
6675  // Graph* g = static_cast<Graph*>(fData->Clone("error"));
6676  // g->SetTitle(fSumTitle);
6677  Int_t xMode = -1;
6678  if (cmn) {
6679  TIter nextC(&fCommon);
6680  HolderCommon* hc = 0;
6681  while ((hc = static_cast<HolderCommon*>(nextC()))) {
6682  if (hc->TestBit(kUsedBit)) continue;
6683  if (xMode < 0) xMode = hc->XMode(fSumOption);
6684  hc->SumPointError(y, xMode, false, quad, exl, exh, wyl, wyh);
6685  if (hc->TestBit(kOnlyWeightBit)) continue;
6686  hc->SumPointError(y, xMode, false, quad, exl, exh, eyl, eyh);
6687  }
6688  }
6689  TIter nextP(&fPoint2Point);
6690  HolderP2P* hp = 0;
6691  while ((hp = static_cast<HolderP2P*>(nextP()))) {
6692  if (hp->TestBit(kUsedBit)) continue;
6693  if (xMode < 0) xMode = hp->XMode(fSumOption);
6694  hp->SumPointError(i, xMode, false, quad, exl, exh, wyl, wyh);
6695  if (hp->TestBit(kOnlyWeightBit)) continue;
6696  hp->SumPointError(i, xMode, false, quad, exl, exh, eyl, eyh);
6697  }
6698  if (quad && !nosqrt) {
6699  eyl = TMath::Sqrt(eyl);
6700  eyh = TMath::Sqrt(eyh);
6701  wyl = TMath::Sqrt(wyl);
6702  wyh = TMath::Sqrt(wyh);
6703  }
6704  return y;
6705  }
6749  static Double_t Round(Double_t v, Int_t p, Int_t& rexpo)
6750  {
6751  if (v == 0) {
6752  // ret.Form("%.*f", p-1, 0);
6753  rexpo = 0;
6754  return 0; // ret.Data();
6755  }
6756 
6757  // Get the sign, take absolute value, get exponent, get scalar, scale
6758  Bool_t neg = (v < 0);
6759  Double_t tmp = TMath::Abs(v);
6760  Int_t expo = Int_t(TMath::Log10(tmp));
6761  Double_t tens = TMath::Power(10, expo - p + 1);
6762  Double_t n = RoundN(tens, tmp); // TMath::Ceil(tmp/tens);
6763 
6764  if (n < TMath::Power(10, p-1)) {
6765  // If scaled is less than 10 to our precision, take off one digit
6766  expo--;
6767  tens = TMath::Power(10, expo - p + 1);
6768  n = RoundN(tens,tmp); // TMath::Ceil(tmp/tens);
6769  }
6770 
6771  if (TMath::Abs((n+1) * tens) <= TMath::Abs(n * tens))
6772  // If the value is not near original, add one
6773  n++;
6774 
6775  if (n >= TMath::Power(10,p)) {
6776  // If value is larger than 10 to our precision, scale down one
6777  n /= 10;
6778  expo++;
6779  }
6780 
6781  rexpo = expo-p+1;
6782  return (neg ? -1 : 1)*n;
6783  }
6784 protected:
6785  static void SwapPoints(Graph* g, Int_t i, Int_t j, Bool_t reflect)
6786  {
6787  Double_t s = (reflect ? -1 : 1);
6788  Double_t tmpX = g->GetX()[i];
6789  Double_t tmpY = g->GetY()[i];
6790  Double_t tmpExl = g->GetErrorXlow(i);
6791  Double_t tmpExh = g->GetErrorXhigh(i);
6792  Double_t tmpEyl = g->GetErrorYlow(i);
6793  Double_t tmpEyh = g->GetErrorYhigh(i);
6794  if (reflect) {
6795  // Swap low and high X error when reflecting
6796  Double_t tmp = tmpExl;
6797  tmpExl = tmpExh;
6798  tmpExh = tmp;
6799  }
6800  g->SetPoint(i, s*g->GetX()[j], g->GetY()[j]);
6801  g->SetPointError(i,
6802  g->GetErrorXlow(j), g->GetErrorXhigh(j),
6803  g->GetErrorYlow(j), g->GetErrorYhigh(j));
6804  g->SetPoint(j, s*tmpX, tmpY);
6805  g->SetPointError(j, tmpExl, tmpExh, tmpEyl, tmpEyh);
6806  }
6811  static const char* FormatKey(const char* key)
6812  {
6813  TString tmp(Form("*%s:", key));
6814  return Form("%-12s", tmp.Data());
6815  }
6824  static const char* ExtractField(const TString& value, Int_t idx)
6825  {
6826  static TString val;
6827  val = "";
6828  TObjArray* tokens = value.Tokenize(":");
6829  Int_t iVal = TMath::Min(Int_t(idx),tokens->GetEntriesFast()-1);
6830  if (iVal < 0) return val;
6831  TString tmp = tokens->At(iVal)->GetName();
6832  val = tmp.Strip(TString::kBoth, ' ');
6833  // tokens->ls();
6834  // ::Info("", "Selecting %d: %s", iVal, val.Data());
6835  delete tokens;
6836  return val.Data();
6837  }
6846  static void StoreQual(TList& quals, Int_t idx,
6847  const char* name, const char* val)
6848  {
6849  TObject* oqv = quals.FindObject(name);
6850  TList* qv = 0;
6851  if (!oqv) {
6852  // ::Info("StoreQual", "Creating list for qualifier %s", name);
6853  qv = new TList;
6854  qv->SetOwner();
6855  qv->SetName(name);
6856  quals.Add(qv);
6857  }
6858  else qv = static_cast<TList*>(oqv);
6859  // ::Info("StoreQual", "Storing %d value %s=%s", idx, name, val);
6860  qv->AddAt(new TObjString(val), idx);
6861  }
6869  static void StoreQual(TList& quals, Int_t idx, TObject* q)
6870  {
6871  StoreQual(quals, idx, q->GetName(), q->GetTitle());
6872  }
6883  std::ostream& ExportKey(std::ostream& out,
6884  const char* which,
6885  Bool_t alsoKey = true,
6886  const char* fill = "<please fill in>") const
6887  {
6888  Bool_t gotit = false;
6889  if (fMap) {
6890  TIter nextK(fMap);
6891  TObject* obj = 0;
6892  while ((obj = nextK())) {
6893  TString k(obj->GetName());
6894  if (!k.EqualTo(which)) continue;
6895  gotit = true;
6896  if (alsoKey) out << FormatKey(obj->GetName());
6897  out << obj->GetTitle() << std::endl;
6898  }
6899  }
6900  if (!gotit) {
6901  if (alsoKey)
6902  out << FormatKey(which);
6903  out << fill << std::endl;
6904  }
6905  return out;
6906  }
6917  std::ostream& ExportHeader(std::ostream& out,
6918  Bool_t alsoTop=false,
6919  Bool_t alsoComment=false,
6920  Int_t nsign=-1) const
6921  {
6922  if (alsoComment) {
6923  out << "# Generated by GraphSysErr\n"
6924  << "# See also\n"
6925  << "# \n"
6926  << "# http://hepdata.cedar.ac.uk/resource/sample.input\n"
6927  << "# http://hepdata.cedar.ac.uk/submittingdata\n"
6928  << "# \n"
6929  << "# Please fill in missing fields\n";
6930  }
6931  if (alsoTop) {
6932  TDatime now;
6933  UserGroup_t* u = gSystem->GetUserInfo();
6934  TString ref = GetKey("reference");
6935  const char* lab = GetKey("laboratory");
6936  const char* accel = GetKey("accelerator");
6937  const char* exper = GetKey("detector");
6938  const char* abs = GetKey("abstract");
6939  const char* months[] = {"JAN","FEB","MAR","APR",
6940  "MAY","JUN","JUL","AUG",
6941  "SEP","OCT","NOV","DEC"};
6942  if (ref.IsNull())
6943  out << "*reference: <journal/archive ref> : <year>" << std::endl;
6944  else {
6945  ExportKey(out, "reference", true, "<reference>");
6946  }
6947  ExportKey(out, "author", true, "<Last name of first author in CAPS>");
6948  ExportKey(out, "doi", true, "<DOI number>");
6949  ExportKey(out, "inspireId", true, "<inSpire ID>");
6950  ExportKey(out, "cdsId", true, "<CDS ID>");
6951  ExportKey(out, "durhamId", true, "<fill in on submission>");
6952  ExportKey(out, "title", true, "<Title of paper>");
6953 
6954  out << FormatKey("status") << "Encoded "
6955  << std::setw(2) << now.GetDay() <<' '<< months[now.GetMonth()-1]<<' '
6956  << std::setw(4) << now.GetYear() << " by "<< u->fRealName.Data()<<"\n"
6957  << FormatKey("experiment") << (lab ? lab : "<lab>") << '-'
6958  << (accel ? accel : "<accelerator>") << '-'
6959  << (exper ? exper : "<experiment>") << '\n'
6960  << FormatKey("detector") << (exper ? exper : "<experiment>") << '\n'
6961  << FormatKey("comment") << (lab ? lab : "<lab>") << '-'
6962  << (accel ? accel : "<accelerator>") << ". "
6963  << (abs ? abs : "<abstract>") << "\n"
6964  << std::endl;
6965  }
6966  out << "# Start of dataset\n";
6967  if (nsign >= 0) out << "# Using " << nsign << " significant digits\n";
6968  out << FormatKey("dataset") << std::endl;
6969  // out << FormatKey("dscomment") << GetTitle() << std::endl;
6970  const char* fields[] = { "location",
6971  "reackey",
6972  "obskey",
6973  "dscomment",
6974  // "qual",
6975  0 };
6976  const char** pfld = fields;
6977  while (*pfld) {
6978  ExportKey(out, *pfld, true);
6979  pfld++;
6980  }
6981  return out;
6982  }
6995  static Int_t ExportError(std::ostream& o,
6996  Double_t low,
6997  Double_t high,
6998  Bool_t nopm,
6999  Bool_t rel,
7000  Int_t nsign)
7001  {
7002  Double_t l = low;
7003  Double_t h = high;
7004  Int_t p = o.precision();
7005  Int_t e = 0;
7006  if (nsign >= 0) {
7007  Int_t pp = nsign;
7008  Int_t le, he;
7009  Double_t lm = Round(l, nsign, le);
7010  Double_t hm = Round(h, nsign, he);
7011  e = TMath::Min(le,he);
7012  l = lm*TMath::Power(10,le);
7013  h = hm*TMath::Power(10,he);
7014  o.precision(pp);
7015  }
7016  // ::Info("ExportError","Before (%2d) -%8f.+%8f After -%8f,+%8f",
7017  // nsign, low,high,l,h);
7018  if (TMath::Abs(h-l) < 1e-12)
7019  o << (nopm ? "" : "+- ") << TMath::Abs(l) << (rel ? " PCT" : "");
7020  else
7021  o << "+" << h << ",-" << l << (rel ? " PCT" : "");
7022  o.precision(p);
7023  return e;
7024  }
7036  std::ostream& ExportPoint(std::ostream& out,
7037  Int_t i,
7038  Bool_t alsoX=true,
7039  Bool_t sysName=true,
7040  Int_t nsign=0) const
7041  {
7042  Int_t p = out.precision();
7043  // out.setf(std::ios::fixed);
7044  if (alsoX) {
7045  Double_t x = GetX(i);
7046  Double_t exl = fData->GetErrorXlow(i);
7047  Double_t exh = fData->GetErrorXhigh(i);
7048  Int_t pp = p;
7049  Int_t ppp = p;
7050  if (nsign > 0) {
7051  Int_t pxl, pxh;
7052  Double_t mxl = Round(exl,nsign,pxl);
7053  Double_t mxh = Round(exh,nsign,pxh);
7054  exl = mxl*TMath::Power(10, pxl);
7055  exh = mxh*TMath::Power(10, pxh);
7056  Int_t mp = TMath::Min(pxl,pxh);
7057  pp = nsign;
7058  ppp = TMath::Ceil(TMath::Log10(TMath::Abs(x)))-mp;
7059  // Info("","pxl=%d pxh=%d mp=%d ppp=%d", pxl, pxh, mp, ppp);
7060  }
7061 
7062  if (TMath::Abs(exl) < 1e-10 && TMath::Abs(exh) < 1e-10) {
7063  if (nsign <= 0) out << ' ' << x;
7064  else {
7065  Int_t px;
7066  Double_t mx = Round(x, nsign, px);
7067  Double_t xx = mx*TMath::Power(10,px);
7068  // Printf("round %f -> %f (%d) -> %.*g", x, mx, px, nsign, xx);
7069  out.precision(nsign);
7070  out << ' ' << xx;
7071  }
7072  }
7073  else {
7074  if (TMath::Abs(exh-exl) < 1e-10) {
7075  out.precision(ppp);
7076  out << ' ' << x-exl << " TO " << x+exh;
7077  }
7078  else {
7079  out.precision(ppp);
7080  out << ' ' << x << ' ';
7081  out.precision(pp);
7082  out << '+' << exh << ",-" << exl;
7083  }
7084  }
7085  out.precision(p);
7086  out << "; ";
7087  }
7088  // out.unsetf(std::ios::fixed);
7089 
7090  Bool_t statRel = fStatRelative;
7091  Double_t y = GetY(i);
7092  Double_t fy = (statRel ? (y == 0 ? 0 : 100./y) : 1);
7093  Double_t eyl = GetStatErrorDown(i) * fy;
7094  Double_t eyh = GetStatErrorUp(i) * fy;
7095  // out << y << ' ';
7096  std::stringstream tmp;
7097  // tmp.setf(std::ios::fixed);
7098  Int_t le = ExportError(tmp, eyl, eyh, false, statRel, nsign);
7099 
7100  if (fPoint2Point.GetEntries() > 0) {
7101  tmp << " (" << std::flush;
7102  TIter nextP(&fPoint2Point);
7103  HolderP2P* holderP2P = 0;
7104  Bool_t first = true;
7105  while ((holderP2P = static_cast<HolderP2P*>(nextP()))) {
7106  Bool_t rel = holderP2P->IsRelative();
7107  fy = (rel ? (y == 0 ? 0 : 100./y) : 1);
7108  if (!first) tmp << ',';
7109  Double_t esh = holderP2P->GetYUp(i)*fy;
7110  Double_t esl = holderP2P->GetYDown(i)*fy;
7111  tmp << "DSYS=";
7112  le = TMath::Min(le,ExportError(tmp, esl, esh, true, rel, nsign));
7113  if (sysName)
7114  tmp << ':' << holderP2P->GetTitle() << std::flush;
7115  first = false;
7116  }
7117  tmp << ')';
7118  }
7119  if (nsign) {
7120  Int_t nY = TMath::Ceil(TMath::Log10(TMath::Abs(y)))-le;
7121  out.precision(nY);
7122  out << y << ' ';
7123  out.precision(p);
7124  }
7125  else
7126  out << y << ' ';
7127  out << tmp.str() << ';';
7128  return out;
7129  }
7148  static Double_t RoundN(Double_t tens, Double_t tmp)
7149  {
7150  // Add small number to not loose least digit if a 1
7151  Double_t n = TMath::Floor(100*tmp/tens+.00001);
7152  Int_t next = int(n/10) % 10;
7153  if (next > 5)
7154  // Round up
7155  return TMath::Ceil(n/100);
7156  if (next < 5)
7157  // Round down
7158  return TMath::Floor(n/100);
7159  Int_t nnext = int(n) % 10;
7160  if (nnext != 0)
7161  // Round up
7162  return TMath::Ceil(n/100);
7163  Int_t last = int(n/100) % 10;
7164  if (last % 2 == 0)
7165  // Round down
7166  return TMath::Floor(n/100);
7167  // Round to nearest even
7168  return TMath::Ceil(n/100);
7169  }
7179  static TString& Token(TObjArray* c, UShort_t idx, Bool_t verbose=true)
7180  {
7181  static TString empty("");
7182  if (idx >= c->GetEntriesFast()) {
7183  if (verbose) {
7184  ::Error("Token", "Token %d does not exist in array", idx);
7185  c->ls();
7186  }
7187  return empty;
7188  }
7189  return static_cast<TObjString*>(c->At(idx))->String();
7190  }
7201  static Bool_t ImportError(const TString& s,
7202  Double_t& el, Double_t& eh, Bool_t& rel)
7203  {
7204  TString tmp(s);
7205  // ::Info("ImportError", "Parsing '%s' for errors", s.Data());
7206  // Allow for PCT (any case) and % as relative markers
7207  if (tmp.Contains("PCT", TString::kIgnoreCase) || tmp.Contains("%"))
7208  rel = true;
7209  if (tmp.BeginsWith("+-") || tmp.BeginsWith("-+")) {
7210  el = eh = tmp.Atof();
7211  // ::Info("ImportError", "Got symmetric %f,%f", el, eh);
7212  return true;
7213  }
7214  TObjArray* tokens = tmp.Tokenize(",");
7215  for (Int_t i = 0; i < tokens->GetEntriesFast(); i++) {
7216  TString& t = Token(tokens, i);
7217  if (t.IsNull()) continue;
7218  t.Remove(TString::kBoth, ' ');
7219  Double_t v = t.Atof();
7220  // ::Info("ImportError", "Got token '%s' -> %f", t.Data(), v);
7221  if (t[0] != '-' && t[0] != '+') {
7222  // ::Info("ImportError", "First char not + or -: %c", t.Data()[0]);
7223  el = eh = TMath::Abs(v);
7224  }
7225  else if (v < 0) {
7226  // ::Info("ImportError", "Value %f < 0 -> setting el=%f", v, -v);
7227  if (el > 0) {
7228  Double_t max = TMath::Max(el, -v);
7229  ::Warning("ImportError",
7230  "Lower error already set to %f, chosing the maximum of "
7231  "(%f,%f) -> %f", el, el, -v, max);
7232  v = -max;
7233  }
7234  el = -v;
7235  }
7236  else {
7237  // ::Info("ImportError", "Value %f >= 0 -> setting eh=%f", v, v);
7238  if (eh > 0) {
7239  Double_t max = TMath::Max(eh, v);
7240  ::Warning("ImportError",
7241  "Lower error already set to %f, chosing the maximum of "
7242  "(%f,%f) -> %f", eh, eh, v, max);
7243  v = max;
7244  }
7245  eh = v;
7246  }
7247  }
7248  // ::Info("ImportError", "Got a-symmetric %f,%f", el, eh);
7249  tokens->Delete();
7250  delete tokens;
7251  return true;
7252  }
7253 
7265  static Bool_t ImportPoint(const TString& s,
7266  Double_t& v,
7267  Double_t& el,
7268  Double_t& eh,
7269  Bool_t& rel)
7270  {
7271  // ::Info("ImportPoint", "s=%s", s.Data());
7272  rel = s.Contains("PCT");
7273  TObjArray* tok = s.Tokenize(" ,");
7274  if (tok->GetEntriesFast() >= 4 &&
7275  Token(tok,2).EqualTo("TO", TString::kIgnoreCase) &&
7276  Token(tok,1).Contains("BIN", TString::kIgnoreCase)) {
7277  v = Token(tok,0).Atof();
7278  TString tmp = Token(tok,1);
7279  TString s1 = tmp.Strip(TString::kBoth, '(');
7280  s1.ReplaceAll("BIN=","");
7281  Double_t v1 = s1.Atof();
7282  tmp = Token(tok,3);
7283  s1 = tmp.Strip(TString::kBoth, ')');
7284  Double_t v2 = s1.Atof();
7285  el = (v-v1);
7286  eh = (v2-v);
7287  }
7288  else if (tok->GetEntriesFast() >= 3 &&
7289  Token(tok,1).EqualTo("TO", TString::kIgnoreCase)) {
7290  Double_t v1 = Token(tok,0).Atof();
7291  Double_t v2 = Token(tok,2).Atof();
7292  v = (v1+v2)/2;
7293  el = (v-v1);
7294  eh = (v2-v);
7295  }
7296  else {
7297  v = Token(tok,0).Atof();
7298  if (tok->GetEntriesFast() >= 2 &&
7299  (Token(tok,1).BeginsWith("+-") ||
7300  Token(tok,1).BeginsWith("-+"))) {
7301  TString t = Token(tok,1);
7302  if (t.Length() > 2) {
7303  t.Remove(0,2);
7304  // ::Info("ImportPoint", "Extract from %s", t.Data());
7305  el = eh = t.Atof();
7306  }
7307  else {
7308  // ::Info("ImportPoint","Extract from next %s",Token(tok, 2).Data());
7309  el = eh = Token(tok,2).Atof();
7310  }
7311  }
7312  else {
7313  for (Int_t i = 1; i < tok->GetEntriesFast(); i++) {
7314  TString t = Token(tok, i);
7315  Char_t m = t[0];
7316  // ::Info("ImportPoint", "%s -> %c", t.Data(), m);
7317  if (m != '+' && m != '-') continue;
7318  if (t.Length() == 1) t = Token(tok, ++i); // Take next
7319  else t.Remove(0,1); // Remove sign
7320  // ::Info("ImportPoint", "%s -> %f", t.Data(), t.Atof());
7321  if (m == '-') el = t.Atof();
7322  else eh = t.Atof();
7323  }
7324  }
7325  }
7326  tok->Delete();
7327  delete tok;
7328  return true;
7329  }
7330  /* @} */
7341  void SqrtPoint(Graph* g, Int_t i)
7342  {
7343  if (kVerbose & kDraw)
7344  printf("%20s Point %3d %f,%f -> ","",
7345  i,g->GetEYlow()[i],g->GetEYhigh()[i]);
7346 
7347  g->GetEYlow()[i] = TMath::Sqrt(g->GetEYlow()[i]);
7348  g->GetEYhigh()[i] = TMath::Sqrt(g->GetEYhigh()[i]);
7349 
7350  if (kVerbose & kDraw)
7351  Printf("%f,%f", g->GetEYlow()[i], g->GetEYhigh()[i]);
7352  }
7358  void SqrtGraph(Graph* g)
7359  {
7360  for (Int_t i = 0; i < g->GetN(); i++)
7361  SqrtPoint(g, i);
7362  }
7363  /* @} */
7373  const char* FormatOption(UInt_t opt)
7374  {
7375  static TString ret;
7376  ret = "";
7377  switch (opt) {
7378  case kNormal: ret = ""; break; // Line with ticks
7379  case kNoTick: ret = "Z0"; break; // Line with no ticks
7380  case kArrow: ret = ">0"; break; // Linw with arrows
7381  case kRect: ret = "20"; break; // Rectangle w/fill
7382  case kBox: ret = "50"; break; // Rectangle w/fill & line
7383  case kFill: ret = "30"; break; // Filled area
7384  case kCurve: ret = "40"; break; // Filled smoothed area
7385  case kHat: ret = "[]0"; break; // Hats
7386  case kBar: ret = "||0"; break; // A bar
7387  case kNone: ret = "XP"; break; // No errors
7388  case kLine: ret = "CX"; break;
7389  case kConnect: ret = "LX"; break;
7390  }
7391  return ret.Data();
7392  }
7399  {
7400  fData = new Graph(n);
7401  fData->SetName(Form("%s_data", fName.Data()));
7402  fData->SetTitle(GetTitle());
7403  }
7414  Int_t FindPoint(Double_t x, Int_t& i1, Int_t& i2, Double_t fac=1) const
7415  {
7416  i1 = -1;
7417  i2 = -1;
7418  Int_t n = GetN();
7419  if (x < GetX(0)-fac*GetErrorXLeft(0) ||
7420  x > GetX(n-1)+fac*GetErrorXRight(n-1)) {
7421  // If we're out side the coverage, set return to -1
7422  // Info("FindPoint", "X=%f Out of range [%f,%f]",
7423  // x,GetX(0),GetX(GetN()-1));
7424  i1 = i2 = -1;
7425  return -2;
7426  }
7427  const Double_t tol = 1e-9;
7428  for (Int_t i=0; i < n; i++) {
7429  // if (TMath::Abs(GetX(i)-x) < tol) {
7430  if (x <= GetX(i)+fac*GetErrorXRight(i) &&
7431  x >= GetX(i)-fac*GetErrorXLeft(i)) {
7432  // if ((GetX(i)+GetErrorXRight(i)-x)<tol &&
7433  // (x - GetX(i) - GetErrorXLeft(i) <tol)) {
7434  // Found matching point
7435  i2 = i1 = i;
7436  break;
7437  }
7438  // Info("FindPoint", "Not @ %d (%f not in %f -%f +%f)",
7439  // i, x, GetX(i), GetErrorXLeft(i), GetErrorXRight(i));
7440  if (x > GetX(i)) {
7441  // X still larger
7442  i1 = i;
7443  }
7444  else {
7445  // X is lower
7446  i2 = i;
7447  break;
7448  }
7449  }
7450  if (i1 == i2)
7451  // If equal, return the point
7452  return i1;
7453  if (TMath::Abs((GetX(i1)+fac*GetErrorXRight(i1)) -
7454  (GetX(i2)-fac*GetErrorXLeft(i2))) > 10*tol) {
7455  // If the two found points are not adjecent, we are in a hole
7456  // and we indicate that.
7457  // Info("FindPoint","Found hole at i1=%d(%f+%f) i2=%d(%f-%f) (%f)",
7458  // i1, GetX(i1), GetErrorXRight(i1),
7459  // i2, GetX(i2), GetErrorXLeft(i2), x);
7460  // Print("xy");
7461  i1 = i2 = -1;
7462  return -2;
7463  }
7464  if (i2 > 0) {
7465  // Otherwise return negative value to indicate we should look at
7466  // two points.
7467  // Info("FindPoint","Found i1=%d i2=%d %f in [%f,%f]", i1, i2,
7468  // x, GetX(i1)-GetErrorXLeft(i1), GetX(i2)+GetErrorXRight(i2));
7469  return -1;
7470  }
7471  // Return negative two to indicate nothing was found
7472  return -2;
7473 
7474  }
7482  TMultiGraph* MakeMulti(Option_t* option)
7483  {
7484  if (!fData) {
7485  Warning("MakeMulti", "No data defined");
7486  return 0;
7487  }
7488 
7489  // --- Process options -------------------------------------------
7490  TString opt(option);
7491  opt.ToUpper();
7492 
7493  Bool_t cmn = opt.Contains("COMMON");
7494  Bool_t combine = opt.Contains("COMBINED");
7495  Bool_t stat = opt.Contains("STAT");
7496  Bool_t quad = opt.Contains("QUAD");
7497  combine = !opt.Contains("STACK");
7498  quad = !opt.Contains("DIRECT");
7499  Bool_t split = opt.Contains("SPLIT");
7500  Int_t xpos = (opt.Contains("WEST") ? -1 : 1);
7501  Int_t ypos = (opt.Contains("MIN") ? -1 :
7502  opt.Contains("MAX") ? 1 : 0);
7503 
7504  // --- Some general information ----------------------------------
7505  Int_t n = fData->GetN();
7506  Double_t dx = TMath::Max(gStyle->GetErrorX(),0.0001F);
7507  if (n <= 0) {
7508  Warning("MakeMulti", "Nothing to do (n=%d)", n);
7509  return 0;
7510  }
7511 
7512  // --- Find location for common errors ---------------------------
7513  Double_t xBase = 0;
7514  if (opt.Contains("XBASE=")) {
7515  Int_t idx = opt.Index("XBASE=")+6;
7516  TString tmp = opt(idx,opt.Length()-idx);
7517  xBase = tmp.Atof();
7518  }
7519  else {
7520  xBase = ((xpos < 0 ?
7521  fData->GetX()[0] - fData->GetEXlow()[0]:
7522  fData->GetX()[n-1] + fData->GetEXhigh()[n-1])
7523  + xpos * 1.2 * dx);
7524  }
7525  Double_t yBase = fData->GetY()[0];
7526  for (Int_t i = 0; i < n; i++) {
7527  Double_t y = fData->GetY()[i];
7528  yBase = (ypos < 0 ? TMath::Min(yBase, y) :
7529  ypos > 0 ? TMath::Max(yBase, y) : yBase + y);
7530  }
7531  if (ypos == 0) yBase /= n;
7532 
7533  // --- Create stack and temp collection --------------------------
7534  TList drawn;
7535  drawn.SetOwner(false);
7536 
7537  Graph* prev = 0;
7538  if (!combine) {
7539  // If we're not combining, we basically copy the graphs to the
7540  // draw list, but modify the errors to show progresive larger
7541  // errors.
7542  prev = fData;
7543  if (cmn) {
7544  // If we should add common errors, do that first
7545  TIter next(&fCommon);
7546  Holder* h = 0;
7547  while ((h = static_cast<Holder*>(next()))) {
7548  Graph* g = h->StackError(prev, (!stat ? prev == fData : false), quad);
7549  if (!g) continue;
7550 
7551  if (quad) SqrtGraph(g);
7552  prev = g;
7553  drawn.AddFirst(g, FormatOption(h->GetDOption()));
7554  }
7555  }
7556  TIter next(&fPoint2Point);
7557  Holder* h = 0;
7558  while ((h = static_cast<Holder*>(next()))) {
7559  Graph* g = h->StackError(prev, (!stat ? prev == fData : false), quad);
7560  if (!g) continue;
7561 
7562  if (quad) SqrtGraph(g);
7563  prev = g;
7564  drawn.AddFirst(g, FormatOption(h->GetDOption()));
7565  }
7566  }
7567  else {
7568  // Otherwise, we are adding all selected errors togethter
7569  Graph* g = static_cast<Graph*>(fData->Clone("error"));
7570  g->SetTitle(fSumTitle);
7571  Holder* h = 0;
7572  for (Int_t i = 0; i < n; i++) {
7573  if (!stat) {
7574  g->SetPointError(i, g->GetEXlow()[i],
7575  g->GetEXhigh()[i],0,0);
7576  }
7577  else if (quad) {
7578  // g->GetEXlow()[i] *= g->GetEXlow()[i];
7579  g->GetEYlow()[i] *= g->GetEYlow()[i];
7580  // g->GetEXhigh()[i] *= g->GetEXhigh()[i];
7581  g->GetEYhigh()[i] *= g->GetEYhigh()[i];
7582  }
7583 
7584 
7585  if (cmn) {
7586  TIter nextC(&fCommon);
7587  while ((h = static_cast<Holder*>(nextC()))) {
7588  h->SumError(g, i, false, quad, fSumOption);
7589 
7590  }
7591  }
7592  TIter nextP(&fPoint2Point);
7593  while ((h = static_cast<Holder*>(nextP())))
7594  h->SumError(g, i, false, quad, fSumOption);
7595 
7596  // Printf("Point # %d", i);
7597  if (quad) SqrtPoint(g, i);
7598  }
7599  fSumLine.Copy(*g);
7600  fSumFill.Copy(*g);
7601  g->SetMarkerStyle(0);
7602  g->SetMarkerSize(0);
7603  g->SetMarkerColor(0);
7604 
7605  drawn.AddFirst(g, FormatOption(fSumOption));
7606  }
7607  // --- Show common errors on the side, if requested --------------
7608  if (!cmn) {
7609  TIter nextC(&fCommon);
7610  HolderCommon* hc = 0;
7611  if (!combine) {
7612  prev = 0;
7613  while ((hc = static_cast<HolderCommon*>(nextC()))) {
7614  Graph* g = hc->BarError(prev, quad && !split, xBase, yBase);
7615  if (!g) {
7616  Warning("MakeMulti", "Got no graph for common error %s",
7617  hc->GetTitle());
7618  continue;
7619  }
7620 
7621  if (quad && !split) SqrtGraph(g);
7622 
7623  if (split) xBase += xpos * 1.2 * dx;
7624  if (!prev)
7625  drawn.AddLast(g, FormatOption(hc->GetDOption()));
7626  else {
7627  // Here, we use the underling linked list of the list so
7628  // that we can add an object before another including the
7629  // possible options
7630  TObjLink* f = drawn.FirstLink();
7631  TObjLink* p = 0;
7632  while (f && f->GetObject()) {
7633