AliPhysics  f05a842 (f05a842)
 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 #else
29 class TH1D;
30 class TH2D;
31 class TH1;
32 class TH2;
33 class TH3;
34 class TAxis;
35 class TList;
36 class TBrowser;
37 class THStack;
38 class AliAODTracklet;
39 class TCanvas; // Autoload
40 class TVirtualPad;
41 #endif
42 
43 //====================================================================
44 #ifndef __CINT__
45 namespace {
46  template <typename T>
47  T GetPar(const char* name, TCollection* l)
48  {
49  TObject* o = l->FindObject(name);
50  if (!o) {
51  Warning("GetPar", "Didn't find parameter %s in %s",
52  name, l->GetName());
53  return T();
54  }
55  TClass* cls = TParameter<T>::Class();
56  if (!o->IsA()->InheritsFrom(cls)) {
57  Warning("GetPar", "Object %s is a %s, not a %s",
58  name, o->ClassName(), cls->GetName());
59  return T();
60  }
61  TParameter<T>* p = static_cast<TParameter<T>*>(o);
62  return p->GetVal();
63  }
64 }
65 #endif
66 
67 //====================================================================
72 {
73 public:
77  enum ECalc {
120  };
126  : TNamed(),
127  fCalc(kProduct),
128  fMask(0xFF),
129  fVeto(0x0)
130  {}
137  AliTrackletBaseWeights(const char* name,
138  const char* title="Sim. tracklet weights")
139  : TNamed(name,title),
140  fCalc(kProduct),
141  fMask(0xFF),
142  fVeto(0x0),
143  fInverse(false)
144  {}
151  : TNamed(o),
152  fCalc (o.fCalc),
153  fMask (o.fMask),
154  fVeto (o.fVeto),
155  fInverse(o.fInverse)
156  {}
161 
170  {
171  if (&o == this) return *this;
172  TNamed::operator=(o);
173  fCalc = o.fCalc;
174  fMask = o.fMask;
175  fVeto = o.fVeto;
176  fInverse = o.fInverse;
177  return *this;
178  }
226  void SetCalc(UChar_t mode=kProduct) { fCalc = mode; }
232  void SetMask(UChar_t mask) { fMask = mask; }
238  void SetVeto(UChar_t veto) { fVeto = veto; }
246  void SetInverse(Bool_t inv) { fInverse = inv; }
254  Bool_t CheckTracklet(const AliAODTracklet* tracklet) const
255  {
256  UChar_t flags = tracklet->GetFlags();
257  if (fMask != 0xFF && (fMask & flags) == 0) {
258  // Info("LookupWeight", "Tracklet 0x%02x does not fullfill mask 0x%02x",
259  // flags, fMask);
260  return false;
261  }
262  if ((fVeto & flags) != 0) {
263  // Info("LookupWeight", "Tracklet 0x%02x vetoed by 0x%02x",flags, fVeto);
264  return false;
265  }
266  return true;
267  }
279  Double_t cent,
280  Double_t ipz,
281  TH2* corr=0) const
282  {
283  if (!CheckTracklet(tracklet)) return 1;
284  Double_t w = CalcWeight(tracklet, cent, ipz, corr);
285  if (fInverse) w = 1/w;
286  return w;
287  }
299  virtual Double_t CalcWeight(AliAODTracklet* tracklet,
300  Double_t cent,
301  Double_t ipz,
302  TH2* corr) const = 0;
308  virtual TCollection* Store(TCollection* parent)
309  {
310  TList* top = new TList;
311  top->SetName(GetName());
312  top->SetOwner(true);
313  parent->Add(top);
314  top->Add(new TParameter<int>("mask", fMask, 'f'));
315  top->Add(new TParameter<int>("veto", fVeto, 'f'));
316  top->Add(new TParameter<int>("calc", fCalc, 'f'));
317  return top;
318  }
327  {
328  TCollection* top = static_cast<TCollection*>(in->FindObject(GetName()));
329  if (!top) {
330  Warning("Retrieve",
331  "Collection %s not found in %s", GetName(), in->GetName());
332  in->ls();
333  return 0;
334  }
335  fCalc = GetPar<int>("calc", top);
336  fMask = GetPar<int>("mask", top);
337  fVeto = GetPar<int>("veto", top);
338  return top;
339  }
343  Bool_t IsFolder() const { return true; }
349  virtual void Print(Option_t* option="") const //*MENU*
350  {
351  gROOT->IndentLevel();
352  Printf("%s : %s", ClassName(), GetName());
353  Printf(" Weight calculation: %s",
354  fCalc == kProduct ? "product" :
355  fCalc == kSquare ? "square" :
356  fCalc == kSum ? "sum" : "average");
357  Printf(" Tracklet mask: 0x%02x", fMask);
358  Printf(" Tracklet veto: 0x%02x", fVeto);
359  }
360  UChar_t fCalc; // Whether the square of the weight is calculated
361  UChar_t fMask; // Which partiles to take
362  UChar_t fVeto; // Which particles not to take
363  Bool_t fInverse; // If true, do 1/w
364  ClassDef(AliTrackletBaseWeights,1); // Base class of weights
365 };
366 
367 //====================================================================
374 {
375 public:
379  enum EMode {
381  kUp = (0x1) << 14,
383  kDown = (0x2) << 14,
385  kDisabled = (0x4) << 14
386  };
387 
389  typedef std::map<short,TH1D*> PdgMap;
396  fPt(0),
397  fAbundance(),
398  fStrangeness()
399  {}
406  AliTrackletPtPidStrWeights(const char* name,
407  const char* title="Sim. tracklet weights");
426 
437  virtual Double_t CalcWeight(AliAODTracklet* tracklet,
438  Double_t cent,
439  Double_t ipz,
440  TH2* corr=0) const;
451  {
452  AddPdgWeight(fAbundance, pdg, h, mode);
453  }
464  {
465  AddPdgWeight(fStrangeness, pdg, h, mode);
466  }
475  Bool_t SetPtWeight(const TH2D* h, UShort_t mode=0);
489  {
490  SetPdgMode(fAbundance, pdg, mode);
491  }
504  {
505  SetPdgMode(fStrangeness, pdg, mode);
506  }
517  void SetPtMode(UShort_t mode);
523  void Draw(Option_t* option=""); //*MENU*
529  void Print(Option_t* option="") const; //*MENU*
556  {
557  return GetPdgWeight(fAbundance, apdg, cent);
558  }
569  {
570  return GetPdgWeight(fStrangeness, apdg, cent);
571  }
572 
573 protected:
585  virtual Double_t CalcWeight(Double_t pT, Short_t pdg, Double_t cent) const;
594  TH1D* GetPdgHist(const PdgMap& m, Short_t pdg) const;
604  virtual Double_t GetPdgWeight(const PdgMap& m,
605  UShort_t apdg,
606  Double_t cent) const;
617  Bool_t AddPdgWeight(PdgMap& map, Short_t pdg, const TH1D* w, UShort_t mode=0);
629  void SetPdgMode(PdgMap& map, Short_t pdg, UShort_t mode);
637  void StoreMap(TCollection* parent, const char* name, PdgMap& m);
647  Bool_t RetrieveMap(TCollection* parent, const char* name, PdgMap& m);
653  void ModStack(THStack* stack);
661  void PrintMap(const PdgMap& m, const char* name, Option_t* options="") const;
668  void PrintHist(const char* tag, TH1* h) const;
669 
670  TH2D* fPt; // Weight by centrality and pT
671  PdgMap fAbundance; // Map for abundance weight
672  PdgMap fStrangeness; // Map for strangeness weight
673  ClassDef(AliTrackletPtPidStrWeights,3); // Weighs for tracklet analysis
674 };
675 
676 //____________________________________________________________________
678  const char* title)
679  : AliTrackletBaseWeights(name, title),
680  fPt(0),
681  fAbundance(),
682  fStrangeness()
683 {}
684 
685 //____________________________________________________________________
689  fPt(0),
690  fAbundance(),
691  fStrangeness()
692 {
693  UInt_t mask = kUp|kDown|kDisabled;
694  SetPtWeight(o.fPt, o.fPt->TestBits(mask));
695  for (PdgMap::const_iterator i = o.fAbundance.begin();
696  i != o.fAbundance.end(); ++i)
697  AddAbundanceWeight(i->first, i->second, i->second->TestBits(mask));
698  for (PdgMap::const_iterator i = o.fStrangeness.begin();
699  i != o.fStrangeness.end(); ++i)
700  AddStrangenessWeight(i->first, i->second, i->second->TestBits(mask));
701 }
702 //____________________________________________________________________
705 {
706  if (&o == this) return *this;
708  fName = o.fName;
709  fTitle = o.fTitle;
710  UInt_t mask = kUp|kDown|kDisabled;
711  SetPtWeight(o.fPt, o.fPt->TestBits(mask));
712  for (PdgMap::const_iterator i = o.fAbundance.begin();
713  i != o.fAbundance.end(); ++i)
714  AddAbundanceWeight(i->first, i->second, i->second->TestBits(mask));
715  for (PdgMap::const_iterator i = o.fStrangeness.begin();
716  i != o.fStrangeness.end(); ++i)
717  AddStrangenessWeight(i->first, i->second, i->second->TestBits(mask));
718  return *this;
719 }
720 //____________________________________________________________________
722  Short_t pdg,
723  const TH1D* w,
724  UShort_t mode)
725 {
726  UShort_t apdg = TMath::Abs(pdg);
727  TH1D* copy = static_cast<TH1D*>(w->Clone(Form("w%d", apdg)));
728  copy->SetDirectory(0);
729  copy->SetXTitle("Centrality [%]");
730  copy->SetBit(mode);
731  PdgMap::iterator i = m.find(apdg);
732  if (i != m.end()) {
733  Warning("AddPdgWeight", "Replacing weight for %d", pdg);
734  delete i->second;
735  i->second = 0;
736  m.erase(i);
737  }
738  m[apdg] = copy;
739  return true;
740 }
741 
742 //____________________________________________________________________
744 {
745  if (fPt) {
746  delete fPt;
747  fPt = 0;
748  }
749  if (!h) return false;
750  fPt = static_cast<TH2D*>(h->Clone("centPt"));
751  fPt->SetDirectory(0);
752  fPt->SetXTitle("Centrality [%]");
753  fPt->SetYTitle("#it{p}_{T} [GeV]");
754  fPt->SetBit(mode);
755  return true;
756 }
757 //____________________________________________________________________
759 {
760  if (fPt) fPt->SetBit(mode);
761 }
762 //____________________________________________________________________
764  Short_t pdg,
765  UShort_t mode)
766 {
767  TH1* h = GetPdgHist(m, pdg);
768  if (h) h->SetBit(mode);
769 }
770 //____________________________________________________________________
772 {
773  UShort_t apdg = TMath::Abs(pdg);
774  PdgMap::const_iterator i = m.find(apdg);
775  if (i == m.end()) return 0;
776  return i->second;
777 }
778 
779 //____________________________________________________________________
781  UShort_t apdg,
782  Double_t cent) const
783 {
784  if (m.size() < 1) return 1;
785  TH1D* h = GetPdgHist(m, apdg);
786  if (!h || h->TestBit(kDisabled)) return 1;
787 
788 #if 1
789  Int_t bC = h->GetXaxis()->FindBin(cent);
790  if (bC < 1 || bC > h->GetNbinsX()) return 1;
791 #else
792  Int_t bC = 1;
793 #endif
794  Double_t add = (h->TestBit(kUp) ? +1 :
795  h->TestBit(kDown) ? -1 : 0) * h->GetBinError(bC);
796  return h->GetBinContent(bC) + add;
797 }
798 //____________________________________________________________________
799 Double_t
801  Short_t pdg,
802  Double_t cent) const
803 {
804  Double_t w = 1;
805  if (fPt && !fPt->TestBit(kDisabled)) {
806  Int_t bC = fPt->GetXaxis()->FindBin(cent);
807  Int_t bpT = fPt->GetYaxis()->FindBin(pT);
808  if (bC >= 1 && bC <= fPt->GetXaxis()->GetNbins() &&
809  bpT >= 1 && bpT <= fPt->GetYaxis()->GetNbins()) {
810  Double_t fac = (pT >= 0.05 ? 1 :
811  (fPt->TestBit(kUp) ? 1.3 :
812  fPt->TestBit(kDown) ? 0.7 : 1));
813  w *= fac*fPt->GetBinContent(bC, bpT);
814  }
815  }
816  UShort_t apdg = TMath::Abs(pdg);
817 
818  w *= GetPdgWeight(fAbundance, apdg, cent);
819  w *= GetPdgWeight(fStrangeness, apdg, cent);
820  // Printf("Weight of pT=%6.3f pdg=%5d cent=%5.1f -> %f", pT, pdg, cent, w);
821  return w;
822 }
823 
824 //____________________________________________________________________
825 Double_t
827  Double_t cent,
828  Double_t ipz,
829  TH2* corr) const
830 {
831 #if 0
832  if (!tracklet->IsSimulated()) {
833  Warning("LookupWeight", "Not a simulated tracklet");
834  return 1;
835  }
836 #endif
837  Double_t w1 = 1, w2 = 1;
838 
839  // AliAODMCTracklet* mc = static_cast<AliAODMCTracklet*>(tracklet);
840  AliAODTracklet* mc = tracklet;
841  Short_t pdg1 = mc->GetParentPdg();
842  Short_t pdg2 = mc->GetParentPdg(true);
843  if (pdg1>0) w1 = CalcWeight(mc->GetParentPt(), pdg1,cent);
844  if (pdg2>0) w2 = CalcWeight(mc->GetParentPt(true),pdg2,cent);
845  else if (fCalc!=kProduct) w2 = w1;
846 
847  if (corr && mc->IsMeasured()) corr->Fill(w1, w2);
848 
849  switch (fCalc) {
850  case kProduct: return w1 * w2;
851  case kSquare: return TMath::Sqrt(w1 * w2);
852  case kSum: return 1+(w1-1)+(w2-1);
853  case kAverage: return 1+((w1-1)+(w2-1))/2;
854  }
855  return 1;
856 }
857 
858 //____________________________________________________________________
860 {
861  if (!stack || !stack->GetHists()) return;
862  TIter next(stack->GetHists());
863  TH1* hist = 0;
864  Color_t colors[] = { kPink+2, kRed+2, kOrange+2, kYellow+2,
865  kSpring+2, kGreen+2, kTeal+2, kCyan+2,
866  kBlue+2, kViolet+2 };
867  Int_t idx = 0;
868  while ((hist = static_cast<TH1*>(next()))) {
869  Color_t c = colors[idx % 10];
870  idx++;
871  hist->SetLineColor(c);
872  hist->SetMarkerColor(c);
873  hist->SetMarkerStyle(20+idx%20);
874  }
875 }
876 
877 //____________________________________________________________________
880 {
882  if (!top) return 0;
883 
884  if (fPt) {
885  TH2* copy = static_cast<TH2*>(fPt->Clone("centPt"));
886  copy->SetDirectory(0);
887  copy->SetBinContent(0,0,1); // For counting merges
888  top->Add(copy);
889  }
890  StoreMap(top, "abundance", fAbundance);
891  StoreMap(top, "strangeness", fStrangeness);
892  return top;
893 }
894 
895 //____________________________________________________________________
896 void
898  const char* name,
899  PdgMap& m)
900 {
901  TList* top = new TList;
902  top->SetName(name);
903  top->SetOwner(true);
904  parent->Add(top);
905  for (PdgMap::const_iterator i = m.begin(); i != m.end(); ++i) {
906  TH1* copy = static_cast<TH1*>(i->second->Clone(Form("w%d",i->first)));
907  copy->SetDirectory(0);
908  copy->SetBinContent(0,1); // For counting merges
909  top->Add(copy);
910  }
911 }
912 
913 
914 //____________________________________________________________________
917 {
919  if (!top) return 0;
920 
921  fPt = static_cast<TH2D*>(top->FindObject("centPt"));
922  if (!fPt)
923  Warning("Retrieve","centPt histogram not found in %s", GetName());
924  else {
925  fPt->SetDirectory(0);
926  Double_t scale = fPt->GetBinContent(0,0);
927  fPt->Scale(1/scale); // Counting merges
928  // fPt->SetBinContent(0,0,1); // Zero merger count
929  }
930  if (!RetrieveMap(top, "abundance", fAbundance)) return 0;
931  if (!RetrieveMap(top, "strangeness", fStrangeness)) return 0;
932 
933  return top;
934 }
935 //____________________________________________________________________
936 Bool_t
938  const char* name,
939  PdgMap& m)
940 {
941  m.clear();
942  TList* top = static_cast<TList*>(parent->FindObject(name));
943  if (!top) {
944  Warning("RetrieveMap",
945  "Collection %s not found in %s", name, parent->GetName());
946  return true;
947  }
948  TIter next(top);
949  TObject* o = 0;
950  while ((o = next())) {
951  if (!o->IsA()->InheritsFrom(TH1D::Class())) continue;
952  TH1D* copy = static_cast<TH1D*>(o);
953  copy->SetDirectory(0);
954  Double_t scale = copy->GetBinContent(0);
955  copy->Scale(1/scale); // Counting merges
956  // copy->SetBinContent(0,1); // Zero merger count
957  TString nme(copy->GetName());
958  nme.Remove(0,1);
959  Int_t apdg = nme.Atoi();
960  AddPdgWeight(m, apdg, copy, copy->TestBits(kUp|kDown|kDisabled));
961  }
962  return true;
963 }
964 
965 //____________________________________________________________________
966 void
968 {
969  TVirtualPad* master = TVirtualPad::Pad();
970  if (!master) {
971  Warning("Draw", "No current pad to draw in");
972  return;
973  }
974 
975  Int_t nPad = 3;
976  Int_t iPad = 0;
977  if (fAbundance.size() <= 0) nPad--;
978  if (fStrangeness.size() <= 0) nPad--;
979  TString opt(option); opt.ToUpper();
980 
981  master->Divide(nPad,1);
982  if (fPt) {
983  THStack* stack = new THStack(fPt, (opt.Contains("C") ? "x" : "y"));
984  ModStack(stack);
985  stack->SetTitle("#it{p}_{T} weights");
986  master->cd(++iPad);
987  stack->Draw("nostack");
988  stack->GetHistogram()->SetXTitle((opt.Contains("C") ?
989  "Centrality [%]" : "#it{p}_{T}"));
990  master->GetPad(iPad)->BuildLegend();
991  master->GetPad(iPad)->Modified();
992  }
993 
994  if (fAbundance.size() > 0) {
995  THStack* ha = new THStack("abundance", "Abundance weights");
996  for (PdgMap::const_iterator i = fAbundance.begin();
997  i != fAbundance.end(); ++i) {
998  ha->Add(i->second);
999  }
1000  ModStack(ha);
1001  master->cd(++iPad);
1002  ha->Draw("nostack");
1003  ha->GetHistogram()->SetXTitle("Centrality [%]");
1004  master->GetPad(iPad)->BuildLegend();
1005  }
1006 
1007  if (fStrangeness.size() > 0) {
1008  THStack* hs = new THStack("strangeness", "Strangeness weights");
1009  for (PdgMap::const_iterator i = fStrangeness.begin();
1010  i != fStrangeness.end(); ++i) {
1011  hs->Add(i->second);
1012  }
1013  ModStack(hs);
1014  master->cd(++iPad);
1015  hs->Draw("nostack");
1016  hs->GetHistogram()->SetXTitle("Centrality [%]");
1017  master->GetPad(iPad)->BuildLegend();
1018  }
1019  master->Modified();
1020 }
1021 
1022 //____________________________________________________________________
1023 void
1025 {
1027  gROOT->IncreaseDirLevel();
1028  PrintHist("pT", fPt);
1029  PrintMap(fAbundance, "Abundance", option);
1030  PrintMap(fStrangeness, "Strangeness", option);
1031 
1032  gROOT->DecreaseDirLevel();
1033 }
1034 //____________________________________________________________________
1035 void
1036 AliTrackletPtPidStrWeights::PrintHist(const char* tag, TH1* h) const
1037 {
1038  if (!h) return;
1039  gROOT->IndentLevel();
1040  Printf("%10s (%c): %p %s/%s", tag,
1041  h->TestBit(kDisabled) ? '0' :
1042  h->TestBit(kUp) ? '+' :
1043  h->TestBit(kDown) ? '-' : '=',
1044  h, h->GetName(), h->GetTitle());
1045 
1046 }
1047 //____________________________________________________________________
1048 void
1050  const char* name,
1051  Option_t* option) const
1052 {
1053  gROOT->IndentLevel();
1054  Printf("Map of PDG codes: %s", name);
1055  gROOT->IncreaseDirLevel();
1056  for (PdgMap::const_iterator i = m.begin(); i != m.end(); ++i)
1057  PrintHist(Form("%10d", i->first), i->second);
1058 
1059  gROOT->DecreaseDirLevel();
1060 }
1061 
1062 
1063 //====================================================================
1068 {
1069 public:
1076  fHistos(),
1077  fCentAxis()
1078  {}
1085  AliTrackletDeltaWeights(const char* name,
1086  const char* title="Sim. tracklet weights");
1097 
1111  void SetCentAxis(const TAxis& axis);
1118  void SetCentAxis(Int_t n, const Double_t* bins);
1126  void SetCentAxis(Int_t n, Double_t low, Double_t high);
1135  Bool_t SetHisto(Int_t bin, TH3* h);
1143  TH3* FindHisto(Double_t cent) const;
1154  static Int_t FindBin(Double_t value, const TAxis* axis);
1165  virtual Double_t CalcWeight(AliAODTracklet* tracklet,
1166  Double_t cent,
1167  Double_t ipZ,
1168  TH2* corr=0) const;
1169 
1179  TH1* Project(TH3* h, char which, const char* name);
1185  void Draw(Option_t* option=""); //*MENU*
1191  void Print(Option_t* option="") const; //*MENU*
1201  void DrawStack(TVirtualPad* pad,
1202  THStack* stack,
1203  const char* xtitle,
1204  Double_t min=0,
1205  Double_t max=2.5);
1206  void DrawOne(Double_t cent=2.5,
1207  Double_t eta=0,
1208  Double_t ipz=0) const; //*MENU*
1213  ClassDef(AliTrackletDeltaWeights,1); // Base class of weights
1214 };
1215 
1216 //____________________________________________________________________
1218  const char* title)
1219  : AliTrackletBaseWeights(name, title),
1220  fHistos(),
1221  fCentAxis()
1222 {
1223  fHistos.SetOwner(true);
1224 }
1225 //____________________________________________________________________
1228  fHistos(),
1229  fCentAxis(o.fCentAxis)
1230 {
1231  fHistos.SetOwner(true);
1232  for (Int_t i = 0; i < o.fHistos.GetEntriesFast(); i++) {
1233  TH3* h = static_cast<TH3*>(o.fHistos.At(i));
1234  if (!h) continue;
1235  h = static_cast<TH3*>(h->Clone());
1236  h->SetDirectory(0);
1237  fHistos.Add(h);
1238  }
1239 }
1240 //____________________________________________________________________
1243 {
1244  if (&o == this) return *this;
1246  fCentAxis = o.fCentAxis;
1247  for (Int_t i = 0; i < o.fHistos.GetEntriesFast(); i++) {
1248  TH3* h = static_cast<TH3*>(o.fHistos.At(i));
1249  if (!h) continue;
1250  h = static_cast<TH3*>(h->Clone());
1251  h->SetDirectory(0);
1252  fHistos.Add(h);
1253  }
1254  return *this;
1255 }
1256 
1257 //____________________________________________________________________
1258 void
1260 {
1261  if (a.GetNbins() && a.GetXbins()->GetArray())
1262  SetCentAxis(a.GetNbins(), a.GetXbins()->GetArray());
1263  else
1264  SetCentAxis(a.GetNbins(), a.GetXmin(), a.GetXmax());
1265 }
1266 
1267 //____________________________________________________________________
1268 void
1270 {
1271  fCentAxis.Set(n, bins);
1272 }
1273 //____________________________________________________________________
1274 void
1276 {
1277  fCentAxis.Set(n, l, h);
1278 }
1279 
1280 //____________________________________________________________________
1281 Bool_t
1283 {
1284  if (bin < 1 || bin > fCentAxis.GetNbins()) {
1285  Warning("SetHisto", "Centrality bin %d out of range [%d,%d]",
1286  bin, 1, fCentAxis.GetNbins());
1287  return false;
1288  }
1289  Double_t c1 = fCentAxis.GetBinLowEdge(bin);
1290  Double_t c2 = fCentAxis.GetBinUpEdge(bin);
1291  TString name;
1292  name.Form("cent%03dd%02d_%03dd%02d",
1293  Int_t(c1), Int_t(c1 * 100) % 100,
1294  Int_t(c2), Int_t(c2 * 100) % 100);
1295  TH3* cpy = static_cast<TH3*>(h->Clone(name));
1296  cpy->SetDirectory(0);
1297  cpy->SetTitle(Form("%6.2f%% - %6.2f%%", c1, c2));
1298  fHistos.AddAtAndExpand(cpy, bin-1);
1299  return true;
1300 }
1301 
1302 //____________________________________________________________________
1303 TH3*
1305 {
1306  Int_t bin = FindBin(cent, &fCentAxis);
1307  if (bin >= fHistos.GetEntriesFast()) return 0;
1308  return static_cast<TH3*>(fHistos.At(bin-1));
1309 }
1310 
1311 //____________________________________________________________________
1312 Int_t
1314  const TAxis* axis)
1315 {
1316  Int_t bin = const_cast<TAxis*>(axis)->FindBin(value);
1317  if (bin < 1) bin = 1;
1318  else if (bin > axis->GetNbins()) bin = axis->GetNbins();
1319  return bin;
1320 }
1321 
1322 
1323 //____________________________________________________________________
1324 Double_t
1326  Double_t cent,
1327  Double_t ipZ,
1328  TH2* corr) const
1329 {
1330  TH3* h = FindHisto(cent);
1331  if (!h) return 1;
1332 
1333  Double_t eta = tracklet->GetEta();
1334  Double_t delta = tracklet->GetDelta();
1335  Int_t etaBin = FindBin(eta, h->GetXaxis());
1336  Int_t deltaBin = FindBin(delta, h->GetYaxis());
1337  Int_t ipzBin = FindBin(ipZ, h->GetZaxis());
1338 
1339  return h->GetBinContent(etaBin, deltaBin, ipzBin);
1340 }
1341 
1342 //____________________________________________________________________
1343 TH1*
1344 AliTrackletDeltaWeights::Project(TH3* h, char which, const char* name)
1345 {
1346  TAxis* a0 = 0;
1347  TAxis* a1 = 0;
1348  TAxis* a2 = 0;
1349  typedef
1350  TH1D* (TH3::*ProjFunc)(const char*,Int_t,Int_t,Int_t,Int_t,Option_t*) const;
1351  ProjFunc func = 0;
1352  switch (which) {
1353  case 'x':
1354  a0 = h->GetXaxis();
1355  a1 = h->GetYaxis();
1356  a2 = h->GetZaxis();
1357  func = &TH3::ProjectionX;
1358  break;
1359  case 'y':
1360  a0 = h->GetYaxis();
1361  a1 = h->GetXaxis();
1362  a2 = h->GetZaxis();
1363  func = &TH3::ProjectionY;
1364  break;
1365  case 'z':
1366  a0 = h->GetZaxis();
1367  a1 = h->GetXaxis();
1368  a2 = h->GetYaxis();
1369  func = &TH3::ProjectionZ;
1370  break;
1371  }
1372 #if 0
1373  TH1* ret = (h->*func)(Form("%s_%s",h->GetName(),name),
1374  1, a1->GetNbins(), 1, a2->GetNbins(),"");
1375 #endif
1376  TString n; n.Form("%s_%s",h->GetName(),name);
1377  TString t(h->GetTitle());
1378  Int_t nx = a0->GetNbins();
1379  TH1* ret = 0;
1380  if (a0->GetXbins() && a0->GetXbins()->GetArray())
1381  ret = new TH1D(n,t,nx,a0->GetXbins()->GetArray());
1382  else
1383  ret = new TH1D(n,t,nx,a0->GetXmin(),a0->GetXmax());
1384  ret->Sumw2();
1385  ret->SetXTitle(a0->GetTitle());
1386  static_cast<const TAttAxis*>(a0)->Copy(*(ret->GetXaxis()));
1387  ret->SetDirectory(0);
1388  ret->SetLineColor(h->GetMarkerColor());
1389  ret->SetMarkerColor(h->GetMarkerColor());
1390  ret->SetFillColor(kWhite);// color);
1391  ret->SetFillStyle(0);
1392  ret->SetMarkerStyle(20);
1393 
1394  for (Int_t i0 = 1; i0 <= ret->GetNbinsX(); i0++) {
1395  Int_t count = 0;
1396  Double_t sum = 0;
1397  Double_t sumE2 = 0;
1398  for (Int_t i1 = 1; i1 <= a1->GetNbins(); i1++) {
1399  for (Int_t i2 = 1; i2 <= a2->GetNbins(); i2++) {
1400  Int_t bin = 0;
1401  switch (which) {
1402  case 'x': bin = h->GetBin(i0, i1, i1); break;
1403  case 'y': bin = h->GetBin(i1, i0, i2); break;
1404  case 'z': bin = h->GetBin(i1, i2, i0); break;
1405  }
1406  Double_t c = h->GetBinContent(bin);
1407  Double_t e = h->GetBinError (bin);
1408  if (c < 1e-3 || e < 1e-6) continue;
1409  count++;
1410  sum += c;
1411  sumE2 += e*e;
1412  }
1413  }
1414  // ret->SetBinContent(i0, ret->GetBinContent(i0)/count);
1415  // ret->SetBinError (i0, ret->GetBinError (i0)/count);
1416  ret->SetBinContent(i0, sum/count);
1417  ret->SetBinError (i0, TMath::Sqrt(sumE2)/count);
1418  }
1419 
1420 
1421  // ret->Scale(1./(a1->GetNbins()*a2->GetNbins()));
1422  return ret;
1423 }
1424 
1425 //____________________________________________________________________
1426 void
1428 {
1430  gROOT->IncreaseDirLevel();
1431  gROOT->IndentLevel();
1432  Printf("%d centrality bins", fCentAxis.GetNbins());
1433  for (Int_t i = 0; i < fHistos.GetEntriesFast(); i++) {
1434  TH3* h = static_cast<TH3*>(fHistos.At(i));
1435  if (!h) continue;
1436  h->Print(option);
1437  }
1438  gROOT->DecreaseDirLevel();
1439 }
1440 //____________________________________________________________________
1441 void
1443  Double_t eta,
1444  Double_t ipz) const
1445 {
1446  TH3* h3 = FindHisto(cent);
1447  if (!h3) return;
1448 
1449  Int_t etaBin = FindBin(eta, h3->GetXaxis());
1450  Int_t ipzBin = FindBin(ipz, h3->GetZaxis());
1451 
1452  TH1* h1 = h3->ProjectionY("tmp", etaBin, etaBin, ipzBin, ipzBin);
1453  h1->SetDirectory(0);
1454  h1->Draw();
1455 }
1456 
1457 //____________________________________________________________________
1458 void
1460  THStack* stack,
1461  const char* xtitle,
1462  Double_t min,
1463  Double_t max)
1464 {
1465  stack->SetMinimum(min);
1466  stack->SetMaximum(max);
1467  pad->SetTicks();
1468  pad->SetGridy();
1469  pad->cd();
1470  stack->Draw("nostack");
1471  TH1* h = stack->GetHistogram();
1472  h->SetXTitle(xtitle);
1473  h->SetYTitle(stack->GetTitle());
1474  h->GetXaxis()->SetNdivisions(205);
1475  h->GetYaxis()->SetNdivisions(205);
1476  pad->Modified();
1477 }
1478 
1479 //____________________________________________________________________
1480 void
1482 {
1483  TVirtualPad* master = TVirtualPad::Pad();
1484  if (!master) {
1485  Warning("Draw", "No current pad to draw in");
1486  return;
1487  }
1488 
1489  Int_t nPad = 3;
1490  Int_t iPad = 0;
1491  TString opt(option); opt.ToUpper();
1492 
1493  master->SetTopMargin(0.01);
1494  master->SetRightMargin(0.01);
1495  if (opt.Contains("V"))
1496  master->Divide(1,nPad);
1497  else
1498  master->Divide(nPad,1, 0, 0);
1499 
1500  THStack* stackEta = new THStack("eta", "#LTw#GT_{#Delta,IP_{z}}");
1501  THStack* stackDelta = new THStack("delta", "#LTw#GT_{#eta,IP_{z}}");
1502  THStack* stackIPz = new THStack("ipz", "#LTw#GT_{#eta,#Delta}");
1503 
1504  for (Int_t i = 0; i < fHistos.GetEntriesFast(); i++) {
1505  TH3* h = static_cast<TH3*>(fHistos.At(i));
1506  if (!h) continue;
1507 
1508  TH1* etaProj = Project(h, 'x', "eta");
1509  TH1* deltaProj = Project(h, 'y', "delta");
1510  TH1* ipzProj = Project(h, 'z', "ipz");
1511  stackEta ->Add(etaProj);
1512  stackDelta->Add(deltaProj);
1513  stackIPz ->Add(ipzProj);
1514  }
1515 
1516  DrawStack(master->cd(1), stackEta, "#eta");
1517  DrawStack(master->cd(2), stackDelta, "#Delta");
1518  DrawStack(master->cd(3), stackIPz, "IP_{z}");
1519 
1520  master->Modified();
1521 }
1522 
1523 #endif
1524 //____________________________________________________________________
1525 //
1526 // EOF
1527 //
1528 
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)
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:26
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 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
AliTrackletDeltaWeights & operator=(const AliTrackletDeltaWeights &o)
ClassDef(AliTrackletPtPidStrWeights, 3)
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:40
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
TH3 * FindHisto(Double_t cent) 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
Double_t LookupWeight(AliAODTracklet *tracklet, Double_t cent, Double_t ipz, TH2 *corr=0) const
UChar_t GetFlags() 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
TH1 * Project(TH3 *h, char which, const char *name)
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
bool Bool_t
Definition: External.C:53
Double_t GetAbundanceWeight(UShort_t apdg, Double_t cent) const
AliTrackletPtPidStrWeights & operator=(const AliTrackletPtPidStrWeights &o)
Bool_t CheckTracklet(const AliAODTracklet *tracklet) const
void DrawStack(TVirtualPad *pad, THStack *stack, const char *xtitle, Double_t min=0, Double_t max=2.5)
Definition: External.C:196
ClassDef(AliTrackletBaseWeights, 1)
Bool_t SetHisto(Int_t bin, TH3 *h)