AliPhysics  05c4c93 (05c4c93)
AliTrackletdNdeta2.C
Go to the documentation of this file.
1 
11 #ifndef ALITRACKLETDNDETA_H
12 #define ALITRACKLETDNDETA_H
13 #include <AliTrackletAODUtils.C>
14 #ifndef __CINT__
15 #include <THStack.h>
16 #include <TFile.h>
17 #include <TLatex.h>
18 #include <TCanvas.h>
19 #include <TLegend.h>
20 #include <TError.h>
21 #include <TParameter.h>
22 #include <TMath.h>
23 #include <TLine.h>
24 #include <TLegendEntry.h>
25 #include <TStyle.h>
26 #include <TBrowser.h>
27 #include <TGraphAsymmErrors.h>
28 #include <TGraphErrors.h>
29 #include <TProfile.h>
30 #include <TFitResult.h>
31 #include <TF1.h>
32 #include <TSystem.h>
33 #include <TFractionFitter.h>
34 #else
35 class TPad;
36 class TLatex;
37 class TObject;
38 class TSeqCollection;
39 class TH1;
40 class TH2;
41 class TF1;
42 class TFitResultPtr;
43 class THStack;
44 class TCanvas;
45 class TVirtualPad;
46 class TFile;
47 class TAxis;
48 class TLegend;
49 class TDirectory;
50 #endif
51 
52 //====================================================================
59 {
65  enum {
67  kGeneral = 0x0001,
69  kParameters = 0x0002,
71  kWeights = 0x0004,
73  kdNdetas = 0x0008,
75  kSpecies = 0x0010,
77  kDeltas = 0x0020,
79  kDetails = 0x0040,
81  kPDF = 0x0100,
82  kPNG = 0x0100,
84  kPause = 0x0200,
86  kLandscape = 0x0400,
88  kAltMarker = 0x0800,
90  kDefaultViz = 0x03ff
91  };
96  enum {
98  kUnitScale = 0x0001,
100  kAverageScale = 0x0002,
102  kEtaScale = 0x0004,
104  kFullScale = 0x0008,
106  kDTFS = 0x0010,
108  kFake = 0x0020,
110  kCentAverage = 0x0040,
112  kClosure = 0x01000,
114  kDebug = 0x100,
116  kDefaultProc = 0x0002
117  };
118  enum {
119  kTopBackground = kAzure-8
120  };
121  //==================================================================
146  //==================================================================
148  TCanvas* fCanvas;
150  TPad* fTop;
152  TPad* fBody;
160  TLatex* fHeader;
162 
167 
168 
180  void Run(UInt_t proc = kDefaultProc,
181  UInt_t viz = kDefaultViz,
182  UShort_t maxBins = 9,
183  const char* dataName = "data.root",
184  const char* simName = "sim.root",
185  const char* output = 0,
186  Double_t fudge = 1);
187  //__________________________________________________________________
202  Bool_t Process(Container* realTop,
203  Container* simTop,
204  TDirectory* outTop,
205  Int_t maxBins);
220  Double_t c2,
221  Container* realTop,
222  Container* simTop,
223  Double_t realEff,
224  Double_t simEff,
225  TDirectory* outTop);
242  Double_t c2,
243  Container* realCont,
244  Container* simCont,
245  Double_t realEff,
246  Double_t simEff,
247  TDirectory* outTop,
248  TDirectory* outDir,
249  Int_t dimen);
250  /* @} */
251  //__________________________________________________________________
266  Bool_t CalculateSEF(Container* simCont,
267  const TAxis* centAxis,
268  TDirectory* out);
269  /* @} */
270  //__________________________________________________________________
285  Bool_t Deltas(Container* realCont,
286  Container* simCont,
287  TDirectory* outParent,
288  Int_t dim);
308  Bool_t Deltas0D(Container* realCont,
309  Container* simCont,
310  TDirectory* outParent);
329  Bool_t Deltas1D(Container* realCont,
330  Container* simCont,
331  TDirectory* outParent);
350  Bool_t Deltas2D(Container* realCont,
351  Container* simCont,
352  TDirectory* outParent);
373  Bool_t Deltas3D(Container* realCont,
374  Container* simCont,
375  TDirectory* outParent);
390  void WriteDeltas(TDirectory* outDir,
391  TH1* realDeltaM, TH1* realDeltaI,
392  TH1* simDeltaM, TH1* simDeltaI,
393  TH1* simDeltaC, TH1* simDeltaP,
394  TH1* simDeltaS, TH1* simDeltaD,
395  TH1* fit);
415  TH1* FractionFit(TDirectory* outDir,
416  TH1* realDeltaM,
417  TH1* simDeltaC,
418  TH1* simDeltaP,
419  TH1* simDeltaS);
420  /* @} */
421 
456  TH1* Results(Container* realCont,
457  Container* simCont,
458  Double_t realEff,
459  Double_t simEff,
460  TDirectory* outParent,
461  Int_t deltaDimen);
471  void FlattenCentrality(TH2* h, Container* c);
472 
473  /* @} */
474  //____________________________________________________________________
483  void ClearCanvas();
489  void CreateCanvas(const TString& outputName);
494  void CloseCanvas();
502  void PrintCanvas(const char* title,
503  const char* shortTitle="page",
504  Float_t size=.7);
515  TLegend* DrawInPad(TVirtualPad* c, Int_t pad, TObject* o, Option_t* opt);
526  void ModLegend(TVirtualPad* p, TLegend* l,
527  Double_t x1, Double_t y1,
528  Double_t x2, Double_t y2);
541  THStack* Make2Stack(const char* name,
542  const char* title,
543  Container* realList,
544  Container* simList,
545  Option_t* dataOpt="",
546  Option_t* simOpt="");
547  /* @} */
548  //____________________________________________________________________
565  Bool_t Visualize(Container* realSums,
566  Container* simSums,
567  Container* realRess,
568  Container* simRess,
569  TDirectory* outTop,
570  Int_t maxBins);
577  void VisualizeGeneral(Container* realList, Container* simList);
583  void VisualizeWeights(Container* simList);
590  void VisualizeFinal(TDirectory* outDir, Int_t i);
597  void VisualizeClosure(TDirectory* outDir, Int_t i);
609  Double_t c2,
610  Container* simList,
611  TDirectory* outTop);
619  Bool_t VisualizeSpecies(Container* simCont);
628  Bool_t VisualizeSpeciesDelta(Container* simCont, TDirectory* outDir);
636  Bool_t VisualizePrimary(Container* simCont);
645  Bool_t VisualizeDelta(TDirectory* outTop, Int_t dimen);
654  Bool_t VisualizeDetails(TDirectory* outTop, Int_t dimen);
663  Bool_t VisualizeResult(TDirectory* outTop, Int_t dimen);
664  /* @} */
665  //____________________________________________________________________
677  void VisualizeParam(const char* name, Double_t& y, const char* val);
685  void VisualizeParam(const char* name, Double_t& y, Double_t val);
693  void VisualizeParam(const char* name, Double_t& y, Int_t val);
701  void VisualizeParam(const char* name, Double_t& y, Bool_t val);
708  void VisualizeParams(Container* pars, const char* title);
715  void VisualizeParams(Container* realSums, Container* simSums);
716  /* @} */
717 
718  //____________________________________________________________________
736  TH1* SetAttr(TH1* h,
737  Color_t color,
738  Style_t marker=20,
739  Double_t size=1.,
740  Style_t fill=0,
741  Style_t line=1,
742  Width_t width=1);
743  static void StackMinMax(THStack* stack, Double_t& min, Double_t& max);
749  const char* ObsTitle() const { return "d#it{N}_{ch}/d#eta"; }
758  Double_t MeanY(TH1* h, Double_t& e);
767  Double_t MeanZ(TH2* h, Double_t& e);
768  /* @} */
769 };
770 
771 //====================================================================
772 namespace {
776  struct SuppressGuard2
777  {
779  Int_t save = 0;
785  SuppressGuard2(Int_t lvl=2000)
786  {
787  save = gErrorIgnoreLevel;
788  gErrorIgnoreLevel = lvl;
789  }
793  ~SuppressGuard2()
794  {
795  gErrorIgnoreLevel = save;
796  }
797  };
798 }
799 
800 //====================================================================
803  fProc(0),
804  fViz(0),
805  fDeltaCut(0),
806  fTailDelta(0),
807  fTailMax(0),
808  fMaxDelta(0),
809  fMinK(.7),
810  fMaxK(2),
811  fMinAlpha(0),
812  fMaxAlpha(2.5),
813  fFudge(0),
814  fSEF(1),
815  fCanvas(0),
816  fTop(0),
817  fBody(0),
818  fRealIsSim(false)
819 {}
820 
821 //====================================================================
823  UInt_t viz,
824  UShort_t maxBins,
825  const char* dataName,
826  const char* simName,
827  const char* outName,
828  Double_t fudge)
829 {
830  // Store options
831  fProc = proc;
832  fViz = viz;
833  fFudge = fudge;
834 
835  Printf("***************************************************\n"
836  " Data file: %s\n"
837  " Simulation file: %s\n"
838  " Output file: %s\n"
839  " Processing: 0x%x\n"
840  " Visualize: 0x%x\n"
841  " Max(Nbins): %d\n"
842  " Fudge: %f\n",
843  dataName, simName, outName, proc, viz, maxBins, fudge);
844 
845 
846  // Open the input files
847  TFile* dataFile = 0;
848  TFile* simFile = 0;
849  if (!(dataFile = OpenFile(dataName))) return;
850  if (!(simFile = OpenFile(simName))) return;
851 
852  // Get some top-level contianers
853  fRealIsSim = false;
854  const char* base = "MidRapidity";
855  Container* realSums = GetC(dataFile, Form("%sSums", base));
856  Container* realRess = GetC(dataFile, Form("%sResults", base));
857  Container* simSums = GetC(simFile, Form("%sMCSums", base));
858  Container* simRess = GetC(simFile, Form("%sMCResults", base));
859  if (!realSums || !realRess) {
860  realSums = GetC(dataFile, Form("%sMCSums", base));
861  realRess = GetC(dataFile, Form("%sMCResults", base));
862  if (realSums && realRess)
863  Warning("Run","\n"
864  "*********************************************\n"
865  "* Testing MC vs. MC: *\n"
866  "* 'Data' file: %23s *\n"
867  "* Simulation file: %23s *\n"
868  "*********************************************\n",
869  dataName, simName);
870  fRealIsSim = true;
871  }
872  if (!realSums || !realRess || !simSums || !simRess) return;
873 
874  // Get parameters from the real data file
875  Container* params = GetC(realSums, "parameters");
876  fDeltaCut = GetD(params, "DeltaCut");
877  fTailDelta = GetD(params, "TailDelta");
878  fTailMax = GetD(params, "TailMax");
879  fMaxDelta = GetD(params, "MaxDelta");
880 
881  // Create output file name
882  TString outBase(outName);
883  if (outBase.IsNull()) outBase.Form("MiddNdeta_0x%04x", fProc);
884  if (outBase.EndsWith(".root")) outBase.ReplaceAll(".root", "");
885 
886  // Make output directory
887  gSystem->mkdir(outBase);
888 
889  // Open the output file
890  TFile* out = TFile::Open(Form("%s/result.root", outBase.Data()), "RECREATE");
891 
892 
893  if (!Process(realRess, simRess, out, maxBins)) return;
894 
895  out->Write();
896 
897  if (fViz == 0) return; // Do not visualize anything
898 
899  Visualize(realSums, simSums, realRess, simRess, out, maxBins);
900 }
901 
902 //====================================================================
904  Container* simTop,
905  TDirectory* outDir,
906  Int_t maxBins)
907 {
908  DebugGuard g(fProc&kDebug,1,"Processing");
909  TH1* realCent = CopyH1(realTop, "cent", "realCent");
910  TH1* simCent = CopyH1(simTop, "cent", "simCent");
911  TH1* realIPz = GetH1(realTop, "ipz");
912  TH1* simIPz = GetH1(simTop, "ipz");
913 
914  // Check consistency of found histograms
915  if (!CheckConsistency(realCent, simCent)) {
916  Warning("Process", "Centrality bins are incompatible, giving up");
917  return false;
918  }
919  if (!CheckConsistency(realIPz, simIPz)) {
920  Warning("Process", "IPz bins are incompatible, giving up");
921  return false;
922  }
923 
924  // Check if we're doing a closure test
925  if (fProc & kClosure) realTop = simTop;
926 
927  PrintAxis(*realCent->GetXaxis(), 2, "Real centrality");
928  PrintAxis(*simCent ->GetXaxis(), 2, "Simulated centrality");
929 
930  THStack* mids = 0;
931  TH1* publ = 0;
932  if (realCent->GetNbinsX() > 1) {
933  Printf("Creating stack for mid rapidity results");
934  mids = new THStack("mids", "");
935  Int_t nbin = 9;
936  Double_t bins[] = { 0., 5., 10., 20., 30., 40., 50., 60., 70., 80., };
937  Double_t vals[] = { 1948, 1590, 1180, 786, 512, 318, 183, 96.3, 44.9, };
938  Double_t errs[] = { 38, 32, 31, 20, 15, 12, 8, 5.8, 3.4, };
939  publ = new TH1D("published", "PRL116,222302 (2016)",
940  nbin, bins);
941  publ->SetMarkerStyle(24);
942  publ->SetMarkerColor(kBlack);
943  publ->SetMarkerSize(1.3);
944  for (Int_t i = 0; i < nbin; i++) {
945  publ->SetBinContent(i+1, vals[i]);
946  publ->SetBinError (i+1, errs[i]);
947  }
948  mids->Add(publ);
949  }
950 
951  // Make histogram for mid-rapidiy results
952  for (Int_t d = 0; d < 4; d++) {
953  if ((fProc & (1 << d)) == 0) continue;
954  Printf("Making directory for final result");
955  TDirectory* dd = outDir->mkdir(Form("final%dd", d));
956  if (realCent->GetNbinsX() > 1) {
957  TH1* mid = Make1D(0,"mid",Form("%s|_{|#eta|<0.5}", ObsTitle()),
958  2+d, 20, *(realCent->GetXaxis()));
959  mid->SetDirectory(dd);
960  mid->SetXTitle("Centrality [%]");
961  mid->SetYTitle(mid->GetTitle());
962  mids->Add(mid);
963  }
964  THStack* full = new THStack("full","");
965  dd->cd();
966  full->Write();
967 
968  if (fRealIsSim) {
969  THStack* clos = new THStack("closure","");
970  clos->Write();
971  }
972  outDir->cd();
973  }
974 
975  // Write centralities to file
976  realCent->Write();
977  simCent ->Write();
978 
979  // Calculate strangeness enhancement factor
980  CalculateSEF(simTop,simCent->GetXaxis(),outDir);
981 
982  // Get trigger efficiencies
983  Double_t realEff = GetD(realTop, "triggerEfficiency", 1);
984  Double_t simEff = GetD(simTop, "triggerEfficiency", 1);
985  Printf("Trigger efficiencies: real=%6.4f sim=%6.4f", realEff, simEff);
986 
987  // Process MB bin
988  if (realCent->GetNbinsX() <= 1)
989  ProcessBin(0, 0, realTop, simTop, realEff, simEff, outDir);
990 
991  // Loop over defined centrality bins
992  for (Int_t i = 1; i <= realCent->GetNbinsX() && i <= maxBins ; i++) {
993  Double_t c1 = realCent->GetXaxis()->GetBinLowEdge(i);
994  Double_t c2 = realCent->GetXaxis()->GetBinUpEdge (i);
995 
996  ProcessBin(c1, c2, realTop, simTop, 1, 1, outDir);
997  }
998  // Close output file
999  outDir->cd();
1000  if (mids) {
1001  mids->Write();
1002 
1003  TH1* mid = static_cast<TH1*>(mids->GetHists()->At(1));
1004  for (Int_t i = 1; i <= mid->GetNbinsX(); i++) {
1005  Double_t c1 = mid->GetXaxis()->GetBinLowEdge(i);
1006  Double_t c2 = mid->GetXaxis()->GetBinUpEdge(i);
1007  Int_t j = publ->GetXaxis()->FindBin((c1+c2)/2);
1008  if (j < 1 || j > publ->GetNbinsX()) continue;
1009  Double_t vh = mid ->GetBinContent(i);
1010  Double_t eh = mid ->GetBinError (i);
1011  Double_t vp = publ->GetBinContent(j);
1012  Double_t ep = publ->GetBinError (j);
1013  Double_t er, r = RatioE(vh,eh, vp, ep, er);
1014  Printf("%5.1f - %5.1f%%: "
1015  "Here %6.1f +/- %4.1f "
1016  "Published %6.1f +/- %4.1f "
1017  "Ratio %5.3f +/- %5.3f",
1018  c1, c2, vh, eh, vp, ep, r, er);
1019  }
1020  }
1021 
1022  return true;
1023 }
1024 
1025 
1026 //____________________________________________________________________
1028  Double_t c2,
1029  Container* realTop,
1030  Container* simTop,
1031  Double_t realEff,
1032  Double_t simEff,
1033  TDirectory* outTop)
1034 {
1035  // Form the folder name
1036  TString centName(CentName(c1,c2));
1037  if (TMath::Abs(c1 - c2) < 1e-6) centName = "all";
1038  DebugGuard g(fProc&kDebug,1,"Processing bin %s", centName.Data());
1039 
1040  // Get centrality containers
1041  Container* realCont = GetC(realTop, centName);
1042  Container* simCont = GetC(simTop, centName);
1043  if (!realCont || !simCont) return false;
1044 
1045  // Test if we have any data at all
1046  TH1* realDist = GetH1(GetC(realCont, "measured"), "delta");
1047  TH1* simDist = GetH1(GetC(realCont, "measured"), "delta");
1048  if (realDist->GetEntries() <= 0 ||
1049  simDist ->GetEntries() <= 0) {
1050  Warning("ProcessBin", "No entries for bin %s", centName.Data());
1051  return false;
1052  }
1053  TDirectory* outDir = outTop->mkdir(centName);
1054 
1055  printf("%5.1f - %5.1f%%:", c1, c2);
1056  for (Int_t i = 0; i < 4; i++) {
1057  if ((fProc & (1 << i)) == 0) continue;
1058  if (!ProcessBin(c1, c2, realCont, simCont, realEff, simEff,
1059  outTop, outDir, i))
1060  return false;
1061  }
1062  printf("\n");
1063  return true;
1064 }
1065 //____________________________________________________________________
1067  Double_t c2,
1068  Container* realCont,
1069  Container* simCont,
1070  Double_t realEff,
1071  Double_t simEff,
1072  TDirectory* outTop,
1073  TDirectory* outDir,
1074  Int_t dimen)
1075 {
1076  if (!outDir) {
1077  Warning("ProcessBin", "No directory passed for %s and %s",
1078  realCont->GetName(), simCont->GetName());
1079  return false;
1080  }
1081  DebugGuard g(fProc&kDebug,1,"Processing bin %s", realCont->GetName());
1082  if (!Deltas(realCont, simCont, outDir, dimen)) return false;
1083  TH1* dndeta = Results(realCont, simCont, realEff, simEff, outDir, dimen);
1084  if (!dndeta) {
1085  Warning("ProcessBin", "Failed on Deltas for %f - %f", c1, c2);
1086  return false;
1087  }
1088 
1089  TDirectory* final = GetT(outTop,Form("final%dd", dimen));
1090  if (!final) {
1091  Warning("ProcessBin", "Failed on results for %f - %f", c1, c2);
1092  return false;
1093  }
1094 
1095  TH1* mid = GetH1(final, "mid");
1096  THStack* full = GetHS(final, "full");
1097  if (!final) {
1098  Warning("ProcessBin", "Missing full (%p)", full);
1099  return false;
1100  }
1101 
1102  Int_t b = 10;
1103  if (mid) {
1104  TF1* f = static_cast<TF1*>(dndeta->GetListOfFunctions()->At(0));
1105  if (!f) {
1106  Warning("ProcessBin", "No fit found on %s", dndeta->GetTitle());
1107  return false;
1108  }
1109 
1110  Double_t c = (c1+c2)/2;
1111  b = mid->GetXaxis()->FindBin(c);
1112  if (b < 1 || b > mid->GetNbinsX()) {
1113  Warning("ProcessBin", "Centrality %f - %f out of range", c1, c2);
1114  return false;
1115  }
1116  mid->SetBinContent(b, f->GetParameter(0));
1117  mid->SetBinError (b, f->GetParError (0));
1118  }
1119 
1120  Color_t tc = CentColor(b);
1121  TH1* copy = static_cast<TH1*>(dndeta->Clone(outDir->GetName()));
1122  copy->SetDirectory(final);
1123  copy->GetListOfFunctions()->Clear();
1124  copy->SetTitle(Form("%5.1f#minus%5.1f%%", c1, c2));
1125  SetAttr(copy, tc);
1126  full->Add(copy);
1127  final->cd();
1128  full->Write(full->GetName(), TObject::kOverwrite);
1129 
1130  THStack* clos = GetHS(final, "closure", false);
1131  TH1* clss = GetH1(GetT(outDir, Form("results%dd",dimen)),
1132  "closure", false);
1133  if (clos && clss) {
1134  copy = static_cast<TH1*>(clss->Clone(outDir->GetName()));
1135  copy->SetDirectory(0);
1136  copy->GetListOfFunctions()->Clear();
1137  copy->SetTitle(Form("%5.1f#minus%5.1f%%", c1, c2));
1138  SetAttr(copy, tc);
1139  clos->Add(copy);
1140  clos->Write(clos->GetName(), TObject::kOverwrite);
1141  }
1142 
1143 
1144  return true;
1145 }
1146 
1147 //====================================================================
1149  const TAxis* centAxis,
1150  TDirectory* out)
1151 {
1152  DebugGuard g(fProc&kDebug,1,"Calculating strangeness enhancement factor");
1153  Double_t c1 = centAxis->GetBinLowEdge(1);
1154  Double_t c2 = centAxis->GetBinUpEdge(1);
1155  TString centName(CentName(c1,c2));
1156  if (TMath::Abs(c1 - c2) < 1e-6) centName = "all";
1157 
1158  // Find the most central bin
1159  Container* centBin = GetC(simCont, centName);
1160  if (!centBin) return false;
1161 
1162  Container* generated = GetC(centBin, "generated");
1163  if (!generated) return false;
1164 
1165  Container* mix = GetC(generated, "mix");
1166  if (!mix) return false;
1167 
1168  THStack* toPion = GetHS(mix, "toPion");
1169  if (!toPion) return false;
1170 
1171  TIter next(toPion->GetHists());
1172  TH1* hist = 0;
1173  Double_t sum = 0;
1174  Double_t sumw = 0;
1175  while ((hist = static_cast<TH1*>(next()))) {
1176  Double_t r2760 = hist->GetBinContent(0);
1177  Double_t e2760 = 0.07*r2760; // Fixed 7% error
1178 
1179  TF1* f = new TF1("f", "pol0", -.5, +.5);
1180  hist->Fit(f, "QN", "", -.5, +.5);
1181  Double_t rHere = f->GetParameter(0);
1182  Double_t eHere = f->GetParError (0);
1183 
1184  Double_t eCh, rCh = RatioE(r2760, e2760, rHere, eHere, eCh);
1185 #if 0
1186  Printf("%20s: @ 2.76TeV=%6.4f+/-%6.4f here=%6.4f+/-%6.4f -> %6.4f+/-%6.4f",
1187  hist->GetTitle(), r2760, e2760, rHere, eHere, rCh, eCh);
1188 #endif
1189  delete f;
1190 
1191  sum += rHere*rCh;
1192  sumw += rHere;
1193  }
1194  Double_t avg = 1;
1195  if (sumw>0) avg = sum / sumw;
1196  Printf(" Weighted average of factor: %6.4f", avg);
1197  Printf(" Preset: %6.4f", fSEF);
1198  if (fSEF == 1) fSEF = avg;
1199  Printf( "Strangeness enhancement factor set to %6.4f", fSEF);
1200 
1201  return true;
1202 }
1203 
1204 //====================================================================
1206  Container* simCont,
1207  TDirectory* outParent,
1208  Int_t dim)
1209 {
1210  DebugGuard g(fProc&kDebug,1,"Doing Delta calculations");
1211  if (!outParent) {
1212  Warning("Deltas", "No directory passed!");
1213  return false;
1214  }
1215  switch (dim) {
1216  case 0: return Deltas0D(realCont, simCont, outParent);
1217  case 1: return Deltas1D(realCont, simCont, outParent);
1218  case 2: return Deltas2D(realCont, simCont, outParent);
1219  case 3: return Deltas3D(realCont, simCont, outParent);
1220  }
1221  return false;
1222 }
1223 //____________________________________________________________________
1225  Container* simCont,
1226  TDirectory* outParent)
1227 {
1228  DebugGuard g(fProc&kDebug,1,"Doing 0D Delta calculations");
1229  // Make an output directory
1230  TDirectory* outDir = outParent->mkdir("delta0d");
1231 
1232  // Get the real and simulated measured folders
1233  Container* realMeas = GetC(realCont, "measured");
1234  Container* simMeas = GetC(simCont, "measured");
1235  if (!realMeas || !simMeas) return false;
1236 
1237  // Create a flat 2D scaling histogram for later user
1238  TH2* h = CopyH2(realMeas, "etaIPz", "scale");
1239  TH1* hp = h->ProjectionX("scaleProj");
1240  hp->Reset();
1241  hp->SetMinimum(fMinK);
1242  hp->SetMaximum(fMaxK);
1243  hp->SetYTitle("k");
1244  h->SetZTitle("k");
1245  h->SetTitle("k=1");
1246  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
1247  hp->SetBinContent(i, 1);
1248  hp->SetBinError (i, 0);
1249  for (Int_t j = 1; j <= h->GetNbinsY(); j++) {
1250  if (h->GetBinContent(i,j) < 1e-6) {
1251  h->SetBinContent(i,j,0);
1252  h->SetBinError (i,j,0);
1253  }
1254  else {
1255  h->SetBinContent(i,j,1);
1256  h->SetBinError (i,j,0);
1257  }
1258  }
1259  }
1260  h->SetDirectory(outDir);
1261  h->SetMinimum(fMinK);
1262  h->SetMaximum(fMaxK);
1263  hp->SetDirectory(outDir);
1264 
1265  // Get the raw projected Delta distributions for each component
1266  TH1* realDeltaM = CopyH1(realMeas, "delta","realDeltaM");
1267  TH1* realDeltaI = CopyH1(GetC(realCont,"injected"), "delta","realDeltaI");
1268  TH1* simDeltaM = CopyH1(simMeas, "delta","simDeltaM");
1269  TH1* simDeltaI = CopyH1(GetC(simCont, "injected"), "delta","simDeltaI");
1270  TH1* simDeltaC = CopyH1(GetC(simCont, "combinatorics"),"delta","simDeltaC");
1271  TH1* simDeltaP = CopyH1(GetC(simCont, "primaries"), "delta","simDeltaP");
1272  TH1* simDeltaS = CopyH1(GetC(simCont, "secondaries"), "delta","simDeltaS");
1273  TH1* simDeltaD = CopyH1(GetC(simCont, "distinct"), "delta","simDeltaD");
1274 
1275  // Get integrated scaling factor for injections, and scale the
1276  // injection distributions by that
1277  Double_t realScaleI = GetD(GetC(realCont,"injected"), "scale");
1278  Double_t realScaleIE = GetD(GetC(realCont,"injected"), "scaleError");
1279  Double_t simScaleI = GetD(GetC(simCont, "injected"), "scale");
1280  Double_t simScaleIE = GetD(GetC(simCont, "injected"), "scaleError");
1281  Scale(realDeltaI, realScaleI, realScaleIE);
1282  Scale(simDeltaI, simScaleI, simScaleIE);
1283  realDeltaI->SetTitle(Form("k_{I}#times%s",realScaleI,realDeltaI->GetTitle()));
1284  simDeltaI ->SetTitle(Form("k_{I'}#times%s",simScaleI,simDeltaI ->GetTitle()));
1285 
1286  TH1* fit = FractionFit(outDir, realDeltaM, simDeltaC, simDeltaP, simDeltaS);
1287  WriteDeltas(outDir, realDeltaM, realDeltaI, simDeltaM, simDeltaI,
1288  simDeltaC, simDeltaP, simDeltaS, simDeltaD, fit);
1289 
1290  outParent->cd();
1291  return true;
1292 }
1293 
1294 //____________________________________________________________________
1296  Container* simCont,
1297  TDirectory* outParent)
1298 {
1299  DebugGuard g(fProc&kDebug,1,"Doing 1D Delta calculations");
1300  // Make an output directory
1301  if (!outParent) {
1302  Warning("Deltas1D", "No directory passed!");
1303  return false;
1304  }
1305  TDirectory* outDir = outParent->mkdir("delta1d");
1306 
1307  // Get the real and simulated measured folders
1308  Container* realMeas = GetC(realCont, "measured");
1309  Container* simMeas = GetC(simCont, "measured");
1310  if (!realMeas || !simMeas) return false;
1311 
1312  // Get the integrated tails of the real and simulated observed
1313  // distribtutions, and calcualte scaling factor.
1314  Double_t realTail = GetD(realMeas, "deltaTail");
1315  Double_t realTailE = GetD(realMeas, "deltaTailError");
1316  Double_t simTail = GetD(simMeas, "deltaTail");
1317  Double_t simTailE = GetD(simMeas, "deltaTailError");
1318  Double_t scaleE = 0;
1319  Double_t scale = RatioE(realTail, realTailE, simTail, simTailE, scaleE);
1320 
1321  // Create a flat 2D scaling histogram for later user
1322  TH2* h = CopyH2(realMeas, "etaIPz", "scale");
1323  TH1* hp = h->ProjectionX("scaleProj");
1324  hp->Reset();
1325  hp->SetMinimum(fMinK);
1326  hp->SetMaximum(fMaxK);
1327  hp->SetYTitle("k");
1328  h->SetZTitle("k");
1329  h->SetTitle(Form("k=%5.3f#pm%5.3f", scale, scaleE));
1330  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
1331  hp->SetBinContent(i, scale);
1332  hp->SetBinError (i, scaleE);
1333  for (Int_t j = 1; j <= h->GetNbinsY(); j++) {
1334  if (h->GetBinContent(i,j) < 1e-6) {
1335  h->SetBinContent(i,j,0);
1336  h->SetBinError (i,j,0);
1337  }
1338  else {
1339  h->SetBinContent(i,j,scale);
1340  h->SetBinError (i,j,scaleE);
1341  }
1342  }
1343  }
1344  h->SetDirectory(outDir);
1345  h->SetMinimum(fMinK);
1346  h->SetMaximum(fMaxK);
1347  hp->SetDirectory(outDir);
1348 
1349  // Get the raw projected Delta distributions for each component
1350  TH1* realDeltaM = CopyH1(realMeas, "delta","realDeltaM");
1351  TH1* realDeltaI = CopyH1(GetC(realCont,"injected"), "delta","realDeltaI");
1352  TH1* simDeltaM = CopyH1(simMeas, "delta","simDeltaM");
1353  TH1* simDeltaI = CopyH1(GetC(simCont, "injected"), "delta","simDeltaI");
1354  TH1* simDeltaC = CopyH1(GetC(simCont, "combinatorics"),"delta","simDeltaC");
1355  TH1* simDeltaP = CopyH1(GetC(simCont, "primaries"), "delta","simDeltaP");
1356  TH1* simDeltaS = CopyH1(GetC(simCont, "secondaries"), "delta","simDeltaS");
1357  TH1* simDeltaD = CopyH1(GetC(simCont, "distinct"), "delta","simDeltaD");
1358 
1359  // Get integrated scaling factor for injections, and scale the
1360  // injection distributions by that
1361  Double_t realScaleI = GetD(GetC(realCont,"injected"), "scale");
1362  Double_t realScaleIE = GetD(GetC(realCont,"injected"), "scaleError");
1363  Double_t simScaleI = GetD(GetC(simCont, "injected"), "scale");
1364  Double_t simScaleIE = GetD(GetC(simCont, "injected"), "scaleError");
1365  Scale(realDeltaI, realScaleI, realScaleIE);
1366  Scale(simDeltaI, simScaleI, simScaleIE);
1367  realDeltaI->SetTitle(Form("k_{I}#times%s", realDeltaI->GetTitle()));
1368  simDeltaI ->SetTitle(Form("k_{I'}#times%s",simDeltaI ->GetTitle()));
1369 
1370  TH1* toScale[] = { simDeltaM,simDeltaI,simDeltaC,
1371  simDeltaP,simDeltaS,simDeltaD, 0};
1372  TH1** pScale = toScale;
1373  while ((*pScale)) {
1374  Scale(*pScale, scale, scaleE);
1375  (*pScale)->SetTitle(Form("k_{M}#times%s", (*pScale)->GetTitle()));
1376  pScale++;
1377  }
1378 
1379  TH1* fit = FractionFit(outDir, realDeltaM, simDeltaC, simDeltaP, simDeltaS);
1380  WriteDeltas(outDir, realDeltaM, realDeltaI, simDeltaM, simDeltaI,
1381  simDeltaC, simDeltaP, simDeltaS, simDeltaD, fit);
1382 
1383  outParent->cd();
1384  return true;
1385 }
1386 
1387 //____________________________________________________________________
1389  Container* simCont,
1390  TDirectory* outParent)
1391 {
1392  DebugGuard g(fProc&kDebug,1,"Doing 2D Delta calculations");
1393  // Make an output directory
1394  TDirectory* outDir = outParent->mkdir("delta2d");
1395 
1396  // Get the real and simulated measured folders
1397  Container* realMeas = GetC(realCont, "measured");
1398  Container* simMeas = GetC(simCont, "measured");
1399  if (!realMeas || !simMeas) return false;
1400 
1401  // Get the eta-differential tails of the real and simulated observed
1402  // distribtutions, and calcualte scaling factor.
1403  TH1* realTail = GetH1(realMeas, "etaDeltaTail");
1404  TH1* simTail = GetH1(simMeas, "etaDeltaTail");
1405  TH1* scale = static_cast<TH1*>(realTail->Clone("scaleProj"));
1406  scale->Divide(simTail);
1407  scale->SetYTitle("k");
1408  Double_t sE, s = MeanY(scale, sE);
1409  TGraphErrors* avgScale = new TGraphErrors(2);
1410  avgScale->SetLineStyle(2);
1411  avgScale->SetLineColor(kBlack);
1412  avgScale->SetFillColor(kYellow);
1413  avgScale->SetFillStyle(3002);
1414  avgScale->SetPoint (0, scale->GetXaxis()->GetXmin(), s);
1415  avgScale->SetPoint (1, scale->GetXaxis()->GetXmax(), s);
1416  avgScale->SetPointError(0, 0, sE);
1417  avgScale->SetPointError(1, 0, sE);
1418  scale->GetListOfFunctions()->Add(avgScale, "le3");
1419 
1420  // Create a flat 2D scaling histogram for later user
1421  TH2* h = CopyH2(realMeas, "etaIPz", "scale");
1422  h->SetZTitle("k");
1423  h->SetTitle(Form("#LTk(#eta)#GT_{#eta}=%5.3f#pm%5.3f", s, sE));
1424  for (Int_t i = 1; i <= h->GetNbinsX(); i++) { // eta
1425  for (Int_t j = 1; j <= h->GetNbinsY(); j++) { // ipz
1426  if (h->GetBinContent(i,j) < 1e-6) {
1427  h->SetBinContent(i,j,0);
1428  h->SetBinError (i,j,0);
1429  }
1430  else {
1431  h->SetBinContent(i,j,scale->GetBinContent(i));
1432  h->SetBinError (i,j,scale->GetBinError (i));
1433  }
1434  }
1435  }
1436  h->SetMinimum(fMinK);
1437  h->SetMaximum(fMaxK);
1438  h->SetDirectory(outDir);
1439  scale->SetMinimum(fMinK);
1440  scale->SetMaximum(fMaxK);
1441  scale->SetDirectory(outDir);
1442 
1443  // Get the raw projected Delta distributions for each component
1444  TH2* r2DeltaM = CopyH2(realMeas, "etaDelta","r2DeltaM");
1445  TH2* r2DeltaI = CopyH2(GetC(realCont,"injected"), "etaDelta","r2DeltaI");
1446  TH2* s2DeltaM = CopyH2(simMeas, "etaDelta","s2DeltaM");
1447  TH2* s2DeltaI = CopyH2(GetC(simCont, "injected"), "etaDelta","s2DeltaI");
1448  TH2* s2DeltaC = CopyH2(GetC(simCont, "combinatorics"),"etaDelta","s2DeltaC");
1449  TH2* s2DeltaP = CopyH2(GetC(simCont, "primaries"), "etaDelta","s2DeltaP");
1450  TH2* s2DeltaS = CopyH2(GetC(simCont, "secondaries"), "etaDelta","s2DeltaS");
1451  TH2* s2DeltaD = CopyH2(GetC(simCont, "distinct"), "etaDelta","s2DeltaD");
1452 
1453  // Get eta-differential scaling factor for injections, and scale the
1454  // injection distributions by that
1455  TH1* rScaleI = GetH1(GetC(realCont,"injected"), "etaScale");
1456  TH1* sScaleI = GetH1(GetC(simCont, "injected"), "etaScale");
1457  Double_t rIE, rI = MeanY(rScaleI, rIE);
1458  Double_t sIE, sI = MeanY(sScaleI, sIE);
1459  Scale(r2DeltaI, rScaleI);
1460  Scale(s2DeltaI, sScaleI);
1461  r2DeltaI ->SetTitle(Form("#LTk_{I}#GT_{#eta}#times%s",
1462  r2DeltaI->GetTitle()));
1463  s2DeltaI ->SetTitle(Form("#LTk_{I'}#GT_{#eta}#times%s",
1464  s2DeltaI->GetTitle()));
1465 
1466  TH2* toScale[] = { s2DeltaM,s2DeltaI,s2DeltaC,s2DeltaP,s2DeltaS,s2DeltaD,0 };
1467  TH2** pScale = toScale;
1468  while ((*pScale)) {
1469  Scale(*pScale, scale);
1470  (*pScale)->SetTitle(Form("#LTk_{M}#GT_{#eta}#times%s",
1471  (*pScale)->GetTitle()));
1472  pScale++;
1473  }
1474 
1475  TH1* rDeltaM = ProjectDelta(r2DeltaM);
1476  TH1* rDeltaI = ProjectDelta(r2DeltaI);
1477  TH1* sDeltaM = ProjectDelta(s2DeltaM);
1478  TH1* sDeltaI = ProjectDelta(s2DeltaI);
1479  TH1* sDeltaC = ProjectDelta(s2DeltaC);
1480  TH1* sDeltaP = ProjectDelta(s2DeltaP);
1481  TH1* sDeltaS = ProjectDelta(s2DeltaS);
1482  TH1* sDeltaD = ProjectDelta(s2DeltaD);
1483  rDeltaM->SetTitle(r2DeltaM->GetTitle()); rDeltaM->SetName("realDeltaM");
1484  rDeltaI->SetTitle(r2DeltaI->GetTitle()); rDeltaI->SetName("realDeltaI");
1485  sDeltaM->SetTitle(s2DeltaM->GetTitle()); sDeltaM->SetName("simDeltaM");
1486  sDeltaI->SetTitle(s2DeltaI->GetTitle()); sDeltaI->SetName("simDeltaI");
1487  sDeltaC->SetTitle(s2DeltaC->GetTitle()); sDeltaC->SetName("simDeltaC");
1488  sDeltaP->SetTitle(s2DeltaP->GetTitle()); sDeltaP->SetName("simDeltaP");
1489  sDeltaS->SetTitle(s2DeltaS->GetTitle()); sDeltaS->SetName("simDeltaS");
1490  sDeltaD->SetTitle(s2DeltaD->GetTitle()); sDeltaD->SetName("simDeltaD");
1491 
1492  TH1* f2 = FractionFit(outDir, r2DeltaM, s2DeltaC, s2DeltaP, s2DeltaS);
1493  TH1* fit = ProjectDelta(static_cast<TH2*>(f2));
1494  WriteDeltas(outDir,rDeltaM,rDeltaI,
1495  sDeltaM,sDeltaI,sDeltaC,
1496  sDeltaP,sDeltaS,sDeltaD,
1497  fit);
1498 
1499  TDirectory* full = outDir->mkdir("full");
1500  r2DeltaM->SetDirectory(full); r2DeltaM->SetName("realDeltaM");
1501  r2DeltaI->SetDirectory(full); r2DeltaI->SetName("realDeltaI");
1502  s2DeltaM->SetDirectory(full); s2DeltaM->SetName("simDeltaM");
1503  s2DeltaI->SetDirectory(full); s2DeltaI->SetName("simDeltaI");
1504  s2DeltaC->SetDirectory(full); s2DeltaC->SetName("simDeltaC");
1505  s2DeltaP->SetDirectory(full); s2DeltaP->SetName("simDeltaP");
1506  s2DeltaS->SetDirectory(full); s2DeltaS->SetName("simDeltaS");
1507  s2DeltaD->SetDirectory(full); s2DeltaD->SetName("simDeltaD");
1508 
1509  outParent->cd();
1510  return true;
1511 }
1512 
1513 //____________________________________________________________________
1515  Container* simCont,
1516  TDirectory* outParent)
1517 {
1518  DebugGuard g(fProc&kDebug,1,"Doing 3D Delta calculations");
1519  // Make an output directory
1520  TDirectory* outDir = outParent->mkdir("delta3d");
1521 
1522  // Get the real and simulated measured folders
1523  Container* realMeas = GetC(realCont, "measured");
1524  Container* simMeas = GetC(simCont, "measured");
1525  if (!realMeas || !simMeas) return false;
1526 
1527  // Get the eta-differential tails of the real and simulated observed
1528  // distribtutions, and calcualte scaling factor.
1529  TH2* realTail = GetH2(realMeas, "etaIPzDeltaTail");
1530  TH2* simTail = GetH2(simMeas, "etaIPzDeltaTail");
1531  TH2* scale = static_cast<TH2*>(realTail->Clone("scale"));
1532  scale->SetDirectory(0);
1533  scale->Divide(simTail);
1534  Double_t sE, s = MeanZ(scale, sE);
1535  scale->SetZTitle("k");
1536  scale->SetTitle(Form("#LTk(#eta)#GT_{#eta,IP_{#it{z}}}=%5.3f#pm%5.3f",
1537  s, sE));
1538  scale->SetDirectory(outDir);
1539  scale->SetMinimum(fMinK);
1540  scale->SetMaximum(fMaxK);
1541  TH1* etaScale = AverageOverIPz(scale, "scaleProj", 1, 0, 0);
1542  etaScale->SetYTitle("#LTk#GT_{IP_{z}}");
1543  etaScale->SetDirectory(outDir);
1544  etaScale->SetMinimum(fMinK);
1545  etaScale->SetMaximum(fMaxK);
1546  TGraphErrors* avgScale = new TGraphErrors(2);
1547  avgScale->SetLineStyle(2);
1548  avgScale->SetLineColor(kBlack);
1549  avgScale->SetFillColor(kYellow);
1550  avgScale->SetFillStyle(3002);
1551  avgScale->SetPoint(0, etaScale->GetXaxis()->GetXmin(), s);
1552  avgScale->SetPoint(1, etaScale->GetXaxis()->GetXmax(), s);
1553  avgScale->SetPointError(0, 0, sE);
1554  avgScale->SetPointError(1, 0, sE);
1555  etaScale->GetListOfFunctions()->Add(avgScale, "le3");
1556 
1557  // Get the raw projected Delta distributions for each component
1558  const char* nm = "etaDeltaIPz";
1559  TH3* r3DeltaM = CopyH3(realMeas, nm,"r3DeltaM");
1560  TH3* r3DeltaI = CopyH3(GetC(realCont,"injected"), nm,"r3DeltaI");
1561  TH3* s3DeltaM = CopyH3(simMeas, nm,"s3DeltaM");
1562  TH3* s3DeltaI = CopyH3(GetC(simCont, "injected"), nm,"s3DeltaI");
1563  TH3* s3DeltaC = CopyH3(GetC(simCont, "combinatorics"),nm,"s3DeltaC");
1564  TH3* s3DeltaP = CopyH3(GetC(simCont, "primaries"), nm,"s3DeltaP");
1565  TH3* s3DeltaS = CopyH3(GetC(simCont, "secondaries"), nm,"s3DeltaS");
1566  TH3* s3DeltaD = CopyH3(GetC(simCont, "distinct"), nm,"s3DeltaD");
1567 
1568  // Get eta-differential scaling factor for injections, and scale the
1569  // injection distributions by that
1570  TH2* rScaleI = GetH2(GetC(realCont,"injected"), "etaIPzScale");
1571  TH2* sScaleI = GetH2(GetC(simCont, "injected"), "etaIPzScale");
1572  Double_t rIE, rI = MeanZ(rScaleI, rIE);
1573  Double_t sIE, sI = MeanZ(sScaleI, sIE);
1574  ScaleDelta(r3DeltaI, rScaleI);
1575  ScaleDelta(s3DeltaI, sScaleI);
1576  r3DeltaI ->SetTitle(Form("#LTk_{I}#GT_{#eta,IP_{#it{z}}}#times%s",
1577  r3DeltaI->GetTitle()));
1578  s3DeltaI ->SetTitle(Form("#LTk_{I'}#GT_{#eta,IP_{#it{z}}}#times%s",
1579  s3DeltaI->GetTitle()));
1580 #if 0
1581  TH2* scale2 = static_cast<TH2*>(scale->Clone("scaleMain"));
1582  scale2->SetDirectory(0);
1583  scale2->Reset();
1584  Int_t sigBin = r3DeltaM->GetYaxis()->FindBin(1.5);
1585  for (Int_t i = 1; i <= r3DeltaM->GetNbinsX(); i++) {
1586  for (Int_t j = 1; j <= r3DeltaM->GetNbinsZ(); j++) {
1587  // Integrate over Delta
1588  Double_t rintg = 0, reintg = 0;
1589  Double_t sintg = 0, seintg = 0;
1590  rintg = r3DeltaM->IntegralAndError(i,i,1,sigBin,j,j,reintg);
1591  sintg = s3DeltaM->IntegralAndError(i,i,1,sigBin,j,j,seintg);
1592  Double_t re, r = RatioE(rintg, reintg, sintg, seintg, re);
1593 
1594  scale2->SetBinContent(i, j, r);
1595  scale2->SetBinError (i, j, re);
1596  }
1597  }
1598  Double_t rS2, rS = MeanZ(scale2, rS2);
1599  Printf("Scalar of Inj %6.3f +/- %6.3f", rS, rS2);
1600 #endif
1601 
1602  TH3* toScale[] = { s3DeltaM,s3DeltaI,s3DeltaC,s3DeltaP,s3DeltaS,s3DeltaD,0};
1603  TH3** pScale = toScale;
1604  while ((*pScale)) {
1605  ScaleDelta(*pScale, scale);
1606  (*pScale)->SetTitle(Form("#LTk_{M}#GT_{#eta,IP_{#it{z}}}#times%s",
1607  (*pScale)->GetTitle()));
1608  pScale++;
1609  }
1610 
1611  TH1* rDeltaM = ProjectDeltaFull(r3DeltaM);
1612  TH1* rDeltaI = ProjectDeltaFull(r3DeltaI);
1613  TH1* sDeltaM = ProjectDeltaFull(s3DeltaM);
1614  TH1* sDeltaI = ProjectDeltaFull(s3DeltaI);
1615  TH1* sDeltaC = ProjectDeltaFull(s3DeltaC);
1616  TH1* sDeltaP = ProjectDeltaFull(s3DeltaP);
1617  TH1* sDeltaS = ProjectDeltaFull(s3DeltaS);
1618  TH1* sDeltaD = ProjectDeltaFull(s3DeltaD);
1619  rDeltaM->SetTitle(r3DeltaM->GetTitle()); rDeltaM->SetName("realDeltaM");
1620  rDeltaI->SetTitle(r3DeltaI->GetTitle()); rDeltaI->SetName("realDeltaI");
1621  sDeltaM->SetTitle(s3DeltaM->GetTitle()); sDeltaM->SetName("simDeltaM");
1622  sDeltaI->SetTitle(s3DeltaI->GetTitle()); sDeltaI->SetName("simDeltaI");
1623  sDeltaC->SetTitle(s3DeltaC->GetTitle()); sDeltaC->SetName("simDeltaC");
1624  sDeltaP->SetTitle(s3DeltaP->GetTitle()); sDeltaP->SetName("simDeltaP");
1625  sDeltaS->SetTitle(s3DeltaS->GetTitle()); sDeltaS->SetName("simDeltaS");
1626  sDeltaD->SetTitle(s3DeltaD->GetTitle()); sDeltaD->SetName("simDeltaD");
1627 
1628  TH1* f3 = FractionFit(outDir, r3DeltaM, s3DeltaC, s3DeltaP, s3DeltaS);
1629  TH3* ff3 = static_cast<TH3*>(f3);
1630 
1631  TDirectory* full = outDir->mkdir("full");
1632  r3DeltaM->SetDirectory(full); r3DeltaM->SetName("realDeltaM");
1633  r3DeltaI->SetDirectory(full); r3DeltaI->SetName("realDeltaI");
1634  s3DeltaM->SetDirectory(full); s3DeltaM->SetName("simDeltaM");
1635  s3DeltaI->SetDirectory(full); s3DeltaI->SetName("simDeltaI");
1636  s3DeltaC->SetDirectory(full); s3DeltaC->SetName("simDeltaC");
1637  s3DeltaP->SetDirectory(full); s3DeltaP->SetName("simDeltaP");
1638  s3DeltaS->SetDirectory(full); s3DeltaS->SetName("simDeltaS");
1639  s3DeltaD->SetDirectory(full); s3DeltaD->SetName("simDeltaD");
1640  TH1* fit = 0;
1641  if (ff3){
1642  ff3->SetDirectory(full);
1643  ff3->SetName("simDeltaF");
1644  TH2* fetaDelta = static_cast<TH2*>(ff3->Project3D("yx e"));
1645  fetaDelta->SetName("simDeltaFF");
1646  fetaDelta->SetDirectory(full);
1647  fetaDelta->Scale(1./ff3->GetNbinsZ());
1648  outDir->cd();
1649  fit = fetaDelta->ProjectionY("simDeltaF");
1650  fit->SetTitle("#Delta_{F}");
1651  fit->SetDirectory(outDir);
1652  fit->Scale(1. / fetaDelta->GetNbinsX());
1653  // delete fetaDelta;
1654  }
1655  WriteDeltas(outDir,rDeltaM,rDeltaI,
1656  sDeltaM,sDeltaI,sDeltaC,
1657  sDeltaP,sDeltaS,sDeltaD,
1658  fit);
1659 
1660  outParent->cd();
1661  return true;
1662 }
1663 
1664 //____________________________________________________________________
1666  TH1* rDeltaM,
1667  TH1* sDeltaC,
1668  TH1* sDeltaP,
1669  TH1* sDeltaS)
1670 {
1671  DebugGuard g(fProc&kDebug,1,"Fraction fitting in %s", outDir->GetName());
1672  // We don't do this, as it doesn't seem to do much.
1673  return 0;
1674  if (!rDeltaM || !sDeltaC || !sDeltaP || !sDeltaS) {
1675  Warning("FractionFit", "Missing M=%p, C'=%s, P'=%s, or S'=%p",
1676  rDeltaM, sDeltaC, sDeltaP, sDeltaS);
1677  return 0;
1678  }
1679  TDirectory* savDir = gDirectory;
1680  gROOT->cd();
1681  Double_t intg = rDeltaM->Integral();
1682  Double_t mintg = sDeltaP->Integral()+sDeltaS->Integral()+sDeltaC->Integral();
1683  TH1* dat = static_cast<TH1*>(rDeltaM->Clone("tmpM"));
1684  dat->SetDirectory(0);
1685  dat->Scale(1./intg);
1686  TH1* sig = static_cast<TH1*>(sDeltaP->Clone("tmpPS"));
1687  sig->SetDirectory(0);
1688  sig->Add(sDeltaS);
1689  sig->Scale(1./sig->Integral()); // mintg);
1690  TH1* bg = static_cast<TH1*>(sDeltaC->Clone("tmpC"));
1691  bg->SetDirectory(0);
1692  bg->Scale(1./bg->Integral()); // mintg);
1693 
1694  TObjArray mc;
1695  mc.SetOwner();
1696  mc.Add(sig);
1697  mc.Add(bg);
1698  // mc.Add(arr[3]);
1699  TFractionFitter f(dat, &mc, "Q");
1700  Int_t status = f.Fit();
1701  savDir->cd();
1702  if (status != 0) {
1703  Warning("FractionFit", "Fit failed w/status=%d", status);
1704  return 0;
1705  }
1706  Printf("\nTemplate fits");
1707  for (Int_t i = 0; i < 2; i++) {
1708  Double_t v, e;
1709  f.GetResult(i, v, e);
1710  Printf("%30s=%6.4f +/- %6.4f",
1711  mc.At(i)->GetName(), e, v);
1712  }
1713  TH1* ret = f.GetPlot();
1714  ret->Scale(intg);
1715  delete dat;
1716  return ret;
1717 }
1718 
1719 //____________________________________________________________________
1720 void AliTrackletdNdeta2::WriteDeltas(TDirectory* outDir,
1721  TH1* rDeltaM, TH1* rDeltaI,
1722  TH1* sDeltaM, TH1* sDeltaI,
1723  TH1* sDeltaC, TH1* sDeltaP,
1724  TH1* sDeltaS, TH1* sDeltaD,
1725  TH1* fit)
1726 {
1727  DebugGuard g(fProc&kDebug,1,"Writing Deltas to %s", outDir->GetName());
1728 
1729  THStack* all = new THStack("all", "");
1730  SetAttr(rDeltaM, kRed+2, 20, 1.0);
1731  SetAttr(rDeltaI, kOrange+2, 21, 1.0);
1732  rDeltaM->SetDirectory(outDir);
1733  rDeltaI->SetDirectory(outDir);
1734  all->Add(rDeltaM);
1735  all->Add(rDeltaI);
1736 
1737  TH1* toScale[] = {sDeltaM,sDeltaI,sDeltaC, sDeltaP,sDeltaS,sDeltaD,fit,0};
1738  Color_t toColor[] = {kRed, kOrange,kMagenta,kGreen, kBlue, kPink, kBlack};
1739  Style_t toStyle[] = {24, 25, 30, 26, 32, 30, 24 };
1740  TH1** pScale = toScale;
1741  Color_t* pColor = toColor;
1742  Style_t* pStyle = toStyle;
1743  while ((*pScale)) {
1744  (*pScale)->SetDirectory(outDir);
1745  all->Add((*pScale));
1746  SetAttr(*pScale, (*pColor)+2, *pStyle, 1.2);
1747  pScale++;
1748  pColor++;
1749  pStyle++;
1750  }
1751  outDir->cd();
1752  all->Write();
1753 
1754  THStack* ratios = new THStack("ratios", "");
1755  TH1* ratioM = static_cast<TH1*>(sDeltaM->Clone("ratioM"));
1756  ratioM->SetTitle("#Delta_{M'}/#Delta_{M}");
1757  ratioM->Divide(rDeltaM);
1758  ratioM->SetDirectory(outDir);
1759  ratios->Add(ratioM);
1760 
1761  TH1* ratioI = static_cast<TH1*>(sDeltaI->Clone("ratioI"));
1762  ratioI->SetTitle("#Delta_{I'}/#Delta_{I}");
1763  ratioI->Divide(rDeltaI);
1764  ratioI->SetDirectory(outDir);
1765  ratios->Add(ratioI);
1766 
1767  TH1* ratioIC = static_cast<TH1*>(sDeltaC->Clone("ratioIC"));
1768  ratioIC->SetTitle("#Delta_{C'}/#Delta_{I}");
1769  ratioIC->Divide(rDeltaI);
1770  ratioIC->SetDirectory(outDir);
1771  ratios->Add(ratioIC);
1772 
1773  if (!fit) { ratios->Write(); return; }
1774 
1775  TH1* ratioF = static_cast<TH1*>(fit->Clone("ratioF"));
1776  ratioF->SetTitle("#Delta_{fit}/#Delta_{M}");
1777  ratioF->Divide(rDeltaM);
1778  ratioF->SetDirectory(outDir);
1779  ratios->Add(ratioF);
1780 
1781  ratios->Write();
1782 }
1783 
1784 
1785 //====================================================================
1787  Container* simCont,
1788  Double_t realEff,
1789  Double_t simEff,
1790  TDirectory* outParent,
1791  Int_t deltaDimen)
1792 {
1793  DebugGuard g(fProc&kDebug,1,"Calculating results for %s",realCont->GetName());
1794 
1795  TDirectory* outDir = outParent->mkdir(Form("results%dd", deltaDimen));
1796  TDirectory* delDir = outParent->GetDirectory(Form("delta%dd", deltaDimen));
1797 
1798  Bool_t fk = (fProc & kFake);
1799  Bool_t sc = (fProc & kDTFS);
1800  Bool_t rs = fRealIsSim;
1801 
1802  TH2* scale = static_cast<TH2*>(delDir->Get("scale"));
1803  TH2* realM = CopyH2(GetC(realCont,"measured"), "etaIPz", "realM");
1804  TH2* realS = CopyH2(GetC(realCont,"measured"), "etaIPz", "realS");
1805  TH2* realC = fk ? CopyH2(GetC(realCont,"measured"), "etaIPz", "realC"):0;
1806  TH2* simM = CopyH2(GetC(simCont, "measured"), "etaIPz", "simM");
1807  TH2* simC = fk ? CopyH2(GetC(simCont, "combinatorics"),"etaIPz", "simC") :0;
1808  TH2* simG = CopyH2(GetC(simCont, "generated"), "etaIPz", "simG");
1809  TH2* simS = CopyH2(GetC(simCont, "measured"), "etaIPz", "simS");
1810  TH2* simA = CopyH2(GetC(simCont, "generated"), "etaIPz", "simA");
1811  TH2* simB = fk ? CopyH2(GetC(simCont, "combinatorics"),"etaIPz", "simB") :0;
1812  TH2* simT = sc ? CopyH2(GetC(simCont, "secondaries"), "dtfs", "simT") :0;
1813  TH2* realG = rs ? CopyH2(GetC(realCont,"generated"), "etaIPz", "realG"):0;
1814  TH1* realZ = CopyH1(realCont, "ipz", "realZ");
1815  TH1* simZ = CopyH1(simCont, "ipz", "simZ");
1816 
1817  Double_t realEt = GetD(realCont, "triggerEfficiency");
1818  Double_t simEt = GetD(simCont, "triggerEfficiency");
1819  Double_t realEz = GetD(realCont, "ipEfficiency");
1820  Double_t simEz = GetD(simCont, "ipEfficiency");
1821 
1822  if (scale->GetEntries() <= 0) {
1823  Warning("Results", "Scale has no entries");
1824  }
1825 
1826  // Scale combinatorial background to measured to get beta
1827  if (simB) {
1828  simB->Divide(simM);
1829  simB->SetTitle("#beta'");
1830  simB->SetZTitle("#beta'");
1831  }
1832 
1833  // Substract combinatorial background from measured to get signal
1834  if (simC) {
1835  simS->Add(simC, -1);
1836  simS->SetTitle("S'");
1837  simS->SetZTitle("S'");
1838  }
1839 
1840  // Possibly subtract secondaries from strange
1841  if (simT) {
1842  simT->SetTitle("T'");
1843  simT->SetZTitle("T'");
1844  simS->Add(simT, -1);
1845  }
1846 
1847  // Scale MC truth primaries by signal to get correction
1848  simA->Divide(simS);
1849  simA->SetTitle("A'");
1850  simA->SetZTitle("#alpha'");
1851 
1852  // Copy simulated beta to real beta, and scale by scalar
1853  TH2* realB = 0;
1854  if (simB) {
1855  realB = static_cast<TH2*>(simB->Clone("realB"));
1856  realB->SetDirectory(0);
1857  realB->SetTitle("#beta");
1858  realB->SetZTitle("#beta");
1859  realB->Multiply(scale);
1860  Scale(realB, fFudge, 0);
1861  }
1862 
1863  // Multiply real beta onto real measured to get background
1864  if (realC) {
1865  realC->Multiply(realB);
1866  realC->SetTitle("C");
1867  realC->SetZTitle("C");
1868  }
1869 
1870  // Substract the real background off the real measured
1871  if (realC) {
1872  realS->Add(realC, -1);
1873  realS->SetTitle("S");
1874  realS->SetZTitle("S");
1875  }
1876 
1877  // Possibly subtract secondaries from strange
1878  TH2* realT = 0;
1879  if (simT) {
1880  realT = static_cast<TH2*>(simT->Clone("realT"));
1881  realT->SetDirectory(0);
1882  realT->SetTitle("T");
1883  realT->SetZTitle("T");
1884  realT->Scale(fSEF);
1885  realS->Add(realT, -1);
1886  }
1887 
1888  // Make a fiducial distribution, and coerce the others to fit this
1889  TH2* fiducial = static_cast<TH2*>(simA->Clone("fiducial"));
1890  fiducial->SetTitle("F");
1891  fiducial->SetZTitle("F");
1892  for (Int_t i = 1; i <= fiducial->GetNbinsX(); i++) {
1893  for (Int_t j = 1; j <= fiducial->GetNbinsY(); j++) {
1894  Double_t c = fiducial->GetBinContent(i,j);
1895  fiducial->SetBinContent(i,j,c > fMinAlpha && c <= fMaxAlpha);
1896  fiducial->SetBinError(i,j,0);
1897  }
1898  }
1899  /* */ realM->Multiply(fiducial);
1900  if (realC) realC->Multiply(fiducial);
1901  /* */ realS->Multiply(fiducial);
1902  if (realB) realB->Multiply(fiducial);
1903  /* */ simM ->Multiply(fiducial);
1904  if (simC) simC ->Multiply(fiducial);
1905  /* */ simS ->Multiply(fiducial);
1906  if (simB) simB ->Multiply(fiducial);
1907  /* */ simA ->Multiply(fiducial);
1908  if (simT) simT ->Multiply(fiducial);
1909  if (realT) realT->Multiply(fiducial);
1910 
1911  // We can now make our result
1912  TH2* result = static_cast<TH2*>(realS->Clone("result"));
1913  result->Multiply(simA);
1914  result->SetTitle("R");
1915  result->SetZTitle("R");
1916 
1917  FlattenCentrality(result,realCont);
1918 
1919  // Output directory for full stuff
1920  TDirectory* full = outDir->mkdir("full");
1921 
1922  // Make a stack of 1D projections
1923  struct Rec {
1924  TH2* h; // histogram to average
1925  TH1* s; // IPz distribution
1926  Double_t e; // IP efficiency
1927  TH2* c; // Mask for bins to use
1928  TH1* p; // The result
1929  Style_t sty; // Marker style
1930  Color_t col; // Color
1931  Float_t siz; // Marker size
1932  const char* tit; // Title
1933  };
1934  Rec sT = { simT, simZ, simEz, 0, 0,30,kSpring+2, 1.4,"strangeness"};
1935  Rec sC = { simC, simZ, simEz, result,0,32,kMagenta+2,1.4,"background"};
1936  Rec sS = { simS, simZ, simEz, result,0,27,kGreen+2, 1.8,"signal" };
1937  Rec sM = { simM, simZ, simEz, result,0,26,kBlue+2, 1.4,"measured" };
1938  Rec sG = { simG, simZ, simEz, 0, 0,24,kRed+2, 1.4,"generated" };
1939  Rec sB = { simB, simZ, simEz, 0, 0,28,kPink+2, 1.4,"#beta" };
1940  Rec rC = { realC, realZ,realEz,result,0,23,kMagenta+2,1.2,"background"};
1941  Rec rT = { realT, realZ,realEz,0, 0,29,kSpring+2, 1.4,"strangeness"};
1942  Rec rS = { realS, realZ,realEz,result,0,33,kGreen+2, 1.6,"signal" };
1943  Rec rM = { realM, realZ,realEz,result,0,22,kBlue+2, 1.2,"measured" };
1944  Rec rR = { result, realZ,realEz,0, 0,20,kRed+2, 1.3,ObsTitle() };
1945  Rec rB = { realB, realZ,realEz,0, 0,34,kPink+2, 1.4,"#beta" };
1946  Rec sA = { simA, simZ, simEz, 0, 0,30,kSpring+2, 1.4,"#alpha" };
1947  Rec sF = { fiducial,simZ, simEz, 0, 0,31,kSpring+2, 1.4,"fiducial" };
1948  Rec rG = { realG, realZ,realEz,0, 0,24,kBlack, 1.4,"truth" };
1949  Rec* recs[]={ &rR, &sG, &rS, &sS, &rM, &sM, &rC, &sC,
1950  &rB, &sB, &sA, &sF, &rG, &rT, &sT, 0 };
1951  Rec** ptr = recs;
1952 
1953  TH1* dndeta = 0;
1954  THStack* all = new THStack("all", "");
1955  if (fViz & kAltMarker) {
1956  rR.sty = 21;
1957  rR.siz = 1.2;
1958  sG.sty = 25;
1959  sG.siz = 1.3;
1960  rG.sty = 25;
1961  rG.siz = 1.3;
1962  }
1963  while ((*ptr)) {
1964  Rec* src = *ptr;
1965  if (!src->h) { ptr++; continue; }
1966  src->h->SetDirectory(full);
1967  src->p = AverageOverIPz(src->h, src->h->GetName(), 1,
1968  src->s, src->e, src->c);
1969  src->p->SetYTitle(src->h->GetZaxis()->GetTitle());
1970  src->p->SetTitle(Form("%s - %s", src->h->GetTitle(), src->tit));
1971  src->p->SetDirectory(outDir);
1972  if (src->e > 1e-3) src->p->Scale(src->e);
1973  if (src->h != simB && src->h != realB &&
1974  src->h != simA && src->h != fiducial) all->Add(src->p);
1975  SetAttr(src->p, src->col, src->sty, src->siz);
1976  if (src->h == result) {
1977  dndeta = src->p;
1978  dndeta->SetYTitle(ObsTitle());
1979  }
1980  ptr++;
1981  }
1982  // Show example calculation
1983  Int_t mi = result->GetXaxis()->FindBin(0.);
1984  Int_t mj = result->GetYaxis()->FindBin(0.);
1985  TString simST;
1986  TString realST;
1987  if (simT) {
1988  simST.Form("-%4.1f", simT->GetBinContent(mi,mj));
1989  realST.Form("-%4.1f", realT->GetBinContent(mi,mj));
1990  }
1991  printf("R=G'/[(1-beta')M'](1-beta)M="
1992  "%6.1f /((1-%6.3f)*%6.1f%s)*(1-%6.3f)*(%6.1f%s)="
1993  "%6.3f * %6.3f="
1994  "%6.1f [%6.1f]",
1995  /* */ sG.h->GetBinContent(mi,mj),
1996  sB.h ? sB.h->GetBinContent(mi,mj) : 0,
1997  /* */ sM.h->GetBinContent(mi,mj),
1998  /* */ simST.Data(),
1999  rB.h ? rB.h->GetBinContent(mi,mj) : 0,
2000  /* */ rM.h->GetBinContent(mi,mj),
2001  /* */ realST.Data(),
2002  /* */ sA.h->GetBinContent(mi,mj),
2003  /* */ rS.h->GetBinContent(mi,mj),
2004  /* */ rR.h->GetBinContent(mi,mj),
2005  /* */ rG.h ? rG.h->GetBinContent(mi,mj) : -1);
2006 
2007  if (rG.p) {
2008  TH1* ratio = RatioH(dndeta, rG.p, "closure");
2009  ratio->SetYTitle("Closure test");
2010  ratio->SetDirectory(outDir);
2011  }
2012 
2013 
2014  THStack* ratios = new THStack("ratios", "");
2015  TH1* scaleC = 0;
2016  ratios->Add(RatioH(rM.p, sM.p, "rMeaured"));
2017  ratios->Add(RatioH(rC.p, sC.p, "rBackground"));
2018  ratios->Add(RatioH(rS.p, sS.p, "rSignal"));
2019  ratios->Add(RatioH(rR.p, sG.p, "rResult"));
2020  ratios->Add((scaleC = AverageOverIPz(scale,scale->GetName(),1,realZ,0,0)));
2021  scaleC->SetMarkerColor(kBlack);
2022  scaleC->SetMarkerColor(kGray);
2023  scaleC->SetMarkerStyle(31);
2024  scaleC->SetLineStyle(2);
2025  scaleC->SetLineColor(kGray);
2026 
2027  TF1* tmp = new TF1("mid", "pol0", -.5, +.5);
2028  dndeta->Fit(tmp, "Q0R+");
2029  TLatex* ltx = new TLatex(0,tmp->GetParameter(0)/2,
2030  Form("%s|_{|#eta|<0.5}=%.1f#pm%.1f",
2031  ObsTitle(),
2032  tmp->GetParameter(0),
2033  tmp->GetParError(0)));
2034 #if 0
2035  printf(" %dD: %6.1f +/- %4.1f (%4.2f)",
2036  deltaDimen, tmp->GetParameter(0), tmp->GetParError(0),
2037  tmp->GetChisquare()/tmp->GetNDF());
2038 #endif
2039 
2040  ltx->SetTextAlign(22);
2041  ltx->SetTextFont(42);
2042  dndeta->GetListOfFunctions()->Add(ltx);
2043 
2044  outDir->cd();
2045  all->Write();
2046  ratios->Write();
2047  if (realB) realB ->SetDirectory(full);
2048  if (simB) simB ->SetDirectory(full);
2049  /* */ simA ->SetDirectory(full);
2050  /* */ fiducial->SetDirectory(full);
2051 
2052  outParent->cd();
2053  dndeta->SetBinContent(0, realEff);
2054  dndeta->SetBinContent(dndeta->GetNbinsX()+1,simEff);
2055 
2056  return dndeta;
2057 }
2058 
2059 //____________________________________________________________________
2061 {
2062  if (!(fProc & kCentAverage)) return;
2063  DebugGuard g(fProc&kDebug,1,"Correcting for non-flat centrality in bin");
2064 
2065  TH1* centIpz = GetH1(c, "centIpz");
2066  if (!centIpz) {
2067  Warning("FlattenCentrality", "No average cent vs IPz data");
2068  return;
2069  }
2070  Double_t meanC = centIpz->GetMean(2);
2071  if (meanC < 0 || meanC > 100) {
2072  Warning("FlattenCentrality", "Average centrality %f funny", meanC);
2073  return;
2074  }
2075  for (Int_t iz = 1; iz <= centIpz->GetXaxis()->GetNbins(); iz++) {
2076  Double_t cc = centIpz->GetBinContent(iz);
2077  Double_t rc = meanC / cc;
2078 #if 0
2079  DebugGuard g2(fProc&kDebug,1,"<c>=%7.4f%% <ci>=%7.4f%% -> %6.4f",
2080  meanC, cc, rc);
2081 #endif
2082  for (Int_t ieta = 1; ieta <= h->GetXaxis()->GetNbins(); ieta++) {
2083  h->SetBinContent(ieta, iz, rc * h->GetBinContent(ieta, iz));
2084  h->SetBinError (ieta, iz, rc * h->GetBinError (ieta, iz));
2085  }
2086  }
2087 }
2088 
2089 
2090 //====================================================================
2092 {
2093  fTop->Clear();
2094  fTop->SetNumber(1);
2095  fTop->SetFillColor(kGray); // kTopBackground);
2096  fTop->SetFillStyle(1001);
2097  fTop->SetBorderSize(0);
2098  fTop->SetBorderMode(0);
2099 
2100  fBody->Clear();
2101  fBody->SetNumber(2);
2102  fBody->SetFillColor(0);
2103  fBody->SetFillStyle(0);
2104  fBody->SetBorderSize(0);
2105  fBody->SetBorderMode(0);
2106  fBody->SetTopMargin(0.01);
2107  fBody->SetRightMargin(0.01);
2108  fBody->SetBottomMargin(0.10);
2109  fBody->SetLeftMargin(0.10);
2110  fBody->SetTicks();
2111 
2112  fCanvas->cd();
2113 }
2114 //____________________________________________________________________
2116 {
2117  gStyle->SetOptStat(0);
2118  gStyle->SetOptTitle(0);
2119 
2120  Int_t h = 1000;
2121  Int_t w = h / TMath::Sqrt(2);
2122  if (fViz & kLandscape) {
2123  Int_t t = h;
2124  h = w;
2125  w = t;
2126  }
2127 
2128  fCanvas = new TCanvas("c",outputName,w,h);
2129  fCanvas->SetFillColor(0);
2130  fCanvas->SetBorderSize(0);
2131  fCanvas->SetBorderMode(0);
2132 
2133  if (fViz & kPDF) {
2134  SuppressGuard2 g;
2135  fCanvas->Print(Form("%s/summary.pdf[", outputName.Data()),
2136  Form("pdf %s", (fViz & kLandscape) ? "Landscape" : ""));
2137  }
2138  // if (fViz & kPNG) {
2139  // }
2140  fCanvas->SetLeftMargin (0.10);
2141  fCanvas->SetRightMargin (0.05);
2142  fCanvas->SetTopMargin (0.05);
2143  fCanvas->SetBottomMargin(0.10);
2144 
2145  Float_t dy = 0.05;
2146  fTop = new TPad("top","Top",0,1-dy,1,1,0,0);
2147  // fTop->SetNumber(1);
2148  // fTop->SetFillColor(kTopBackground);
2149  // fTop->SetBorderSize(0);
2150  // fTop->SetBorderMode(0);
2151  fCanvas->cd();
2152  fTop->Draw();
2153 
2154  fBody = new TPad("body","Body",0,0,1,1-dy,0,0);
2155  fBody->SetNumber(2);
2156  // fBody->SetFillColor(0);
2157  // fBody->SetFillStyle(0);
2158  // fBody->SetBorderSize(0);
2159  // fBody->SetBorderMode(0);
2160  fCanvas->cd();
2161  fBody->Draw();
2162 
2163  ClearCanvas();
2164 
2165  fHeader = new TLatex(.5, .5, "Title");
2166  fHeader->SetNDC();
2167  fHeader->SetTextAlign(22);
2168  // fHeader->SetTextColor(kWhite);
2169  fHeader->SetTextFont(62);
2170  fHeader->SetTextSize(0.7);
2171 
2172  fCanvas->cd();
2173 }
2174 //____________________________________________________________________
2176 {
2177  if ((fViz & kPDF) && fCanvas) {
2178  SuppressGuard2 g;
2179  fCanvas->Print(Form("%s/summary.pdf]", fCanvas->GetTitle()),
2180  Form("pdf %s Title:%s",
2181  (fViz & kLandscape) ? "Landscape" : "",
2182  fLastTitle.Data()));
2183  Printf("PDF %s written", fCanvas->GetTitle());
2184  }
2185  if (fCanvas) fCanvas->Close();
2186  fCanvas = 0;
2187 }
2188 
2189 //____________________________________________________________________
2191  const char* shortTitle,
2192  Float_t size)
2193 {
2194  if (fTop) {
2195  fTop->cd();
2196  fHeader->SetTextSize(size);
2197  fHeader->DrawLatex(.5,.5,title);
2198  }
2199  fCanvas->Modified();
2200  fCanvas->Update();
2201  fCanvas->cd();
2202 
2203  if (fViz & kPDF) {
2204  TString tit;
2205  tit.Form("pdf %s Title:%s", (fViz & kLandscape) ? "Landscape" : "",
2206  title);
2207  // Suppress prints
2208  SuppressGuard2 g;
2209  fCanvas->Print(Form("%s/summary.pdf",fCanvas->GetTitle()), tit);
2210  }
2211  static Int_t cnt = 1;
2212  if (fViz & kPNG) {
2213  SuppressGuard2 g;
2214  fCanvas->Print(Form("%s/%03d_%s.png",fCanvas->GetTitle(),cnt,shortTitle));
2215  cnt++;
2216  }
2217  fLastTitle = title;
2218  if (fViz & kPause) fCanvas->WaitPrimitive();
2219 }
2220 //____________________________________________________________________
2221 TLegend* AliTrackletdNdeta2::DrawInPad(TVirtualPad* c,
2222  Int_t pad,
2223  TObject* o,
2224  Option_t* opt)
2225 {
2226  if (!o) {
2227  // Warning("", "Nothing to draw in pad %d", pad);
2228  return 0;
2229  }
2230  TLegend* l = 0;
2231  TVirtualPad* p = c->cd(pad);
2232  if (!p) {
2233  Warning("DrawInPad", "Sub-pad %d does not exist", pad);
2234  return 0;
2235  }
2236  TString option(opt);
2237  option.ToLower();
2238  if (option.Contains("logx")) { p->SetLogx(); option.ReplaceAll("logx",""); }
2239  if (option.Contains("logy")) { p->SetLogy(); option.ReplaceAll("logy",""); }
2240  if (option.Contains("logz")) { p->SetLogz(); option.ReplaceAll("logz",""); }
2241  if (option.Contains("grid")) { p->SetGridx(); p->SetGridy();
2242  option.ReplaceAll("grid",""); }
2243  Int_t leg = 0;
2244  if (option.Contains("leg3")) { leg = 3; option.ReplaceAll("leg3",""); }
2245  if (option.Contains("leg2")) { leg = 2; option.ReplaceAll("leg2",""); }
2246  if (option.Contains("leg")) { leg = 1; option.ReplaceAll("leg",""); }
2247  // Printf("Drawing %p %s with %s", o, o->GetName(), option.Data());
2248  o->Draw(option);
2249  if (leg) {
2250  l = p->BuildLegend(0.5, 0.73, .98, .98);
2251  l->SetNColumns(leg);
2252  TObject* frame = 0;
2253  TIter next(l->GetListOfPrimitives());
2254  TLegendEntry* ent = 0;
2255  while ((ent = static_cast<TLegendEntry*>(next()))) {
2256  if (TString(ent->GetLabel()).EqualTo("frame")) frame = ent;
2257  }
2258  if (frame) l->GetListOfPrimitives()->Remove(frame);
2259  // l->GetListOfPrimitives()->Print();
2260  }
2261  p->Modified();
2262  // p->Update();
2263  // p->cd();
2264  return l;
2265 }
2266 //____________________________________________________________________
2267 void AliTrackletdNdeta2::ModLegend(TVirtualPad* p, TLegend* l,
2268  Double_t x1, Double_t y1,
2269  Double_t x2, Double_t y2)
2270 {
2271  if (!p || !l) return;
2272  Double_t px1 = p->GetX1();
2273  Double_t px2 = p->GetX2();
2274  Double_t py1 = p->GetY1();
2275  Double_t py2 = p->GetY2();
2276  l->SetX1(px1+(px2-px1)*x1);
2277  l->SetX2(px1+(px2-px1)*x2);
2278  l->SetY1(py1+(py2-py1)*y1);
2279  l->SetY2(py1+(py2-py1)*y2);
2280  p->Modified();
2281 }
2282 //____________________________________________________________________
2283 THStack* AliTrackletdNdeta2::Make2Stack(const char* name,
2284  const char* title,
2285  Container* realList,
2286  Container* simList,
2287  Option_t* realOpt,
2288  Option_t* simOpt)
2289 {
2290  TString nme(name);
2291  TH1* real = CopyH1(realList, name, Form("real%s",name));
2292  TH1* sim = CopyH1(simList, name, Form("sim%s",name));
2293  if (real->GetNbinsX() <= 1 && sim->GetNbinsX() <= 1) return 0;
2294  THStack* stack = new THStack(name, title);
2295  real->SetMarkerStyle(20);
2296  sim ->SetMarkerStyle(24);
2297  real->SetFillStyle(3004);
2298  sim ->SetFillStyle(3005);
2299  real->SetBarWidth(0.4);
2300  sim ->SetBarWidth(0.4);
2301  real->SetBarOffset(0.1);
2302  sim ->SetBarOffset(0.5);
2303  TString dtit(real->GetTitle());
2304  if (dtit.Contains("\\")) dtit.Form("%s\\hbox{ - real}", real->GetTitle());
2305  else dtit.Form("%s - real", real->GetTitle());
2306  real->SetTitle(dtit);
2307  TString stit(sim->GetTitle());
2308  if (stit.Contains("\\")) stit.Form("%s\\hbox{ - sim.}", sim->GetTitle());
2309  else stit.Form("%s - sim.", sim->GetTitle());
2310  sim->SetTitle(stit);
2311  stack->Add(real, realOpt);
2312  stack->Add(sim, simOpt);
2313  return stack;
2314 }
2315 
2316 
2317 //====================================================================
2319  Container* simSums,
2320  Container* realRess,
2321  Container* simRess,
2322  TDirectory* outDir,
2323  Int_t maxBins)
2324 {
2325  DebugGuard g(fProc&kDebug,1,"Visualizing");
2326  // --- Visualization -----------------------------------------------
2327  TH1* realCent = GetH1(outDir, "realCent");
2328  if (!realCent) {
2329  Warning("Visualize", "realCent histogram not found");
2330  return false;
2331  }
2332  TString outName(gSystem->DirName(outDir->GetName()));
2333  outName.ReplaceAll(".root", "");
2334  CreateCanvas(outName);
2335  if (fViz & kParameters) VisualizeParams(realSums, simSums);
2336  if (fViz & kGeneral) VisualizeGeneral(realRess, simRess);
2337  if (fViz & kWeights) VisualizeWeights(simRess);
2338 
2339  if (fViz & kdNdetas) {
2340  for (Int_t i = 0; i < 4; i++) {
2341  if ((fProc & (1 << i)) == 0) continue;
2342  VisualizeFinal(outDir, i);
2343  VisualizeClosure(outDir, i);
2344  }
2345  }
2346  for (Int_t i = 1; i <= realCent->GetNbinsX() && i <= maxBins ; i++) {
2347  Double_t c1 = realCent->GetXaxis()->GetBinLowEdge(i);
2348  Double_t c2 = realCent->GetXaxis()->GetBinUpEdge (i);
2349 
2350  VisualizeBin(c1, c2, simRess, outDir);
2351  }
2352  CloseCanvas();
2353 }
2354 
2355 //====================================================================
2357  Container* simList)
2358 {
2359  DebugGuard g(fProc&kDebug,1,"Visualizing general things");
2360  THStack* ipz = Make2Stack("ipz", "IP_{#it{z}}", realList,simList);
2361  THStack* cent = Make2Stack("cent", "Centrality [%]", realList,simList);
2362  THStack* status = Make2Stack("status","Task status", realList,simList,
2363  "B text90", "B text90");
2364  THStack* centTr = Make2Stack("centTracklets", "#LTtracklets#GT",
2365  realList, simList, "E", "E");
2366  ClearCanvas();
2367  fBody->Divide(1,3);
2368  for (Int_t i = 1; i <= 3; i++) {
2369  if (i < 3) fBody->GetPad(i)->SetRightMargin(0.01);
2370  fBody->GetPad(i)->SetTopMargin(0.01);
2371  }
2372  TLegend* l = 0;
2373  TVirtualPad* q = fBody->GetPad(1);
2374  q->Divide(2,1);
2375  l = DrawInPad(q,1,ipz, "nostack leg");
2376  ModLegend(q->GetPad(1),l,.4,.1,.75,.4);
2377  l = DrawInPad(q,2,cent, "nostack leg");
2378  ModLegend(q->GetPad(2),l,.6,.1,.99,.4);
2379  q = fBody->GetPad(2);
2380  q->Divide(2,1);
2381  l = DrawInPad(q,1,status, "nostack hist text90 leg");
2382  ModLegend(q->GetPad(1),l,.5,.7,.99,.99);
2383  l = DrawInPad(q,2,centTr, "nostack leg");
2384  ModLegend(q->GetPad(2),l,.5,.7,.99,.99);
2385 
2386  TH2* real = GetH2(realList, "etaPhi");
2387  TH2* sim = GetH2(simList, "etaPhi");
2388  if (sim) {
2389  sim->SetMarkerColor(kBlack);
2390  sim->SetMarkerStyle(0);
2391  sim->SetMarkerSize(1);
2392  sim->SetLineColor(kBlack);
2393  sim->SetFillColor(kBlack);
2394  sim->SetFillStyle(0);
2395  sim->SetName("simEtaPhi");
2396  }
2397  DrawInPad(fBody, 3, real, "colz");
2398  DrawInPad(fBody, 3, sim, "box same");
2399 
2400  PrintCanvas("General information","general");
2401 }
2402 
2403 namespace {
2404  void SetCentColors(THStack* s, TH1* dist=0)
2405  {
2406  if (!s || !s->GetHists()) return;
2407 
2408  TIter next(s->GetHists());
2409  TH1* h = 0;
2410  Int_t i = 0;
2411  Double_t min = +10000;
2412  Double_t max = -10000;
2413  while ((h = static_cast<TH1*>(next()))) {
2414  Color_t c = AliTrackletAODUtils::CentColor(i);
2415  h->SetMarkerColor(c);
2416  h->SetFillColor(c);
2417  h->SetLineColor(c);
2418  h->SetMarkerStyle(20+(i%4));
2419  h->SetDirectory(0);
2420  min = TMath::Min(h->GetMinimum(), min);
2421  max = TMath::Max(h->GetMaximum(), max);
2422  i++;
2423  if (!dist) continue;
2424  h->SetTitle(Form("%5.1f-%5.1f%%",
2425  dist->GetXaxis()->GetBinLowEdge(i),
2426  dist->GetXaxis()->GetBinUpEdge(i)));
2427 
2428  }
2429  s->SetMinimum(min*.9);
2430  s->SetMaximum(max*1.1);
2431  }
2432  THStack* GetPdgStack(AliTrackletAODUtils::Container* w, const char* name)
2433  {
2435  if (!c) return 0;
2436 
2437  THStack* s = new THStack(name, "");
2438  TH1* h = 0;
2439  TIter n(c);
2440  while ((h = static_cast<TH1*>(n()))) s->Add(h);
2441 
2442  if (!s->GetHists() || s->GetHists()->GetEntries() <= 0) {
2443  delete s;
2444  s = 0;
2445  }
2446  return s;
2447  }
2448 }
2449 //____________________________________________________________________
2451 {
2452  DebugGuard g(fProc&kDebug,1,"Visualizing weights");
2453  TH1* c = GetH1(simList, "cent");
2454  Container* w = GetC(simList,"weights", false);
2455  if (!w) return;
2456 
2457  Float_t right = .75;
2458  ClearCanvas();
2459  fBody->SetTopMargin(0.01);
2460  fBody->Divide(1,3);
2461  TVirtualPad* pp[] = { fBody->GetPad(1), fBody->GetPad(2) };
2462  for (Int_t i = 0; i < 2; i++) {
2463  pp[i]->SetRightMargin(0.01);
2464  pp[i]->SetTopMargin(0.01);
2465  pp[i]->SetPad(pp[i]->GetXlowNDC(), pp[i]->GetYlowNDC(),
2466  right, pp[i]->GetYlowNDC()+pp[i]->GetHNDC());
2467  pp[i]->Modified();
2468  }
2469  TH2* hp = GetH2(w, "centPt");
2470  THStack* ef = new THStack(GetP2(simList,"etaWeight"),"x","effWeights","");
2471  THStack* ab = GetPdgStack(w, "abundance");
2472  THStack* st = GetPdgStack(w, "strangeness");
2473  SetCentColors(ef, c);
2474  SetCentColors(ab);
2475  SetCentColors(st);
2476 
2477  Double_t eMin = +1e9;
2478  Double_t eMax = -1e9;
2479  TIter next(ef->GetHists());
2480  TH1* h = 0;
2481  while ((h = static_cast<TH1*>(next()))) {
2482  eMin = TMath::Min(h->GetMinimum(), eMin);
2483  eMax = TMath::Max(h->GetMaximum(), eMax);
2484  }
2485  ef->SetMinimum(eMax);
2486  ef->SetMinimum(eMin);
2487  // ef->SetMaximum(1.02);
2488  TLegend* l = DrawInPad(fBody, 1, ef, "nostack p leg");
2489  ef->GetHistogram()->SetYTitle("Average weight");
2490  ef->GetHistogram()->SetXTitle("#eta");
2491  fBody->GetPad(1)->GetListOfPrimitives()->Remove(l);
2492  fBody->GetPad(1)->Modified();
2493 
2494  if (hp) {
2495  if (hp->GetNbinsY() > 1) {
2496  THStack* pt = new THStack(hp, "y","pt","");
2497  SetCentColors(pt);
2498  DrawInPad(fBody, 2, pt, "nostack");
2499  pt->GetHistogram()->SetYTitle("Weight");
2500  pt->GetHistogram()->SetXTitle("#it{p}_{T}");
2501  }
2502  else {
2503  TArrayD bins(hp->GetNbinsX()+1);
2504  bins[0] = hp->GetXaxis()->GetBinLowEdge(1);
2505  for (Int_t i = 1; i <= hp->GetNbinsX(); i++)
2506  bins[i] = hp->GetXaxis()->GetBinUpEdge(1);
2507  TH1* pt = new TH1D("pt","", bins.GetSize()-1,bins.GetArray());
2508  for (Int_t i = 1; i <= hp->GetNbinsX(); i++) {
2509  pt->SetBinContent(i, hp->GetBinContent(i,1));
2510  pt->SetBinError (i, hp->GetBinError (i,1));
2511  }
2512  pt->SetYTitle("Weight");
2513  pt->SetXTitle("Centrality [%]");
2514  pt->SetMarkerStyle(2);
2515  pt->SetMarkerColor(kRed+2);
2516  DrawInPad(fBody, 2, pt, "");
2517  }
2518  }
2519 
2520  fBody->cd();
2521  if (l) {
2522  l->Draw();
2523  ModLegend(fBody, l, right, pp[1]->GetYlowNDC(),
2524  .99, 1-fBody->GetTopMargin());
2525  fBody->Modified();
2526  }
2527 
2528  TVirtualPad* p3 = fBody->GetPad(3);
2529  p3->SetTopMargin(0.01);
2530  p3->SetRightMargin(0.01);
2531  p3->Divide(2,1);
2532  p3->GetPad(1)->SetRightMargin(0.01);
2533  p3->GetPad(2)->SetRightMargin(0.01);
2534  p3->GetPad(1)->SetTopMargin(0.01);
2535  p3->GetPad(2)->SetTopMargin(0.01);
2536 
2537  DrawInPad(p3, 1, ab, "nostack leg");
2538  DrawInPad(p3, 2, st, "nostack leg");
2539 
2540  if (ab) {
2541  ab->GetHistogram()->SetYTitle("Weight");
2542  ab->GetHistogram()->SetXTitle("Centrality [%]");
2543  }
2544  if (st) {
2545  st->GetHistogram()->SetYTitle("Weight");
2546  st->GetHistogram()->SetXTitle("Centrality [%]");
2547  }
2548  p3->GetPad(1)->Modified();
2549  p3->GetPad(2)->Modified();
2550 
2551  PrintCanvas("Simulation weights","weights");
2552 }
2553 
2554 //____________________________________________________________________
2556  Double_t& min,
2557  Double_t& max)
2558 {
2559  if (!stack->GetHists()) return;
2560  min = +1e6;
2561  max = -1e6;
2562  TIter next(stack->GetHists());
2563  TH1* hist = 0;
2564  while ((hist = static_cast<TH1*>(next()))) {
2565  if (hist->GetEntries() <= 0) continue;
2566  min = TMath::Min(min, hist->GetMinimum());
2567  max = TMath::Max(max, hist->GetMaximum());
2568  }
2569  stack->SetMinimum(min);
2570  stack->SetMaximum(max);
2571 }
2572 
2573 
2574 
2575 //____________________________________________________________________
2576 void AliTrackletdNdeta2::VisualizeFinal(TDirectory* outDir, Int_t i)
2577 {
2578  DebugGuard g(fProc&kDebug,1,"Visualizing final @ %d", i);
2579  TDirectory* dd = outDir->GetDirectory(Form("final%dd", i));
2580  if (!dd) return;
2581 
2582  THStack* all = GetHS(dd, "full");
2583  TH1* mid = GetH1(dd, "mid");
2584  TH1* pub = GetH1(outDir, "published");
2585  if (!mid) return;
2586 
2587  Double_t min, max;
2588  StackMinMax(all, min, max);
2589  if (all->GetHists() && all->GetHists()->GetEntries() > 1 && min <= 0) min=2;
2590  if (mid->GetMinimum() <= 0) mid->SetMinimum(min);
2591  if (fProc&kDebug) all->GetHists()->Print();
2592 
2593  max = TMath::Max(max,mid->GetMaximum());
2594  min = TMath::Min(min,mid->GetMinimum());
2595  all->SetMinimum(.9*min);
2596  all->SetMaximum(1.2*max);
2597  mid->SetMinimum(.9*min);
2598  mid->SetMaximum(1.2*max);
2599  mid->SetLineColor(kBlack);
2600  mid->SetFillColor(kBlack);
2601  mid->SetMarkerColor(kBlack);
2602  ClearCanvas();
2603 
2604  TPad* p1 = new TPad("p1","p1",0,0,.4,1);
2605  p1->SetTopMargin(0.01);
2606  p1->SetRightMargin(0.0);
2607  p1->SetLeftMargin(0.12);
2608  p1->SetBottomMargin(0.15);
2609  p1->SetTicks();
2610  fBody->cd();
2611  p1->Draw();
2612  p1->SetNumber(2);
2613 
2614  Double_t right = .7;
2615  TPad* p2 = new TPad("p2","p2",.4,0,1,1);
2616  p2->SetTopMargin(0.01);
2617  p2->SetRightMargin(1-right);
2618  p2->SetLeftMargin(0.0);
2619  p2->SetBottomMargin(0.15);
2620  p2->SetTicks();
2621  fBody->cd();
2622  p2->Draw();
2623  p2->SetNumber(2);
2624 
2625 
2626  DrawInPad(p1,0, mid, "logy grid");
2627  TLegend* l = DrawInPad(p1,0, pub, "logy grid same leg");
2628  ModLegend(p1, l, .4, .90, .99, .99);
2629  l->SetMargin(0.2);
2630  l->SetEntrySeparation(0.1);
2631  l->SetTextSize(0.04);
2632 
2633  l = DrawInPad(p2,0, all, "nostack logy grid leg");
2634  if (all && all->GetHistogram()) all->GetHistogram()->SetXTitle("#eta");
2635 
2636  ModLegend(p2, l, right, .15, .99, .99);
2637  l->SetMargin(0.2);
2638  l->SetEntrySeparation(0.1);
2639  l->SetTextSize(0.04);
2640  p1->Modified();
2641  p2->Modified();
2642  fBody->Modified();
2643 
2644  const char* what = (i == 3 ? "d^{3}N/(d#Deltad#etadIP_{z})" :
2645  i == 2 ? "d^{2}N/(d#Deltad#eta)" :
2646  i == 1 ? "dN/d#Delta" : "dN/d#Delta (k#equiv1)");
2647  PrintCanvas(Form("Results #topbar %s", what), "results");
2648 }
2649 
2650 //____________________________________________________________________
2651 void AliTrackletdNdeta2::VisualizeClosure(TDirectory* outDir, Int_t i)
2652 {
2653  DebugGuard g(fProc&kDebug,1,"Visualizing closure @ %d", i);
2654  TDirectory* dd = outDir->GetDirectory(Form("final%dd", i));
2655  if (!dd) return;
2656 
2657  THStack* all = GetHS(dd, "closure", false);
2658  if (!all) return;
2659 
2660  Double_t max = all->GetMaximum("nostack");
2661  Double_t min = all->GetMinimum("nostack");
2662  all->SetMinimum(.95*min);
2663  all->SetMaximum(1.05*max);
2664 
2665  ClearCanvas();
2666  Double_t right = .3;
2667  fBody->SetRightMargin(right);
2668  TLegend* l = DrawInPad(fBody,0, all, "grid nostack leg");
2669  ModLegend(fBody, l, 1-right+.01, fBody->GetBottomMargin(), .99,
2670  1-fBody->GetTopMargin());
2671  l->SetMargin(0.2);
2672  l->SetEntrySeparation(0.1);
2673  l->SetTextSize(0.04);
2674  l->SetBorderSize(0);
2675  fBody->Modified();
2676 
2677  const char* what = (i == 3 ? "d^{3}N/(d#Deltad#etadIP_{z})" :
2678  i == 2 ? "d^{2}N/(d#Deltad#eta)" :
2679  i == 1 ? "dN/d#Delta" : "dN/d#Delta (k#equiv1)");
2680  PrintCanvas(Form("Closure #topbar %s", what), "closure");
2681 }
2682 
2683 //====================================================================
2685  Double_t& y,
2686  const char* val)
2687 {
2688  TLatex* ln = new TLatex(.49, y, name);
2689  ln->SetTextAlign(31);
2690  ln->SetTextSize(0.02/gPad->GetHNDC());
2691  ln->SetNDC();
2692  ln->Draw();
2693  TLatex* lv = new TLatex(.51, y, val);
2694  lv->SetTextAlign(11);
2695  lv->SetTextSize(0.02/gPad->GetHNDC());
2696  lv->SetNDC();
2697  lv->Draw();
2698  y -= 0.025/gPad->GetHNDC();
2699 }
2700 //____________________________________________________________________
2702  Double_t& y,
2703  Double_t val)
2704 {
2705  VisualizeParam(name, y, Form("%f", val));
2706 }
2707 //____________________________________________________________________
2709  Double_t& y,
2710  Int_t val)
2711 {
2712  VisualizeParam(name, y, Form("%d", val));
2713 }
2714 //____________________________________________________________________
2716  Double_t& y,
2717  Bool_t val)
2718 {
2719  VisualizeParam(name, y, val ? "yes" : "no");
2720 }
2721 //____________________________________________________________________
2723  const char* title)
2724 {
2725  // c->Clear();
2726  Double_t y = .9;
2727  TLatex* latex = new TLatex(.5, y, title);
2728  latex->SetTextAlign(21);
2729  latex->SetTextSize(0.023/gPad->GetHNDC());
2730  latex->SetNDC();
2731  latex->Draw();
2732  y -= 0.028/gPad->GetHNDC();
2733  if (!pars) return;
2734  VisualizeParam("#delta#phi shift", y, GetD(pars, "DPhiShift"));
2735  VisualizeParam("Shifted #delta#phi cut",y, GetD(pars, "ShiftedDPhiCut"));
2736  VisualizeParam("#Delta cut", y, GetD(pars, "DeltaCut"));
2737  VisualizeParam("max#Delta", y, GetD(pars, "MaxDelta"));
2738  VisualizeParam("min#Delta_{tail}", y, GetD(pars, "TailDelta"));
2739  VisualizeParam("max#Delta_{tail}", y, GetD(pars, "TailMax"));
2740  VisualizeParam("abs.min#it{c}", y, GetD(pars, "AbsMinCent"));
2741 }
2742 //____________________________________________________________________
2744  Container* simSums)
2745 {
2746  DebugGuard g(fProc&kDebug,1,"Visualizing parameters");
2747  ClearCanvas();
2748  fBody->Divide((fViz & kLandscape ? 3 : 1),
2749  (fViz & kLandscape ? 1 : 3), 0, 0);
2750 
2751  TVirtualPad* p1 = fBody->GetPad(1);
2752  TVirtualPad* p2 = fBody->GetPad(2);
2753  TVirtualPad* p3 = fBody->GetPad(3);
2754 #if 0
2755  if (!(fViz & kLandscape)) {
2756  Double_t yr = .2;
2757  p1->SetPad(0,1-yr,1,1);
2758  p2->SetPad(0,(1-yr)/2,1,1-yr);
2759  p3->SetPad(0,0,1,(1-yr)/2);
2760  }
2761 #endif
2762 
2763  // Post-processing stuff
2764  fBody->cd(1);
2765  Double_t y = .80;
2766  TLatex* latex = new TLatex(.5, y, "Post-processing");
2767  latex->SetTextAlign(21);
2768  latex->SetTextSize(0.023/p1->GetHNDC());
2769  latex->SetNDC();
2770  latex->Draw();
2771  y -= 0.028/p1->GetHNDC();
2772  TString scaleM("");
2773  if (fProc & kUnitScale) scaleM.Append(" unit");
2774  if (fProc & kAverageScale) scaleM.Append(" const");
2775  if (fProc & kEtaScale) scaleM.Append(" dN/d#eta");
2776  if (fProc & kFullScale) scaleM.Append(" d^{2}N/(d#etadIP_{z})");
2777  TString tmp = scaleM.Strip(TString::kBoth, ' ');
2778  VisualizeParam("Scaling of comb. bg.", y, tmp);
2779  VisualizeParam("min#alpha", y, fMinAlpha);
2780  VisualizeParam("max#alpha", y, fMaxAlpha);
2781  VisualizeParam("min#it{k}", y, fMinK);
2782  VisualizeParam("max#it{k}", y, fMaxK);
2783  VisualizeParam("Flatten c", y, (fProc & kCentAverage) != 0);
2784  if (fFudge != 1) VisualizeParam("Fudge", y, fFudge);
2785 
2786  // From tasks
2787  fBody->cd(2);
2788  VisualizeParams(GetC(realSums, "parameters"), "Real data");
2789  fBody->cd(3);
2790  VisualizeParams(GetC(simSums, "parameters"), "Simulated data");
2791  PrintCanvas("Parameters", "parameters");
2792 }
2793 
2794 //====================================================================
2796  Double_t c2,
2797  Container* simList,
2798  TDirectory* outTop)
2799 {
2800  // Form the folder name
2801  TString centName(CentName(c1,c2));
2802  if (TMath::Abs(c1 - c2) < 1e-6) centName = "all";
2803  DebugGuard g(fProc&kDebug,1,"Visualizing bin %s", centName.Data());
2804 
2805  fLastBin.Form("%.1f#minus%.1f%%", c1, c2);
2806  fLastShort = centName;
2807 
2808  TDirectory* outDir = outTop->GetDirectory(centName);
2809  if (!outDir) {
2810  Warning("VisualizeBin", "Directory %s not found in %s",
2811  centName.Data(), outTop->GetName());
2812  return false;
2813  }
2814 
2815  Printf("%5.1f - %5.1f%%", c1, c2);
2816  if (fViz & kSpecies) {
2817  VisualizeSpecies(GetC(simList, centName));
2818  if (fViz & kDeltas)
2819  VisualizeSpeciesDelta(GetC(simList, centName), outDir);
2820  VisualizePrimary(GetC(GetC(simList, centName), "generated"));
2821  }
2822  if (fViz & kdNdetas) {
2823  for (Int_t i = 0; i < 4; i++) {
2824  if ((fProc & (1 << i)) == 0) continue;
2825  if (fViz & kDeltas) VisualizeDelta(outDir, i);
2826  if (fViz & kDetails) VisualizeDetails(outDir, i);
2827  if (fViz & kdNdetas) VisualizeResult(outDir, i);
2828  }
2829  }
2830  return true;
2831 }
2832 //____________________________________________________________________
2834 {
2835  DebugGuard g(fProc&kDebug,1,"Visualizing species");
2836  if (!simCont) return true;
2837  ClearCanvas();
2838  Container* meas = GetC(simCont, "measured");
2839  Container* comb = GetC(simCont, "combinatorics");
2840  Container* prim = GetC(simCont, "primaries");
2841  Container* seco = GetC(simCont, "secondaries");
2842  Container* cont[] = { meas, comb, prim, seco };
2843  const char* tit[] = { "M' by primary mother's specie",
2844  "C' by primary mother's specie",
2845  "P' by species",
2846  "S' by primary mother's specie" };
2847  fBody->SetTopMargin(0.01);
2848  fBody->SetRightMargin(0.01);
2849  fBody->Divide(2,2,0.01,0.01);
2850 
2851  for (Int_t i = 0; i < 4; i++) {
2852  if (!cont[i]) continue;
2853 
2854  Container* species = GetC(cont[i], "types");
2855  if (!species) continue;
2856 
2857  THStack* all = GetHS(species, "all");
2858  THStack* toPion = GetHS(species, "toPion");
2859 
2860 
2861  TVirtualPad* p = fBody->GetPad(i+1);
2862  p->SetTopMargin(0.10);
2863  p->SetRightMargin(0.15);
2864  p->SetBottomMargin(0.15);
2865  p->cd();
2866  TLatex* ltx = new TLatex(.5, .99, tit[i]);
2867  ltx->SetTextAlign(23);
2868  ltx->SetTextSize(0.04);
2869  ltx->SetTextFont(42);
2870  ltx->Draw();
2871 
2872  p->Divide(1,2,0,0);
2873  TLegend* l = DrawInPad(p, 1, all, "nostack grid leg");
2874  DrawInPad(p, 2, toPion, "nostack grid");
2875  all->GetHistogram()->SetYTitle("dN_{X}/d#eta");
2876  all->GetHistogram()->SetXTitle("#eta");
2877  all->GetHistogram()->GetYaxis()->SetTitleSize(0.08);
2878  all->GetHistogram()->GetYaxis()->SetTitleOffset(0.6);
2879  all->GetHistogram()->GetYaxis()->SetLabelSize(0.08);
2880  toPion->GetHistogram()->SetYTitle("Relative to #pi^{#pm} mothers");
2881  toPion->GetHistogram()->SetXTitle("#eta");
2882  toPion->GetHistogram()->GetYaxis()->SetTitleSize(0.08);
2883  toPion->GetHistogram()->GetYaxis()->SetTitleOffset(0.6);
2884  toPion->GetHistogram()->GetYaxis()->SetLabelSize(0.08);
2885  toPion->GetHistogram()->GetXaxis()->SetTitleSize(0.08);
2886  toPion->GetHistogram()->GetXaxis()->SetTitleOffset(0.6);
2887  toPion->GetHistogram()->GetXaxis()->SetLabelSize(0.08);
2888  p->GetPad(1)->GetListOfPrimitives()->Remove(l);
2889  p->GetPad(1)->Modified();
2890  p->GetPad(2)->Modified();
2891  p->cd();
2892  l->Draw();
2893 
2894  ModLegend(p, l, .85, p->GetBottomMargin(), .99, 1-p->GetTopMargin());
2895  p->Modified();
2896  }
2897 
2898  PrintCanvas(Form("Species #topbar %s", fLastBin.Data()),
2899  Form("%s_species", fLastShort.Data()));
2900  return true;
2901 }
2902 
2903 //____________________________________________________________________
2905  TDirectory* outDir)
2906 {
2907  DebugGuard g(fProc&kDebug,1,"Visualizing species deltas");
2908  if (!simCont) return true;
2909  ClearCanvas();
2910  Container* comb = GetC(GetC(simCont, "combinatorics"), "specie", false);
2911  Container* prim = GetC(GetC(simCont, "primaries"), "specie", false);
2912  Container* seco = GetC(GetC(simCont, "secondaries"), "specie", false);
2913  Container* cont[] = { comb, prim, seco, 0 };
2914  if (!comb || !prim || !seco) return true;
2915 
2916  const char* tit[] = { "C' by primary mother's specie",
2917  "P' by species",
2918  "S' by primary mother's specie" };
2919 
2920  // --- Draw the delta distributions --------------------------------
2921  THStack* hp = 0;
2922  THStack* hs = 0;
2923  THStack* hc = 0;
2924  Double_t rr = .8;
2925  for (Int_t i = 0; i < 2; i++) {
2926  TString sub = (i == 0 ? "mid" : "fwd");
2927  TPad* pad = new TPad("pad","pad",i*rr/2,0,(i+1)*rr/2,1);
2928  pad->SetTopMargin(0.10);
2929  pad->SetRightMargin(0.01);
2930  pad->SetNumber(i+1);
2931  fBody->cd();
2932  pad->Draw();
2933  pad->Divide(1,3,0,0);
2934  pad->GetPad(1)->SetRightMargin(0.01);
2935  pad->GetPad(2)->SetRightMargin(0.01);
2936  pad->GetPad(3)->SetRightMargin(0.01);
2937  DrawInPad(pad, 1, hp = GetHS(GetC(prim, sub), "all"),
2938  "nostack grid logy logx");
2939  DrawInPad(pad, 2, hs = GetHS(GetC(seco, sub), "all"),
2940  "nostack grid logy logx");
2941  DrawInPad(pad, 3, hc = GetHS(GetC(comb, sub), "all"),
2942  "nostack grid logy logx");
2943  hp->GetHistogram()->SetYTitle("Primaries");
2944  hs->GetHistogram()->SetYTitle("Secondaries");
2945  hc->GetHistogram()->SetYTitle("Combinatorial");
2946  hp->GetHistogram()->SetXTitle("#Delta");
2947  hs->GetHistogram()->SetXTitle("#Delta");
2948  hc->GetHistogram()->SetXTitle("#Delta");
2949  TLatex* txt = new TLatex(pad->GetLeftMargin()+
2950  (1-pad->GetLeftMargin()-pad->GetRightMargin())/2,
2951  1-pad->GetTopMargin()+.01, hp->GetTitle());
2952  txt->SetNDC();
2953  txt->SetTextAlign(21);
2954  txt->SetTextSize(0.07);
2955  DrawInPad(pad,0,txt,"");
2956  }
2957  // --- Make a legend -----------------------------------------------
2958  const TAxis& pdgs = PdgAxis();
2959  TLegend* l = new TLegend(rr, fBody->GetBottomMargin(),
2960  .99, 1-fBody->GetTopMargin());
2961  l->SetFillStyle(0);
2962  l->SetBorderSize(0);
2963  TLegendEntry* e = 0;
2964  for (Int_t i = 0; i < 3; i++) {
2965  THStack* tmp = i == 0 ? hs : i == 1 ? hc : hp;
2966  TIter next(tmp->GetHists());
2967  TH1* h = 0;
2968  while ((h = static_cast<TH1*>(next()))) {
2969  Int_t bin = h->GetUniqueID();
2970  TString lbl = pdgs.GetBinLabel(bin);
2971  Int_t pdg = lbl.Atoi();
2972  TString nme;
2973  Color_t col;
2974  Style_t sty;
2975  PdgAttr(pdg, nme, col, sty);
2976  if (nme.EqualTo("0")) {
2977  nme = "All";
2978  sty = 24;
2979  col = kBlack;
2980  h->SetMarkerStyle(sty);
2981  }
2982  nme.Append("#rightarrowX_{ch}");
2983  TIter nextE(l->GetListOfPrimitives());
2984  while ((e = static_cast<TLegendEntry*>(nextE()))) {
2985  if (nme.EqualTo(e->GetLabel())) break;
2986  }
2987  if (e) continue;
2988  e = l->AddEntry(lbl,nme,"p");
2989  e->SetMarkerStyle(sty);
2990  e->SetMarkerColor(col);
2991  e->SetLineColor (col);
2992  }
2993  }
2994  fBody->cd();
2995  l->Draw();
2996 
2997  PrintCanvas(Form("Species #Delta #topbar %s", fLastBin.Data()),
2998  Form("%s_species_delta", fLastShort.Data()));
2999 
3000  // --- Draw relative contributions to signal/background ------------
3001  ClearCanvas();
3002  fBody->SetTopMargin(0.2);
3003  fBody->Divide(2,3,0,0);
3004  gStyle->SetPaintTextFormat("9.7f");
3005  Bool_t drawLog = true;
3006  for (Int_t i = 0; i < 2; i++) {
3007  TString sub = (i == 0 ? "mid" : "fwd");
3008  for (Int_t j = 0; j < 3; j++) {
3009  Container* par = (j == 0 ? prim :
3010  j == 1 ? seco : comb);
3011  TString parT = (j == 0 ? "Primaries" :
3012  j == 1 ? "Secondaries" : "Combinatorics");
3013  Int_t padNo = j * 2 + i + 1;
3014  TVirtualPad* pad = fBody->cd(padNo);
3015  if ((padNo % 2) == 0) pad->SetRightMargin(0.01);
3016  pad->SetTicks();
3017  pad->SetGridx();
3018  pad->SetGridy();
3019  THStack* rhs = GetHS(GetC(par, sub), "ratios");
3020  TObjLink* ptr = rhs->GetHists()->FirstLink();
3021  while (ptr) {
3022  ptr->SetOption("hist bar0");
3023  TH1* h = static_cast<TH1*>(ptr->GetObject());
3024  h->SetMarkerSize(2);
3025  ptr = ptr->Next();
3026  }
3027  rhs->SetMaximum(drawLog ? 2 : 1.1);
3028  rhs->SetMinimum(drawLog ? 1e-8 : 0);
3029  // Printf("Draw %s/%s in %d", sub.Data(), parT.Data(), padNo);
3030  DrawInPad(fBody, padNo, hp = rhs, Form("nostack %s",
3031  drawLog ? "logy" : ""));
3032  hp->GetHistogram()->SetYTitle(parT);
3033  }
3034  }
3035  // --- Make a legend -----------------------------------------------
3036  l = new TLegend(fBody->GetLeftMargin(),
3037  1-fBody->GetTopMargin(),
3038  .99, .99);
3039  l->SetNColumns(2);
3040  l->SetFillStyle(0);
3041  l->SetBorderSize(0);
3042  e = l->AddEntry("dummy", "Signal #minus #Delta<1.5", "f");
3043  e->SetFillColor(kGreen+1);
3044  e->SetFillStyle(1001);
3045  e = l->AddEntry("dummy", "Background #minus 5<#Delta<25", "f");
3046  e->SetFillColor(kRed+1);
3047  e->SetFillStyle(1001);
3048  fBody->cd();
3049  l->Draw();
3050 
3051  PrintCanvas(Form("Species contribution #topbar %s", fLastBin.Data()),
3052  Form("%s_species_contrib", fLastShort.Data()));
3053 
3054 
3055  // --- Draw selected particles -------------------------------------
3056  TH1* lt1p = GetH1(GetC(prim, "mid"), "totalIntegrals");
3057  TH1* lt1s = GetH1(GetC(seco, "mid"), "totalIntegrals");
3058  TH1* lt1c = GetH1(GetC(comb, "mid"), "totalIntegrals");
3059  TH1* gt1p = GetH1(GetC(prim, "fwd"), "totalIntegrals");
3060  TH1* gt1s = GetH1(GetC(seco, "fwd"), "totalIntegrals");
3061  TH1* gt1c = GetH1(GetC(comb, "fwd"), "totalIntegrals");
3062  Double_t totalA = (lt1p->GetBinContent(1)+
3063  lt1s->GetBinContent(1)+
3064  lt1c->GetBinContent(1)+
3065  gt1p->GetBinContent(1)+
3066  gt1s->GetBinContent(1)+
3067  gt1c->GetBinContent(1));
3068 
3069  // --- Draw relative contributions to signal/background ------------
3070  ClearCanvas();
3071  fBody->SetTopMargin(0.2);
3072  fBody->SetRightMargin(0.01);
3073  fBody->Divide(2,3,0,0);
3074  gStyle->SetPaintTextFormat("6.2f");
3075  drawLog = false;
3076  for (Int_t i = 0; i < 2; i++) {
3077  TString sub = (i == 0 ? "mid" : "fwd");
3078  TString bin = (i == 0 ? "|#eta|<1" : "|#eta|>1");
3079  for (Int_t j = 0; j < 3; j++) {
3080  Container* par = (j == 0 ? prim :
3081  j == 1 ? seco : comb);
3082  TString parT = (j == 0 ? "Primaries" :
3083  j == 1 ? "Secondaries" : "Combinatorics");
3084  Int_t padNo = j * 2 + i + 1;
3085  TVirtualPad* pad = fBody->cd(padNo);
3086  if ((padNo % 2) == 0) pad->SetRightMargin(0.10);
3087  else pad->SetLeftMargin(0.20);
3088  pad->SetTicks();
3089  pad->SetGridx();
3090  pad->SetGridy();
3091  TH1* itg = GetH1(GetC(par, sub), "totalIntegrals");
3092  THStack* rhs = GetHS(GetC(par, sub), "rows");
3093  TObjLink* ptr = rhs->GetHists()->FirstLink();
3094  while (ptr) {
3095  ptr->SetOption("hist bar0 text30");
3096  TH1* h = static_cast<TH1*>(ptr->GetObject());
3097  h->SetMarkerSize(3);
3098  ptr = ptr->Next();
3099  }
3100  rhs->SetMaximum(drawLog ? 200 : 110);
3101  rhs->SetMinimum(drawLog ? 1e-3 : 0);
3102  // Printf("Draw %s/%s in %d", sub.Data(), parT.Data(), padNo);
3103  DrawInPad(fBody, padNo, rhs, Form("nostack %s",
3104  drawLog ? "logy" : ""));
3105  rhs->GetHistogram()->SetYTitle(parT);
3106  rhs->GetHistogram()->GetXaxis()->SetLabelSize(0.03/pad->GetHNDC());
3107  rhs->GetHistogram()->GetYaxis()->SetLabelSize(0.03/pad->GetHNDC());
3108  rhs->GetHistogram()->GetYaxis()->SetTitleSize(0.03/pad->GetHNDC());
3109  rhs->GetHistogram()->GetYaxis()->SetNdivisions(207);
3110  rhs->GetHistogram()->GetYaxis()->SetTitleOffset(1);
3111 
3112  pad->cd();
3113  Double_t total = itg->GetBinContent(1);
3114  TLatex* txt = new TLatex(pad->GetLeftMargin()+.05,
3115  .99, Form("%5.2f%% of all",
3116  100*total/totalA));
3117  txt->SetTextAlign(13);
3118  txt->SetNDC();
3119  txt->SetTextSize(0.06);
3120  txt->Draw();
3121 
3122  txt = new TLatex(pad->GetLeftMargin()+.05,
3123  .99-txt->GetTextSize(),
3124  Form("Signal %5.2f%% of all",
3125  100*itg->GetBinContent(2)/totalA));
3126  txt->SetTextAlign(13);
3127  txt->SetNDC();
3128  txt->SetTextSize(0.06);
3129  txt->SetTextColor(kGreen+1);
3130  txt->Draw();
3131 
3132  txt = new TLatex(pad->GetLeftMargin()+.05,
3133  .99-2*txt->GetTextSize(),
3134  Form("Background %5.2f%% of all",
3135  100*itg->GetBinContent(3)/totalA));
3136  txt->SetTextAlign(13);
3137  txt->SetNDC();
3138  txt->SetTextSize(0.06);
3139  txt->SetTextColor(kRed+1);
3140  txt->Draw();
3141  }
3142  }
3143  // --- Make a legend -----------------------------------------------
3144  l = new TLegend(fBody->GetLeftMargin(),
3145  1-fBody->GetTopMargin(),
3146  .99, .99);
3147  l->SetNColumns(2);
3148  l->SetFillStyle(0);
3149  l->SetBorderSize(0);
3150  e = l->AddEntry("dummy", "Signal #minus #Delta<1.5", "f");
3151  e->SetFillColor(kGreen+1);
3152  e->SetFillStyle(1001);
3153  e = l->AddEntry("dummy", "Background #minus 5<#Delta<25", "f");
3154  e->SetFillColor(kRed+1);
3155  e->SetFillStyle(1001);
3156  fBody->cd();
3157  l->Draw();
3158 
3159  PrintCanvas(Form("Species strangeness #topbar %s", fLastBin.Data()),
3160  Form("%s_species_strange", fLastShort.Data()));
3161 
3162 
3163 
3164 }
3165 
3166 
3167 //____________________________________________________________________
3169 {
3170  DebugGuard g(fProc&kDebug,1,"Visualizing primaries");
3171  if (!simCont) return true;
3172  ClearCanvas();
3173 
3174  THStack* all = GetHS(GetC(simCont,"mix"),"all");
3175  THStack* toPion = GetHS(GetC(simCont,"mix"),"toPion");
3176  THStack* toAll = GetHS(GetC(simCont,"mix"),"toAll");
3177  TH2* etaPt = GetH2(simCont, "etaPt");
3178  Int_t etaM = etaPt->GetXaxis()->FindBin(-.5);
3179  Int_t etaP = etaPt->GetXaxis()->FindBin(+.5);
3180  TH1* pt = etaPt->ProjectionY("pt", etaM, etaP);
3181  pt->GetYaxis()->SetTitle("d#it{N}/d#it{p}_{T}");
3182  pt->GetYaxis()->SetTitleSize(0.08);
3183  pt->GetYaxis()->SetLabelSize(0.08);
3184  pt->GetYaxis()->SetTitleOffset(0.6);
3185  pt->GetXaxis()->SetTitleSize(0.08);
3186  pt->GetXaxis()->SetLabelSize(0.08);
3187  pt->GetXaxis()->SetTitleOffset(0.6);
3188 
3189  if (false) {
3190  TIter next(toPion->GetHists());
3191  TH1* rat = 0;
3192  while ((rat = static_cast<TH1*>(next()))) {
3193  TGraphErrors* g =
3194  static_cast<TGraphErrors*>(rat->GetListOfFunctions()
3195  ->FindObject(Form("%s_2760",
3196  rat->GetName())));
3197  TF1* f = new TF1("fit", "pol0", -.5, +.5);
3198  rat->Fit(f, "Q0R+", "", -.5, +.5);
3199  Double_t re, r = RatioE(f->GetParameter(0), f->GetParError(0),
3200  g->GetY()[0], g->GetEY()[0], re);
3201  Printf("%10s: 2760: %6.4f +/- %6.4f Here: %6.4f +/- %6.4f "
3202  "Ratio: %6.4f +/- %6.4f",
3203  rat->GetName(), g->GetY()[0], g->GetEY()[0],
3204  f->GetParameter(0), f->GetParError(0), r, re);
3205 
3206  f->SetLineColor(rat->GetLineColor());
3207  f->SetLineStyle(7);
3208  }
3209  }
3210 
3211  TPad* p1 = new TPad("p1","p1",0,.3,1,1);
3212  p1->SetTopMargin(.01);
3213  p1->SetRightMargin(.01);
3214  fBody->cd();
3215  p1->Draw();
3216  p1->SetNumber(1);
3217  p1->Divide(1,3,0,0);
3218  p1->GetPad(1)->SetRightMargin(0.2);
3219  p1->GetPad(2)->SetRightMargin(0.2);
3220  p1->GetPad(3)->SetRightMargin(0.2);
3221  TLegend* l = DrawInPad(p1,1,all, "leg2 nostack logy");
3222  all->GetHistogram()->GetYaxis()->SetTitle("d#it{N}_{X}/d#eta");
3223  all->GetHistogram()->GetYaxis()->SetTitleSize(0.08);
3224  all->GetHistogram()->GetYaxis()->SetLabelSize(0.08);
3225  all->GetHistogram()->GetYaxis()->SetTitleOffset(0.6);
3226  ModLegend(p1->GetPad(1), l,
3227  1-p1->GetPad(1)->GetRightMargin(),
3228  p1->GetPad(1)->GetBottomMargin(),
3229  1-p1->GetPad(1)->GetTopMargin(),
3230  .99);
3231  l->SetBorderSize(0);
3232 
3233  l = DrawInPad(p1,2,toAll, "nostack leg2 logy");
3234  toAll->GetHistogram()->GetYaxis()->SetTitle("Ratio to all");
3235  toAll->GetHistogram()->GetYaxis()->SetTitleSize(0.08);
3236  toAll->GetHistogram()->GetYaxis()->SetLabelSize(0.08);
3237  toAll->GetHistogram()->GetYaxis()->SetTitleOffset(0.6);
3238  toAll->GetHistogram()->GetXaxis()->SetTitle("#eta");
3239  toAll->GetHistogram()->GetXaxis()->SetTitleSize(0.08);
3240  toAll->GetHistogram()->GetXaxis()->SetLabelSize(0.08);
3241  toAll->GetHistogram()->GetXaxis()->SetTitleOffset(0.6);
3242  ModLegend(p1->GetPad(2), l,
3243  1-p1->GetPad(2)->GetRightMargin(),
3244  p1->GetPad(2)->GetBottomMargin(),
3245  1-p1->GetPad(2)->GetTopMargin(),
3246  .99);
3247  l->SetBorderSize(0);
3248 
3249  l = DrawInPad(p1,3,toPion, "nostack leg");
3250  toPion->GetHistogram()->GetYaxis()->SetTitle("Ratio to #pi");
3251  toPion->GetHistogram()->GetYaxis()->SetTitleSize(0.08);
3252  toPion->GetHistogram()->GetYaxis()->SetLabelSize(0.08);
3253  toPion->GetHistogram()->GetYaxis()->SetTitleOffset(0.6);
3254  toPion->GetHistogram()->GetXaxis()->SetTitle("#eta");
3255  toPion->GetHistogram()->GetXaxis()->SetTitleSize(0.08);
3256  toPion->GetHistogram()->GetXaxis()->SetLabelSize(0.08);
3257  toPion->GetHistogram()->GetXaxis()->SetTitleOffset(0.6);
3258  ModLegend(p1->GetPad(3), l,
3259  1-p1->GetPad(3)->GetRightMargin(),
3260  p1->GetPad(3)->GetBottomMargin(),
3261  1-p1->GetPad(3)->GetTopMargin(),
3262  .99);
3263  l->SetBorderSize(0);
3264 
3265  p1->Modified();
3266  p1->Update();
3267  p1->cd();
3268 
3269  TPad* p2 = new TPad("p2","p2",0,0,1,.3);
3270  p2->SetTopMargin(.01);
3271  p2->SetRightMargin(.01);
3272  p2->SetBottomMargin(0.15);
3273  fBody->cd();
3274  p2->Draw();
3275  p2->SetNumber(2);
3276  DrawInPad(p2,0,pt, "logx logy");
3277  p2->Modified();
3278  p2->Update();
3279  p2->cd();
3280 
3281  fBody->Modified();
3282 
3283  PrintCanvas(Form("Primary species #topbar %s", fLastBin.Data()),
3284  Form("%s_primary_species", fLastShort.Data()));
3285 
3286  return true;
3287 }
3288 
3289 //____________________________________________________________________
3291  Int_t dimen)
3292 {
3293  DebugGuard g(fProc&kDebug,1,"Visualizing delta in %s/%dD",
3294  outTop->GetName(), dimen);
3295  TDirectory* outDir = outTop->GetDirectory(Form("delta%dd", dimen));
3296  if (!outDir) {
3297  Warning("VisualizeDelta", "Directory detla%dd not found in %s",
3298  dimen, outTop->GetName());
3299  return false;
3300  }
3301 
3302  ClearCanvas();
3303  fBody->cd();
3304  TPad* pq = new TPad("p1","p1",0, .3, 1, 1);
3305  pq->SetNumber(1);
3306  pq->Draw();
3307 
3308  TVirtualPad* q = fBody->cd(1);
3309  q->SetTopMargin(0.01);
3310  q->SetRightMargin(0.01);
3311  q->Divide(1,2,0,0);
3312  q->GetPad(1)->SetRightMargin(0.15);
3313  q->GetPad(2)->SetRightMargin(0.15);
3314  TVirtualPad* qq = q->GetPad(1);
3315  THStack* all = GetHS(outDir,"all");
3316  TLegend* l = DrawInPad(q,1,all,"nostack logx logy grid leg");
3317  l->SetBorderSize(0);
3318  l->SetFillStyle(0);
3319  l->SetMargin(0.2);
3320  l->SetEntrySeparation(0.1);
3321  all->GetHistogram()->GetYaxis()->SetTitle("d#it{N}/d#Delta");
3322  all->GetHistogram()->GetYaxis()->SetLabelSize(0.06);
3323  all->GetHistogram()->GetYaxis()->SetTitleSize(0.06);
3324  all->GetHistogram()->GetYaxis()->SetTitleOffset(0.6);
3325  ModLegend(qq, l, 1-qq->GetRightMargin(), qq->GetBottomMargin(), .99,
3326  1-qq->GetTopMargin()-.01);
3327 
3328 
3329  THStack* ratios = GetHS(outDir,"ratios");
3330  ratios->SetMinimum(.6);
3331  ratios->SetMaximum(1.4);
3332  qq = q->GetPad(2);
3333  l = DrawInPad(q,2,ratios,"nostack logx grid leg");
3334  l->SetBorderSize(0);
3335  l->SetFillStyle(0);
3336  l->SetMargin(0.2);
3337  l->SetEntrySeparation(0.1);
3338  ratios->GetHistogram()->GetXaxis()->SetTitle("#Delta");
3339  ratios->GetHistogram()->GetYaxis()->SetTitle("Ratio");
3340  ratios->GetHistogram()->GetYaxis()->SetLabelSize(0.06);
3341  ratios->GetHistogram()->GetYaxis()->SetTitleSize(0.06);
3342  ratios->GetHistogram()->GetYaxis()->SetTitleOffset(0.6);
3343  ratios->GetHistogram()->GetXaxis()->SetLabelSize(0.06);
3344  ratios->GetHistogram()->GetXaxis()->SetTitleSize(0.06);
3345  ratios->GetHistogram()->GetXaxis()->SetTitleOffset(0.6);
3346 
3347  ModLegend(qq, l, 1-qq->GetRightMargin(), qq->GetBottomMargin(), .99,
3348  1-qq->GetTopMargin()-.01);
3349 
3350  fBody->cd();
3351  pq = new TPad("p2","p2",0, 0, 1, .3);
3352  pq->SetBottomMargin(0.25);
3353  pq->SetNumber(2);
3354  pq->Draw();
3355 
3356  q = fBody->cd(2);
3357  q->SetTopMargin(0.01);
3358  q->SetRightMargin(0.01);
3359  q->Divide(1,2,0,0);
3360  q->GetPad(1)->SetRightMargin(0.10);
3361  q->GetPad(2)->SetRightMargin(0.10);
3362  TH2* scale = GetH2(outDir,"scale");
3363  TH1* scaleProj = GetH1(outDir,"scaleProj");
3364  DrawInPad(q,1,scale, "colz");
3365  DrawInPad(q,2,scaleProj,"");
3366  scale->SetYTitle("IP_{#it{z}}");
3367  scaleProj->SetYTitle("#it{k}");
3368  scaleProj->SetXTitle("#eta");
3369  scale->GetYaxis()->SetLabelSize(0.12);
3370  scale->GetYaxis()->SetTitleSize(0.12);
3371  scale->GetYaxis()->SetTitleOffset(0.4);
3372  scaleProj->GetYaxis()->SetLabelSize(0.12);
3373  scaleProj->GetYaxis()->SetTitleSize(0.12);
3374  scaleProj->GetYaxis()->SetTitleOffset(0.4);
3375  scaleProj->GetYaxis()->SetNdivisions(207);
3376  scaleProj->GetXaxis()->SetLabelSize(0.12);
3377  scaleProj->GetXaxis()->SetTitleSize(0.12);
3378  scaleProj->GetXaxis()->SetTitleOffset(0.6);
3379 
3380  const char* what = (dimen == 3 ? "d^{3}N/(d#Deltad#etadIP_{z})" :
3381  dimen == 2 ? "d^{2}N/(d#Deltad#eta)" :
3382  dimen == 1 ? "dN/d#Delta" : "dN/d#Delta (k#equiv1)");
3383  PrintCanvas(Form("%s #topbar %s", what, fLastBin.Data()),
3384  Form("%s_deltas", fLastShort.Data()));
3385 
3386  return true;
3387 }
3388 
3389 //____________________________________________________________________
3391  Int_t dimen)
3392 {
3393  DebugGuard g(fProc&kDebug,1,"Visualize details for %s %dD",
3394  outTop->GetName(), dimen);
3395  TDirectory* outDir = outTop->GetDirectory(Form("results%dd", dimen));
3396  if (!outDir) {
3397  Warning("VisualizeDelta", "Directory results%dd not found in %s",
3398  dimen, outTop->GetName());
3399  return false;
3400  }
3401  Double_t tbase = 0.015;
3402  ClearCanvas();
3403  fBody->cd();
3404  fBody->SetTopMargin(0.01);
3405  fBody->SetRightMargin(0.01);
3406  fBody->SetBottomMargin(0.15);
3407  fBody->SetLeftMargin(0.15);
3408  fBody->Divide(2,6,0,0);
3409 
3410  const char* names[] = { "realM", "simM",
3411  "realC", "simC",
3412  "realB", "simB",
3413  "realS", "simS",
3414  "simA", "fiducial",
3415  "result", "simG" };
3416  const char* titles[] = { "M", "M'",
3417  "C", "C'",
3418  "#beta", "#beta'",
3419  "S", "S'",
3420  "#alpha", "F'",
3421  "R", "G'" };
3422 
3423  for (Int_t i = 0; i < 12; i++) {
3424  const char* name = names[i];
3425  const char* title = titles[i];
3426  TVirtualPad* pad = fBody->GetPad(i+1);
3427  pad->SetLeftMargin(0.15);
3428  pad->SetRightMargin(0);
3429  pad->Divide(1,2,0,0);
3430  TVirtualPad* q2 = pad->GetPad(1); q2->SetRightMargin(0.15);
3431  TVirtualPad* q1 = pad->GetPad(2); q1->SetRightMargin(0.15);
3432 
3433  TH2* h2 = GetH2(outDir, Form("full/%s", name));
3434  TH1* h1 = GetH1(outDir, name);
3435  if (!h2 || !h1) {
3436  Warning("VisualizeDetails", "Didn't find full/%s (%p) or %s (%p)",
3437  name, name);
3438  continue;
3439  }
3440  h2 = static_cast<TH2*>(h2->Clone());
3441  h1 = static_cast<TH1*>(h1->Clone());
3442  h2->SetDirectory(0);
3443  h1->SetDirectory(0);
3444  h2->SetXTitle("#eta");
3445  h1->SetXTitle("#eta");
3446  h2->SetYTitle(title);
3447  h1->SetYTitle("");
3448  TAxis* axis[] = { h2->GetXaxis(), h2->GetYaxis(), h2->GetZaxis(),
3449  h1->GetXaxis(), h1->GetYaxis(), 0 };
3450  TVirtualPad* pads[] = { q2, q2, q2, q1, q1, 0 };
3451  TAxis** pa = axis;
3452  TVirtualPad** pq = pads;
3453  while ((*pa)) {
3454  (*pa)->SetTitleSize(tbase/pad->GetHNDC()/(*pq)->GetHNDC());
3455  (*pa)->SetLabelSize(tbase/pad->GetHNDC()/(*pq)->GetHNDC());
3456  (*pa)->SetTitleOffset(0.4);
3457  (*pa)->SetNdivisions(207);
3458  pa++;
3459  pq++;
3460  }
3461  DrawInPad(pad, 1, h2, "colz");
3462  DrawInPad(pad, 2, h1, "");
3463  pad->Modified();
3464  }
3465  fBody->Modified();
3466 
3467  PrintCanvas(Form("Details #topbar %s", fLastBin.Data()),
3468  Form("%s_details", fLastShort.Data()));
3469  return true;
3470 }
3471 
3472 //____________________________________________________________________
3474  Int_t dimen)
3475 {
3476  DebugGuard g(fProc&kDebug,1,"Visualize results for %s %dD",
3477  outTop->GetName(), dimen);
3478  TDirectory* outDir = outTop->GetDirectory(Form("results%dd", dimen));
3479  if (!outDir) {
3480  Warning("VisualizeDelta", "Directory results%dd not found in %s",
3481  dimen, outTop->GetName());
3482  return false;
3483  }
3484 
3485  ClearCanvas();
3486  Double_t tbase = 0.03;
3487  Double_t yr = .2;
3488  Double_t yf = .2;
3489  TPad* p1 = new TPad("p1","p1",0,yf,1,1);
3490  p1->SetTopMargin(yr);
3491  p1->SetRightMargin(0.01);
3492  p1->SetLeftMargin(0.15);
3493  p1->SetBottomMargin(0);
3494  p1->SetTicks();
3495  fBody->cd();
3496  p1->Draw();
3497  p1->SetNumber(1);
3498 
3499  TPad* p2 = new TPad("p2","p2",0,0,1,yf);
3500  p2->SetTopMargin(0.01);
3501  p2->SetRightMargin(.01);
3502  p2->SetLeftMargin(0.15);
3503  p2->SetBottomMargin(0.20);
3504  p2->SetTicks();
3505  fBody->cd();
3506  p2->Draw();
3507  p2->SetNumber(2);
3508 
3509  THStack* all = GetHS(outDir, "all");
3510  TLegend* l = DrawInPad(fBody,1,all, "nostack leg2");
3511  all->GetHistogram()->SetXTitle("#eta");
3512  all->GetHistogram()->SetYTitle(ObsTitle());
3513  all->GetHistogram()->GetYaxis()->SetTitleOffset(1.7);
3514  all->GetHistogram()->GetYaxis()->SetTitleSize(tbase/(1-yf));
3515  all->GetHistogram()->GetYaxis()->SetLabelSize(tbase/(1-yf));
3516  all->GetHistogram()->GetYaxis()->SetNdivisions(205);
3517  l->SetBorderSize(0);
3518  l->SetFillStyle(0);
3519  l->SetMargin(0.12);
3520  l->SetEntrySeparation(0.1);
3521  l->SetHeader("R=G'/(M'-C')#times(M-C) #kern[2]{ } [C=kC'/M'#timesM]");
3522  // l->SetTextSize(0.04);
3523  // l->SetTextAlign(12);
3524  ModLegend(p1, l, p1->GetLeftMargin()-.01, 1-yr,
3525  1-p1->GetRightMargin(), .99);
3526 
3527  THStack* ratios = GetHS(outDir, "ratios");
3528  FixMinMax(ratios, true);
3529  DrawInPad(fBody,2,ratios, "nostack");
3530  ratios->GetHistogram()->SetXTitle("#eta");
3531  ratios->GetHistogram()->SetYTitle("Ratios");
3532  ratios->GetHistogram()->GetYaxis()->SetTitleOffset(.45);
3533  ratios->GetHistogram()->GetXaxis()->SetTitleOffset(.7);
3534  ratios->GetHistogram()->GetYaxis()->SetTitleSize(tbase/yr);
3535  ratios->GetHistogram()->GetYaxis()->SetLabelSize(tbase/yr);
3536  ratios->GetHistogram()->GetXaxis()->SetTitleSize(tbase/yr);
3537  ratios->GetHistogram()->GetXaxis()->SetLabelSize(tbase/yr);
3538  ratios->GetHistogram()->GetYaxis()->SetNdivisions(205);
3539 
3540  fBody->Modified();
3541  // fBody->GetListOfPrimitives()->Print();
3542 
3543  const char* what = (dimen == 3 ? "k(#eta,IP_{z})" :
3544  dimen == 2 ? "k(#eta)" :
3545  dimen == 1 ? "k=const." : "k#equiv1");
3546  PrintCanvas(Form("%s [%s] #topbar %s", ObsTitle(), what, fLastBin.Data()),
3547  Form("%s_summary", fLastShort.Data()));
3548 
3549  return true;
3550 }
3551 
3552 //====================================================================
3554  Color_t color,
3555  Style_t marker,
3556  Double_t size,
3557  Style_t fill,
3558  Style_t line,
3559  Width_t width)
3560 {
3561  h->SetMarkerColor(color);
3562  h->SetMarkerStyle(marker);
3563  h->SetMarkerSize(size);
3564  h->SetFillColor(color);
3565  h->SetFillStyle(fill);
3566  h->SetLineColor(color);
3567  h->SetLineStyle(line);
3568  h->SetLineWidth(width);
3569  h->GetXaxis()->SetNdivisions(210);
3570  h->GetYaxis()->SetNdivisions(210);
3571 }
3572 
3573 //____________________________________________________________________
3575 {
3576  Double_t sum = 0;
3577  Double_t sumw = 0;
3578  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
3579  Double_t c = h->GetBinContent(i);
3580  Double_t e = h->GetBinError (i);
3581  if (c < 1e-6 || e/c > 0.1) continue;
3582  Double_t w = 1/e/e;
3583  sum += w * c;
3584  sumw += w;
3585  }
3586  e = TMath::Sqrt(1/sumw);
3587  return sum / sumw;
3588 }
3589 //____________________________________________________________________
3591 {
3592  Double_t sum = 0;
3593  Double_t sumw = 0;
3594  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
3595  for (Int_t j = 1; j <= h->GetNbinsY(); j++) {
3596  Double_t c = h->GetBinContent(i,j);
3597  Double_t e = h->GetBinError (i,j);
3598  if (c < 1e-6 || e/c > 0.1) continue;
3599  Double_t w = 1/e/e;
3600  sum += w * c;
3601  sumw += w;
3602  }
3603  }
3604  e = TMath::Sqrt(1/sumw);
3605  return sum / sumw;
3606 }
3607 
3608 
3609 
3610 
3611 #endif
3612 //____________________________________________________________________
3613 //
3614 // EOF
3615 //
Int_t color[]
print message on plot with ok/not ok
Int_t pdg
Bool_t Deltas1D(Container *realCont, Container *simCont, TDirectory *outParent)
static const TAxis & PdgAxis()
static void PdgAttr(Int_t pdg, TString &nme, Color_t &c, Style_t &s)
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
const char * title
Definition: MakeQAPdf.C:27
void Run(UInt_t proc=kDefaultProc, UInt_t viz=kDefaultViz, UShort_t maxBins=9, const char *dataName="data.root", const char *simName="sim.root", const char *output=0, Double_t fudge=1)
void VisualizeParam(const char *name, Double_t &y, const char *val)
Definition: External.C:244
static void PrintAxis(const TAxis &axis, Int_t nSig=2, const char *alt=0)
void VisualizeGeneral(Container *realList, Container *simList)
void PrintCanvas(const char *title, const char *shortTitle="page", Float_t size=.7)
Bool_t VisualizeDetails(TDirectory *outTop, Int_t dimen)
Bool_t mix
Definition: InvMassFit.C:51
Bool_t VisualizeDelta(TDirectory *outTop, Int_t dimen)
static TFile * OpenFile(const char *filename)
void FlattenCentrality(TH2 *h, Container *c)
TSystem * gSystem
Bool_t Visualize(Container *realSums, Container *simSums, Container *realRess, Container *simRess, TDirectory *outTop, Int_t maxBins)
static TH1 * GetH1(Container *parent, const char *name, Bool_t verb=true)
static const char * CentName(Double_t c1, Double_t c2)
TH1 * FractionFit(TDirectory *outDir, TH1 *realDeltaM, TH1 *simDeltaC, TH1 *simDeltaP, TH1 *simDeltaS)
static TH3 * ScaleDelta(TH3 *h, TH2 *scale)
THStack * Make2Stack(const char *name, const char *title, Container *realList, Container *simList, Option_t *dataOpt="", Option_t *simOpt="")
TCanvas * c
Definition: TestFitELoss.C:172
Bool_t Deltas2D(Container *realCont, Container *simCont, TDirectory *outParent)
TLegend * DrawInPad(TVirtualPad *c, Int_t pad, TObject *o, Option_t *opt)
void VisualizeParams(Container *pars, const char *title)
AliStack * stack
static TH1 * ProjectDelta(TH2 *h)
Bool_t VisualizePrimary(Container *simCont)
static Double_t GetD(Container *parent, const char *name, Double_t def=-1, Bool_t verb=true)
Bool_t Deltas(Container *realCont, Container *simCont, TDirectory *outParent, Int_t dim)
Bool_t VisualizeBin(Double_t c1, Double_t c2, Container *simList, TDirectory *outTop)
static Bool_t CheckConsistency(const TH1 *h1, const TH1 *h2)
void WriteDeltas(TDirectory *outDir, TH1 *realDeltaM, TH1 *realDeltaI, TH1 *simDeltaM, TH1 *simDeltaI, TH1 *simDeltaC, TH1 *simDeltaP, TH1 *simDeltaS, TH1 *simDeltaD, TH1 *fit)
static TH1 * Scale(TH1 *h, Double_t x, Double_t xe)
void VisualizeFinal(TDirectory *outDir, Int_t i)
Utilities for midrapidity analysis.
void ModLegend(TVirtualPad *p, TLegend *l, Double_t x1, Double_t y1, Double_t x2, Double_t y2)
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)
Double_t MeanZ(TH2 *h, Double_t &e)
static TH1 * ProjectDeltaFull(TH3 *h)
unsigned int UInt_t
Definition: External.C:33
Bool_t VisualizeResult(TDirectory *outTop, Int_t dimen)
float Float_t
Definition: External.C:68
const char * ObsTitle() const
Double_t MeanY(TH1 *h, Double_t &e)
Definition: External.C:212
TH1 * Results(Container *realCont, Container *simCont, Double_t realEff, Double_t simEff, TDirectory *outParent, Int_t deltaDimen)
Bool_t Process(Container *realTop, Container *simTop, TDirectory *outTop, Int_t maxBins)
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)
static TH1 * CopyH1(Container *parent, const char *name, const char *newName=0, Bool_t verb=true)
static TH1 * AverageOverIPz(TH2 *h, const char *name, UShort_t mode, TH1 *ipz, Double_t ipEff=0, TH2 *mask=0, Bool_t verb=false)
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)
Bool_t ProcessBin(Double_t c1, Double_t c2, Container *realTop, Container *simTop, Double_t realEff, Double_t simEff, TDirectory *outTop)
Bool_t CalculateSEF(Container *simCont, const TAxis *centAxis, TDirectory *out)
static THStack * GetHS(Container *parent, const char *name, Bool_t verb=true)
void CreateCanvas(const TString &outputName)
Bool_t VisualizeSpeciesDelta(Container *simCont, TDirectory *outDir)
Definition: External.C:220
static Container * GetC(Container *parent, const char *name, Bool_t verb=true)
static void StackMinMax(THStack *stack, Double_t &min, Double_t &max)
unsigned short UShort_t
Definition: External.C:28
TList * ef
Definition: TestFitELoss.C:145
const char Option_t
Definition: External.C:48
AliTrackletAODUtils::Container Container
void VisualizeWeights(Container *simList)
static TH3 * CopyH3(Container *parent, const char *name, const char *newName=0, Bool_t verb=true)
bool Bool_t
Definition: External.C:53
Bool_t Deltas3D(Container *realCont, Container *simCont, TDirectory *outParent)
void VisualizeClosure(TDirectory *outDir, Int_t i)
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)
Bool_t Deltas0D(Container *realCont, Container *simCont, TDirectory *outParent)
TH1 * SetAttr(TH1 *h, Color_t color, Style_t marker=20, Double_t size=1., Style_t fill=0, Style_t line=1, Width_t width=1)
Definition: External.C:196
Bool_t VisualizeSpecies(Container *simCont)