AliPhysics  ebc57ae (ebc57ae)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliTrackletWeights.C
Go to the documentation of this file.
1 
11 #ifndef ALITRACKLETWEIGHTS_C
12 #define ALITRACKLETWEIGHTS_C
13 #include <TNamed.h>
14 #include <map>
15 #ifndef __CINT__
16 # include <TH1.h>
17 # include <TH2.h>
18 # include <TH3.h>
19 # include <TAxis.h>
20 # include <TList.h>
21 # include <TBrowser.h>
22 # include <TCanvas.h>
23 # include <THStack.h>
24 # include <TROOT.h>
25 # include <TClass.h>
26 # include "AliAODTracklet.C"
27 # include <TParameter.h>
28 # include <TParticle.h>
29 #else
30 class TH1D;
31 class TH2D;
32 class TH1;
33 class TH2;
34 class TH3;
35 class TAxis;
36 class TList;
37 class TBrowser;
38 class THStack;
39 class AliAODTracklet;
40 class TCanvas; // Autoload
41 class TVirtualPad;
42 class TParticle;
43 #endif
44 
45 //====================================================================
46 #ifndef __CINT__
47 namespace {
48  template <typename T>
49  T GetPar(const char* name, TCollection* l)
50  {
51  TObject* o = l->FindObject(name);
52  if (!o) {
53  Warning("GetPar", "Didn't find parameter %s in %s",
54  name, l->GetName());
55  return T();
56  }
57  TClass* cls = TParameter<T>::Class();
58  if (!o->IsA()->InheritsFrom(cls)) {
59  Warning("GetPar", "Object %s is a %s, not a %s",
60  name, o->ClassName(), cls->GetName());
61  return T();
62  }
63  TParameter<T>* p = static_cast<TParameter<T>*>(o);
64  return p->GetVal();
65  }
66 }
67 #endif
68 
69 //====================================================================
74 {
75 public:
79  enum ECalc {
122  };
128  : TNamed(),
129  fCalc(kProduct),
130  fMask(0xFF),
131  fVeto(0x0),
132  fInverse(false),
133  fDebug(0)
134  {}
141  AliTrackletBaseWeights(const char* name,
142  const char* title="Sim. tracklet weights")
143  : TNamed(name,title),
144  fCalc(kProduct),
145  fMask(0xFF),
146  fVeto(0x0),
147  fInverse(false),
148  fDebug(0)
149  {}
156  : TNamed(o),
157  fCalc (o.fCalc),
158  fMask (o.fMask),
159  fVeto (o.fVeto),
160  fInverse(o.fInverse),
161  fDebug (o.fDebug)
162  {}
167 
176  {
177  if (&o == this) return *this;
178  TNamed::operator=(o);
179  fCalc = o.fCalc;
180  fMask = o.fMask;
181  fVeto = o.fVeto;
182  fInverse = o.fInverse;
183  fDebug = o.fDebug;
184  return *this;
185  }
233  void SetCalc(UChar_t mode=kProduct) { fCalc = mode; }
239  void SetMask(UChar_t mask) { fMask = mask; }
245  void SetVeto(UChar_t veto) { fVeto = veto; }
253  void SetInverse(Bool_t inv) { fInverse = inv; }
259  void SetDebug(Int_t lvl) { fDebug = lvl; }
267  Bool_t CheckTracklet(const AliAODTracklet* tracklet) const
268  {
269  UChar_t flags = tracklet->GetFlags();
270  if (fMask != 0xFF && (fMask & flags) == 0) {
271  if (fDebug > 3)
272  Info("LookupWeight", "Tracklet 0x%02x does not fullfill mask 0x%02x",
273  flags, fMask);
274  return false;
275  }
276  if ((fVeto & flags) != 0) {
277  if (fDebug > 3)
278  Info("LookupWeight", "Tracklet 0x%02x vetoed by 0x%02x",flags, fVeto);
279  return false;
280  }
281  return true;
282  }
294  Double_t cent,
295  Double_t ipz,
296  TH2* corr=0) const
297  {
298  if (!CheckTracklet(tracklet)) return 1;
299  Double_t w = CalcWeight(tracklet, cent, ipz, corr);
300  if (fInverse) w = 1/w;
301  return w;
302  }
312  Double_t LookupWeight(TParticle* particle,
313  Double_t cent,
314  Double_t ipz) const
315  {
316  Double_t w = CalcWeight(particle, cent, ipz);
317  if (fInverse) w = 1/w;
318  return w;
319  }
331  virtual Double_t CalcWeight(AliAODTracklet* tracklet,
332  Double_t cent,
333  Double_t ipz,
334  TH2* corr) const = 0;
344  virtual Double_t CalcWeight(TParticle* particle,
345  Double_t cent,
346  Double_t ipZ) const = 0;
352  virtual TCollection* Store(TCollection* parent)
353  {
354  TList* top = new TList;
355  top->SetName(GetName());
356  top->SetOwner(true);
357  parent->Add(top);
358  top->Add(new TParameter<int> ("mask", fMask, 'f'));
359  top->Add(new TParameter<int> ("veto", fVeto, 'f'));
360  top->Add(new TParameter<int> ("calc", fCalc, 'f'));
361  top->Add(new TParameter<bool>("inv", fInverse, 'f'));
362  return top;
363  }
372  {
373  TCollection* top = static_cast<TCollection*>(in->FindObject(GetName()));
374  if (!top) {
375  Warning("Retrieve",
376  "Collection %s not found in %s", GetName(), in->GetName());
377  in->ls();
378  return 0;
379  }
380  fCalc = GetPar<int>("calc", top);
381  fMask = GetPar<int>("mask", top);
382  fVeto = GetPar<int>("veto", top);
383  return top;
384  }
388  Bool_t IsFolder() const { return true; }
394  virtual void Print(Option_t* option="") const //*MENU*
395  {
396  gROOT->IndentLevel();
397  Printf("%s : %s", ClassName(), GetName());
398  Printf(" Weight calculation: %s",
399  fCalc == kProduct ? "product" :
400  fCalc == kSquare ? "square" :
401  fCalc == kSum ? "sum" : "average");
402  Printf(" Tracklet mask: 0x%02x", fMask);
403  Printf(" Tracklet veto: 0x%02x", fVeto);
404  Printf(" Take inverse: %s", fInverse ? "yes" : "no");
405  }
406  UChar_t fCalc; // Whether the square of the weight is calculated
407  UChar_t fMask; // Which partiles to take
408  UChar_t fVeto; // Which particles not to take
409  Bool_t fInverse; // If true, do 1/w
410  Int_t fDebug; // Debug level
411  ClassDef(AliTrackletBaseWeights,1); // Base class of weights
412 };
413 
414 //====================================================================
421 {
422 public:
426  enum EMode {
428  kUp = (0x1) << 14,
430  kDown = (0x2) << 14,
432  kDisabled = (0x4) << 14
433  };
434 
436  typedef std::map<short,TH1D*> PdgMap;
443  fPt(0),
444  fAbundance(),
445  fStrangeness()
446  {}
453  AliTrackletPtPidStrWeights(const char* name,
454  const char* title="Sim. tracklet weights");
473 
484  virtual Double_t CalcWeight(AliAODTracklet* tracklet,
485  Double_t cent,
486  Double_t ipz,
487  TH2* corr=0) const;
488  virtual Double_t CalcWeight(TParticle* particle,
489  Double_t cent,
490  Double_t ipZ) const;
501  {
502  AddPdgWeight(fAbundance, pdg, h, mode);
503  }
514  {
515  AddPdgWeight(fStrangeness, pdg, h, mode);
516  }
525  Bool_t SetPtWeight(const TH2D* h, UShort_t mode=0);
539  {
540  SetPdgMode(fAbundance, pdg, mode);
541  }
554  {
555  SetPdgMode(fStrangeness, pdg, mode);
556  }
567  void SetPtMode(UShort_t mode);
573  void Draw(Option_t* option=""); //*MENU*
579  void Print(Option_t* option="") const; //*MENU*
606  {
607  return GetPdgWeight(fAbundance, apdg, cent);
608  }
619  {
620  return GetPdgWeight(fStrangeness, apdg, cent);
621  }
622 
623 protected:
635  virtual Double_t CalcWeight(Double_t pT, Short_t pdg, Double_t cent) const;
644  TH1D* GetPdgHist(const PdgMap& m, Short_t pdg) const;
654  virtual Double_t GetPdgWeight(const PdgMap& m,
655  UShort_t apdg,
656  Double_t cent) const;
667  Bool_t AddPdgWeight(PdgMap& map, Short_t pdg, const TH1D* w, UShort_t mode=0);
679  void SetPdgMode(PdgMap& map, Short_t pdg, UShort_t mode);
687  void StoreMap(TCollection* parent, const char* name, PdgMap& m);
697  Bool_t RetrieveMap(TCollection* parent, const char* name, PdgMap& m);
703  void ModStack(THStack* stack);
711  void PrintMap(const PdgMap& m, const char* name, Option_t* options="") const;
718  void PrintHist(const char* tag, TH1* h) const;
719 
720  TH2D* fPt; // Weight by centrality and pT
721  PdgMap fAbundance; // Map for abundance weight
722  PdgMap fStrangeness; // Map for strangeness weight
723  ClassDef(AliTrackletPtPidStrWeights,3); // Weighs for tracklet analysis
724 };
725 
726 //____________________________________________________________________
728  const char* title)
729  : AliTrackletBaseWeights(name, title),
730  fPt(0),
731  fAbundance(),
732  fStrangeness()
733 {}
734 
735 //____________________________________________________________________
739  fPt(0),
740  fAbundance(),
741  fStrangeness()
742 {
743  UInt_t mask = kUp|kDown|kDisabled;
744  SetPtWeight(o.fPt, o.fPt->TestBits(mask));
745  for (PdgMap::const_iterator i = o.fAbundance.begin();
746  i != o.fAbundance.end(); ++i)
747  AddAbundanceWeight(i->first, i->second, i->second->TestBits(mask));
748  for (PdgMap::const_iterator i = o.fStrangeness.begin();
749  i != o.fStrangeness.end(); ++i)
750  AddStrangenessWeight(i->first, i->second, i->second->TestBits(mask));
751 }
752 //____________________________________________________________________
755 {
756  if (&o == this) return *this;
758  fName = o.fName;
759  fTitle = o.fTitle;
760  UInt_t mask = kUp|kDown|kDisabled;
761  SetPtWeight(o.fPt, o.fPt->TestBits(mask));
762  for (PdgMap::const_iterator i = o.fAbundance.begin();
763  i != o.fAbundance.end(); ++i)
764  AddAbundanceWeight(i->first, i->second, i->second->TestBits(mask));
765  for (PdgMap::const_iterator i = o.fStrangeness.begin();
766  i != o.fStrangeness.end(); ++i)
767  AddStrangenessWeight(i->first, i->second, i->second->TestBits(mask));
768  return *this;
769 }
770 //____________________________________________________________________
772  Short_t pdg,
773  const TH1D* w,
774  UShort_t mode)
775 {
776  UShort_t apdg = TMath::Abs(pdg);
777  TH1D* copy = static_cast<TH1D*>(w->Clone(Form("w%d", apdg)));
778  copy->SetDirectory(0);
779  copy->SetXTitle("Centrality [%]");
780  copy->SetBit(mode);
781  PdgMap::iterator i = m.find(apdg);
782  if (i != m.end()) {
783  Warning("AddPdgWeight", "Replacing weight for %d", pdg);
784  delete i->second;
785  i->second = 0;
786  m.erase(i);
787  }
788  m[apdg] = copy;
789  return true;
790 }
791 
792 //____________________________________________________________________
794 {
795  if (fPt) {
796  delete fPt;
797  fPt = 0;
798  }
799  if (!h) return false;
800  fPt = static_cast<TH2D*>(h->Clone("centPt"));
801  fPt->SetDirectory(0);
802  fPt->SetXTitle("Centrality [%]");
803  fPt->SetYTitle("#it{p}_{T} [GeV]");
804  fPt->SetBit(mode);
805  return true;
806 }
807 //____________________________________________________________________
809 {
810  if (fPt) fPt->SetBit(mode);
811 }
812 //____________________________________________________________________
814  Short_t pdg,
815  UShort_t mode)
816 {
817  TH1* h = GetPdgHist(m, pdg);
818  if (h) h->SetBit(mode);
819 }
820 //____________________________________________________________________
822 {
823  UShort_t apdg = TMath::Abs(pdg);
824  PdgMap::const_iterator i = m.find(apdg);
825  if (i == m.end()) return 0;
826  return i->second;
827 }
828 
829 //____________________________________________________________________
831  UShort_t apdg,
832  Double_t cent) const
833 {
834  if (m.size() < 1) return 1;
835  TH1D* h = GetPdgHist(m, apdg);
836  if (!h || h->TestBit(kDisabled)) return 1;
837 
838 #if 1
839  Int_t bC = h->GetXaxis()->FindBin(cent);
840  if (bC < 1 || bC > h->GetNbinsX()) return 1;
841 #else
842  Int_t bC = 1;
843 #endif
844  Double_t add = (h->TestBit(kUp) ? +1 :
845  h->TestBit(kDown) ? -1 : 0) * h->GetBinError(bC);
846  return h->GetBinContent(bC) + add;
847 }
848 //____________________________________________________________________
849 Double_t
851  Short_t pdg,
852  Double_t cent) const
853 {
854  Double_t w = 1;
855  if (fPt && !fPt->TestBit(kDisabled)) {
856  Int_t bC = fPt->GetXaxis()->FindBin(cent);
857  Int_t bpT = fPt->GetYaxis()->FindBin(pT);
858  if (bC >= 1 && bC <= fPt->GetXaxis()->GetNbins() &&
859  bpT >= 1 && bpT <= fPt->GetYaxis()->GetNbins()) {
860  Double_t fac = (pT >= 0.05 ? 1 :
861  (fPt->TestBit(kUp) ? 1.3 :
862  fPt->TestBit(kDown) ? 0.7 : 1));
863  Double_t ptW = fPt->GetBinContent(bC, bpT);
864  w *= fac*ptW;
865  if (fDebug > 3)
866  Info("CalcWeight", "pT=%5.2f -> %4.2f * %6.4f = %6.4f",
867  pT, fac, ptW, fac*ptW);
868  }
869  }
870  UShort_t apdg = TMath::Abs(pdg);
871  Double_t aW = GetPdgWeight(fAbundance, apdg, cent);
872  Double_t sW = GetPdgWeight(fStrangeness, apdg, cent);
873  w *= aW * sW;
874  if (fDebug > 3)
875  Info("CalcWeight","pdg=%4d -> %6.4f * %6.4f = %6.4f -> %6.4f",
876  pdg, aW, sW, aW*sW, w);
877  // Printf("Weight of pT=%6.3f pdg=%5d cent=%5.1f -> %f", pT, pdg, cent, w);
878  return w;
879 }
880 
881 //____________________________________________________________________
882 Double_t
884  Double_t cent,
885  Double_t ipz) const
886 {
887  Int_t pdg = particle->GetPdgCode();
888  Double_t pT = particle->Pt();
889  return CalcWeight(pT, pdg, cent);
890 }
891 
892 //____________________________________________________________________
893 Double_t
895  Double_t cent,
896  Double_t ipz,
897  TH2* corr) const
898 {
899 #if 0
900  if (!tracklet->IsSimulated()) {
901  Warning("CalcWeight", "Not a simulated tracklet");
902  return 1;
903  }
904 #endif
905  Double_t w1 = 1, w2 = 1;
906 
907  // AliAODMCTracklet* mc = static_cast<AliAODMCTracklet*>(tracklet);
908  AliAODTracklet* mc = tracklet;
909  Short_t pdg1 = mc->GetParentPdg();
910  Short_t pdg2 = mc->GetParentPdg(true);
911  if (pdg1!=0) w1 = CalcWeight(mc->GetParentPt(), pdg1,cent);
912  if (pdg2!=0) w2 = CalcWeight(mc->GetParentPt(true),pdg2,cent);
913  else if (fCalc!=kProduct) w2 = w1;
914 
915  if (corr && mc->IsMeasured()) corr->Fill(w1, w2);
916 
917  Double_t ret = 1;
918  const char* cm = "?";
919  switch (fCalc) {
920  case kProduct: ret = w1 * w2; cm="product"; break;
921  case kSquare: ret = TMath::Sqrt(w1 * w2); cm="square"; break;
922  case kSum: ret = 1+(w1-1)+(w2-1); cm="sum"; break;
923  case kAverage: ret = 1+((w1-1)+(w2-1))/2; cm="average"; break;
924  }
925 
926  if (fDebug > 1)
927  Printf("pdg1=%5d -> %6.4f pdg2=%5d -> %6.4f (%10s) -> %6.4f",
928  pdg1, w1, pdg2, w2, cm, ret);
929  return ret;
930 }
931 
932 //____________________________________________________________________
934 {
935  if (!stack || !stack->GetHists()) return;
936  TIter next(stack->GetHists());
937  TH1* hist = 0;
938  Color_t colors[] = { kPink+2, kRed+2, kOrange+2, kYellow+2,
939  kSpring+2, kGreen+2, kTeal+2, kCyan+2,
940  kBlue+2, kViolet+2 };
941  Int_t idx = 0;
942  while ((hist = static_cast<TH1*>(next()))) {
943  Color_t c = colors[idx % 10];
944  idx++;
945  hist->SetLineColor(c);
946  hist->SetMarkerColor(c);
947  hist->SetMarkerStyle(20+idx%20);
948  }
949 }
950 
951 //____________________________________________________________________
954 {
956  if (!top) return 0;
957 
958  if (fPt) {
959  TH2* copy = static_cast<TH2*>(fPt->Clone("centPt"));
960  copy->SetDirectory(0);
961  copy->SetBinContent(0,0,1); // For counting merges
962  top->Add(copy);
963  }
964  StoreMap(top, "abundance", fAbundance);
965  StoreMap(top, "strangeness", fStrangeness);
966  return top;
967 }
968 
969 //____________________________________________________________________
970 void
972  const char* name,
973  PdgMap& m)
974 {
975  TList* top = new TList;
976  top->SetName(name);
977  top->SetOwner(true);
978  parent->Add(top);
979  for (PdgMap::const_iterator i = m.begin(); i != m.end(); ++i) {
980  TH1* copy = static_cast<TH1*>(i->second->Clone(Form("w%d",i->first)));
981  copy->SetDirectory(0);
982  copy->SetBinContent(0,1); // For counting merges
983  top->Add(copy);
984  }
985 }
986 
987 
988 //____________________________________________________________________
991 {
993  if (!top) return 0;
994 
995  fPt = static_cast<TH2D*>(top->FindObject("centPt"));
996  if (!fPt)
997  Warning("Retrieve","centPt histogram not found in %s", GetName());
998  else {
999  fPt->SetDirectory(0);
1000  Double_t scale = fPt->GetBinContent(0,0);
1001  fPt->Scale(1/scale); // Counting merges
1002  // fPt->SetBinContent(0,0,1); // Zero merger count
1003  }
1004  if (!RetrieveMap(top, "abundance", fAbundance)) return 0;
1005  if (!RetrieveMap(top, "strangeness", fStrangeness)) return 0;
1006 
1007  return top;
1008 }
1009 //____________________________________________________________________
1010 Bool_t
1012  const char* name,
1013  PdgMap& m)
1014 {
1015  m.clear();
1016  TList* top = static_cast<TList*>(parent->FindObject(name));
1017  if (!top) {
1018  Warning("RetrieveMap",
1019  "Collection %s not found in %s", name, parent->GetName());
1020  return true;
1021  }
1022  TIter next(top);
1023  TObject* o = 0;
1024  while ((o = next())) {
1025  if (!o->IsA()->InheritsFrom(TH1D::Class())) continue;
1026  TH1D* copy = static_cast<TH1D*>(o);
1027  copy->SetDirectory(0);
1028  Double_t scale = copy->GetBinContent(0);
1029  copy->Scale(1/scale); // Counting merges
1030  // copy->SetBinContent(0,1); // Zero merger count
1031  TString nme(copy->GetName());
1032  nme.Remove(0,1);
1033  Int_t apdg = nme.Atoi();
1034  AddPdgWeight(m, apdg, copy, copy->TestBits(kUp|kDown|kDisabled));
1035  }
1036  return true;
1037 }
1038 
1039 //____________________________________________________________________
1040 void
1042 {
1043  TVirtualPad* master = TVirtualPad::Pad();
1044  if (!master) {
1045  Warning("Draw", "No current pad to draw in");
1046  return;
1047  }
1048 
1049  Int_t nPad = 3;
1050  Int_t iPad = 0;
1051  if (fAbundance.size() <= 0) nPad--;
1052  if (fStrangeness.size() <= 0) nPad--;
1053  TString opt(option); opt.ToUpper();
1054 
1055  master->Divide(nPad,1);
1056  if (fPt) {
1057  THStack* stack = new THStack(fPt, (opt.Contains("C") ? "x" : "y"));
1058  ModStack(stack);
1059  stack->SetTitle("#it{p}_{T} weights");
1060  master->cd(++iPad);
1061  stack->Draw("nostack");
1062  stack->GetHistogram()->SetXTitle((opt.Contains("C") ?
1063  "Centrality [%]" : "#it{p}_{T}"));
1064  master->GetPad(iPad)->BuildLegend();
1065  master->GetPad(iPad)->Modified();
1066  }
1067 
1068  if (fAbundance.size() > 0) {
1069  THStack* ha = new THStack("abundance", "Abundance weights");
1070  for (PdgMap::const_iterator i = fAbundance.begin();
1071  i != fAbundance.end(); ++i) {
1072  ha->Add(i->second);
1073  }
1074  ModStack(ha);
1075  master->cd(++iPad);
1076  ha->Draw("nostack");
1077  ha->GetHistogram()->SetXTitle("Centrality [%]");
1078  master->GetPad(iPad)->BuildLegend();
1079  }
1080 
1081  if (fStrangeness.size() > 0) {
1082  THStack* hs = new THStack("strangeness", "Strangeness weights");
1083  for (PdgMap::const_iterator i = fStrangeness.begin();
1084  i != fStrangeness.end(); ++i) {
1085  hs->Add(i->second);
1086  }
1087  ModStack(hs);
1088  master->cd(++iPad);
1089  hs->Draw("nostack");
1090  hs->GetHistogram()->SetXTitle("Centrality [%]");
1091  master->GetPad(iPad)->BuildLegend();
1092  }
1093  master->Modified();
1094 }
1095 
1096 //____________________________________________________________________
1097 void
1099 {
1101  gROOT->IncreaseDirLevel();
1102  PrintHist("pT", fPt);
1103  PrintMap(fAbundance, "Abundance", option);
1104  PrintMap(fStrangeness, "Strangeness", option);
1105 
1106  gROOT->DecreaseDirLevel();
1107 }
1108 //____________________________________________________________________
1109 void
1110 AliTrackletPtPidStrWeights::PrintHist(const char* tag, TH1* h) const
1111 {
1112  if (!h) return;
1113  gROOT->IndentLevel();
1114  Printf("%10s (%c): %p %s/%s", tag,
1115  h->TestBit(kDisabled) ? '0' :
1116  h->TestBit(kUp) ? '+' :
1117  h->TestBit(kDown) ? '-' : '=',
1118  h, h->GetName(), h->GetTitle());
1119 
1120 }
1121 //____________________________________________________________________
1122 void
1124  const char* name,
1125  Option_t* option) const
1126 {
1127  gROOT->IndentLevel();
1128  Printf("Map of PDG codes: %s", name);
1129  gROOT->IncreaseDirLevel();
1130  for (PdgMap::const_iterator i = m.begin(); i != m.end(); ++i)
1131  PrintHist(Form("%10d", i->first), i->second);
1132 
1133  gROOT->DecreaseDirLevel();
1134 }
1135 
1136 
1137 //====================================================================
1142 {
1143 public:
1150  fHistos(),
1151  fCentAxis()
1152  {}
1159  AliTrackletDeltaWeights(const char* name,
1160  const char* title="Sim. tracklet weights");
1171 
1185  void SetCentAxis(const TAxis& axis);
1192  void SetCentAxis(Int_t n, const Double_t* bins);
1200  void SetCentAxis(Int_t n, Double_t low, Double_t high);
1209  Bool_t SetHisto(Int_t bin, TH1* h);
1217  TH1* FindHisto(Double_t cent) const;
1228  static Int_t FindBin(Double_t value, const TAxis* axis);
1239  virtual Double_t CalcWeight(AliAODTracklet* tracklet,
1240  Double_t cent,
1241  Double_t ipZ,
1242  TH2* corr=0) const;
1243  virtual Double_t CalcWeight(TParticle* particle,
1244  Double_t cent,
1245  Double_t ipZ) const;
1255  TH1* Project(TH1* h, char which, const char* name);
1261  void Draw(Option_t* option=""); //*MENU*
1267  void Print(Option_t* option="") const; //*MENU*
1278  void DrawStack(TVirtualPad* pad,
1279  Int_t nPad,
1280  THStack* stack,
1281  const char* xtitle,
1282  Double_t min=0,
1283  Double_t max=2.5);
1284  void DrawOne(Double_t cent=2.5,
1285  Double_t eta=0,
1286  Double_t ipz=0) const; //*MENU*
1291  ClassDef(AliTrackletDeltaWeights,1); // Base class of weights
1292 };
1293 
1294 //____________________________________________________________________
1296  const char* title)
1297  : AliTrackletBaseWeights(name, title),
1298  fHistos(),
1299  fCentAxis()
1300 {
1301  fHistos.SetOwner(true);
1302 }
1303 //____________________________________________________________________
1306  fHistos(),
1307  fCentAxis(o.fCentAxis)
1308 {
1309  fHistos.SetOwner(true);
1310  for (Int_t i = 0; i < o.fHistos.GetEntriesFast(); i++) {
1311  TH1* h = static_cast<TH1*>(o.fHistos.At(i));
1312  if (!h) continue;
1313  h = static_cast<TH1*>(h->Clone());
1314  h->SetDirectory(0);
1315  fHistos.Add(h);
1316  }
1317 }
1318 //____________________________________________________________________
1321 {
1322  if (&o == this) return *this;
1324  fCentAxis = o.fCentAxis;
1325  for (Int_t i = 0; i < o.fHistos.GetEntriesFast(); i++) {
1326  TH1* h = static_cast<TH1*>(o.fHistos.At(i));
1327  if (!h) continue;
1328  h = static_cast<TH1*>(h->Clone());
1329  h->SetDirectory(0);
1330  fHistos.Add(h);
1331  }
1332  return *this;
1333 }
1334 
1335 //____________________________________________________________________
1336 void
1338 {
1339  if (a.GetNbins() && a.GetXbins()->GetArray())
1340  SetCentAxis(a.GetNbins(), a.GetXbins()->GetArray());
1341  else
1342  SetCentAxis(a.GetNbins(), a.GetXmin(), a.GetXmax());
1343 }
1344 
1345 //____________________________________________________________________
1346 void
1348 {
1349  fCentAxis.Set(n, bins);
1350 }
1351 //____________________________________________________________________
1352 void
1354 {
1355  fCentAxis.Set(n, l, h);
1356 }
1357 
1358 //____________________________________________________________________
1359 Bool_t
1361 {
1362  if (bin < 1 || bin > fCentAxis.GetNbins()) {
1363  Warning("SetHisto", "Centrality bin %d out of range [%d,%d]",
1364  bin, 1, fCentAxis.GetNbins());
1365  return false;
1366  }
1367  Double_t c1 = fCentAxis.GetBinLowEdge(bin);
1368  Double_t c2 = fCentAxis.GetBinUpEdge(bin);
1369  TString name;
1370  name.Form("cent%03dd%02d_%03dd%02d",
1371  Int_t(c1), Int_t(c1 * 100) % 100,
1372  Int_t(c2), Int_t(c2 * 100) % 100);
1373  TH1* cpy = static_cast<TH1*>(h->Clone(name));
1374  cpy->SetDirectory(0);
1375  cpy->SetTitle(Form("%6.2f%% - %6.2f%%", c1, c2));
1376  fHistos.AddAtAndExpand(cpy, bin-1);
1377  return true;
1378 }
1379 
1380 //____________________________________________________________________
1381 TH1*
1383 {
1384  Int_t bin = FindBin(cent, &fCentAxis);
1385  if (bin >= fHistos.GetEntriesFast()) return 0;
1386  return static_cast<TH1*>(fHistos.At(bin-1));
1387 }
1388 
1389 //____________________________________________________________________
1390 Int_t
1392  const TAxis* axis)
1393 {
1394  if (!axis) return 1;
1395  Int_t bin = const_cast<TAxis*>(axis)->FindBin(value);
1396  if (bin < 1) bin = 1;
1397  else if (bin > axis->GetNbins()) bin = axis->GetNbins();
1398  return bin;
1399 }
1400 
1401 //____________________________________________________________________
1402 Double_t
1404  Double_t cent,
1405  Double_t ipZ) const
1406 {
1407  // Perhaps this method doesn't really make much sense
1408  TH1* h = FindHisto(cent);
1409  if (!h) return 1;
1410 
1411  Double_t eta = particle->Eta();
1412  Double_t delta = 0;
1413  Int_t etaBin = FindBin(eta, h->GetXaxis());
1414  Int_t deltaBin = FindBin(delta, h->GetYaxis());
1415  Int_t ipzBin = FindBin(ipZ, h->GetZaxis());
1416 
1417  return h->GetBinContent(etaBin, deltaBin, ipzBin);
1418 }
1419 
1420 //____________________________________________________________________
1421 Double_t
1423  Double_t cent,
1424  Double_t ipZ,
1425  TH2* corr) const
1426 {
1427  if (tracklet->IsGenerated()) return 1;
1428  TH1* h = FindHisto(cent);
1429  if (!h) return 1;
1430 
1431  Double_t eta = tracklet->GetEta();
1432  Double_t delta = tracklet->GetDelta();
1433  Int_t etaBin = FindBin(eta, h->GetXaxis());
1434  Int_t deltaBin = FindBin(delta, h->GetYaxis());
1435  Int_t ipzBin = FindBin(ipZ, h->GetZaxis());
1436 
1437  return h->GetBinContent(etaBin, deltaBin, ipzBin);
1438 }
1439 
1440 //____________________________________________________________________
1441 TH1*
1442 AliTrackletDeltaWeights::Project(TH1* h, char which, const char* name)
1443 {
1444  TString n; n.Form("%s_%s",h->GetName(),name);
1445  TString t(h->GetTitle());
1446  TH1* ret = 0;
1447  TAxis* a0 = 0;
1448  TAxis* a1 = 0;
1449  TAxis* a2 = 0;
1450  typedef
1451  TH1D* (TH3::*ProjFunc)(const char*,Int_t,Int_t,Int_t,Int_t,Option_t*) const;
1452  ProjFunc func = 0;
1453  switch (which) {
1454  case 'x':
1455  a0 = h->GetXaxis();
1456  a1 = h->GetYaxis();
1457  a2 = h->GetZaxis();
1458  func = &TH3::ProjectionX;
1459  break;
1460  case 'y':
1461  a0 = h->GetYaxis();
1462  a1 = h->GetXaxis();
1463  a2 = h->GetZaxis();
1464  func = &TH3::ProjectionY;
1465  break;
1466  case 'z':
1467  a0 = h->GetZaxis();
1468  a1 = h->GetXaxis();
1469  a2 = h->GetYaxis();
1470  func = &TH3::ProjectionZ;
1471  break;
1472  }
1473  Int_t nx = a0->GetNbins();
1474  if (nx <= 1) return 0;
1475  if (a0->GetXbins() && a0->GetXbins()->GetArray())
1476  ret = new TH1D(n,t,nx,a0->GetXbins()->GetArray());
1477  else
1478  ret = new TH1D(n,t,nx,a0->GetXmin(),a0->GetXmax());
1479  ret->Sumw2();
1480  ret->SetXTitle(a0->GetTitle());
1481  static_cast<const TAttAxis*>(a0)->Copy(*(ret->GetXaxis()));
1482  ret->SetDirectory(0);
1483  ret->SetLineColor(h->GetMarkerColor());
1484  ret->SetMarkerColor(h->GetMarkerColor());
1485  ret->SetFillColor(kWhite);// color);
1486  ret->SetFillStyle(0);
1487  ret->SetMarkerStyle(20);
1488 
1489  for (Int_t i0 = 1; i0 <= ret->GetNbinsX(); i0++) {
1490  Int_t count = 0;
1491  Double_t sum = 0;
1492  Double_t sumE2 = 0;
1493  for (Int_t i1 = 1; i1 <= a1->GetNbins(); i1++) {
1494  for (Int_t i2 = 1; i2 <= a2->GetNbins(); i2++) {
1495  Int_t bin = 0;
1496  switch (which) {
1497  case 'x': bin = h->GetBin(i0, i1, i1); break;
1498  case 'y': bin = h->GetBin(i1, i0, i2); break;
1499  case 'z': bin = h->GetBin(i1, i2, i0); break;
1500  }
1501  Double_t c = h->GetBinContent(bin);
1502  Double_t e = h->GetBinError (bin);
1503  if (c < 1e-3 || e < 1e-6) continue;
1504  count++;
1505  sum += c;
1506  sumE2 += e*e;
1507  }
1508  }
1509  // ret->SetBinContent(i0, ret->GetBinContent(i0)/count);
1510  // ret->SetBinError (i0, ret->GetBinError (i0)/count);
1511  ret->SetBinContent(i0, sum/count);
1512  ret->SetBinError (i0, TMath::Sqrt(sumE2)/count);
1513  }
1514 
1515 
1516  // ret->Scale(1./(a1->GetNbins()*a2->GetNbins()));
1517  return ret;
1518 }
1519 
1520 //____________________________________________________________________
1521 void
1523 {
1525  gROOT->IncreaseDirLevel();
1526  gROOT->IndentLevel();
1527  Printf("%d centrality bins", fCentAxis.GetNbins());
1528  for (Int_t i = 0; i < fHistos.GetEntriesFast(); i++) {
1529  TH1* h = static_cast<TH1*>(fHistos.At(i));
1530  if (!h) continue;
1531  h->Print(option);
1532  }
1533  gROOT->DecreaseDirLevel();
1534 }
1535 //____________________________________________________________________
1536 void
1538  Double_t eta,
1539  Double_t ipz) const
1540 {
1541  TH1* h = FindHisto(cent);
1542  if (!h) return;
1543 
1544  Int_t etaBin = FindBin(eta, h->GetXaxis());
1545  Int_t ipzBin = FindBin(ipz, h->GetZaxis());
1546 
1547  TH1* tmp = 0;
1548  if (h->IsA()->InheritsFrom(TH2::Class()))
1549  tmp = static_cast<TH2*>(h)->ProjectionY("tmp", etaBin, etaBin);
1550  if (h->IsA()->InheritsFrom(TH3::Class()))
1551  tmp = static_cast<TH3*>(h)->ProjectionY("tmp", etaBin, etaBin);
1552  else
1553  tmp = static_cast<TH1*>(h->Clone("tmp"));
1554  if (!tmp) return;
1555  tmp->SetDirectory(0);
1556  tmp->Draw();
1557 }
1558 
1559 //____________________________________________________________________
1560 void
1562  Int_t nPad,
1563  THStack* stack,
1564  const char* xtitle,
1565  Double_t min,
1566  Double_t max)
1567 {
1568  if (!stack->GetHists()) return;
1569  stack->SetMinimum(min);
1570  stack->SetMaximum(max);
1571  pad->SetTicks();
1572  pad->SetGridy();
1573  pad->cd();
1574  stack->Draw("nostack");
1575  TH1* h = stack->GetHistogram();
1576  h->SetXTitle(xtitle);
1577  h->SetYTitle(stack->GetTitle());
1578  h->GetXaxis()->SetNdivisions(205);
1579  h->GetYaxis()->SetNdivisions(205);
1580  h->GetXaxis()->SetLabelOffset(0.005*pad->GetWNDC());
1581  h->GetXaxis()->SetTitleOffset(1.5*pad->GetWNDC());
1582  h->GetXaxis()->SetLabelSize(0.03/pad->GetWNDC()/*nPad*/);
1583  h->GetYaxis()->SetLabelSize(0.03/pad->GetWNDC()/*nPad*/);
1584  h->GetXaxis()->SetTitleSize(0.03/pad->GetWNDC()/*nPad*/);
1585  h->GetYaxis()->SetTitleSize(0.03/pad->GetWNDC()/*nPad*/);
1586  pad->Modified();
1587 }
1588 
1589 //____________________________________________________________________
1590 void
1592 {
1593  TVirtualPad* master = TVirtualPad::Pad();
1594  if (!master) {
1595  Warning("Draw", "No current pad to draw in");
1596  return;
1597  }
1598 
1599  TString opt(option); opt.ToUpper();
1600  THStack* stackEta = new THStack("eta", "#LTw#GT_{#Delta,IP_{z}}");
1601  THStack* stackDelta = new THStack("delta", "#LTw#GT_{#eta,IP_{z}}");
1602  THStack* stackIPz = new THStack("ipz", "#LTw#GT_{#eta,#Delta}");
1603 
1604  for (Int_t i = 0; i < fHistos.GetEntriesFast(); i++) {
1605  TH3* h = static_cast<TH3*>(fHistos.At(i));
1606  if (!h) continue;
1607 
1608  TH1* etaProj = Project(h, 'x', "eta");
1609  TH1* deltaProj = Project(h, 'y', "delta");
1610  TH1* ipzProj = Project(h, 'z', "ipz");
1611  stackEta ->Add(etaProj);
1612  stackDelta->Add(deltaProj);
1613  stackIPz ->Add(ipzProj);
1614  }
1615 
1616  Int_t nPad = 0;
1617  if (stackEta ->GetHists()) nPad++;
1618  if (stackDelta->GetHists()) nPad++;
1619  if (stackIPz ->GetHists()) nPad++;
1620 
1621  master->SetTopMargin(0.01);
1622  master->SetRightMargin(0.01);
1623  master->SetLeftMargin(0.07*nPad);
1624  if (opt.Contains("V"))
1625  master->Divide(1,nPad);
1626  else
1627  master->Divide(nPad,1, 0, 0);
1628 
1629  if (nPad >= 1) DrawStack(master->cd(1), nPad, stackEta, "#eta");
1630  if (nPad >= 2) DrawStack(master->cd(2), nPad, stackDelta, "#Delta");
1631  if (nPad >= 3) DrawStack(master->cd(3), nPad, stackIPz, "IP_{z}");
1632  master->GetPad(nPad)->SetRightMargin(0.01);
1633 
1634  master->Modified();
1635 }
1636 
1637 #endif
1638 //____________________________________________________________________
1639 //
1640 // EOF
1641 //
1642 
ClassDef(AliTrackletDeltaWeights, 1)
virtual Double_t CalcWeight(AliAODTracklet *tracklet, Double_t cent, Double_t ipz, TH2 *corr=0) const
Int_t pdg
void SetPtMode(UShort_t mode)
Double_t LookupWeight(TParticle *particle, Double_t cent, Double_t ipz) const
static Int_t FindBin(Double_t value, const TAxis *axis)
Real_t GetEta() const
double Double_t
Definition: External.C:58
Bool_t AddAbundanceWeight(Short_t pdg, const TH1D *h, UShort_t mode=0)
TCollection * Store(TCollection *out)
const char * title
Definition: MakeQAPdf.C:27
virtual TCollection * Store(TCollection *parent)
AliTrackletBaseWeights & operator=(const AliTrackletBaseWeights &o)
void Draw(Option_t *option="")
Definition: External.C:244
AliTrackletBaseWeights(const char *name, const char *title="Sim. tracklet weights")
void SetCalc(UChar_t mode=kProduct)
Bool_t AddStrangenessWeight(Short_t pdg, const TH1D *h, UShort_t mode=0)
void SetStrangenessMode(Short_t pdg, UShort_t mode)
void StoreMap(TCollection *parent, const char *name, PdgMap &m)
Bool_t SetHisto(Int_t bin, TH1 *h)
Bool_t RetrieveMap(TCollection *parent, const char *name, PdgMap &m)
void SetInverse(Bool_t inv)
TCanvas * c
Definition: TestFitELoss.C:172
virtual Double_t CalcWeight(AliAODTracklet *tracklet, Double_t cent, Double_t ipz, TH2 *corr) const =0
void DrawStack(TVirtualPad *pad, Int_t nPad, THStack *stack, const char *xtitle, Double_t min=0, Double_t max=2.5)
AliTrackletDeltaWeights & operator=(const AliTrackletDeltaWeights &o)
ClassDef(AliTrackletPtPidStrWeights, 3)
AliStack * stack
void Print(Option_t *option="") const
void Draw(Option_t *option="")
virtual Real_t GetParentPt(Bool_t second=false) const
void SetMask(UChar_t mask)
UShort_t T(UShort_t m, UShort_t t)
Definition: RingBits.C:60
virtual TCollection * Retrieve(TCollection *in)
Bool_t IsSimulated() const
void ModStack(THStack *stack)
Tracklet AOD object.
void SetCentAxis(const TAxis &axis)
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
TCollection * Retrieve(TCollection *in)
AliTrackletBaseWeights(const AliTrackletBaseWeights &o)
Definition: External.C:228
Definition: External.C:212
Bool_t IsMeasured() const
virtual Short_t GetParentPdg(Bool_t second=false) const
virtual void Print(Option_t *option="") const
virtual Double_t CalcWeight(AliAODTracklet *tracklet, Double_t cent, Double_t ipZ, TH2 *corr=0) const
Int_t mode
Definition: anaM.C:41
Bool_t SetPtWeight(const TH2D *h, UShort_t mode=0)
Double_t GetStrangenessWeight(UShort_t apdg, Double_t cent) const
TH1D * GetPdgHist(const PdgMap &m, Short_t pdg) const
void SetVeto(UChar_t veto)
Bool_t AddPdgWeight(PdgMap &map, Short_t pdg, const TH1D *w, UShort_t mode=0)
short Short_t
Definition: External.C:23
Int_t colors[nPtBins]
Double_t LookupWeight(AliAODTracklet *tracklet, Double_t cent, Double_t ipz, TH2 *corr=0) const
UChar_t GetFlags() const
Bool_t IsGenerated() const
void SetPdgMode(PdgMap &map, Short_t pdg, UShort_t mode)
std::map< short, TH1D * > PdgMap
void DrawOne(Double_t cent=2.5, Double_t eta=0, Double_t ipz=0) const
Definition: External.C:220
void PrintHist(const char *tag, TH1 *h) const
void SetAbundanceMode(Short_t pdg, UShort_t mode)
void Print(Option_t *option="") const
virtual Double_t GetPdgWeight(const PdgMap &m, UShort_t apdg, Double_t cent) const
unsigned short UShort_t
Definition: External.C:28
const char Option_t
Definition: External.C:48
void PrintMap(const PdgMap &m, const char *name, Option_t *options="") const
Real_t GetDelta() const
TList * GetHists(UShort_t flags, const char *var, const char *stackName="result", const char *sub="")
Definition: Compare.C:81
bool Bool_t
Definition: External.C:53
TH1 * FindHisto(Double_t cent) const
Double_t GetAbundanceWeight(UShort_t apdg, Double_t cent) const
AliTrackletPtPidStrWeights & operator=(const AliTrackletPtPidStrWeights &o)
Bool_t CheckTracklet(const AliAODTracklet *tracklet) const
TH1 * Project(TH1 *h, char which, const char *name)
Definition: External.C:196
ClassDef(AliTrackletBaseWeights, 1)