AliPhysics  d497547 (d497547)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliTrackletdNdetaUtils.C
Go to the documentation of this file.
1 #ifndef ALITRACKLETDNDETAUTILS_H
2 #define ALITRACKLETDNDETAUTILS_H
3 #ifndef __CINT__
4 # include <TH2.h>
5 # include <TH1.h>
6 # include <TList.h>
7 # include <TParameter.h>
8 # include <TError.h>
9 # include <TMath.h>
10 # include <TDirectory.h>
11 # include <THashList.h>
12 #else
13 class TList;
14 class TH1;
15 class TH2;
16 class TAxis;
17 class TDirectory;
18 #endif
19 
25 {
26 public:
27  typedef TList Container;
36  //_________________________________________________________________________
50  static Bool_t CheckAxisNBins(const char* which,
51  const TAxis* a1,
52  const TAxis* a2);
62  static Bool_t CheckAxisLimits(const char* which,
63  const TAxis* a1,
64  const TAxis* a2);
75  static Bool_t CheckAxisBins(const char* which,
76  const TAxis* a1,
77  const TAxis* a2);
87  static Bool_t CheckAxisLabels(const char* which,
88  const TAxis* a1,
89  const TAxis* a2);
100  static Bool_t CheckAxis(const char* which,
101  const TAxis* a1,
102  const TAxis* a2,
103  Bool_t alsoLbls);
113  static Bool_t CheckConsistency(const TH1* h1, const TH1* h2);
114  /* @} */
115  //__________________________________________________________________
129  static TObject* GetO(Container* parent,const char* name,TClass* cls);
139  static TObject* GetO(TDirectory* parent,const char* name,TClass* cls);
148  static TH1* GetH1(Container* parent, const char* name);
157  static TH2* GetH2(Container* parent, const char* name);
166  static Container* GetC(Container* parent, const char* name);
175  static Container* GetC(TDirectory* parent, const char* name);
185  static Double_t GetD(Container* parent, const char* name, Double_t def=-1);
195  static Int_t GetI(Container* parent, const char* name, Int_t def=-1);
205  static Int_t GetB(Container* parent, const char* name, Bool_t def=false);
215  static TH1* CopyH1(Container* parent, const char* name,const char* newName=0);
225  static TH2* CopyH2(Container* parent, const char* name,const char* newName=0);
238  static TH1* Make1D(Container& c,
239  const TString& name,
240  const TString& title,
241  Color_t color,
242  Style_t style,
243  const TAxis& xAxis);
257  static TH2* Make2D(Container& c,
258  const TString& name,
259  const TString& title,
260  Color_t color,
261  Style_t style,
262  const TAxis& xAxis,
263  const TAxis& yAxis);
264 
274  static TH1* Scale(TH1* h, Double_t x, Double_t xe);
284  static TH2* Scale(TH2* h, Double_t x, Double_t xe);
294  static TH2* Scale(TH2* h, TH1* s);
301  static void FixAxis(TAxis& axis, const char* title=0);
309  static void ScaleAxis(TAxis& axis, Double_t fact=1);
317  static void SetAxis(TAxis& axis, Int_t n, Double_t* borders);
326  static void SetAxis(TAxis& axis, const TString& spec, const char* sep=":,");
335  static void SetAxis(TAxis& axis, Int_t n, Double_t l, Double_t h);
343  static void SetAxis(TAxis& axis, Int_t n, Double_t m);
351  static void PrintAxis(const TAxis& axis, Int_t nSig=2, const char* alt=0);
362  static TH2* ScaleToIPz(TH2* h, TH1* ipZ, Bool_t full=false);
383  static TH1* AverageOverIPz(TH2* h,
384  const char* name,
385  UShort_t mode,
386  TH1* ipz,
387  TH2* mask=0);
396  static TObject* CloneAndAdd(Container* c, TObject* o);
409  static Double_t Integrate(TH1* h, Double_t min, Double_t max, Double_t& err);
421  static Double_t RatioE(Double_t n, Double_t en,
422  Double_t d, Double_t ed,
423  Double_t& er);
435  static Double_t RatioE2(Double_t n, Double_t e2n,
436  Double_t d, Double_t e2d,
437  Double_t& e2r);
438  /* @} */
439 protected:
440  ClassDef(AliTrackletdNdetaUtils,1); // Utilities
441 };
442 
443 //====================================================================
444 // Utilities
445 //____________________________________________________________________
447  const TAxis* a1,
448  const TAxis* a2)
449 {
450  if (a1->GetNbins() != a2->GetNbins()) {
451  ::Warning("CheckAxisNBins", "Incompatible number %s bins: %d vs %d",
452  which, a1->GetNbins(), a2->GetNbins());
453  return false;
454  }
455  return true;
456 }
457 //____________________________________________________________________
459  const TAxis* a1,
460  const TAxis* a2)
461 {
462  if (!TMath::AreEqualRel(a1->GetXmin(), a2->GetXmin(),1.E-12) ||
463  !TMath::AreEqualRel(a1->GetXmax(), a2->GetXmax(),1.E-12)) {
464  Warning("CheckAxisLimits",
465  "Limits of %s axis incompatible [%f,%f] vs [%f,%f]", which,
466  a1->GetXmin(), a1->GetXmax(), a2->GetXmin(), a2->GetXmax());
467  return false;
468  }
469  return true;
470 }
471 //____________________________________________________________________
473  const TAxis* a1,
474  const TAxis* a2)
475 {
476  const TArrayD * h1Array = a1->GetXbins();
477  const TArrayD * h2Array = a2->GetXbins();
478  Int_t fN = h1Array->fN;
479  if ( fN == 0 ) return true;
480  if (h2Array->fN != fN) {
481  // Redundant?
482  Warning("CheckAxisBins", "Not equal number of %s bin limits: %d vs %d",
483  which, fN, h2Array->fN);
484  return false;
485  }
486  else {
487  for (int i = 0; i < fN; ++i) {
488  if (!TMath::AreEqualRel(h1Array->GetAt(i),h2Array->GetAt(i),1E-10)) {
489  Warning("CheckAxisBins",
490  "%s limit # %3d incompatible: %f vs %f",
491  which, i, h1Array->GetAt(i),h2Array->GetAt(i));
492  return false;
493  }
494  }
495  }
496  return true;
497 }
498 //____________________________________________________________________
500  const TAxis* a1,
501  const TAxis* a2)
502 {
503  // check that axis have same labels
504  THashList *l1 = (const_cast<TAxis*>(a1))->GetLabels();
505  THashList *l2 = (const_cast<TAxis*>(a2))->GetLabels();
506 
507  if (!l1 && !l2) return true;
508  if (!l1 || !l2) {
509  Warning("CheckAxisLabels", "Difference in %s labels: %p vs %p",
510  which, l1, l2);
511  return false;
512  }
513  // check now labels sizes are the same
514  if (l1->GetSize() != l2->GetSize()) {
515  Warning("CheckAxisLabels", "Different number of %s labels: %d vs %d",
516  which, l1->GetSize(), l2->GetSize());
517  return false;
518  }
519  for (int i = 1; i <= a1->GetNbins(); ++i) {
520  TString label1 = a1->GetBinLabel(i);
521  TString label2 = a2->GetBinLabel(i);
522  if (label1 != label2) {
523  Warning("CheckAxisLabels", "%s label # %d not the same: '%s' vs '%s'",
524  which, i, label1.Data(), label2.Data());
525  return false;
526  }
527  }
528 
529  return true;
530 }
531 //____________________________________________________________________
533  const TAxis* a1,
534  const TAxis* a2,
535  Bool_t alsoLbls)
536 {
537  if (!CheckAxisNBins (which, a1, a2)) return false;
538  if (!CheckAxisLimits(which, a1, a2)) return false;
539  if (!CheckAxisBins (which, a1, a2)) return false;
540  if (alsoLbls && !CheckAxisLabels(which, a1, a2)) return false;
541  return true;
542 }
543 //____________________________________________________________________
545 {
546  // Check histogram compatibility
547  if (h1 == h2) return true;
548 
549  if (h1->GetDimension() != h2->GetDimension() ) {
550  Warning("CheckConsistency",
551  "%s and %s have different dimensions %d vs %d",
552  h1->GetName(), h2->GetName(),
553  h1->GetDimension(), h2->GetDimension());
554  return false;
555  }
556  Int_t dim = h1->GetDimension();
557 
558  Bool_t alsoLbls = (h1->GetEntries() != 0 && h2->GetEntries() != 0);
559  if (!CheckAxis("X", h1->GetXaxis(), h2->GetXaxis(), alsoLbls)) return false;
560  if (dim > 1 &&
561  !CheckAxis("Y", h1->GetYaxis(), h2->GetYaxis(), alsoLbls)) return false;
562  if (dim > 2 &&
563  !CheckAxis("Z", h1->GetZaxis(), h2->GetZaxis(), alsoLbls)) return false;
564 
565  return true;
566 }
567 
568 //____________________________________________________________________
570 {
571  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
572  Double_t dc = s->GetBinContent(i);
573  Double_t de = s->GetBinError (i);
574  Double_t dr = (dc > 1e-6 ? 1/dc : 0);
575  Double_t dq = dr * de;
576  for (Int_t j = 1; j <= h->GetNbinsY(); j++) {
577  Double_t nc = h->GetBinContent(i,j);
578  Double_t ne = h->GetBinError (i,j);
579  Double_t ns = (nc > 0 ? ne/nc : 0);
580  Double_t sc = dr * nc;
581  Double_t se = sc*TMath::Sqrt(ns*ns+dq*dq);
582  // Printf("Setting bin %2d,%2d to %f +/- %f", sc, se);
583  h->SetBinContent(i,j,sc);
584  h->SetBinError (i,j,se);
585  }
586  }
587  return h;
588 }
589 
590 //____________________________________________________________________
592 {
593  Double_t rr = xe/x;
594  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
595  for (Int_t j = 1; j <= h->GetNbinsY(); j++) {
596  Double_t c = h->GetBinContent(i,j);
597  Double_t e = h->GetBinError (i,j);
598  Double_t s = (c > 0 ? e/c : 0);
599  Double_t sc = x * c;
600  Double_t se = sc*TMath::Sqrt(s*s+rr*rr);
601  // Printf("Setting bin %2d,%2d to %f +/- %f", sc, se);
602  h->SetBinContent(i,j,sc);
603  h->SetBinError (i,j,se);
604  }
605  }
606  return h;
607 }
608 
609 //____________________________________________________________________
611 {
612  Double_t rr = xe/x;
613  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
614  Double_t c = h->GetBinContent(i);
615  Double_t e = h->GetBinError (i);
616  Double_t s = (c > 0 ? e/c : 0);
617  Double_t sc = x * c;
618  Double_t se = sc*TMath::Sqrt(s*s+rr*rr);
619  // Printf("Setting bin %2d,%2d to %f +/- %f", sc, se);
620  h->SetBinContent(i,sc);
621  h->SetBinError (i,se);
622  }
623  return h;
624 }
625 
626 //____________________________________________________________________
628  const char* name,
629  TClass* cls)
630 {
631  if (!parent) {
632  ::Warning("GetO", "No parent container passed");
633  return 0;
634  }
635  TObject* o = parent->FindObject(name);
636  if (!o) {
637  ::Warning("GetO", "Object \"%s\" not found in \"%s\"",
638  name, parent->GetName());
639  // parent->ls();
640  return 0;
641  }
642  if (!cls) return o;
643  if (!o->IsA()->InheritsFrom(cls)) {
644  ::Warning("GetO", "\"%s\" is an object of class %s, not %s",
645  name, o->ClassName(), cls->GetName());
646  return 0;
647  }
648  return o;
649 }
650 //____________________________________________________________________
652  const char* name,
653  TClass* cls)
654 {
655  if (!parent) {
656  ::Warning("GetO", "No parent directory passed");
657  return 0;
658  }
659  TObject* o = parent->Get(name);
660  if (!o) {
661  ::Warning("GetO", "Object \"%s\" not found in \"%s\"",
662  name, parent->GetName());
663  parent->ls();
664  return 0;
665  }
666  if (!cls) return o;
667  if (!o->IsA()->InheritsFrom(cls)) {
668  ::Warning("GetO", "\"%s\" is an object of class %s, not %s",
669  name, o->ClassName(), cls->GetName());
670  return 0;
671  }
672  return o;
673 }
674 //____________________________________________________________________
675 TH1* AliTrackletdNdetaUtils::GetH1(Container* parent, const char* name)
676 {
677  return static_cast<TH1*>(GetO(parent, name, TH1::Class()));
678 }
679 //____________________________________________________________________
680 TH2* AliTrackletdNdetaUtils::GetH2(Container* parent, const char* name)
681 {
682  return static_cast<TH2*>(GetO(parent, name, TH2::Class()));
683 }
684 //____________________________________________________________________
686 AliTrackletdNdetaUtils::GetC(Container* parent, const char* name)
687 {
688  return static_cast<Container*>(GetO(parent, name, Container::Class()));
689 }
690 //____________________________________________________________________
692 AliTrackletdNdetaUtils::GetC(TDirectory* parent, const char* name)
693 {
694  return static_cast<Container*>(GetO(parent, name, Container::Class()));
695 }
696 //____________________________________________________________________
698  const char* name,
699  Double_t def)
700 {
701  TParameter<double>* p =
702  static_cast<TParameter<double>*>(GetO(parent, name,
704  if (!p) return def;
705  return p->GetVal();
706 }
707 //____________________________________________________________________
709  const char* name,
710  Int_t def)
711 {
712  TParameter<int>* p =
713  static_cast<TParameter<int>*>(GetO(parent, name, TParameter<int>::Class()));
714  if (!p) return def;
715  return p->GetVal();
716 }
717 //____________________________________________________________________
719  const char* name,
720  Bool_t def)
721 {
722  TParameter<bool>* p =
723  static_cast<TParameter<bool>*>(GetO(parent, name,
725  if (!p) return def;
726  return p->GetVal();
727 }
728 //____________________________________________________________________
730  const char* name,
731  const char* newName)
732 {
733  TH1* orig = GetH1(parent, name);
734  if (!orig) return 0;
735  TH1* ret = static_cast<TH1*>(orig->Clone(newName ? newName : name));
736  ret->SetDirectory(0); // Release from file container
737  return ret;
738 }
739 //____________________________________________________________________
741  const char* name,
742  const char* newName)
743 {
744  TH2* orig = GetH2(parent, name);
745  if (!orig) return 0;
746  TH2* ret = static_cast<TH2*>(orig->Clone(newName ? newName : name));
747  ret->SetDirectory(0); // Release from file container
748  return ret;
749 }
750 //____________________________________________________________________
752  const TString& name,
753  const TString& title,
754  Color_t color,
755  Style_t style,
756  const TAxis& xAxis)
757 {
758  TString n = name;
759  TString t = title;
760  TH1* ret = 0;
761  Int_t nx = xAxis.GetNbins();
762  if (t.IsNull()) t = xAxis.GetTitle();
763  if (xAxis.GetXbins() && xAxis.GetXbins()->GetArray())
764  ret = new TH1D(n,t,nx,xAxis.GetXbins()->GetArray());
765  else
766  ret = new TH1D(n,t,nx,xAxis.GetXmin(),xAxis.GetXmax());
767  ret->Sumw2();
768  ret->SetXTitle(xAxis.GetTitle());
769  static_cast<const TAttAxis&>(xAxis).Copy(*(ret->GetXaxis()));
770  ret->SetDirectory(0);
771  ret->SetLineColor(color);
772  ret->SetMarkerColor(color);
773  ret->SetFillColor(kWhite);// color);
774  ret->SetFillStyle(0);
775  ret->SetMarkerStyle(style);
776  if (const_cast<TAxis&>(xAxis).GetLabels()) {
777  for (Int_t i = 1; i <= xAxis.GetNbins(); i++)
778  ret->GetXaxis()->SetBinLabel(i, xAxis.GetBinLabel(i));
779  }
780  switch (style) {
781  case 27:
782  case 33:
783  ret->SetMarkerSize(1.4);
784  break;
785  case 29:
786  case 30:
787  ret->SetMarkerSize(1.2);
788  break;
789  }
790  c.Add(ret);
791  return ret;
792 }
793 //____________________________________________________________________
795  const TString& name,
796  const TString& title,
797  Color_t color,
798  Style_t style,
799  const TAxis& xAxis,
800  const TAxis& yAxis)
801 {
802  TString n = name;
803  TString t = title;
804  TH2* ret = 0;
805  Int_t nx = xAxis.GetNbins();
806  Int_t ny = yAxis.GetNbins();
807  const Double_t* xb = (xAxis.GetXbins() && xAxis.GetXbins()->GetArray() ?
808  xAxis.GetXbins()->GetArray() : 0);
809  const Double_t* yb = (yAxis.GetXbins() && yAxis.GetXbins()->GetArray() ?
810  yAxis.GetXbins()->GetArray() : 0);
811  if (t.IsNull())
812  t.Form("%s vs %s", yAxis.GetTitle(), xAxis.GetTitle());
813  if (xb) {
814  if (yb) ret = new TH2D(n,t,nx,xb,ny,yb);
815  else ret = new TH2D(n,t,
816  nx,xb,
817  ny,yAxis.GetXmin(),yAxis.GetXmax());
818  }
819  else {
820  if (yb) ret = new TH2D(n,t,
821  nx,xAxis.GetXmin(),xAxis.GetXmax(),
822  ny,yb);
823  else ret = new TH2D(n,t,
824  nx,xAxis.GetXmin(),xAxis.GetXmax(),
825  ny,yAxis.GetXmin(),yAxis.GetXmax());
826  }
827  ret->Sumw2();
828  ret->SetXTitle(xAxis.GetTitle());
829  ret->SetYTitle(yAxis.GetTitle());
830  ret->SetLineColor(color);
831  ret->SetMarkerColor(color);
832  ret->SetFillColor(color);
833  ret->SetMarkerStyle(style);
834  static_cast<const TAttAxis&>(xAxis).Copy(*(ret->GetXaxis()));
835  static_cast<const TAttAxis&>(yAxis).Copy(*(ret->GetYaxis()));
836  ret->SetDirectory(0);
837  if (const_cast<TAxis&>(xAxis).GetLabels()) {
838  for (Int_t i = 1; i <= xAxis.GetNbins(); i++)
839  ret->GetXaxis()->SetBinLabel(i, xAxis.GetBinLabel(i));
840  }
841  if (const_cast<TAxis&>(yAxis).GetLabels()) {
842  for (Int_t i = 1; i <= yAxis.GetNbins(); i++)
843  ret->GetYaxis()->SetBinLabel(i, yAxis.GetBinLabel(i));
844  }
845  c.Add(ret);
846  return ret;
847 }
848 //____________________________________________________________________
850 {
851  if (title && title[0] != '\0') axis.SetTitle(title);
852  axis. SetNdivisions(210);
853  axis. SetLabelFont(42);
854  axis. SetLabelSize(0.03);
855  axis. SetLabelOffset(0.005);
856  axis. SetLabelColor(kBlack);
857  axis. SetTitleOffset(1);
858  axis. SetTitleFont(42);
859  axis. SetTitleSize(0.04);
860  axis. SetTitleColor(kBlack);
861  axis. SetTickLength(0.03);
862  axis. SetAxisColor(kBlack);
863 }
864 //____________________________________________________________________
866  Double_t fact)
867 {
868  if (ret.GetXbins()->GetArray()) {
869  TArrayD bins(*ret.GetXbins());
870  for (Int_t i = 0; i < bins.GetSize(); i++) bins[i] *= fact;
871  ret.Set(ret.GetNbins(), bins.GetArray());
872  }
873  else {
874  ret.Set(ret.GetNbins(), fact*ret.GetXmin(), fact*ret.GetXmax());
875  }
876  FixAxis(ret);
877 }
878 //____________________________________________________________________
880 {
881  axis.Set(n, borders);
882  FixAxis(axis);
883 }
884 //____________________________________________________________________
886  const TString& spec,
887  const char* sep)
888 {
889  TString s(spec);
890  Bool_t isRange = false, isUnit = false;
891  if (s[0] == 'r' || s[0] == 'R') {
892  isRange = true;
893  s.Remove(0,1);
894  }
895  if (s[0] == 'u' || s[0] == 'U') {
896  isUnit = true;
897  s.Remove(0,1);
898  }
899  TObjArray* tokens = s.Tokenize(sep);
900  TArrayD bins(tokens->GetEntries());
901  TObjString* token = 0;
902  TIter next(tokens);
903  Int_t i = 0;
904  while ((token = static_cast<TObjString*>(next()))) {
905  Double_t v = token->String().Atof();
906  bins[i] = v;
907  i++;
908  }
909  delete tokens;
910  if (isUnit) {
911  if (bins.GetSize() > 1)
912  SetAxis(axis, Int_t(bins[1]-bins[0]), bins[0], bins[1]);
913  else
914  SetAxis(axis, 2*Int_t(bins[0]), bins[0]);
915  }
916  else if (isRange) {
917  Int_t nBins = Int_t(bins[0]);
918  if (bins.GetSize() > 2)
919  SetAxis(axis, nBins, bins[1], bins[2]);
920  else
921  SetAxis(axis, nBins, bins[1]);
922  }
923  else
924  SetAxis(axis, bins.GetSize()-1,bins.GetArray());
925 }
926 //____________________________________________________________________
928  Int_t n,
929  Double_t l,
930  Double_t h)
931 {
932  axis.Set(n, l, h);
933  FixAxis(axis);
934 }
935 //____________________________________________________________________
937 {
938  SetAxis(axis, n, -TMath::Abs(m), +TMath::Abs(m));
939 }
940 //____________________________________________________________________
942  Int_t nSig,
943  const char* alt)
944 {
945  printf(" %17s axis: ", alt ? alt : axis.GetTitle());
946  if (axis.GetXbins() && axis.GetXbins()->GetArray()) {
947  printf("%.*f", nSig, axis.GetBinLowEdge(1));
948  for (Int_t i = 1; i <= axis.GetNbins(); i++)
949  printf(":%.*f", nSig, axis.GetBinUpEdge(i));
950  }
951  else
952  printf("%d bins between %.*f and %.*f",
953  axis.GetNbins(), nSig, axis.GetXmin(),nSig,axis.GetXmax());
954  printf("\n");
955 }
956 //____________________________________________________________________
958 {
959  if (!h) {
960  ::Warning("ScaleToIPz","Nothing to scale");
961  return 0;
962  }
963  if (!ipZ) {
964  ::Warning("ScaleToIPz","Nothing to scale by");
965  return 0;
966  }
967  TH2* ret = static_cast<TH2*>(h->Clone());
968  ret->SetDirectory(0);
969  if (!ipZ) return ret;
970  for (Int_t iy = 1; iy <= ret->GetNbinsY(); iy++) {
971  Double_t z = ret->GetYaxis()->GetBinCenter(iy);
972  Int_t bin = ipZ->GetXaxis()->FindBin(z);
973  Double_t nEv = ipZ->GetBinContent(bin);
974  Double_t eEv = ipZ->GetBinError (bin);
975  Double_t esc = (nEv > 0 ? 1./nEv : 0);
976  Double_t rE2 = esc*esc*eEv*eEv;
977  for (Int_t ix = 1; ix <= ret->GetNbinsX(); ix++) {
978  Double_t c = ret->GetBinContent(ix,iy);
979  Double_t e = ret->GetBinError (ix,iy);
980  Double_t r = (c > 0 ? e/c : 0);
981  // Scale by number of events, and error propagate
982  Double_t sc = c * esc;
983  Double_t se = 0;
984  if (full) se = sc * TMath::Sqrt(r*r+rE2);
985  else se = e * esc;
986  Double_t scl = 1 / ret->GetXaxis()->GetBinWidth(ix);
987  ret->SetBinContent(ix, iy, scl*sc);
988  ret->SetBinError (ix, iy, scl*se);
989  }
990  }
991  return ret;
992 }
993 //____________________________________________________________________
995  const char* name,
996  UShort_t mode,
997  TH1* ipz,
998  TH2* other)
999 {
1000  if (!h) return 0;
1001 
1002  Int_t nIPz = h->GetNbinsY();
1003  Int_t nEta = h->GetNbinsX();
1004  TH1* p = h->ProjectionX(name,1,nIPz,"e");
1005  TH2* mask = (other ? other : h);
1006  p->SetDirectory(0);
1007  p->SetFillColor(0);
1008  p->SetFillStyle(0);
1009  p->SetYTitle(Form("#LT%s#GT", h->GetYaxis()->GetTitle()));
1010  p->Reset();
1011  for (Int_t etaBin = 1; etaBin <= nEta; etaBin++) {
1012  TArrayD hv(nIPz);
1013  TArrayD he(nIPz);
1014  TArrayD hr(nIPz);
1015  TArrayI hb(nIPz);
1016  Int_t j = 0;
1017  for (Int_t ipzBin = 1; ipzBin <= nIPz; ipzBin++) {
1018  Double_t bc = mask->GetBinContent(etaBin, ipzBin);
1019  if (bc < 1e-9) continue; // Low value
1020  Double_t be = mask->GetBinError (etaBin, ipzBin);
1021  if (TMath::IsNaN(be)) continue; // Bad error value
1022  hv[j] = bc;
1023  he[j] = be;
1024  hr[j] = be/bc;
1025  hb[j] = ipzBin;
1026  j++;
1027  }
1028  // Now we have all interesting values. Sort the relative error
1029  // array to get the most significant first. Note, we sort on the
1030  // mask errors which may not be the same as the histogram errors.
1031  TArrayI idx(nIPz);
1032  TMath::Sort(j, hr.fArray, idx.fArray, false);
1033  Double_t nsum = 0; // Running weighted sum
1034  Double_t nsume = 0; // Running weighted sum error
1035  Double_t dsum = 0;
1036  Double_t dsume = 0;
1037  Int_t n = 0;
1038  Double_t rat = 1e99;
1039 
1040  Int_t k = 0;
1041  for (k = 0; k < j; k++) {
1042  Int_t l = idx[k]; // Sorted index - ipzBin
1043  Int_t ipzBin = hb[l];
1044  Double_t hvv = hv[l];
1045  Double_t hee = he[l];
1046  Double_t x = TMath::Sqrt(nsume+hee*hee)/(nsum+hvv);
1047  if (x > rat) {
1048  continue; // Ignore - does not help
1049  }
1050 
1051  Double_t by = mask->GetYaxis()->GetBinCenter(ipzBin);
1052  Int_t ib = ipz ? ipz->FindBin(by) : 0;
1053  rat = x;
1054  nsum += h->GetBinContent(etaBin, ipzBin);
1055  nsume += TMath::Power(h->GetBinError(etaBin, ipzBin), 2);
1056  // If we do not have the vertex distribution, then just count
1057  // number of observations.
1058  dsum += !ipz ? 1 : ipz->GetBinContent(ib);
1059  dsume += !ipz ? 0 : TMath::Power(ipz->GetBinError(ib),2);
1060  n++;
1061  }
1062  if (k == 0 || n == 0) {
1063  ::Warning("Average", "Eta bin # %3d has no data",etaBin);
1064  continue; // This eta bin empty!
1065  }
1066  Double_t norm = (mode > 0 ? n : dsum);
1067  Double_t rne = nsume/nsum/nsum;
1068  Double_t rde = dsume/dsum/dsum;
1069  Double_t av = nsum/norm;
1070  Double_t ave = 0;
1071  if (mode==2) ave = TMath::Sqrt(nsume)/n;
1072  else ave = av*TMath::Sqrt(rne+rde);
1073  p->SetBinContent(etaBin, av);
1074  p->SetBinError (etaBin, ave);
1075  }
1076  if (mode == 0) p->Scale(1, "width");
1077  return p;
1078 }
1079 
1080 //____________________________________________________________________
1082 {
1083  if (!o) {
1084  ::Warning("CloneAndAdd", "Nothing to clone");
1085  return 0;
1086  }
1087  TObject* copy = o->Clone();
1088  if (copy->IsA()->InheritsFrom(TH1::Class()))
1089  // Release from underlying directory
1090  static_cast<TH1*>(copy)->SetDirectory(0);
1091  if (c)
1092  // Add to output container
1093  c->Add(copy);
1094  return copy;
1095 }
1096 //____________________________________________________________________
1098  Double_t min,
1099  Double_t max,
1100  Double_t& err)
1101 {
1102  const Double_t epsilon = 1e-6;
1103  Int_t bMin = h->FindBin(min+epsilon);
1104  Int_t bMax = h->FindBin(max-epsilon);
1105  if (bMin < 1) bMin = 0; // Include underflow bin
1106  Double_t val = h->IntegralAndError(bMin, bMax, err);
1107  // For a-symmetric histograms, return
1108  if (TMath::Abs(h->GetXaxis()->GetXmin()+h->GetXaxis()->GetXmax())>=epsilon)
1109  return val;
1110 
1111  // Otherwise, also integrate negative side
1112  Double_t err2;
1113  bMin = h->FindBin(-min+epsilon);
1114  bMax = h->FindBin(-max-epsilon);
1115  val += h->IntegralAndError(bMin, bMax, err2);
1116  err = TMath::Sqrt(err*err+err2*err2);
1117  return val;
1118 }
1119 //____________________________________________________________________
1121  Double_t d, Double_t ed,
1122  Double_t& er)
1123 {
1124  Double_t r = 0;
1125  er = 0;
1126  if (TMath::Abs(n) < 1e-16 || TMath::Abs(d) < 1e-16) return 0;
1127  r = n/d;
1128  er = TMath::Sqrt(en*en/n/n + ed*ed/d/d);
1129  return r;
1130 }
1131 //____________________________________________________________________
1133  Double_t d, Double_t e2d,
1134  Double_t& e2r)
1135 {
1136  Double_t r = 0;
1137  e2r = 0;
1138  if (TMath::Abs(n) < 1e-16 || TMath::Abs(d) < 1e-16) return 0;
1139  r = n/d;
1140  e2r = (e2n/n/n + e2d/d/d);
1141  return r;
1142 }
1143 
1144 #endif
1145 // Local Variables:
1146 // mode: C++
1147 // End:
1148 
Int_t color[]
option to what and if export to output file
static TH2 * CopyH2(Container *parent, const char *name, const char *newName=0)
double Double_t
Definition: External.C:58
const char * title
Definition: MakeQAPdf.C:27
static Bool_t CheckAxisNBins(const char *which, const TAxis *a1, const TAxis *a2)
static Bool_t CheckAxis(const char *which, const TAxis *a1, const TAxis *a2, Bool_t alsoLbls)
static TH1 * AverageOverIPz(TH2 *h, const char *name, UShort_t mode, TH1 *ipz, TH2 *mask=0)
static TH1 * Make1D(Container &c, const TString &name, const TString &title, Color_t color, Style_t style, const TAxis &xAxis)
static void SetAxis(TAxis &axis, Int_t n, Double_t *borders)
static Int_t GetI(Container *parent, const char *name, Int_t def=-1)
static Double_t RatioE(Double_t n, Double_t en, Double_t d, Double_t ed, Double_t &er)
static Double_t RatioE2(Double_t n, Double_t e2n, Double_t d, Double_t e2d, Double_t &e2r)
static Bool_t CheckAxisLabels(const char *which, const TAxis *a1, const TAxis *a2)
TCanvas * c
Definition: TestFitELoss.C:172
static Double_t Integrate(TH1 *h, Double_t min, Double_t max, Double_t &err)
static TH2 * ScaleToIPz(TH2 *h, TH1 *ipZ, Bool_t full=false)
static void PrintAxis(const TAxis &axis, Int_t nSig=2, const char *alt=0)
static Container * GetC(Container *parent, const char *name)
static TH1 * CopyH1(Container *parent, const char *name, const char *newName=0)
static Int_t GetB(Container *parent, const char *name, Bool_t def=false)
int Int_t
Definition: External.C:63
static Bool_t CheckAxisLimits(const char *which, const TAxis *a1, const TAxis *a2)
static void FixAxis(TAxis &axis, const char *title=0)
static void ScaleAxis(TAxis &axis, Double_t fact=1)
Definition: External.C:228
Definition: External.C:212
static TH1 * Scale(TH1 *h, Double_t x, Double_t xe)
Int_t mode
Definition: anaM.C:41
static Double_t GetD(Container *parent, const char *name, Double_t def=-1)
Definition: External.C:220
static TH1 * GetH1(Container *parent, const char *name)
unsigned short UShort_t
Definition: External.C:28
bool Bool_t
Definition: External.C:53
static Bool_t CheckConsistency(const TH1 *h1, const TH1 *h2)
static TH2 * Make2D(Container &c, const TString &name, const TString &title, Color_t color, Style_t style, const TAxis &xAxis, const TAxis &yAxis)
static TObject * CloneAndAdd(Container *c, TObject *o)
static TObject * GetO(Container *parent, const char *name, TClass *cls)
static Bool_t CheckAxisBins(const char *which, const TAxis *a1, const TAxis *a2)
Definition: External.C:196
static TH2 * GetH2(Container *parent, const char *name)