AliPhysics  a4b41ad (a4b41ad)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliTrackletAODUtils.C
Go to the documentation of this file.
1 
12 #ifndef ALITRACKLETAODUTILS_H
13 #define ALITRACKLETAODUTILS_H
14 #ifndef __CINT__
15 # include <TH1.h>
16 # include <TH2.h>
17 # include <TH3.h>
18 # include <TList.h>
19 # include <TParameter.h>
20 # include <TError.h>
21 # include <TMath.h>
22 # include <TFile.h>
23 # include <TDirectory.h>
24 # include <THashList.h>
25 # include <TProfile.h>
26 # include <TProfile2D.h>
27 # include <TDatabasePDG.h>
28 # include <THStack.h>
29 #else
30 class TList;
31 class TH1;
32 class TH2;
33 class TH3;
34 class TProfile;
35 class TProfile2D;
36 class TAxis;
37 class TDirectory;
38 class TFile;
39 class TDatabasePDG; // Auto-load
40 class THStack;
41 #endif
42 
56 {
57 public:
58  typedef TList Container;
66  virtual ~AliTrackletAODUtils() {}
67  //_________________________________________________________________________
81  static Bool_t CheckAxisNBins(const char* which,
82  const TAxis* a1,
83  const TAxis* a2);
93  static Bool_t CheckAxisLimits(const char* which,
94  const TAxis* a1,
95  const TAxis* a2);
106  static Bool_t CheckAxisBins(const char* which,
107  const TAxis* a1,
108  const TAxis* a2);
118  static Bool_t CheckAxisLabels(const char* which,
119  const TAxis* a1,
120  const TAxis* a2);
131  static Bool_t CheckAxis(const char* which,
132  const TAxis* a1,
133  const TAxis* a2,
134  Bool_t alsoLbls);
144  static Bool_t CheckConsistency(const TH1* h1, const TH1* h2);
145  /* @} */
146  //__________________________________________________________________
161  static TObject* GetO(Container* parent,
162  const char* name,
163  TClass* cls=0,
164  Bool_t verb=true);
175  static TObject* GetO(TDirectory* parent,
176  const char* name,
177  TClass* cls=0,
178  Bool_t verb=true);
188  static TH1* GetH1(Container* parent, const char* name, Bool_t verb=true);
198  static TH1* GetH1(TDirectory* parent, const char* name, Bool_t verb=true);
208  static TProfile* GetP1(Container* parent, const char* name, Bool_t verb=true);
218  static TH2* GetH2(Container* parent, const char* name, Bool_t verb=true);
228  static TH2* GetH2(TDirectory* parent, const char* name, Bool_t verb=true);
238  static TH3* GetH3(Container* parent, const char* name, Bool_t verb=true);
248  static TH3* GetH3(TDirectory* parent, const char* name, Bool_t verb=true);
258  static TProfile2D* GetP2(Container* parent, const char* name,
259  Bool_t verb=true);
269  static TProfile* GetP(Container* parent, const char* name, Bool_t verb=true)
270  {
271  return GetP1(parent, name, verb);
272  }
282  static THStack* GetHS(Container* parent, const char* name,Bool_t verb=true);
292  static THStack* GetHS(TDirectory* parent, const char* name,Bool_t verb=true);
302  static Container* GetC(Container* parent, const char* name, Bool_t verb=true);
312  static Container* GetC(TDirectory* parent, const char* name,Bool_t verb=true);
322  static TDirectory* GetT(TDirectory* parent,const char* name,Bool_t verb=true);
333  static Double_t GetD(Container* parent, const char* name,
334  Double_t def=-1, Bool_t verb=true);
345  static Int_t GetI(Container* parent, const char* name,
346  Int_t def=-1, Bool_t verb=true);
357  static Int_t GetB(Container* parent, const char* name,
358  Bool_t def=false, Bool_t verb=true);
369  static TH1* CopyH1(Container* parent,
370  const char* name,
371  const char* newName=0,
372  Bool_t verb=true);
383  static TH2* CopyH2(Container* parent,
384  const char* name,
385  const char* newName=0,
386  Bool_t verb=true);
397  static TH3* CopyH3(Container* parent,
398  const char* name,
399  const char* newName=0,
400  Bool_t verb=true);
407  static void CopyAttr(const TH1* src, TH1* dest);
420  static TH1* Make1D(Container* c,
421  const TString& name,
422  const TString& title,
423  Color_t color,
424  Style_t style,
425  const TAxis& xAxis);
438  static TProfile* Make1P(Container* c,
439  const TString& name,
440  const TString& title,
441  Color_t color,
442  Style_t style,
443  const TAxis& xAxis);
457  static TH2* Make2D(Container* c,
458  const TString& name,
459  const TString& title,
460  Color_t color,
461  Style_t style,
462  const TAxis& xAxis,
463  const TAxis& yAxis);
478  static TH3* Make3D(Container* c,
479  const TString& name,
480  const TString& title,
481  Color_t color,
482  Style_t style,
483  const TAxis& xAxis,
484  const TAxis& yAxis,
485  const TAxis& zAxis);
499  static TProfile2D* Make2P(Container* c,
500  const TString& name,
501  const TString& title,
502  Color_t color,
503  Style_t style,
504  const TAxis& xAxis,
505  const TAxis& yAxis);
506 
516  static TH1* Scale(TH1* h, Double_t x, Double_t xe);
526  static TH2* Scale(TH2* h, Double_t x, Double_t xe);
536  static TH2* Scale(TH2* h, TH1* s);
543  static void FixAxis(TAxis& axis, const char* title=0);
551  static void ScaleAxis(TAxis& ret, Double_t fact=1);
559  static void SetAxis(TAxis& axis, Int_t n, Double_t* borders);
568  static void SetAxis(TAxis& axis, const TString& spec, const char* sep=":,");
577  static void SetAxis(TAxis& axis, Int_t n, Double_t l, Double_t h);
585  static void SetAxis(TAxis& axis, Int_t n, Double_t m);
593  static void PrintAxis(const TAxis& axis, Int_t nSig=2, const char* alt=0);
604  static TH2* ScaleToIPz(TH2* h, TH1* ipZ, Bool_t full=false);
605  static TH3* ScaleToIPz(TH3* h, TH1* ipZ, Bool_t full=false);
606  static TH3* ScaleDelta(TH3* h, TH2* scale);
631  static TH2* ProjectEtaDelta(TH3* h);
647  static TH1* ProjectDelta(TH2* h);
663  static TH1* ProjectDeltaFull(TH3* h);
685  static TH1* AverageOverIPz(TH2* h,
686  const char* name,
687  UShort_t mode,
688  TH1* ipz,
689  TH2* mask=0,
690  Bool_t verb=true);
699  static TObject* CloneAndAdd(Container* c, TObject* o);
712  static Double_t Integrate(TH1* h, Double_t min, Double_t max, Double_t& err);
724  static Double_t RatioE(Double_t n, Double_t en,
725  Double_t d, Double_t ed,
726  Double_t& er);
738  static Double_t RatioE2(Double_t n, Double_t e2n,
739  Double_t d, Double_t e2d,
740  Double_t& e2r);
750  static TH1* RatioH(const TH1* num, const TH1* denom, const char* name=0);
757  static void FixMinMax(TH1* h, Bool_t ignoreZero=true);
764  static void FixMinMax(THStack* h, Bool_t ignoreZero=true);
765 
766  /* @} */
775  static void PdgAttr(Int_t pdg, TString& nme, Color_t& c, Style_t& s);
780  static Int_t* PdgArray(Int_t& size);
786  static const TAxis& PdgAxis();
794  static Int_t PdgBin(Int_t pdg);
795  /* @} */
804  static const char* CentName(Double_t c1, Double_t c2);
814  static Color_t CentColor(const TAxis& axis,
815  Double_t c1,
816  Double_t c2);
824  static Color_t CentColor(Int_t bin);
832  static TFile* OpenFile(const char* filename);
833 protected:
834  ClassDef(AliTrackletAODUtils,1); // Utilities
835 };
836 
837 //====================================================================
838 // Utilities
839 //____________________________________________________________________
841  const TAxis* a1,
842  const TAxis* a2)
843 {
844  if (a1->GetNbins() != a2->GetNbins()) {
845  ::Warning("CheckAxisNBins", "Incompatible number %s bins: %d vs %d",
846  which, a1->GetNbins(), a2->GetNbins());
847  return false;
848  }
849  return true;
850 }
851 //____________________________________________________________________
853  const TAxis* a1,
854  const TAxis* a2)
855 {
856  if (TMath::AreEqualAbs(a1->GetXmin(), a2->GetXmin(), 1.e-3) &&
857  TMath::AreEqualAbs(a1->GetXmax(), a2->GetXmax(), 1.e-3)) return true;
858  if (!TMath::AreEqualRel(a1->GetXmin(), a2->GetXmin(),1.E-12) ||
859  !TMath::AreEqualRel(a1->GetXmax(), a2->GetXmax(),1.E-12)) {
860  Warning("CheckAxisLimits",
861  "Limits of %s axis incompatible [%f,%f] vs [%f,%f]", which,
862  a1->GetXmin(), a1->GetXmax(), a2->GetXmin(), a2->GetXmax());
863  return false;
864  }
865  return true;
866 }
867 //____________________________________________________________________
869  const TAxis* a1,
870  const TAxis* a2)
871 {
872  const TArrayD * h1Array = a1->GetXbins();
873  const TArrayD * h2Array = a2->GetXbins();
874  Int_t fN = h1Array->fN;
875  if ( fN == 0 ) return true;
876  if (h2Array->fN != fN) {
877  // Redundant?
878  Warning("CheckAxisBins", "Not equal number of %s bin limits: %d vs %d",
879  which, fN, h2Array->fN);
880  return false;
881  }
882  else {
883  for (int i = 0; i < fN; ++i) {
884  if (TMath::AreEqualAbs(h1Array->GetAt(i),
885  h2Array->GetAt(i), 1.e-3)) continue;
886  if (!TMath::AreEqualRel(h1Array->GetAt(i),h2Array->GetAt(i),1E-10)) {
887  Warning("CheckAxisBins",
888  "%s limit # %3d incompatible: %f vs %f",
889  which, i, h1Array->GetAt(i),h2Array->GetAt(i));
890  return false;
891  }
892  }
893  }
894  return true;
895 }
896 //____________________________________________________________________
898  const TAxis* a1,
899  const TAxis* a2)
900 {
901  // check that axis have same labels
902  THashList *l1 = (const_cast<TAxis*>(a1))->GetLabels();
903  THashList *l2 = (const_cast<TAxis*>(a2))->GetLabels();
904 
905  if (!l1 && !l2) return true;
906  if (!l1 || !l2) {
907  Warning("CheckAxisLabels", "Difference in %s labels: %p vs %p",
908  which, l1, l2);
909  return false;
910  }
911  // check now labels sizes are the same
912  if (l1->GetSize() != l2->GetSize()) {
913  Warning("CheckAxisLabels", "Different number of %s labels: %d vs %d",
914  which, l1->GetSize(), l2->GetSize());
915  return false;
916  }
917  for (int i = 1; i <= a1->GetNbins(); ++i) {
918  TString label1 = a1->GetBinLabel(i);
919  TString label2 = a2->GetBinLabel(i);
920  if (label1 != label2) {
921  Warning("CheckAxisLabels", "%s label # %d not the same: '%s' vs '%s'",
922  which, i, label1.Data(), label2.Data());
923  return false;
924  }
925  }
926 
927  return true;
928 }
929 //____________________________________________________________________
931  const TAxis* a1,
932  const TAxis* a2,
933  Bool_t alsoLbls)
934 {
935  if (a1->GetNbins() == 0 && a2->GetNbins() == 0) return true;
936  if (!CheckAxisNBins (which, a1, a2)) return false;
937  if (!CheckAxisLimits(which, a1, a2)) return false;
938  if (!CheckAxisBins (which, a1, a2)) return false;
939  if (alsoLbls && !CheckAxisLabels(which, a1, a2)) return false;
940  return true;
941 }
942 //____________________________________________________________________
944 {
945  // Check histogram compatibility
946  if (h1 == h2) return true;
947 
948  if (h1->GetDimension() != h2->GetDimension() ) {
949  Warning("CheckConsistency",
950  "%s and %s have different dimensions %d vs %d",
951  h1->GetName(), h2->GetName(),
952  h1->GetDimension(), h2->GetDimension());
953  return false;
954  }
955  Int_t dim = h1->GetDimension();
956 
957  Bool_t alsoLbls = (h1->GetEntries() != 0 && h2->GetEntries() != 0);
958  if (!CheckAxis("X", h1->GetXaxis(), h2->GetXaxis(), alsoLbls)) return false;
959  if (dim > 1 &&
960  !CheckAxis("Y", h1->GetYaxis(), h2->GetYaxis(), alsoLbls)) return false;
961  if (dim > 2 &&
962  !CheckAxis("Z", h1->GetZaxis(), h2->GetZaxis(), alsoLbls)) return false;
963 
964  return true;
965 }
966 
967 //____________________________________________________________________
969 {
970  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
971  Double_t dc = s->GetBinContent(i);
972  Double_t de = s->GetBinError (i);
973  Double_t dr = (dc > 1e-6 ? de/dc : 0);
974  for (Int_t j = 1; j <= h->GetNbinsY(); j++) {
975  Double_t nc = h->GetBinContent(i,j);
976  Double_t ne = h->GetBinError (i,j);
977  Double_t ns = (nc > 0 ? ne/nc : 0);
978  Double_t sc = dc * nc;
979  Double_t se = sc*TMath::Sqrt(ns*ns+dr*dr);
980  // Printf("Setting bin %2d,%2d to %f +/- %f", sc, se);
981  h->SetBinContent(i,j,sc);
982  h->SetBinError (i,j,se);
983  }
984  }
985  return h;
986 }
987 
988 //____________________________________________________________________
990 {
991  Double_t rr = xe/x;
992  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
993  for (Int_t j = 1; j <= h->GetNbinsY(); j++) {
994  Double_t c = h->GetBinContent(i,j);
995  Double_t e = h->GetBinError (i,j);
996  Double_t s = (c > 0 ? e/c : 0);
997  Double_t sc = x * c;
998  Double_t se = sc*TMath::Sqrt(s*s+rr*rr);
999  // Printf("Setting bin %2d,%2d to %f +/- %f", sc, se);
1000  h->SetBinContent(i,j,sc);
1001  h->SetBinError (i,j,se);
1002  }
1003  }
1004  return h;
1005 }
1006 
1007 //____________________________________________________________________
1009 {
1010  Double_t rr = xe/x;
1011  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
1012  Double_t c = h->GetBinContent(i);
1013  Double_t e = h->GetBinError (i);
1014  Double_t s = (c > 0 ? e/c : 0);
1015  Double_t sc = x * c;
1016  Double_t se = sc*TMath::Sqrt(s*s+rr*rr);
1017  // Printf("Setting bin %2d,%2d to %f +/- %f", sc, se);
1018  h->SetBinContent(i,sc);
1019  h->SetBinError (i,se);
1020  }
1021  return h;
1022 }
1023 
1024 //____________________________________________________________________
1026  const char* name,
1027  TClass* cls,
1028  Bool_t verb)
1029 {
1030  if (!parent) {
1031  if (verb) ::Warning("GetO", "No parent container passed");
1032  return 0;
1033  }
1034  TObject* o = parent->FindObject(name);
1035  if (!o) {
1036  if (verb) ::Warning("GetO", "Object \"%s\" not found in \"%s\"",
1037  name, parent->GetName());
1038  // parent->ls();
1039  return 0;
1040  }
1041  if (!cls) return o;
1042  if (!o->IsA()->InheritsFrom(cls)) {
1043  if (verb) ::Warning("GetO", "\"%s\" is an object of class %s, not %s",
1044  name, o->ClassName(), cls->GetName());
1045  return 0;
1046  }
1047  return o;
1048 }
1049 //____________________________________________________________________
1051  const char* name,
1052  TClass* cls,
1053  Bool_t verb)
1054 {
1055  if (!parent) {
1056  if (!verb) ::Warning("GetO", "No parent directory passed");
1057  return 0;
1058  }
1059  TObject* o = parent->Get(name);
1060  if (!o) {
1061  if (verb) ::Warning("GetO", "Object \"%s\" not found in \"%s\"",
1062  name, parent->GetName());
1063  // parent->ls();
1064  return 0;
1065  }
1066  if (!cls) return o;
1067  if (!o->IsA()->InheritsFrom(cls)) {
1068  if (verb) ::Warning("GetO", "\"%s\" is an object of class %s, not %s",
1069  name, o->ClassName(), cls->GetName());
1070  return 0;
1071  }
1072  return o;
1073 }
1074 //____________________________________________________________________
1075 TH1* AliTrackletAODUtils::GetH1(Container* parent, const char* name, Bool_t v)
1076 {
1077  return static_cast<TH1*>(GetO(parent, name, TH1::Class(), v));
1078 }
1079 //____________________________________________________________________
1080 TH1* AliTrackletAODUtils::GetH1(TDirectory* parent, const char* name, Bool_t v)
1081 {
1082  return static_cast<TH1*>(GetO(parent, name, TH1::Class(), v));
1083 }
1084 //____________________________________________________________________
1085 TH2* AliTrackletAODUtils::GetH2(Container* parent, const char* name, Bool_t v)
1086 {
1087  return static_cast<TH2*>(GetO(parent, name, TH2::Class(), v));
1088 }
1089 //____________________________________________________________________
1090 TH2* AliTrackletAODUtils::GetH2(TDirectory* parent, const char* name, Bool_t v)
1091 {
1092  return static_cast<TH2*>(GetO(parent, name, TH2::Class(), v));
1093 }
1094 //____________________________________________________________________
1095 TH3* AliTrackletAODUtils::GetH3(Container* parent, const char* name, Bool_t v)
1096 {
1097  return static_cast<TH3*>(GetO(parent, name, TH3::Class(), v));
1098 }
1099 //____________________________________________________________________
1100 TH3* AliTrackletAODUtils::GetH3(TDirectory* parent, const char* name, Bool_t v)
1101 {
1102  return static_cast<TH3*>(GetO(parent, name, TH3::Class(), v));
1103 }
1104 //____________________________________________________________________
1105 TProfile* AliTrackletAODUtils::GetP1(Container* parent,const char* name,
1106  Bool_t v)
1107 {
1108  return static_cast<TProfile*>(GetO(parent, name, TProfile::Class(), v));
1109 }
1110 //____________________________________________________________________
1111 TProfile2D* AliTrackletAODUtils::GetP2(Container* parent, const char* name,
1112  Bool_t v)
1113 {
1114  return static_cast<TProfile2D*>(GetO(parent, name, TProfile2D::Class(), v));
1115 }
1116 //____________________________________________________________________
1117 THStack*
1118 AliTrackletAODUtils::GetHS(TDirectory* parent, const char* name, Bool_t v)
1119 {
1120  return static_cast<THStack*>(GetO(parent, name, THStack::Class(), v));
1121 }
1122 //____________________________________________________________________
1123 THStack*
1124 AliTrackletAODUtils::GetHS(Container* parent, const char* name, Bool_t v)
1125 {
1126  return static_cast<THStack*>(GetO(parent, name, THStack::Class(), v));
1127 }
1128 //____________________________________________________________________
1130 AliTrackletAODUtils::GetC(Container* parent, const char* name, Bool_t v)
1131 {
1132  return static_cast<Container*>(GetO(parent, name, Container::Class(), v));
1133 }
1134 //____________________________________________________________________
1136 AliTrackletAODUtils::GetC(TDirectory* parent, const char* name, Bool_t v)
1137 {
1138  return static_cast<Container*>(GetO(parent, name, Container::Class(), v));
1139 }
1140 //____________________________________________________________________
1141 TDirectory*
1142 AliTrackletAODUtils::GetT(TDirectory* parent, const char* name, Bool_t v)
1143 {
1144  TDirectory* d = parent->GetDirectory(name);
1145  if (!d) {
1146  if (v) ::Warning("GetO", "Directory \"%s\" not found in \"%s\"",
1147  name, parent->GetName());
1148  return 0;
1149  }
1150  return d;
1151 }
1152 //____________________________________________________________________
1154  const char* name,
1155  Double_t def,
1156  Bool_t v)
1157 {
1158  TParameter<double>* p =
1159  static_cast<TParameter<double>*>(GetO(parent, name,
1161  if (!p) return def;
1162  return p->GetVal();
1163 }
1164 //____________________________________________________________________
1166  const char* name,
1167  Int_t def,
1168  Bool_t v)
1169 {
1170  TParameter<int>* p =
1171  static_cast<TParameter<int>*>(GetO(parent, name,
1172  TParameter<int>::Class(),v));
1173  if (!p) return def;
1174  return p->GetVal();
1175 }
1176 //____________________________________________________________________
1178  const char* name,
1179  Bool_t def,
1180  Bool_t v)
1181 {
1182  TParameter<bool>* p =
1183  static_cast<TParameter<bool>*>(GetO(parent, name,
1184  TParameter<bool>::Class(), v));
1185  if (!p) return def;
1186  return p->GetVal();
1187 }
1188 //____________________________________________________________________
1190  const char* name,
1191  const char* newName,
1192  Bool_t v)
1193 {
1194  TH1* orig = GetH1(parent, name, v);
1195  if (!orig) return 0;
1196  TH1* ret = static_cast<TH1*>(orig->Clone(newName ? newName : name));
1197  ret->SetDirectory(0); // Release from file container
1198  return ret;
1199 }
1200 //____________________________________________________________________
1202  const char* name,
1203  const char* newName,
1204  Bool_t v)
1205 {
1206  TH2* orig = GetH2(parent, name, v);
1207  if (!orig) return 0;
1208  TH2* ret = static_cast<TH2*>(orig->Clone(newName ? newName : name));
1209  ret->SetDirectory(0); // Release from file container
1210  return ret;
1211 }
1212 //____________________________________________________________________
1214  const char* name,
1215  const char* newName,
1216  Bool_t v)
1217 {
1218  TH3* orig = GetH3(parent, name, v);
1219  if (!orig) return 0;
1220  TH3* ret = static_cast<TH3*>(orig->Clone(newName ? newName : name));
1221  ret->SetDirectory(0); // Release from file container
1222  return ret;
1223 }
1224 
1225 //____________________________________________________________________
1226 void AliTrackletAODUtils::CopyAttr(const TH1* src, TH1* dest)
1227 {
1228  if (!src || !dest) return;
1229  dest->SetMarkerStyle(src->GetMarkerStyle());
1230  dest->SetMarkerColor(src->GetMarkerColor());
1231  dest->SetMarkerSize (src->GetMarkerSize());
1232  dest->SetLineStyle (src->GetLineStyle());
1233  dest->SetLineColor (src->GetLineColor());
1234  dest->SetLineWidth (src->GetLineWidth());
1235  dest->SetFillStyle (src->GetFillStyle());
1236  dest->SetFillColor (src->GetFillColor());
1237 }
1238 
1239 
1240 //____________________________________________________________________
1242  const TString& name,
1243  const TString& title,
1244  Color_t color,
1245  Style_t style,
1246  const TAxis& xAxis)
1247 {
1248  TString n = name;
1249  TString t = title;
1250  TH1* ret = 0;
1251  Int_t nx = xAxis.GetNbins();
1252  if (t.IsNull()) t = xAxis.GetTitle();
1253  if (xAxis.GetXbins() && xAxis.GetXbins()->GetArray())
1254  ret = new TH1D(n,t,nx,xAxis.GetXbins()->GetArray());
1255  else
1256  ret = new TH1D(n,t,nx,xAxis.GetXmin(),xAxis.GetXmax());
1257  ret->Sumw2();
1258  ret->SetXTitle(xAxis.GetTitle());
1259  static_cast<const TAttAxis&>(xAxis).Copy(*(ret->GetXaxis()));
1260  ret->SetDirectory(0);
1261  ret->SetLineColor(color);
1262  ret->SetMarkerColor(color);
1263  ret->SetFillColor(kWhite);// color);
1264  ret->SetFillStyle(0);
1265  ret->SetMarkerStyle(style);
1266  if (const_cast<TAxis&>(xAxis).GetLabels()) {
1267  for (Int_t i = 1; i <= xAxis.GetNbins(); i++)
1268  ret->GetXaxis()->SetBinLabel(i, xAxis.GetBinLabel(i));
1269  }
1270  switch (style) {
1271  case 27:
1272  case 33:
1273  ret->SetMarkerSize(1.4);
1274  break;
1275  case 29:
1276  case 30:
1277  ret->SetMarkerSize(1.2);
1278  break;
1279  }
1280  if (c) c->Add(ret);
1281  return ret;
1282 }
1283 //____________________________________________________________________
1285  const TString& name,
1286  const TString& title,
1287  Color_t color,
1288  Style_t style,
1289  const TAxis& xAxis)
1290 {
1291  TString n = name;
1292  TString t = title;
1293  TProfile* ret = 0;
1294  Int_t nx = xAxis.GetNbins();
1295  if (t.IsNull()) t = xAxis.GetTitle();
1296  if (xAxis.GetXbins() && xAxis.GetXbins()->GetArray())
1297  ret = new TProfile(n,t,nx,xAxis.GetXbins()->GetArray());
1298  else
1299  ret = new TProfile(n,t,nx,xAxis.GetXmin(),xAxis.GetXmax());
1300  ret->Sumw2();
1301  ret->SetXTitle(xAxis.GetTitle());
1302  static_cast<const TAttAxis&>(xAxis).Copy(*(ret->GetXaxis()));
1303  ret->SetDirectory(0);
1304  ret->SetLineColor(color);
1305  ret->SetMarkerColor(color);
1306  ret->SetFillColor(kWhite);// color);
1307  ret->SetFillStyle(0);
1308  ret->SetMarkerStyle(style);
1309  if (const_cast<TAxis&>(xAxis).GetLabels()) {
1310  for (Int_t i = 1; i <= xAxis.GetNbins(); i++)
1311  ret->GetXaxis()->SetBinLabel(i, xAxis.GetBinLabel(i));
1312  }
1313  switch (style) {
1314  case 27:
1315  case 33:
1316  ret->SetMarkerSize(1.4);
1317  break;
1318  case 29:
1319  case 30:
1320  ret->SetMarkerSize(1.2);
1321  break;
1322  }
1323  if (c) c->Add(ret);
1324  return ret;
1325 }
1326 //____________________________________________________________________
1328  const TString& name,
1329  const TString& title,
1330  Color_t color,
1331  Style_t style,
1332  const TAxis& xAxis,
1333  const TAxis& yAxis)
1334 {
1335  TString n = name;
1336  TString t = title;
1337  TH2* ret = 0;
1338  Int_t nx = xAxis.GetNbins();
1339  Int_t ny = yAxis.GetNbins();
1340  const Double_t* xb = (xAxis.GetXbins() && xAxis.GetXbins()->GetArray() ?
1341  xAxis.GetXbins()->GetArray() : 0);
1342  const Double_t* yb = (yAxis.GetXbins() && yAxis.GetXbins()->GetArray() ?
1343  yAxis.GetXbins()->GetArray() : 0);
1344  if (t.IsNull())
1345  t.Form("%s vs %s", yAxis.GetTitle(), xAxis.GetTitle());
1346  if (xb) {
1347  if (yb) ret = new TH2D(n,t,nx,xb,ny,yb);
1348  else ret = new TH2D(n,t,
1349  nx,xb,
1350  ny,yAxis.GetXmin(),yAxis.GetXmax());
1351  }
1352  else {
1353  if (yb) ret = new TH2D(n,t,
1354  nx,xAxis.GetXmin(),xAxis.GetXmax(),
1355  ny,yb);
1356  else ret = new TH2D(n,t,
1357  nx,xAxis.GetXmin(),xAxis.GetXmax(),
1358  ny,yAxis.GetXmin(),yAxis.GetXmax());
1359  }
1360  ret->Sumw2();
1361  ret->SetXTitle(xAxis.GetTitle());
1362  ret->SetYTitle(yAxis.GetTitle());
1363  ret->SetLineColor(color);
1364  ret->SetMarkerColor(color);
1365  ret->SetFillColor(color);
1366  ret->SetMarkerStyle(style);
1367  static_cast<const TAttAxis&>(xAxis).Copy(*(ret->GetXaxis()));
1368  static_cast<const TAttAxis&>(yAxis).Copy(*(ret->GetYaxis()));
1369  ret->SetDirectory(0);
1370  if (const_cast<TAxis&>(xAxis).GetLabels()) {
1371  for (Int_t i = 1; i <= xAxis.GetNbins(); i++)
1372  ret->GetXaxis()->SetBinLabel(i, xAxis.GetBinLabel(i));
1373  }
1374  if (const_cast<TAxis&>(yAxis).GetLabels()) {
1375  for (Int_t i = 1; i <= yAxis.GetNbins(); i++)
1376  ret->GetYaxis()->SetBinLabel(i, yAxis.GetBinLabel(i));
1377  }
1378  if (c) c->Add(ret);
1379  return ret;
1380 }
1381 //____________________________________________________________________
1383  const TString& name,
1384  const TString& title,
1385  Color_t color,
1386  Style_t style,
1387  const TAxis& xAxis,
1388  const TAxis& yAxis,
1389  const TAxis& zAxis)
1390 {
1391  TString n = name;
1392  TString t = title;
1393  TH3* ret = 0;
1394  Int_t nx = xAxis.GetNbins();
1395  Int_t ny = yAxis.GetNbins();
1396  Int_t nz = zAxis.GetNbins();
1397  Double_t* xb = (xAxis.GetXbins() && xAxis.GetXbins()->GetArray() ?
1398  const_cast<Double_t*>(xAxis.GetXbins()->GetArray()) : 0);
1399  Double_t* yb = (yAxis.GetXbins() && yAxis.GetXbins()->GetArray() ?
1400  const_cast<Double_t*>(yAxis.GetXbins()->GetArray()) : 0);
1401  Double_t* zb = (zAxis.GetXbins() && zAxis.GetXbins()->GetArray() ?
1402  const_cast<Double_t*>(zAxis.GetXbins()->GetArray()) : 0);
1403  if (t.IsNull())
1404  t.Form("%s vs %s vs %s",
1405  zAxis.GetTitle(), yAxis.GetTitle(), xAxis.GetTitle());
1406  if (xb || yb || zb) {
1407  // One or more axis are defined as arrays. Make sure the rest are
1408  // also arrays.
1409  if (xb) {
1410  xb = new Double_t[nx+1];
1411  Double_t dx = (xAxis.GetXmax()-xAxis.GetXmin())/nx;
1412  xb[0] = xAxis.GetXmin();
1413  for (Int_t i = 1; i <= nx; i++) xb[i] = xb[i-1]+dx;
1414  }
1415  if (yb) {
1416  yb = new Double_t[ny+1];
1417  Double_t dy = (yAxis.GetXmax()-yAxis.GetXmin())/ny;
1418  yb[0] = yAxis.GetXmin();
1419  for (Int_t i = 1; i <= ny; i++) yb[i] = yb[i-1]+dy;
1420  }
1421  if (zb) {
1422  zb = new Double_t[nz+1];
1423  Double_t dz = (zAxis.GetXmax()-zAxis.GetXmin())/nz;
1424  zb[0] = zAxis.GetXmin();
1425  for (Int_t i = 1; i <= nz; i++) zb[i] = zb[i-1]+dz;
1426  }
1427  ret = new TH3D(n,t,nx,xb,ny,yb,nz,zb);
1428  }
1429  else {
1430  ret = new TH3D(n,t,
1431  nx, xAxis.GetXmin(), xAxis.GetXmax(),
1432  ny, yAxis.GetXmin(), yAxis.GetXmax(),
1433  nz, zAxis.GetXmin(), zAxis.GetXmax());
1434  }
1435  ret->Sumw2();
1436  ret->SetXTitle(xAxis.GetTitle());
1437  ret->SetYTitle(yAxis.GetTitle());
1438  ret->SetZTitle(zAxis.GetTitle());
1439  ret->SetLineColor(color);
1440  ret->SetMarkerColor(color);
1441  ret->SetFillColor(color);
1442  ret->SetMarkerStyle(style);
1443  static_cast<const TAttAxis&>(xAxis).Copy(*(ret->GetXaxis()));
1444  static_cast<const TAttAxis&>(yAxis).Copy(*(ret->GetYaxis()));
1445  static_cast<const TAttAxis&>(zAxis).Copy(*(ret->GetZaxis()));
1446  ret->SetDirectory(0);
1447  if (const_cast<TAxis&>(xAxis).GetLabels()) {
1448  for (Int_t i = 1; i <= xAxis.GetNbins(); i++)
1449  ret->GetXaxis()->SetBinLabel(i, xAxis.GetBinLabel(i));
1450  }
1451  if (const_cast<TAxis&>(yAxis).GetLabels()) {
1452  for (Int_t i = 1; i <= yAxis.GetNbins(); i++)
1453  ret->GetYaxis()->SetBinLabel(i, yAxis.GetBinLabel(i));
1454  }
1455  if (const_cast<TAxis&>(zAxis).GetLabels()) {
1456  for (Int_t i = 1; i <= zAxis.GetNbins(); i++)
1457  ret->GetZaxis()->SetBinLabel(i, zAxis.GetBinLabel(i));
1458  }
1459  if (c) c->Add(ret);
1460  return ret;
1461 }
1462 //____________________________________________________________________
1464  const TString& name,
1465  const TString& title,
1466  Color_t color,
1467  Style_t style,
1468  const TAxis& xAxis,
1469  const TAxis& yAxis)
1470 {
1471  TString n = name;
1472  TString t = title;
1473  TProfile2D* ret = 0;
1474  Int_t nx = xAxis.GetNbins();
1475  Int_t ny = yAxis.GetNbins();
1476  const Double_t* xb = (xAxis.GetXbins() && xAxis.GetXbins()->GetArray() ?
1477  xAxis.GetXbins()->GetArray() : 0);
1478  const Double_t* yb = (yAxis.GetXbins() && yAxis.GetXbins()->GetArray() ?
1479  yAxis.GetXbins()->GetArray() : 0);
1480  if (t.IsNull())
1481  t.Form("%s vs %s", yAxis.GetTitle(), xAxis.GetTitle());
1482  if (xb) {
1483  if (yb) ret = new TProfile2D(n,t,nx,xb,ny,yb);
1484  else ret = new TProfile2D(n,t,
1485  nx,xb,
1486  ny,yAxis.GetXmin(),yAxis.GetXmax());
1487  }
1488  else {
1489  if (yb) ret = new TProfile2D(n,t,
1490  nx,xAxis.GetXmin(),xAxis.GetXmax(),
1491  ny,yb);
1492  else ret = new TProfile2D(n,t,
1493  nx,xAxis.GetXmin(),xAxis.GetXmax(),
1494  ny,yAxis.GetXmin(),yAxis.GetXmax());
1495  }
1496  ret->Sumw2();
1497  ret->SetXTitle(xAxis.GetTitle());
1498  ret->SetYTitle(yAxis.GetTitle());
1499  ret->SetLineColor(color);
1500  ret->SetMarkerColor(color);
1501  ret->SetFillColor(color);
1502  ret->SetMarkerStyle(style);
1503  static_cast<const TAttAxis&>(xAxis).Copy(*(ret->GetXaxis()));
1504  static_cast<const TAttAxis&>(yAxis).Copy(*(ret->GetYaxis()));
1505  ret->SetDirectory(0);
1506  if (const_cast<TAxis&>(xAxis).GetLabels()) {
1507  for (Int_t i = 1; i <= xAxis.GetNbins(); i++)
1508  ret->GetXaxis()->SetBinLabel(i, xAxis.GetBinLabel(i));
1509  }
1510  if (const_cast<TAxis&>(yAxis).GetLabels()) {
1511  for (Int_t i = 1; i <= yAxis.GetNbins(); i++)
1512  ret->GetYaxis()->SetBinLabel(i, yAxis.GetBinLabel(i));
1513  }
1514  if (c) c->Add(ret);
1515  return ret;
1516 }
1517 //____________________________________________________________________
1518 void AliTrackletAODUtils::FixAxis(TAxis& axis, const char* title)
1519 {
1520  if (title && title[0] != '\0') axis.SetTitle(title);
1521  axis. SetNdivisions(210);
1522  axis. SetLabelFont(42);
1523  axis. SetLabelSize(0.03);
1524  axis. SetLabelOffset(0.005);
1525  axis. SetLabelColor(kBlack);
1526  axis. SetTitleOffset(1);
1527  axis. SetTitleFont(42);
1528  axis. SetTitleSize(0.04);
1529  axis. SetTitleColor(kBlack);
1530  axis. SetTickLength(0.03);
1531  axis. SetAxisColor(kBlack);
1532 }
1533 //____________________________________________________________________
1535  Double_t fact)
1536 {
1537  if (ret.GetXbins()->GetArray()) {
1538  TArrayD bins(*ret.GetXbins());
1539  for (Int_t i = 0; i < bins.GetSize(); i++) bins[i] *= fact;
1540  ret.Set(ret.GetNbins(), bins.GetArray());
1541  }
1542  else {
1543  ret.Set(ret.GetNbins(), fact*ret.GetXmin(), fact*ret.GetXmax());
1544  }
1545  FixAxis(ret);
1546 }
1547 //____________________________________________________________________
1549 {
1550  axis.Set(n, borders);
1551  FixAxis(axis);
1552 }
1553 //____________________________________________________________________
1555  const TString& spec,
1556  const char* sep)
1557 {
1558  TString s(spec);
1559  Bool_t isRange = false, isUnit = false;
1560  if (s[0] == 'r' || s[0] == 'R') {
1561  isRange = true;
1562  s.Remove(0,1);
1563  }
1564  if (s[0] == 'u' || s[0] == 'U') {
1565  isUnit = true;
1566  s.Remove(0,1);
1567  }
1568  TObjArray* tokens = s.Tokenize(sep);
1569  TArrayD bins(tokens->GetEntries());
1570  TObjString* token = 0;
1571  TIter next(tokens);
1572  Int_t i = 0;
1573  while ((token = static_cast<TObjString*>(next()))) {
1574  Double_t v = token->String().Atof();
1575  bins[i] = v;
1576  i++;
1577  }
1578  delete tokens;
1579  if (isUnit) {
1580  if (bins.GetSize() > 1)
1581  SetAxis(axis, Int_t(bins[1]-bins[0]), bins[0], bins[1]);
1582  else
1583  SetAxis(axis, 2*Int_t(bins[0]), bins[0]);
1584  }
1585  else if (isRange) {
1586  Int_t nBins = Int_t(bins[0]);
1587  if (bins.GetSize() > 2)
1588  SetAxis(axis, nBins, bins[1], bins[2]);
1589  else
1590  SetAxis(axis, nBins, bins[1]);
1591  }
1592  else
1593  SetAxis(axis, bins.GetSize()-1,bins.GetArray());
1594 }
1595 //____________________________________________________________________
1597  Int_t n,
1598  Double_t l,
1599  Double_t h)
1600 {
1601  axis.Set(n, l, h);
1602  FixAxis(axis);
1603 }
1604 //____________________________________________________________________
1606 {
1607  SetAxis(axis, n, -TMath::Abs(m), +TMath::Abs(m));
1608 }
1609 //____________________________________________________________________
1611  Int_t nSig,
1612  const char* alt)
1613 {
1614  printf(" %17s axis: ", alt ? alt : axis.GetTitle());
1615  if (axis.GetXbins() && axis.GetXbins()->GetArray()) {
1616  printf("%.*f", nSig, axis.GetBinLowEdge(1));
1617  for (Int_t i = 1; i <= axis.GetNbins(); i++)
1618  printf(":%.*f", nSig, axis.GetBinUpEdge(i));
1619  }
1620  else
1621  printf("%d bins between %.*f and %.*f",
1622  axis.GetNbins(), nSig, axis.GetXmin(),nSig,axis.GetXmax());
1623  printf("\n");
1624 }
1625 //____________________________________________________________________
1627 {
1628  if (!h) {
1629  ::Warning("ScaleToIPz","Nothing to scale");
1630  return 0;
1631  }
1632  if (!ipZ) {
1633  ::Warning("ScaleToIPz","Nothing to scale by");
1634  return 0;
1635  }
1636  TH2* ret = static_cast<TH2*>(h->Clone());
1637  ret->SetDirectory(0);
1638  if (!ipZ) return ret;
1639  for (Int_t iy = 1; iy <= ret->GetNbinsY(); iy++) {
1640  Double_t z = ret->GetYaxis()->GetBinCenter(iy);
1641  Int_t bin = ipZ->GetXaxis()->FindBin(z);
1642  Double_t nEv = ipZ->GetBinContent(bin);
1643  Double_t eEv = ipZ->GetBinError (bin);
1644  Double_t esc = (nEv > 0 ? 1./nEv : 0);
1645  Double_t rE2 = esc*esc*eEv*eEv;
1646  for (Int_t ix = 1; ix <= ret->GetNbinsX(); ix++) {
1647  Double_t c = ret->GetBinContent(ix,iy);
1648  Double_t e = ret->GetBinError (ix,iy);
1649  Double_t r = (c > 0 ? e/c : 0);
1650  // Scale by number of events, and error propagate
1651  Double_t sc = c * esc;
1652  Double_t se = 0;
1653  if (full) se = sc * TMath::Sqrt(r*r+rE2);
1654  else se = e * esc;
1655  Double_t scl = 1 / ret->GetXaxis()->GetBinWidth(ix);
1656  ret->SetBinContent(ix, iy, scl*sc);
1657  ret->SetBinError (ix, iy, scl*se);
1658  }
1659  }
1660  return ret;
1661 }
1662 //____________________________________________________________________
1664 {
1665  if (!h) {
1666  ::Warning("ScaleToIPz","Nothing to scale");
1667  return 0;
1668  }
1669  if (!ipZ) {
1670  ::Warning("ScaleToIPz","Nothing to scale by");
1671  return 0;
1672  }
1673  TH3* ret = static_cast<TH3*>(h->Clone());
1674  ret->SetDirectory(0);
1675  if (!ipZ) return ret;
1676  for (Int_t iz = 1; iz <= ret->GetNbinsZ(); iz++) {
1677  Double_t z = ret->GetZaxis()->GetBinCenter(iz);
1678  Int_t bin = ipZ->GetXaxis()->FindBin(z);
1679  Double_t nEv = ipZ->GetBinContent(bin);
1680  Double_t eEv = ipZ->GetBinError (bin);
1681  Double_t esc = (nEv > 0 ? 1./nEv : 0);
1682  Double_t rE2 = esc*esc*eEv*eEv;
1683  for (Int_t iy = 1; iy <= ret->GetNbinsY(); iy++) {
1684  for (Int_t ix = 1; ix <= ret->GetNbinsX(); ix++) {
1685  Double_t c = ret->GetBinContent(ix,iy,iz);
1686  Double_t e = ret->GetBinError (ix,iy,iz);
1687  Double_t r = (c > 0 ? e/c : 0);
1688  // Scale by number of events, and error propagate
1689  Double_t sc = c * esc;
1690  Double_t se = 0;
1691  if (full) se = sc * TMath::Sqrt(r*r+rE2);
1692  else se = e * esc;
1693  Double_t scl = 1 / ret->GetXaxis()->GetBinWidth(ix);
1694  ret->SetBinContent(ix, iy, iz, scl*sc);
1695  ret->SetBinError (ix, iy, iz, scl*se);
1696  }
1697  }
1698  }
1699  return ret;
1700 }
1701 //____________________________________________________________________
1703 {
1704  for (Int_t i = 1; i <= h->GetNbinsX(); i++) { // eta
1705  for (Int_t j = 1; j <= h->GetNbinsZ(); j++) { // IPz
1706  Double_t scale = etaIPzScale->GetBinContent(i,j);
1707  Double_t scaleE = etaIPzScale->GetBinError (i,j);
1708  Double_t q = (scale > 0 ? scaleE / scale : 0);
1709  for (Int_t k = 1; k <= h->GetNbinsY(); k++) { // delta
1710  Double_t c = h->GetBinContent(i,k,j);
1711  Double_t e = h->GetBinError (i,k,j);
1712  Double_t r = (c > 0 ? e / c : 0);
1713  Double_t w = c * scale;
1714  Double_t v = w * TMath::Sqrt(r*r + q*q);
1715  h->SetBinContent(i,k,j,w);
1716  h->SetBinError (i,k,j,v);
1717 #if 0
1718  Printf("%2d,%3d,%2d=(%9g+/-%9g)*(%9g+/-%9g)=(%9g+/-%9g)",
1719  i,k,j,c,e,scale,scaleE,w,v);
1720 #endif
1721  }
1722  }
1723  }
1724  return h;
1725 }
1726 //____________________________________________________________________
1727 // @cond
1729 {
1730  TH2* etaDelta = static_cast<TH2*>(h->Project3D("yx e"));
1731  etaDelta->SetName("etaDelta");
1732  etaDelta->SetTitle(h->GetTitle());
1733  etaDelta->SetDirectory(0);
1734  etaDelta->SetZTitle("d^{2}#it{N}/(d#etad#Delta)");
1735 #if 1
1736  // Reset content of projected histogram and calculate averages
1737  etaDelta->Reset();
1738  for (Int_t i = 1; i <= h->GetNbinsX(); i++) { // Loop over eta
1739  for (Int_t j = 1; j <= h->GetNbinsY(); j++) { // Loop over Delta
1740  Double_t sum = 0;
1741  Double_t sumw = 0;
1742  Int_t cnt = 0;
1743  for (Int_t k = 1; k <= h->GetNbinsZ(); k++) { // Loop over IPz
1744  Double_t c = h->GetBinContent(i,j,k);
1745  Double_t e = h->GetBinError (i,j,k);
1746  if (c < 1e-6 || e/c > 1) continue;
1747 #if 0
1748  Double_t w = 1/e/e;
1749  sum += w * c;
1750  sumw += w;
1751 #else
1752  sum += c;
1753  sumw += e*e;
1754 #endif
1755  cnt += 1;
1756  }
1757  if (sumw < 1e-6) continue;
1758 #if 0
1759  etaDelta->SetBinContent(i,j,sum/sumw);
1760  etaDelta->SetBinError (i,j,TMath::Sqrt(1/sumw));
1761 #else
1762  etaDelta->SetBinContent(i,j,sum/cnt);
1763  etaDelta->SetBinError (i,j,TMath::Sqrt(sumw)/cnt);
1764 #endif
1765  }
1766  }
1767 #endif
1768  return etaDelta;
1769 }
1770 // @endcond
1771 //____________________________________________________________________
1773 {
1774  if (!h) {
1775  Warning("ProjectDelta", "No histogram passed");
1776  return 0;
1777  }
1778  TH1* delta = h->ProjectionY("delta");
1779  delta->SetDirectory(0);
1780  delta->SetTitle(h->GetTitle());
1781  delta->SetYTitle("d#it{N}/d#Delta");
1782  delta->Scale(1./h->GetNbinsX());
1783  return delta;
1784 }
1785 //____________________________________________________________________
1787 {
1788  if (!h) {
1789  Warning("ProjectDelta", "No histogram passed");
1790  return 0;
1791  }
1792  TH2* tmp = ProjectEtaDelta(h);
1793  TH1* delta = ProjectDelta(tmp);
1794  delta->SetDirectory(0);
1795  tmp->SetDirectory(0);
1796  delete tmp;
1797  return delta;
1798 }
1799 
1800 //____________________________________________________________________
1802  const char* name,
1803  UShort_t mode,
1804  TH1* ipz,
1805  TH2* other,
1806  Bool_t verb)
1807 {
1808  if (!h) return 0;
1809 
1810  Int_t nIPz = h->GetNbinsY();
1811  Int_t nEta = h->GetNbinsX();
1812  TH1* p = h->ProjectionX(name,1,nIPz,"e");
1813  TH2* mask = (other ? other : h);
1814  p->SetDirectory(0);
1815  p->SetFillColor(0);
1816  p->SetFillStyle(0);
1817  p->SetYTitle(Form("#LT%s#GT", h->GetYaxis()->GetTitle()));
1818  p->Reset();
1819  for (Int_t etaBin = 1; etaBin <= nEta; etaBin++) {
1820  TArrayD hv(nIPz);
1821  TArrayD he(nIPz);
1822  TArrayD hr(nIPz);
1823  TArrayI hb(nIPz);
1824  Int_t j = 0;
1825  for (Int_t ipzBin = 1; ipzBin <= nIPz; ipzBin++) {
1826  Double_t bc = mask->GetBinContent(etaBin, ipzBin);
1827  if (bc < 1e-9) continue; // Low value
1828  Double_t be = mask->GetBinError (etaBin, ipzBin);
1829  if (TMath::IsNaN(be)) continue; // Bad error value
1830  hv[j] = bc;
1831  he[j] = be;
1832  hr[j] = be/bc;
1833  hb[j] = ipzBin;
1834  j++;
1835  }
1836  // Now we have all interesting values. Sort the relative error
1837  // array to get the most significant first. Note, we sort on the
1838  // mask errors which may not be the same as the histogram errors.
1839  TArrayI idx(nIPz);
1840  TMath::Sort(j, hr.fArray, idx.fArray, false);
1841  Double_t nsum = 0; // Running weighted sum
1842  Double_t nsume = 0; // Running weighted sum error
1843  Double_t dsum = 0;
1844  Double_t dsume = 0;
1845  Int_t n = 0;
1846  Double_t rat = 1e99;
1847 
1848  Int_t k = 0;
1849  for (k = 0; k < j; k++) {
1850  Int_t l = idx[k]; // Sorted index - ipzBin
1851  Int_t ipzBin = hb[l];
1852  Double_t hvv = hv[l];
1853  Double_t hee = he[l];
1854  Double_t x = TMath::Sqrt(nsume+hee*hee)/(nsum+hvv);
1855  if (x > rat) {
1856  continue; // Ignore - does not help
1857  }
1858 
1859  Double_t by = mask->GetYaxis()->GetBinCenter(ipzBin);
1860  Int_t ib = ipz ? ipz->FindBin(by) : 0;
1861  rat = x;
1862  nsum += h->GetBinContent(etaBin, ipzBin);
1863  nsume += TMath::Power(h->GetBinError(etaBin, ipzBin), 2);
1864  // If we do not have the vertex distribution, then just count
1865  // number of observations.
1866  dsum += !ipz ? 1 : ipz->GetBinContent(ib);
1867  dsume += !ipz ? 0 : TMath::Power(ipz->GetBinError(ib),2);
1868  n++;
1869  }
1870  if (k == 0 || n == 0) {
1871  if (verb)
1872  ::Warning("Average", "Eta bin # %3d of %10s has no data",
1873  etaBin, h->GetName());
1874  continue; // This eta bin empty!
1875  }
1876  Double_t norm = (mode > 0 ? n : dsum);
1877  Double_t rne = nsume/nsum/nsum;
1878  Double_t rde = dsume/dsum/dsum;
1879  Double_t av = nsum/norm;
1880  Double_t ave = 0;
1881  if (mode==2) ave = TMath::Sqrt(nsume)/n;
1882  else ave = av*TMath::Sqrt(rne+rde);
1883  p->SetBinContent(etaBin, av);
1884  p->SetBinError (etaBin, ave);
1885  }
1886  if (mode == 0) p->Scale(1, "width");
1887  return p;
1888 }
1889 
1890 //____________________________________________________________________
1892 {
1893  if (!o) {
1894  ::Warning("CloneAndAdd", "Nothing to clone");
1895  return 0;
1896  }
1897  TObject* copy = o->Clone();
1898  if (copy->IsA()->InheritsFrom(TH1::Class()))
1899  // Release from underlying directory
1900  static_cast<TH1*>(copy)->SetDirectory(0);
1901  if (c)
1902  // Add to output container
1903  c->Add(copy);
1904  return copy;
1905 }
1906 //____________________________________________________________________
1908  Double_t min,
1909  Double_t max,
1910  Double_t& err)
1911 {
1912  const Double_t epsilon = 1e-6;
1913  Int_t bMin = h->FindBin(min+epsilon);
1914  Int_t bMax = h->FindBin(max-epsilon);
1915  if (bMin < 1) bMin = 0; // Include underflow bin
1916  Double_t val = h->IntegralAndError(bMin, bMax, err);
1917  // For a-symmetric histograms, return
1918  if (TMath::Abs(h->GetXaxis()->GetXmin()+h->GetXaxis()->GetXmax())>=epsilon)
1919  return val;
1920 
1921  // Otherwise, also integrate negative side
1922  Double_t err2;
1923  bMin = h->FindBin(-min+epsilon);
1924  bMax = h->FindBin(-max-epsilon);
1925  val += h->IntegralAndError(bMin, bMax, err2);
1926  err = TMath::Sqrt(err*err+err2*err2);
1927  return val;
1928 }
1929 //____________________________________________________________________
1931  Double_t d, Double_t ed,
1932  Double_t& er)
1933 {
1934  Double_t r = 0;
1935  er = 0;
1936  if (TMath::Abs(n) < 1e-16 || TMath::Abs(d) < 1e-16) return 0;
1937  r = n/d;
1938  er = TMath::Sqrt(en*en/n/n + ed*ed/d/d);
1939  return r;
1940 }
1941 //____________________________________________________________________
1943  Double_t d, Double_t e2d,
1944  Double_t& e2r)
1945 {
1946  Double_t r = 0;
1947  e2r = 0;
1948  if (TMath::Abs(n) < 1e-16 || TMath::Abs(d) < 1e-16) return 0;
1949  r = n/d;
1950  e2r = (e2n/n/n + e2d/d/d);
1951  return r;
1952 }
1953 //____________________________________________________________________
1955  const TH1* denom,
1956  const char* name)
1957 {
1958  if (!num || !denom) {
1959  ::Warning("RatioH", "Numerator (%p) or denominator (%p) missing",
1960  num, denom);
1961  return 0;
1962  }
1963  if (!CheckConsistency(num, denom)) return 0;
1964 
1965  TH1* ratio = static_cast<TH1*>(num->Clone(name ? name : num->GetName()));
1966  ratio->SetDirectory(0);
1967  ratio->SetMinimum(-1111);
1968  ratio->SetMaximum(-1111);
1969  ratio->GetListOfFunctions()->Clear();
1970  ratio->Divide(denom);
1971  ratio->SetYTitle(Form("%s / %s",
1972  num ->GetYaxis()->GetTitle(),
1973  denom->GetYaxis()->GetTitle()));
1974  ratio->SetTitle(Form("%s / %s", num ->GetTitle(), denom->GetTitle()));
1975  FixMinMax(ratio, true);
1976 
1977  return ratio;
1978 }
1979 //____________________________________________________________________
1981 {
1982  Double_t min = +1e99;
1983  Double_t max = -1e99;
1984  for (Int_t ix = 1; ix <= h->GetNbinsX(); ix++) {
1985  for (Int_t iy = 1; iy <= h->GetNbinsY(); iy++) {
1986  for (Int_t iz = 1; iz <= h->GetNbinsZ(); iz++) {
1987  Int_t bin = h->GetBin(ix,iy,iz);
1988  Double_t c = h->GetBinContent(bin);
1989  Double_t e = h->GetBinError (bin);
1990  if (c == 0 && ignoreZero) continue;
1991  min = TMath::Min(min, c-e);
1992  max = TMath::Max(max, c+e);
1993  }
1994  }
1995  }
1996  h->SetMinimum(min);
1997  h->SetMaximum(max);
1998 }
1999 //____________________________________________________________________
2000 void AliTrackletAODUtils::FixMinMax(THStack* stack, Bool_t ignoreZero)
2001 {
2002  Double_t min = +1e99;
2003  Double_t max = -1e99;
2004  TIter next(stack->GetHists());
2005  TH1* h = 0;
2006  while ((h = static_cast<TH1*>(next()))) {
2007  FixMinMax(h, ignoreZero);
2008  min = TMath::Min(min, h->GetMinimum());
2009  max = TMath::Max(max, h->GetMaximum());
2010  }
2011  stack->SetMinimum(min);
2012  stack->SetMaximum(max);
2013 }
2014 
2015 //____________________________________________________________________
2017 {
2018  static Int_t codes[] = {
2019  211, // #pi^{+}
2020  2212, // p
2021  321, // K^{+}
2022  323, // K^{*+}
2023  11, // e^{-}
2024  13, // #mu^{-}
2025  213, // #rho^{+}
2026  411, // D^{+}
2027  413, // D^{*+}
2028  431, // D_{s}^{+}
2029  433, // D_{s}^{*+}
2030  1114, // #Delta^{-}
2031  2214, // #Delta^{+}
2032  2224, // #Delta^{++}
2033  3112, // #Sigma^{-}
2034  3222, // #Sigma^{+}
2035  3114, // #Sigma^{*-}
2036  3224, // #Sigma^{*+}
2037  4214, // #Sigma^{*+}_{c}
2038  4224, // #Sigma^{*++}_{c}
2039  3312, // #Xi^{-}
2040  3314, // #Xi^{*-}
2041  4122, // #Lambda^{+}_{c}
2042  2112, // n
2043  2114, // #Delta^{0}
2044  22, // #gamma
2045  310, // K^{0}_{S}
2046  130, // K^{0}_{L}
2047  311, // K^{0}
2048  313, // K^{*}
2049  221, // #eta
2050  111, // #pi^{0}
2051  113, // #rho^{0}
2052  333, // #varphi
2053  331, // #eta'
2054  223, // #omega
2055  3122, // #Lambda
2056  3212, // #Sigma^{0}
2057  4114, // #Sigma^{*0}_{c}
2058  3214, // #Sigma^{*0}
2059  421, // D^{0}
2060  423, // D^{*0}
2061  3322, // #Xi_{0}
2062  3324, // #Xi^{*0}
2063  4132, // #Xi^{0}_{c}
2064  4314, // #Xi^{*0}_{c}
2065  1000000000 // Nuclei
2066  };
2067  static Int_t nCodes = sizeof(codes) / sizeof(Int_t);
2068  static Int_t* sorted = 0;
2069  size = nCodes;
2070  if (sorted) return sorted;
2071 
2072  sorted = new Int_t[nCodes];
2073  Int_t* idx = new Int_t[nCodes];
2074  TMath::Sort(nCodes, codes, idx, false);
2075  for (Int_t i = 0; i < nCodes; i++) sorted[i] = codes[idx[i]];
2076  delete [] idx;
2077  return sorted;
2078 }
2079 
2080 //____________________________________________________________________
2082  Color_t& c, Style_t& s)
2083 {
2084  // Leptons are black
2085  // Non-strange, non-charm meson are red
2086  // Non-strange, non-charm baryons are magenta
2087  // Strange-mesons are blue
2088  // Strange-baryons are cyan
2089  // Charmed mesons are green
2090  // Charmed baryons are yellow
2091  switch (pdg) {
2092  case -1: c=kPink +7; s=20; nme="?"; break;
2093  case 11: c=kBlack +0; s=20; nme="e^{-}"; break;
2094  case 13: c=kBlack +0; s=21; nme="#mu^{-}"; break;
2095  case 22: c=kBlack +0; s=22; nme="#gamma"; break;
2096  case 111: c=kRed +2; s=20; nme="#pi^{0}"; break;
2097  case 211: c=kRed +2; s=21; nme="#pi^{+}"; break;
2098  case 213: c=kRed +2; s=22; nme="#rho^{+}"; break;
2099  case 113: c=kRed +2; s=23; nme="#rho^{0}"; break;
2100  case 223: c=kRed +2; s=24; nme="#omega"; break;
2101  case 321: c=kBlue +2; s=20; nme="K^{+}"; break;
2102  case 323: c=kBlue +2; s=21; nme="K^{*+}"; break;
2103  case 310: c=kBlue +2; s=22; nme="K^{0}_{S}"; break;
2104  case 130: c=kBlue +2; s=23; nme="K^{0}_{L}"; break;
2105  case 311: c=kBlue +2; s=24; nme="K^{0}"; break;
2106  case 313: c=kBlue +2; s=25; nme="K^{*}"; break;
2107  case 221: c=kBlue +2; s=26; nme="#eta"; break;
2108  case 333: c=kBlue +2; s=27; nme="#varphi"; break;
2109  case 331: c=kBlue +2; s=28; nme="#eta'"; break;
2110  case 411: c=kGreen +2; s=20; nme="D^{+}"; break;
2111  case 413: c=kGreen +2; s=21; nme="D^{*+}"; break;
2112  case 421: c=kGreen +2; s=22; nme="D^{0}"; break;
2113  case 423: c=kGreen +2; s=23; nme="D^{*0}"; break;
2114  case 431: c=kGreen +2; s=24; nme="D_{s}^{+}"; break;
2115  case 433: c=kGreen +2; s=25; nme="D_{s}^{*+}"; break;
2116  case 2212: c=kMagenta+2; s=20; nme="p"; break;
2117  case 2112: c=kMagenta+2; s=21; nme="n"; break;
2118  case 2114: c=kMagenta+2; s=22; nme="#Delta^{0}"; break;
2119  case 1114: c=kMagenta+2; s=23; nme="#Delta^{-}"; break;
2120  case 2214: c=kMagenta+2; s=24; nme="#Delta^{+}"; break;
2121  case 2224: c=kMagenta+2; s=25; nme="#Delta^{++}"; break;
2122  case 3112: c=kCyan +2; s=20; nme="#Sigma^{-}"; break;
2123  case 3222: c=kCyan +2; s=21; nme="#Sigma^{+}"; break;
2124  case 3114: c=kCyan +2; s=22; nme="#Sigma^{*-}"; break;
2125  case 3224: c=kCyan +2; s=23; nme="#Sigma^{*+}"; break;
2126  case 3312: c=kCyan +2; s=24; nme="#Xi^{-}"; break;
2127  case 3314: c=kCyan +2; s=25; nme="#Xi^{*-}"; break;
2128  case 3122: c=kCyan +2; s=26; nme="#Lambda"; break;
2129  case 3212: c=kCyan +2; s=27; nme="#Sigma^{0}"; break;
2130  case 3214: c=kCyan +2; s=28; nme="#Sigma^{*0}"; break;
2131  case 3322: c=kCyan +2; s=29; nme="#Xi^{0}"; break;
2132  case 3324: c=kCyan +2; s=30; nme="#Xi^{*0}"; break;
2133  case 4214: c=kYellow +2; s=20; nme="#Sigma^{*+}_{c}"; break;
2134  case 4224: c=kYellow +2; s=21; nme="#Sigma^{*++}_{c}";break;
2135  case 4122: c=kYellow +2; s=22; nme="#Lambda^{+}_{c}"; break;
2136  case 4114: c=kYellow +2; s=23; nme="#Sigma^{*0}_{c}"; break;
2137  case 4132: c=kYellow +2; s=24; nme="#Xi^{0}_{c}"; break;
2138  case 4314: c=kYellow +2; s=25; nme="#Xi^{*0}_{c}"; break;
2139  case 1000000000:c=kPink +2; s=20; nme="Nuclei"; break;
2140  default: c=kGray; s=1; nme.Form("%d",pdg); break;
2141  };
2142 }
2143 
2144 //____________________________________________________________________
2146 {
2147  Int_t size;
2148  Int_t* array = PdgArray(size);
2149  Int_t apdg = TMath::Abs(pdg);
2150  Int_t idx = TMath::BinarySearch(size, array, apdg);
2151  if (idx == size - 1) return idx+1;
2152  if (array[idx] != apdg) return size+1;
2153  return idx+1;
2154 }
2155 
2156 //____________________________________________________________________
2158 {
2159  static TAxis* axis = 0;
2160  if (axis) return *axis;
2161 
2162  Int_t size;
2163  Int_t* array = PdgArray(size);
2164  axis = new TAxis(size+1, +.5, size+1.5);
2165  // TDatabasePDG* pdgDb = TDatabasePDG::Instance();
2166  for (Int_t i = 1; i <= size; i++) {
2167  Int_t pdg = array[i-1];
2168  axis->SetBinLabel(i, Form("%d", pdg));
2169  // TString name = "?";
2170  // TParticlePDG* type = pdgDb->GetParticle(pdg);
2171  // if (type) name = type->GetName();
2172  // axis->SetBinLabel(i, name.Data());
2173  }
2174  axis->SetBinLabel(size+1, "-1");
2175  FixAxis(*axis);
2176  return *axis;
2177 }
2178 //____________________________________________________________________
2180 {
2181  TFile* file = TFile::Open(filename, "READ");
2182  if (!file) {
2183  ::Warning("OpenFile", "Failed to open \"%s\"", filename);
2184  return 0;
2185  }
2186  return file;
2187 }
2188 //____________________________________________________________________
2190 {
2191  static TString tmp;
2192  tmp.Form("cent%03dd%02d_%03dd%02d",
2193  Int_t(c1), Int_t(c1*100)%100,
2194  Int_t(c2), Int_t(c2*100)%100);
2195  return tmp.Data();
2196 }
2197 //____________________________________________________________________
2199 {
2200  const Color_t cc[] = { kMagenta+2, // 0
2201  kBlue+2, // 1
2202  kAzure-1, // 2 // 10,
2203  kCyan+2, // 3
2204  kGreen+1, // 4
2205  kSpring+5, // 5 //+10,
2206  kYellow+1, // 6
2207  kOrange+5, // 7 //+10,
2208  kRed+1, // 8
2209  kPink+5, // 9 //+10,
2210  kBlack }; // 10
2211  Color_t tc = cc[bin % 10];
2212  return tc;
2213 }
2214 //____________________________________________________________________
2216  Double_t c1,
2217  Double_t c2)
2218 {
2219  Double_t avgC = (c1+c2)/2;
2220  Int_t binC = const_cast<TAxis&>(axis).FindBin(avgC);
2221  return CentColor(binC);
2222 }
2223 
2224 
2225 
2226 #endif
2227 // Local Variables:
2228 // mode: C++
2229 // End:
2230 
Int_t color[]
static TProfile * GetP1(Container *parent, const char *name, Bool_t verb=true)
static Double_t RatioE2(Double_t n, Double_t e2n, Double_t d, Double_t e2d, Double_t &e2r)
Int_t pdg
const char * filename
Definition: TestFCM.C:1
static const TAxis & PdgAxis()
static void PdgAttr(Int_t pdg, TString &nme, Color_t &c, Style_t &s)
static TH2 * Make2D(Container *c, const TString &name, const TString &title, Color_t color, Style_t style, const TAxis &xAxis, const TAxis &yAxis)
static TDirectory * GetT(TDirectory *parent, const char *name, Bool_t verb=true)
const Color_t cc[]
Definition: DrawKs.C:1
double Double_t
Definition: External.C:58
static void ScaleAxis(TAxis &ret, Double_t fact=1)
const char * title
Definition: MakeQAPdf.C:26
Definition: External.C:244
static void PrintAxis(const TAxis &axis, Int_t nSig=2, const char *alt=0)
static void SetAxis(TAxis &axis, Int_t n, Double_t *borders)
static TFile * OpenFile(const char *filename)
static TProfile2D * Make2P(Container *c, const TString &name, const TString &title, Color_t color, Style_t style, const TAxis &xAxis, const TAxis &yAxis)
static Bool_t CheckAxisBins(const char *which, const TAxis *a1, const TAxis *a2)
static TH1 * GetH1(Container *parent, const char *name, Bool_t verb=true)
static const char * CentName(Double_t c1, Double_t c2)
static Bool_t CheckAxisLabels(const char *which, const TAxis *a1, const TAxis *a2)
static TH3 * ScaleDelta(TH3 *h, TH2 *scale)
TCanvas * c
Definition: TestFitELoss.C:172
static Bool_t CheckAxis(const char *which, const TAxis *a1, const TAxis *a2, Bool_t alsoLbls)
static Int_t GetI(Container *parent, const char *name, Int_t def=-1, Bool_t verb=true)
static TH2 * ScaleToIPz(TH2 *h, TH1 *ipZ, Bool_t full=false)
static TH2 * ProjectEtaDelta(TH3 *h)
static TH1 * AverageOverIPz(TH2 *h, const char *name, UShort_t mode, TH1 *ipz, TH2 *mask=0, Bool_t verb=true)
static TH1 * ProjectDelta(TH2 *h)
static Double_t GetD(Container *parent, const char *name, Double_t def=-1, Bool_t verb=true)
static TH3 * Make3D(Container *c, const TString &name, const TString &title, Color_t color, Style_t style, const TAxis &xAxis, const TAxis &yAxis, const TAxis &zAxis)
static Bool_t CheckConsistency(const TH1 *h1, const TH1 *h2)
static TH1 * Scale(TH1 *h, Double_t x, Double_t xe)
static Int_t GetB(Container *parent, const char *name, Bool_t def=false, Bool_t verb=true)
static void FixMinMax(TH1 *h, Bool_t ignoreZero=true)
int Int_t
Definition: External.C:63
static TProfile2D * GetP2(Container *parent, const char *name, Bool_t verb=true)
static TH1 * ProjectDeltaFull(TH3 *h)
static Int_t * PdgArray(Int_t &size)
static Double_t Integrate(TH1 *h, Double_t min, Double_t max, Double_t &err)
Definition: External.C:252
static Bool_t CheckAxisNBins(const char *which, const TAxis *a1, const TAxis *a2)
Definition: External.C:228
Definition: External.C:212
ClassDef(AliTrackletAODUtils, 1)
static Bool_t CheckAxisLimits(const char *which, const TAxis *a1, const TAxis *a2)
static Color_t CentColor(const TAxis &axis, Double_t c1, Double_t c2)
static TH2 * CopyH2(Container *parent, const char *name, const char *newName=0, Bool_t verb=true)
Int_t mode
Definition: anaM.C:40
static TH1 * CopyH1(Container *parent, const char *name, const char *newName=0, Bool_t verb=true)
static TH3 * GetH3(Container *parent, const char *name, Bool_t verb=true)
static Double_t RatioE(Double_t n, Double_t en, Double_t d, Double_t ed, Double_t &er)
static TH2 * GetH2(Container *parent, const char *name, Bool_t verb=true)
static TObject * GetO(Container *parent, const char *name, TClass *cls=0, Bool_t verb=true)
static THStack * GetHS(Container *parent, const char *name, Bool_t verb=true)
Definition: External.C:220
TFile * file
static Container * GetC(Container *parent, const char *name, Bool_t verb=true)
static Int_t PdgBin(Int_t pdg)
unsigned short UShort_t
Definition: External.C:28
static TH3 * CopyH3(Container *parent, const char *name, const char *newName=0, Bool_t verb=true)
bool Bool_t
Definition: External.C:53
static TObject * CloneAndAdd(Container *c, TObject *o)
static TH1 * RatioH(const TH1 *num, const TH1 *denom, const char *name=0)
static TH1 * Make1D(Container *c, const TString &name, const TString &title, Color_t color, Style_t style, const TAxis &xAxis)
static TProfile * Make1P(Container *c, const TString &name, const TString &title, Color_t color, Style_t style, const TAxis &xAxis)
static void FixAxis(TAxis &axis, const char *title=0)
Definition: External.C:196
static TProfile * GetP(Container *parent, const char *name, Bool_t verb=true)
static void CopyAttr(const TH1 *src, TH1 *dest)