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