AliPhysics  026afea (026afea)
 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 = 0x0010,
75  kDeltas = 0x0020,
77  kBackgrounds = 0x0040,
79  kAlphas = 0x0080,
81  kPDF = 0x0100,
82  kPNG = 0x0100,
84  kPause = 0x0200,
86  kLandscape = 0x0400,
88  kAltMarker = 0x0800,
90  kDefaultViz = 0x07ff
91  };
96  enum {
98  kUnitScale = 0x0001,
100  kAverageScale = 0x0002,
102  kEtaScale = 0x0004,
104  kFullScale = 0x0008,
106  kClosure = 0x01000,
108  kDefaultProc = 0x00083
109  };
110  enum {
111  kTopBackground = kAzure-8
112  };
113  //==================================================================
136  //==================================================================
138  TCanvas* fCanvas;
140  TPad* fTop;
142  TPad* fBody;
150  TLatex* fHeader;
152 
157 
158 
170  void Run(UInt_t proc = kDefaultProc,
171  UInt_t viz = kDefaultViz,
172  UShort_t maxBins = 9,
173  const char* dataName = "data.root",
174  const char* simName = "sim.root",
175  const char* output = 0,
176  Double_t fudge = 1);
187  Bool_t Process(Container* realTop,
188  Container* simTop,
189  TDirectory* outTop,
190  Int_t maxBins);
203  Double_t c2,
204  Container* realTop,
205  Container* simTop,
206  TDirectory* outTop);
221  Double_t c2,
222  Container* realCont,
223  Container* simCont,
224  TDirectory* outTop,
225  TDirectory* outDir,
226  Int_t dimen);
227 //__________________________________________________________________
242  Bool_t Deltas(Container* realCont,
243  Container* simCont,
244  TDirectory* outParent,
245  Int_t dim);
265  Bool_t Deltas0D(Container* realCont,
266  Container* simCont,
267  TDirectory* outParent);
286  Bool_t Deltas1D(Container* realCont,
287  Container* simCont,
288  TDirectory* outParent);
307  Bool_t Deltas2D(Container* realCont,
308  Container* simCont,
309  TDirectory* outParent);
330  Bool_t Deltas3D(Container* realCont,
331  Container* simCont,
332  TDirectory* outParent);
347  void WriteDeltas(TDirectory* outDir,
348  TH1* realDeltaM, TH1* realDeltaI,
349  TH1* simDeltaM, TH1* simDeltaI,
350  TH1* simDeltaC, TH1* simDeltaP,
351  TH1* simDeltaS, TH1* simDeltaD,
352  TH1* fit);
372  TH1* FractionFit(TDirectory* outDir,
373  TH1* realDeltaM,
374  TH1* simDeltaC,
375  TH1* simDeltaP,
376  TH1* simDeltaS);
377  /* @} */
378 
410  TH1* Results(Container* realCont,
411  Container* simCont,
412  TDirectory* outParent,
413  Int_t deltaDimen);
414 
415  /* @} */
416  //____________________________________________________________________
425  void ClearCanvas();
431  void CreateCanvas(const TString& outputName);
436  void CloseCanvas();
444  void PrintCanvas(const char* title,
445  const char* shortTitle="page",
446  Float_t size=.7);
457  TLegend* DrawInPad(TVirtualPad* c, Int_t pad, TObject* o, Option_t* opt);
468  void ModLegend(TVirtualPad* p, TLegend* l,
469  Double_t x1, Double_t y1,
470  Double_t x2, Double_t y2);
483  THStack* Make2Stack(const char* name,
484  const char* title,
485  Container* realList,
486  Container* simList,
487  Option_t* dataOpt="",
488  Option_t* simOpt="");
489  /* @} */
490  //____________________________________________________________________
507  Bool_t Visualize(Container* realSums,
508  Container* simSums,
509  Container* realRess,
510  Container* simRess,
511  TDirectory* outTop,
512  Int_t maxBins);
519  void VisualizeGeneral(Container* realList, Container* simList);
525  void VisualizeWeights(Container* simList);
532  void VisualizeFinal(TDirectory* outDir, Int_t i);
539  void VisualizeClosure(TDirectory* outDir, Int_t i);
551  Double_t c2,
552  Container* simList,
553  TDirectory* outTop);
578  Bool_t VisualizeDelta(TDirectory* outTop, Int_t dimen);
587  Bool_t VisualizeDetails(TDirectory* outTop, Int_t dimen);
596  Bool_t VisualizeResult(TDirectory* outTop, Int_t dimen);
597  /* @} */
598  //____________________________________________________________________
610  void VisualizeParam(const char* name, Double_t& y, const char* val);
618  void VisualizeParam(const char* name, Double_t& y, Double_t val);
626  void VisualizeParam(const char* name, Double_t& y, Int_t val);
634  void VisualizeParam(const char* name, Double_t& y, Bool_t val);
641  void VisualizeParams(Container* pars, const char* title);
648  void VisualizeParams(Container* realSums, Container* simSums);
649  /* @} */
650 
651  //____________________________________________________________________
663  TFile* OpenFile(const char* filename);
677  TH1* SetAttr(TH1* h,
678  Color_t color,
679  Style_t marker=20,
680  Double_t size=1.,
681  Style_t fill=0,
682  Style_t line=1,
683  Width_t width=1);
689  const char* ObsTitle() const { return "d#it{N}_{ch}/d#eta"; }
698  const char* CentName(Double_t c1, Double_t c2);
707  Double_t MeanY(TH1* h, Double_t& e);
716  Double_t MeanZ(TH2* h, Double_t& e);
717  /* @} */
718 };
719 
720 //====================================================================
721 namespace {
725  struct SuppressGuard2
726  {
728  Int_t save = 0;
734  SuppressGuard2(Int_t lvl=2000)
735  {
736  save = gErrorIgnoreLevel;
737  gErrorIgnoreLevel = lvl;
738  }
742  ~SuppressGuard2()
743  {
744  gErrorIgnoreLevel = save;
745  }
746  };
747 }
748 
749 //====================================================================
752  fProc(0),
753  fViz(0),
754  fDeltaCut(0),
755  fTailDelta(0),
756  fTailMax(0),
757  fMaxDelta(0),
758  fMinK(.7),
759  fMaxK(2),
760  fMinAlpha(0),
761  fMaxAlpha(2.5),
762  fFudge(0),
763  fCanvas(0),
764  fTop(0),
765  fBody(0),
766  fRealIsSim(false)
767 {}
768 
769 //====================================================================
771  UInt_t viz,
772  UShort_t maxBins,
773  const char* dataName,
774  const char* simName,
775  const char* outName,
776  Double_t fudge)
777 {
778  // Store options
779  fProc = proc;
780  fViz = viz;
781  fFudge = fudge;
782 
783  Printf("***************************************************\n"
784  " Data file: %s\n"
785  " Simulation file: %s\n"
786  " Output file: %s\n",
787  dataName, simName, outName);
788 
789 
790  // Open the input files
791  TFile* dataFile = 0;
792  TFile* simFile = 0;
793  if (!(dataFile = OpenFile(dataName))) return;
794  if (!(simFile = OpenFile(simName))) return;
795 
796  // Get some top-level contianers
797  fRealIsSim = false;
798  const char* base = "MidRapidity";
799  Container* realSums = GetC(dataFile, Form("%sSums", base));
800  Container* realRess = GetC(dataFile, Form("%sResults", base));
801  Container* simSums = GetC(simFile, Form("%sMCSums", base));
802  Container* simRess = GetC(simFile, Form("%sMCResults", base));
803  if (!realSums || !realRess) {
804  realSums = GetC(dataFile, Form("%sMCSums", base));
805  realRess = GetC(dataFile, Form("%sMCResults", base));
806  if (realSums && realRess)
807  Warning("Run","\n"
808  "*********************************************\n"
809  "* Testing MC vs. MC: *\n"
810  "* 'Data' file: %23s *\n"
811  "* Simulation file: %23s *\n"
812  "*********************************************\n",
813  dataName, simName);
814  fRealIsSim = true;
815  }
816  if (!realSums || !realRess || !simSums || !simRess) return;
817 
818  // Get parameters from the real data file
819  Container* params = GetC(realSums, "parameters");
820  fDeltaCut = GetD(params, "DeltaCut");
821  fTailDelta = GetD(params, "TailDelta");
822  fTailMax = GetD(params, "TailMax");
823  fMaxDelta = GetD(params, "MaxDelta");
824 
825  // Create output file name
826  TString outBase(outName);
827  if (outBase.IsNull()) outBase.Form("MiddNdeta_0x%04x", fProc);
828  if (outBase.EndsWith(".root")) outBase.ReplaceAll(".root", "");
829 
830  // Make output directory
831  gSystem->mkdir(outBase);
832 
833  // Open the output file
834  TFile* out = TFile::Open(Form("%s/result.root", outBase.Data()), "RECREATE");
835 
836 
837  if (!Process(realRess, simRess, out, maxBins)) return;
838 
839  out->Write();
840 
841  Visualize(realSums, simSums, realRess, simRess, out, maxBins);
842 }
843 
844 //====================================================================
846  Container* simTop,
847  TDirectory* outDir,
848  Int_t maxBins)
849 {
850  TH1* realCent = CopyH1(realTop, "cent", "realCent");
851  TH1* simCent = CopyH1(simTop, "cent", "simCent");
852  TH1* realIPz = GetH1(realTop, "ipz");
853  TH1* simIPz = GetH1(simTop, "ipz");
854 
855  // Check consistency of found histograms
856  if (!CheckConsistency(realCent, simCent)) {
857  Warning("Post", "Centrality bins are incompatible, giving up");
858  return false;
859  }
860  if (!CheckConsistency(realIPz, simIPz)) {
861  Warning("Post", "IPz bins are incompatible, giving up");
862  return false;
863  }
864 
865  // Check if we're doing a closure test
866  if (fProc & kClosure) realTop = simTop;
867 
868 
869  THStack* mids = new THStack("mids", "");
870  TH1* publ = 0;
871  {
872  Int_t nbin = 9;
873  Double_t bins[] = { 0., 5., 10., 20., 30., 40., 50., 60., 70., 80., };
874  Double_t vals[] = { 1948, 1590, 1180, 786, 512, 318, 183, 96.3, 44.9, };
875  Double_t errs[] = { 38, 32, 31, 20, 15, 12, 8, 5.8, 3.4, };
876  publ = new TH1D("published", "PRL116,222302 (2016)",
877  nbin, bins);
878  publ->SetMarkerStyle(24);
879  publ->SetMarkerColor(kBlack);
880  publ->SetMarkerSize(1.3);
881  for (Int_t i = 0; i < nbin; i++) {
882  publ->SetBinContent(i+1, vals[i]);
883  publ->SetBinError (i+1, errs[i]);
884  }
885  mids->Add(publ);
886  }
887 
888  // Make histogram for mid-rapidiy results
889  for (Int_t d = 0; d < 4; d++) {
890  if ((fProc & (1 << d)) == 0) continue;
891  TDirectory* dd = outDir->mkdir(Form("final%dd", d));
892  TH1* mid = Make1D(0,"mid",Form("%s|_{|#eta|<0.5}", ObsTitle()),
893  2+d, 20, *(realCent->GetXaxis()));
894  mid->SetDirectory(dd);
895  mid->SetXTitle("Centrality [%]");
896  mid->SetYTitle(mid->GetTitle());
897  mids->Add(mid);
898  THStack* full = new THStack("full","");
899  dd->cd();
900  full->Write();
901 
902  if (fRealIsSim) {
903  THStack* clos = new THStack("closure","");
904  clos->Write();
905  }
906  outDir->cd();
907  }
908 
909  // Write centralities to file
910  realCent->Write();
911  simCent ->Write();
912 
913  // Loop over defined centrality bins
914  for (Int_t i = 1; i <= realCent->GetNbinsX() && i <= maxBins ; i++) {
915  Double_t c1 = realCent->GetXaxis()->GetBinLowEdge(i);
916  Double_t c2 = realCent->GetXaxis()->GetBinUpEdge (i);
917 
918  ProcessBin(c1, c2, realTop, simTop, outDir);
919  }
920  // Close output file
921  outDir->cd();
922  mids->Write();
923 
924  TH1* mid = static_cast<TH1*>(mids->GetHists()->At(1));
925  for (Int_t i = 1; i <= mid->GetNbinsX(); i++) {
926  Double_t c1 = mid->GetXaxis()->GetBinLowEdge(i);
927  Double_t c2 = mid->GetXaxis()->GetBinUpEdge(i);
928  Int_t j = publ->GetXaxis()->FindBin((c1+c2)/2);
929  if (j < 1 || j > publ->GetNbinsX()) continue;
930  Double_t vh = mid ->GetBinContent(i);
931  Double_t eh = mid ->GetBinError (i);
932  Double_t vp = publ->GetBinContent(j);
933  Double_t ep = publ->GetBinError (j);
934  Double_t er, r = RatioE(vh,eh, vp, ep, er);
935  Printf("%5.1f - %5.1f%%: "
936  "Here %6.1f +/- %4.1f "
937  "Published %6.1f +/- %4.1f "
938  "Ratio %5.3f +/- %5.3f",
939  c1, c2, vh, eh, vp, ep, r, er);
940  }
941 
942  return true;
943 }
944 
945 
946 //____________________________________________________________________
948  Double_t c2,
949  Container* realTop,
950  Container* simTop,
951  TDirectory* outTop)
952 {
953  // Form the folder name
954  TString centName(CentName(c1,c2));
955 
956  // Get centrality containers
957  Container* realCont = GetC(realTop, centName);
958  Container* simCont = GetC(simTop, centName);
959  if (!realCont || !simCont) return false;
960 
961  TDirectory* outDir = outTop->mkdir(centName);
962 
963  printf("%5.1f - %5.1f%%:", c1, c2);
964  for (Int_t i = 0; i < 4; i++) {
965  if ((fProc & (1 << i)) == 0) continue;
966  if (!ProcessBin(c1, c2, realCont, simCont, outTop, outDir, i))
967  return false;
968  }
969  printf("\n");
970  return true;
971 }
972 //____________________________________________________________________
974  Double_t c2,
975  Container* realCont,
976  Container* simCont,
977  TDirectory* outTop,
978  TDirectory* outDir,
979  Int_t dimen)
980 {
981  if (!Deltas(realCont, simCont, outDir, dimen)) return false;
982  TH1* dndeta = Results(realCont, simCont, outDir, dimen);
983  if (!dndeta) {
984  Warning("ProcessBin", "Failed on Deltas for %f - %f", c1, c2);
985  return false;
986  }
987 
988  TDirectory* final = GetT(outTop,Form("final%dd", dimen));
989  if (!final) {
990  Warning("ProcessBin", "Failed on results for %f - %f", c1, c2);
991  return false;
992  }
993 
994  TH1* mid = static_cast<TH1*> (GetO(final, "mid"));
995  THStack* full = static_cast<THStack*>(GetO(final, "full"));
996  if (!mid || !final) {
997  Warning("ProcessBin", "Missing one of mid (%p) or full (%p)", mid, full);
998  return false;
999  }
1000 
1001  TF1* f = static_cast<TF1*>(dndeta->GetListOfFunctions()->At(0));
1002  if (!f) {
1003  Warning("ProcessBin", "No fit found on %s", dndeta->GetTitle());
1004  return false;
1005  }
1006 
1007  Double_t c = (c1+c2)/2;
1008  Int_t b = mid->GetXaxis()->FindBin(c);
1009  if (b < 1 || b > mid->GetNbinsX()) {
1010  Warning("ProcessBin", "Centrality %f - %f out of range", c1, c2);
1011  return false;
1012  }
1013 
1014  mid->SetBinContent(b, f->GetParameter(0));
1015  mid->SetBinError (b, f->GetParError (0));
1016 
1017  const Color_t cc[] = { kMagenta+2, // 0
1018  kBlue+2, // 1
1019  kAzure-1, // 2 // 10,
1020  kCyan+2, // 3
1021  kGreen+1, // 4
1022  kSpring+5, // 5 //+10,
1023  kYellow+1, // 6
1024  kOrange+5, // 7 //+10,
1025  kRed+1, // 8
1026  kPink+5, // 9 //+10,
1027  kBlack }; // 10
1028  Color_t tc = cc[b % 10];
1029  TH1* copy = static_cast<TH1*>(dndeta->Clone(outDir->GetName()));
1030  copy->SetDirectory(final);
1031  copy->GetListOfFunctions()->Clear();
1032  copy->SetTitle(Form("%5.1f#minus%5.1f%%", c1, c2));
1033  SetAttr(copy, tc);
1034  full->Add(copy);
1035  final->cd();
1036  full->Write(full->GetName(), TObject::kOverwrite);
1037 
1038  THStack* clos = static_cast<THStack*>(GetO(final, "closure"));
1039  TH1* clss = GetH1(GetT(outDir, Form("results%dd",dimen)), "closure");
1040  if (clos && clss) {
1041  copy = static_cast<TH1*>(clss->Clone(outDir->GetName()));
1042  copy->SetDirectory(0);
1043  copy->GetListOfFunctions()->Clear();
1044  copy->SetTitle(Form("%5.1f#minus%5.1f%%", c1, c2));
1045  SetAttr(copy, tc);
1046  clos->Add(copy);
1047  clos->Write(clos->GetName(), TObject::kOverwrite);
1048  }
1049 
1050 
1051  return true;
1052 }
1053 
1054 
1055 //====================================================================
1057  Container* simCont,
1058  TDirectory* outParent,
1059  Int_t dim)
1060 {
1061  switch (dim) {
1062  case 0: return Deltas0D(realCont, simCont, outParent);
1063  case 1: return Deltas1D(realCont, simCont, outParent);
1064  case 2: return Deltas2D(realCont, simCont, outParent);
1065  case 3: return Deltas3D(realCont, simCont, outParent);
1066  }
1067  return false;
1068 }
1069 //____________________________________________________________________
1071  Container* simCont,
1072  TDirectory* outParent)
1073 {
1074  // Make an output directory
1075  TDirectory* outDir = outParent->mkdir("delta0d");
1076 
1077  // Get the real and simulated measured folders
1078  Container* realMeas = GetC(realCont, "measured");
1079  Container* simMeas = GetC(simCont, "measured");
1080  if (!realMeas || !simMeas) return false;
1081 
1082  // Create a flat 2D scaling histogram for later user
1083  TH2* h = CopyH2(realMeas, "etaIPz", "scale");
1084  TH1* hp = h->ProjectionX("scaleProj");
1085  hp->Reset();
1086  hp->SetMinimum(fMinK);
1087  hp->SetMaximum(fMaxK);
1088  hp->SetYTitle("k");
1089  h->SetZTitle("k");
1090  h->SetTitle("k=1");
1091  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
1092  hp->SetBinContent(i, 1);
1093  hp->SetBinError (i, 0);
1094  for (Int_t j = 1; j <= h->GetNbinsY(); j++) {
1095  if (h->GetBinContent(i,j) < 1e-6) {
1096  h->SetBinContent(i,j,0);
1097  h->SetBinError (i,j,0);
1098  }
1099  else {
1100  h->SetBinContent(i,j,1);
1101  h->SetBinError (i,j,0);
1102  }
1103  }
1104  }
1105  h->SetDirectory(outDir);
1106  h->SetMinimum(fMinK);
1107  h->SetMaximum(fMaxK);
1108  hp->SetDirectory(outDir);
1109 
1110  // Get the raw projected Delta distributions for each component
1111  TH1* realDeltaM = CopyH1(realMeas, "delta","realDeltaM");
1112  TH1* realDeltaI = CopyH1(GetC(realCont,"injected"), "delta","realDeltaI");
1113  TH1* simDeltaM = CopyH1(simMeas, "delta","simDeltaM");
1114  TH1* simDeltaI = CopyH1(GetC(simCont, "injected"), "delta","simDeltaI");
1115  TH1* simDeltaC = CopyH1(GetC(simCont, "combinatorics"),"delta","simDeltaC");
1116  TH1* simDeltaP = CopyH1(GetC(simCont, "primaries"), "delta","simDeltaP");
1117  TH1* simDeltaS = CopyH1(GetC(simCont, "secondaries"), "delta","simDeltaS");
1118  TH1* simDeltaD = CopyH1(GetC(simCont, "distinct"), "delta","simDeltaD");
1119 
1120  // Get integrated scaling factor for injections, and scale the
1121  // injection distributions by that
1122  Double_t realScaleI = GetD(GetC(realCont,"injected"), "scale");
1123  Double_t realScaleIE = GetD(GetC(realCont,"injected"), "scaleError");
1124  Double_t simScaleI = GetD(GetC(simCont, "injected"), "scale");
1125  Double_t simScaleIE = GetD(GetC(simCont, "injected"), "scaleError");
1126  Scale(realDeltaI, realScaleI, realScaleIE);
1127  Scale(simDeltaI, simScaleI, simScaleIE);
1128  realDeltaI->SetTitle(Form("k_{I}#times%s",realScaleI,realDeltaI->GetTitle()));
1129  simDeltaI ->SetTitle(Form("k_{I'}#times%s",simScaleI,simDeltaI ->GetTitle()));
1130 
1131  TH1* fit = FractionFit(outDir, realDeltaM, simDeltaC, simDeltaP, simDeltaS);
1132  WriteDeltas(outDir, realDeltaM, realDeltaI, simDeltaM, simDeltaI,
1133  simDeltaC, simDeltaP, simDeltaS, simDeltaD, fit);
1134 
1135  outParent->cd();
1136  return true;
1137 }
1138 
1139 //____________________________________________________________________
1141  Container* simCont,
1142  TDirectory* outParent)
1143 {
1144  // Make an output directory
1145  TDirectory* outDir = outParent->mkdir("delta1d");
1146 
1147  // Get the real and simulated measured folders
1148  Container* realMeas = GetC(realCont, "measured");
1149  Container* simMeas = GetC(simCont, "measured");
1150  if (!realMeas || !simMeas) return false;
1151 
1152  // Get the integrated tails of the real and simulated observed
1153  // distribtutions, and calcualte scaling factor.
1154  Double_t realTail = GetD(realMeas, "deltaTail");
1155  Double_t realTailE = GetD(realMeas, "deltaTailError");
1156  Double_t simTail = GetD(simMeas, "deltaTail");
1157  Double_t simTailE = GetD(simMeas, "deltaTailError");
1158  Double_t scaleE = 0;
1159  Double_t scale = RatioE(realTail, realTailE, simTail, simTailE, scaleE);
1160 
1161  // Create a flat 2D scaling histogram for later user
1162  TH2* h = CopyH2(realMeas, "etaIPz", "scale");
1163  TH1* hp = h->ProjectionX("scaleProj");
1164  hp->Reset();
1165  hp->SetMinimum(fMinK);
1166  hp->SetMaximum(fMaxK);
1167  hp->SetYTitle("k");
1168  h->SetZTitle("k");
1169  h->SetTitle(Form("k=%5.3f#pm%5.3f", scale, scaleE));
1170  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
1171  hp->SetBinContent(i, scale);
1172  hp->SetBinError (i, scaleE);
1173  for (Int_t j = 1; j <= h->GetNbinsY(); j++) {
1174  if (h->GetBinContent(i,j) < 1e-6) {
1175  h->SetBinContent(i,j,0);
1176  h->SetBinError (i,j,0);
1177  }
1178  else {
1179  h->SetBinContent(i,j,scale);
1180  h->SetBinError (i,j,scaleE);
1181  }
1182  }
1183  }
1184  h->SetDirectory(outDir);
1185  h->SetMinimum(fMinK);
1186  h->SetMaximum(fMaxK);
1187  hp->SetDirectory(outDir);
1188 
1189  // Get the raw projected Delta distributions for each component
1190  TH1* realDeltaM = CopyH1(realMeas, "delta","realDeltaM");
1191  TH1* realDeltaI = CopyH1(GetC(realCont,"injected"), "delta","realDeltaI");
1192  TH1* simDeltaM = CopyH1(simMeas, "delta","simDeltaM");
1193  TH1* simDeltaI = CopyH1(GetC(simCont, "injected"), "delta","simDeltaI");
1194  TH1* simDeltaC = CopyH1(GetC(simCont, "combinatorics"),"delta","simDeltaC");
1195  TH1* simDeltaP = CopyH1(GetC(simCont, "primaries"), "delta","simDeltaP");
1196  TH1* simDeltaS = CopyH1(GetC(simCont, "secondaries"), "delta","simDeltaS");
1197  TH1* simDeltaD = CopyH1(GetC(simCont, "distinct"), "delta","simDeltaD");
1198 
1199  // Get integrated scaling factor for injections, and scale the
1200  // injection distributions by that
1201  Double_t realScaleI = GetD(GetC(realCont,"injected"), "scale");
1202  Double_t realScaleIE = GetD(GetC(realCont,"injected"), "scaleError");
1203  Double_t simScaleI = GetD(GetC(simCont, "injected"), "scale");
1204  Double_t simScaleIE = GetD(GetC(simCont, "injected"), "scaleError");
1205  Scale(realDeltaI, realScaleI, realScaleIE);
1206  Scale(simDeltaI, simScaleI, simScaleIE);
1207  realDeltaI->SetTitle(Form("k_{I}#times%s", realDeltaI->GetTitle()));
1208  simDeltaI ->SetTitle(Form("k_{I'}#times%s",simDeltaI ->GetTitle()));
1209 
1210  TH1* toScale[] = { simDeltaM,simDeltaI,simDeltaC,
1211  simDeltaP,simDeltaS,simDeltaD, 0};
1212  TH1** pScale = toScale;
1213  while ((*pScale)) {
1214  Scale(*pScale, scale, scaleE);
1215  (*pScale)->SetTitle(Form("k_{M}#times%s", (*pScale)->GetTitle()));
1216  pScale++;
1217  }
1218 
1219  TH1* fit = FractionFit(outDir, realDeltaM, simDeltaC, simDeltaP, simDeltaS);
1220  WriteDeltas(outDir, realDeltaM, realDeltaI, simDeltaM, simDeltaI,
1221  simDeltaC, simDeltaP, simDeltaS, simDeltaD, fit);
1222 
1223  outParent->cd();
1224  return true;
1225 }
1226 
1227 //____________________________________________________________________
1229  Container* simCont,
1230  TDirectory* outParent)
1231 {
1232  // Make an output directory
1233  TDirectory* outDir = outParent->mkdir("delta2d");
1234 
1235  // Get the real and simulated measured folders
1236  Container* realMeas = GetC(realCont, "measured");
1237  Container* simMeas = GetC(simCont, "measured");
1238  if (!realMeas || !simMeas) return false;
1239 
1240  // Get the eta-differential tails of the real and simulated observed
1241  // distribtutions, and calcualte scaling factor.
1242  TH1* realTail = GetH1(realMeas, "etaDeltaTail");
1243  TH1* simTail = GetH1(simMeas, "etaDeltaTail");
1244  TH1* scale = static_cast<TH1*>(realTail->Clone("scaleProj"));
1245  scale->Divide(simTail);
1246  scale->SetYTitle("k");
1247  Double_t sE, s = MeanY(scale, sE);
1248  TGraphErrors* g = new TGraphErrors(2);
1249  g->SetLineStyle(2);
1250  g->SetLineColor(kBlack);
1251  g->SetFillColor(kYellow);
1252  g->SetFillStyle(3002);
1253  g->SetPoint(0, scale->GetXaxis()->GetXmin(), s); g->SetPointError(0, 0, sE);
1254  g->SetPoint(1, scale->GetXaxis()->GetXmax(), s); g->SetPointError(1, 0, sE);
1255  scale->GetListOfFunctions()->Add(g, "le3");
1256 
1257  // Create a flat 2D scaling histogram for later user
1258  TH2* h = CopyH2(realMeas, "etaIPz", "scale");
1259  h->SetZTitle("k");
1260  h->SetTitle(Form("#LTk(#eta)#GT_{#eta}=%5.3f#pm%5.3f", s, sE));
1261  for (Int_t i = 1; i <= h->GetNbinsX(); i++) { // eta
1262  for (Int_t j = 1; j <= h->GetNbinsY(); j++) { // ipz
1263  if (h->GetBinContent(i,j) < 1e-6) {
1264  h->SetBinContent(i,j,0);
1265  h->SetBinError (i,j,0);
1266  }
1267  else {
1268  h->SetBinContent(i,j,scale->GetBinContent(i));
1269  h->SetBinError (i,j,scale->GetBinError (i));
1270  }
1271  }
1272  }
1273  h->SetMinimum(fMinK);
1274  h->SetMaximum(fMaxK);
1275  h->SetDirectory(outDir);
1276  scale->SetMinimum(fMinK);
1277  scale->SetMaximum(fMaxK);
1278  scale->SetDirectory(outDir);
1279 
1280  // Get the raw projected Delta distributions for each component
1281  TH2* r2DeltaM = CopyH2(realMeas, "etaDelta","r2DeltaM");
1282  TH2* r2DeltaI = CopyH2(GetC(realCont,"injected"), "etaDelta","r2DeltaI");
1283  TH2* s2DeltaM = CopyH2(simMeas, "etaDelta","s2DeltaM");
1284  TH2* s2DeltaI = CopyH2(GetC(simCont, "injected"), "etaDelta","s2DeltaI");
1285  TH2* s2DeltaC = CopyH2(GetC(simCont, "combinatorics"),"etaDelta","s2DeltaC");
1286  TH2* s2DeltaP = CopyH2(GetC(simCont, "primaries"), "etaDelta","s2DeltaP");
1287  TH2* s2DeltaS = CopyH2(GetC(simCont, "secondaries"), "etaDelta","s2DeltaS");
1288  TH2* s2DeltaD = CopyH2(GetC(simCont, "distinct"), "etaDelta","s2DeltaD");
1289 
1290  // Get eta-differential scaling factor for injections, and scale the
1291  // injection distributions by that
1292  TH1* rScaleI = GetH1(GetC(realCont,"injected"), "etaScale");
1293  TH1* sScaleI = GetH1(GetC(simCont, "injected"), "etaScale");
1294  Double_t rIE, rI = MeanY(rScaleI, rIE);
1295  Double_t sIE, sI = MeanY(sScaleI, sIE);
1296  Scale(r2DeltaI, rScaleI);
1297  Scale(s2DeltaI, sScaleI);
1298  r2DeltaI ->SetTitle(Form("#LTk_{I}#GT_{#eta}#times%s",
1299  r2DeltaI->GetTitle()));
1300  s2DeltaI ->SetTitle(Form("#LTk_{I'}#GT_{#eta}#times%s",
1301  s2DeltaI->GetTitle()));
1302 
1303  TH2* toScale[] = { s2DeltaM,s2DeltaI,s2DeltaC,s2DeltaP,s2DeltaS,s2DeltaD,0 };
1304  TH2** pScale = toScale;
1305  while ((*pScale)) {
1306  Scale(*pScale, scale);
1307  (*pScale)->SetTitle(Form("#LTk_{M}#GT_{#eta}#times%s",
1308  (*pScale)->GetTitle()));
1309  pScale++;
1310  }
1311 
1312  TH1* rDeltaM = ProjectDelta(r2DeltaM);
1313  TH1* rDeltaI = ProjectDelta(r2DeltaI);
1314  TH1* sDeltaM = ProjectDelta(s2DeltaM);
1315  TH1* sDeltaI = ProjectDelta(s2DeltaI);
1316  TH1* sDeltaC = ProjectDelta(s2DeltaC);
1317  TH1* sDeltaP = ProjectDelta(s2DeltaP);
1318  TH1* sDeltaS = ProjectDelta(s2DeltaS);
1319  TH1* sDeltaD = ProjectDelta(s2DeltaD);
1320  rDeltaM->SetTitle(r2DeltaM->GetTitle()); rDeltaM->SetName("realDeltaM");
1321  rDeltaI->SetTitle(r2DeltaI->GetTitle()); rDeltaI->SetName("realDeltaI");
1322  sDeltaM->SetTitle(s2DeltaM->GetTitle()); sDeltaM->SetName("simDeltaM");
1323  sDeltaI->SetTitle(s2DeltaI->GetTitle()); sDeltaI->SetName("simDeltaI");
1324  sDeltaC->SetTitle(s2DeltaC->GetTitle()); sDeltaC->SetName("simDeltaC");
1325  sDeltaP->SetTitle(s2DeltaP->GetTitle()); sDeltaP->SetName("simDeltaP");
1326  sDeltaS->SetTitle(s2DeltaS->GetTitle()); sDeltaS->SetName("simDeltaS");
1327  sDeltaD->SetTitle(s2DeltaD->GetTitle()); sDeltaD->SetName("simDeltaD");
1328 
1329  TH1* f2 = FractionFit(outDir, r2DeltaM, s2DeltaC, s2DeltaP, s2DeltaS);
1330  TH1* fit = ProjectDelta(static_cast<TH2*>(f2));
1331  WriteDeltas(outDir,rDeltaM,rDeltaI,
1332  sDeltaM,sDeltaI,sDeltaC,
1333  sDeltaP,sDeltaS,sDeltaD,
1334  fit);
1335 
1336  TDirectory* full = outDir->mkdir("full");
1337  r2DeltaM->SetDirectory(full); r2DeltaM->SetName("realDeltaM");
1338  r2DeltaI->SetDirectory(full); r2DeltaI->SetName("realDeltaI");
1339  s2DeltaM->SetDirectory(full); s2DeltaM->SetName("simDeltaM");
1340  s2DeltaI->SetDirectory(full); s2DeltaI->SetName("simDeltaI");
1341  s2DeltaC->SetDirectory(full); s2DeltaC->SetName("simDeltaC");
1342  s2DeltaP->SetDirectory(full); s2DeltaP->SetName("simDeltaP");
1343  s2DeltaS->SetDirectory(full); s2DeltaS->SetName("simDeltaS");
1344  s2DeltaD->SetDirectory(full); s2DeltaD->SetName("simDeltaD");
1345 
1346  outParent->cd();
1347  return true;
1348 }
1349 
1350 //____________________________________________________________________
1352  Container* simCont,
1353  TDirectory* outParent)
1354 {
1355  // Make an output directory
1356  TDirectory* outDir = outParent->mkdir("delta3d");
1357 
1358  // Get the real and simulated measured folders
1359  Container* realMeas = GetC(realCont, "measured");
1360  Container* simMeas = GetC(simCont, "measured");
1361  if (!realMeas || !simMeas) return false;
1362 
1363  // Get the eta-differential tails of the real and simulated observed
1364  // distribtutions, and calcualte scaling factor.
1365  TH2* realTail = GetH2(realMeas, "etaIPzDeltaTail");
1366  TH2* simTail = GetH2(simMeas, "etaIPzDeltaTail");
1367  TH2* scale = static_cast<TH2*>(realTail->Clone("scale"));
1368  scale->SetDirectory(0);
1369  scale->Divide(simTail);
1370  Double_t sE, s = MeanZ(scale, sE);
1371  scale->SetZTitle("k");
1372  scale->SetTitle(Form("#LTk(#eta)#GT_{#eta,IP_{#it{z}}}=%5.3f#pm%5.3f",
1373  s, sE));
1374  scale->SetDirectory(outDir);
1375  scale->SetMinimum(fMinK);
1376  scale->SetMaximum(fMaxK);
1377  TH1* etaScale = AverageOverIPz(scale, "scaleProj", 1, 0, 0);
1378  etaScale->SetYTitle("#LTk#GT_{IP_{z}}");
1379  etaScale->SetDirectory(outDir);
1380  etaScale->SetMinimum(fMinK);
1381  etaScale->SetMaximum(fMaxK);
1382  TGraphErrors* g = new TGraphErrors(2);
1383  g->SetLineStyle(2);
1384  g->SetLineColor(kBlack);
1385  g->SetFillColor(kYellow);
1386  g->SetFillStyle(3002);
1387  g->SetPoint(0, etaScale->GetXaxis()->GetXmin(), s);
1388  g->SetPoint(1, etaScale->GetXaxis()->GetXmax(), s);
1389  g->SetPointError(0, 0, sE);
1390  g->SetPointError(1, 0, sE);
1391  etaScale->GetListOfFunctions()->Add(g, "le3");
1392 
1393  // Get the raw projected Delta distributions for each component
1394  const char* nm = "etaDeltaIPz";
1395  TH3* r3DeltaM = CopyH3(realMeas, nm,"r3DeltaM");
1396  TH3* r3DeltaI = CopyH3(GetC(realCont,"injected"), nm,"r3DeltaI");
1397  TH3* s3DeltaM = CopyH3(simMeas, nm,"s3DeltaM");
1398  TH3* s3DeltaI = CopyH3(GetC(simCont, "injected"), nm,"s3DeltaI");
1399  TH3* s3DeltaC = CopyH3(GetC(simCont, "combinatorics"),nm,"s3DeltaC");
1400  TH3* s3DeltaP = CopyH3(GetC(simCont, "primaries"), nm,"s3DeltaP");
1401  TH3* s3DeltaS = CopyH3(GetC(simCont, "secondaries"), nm,"s3DeltaS");
1402  TH3* s3DeltaD = CopyH3(GetC(simCont, "distinct"), nm,"s3DeltaD");
1403 
1404  // Get eta-differential scaling factor for injections, and scale the
1405  // injection distributions by that
1406  TH2* rScaleI = GetH2(GetC(realCont,"injected"), "etaIPzScale");
1407  TH2* sScaleI = GetH2(GetC(simCont, "injected"), "etaIPzScale");
1408  Double_t rIE, rI = MeanZ(rScaleI, rIE);
1409  Double_t sIE, sI = MeanZ(sScaleI, sIE);
1410  ScaleDelta(r3DeltaI, rScaleI);
1411  ScaleDelta(s3DeltaI, sScaleI);
1412  r3DeltaI ->SetTitle(Form("#LTk_{I}#GT_{#eta,IP_{#it{z}}}#times%s",
1413  r3DeltaI->GetTitle()));
1414  s3DeltaI ->SetTitle(Form("#LTk_{I'}#GT_{#eta,IP_{#it{z}}}#times%s",
1415  s3DeltaI->GetTitle()));
1416 #if 0
1417  TH2* scale2 = static_cast<TH2*>(scale->Clone("scaleMain"));
1418  scale2->SetDirectory(0);
1419  scale2->Reset();
1420  Int_t sigBin = r3DeltaM->GetYaxis()->FindBin(1.5);
1421  for (Int_t i = 1; i <= r3DeltaM->GetNbinsX(); i++) {
1422  for (Int_t j = 1; j <= r3DeltaM->GetNbinsZ(); j++) {
1423  // Integrate over Delta
1424  Double_t rintg = 0, reintg = 0;
1425  Double_t sintg = 0, seintg = 0;
1426  rintg = r3DeltaM->IntegralAndError(i,i,1,sigBin,j,j,reintg);
1427  sintg = s3DeltaM->IntegralAndError(i,i,1,sigBin,j,j,seintg);
1428  Double_t re, r = RatioE(rintg, reintg, sintg, seintg, re);
1429 
1430  scale2->SetBinContent(i, j, r);
1431  scale2->SetBinError (i, j, re);
1432  }
1433  }
1434  Double_t rS2, rS = MeanZ(scale2, rS2);
1435  Printf("Scalar of Inj %6.3f +/- %6.3f", rS, rS2);
1436 #endif
1437 
1438  TH3* toScale[] = { s3DeltaM,s3DeltaI,s3DeltaC,s3DeltaP,s3DeltaS,s3DeltaD,0};
1439  TH3** pScale = toScale;
1440  while ((*pScale)) {
1441  ScaleDelta(*pScale, scale);
1442  (*pScale)->SetTitle(Form("#LTk_{M}#GT_{#eta,IP_{#it{z}}}#times%s",
1443  (*pScale)->GetTitle()));
1444  pScale++;
1445  }
1446 
1447  TH1* rDeltaM = ProjectDeltaFull(r3DeltaM);
1448  TH1* rDeltaI = ProjectDeltaFull(r3DeltaI);
1449  TH1* sDeltaM = ProjectDeltaFull(s3DeltaM);
1450  TH1* sDeltaI = ProjectDeltaFull(s3DeltaI);
1451  TH1* sDeltaC = ProjectDeltaFull(s3DeltaC);
1452  TH1* sDeltaP = ProjectDeltaFull(s3DeltaP);
1453  TH1* sDeltaS = ProjectDeltaFull(s3DeltaS);
1454  TH1* sDeltaD = ProjectDeltaFull(s3DeltaD);
1455  rDeltaM->SetTitle(r3DeltaM->GetTitle()); rDeltaM->SetName("realDeltaM");
1456  rDeltaI->SetTitle(r3DeltaI->GetTitle()); rDeltaI->SetName("realDeltaI");
1457  sDeltaM->SetTitle(s3DeltaM->GetTitle()); sDeltaM->SetName("simDeltaM");
1458  sDeltaI->SetTitle(s3DeltaI->GetTitle()); sDeltaI->SetName("simDeltaI");
1459  sDeltaC->SetTitle(s3DeltaC->GetTitle()); sDeltaC->SetName("simDeltaC");
1460  sDeltaP->SetTitle(s3DeltaP->GetTitle()); sDeltaP->SetName("simDeltaP");
1461  sDeltaS->SetTitle(s3DeltaS->GetTitle()); sDeltaS->SetName("simDeltaS");
1462  sDeltaD->SetTitle(s3DeltaD->GetTitle()); sDeltaD->SetName("simDeltaD");
1463 
1464  TH1* f3 = FractionFit(outDir, r3DeltaM, s3DeltaC, s3DeltaP, s3DeltaS);
1465  TH1* fit = ProjectDeltaFull(static_cast<TH3*>(f3));
1466  WriteDeltas(outDir,rDeltaM,rDeltaI,
1467  sDeltaM,sDeltaI,sDeltaC,
1468  sDeltaP,sDeltaS,sDeltaD,
1469  fit);
1470 
1471  TDirectory* full = outDir->mkdir("full");
1472  r3DeltaM->SetDirectory(full); r3DeltaM->SetName("realDeltaM");
1473  r3DeltaI->SetDirectory(full); r3DeltaI->SetName("realDeltaI");
1474  s3DeltaM->SetDirectory(full); s3DeltaM->SetName("simDeltaM");
1475  s3DeltaI->SetDirectory(full); s3DeltaI->SetName("simDeltaI");
1476  s3DeltaC->SetDirectory(full); s3DeltaC->SetName("simDeltaC");
1477  s3DeltaP->SetDirectory(full); s3DeltaP->SetName("simDeltaP");
1478  s3DeltaS->SetDirectory(full); s3DeltaS->SetName("simDeltaS");
1479  s3DeltaD->SetDirectory(full); s3DeltaD->SetName("simDeltaD");
1480 
1481  outParent->cd();
1482  return true;
1483 }
1484 
1485 //____________________________________________________________________
1487  TH1* rDeltaM,
1488  TH1* sDeltaC,
1489  TH1* sDeltaP,
1490  TH1* sDeltaS)
1491 {
1492  return 0;
1493  if (!rDeltaM || !sDeltaC || !sDeltaP || !sDeltaS) {
1494  Warning("FractionFit", "Missing M=%p, C'=%s, P'=%s, or S'=%p",
1495  rDeltaM, sDeltaC, sDeltaP, sDeltaS);
1496  return 0;
1497  }
1498  TH1* arr[] = { rDeltaM, sDeltaC, sDeltaP, sDeltaS, 0 };
1499  TH1** ptr = arr;
1500  while (*ptr) {
1501  TH1* tmp = static_cast<TH1*>((*ptr)->Clone());
1502  tmp->SetDirectory(0);
1503  Double_t intg = tmp->Integral();
1504  tmp->Scale(1./intg);
1505  *ptr = tmp;
1506  ptr++;
1507  }
1508  TObjArray mc;
1509  mc.SetOwner();
1510  mc.Add(arr[1]);
1511  mc.Add(arr[2]);
1512  mc.Add(arr[3]);
1513  TFractionFitter f(arr[0], &mc, "Q");
1514  Int_t status = f.Fit();
1515  if (status != 0) {
1516  Warning("FractionFit", "Fit failed w/status=%d", status);
1517  return 0;
1518  }
1519  Printf("\nTemplate fits");
1520  for (Int_t i = 0; i < 3; i++) {
1521  Double_t v, e;
1522  f.GetResult(i, v, e);
1523  Printf("%30s=%6.4f +/- %6.4f",
1524  mc.At(i)->GetName(), e, v);
1525  }
1526  TH1* ret = f.GetPlot();
1527 
1528  if (ret) {
1529  ret->Scale(rDeltaM->Integral());
1530  }
1531  return ret;
1532 }
1533 
1534 //____________________________________________________________________
1535 void AliTrackletdNdeta2::WriteDeltas(TDirectory* outDir,
1536  TH1* rDeltaM, TH1* rDeltaI,
1537  TH1* sDeltaM, TH1* sDeltaI,
1538  TH1* sDeltaC, TH1* sDeltaP,
1539  TH1* sDeltaS, TH1* sDeltaD,
1540  TH1* fit)
1541 {
1542  THStack* all = new THStack("all", "");
1543  SetAttr(rDeltaM, kRed+2, 20, 1.0);
1544  SetAttr(rDeltaI, kOrange+2, 21, 1.0);
1545  rDeltaM->SetDirectory(outDir);
1546  rDeltaI->SetDirectory(outDir);
1547  all->Add(rDeltaM);
1548  all->Add(rDeltaI);
1549 
1550  TH1* toScale[] = { sDeltaM,sDeltaI,sDeltaC, sDeltaP,sDeltaS,sDeltaD,fit,0};
1551  Color_t toColor[] = { kRed, kOrange,kMagenta,kGreen, kBlue, kPink, kBlack };
1552  Style_t toStyle[] = { 24, 25, 30, 26, 32, 30, 24 };
1553  TH1** pScale = toScale;
1554  Color_t* pColor = toColor;
1555  Style_t* pStyle = toStyle;
1556  while ((*pScale)) {
1557  (*pScale)->SetDirectory(outDir);
1558  all->Add((*pScale));
1559  SetAttr(*pScale, (*pColor)+2, *pStyle, 1.2);
1560  pScale++;
1561  pColor++;
1562  pStyle++;
1563  }
1564  outDir->cd();
1565  all->Write();
1566 
1567  THStack* ratios = new THStack("ratios", "");
1568  TH1* ratioM = static_cast<TH1*>(sDeltaM->Clone("ratioM"));
1569  ratioM->SetTitle("#Delta_{M'}/#Delta_{M}");
1570  ratioM->Divide(rDeltaM);
1571  ratioM->SetDirectory(outDir);
1572  ratios->Add(ratioM);
1573 
1574  TH1* ratioI = static_cast<TH1*>(sDeltaI->Clone("ratioI"));
1575  ratioI->SetTitle("#Delta_{I'}/#Delta_{I}");
1576  ratioI->Divide(rDeltaI);
1577  ratioI->SetDirectory(outDir);
1578  ratios->Add(ratioI);
1579 
1580  TH1* ratioIC = static_cast<TH1*>(sDeltaC->Clone("ratioIC"));
1581  ratioIC->SetTitle("#Delta_{C'}/#Delta_{I}");
1582  ratioIC->Divide(rDeltaI);
1583  ratioIC->SetDirectory(outDir);
1584  ratios->Add(ratioIC);
1585 
1586  TH1* ratioF = static_cast<TH1*>(sDeltaC->Clone("ratioF"));
1587  ratioF->SetTitle("#Delta_{fit}/#Delta_{M}");
1588  ratioF->Divide(rDeltaM);
1589  ratioF->SetDirectory(outDir);
1590  ratios->Add(ratioF);
1591 
1592  ratios->Write();
1593 }
1594 
1595 
1596 //====================================================================
1598  Container* simCont,
1599  TDirectory* outParent,
1600  Int_t deltaDimen)
1601 {
1602  TDirectory* outDir = outParent->mkdir(Form("results%dd", deltaDimen));
1603  TDirectory* delDir = outParent->GetDirectory(Form("delta%dd", deltaDimen));
1604 
1605  TH2* scale = static_cast<TH2*>(delDir->Get("scale"));
1606  TH2* realM = CopyH2(GetC(realCont, "measured"), "etaIPz", "realM");
1607  TH2* realS = CopyH2(GetC(realCont, "measured"), "etaIPz", "realS");
1608  TH2* realC = CopyH2(GetC(realCont, "measured"), "etaIPz", "realC");
1609  TH2* simM = CopyH2(GetC(simCont, "measured"), "etaIPz", "simM");
1610  TH2* simC = CopyH2(GetC(simCont, "combinatorics"), "etaIPz", "simC");
1611  TH2* simG = CopyH2(GetC(simCont, "generated"), "etaIPz", "simG");
1612  TH2* simS = CopyH2(GetC(simCont, "measured"), "etaIPz", "simS");
1613  TH2* simA = CopyH2(GetC(simCont, "generated"), "etaIPz", "simA");
1614  TH2* simB = CopyH2(GetC(simCont, "combinatorics"), "etaIPz", "simB");
1615  TH2* realG = (!fRealIsSim ? 0 :
1616  CopyH2(GetC(realCont,"generated"), "etaIPz", "realG"));
1617  TH1* realZ = CopyH1(realCont, "ipz", "realZ");
1618  TH1* simZ = CopyH1(simCont, "ipz", "simZ");
1619 
1620  // Scale combinatorial background to measured to get beta
1621  simB->Divide(simM);
1622  simB->SetTitle("#beta'");
1623  simB->SetZTitle("#beta'");
1624 
1625  // Copy simulated beta to real beta, and scale by scalar
1626  TH2* realB = static_cast<TH2*>(simB->Clone("realB"));
1627  realB->SetDirectory(0);
1628  realB->SetTitle("#beta");
1629  realB->SetZTitle("#beta");
1630  realB->Multiply(scale);
1631  Scale(realB, fFudge, 0);
1632 
1633  // Multiply real beta onto real measured to get background
1634  realC->Multiply(realB);
1635  realC->SetTitle("C");
1636  realC->SetZTitle("C");
1637 
1638  // Substract the real background off the real measured
1639  realS->Add(realC, -1);
1640  realS->SetTitle("S");
1641  realS->SetZTitle("S");
1642 
1643  // Substract combinatorial background from measured to get signal
1644  simS->Add(simC, -1);
1645  simS->SetTitle("S'");
1646  simS->SetZTitle("S'");
1647 
1648  // Scale MC truth primaries by signal to get correction
1649  simA->Divide(simS);
1650  simA->SetTitle("A'");
1651  simA->SetZTitle("#alpha'");
1652 
1653  // Make a fiducial distribution, and coerce the others to fit this
1654  TH2* fiducial = static_cast<TH2*>(simA->Clone("fiducial"));
1655  fiducial->SetTitle("F");
1656  fiducial->SetZTitle("F");
1657  for (Int_t i = 1; i <= fiducial->GetNbinsX(); i++) {
1658  for (Int_t j = 1; j <= fiducial->GetNbinsY(); j++) {
1659  Double_t c = fiducial->GetBinContent(i,j);
1660  fiducial->SetBinContent(i,j,c > fMinAlpha && c <= fMaxAlpha);
1661  fiducial->SetBinError(i,j,0);
1662  }
1663  }
1664 #if 1
1665  realM->Multiply(fiducial);
1666  realC->Multiply(fiducial);
1667  realS->Multiply(fiducial);
1668  realB->Multiply(fiducial);
1669  simM ->Multiply(fiducial);
1670  simC ->Multiply(fiducial);
1671  simS ->Multiply(fiducial);
1672  simB ->Multiply(fiducial);
1673  simA ->Multiply(fiducial);
1674 #endif
1675 
1676  // We can now make our result
1677  TH2* result = static_cast<TH2*>(realS->Clone("result"));
1678  result->Multiply(simA);
1679  result->SetTitle("R");
1680  result->SetZTitle("R");
1681 
1682  // Output directory for full stuff
1683  TDirectory* full = outDir->mkdir("full");
1684 
1685  // Make a stack of 1D projections
1686  struct Rec {
1687  TH2* h;
1688  TH1* s;
1689  TH2* c;
1690  TH1* p;
1691  Style_t sty;
1692  Color_t col;
1693  Float_t siz;
1694  const char* tit;
1695  };
1696  Rec sC = { simC, simZ, result, 0, 32, kMagenta+2, 1.4, "background" };
1697  Rec sS = { simS, simZ, result, 0, 27, kGreen+2, 1.8, "signal" };
1698  Rec sM = { simM, simZ, result, 0, 26, kBlue+2, 1.4, "measured" };
1699  Rec sG = { simG, simZ, 0, 0, 24, kRed+2, 1.4, "generated" };
1700  Rec sB = { simB, simZ, 0, 0, 28, kPink+2, 1.4, "#beta" };
1701  Rec rC = { realC, realZ, result, 0, 23, kMagenta+2, 1.2, "background" };
1702  Rec rS = { realS, realZ, result, 0, 33, kGreen+2, 1.6, "signal" };
1703  Rec rM = { realM, realZ, result, 0, 22, kBlue+2, 1.2, "measured" };
1704  Rec rR = { result,realZ, 0, 0, 20, kRed+2, 1.3, ObsTitle() };
1705  Rec rB = { realB, realZ, 0, 0, 34, kPink+2, 1.4, "#beta" };
1706  Rec sA = { simA, simZ, 0, 0, 30, kSpring+2, 1.4, "#alpha" };
1707  Rec sF = { fiducial,simZ,0, 0, 31, kSpring+2, 1.4, "fiducial" };
1708  Rec rG = { realG, realZ, 0, 0, 24, kBlack, 1.4, "truth" };
1709  Rec* recs[]={ &rR, &sG, &rS, &sS, &rM, &sM, &rC, &sC,
1710  &rB, &sB, &sA, &sF, &rG, 0 };
1711  Rec** ptr = recs;
1712  TH1* dndeta = 0;
1713  THStack* all = new THStack("all", "");
1714  if (fViz & kAltMarker) {
1715  rR.sty = 21;
1716  rR.siz = 1.2;
1717  sG.sty = 25;
1718  sG.siz = 1.3;
1719  rG.sty = 25;
1720  rG.siz = 1.3;
1721  }
1722  while ((*ptr)) {
1723  Rec* src = *ptr;
1724  if (!src->h) { ptr++; continue; }
1725  src->h->SetDirectory(full);
1726  src->p = AverageOverIPz(src->h, src->h->GetName(), 1,
1727  src->s, src->c);
1728  src->p->SetYTitle(src->h->GetZaxis()->GetTitle());
1729  src->p->SetTitle(Form("%s - %s", src->h->GetTitle(), src->tit));
1730  src->p->SetDirectory(outDir);
1731  if (src->h != simB && src->h != realB &&
1732  src->h != simA && src->h != fiducial) all->Add(src->p);
1733  SetAttr(src->p, src->col, src->sty, src->siz);
1734  if (src->h == result) {
1735  dndeta = src->p;
1736  dndeta->SetYTitle(ObsTitle());
1737  }
1738  ptr++;
1739  }
1740  // Show example calculation
1741  Int_t mi = result->GetXaxis()->FindBin(0.);
1742  Int_t mj = result->GetYaxis()->FindBin(0.);
1743  printf("R=G'/[(1-beta')M'](1-beta)M="
1744  "%6.1f /((1-%6.3f)*%6.1f)*(1-%6.3f)*%6.1f="
1745  "%6.3f * %6.3f * %6.1f="
1746  "%6.1f [%6.1f]",
1747  sG.h->GetBinContent(mi,mj), // p->GetBinContent(mi),
1748  sB.h->GetBinContent(mi,mj), // p->GetBinContent(mi),
1749  sM.h->GetBinContent(mi,mj), // p->GetBinContent(mi),
1750  rB.h->GetBinContent(mi,mj), // p->GetBinContent(mi),
1751  rM.h->GetBinContent(mi,mj), // p->GetBinContent(mi),
1752  sA.h->GetBinContent(mi,mj), // p->GetBinContent(mi),
1753  1-rB.h->GetBinContent(mi,mj), // p->GetBinContent(mi),
1754  rM.h->GetBinContent(mi,mj), // p->GetBinContent(mi),
1755  rR.h->GetBinContent(mi,mj), // p->GetBinContent(mb));
1756  rG.h ? rG.h->GetBinContent(mi,mj) : -1);
1757 
1758  if (rG.p) {
1759  TH1* ratio = RatioH(dndeta, rG.p, "closure");
1760  ratio->SetYTitle("Closure test");
1761  ratio->SetDirectory(outDir);
1762  }
1763 
1764  THStack* ratios = new THStack("ratios", "");
1765  ratios->Add(RatioH(rM.p, sM.p, "rMeaured"));
1766  ratios->Add(RatioH(rC.p, sC.p, "rBackground"));
1767  ratios->Add(RatioH(rS.p, sS.p, "rSignal"));
1768  ratios->Add(RatioH(rR.p, sG.p, "rResult"));
1769  ratios->Add(AverageOverIPz(scale, scale->GetName(), 1, realZ, 0));
1770 
1771  TF1* tmp = new TF1("mid", "pol0", -.5, +.5);
1772  dndeta->Fit(tmp, "Q0R+");
1773  TLatex* ltx = new TLatex(0,tmp->GetParameter(0)/2,
1774  Form("%s|_{|#eta|<0.5}=%.1f#pm%.1f",
1775  ObsTitle(),
1776  tmp->GetParameter(0),
1777  tmp->GetParError(0)));
1778 #if 0
1779  printf(" %dD: %6.1f +/- %4.1f (%4.2f)",
1780  deltaDimen, tmp->GetParameter(0), tmp->GetParError(0),
1781  tmp->GetChisquare()/tmp->GetNDF());
1782 #endif
1783 
1784  ltx->SetTextAlign(22);
1785  ltx->SetTextFont(42);
1786  dndeta->GetListOfFunctions()->Add(ltx);
1787 
1788  outDir->cd();
1789  all->Write();
1790  ratios->Write();
1791  realB ->SetDirectory(full);
1792  simB ->SetDirectory(full);
1793  simA ->SetDirectory(full);
1794  fiducial->SetDirectory(full);
1795 
1796  outParent->cd();
1797 
1798  return dndeta;
1799 }
1800 
1801 
1802 //====================================================================
1804 {
1805  fTop->Clear();
1806  fTop->SetNumber(1);
1807  fTop->SetFillColor(kGray); // kTopBackground);
1808  fTop->SetFillStyle(1001);
1809  fTop->SetBorderSize(0);
1810  fTop->SetBorderMode(0);
1811 
1812  fBody->Clear();
1813  fBody->SetNumber(2);
1814  fBody->SetFillColor(0);
1815  fBody->SetFillStyle(0);
1816  fBody->SetBorderSize(0);
1817  fBody->SetBorderMode(0);
1818  fBody->SetTopMargin(0.01);
1819  fBody->SetRightMargin(0.01);
1820  fBody->SetBottomMargin(0.10);
1821  fBody->SetLeftMargin(0.10);
1822  fBody->SetTicks();
1823 
1824  fCanvas->cd();
1825 }
1826 //____________________________________________________________________
1828 {
1829  gStyle->SetOptStat(0);
1830  gStyle->SetOptTitle(0);
1831 
1832  Int_t h = 1000;
1833  Int_t w = h / TMath::Sqrt(2);
1834  if (fViz & kLandscape) {
1835  Int_t t = h;
1836  h = w;
1837  w = t;
1838  }
1839 
1840  fCanvas = new TCanvas("c",outputName,w,h);
1841  fCanvas->SetFillColor(0);
1842  fCanvas->SetBorderSize(0);
1843  fCanvas->SetBorderMode(0);
1844 
1845  if (fViz & kPDF) {
1846  SuppressGuard2 g;
1847  fCanvas->Print(Form("%s/summary.pdf[", outputName.Data()),
1848  Form("pdf %s", (fViz & kLandscape) ? "Landscape" : ""));
1849  }
1850  // if (fViz & kPNG) {
1851  // }
1852  fCanvas->SetLeftMargin (0.10);
1853  fCanvas->SetRightMargin (0.05);
1854  fCanvas->SetTopMargin (0.05);
1855  fCanvas->SetBottomMargin(0.10);
1856 
1857  Float_t dy = 0.05;
1858  fTop = new TPad("top","Top",0,1-dy,1,1,0,0);
1859  // fTop->SetNumber(1);
1860  // fTop->SetFillColor(kTopBackground);
1861  // fTop->SetBorderSize(0);
1862  // fTop->SetBorderMode(0);
1863  fCanvas->cd();
1864  fTop->Draw();
1865 
1866  fBody = new TPad("body","Body",0,0,1,1-dy,0,0);
1867  fBody->SetNumber(2);
1868  // fBody->SetFillColor(0);
1869  // fBody->SetFillStyle(0);
1870  // fBody->SetBorderSize(0);
1871  // fBody->SetBorderMode(0);
1872  fCanvas->cd();
1873  fBody->Draw();
1874 
1875  ClearCanvas();
1876 
1877  fHeader = new TLatex(.5, .5, "Title");
1878  fHeader->SetNDC();
1879  fHeader->SetTextAlign(22);
1880  // fHeader->SetTextColor(kWhite);
1881  fHeader->SetTextFont(62);
1882  fHeader->SetTextSize(0.7);
1883 
1884  fCanvas->cd();
1885 }
1886 //____________________________________________________________________
1888 {
1889  if ((fViz & kPDF) && fCanvas) {
1890  SuppressGuard2 g;
1891  fCanvas->Print(Form("%s/summary.pdf]", fCanvas->GetTitle()),
1892  Form("pdf %s Title:%s",
1893  (fViz & kLandscape) ? "Landscape" : "",
1894  fLastTitle.Data()));
1895  Printf("PDF %s written", fCanvas->GetTitle());
1896  }
1897  if (fCanvas) fCanvas->Close();
1898  fCanvas = 0;
1899 }
1900 
1901 //____________________________________________________________________
1903  const char* shortTitle,
1904  Float_t size)
1905 {
1906  if (fTop) {
1907  fTop->cd();
1908  fHeader->SetTextSize(size);
1909  fHeader->DrawLatex(.5,.5,title);
1910  }
1911  fCanvas->Modified();
1912  fCanvas->Update();
1913  fCanvas->cd();
1914 
1915  if (fViz & kPDF) {
1916  TString tit;
1917  tit.Form("pdf %s Title:%s", (fViz & kLandscape) ? "Landscape" : "",
1918  title);
1919  // Suppress prints
1920  SuppressGuard2 g;
1921  fCanvas->Print(Form("%s/summary.pdf",fCanvas->GetTitle()), tit);
1922  }
1923  static Int_t cnt = 1;
1924  if (fViz & kPNG) {
1925  SuppressGuard2 g;
1926  fCanvas->Print(Form("%s/%03d_%s.png",fCanvas->GetTitle(),cnt,shortTitle));
1927  cnt++;
1928  }
1929  fLastTitle = title;
1930  if (fViz & kPause) fCanvas->WaitPrimitive();
1931 }
1932 //____________________________________________________________________
1933 TLegend* AliTrackletdNdeta2::DrawInPad(TVirtualPad* c,
1934  Int_t pad,
1935  TObject* o,
1936  Option_t* opt)
1937 {
1938  if (!o) {
1939  // Warning("", "Nothing to draw in pad %d", pad);
1940  return 0;
1941  }
1942  TLegend* l = 0;
1943  TVirtualPad* p = c->cd(pad);
1944  if (!p) {
1945  Warning("DrawInPad", "Sub-pad %d does not exist", pad);
1946  return 0;
1947  }
1948  TString option(opt);
1949  option.ToLower();
1950  if (option.Contains("logx")) { p->SetLogx(); option.ReplaceAll("logx",""); }
1951  if (option.Contains("logy")) { p->SetLogy(); option.ReplaceAll("logy",""); }
1952  if (option.Contains("logz")) { p->SetLogz(); option.ReplaceAll("logz",""); }
1953  if (option.Contains("grid")) { p->SetGridx(); p->SetGridy();
1954  option.ReplaceAll("grid",""); }
1955  Int_t leg = 0;
1956  if (option.Contains("leg3")) { leg = 3; option.ReplaceAll("leg3",""); }
1957  if (option.Contains("leg2")) { leg = 2; option.ReplaceAll("leg2",""); }
1958  if (option.Contains("leg")) { leg = 1; option.ReplaceAll("leg",""); }
1959  // Printf("Drawing %p %s with %s", o, o->GetName(), option.Data());
1960  o->Draw(option);
1961  if (leg) {
1962  l = p->BuildLegend(0.5, 0.73, .98, .98);
1963  l->SetNColumns(leg);
1964  TObject* frame = 0;
1965  TIter next(l->GetListOfPrimitives());
1966  TLegendEntry* ent = 0;
1967  while ((ent = static_cast<TLegendEntry*>(next()))) {
1968  if (TString(ent->GetLabel()).EqualTo("frame")) frame = ent;
1969  }
1970  if (frame) l->GetListOfPrimitives()->Remove(frame);
1971  // l->GetListOfPrimitives()->Print();
1972  }
1973  p->Modified();
1974  // p->Update();
1975  // p->cd();
1976  return l;
1977 }
1978 //____________________________________________________________________
1979 void AliTrackletdNdeta2::ModLegend(TVirtualPad* p, TLegend* l,
1980  Double_t x1, Double_t y1,
1981  Double_t x2, Double_t y2)
1982 {
1983  Double_t px1 = p->GetX1();
1984  Double_t px2 = p->GetX2();
1985  Double_t py1 = p->GetY1();
1986  Double_t py2 = p->GetY2();
1987  l->SetX1(px1+(px2-px1)*x1);
1988  l->SetX2(px1+(px2-px1)*x2);
1989  l->SetY1(py1+(py2-py1)*y1);
1990  l->SetY2(py1+(py2-py1)*y2);
1991  p->Modified();
1992 }
1993 //____________________________________________________________________
1994 THStack* AliTrackletdNdeta2::Make2Stack(const char* name,
1995  const char* title,
1996  Container* realList,
1997  Container* simList,
1998  Option_t* realOpt,
1999  Option_t* simOpt)
2000 {
2001  TString nme(name);
2002  THStack* stack = new THStack(name, title);
2003  TH1* real = CopyH1(realList, name, Form("real%s",name));
2004  TH1* sim = CopyH1(simList, name, Form("sim%s",name));
2005  real->SetMarkerStyle(20);
2006  sim ->SetMarkerStyle(24);
2007  real->SetFillStyle(3004);
2008  sim ->SetFillStyle(3005);
2009  real->SetBarWidth(0.4);
2010  sim ->SetBarWidth(0.4);
2011  real->SetBarOffset(0.1);
2012  sim ->SetBarOffset(0.5);
2013  TString dtit(real->GetTitle());
2014  if (dtit.Contains("\\")) dtit.Form("%s\\hbox{ - real}", real->GetTitle());
2015  else dtit.Form("%s - real", real->GetTitle());
2016  real->SetTitle(dtit);
2017  TString stit(sim->GetTitle());
2018  if (stit.Contains("\\")) stit.Form("%s\\hbox{ - sim.}", sim->GetTitle());
2019  else stit.Form("%s - sim.", sim->GetTitle());
2020  sim->SetTitle(stit);
2021  stack->Add(real, realOpt);
2022  stack->Add(sim, simOpt);
2023  return stack;
2024 }
2025 
2026 
2027 //====================================================================
2029  Container* simSums,
2030  Container* realRess,
2031  Container* simRess,
2032  TDirectory* outDir,
2033  Int_t maxBins)
2034 {
2035  // --- Visualization -----------------------------------------------
2036  TH1* realCent = static_cast<TH1*>(GetO(outDir, "realCent"));
2037  if (!realCent) {
2038  Warning("Visualize", "realCent histogram not found");
2039  return false;
2040  }
2041  TString outName(gSystem->DirName(outDir->GetName()));
2042  outName.ReplaceAll(".root", "");
2043  CreateCanvas(outName);
2044  VisualizeParams(realSums, simSums);
2045  VisualizeGeneral(realRess, simRess);
2046  VisualizeWeights(simRess);
2047  for (Int_t i = 0; i < 4; i++) {
2048  if ((fProc & (1 << i)) == 0) continue;
2049  VisualizeFinal(outDir, i);
2050  VisualizeClosure(outDir, i);
2051  }
2052  for (Int_t i = 1; i <= realCent->GetNbinsX() && i <= maxBins ; i++) {
2053  Double_t c1 = realCent->GetXaxis()->GetBinLowEdge(i);
2054  Double_t c2 = realCent->GetXaxis()->GetBinUpEdge (i);
2055 
2056  VisualizeBin(c1, c2, simRess, outDir);
2057  }
2058  CloseCanvas();
2059 }
2060 
2061 //====================================================================
2063  Container* simList)
2064 {
2065  THStack* ipz = Make2Stack("ipz", "IP_{#it{z}}", realList,simList);
2066  THStack* cent = Make2Stack("cent", "Centrality [%]", realList,simList);
2067  THStack* status = Make2Stack("status","Task status", realList,simList,
2068  "B text90", "B text90");
2069  THStack* centTr = Make2Stack("centTracklets", "#LTtracklets#GT",
2070  realList, simList, "E", "E");
2071  ClearCanvas();
2072  fBody->Divide(1,3);
2073  for (Int_t i = 1; i <= 3; i++) {
2074  if (i < 3) fBody->GetPad(i)->SetRightMargin(0.01);
2075  fBody->GetPad(i)->SetTopMargin(0.01);
2076  }
2077  TLegend* l = 0;
2078  TVirtualPad* q = fBody->GetPad(1);
2079  q->Divide(2,1);
2080  l = DrawInPad(q,1,ipz, "nostack leg");
2081  ModLegend(q->GetPad(1),l,.4,.1,.75,.4);
2082  l = DrawInPad(q,2,cent, "nostack leg");
2083  ModLegend(q->GetPad(2),l,.6,.1,.99,.4);
2084  q = fBody->GetPad(2);
2085  q->Divide(2,1);
2086  l = DrawInPad(q,1,status, "nostack hist text90 leg");
2087  ModLegend(q->GetPad(1),l,.5,.7,.99,.99);
2088  l = DrawInPad(q,2,centTr, "nostack leg");
2089  ModLegend(q->GetPad(2),l,.5,.7,.99,.99);
2090 
2091  TH2* real = GetH2(realList, "etaPhi");
2092  TH2* sim = GetH2(simList, "etaPhi");
2093  if (sim) {
2094  sim->SetMarkerColor(kBlack);
2095  sim->SetMarkerStyle(0);
2096  sim->SetMarkerSize(1);
2097  sim->SetLineColor(kBlack);
2098  sim->SetFillColor(kBlack);
2099  sim->SetFillStyle(0);
2100  sim->SetName("simEtaPhi");
2101  }
2102  DrawInPad(fBody, 3, real, "colz");
2103  DrawInPad(fBody, 3, sim, "box same");
2104 
2105  PrintCanvas("General information","general");
2106 }
2107 
2108 namespace {
2109  void SetCentColors(THStack* s, TH1* dist=0)
2110  {
2111  if (!s || !s->GetHists()) return;
2112 
2113  const Color_t cc[] = { kMagenta+2, // 0
2114  kBlue+2, // 1
2115  kAzure-1, // 2 // 10,
2116  kCyan+2, // 3
2117  kGreen+1, // 4
2118  kSpring+5, // 5 //+10,
2119  kYellow+1, // 6
2120  kOrange+5, // 7 //+10,
2121  kRed+1, // 8
2122  kPink+5, // 9 //+10,
2123  kBlack }; // 10
2124  TIter next(s->GetHists());
2125  TH1* h = 0;
2126  Int_t i = 0;
2127  Double_t min = +10000;
2128  Double_t max = -10000;
2129  while ((h = static_cast<TH1*>(next()))) {
2130  Color_t c = cc[i % 10];
2131  h->SetMarkerColor(c);
2132  h->SetFillColor(c);
2133  h->SetLineColor(c);
2134  h->SetMarkerStyle(20+(i%4));
2135  h->SetDirectory(0);
2136  min = TMath::Min(h->GetMinimum(), min);
2137  max = TMath::Max(h->GetMaximum(), max);
2138  i++;
2139  if (!dist) continue;
2140  h->SetTitle(Form("%5.1f-%5.1f%%",
2141  dist->GetXaxis()->GetBinLowEdge(i),
2142  dist->GetXaxis()->GetBinUpEdge(i)));
2143 
2144  }
2145  s->SetMinimum(min*.9);
2146  s->SetMaximum(max*1.1);
2147  }
2148  THStack* GetPdgStack(AliTrackletAODUtils::Container* w, const char* name)
2149  {
2151  if (!c) return 0;
2152 
2153  THStack* s = new THStack(name, "");
2154  TH1* h = 0;
2155  TIter n(c);
2156  while ((h = static_cast<TH1*>(n()))) s->Add(h);
2157 
2158  if (!s->GetHists() || s->GetHists()->GetEntries() <= 0) {
2159  delete s;
2160  s = 0;
2161  }
2162  return s;
2163  }
2164 }
2165 //____________________________________________________________________
2167 {
2168  TH1* c = GetH1(simList, "cent");
2169  Container* w = GetC(simList,"weights", false);
2170  if (!w) return;
2171 
2172  Float_t right = .75;
2173  ClearCanvas();
2174  fBody->SetTopMargin(0.01);
2175  fBody->Divide(1,3);
2176  TVirtualPad* pp[] = { fBody->GetPad(1), fBody->GetPad(2) };
2177  for (Int_t i = 0; i < 2; i++) {
2178  pp[i]->SetRightMargin(0.01);
2179  pp[i]->SetTopMargin(0.01);
2180  pp[i]->SetPad(pp[i]->GetXlowNDC(), pp[i]->GetYlowNDC(),
2181  right, pp[i]->GetYlowNDC()+pp[i]->GetHNDC());
2182  pp[i]->Modified();
2183  }
2184  TH2* hp = GetH2(w, "centPt");
2185  THStack* ef = new THStack(GetP2(simList,"etaWeight"),"x","effWeights","");
2186  THStack* ab = GetPdgStack(w, "abundance");
2187  THStack* st = GetPdgStack(w, "strangeness");
2188  SetCentColors(ef, c);
2189  SetCentColors(ab);
2190  SetCentColors(st);
2191 
2192  Double_t eMin = +1e9;
2193  Double_t eMax = -1e9;
2194  TIter next(ef->GetHists());
2195  TH1* h = 0;
2196  while ((h = static_cast<TH1*>(next()))) {
2197  eMin = TMath::Min(h->GetMinimum(), eMin);
2198  eMax = TMath::Max(h->GetMaximum(), eMax);
2199  }
2200  ef->SetMinimum(eMax);
2201  ef->SetMinimum(eMin);
2202  // ef->SetMaximum(1.02);
2203  TLegend* l = DrawInPad(fBody, 1, ef, "nostack p leg");
2204  ef->GetHistogram()->SetYTitle("Average weight");
2205  ef->GetHistogram()->SetXTitle("#eta");
2206  fBody->GetPad(1)->GetListOfPrimitives()->Remove(l);
2207  fBody->GetPad(1)->Modified();
2208 
2209  if (hp) {
2210  if (hp->GetNbinsY() > 1) {
2211  THStack* pt = new THStack(hp, "y","pt","");
2212  SetCentColors(pt);
2213  DrawInPad(fBody, 2, pt, "nostack");
2214  pt->GetHistogram()->SetYTitle("Weight");
2215  pt->GetHistogram()->SetXTitle("#it{p}_{T}");
2216  }
2217  else {
2218  TArrayD bins(hp->GetNbinsX()+1);
2219  bins[0] = hp->GetXaxis()->GetBinLowEdge(1);
2220  for (Int_t i = 1; i <= hp->GetNbinsX(); i++)
2221  bins[i] = hp->GetXaxis()->GetBinUpEdge(1);
2222  TH1* pt = new TH1D("pt","", bins.GetSize()-1,bins.GetArray());
2223  for (Int_t i = 1; i <= hp->GetNbinsX(); i++) {
2224  pt->SetBinContent(i, hp->GetBinContent(i,1));
2225  pt->SetBinError (i, hp->GetBinError (i,1));
2226  }
2227  pt->SetYTitle("Weight");
2228  pt->SetXTitle("Centrality [%]");
2229  pt->SetMarkerStyle(2);
2230  pt->SetMarkerColor(kRed+2);
2231  DrawInPad(fBody, 2, pt, "");
2232  }
2233  }
2234 
2235  fBody->cd();
2236  if (l) {
2237  l->Draw();
2238  ModLegend(fBody, l, right, pp[1]->GetYlowNDC(),
2239  .99, 1-fBody->GetTopMargin());
2240  fBody->Modified();
2241  }
2242 
2243  TVirtualPad* p3 = fBody->GetPad(3);
2244  p3->SetTopMargin(0.01);
2245  p3->SetRightMargin(0.01);
2246  p3->Divide(2,1);
2247  p3->GetPad(1)->SetRightMargin(0.01);
2248  p3->GetPad(2)->SetRightMargin(0.01);
2249  p3->GetPad(1)->SetTopMargin(0.01);
2250  p3->GetPad(2)->SetTopMargin(0.01);
2251 
2252  DrawInPad(p3, 1, ab, "nostack leg");
2253  DrawInPad(p3, 2, st, "nostack leg");
2254 
2255  if (ab) {
2256  ab->GetHistogram()->SetYTitle("Weight");
2257  ab->GetHistogram()->SetXTitle("Centrality [%]");
2258  }
2259  if (st) {
2260  st->GetHistogram()->SetYTitle("Weight");
2261  st->GetHistogram()->SetXTitle("Centrality [%]");
2262  }
2263  p3->GetPad(1)->Modified();
2264  p3->GetPad(2)->Modified();
2265 
2266  PrintCanvas("Simulation weights","weights");
2267 }
2268 
2269 //____________________________________________________________________
2270 void AliTrackletdNdeta2::VisualizeFinal(TDirectory* outDir, Int_t i)
2271 {
2272  TDirectory* dd = outDir->GetDirectory(Form("final%dd", i));
2273  if (!dd) return;
2274 
2275  THStack* all = static_cast<THStack*>(GetO(dd, "full"));
2276  TH1* mid = static_cast<TH1*> (GetO(dd, "mid"));
2277  TH1* pub = static_cast<TH1*> (GetO(outDir, "published"));
2278  if (mid->GetMinimum() <= 0) mid->SetMinimum(all->GetMinimum("nostack"));
2279  Double_t max = TMath::Max(all->GetMaximum("nostack"),mid->GetMaximum());
2280  Double_t min = TMath::Min(all->GetMinimum("nostack"),mid->GetMinimum());
2281  all->SetMinimum(.9*min);
2282  mid->SetMinimum(.9*min);
2283  all->SetMaximum(1.2*max);
2284  mid->SetMaximum(1.2*max);
2285  mid->SetLineColor(kBlack);
2286  mid->SetFillColor(kBlack);
2287  mid->SetMarkerColor(kBlack);
2288 
2289  ClearCanvas();
2290 
2291  TPad* p1 = new TPad("p1","p1",0,0,.4,1);
2292  p1->SetTopMargin(0.01);
2293  p1->SetRightMargin(0.0);
2294  p1->SetLeftMargin(0.12);
2295  p1->SetBottomMargin(0.15);
2296  p1->SetTicks();
2297  fBody->cd();
2298  p1->Draw();
2299  p1->SetNumber(2);
2300 
2301  Double_t right = .7;
2302  TPad* p2 = new TPad("p2","p2",.4,0,1,1);
2303  p2->SetTopMargin(0.01);
2304  p2->SetRightMargin(1-right);
2305  p2->SetLeftMargin(0.0);
2306  p2->SetBottomMargin(0.15);
2307  p2->SetTicks();
2308  fBody->cd();
2309  p2->Draw();
2310  p2->SetNumber(2);
2311 
2312 
2313  DrawInPad(p1,0, mid, "logy grid");
2314  TLegend* l = DrawInPad(p1,0, pub, "logy grid same leg");
2315  ModLegend(p1, l, .4, .90, .99, .99);
2316  l->SetMargin(0.2);
2317  l->SetEntrySeparation(0.1);
2318  l->SetTextSize(0.04);
2319 
2320  l = DrawInPad(p2,0, all, "nostack logy grid leg");
2321  all->GetHistogram()->SetXTitle("#eta");
2322 
2323  ModLegend(p2, l, right, .15, .99, .99);
2324  l->SetMargin(0.2);
2325  l->SetEntrySeparation(0.1);
2326  l->SetTextSize(0.04);
2327  p1->Modified();
2328  p2->Modified();
2329  fBody->Modified();
2330 
2331  const char* what = (i == 3 ? "d^{3}N/(d#Deltad#etadIP_{z})" :
2332  i == 2 ? "d^{2}N/(d#Deltad#eta)" :
2333  i == 1 ? "dN/d#Delta" : "dN/d#Delta (k#equiv1)");
2334  PrintCanvas(Form("Results #topbar %s", what), "results");
2335 }
2336 
2337 //____________________________________________________________________
2338 void AliTrackletdNdeta2::VisualizeClosure(TDirectory* outDir, Int_t i)
2339 {
2340  TDirectory* dd = outDir->GetDirectory(Form("final%dd", i));
2341  if (!dd) return;
2342 
2343  THStack* all = static_cast<THStack*>(GetO(dd, "closure"));
2344  if (!all) return;
2345 
2346  Double_t max = all->GetMaximum("nostack");
2347  Double_t min = all->GetMinimum("nostack");
2348  all->SetMinimum(.95*min);
2349  all->SetMaximum(1.05*max);
2350 
2351  ClearCanvas();
2352  Double_t right = .3;
2353  fBody->SetRightMargin(right);
2354  TLegend* l = DrawInPad(fBody,0, all, "grid nostack leg");
2355  ModLegend(fBody, l, 1-right+.01, fBody->GetBottomMargin(), .99,
2356  1-fBody->GetTopMargin());
2357  l->SetMargin(0.2);
2358  l->SetEntrySeparation(0.1);
2359  l->SetTextSize(0.04);
2360  l->SetBorderSize(0);
2361  fBody->Modified();
2362 
2363  const char* what = (i == 3 ? "d^{3}N/(d#Deltad#etadIP_{z})" :
2364  i == 2 ? "d^{2}N/(d#Deltad#eta)" :
2365  i == 1 ? "dN/d#Delta" : "dN/d#Delta (k#equiv1)");
2366  PrintCanvas(Form("Closure #topbar %s", what), "closure");
2367 }
2368 
2369 //====================================================================
2371  Double_t& y,
2372  const char* val)
2373 {
2374  TLatex* ln = new TLatex(.49, y, name);
2375  ln->SetTextAlign(31);
2376  ln->SetTextSize(0.02/gPad->GetHNDC());
2377  ln->SetNDC();
2378  ln->Draw();
2379  TLatex* lv = new TLatex(.51, y, val);
2380  lv->SetTextAlign(11);
2381  lv->SetTextSize(0.02/gPad->GetHNDC());
2382  lv->SetNDC();
2383  lv->Draw();
2384  y -= 0.025/gPad->GetHNDC();
2385 }
2386 //____________________________________________________________________
2388  Double_t& y,
2389  Double_t val)
2390 {
2391  VisualizeParam(name, y, Form("%f", val));
2392 }
2393 //____________________________________________________________________
2395  Double_t& y,
2396  Int_t val)
2397 {
2398  VisualizeParam(name, y, Form("%d", val));
2399 }
2400 //____________________________________________________________________
2402  Double_t& y,
2403  Bool_t val)
2404 {
2405  VisualizeParam(name, y, val ? "yes" : "no");
2406 }
2407 //____________________________________________________________________
2409  const char* title)
2410 {
2411  // c->Clear();
2412  Double_t y = .9;
2413  TLatex* latex = new TLatex(.5, y, title);
2414  latex->SetTextAlign(21);
2415  latex->SetTextSize(0.023/gPad->GetHNDC());
2416  latex->SetNDC();
2417  latex->Draw();
2418  y -= 0.028/gPad->GetHNDC();
2419  if (!pars) return;
2420  VisualizeParam("#delta#phi shift", y, GetD(pars, "DPhiShift"));
2421  VisualizeParam("Shifted #delta#phi cut",y, GetD(pars, "ShiftedDPhiCut"));
2422  VisualizeParam("#Delta cut", y, GetD(pars, "DeltaCut"));
2423  VisualizeParam("max#Delta", y, GetD(pars, "MaxDelta"));
2424  VisualizeParam("min#Delta_{tail}", y, GetD(pars, "TailDelta"));
2425  VisualizeParam("max#Delta_{tail}", y, GetD(pars, "TailMax"));
2426  VisualizeParam("abs.min#it{c}", y, GetD(pars, "AbsMinCent"));
2427 }
2428 //____________________________________________________________________
2430  Container* simSums)
2431 {
2432  ClearCanvas();
2433  fBody->Divide((fViz & kLandscape ? 3 : 1),
2434  (fViz & kLandscape ? 1 : 3), 0, 0);
2435 
2436  TVirtualPad* p1 = fBody->GetPad(1);
2437  TVirtualPad* p2 = fBody->GetPad(2);
2438  TVirtualPad* p3 = fBody->GetPad(3);
2439 #if 0
2440  if (!(fViz & kLandscape)) {
2441  Double_t yr = .2;
2442  p1->SetPad(0,1-yr,1,1);
2443  p2->SetPad(0,(1-yr)/2,1,1-yr);
2444  p3->SetPad(0,0,1,(1-yr)/2);
2445  }
2446 #endif
2447 
2448  // Post-processing stuff
2449  fBody->cd(1);
2450  Double_t y = .80;
2451  TLatex* latex = new TLatex(.5, y, "Post-processing");
2452  latex->SetTextAlign(21);
2453  latex->SetTextSize(0.023/p1->GetHNDC());
2454  latex->SetNDC();
2455  latex->Draw();
2456  y -= 0.028/p1->GetHNDC();
2457  TString scaleM("");
2458  if (fProc & kUnitScale) scaleM.Append(" unit");
2459  if (fProc & kAverageScale) scaleM.Append(" const");
2460  if (fProc & kEtaScale) scaleM.Append(" dN/d#eta");
2461  if (fProc & kFullScale) scaleM.Append(" d^{2}N/(d#etadIP_{z})");
2462  TString tmp = scaleM.Strip(TString::kBoth, ' ');
2463  VisualizeParam("Scaling of comb. bg.", y, tmp);
2464  VisualizeParam("min#alpha", y, fMinAlpha);
2465  VisualizeParam("max#alpha", y, fMaxAlpha);
2466  VisualizeParam("min#it{k}", y, fMinK);
2467  VisualizeParam("max#it{k}", y, fMaxK);
2468  if (fFudge != 1) VisualizeParam("Fudge", y, fFudge);
2469 
2470  // From tasks
2471  fBody->cd(2);
2472  VisualizeParams(GetC(realSums, "parameters"), "Real data");
2473  fBody->cd(3);
2474  VisualizeParams(GetC(simSums, "parameters"), "Simulated data");
2475  PrintCanvas("Parameters", "parameters");
2476 }
2477 
2478 //====================================================================
2480  Double_t c2,
2481  Container* simList,
2482  TDirectory* outTop)
2483 {
2484  // Form the folder name
2485  TString centName(CentName(c1,c2));
2486  fLastBin.Form("%.1f#minus%.1f%%", c1, c2);
2487  fLastShort = centName;
2488 
2489  TDirectory* outDir = outTop->GetDirectory(centName);
2490  if (!outDir) {
2491  Warning("VisualizeBin", "Directory %s not found in %s",
2492  centName.Data(), outTop->GetName());
2493  return false;
2494  }
2495 
2496  Printf("%5.1f - %5.1f%%", c1, c2);
2497  VisualizeSpecies(GetC(simList, centName));
2498  VisualizePrimary(GetC(GetC(simList, centName), "generated"));
2499  for (Int_t i = 0; i < 4; i++) {
2500  if ((fProc & (1 << i)) == 0) continue;
2501  VisualizeDelta(outDir, i);
2502  VisualizeDetails(outDir, i);
2503  VisualizeResult(outDir, i);
2504  }
2505 
2506  return true;
2507 }
2508 //____________________________________________________________________
2510 {
2511  if (!simCont) return true;
2512  ClearCanvas();
2513  Container* meas = GetC(simCont, "measured");
2514  Container* comb = GetC(simCont, "combinatorics");
2515  Container* prim = GetC(simCont, "primaries");
2516  Container* seco = GetC(simCont, "secondaries");
2517  Container* cont[] = { meas, comb, prim, seco };
2518  const char* tit[] = { "M' by primary mother's specie",
2519  "C' by primary mother's specie",
2520  "P' by species",
2521  "S' by primary mother's specie" };
2522  fBody->SetTopMargin(0.01);
2523  fBody->SetRightMargin(0.01);
2524  fBody->Divide(2,2,0.01,0.01);
2525 
2526  for (Int_t i = 0; i < 4; i++) {
2527  if (!cont[i]) continue;
2528 
2529  Container* species = GetC(cont[i], "types");
2530  if (!species) continue;
2531 
2532  THStack* all = static_cast<THStack*>(GetO(species, "all"));
2533  THStack* toPion = static_cast<THStack*>(GetO(species, "toPion"));
2534 
2535 
2536  TVirtualPad* p = fBody->GetPad(i+1);
2537  p->SetTopMargin(0.10);
2538  p->SetRightMargin(0.15);
2539  p->SetBottomMargin(0.15);
2540  p->cd();
2541  TLatex* ltx = new TLatex(.5, .99, tit[i]);
2542  ltx->SetTextAlign(23);
2543  ltx->SetTextSize(0.04);
2544  ltx->SetTextFont(42);
2545  ltx->Draw();
2546 
2547  p->Divide(1,2,0,0);
2548  TLegend* l = DrawInPad(p, 1, all, "nostack grid leg");
2549  DrawInPad(p, 2, toPion, "nostack grid");
2550  all->GetHistogram()->SetYTitle("dN_{X}/d#eta");
2551  all->GetHistogram()->SetXTitle("#eta");
2552  all->GetHistogram()->GetYaxis()->SetTitleSize(0.08);
2553  all->GetHistogram()->GetYaxis()->SetTitleOffset(0.6);
2554  all->GetHistogram()->GetYaxis()->SetLabelSize(0.08);
2555  toPion->GetHistogram()->SetYTitle("Relative to #pi^{#pm} mothers");
2556  toPion->GetHistogram()->SetXTitle("#eta");
2557  toPion->GetHistogram()->GetYaxis()->SetTitleSize(0.08);
2558  toPion->GetHistogram()->GetYaxis()->SetTitleOffset(0.6);
2559  toPion->GetHistogram()->GetYaxis()->SetLabelSize(0.08);
2560  toPion->GetHistogram()->GetXaxis()->SetTitleSize(0.08);
2561  toPion->GetHistogram()->GetXaxis()->SetTitleOffset(0.6);
2562  toPion->GetHistogram()->GetXaxis()->SetLabelSize(0.08);
2563  p->GetPad(1)->GetListOfPrimitives()->Remove(l);
2564  p->GetPad(1)->Modified();
2565  p->GetPad(2)->Modified();
2566  p->cd();
2567  l->Draw();
2568 
2569  ModLegend(p, l, .85, p->GetBottomMargin(), .99, 1-p->GetTopMargin());
2570  p->Modified();
2571  }
2572 
2573  PrintCanvas(Form("Species #topbar %s", fLastBin.Data()),
2574  Form("%s_species", fLastShort.Data()));
2575  return true;
2576 }
2577 
2578 //____________________________________________________________________
2580 {
2581  if (!simCont) return true;
2582  ClearCanvas();
2583 
2584  THStack* all = static_cast<THStack*>(GetO(GetC(simCont,"mix"),"all"));
2585  THStack* toPion = static_cast<THStack*>(GetO(GetC(simCont,"mix"),"toPion"));
2586  THStack* toAll = static_cast<THStack*>(GetO(GetC(simCont,"mix"),"toAll"));
2587  TH2* etaPt = GetH2(simCont, "etaPt");
2588  Int_t etaM = etaPt->GetXaxis()->FindBin(-.5);
2589  Int_t etaP = etaPt->GetXaxis()->FindBin(+.5);
2590  TH1* pt = etaPt->ProjectionY("pt", etaM, etaP);
2591  pt->GetYaxis()->SetTitle("d#it{N}/d#it{p}_{T}");
2592  pt->GetYaxis()->SetTitleSize(0.08);
2593  pt->GetYaxis()->SetLabelSize(0.08);
2594  pt->GetYaxis()->SetTitleOffset(0.6);
2595  pt->GetXaxis()->SetTitleSize(0.08);
2596  pt->GetXaxis()->SetLabelSize(0.08);
2597  pt->GetXaxis()->SetTitleOffset(0.6);
2598 
2599  if (false) {
2600  TIter next(toPion->GetHists());
2601  TH1* rat = 0;
2602  while ((rat = static_cast<TH1*>(next()))) {
2603  TGraphErrors* g =
2604  static_cast<TGraphErrors*>(rat->GetListOfFunctions()
2605  ->FindObject(Form("%s_2760",
2606  rat->GetName())));
2607  TF1* f = new TF1("fit", "pol0", -.5, +.5);
2608  rat->Fit(f, "Q0R+", "", -.5, +.5);
2609  Double_t re, r = RatioE(f->GetParameter(0), f->GetParError(0),
2610  g->GetY()[0], g->GetEY()[0], re);
2611  Printf("%10s: 2760: %6.4f +/- %6.4f Here: %6.4f +/- %6.4f "
2612  "Ratio: %6.4f +/- %6.4f",
2613  rat->GetName(), g->GetY()[0], g->GetEY()[0],
2614  f->GetParameter(0), f->GetParError(0), r, re);
2615 
2616  f->SetLineColor(rat->GetLineColor());
2617  f->SetLineStyle(7);
2618  }
2619  }
2620 
2621  TPad* p1 = new TPad("p1","p1",0,.3,1,1);
2622  p1->SetTopMargin(.01);
2623  p1->SetRightMargin(.01);
2624  fBody->cd();
2625  p1->Draw();
2626  p1->SetNumber(1);
2627  p1->Divide(1,3,0,0);
2628  p1->GetPad(1)->SetRightMargin(0.2);
2629  p1->GetPad(2)->SetRightMargin(0.2);
2630  p1->GetPad(3)->SetRightMargin(0.2);
2631  TLegend* l = DrawInPad(p1,1,all, "leg2 nostack logy");
2632  all->GetHistogram()->GetYaxis()->SetTitle("d#it{N}_{X}/d#eta");
2633  all->GetHistogram()->GetYaxis()->SetTitleSize(0.08);
2634  all->GetHistogram()->GetYaxis()->SetLabelSize(0.08);
2635  all->GetHistogram()->GetYaxis()->SetTitleOffset(0.6);
2636  ModLegend(p1->GetPad(1), l,
2637  1-p1->GetPad(1)->GetRightMargin(),
2638  p1->GetPad(1)->GetBottomMargin(),
2639  1-p1->GetPad(1)->GetTopMargin(),
2640  .99);
2641  l->SetBorderSize(0);
2642 
2643  l = DrawInPad(p1,2,toAll, "nostack leg2 logy");
2644  toAll->GetHistogram()->GetYaxis()->SetTitle("Ratio to all");
2645  toAll->GetHistogram()->GetYaxis()->SetTitleSize(0.08);
2646  toAll->GetHistogram()->GetYaxis()->SetLabelSize(0.08);
2647  toAll->GetHistogram()->GetYaxis()->SetTitleOffset(0.6);
2648  toAll->GetHistogram()->GetXaxis()->SetTitle("#eta");
2649  toAll->GetHistogram()->GetXaxis()->SetTitleSize(0.08);
2650  toAll->GetHistogram()->GetXaxis()->SetLabelSize(0.08);
2651  toAll->GetHistogram()->GetXaxis()->SetTitleOffset(0.6);
2652  ModLegend(p1->GetPad(2), l,
2653  1-p1->GetPad(2)->GetRightMargin(),
2654  p1->GetPad(2)->GetBottomMargin(),
2655  1-p1->GetPad(2)->GetTopMargin(),
2656  .99);
2657  l->SetBorderSize(0);
2658 
2659  l = DrawInPad(p1,3,toPion, "nostack leg");
2660  toPion->GetHistogram()->GetYaxis()->SetTitle("Ratio to #pi");
2661  toPion->GetHistogram()->GetYaxis()->SetTitleSize(0.08);
2662  toPion->GetHistogram()->GetYaxis()->SetLabelSize(0.08);
2663  toPion->GetHistogram()->GetYaxis()->SetTitleOffset(0.6);
2664  toPion->GetHistogram()->GetXaxis()->SetTitle("#eta");
2665  toPion->GetHistogram()->GetXaxis()->SetTitleSize(0.08);
2666  toPion->GetHistogram()->GetXaxis()->SetLabelSize(0.08);
2667  toPion->GetHistogram()->GetXaxis()->SetTitleOffset(0.6);
2668  ModLegend(p1->GetPad(3), l,
2669  1-p1->GetPad(3)->GetRightMargin(),
2670  p1->GetPad(3)->GetBottomMargin(),
2671  1-p1->GetPad(3)->GetTopMargin(),
2672  .99);
2673  l->SetBorderSize(0);
2674 
2675  p1->Modified();
2676  p1->Update();
2677  p1->cd();
2678 
2679  TPad* p2 = new TPad("p2","p2",0,0,1,.3);
2680  p2->SetTopMargin(.01);
2681  p2->SetRightMargin(.01);
2682  p2->SetBottomMargin(0.15);
2683  fBody->cd();
2684  p2->Draw();
2685  p2->SetNumber(2);
2686  DrawInPad(p2,0,pt, "logx logy");
2687  p2->Modified();
2688  p2->Update();
2689  p2->cd();
2690 
2691  fBody->Modified();
2692 
2693  PrintCanvas(Form("Primary species #topbar %s", fLastBin.Data()),
2694  Form("%s_primary_species", fLastShort.Data()));
2695 
2696  return true;
2697 }
2698 
2699 //____________________________________________________________________
2701  Int_t dimen)
2702 {
2703  TDirectory* outDir = outTop->GetDirectory(Form("delta%dd", dimen));
2704  if (!outDir) {
2705  Warning("VisualizeDelta", "Directory detla%dd not found in %s",
2706  dimen, outTop->GetName());
2707  return false;
2708  }
2709 
2710  ClearCanvas();
2711  fBody->cd();
2712  TPad* pq = new TPad("p1","p1",0, .3, 1, 1);
2713  pq->SetNumber(1);
2714  pq->Draw();
2715 
2716  TVirtualPad* q = fBody->cd(1);
2717  q->SetTopMargin(0.01);
2718  q->SetRightMargin(0.01);
2719  q->Divide(1,2,0,0);
2720  q->GetPad(1)->SetRightMargin(0.15);
2721  q->GetPad(2)->SetRightMargin(0.15);
2722  TVirtualPad* qq = q->GetPad(1);
2723  THStack* all = static_cast<THStack*>(GetO(outDir,"all"));
2724  TLegend* l = DrawInPad(q,1,all,"nostack logx logy grid leg");
2725  l->SetBorderSize(0);
2726  l->SetFillStyle(0);
2727  l->SetMargin(0.2);
2728  l->SetEntrySeparation(0.1);
2729  all->GetHistogram()->GetYaxis()->SetTitle("d#it{N}/d#Delta");
2730  all->GetHistogram()->GetYaxis()->SetLabelSize(0.06);
2731  all->GetHistogram()->GetYaxis()->SetTitleSize(0.06);
2732  all->GetHistogram()->GetYaxis()->SetTitleOffset(0.6);
2733  ModLegend(qq, l, 1-qq->GetRightMargin(), qq->GetBottomMargin(), .99,
2734  1-qq->GetTopMargin()-.01);
2735 
2736 
2737  THStack* ratios = static_cast<THStack*>(GetO(outDir,"ratios"));
2738  ratios->SetMinimum(.6);
2739  ratios->SetMaximum(1.4);
2740  qq = q->GetPad(2);
2741  l = DrawInPad(q,2,ratios,"nostack logx grid leg");
2742  l->SetBorderSize(0);
2743  l->SetFillStyle(0);
2744  l->SetMargin(0.2);
2745  l->SetEntrySeparation(0.1);
2746  ratios->GetHistogram()->GetXaxis()->SetTitle("#Delta");
2747  ratios->GetHistogram()->GetYaxis()->SetTitle("Ratio");
2748  ratios->GetHistogram()->GetYaxis()->SetLabelSize(0.06);
2749  ratios->GetHistogram()->GetYaxis()->SetTitleSize(0.06);
2750  ratios->GetHistogram()->GetYaxis()->SetTitleOffset(0.6);
2751  ratios->GetHistogram()->GetXaxis()->SetLabelSize(0.06);
2752  ratios->GetHistogram()->GetXaxis()->SetTitleSize(0.06);
2753  ratios->GetHistogram()->GetXaxis()->SetTitleOffset(0.6);
2754 
2755  ModLegend(qq, l, 1-qq->GetRightMargin(), qq->GetBottomMargin(), .99,
2756  1-qq->GetTopMargin()-.01);
2757 
2758  fBody->cd();
2759  pq = new TPad("p2","p2",0, 0, 1, .3);
2760  pq->SetBottomMargin(0.25);
2761  pq->SetNumber(2);
2762  pq->Draw();
2763 
2764  q = fBody->cd(2);
2765  q->SetTopMargin(0.01);
2766  q->SetRightMargin(0.01);
2767  q->Divide(1,2,0,0);
2768  q->GetPad(1)->SetRightMargin(0.10);
2769  q->GetPad(2)->SetRightMargin(0.10);
2770  TH2* scale = static_cast<TH2*>(GetO(outDir,"scale"));
2771  TH1* scaleProj = static_cast<TH2*>(GetO(outDir,"scaleProj"));
2772  DrawInPad(q,1,scale, "colz");
2773  DrawInPad(q,2,scaleProj,"");
2774  scale->SetYTitle("IP_{#it{z}}");
2775  scaleProj->SetYTitle("#it{k}");
2776  scaleProj->SetXTitle("#eta");
2777  scale->GetYaxis()->SetLabelSize(0.12);
2778  scale->GetYaxis()->SetTitleSize(0.12);
2779  scale->GetYaxis()->SetTitleOffset(0.4);
2780  scaleProj->GetYaxis()->SetLabelSize(0.12);
2781  scaleProj->GetYaxis()->SetTitleSize(0.12);
2782  scaleProj->GetYaxis()->SetTitleOffset(0.4);
2783  scaleProj->GetYaxis()->SetNdivisions(207);
2784  scaleProj->GetXaxis()->SetLabelSize(0.12);
2785  scaleProj->GetXaxis()->SetTitleSize(0.12);
2786  scaleProj->GetXaxis()->SetTitleOffset(0.6);
2787 
2788  const char* what = (dimen == 3 ? "d^{3}N/(d#Deltad#etadIP_{z})" :
2789  dimen == 2 ? "d^{2}N/(d#Deltad#eta)" :
2790  dimen == 1 ? "dN/d#Delta" : "dN/d#Delta (k#equiv1)");
2791  PrintCanvas(Form("%s #topbar %s", what, fLastBin.Data()),
2792  Form("%s_deltas", fLastShort.Data()));
2793 
2794  return true;
2795 }
2796 
2797 //____________________________________________________________________
2799  Int_t dimen)
2800 {
2801  TDirectory* outDir = outTop->GetDirectory(Form("results%dd", dimen));
2802  if (!outDir) {
2803  Warning("VisualizeDelta", "Directory results%dd not found in %s",
2804  dimen, outTop->GetName());
2805  return false;
2806  }
2807  Double_t tbase = 0.015;
2808  ClearCanvas();
2809  fBody->cd();
2810  fBody->SetTopMargin(0.01);
2811  fBody->SetRightMargin(0.01);
2812  fBody->SetBottomMargin(0.15);
2813  fBody->SetLeftMargin(0.15);
2814  fBody->Divide(2,6,0,0);
2815 
2816  const char* names[] = { "realM", "simM",
2817  "realC", "simC",
2818  "realB", "simB",
2819  "realS", "simS",
2820  "simA", "fiducial",
2821  "result", "simG" };
2822  const char* titles[] = { "M", "M'",
2823  "C", "C'",
2824  "#beta", "#beta'",
2825  "S", "S'",
2826  "#alpha", "F'",
2827  "R", "G'" };
2828 
2829  for (Int_t i = 0; i < 12; i++) {
2830  const char* name = names[i];
2831  const char* title = titles[i];
2832  TVirtualPad* pad = fBody->GetPad(i+1);
2833  pad->SetLeftMargin(0.15);
2834  pad->SetRightMargin(0);
2835  pad->Divide(1,2,0,0);
2836  TVirtualPad* q2 = pad->GetPad(1); q2->SetRightMargin(0.15);
2837  TVirtualPad* q1 = pad->GetPad(2); q1->SetRightMargin(0.15);
2838 
2839  TH2* h2 = static_cast<TH2*>(outDir->Get(Form("full/%s", name)));
2840  TH1* h1 = static_cast<TH1*>(outDir->Get(name));
2841  if (!h2 || !h1) {
2842  Warning("VisualizeDetails", "Didn't find full/%s (%p) or %s (%p)",
2843  name, name);
2844  continue;
2845  }
2846  h2 = static_cast<TH2*>(h2->Clone());
2847  h1 = static_cast<TH1*>(h1->Clone());
2848  h2->SetDirectory(0);
2849  h1->SetDirectory(0);
2850  h2->SetXTitle("#eta");
2851  h1->SetXTitle("#eta");
2852  h2->SetYTitle(title);
2853  h1->SetYTitle("");
2854  TAxis* axis[] = { h2->GetXaxis(), h2->GetYaxis(), h2->GetZaxis(),
2855  h1->GetXaxis(), h1->GetYaxis(), 0 };
2856  TVirtualPad* pads[] = { q2, q2, q2, q1, q1, 0 };
2857  TAxis** pa = axis;
2858  TVirtualPad** pq = pads;
2859  while ((*pa)) {
2860  (*pa)->SetTitleSize(tbase/pad->GetHNDC()/(*pq)->GetHNDC());
2861  (*pa)->SetLabelSize(tbase/pad->GetHNDC()/(*pq)->GetHNDC());
2862  (*pa)->SetTitleOffset(0.4);
2863  (*pa)->SetNdivisions(207);
2864  pa++;
2865  pq++;
2866  }
2867  DrawInPad(pad, 1, h2, "colz");
2868  DrawInPad(pad, 2, h1, "");
2869  pad->Modified();
2870  }
2871  fBody->Modified();
2872 
2873  PrintCanvas(Form("Details #topbar %s", fLastBin.Data()),
2874  Form("%s_details", fLastShort.Data()));
2875  return true;
2876 }
2877 
2878 //____________________________________________________________________
2880  Int_t dimen)
2881 {
2882  TDirectory* outDir = outTop->GetDirectory(Form("results%dd", dimen));
2883  if (!outDir) {
2884  Warning("VisualizeDelta", "Directory results%dd not found in %s",
2885  dimen, outTop->GetName());
2886  return false;
2887  }
2888 
2889  ClearCanvas();
2890  Double_t tbase = 0.03;
2891  Double_t yr = .2;
2892  Double_t yf = .2;
2893  TPad* p1 = new TPad("p1","p1",0,yf,1,1);
2894  p1->SetTopMargin(yr);
2895  p1->SetRightMargin(0.01);
2896  p1->SetLeftMargin(0.15);
2897  p1->SetBottomMargin(0);
2898  p1->SetTicks();
2899  fBody->cd();
2900  p1->Draw();
2901  p1->SetNumber(1);
2902 
2903  TPad* p2 = new TPad("p2","p2",0,0,1,yf);
2904  p2->SetTopMargin(0.01);
2905  p2->SetRightMargin(.01);
2906  p2->SetLeftMargin(0.15);
2907  p2->SetBottomMargin(0.20);
2908  p2->SetTicks();
2909  fBody->cd();
2910  p2->Draw();
2911  p2->SetNumber(2);
2912 
2913  THStack* all = static_cast<THStack*>(GetO(outDir, "all"));
2914  TLegend* l = DrawInPad(fBody,1,all, "nostack leg2");
2915  all->GetHistogram()->SetXTitle("#eta");
2916  all->GetHistogram()->SetYTitle(ObsTitle());
2917  all->GetHistogram()->GetYaxis()->SetTitleOffset(1.7);
2918  all->GetHistogram()->GetYaxis()->SetTitleSize(tbase/(1-yf));
2919  all->GetHistogram()->GetYaxis()->SetLabelSize(tbase/(1-yf));
2920  all->GetHistogram()->GetYaxis()->SetNdivisions(205);
2921  l->SetBorderSize(0);
2922  l->SetFillStyle(0);
2923  l->SetMargin(0.12);
2924  l->SetEntrySeparation(0.1);
2925  l->SetHeader("R=G'/(M'-C')#times(M-C) #kern[2]{ } [C=kC'/M'#timesM]");
2926  // l->SetTextSize(0.04);
2927  // l->SetTextAlign(12);
2928  ModLegend(p1, l, p1->GetLeftMargin()-.01, 1-yr,
2929  1-p1->GetRightMargin(), .99);
2930 
2931  THStack* ratios = static_cast<THStack*>(GetO(outDir, "ratios"));
2932  FixMinMax(ratios, true);
2933  DrawInPad(fBody,2,ratios, "nostack");
2934  ratios->GetHistogram()->SetXTitle("#eta");
2935  ratios->GetHistogram()->SetYTitle("Ratios");
2936  ratios->GetHistogram()->GetYaxis()->SetTitleOffset(.45);
2937  ratios->GetHistogram()->GetXaxis()->SetTitleOffset(.7);
2938  ratios->GetHistogram()->GetYaxis()->SetTitleSize(tbase/yr);
2939  ratios->GetHistogram()->GetYaxis()->SetLabelSize(tbase/yr);
2940  ratios->GetHistogram()->GetXaxis()->SetTitleSize(tbase/yr);
2941  ratios->GetHistogram()->GetXaxis()->SetLabelSize(tbase/yr);
2942  ratios->GetHistogram()->GetYaxis()->SetNdivisions(205);
2943 
2944  fBody->Modified();
2945  // fBody->GetListOfPrimitives()->Print();
2946 
2947  const char* what = (dimen == 3 ? "k(#eta,IP_{z})" :
2948  dimen == 2 ? "k(#eta)" :
2949  dimen == 1 ? "k=const." : "k#equiv1");
2950  PrintCanvas(Form("%s [%s] #topbar %s", ObsTitle(), what, fLastBin.Data()),
2951  Form("%s_summary", fLastShort.Data()));
2952 
2953  return true;
2954 }
2955 
2956 //====================================================================
2958 {
2959  TFile* file = TFile::Open(filename, "READ");
2960  if (!file) {
2961  Warning("OpenFile", "Failed to open \"%s\"", filename);
2962  return 0;
2963  }
2964  return file;
2965 }
2966 //____________________________________________________________________
2968 {
2969  static TString tmp;
2970  tmp.Form("cent%06.2f_%06.2f", c1, c2);
2971  tmp.ReplaceAll(".", "d");
2972  return tmp.Data();
2973 }
2974 
2975 //____________________________________________________________________
2977  Color_t color,
2978  Style_t marker,
2979  Double_t size,
2980  Style_t fill,
2981  Style_t line,
2982  Width_t width)
2983 {
2984  h->SetMarkerColor(color);
2985  h->SetMarkerStyle(marker);
2986  h->SetMarkerSize(size);
2987  h->SetFillColor(color);
2988  h->SetFillStyle(fill);
2989  h->SetLineColor(color);
2990  h->SetLineStyle(line);
2991  h->SetLineWidth(width);
2992  h->GetXaxis()->SetNdivisions(210);
2993  h->GetYaxis()->SetNdivisions(210);
2994 }
2995 
2996 //____________________________________________________________________
2998 {
2999  Double_t sum = 0;
3000  Double_t sumw = 0;
3001  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
3002  Double_t c = h->GetBinContent(i);
3003  Double_t e = h->GetBinError (i);
3004  if (c < 1e-6 || e/c > 0.1) continue;
3005  Double_t w = 1/e/e;
3006  sum += w * c;
3007  sumw += w;
3008  }
3009  e = TMath::Sqrt(1/sumw);
3010  return sum / sumw;
3011 }
3012 //____________________________________________________________________
3014 {
3015  Double_t sum = 0;
3016  Double_t sumw = 0;
3017  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
3018  for (Int_t j = 1; j <= h->GetNbinsY(); j++) {
3019  Double_t c = h->GetBinContent(i,j);
3020  Double_t e = h->GetBinError (i,j);
3021  if (c < 1e-6 || e/c > 0.1) continue;
3022  Double_t w = 1/e/e;
3023  sum += w * c;
3024  sumw += w;
3025  }
3026  }
3027  e = TMath::Sqrt(1/sumw);
3028  return sum / sumw;
3029 }
3030 
3031 
3032 
3033 
3034 #endif
3035 //____________________________________________________________________
3036 //
3037 // EOF
3038 //
Int_t color[]
Bool_t ProcessBin(Double_t c1, Double_t c2, Container *realTop, Container *simTop, TDirectory *outTop)
Bool_t Deltas1D(Container *realCont, Container *simCont, TDirectory *outParent)
const char * filename
Definition: TestFCM.C:1
static TDirectory * GetT(TDirectory *parent, const char *name, Bool_t verb=true)
const Color_t cc[]
Definition: DrawKs.C:1
double Double_t
Definition: External.C:58
const char * title
Definition: MakeQAPdf.C: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
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)
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)
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)
TFile * OpenFile(const char *filename)
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 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)
static TObject * GetO(Container *parent, const char *name, TClass *cls=0, Bool_t verb=true)
void CreateCanvas(const TString &outputName)
Definition: External.C:220
TFile * file
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)
const char * CentName(Double_t c1, Double_t c2)
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)