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