AliPhysics  781d0c7 (781d0c7)
 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 
256 {
257 public:
258  typedef TList Container;
266  virtual ~AliTrackletAODUtils() {}
267  //_________________________________________________________________________
281  static Bool_t CheckAxisNBins(const char* which,
282  const TAxis* a1,
283  const TAxis* a2);
293  static Bool_t CheckAxisLimits(const char* which,
294  const TAxis* a1,
295  const TAxis* a2);
306  static Bool_t CheckAxisBins(const char* which,
307  const TAxis* a1,
308  const TAxis* a2);
318  static Bool_t CheckAxisLabels(const char* which,
319  const TAxis* a1,
320  const TAxis* a2);
331  static Bool_t CheckAxis(const char* which,
332  const TAxis* a1,
333  const TAxis* a2,
334  Bool_t alsoLbls);
344  static Bool_t CheckConsistency(const TH1* h1, const TH1* h2);
345  /* @} */
346  //__________________________________________________________________
361  static TObject* GetO(Container* parent,
362  const char* name,
363  TClass* cls=0,
364  Bool_t verb=true);
375  static TObject* GetO(TDirectory* parent,
376  const char* name,
377  TClass* cls=0,
378  Bool_t verb=true);
388  static TH1* GetH1(Container* parent, const char* name, Bool_t verb=true);
398  static TH1* GetH1(TDirectory* parent, const char* name, Bool_t verb=true);
408  static TProfile* GetP1(Container* parent, const char* name, Bool_t verb=true);
418  static TH2* GetH2(Container* parent, const char* name, Bool_t verb=true);
428  static TH2* GetH2(TDirectory* parent, const char* name, Bool_t verb=true);
438  static TH3* GetH3(Container* parent, const char* name, Bool_t verb=true);
448  static TH3* GetH3(TDirectory* parent, const char* name, Bool_t verb=true);
458  static TProfile2D* GetP2(Container* parent, const char* name,
459  Bool_t verb=true);
469  static TProfile* GetP(Container* parent, const char* name, Bool_t verb=true)
470  {
471  return GetP1(parent, name, verb);
472  }
482  static THStack* GetHS(Container* parent, const char* name,Bool_t verb=true);
492  static THStack* GetHS(TDirectory* parent, const char* name,Bool_t verb=true);
502  static Container* GetC(Container* parent, const char* name, Bool_t verb=true);
512  static Container* GetC(TDirectory* parent, const char* name,Bool_t verb=true);
522  static TDirectory* GetT(TDirectory* parent,const char* name,Bool_t verb=true);
533  static Double_t GetD(Container* parent, const char* name,
534  Double_t def=-1, Bool_t verb=true);
545  static Int_t GetI(Container* parent, const char* name,
546  Int_t def=-1, Bool_t verb=true);
557  static Int_t GetB(Container* parent, const char* name,
558  Bool_t def=false, Bool_t verb=true);
569  static TH1* CopyH1(Container* parent,
570  const char* name,
571  const char* newName=0,
572  Bool_t verb=true);
583  static TH2* CopyH2(Container* parent,
584  const char* name,
585  const char* newName=0,
586  Bool_t verb=true);
597  static TH3* CopyH3(Container* parent,
598  const char* name,
599  const char* newName=0,
600  Bool_t verb=true);
607  static void CopyAttr(const TH1* src, TH1* dest);
620  static TH1* Make1D(Container* c,
621  const TString& name,
622  const TString& title,
623  Color_t color,
624  Style_t style,
625  const TAxis& xAxis);
638  static TProfile* Make1P(Container* c,
639  const TString& name,
640  const TString& title,
641  Color_t color,
642  Style_t style,
643  const TAxis& xAxis);
657  static TH2* Make2D(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);
678  static TH3* Make3D(Container* c,
679  const TString& name,
680  const TString& title,
681  Color_t color,
682  Style_t style,
683  const TAxis& xAxis,
684  const TAxis& yAxis,
685  const TAxis& zAxis);
699  static TProfile2D* Make2P(Container* c,
700  const TString& name,
701  const TString& title,
702  Color_t color,
703  Style_t style,
704  const TAxis& xAxis,
705  const TAxis& yAxis);
706 
716  static TH1* Scale(TH1* h, Double_t x, Double_t xe);
726  static TH2* Scale(TH2* h, Double_t x, Double_t xe);
736  static TH2* Scale(TH2* h, TH1* s);
743  static void FixAxis(TAxis& axis, const char* title=0);
751  static void ScaleAxis(TAxis& ret, Double_t fact=1);
759  static void SetAxis(TAxis& axis, Int_t n, Double_t* borders);
768  static void SetAxis(TAxis& axis, const TString& spec, const char* sep=":,");
777  static void SetAxis(TAxis& axis, Int_t n, Double_t l, Double_t h);
785  static void SetAxis(TAxis& axis, Int_t n, Double_t m);
793  static void PrintAxis(const TAxis& axis, Int_t nSig=2, const char* alt=0);
804  static TH2* ScaleToIPz(TH2* h, TH1* ipZ, Bool_t full=false);
805  static TH3* ScaleToIPz(TH3* h, TH1* ipZ, Bool_t full=false);
806  static TH3* ScaleDelta(TH3* h, TH2* scale);
831  static TH2* ProjectEtaDelta(TH3* h);
847  static TH1* ProjectDelta(TH2* h);
863  static TH1* ProjectDeltaFull(TH3* h);
885  static TH1* AverageOverIPz(TH2* h,
886  const char* name,
887  UShort_t mode,
888  TH1* ipz,
889  TH2* mask=0,
890  Bool_t verb=true);
899  static TObject* CloneAndAdd(Container* c, TObject* o);
912  static Double_t Integrate(TH1* h, Double_t min, Double_t max, Double_t& err);
924  static Double_t RatioE(Double_t n, Double_t en,
925  Double_t d, Double_t ed,
926  Double_t& er);
938  static Double_t RatioE2(Double_t n, Double_t e2n,
939  Double_t d, Double_t e2d,
940  Double_t& e2r);
950  static TH1* RatioH(const TH1* num, const TH1* denom, const char* name=0);
957  static void FixMinMax(TH1* h, Bool_t ignoreZero=true);
964  static void FixMinMax(THStack* h, Bool_t ignoreZero=true);
965 
966  /* @} */
975  static void PdgAttr(Int_t pdg, TString& nme, Color_t& c, Style_t& s);
980  static Int_t* PdgArray(Int_t& size);
986  static const TAxis& PdgAxis();
994  static Int_t PdgBin(Int_t pdg);
995  /* @} */
1004  static const char* CentName(Double_t c1, Double_t c2);
1014  static Color_t CentColor(const TAxis& axis,
1015  Double_t c1,
1016  Double_t c2);
1024  static Color_t CentColor(Int_t bin);
1032  static TFile* OpenFile(const char* filename);
1033 protected:
1034  ClassDef(AliTrackletAODUtils,1); // Utilities
1035 };
1036 
1037 //====================================================================
1038 // Utilities
1039 //____________________________________________________________________
1041  const TAxis* a1,
1042  const TAxis* a2)
1043 {
1044  if (a1->GetNbins() != a2->GetNbins()) {
1045  ::Warning("CheckAxisNBins", "Incompatible number %s bins: %d vs %d",
1046  which, a1->GetNbins(), a2->GetNbins());
1047  return false;
1048  }
1049  return true;
1050 }
1051 //____________________________________________________________________
1053  const TAxis* a1,
1054  const TAxis* a2)
1055 {
1056  if (TMath::AreEqualAbs(a1->GetXmin(), a2->GetXmin(), 1.e-3) &&
1057  TMath::AreEqualAbs(a1->GetXmax(), a2->GetXmax(), 1.e-3)) return true;
1058  if (!TMath::AreEqualRel(a1->GetXmin(), a2->GetXmin(),1.E-12) ||
1059  !TMath::AreEqualRel(a1->GetXmax(), a2->GetXmax(),1.E-12)) {
1060  Warning("CheckAxisLimits",
1061  "Limits of %s axis incompatible [%f,%f] vs [%f,%f]", which,
1062  a1->GetXmin(), a1->GetXmax(), a2->GetXmin(), a2->GetXmax());
1063  return false;
1064  }
1065  return true;
1066 }
1067 //____________________________________________________________________
1069  const TAxis* a1,
1070  const TAxis* a2)
1071 {
1072  const TArrayD * h1Array = a1->GetXbins();
1073  const TArrayD * h2Array = a2->GetXbins();
1074  Int_t fN = h1Array->fN;
1075  if ( fN == 0 ) return true;
1076  if (h2Array->fN != fN) {
1077  // Redundant?
1078  Warning("CheckAxisBins", "Not equal number of %s bin limits: %d vs %d",
1079  which, fN, h2Array->fN);
1080  return false;
1081  }
1082  else {
1083  for (int i = 0; i < fN; ++i) {
1084  if (TMath::AreEqualAbs(h1Array->GetAt(i),
1085  h2Array->GetAt(i), 1.e-3)) continue;
1086  if (!TMath::AreEqualRel(h1Array->GetAt(i),h2Array->GetAt(i),1E-10)) {
1087  Warning("CheckAxisBins",
1088  "%s limit # %3d incompatible: %f vs %f",
1089  which, i, h1Array->GetAt(i),h2Array->GetAt(i));
1090  return false;
1091  }
1092  }
1093  }
1094  return true;
1095 }
1096 //____________________________________________________________________
1098  const TAxis* a1,
1099  const TAxis* a2)
1100 {
1101  // check that axis have same labels
1102  THashList *l1 = (const_cast<TAxis*>(a1))->GetLabels();
1103  THashList *l2 = (const_cast<TAxis*>(a2))->GetLabels();
1104 
1105  if (!l1 && !l2) return true;
1106  if (!l1 || !l2) {
1107  Warning("CheckAxisLabels", "Difference in %s labels: %p vs %p",
1108  which, l1, l2);
1109  return false;
1110  }
1111  // check now labels sizes are the same
1112  if (l1->GetSize() != l2->GetSize()) {
1113  Warning("CheckAxisLabels", "Different number of %s labels: %d vs %d",
1114  which, l1->GetSize(), l2->GetSize());
1115  return false;
1116  }
1117  for (int i = 1; i <= a1->GetNbins(); ++i) {
1118  TString label1 = a1->GetBinLabel(i);
1119  TString label2 = a2->GetBinLabel(i);
1120  if (label1 != label2) {
1121  Warning("CheckAxisLabels", "%s label # %d not the same: '%s' vs '%s'",
1122  which, i, label1.Data(), label2.Data());
1123  return false;
1124  }
1125  }
1126 
1127  return true;
1128 }
1129 //____________________________________________________________________
1131  const TAxis* a1,
1132  const TAxis* a2,
1133  Bool_t alsoLbls)
1134 {
1135  if (a1->GetNbins() == 0 && a2->GetNbins() == 0) return true;
1136  if (!CheckAxisNBins (which, a1, a2)) return false;
1137  if (!CheckAxisLimits(which, a1, a2)) return false;
1138  if (!CheckAxisBins (which, a1, a2)) return false;
1139  if (alsoLbls && !CheckAxisLabels(which, a1, a2)) return false;
1140  return true;
1141 }
1142 //____________________________________________________________________
1144 {
1145  // Check histogram compatibility
1146  if (h1 == h2) return true;
1147 
1148  if (h1->GetDimension() != h2->GetDimension() ) {
1149  Warning("CheckConsistency",
1150  "%s and %s have different dimensions %d vs %d",
1151  h1->GetName(), h2->GetName(),
1152  h1->GetDimension(), h2->GetDimension());
1153  return false;
1154  }
1155  Int_t dim = h1->GetDimension();
1156 
1157  Bool_t alsoLbls = (h1->GetEntries() != 0 && h2->GetEntries() != 0);
1158  if (!CheckAxis("X", h1->GetXaxis(), h2->GetXaxis(), alsoLbls)) return false;
1159  if (dim > 1 &&
1160  !CheckAxis("Y", h1->GetYaxis(), h2->GetYaxis(), alsoLbls)) return false;
1161  if (dim > 2 &&
1162  !CheckAxis("Z", h1->GetZaxis(), h2->GetZaxis(), alsoLbls)) return false;
1163 
1164  return true;
1165 }
1166 
1167 //____________________________________________________________________
1169 {
1170  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
1171  Double_t dc = s->GetBinContent(i);
1172  Double_t de = s->GetBinError (i);
1173  Double_t dr = (dc > 1e-6 ? de/dc : 0);
1174  for (Int_t j = 1; j <= h->GetNbinsY(); j++) {
1175  Double_t nc = h->GetBinContent(i,j);
1176  Double_t ne = h->GetBinError (i,j);
1177  Double_t ns = (nc > 0 ? ne/nc : 0);
1178  Double_t sc = dc * nc;
1179  Double_t se = sc*TMath::Sqrt(ns*ns+dr*dr);
1180  // Printf("Setting bin %2d,%2d to %f +/- %f", sc, se);
1181  h->SetBinContent(i,j,sc);
1182  h->SetBinError (i,j,se);
1183  }
1184  }
1185  return h;
1186 }
1187 
1188 //____________________________________________________________________
1190 {
1191  Double_t rr = xe/x;
1192  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
1193  for (Int_t j = 1; j <= h->GetNbinsY(); j++) {
1194  Double_t c = h->GetBinContent(i,j);
1195  Double_t e = h->GetBinError (i,j);
1196  Double_t s = (c > 0 ? e/c : 0);
1197  Double_t sc = x * c;
1198  Double_t se = sc*TMath::Sqrt(s*s+rr*rr);
1199  // Printf("Setting bin %2d,%2d to %f +/- %f", sc, se);
1200  h->SetBinContent(i,j,sc);
1201  h->SetBinError (i,j,se);
1202  }
1203  }
1204  return h;
1205 }
1206 
1207 //____________________________________________________________________
1209 {
1210  Double_t rr = xe/x;
1211  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
1212  Double_t c = h->GetBinContent(i);
1213  Double_t e = h->GetBinError (i);
1214  Double_t s = (c > 0 ? e/c : 0);
1215  Double_t sc = x * c;
1216  Double_t se = sc*TMath::Sqrt(s*s+rr*rr);
1217  // Printf("Setting bin %2d,%2d to %f +/- %f", sc, se);
1218  h->SetBinContent(i,sc);
1219  h->SetBinError (i,se);
1220  }
1221  return h;
1222 }
1223 
1224 //____________________________________________________________________
1226  const char* name,
1227  TClass* cls,
1228  Bool_t verb)
1229 {
1230  if (!parent) {
1231  if (verb) ::Warning("GetO", "No parent container passed");
1232  return 0;
1233  }
1234  TObject* o = parent->FindObject(name);
1235  if (!o) {
1236  if (verb) ::Warning("GetO", "Object \"%s\" not found in \"%s\"",
1237  name, parent->GetName());
1238  // parent->ls();
1239  return 0;
1240  }
1241  if (!cls) return o;
1242  if (!o->IsA()->InheritsFrom(cls)) {
1243  if (verb) ::Warning("GetO", "\"%s\" is an object of class %s, not %s",
1244  name, o->ClassName(), cls->GetName());
1245  return 0;
1246  }
1247  return o;
1248 }
1249 //____________________________________________________________________
1251  const char* name,
1252  TClass* cls,
1253  Bool_t verb)
1254 {
1255  if (!parent) {
1256  if (!verb) ::Warning("GetO", "No parent directory passed");
1257  return 0;
1258  }
1259  TObject* o = parent->Get(name);
1260  if (!o) {
1261  if (verb) ::Warning("GetO", "Object \"%s\" not found in \"%s\"",
1262  name, parent->GetName());
1263  // parent->ls();
1264  return 0;
1265  }
1266  if (!cls) return o;
1267  if (!o->IsA()->InheritsFrom(cls)) {
1268  if (verb) ::Warning("GetO", "\"%s\" is an object of class %s, not %s",
1269  name, o->ClassName(), cls->GetName());
1270  return 0;
1271  }
1272  return o;
1273 }
1274 //____________________________________________________________________
1275 TH1* AliTrackletAODUtils::GetH1(Container* parent, const char* name, Bool_t v)
1276 {
1277  return static_cast<TH1*>(GetO(parent, name, TH1::Class(), v));
1278 }
1279 //____________________________________________________________________
1280 TH1* AliTrackletAODUtils::GetH1(TDirectory* parent, const char* name, Bool_t v)
1281 {
1282  return static_cast<TH1*>(GetO(parent, name, TH1::Class(), v));
1283 }
1284 //____________________________________________________________________
1285 TH2* AliTrackletAODUtils::GetH2(Container* parent, const char* name, Bool_t v)
1286 {
1287  return static_cast<TH2*>(GetO(parent, name, TH2::Class(), v));
1288 }
1289 //____________________________________________________________________
1290 TH2* AliTrackletAODUtils::GetH2(TDirectory* parent, const char* name, Bool_t v)
1291 {
1292  return static_cast<TH2*>(GetO(parent, name, TH2::Class(), v));
1293 }
1294 //____________________________________________________________________
1295 TH3* AliTrackletAODUtils::GetH3(Container* parent, const char* name, Bool_t v)
1296 {
1297  return static_cast<TH3*>(GetO(parent, name, TH3::Class(), v));
1298 }
1299 //____________________________________________________________________
1300 TH3* AliTrackletAODUtils::GetH3(TDirectory* parent, const char* name, Bool_t v)
1301 {
1302  return static_cast<TH3*>(GetO(parent, name, TH3::Class(), v));
1303 }
1304 //____________________________________________________________________
1305 TProfile* AliTrackletAODUtils::GetP1(Container* parent,const char* name,
1306  Bool_t v)
1307 {
1308  return static_cast<TProfile*>(GetO(parent, name, TProfile::Class(), v));
1309 }
1310 //____________________________________________________________________
1311 TProfile2D* AliTrackletAODUtils::GetP2(Container* parent, const char* name,
1312  Bool_t v)
1313 {
1314  return static_cast<TProfile2D*>(GetO(parent, name, TProfile2D::Class(), v));
1315 }
1316 //____________________________________________________________________
1317 THStack*
1318 AliTrackletAODUtils::GetHS(TDirectory* parent, const char* name, Bool_t v)
1319 {
1320  return static_cast<THStack*>(GetO(parent, name, THStack::Class(), v));
1321 }
1322 //____________________________________________________________________
1323 THStack*
1324 AliTrackletAODUtils::GetHS(Container* parent, const char* name, Bool_t v)
1325 {
1326  return static_cast<THStack*>(GetO(parent, name, THStack::Class(), v));
1327 }
1328 //____________________________________________________________________
1330 AliTrackletAODUtils::GetC(Container* parent, const char* name, Bool_t v)
1331 {
1332  return static_cast<Container*>(GetO(parent, name, Container::Class(), v));
1333 }
1334 //____________________________________________________________________
1336 AliTrackletAODUtils::GetC(TDirectory* parent, const char* name, Bool_t v)
1337 {
1338  return static_cast<Container*>(GetO(parent, name, Container::Class(), v));
1339 }
1340 //____________________________________________________________________
1341 TDirectory*
1342 AliTrackletAODUtils::GetT(TDirectory* parent, const char* name, Bool_t v)
1343 {
1344  TDirectory* d = parent->GetDirectory(name);
1345  if (!d) {
1346  if (v) ::Warning("GetO", "Directory \"%s\" not found in \"%s\"",
1347  name, parent->GetName());
1348  return 0;
1349  }
1350  return d;
1351 }
1352 //____________________________________________________________________
1354  const char* name,
1355  Double_t def,
1356  Bool_t v)
1357 {
1358  TParameter<double>* p =
1359  static_cast<TParameter<double>*>(GetO(parent, name,
1361  if (!p) return def;
1362  return p->GetVal();
1363 }
1364 //____________________________________________________________________
1366  const char* name,
1367  Int_t def,
1368  Bool_t v)
1369 {
1370  TParameter<int>* p =
1371  static_cast<TParameter<int>*>(GetO(parent, name,
1372  TParameter<int>::Class(),v));
1373  if (!p) return def;
1374  return p->GetVal();
1375 }
1376 //____________________________________________________________________
1378  const char* name,
1379  Bool_t def,
1380  Bool_t v)
1381 {
1382  TParameter<bool>* p =
1383  static_cast<TParameter<bool>*>(GetO(parent, name,
1384  TParameter<bool>::Class(), v));
1385  if (!p) return def;
1386  return p->GetVal();
1387 }
1388 //____________________________________________________________________
1390  const char* name,
1391  const char* newName,
1392  Bool_t v)
1393 {
1394  TH1* orig = GetH1(parent, name, v);
1395  if (!orig) return 0;
1396  TH1* ret = static_cast<TH1*>(orig->Clone(newName ? newName : name));
1397  ret->SetDirectory(0); // Release from file container
1398  return ret;
1399 }
1400 //____________________________________________________________________
1402  const char* name,
1403  const char* newName,
1404  Bool_t v)
1405 {
1406  TH2* orig = GetH2(parent, name, v);
1407  if (!orig) return 0;
1408  TH2* ret = static_cast<TH2*>(orig->Clone(newName ? newName : name));
1409  ret->SetDirectory(0); // Release from file container
1410  return ret;
1411 }
1412 //____________________________________________________________________
1414  const char* name,
1415  const char* newName,
1416  Bool_t v)
1417 {
1418  TH3* orig = GetH3(parent, name, v);
1419  if (!orig) return 0;
1420  TH3* ret = static_cast<TH3*>(orig->Clone(newName ? newName : name));
1421  ret->SetDirectory(0); // Release from file container
1422  return ret;
1423 }
1424 
1425 //____________________________________________________________________
1426 void AliTrackletAODUtils::CopyAttr(const TH1* src, TH1* dest)
1427 {
1428  if (!src || !dest) return;
1429  dest->SetMarkerStyle(src->GetMarkerStyle());
1430  dest->SetMarkerColor(src->GetMarkerColor());
1431  dest->SetMarkerSize (src->GetMarkerSize());
1432  dest->SetLineStyle (src->GetLineStyle());
1433  dest->SetLineColor (src->GetLineColor());
1434  dest->SetLineWidth (src->GetLineWidth());
1435  dest->SetFillStyle (src->GetFillStyle());
1436  dest->SetFillColor (src->GetFillColor());
1437 }
1438 
1439 
1440 //____________________________________________________________________
1442  const TString& name,
1443  const TString& title,
1444  Color_t color,
1445  Style_t style,
1446  const TAxis& xAxis)
1447 {
1448  TString n = name;
1449  TString t = title;
1450  TH1* ret = 0;
1451  Int_t nx = xAxis.GetNbins();
1452  if (t.IsNull()) t = xAxis.GetTitle();
1453  if (xAxis.GetXbins() && xAxis.GetXbins()->GetArray())
1454  ret = new TH1D(n,t,nx,xAxis.GetXbins()->GetArray());
1455  else
1456  ret = new TH1D(n,t,nx,xAxis.GetXmin(),xAxis.GetXmax());
1457  ret->Sumw2();
1458  ret->SetXTitle(xAxis.GetTitle());
1459  static_cast<const TAttAxis&>(xAxis).Copy(*(ret->GetXaxis()));
1460  ret->SetDirectory(0);
1461  ret->SetLineColor(color);
1462  ret->SetMarkerColor(color);
1463  ret->SetFillColor(kWhite);// color);
1464  ret->SetFillStyle(0);
1465  ret->SetMarkerStyle(style);
1466  if (const_cast<TAxis&>(xAxis).GetLabels()) {
1467  for (Int_t i = 1; i <= xAxis.GetNbins(); i++)
1468  ret->GetXaxis()->SetBinLabel(i, xAxis.GetBinLabel(i));
1469  }
1470  switch (style) {
1471  case 27:
1472  case 33:
1473  ret->SetMarkerSize(1.4);
1474  break;
1475  case 29:
1476  case 30:
1477  ret->SetMarkerSize(1.2);
1478  break;
1479  }
1480  if (c) c->Add(ret);
1481  return ret;
1482 }
1483 //____________________________________________________________________
1485  const TString& name,
1486  const TString& title,
1487  Color_t color,
1488  Style_t style,
1489  const TAxis& xAxis)
1490 {
1491  TString n = name;
1492  TString t = title;
1493  TProfile* ret = 0;
1494  Int_t nx = xAxis.GetNbins();
1495  if (t.IsNull()) t = xAxis.GetTitle();
1496  if (xAxis.GetXbins() && xAxis.GetXbins()->GetArray())
1497  ret = new TProfile(n,t,nx,xAxis.GetXbins()->GetArray());
1498  else
1499  ret = new TProfile(n,t,nx,xAxis.GetXmin(),xAxis.GetXmax());
1500  ret->Sumw2();
1501  ret->SetXTitle(xAxis.GetTitle());
1502  static_cast<const TAttAxis&>(xAxis).Copy(*(ret->GetXaxis()));
1503  ret->SetDirectory(0);
1504  ret->SetLineColor(color);
1505  ret->SetMarkerColor(color);
1506  ret->SetFillColor(kWhite);// color);
1507  ret->SetFillStyle(0);
1508  ret->SetMarkerStyle(style);
1509  if (const_cast<TAxis&>(xAxis).GetLabels()) {
1510  for (Int_t i = 1; i <= xAxis.GetNbins(); i++)
1511  ret->GetXaxis()->SetBinLabel(i, xAxis.GetBinLabel(i));
1512  }
1513  switch (style) {
1514  case 27:
1515  case 33:
1516  ret->SetMarkerSize(1.4);
1517  break;
1518  case 29:
1519  case 30:
1520  ret->SetMarkerSize(1.2);
1521  break;
1522  }
1523  if (c) c->Add(ret);
1524  return ret;
1525 }
1526 //____________________________________________________________________
1528  const TString& name,
1529  const TString& title,
1530  Color_t color,
1531  Style_t style,
1532  const TAxis& xAxis,
1533  const TAxis& yAxis)
1534 {
1535  TString n = name;
1536  TString t = title;
1537  TH2* ret = 0;
1538  Int_t nx = xAxis.GetNbins();
1539  Int_t ny = yAxis.GetNbins();
1540  const Double_t* xb = (xAxis.GetXbins() && xAxis.GetXbins()->GetArray() ?
1541  xAxis.GetXbins()->GetArray() : 0);
1542  const Double_t* yb = (yAxis.GetXbins() && yAxis.GetXbins()->GetArray() ?
1543  yAxis.GetXbins()->GetArray() : 0);
1544  if (t.IsNull())
1545  t.Form("%s vs %s", yAxis.GetTitle(), xAxis.GetTitle());
1546  if (xb) {
1547  if (yb) ret = new TH2D(n,t,nx,xb,ny,yb);
1548  else ret = new TH2D(n,t,
1549  nx,xb,
1550  ny,yAxis.GetXmin(),yAxis.GetXmax());
1551  }
1552  else {
1553  if (yb) ret = new TH2D(n,t,
1554  nx,xAxis.GetXmin(),xAxis.GetXmax(),
1555  ny,yb);
1556  else ret = new TH2D(n,t,
1557  nx,xAxis.GetXmin(),xAxis.GetXmax(),
1558  ny,yAxis.GetXmin(),yAxis.GetXmax());
1559  }
1560  ret->Sumw2();
1561  ret->SetXTitle(xAxis.GetTitle());
1562  ret->SetYTitle(yAxis.GetTitle());
1563  ret->SetLineColor(color);
1564  ret->SetMarkerColor(color);
1565  ret->SetFillColor(color);
1566  ret->SetMarkerStyle(style);
1567  static_cast<const TAttAxis&>(xAxis).Copy(*(ret->GetXaxis()));
1568  static_cast<const TAttAxis&>(yAxis).Copy(*(ret->GetYaxis()));
1569  ret->SetDirectory(0);
1570  if (const_cast<TAxis&>(xAxis).GetLabels()) {
1571  for (Int_t i = 1; i <= xAxis.GetNbins(); i++)
1572  ret->GetXaxis()->SetBinLabel(i, xAxis.GetBinLabel(i));
1573  }
1574  if (const_cast<TAxis&>(yAxis).GetLabels()) {
1575  for (Int_t i = 1; i <= yAxis.GetNbins(); i++)
1576  ret->GetYaxis()->SetBinLabel(i, yAxis.GetBinLabel(i));
1577  }
1578  if (c) c->Add(ret);
1579  return ret;
1580 }
1581 //____________________________________________________________________
1583  const TString& name,
1584  const TString& title,
1585  Color_t color,
1586  Style_t style,
1587  const TAxis& xAxis,
1588  const TAxis& yAxis,
1589  const TAxis& zAxis)
1590 {
1591  TString n = name;
1592  TString t = title;
1593  TH3* ret = 0;
1594  Int_t nx = xAxis.GetNbins();
1595  Int_t ny = yAxis.GetNbins();
1596  Int_t nz = zAxis.GetNbins();
1597  Double_t* xb = (xAxis.GetXbins() && xAxis.GetXbins()->GetArray() ?
1598  const_cast<Double_t*>(xAxis.GetXbins()->GetArray()) : 0);
1599  Double_t* yb = (yAxis.GetXbins() && yAxis.GetXbins()->GetArray() ?
1600  const_cast<Double_t*>(yAxis.GetXbins()->GetArray()) : 0);
1601  Double_t* zb = (zAxis.GetXbins() && zAxis.GetXbins()->GetArray() ?
1602  const_cast<Double_t*>(zAxis.GetXbins()->GetArray()) : 0);
1603  if (t.IsNull())
1604  t.Form("%s vs %s vs %s",
1605  zAxis.GetTitle(), yAxis.GetTitle(), xAxis.GetTitle());
1606  if (xb || yb || zb) {
1607  // One or more axis are defined as arrays. Make sure the rest are
1608  // also arrays.
1609  if (xb) {
1610  xb = new Double_t[nx+1];
1611  Double_t dx = (xAxis.GetXmax()-xAxis.GetXmin())/nx;
1612  xb[0] = xAxis.GetXmin();
1613  for (Int_t i = 1; i <= nx; i++) xb[i] = xb[i-1]+dx;
1614  }
1615  if (yb) {
1616  yb = new Double_t[ny+1];
1617  Double_t dy = (yAxis.GetXmax()-yAxis.GetXmin())/ny;
1618  yb[0] = yAxis.GetXmin();
1619  for (Int_t i = 1; i <= ny; i++) yb[i] = yb[i-1]+dy;
1620  }
1621  if (zb) {
1622  zb = new Double_t[nz+1];
1623  Double_t dz = (zAxis.GetXmax()-zAxis.GetXmin())/nz;
1624  zb[0] = zAxis.GetXmin();
1625  for (Int_t i = 1; i <= nz; i++) zb[i] = zb[i-1]+dz;
1626  }
1627  ret = new TH3D(n,t,nx,xb,ny,yb,nz,zb);
1628  }
1629  else {
1630  ret = new TH3D(n,t,
1631  nx, xAxis.GetXmin(), xAxis.GetXmax(),
1632  ny, yAxis.GetXmin(), yAxis.GetXmax(),
1633  nz, zAxis.GetXmin(), zAxis.GetXmax());
1634  }
1635  ret->Sumw2();
1636  ret->SetXTitle(xAxis.GetTitle());
1637  ret->SetYTitle(yAxis.GetTitle());
1638  ret->SetZTitle(zAxis.GetTitle());
1639  ret->SetLineColor(color);
1640  ret->SetMarkerColor(color);
1641  ret->SetFillColor(color);
1642  ret->SetMarkerStyle(style);
1643  static_cast<const TAttAxis&>(xAxis).Copy(*(ret->GetXaxis()));
1644  static_cast<const TAttAxis&>(yAxis).Copy(*(ret->GetYaxis()));
1645  static_cast<const TAttAxis&>(zAxis).Copy(*(ret->GetZaxis()));
1646  ret->SetDirectory(0);
1647  if (const_cast<TAxis&>(xAxis).GetLabels()) {
1648  for (Int_t i = 1; i <= xAxis.GetNbins(); i++)
1649  ret->GetXaxis()->SetBinLabel(i, xAxis.GetBinLabel(i));
1650  }
1651  if (const_cast<TAxis&>(yAxis).GetLabels()) {
1652  for (Int_t i = 1; i <= yAxis.GetNbins(); i++)
1653  ret->GetYaxis()->SetBinLabel(i, yAxis.GetBinLabel(i));
1654  }
1655  if (const_cast<TAxis&>(zAxis).GetLabels()) {
1656  for (Int_t i = 1; i <= zAxis.GetNbins(); i++)
1657  ret->GetZaxis()->SetBinLabel(i, zAxis.GetBinLabel(i));
1658  }
1659  if (c) c->Add(ret);
1660  return ret;
1661 }
1662 //____________________________________________________________________
1664  const TString& name,
1665  const TString& title,
1666  Color_t color,
1667  Style_t style,
1668  const TAxis& xAxis,
1669  const TAxis& yAxis)
1670 {
1671  TString n = name;
1672  TString t = title;
1673  TProfile2D* ret = 0;
1674  Int_t nx = xAxis.GetNbins();
1675  Int_t ny = yAxis.GetNbins();
1676  const Double_t* xb = (xAxis.GetXbins() && xAxis.GetXbins()->GetArray() ?
1677  xAxis.GetXbins()->GetArray() : 0);
1678  const Double_t* yb = (yAxis.GetXbins() && yAxis.GetXbins()->GetArray() ?
1679  yAxis.GetXbins()->GetArray() : 0);
1680  if (t.IsNull())
1681  t.Form("%s vs %s", yAxis.GetTitle(), xAxis.GetTitle());
1682  if (xb) {
1683  if (yb) ret = new TProfile2D(n,t,nx,xb,ny,yb);
1684  else ret = new TProfile2D(n,t,
1685  nx,xb,
1686  ny,yAxis.GetXmin(),yAxis.GetXmax());
1687  }
1688  else {
1689  if (yb) ret = new TProfile2D(n,t,
1690  nx,xAxis.GetXmin(),xAxis.GetXmax(),
1691  ny,yb);
1692  else ret = new TProfile2D(n,t,
1693  nx,xAxis.GetXmin(),xAxis.GetXmax(),
1694  ny,yAxis.GetXmin(),yAxis.GetXmax());
1695  }
1696  ret->Sumw2();
1697  ret->SetXTitle(xAxis.GetTitle());
1698  ret->SetYTitle(yAxis.GetTitle());
1699  ret->SetLineColor(color);
1700  ret->SetMarkerColor(color);
1701  ret->SetFillColor(color);
1702  ret->SetMarkerStyle(style);
1703  static_cast<const TAttAxis&>(xAxis).Copy(*(ret->GetXaxis()));
1704  static_cast<const TAttAxis&>(yAxis).Copy(*(ret->GetYaxis()));
1705  ret->SetDirectory(0);
1706  if (const_cast<TAxis&>(xAxis).GetLabels()) {
1707  for (Int_t i = 1; i <= xAxis.GetNbins(); i++)
1708  ret->GetXaxis()->SetBinLabel(i, xAxis.GetBinLabel(i));
1709  }
1710  if (const_cast<TAxis&>(yAxis).GetLabels()) {
1711  for (Int_t i = 1; i <= yAxis.GetNbins(); i++)
1712  ret->GetYaxis()->SetBinLabel(i, yAxis.GetBinLabel(i));
1713  }
1714  if (c) c->Add(ret);
1715  return ret;
1716 }
1717 //____________________________________________________________________
1718 void AliTrackletAODUtils::FixAxis(TAxis& axis, const char* title)
1719 {
1720  if (title && title[0] != '\0') axis.SetTitle(title);
1721  axis. SetNdivisions(210);
1722  axis. SetLabelFont(42);
1723  axis. SetLabelSize(0.03);
1724  axis. SetLabelOffset(0.005);
1725  axis. SetLabelColor(kBlack);
1726  axis. SetTitleOffset(1);
1727  axis. SetTitleFont(42);
1728  axis. SetTitleSize(0.04);
1729  axis. SetTitleColor(kBlack);
1730  axis. SetTickLength(0.03);
1731  axis. SetAxisColor(kBlack);
1732 }
1733 //____________________________________________________________________
1735  Double_t fact)
1736 {
1737  if (ret.GetXbins()->GetArray()) {
1738  TArrayD bins(*ret.GetXbins());
1739  for (Int_t i = 0; i < bins.GetSize(); i++) bins[i] *= fact;
1740  ret.Set(ret.GetNbins(), bins.GetArray());
1741  }
1742  else {
1743  ret.Set(ret.GetNbins(), fact*ret.GetXmin(), fact*ret.GetXmax());
1744  }
1745  FixAxis(ret);
1746 }
1747 //____________________________________________________________________
1749 {
1750  axis.Set(n, borders);
1751  FixAxis(axis);
1752 }
1753 //____________________________________________________________________
1755  const TString& spec,
1756  const char* sep)
1757 {
1758  TString s(spec);
1759  Bool_t isRange = false, isUnit = false;
1760  if (s[0] == 'r' || s[0] == 'R') {
1761  isRange = true;
1762  s.Remove(0,1);
1763  }
1764  if (s[0] == 'u' || s[0] == 'U') {
1765  isUnit = true;
1766  s.Remove(0,1);
1767  }
1768  TObjArray* tokens = s.Tokenize(sep);
1769  TArrayD bins(tokens->GetEntries());
1770  TObjString* token = 0;
1771  TIter next(tokens);
1772  Int_t i = 0;
1773  while ((token = static_cast<TObjString*>(next()))) {
1774  Double_t v = token->String().Atof();
1775  bins[i] = v;
1776  i++;
1777  }
1778  delete tokens;
1779  if (isUnit) {
1780  if (bins.GetSize() > 1)
1781  SetAxis(axis, Int_t(bins[1]-bins[0]), bins[0], bins[1]);
1782  else
1783  SetAxis(axis, 2*Int_t(bins[0]), bins[0]);
1784  }
1785  else if (isRange) {
1786  Int_t nBins = Int_t(bins[0]);
1787  if (bins.GetSize() > 2)
1788  SetAxis(axis, nBins, bins[1], bins[2]);
1789  else
1790  SetAxis(axis, nBins, bins[1]);
1791  }
1792  else
1793  SetAxis(axis, bins.GetSize()-1,bins.GetArray());
1794 }
1795 //____________________________________________________________________
1797  Int_t n,
1798  Double_t l,
1799  Double_t h)
1800 {
1801  axis.Set(n, l, h);
1802  FixAxis(axis);
1803 }
1804 //____________________________________________________________________
1806 {
1807  SetAxis(axis, n, -TMath::Abs(m), +TMath::Abs(m));
1808 }
1809 //____________________________________________________________________
1811  Int_t nSig,
1812  const char* alt)
1813 {
1814  printf(" %17s axis: ", alt ? alt : axis.GetTitle());
1815  if (axis.GetXbins() && axis.GetXbins()->GetArray()) {
1816  printf("%.*f", nSig, axis.GetBinLowEdge(1));
1817  for (Int_t i = 1; i <= axis.GetNbins(); i++)
1818  printf(":%.*f", nSig, axis.GetBinUpEdge(i));
1819  }
1820  else
1821  printf("%d bins between %.*f and %.*f",
1822  axis.GetNbins(), nSig, axis.GetXmin(),nSig,axis.GetXmax());
1823  printf("\n");
1824 }
1825 //____________________________________________________________________
1827 {
1828  if (!h) {
1829  ::Warning("ScaleToIPz","Nothing to scale");
1830  return 0;
1831  }
1832  if (!ipZ) {
1833  ::Warning("ScaleToIPz","Nothing to scale by");
1834  return 0;
1835  }
1836  TH2* ret = static_cast<TH2*>(h->Clone());
1837  ret->SetDirectory(0);
1838  if (!ipZ) return ret;
1839  for (Int_t iy = 1; iy <= ret->GetNbinsY(); iy++) {
1840  Double_t z = ret->GetYaxis()->GetBinCenter(iy);
1841  Int_t bin = ipZ->GetXaxis()->FindBin(z);
1842  Double_t nEv = ipZ->GetBinContent(bin);
1843  Double_t eEv = ipZ->GetBinError (bin);
1844  Double_t esc = (nEv > 0 ? 1./nEv : 0);
1845  Double_t rE2 = esc*esc*eEv*eEv;
1846  for (Int_t ix = 1; ix <= ret->GetNbinsX(); ix++) {
1847  Double_t c = ret->GetBinContent(ix,iy);
1848  Double_t e = ret->GetBinError (ix,iy);
1849  Double_t r = (c > 0 ? e/c : 0);
1850  // Scale by number of events, and error propagate
1851  Double_t sc = c * esc;
1852  Double_t se = 0;
1853  if (full) se = sc * TMath::Sqrt(r*r+rE2);
1854  else se = e * esc;
1855  Double_t scl = 1 / ret->GetXaxis()->GetBinWidth(ix);
1856  ret->SetBinContent(ix, iy, scl*sc);
1857  ret->SetBinError (ix, iy, scl*se);
1858  }
1859  }
1860  return ret;
1861 }
1862 //____________________________________________________________________
1864 {
1865  if (!h) {
1866  ::Warning("ScaleToIPz","Nothing to scale");
1867  return 0;
1868  }
1869  if (!ipZ) {
1870  ::Warning("ScaleToIPz","Nothing to scale by");
1871  return 0;
1872  }
1873  TH3* ret = static_cast<TH3*>(h->Clone());
1874  ret->SetDirectory(0);
1875  if (!ipZ) return ret;
1876  for (Int_t iz = 1; iz <= ret->GetNbinsZ(); iz++) {
1877  Double_t z = ret->GetZaxis()->GetBinCenter(iz);
1878  Int_t bin = ipZ->GetXaxis()->FindBin(z);
1879  Double_t nEv = ipZ->GetBinContent(bin);
1880  Double_t eEv = ipZ->GetBinError (bin);
1881  Double_t esc = (nEv > 0 ? 1./nEv : 0);
1882  Double_t rE2 = esc*esc*eEv*eEv;
1883  for (Int_t iy = 1; iy <= ret->GetNbinsY(); iy++) {
1884  for (Int_t ix = 1; ix <= ret->GetNbinsX(); ix++) {
1885  Double_t c = ret->GetBinContent(ix,iy,iz);
1886  Double_t e = ret->GetBinError (ix,iy,iz);
1887  Double_t r = (c > 0 ? e/c : 0);
1888  // Scale by number of events, and error propagate
1889  Double_t sc = c * esc;
1890  Double_t se = 0;
1891  if (full) se = sc * TMath::Sqrt(r*r+rE2);
1892  else se = e * esc;
1893  Double_t scl = 1 / ret->GetXaxis()->GetBinWidth(ix);
1894  ret->SetBinContent(ix, iy, iz, scl*sc);
1895  ret->SetBinError (ix, iy, iz, scl*se);
1896  }
1897  }
1898  }
1899  return ret;
1900 }
1901 //____________________________________________________________________
1903 {
1904  for (Int_t i = 1; i <= h->GetNbinsX(); i++) { // eta
1905  for (Int_t j = 1; j <= h->GetNbinsZ(); j++) { // IPz
1906  Double_t scale = etaIPzScale->GetBinContent(i,j);
1907  Double_t scaleE = etaIPzScale->GetBinError (i,j);
1908  Double_t q = (scale > 0 ? scaleE / scale : 0);
1909  for (Int_t k = 1; k <= h->GetNbinsY(); k++) { // delta
1910  Double_t c = h->GetBinContent(i,k,j);
1911  Double_t e = h->GetBinError (i,k,j);
1912  Double_t r = (c > 0 ? e / c : 0);
1913  Double_t w = c * scale;
1914  Double_t v = w * TMath::Sqrt(r*r + q*q);
1915  h->SetBinContent(i,k,j,w);
1916  h->SetBinError (i,k,j,v);
1917 #if 0
1918  Printf("%2d,%3d,%2d=(%9g+/-%9g)*(%9g+/-%9g)=(%9g+/-%9g)",
1919  i,k,j,c,e,scale,scaleE,w,v);
1920 #endif
1921  }
1922  }
1923  }
1924  return h;
1925 }
1926 //____________________________________________________________________
1927 // @cond
1929 {
1930  TH2* etaDelta = static_cast<TH2*>(h->Project3D("yx e"));
1931  etaDelta->SetName("etaDelta");
1932  etaDelta->SetTitle(h->GetTitle());
1933  etaDelta->SetDirectory(0);
1934  etaDelta->SetZTitle("d^{2}#it{N}/(d#etad#Delta)");
1935 #if 1
1936  // Reset content of projected histogram and calculate averages
1937  etaDelta->Reset();
1938  for (Int_t i = 1; i <= h->GetNbinsX(); i++) { // Loop over eta
1939  for (Int_t j = 1; j <= h->GetNbinsY(); j++) { // Loop over Delta
1940  Double_t sum = 0;
1941  Double_t sumw = 0;
1942  Int_t cnt = 0;
1943  for (Int_t k = 1; k <= h->GetNbinsZ(); k++) { // Loop over IPz
1944  Double_t c = h->GetBinContent(i,j,k);
1945  Double_t e = h->GetBinError (i,j,k);
1946  if (c < 1e-6 || e/c > 1) continue;
1947 #if 0
1948  Double_t w = 1/e/e;
1949  sum += w * c;
1950  sumw += w;
1951 #else
1952  sum += c;
1953  sumw += e*e;
1954 #endif
1955  cnt += 1;
1956  }
1957  if (sumw < 1e-6) continue;
1958 #if 0
1959  etaDelta->SetBinContent(i,j,sum/sumw);
1960  etaDelta->SetBinError (i,j,TMath::Sqrt(1/sumw));
1961 #else
1962  etaDelta->SetBinContent(i,j,sum/cnt);
1963  etaDelta->SetBinError (i,j,TMath::Sqrt(sumw)/cnt);
1964 #endif
1965  }
1966  }
1967 #endif
1968  return etaDelta;
1969 }
1970 // @endcond
1971 //____________________________________________________________________
1973 {
1974  if (!h) {
1975  Warning("ProjectDelta", "No histogram passed");
1976  return 0;
1977  }
1978  TH1* delta = h->ProjectionY("delta");
1979  delta->SetDirectory(0);
1980  delta->SetTitle(h->GetTitle());
1981  delta->SetYTitle("d#it{N}/d#Delta");
1982  delta->Scale(1./h->GetNbinsX());
1983  return delta;
1984 }
1985 //____________________________________________________________________
1987 {
1988  if (!h) {
1989  Warning("ProjectDelta", "No histogram passed");
1990  return 0;
1991  }
1992  TH2* tmp = ProjectEtaDelta(h);
1993  TH1* delta = ProjectDelta(tmp);
1994  delta->SetDirectory(0);
1995  tmp->SetDirectory(0);
1996  delete tmp;
1997  return delta;
1998 }
1999 
2000 //____________________________________________________________________
2002  const char* name,
2003  UShort_t mode,
2004  TH1* ipz,
2005  TH2* other,
2006  Bool_t verb)
2007 {
2008  if (!h) return 0;
2009 
2010  Int_t nIPz = h->GetNbinsY();
2011  Int_t nEta = h->GetNbinsX();
2012  TH1* p = h->ProjectionX(name,1,nIPz,"e");
2013  TH2* mask = (other ? other : h);
2014  p->SetDirectory(0);
2015  p->SetFillColor(0);
2016  p->SetFillStyle(0);
2017  p->SetYTitle(Form("#LT%s#GT", h->GetYaxis()->GetTitle()));
2018  p->Reset();
2019  for (Int_t etaBin = 1; etaBin <= nEta; etaBin++) {
2020  TArrayD hv(nIPz);
2021  TArrayD he(nIPz);
2022  TArrayD hr(nIPz);
2023  TArrayI hb(nIPz);
2024  Int_t j = 0;
2025  for (Int_t ipzBin = 1; ipzBin <= nIPz; ipzBin++) {
2026  Double_t bc = mask->GetBinContent(etaBin, ipzBin);
2027  if (bc < 1e-9) continue; // Low value
2028  Double_t be = mask->GetBinError (etaBin, ipzBin);
2029  if (TMath::IsNaN(be)) continue; // Bad error value
2030  hv[j] = bc;
2031  he[j] = be;
2032  hr[j] = be/bc;
2033  hb[j] = ipzBin;
2034  j++;
2035  }
2036  // Now we have all interesting values. Sort the relative error
2037  // array to get the most significant first. Note, we sort on the
2038  // mask errors which may not be the same as the histogram errors.
2039  TArrayI idx(nIPz);
2040  TMath::Sort(j, hr.fArray, idx.fArray, false);
2041  Double_t nsum = 0; // Running weighted sum
2042  Double_t nsume = 0; // Running weighted sum error
2043  Double_t dsum = 0;
2044  Double_t dsume = 0;
2045  Int_t n = 0;
2046  Double_t rat = 1e99;
2047 
2048  Int_t k = 0;
2049  for (k = 0; k < j; k++) {
2050  Int_t l = idx[k]; // Sorted index - ipzBin
2051  Int_t ipzBin = hb[l];
2052  Double_t hvv = hv[l];
2053  Double_t hee = he[l];
2054  Double_t x = TMath::Sqrt(nsume+hee*hee)/(nsum+hvv);
2055  if (x > rat) {
2056  continue; // Ignore - does not help
2057  }
2058 
2059  Double_t by = mask->GetYaxis()->GetBinCenter(ipzBin);
2060  Int_t ib = ipz ? ipz->FindBin(by) : 0;
2061  rat = x;
2062  nsum += h->GetBinContent(etaBin, ipzBin);
2063  nsume += TMath::Power(h->GetBinError(etaBin, ipzBin), 2);
2064  // If we do not have the vertex distribution, then just count
2065  // number of observations.
2066  dsum += !ipz ? 1 : ipz->GetBinContent(ib);
2067  dsume += !ipz ? 0 : TMath::Power(ipz->GetBinError(ib),2);
2068  n++;
2069  }
2070  if (k == 0 || n == 0) {
2071  if (verb)
2072  ::Warning("Average", "Eta bin # %3d has no data",etaBin);
2073  continue; // This eta bin empty!
2074  }
2075  Double_t norm = (mode > 0 ? n : dsum);
2076  Double_t rne = nsume/nsum/nsum;
2077  Double_t rde = dsume/dsum/dsum;
2078  Double_t av = nsum/norm;
2079  Double_t ave = 0;
2080  if (mode==2) ave = TMath::Sqrt(nsume)/n;
2081  else ave = av*TMath::Sqrt(rne+rde);
2082  p->SetBinContent(etaBin, av);
2083  p->SetBinError (etaBin, ave);
2084  }
2085  if (mode == 0) p->Scale(1, "width");
2086  return p;
2087 }
2088 
2089 //____________________________________________________________________
2091 {
2092  if (!o) {
2093  ::Warning("CloneAndAdd", "Nothing to clone");
2094  return 0;
2095  }
2096  TObject* copy = o->Clone();
2097  if (copy->IsA()->InheritsFrom(TH1::Class()))
2098  // Release from underlying directory
2099  static_cast<TH1*>(copy)->SetDirectory(0);
2100  if (c)
2101  // Add to output container
2102  c->Add(copy);
2103  return copy;
2104 }
2105 //____________________________________________________________________
2107  Double_t min,
2108  Double_t max,
2109  Double_t& err)
2110 {
2111  const Double_t epsilon = 1e-6;
2112  Int_t bMin = h->FindBin(min+epsilon);
2113  Int_t bMax = h->FindBin(max-epsilon);
2114  if (bMin < 1) bMin = 0; // Include underflow bin
2115  Double_t val = h->IntegralAndError(bMin, bMax, err);
2116  // For a-symmetric histograms, return
2117  if (TMath::Abs(h->GetXaxis()->GetXmin()+h->GetXaxis()->GetXmax())>=epsilon)
2118  return val;
2119 
2120  // Otherwise, also integrate negative side
2121  Double_t err2;
2122  bMin = h->FindBin(-min+epsilon);
2123  bMax = h->FindBin(-max-epsilon);
2124  val += h->IntegralAndError(bMin, bMax, err2);
2125  err = TMath::Sqrt(err*err+err2*err2);
2126  return val;
2127 }
2128 //____________________________________________________________________
2130  Double_t d, Double_t ed,
2131  Double_t& er)
2132 {
2133  Double_t r = 0;
2134  er = 0;
2135  if (TMath::Abs(n) < 1e-16 || TMath::Abs(d) < 1e-16) return 0;
2136  r = n/d;
2137  er = TMath::Sqrt(en*en/n/n + ed*ed/d/d);
2138  return r;
2139 }
2140 //____________________________________________________________________
2142  Double_t d, Double_t e2d,
2143  Double_t& e2r)
2144 {
2145  Double_t r = 0;
2146  e2r = 0;
2147  if (TMath::Abs(n) < 1e-16 || TMath::Abs(d) < 1e-16) return 0;
2148  r = n/d;
2149  e2r = (e2n/n/n + e2d/d/d);
2150  return r;
2151 }
2152 //____________________________________________________________________
2154  const TH1* denom,
2155  const char* name)
2156 {
2157  if (!num || !denom) {
2158  ::Warning("RatioH", "Numerator (%p) or denominator (%p) missing",
2159  num, denom);
2160  return 0;
2161  }
2162  if (!CheckConsistency(num, denom)) return 0;
2163 
2164  TH1* ratio = static_cast<TH1*>(num->Clone(name ? name : num->GetName()));
2165  ratio->SetDirectory(0);
2166  ratio->SetMinimum(-1111);
2167  ratio->SetMaximum(-1111);
2168  ratio->GetListOfFunctions()->Clear();
2169  ratio->Divide(denom);
2170  ratio->SetYTitle(Form("%s / %s",
2171  num ->GetYaxis()->GetTitle(),
2172  denom->GetYaxis()->GetTitle()));
2173  ratio->SetTitle(Form("%s / %s", num ->GetTitle(), denom->GetTitle()));
2174  FixMinMax(ratio, true);
2175 
2176  return ratio;
2177 }
2178 //____________________________________________________________________
2180 {
2181  Double_t min = +1e99;
2182  Double_t max = -1e99;
2183  for (Int_t ix = 1; ix <= h->GetNbinsX(); ix++) {
2184  for (Int_t iy = 1; iy <= h->GetNbinsY(); iy++) {
2185  for (Int_t iz = 1; iz <= h->GetNbinsZ(); iz++) {
2186  Int_t bin = h->GetBin(ix,iy,iz);
2187  Double_t c = h->GetBinContent(bin);
2188  Double_t e = h->GetBinError (bin);
2189  if (c == 0 && ignoreZero) continue;
2190  min = TMath::Min(min, c-e);
2191  max = TMath::Max(max, c+e);
2192  }
2193  }
2194  }
2195  h->SetMinimum(min);
2196  h->SetMaximum(max);
2197 }
2198 //____________________________________________________________________
2199 void AliTrackletAODUtils::FixMinMax(THStack* stack, Bool_t ignoreZero)
2200 {
2201  Double_t min = +1e99;
2202  Double_t max = -1e99;
2203  TIter next(stack->GetHists());
2204  TH1* h = 0;
2205  while ((h = static_cast<TH1*>(next()))) {
2206  FixMinMax(h, ignoreZero);
2207  min = TMath::Min(min, h->GetMinimum());
2208  max = TMath::Max(max, h->GetMaximum());
2209  }
2210  stack->SetMinimum(min);
2211  stack->SetMaximum(max);
2212 }
2213 
2214 //____________________________________________________________________
2216 {
2217  static Int_t codes[] = {
2218  211, // #pi^{+}
2219  2212, // p
2220  321, // K^{+}
2221  323, // K^{*+}
2222  11, // e^{-}
2223  13, // #mu^{-}
2224  213, // #rho^{+}
2225  411, // D^{+}
2226  413, // D^{*+}
2227  431, // D_{s}^{+}
2228  433, // D_{s}^{*+}
2229  1114, // #Delta^{-}
2230  2214, // #Delta^{+}
2231  2224, // #Delta^{++}
2232  3112, // #Sigma^{-}
2233  3222, // #Sigma^{+}
2234  3114, // #Sigma^{*-}
2235  3224, // #Sigma^{*+}
2236  4214, // #Sigma^{*+}_{c}
2237  4224, // #Sigma^{*++}_{c}
2238  3312, // #Xi^{-}
2239  3314, // #Xi^{*-}
2240  4122, // #Lambda^{+}_{c}
2241  2112, // n
2242  2114, // #Delta^{0}
2243  22, // #gamma
2244  310, // K^{0}_{S}
2245  130, // K^{0}_{L}
2246  311, // K^{0}
2247  313, // K^{*}
2248  221, // #eta
2249  111, // #pi^{0}
2250  113, // #rho^{0}
2251  333, // #varphi
2252  331, // #eta'
2253  223, // #omega
2254  3122, // #Lambda
2255  3212, // #Sigma^{0}
2256  4114, // #Sigma^{*0}_{c}
2257  3214, // #Sigma^{*0}
2258  421, // D^{0}
2259  423, // D^{*0}
2260  3322, // #Xi_{0}
2261  3324, // #Xi^{*0}
2262  4132, // #Xi^{0}_{c}
2263  4314, // #Xi^{*0}_{c}
2264  1000000000 // Nuclei
2265  };
2266  static Int_t nCodes = sizeof(codes) / sizeof(Int_t);
2267  static Int_t* sorted = 0;
2268  size = nCodes;
2269  if (sorted) return sorted;
2270 
2271  sorted = new Int_t[nCodes];
2272  Int_t* idx = new Int_t[nCodes];
2273  TMath::Sort(nCodes, codes, idx, false);
2274  for (Int_t i = 0; i < nCodes; i++) sorted[i] = codes[idx[i]];
2275  delete [] idx;
2276  return sorted;
2277 }
2278 
2279 //____________________________________________________________________
2281  Color_t& c, Style_t& s)
2282 {
2283  // Leptons are black
2284  // Non-strange, non-charm meson are red
2285  // Non-strange, non-charm baryons are magenta
2286  // Strange-mesons are blue
2287  // Strange-baryons are cyan
2288  // Charmed mesons are green
2289  // Charmed baryons are yellow
2290  switch (pdg) {
2291  case -1: c=kPink +7; s=20; nme="?"; break;
2292  case 11: c=kBlack +0; s=20; nme="e^{-}"; break;
2293  case 13: c=kBlack +0; s=21; nme="#mu^{-}"; break;
2294  case 22: c=kBlack +0; s=22; nme="#gamma"; break;
2295  case 111: c=kRed +2; s=20; nme="#pi^{0}"; break;
2296  case 211: c=kRed +2; s=21; nme="#pi^{+}"; break;
2297  case 213: c=kRed +2; s=22; nme="#rho^{+}"; break;
2298  case 113: c=kRed +2; s=23; nme="#rho^{0}"; break;
2299  case 223: c=kRed +2; s=24; nme="#omega"; break;
2300  case 321: c=kBlue +2; s=20; nme="K^{+}"; break;
2301  case 323: c=kBlue +2; s=21; nme="K^{*+}"; break;
2302  case 310: c=kBlue +2; s=22; nme="K^{0}_{S}"; break;
2303  case 130: c=kBlue +2; s=23; nme="K^{0}_{L}"; break;
2304  case 311: c=kBlue +2; s=24; nme="K^{0}"; break;
2305  case 313: c=kBlue +2; s=25; nme="K^{*}"; break;
2306  case 221: c=kBlue +2; s=26; nme="#eta"; break;
2307  case 333: c=kBlue +2; s=27; nme="#varphi"; break;
2308  case 331: c=kBlue +2; s=28; nme="#eta'"; break;
2309  case 411: c=kGreen +2; s=20; nme="D^{+}"; break;
2310  case 413: c=kGreen +2; s=21; nme="D^{*+}"; break;
2311  case 421: c=kGreen +2; s=22; nme="D^{0}"; break;
2312  case 423: c=kGreen +2; s=23; nme="D^{*0}"; break;
2313  case 431: c=kGreen +2; s=24; nme="D_{s}^{+}"; break;
2314  case 433: c=kGreen +2; s=25; nme="D_{s}^{*+}"; break;
2315  case 2212: c=kMagenta+2; s=20; nme="p"; break;
2316  case 2112: c=kMagenta+2; s=21; nme="n"; break;
2317  case 2114: c=kMagenta+2; s=22; nme="#Delta^{0}"; break;
2318  case 1114: c=kMagenta+2; s=23; nme="#Delta^{-}"; break;
2319  case 2214: c=kMagenta+2; s=24; nme="#Delta^{+}"; break;
2320  case 2224: c=kMagenta+2; s=25; nme="#Delta^{++}"; break;
2321  case 3112: c=kCyan +2; s=20; nme="#Sigma^{-}"; break;
2322  case 3222: c=kCyan +2; s=21; nme="#Sigma^{+}"; break;
2323  case 3114: c=kCyan +2; s=22; nme="#Sigma^{*-}"; break;
2324  case 3224: c=kCyan +2; s=23; nme="#Sigma^{*+}"; break;
2325  case 3312: c=kCyan +2; s=24; nme="#Xi^{-}"; break;
2326  case 3314: c=kCyan +2; s=25; nme="#Xi^{*-}"; break;
2327  case 3122: c=kCyan +2; s=26; nme="#Lambda"; break;
2328  case 3212: c=kCyan +2; s=27; nme="#Sigma^{0}"; break;
2329  case 3214: c=kCyan +2; s=28; nme="#Sigma^{*0}"; break;
2330  case 3322: c=kCyan +2; s=29; nme="#Xi^{0}"; break;
2331  case 3324: c=kCyan +2; s=30; nme="#Xi^{*0}"; break;
2332  case 4214: c=kYellow +2; s=20; nme="#Sigma^{*+}_{c}"; break;
2333  case 4224: c=kYellow +2; s=21; nme="#Sigma^{*++}_{c}";break;
2334  case 4122: c=kYellow +2; s=22; nme="#Lambda^{+}_{c}"; break;
2335  case 4114: c=kYellow +2; s=23; nme="#Sigma^{*0}_{c}"; break;
2336  case 4132: c=kYellow +2; s=24; nme="#Xi^{0}_{c}"; break;
2337  case 4314: c=kYellow +2; s=25; nme="#Xi^{*0}_{c}"; break;
2338  case 1000000000:c=kPink +2; s=20; nme="Nuclei"; break;
2339  default: c=kGray; s=1; nme.Form("%d",pdg); break;
2340  };
2341 }
2342 
2343 //____________________________________________________________________
2345 {
2346  Int_t size;
2347  Int_t* array = PdgArray(size);
2348  Int_t apdg = TMath::Abs(pdg);
2349  Int_t idx = TMath::BinarySearch(size, array, apdg);
2350  if (idx == size - 1) return idx+1;
2351  if (array[idx] != apdg) return size+1;
2352  return idx+1;
2353 }
2354 
2355 //____________________________________________________________________
2357 {
2358  static TAxis* axis = 0;
2359  if (axis) return *axis;
2360 
2361  Int_t size;
2362  Int_t* array = PdgArray(size);
2363  axis = new TAxis(size+1, +.5, size+1.5);
2364  // TDatabasePDG* pdgDb = TDatabasePDG::Instance();
2365  for (Int_t i = 1; i <= size; i++) {
2366  Int_t pdg = array[i-1];
2367  axis->SetBinLabel(i, Form("%d", pdg));
2368  // TString name = "?";
2369  // TParticlePDG* type = pdgDb->GetParticle(pdg);
2370  // if (type) name = type->GetName();
2371  // axis->SetBinLabel(i, name.Data());
2372  }
2373  axis->SetBinLabel(size+1, "-1");
2374  FixAxis(*axis);
2375  return *axis;
2376 }
2377 //____________________________________________________________________
2379 {
2380  TFile* file = TFile::Open(filename, "READ");
2381  if (!file) {
2382  ::Warning("OpenFile", "Failed to open \"%s\"", filename);
2383  return 0;
2384  }
2385  return file;
2386 }
2387 //____________________________________________________________________
2389 {
2390  static TString tmp;
2391  tmp.Form("cent%03dd%02d_%03dd%02d",
2392  Int_t(c1), Int_t(c1*100)%100,
2393  Int_t(c2), Int_t(c2*100)%100);
2394  return tmp.Data();
2395 }
2396 //____________________________________________________________________
2398 {
2399  const Color_t cc[] = { kMagenta+2, // 0
2400  kBlue+2, // 1
2401  kAzure-1, // 2 // 10,
2402  kCyan+2, // 3
2403  kGreen+1, // 4
2404  kSpring+5, // 5 //+10,
2405  kYellow+1, // 6
2406  kOrange+5, // 7 //+10,
2407  kRed+1, // 8
2408  kPink+5, // 9 //+10,
2409  kBlack }; // 10
2410  Color_t tc = cc[bin % 10];
2411  return tc;
2412 }
2413 //____________________________________________________________________
2415  Double_t c1,
2416  Double_t c2)
2417 {
2418  Double_t avgC = (c1+c2)/2;
2419  Int_t binC = const_cast<TAxis&>(axis).FindBin(avgC);
2420  return CentColor(binC);
2421 }
2422 
2423 
2424 
2425 #endif
2426 // Local Variables:
2427 // mode: C++
2428 // End:
2429 
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)