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