AliPhysics  ef3b16e (ef3b16e)
DeltaExpectations.C
Go to the documentation of this file.
1 
12 #ifndef __CINT__
13 # include <vector>
14 # include <TFile.h>
15 # include <TError.h>
16 # include <TCollection.h>
17 # include <TString.h>
18 # include <TH1.h>
19 # include <TH2.h>
20 # include <THStack.h>
21 # include <TClass.h>
22 # include <TBrowser.h>
23 # include <TCanvas.h>
24 # include <TLegend.h>
25 # include <TLegendEntry.h>
26 # include <TStyle.h>
27 #else
28 class TFile;
29 class TClass;
30 class TH1;
31 class TH2;
32 class TDirectory;
33 class THStack;
34 class TCollection;
35 class TBrowser;
36 class TCanvas;
37 class TVirtualPad;
38 #endif
39 
45 const char* parFmt = "%-10s | %6.3f |";
46 const char* rowFmt = " %5.2f | %5.2f | %5.2f | %5.2f |";
47 const char* eRowFmt = " %5.2f+/-%5.2f | %5.2f+/-%5.2f | "
48  "%5.2f+/-%5.2f | %5.2f+/-%5.2f |";
49 struct Utilities
50 {
59  static TObject* ChkC(TObject* o, TClass* c)
60  {
61  if (!o) return 0;
62  if (!c) return o;
63  if (!o->IsA()->InheritsFrom(c)) {
64  ::Warning("ChkC", "Object %s is not a %s but a %s",
65  o->GetName(), c->GetName(), o->ClassName());
66  return 0;
67  }
68  return o;
69  }
79  static TObject* ChkO(TObject* src, TObject* o, const char* name)
80  {
81  if (!o) {
82  ::Warning("ChkO", "Object %s not found in %s", name, src->GetName());
83  src->ls();
84  return 0;
85  }
86  return o;
87  }
97  static TObject* GetO(TDirectory* d, const char* name, TClass* cls)
98  {
99  if (!d) {
100  ::Warning("GetO", "No directory passed for %s", name);
101  return 0;
102  }
103  TObject* o = d->Get(name);
104  return ChkC(ChkO(d, o, name), cls);
105  }
115  static TObject* GetO(TCollection* d, const char* name, TClass* cls)
116  {
117  if (!d) {
118  ::Warning("GetO", "No collection passed for %s", name);
119  return 0;
120  }
121  TObject* o = d->FindObject(name);
122  return ChkC(ChkO(d, o, name), cls);
123  }
132  static TCollection* GetC(TDirectory* d, const char* name)
133  {
134  return static_cast<TCollection*>(GetO(d, name, TCollection::Class()));
135  }
144  static TCollection* GetC(TCollection* d, const char* name)
145  {
146  return static_cast<TCollection*>(GetO(d, name, TCollection::Class()));
147  }
156  static TH1* GetH1(TCollection* d, const char* name)
157  {
158  return static_cast<TH1*>(GetO(d, name, TH1::Class()));
159  }
168  static TH2* GetH2(TCollection* d, const char* name)
169  {
170  return static_cast<TH2*>(GetO(d, name, TH2::Class()));
171  }
180  static THStack* GetHS(TCollection* d, const char* name)
181  {
182  return static_cast<THStack*>(GetO(d, name, THStack::Class()));
183  }
184 };
185 
192 {
193  enum EColumn {
197  kCb
198  };
199  enum ERow {
203  kAsIs
204  };
205  enum EKind {
206  kK0S = 1,
211  kOther
212  };
213  struct Row
214  {
231  {
232  return (i == kPs ? Ps :
233  i == kSs ? Ss :
234  i == kCs ? Cs :
235  i == kCb ? Cb : 0);
236  }
245  {
246  return (i == kPs ? ePs :
247  i == kSs ? eSs :
248  i == kCs ? eCs :
249  i == kCb ? eCb : 0);
250  }
257  void Scale(EColumn i, Double_t s)
258  {
259  switch (i) {
260  case kPs: Ps *= s; ePs *= s; break;
261  case kSs: Ss *= s; eSs *= s; break;
262  case kCs: Cs *= s; eCs *= s; break;
263  case kCb: Cb *= s; eCb *= s; break;
264  }
265  }
270  void Print(Bool_t err=false) const
271  {
272  if (!err) printf(rowFmt, Ps, Ss, Cs, Cb);
273  else printf(eRowFmt, Ps, ePs, Ss, eSs, Cs, eCs, Cb, eCb);
274  }
275  };
281  struct Particle
282  {
297  {
298  return (i == kReduced ? reduced :
299  i == kExpectation ? expected :
300  i == kReweighed ? reweighed :
301  asis);
302  }
307  void Print(Bool_t err=false)
308  {
309  printf(parFmt, n.Data(), w);
310  reduced. Print(err);
311  expected. Print(err);
312  reweighed.Print(err);
313  asis. Print(err);
314  printf("\n");
315  }
329  EKind kind,
330  TCollection* bin,
331  const TString& sub,
332  const TString& oth,
333  Bool_t abso=true)
334  {
335  TString pn = "primaries";
336  TString sn = "secondaries";
337  TString cn = "combinatorics";
338  TCollection* pc = GetC(GetC(GetC(bin,pn),"specie"),sub);
339  TCollection* sc = GetC(GetC(GetC(bin,sn),"specie"),sub);
340  TCollection* cc = GetC(GetC(GetC(bin,cn),"specie"),sub);
341  TCollection* po = GetC(GetC(GetC(bin,pn),"specie"),oth);
342  TCollection* so = GetC(GetC(GetC(bin,sn),"specie"),oth);
343  TCollection* co = GetC(GetC(GetC(bin,cn),"specie"),oth);
344  THStack* hp = GetHS(pc,"rows");
345  THStack* hs = GetHS(sc,"rows");
346  THStack* hc = GetHS(cc,"rows");
347  TH1* rps = GetH1(hp->GetHists(), "rowSig");
348  TH1* rss = GetH1(hs->GetHists(), "rowSig");
349  TH1* rcs = GetH1(hc->GetHists(), "rowSig");
350  TH1* rcb = GetH1(hc->GetHists(), "rowBg");
351  TH1* pi = GetH1(pc, "totalIntegrals");
352  TH1* si = GetH1(sc, "totalIntegrals");
353  TH1* ci = GetH1(cc, "totalIntegrals");
354  TH1* pj = GetH1(po, "totalIntegrals");
355  TH1* sj = GetH1(so, "totalIntegrals");
356  TH1* cj = GetH1(co, "totalIntegrals");
357  Double_t ww = (row == kExpectation ? w : 1);
358  Int_t bn = kind;
359  Row& ret = Get(row);
360  if (n.IsNull()) n = rps->GetXaxis()->GetBinLabel(bn);
361  if (abso) {
362  ret.Ps = ww*pi->GetBinContent(2)*rps->GetBinContent(bn);
363  ret.Ss = ww*si->GetBinContent(2)*rss->GetBinContent(bn);
364  ret.Cs = ww*ci->GetBinContent(2)*rcs->GetBinContent(bn);
365  ret.Cb = ww*ci->GetBinContent(3)*rcb->GetBinContent(bn);
366  ret.ePs = ww*pi->GetBinContent(2)*rps->GetBinError (bn);
367  ret.eSs = ww*si->GetBinContent(2)*rss->GetBinError (bn);
368  ret.eCs = ww*ci->GetBinContent(2)*rcs->GetBinError (bn);
369  ret.eCb = ww*ci->GetBinContent(3)*rcb->GetBinError (bn);
370  }
371  else {
372  ret.Ps = ww*rps->GetBinContent(bn);
373  ret.Ss = ww*rss->GetBinContent(bn);
374  ret.Cs = ww*rcs->GetBinContent(bn);
375  ret.Cb = ww*rcb->GetBinContent(bn);
376  ret.ePs = ww*rps->GetBinError (bn);
377  ret.eSs = ww*rss->GetBinError (bn);
378  ret.eCs = ww*rcs->GetBinError (bn);
379  ret.eCb = ww*rcb->GetBinError (bn);
380  }
381  return this;
382  }
383  };
384 #ifndef __CINT__
385  typedef std::vector<Particle*> ParticleList;
386 #else
387  typedef std::vector<void*> ParticleList;
388 #endif
389  ParticleList particles;
390  Double_t fP; // Ratio of primaries to all tracklets
391  Double_t fS; // Ratio of secondaries to all tracklets
392  Double_t fC; // Ratio of combinatorics to all tracklets
393  Double_t fPs; // Ratio of signal primaries to all tracklets
394  Double_t fSs; // Ratio of signal secondaries to all tracklets
395  Double_t fCs; // Ratio of signal combinatorics to all tracklets
400  void Print(Bool_t err=false)
401  {
402  printf("%-18s |", "Sample");
403  if (err) {
404  printf(" %-29s |", "Reduced");
405  printf(" %-29s |", "Expected");
406  printf(" %-29s |", "Reweighed");
407  printf(" %-29s |", "As-is");
408  }
409  else {
410  printf(" %-61s |", "Reduced");
411  printf(" %-61s |", "Expected");
412  printf(" %-61s |", "Reweighed");
413  printf(" %-61s |", "As-is");
414  }
415  printf("\n");
416  TString head(parFmt);
417  for (Int_t i = 0; i < 4; i++) head.Append(rowFmt);
418  head.ReplaceAll("%6.3f", "%6s");
419  head.ReplaceAll("%5.2f", err ? "%13s" : "%5s");
420  TString lne(head);
421  lne.ReplaceAll("|", "+");
422  lne.ReplaceAll("%5s", "-----");
423  lne.ReplaceAll("%6s", "------");
424  lne.ReplaceAll("%10s", "----------");
425  lne.ReplaceAll("%13s", "-------------");
426  lne.ReplaceAll(" ", "-");
427  Printf(head.Data(), "Type", "W",
428  "Ps", "Ss", "Cs", "Cb",
429  "Ps", "Ss", "Cs", "Cb",
430  "Ps", "Ss", "Cs", "Cb",
431  "Ps", "Ss", "Cs", "Cb");
432  Printf(lne);
433  for (ParticleList::const_iterator i = particles.begin();
434  i != particles.end(); ++i)
435  (*i)->Print(err);
436  Printf(lne);
437  }
443  {
444  for (Int_t i = 0; i < 4; i++) {
445  EColumn c = (i == 0 ? kPs :
446  i == 1 ? kSs :
447  i == 2 ? kCs : kCb);
448  Double_t sum = 0;
449  for (ParticleList::const_iterator j = particles.begin();
450  j != particles.end(); ++j)
451  sum += (*j)->expected.Get(c);
452  for (ParticleList::const_iterator j = particles.begin();
453  j != particles.end(); ++j)
454  (*j)->expected.Scale(c,100/sum);
455  }
456  }
457 
458 
475  EKind kind,
476  TCollection* reduced,
477  TCollection* reweighed,
478  TCollection* asis,
479  const TString& sub,
480  const TString& oth,
481  Bool_t abso)
482  {
483  Particle* ret = new Particle;
484  ret->w = w;
485  ret->Read(kReduced, kind, reduced, sub, oth, abso);
486  ret->Read(kExpectation, kind, reduced, sub, oth, abso);
487  ret->Read(kReweighed, kind, reweighed, sub, oth, abso);
488  ret->Read(kAsIs, kind, asis, sub, oth, abso);
489  particles.push_back(ret);
490  return ret;
491  }
501  THStack* MakeStack(EColumn column,
502  const char* name,
503  const char* title,
504  Int_t mode)
505  {
506  const char* names[] = { "red", "exp", "rew", "asi" };
507  const char* titles[] = { "Reduced", "Expected", "Reweighed", "As-is" };
508  const Color_t colors[] = { kRed+1, kMagenta+1, kBlue+1, kGreen+1 };
509  Int_t nPart = particles.size();
510  Double_t bW = 0.15;
511  Double_t bS = 0.05;
512  Double_t bT = 4*bW+3*bS;
513  Double_t bO = (1-bT)/2;
514  THStack* stack = new THStack(name, title);
515  for (Int_t i = 0; i < 4; i++) {
516  ERow q = (i == 0 ? kReduced :
517  i == 1 ? kExpectation :
518  i == 2 ? kReweighed : kAsIs);
519  TH1* h = new TH1D(names[i],titles[i],nPart,0,nPart);
520  h->SetFillColor(colors[i]);
521  h->SetLineColor(colors[i]);
522  h->SetMarkerColor(colors[i]);
523  h->SetBarOffset(bO+i*(bW+bS));
524  h->SetMarkerSize(mode == 2 ? 3.5 : 2);
525  h->SetBarWidth (bW);
526  h->SetDirectory(0);
527  h->SetYTitle("Fraction [%]");
528  stack->Add(h, "bar0 text90");
529 
530  for (Int_t j = 0; j < nPart; j++) {
531  Particle* p = particles[j];
532  Row& r = p->Get(q);
533  h->GetXaxis()->SetBinLabel(j+1,p->n);
534  h->SetBinContent (j+1,r.Get (column));
535  h->SetBinError (j+1,r.GetE(column));
536  }
537  }
538  return stack;
539  }
549  void DrawStack(TVirtualPad* mother,
550  EColumn column,
551  THStack* stack,
552  Bool_t logy,
553  Bool_t abso,
554  Int_t mode)
555  {
556  TVirtualPad* p = mother->cd(column+1);
557  p->SetGridx();
558  p->SetGridy();
559  p->SetTicks();
560  if (logy) p->SetLogy();
561  if (mode==2 || (p->GetNumber() % (mode==1 ? 4 : 2)) == 0)
562  p->SetRightMargin(0.01);
563  if (!abso) {
564  stack->SetMinimum(0.003);
565  stack->SetMaximum(logy ? (mode == 2 ? 900 : 400) : 110);
566  }
567  gStyle->SetTitleFontSize(0.02/p->GetHNDC());
568  stack->Draw("nostack hist");
569  stack->GetHistogram()->GetYaxis()->SetNdivisions(210);
570  stack->GetHistogram()->GetYaxis()->SetLabelSize(0.03/p->GetHNDC());
571  stack->GetHistogram()->GetYaxis()->SetTitleSize(0.03/p->GetHNDC());
572  stack->GetHistogram()->GetXaxis()->SetLabelSize(0.03/p->GetHNDC());
573  TString tit = abso ? "#it{N}_{tracklet}" : "Fraction [%]";
574  if ((mode == 2 || mode == 0) && p->GetNumber() != 1) tit="";
575  stack->GetHistogram()->GetYaxis()->SetTitle(tit);
576  stack->GetHistogram()->GetYaxis()->SetTitleOffset(mode == 2 ? 0.5 : 1);
577  }
578 
586  void Run(Double_t c1=0, Double_t c2=0, Bool_t mid=true)
587  {
588  Bool_t abso = false;
589  Bool_t logy = true;
590  Int_t mode = 2; // 0: square, 1: landscape, 2: portrait
591 
592  TString bin("all");
593  if (c2 > c1)
594  bin.Form("cent%03dd%02d_%03dd%02d",
595  Int_t(c1), Int_t(100*c1) % 100,
596  Int_t(c2), Int_t(100*c2) % 100);
597 
598  TFile* redFile = TFile::Open("stk.root", "READ");
599  TFile* rweFile = TFile::Open("wstk.root", "READ");
600  TFile* trhFile = TFile::Open("hijing.root", "READ");
601  if (!redFile || !rweFile || !trhFile) return;
602  TString sub = (mid ? "mid" : "fwd");
603  TString oth = (mid ? "fwd" : "mid");
604  TCollection* redTop = GetC(GetC(redFile,"MidRapidityMCResults"),bin.Data());
605  TCollection* rweTop = GetC(GetC(rweFile,"MidRapidityMCResults"),bin.Data());
606  TCollection* trhTop = GetC(GetC(trhFile,"MidRapidityMCResults"),bin.Data());
607 
608  Create(1.52233, kK0S, redTop, rweTop, trhTop, sub, oth, abso);
609  Create(1.41178, kKpm, redTop, rweTop, trhTop, sub, oth, abso);
610  Create(2.75002, kLambda, redTop, rweTop, trhTop, sub, oth, abso);
611  Create(3.24110, kXi, redTop, rweTop, trhTop, sub, oth, abso);
612  Create(2.75002, kSigma, redTop, rweTop, trhTop, sub, oth, abso);
613  Create(1, kOther, redTop, rweTop, trhTop, sub, oth, abso);
614 
615  if (!abso) ScaleExpected();
616  Print(true);
617 
618  THStack* Ps = MakeStack(kPs, "Ps", "Signal primaries", mode);
619  THStack* Ss = MakeStack(kSs, "Ss", "Signal secondaries", mode);
620  THStack* Cs = MakeStack(kCs, "Cs", "Signal fake", mode);
621  THStack* Cb = MakeStack(kCb, "Cb", "Background fake", mode);
622 
623  gStyle->SetPaintTextFormat("5.2f%%");
624 
625  Int_t cw = (mode == 1 ? 1600 : mode == 2 ? 800 : 1000);
626  Int_t ch = (mode == 1 ? cw/2 : mode == 2 ? 1.2*cw : cw);
627  TCanvas* c = new TCanvas(Form("deltaExpectation_%s_%s",
628  bin.Data(), mid?"mid":"fwd"),
629  "DeltaCanvas",cw,ch);
630  c->SetTopMargin(0.01);
631  c->SetRightMargin(0.01);
632  c->SetLeftMargin(mode == 1 ? 0.09 : 0.13);
633  switch (mode) {
634  case 1: c->Divide(4,1,0,0); break;
635  case 2: c->Divide(1,4,0,0); break;
636  default: c->Divide(2,2,0,0); break;
637  }
638 
639  DrawStack(c, kPs, Ps, logy, abso, mode);
640  DrawStack(c, kSs, Ss, logy, abso, mode);
641  DrawStack(c, kCs, Cs, logy, abso, mode);
642  DrawStack(c, kCb, Cb, logy, abso, mode);
643 
644  TVirtualPad* p = c->cd(1);
645  TString t = mid ? "|#eta|<1" : "|#eta|>1";
646  if (c1 >= c2) t.Append(" MB");
647  else t.Append(Form(" %4.1f#minus%4.1f%%", c1, c2));
648  TLegend* l = (logy
649  ? p->BuildLegend(0.45,(mode==2?0.3:0.5),0.95,0.9, t)
650  : p->BuildLegend(0.16,0.4,0.6,0.8, t));
651  l->SetFillStyle(0);
652  l->SetBorderSize(0);
653 
654  TFile* output = TFile::Open(Form("%s.root", c->GetName()),"RECREATE");
655  output->cd();
656  THStack* stacks[] = { Ps, Ss, Cs, Cb, 0 };
657  THStack** ptr = stacks;
658  while (*ptr) {
659  THStack* stack = *ptr;
660  TString name = stack->GetTitle(); name.ReplaceAll(" ","");
661  TDirectory* dir = output->mkdir(name);
662  dir->cd();
663  stack->Write("all");
664  TIter next(stack->GetHists());
665  TH1* hist = 0;
666  while ((hist = static_cast<TH1*>(next()))) {
667  name = hist->GetTitle(); name.ReplaceAll("-", "");
668  hist->Write(name);
669  }
670  ptr++;
671  output->cd();
672  }
673  output->Write();
674 
675 
676  c->SaveAs(Form("%s.png",c->GetName()));
677  }
678 };
679 
691 {
692  Bool_t mid=true;
694  cm.Run(c1, c2, mid);
696  cf.Run(c1, c2, !mid);
697 }
698 //
699 // EOF
700 //
Particle * Read(ERow row, EKind kind, TCollection *bin, const TString &sub, const TString &oth, Bool_t abso=true)
static TObject * ChkO(TObject *src, TObject *o, const char *name)
static TObject * GetO(TDirectory *d, const char *name, TClass *cls)
void Print(std::ostream &o, const char *name, Double_t dT, Double_t dVM, Double_t alldT, Double_t alldVM)
Definition: PlotSysInfo.C:121
const Color_t cc[]
Definition: DrawKs.C:1
double Double_t
Definition: External.C:58
const char * title
Definition: MakeQAPdf.C:27
static TH1 * GetH1(TCollection *d, const char *name)
static TCollection * GetC(TCollection *d, const char *name)
const char * parFmt
static TCollection * GetC(TDirectory *d, const char *name)
TCanvas * c
Definition: TestFitELoss.C:172
void Print(Bool_t err=false)
void Scale(EColumn i, Double_t s)
static TObject * GetO(TCollection *d, const char *name, TClass *cls)
void DeltaExpectations(Double_t c1=0, Double_t c2=5)
AliStack * stack
void Run(Double_t c1=0, Double_t c2=0, Bool_t mid=true)
void DrawStack(TVirtualPad *p, THStack *s, Double_t min, Double_t max)
Definition: DrawDeltas3.C:48
static TH2 * GetH2(TCollection *d, const char *name)
void Print(Bool_t err=false) const
ParticleList particles
int Int_t
Definition: External.C:63
Double_t Get(EColumn i) const
Double_t GetE(EColumn i) const
static TObject * ChkC(TObject *o, TClass *c)
Definition: External.C:212
static THStack * GetHS(TCollection *d, const char *name)
Int_t mode
Definition: anaM.C:41
void Print(Bool_t err=false)
Particle * Create(Double_t w, EKind kind, TCollection *reduced, TCollection *reweighed, TCollection *asis, const TString &sub, const TString &oth, Bool_t abso)
Int_t colors[nPtBins]
const char * eRowFmt
const char * rowFmt
THStack * MakeStack(EColumn column, const char *name, const char *title, Int_t mode)
Definition: External.C:220
std::vector< Particle * > ParticleList
const Double_t pi
bool Bool_t
Definition: External.C:53
void DrawStack(TVirtualPad *mother, EColumn column, THStack *stack, Bool_t logy, Bool_t abso, Int_t mode)
THStack * MakeStack(Bool_t isNSD, Bool_t showClusters, Double_t strangeCorr, Double_t maxFactor=1.1)
Definition: ExtractData.C:204
Definition: External.C:196
TDirectoryFile * dir