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