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