AliPhysics  026afea (026afea)
 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 <TDirectory.h>
23 # include <THashList.h>
24 # include <TProfile.h>
25 # include <TProfile2D.h>
26 # include <TDatabasePDG.h>
27 # include <THStack.h>
28 #else
29 class TList;
30 class TH1;
31 class TH2;
32 class TH3;
33 class TProfile;
34 class TProfile2D;
35 class TAxis;
36 class TDirectory;
37 class TDatabasePDG; // Auto-load
38 class THStack;
39 #endif
40 
254 {
255 public:
256  typedef TList Container;
264  virtual ~AliTrackletAODUtils() {}
265  //_________________________________________________________________________
279  static Bool_t CheckAxisNBins(const char* which,
280  const TAxis* a1,
281  const TAxis* a2);
291  static Bool_t CheckAxisLimits(const char* which,
292  const TAxis* a1,
293  const TAxis* a2);
304  static Bool_t CheckAxisBins(const char* which,
305  const TAxis* a1,
306  const TAxis* a2);
316  static Bool_t CheckAxisLabels(const char* which,
317  const TAxis* a1,
318  const TAxis* a2);
329  static Bool_t CheckAxis(const char* which,
330  const TAxis* a1,
331  const TAxis* a2,
332  Bool_t alsoLbls);
342  static Bool_t CheckConsistency(const TH1* h1, const TH1* h2);
343  /* @} */
344  //__________________________________________________________________
359  static TObject* GetO(Container* parent,
360  const char* name,
361  TClass* cls=0,
362  Bool_t verb=true);
373  static TObject* GetO(TDirectory* parent,
374  const char* name,
375  TClass* cls=0,
376  Bool_t verb=true);
386  static TH1* GetH1(Container* parent, const char* name, Bool_t verb=true);
396  static TH1* GetH1(TDirectory* parent, const char* name, Bool_t verb=true);
406  static TProfile* GetP1(Container* parent, const char* name, Bool_t verb=true);
416  static TH2* GetH2(Container* parent, const char* name, Bool_t verb=true);
426  static TH3* GetH3(Container* parent, const char* name, Bool_t verb=true);
436  static TProfile2D* GetP2(Container* parent, const char* name,
437  Bool_t verb=true);
447  static TProfile* GetP(Container* parent, const char* name, Bool_t verb=true)
448  {
449  return GetP1(parent, name, verb);
450  }
460  static Container* GetC(Container* parent, const char* name, Bool_t verb=true);
470  static Container* GetC(TDirectory* parent, const char* name,Bool_t verb=true);
480  static TDirectory* GetT(TDirectory* parent, const char* name,Bool_t verb=true);
491  static Double_t GetD(Container* parent, const char* name,
492  Double_t def=-1, Bool_t verb=true);
503  static Int_t GetI(Container* parent, const char* name,
504  Int_t def=-1, Bool_t verb=true);
515  static Int_t GetB(Container* parent, const char* name,
516  Bool_t def=false, Bool_t verb=true);
527  static TH1* CopyH1(Container* parent,
528  const char* name,
529  const char* newName=0,
530  Bool_t verb=true);
541  static TH2* CopyH2(Container* parent,
542  const char* name,
543  const char* newName=0,
544  Bool_t verb=true);
555  static TH3* CopyH3(Container* parent,
556  const char* name,
557  const char* newName=0,
558  Bool_t verb=true);
565  static void CopyAttr(const TH1* src, TH1* dest);
578  static TH1* Make1D(Container* c,
579  const TString& name,
580  const TString& title,
581  Color_t color,
582  Style_t style,
583  const TAxis& xAxis);
596  static TProfile* Make1P(Container* c,
597  const TString& name,
598  const TString& title,
599  Color_t color,
600  Style_t style,
601  const TAxis& xAxis);
615  static TH2* Make2D(Container* c,
616  const TString& name,
617  const TString& title,
618  Color_t color,
619  Style_t style,
620  const TAxis& xAxis,
621  const TAxis& yAxis);
636  static TH3* Make3D(Container* c,
637  const TString& name,
638  const TString& title,
639  Color_t color,
640  Style_t style,
641  const TAxis& xAxis,
642  const TAxis& yAxis,
643  const TAxis& zAxis);
657  static TProfile2D* Make2P(Container* c,
658  const TString& name,
659  const TString& title,
660  Color_t color,
661  Style_t style,
662  const TAxis& xAxis,
663  const TAxis& yAxis);
664 
674  static TH1* Scale(TH1* h, Double_t x, Double_t xe);
684  static TH2* Scale(TH2* h, Double_t x, Double_t xe);
694  static TH2* Scale(TH2* h, TH1* s);
701  static void FixAxis(TAxis& axis, const char* title=0);
709  static void ScaleAxis(TAxis& ret, Double_t fact=1);
717  static void SetAxis(TAxis& axis, Int_t n, Double_t* borders);
726  static void SetAxis(TAxis& axis, const TString& spec, const char* sep=":,");
735  static void SetAxis(TAxis& axis, Int_t n, Double_t l, Double_t h);
743  static void SetAxis(TAxis& axis, Int_t n, Double_t m);
751  static void PrintAxis(const TAxis& axis, Int_t nSig=2, const char* alt=0);
762  static TH2* ScaleToIPz(TH2* h, TH1* ipZ, Bool_t full=false);
763  static TH3* ScaleToIPz(TH3* h, TH1* ipZ, Bool_t full=false);
764  static TH3* ScaleDelta(TH3* h, TH2* scale);
789  static TH2* ProjectEtaDelta(TH3* h);
805  static TH1* ProjectDelta(TH2* h);
821  static TH1* ProjectDeltaFull(TH3* h);
843  static TH1* AverageOverIPz(TH2* h,
844  const char* name,
845  UShort_t mode,
846  TH1* ipz,
847  TH2* mask=0,
848  Bool_t verb=true);
857  static TObject* CloneAndAdd(Container* c, TObject* o);
870  static Double_t Integrate(TH1* h, Double_t min, Double_t max, Double_t& err);
882  static Double_t RatioE(Double_t n, Double_t en,
883  Double_t d, Double_t ed,
884  Double_t& er);
896  static Double_t RatioE2(Double_t n, Double_t e2n,
897  Double_t d, Double_t e2d,
898  Double_t& e2r);
908  static TH1* RatioH(const TH1* num, const TH1* denom, const char* name=0);
915  static void FixMinMax(TH1* h, Bool_t ignoreZero=true);
922  static void FixMinMax(THStack* h, Bool_t ignoreZero=true);
923 
924  /* @} */
933  static void PdgAttr(Int_t pdg, TString& nme, Color_t& c, Style_t& s);
938  static Int_t* PdgArray(Int_t& size);
944  static const TAxis& PdgAxis();
952  static Int_t PdgBin(Int_t pdg);
953  /* @} */
954 protected:
955  ClassDef(AliTrackletAODUtils,1); // Utilities
956 };
957 
958 //====================================================================
959 // Utilities
960 //____________________________________________________________________
962  const TAxis* a1,
963  const TAxis* a2)
964 {
965  if (a1->GetNbins() != a2->GetNbins()) {
966  ::Warning("CheckAxisNBins", "Incompatible number %s bins: %d vs %d",
967  which, a1->GetNbins(), a2->GetNbins());
968  return false;
969  }
970  return true;
971 }
972 //____________________________________________________________________
974  const TAxis* a1,
975  const TAxis* a2)
976 {
977  if (!TMath::AreEqualRel(a1->GetXmin(), a2->GetXmin(),1.E-12) ||
978  !TMath::AreEqualRel(a1->GetXmax(), a2->GetXmax(),1.E-12)) {
979  Warning("CheckAxisLimits",
980  "Limits of %s axis incompatible [%f,%f] vs [%f,%f]", which,
981  a1->GetXmin(), a1->GetXmax(), a2->GetXmin(), a2->GetXmax());
982  return false;
983  }
984  return true;
985 }
986 //____________________________________________________________________
988  const TAxis* a1,
989  const TAxis* a2)
990 {
991  const TArrayD * h1Array = a1->GetXbins();
992  const TArrayD * h2Array = a2->GetXbins();
993  Int_t fN = h1Array->fN;
994  if ( fN == 0 ) return true;
995  if (h2Array->fN != fN) {
996  // Redundant?
997  Warning("CheckAxisBins", "Not equal number of %s bin limits: %d vs %d",
998  which, fN, h2Array->fN);
999  return false;
1000  }
1001  else {
1002  for (int i = 0; i < fN; ++i) {
1003  if (!TMath::AreEqualRel(h1Array->GetAt(i),h2Array->GetAt(i),1E-10)) {
1004  Warning("CheckAxisBins",
1005  "%s limit # %3d incompatible: %f vs %f",
1006  which, i, h1Array->GetAt(i),h2Array->GetAt(i));
1007  return false;
1008  }
1009  }
1010  }
1011  return true;
1012 }
1013 //____________________________________________________________________
1015  const TAxis* a1,
1016  const TAxis* a2)
1017 {
1018  // check that axis have same labels
1019  THashList *l1 = (const_cast<TAxis*>(a1))->GetLabels();
1020  THashList *l2 = (const_cast<TAxis*>(a2))->GetLabels();
1021 
1022  if (!l1 && !l2) return true;
1023  if (!l1 || !l2) {
1024  Warning("CheckAxisLabels", "Difference in %s labels: %p vs %p",
1025  which, l1, l2);
1026  return false;
1027  }
1028  // check now labels sizes are the same
1029  if (l1->GetSize() != l2->GetSize()) {
1030  Warning("CheckAxisLabels", "Different number of %s labels: %d vs %d",
1031  which, l1->GetSize(), l2->GetSize());
1032  return false;
1033  }
1034  for (int i = 1; i <= a1->GetNbins(); ++i) {
1035  TString label1 = a1->GetBinLabel(i);
1036  TString label2 = a2->GetBinLabel(i);
1037  if (label1 != label2) {
1038  Warning("CheckAxisLabels", "%s label # %d not the same: '%s' vs '%s'",
1039  which, i, label1.Data(), label2.Data());
1040  return false;
1041  }
1042  }
1043 
1044  return true;
1045 }
1046 //____________________________________________________________________
1048  const TAxis* a1,
1049  const TAxis* a2,
1050  Bool_t alsoLbls)
1051 {
1052  if (!CheckAxisNBins (which, a1, a2)) return false;
1053  if (!CheckAxisLimits(which, a1, a2)) return false;
1054  if (!CheckAxisBins (which, a1, a2)) return false;
1055  if (alsoLbls && !CheckAxisLabels(which, a1, a2)) return false;
1056  return true;
1057 }
1058 //____________________________________________________________________
1060 {
1061  // Check histogram compatibility
1062  if (h1 == h2) return true;
1063 
1064  if (h1->GetDimension() != h2->GetDimension() ) {
1065  Warning("CheckConsistency",
1066  "%s and %s have different dimensions %d vs %d",
1067  h1->GetName(), h2->GetName(),
1068  h1->GetDimension(), h2->GetDimension());
1069  return false;
1070  }
1071  Int_t dim = h1->GetDimension();
1072 
1073  Bool_t alsoLbls = (h1->GetEntries() != 0 && h2->GetEntries() != 0);
1074  if (!CheckAxis("X", h1->GetXaxis(), h2->GetXaxis(), alsoLbls)) return false;
1075  if (dim > 1 &&
1076  !CheckAxis("Y", h1->GetYaxis(), h2->GetYaxis(), alsoLbls)) return false;
1077  if (dim > 2 &&
1078  !CheckAxis("Z", h1->GetZaxis(), h2->GetZaxis(), alsoLbls)) return false;
1079 
1080  return true;
1081 }
1082 
1083 //____________________________________________________________________
1085 {
1086  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
1087  Double_t dc = s->GetBinContent(i);
1088  Double_t de = s->GetBinError (i);
1089  Double_t dr = (dc > 1e-6 ? de/dc : 0);
1090  for (Int_t j = 1; j <= h->GetNbinsY(); j++) {
1091  Double_t nc = h->GetBinContent(i,j);
1092  Double_t ne = h->GetBinError (i,j);
1093  Double_t ns = (nc > 0 ? ne/nc : 0);
1094  Double_t sc = dc * nc;
1095  Double_t se = sc*TMath::Sqrt(ns*ns+dr*dr);
1096  // Printf("Setting bin %2d,%2d to %f +/- %f", sc, se);
1097  h->SetBinContent(i,j,sc);
1098  h->SetBinError (i,j,se);
1099  }
1100  }
1101  return h;
1102 }
1103 
1104 //____________________________________________________________________
1106 {
1107  Double_t rr = xe/x;
1108  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
1109  for (Int_t j = 1; j <= h->GetNbinsY(); j++) {
1110  Double_t c = h->GetBinContent(i,j);
1111  Double_t e = h->GetBinError (i,j);
1112  Double_t s = (c > 0 ? e/c : 0);
1113  Double_t sc = x * c;
1114  Double_t se = sc*TMath::Sqrt(s*s+rr*rr);
1115  // Printf("Setting bin %2d,%2d to %f +/- %f", sc, se);
1116  h->SetBinContent(i,j,sc);
1117  h->SetBinError (i,j,se);
1118  }
1119  }
1120  return h;
1121 }
1122 
1123 //____________________________________________________________________
1125 {
1126  Double_t rr = xe/x;
1127  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
1128  Double_t c = h->GetBinContent(i);
1129  Double_t e = h->GetBinError (i);
1130  Double_t s = (c > 0 ? e/c : 0);
1131  Double_t sc = x * c;
1132  Double_t se = sc*TMath::Sqrt(s*s+rr*rr);
1133  // Printf("Setting bin %2d,%2d to %f +/- %f", sc, se);
1134  h->SetBinContent(i,sc);
1135  h->SetBinError (i,se);
1136  }
1137  return h;
1138 }
1139 
1140 //____________________________________________________________________
1142  const char* name,
1143  TClass* cls,
1144  Bool_t verb)
1145 {
1146  if (!parent) {
1147  if (verb) ::Warning("GetO", "No parent container passed");
1148  return 0;
1149  }
1150  TObject* o = parent->FindObject(name);
1151  if (!o) {
1152  if (verb) ::Warning("GetO", "Object \"%s\" not found in \"%s\"",
1153  name, parent->GetName());
1154  // parent->ls();
1155  return 0;
1156  }
1157  if (!cls) return o;
1158  if (!o->IsA()->InheritsFrom(cls)) {
1159  if (verb) ::Warning("GetO", "\"%s\" is an object of class %s, not %s",
1160  name, o->ClassName(), cls->GetName());
1161  return 0;
1162  }
1163  return o;
1164 }
1165 //____________________________________________________________________
1167  const char* name,
1168  TClass* cls,
1169  Bool_t verb)
1170 {
1171  if (!parent) {
1172  if (!verb) ::Warning("GetO", "No parent directory passed");
1173  return 0;
1174  }
1175  TObject* o = parent->Get(name);
1176  if (!o) {
1177  if (verb) ::Warning("GetO", "Object \"%s\" not found in \"%s\"",
1178  name, parent->GetName());
1179  parent->ls();
1180  return 0;
1181  }
1182  if (!cls) return o;
1183  if (!o->IsA()->InheritsFrom(cls)) {
1184  if (verb) ::Warning("GetO", "\"%s\" is an object of class %s, not %s",
1185  name, o->ClassName(), cls->GetName());
1186  return 0;
1187  }
1188  return o;
1189 }
1190 //____________________________________________________________________
1191 TH1* AliTrackletAODUtils::GetH1(Container* parent, const char* name, Bool_t v)
1192 {
1193  return static_cast<TH1*>(GetO(parent, name, TH1::Class(), v));
1194 }
1195 //____________________________________________________________________
1196 TH1* AliTrackletAODUtils::GetH1(TDirectory* parent, const char* name, Bool_t v)
1197 {
1198  return static_cast<TH1*>(GetO(parent, name, TH1::Class(), v));
1199 }
1200 //____________________________________________________________________
1201 TH2* AliTrackletAODUtils::GetH2(Container* parent, const char* name, Bool_t v)
1202 {
1203  return static_cast<TH2*>(GetO(parent, name, TH2::Class(), v));
1204 }
1205 //____________________________________________________________________
1206 TH3* AliTrackletAODUtils::GetH3(Container* parent, const char* name, Bool_t v)
1207 {
1208  return static_cast<TH3*>(GetO(parent, name, TH3::Class(), v));
1209 }
1210 //____________________________________________________________________
1211 TProfile* AliTrackletAODUtils::GetP1(Container* parent,const char* name,
1212  Bool_t v)
1213 {
1214  return static_cast<TProfile*>(GetO(parent, name, TProfile::Class(), v));
1215 }
1216 //____________________________________________________________________
1217 TProfile2D* AliTrackletAODUtils::GetP2(Container* parent, const char* name,
1218  Bool_t v)
1219 {
1220  return static_cast<TProfile2D*>(GetO(parent, name, TProfile2D::Class(), v));
1221 }
1222 //____________________________________________________________________
1224 AliTrackletAODUtils::GetC(Container* parent, const char* name, Bool_t v)
1225 {
1226  return static_cast<Container*>(GetO(parent, name, Container::Class(), v));
1227 }
1228 //____________________________________________________________________
1230 AliTrackletAODUtils::GetC(TDirectory* parent, const char* name, Bool_t v)
1231 {
1232  return static_cast<Container*>(GetO(parent, name, Container::Class(), v));
1233 }
1234 //____________________________________________________________________
1235 TDirectory*
1236 AliTrackletAODUtils::GetT(TDirectory* parent, const char* name, Bool_t v)
1237 {
1238  TDirectory* d = parent->GetDirectory(name);
1239  if (!d) {
1240  if (v) ::Warning("GetO", "Directory \"%s\" not found in \"%s\"",
1241  name, parent->GetName());
1242  return 0;
1243  }
1244  return d;
1245 }
1246 //____________________________________________________________________
1248  const char* name,
1249  Double_t def,
1250  Bool_t v)
1251 {
1252  TParameter<double>* p =
1253  static_cast<TParameter<double>*>(GetO(parent, name,
1255  if (!p) return def;
1256  return p->GetVal();
1257 }
1258 //____________________________________________________________________
1260  const char* name,
1261  Int_t def,
1262  Bool_t v)
1263 {
1264  TParameter<int>* p =
1265  static_cast<TParameter<int>*>(GetO(parent, name,
1266  TParameter<int>::Class(),v));
1267  if (!p) return def;
1268  return p->GetVal();
1269 }
1270 //____________________________________________________________________
1272  const char* name,
1273  Bool_t def,
1274  Bool_t v)
1275 {
1276  TParameter<bool>* p =
1277  static_cast<TParameter<bool>*>(GetO(parent, name,
1278  TParameter<bool>::Class(), v));
1279  if (!p) return def;
1280  return p->GetVal();
1281 }
1282 //____________________________________________________________________
1284  const char* name,
1285  const char* newName,
1286  Bool_t v)
1287 {
1288  TH1* orig = GetH1(parent, name, v);
1289  if (!orig) return 0;
1290  TH1* ret = static_cast<TH1*>(orig->Clone(newName ? newName : name));
1291  ret->SetDirectory(0); // Release from file container
1292  return ret;
1293 }
1294 //____________________________________________________________________
1296  const char* name,
1297  const char* newName,
1298  Bool_t v)
1299 {
1300  TH2* orig = GetH2(parent, name, v);
1301  if (!orig) return 0;
1302  TH2* ret = static_cast<TH2*>(orig->Clone(newName ? newName : name));
1303  ret->SetDirectory(0); // Release from file container
1304  return ret;
1305 }
1306 //____________________________________________________________________
1308  const char* name,
1309  const char* newName,
1310  Bool_t v)
1311 {
1312  TH3* orig = GetH3(parent, name, v);
1313  if (!orig) return 0;
1314  TH3* ret = static_cast<TH3*>(orig->Clone(newName ? newName : name));
1315  ret->SetDirectory(0); // Release from file container
1316  return ret;
1317 }
1318 
1319 //____________________________________________________________________
1320 void AliTrackletAODUtils::CopyAttr(const TH1* src, TH1* dest)
1321 {
1322  if (!src || !dest) return;
1323  dest->SetMarkerStyle(src->GetMarkerStyle());
1324  dest->SetMarkerColor(src->GetMarkerColor());
1325  dest->SetMarkerSize (src->GetMarkerSize());
1326  dest->SetLineStyle (src->GetLineStyle());
1327  dest->SetLineColor (src->GetLineColor());
1328  dest->SetLineWidth (src->GetLineWidth());
1329  dest->SetFillStyle (src->GetFillStyle());
1330  dest->SetFillColor (src->GetFillColor());
1331 }
1332 
1333 
1334 //____________________________________________________________________
1336  const TString& name,
1337  const TString& title,
1338  Color_t color,
1339  Style_t style,
1340  const TAxis& xAxis)
1341 {
1342  TString n = name;
1343  TString t = title;
1344  TH1* ret = 0;
1345  Int_t nx = xAxis.GetNbins();
1346  if (t.IsNull()) t = xAxis.GetTitle();
1347  if (xAxis.GetXbins() && xAxis.GetXbins()->GetArray())
1348  ret = new TH1D(n,t,nx,xAxis.GetXbins()->GetArray());
1349  else
1350  ret = new TH1D(n,t,nx,xAxis.GetXmin(),xAxis.GetXmax());
1351  ret->Sumw2();
1352  ret->SetXTitle(xAxis.GetTitle());
1353  static_cast<const TAttAxis&>(xAxis).Copy(*(ret->GetXaxis()));
1354  ret->SetDirectory(0);
1355  ret->SetLineColor(color);
1356  ret->SetMarkerColor(color);
1357  ret->SetFillColor(kWhite);// color);
1358  ret->SetFillStyle(0);
1359  ret->SetMarkerStyle(style);
1360  if (const_cast<TAxis&>(xAxis).GetLabels()) {
1361  for (Int_t i = 1; i <= xAxis.GetNbins(); i++)
1362  ret->GetXaxis()->SetBinLabel(i, xAxis.GetBinLabel(i));
1363  }
1364  switch (style) {
1365  case 27:
1366  case 33:
1367  ret->SetMarkerSize(1.4);
1368  break;
1369  case 29:
1370  case 30:
1371  ret->SetMarkerSize(1.2);
1372  break;
1373  }
1374  if (c) c->Add(ret);
1375  return ret;
1376 }
1377 //____________________________________________________________________
1379  const TString& name,
1380  const TString& title,
1381  Color_t color,
1382  Style_t style,
1383  const TAxis& xAxis)
1384 {
1385  TString n = name;
1386  TString t = title;
1387  TProfile* ret = 0;
1388  Int_t nx = xAxis.GetNbins();
1389  if (t.IsNull()) t = xAxis.GetTitle();
1390  if (xAxis.GetXbins() && xAxis.GetXbins()->GetArray())
1391  ret = new TProfile(n,t,nx,xAxis.GetXbins()->GetArray());
1392  else
1393  ret = new TProfile(n,t,nx,xAxis.GetXmin(),xAxis.GetXmax());
1394  ret->Sumw2();
1395  ret->SetXTitle(xAxis.GetTitle());
1396  static_cast<const TAttAxis&>(xAxis).Copy(*(ret->GetXaxis()));
1397  ret->SetDirectory(0);
1398  ret->SetLineColor(color);
1399  ret->SetMarkerColor(color);
1400  ret->SetFillColor(kWhite);// color);
1401  ret->SetFillStyle(0);
1402  ret->SetMarkerStyle(style);
1403  if (const_cast<TAxis&>(xAxis).GetLabels()) {
1404  for (Int_t i = 1; i <= xAxis.GetNbins(); i++)
1405  ret->GetXaxis()->SetBinLabel(i, xAxis.GetBinLabel(i));
1406  }
1407  switch (style) {
1408  case 27:
1409  case 33:
1410  ret->SetMarkerSize(1.4);
1411  break;
1412  case 29:
1413  case 30:
1414  ret->SetMarkerSize(1.2);
1415  break;
1416  }
1417  if (c) c->Add(ret);
1418  return ret;
1419 }
1420 //____________________________________________________________________
1422  const TString& name,
1423  const TString& title,
1424  Color_t color,
1425  Style_t style,
1426  const TAxis& xAxis,
1427  const TAxis& yAxis)
1428 {
1429  TString n = name;
1430  TString t = title;
1431  TH2* ret = 0;
1432  Int_t nx = xAxis.GetNbins();
1433  Int_t ny = yAxis.GetNbins();
1434  const Double_t* xb = (xAxis.GetXbins() && xAxis.GetXbins()->GetArray() ?
1435  xAxis.GetXbins()->GetArray() : 0);
1436  const Double_t* yb = (yAxis.GetXbins() && yAxis.GetXbins()->GetArray() ?
1437  yAxis.GetXbins()->GetArray() : 0);
1438  if (t.IsNull())
1439  t.Form("%s vs %s", yAxis.GetTitle(), xAxis.GetTitle());
1440  if (xb) {
1441  if (yb) ret = new TH2D(n,t,nx,xb,ny,yb);
1442  else ret = new TH2D(n,t,
1443  nx,xb,
1444  ny,yAxis.GetXmin(),yAxis.GetXmax());
1445  }
1446  else {
1447  if (yb) ret = new TH2D(n,t,
1448  nx,xAxis.GetXmin(),xAxis.GetXmax(),
1449  ny,yb);
1450  else ret = new TH2D(n,t,
1451  nx,xAxis.GetXmin(),xAxis.GetXmax(),
1452  ny,yAxis.GetXmin(),yAxis.GetXmax());
1453  }
1454  ret->Sumw2();
1455  ret->SetXTitle(xAxis.GetTitle());
1456  ret->SetYTitle(yAxis.GetTitle());
1457  ret->SetLineColor(color);
1458  ret->SetMarkerColor(color);
1459  ret->SetFillColor(color);
1460  ret->SetMarkerStyle(style);
1461  static_cast<const TAttAxis&>(xAxis).Copy(*(ret->GetXaxis()));
1462  static_cast<const TAttAxis&>(yAxis).Copy(*(ret->GetYaxis()));
1463  ret->SetDirectory(0);
1464  if (const_cast<TAxis&>(xAxis).GetLabels()) {
1465  for (Int_t i = 1; i <= xAxis.GetNbins(); i++)
1466  ret->GetXaxis()->SetBinLabel(i, xAxis.GetBinLabel(i));
1467  }
1468  if (const_cast<TAxis&>(yAxis).GetLabels()) {
1469  for (Int_t i = 1; i <= yAxis.GetNbins(); i++)
1470  ret->GetYaxis()->SetBinLabel(i, yAxis.GetBinLabel(i));
1471  }
1472  if (c) c->Add(ret);
1473  return ret;
1474 }
1475 //____________________________________________________________________
1477  const TString& name,
1478  const TString& title,
1479  Color_t color,
1480  Style_t style,
1481  const TAxis& xAxis,
1482  const TAxis& yAxis,
1483  const TAxis& zAxis)
1484 {
1485  TString n = name;
1486  TString t = title;
1487  TH3* ret = 0;
1488  Int_t nx = xAxis.GetNbins();
1489  Int_t ny = yAxis.GetNbins();
1490  Int_t nz = zAxis.GetNbins();
1491  Double_t* xb = (xAxis.GetXbins() && xAxis.GetXbins()->GetArray() ?
1492  const_cast<Double_t*>(xAxis.GetXbins()->GetArray()) : 0);
1493  Double_t* yb = (yAxis.GetXbins() && yAxis.GetXbins()->GetArray() ?
1494  const_cast<Double_t*>(yAxis.GetXbins()->GetArray()) : 0);
1495  Double_t* zb = (zAxis.GetXbins() && zAxis.GetXbins()->GetArray() ?
1496  const_cast<Double_t*>(zAxis.GetXbins()->GetArray()) : 0);
1497  if (t.IsNull())
1498  t.Form("%s vs %s vs %s",
1499  zAxis.GetTitle(), yAxis.GetTitle(), xAxis.GetTitle());
1500  if (xb || yb || zb) {
1501  // One or more axis are defined as arrays. Make sure the rest are
1502  // also arrays.
1503  if (xb) {
1504  xb = new Double_t[nx+1];
1505  Double_t dx = (xAxis.GetXmax()-xAxis.GetXmin())/nx;
1506  xb[0] = xAxis.GetXmin();
1507  for (Int_t i = 1; i <= nx; i++) xb[i] = xb[i-1]+dx;
1508  }
1509  if (yb) {
1510  yb = new Double_t[ny+1];
1511  Double_t dy = (yAxis.GetXmax()-yAxis.GetXmin())/ny;
1512  yb[0] = yAxis.GetXmin();
1513  for (Int_t i = 1; i <= ny; i++) yb[i] = yb[i-1]+dy;
1514  }
1515  if (zb) {
1516  zb = new Double_t[nz+1];
1517  Double_t dz = (zAxis.GetXmax()-zAxis.GetXmin())/nz;
1518  zb[0] = zAxis.GetXmin();
1519  for (Int_t i = 1; i <= nz; i++) zb[i] = zb[i-1]+dz;
1520  }
1521  ret = new TH3D(n,t,nx,xb,ny,yb,nz,zb);
1522  }
1523  else {
1524  ret = new TH3D(n,t,
1525  nx, xAxis.GetXmin(), xAxis.GetXmax(),
1526  ny, yAxis.GetXmin(), yAxis.GetXmax(),
1527  nz, zAxis.GetXmin(), zAxis.GetXmax());
1528  }
1529  ret->Sumw2();
1530  ret->SetXTitle(xAxis.GetTitle());
1531  ret->SetYTitle(yAxis.GetTitle());
1532  ret->SetZTitle(zAxis.GetTitle());
1533  ret->SetLineColor(color);
1534  ret->SetMarkerColor(color);
1535  ret->SetFillColor(color);
1536  ret->SetMarkerStyle(style);
1537  static_cast<const TAttAxis&>(xAxis).Copy(*(ret->GetXaxis()));
1538  static_cast<const TAttAxis&>(yAxis).Copy(*(ret->GetYaxis()));
1539  static_cast<const TAttAxis&>(zAxis).Copy(*(ret->GetZaxis()));
1540  ret->SetDirectory(0);
1541  if (const_cast<TAxis&>(xAxis).GetLabels()) {
1542  for (Int_t i = 1; i <= xAxis.GetNbins(); i++)
1543  ret->GetXaxis()->SetBinLabel(i, xAxis.GetBinLabel(i));
1544  }
1545  if (const_cast<TAxis&>(yAxis).GetLabels()) {
1546  for (Int_t i = 1; i <= yAxis.GetNbins(); i++)
1547  ret->GetYaxis()->SetBinLabel(i, yAxis.GetBinLabel(i));
1548  }
1549  if (const_cast<TAxis&>(zAxis).GetLabels()) {
1550  for (Int_t i = 1; i <= zAxis.GetNbins(); i++)
1551  ret->GetZaxis()->SetBinLabel(i, zAxis.GetBinLabel(i));
1552  }
1553  if (c) c->Add(ret);
1554  return ret;
1555 }
1556 //____________________________________________________________________
1558  const TString& name,
1559  const TString& title,
1560  Color_t color,
1561  Style_t style,
1562  const TAxis& xAxis,
1563  const TAxis& yAxis)
1564 {
1565  TString n = name;
1566  TString t = title;
1567  TProfile2D* ret = 0;
1568  Int_t nx = xAxis.GetNbins();
1569  Int_t ny = yAxis.GetNbins();
1570  const Double_t* xb = (xAxis.GetXbins() && xAxis.GetXbins()->GetArray() ?
1571  xAxis.GetXbins()->GetArray() : 0);
1572  const Double_t* yb = (yAxis.GetXbins() && yAxis.GetXbins()->GetArray() ?
1573  yAxis.GetXbins()->GetArray() : 0);
1574  if (t.IsNull())
1575  t.Form("%s vs %s", yAxis.GetTitle(), xAxis.GetTitle());
1576  if (xb) {
1577  if (yb) ret = new TProfile2D(n,t,nx,xb,ny,yb);
1578  else ret = new TProfile2D(n,t,
1579  nx,xb,
1580  ny,yAxis.GetXmin(),yAxis.GetXmax());
1581  }
1582  else {
1583  if (yb) ret = new TProfile2D(n,t,
1584  nx,xAxis.GetXmin(),xAxis.GetXmax(),
1585  ny,yb);
1586  else ret = new TProfile2D(n,t,
1587  nx,xAxis.GetXmin(),xAxis.GetXmax(),
1588  ny,yAxis.GetXmin(),yAxis.GetXmax());
1589  }
1590  ret->Sumw2();
1591  ret->SetXTitle(xAxis.GetTitle());
1592  ret->SetYTitle(yAxis.GetTitle());
1593  ret->SetLineColor(color);
1594  ret->SetMarkerColor(color);
1595  ret->SetFillColor(color);
1596  ret->SetMarkerStyle(style);
1597  static_cast<const TAttAxis&>(xAxis).Copy(*(ret->GetXaxis()));
1598  static_cast<const TAttAxis&>(yAxis).Copy(*(ret->GetYaxis()));
1599  ret->SetDirectory(0);
1600  if (const_cast<TAxis&>(xAxis).GetLabels()) {
1601  for (Int_t i = 1; i <= xAxis.GetNbins(); i++)
1602  ret->GetXaxis()->SetBinLabel(i, xAxis.GetBinLabel(i));
1603  }
1604  if (const_cast<TAxis&>(yAxis).GetLabels()) {
1605  for (Int_t i = 1; i <= yAxis.GetNbins(); i++)
1606  ret->GetYaxis()->SetBinLabel(i, yAxis.GetBinLabel(i));
1607  }
1608  if (c) c->Add(ret);
1609  return ret;
1610 }
1611 //____________________________________________________________________
1612 void AliTrackletAODUtils::FixAxis(TAxis& axis, const char* title)
1613 {
1614  if (title && title[0] != '\0') axis.SetTitle(title);
1615  axis. SetNdivisions(210);
1616  axis. SetLabelFont(42);
1617  axis. SetLabelSize(0.03);
1618  axis. SetLabelOffset(0.005);
1619  axis. SetLabelColor(kBlack);
1620  axis. SetTitleOffset(1);
1621  axis. SetTitleFont(42);
1622  axis. SetTitleSize(0.04);
1623  axis. SetTitleColor(kBlack);
1624  axis. SetTickLength(0.03);
1625  axis. SetAxisColor(kBlack);
1626 }
1627 //____________________________________________________________________
1629  Double_t fact)
1630 {
1631  if (ret.GetXbins()->GetArray()) {
1632  TArrayD bins(*ret.GetXbins());
1633  for (Int_t i = 0; i < bins.GetSize(); i++) bins[i] *= fact;
1634  ret.Set(ret.GetNbins(), bins.GetArray());
1635  }
1636  else {
1637  ret.Set(ret.GetNbins(), fact*ret.GetXmin(), fact*ret.GetXmax());
1638  }
1639  FixAxis(ret);
1640 }
1641 //____________________________________________________________________
1643 {
1644  axis.Set(n, borders);
1645  FixAxis(axis);
1646 }
1647 //____________________________________________________________________
1649  const TString& spec,
1650  const char* sep)
1651 {
1652  TString s(spec);
1653  Bool_t isRange = false, isUnit = false;
1654  if (s[0] == 'r' || s[0] == 'R') {
1655  isRange = true;
1656  s.Remove(0,1);
1657  }
1658  if (s[0] == 'u' || s[0] == 'U') {
1659  isUnit = true;
1660  s.Remove(0,1);
1661  }
1662  TObjArray* tokens = s.Tokenize(sep);
1663  TArrayD bins(tokens->GetEntries());
1664  TObjString* token = 0;
1665  TIter next(tokens);
1666  Int_t i = 0;
1667  while ((token = static_cast<TObjString*>(next()))) {
1668  Double_t v = token->String().Atof();
1669  bins[i] = v;
1670  i++;
1671  }
1672  delete tokens;
1673  if (isUnit) {
1674  if (bins.GetSize() > 1)
1675  SetAxis(axis, Int_t(bins[1]-bins[0]), bins[0], bins[1]);
1676  else
1677  SetAxis(axis, 2*Int_t(bins[0]), bins[0]);
1678  }
1679  else if (isRange) {
1680  Int_t nBins = Int_t(bins[0]);
1681  if (bins.GetSize() > 2)
1682  SetAxis(axis, nBins, bins[1], bins[2]);
1683  else
1684  SetAxis(axis, nBins, bins[1]);
1685  }
1686  else
1687  SetAxis(axis, bins.GetSize()-1,bins.GetArray());
1688 }
1689 //____________________________________________________________________
1691  Int_t n,
1692  Double_t l,
1693  Double_t h)
1694 {
1695  axis.Set(n, l, h);
1696  FixAxis(axis);
1697 }
1698 //____________________________________________________________________
1700 {
1701  SetAxis(axis, n, -TMath::Abs(m), +TMath::Abs(m));
1702 }
1703 //____________________________________________________________________
1705  Int_t nSig,
1706  const char* alt)
1707 {
1708  printf(" %17s axis: ", alt ? alt : axis.GetTitle());
1709  if (axis.GetXbins() && axis.GetXbins()->GetArray()) {
1710  printf("%.*f", nSig, axis.GetBinLowEdge(1));
1711  for (Int_t i = 1; i <= axis.GetNbins(); i++)
1712  printf(":%.*f", nSig, axis.GetBinUpEdge(i));
1713  }
1714  else
1715  printf("%d bins between %.*f and %.*f",
1716  axis.GetNbins(), nSig, axis.GetXmin(),nSig,axis.GetXmax());
1717  printf("\n");
1718 }
1719 //____________________________________________________________________
1721 {
1722  if (!h) {
1723  ::Warning("ScaleToIPz","Nothing to scale");
1724  return 0;
1725  }
1726  if (!ipZ) {
1727  ::Warning("ScaleToIPz","Nothing to scale by");
1728  return 0;
1729  }
1730  TH2* ret = static_cast<TH2*>(h->Clone());
1731  ret->SetDirectory(0);
1732  if (!ipZ) return ret;
1733  for (Int_t iy = 1; iy <= ret->GetNbinsY(); iy++) {
1734  Double_t z = ret->GetYaxis()->GetBinCenter(iy);
1735  Int_t bin = ipZ->GetXaxis()->FindBin(z);
1736  Double_t nEv = ipZ->GetBinContent(bin);
1737  Double_t eEv = ipZ->GetBinError (bin);
1738  Double_t esc = (nEv > 0 ? 1./nEv : 0);
1739  Double_t rE2 = esc*esc*eEv*eEv;
1740  for (Int_t ix = 1; ix <= ret->GetNbinsX(); ix++) {
1741  Double_t c = ret->GetBinContent(ix,iy);
1742  Double_t e = ret->GetBinError (ix,iy);
1743  Double_t r = (c > 0 ? e/c : 0);
1744  // Scale by number of events, and error propagate
1745  Double_t sc = c * esc;
1746  Double_t se = 0;
1747  if (full) se = sc * TMath::Sqrt(r*r+rE2);
1748  else se = e * esc;
1749  Double_t scl = 1 / ret->GetXaxis()->GetBinWidth(ix);
1750  ret->SetBinContent(ix, iy, scl*sc);
1751  ret->SetBinError (ix, iy, scl*se);
1752  }
1753  }
1754  return ret;
1755 }
1756 //____________________________________________________________________
1758 {
1759  if (!h) {
1760  ::Warning("ScaleToIPz","Nothing to scale");
1761  return 0;
1762  }
1763  if (!ipZ) {
1764  ::Warning("ScaleToIPz","Nothing to scale by");
1765  return 0;
1766  }
1767  TH3* ret = static_cast<TH3*>(h->Clone());
1768  ret->SetDirectory(0);
1769  if (!ipZ) return ret;
1770  for (Int_t iz = 1; iz <= ret->GetNbinsZ(); iz++) {
1771  Double_t z = ret->GetZaxis()->GetBinCenter(iz);
1772  Int_t bin = ipZ->GetXaxis()->FindBin(z);
1773  Double_t nEv = ipZ->GetBinContent(bin);
1774  Double_t eEv = ipZ->GetBinError (bin);
1775  Double_t esc = (nEv > 0 ? 1./nEv : 0);
1776  Double_t rE2 = esc*esc*eEv*eEv;
1777  for (Int_t iy = 1; iy <= ret->GetNbinsY(); iy++) {
1778  for (Int_t ix = 1; ix <= ret->GetNbinsX(); ix++) {
1779  Double_t c = ret->GetBinContent(ix,iy,iz);
1780  Double_t e = ret->GetBinError (ix,iy,iz);
1781  Double_t r = (c > 0 ? e/c : 0);
1782  // Scale by number of events, and error propagate
1783  Double_t sc = c * esc;
1784  Double_t se = 0;
1785  if (full) se = sc * TMath::Sqrt(r*r+rE2);
1786  else se = e * esc;
1787  Double_t scl = 1 / ret->GetXaxis()->GetBinWidth(ix);
1788  ret->SetBinContent(ix, iy, iz, scl*sc);
1789  ret->SetBinError (ix, iy, iz, scl*se);
1790  }
1791  }
1792  }
1793  return ret;
1794 }
1795 //____________________________________________________________________
1797 {
1798  for (Int_t i = 1; i <= h->GetNbinsX(); i++) { // eta
1799  for (Int_t j = 1; j <= h->GetNbinsZ(); j++) { // IPz
1800  Double_t scale = etaIPzScale->GetBinContent(i,j);
1801  Double_t scaleE = etaIPzScale->GetBinError (i,j);
1802  Double_t q = (scale > 0 ? scaleE / scale : 0);
1803  for (Int_t k = 1; k <= h->GetNbinsY(); k++) { // delta
1804  Double_t c = h->GetBinContent(i,k,j);
1805  Double_t e = h->GetBinError (i,k,j);
1806  Double_t r = (c > 0 ? e / c : 0);
1807  Double_t w = c * scale;
1808  Double_t v = w * TMath::Sqrt(r*r + q*q);
1809  h->SetBinContent(i,k,j,w);
1810  h->SetBinError (i,k,j,v);
1811 #if 0
1812  Printf("%2d,%3d,%2d=(%9g+/-%9g)*(%9g+/-%9g)=(%9g+/-%9g)",
1813  i,k,j,c,e,scale,scaleE,w,v);
1814 #endif
1815  }
1816  }
1817  }
1818  return h;
1819 }
1820 //____________________________________________________________________
1821 // @cond
1822 TH2* AliTrackletAODUtils::ProjectEtaDelta(TH3* h)
1823 {
1824  TH2* etaDelta = static_cast<TH2*>(h->Project3D("yx e"));
1825  etaDelta->SetName("etaDelta");
1826  etaDelta->SetTitle(h->GetTitle());
1827  etaDelta->SetDirectory(0);
1828  etaDelta->SetZTitle("d^{2}#it{N}/(d#etad#Delta)");
1829 #if 1
1830  // Reset content of projected histogram and calculate averages
1831  etaDelta->Reset();
1832  for (Int_t i = 1; i <= h->GetNbinsX(); i++) { // Loop over eta
1833  for (Int_t j = 1; j <= h->GetNbinsY(); j++) { // Loop over Delta
1834  Double_t sum = 0;
1835  Double_t sumw = 0;
1836  Int_t cnt = 0;
1837  for (Int_t k = 1; k <= h->GetNbinsZ(); k++) { // Loop over IPz
1838  Double_t c = h->GetBinContent(i,j,k);
1839  Double_t e = h->GetBinError (i,j,k);
1840  if (c < 1e-6 || e/c > 1) continue;
1841 #if 0
1842  Double_t w = 1/e/e;
1843  sum += w * c;
1844  sumw += w;
1845 #else
1846  sum += c;
1847  sumw += e*e;
1848 #endif
1849  cnt += 1;
1850  }
1851  if (sumw < 1e-6) continue;
1852 #if 0
1853  etaDelta->SetBinContent(i,j,sum/sumw);
1854  etaDelta->SetBinError (i,j,TMath::Sqrt(1/sumw));
1855 #else
1856  etaDelta->SetBinContent(i,j,sum/cnt);
1857  etaDelta->SetBinError (i,j,TMath::Sqrt(sumw)/cnt);
1858 #endif
1859  }
1860  }
1861 #endif
1862  return etaDelta;
1863 }
1864 // @endcond
1865 //____________________________________________________________________
1867 {
1868  if (!h) {
1869  Warning("ProjectDelta", "No histogram passed");
1870  return 0;
1871  }
1872  TH1* delta = h->ProjectionY("delta");
1873  delta->SetDirectory(0);
1874  delta->SetTitle(h->GetTitle());
1875  delta->SetYTitle("d#it{N}/d#Delta");
1876  delta->Scale(1./h->GetNbinsX());
1877  return delta;
1878 }
1879 //____________________________________________________________________
1881 {
1882  if (!h) {
1883  Warning("ProjectDelta", "No histogram passed");
1884  return 0;
1885  }
1886  TH2* tmp = ProjectEtaDelta(h);
1887  TH1* delta = ProjectDelta(tmp);
1888  delta->SetDirectory(0);
1889  tmp->SetDirectory(0);
1890  delete tmp;
1891  return delta;
1892 }
1893 
1894 //____________________________________________________________________
1896  const char* name,
1897  UShort_t mode,
1898  TH1* ipz,
1899  TH2* other,
1900  Bool_t verb)
1901 {
1902  if (!h) return 0;
1903 
1904  Int_t nIPz = h->GetNbinsY();
1905  Int_t nEta = h->GetNbinsX();
1906  TH1* p = h->ProjectionX(name,1,nIPz,"e");
1907  TH2* mask = (other ? other : h);
1908  p->SetDirectory(0);
1909  p->SetFillColor(0);
1910  p->SetFillStyle(0);
1911  p->SetYTitle(Form("#LT%s#GT", h->GetYaxis()->GetTitle()));
1912  p->Reset();
1913  for (Int_t etaBin = 1; etaBin <= nEta; etaBin++) {
1914  TArrayD hv(nIPz);
1915  TArrayD he(nIPz);
1916  TArrayD hr(nIPz);
1917  TArrayI hb(nIPz);
1918  Int_t j = 0;
1919  for (Int_t ipzBin = 1; ipzBin <= nIPz; ipzBin++) {
1920  Double_t bc = mask->GetBinContent(etaBin, ipzBin);
1921  if (bc < 1e-9) continue; // Low value
1922  Double_t be = mask->GetBinError (etaBin, ipzBin);
1923  if (TMath::IsNaN(be)) continue; // Bad error value
1924  hv[j] = bc;
1925  he[j] = be;
1926  hr[j] = be/bc;
1927  hb[j] = ipzBin;
1928  j++;
1929  }
1930  // Now we have all interesting values. Sort the relative error
1931  // array to get the most significant first. Note, we sort on the
1932  // mask errors which may not be the same as the histogram errors.
1933  TArrayI idx(nIPz);
1934  TMath::Sort(j, hr.fArray, idx.fArray, false);
1935  Double_t nsum = 0; // Running weighted sum
1936  Double_t nsume = 0; // Running weighted sum error
1937  Double_t dsum = 0;
1938  Double_t dsume = 0;
1939  Int_t n = 0;
1940  Double_t rat = 1e99;
1941 
1942  Int_t k = 0;
1943  for (k = 0; k < j; k++) {
1944  Int_t l = idx[k]; // Sorted index - ipzBin
1945  Int_t ipzBin = hb[l];
1946  Double_t hvv = hv[l];
1947  Double_t hee = he[l];
1948  Double_t x = TMath::Sqrt(nsume+hee*hee)/(nsum+hvv);
1949  if (x > rat) {
1950  continue; // Ignore - does not help
1951  }
1952 
1953  Double_t by = mask->GetYaxis()->GetBinCenter(ipzBin);
1954  Int_t ib = ipz ? ipz->FindBin(by) : 0;
1955  rat = x;
1956  nsum += h->GetBinContent(etaBin, ipzBin);
1957  nsume += TMath::Power(h->GetBinError(etaBin, ipzBin), 2);
1958  // If we do not have the vertex distribution, then just count
1959  // number of observations.
1960  dsum += !ipz ? 1 : ipz->GetBinContent(ib);
1961  dsume += !ipz ? 0 : TMath::Power(ipz->GetBinError(ib),2);
1962  n++;
1963  }
1964  if (k == 0 || n == 0) {
1965  if (verb)
1966  ::Warning("Average", "Eta bin # %3d has no data",etaBin);
1967  continue; // This eta bin empty!
1968  }
1969  Double_t norm = (mode > 0 ? n : dsum);
1970  Double_t rne = nsume/nsum/nsum;
1971  Double_t rde = dsume/dsum/dsum;
1972  Double_t av = nsum/norm;
1973  Double_t ave = 0;
1974  if (mode==2) ave = TMath::Sqrt(nsume)/n;
1975  else ave = av*TMath::Sqrt(rne+rde);
1976  p->SetBinContent(etaBin, av);
1977  p->SetBinError (etaBin, ave);
1978  }
1979  if (mode == 0) p->Scale(1, "width");
1980  return p;
1981 }
1982 
1983 //____________________________________________________________________
1985 {
1986  if (!o) {
1987  ::Warning("CloneAndAdd", "Nothing to clone");
1988  return 0;
1989  }
1990  TObject* copy = o->Clone();
1991  if (copy->IsA()->InheritsFrom(TH1::Class()))
1992  // Release from underlying directory
1993  static_cast<TH1*>(copy)->SetDirectory(0);
1994  if (c)
1995  // Add to output container
1996  c->Add(copy);
1997  return copy;
1998 }
1999 //____________________________________________________________________
2001  Double_t min,
2002  Double_t max,
2003  Double_t& err)
2004 {
2005  const Double_t epsilon = 1e-6;
2006  Int_t bMin = h->FindBin(min+epsilon);
2007  Int_t bMax = h->FindBin(max-epsilon);
2008  if (bMin < 1) bMin = 0; // Include underflow bin
2009  Double_t val = h->IntegralAndError(bMin, bMax, err);
2010  // For a-symmetric histograms, return
2011  if (TMath::Abs(h->GetXaxis()->GetXmin()+h->GetXaxis()->GetXmax())>=epsilon)
2012  return val;
2013 
2014  // Otherwise, also integrate negative side
2015  Double_t err2;
2016  bMin = h->FindBin(-min+epsilon);
2017  bMax = h->FindBin(-max-epsilon);
2018  val += h->IntegralAndError(bMin, bMax, err2);
2019  err = TMath::Sqrt(err*err+err2*err2);
2020  return val;
2021 }
2022 //____________________________________________________________________
2024  Double_t d, Double_t ed,
2025  Double_t& er)
2026 {
2027  Double_t r = 0;
2028  er = 0;
2029  if (TMath::Abs(n) < 1e-16 || TMath::Abs(d) < 1e-16) return 0;
2030  r = n/d;
2031  er = TMath::Sqrt(en*en/n/n + ed*ed/d/d);
2032  return r;
2033 }
2034 //____________________________________________________________________
2036  Double_t d, Double_t e2d,
2037  Double_t& e2r)
2038 {
2039  Double_t r = 0;
2040  e2r = 0;
2041  if (TMath::Abs(n) < 1e-16 || TMath::Abs(d) < 1e-16) return 0;
2042  r = n/d;
2043  e2r = (e2n/n/n + e2d/d/d);
2044  return r;
2045 }
2046 //____________________________________________________________________
2048  const TH1* denom,
2049  const char* name)
2050 {
2051  if (!num || !denom) {
2052  ::Warning("RatioH", "Numerator (%p) or denominator (%p) missing",
2053  num, denom);
2054  return 0;
2055  }
2056  if (!CheckConsistency(num, denom)) return 0;
2057 
2058  TH1* ratio = static_cast<TH1*>(num->Clone(name ? name : num->GetName()));
2059  ratio->SetDirectory(0);
2060  ratio->SetMinimum(-1111);
2061  ratio->SetMaximum(-1111);
2062  ratio->GetListOfFunctions()->Clear();
2063  ratio->Divide(denom);
2064  ratio->SetYTitle(Form("%s / %s",
2065  num ->GetYaxis()->GetTitle(),
2066  denom->GetYaxis()->GetTitle()));
2067  ratio->SetTitle(Form("%s / %s", num ->GetTitle(), denom->GetTitle()));
2068  FixMinMax(ratio, true);
2069 
2070  return ratio;
2071 }
2072 //____________________________________________________________________
2074 {
2075  Double_t min = +1e99;
2076  Double_t max = -1e99;
2077  for (Int_t ix = 1; ix <= h->GetNbinsX(); ix++) {
2078  for (Int_t iy = 1; iy <= h->GetNbinsY(); iy++) {
2079  for (Int_t iz = 1; iz <= h->GetNbinsZ(); iz++) {
2080  Int_t bin = h->GetBin(ix,iy,iz);
2081  Double_t c = h->GetBinContent(bin);
2082  Double_t e = h->GetBinError (bin);
2083  if (c == 0 && ignoreZero) continue;
2084  min = TMath::Min(min, c-e);
2085  max = TMath::Max(max, c+e);
2086  }
2087  }
2088  }
2089  h->SetMinimum(min);
2090  h->SetMaximum(max);
2091 }
2092 //____________________________________________________________________
2093 void AliTrackletAODUtils::FixMinMax(THStack* stack, Bool_t ignoreZero)
2094 {
2095  Double_t min = +1e99;
2096  Double_t max = -1e99;
2097  TIter next(stack->GetHists());
2098  TH1* h = 0;
2099  while ((h = static_cast<TH1*>(next()))) {
2100  FixMinMax(h, ignoreZero);
2101  min = TMath::Min(min, h->GetMinimum());
2102  max = TMath::Max(max, h->GetMaximum());
2103  }
2104  stack->SetMinimum(min);
2105  stack->SetMaximum(max);
2106 }
2107 
2108 //____________________________________________________________________
2110 {
2111  static Int_t codes[] = {
2112  211, // #pi^{+}
2113  2212, // p
2114  321, // K^{+}
2115  323, // K^{*+}
2116  11, // e^{-}
2117  13, // #mu^{-}
2118  213, // #rho^{+}
2119  411, // D^{+}
2120  413, // D^{*+}
2121  431, // D_{s}^{+}
2122  433, // D_{s}^{*+}
2123  1114, // #Delta^{-}
2124  2214, // #Delta^{+}
2125  2224, // #Delta^{++}
2126  3112, // #Sigma^{-}
2127  3222, // #Sigma^{+}
2128  3114, // #Sigma^{*-}
2129  3224, // #Sigma^{*+}
2130  4214, // #Sigma^{*+}_{c}
2131  4224, // #Sigma^{*++}_{c}
2132  3312, // #Xi^{-}
2133  3314, // #Xi^{*-}
2134  4122, // #Lambda^{+}_{c}
2135  2112, // n
2136  2114, // #Delta^{0}
2137  22, // #gamma
2138  310, // K^{0}_{S}
2139  130, // K^{0}_{L}
2140  311, // K^{0}
2141  313, // K^{*}
2142  221, // #eta
2143  111, // #pi^{0}
2144  113, // #rho^{0}
2145  333, // #varphi
2146  331, // #eta'
2147  223, // #omega
2148  3122, // #Lambda
2149  3212, // #Sigma^{0}
2150  4114, // #Sigma^{*0}_{c}
2151  3214, // #Sigma^{*0}
2152  421, // D^{0}
2153  423, // D^{*0}
2154  3322, // #Xi_{0}
2155  3324, // #Xi^{*0}
2156  4132, // #Xi^{0}_{c}
2157  4314, // #Xi^{*0}_{c}
2158  1000000000 // Nuclei
2159  };
2160  static Int_t nCodes = sizeof(codes) / sizeof(Int_t);
2161  static Int_t* sorted = 0;
2162  size = nCodes;
2163  if (sorted) return sorted;
2164 
2165  sorted = new Int_t[nCodes];
2166  Int_t* idx = new Int_t[nCodes];
2167  TMath::Sort(nCodes, codes, idx, false);
2168  for (Int_t i = 0; i < nCodes; i++) sorted[i] = codes[idx[i]];
2169  delete [] idx;
2170  return sorted;
2171 }
2172 
2173 //____________________________________________________________________
2175  Color_t& c, Style_t& s)
2176 {
2177  // Leptons are black
2178  // Non-strange, non-charm meson are red
2179  // Non-strange, non-charm baryons are magenta
2180  // Strange-mesons are blue
2181  // Strange-baryons are cyan
2182  // Charmed mesons are green
2183  // Charmed baryons are yellow
2184  switch (pdg) {
2185  case -1: c=kPink +7; s=20; nme="?"; break;
2186  case 11: c=kBlack +0; s=20; nme="e^{-}"; break;
2187  case 13: c=kBlack +0; s=21; nme="#mu^{-}"; break;
2188  case 22: c=kBlack +0; s=22; nme="#gamma"; break;
2189  case 111: c=kRed +2; s=20; nme="#pi^{0}"; break;
2190  case 211: c=kRed +2; s=21; nme="#pi^{+}"; break;
2191  case 213: c=kRed +2; s=22; nme="#rho^{+}"; break;
2192  case 113: c=kRed +2; s=23; nme="#rho^{0}"; break;
2193  case 223: c=kRed +2; s=24; nme="#omega"; break;
2194  case 321: c=kBlue +2; s=20; nme="K^{+}"; break;
2195  case 323: c=kBlue +2; s=21; nme="K^{*+}"; break;
2196  case 310: c=kBlue +2; s=22; nme="K^{0}_{S}"; break;
2197  case 130: c=kBlue +2; s=23; nme="K^{0}_{L}"; break;
2198  case 311: c=kBlue +2; s=24; nme="K^{0}"; break;
2199  case 313: c=kBlue +2; s=25; nme="K^{*}"; break;
2200  case 221: c=kBlue +2; s=26; nme="#eta"; break;
2201  case 333: c=kBlue +2; s=27; nme="#varphi"; break;
2202  case 331: c=kBlue +2; s=28; nme="#eta'"; break;
2203  case 411: c=kGreen +2; s=20; nme="D^{+}"; break;
2204  case 413: c=kGreen +2; s=21; nme="D^{*+}"; break;
2205  case 421: c=kGreen +2; s=22; nme="D^{0}"; break;
2206  case 423: c=kGreen +2; s=23; nme="D^{*0}"; break;
2207  case 431: c=kGreen +2; s=24; nme="D_{s}^{+}"; break;
2208  case 433: c=kGreen +2; s=25; nme="D_{s}^{*+}"; break;
2209  case 2212: c=kMagenta+2; s=20; nme="p"; break;
2210  case 2112: c=kMagenta+2; s=21; nme="n"; break;
2211  case 2114: c=kMagenta+2; s=22; nme="#Delta^{0}"; break;
2212  case 1114: c=kMagenta+2; s=23; nme="#Delta^{-}"; break;
2213  case 2214: c=kMagenta+2; s=24; nme="#Delta^{+}"; break;
2214  case 2224: c=kMagenta+2; s=25; nme="#Delta^{++}"; break;
2215  case 3112: c=kCyan +2; s=20; nme="#Sigma^{-}"; break;
2216  case 3222: c=kCyan +2; s=21; nme="#Sigma^{+}"; break;
2217  case 3114: c=kCyan +2; s=22; nme="#Sigma^{*-}"; break;
2218  case 3224: c=kCyan +2; s=23; nme="#Sigma^{*+}"; break;
2219  case 3312: c=kCyan +2; s=24; nme="#Xi^{-}"; break;
2220  case 3314: c=kCyan +2; s=25; nme="#Xi^{*-}"; break;
2221  case 3122: c=kCyan +2; s=26; nme="#Lambda"; break;
2222  case 3212: c=kCyan +2; s=27; nme="#Sigma^{0}"; break;
2223  case 3214: c=kCyan +2; s=28; nme="#Sigma^{*0}"; break;
2224  case 3322: c=kCyan +2; s=29; nme="#Xi^{0}"; break;
2225  case 3324: c=kCyan +2; s=30; nme="#Xi^{*0}"; break;
2226  case 4214: c=kYellow +2; s=20; nme="#Sigma^{*+}_{c}"; break;
2227  case 4224: c=kYellow +2; s=21; nme="#Sigma^{*++}_{c}";break;
2228  case 4122: c=kYellow +2; s=22; nme="#Lambda^{+}_{c}"; break;
2229  case 4114: c=kYellow +2; s=23; nme="#Sigma^{*0}_{c}"; break;
2230  case 4132: c=kYellow +2; s=24; nme="#Xi^{0}_{c}"; break;
2231  case 4314: c=kYellow +2; s=25; nme="#Xi^{*0}_{c}"; break;
2232  case 1000000000:c=kPink +2; s=20; nme="Nuclei"; break;
2233  default: c=kGray; s=1; nme.Form("%d",pdg); break;
2234  };
2235 }
2236 
2237 //____________________________________________________________________
2239 {
2240  Int_t size;
2241  Int_t* array = PdgArray(size);
2242  Int_t apdg = TMath::Abs(pdg);
2243  Int_t idx = TMath::BinarySearch(size, array, apdg);
2244  if (idx == size - 1) return idx+1;
2245  if (array[idx] != apdg) return size+1;
2246  return idx+1;
2247 }
2248 
2249 //____________________________________________________________________
2251 {
2252  static TAxis* axis = 0;
2253  if (axis) return *axis;
2254 
2255  Int_t size;
2256  Int_t* array = PdgArray(size);
2257  axis = new TAxis(size+1, +.5, size+1.5);
2258  // TDatabasePDG* pdgDb = TDatabasePDG::Instance();
2259  for (Int_t i = 1; i <= size; i++) {
2260  Int_t pdg = array[i-1];
2261  axis->SetBinLabel(i, Form("%d", pdg));
2262  // TString name = "?";
2263  // TParticlePDG* type = pdgDb->GetParticle(pdg);
2264  // if (type) name = type->GetName();
2265  // axis->SetBinLabel(i, name.Data());
2266  }
2267  axis->SetBinLabel(size+1, "-1");
2268  FixAxis(*axis);
2269  return *axis;
2270 }
2271 
2272 #endif
2273 // Local Variables:
2274 // mode: C++
2275 // End:
2276 
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
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)
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 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 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 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 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)
Definition: External.C:220
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)